comments

Author: vpuri3

src/pardiso.jl

@@ -1,5 +1,5 @@
 
-## Paradiso
+## Pardiso
 
 import Pardiso
 
@@ -17,10 +17,10 @@ Base.@kwdef struct PardisoJL <: SciMLLinearSolveAlgorithm
 end
 
 MKLPardisoFactorize(;kwargs...) = PardisoJL(;fact_phase=Pardiso.NUM_FACT,
-                                            solve_phase=Pardiso.SOLVE_ITERATIVE_REFINE,
-                                            kwargs...)
+                                             solve_phase=Pardiso.SOLVE_ITERATIVE_REFINE,
+                                             kwargs...)
 MKLPardisoIterate(;kwargs...) = PardisoJL(;solve_phase=Pardiso.NUM_FACT_SOLVE_REFINE,
-                                          kwargs...)
+                                           kwargs...)
 
 # TODO schur complement functionality
 
@@ -50,9 +50,7 @@ function init_cacheval(alg::PardisoJL, cache::LinearCache)
     matrix_type !== nothing && Pardiso.set_matrixtype!(solver, matrix_type)
     cache.verbose && Pardiso.set_msglvl!(solver, Pardiso.MESSAGE_LEVEL_ON)
 
-    """
-    pass in vector of tuples like [(iparm::Int, key::Int) ...]
-    """
+    # pass in vector of tuples like [(iparm::Int, key::Int) ...]
     if iparm !== nothing
         for i in length(iparm)
             Pardiso.set_iparm!(solver, iparm[i]...)
@@ -66,7 +64,6 @@ function init_cacheval(alg::PardisoJL, cache::LinearCache)
     end
 
     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
@@ -76,12 +73,6 @@ function init_cacheval(alg::PardisoJL, cache::LinearCache)
         Pardiso.pardiso(solver, u, A, b)
     end
 
-    """
-    use ipram/dpram to set tolerances???
-    """
-#   abstol = cache.abstol
-#   reltol = cache.reltol
-
     return solver
 end
 

test/runtests.jl

@@ -131,22 +131,35 @@ end
 
     using Pardiso, SparseArrays
 
-    A = sparse([ 1. 0 -2  3
+    A1 = sparse([ 1. 0 -2  3
                  0  5  1  2
                 -2  1  4 -7
                  3  2 -7  5 ])
-    b = rand(4)
-    u = zero(b)
+    b1 = rand(4)
+    prob1 = LinearProblem(A1, b1)
+
+    lambda = 3
+    e = ones(n)
+    e2 = ones(n-1)
+    A2 = spdiagm(-1 => im*e2, 0 => lambda*e, 1 => -im*e2)
+    b2 = rand(n) + im * zeros(n)
+
+    prob2 = LinearProblem(A2, b2)
 
-    prob = LinearProblem(A, b)
     for alg in (
                 PardisoJL(),
                 MKLPardisoFactorize(),
                 MKLPardisoIterate(),
                )
 
-        u = solve(prob, alg; cache_kwargs...)
-        @test A * u ≈ b
+        u = solve(prob1, alg; cache_kwargs...)
+        @test A1 * u ≈ b1
+
+        # common interface doesn't support complex types
+        # https://github.com/SciML/LinearSolve.jl/issues/38
+
+#       u = solve(prob2, alg; cache_kwargs...)
+#       @test A2 * u ≈ b2
     end
 
 end