Pardiso.jl icon indicating copy to clipboard operation
Pardiso.jl copied to clipboard

MKL Pardiso with multiple threads

Open migarstka opened this issue 5 years ago • 6 comments

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

migarstka avatar May 20 '19 08:05 migarstka

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 avatar May 20 '19 09:05 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?

KristofferC avatar Feb 07 '20 18:02 KristofferC

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?

migarstka avatar May 19 '20 10:05 migarstka

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 https://github.com/JuliaSparse/Pardiso.jl/pull/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).

KristofferC avatar May 19 '20 11:05 KristofferC

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

migarstka avatar May 19 '20 11:05 migarstka

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

KristofferC avatar May 19 '20 12:05 KristofferC