one test passing.

Author: vpuri3

Project.toml

@@ -33,6 +33,7 @@ julia = "1.6"
 DiffEqProblemLibrary = "a077e3f3-b75c-5d7f-a0c6-6bc4c8ec64a9"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
 
 [targets]
-test = ["Test", "OrdinaryDiffEq", "DiffEqProblemLibrary"]
+test = ["Test", "OrdinaryDiffEq", "DiffEqProblemLibrary", "Pardiso"]

src/common.jl

@@ -56,9 +56,9 @@ SciMLBase.init(prob::LinearProblem, args...; kwargs...) = SciMLBase.init(prob,no
 
 function SciMLBase.init(prob::LinearProblem, alg::Union{SciMLLinearSolveAlgorithm,Nothing}, args...;
                         alias_A = false, alias_b = false,
-                        abstol=√eps(eltype(prob.A)),
-                        reltol=√eps(eltype(prob.A)),
-                        maxiters=length(prob.b),
+                        abstol=√eps(eltype(prob.A)), # TODO handle when eps()
+                        reltol=√eps(eltype(prob.A)), # not defined on eltype
+                        maxiters=length(prob.b),     # like Int of Complex
                         verbose=false,
                         Pl = nothing,
                         Pr = nothing,

src/pardiso.jl

@@ -5,41 +5,35 @@ import Pardiso
 
 export PardisoJL
 
-struct PardisoJL{A} <: SciMLLinearSolveAlgorithm
-    nthreads::Union{Int, Nothing}
-    solver_type::Union{Int, Pardiso.Solver, Nothing}
-    matrix_type::Union{Int, Pardiso.MatrixType, Nothing}
-    solve_phase::Union{Int, Pardiso.Phase, Nothing}
-    release_phase::Union{Int, Nothing}
-    iparm::Union{A, Nothing}
-    dparm::Union{A, Nothing}
-end
-
-function PardisoJL(solver_type=Pardiso.nothing,)
-
-    return PardisoJL(nthreads, solver_type, matrix_type, solve_phase,
-                     release_phase, iparm, dparm)
+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
+    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
 
 function init_cacheval(alg::PardisoJL, cache::LinearCache)
-    @unpack nthreads, solver_type, matrix_type, iparm, dparm = alg
+    @unpack nprocs, solver_type, matrix_type, iparm, dparm = alg
 
-    solver = Pardiso.PARDISO_LOADED[] ? PardisoSolver() : MKLPardisoSolver()
+    solver = Pardiso.PARDISO_LOADED[] ? Pardiso.PardisoSolver() : Pardiso.MKLPardisoSolver()
 
     Pardiso.pardisoinit(solver) # default initialization
 
-    nthreads    !== nothing && Pardiso.set_nprocs!(ps, nthreads)
+    nprocs      !== nothing && Pardiso.set_nprocs!(ps, nprocs)
     solver_type !== nothing && Pardiso.set_solver!(solver, key)
     matrix_type !== nothing && Pardiso.set_matrixtype!(solver, matrix_type)
     cache.verbose && Pardiso.set_msglvl!(solver, Pardiso.MESSAGE_LEVEL_ON)
 
-    iparm !== nothing && begin # pass in vector of tuples like [(iparm, key)]
+    if iparm !== nothing # pass in vector of tuples like [(iparm, key)]
         for i in length(iparm)
             Pardiso.set_iparm!(solver, iparm[i]...)
         end
     end
 
-    dparm !== nothing && begin
+    if dparm !== nothing
         for i in length(dparm)
             Pardiso.set_dparm!(solver, dparm[i]...)
         end
@@ -49,7 +43,10 @@ function init_cacheval(alg::PardisoJL, cache::LinearCache)
 end
 
 function SciMLBase.solve(cache::LinearCache, alg::PardisoJL; kwargs...)
-    @unpack A, b, u, cacheval = cache
+    @unpack A, b, u = cache
+    if A isa DiffEqArrayOperator
+        A = A.A
+    end
 
     if cache.isfresh
         solver = init_cacheval(alg, cache)
@@ -58,15 +55,15 @@ function SciMLBase.solve(cache::LinearCache, alg::PardisoJL; kwargs...)
 
     abstol = cache.abstol
     reltol = cache.reltol
-    kwargs = (abstol=abstol, reltol=reltol, alg.kwargs...)
+    kwargs = (abstol=abstol, reltol=reltol)
 
     """
     figure out whatever phase is. should set_phase call be in init_cacheval?
     can we use phase to store factorization in cache?
     """
-    Pardiso.set_phase!(cacheval, alg.solve_phase)
-    Pardiso.solve!(cacheval, u, A, b)
-    Pardiso.set_phase!(cacheval, alg.release_phase) # is this necessary?
+    alg.solve_phase !== nothing && Pardiso.set_phase!(cacheval, alg.solve_phase)
+    Pardiso.solve!(cache.cacheval, u, A, b)
+    alg.release_phase !== nothing && Pardiso.set_phase!(cacheval, alg.release_phase) # is this necessary?
 
     return cache.u
 end

test/runtests.jl

@@ -126,6 +126,25 @@ end
     end
 end
 
+@testset "PardisoJL" begin
+    #@test_broken alg = PardisoJL()
+
+    using Pardiso, SparseArrays
+    verbose = true
+
+    A = sparse([ 1. 0 -2  3
+                 0  5  1  2
+                -2  1  4 -7
+                 3  2 -7  5 ])
+    b = rand(4)
+
+    prob = LinearProblem(A, b)
+    alg = PardisoJL()
+
+    u = solve(prob, alg; verbose=true)
+
+end
+
 @testset "Preconditioners" begin
     @testset "scaling_preconditioner" begin
         s = rand()
@@ -174,7 +193,4 @@ end
     end
 end
 
-@testset "PardisoJL" begin
-end
-
 end # testset