Merge pull request #505 from j-fu/rework_iparm_setting

Rework iparm setting
SciML · Jun 5, 2024 · 1f8eb0d · 1f8eb0d
2 parents 270b56d + a339509
commit 1f8eb0d
Show file tree

Hide file tree

Showing 3 changed files with 93 additions and 42 deletions.
diff --git a/ext/LinearSolvePardisoExt.jl b/ext/LinearSolvePardisoExt.jl
@@ -22,23 +22,28 @@ function LinearSolve.init_cacheval(alg::PardisoJL,
         reltol,
         verbose::Bool,
         assumptions::LinearSolve.OperatorAssumptions)
-    @unpack nprocs, solver_type, matrix_type, iparm, dparm = alg
+    @unpack nprocs, solver_type, matrix_type, cache_analysis, iparm, dparm = alg
     A = convert(AbstractMatrix, A)
 
+    transposed_iparm = 1
     solver = if Pardiso.PARDISO_LOADED[]
         solver = Pardiso.PardisoSolver()
+        Pardiso.pardisoinit(solver)
         solver_type !== nothing && Pardiso.set_solver!(solver, solver_type)
 
         solver
     else
         solver = Pardiso.MKLPardisoSolver()
+        Pardiso.pardisoinit(solver)
         nprocs !== nothing && Pardiso.set_nprocs!(solver, nprocs)
 
+        # for mkl 1 means conjugated and 2 means transposed.
+        # https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2024-0/pardiso-iparm-parameter.html#IPARM37
+        transposed_iparm = 2
+
         solver
     end
 
-    Pardiso.pardisoinit(solver) # default initialization
-
     if matrix_type !== nothing
         Pardiso.set_matrixtype!(solver, matrix_type)
     else
@@ -52,22 +57,6 @@ function LinearSolve.init_cacheval(alg::PardisoJL,
     end
     verbose && Pardiso.set_msglvl!(solver, Pardiso.MESSAGE_LEVEL_ON)
 
-    # pass in vector of tuples like [(iparm::Int, key::Int) ...]
-    if iparm !== nothing
-        for i in iparm
-            Pardiso.set_iparm!(solver, i...)
-        end
-    end
-
-    if dparm !== nothing
-        for d in dparm
-            Pardiso.set_dparm!(solver, d...)
-        end
-    end
-
-    # Make sure to say it's transposed because its CSC not CSR
-    Pardiso.set_iparm!(solver, 12, 1)
-
     #=
     Note: It is recommended to use IPARM(11)=1 (scaling) and IPARM(13)=1 (matchings) for
     highly indefinite symmetric matrices e.g. from interior point optimizations or saddle point problems.
@@ -79,23 +68,44 @@ function LinearSolve.init_cacheval(alg::PardisoJL,
     be changed to Pardiso.ANALYSIS_NUM_FACT in the solver loop otherwise instabilities
     occur in the example https://github.com/SciML/OrdinaryDiffEq.jl/issues/1569
     =#
-    Pardiso.set_iparm!(solver, 11, 0)
-    Pardiso.set_iparm!(solver, 13, 0)
-
-    Pardiso.set_phase!(solver, Pardiso.ANALYSIS)
+    if cache_analysis
+        Pardiso.set_iparm!(solver, 11, 0)
+        Pardiso.set_iparm!(solver, 13, 0)
+    end
 
     if alg.solver_type == 1
         # PARDISO uses a numerical factorization A = LU for the first system and
         # applies these exact factors L and U for the next steps in a
         # preconditioned Krylov-Subspace iteration. If the iteration does not
         # converge, the solver will automatically switch back to the numerical factorization.
-        Pardiso.set_iparm!(solver, 3, round(Int, abs(log10(reltol)), RoundDown) * 10 + 1)
+        # Be aware that in the intel docs, iparm indexes are one lower.
+        Pardiso.set_iparm!(solver, 4, round(Int, abs(log10(reltol)), RoundDown) * 10 + 1)
     end
 
-    Pardiso.pardiso(solver,
-        u,
-        SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A)),
-        b)
+    # pass in vector of tuples like [(iparm::Int, key::Int) ...]
+    if iparm !== nothing
+        for i in iparm
+            Pardiso.set_iparm!(solver, i...)
+        end
+    end
+
+    if dparm !== nothing
+        for d in dparm
+            Pardiso.set_dparm!(solver, d...)
+        end
+    end
+
+    # Make sure to say it's transposed because its CSC not CSR
+    # This is also the only value which should not be overwritten by users
+    Pardiso.set_iparm!(solver, 12, transposed_iparm)
+
+    if cache_analysis
+        Pardiso.set_phase!(solver, Pardiso.ANALYSIS)
+        Pardiso.pardiso(solver,
+            u,
+            SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A)),
+            b)
+    end
 
     return solver
 end
@@ -105,13 +115,14 @@ function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::PardisoJL; kwargs
     A = convert(AbstractMatrix, A)
 
     if cache.isfresh
-        Pardiso.set_phase!(cache.cacheval, Pardiso.NUM_FACT)
+        phase = alg.cache_analysis ? Pardiso.NUM_FACT : Pardiso.ANALYSIS_NUM_FACT
+        Pardiso.set_phase!(cache.cacheval, phase)
         Pardiso.pardiso(cache.cacheval, A, eltype(A)[])
         cache.isfresh = false
     end
-
     Pardiso.set_phase!(cache.cacheval, Pardiso.SOLVE_ITERATIVE_REFINE)
     Pardiso.pardiso(cache.cacheval, u, A, b)
+
     return SciMLBase.build_linear_solution(alg, cache.u, nothing, cache)
 end
 

diff --git a/src/extension_algs.jl b/src/extension_algs.jl
@@ -86,6 +86,7 @@ end
 ```julia
 MKLPardisoFactorize(; nprocs::Union{Int, Nothing} = nothing,
     matrix_type = nothing,
+    cache_analysis = false,
     iparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing,
     dparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing)
 ```
@@ -98,7 +99,11 @@ A sparse factorization method using MKL Pardiso.
 
 ## Keyword Arguments
 
-For the definition of the keyword arguments, see the Pardiso.jl documentation.
+Setting `cache_analysis = true` disables Pardiso's scaling and matching defaults
+and caches the result of the initial analysis phase for all further computations
+with this solver.
+
+For the definition of the other keyword arguments, see the Pardiso.jl documentation.
 All values default to `nothing` and the solver internally determines the values
 given the input types, and these keyword arguments are only for overriding the
 default handling process. This should not be required by most users.
@@ -109,6 +114,7 @@ MKLPardisoFactorize(; kwargs...) = PardisoJL(; solver_type = 0, kwargs...)
 ```julia
 MKLPardisoIterate(; nprocs::Union{Int, Nothing} = nothing,
     matrix_type = nothing,
+    cache_analysis = false,
     iparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing,
     dparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing)
 ```
@@ -121,7 +127,11 @@ A mixed factorization+iterative method using MKL Pardiso.
 
 ## Keyword Arguments
 
-For the definition of the keyword arguments, see the Pardiso.jl documentation.
+Setting `cache_analysis = true` disables Pardiso's scaling and matching defaults
+and caches the result of the initial analysis phase for all further computations
+with this solver.
+
+For the definition of the other keyword arguments, see the Pardiso.jl documentation.
 All values default to `nothing` and the solver internally determines the values
 given the input types, and these keyword arguments are only for overriding the
 default handling process. This should not be required by most users.
@@ -133,6 +143,7 @@ MKLPardisoIterate(; kwargs...) = PardisoJL(; solver_type = 1, kwargs...)
 PardisoJL(; nprocs::Union{Int, Nothing} = nothing,
     solver_type = nothing,
     matrix_type = nothing,
+    cache_analysis = false,
     iparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing,
     dparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing)
 ```
@@ -145,6 +156,10 @@ A generic method using MKL Pardiso. Specifying `solver_type` is required.
 
 ## Keyword Arguments
 
+Setting `cache_analysis = true` disables Pardiso's scaling and matching defaults
+and caches the result of the initial analysis phase for all further computations
+with this solver.
+
 For the definition of the keyword arguments, see the Pardiso.jl documentation.
 All values default to `nothing` and the solver internally determines the values
 given the input types, and these keyword arguments are only for overriding the
@@ -154,14 +169,16 @@ struct PardisoJL{T1, T2} <: LinearSolve.SciMLLinearSolveAlgorithm
     nprocs::Union{Int, Nothing}
     solver_type::T1
     matrix_type::T2
+    cache_analysis::Bool
     iparm::Union{Vector{Tuple{Int, Int}}, Nothing}
     dparm::Union{Vector{Tuple{Int, Int}}, Nothing}
 
     function PardisoJL(; nprocs::Union{Int, Nothing} = nothing,
             solver_type = nothing,
             matrix_type = nothing,
+            cache_analysis = false,
             iparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing,
-            dparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing)
+            dparm::Union{Vector{Tuple{Int, Float64}}, Nothing} = nothing)
         ext = Base.get_extension(@__MODULE__, :LinearSolvePardisoExt)
         if ext === nothing
             error("PardisoJL requires that Pardiso is loaded, i.e. `using Pardiso`")
@@ -170,7 +187,8 @@ struct PardisoJL{T1, T2} <: LinearSolve.SciMLLinearSolveAlgorithm
             T2 = typeof(matrix_type)
             @assert T1 <: Union{Int, Nothing, ext.Pardiso.Solver}
             @assert T2 <: Union{Int, Nothing, ext.Pardiso.MatrixType}
-            return new{T1, T2}(nprocs, solver_type, matrix_type, iparm, dparm)
+            return new{T1, T2}(
+                nprocs, solver_type, matrix_type, cache_analysis, iparm, dparm)
         end
     end
 end

diff --git a/test/pardiso/pardiso.jl b/test/pardiso/pardiso.jl
@@ -1,4 +1,4 @@
-using LinearSolve, SparseArrays, Random
+using LinearSolve, SparseArrays, Random, LinearAlgebra
 import Pardiso
 
 A1 = sparse([1.0 0 -2 3
@@ -14,19 +14,22 @@ e = ones(n)
 e2 = ones(n - 1)
 A2 = spdiagm(-1 => im * e2, 0 => lambda * e, 1 => -im * e2)
 b2 = rand(n) + im * zeros(n)
-cache_kwargs = (; verbose = true, abstol = 1e-8, reltol = 1e-8, maxiter = 30)
-
 prob2 = LinearProblem(A2, b2)
 
+cache_kwargs = (; abstol = 1e-8, reltol = 1e-8, maxiter = 30)
+
 for alg in (PardisoJL(), MKLPardisoFactorize(), MKLPardisoIterate())
-    u = solve(prob1, alg; cache_kwargs...).u
-    @test A1 * u ≈ b1
+     u = solve(prob1, alg; cache_kwargs...).u
+     @test A1 * u ≈ b1
 
-    u = solve(prob2, alg; cache_kwargs...).u
-    @test eltype(u) <: Complex
-    @test_broken A2 * u ≈ b2
+     u = solve(prob2, alg; cache_kwargs...).u
+     @test eltype(u) <: Complex
+     @test A2 * u ≈ b2
 end
 
+return
+
+
 Random.seed!(10)
 A = sprand(n, n, 0.8);
 A2 = 2.0 .* A;
@@ -53,6 +56,25 @@ sol33 = solve(linsolve)
 @test sol12.u ≈ sol32.u
 @test sol13.u ≈ sol33.u
 
+
+# Test for problem from #497
+function makeA()
+    n = 60
+    colptr = [1, 4, 7, 11, 15, 17, 22, 26, 30, 34, 38, 40, 46, 50, 54, 58, 62, 64, 70, 74, 78, 82, 86, 88, 94, 98, 102, 106, 110, 112, 118, 122, 126, 130, 134, 136, 142, 146, 150, 154, 158, 160, 166, 170, 174, 178, 182, 184, 190, 194, 198, 202, 206, 208, 214, 218, 222, 224, 226, 228, 232]
+    rowval = [1, 3, 4, 1, 2, 4, 2, 4, 9, 10, 3, 5, 11, 12, 1, 3, 2, 4, 6, 11, 12, 2, 7, 9, 10, 2, 7, 8, 10, 8, 10, 15, 16, 9, 11, 17, 18, 7, 9, 2, 8, 10, 12, 17, 18, 8, 13, 15, 16, 8, 13, 14, 16, 14, 16, 21, 22, 15, 17, 23, 24, 13, 15, 8, 14, 16, 18, 23, 24, 14, 19, 21, 22, 14, 19, 20, 22, 20, 22, 27, 28, 21, 23, 29, 30, 19, 21, 14, 20, 22, 24, 29, 30, 20, 25, 27, 28, 20, 25, 26, 28, 26, 28, 33, 34, 27, 29, 35, 36, 25, 27, 20, 26, 28, 30, 35, 36, 26, 31, 33, 34, 26, 31, 32, 34, 32, 34, 39, 40, 33, 35, 41, 42, 31, 33, 26, 32, 34, 36, 41, 42, 32, 37, 39, 40, 32, 37, 38, 40, 38, 40, 45, 46, 39, 41, 47, 48, 37, 39, 32, 38, 40, 42, 47, 48, 38, 43, 45, 46, 38, 43, 44, 46, 44, 46, 51, 52, 45, 47, 53, 54, 43, 45, 38, 44, 46, 48, 53, 54, 44, 49, 51, 52, 44, 49, 50, 52, 50, 52, 57, 58, 51, 53, 59, 60, 49, 51, 44, 50, 52, 54, 59, 60, 50, 55, 57, 58, 50, 55, 56, 58, 56, 58, 57, 59, 55, 57, 50, 56, 58, 60]
+    nzval = [-0.64, 1.0, -1.0, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -0.03510101010101016, -0.975, -1.0806825309567203, 1.0, -0.95, -0.025, 2.370597639417811, -2.3705976394178108, -11.083604432603583, -0.2770901108150896, 1.0, -0.025, -0.95, -0.3564, -0.64, 1.0, -1.0, 13.792569659442691, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -0.03510101010101016, -0.975, -1.0806825309567203, 1.0, -0.95, -0.025, 2.370597639417811, -2.3705976394178108, 10.698449178570607, -11.083604432603583, -0.2770901108150896, 1.0, -0.025, -0.95, -0.3564, -0.64, 1.0, -1.0, 13.792569659442691, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -0.03510101010101016, -0.975, -1.0806825309567203, 1.0, -0.95, -0.025, 2.370597639417811, -2.3705976394178108, 10.698449178570607, -11.083604432603583, -0.2770901108150896, 1.0, -0.025, -0.95, -0.3564, -0.64, 1.0, -1.0, 13.792569659442691, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -0.03510101010101016, -0.975, -1.0806825309567203, 1.0, -0.95, -0.025, 2.370597639417811, -2.3705976394178108, 10.698449178570607, -11.083604432603583, -0.2770901108150896, 1.0, -0.025, -0.95, -0.3564, -0.64, 1.0, -1.0, 13.792569659442691, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -0.03510101010101016, -0.975, -1.0806825309567203, 1.0, -0.95, -0.025, 2.370597639417811, -2.3705976394178108, 10.698449178570607, -11.083604432603583, -0.2770901108150896, 1.0, -0.025, -0.95, -0.3564, -0.64, 1.0, -1.0, 13.792569659442691, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -0.03510101010101016, -0.975, -1.0806825309567203, 1.0, -0.95, -0.025, 2.370597639417811, -2.3705976394178108, 10.698449178570607, -11.083604432603583, -0.2770901108150896, 1.0, -0.025, -0.95, -0.3564, -0.64, 1.0, -1.0, 13.792569659442691, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -0.03510101010101016, -0.975, -1.0806825309567203, 1.0, -0.95, -0.025, 2.370597639417811, -2.3705976394178108, 10.698449178570607, -11.083604432603583, -0.2770901108150896, 1.0, -0.025, -0.95, -0.3564, -0.64, 1.0, -1.0, 13.792569659442691, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -0.03510101010101016, -0.975, -1.0806825309567203, 1.0, -0.95, -0.025, 2.370597639417811, -2.3705976394178108, 10.698449178570607, -11.083604432603583, -0.2770901108150896, 1.0, -0.025, -0.95, -0.3564, -0.64, 1.0, -1.0, 13.792569659442691, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -0.03510101010101016, -0.975, -1.0806825309567203, 1.0, -0.95, -0.025, 2.370597639417811, -2.3705976394178108, 10.698449178570607, -11.083604432603583, -0.2770901108150896, 1.0, -0.025, -0.95, -0.3564, -0.64, 1.0, -1.0, 13.792569659442691, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -1.0806825309567203, 1.0, 2.370597639417811, -2.3705976394178108, 10.698449178570607, -11.083604432603583, -0.2770901108150896, 1.0]
+    A = SparseMatrixCSC(n, n, colptr, rowval, nzval)
+    return(A)
+end
+
+A=makeA()
+u0=fill(0.1,size(A,2))
+linprob = LinearProblem(A, A*u0)
+u = LinearSolve.solve(linprob, PardisoJL())
+@test norm(u-u0) < 1.0e-14
+
+
+
 # Testing and demonstrating Pardiso.set_iparm! for MKLPardisoSolver
 solver = Pardiso.MKLPardisoSolver()
 iparm = [