OpenMP + MKL does not give matching result between LU inverse and QR inverse

Summary

With OMP_NUM_THREADS and MKL_NUM_THREADS both greater than 1, the LU inverse and QR inverse do not match each other for certain matrix. However, when any of the two is turned to 1, then the inverses match.

Environment

  • Operating System : Linux
  • Architecture : x64/Arm64/PowerPC ...
  • Eigen Version : 3.4.0
  • Compiler Version : Gcc9.3.0
  • Compile Flags : none
  • Vector Extension : none

Minimal Example

I attached my testing files, including

  • matrix.csv: a matrix generated from my application
  • eigen_run.cpp: a script reading the matrix and compute its inverse using different methods
  • flare.cpp: an almost empty file (with only an empty main function that does nothing)
  • CMakeLists.txt: it creates a static library flare using flare.cpp (linked to MKL and OMP), and then connect eigen_run to flare.

Note: In the CMakeLists.txt, the lines of linking OMP (shown below) are commented. And you can uncomment it to reproduce the error:

# Link to OpenMP.
find_package(OpenMP)
if(OpenMP_CXX_FOUND)
    target_link_libraries(flare PUBLIC OpenMP::OpenMP_CXX)
endif()

my_eigen.tar.gz

Steps to reproduce

  1. Make sure OpenMP and MKL are installed and can be found. (Or you can modify the CMakeLists correspondingly)
  2. Uncomment the "Link to OpenMP" block in CMakeLists.txt
  3. Run following
cd my_eigen
mkdir build
cd build
cmake ..
make -j
cp ../matrix.csv .
export OMP_NUM_THREADS=<set to 1 or 2>
export MKL_NUM_THREADS=<set to 1 or 2>
./eigen_run

What is the current bug behavior?

When OMP_NUM_THREADS=2 and MKL_NUM_THREADS=2, the A.inverse() * y does not match qr(A).solve(y). Below is a table of different combinations.

mkl omp match
1 1 yes
1 2 yes
2 1 yes
2 2 no

If I compile the eigen_run with the "Link to OpenMP" block commented out. Then everything works.