MKL Pardiso with multiple threads

[migarstka]

I am encountering a problem when trying to set the number of threads for the MKLPardisoSolver. In this example, I set the number of threads to 2:

using Pardiso, LinearAlgebra, SparseArrays
n = 4
m = 3

A = sparse([ 1. 0 -2  3
             0  5  1  2
            -2  1  4 -7
             3  2 -7  5 ])

B = rand(n,m)

ps = MKLPardisoSolver()

set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
set_nprocs!(ps, 2)

X1 = solve(ps, A, B)
get_nprocs(ps)

However, the last line get_nprocs(ps) returns 1. Running the script, produces the following output:

=== PARDISO is running in In-Core mode, because iparam(60)=0 ===

Percentage of computed non-zeros for LL^T factorization
 25 %  100 % 
*** Error in PARDISO  ( numerical_factorization) error_num= -1
*** Error in PARDISO: zero or negative pivot, A is not SPD-matrix

=== PARDISO: solving a symmetric 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
Scaling is turned ON
Matching is turned ON
Single-level factorization algorithm is turned ON


Summary: ( starting phase is reordering, ending phase is solution )
================

Times:
======
Time spent in calculations of symmetric matrix portrait (fulladj): 0.000002 s
Time spent in reordering of the initial matrix (reorder)         : 0.000068 s
Time spent in symbolic factorization (symbfct)                   : 0.000015 s
Time spent in data preparations for factorization (parlist)      : 0.000000 s
Time spent in copying matrix to internal data structure (A to LU): 0.000000 s
Time spent in factorization step (numfct)                        : 0.000000 s
Time spent in allocation of internal data structures (malloc)    : 0.000083 s
Time spent in additional calculations                            : 0.000061 s
Total time spent                                                 : 0.000229 s

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

< Linear system Ax = b >
             number of equations:           4
             number of non-zeros in A:      9
             number of non-zeros in A (%): 56.250000

             number of right-hand sides:    3

< 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:                    2
             size of largest supernode:               3
             number of non-zeros in L:                12
             number of non-zeros in U:                1
             number of non-zeros in L+U:              13
             gflop   for the numerical factorization: 0.000000


Percentage of computed non-zeros for LL^T factorization
 25 %  100 % 

=== PARDISO: solving a symmetric indefinite system ===
1-based array indexing is turned ON
PARDISO double precision computation is turned ON
METIS algorithm at reorder step is turned ON
Scaling is turned ON
Matching is turned ON
Single-level factorization algorithm is turned ON


Summary: ( starting phase is reordering, ending phase is solution )
================

Times:
======
Time spent in calculations of symmetric matrix portrait (fulladj): 0.000002 s
Time spent in reordering of the initial matrix (reorder)         : 0.000000 s
Time spent in symbolic factorization (symbfct)                   : 0.000015 s
Time spent in data preparations for factorization (parlist)      : 0.000001 s
Time spent in copying matrix to internal data structure (A to LU): 0.000000 s
Time spent in factorization step (numfct)                        : 0.000056 s
Time spent in direct solver at solve step (solve)                : 0.000025 s
Time spent in allocation of internal data structures (malloc)    : 0.000297 s
Time spent in additional calculations                            : 0.000007 s
Total time spent                                                 : 0.000403 s

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

< Linear system Ax = b >
             number of equations:           4
             number of non-zeros in A:      9
             number of non-zeros in A (%): 56.250000

             number of right-hand sides:    3

< 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:                    2
             size of largest supernode:               3
             number of non-zeros in L:                12
             number of non-zeros in U:                1
             number of non-zeros in L+U:              13
             gflop   for the numerical factorization: 0.000000

             gflop/s for the numerical factorization: 0.000321

1

Again it seems like only one thread is used. I am using a Late 2013 MacBook Pro with Julia 1.1. Multicore Pardiso 6.0 works for me

[KristofferC]

The set_nprocs! calls into mkl_domain_set_num_threads with the domain 4 which should be the Pardiso domain. Not sure why that wouldn't work.

[KristofferC]

Something seems off with your MKL, it works for me

julia> ps = MKLPardisoSolver();

julia> set_nprocs!(ps, 2)

julia> get_nprocs(ps)
2

julia> set_nprocs!(ps, 8)

julia> get_nprocs(ps)
8

Perhaps upgrade your MKL?

[migarstka]

Sorry to open this again.
I updated and reinstalled Julia/Pardiso and IntelMKL and still can't use set_nprocs!(). All the unit tests pass except for

Pardiso.MklInt = Int32
Testing DataType[MKLPardisoSolver, PardisoSolver]
Test Summary: | Pass  Total
solving       |  192    192
The factors have 14 nonzero entries.
The matrix has 0 positive and 0 negative eigenvalues.
PARDISO performed 0 iterative refinement steps.
The maximum residual for the solution X is 1.33e-15.
The factors have 13 nonzero entries.
The matrix has 3 positive and 1 negative eigenvalues.
PARDISO performed 0 iterative refinement steps.
The maximum residual for the solution X is 1.98e-14.
The factors have 289 nonzero entries.
PARDISO performed 0 iterative refinement steps.
The maximum residual for the solution is 4.48e-16.
procs: Test Failed at /Users/Micha/.julia/packages/Pardiso/yZsYO/test/runtests.jl:70
  Expression: get_nprocs(ps) == 2
   Evaluated: 1 == 2
Stacktrace:
 [1] top-level scope at /Users/Micha/.julia/packages/Pardiso/yZsYO/test/runtests.jl:70
 [2] top-level scope at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/Test/src/Test.jl:1113
 [3] top-level scope at /Users/Micha/.julia/packages/Pardiso/yZsYO/test/runtests.jl:67
Test Summary: | Pass  Fail  Total
procs         |    1     1      2
ERROR: LoadError: Some tests did not pass: 1 passed, 1 failed, 0 errored, 0 broken.
in expression starting at /Users/Micha/.julia/packages/Pardiso/yZsYO/test/runtests.jl:65
ERROR: Package Pardiso errored during testing

Any idea what else I could try?

[KristofferC]

Does it work if you do not have the non-MKL Pardiso loaded? Remove ENV["JULIA_PARDISO"], run build Pardiso and try re-run the tests.

I had some issues with this when both having MKL and non-MKL Pardiso loaded in the same session (see #67). See the note in the README:

Note: Weird errors and problems with MKL Pardiso has been observed when Pardiso 6.0 is enabled (likely because some library that is needed by Pardiso 6.0 is problematic with MKL). If you want to use MKL Pardiso it is better ot just disable Paridso 6.0 by not setting the environment variable JULIA_PARDISO (and rerunning build Pardiso).

[migarstka]

Thanks a lot. That seems to have solved the problem. That's a weird bug.

Project Pardiso disabled:

julia> Pardiso.show_build_log()
Pardiso library
===============
Looking for libraries with name: libpardiso600-MACOS-X86-64.dylib.
INFO: use the `JULIA_PARDISO` environment variable to set a path to the folder where the Pardiso library is located
Looking in "/Users/Micha/.julia/packages/Pardiso/yZsYO/deps" for libraries
did not find libpardiso, assuming PARDISO 5/6 is not installed

MKL Pardiso
=============
found MKLROOT environment varaible, enabling local MKL

and

ps = MKLPardisoSolver();

julia> set_nprocs!(ps, 2)

julia> get_nprocs(ps)
2

[KristofferC]

We have to open some libraries for non-MKL Pardiso to work and I am guessing one of them causes some issue with MKL...