It has been discussed that Intel MKL can exhibit irreproducible behavior under certain conditions. In fact, this is a known thing and is described by Intel as Conditional Numerical Reproducibility.A colleague has recently encountered this issue while trying to find the eigenvalues of a matrix and he realized that the function eigen in R is not deterministic on the server he was using. In particular, eigenvalues can be slightly different (<1e-12) and the signs of the eigenvectors are non-deterministic.I tested his hypothesis by considering this matrix corr <- structure(c(1, 0.250050163743823, 0.23435347961023, 0.482039807584937,                 0.260244618444847, 0.245304368749426, 0.486808023309575, 0.479348415857738,                 0.250050163743823, 1, 0.24648658148082, 0.246713558484002, 0.249547953702824,                 0.249424786521267, 0.247921981042718, 0.493767069248003, 0.23435347961023,                 0.24648658148082, 1, 0.0140300592972217, 0.229687680789224, 0.251843633087952,                 0.0102624595100619, 0.0160021143836412, 0.482039807584937, 0.246713558484002,                 0.0140300592972217, 1, 0.47661703967983, 0.00364166259778339,                 0.257242621615393, 0.507588614471663, 0.260244618444847, 0.249547953702824,                 0.229687680789224, 0.47661703967983, 1, 0.243407621930236, 0.482342591096063,                 0.473425275160343, 0.245304368749426, 0.249424786521267, 0.251843633087952,                 0.00364166259778339, 0.243407621930236, 1, 0.00266373582273809,                 0.00415353194451873, 0.486808023309575, 0.247921981042718, 0.0102624595100619,                 0.257242621615393, 0.482342591096063, 0.00266373582273809, 1,                 0.505031653312059, 0.479348415857738, 0.493767069248003, 0.0160021143836412,                 0.507588614471663, 0.473425275160343, 0.00415353194451873, 0.505031653312059,                 1), .Dim = structure(c(8L, 8L), .Names = c(rt, rt)), .Dimnames = structure(list(                    rt = c(A, B, C, D, E, F, G, H                    ), rt = c(A, B, C, D, E, F, G,                               H)), .Names = c(rt, rt)))after computing its eigenvectors and eigenvalues I saw that they are indeed not the same:> percDet <- 100 * mean(sapply(1:1000, function(...) {+     identical(eigen(corr, symmetric = TRUE),+               eigen(corr, symmetric = TRUE))}))> message(% eigen determinism: , percDet, %)% eigen determinism: 75.7%> percDet <- 100 * mean(sapply(1:1000, function(...) {+     identical(eigen(corr, symmetric = TRUE)$values,+               eigen(corr, symmetric = TRUE)$values)}))> message(% eigenvalues determinism: , percDet, %)% eigenvalues determinism: 76.3%After some investigation I've found out that the issue was due to the Intel MKL configuration we used for R. Specifically, in our Ansible role we were configuring R as follows:./configure --enable-R-shlib --with-blas=\-Wl,--no-as-needed -L${MKLROOT}/lib/intel64 -L{{ mkl_install_dir }}/compiler/lib/intel64 -lmkl_gf_lp64 -lmkl_core -lmkl_intel_thread -liomp5 -lpthread -lm\Changing -lmkl_gf_lp64 part to -lmkl_intel_lp64 as below ./configure --enable-R-shlib --with-blas=-Wl,--no-as-needed -L${MKLROOT}/lib/intel64 -L{{ mkl_install_dir }}/compiler/lib/intel64 -lmkl_intel_lp64 -lmkl_core -lmkl_intel_thread -liomp5 -lpthread\solved the issue:> percDet <- 100 * mean(sapply(1:1000, function(...) {+ identical(eigen(corr, symmetric = TRUE),+ eigen(corr, symmetric = TRUE))}))> message(% eigen determinism: , percDet, %)% eigen determinism: 100%> percDet <- 100 * mean(sapply(1:1000, function(...) {+ identical(eigen(corr, symmetric = TRUE)$values,+ eigen(corr, symmetric = TRUE)$values)}))> message(% eigenvalues determinism: , percDet, %)% eigenvalues determinism: 100%While I have found a related question Intel MKL - Difference between mkl_intel_lp64 and mkl_gf_lp64 helpful, it doesn't explain why this leads to a reducibility issue when multiple threads are used. So, I'm wondering whether there is a good explanation to this behavior and why changing mkl_gf_lp64 flag to mkl_intel_lp64 resolves the issue.
The difference between mkl_intel_lp64 vs mkl_gf_lp64 in a numerical reproducibility issue with Intel MKL
blas;r;compiling;intel mkl
The various Fortran standards allow a lot of compiler dependent behaviour in terms of function binary interfaces when being called with complicated data types such as Fortran90 style arrays and complex numbers. This means calling code compiled with one compiler from another one is not guaranteed to do what you expect, and can lead to grabbing the wrong bit of memory for a variable, thus giving spurious results.You're using the intel compiler suite to build R, so you need to to link to the intel build of the fortran library for your BLAS to behave correctly. This is what the -lmkl_intel_lp64, embedded inside an instruction to the linker, does (as opposed to your previous version, which was linking against a gfortran version, `libmkl_gf_lp64). Since threading makes it more likely that the wrong memory you grabbed was different between different runs, this showed up as non-deterministic behaviour, but there's no prior reason to thing that a serial run linking against the wrong library was generating the right result either.