Possible threading bug with MKL on macOS using Julia 1.7.0

[BenjaminBorn]

I encountered a possible MKL threading bug in Julia 1.7.0 on macOS but don't know if it is the same as the OpenBLAS threading issue (e.g. #890 ), so please close if this is a duplicate.

The following code is a simplified version of code that generates impulse response functions from a state-space model. It works with 1 BLAS thread but the matrix based iteration x2 explodes with any higher thread number, even though the relevant row of the matrix only has one non-zero entry of 0.9 and should give the same path as x1, which is essentially an AR1. On Linux, the code runs for any number of threads and x1 and x2 follow the same impulse response.

using Plots, MKL, LinearAlgebra
BLAS.set_num_threads(2)

var_index = 150 # index of variable of interest
LOMstate = randn(300, 300) # law of motion
LOMstate[var_index, :] .= 0.0 # set the law of motion for the variable of interest to an AR with coefficient 0.9
LOMstate[var_index, var_index] = 0.9

x0 = zeros(size(LOMstate, 1), 1)
x0[var_index] = 1.0 # shock to the variable of interest

irf_horizon = 40
x1 = x0 * ones(1, irf_horizon + 1)
x2 = x0 * ones(1, irf_horizon + 1)

for t = 1:irf_horizon
    x1[var_index, t+1] = LOMstate[var_index, var_index] * x1[var_index, t] # just iterate over AR1
    x2[:, t+1] = LOMstate * x2[:, t] # iterate full system but variable of interest should be unaffected by rest of system
end

plt1 = plot(x1[var_index, :], label = "IRF AR1 based")
plt1 = plot!(x2[var_index, :], label = "IRF matrix based")

With 1 thread

With 2 or more threads (might look different due to randomLOMstate matrix)

My system:

Julia Version 1.7.0
Commit 3bf9d17731 (2021-11-30 12:12 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin19.5.0)
  CPU: Intel(R) Xeon(R) W-2191B CPU @ 2.30GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, skylake-avx512)
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS = 1 

I'm happy to provide more information.

[BenjaminBorn]

This is still present in 1.7.1 and MKL. OpenBLAS works for any number of threads.

[oscardssmith]

is this our fault or mkl's

[ViralBShah]

This rings a bell. I remember @staticfloat saying this is a conflict between Intel's openmp and GCC's. Perhaps similar to or the same as #845

[ViralBShah]

Out of curiosity, can you try on Linux and Windows?

[BenjaminBorn]

On Linux, the code runs for any number of threads and x1 and x2 follow the same impulse response. I don't have access to a Windows machine, so cannot test.

[ViralBShah]

After using MKL, can you try

MKL.set_threading_layer(MKL.THREADING_GNU)

And see if that works as expected?

[BenjaminBorn]

Unfortunately, MKL.set_threading_layer(MKL.THREADING_GNU) does not solve the issue.

[ViralBShah]

How about THREADING_SEQUENTIAL instead? Maybe make the change in MKL.jl as in https://github.com/JuliaLinearAlgebra/MKL.jl/pull/99/files directly.

[BenjaminBorn]

THREADING_SEQUENTIAL works if I make the change directly in MKL.jl as you suggested. It did not work via MKL.set_threading_layer. Setting THREADING_GNU directly in MKL.jl results in an error (ERROR: LoadError: InitError: return value was -1) during precompilation.

Does THREADING_SEQUENTIAL essentially set the number of BLAS threads to 1? Because that's what I now get from BLAS.get_num_threads().

[ViralBShah]

That is right. GNU does not work if you set it in MKL.jl as I found out too. SEQUENTIAL will set MKL BLAS threads to 1.

[ViralBShah]

We're now defaulting to single threaded MKL on mac because of too many issues. JuliaLinearAlgebra/MKL.jl#98

I'm closing because this is not a Julia bug.