Ordering problem in Julia 1.9 #90

jd-lara · 2023-04-12T19:30:10Z

Since updating to Julia 1.9 I keep running into this error which didn't happen before, this is an MWE based on the examples

using Pardiso
using SparseArrays
using Random
using Printf
using Test

# Script parameters.
# -----------------
verbose = true
n       = 4  # The number of equations.
m       = 3  # The number of right-hand sides.

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

# Generate a random collection of right-hand sides.
B = rand(n,m)

# Initialize the PARDISO internal data structures.
# ps = PardisoSolver()
ps = MKLPardisoSolver()

if verbose
    set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
end

# If we want, we could just solve the system right now.
# Pardiso.jl will automatically detect the correct matrix type,
# solve the system and free the data
X1 = solve(ps, A, B)

Which results in

RROR: Reordering problem.
Stacktrace:
 [1] check_error(ps::MKLPardisoSolver, err::Int32)
   @ Pardiso ~/.julia/packages/Pardiso/rrIt7/src/mkl_pardiso.jl:80
 [2] ccall_pardiso(ps::MKLPardisoSolver, N::Int64, nzval::Vector{Float64}, colptr::Vector{Int64}, rowval::Vector{Int64}, NRHS::Int64, B::Matrix{Float64}, X::Matrix{Float64})
   @ Pardiso ~/.julia/packages/Pardiso/rrIt7/src/mkl_pardiso.jl:73
 [3] pardiso(ps::MKLPardisoSolver, X::Matrix{Float64}, A::SparseMatrixCSC{Float64, Int64}, B::Matrix{Float64})
   @ Pardiso ~/.julia/packages/Pardiso/rrIt7/src/Pardiso.jl:346
 [4] solve!(ps::MKLPardisoSolver, X::Matrix{Float64}, A::SparseMatrixCSC{Float64, Int64}, B::Matrix{Float64}, T::Symbol)
   @ Pardiso ~/.julia/packages/Pardiso/rrIt7/src/Pardiso.jl:260
 [5] solve
   @ ~/.julia/packages/Pardiso/rrIt7/src/Pardiso.jl:222 [inlined]
 [6] solve(ps::MKLPardisoSolver, A::SparseMatrixCSC{Float64, Int64}, B::Matrix{Float64})
   @ Pardiso ~/.julia/packages/Pardiso/rrIt7/src/Pardiso.jl:221
 [7] top-level scope
   @ REPL[15]:1

The text was updated successfully, but these errors were encountered:

jd-lara · 2023-04-12T19:35:48Z

Could be related to #28

jd-lara · 2023-04-12T21:03:12Z

For more context

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

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

=== 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
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.003684 s
Time spent in reordering of the initial matrix (reorder)         : 0.002538 s
Time spent in symbolic factorization (symbfct)                   : 0.002417 s
Time spent in data preparations for factorization (parlist)      : 0.000160 s
Time spent in copying matrix to internal data structure (A to LU): 0.000000 s
Time spent in factorization step (numfct)                        : 0.006884 s
Time spent in direct solver at solve step (solve)                : 0.002029 s
Time spent in allocation of internal data structures (malloc)    : 0.078254 s
Time spent in additional calculations                            : 0.000208 s
Total time spent                                                 : 0.096174 s

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

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

             number of right-hand sides:    4

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

             gflop/s for the numerical factorization: 0.000001

*** Error in PARDISO  (        reordering_phase) error_num= -180
*** error PARDISO: reordering, symbolic factorization

=== 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


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

Times:
======
Time spent in calculations of symmetric matrix portrait (fulladj): 0.000003 s
Time spent in reordering of the initial matrix (reorder)         : 0.000034 s
Time spent in symbolic factorization (symbfct)                   : 0.000277 s
Time spent in allocation of internal data structures (malloc)    : 0.000179 s
Time spent in additional calculations                            : 0.000006 s
Total time spent                                                 : 0.000499 s

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

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

             number of right-hand sides:    4

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

KristofferC · 2023-04-12T21:04:36Z

So you are saying that just toggling between julia 1.8 and 1.9 gives this behavior, even with the same version of Pardiso and its dependencies?

jd-lara · 2023-04-12T21:10:35Z

Output 1.9

               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.9.0-rc2 (2023-04-01)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> using Pardiso

julia> using SparseArrays

julia> using Random

julia> using Printf

julia> using Test

julia> # Script parameters.
       # -----------------
       verbose = true
true

julia> n       = 4  # The number of equations.
4

julia> m       = 3  # The number of right-hand sides.
3

julia> A = sparse([ 1. 0 -2  3
                    0  5  1  2
                   -2  1  4 -7
                    3  2 -7  5 ])
4×4 SparseMatrixCSC{Float64, Int64} with 14 stored entries:
  1.0   ⋅   -2.0   3.0
   ⋅   5.0   1.0   2.0
 -2.0  1.0   4.0  -7.0
  3.0  2.0  -7.0   5.0

julia> # Generate a random collection of right-hand sides.
       B = rand(n,m)
4×3 Matrix{Float64}:
 0.175154  0.918753  0.364071
 0.302539  0.604431  0.237096
 0.460432  0.470347  0.528947
 0.834289  0.783474  0.451318

julia> # Initialize the PARDISO internal data structures.
       # ps = PardisoSolver()
       ps = MKLPardisoSolver()
MKLPardisoSolver:
	Matrix type: Real nonsymmetric
	Phase: Analysis, numerical factorization, solve, iterative refinement

julia> if verbose
           set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
       end
MESSAGE_LEVEL_ON::MessageLevel = 1

julia> # If we want, we could just solve the system right now.
       # Pardiso.jl will automatically detect the correct matrix type,
       # solve the system and free the data
       X1 = solve(ps, A, B)
*** Error in PARDISO  (        reordering_phase) error_num= -180
*** error PARDISO: reordering, symbolic factorization

=== 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


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

Times:
======
Time spent in calculations of symmetric matrix portrait (fulladj): 0.000783 s
Time spent in reordering of the initial matrix (reorder)         : 0.000288 s
Time spent in symbolic factorization (symbfct)                   : 0.001417 s
Time spent in allocation of internal data structures (malloc)    : 0.060202 s
Time spent in additional calculations                            : 0.000020 s
Total time spent                                                 : 0.062710 s

Statistics:
===========
Parallel Direct Factorization is running on 10 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:                    4
             size of largest supernode:               1
             number of non-zeros in L:                4
             number of non-zeros in U:                1
             number of non-zeros in L+U:              5
             gflop   for the numerical factorization: 0.000000

ERROR: Reordering problem.
Stacktrace:
 [1] check_error(ps::MKLPardisoSolver, err::Int32)
   @ Pardiso ~/.julia/packages/Pardiso/rrIt7/src/mkl_pardiso.jl:80
 [2] ccall_pardiso(ps::MKLPardisoSolver, N::Int64, nzval::Vector{Float64}, colptr::Vector{Int64}, rowval::Vector{Int64}, NRHS::Int64, B::Matrix{Float64}, X::Matrix{Float64})
   @ Pardiso ~/.julia/packages/Pardiso/rrIt7/src/mkl_pardiso.jl:73
 [3] pardiso(ps::MKLPardisoSolver, X::Matrix{Float64}, A::SparseMatrixCSC{Float64, Int64}, B::Matrix{Float64})
   @ Pardiso ~/.julia/packages/Pardiso/rrIt7/src/Pardiso.jl:346
 [4] solve!(ps::MKLPardisoSolver, X::Matrix{Float64}, A::SparseMatrixCSC{Float64, Int64}, B::Matrix{Float64}, T::Symbol)
   @ Pardiso ~/.julia/packages/Pardiso/rrIt7/src/Pardiso.jl:260
 [5] solve
   @ ~/.julia/packages/Pardiso/rrIt7/src/Pardiso.jl:222 [inlined]
 [6] solve(ps::MKLPardisoSolver, A::SparseMatrixCSC{Float64, Int64}, B::Matrix{Float64})
   @ Pardiso ~/.julia/packages/Pardiso/rrIt7/src/Pardiso.jl:221
 [7] top-level scope
   @ REPL[13]:4

Julia>

Output 1.8.5

               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.8.5 (2023-01-08)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> using Pardiso

julia> using SparseArrays

julia> using Random

julia> using Printf

julia> using Test

julia> # Script parameters.
       # -----------------
       verbose = true
true

julia> n       = 4  # The number of equations.
4

julia> m       = 3  # The number of right-hand sides.
3

julia> A = sparse([ 1. 0 -2  3
                    0  5  1  2
                   -2  1  4 -7
                    3  2 -7  5 ])
4×4 SparseMatrixCSC{Float64, Int64} with 14 stored entries:
  1.0   ⋅   -2.0   3.0
   ⋅   5.0   1.0   2.0
 -2.0  1.0   4.0  -7.0
  3.0  2.0  -7.0   5.0

julia> # Generate a random collection of right-hand sides.
       B = rand(n,m)
4×3 Matrix{Float64}:
 0.916727  0.301987  0.333427
 0.683362  0.4072    0.863223
 0.056004  0.599119  0.325091
 0.774732  0.33897   0.733145

julia> # Initialize the PARDISO internal data structures.
       # ps = PardisoSolver()
       ps = MKLPardisoSolver()
MKLPardisoSolver:
	Matrix type: Real nonsymmetric
	Phase: Analysis, numerical factorization, solve, iterative refinement

julia> if verbose
           set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
       end
MESSAGE_LEVEL_ON::MessageLevel = 1

julia> # If we want, we could just solve the system right now.
       # Pardiso.jl will automatically detect the correct matrix type,
       # solve the system and free the data
       X1 = solve(ps, A, B)
=== 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
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.000022 s
Time spent in reordering of the initial matrix (reorder)         : 0.000287 s
Time spent in symbolic factorization (symbfct)                   : 0.001157 s
Time spent in data preparations for factorization (parlist)      : 0.000002 s
Time spent in copying matrix to internal data structure (A to LU): 0.000001 s
Time spent in factorization step (numfct)                        : 0.000000 s
Time spent in allocation of internal data structures (malloc)    : 0.051290 s
Time spent in additional calculations                            : 0.001319 s
Total time spent                                                 : 0.054078 s

Statistics:
===========
Parallel Direct Factorization is running on 10 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

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

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
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.000000 s
Time spent in reordering of the initial matrix (reorder)         : 0.000013 s
Time spent in symbolic factorization (symbfct)                   : 0.002292 s
Time spent in data preparations for factorization (parlist)      : 0.000002 s
Time spent in copying matrix to internal data structure (A to LU): 0.000000 s
Time spent in factorization step (numfct)                        : 0.000097 s
Time spent in direct solver at solve step (solve)                : 0.000729 s
Time spent in allocation of internal data structures (malloc)    : 0.000109 s
Time spent in additional calculations                            : 0.000065 s
Total time spent                                                 : 0.003307 s

Statistics:
===========
Parallel Direct Factorization is running on 10 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.000186

4×3 Matrix{Float64}:
 32.7511   18.2624    13.722
 -1.0142   -0.513195  -0.277561
 11.5617    6.40575    4.79004
 -2.90366  -1.71629   -1.26951

Julia>

For reference here are the installed deps

jd-lara · 2023-04-12T21:30:20Z

@KristofferC upon checking it seems that the LibBlasTrampoline versions are different I don't know is that's the problem

mzy2240 · 2023-05-10T20:24:54Z

Not sure whether it is related: JuliaLinearAlgebra/libblastrampoline#116

staticfloat · 2023-05-10T21:20:06Z

Can you print out LinearAlgebra.BLAS.get_config() after the tests fail?

staticfloat · 2023-05-10T21:23:11Z

Also, can you show exactly what version of MKL_jll you are loading on both Julia versions? pkg> status -m should show it.

KristofferC · 2023-05-11T07:18:14Z

I cannot reproduce this with MKL_jll v2022.2.0+0. It solves ok on both 1.8 and 1.9. What OS is this is on?

jd-lara · 2023-05-12T18:54:31Z

@KristofferC this is on MacOS rosetta.

KristofferC · 2023-05-12T18:58:04Z

And exact version of MKL_jll?

jd-lara · 2023-05-12T18:59:02Z

And exact version of MKL_jll?

MKL_jll v2023.1.0+0 and it seems to downgrade Pardiso.jl to 0.5.0...

KristofferC · 2023-05-12T19:05:41Z

We need to bump MKL_jll compatibility.

KristofferC · 2023-05-12T19:12:52Z

I can reproduce this in rosetta. I get it on MKL_jll 2021, 2022 and 2023. I will test on Julia 1.8 now.

And indeed on 1.8 it works..

KristofferC · 2023-05-15T09:54:50Z

Bisected this to JuliaLang/julia#44360

I'm guessing it has something to do with libblastrampoline_jll now being loaded in the sysimage?

staticfloat · 2023-05-16T19:15:13Z

This should be fixed by JuliaLang/julia#49842, you can test it by downloading the x86_64 macOS build here: https://buildkite.com/organizations/julialang/pipelines/julia-master/builds/23999/jobs/018825e2-7226-49a1-9f7e-ecec229c537d/artifacts/018825ee-b753-4a32-ba27-b79c095e79dc

ViralBShah · 2024-02-15T05:41:40Z

Please reopen if still an issue.

jd-lara · 2024-02-25T19:26:11Z

Seems this is fixed in Pardiso 0.5.5 and Julia 1.10

KristofferC mentioned this issue May 15, 2023

bump MKL compat #91

Merged

KristofferC mentioned this issue May 15, 2023

Simplification of libblastrampoline stuff in LinearAlgebra JuliaLang/julia#44360

Merged

staticfloat mentioned this issue May 16, 2023

Add LazilyLoadedLibrary type to Libdl JuliaLang/julia#49843

Closed

ViralBShah closed this as completed Feb 15, 2024

mipals mentioned this issue Feb 27, 2024

Using JuMP causes Pardiso crash #88

Open

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Ordering problem in Julia 1.9 #90

Ordering problem in Julia 1.9 #90

jd-lara commented Apr 12, 2023 •

edited

Loading

jd-lara commented Apr 12, 2023

jd-lara commented Apr 12, 2023

KristofferC commented Apr 12, 2023

jd-lara commented Apr 12, 2023 •

edited

Loading

jd-lara commented Apr 12, 2023

mzy2240 commented May 10, 2023

staticfloat commented May 10, 2023 •

edited

Loading

staticfloat commented May 10, 2023

KristofferC commented May 11, 2023

jd-lara commented May 12, 2023

KristofferC commented May 12, 2023

jd-lara commented May 12, 2023

KristofferC commented May 12, 2023

KristofferC commented May 12, 2023 •

edited

Loading

KristofferC commented May 15, 2023 •

edited

Loading

staticfloat commented May 16, 2023

ViralBShah commented Feb 15, 2024

jd-lara commented Feb 25, 2024

Ordering problem in Julia 1.9 #90

Ordering problem in Julia 1.9 #90

Comments

jd-lara commented Apr 12, 2023 • edited Loading

jd-lara commented Apr 12, 2023

jd-lara commented Apr 12, 2023

KristofferC commented Apr 12, 2023

jd-lara commented Apr 12, 2023 • edited Loading

jd-lara commented Apr 12, 2023

mzy2240 commented May 10, 2023

staticfloat commented May 10, 2023 • edited Loading

staticfloat commented May 10, 2023

KristofferC commented May 11, 2023

jd-lara commented May 12, 2023

KristofferC commented May 12, 2023

jd-lara commented May 12, 2023

KristofferC commented May 12, 2023

KristofferC commented May 12, 2023 • edited Loading

KristofferC commented May 15, 2023 • edited Loading

staticfloat commented May 16, 2023

ViralBShah commented Feb 15, 2024

jd-lara commented Feb 25, 2024

jd-lara commented Apr 12, 2023 •

edited

Loading

jd-lara commented Apr 12, 2023 •

edited

Loading

staticfloat commented May 10, 2023 •

edited

Loading

KristofferC commented May 12, 2023 •

edited

Loading

KristofferC commented May 15, 2023 •

edited

Loading