Follow matrix type branching from MKL Pardiso docs (#114)

Co-authored-by: Kristoffer <kcarlsson89@gmail.com>

Author: briochemc

Project.toml

@@ -17,7 +17,8 @@ julia = "1.6"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Random", "Printf", "Pkg"]
+test = ["StableRNGs", "Test", "Random", "Printf", "Pkg"]

src/Pardiso.jl

@@ -28,7 +28,7 @@ end
 
 MKL_LOAD_FAILED = false
 
-mkl_is_available() = (LOCAL_MKL_FOUND || MKL_jll.is_available()) && !MKL_LOAD_FAILED 
+mkl_is_available() = (LOCAL_MKL_FOUND || MKL_jll.is_available()) && !MKL_LOAD_FAILED
 
 if LinearAlgebra.BLAS.vendor() === :mkl && LinearAlgebra.BlasInt == Int64
     const MklInt = Int64
@@ -128,7 +128,7 @@ function __init__()
     elseif MKL_jll.is_available()
         libmkl_rt[] = MKL_jll.libmkl_rt_path
     end
-    
+
     if !haskey(ENV, "PARDISOLICMESSAGE")
         ENV["PARDISOLICMESSAGE"] = 1
     end
@@ -137,7 +137,7 @@ function __init__()
         @warn "MKLROOT not set, MKL Pardiso solver will not be functional"
     end
 
-    if mkl_is_available() 
+    if mkl_is_available()
         try
             libmklpardiso = Libdl.dlopen(libmkl_rt[])
             mklpardiso_f = Libdl.dlsym(libmklpardiso, "pardiso")
@@ -146,7 +146,7 @@ function __init__()
             MKL_LOAD_FAILED = true
         end
     end
-    
+
     # This is apparently needed for MKL to not get stuck on 1 thread when
     # libpardiso is loaded in the block below...
     if libmkl_rt[] !== ""
@@ -241,6 +241,85 @@ function fix_iparm!(ps::AbstractPardisoSolver, T::Symbol)
     end
 end
 
+# Copied from SparseArrays.jl but with a tweak
+# to check symmetry of the sparsity pattern
+function _is_hermsym(A::SparseMatrixCSC, check::Function)
+    m, n = size(A)
+    if m != n; return false; end
+
+    colptr = SparseArrays.getcolptr(A)
+    rowval = rowvals(A)
+    nzval = nonzeros(A)
+    tracker = copy(SparseArrays.getcolptr(A))
+    for col in axes(A,2)
+        # `tracker` is updated such that, for symmetric matrices,
+        # the loop below starts from an element at or below the
+        # diagonal element of column `col`"
+        for p = tracker[col]:colptr[col+1]-1
+            val = nzval[p]
+            row = rowval[p]
+
+            # Ignore stored zeros
+            if iszero(val)
+                continue
+            end
+
+            # If the matrix was symmetric we should have updated
+            # the tracker to start at the diagonal or below. Here
+            # we are above the diagonal so the matrix can't be symmetric.
+            if row < col
+                return false
+            end
+
+            # Diagonal element
+            if row == col
+                if !check(val, val)
+                    return false
+                end
+            else
+                offset = tracker[row]
+
+                # If the matrix is unsymmetric, there might not exist
+                # a rowval[offset]
+                if offset > length(rowval)
+                    return false
+                end
+
+                row2 = rowval[offset]
+
+                # row2 can be less than col if the tracker didn't
+                # get updated due to stored zeros in previous elements.
+                # We therefore "catch up" here while making sure that
+                # the elements are actually zero.
+                while row2 < col
+                    if _isnotzero(nzval[offset])
+                        return false
+                    end
+                    offset += 1
+                    row2 = rowval[offset]
+                    tracker[row] += 1
+                end
+
+                # Non zero A[i,j] exists but A[j,i] does not exist
+                if row2 > col
+                    return false
+                end
+
+                # A[i,j] and A[j,i] exists
+                if row2 == col
+                    if !check(val, nzval[offset])
+                        return false
+                    end
+                    tracker[row] += 1
+                end
+            end
+        end
+    end
+    return true
+end
+
+isstructurallysymmetric(A::SparseMatrixCSC) = _is_hermsym(A, (x,y) -> true)
+
 function solve!(ps::AbstractPardisoSolver, X::StridedVecOrMat{Tv},
                 A::SparseMatrixCSC{Tv,Ti}, B::StridedVecOrMat{Tv},
                 T::Symbol=:N) where {Ti, Tv <: PardisoNumTypes}
@@ -253,38 +332,72 @@ function solve!(ps::AbstractPardisoSolver, X::StridedVecOrMat{Tv},
     #   - On pos def exception, solve instead with symmetric indefinite.
     # - If complex and symmetric, solve with symmetric complex solver
     # - Else solve as unsymmetric.
-    if ishermitian(A)
-        eltype(A) <: Union{Float32, Float64} ? set_matrixtype!(ps, REAL_SYM_POSDEF) : set_matrixtype!(ps, COMPLEX_HERM_POSDEF)
-        pardisoinit(ps)
-        fix_iparm!(ps, T)
-        eltype(A) <: Union{Float32, ComplexF32} ? set_iparm!(ps, 28, 1) : nothing
-        try
+    if Tv <: Union{Float32, Float64}
+        if issymmetric(A)
+            set_matrixtype!(ps, REAL_SYM_POSDEF)
+            pardisoinit(ps)
+            fix_iparm!(ps, T)
+            Tv == Float32 && set_iparm!(ps, 28, 1)
+            try
+                pardiso(ps, X, get_matrix(ps, A, T), B)
+            catch e
+                set_phase!(ps, RELEASE_ALL)
+                pardiso(ps, X, A, B)
+                set_phase!(ps, ANALYSIS_NUM_FACT_SOLVE_REFINE)
+                if !isa(e, PardisoPosDefException)
+                    rethrow()
+                end
+                set_matrixtype!(ps, REAL_SYM_INDEF)
+                pardisoinit(ps)
+                fix_iparm!(ps, T)
+                Tv == Float32 && set_iparm!(ps, 28, 1)
+                pardiso(ps, X, get_matrix(ps, A, T), B)
+            end
+        else
+            if isstructurallysymmetric(A)
+                set_matrixtype!(ps, REAL_SYM)
+            else
+                set_matrixtype!(ps, REAL_NONSYM)
+            end
+            pardisoinit(ps)
+            fix_iparm!(ps, T)
+            Tv == Float32 && set_iparm!(ps, 28, 1)
             pardiso(ps, X, get_matrix(ps, A, T), B)
-        catch e
-            set_phase!(ps, RELEASE_ALL)
-            pardiso(ps, X, A, B)
-            set_phase!(ps, ANALYSIS_NUM_FACT_SOLVE_REFINE)
-            if !isa(e, PardisoPosDefException)
-                rethrow()
+        end
+    else # Tv <: Union{ComplexF64, ComplexF32}
+        if ishermitian(A)
+            set_matrixtype!(ps, COMPLEX_HERM_POSDEF)
+            pardisoinit(ps)
+            fix_iparm!(ps, T)
+            Tv == ComplexF32 && set_iparm!(ps, 28, 1)
+            try
+                pardiso(ps, X, get_matrix(ps, A, T), B)
+            catch e
+                set_phase!(ps, RELEASE_ALL)
+                pardiso(ps, X, A, B)
+                set_phase!(ps, ANALYSIS_NUM_FACT_SOLVE_REFINE)
+                if !isa(e, PardisoPosDefException)
+                    rethrow()
+                end
+                set_matrixtype!(ps, COMPLEX_HERM_INDEF)
+                pardisoinit(ps)
+                fix_iparm!(ps, T)
+                Tv == ComplexF32 && set_iparm!(ps, 28, 1)
+                pardiso(ps, X, get_matrix(ps, A, T), B)
+            end
+        else
+            if issymmetric(A)
+                set_matrixtype!(ps, COMPLEX_SYM)
+            elseif isstructurallysymmetric(A)
+                set_matrixtype!(ps, COMPLEX_STRUCT_SYM)
+            else
+                set_matrixtype!(ps, COMPLEX_NONSYM)
             end
-            eltype(A) <: Union{Float32, Float64} ? set_matrixtype!(ps, REAL_SYM_INDEF) : set_matrixtype!(ps, COMPLEX_HERM_INDEF)
             pardisoinit(ps)
             fix_iparm!(ps, T)
-            eltype(A) <: Union{Float32, ComplexF32} ? set_iparm!(ps, 28, 1) : nothing
+            Tv == ComplexF32 && set_iparm!(ps, 28, 1)
             pardiso(ps, X, get_matrix(ps, A, T), B)
         end
-    elseif issymmetric(A)
-        set_matrixtype!(ps, COMPLEX_SYM)
-        pardisoinit(ps)
-        fix_iparm!(ps, T)
-        eltype(A) <: Union{Float32, ComplexF32} ? set_iparm!(ps, 28, 1) : nothing
-        pardiso(ps, X, get_matrix(ps, A, T), B)
-    else
-        eltype(A) <: Union{Float32, Float64} ? set_matrixtype!(ps, REAL_NONSYM) : set_matrixtype!(ps, COMPLEX_NONSYM)
-        pardisoinit(ps)
-        fix_iparm!(ps, T)
-        eltype(A) <: Union{Float32, ComplexF32} ? set_iparm!(ps, 28, 1) : nothing
-        pardiso(ps, X, get_matrix(ps, A, T), B)
     end
 
     # Release memory, TODO: We are running the convert on IA and JA here
@@ -344,7 +457,7 @@ function pardiso(ps::AbstractPardisoSolver, X::StridedVecOrMat{Tv}, A::SparseMat
     end
 
     N = size(A, 2)
-    
+
     resize!(ps.perm, size(B, 1))
 
     NRHS = size(B, 2)

src/enums.jl

@@ -12,7 +12,7 @@
 )
 
 Base.isreal(v::MatrixType) = v in (REAL_SYM, REAL_SYM_POSDEF, REAL_SYM_INDEF, REAL_NONSYM)
-LinearAlgebra.issymmetric(v::MatrixType) = v in (REAL_SYM, REAL_SYM_POSDEF, REAL_SYM_INDEF, COMPLEX_STRUCT_SYM,
+LinearAlgebra.issymmetric(v::MatrixType) = v in (REAL_SYM_POSDEF, REAL_SYM_INDEF,
                                         COMPLEX_HERM_POSDEF, COMPLEX_HERM_INDEF, COMPLEX_SYM)
 LinearAlgebra.ishermitian(v::MatrixType) = v in (REAL_SYM_POSDEF, COMPLEX_HERM_POSDEF, COMPLEX_HERM_INDEF)
 isposornegdef(v::MatrixType) = v in (REAL_SYM_POSDEF, REAL_SYM_INDEF, COMPLEX_HERM_POSDEF, COMPLEX_HERM_INDEF)
@@ -30,7 +30,7 @@ const MATRIX_STRING = Dict{MatrixType, String}(
 )
 
 const REAL_MATRIX_TYPES = [REAL_SYM, REAL_SYM_POSDEF, REAL_SYM_INDEF, REAL_NONSYM]
-const COMPLEX_MATRIX_TYPES = [COMPLEX_STRUCT_SYM, COMPLEX_HERM_POSDEF, COMPLEX_HERM_INDEF, COMPLEX_NONSYM]
+const COMPLEX_MATRIX_TYPES = [COMPLEX_STRUCT_SYM, COMPLEX_HERM_POSDEF, COMPLEX_HERM_INDEF, COMPLEX_NONSYM, COMPLEX_SYM]
 
 # Messages
 @enum(MessageLevel::Int32,

test/runtests.jl

@@ -9,12 +9,11 @@ end
 
 using Test
 using Pardiso
+using StableRNGs
 using Random
 using SparseArrays
 using LinearAlgebra
 
-Random.seed!(1234)
-
 available_solvers = empty([Pardiso.AbstractPardisoSolver])
 if Pardiso.mkl_is_available()
     push!(available_solvers, MKLPardisoSolver)
@@ -27,6 +26,8 @@ else
     @warn "Not testing panua Pardiso solver"
 end
 
+const rng = StableRNG(1)
+
 @show Pardiso.MklInt
 
 println("Testing ", available_solvers)
@@ -38,9 +39,10 @@ supported_eltypes(ps::MKLPardisoSolver) = (Float32, ComplexF32, Float64, Complex
 @testset "solving" begin
 for pardiso_type in available_solvers
     ps = pardiso_type()
+    Random.seed!(rng, 1234)
     for T in supported_eltypes(ps)
-        A1 = sparse(rand(T, 10,10))
-        for B in (rand(T, 10, 2), view(rand(T, 10, 4), 1:10, 2:3))
+        A1 = sparse(rand(rng, T, 10,10))
+        for B in (rand(rng, T, 10, 2), view(rand(rng, T, 10, 4), 1:10, 2:3))
             X = similar(B)
             # Test unsymmetric, herm indef, herm posdef and symmetric
             for A in SparseMatrixCSC[A1, A1 + A1', A1'A1 + I, transpose(A1) + A1]
@@ -91,6 +93,7 @@ end
 
 if Pardiso.PARDISO_LOADED[]
 @testset "schur" begin
+    Random.seed!(rng, 1234)
     # reproduce example from Pardiso website
     include("schur_matrix_def.jl")
     @test norm(real(D) - real(C)*rA⁻¹*real(B) - s) < 1e-10*(8)^2
@@ -108,11 +111,11 @@ if Pardiso.PARDISO_LOADED[]
             set_matrixtype!(ps, 13)
         end
         for j ∈ 1:100
-            A = 5I + sprand(T,m,m,p)
+            A = 5I + sprand(rng,T,m,m,p)
             A⁻¹ = inv(Matrix(A))
-            B = sprand(T,m,n,p)
-            C = sprand(T,n,m,p)
-            D = 5I + sprand(T,n,n,p)
+            B = sprand(rng,T,m,n,p)
+            C = sprand(rng,T,n,m,p)
+            D = 5I + sprand(rng,T,n,n,p)
             M = [A B; C D]
 
             # test integer block specification
@@ -137,13 +140,14 @@ end # testset
 end
 
 @testset "error checks" begin
+Random.seed!(rng, 1234)
 for pardiso_type in available_solvers
 
     ps = pardiso_type()
 
-    A = sparse(rand(10,10))
-    B = rand(10, 2)
-    X = rand(10, 2)
+    A = sparse(rand(rng,10,10))
+    B = rand(rng, 10, 2)
+    X = rand(rng, 10, 2)
 
     if pardiso_type == PardisoSolver
         printstats(ps, A, B)
@@ -160,7 +164,7 @@ for pardiso_type in available_solvers
     X = zeros(12, 2)
     @test_throws DimensionMismatch solve!(ps,X, A, B)
 
-    B = rand(12, 2)
+    B = rand(rng, 12, 2)
     @test_throws DimensionMismatch solve(ps, A, B)
 end
 end # testset
@@ -201,9 +205,10 @@ for pardiso_type in available_solvers
 end
 
 @testset "pardiso" begin
+    Random.seed!(rng, 1234)
     for pardiso_type in available_solvers
-        A = sparse(rand(2,2) + im * rand(2,2))
-        b = rand(2)          + im * rand(2)
+        A = sparse(rand(rng,2,2) + im * rand(rng,2,2))
+        b = rand(rng,2)          + im * rand(rng,2)
         ps = pardiso_type()
         set_matrixtype!(ps, Pardiso.COMPLEX_NONSYM)
         x = Pardiso.solve(ps, A, b);