[MKL] Get/Set number of threads - should get/set for BLAS, not for the whole MKL

[ylvain]

If MKL.jl (BLAS/LAPACK) is used with other packages providing access to others MKL features (FFT, VML, etc.) then a call to BLAS.set_num_threads(_) will set the maximum number of threads used by all this packages.

So if you need to manage threading yourself, e.g. for an embarrassingly parallelizable algo working on a lots of small matrix, then you would do a BLAS.set_num_threads(1). But at this point your FFT and VML calls also become single threaded.

MKL has a notion of "domain", see: mkl\include\mkl_types.h

• define MKL_DOMAIN_ALL 0
• define MKL_DOMAIN_BLAS 1
• define MKL_DOMAIN_FFT 2
• define MKL_DOMAIN_VML 3
• define MKL_DOMAIN_PARDISO 4

All service settings functions should work on the domain they belong.
See MKL support-functions

The correct call:

set_num_threads(numthreads::Int) =
    Bool(ccall((:MKL_Domain_Set_Num_Threads, MKL_jll.libmkl_rt),
    Cuint, (Cint,Cint), numthreads, MKL_DOMAIN_BLAS))

Since the method is BLAS.set_num_threads, not VML or FFT...
The current behavior interfere with VML, and as it is with the PR JuliaMath/IntelVectorMath.jl#55.

I think one should special case MKL in
https://github.com/JuliaLinearAlgebra/libblastrampoline/blob/main/src/threading.c
(and maybe elsewhere) since it's more than a BLAS/LAPACK library.

[staticfloat]

Hmmm, I think I agree! I would be happy to accept a PR making this change, to force MKL to use MKL_Domain_Set_Num_Threads instead of MKL_Set_Num_Threads.

[ViralBShah]

You also need to use MKL_Domain_Get_Max_Threads and since the Domain value has to be passed, the internal APIs inside LBT need to be adjusted - I think they expect these APIs roughly to be similar across different BLAS, but now that will change. Should be straightforward to incorporate, but it appears to need a bit more work than just find and replace.

[staticfloat]

Yes, we need to essentially special-case this.

[ViralBShah]

This is pretty difficult to do in the current system. It would be best to not special case this mechanism in libblastrampoline, and instead let the various BLAS packages provide this at the Julia level.