Merge pull request #58 from SciML/complex

Fix default tolerance for complex
SciML · Dec 15, 2021 · e2fb383 · e2fb383
2 parents d8054e0 + 50331e2
commit e2fb383
Show file tree

Hide file tree

Showing 3 changed files with 22 additions and 9 deletions.
diff --git a/src/common.jl b/src/common.jl
@@ -69,10 +69,15 @@ init_cacheval(alg::Union{SciMLLinearSolveAlgorithm,Nothing}, A, b, u) = nothing
 
 SciMLBase.init(prob::LinearProblem, args...; kwargs...) = SciMLBase.init(prob,nothing,args...;kwargs...)
 
+default_tol(::Type{T}) where T = √(eps(T))
+default_tol(::Type{Complex{T}}) where T = √(eps(T))
+default_tol(::Type{<:Rational}) = 0
+default_tol(::Type{<:Integer}) = 0
+
 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)),
+                        abstol=default_tol(eltype(prob.A)),
+                        reltol=default_tol(eltype(prob.A)),
                         maxiters=length(prob.b),
                         verbose=false,
                         Pl = nothing,

diff --git a/src/pardiso.jl b/src/pardiso.jl
@@ -40,7 +40,17 @@ function init_cacheval(alg::PardisoJL, cache::LinearCache)
 
     Pardiso.pardisoinit(solver) # default initialization
 
-    matrix_type !== nothing && Pardiso.set_matrixtype!(solver, matrix_type)
+    if matrix_type !== nothing
+        Pardiso.set_matrixtype!(solver, matrix_type)
+    else
+        if eltype(A) <: Real
+            Pardiso.set_matrixtype!(solver, Pardiso.REAL_NONSYM)
+        elseif eltype(A) <: Complex
+            Pardiso.set_matrixtype!(solver, Pardiso.COMPLEX_NONSYM)
+        else
+            error("Number type not supported by Pardiso")
+        end
+    end
     cache.verbose && Pardiso.set_msglvl!(solver, Pardiso.MESSAGE_LEVEL_ON)
 
     # pass in vector of tuples like [(iparm::Int, key::Int) ...]

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -186,14 +186,12 @@ end
                 MKLPardisoIterate(),
                )
 
-        u = solve(prob1, alg; cache_kwargs...)
+        u = solve(prob1, alg; cache_kwargs...).u
         @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
+        u = solve(prob2, alg; cache_kwargs...).u
+        @test eltype(u) <: Complex
+        @test_broken A2 * u ≈ b2
     end
 
 end