upadtes

Author: vpuri3

src/pardiso.jl

@@ -3,46 +3,57 @@
 
 import Pardiso
 
-export PardisoJL, PardisoJLFactorize, PardisoJLIterate
+export PardisoJL, MKLPardisoFactorize, MKLPardisoIterate
 
 Base.@kwdef struct PardisoJL <: SciMLLinearSolveAlgorithm
     nprocs::Union{Int, Nothing} = nothing
     solver_type::Union{Int, Pardiso.Solver, Nothing} = nothing
     matrix_type::Union{Int, Pardiso.MatrixType, Nothing} = nothing
+    fact_phase::Union{Int, Pardiso.Phase, Nothing} = nothing
     solve_phase::Union{Int, Pardiso.Phase, Nothing} = nothing
     release_phase::Union{Int, Nothing} = nothing
     iparm::Union{Vector{Tuple{Int,Int}}, Nothing} = nothing
     dparm::Union{Vector{Tuple{Int,Int}}, Nothing} = nothing
 end
 
-PardisoJLFactorize(;kwargs...) = PardisoJL(;solver_type=0,
-                                         solve_phase=Pardiso.NUM_FACT,
-                                         kwargs...)
-PardisoJLIterate(;kwargs...) = PardisoJL(;solver_type=1,
-                                         solve_phase=Pardiso.SOLVE_ITERATIVE_REFINE,
-                                         kwargs...)
+MKLPardisoFactorize(;kwargs...) = PardisoJL(;fact_phase=Pardiso.NUM_FACT,
+                                            solve_phase=Pardiso.SOLVE_ITERATIVE_REFINE,
+                                            kwargs...)
+MKLPardisoIterate(;kwargs...) = PardisoJL(;solve_phase=Pardiso.NUM_FACT_SOLVE_REFINE,
+                                          kwargs...)
 
 # TODO schur complement functionality
 
 function init_cacheval(alg::PardisoJL, cache::LinearCache)
-    @unpack nprocs, solver_type, matrix_type, iparm, dparm = alg
+    @unpack nprocs, solver_type, matrix_type, fact_phase, solve_phase, iparm, dparm = alg
+    @unpack A, b, u = cache
+
+    if A isa DiffEqArrayOperator
+        A = A.A
+    end
 
     solver =
     if Pardiso.PARDISO_LOADED[]
-        Pardiso.PardisoSolver()
+        solver = Pardiso.PardisoSolver()
         solver_type !== nothing && Pardiso.set_solver!(solver, solver_type)
+
+        solver
     else
-        Pardiso.MKLPardisoSolver()
+        solver = Pardiso.MKLPardisoSolver()
         nprocs !== nothing && Pardiso.set_nprocs!(solver, nprocs)
+
+        solver
     end
 
     Pardiso.pardisoinit(solver) # default initialization
 
-    @show solver
     matrix_type !== nothing && Pardiso.set_matrixtype!(solver, matrix_type)
     cache.verbose && Pardiso.set_msglvl!(solver, Pardiso.MESSAGE_LEVEL_ON)
 
-    if iparm !== nothing # pass in vector of tuples like [(iparm, key) ...]
+    """
+    pass in vector of tuples like [(iparm::Int, key::Int) ...]
+    """
+    if iparm !== nothing
         for i in length(iparm)
             Pardiso.set_iparm!(solver, iparm[i]...)
         end
@@ -54,7 +65,21 @@ function init_cacheval(alg::PardisoJL, cache::LinearCache)
         end
     end
 
-    Pardiso.set_phase!(cacheval, Pardiso.ANALYSIS)
+    if (fact_phase !== nothing) | (solve_phase !== nothing)
+        # ensure phase is being changed afterwards?
+        Pardiso.set_phase!(solver, Pardiso.ANALYSIS)
+        Pardiso.pardiso(solver, u, A, b)
+    end
+
+    if fact_phase !== nothing
+        Pardiso.set_phase!(solver, fact_phase)
+        Pardiso.pardiso(solver, u, A, b)
+    end
+
+    #   ipram/dpram
+#   abstol = cache.abstol
+#   reltol = cache.reltol
+#   kwargs = (abstol=abstol, reltol=reltol)
 
     return solver
 end
@@ -70,13 +95,9 @@ function SciMLBase.solve(cache::LinearCache, alg::PardisoJL; kwargs...)
         cache = set_cacheval(cache, solver)
     end
 
-    abstol = cache.abstol
-    reltol = cache.reltol
-    kwargs = (abstol=abstol, reltol=reltol)
-
-    alg.solve_phase !== nothing && Pardiso.set_phase!(cacheval, alg.solve_phase)
+    alg.solve_phase !== nothing && Pardiso.set_phase!(cache.cacheval, alg.solve_phase)
     Pardiso.pardiso(cache.cacheval, u, A, b)
-    alg.release_phase !== nothing && Pardiso.set_phase!(cacheval, alg.release_phase)
+    alg.release_phase !== nothing && Pardiso.set_phase!(cache.cacheval, alg.release_phase)
 
     return cache.u
 end

test/runtests.jl

@@ -127,10 +127,9 @@ end
 end
 
 @testset "PardisoJL" begin
-    @test_broken alg = PardisoJL()
+    @test_throws UndefVarError alg = PardisoJL()
 
     using Pardiso, SparseArrays
-    verbose = true
 
     A = sparse([ 1. 0 -2  3
                  0  5  1  2
@@ -140,12 +139,15 @@ end
     u = zero(b)
 
     prob = LinearProblem(A, b)
-    for alg in (PardisoJL(),
-                PardisoJLFactorize(),
-                PardisoJLIterate(), # not with MKLPardisoSolver
+    for alg in (
+                PardisoJL(),
+                MKLPardisoFactorize(),
+                MKLPardisoIterate(),
                )
 
         u = solve(prob, alg; verbose=true)
+
+        @test A * u ≈ b
     end
 
 end