Exposing single precision support from MKL (#108)

* - Exposing single precision support from MKL
 - Adding check for square matrix
 - Changing Project.lic to Panua.lic and name of of file form "project_pardiso.jl" to "panua_pardiso.jl"

* Update runtests.jl

Keeping tests running for 1.6 by always computing the reference solution using Float64/CompelxF64.

* Update runtests.jl

Converting results back to Float32/ComplexF32.

* Adding a diagonal to the `herm posdef` to make it well-condtioned. Testing that the solution actually works rather than comparing with SparseArrays

* Adding exception to catch that iparm[28]=1 for Float32/ComplexF32 and MKLPardisoSolver.

Author: mipals

src/Pardiso.jl

@@ -64,7 +64,7 @@ end
 Base.showerror(io::IO, e::Union{PardisoException,PardisoPosDefException}) = print(io, e.info);
 
 
-const PardisoNumTypes = Union{Float64,ComplexF64}
+const PardisoNumTypes = Union{Float32, ComplexF32, Float64, ComplexF64}
 
 abstract type AbstractPardisoSolver end
 
@@ -173,7 +173,7 @@ function __init__()
     end
 end
 include("enums.jl")
-include("project_pardiso.jl")
+include("panua_pardiso.jl")
 include("mkl_pardiso.jl")
 
 # Getters and setters
@@ -218,9 +218,9 @@ end
 
 function solve(ps::AbstractPardisoSolver, A::SparseMatrixCSC{Tv,Ti},
                B::StridedVecOrMat{Tv}, T::Symbol=:N) where {Ti, Tv <: PardisoNumTypes}
-  X = copy(B)
-  solve!(ps, X, A, B, T)
-  return X
+    X = copy(B)
+    solve!(ps, X, A, B, T)
+    return X
 end
 
 function fix_iparm!(ps::AbstractPardisoSolver, T::Symbol)
@@ -244,6 +244,7 @@ end
 function solve!(ps::AbstractPardisoSolver, X::StridedVecOrMat{Tv},
                 A::SparseMatrixCSC{Tv,Ti}, B::StridedVecOrMat{Tv},
                 T::Symbol=:N) where {Ti, Tv <: PardisoNumTypes}
+    LinearAlgebra.checksquare(A)
     set_phase!(ps, ANALYSIS_NUM_FACT_SOLVE_REFINE)
 
     # This is the heuristics for choosing what matrix type to use
@@ -253,9 +254,10 @@ function solve!(ps::AbstractPardisoSolver, X::StridedVecOrMat{Tv},
     # - If complex and symmetric, solve with symmetric complex solver
     # - Else solve as unsymmetric.
     if ishermitian(A)
-        eltype(A) == Float64 ? set_matrixtype!(ps, REAL_SYM_POSDEF) : set_matrixtype!(ps, COMPLEX_HERM_POSDEF)
+        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
             pardiso(ps, X, get_matrix(ps, A, T), B)
         catch e
@@ -265,20 +267,23 @@ function solve!(ps::AbstractPardisoSolver, X::StridedVecOrMat{Tv},
             if !isa(e, PardisoPosDefException)
                 rethrow()
             end
-            eltype(A) == Float64 ? set_matrixtype!(ps, REAL_SYM_INDEF) : set_matrixtype!(ps, COMPLEX_HERM_INDEF)
+            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
             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) == Float64 ? set_matrixtype!(ps, REAL_NONSYM) : set_matrixtype!(ps, COMPLEX_NONSYM)
+        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
 
@@ -333,6 +338,11 @@ function pardiso(ps::AbstractPardisoSolver, X::StridedVecOrMat{Tv}, A::SparseMat
                                     "has a complex matrix type set: $(get_matrixtype(ps))")))
     end
 
+    if Tv <: Union{Float32, ComplexF32} && typeof(ps) <: MKLPardisoSolver && ps.iparm[28] != 1
+        throw(ErrorException(string("input matrix is Float32/ComplexF32 while MKLPardisoSolver ",
+                                    "have iparm[28]=$(ps.iparm[28]) rather than 1.")))
+    end
+
     N = size(A, 2)
 
     resize!(ps.perm, size(B, 1))
@@ -435,6 +445,7 @@ function pardisogetschur(ps::AbstractPardisoSolver)
 end
 
 function dim_check(X, A, B)
+    LinearAlgebra.checksquare(A)
     size(X) == size(B) || throw(DimensionMismatch(string("solution has $(size(X)), ",
                                                          "RHS has size as $(size(B)).")))
     size(A, 1) == size(B, 1) || throw(DimensionMismatch(string("matrix has $(size(A,1)) ",

src/project_pardiso.jl renamed to src/panua_pardiso.jl

@@ -9,6 +9,8 @@ mutable struct PardisoSolver <: AbstractPardisoSolver
     maxfct::Int32
     mnum::Int32
     perm::Vector{Int32}
+    colptr::Vector{Int32}
+    rowval::Vector{Int32}
 end
 
 function PardisoSolver()
@@ -34,9 +36,11 @@ function PardisoSolver()
     mnum = 1
     maxfct = 1
     perm = Int32[]
+    colptr = Int32[]
+    rowval = Int32[]
 
     ps = PardisoSolver(pt, iparm, dparm, mtype, solver,
-                  phase, msglvl, maxfct, mnum, perm)
+                  phase, msglvl, maxfct, mnum, perm, colptr, rowval)
 
     return ps
 end
@@ -77,9 +81,14 @@ end
 
 @inline function ccall_pardiso(ps::PardisoSolver, N::Integer, nzval::Vector{Tv},
                                    colptr, rowval, NRHS::Integer, B::StridedVecOrMat{Tv}, X::StridedVecOrMat{Tv}) where {Tv}
+    (Tv == Float32 || Tv == ComplexF32) && throw(ArgumentError("Single precision input matrix only supported by MKL."))
+
     N = Int32(N)
-    colptr = convert(Vector{Int32}, colptr)
-    rowval = convert(Vector{Int32}, rowval)
+    # Save new colptr and rowvals if a new analysis phase is run
+    if ps.phase in [ANALYSIS, ANALYSIS_NUM_FACT, ANALYSIS_NUM_FACT_SOLVE_REFINE]
+        ps.colptr = convert(Vector{Int32}, colptr)
+        ps.rowval = convert(Vector{Int32}, rowval)
+    end
     resize!(ps.perm, size(B, 1))
     NRHS = Int32(NRHS)
 
@@ -90,7 +99,7 @@ end
            Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Ptr{Tv}, Ptr{Tv},
            Ptr{Int32}, Ptr{Float64}),
           ps.pt, Ref(ps.maxfct), Ref(Int32(ps.mnum)), Ref(Int32(ps.mtype)), Ref(Int32(ps.phase)),
-          Ref(N), nzval, colptr, rowval, ps.perm,
+          Ref(N), nzval, ps.colptr, ps.rowval, ps.perm,
           Ref(NRHS), ps.iparm, Ref(Int32(ps.msglvl)), B, X,
           ERR, ps.dparm)
     check_error(ps, ERR[])
@@ -99,6 +108,8 @@ end
 
 @inline function ccall_pardiso_get_schur(ps::PardisoSolver, S::Vector{Tv},
                                    IS::Vector{Int32}, JS::Vector{Int32}) where Tv
+    (Tv == Float32 || Tv == ComplexF32) && throw(ArgumentError("Single precision input matrix only supported by MKL."))
+
     ccall(pardiso_get_schur_f[], Cvoid,
           (Ptr{Int}, Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Ptr{Tv},
            Ptr{Int32}, Ptr{Int32}),
@@ -112,8 +123,8 @@ function printstats(ps::PardisoSolver, A::SparseMatrixCSC{Tv, Ti},
                     B::StridedVecOrMat{Tv}) where {Ti,Tv <: PardisoNumTypes}
     N = Int32(size(A, 2))
     AA = A.nzval
-    IA = convert(Vector{Int32}, A.colptr)
-    JA = convert(Vector{Int32}, A.rowval)
+    IA = ps.colptr
+    JA = ps.rowval
     NRHS = Int32(size(B, 2))
     ERR = Ref{Int32}(0)
     if Tv <: Complex
@@ -181,7 +192,7 @@ function check_error(ps::PardisoSolver, err::Integer)
     err != -6  || throw(PardisoException("Preordering failed (matrix types 11, 13 only)."))
     err != -7  || throw(PardisoException("Diagonal matrix problem."))
     err != -8  || throw(PardisoException("32-bit integer overflow problem."))
-    err != -10 || throw(PardisoException("No license file pardiso.lic found."))
+    err != -10 || throw(PardisoException("No license file panua.lic found."))
     err != -11 || throw(PardisoException("License is expired."))
     err != -12 || throw(PardisoException("Wrong username or hostname."))
     err != -100|| throw(PardisoException("Reached maximum number of Krylov-subspace iteration in iterative solver."))

test/runtests.jl

@@ -24,40 +24,43 @@ end
 if Pardiso.PARDISO_LOADED[]
     push!(available_solvers, PardisoSolver)
 else
-    @warn "Not testing project Pardiso solver"
+    @warn "Not testing panua Pardiso solver"
 end
 
 @show Pardiso.MklInt
 
 println("Testing ", available_solvers)
 
+supported_eltypes(ps::PardisoSolver) = (Float64, ComplexF64)
+supported_eltypes(ps::MKLPardisoSolver) = (Float32, ComplexF32, Float64, ComplexF64)
+
 # Test solver + for real and complex data
 @testset "solving" begin
 for pardiso_type in available_solvers
     ps = pardiso_type()
-    for T in (Float64, ComplexF64)
+    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))
             X = similar(B)
             # Test unsymmetric, herm indef, herm posdef and symmetric
-            for A in SparseMatrixCSC[A1, A1 + A1', A1'A1, transpose(A1) + A1]
+            for A in SparseMatrixCSC[A1, A1 + A1', A1'A1 + I, transpose(A1) + A1]
                 solve!(ps, X, A, B)
-                @test X ≈ A\Matrix(B)
+                @test A*X ≈ B
 
                 X = solve(ps, A, B)
-                @test X ≈ A\Matrix(B)
+                @test A*X ≈ B
 
                 solve!(ps, X, A, B, :C)
-                @test X ≈ A'\Matrix(B)
+                @test A'*X ≈ B
 
                 X = solve(ps, A, B, :C)
-                @test X ≈ A'\Matrix(B)
+                @test A'*X ≈ B
 
                 solve!(ps, X, A, B, :T)
-                @test X ≈ copy(transpose(A))\Matrix(B)
+                @test transpose(A)*X ≈ B
 
                 X = solve(ps, A, B, :T)
-                @test X ≈ copy(transpose(A))\Matrix(B)
+                @test transpose(A)*X ≈ B
             end
         end
     end