Merge pull request #88 from SciML/pardiso

Fix Pardiso precaching form

Author: ChrisRackauckas

Project.toml

@@ -38,8 +38,9 @@ julia = "1.6"
 [extras]
 Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Pardiso", "Pkg", "SafeTestsets"]
+test = ["Test", "Pardiso", "Pkg", "Random", "SafeTestsets"]

src/pardiso.jl

@@ -6,17 +6,15 @@ Base.@kwdef struct PardisoJL <: SciMLLinearSolveAlgorithm
     dparm::Union{Vector{Tuple{Int,Int}}, Nothing} = nothing
 end
 
-MKLPardisoFactorize(;kwargs...) = PardisoJL(;kwargs...)
-MKLPardisoIterate(;kwargs...) = PardisoJL(;kwargs...)
+MKLPardisoFactorize(;kwargs...) = PardisoJL(;solver_type = 0,kwargs...)
+MKLPardisoIterate(;kwargs...) = PardisoJL(;solver_type = 1,kwargs...)
+needs_concrete_A(alg::PardisoJL) = true
 
 # TODO schur complement functionality
 
 function init_cacheval(alg::PardisoJL, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)
     @unpack nprocs, solver_type, matrix_type, iparm, dparm = alg
-
-    if A isa DiffEqArrayOperator
-        A = A.A
-    end
+    A = convert(AbstractMatrix,A)
 
     solver =
     if Pardiso.PARDISO_LOADED[]
@@ -59,6 +57,9 @@ function init_cacheval(alg::PardisoJL, A, b, u, Pl, Pr, maxiters, abstol, reltol
         end
     end
 
+    # Make sure to say it's transposed because its CSC not CSR
+    Pardiso.set_iparm!(solver,12, 1)
+
     Pardiso.set_phase!(solver, Pardiso.ANALYSIS)
     Pardiso.pardiso(solver, u, A, b)
 
@@ -67,23 +68,18 @@ end
 
 function SciMLBase.solve(cache::LinearCache, alg::PardisoJL; kwargs...)
     @unpack A, b, u = cache
-    if A isa DiffEqArrayOperator
-        A = A.A
-    end
+    A = convert(AbstractMatrix,A)
 
     if cache.isfresh
         Pardiso.set_phase!(cache.cacheval, Pardiso.NUM_FACT)
-        Pardiso.pardiso(cache.cacheval, cache.u, cache.A, cache.b)
+        Pardiso.pardiso(cache.cacheval, u, A, b)
     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
 
-needsconcreteA(alg::PardisoJL) = true
-
 # Add finalizer to release memory
 # Pardiso.set_phase!(cache.cacheval, Pardiso.RELEASE_ALL)
 

test/basictests.jl

@@ -1,5 +1,6 @@
 using LinearSolve, LinearAlgebra, SparseArrays
 using Test
+import Random
 
 n = 8
 A = Matrix(I,n,n)
@@ -232,4 +233,41 @@ end
     end
 end
 
+@testset "Sparse Precaching" begin
+    n = 4
+    Random.seed!(10)
+    A = sprand(n,n,0.8); A2 = 2.0 .* A
+    b1 = rand(n); b2 = rand(n)
+
+    prob = LinearProblem(copy(A), copy(b1))
+    linsolve = init(prob,UMFPACKFactorization())
+    sol11 = solve(linsolve)
+    linsolve = LinearSolve.set_b(sol11.cache,copy(b2))
+    sol12 = solve(linsolve)
+    linsolve = LinearSolve.set_A(sol12.cache,copy(A2))
+    sol13 = solve(linsolve)
+
+    prob = LinearProblem(copy(A), copy(b1))
+    linsolve = init(prob,KLUFactorization())
+    sol21 = solve(linsolve)
+    linsolve = LinearSolve.set_b(sol21.cache,copy(b2))
+    sol22 = solve(linsolve)
+    linsolve = LinearSolve.set_A(sol22.cache,copy(A2))
+    sol23 = solve(linsolve)
+
+    linsolve = init(prob,MKLPardisoFactorize())
+    sol31 = solve(linsolve)
+    linsolve = LinearSolve.set_b(sol31.cache,copy(b2))
+    sol32 = solve(linsolve)
+    linsolve = LinearSolve.set_A(sol32.cache,copy(A2))
+    sol33 = solve(linsolve)
+
+    @test sol11.u ≈ sol21.u
+    @test sol11.u ≈ sol31.u
+    @test sol12.u ≈ sol22.u
+    @test sol12.u ≈ sol32.u
+    @test sol13.u ≈ sol23.u
+    @test sol13.u ≈ sol33.u
+end
+
 end # testset