Merge branch 'main' into pardiso-vendor

j-fu · j-fu · commit 372ea1efa78c · 2024-06-05T12:29:08.000+02:00
diff --git a/ext/LinearSolvePardisoExt.jl b/ext/LinearSolvePardisoExt.jl
@@ -22,7 +22,7 @@ function LinearSolve.init_cacheval(alg::PardisoJL,
         reltol,
         verbose::Bool,
         assumptions::LinearSolve.OperatorAssumptions)
-    @unpack nprocs, solver_type, matrix_type, iparm, dparm, vendor = alg
+    @unpack nprocs, solver_type, matrix_type, cache_analysis, iparm, dparm, vendor = alg
     A = convert(AbstractMatrix, A)
 
     if isnothing(vendor)
@@ -37,8 +37,12 @@ function LinearSolve.init_cacheval(alg::PardisoJL,
     solver = if vendor == :MKL
         solver = if Pardiso.mkl_is_available()
             solver = Pardiso.MKLPardisoSolver()
+            Pardiso.pardisoinit(solver)
             nprocs !== nothing && Pardiso.set_nprocs!(solver, nprocs)
-            transposed_iparm = 2 # see https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2024-0/pardiso-iparm-parameter.html#IPARM37
+
+            # 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
         else
@@ -47,6 +51,7 @@ function LinearSolve.init_cacheval(alg::PardisoJL,
     elseif vendor == :Panua
         solver = if Pardiso.panua_is_available()
             solver = Pardiso.PardisoSolver()
+            Pardiso.pardisoinit(solver)
             solver_type !== nothing && Pardiso.set_solver!(solver, solver_type)
             
             solver
@@ -57,8 +62,6 @@ function LinearSolve.init_cacheval(alg::PardisoJL,
         error("Pardiso vendor must be either `:MKL` or `:Panua`")
     end
     
-    Pardiso.pardisoinit(solver) # default initialization
-
     if matrix_type !== nothing
         Pardiso.set_matrixtype!(solver, matrix_type)
     else
@@ -72,22 +75,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, transposed_iparm)
-
     #=
     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.
@@ -99,36 +86,58 @@ 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.
+        # 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
+
+>>>>>>> main
     return solver
 end
 
 function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::PardisoJL; kwargs...)
     @unpack A, b, u = cache
     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)
 
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(; vendor=:MKL, solver_type = 0, kwa
 ```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(; vendor=:MKL, solver_type = 1, kwarg
 ```julia
 PanuaPardisoFactorize(; 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)
 ```
@@ -145,6 +156,10 @@ A sparse factorization method using Panua Pardiso.
 
 ## 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
@@ -207,13 +222,15 @@ 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}
     vendor::Union{Symbol,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,
             vendor::Union{Symbol,Nothing}=nothing )
@@ -225,7 +242,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, vendor)
+            return new{T1, T2}(
+                nprocs, solver_type, matrix_type, cache_analysis, iparm, dparm, vendor)
         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
@@ -38,6 +38,9 @@ for alg in algs
     @test A2 * u ≈ b2
 end
 
+return
+
+
 Random.seed!(10)
 A = sprand(n, n, 0.8);
 A2 = 2.0 .* A;
@@ -64,6 +67,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 = [
