Problems when using MKL BLAS

[KristofferC]

Right now, when using Julia with MKL BLAS, PARDISO complains about not finding dgemv_64_.

Related issue: https://github.com/JuliaLang/julia/issues/4923

[KristofferC]

Also, on the computer with the normal blas installation, leaving out this line: https://github.com/KristofferC/Pardiso.jl/blob/master/src/Pardiso.jl#L13 gives:

symbol lookup error:: undefined symbol: dswap_

Changing it to dlopen libopenblas does not help. I thought dswap was a part of blas which should be linked in to julia. Is it some name mangling thing?

@tkelman if you have any spare time, any ideas?

[tkelman]

Is that the right line? export solve, solve! ?

So the issue with #4923 is that on 64 bit systems, Julia made the unfortunate decision to use the 64-bit-integer interface to BLAS, which is nonstandard (there's some precedent with Matlab, but it causes issues there as well). You may want to look into which blas library pardiso was compiled against, and whether it was using the lp64 interface or the ilp64 interface. Those two interfaces should really not be using the same function names, as it causes all sorts of problems at runtime. So for the copy of openblas that we build with Julia when we're using 64 bit ints, we now rename the symbols so there's no possibility of confusion.

If pardiso was built against the more common lp64 blas where integers are 32 bit, then you want to dlopen an lp64 blas library with conventional symbol names - install the best blas library you can get via your package manager and dlopen'ing libblas should work.

[tkelman]

Also relevant here is exactly how are you building Julia, against which blas library, etc. Julia from most package managers should be using 32 bit integers in blas for compatibility with existing distribution blas libraries. A source build or the generic binaries with the bundled openblas will use 64 bit integers. If you build Julia itself using MKL, then I think it can be configured to use either of the interfaces depending on the build-time setting of USE_BLAS64.

[KristofferC]

That was the wrong line, this is what I meant https://github.com/KristofferC/Pardiso.jl/blob/master/src/Pardiso.jl#L16

I would guess Pardiso uses lp64 since the index vectors you pass for the sparse matrices should be in 32 bit.

Thanks for the explanation, I realize now why I needed to load libblas and why julias openblas did not work,

Now I have to figure out why the pardiso solver is super slow...

using Pardiso
ps = PardisoSolver()
for i = [0.1, 0.2, 0.3, 0.4, 0.5]
    n = 3 + i
    A = sparse(eye(trunc(Int, 10^n)) + sprand(trunc(Int,10^n), trunc(Int, 10^n), 0.001))
    x = rand(trunc(Int, 10^n))
    X = similar(x)
    println("-----")
    println("pardiso:")
    @time solve!(ps, X, A, x)
    println("umfpack:")
    @time A\x
    println("-----")
end

-----
pardiso:
elapsed time: 0.012262463 seconds (16 kB allocated)
umfpack:
elapsed time: 0.001619705 seconds (87 kB allocated)
-----
-----
pardiso:
elapsed time: 0.023350006 seconds (22 kB allocated)
umfpack:
elapsed time: 0.002647499 seconds (114 kB allocated)
-----
-----
pardiso:
elapsed time: 0.055432557 seconds (31 kB allocated)
umfpack:
elapsed time: 0.006090313 seconds (153 kB allocated)
-----
-----
pardiso:
elapsed time: 0.110756689 seconds (45 kB allocated)
umfpack:
elapsed time: 0.021798734 seconds (211 kB allocated)
-----
-----
pardiso:
elapsed time: 0.355305781 seconds (64 kB allocated)
umfpack:
elapsed time: 0.061372032 seconds (294 kB allocated)
-----

[tkelman]

Ah, okay, dlopen'ing libblas makes more sense, the system version should be lp64 in all cases that I'm aware of. If the sparse index vectors are 32-bit ints that's probably a good sign that the blas calls are also lp64, but it's not absolutely guaranteed. The external sparse API could potentially differ from the internal dense API that gets used on sub-blocks.

Check what's going on with multithreading, if both pardiso and your blas library (assuming your system blas is something multithreaded like openblas) are trying to use threads at the same time that can cause nasty contention.

[KristofferC]

This is the linking advice which seems to mean lp64?

ifort <source/objects files> -o <executable>-L <Path to directory of PARDISO> -lpardiso500-INTEL1301-X86-64-L <Path to directory of LAPACK/BLAS>-l <fast LAPACK and BLAS libraries>-L/opt/intel/mkl/lib/intel64 -lmkl gf-lmkl core -lm -lgfortran -fopenmp lp64 -lmkl sequential

However, they also say "The sparse direct solvers use 8-byte integer length internally that allows storage of factors with more than 2^32 nonzeros elements".

I had my suspicion about threads fighting being the problem. I will look into it more tomorrow, late here now :)

Thanks again for all your help. Really appreciated.

[KristofferC]

@tkelman I am thinking about how to get MKL Pardiso to be implemented in a good way.

Firstly, do you think it makes sense to have them in the same package or have separate MKLPardiso.jl package. There will be some code duplication but maybe it is the easiest?

If not, the following is a little bit of a problem: Let's say that we detect MKL Blas. That does not automatically mean that the user want's to use the MKL version of Pardiso. Maybe he/she is just wants to use MKL Blas with official Paridso.

So there need to be some way to specify which Pardiso version you want to run. This could be done maybe with two types, one ParidsoSolver as I have now and now PardisoMKLSolver and then you can dispatch on the type for the things that are different between the two versions.

My guess is that there will be quite a lot of

if typeof(ps) == PardisoSolver
    blah
elseif typeof(ps) == PardisoMKLSolver
    bleh
end

but I guess that is unavoidable.

Any comments?

[tkelman]

Good to think about, I'm not sure I have any meaningful answers for you. If a large amount of code will be shared, I think using the same package makes sense. But as you say, auto-detection could get tricky. So different types sounds like a good idea to me. Usually you don't need to do if typeof(ps) == PardisoSolver, you can use dispatch for that. Put common code in separate helper functions, then dispatch on different methods for the code that needs to be different.

[KristofferC]

Ok, I have a branch up now where MKL Pardiso also works. Still having some troubles with the julia that is actually built with MKL Blas.

However, I can't even use Julias inbuilt umfpack functions on that build so it seems kinda nuked..

julia-0.4: symbol lookup error: /home/kristoffer/julia/usr/bin/../lib/libumfpack.so: undefined symbol: dgemv_64_

[tkelman]

Are you switching back and forth between MKL and OpenBLAS within the same clone? You'll need to do make -C deps distclean-suitesparse distclean-arpack any time you try to switch.

[RoyiAvital]

If I use JuliaPro with Intel MKL, will it work out of the box and use the built in MKL Library?