JuliaSparse / Pardiso.jl

Calling the PARDISO library from Julia
Other
100 stars 27 forks source link

Linking with libmkl_rt fails on macOS Sierra + MKL 18 #28

Closed mcovalt closed 4 years ago

mcovalt commented 6 years ago

The currently selected MKL library is not working for my system:

global const libmkl_core = Libdl.dlopen(string(MKLROOT, "/lib/libmkl_rt"), Libdl.RTLD_GLOBAL)

If using this library, I get an Input inconsistent error when running the tests. Changing this library to:

global const libmkl_core = Libdl.dlopen(string(MKLROOT, "/lib/libmkl_intel_lp64"), Libdl.RTLD_GLOBAL)

Manages to fix the problem (for me).

I'm no expert on Intel MKL linking, but from what I thought, libmkl_rt is the "catch all" link that loads all the correct layers (threading, compute, and interface). I'm surprised that just linking the interface layer (libmkl_intel_lp64) allows Pardiso to run nicely, but it does.

KristofferC commented 6 years ago

For me, linking works fine but then, using a simple example, I get:

julia> solve!(ps, X, A, B)
julia(17480,0x7fff8c442340) malloc: *** mach_vm_map(size=18446744059401457664) failed (error code=3)
*** error: can't allocate region
*** set a breakpoint in malloc_error_break to debug
ERROR: Not enough memory.
Stacktrace:
 [1] check_error(::Pardiso.MKLPardisoSolver, ::Int32) at /Users/kristoffer/.julia/v0.6/Pardiso/src/mkl_pardiso.jl:92
 [2] ccall_pardiso(::Pardiso.MKLPardisoSolver, ::Int32, ::Array{Float64,1}, ::Array{Int32,1}, ::Array{Int32,1}, ::Int32, ::Array{Float64,2}, ::Array{Float64,2}) at /Users/kristoffer/.julia/v0.6/Pardiso/src/mkl_pardiso.jl:87
 [3] pardiso(::Pardiso.MKLPardisoSolver, ::Array{Float64,2}, ::SparseMatrixCSC{Float64,Int64}, ::Array{Float64,2}) at /Users/kristoffer/.julia/v0.6/Pardiso/src/Pardiso.jl:264
 [4] solve!(::Pardiso.MKLPardisoSolver, ::Array{Float64,2}, ::SparseMatrixCSC{Float64,Int64}, ::Array{Float64,2}, ::Symbol) at /Users/kristoffer/.julia/v0.6/Pardiso/src/Pardiso.jl:226
 [5] solve!(::Pardiso.MKLPardisoSolver, ::Array{Float64,2}, ::SparseMatrixCSC{Float64,Int64}, ::Array{Float64,2}) at /Users/kristoffer/.julia/v0.6/Pardiso/src/Pardiso.jl:176
KristofferC commented 6 years ago

Running the example at https://github.com/JuliaSparse/Pardiso.jl/blob/master/examples/exampleherm.jl I get:

julia> pardiso(ps, A_pardiso, b)
*** Error in PARDISO  (        reordering_phase) error_num= -180
*** error PARDISO: reordering, symbolic factorization

=== PARDISO: solving a Hermitian positive definite system ===
1-based array indexing is turned ON
PARDISO double precision computation is turned ON
METIS algorithm at reorder step is turned ON

Summary: ( reordering phase )
================

Times:
======
Time spent in calculations of symmetric matrix portrait (fulladj): 0.017118 s
Time spent in reordering of the initial matrix (reorder)         : 0.004487 s
Time spent in symbolic factorization (symbfct)                   : 0.002134 s
Time spent in allocation of internal data structures (malloc)    : 0.066873 s
Time spent in additional calculations                            : 0.000942 s
Total time spent                                                 : 0.091554 s

Statistics:
===========
Parallel Direct Factorization is running on 4 OpenMP

< Linear system Ax = b >
             number of equations:           100
             number of non-zeros in A:      199
             number of non-zeros in A (%): 1.990000

             number of right-hand sides:    1

< Factors L and U >
             number of columns for each panel: 128
             number of independent subgraphs:  0
< Preprocessing with state of the art partitioning metis>
             number of supernodes:                    100
             size of largest supernode:               1
             number of non-zeros in L:                34376
             number of non-zeros in U:                1
             number of non-zeros in L+U:              34377
ERROR: Reordering problem.
Stacktrace:
 [1] check_error(::Pardiso.MKLPardisoSolver, ::Int32) at /Users/kristoffer/.julia/v0.6/Pardiso/src/mkl_pardiso.jl:93
 [2] ccall_pardiso(::Pardiso.MKLPardisoSolver, ::Int32, ::Array{Complex{Float64},1}, ::Array{Int32,1}, ::Array{Int32,1}, ::Int32, ::Array{Complex{Float64},2}, ::Array{Complex{Float64},1}) at /Users/kristoffer/.julia/v0.6/Pardiso/src/mkl_pardiso.jl:87
 [3] pardiso(::Pardiso.MKLPardisoSolver, ::Array{Complex{Float64},1}, ::SparseMatrixCSC{Complex{Float64},Int64}, ::Array{Complex{Float64},2}) at /Users/kristoffer/.julia/v0.6/Pardiso/src/Pardiso.jl:264
 [4] pardiso(::Pardiso.MKLPardisoSolver, ::SparseMatrixCSC{Complex{Float64},Int64}, ::Array{Complex{Float64},2}) at /Users/kristoffer/.julia/v0.6/Pardiso/src/Pardiso.jl:269
KristofferC commented 6 years ago

Using normal pardiso it prints:

julia> pardiso(ps, A_pardiso, b)
================  PARDISO: solving a  Herm. pos. def.  system  ================

Summary PARDISO 5.0.0: ( reorder to reorder )
=======================

Times:
======

      Time fulladj: 0.000499 s
      Time reorder: 0.001446 s
      Time symbfct: 0.002001 s
      Time parlist: 0.000051 s
      Time malloc : 0.000926 s
      Time total  : 0.007176 s total - sum: 0.002253 s

Statistics:
===========
 < Parallel Direct Factorization with #cores: >              1
 <                                and #nodes: >              1
 < Numerical Factorization with Level-3 BLAS performance >

 < Linear system Ax = b>
             #equations:                                     100
             #non-zeros in A:                                199
             non-zeros in A (%):                             1.990000
             #right-hand sides:                              1

< Factors L and U >
             #columns for each panel:                        80
             # of independent subgraphs:                     0
 < preprocessing with state of the art partitioning metis>
             #supernodes:                                    96
             size of largest supernode:                      3
             number of nonzeros in L                         290
             number of nonzeros in U                         1
             number of nonzeros in L+U                       291
             number of perturbed pivots                      0
             number of nodes in solve                        100
             Gflop   for the numerical factorization:        0.000003
KristofferC commented 4 years ago

Ref https://github.com/oxfordcontrol/osqp-python/issues/14, maybe there is some other MKL library loaded? Tried locally on mac and on CI and it seems to work now.