Bug fixes+ type fixes

* Don't check for symmetry in Pardiso
* Set proper transpose flag for Pardisoa
* Fix type handling for Factorization()

Author: j-fu

Project.toml

@@ -1,7 +1,7 @@
 name = "ExtendableSparse"
 uuid = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
 authors = ["Juergen Fuhrmann <juergen.fuhrmann@wias-berlin.de>"]
-version = "0.6.1"
+version = "0.6.2"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"

src/amg.jl

@@ -1,6 +1,6 @@
 mutable struct AMGPreconditioner{Tv, Ti} <: AbstractPreconditioner{Tv,Ti}
     A::ExtendableSparseMatrix{Tv,Ti}
-    fact
+    fact::AlgebraicMultigrid.Preconditioner
     max_levels::Int
     max_coarse::Int
     function AMGPreconditioner{Tv,Ti}(;max_levels=10, max_coarse=10) where {Tv,Ti}
@@ -13,13 +13,13 @@ end
 
 """
 ```
-AMGPreconditioner(;max_levels=10, max_coarse=10)
+AMGPreconditioner(;max_levels=10, max_coarse=10, valuetype=Float64,indextype=Int64)
 AMGPreconditioner(matrix;max_levels=10, max_coarse=10)
 ```
 
 Create the  [`AMGPreconditioner`](@ref) wrapping the Ruge-Stüben AMG preconditioner from [AlgebraicMultigrid.jl](https://github.com/JuliaLinearAlgebra/AlgebraicMultigrid.jl)
 """
-AMGPreconditioner(;kwargs...)=AMGPreconditioner{Float64,Int64}(;kwargs...)
+AMGPreconditioner(;valuetype::Type=Float64, indextype::Type=Int64,kwargs...)=AMGPreconditioner{valuetype,indextype}(;kwargs...)
 
 @eval begin
     @makefrommatrix AMGPreconditioner

src/cholmod_cholesky.jl

@@ -6,12 +6,12 @@ mutable struct CholeskyFactorization{Tv, Ti} <: AbstractLUFactorization{Tv,Ti}
 end
 
 """
-CholeskyFactorization()
+CholeskyFactorization(;valuetype=Float64, indextype=Int64)
 CholeskyFactorization(matrix)
 
 Default Cholesky factorization via cholmod.
 """
-CholeskyFactorization()=CholeskyFactorization{Float64,Int64}(nothing,nothing,0,nothing)
+CholeskyFactorization(;valuetype::Type=Float64, indextype::Type=Int64)=CholeskyFactorization{valuetype,indextype}(nothing,nothing,0,nothing)
 
 
 function update!(cholfact::CholeskyFactorization)

src/factorizations.jl

@@ -117,7 +117,7 @@ LinearAlgebra.ldiv!(fact::AbstractFactorization, v)=ldiv!(fact.fact,v)
 
 macro makefrommatrix(fact)
     return quote
-        $(esc(fact))(A::ExtendableSparseMatrix{Tv,Ti}; kwargs...) where {Tv,Ti} = factorize!($(esc(fact))(),A; kwargs...)
+        $(esc(fact))(A::ExtendableSparseMatrix{Tv,Ti}; kwargs...) where {Tv,Ti} = factorize!($(esc(fact))(;valuetype=Tv, indextype=Ti),A; kwargs...)
         $(esc(fact))(A::SparseMatrixCSC{Tv,Ti}; kwargs...) where {Tv,Ti} = $(esc(fact))(ExtendableSparseMatrix(A); kwargs...)
     end
 end

src/ilu0.jl

@@ -12,13 +12,13 @@ end
 
 """
 ```
-ILU0Preconditioner()
+ILU0Preconditioner(;valuetype=Float64,indextype=Int64)
 ILU0Preconditioner(matrix)
 ```
 
 ILU preconditioner with zero fill-in.
 """
-ILU0Preconditioner()=ILU0Preconditioner{Float64,Int64}()
+ILU0Preconditioner(;valuetype::Type=Float64, indextype::Type=Int64)=ILU0Preconditioner{valuetype,indextype}()
 
 
 function update!(precon::ILU0Preconditioner{Tv,Ti}) where {Tv,Ti}

src/ilut.jl

@@ -1,5 +1,5 @@
 mutable struct ILUTPreconditioner{Tv, Ti} <: AbstractPreconditioner{Tv,Ti}
-    A::ExtendableSparseMatrix
+    A::ExtendableSparseMatrix{Tv,Ti}
     fact::IncompleteLU.ILUFactorization{Tv,Ti}
     droptol::Float64
     function ILUTPreconditioner{Tv,Ti}(;droptol=1.0e-3) where {Tv,Ti}
@@ -11,15 +11,15 @@ end
 
 """
 ```
-ILUTPreconditioner(;droptol=1.0e-3)
+ILUTPreconditioner(;droptol=1.0e-3,valuetype=Float64, indextype=Int64)
 ILUTPreconditioner(matrix; droptol=1.0e-3)
 ```
 
 Create the [`ILUTPreconditioner`](@ref) wrapping the one 
 from [IncompleteLU.jl](https://github.com/haampie/IncompleteLU.jl)
 For using this, you need to issue `using IncompleteLU`.
 """
-ILUTPreconditioner(;kwargs...)=ILUTPreconditioner{Float64,Int64}(;kwargs...)
+ILUTPreconditioner(;valuetype::Type=Float64, indextype::Type=Int64, kwargs...)=ILUTPreconditioner{valuetype,indextype}(;kwargs...)
 
 @eval begin
     @makefrommatrix ILUTPreconditioner

src/jacobi.jl

@@ -11,13 +11,13 @@ end
 
 """
 ```
-JacobiPreconditioner()
+JacobiPreconditioner(;valuetype=Float64, indextype=Int64)
 JacobiPreconditioner(matrix)
 ```
 
 Jacobi preconditioner
 """
-JacobiPreconditioner()=JacobiPreconditioner{Float64,Int64}()
+JacobiPreconditioner(;valuetype::Type=Float64, indextype::Type=Int64)=JacobiPreconditioner{valuetype,indextype}()
 
 function update!(precon::JacobiPreconditioner{Tv,Ti}) where {Tv,Ti}
     flush!(precon.A)

src/parallel_jacobi.jl

@@ -24,13 +24,13 @@ end
 
 """
 ```
-ParallelJacobiPreconditioner()
+ParallelJacobiPreconditioner(;valuetype=Float64, indextype=Int64)
 ParallelJacobiPreconditioner(matrix)
 ```
 
 Parallel Jacobi preconditioner
 """
-ParallelJacobiPreconditioner()=ParallelJacobiPreconditioner{Float64,Int64}()
+ParallelJacobiPreconditioner(;valuetype::Type=Float64, indextype::Type=Int64)=ParallelJacobiPreconditioner{valuetype,indextype}()
 
 function  LinearAlgebra.ldiv!(u::AbstractArray{T,1} where T, precon::ParallelJacobiPreconditioner, v::AbstractArray{T,1} where T)
     invdiag=precon.invdiag

src/pardiso_lu.jl

@@ -11,7 +11,7 @@ PardisoLU{Tv,Ti}() where {Tv,Ti} =PardisoLU{Tv,Ti}(nothing,Pardiso.PardisoSolver
 
 """
 ```
-PardisoLU()
+PardisoLU(;valuetype=Float64, indextype=Int64)
 PardisoLU(matrix)
 ```
 
@@ -27,20 +27,19 @@ plu=PardisoLU()
 Pardiso.set_iparm!(plu.ps,5,13.0)
 ```
 """
-PardisoLU()=PardisoLU{Float64,Int64}(nothing,Pardiso.PardisoSolver(),0)
+PardisoLU(;valuetype::Type=Float64, indextype::Type=Int64)=PardisoLU{valuetype,indextype}(nothing,Pardiso.PardisoSolver(),0)
 
 
 mutable struct MKLPardisoLU{Tv, Ti} <: AbstractPardisoLU{Tv,Ti}
     A::Union{ExtendableSparseMatrix{Tv,Ti},Nothing}
     ps::Pardiso.MKLPardisoSolver
     phash::UInt64
 end
-MKLPardisoLU{Tv,Ti}() where {Tv,Ti} = MKLPardisoLU{Tv,Ti}(nothing,Pardiso.MKLPardisoSolver(),0)
 
 
 """
 ```
-MKLPardisoLU()
+MKLPardisoLU(;valuetype=Float64, indextype=Int64)
 MKLPardisoLU(matrix)
 ```
 
@@ -56,40 +55,29 @@ plu=MKLPardisoLU()
 Pardiso.set_iparm!(plu.ps,5,13.0)
 ```
 """
-MKLPardisoLU()=MKLPardisoLU{Float64,Int64}(nothing,Pardiso.MKLPardisoSolver(),0)
-
-
-function Pardiso.set_matrixtype!(ps, A::ExtendableSparseMatrix)
-    Acsc=A.cscmatrix
-    
-    if eltype(Acsc)==Float64 && issymmetric(Acsc)
-        Pardiso.set_matrixtype!(ps, Pardiso.REAL_SYM)
-    elseif eltype(Acsc)==Float64
-        Pardiso.set_matrixtype!(ps, Pardiso.REAL_NONSYM)
-    elseif eltype(Acsc)==Complex64 && ishermitian(Acsc)
-        Pardiso.set_matrixtype!(ps, Pardiso.COMPLEX_HERM_INDEF)
-    elseif eltype(Acsc)==Complex64
-        Pardiso.set_matrixtype!(ps, Pardiso.COMPLEX_NONYSYM)
-    else
-        error("unable to detect matrix type")
-    end
-end
+MKLPardisoLU(;valuetype::Type=Float64, indextype::Type=Int64)=MKLPardisoLU{valuetype,indextype}(nothing,Pardiso.MKLPardisoSolver(),0)
 
 
 function update!(lufact::AbstractPardisoLU{Tv,Ti}) where {Tv, Ti}
     ps=lufact.ps
     flush!(lufact.A)
     Acsc=lufact.A.cscmatrix
-    if lufact.phash!=lufact.A.phash
+    if lufact.phash!=lufact.A.phash 
         Pardiso.pardisoinit(ps)
         Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
         Pardiso.pardiso(ps, Tv[], Acsc, Tv[])
-        Pardiso.set_matrixtype!(ps,lufact.A)
+
+        if Tv<:Complex
+            Pardiso.set_matrixtype!(ps,Pardiso.COMPLEX_NONSYM)
+        else
+            Pardiso.set_matrixtype!(ps,Pardiso.REAL_NONSYM)
+        end
         Pardiso.set_phase!(ps, Pardiso.ANALYSIS_NUM_FACT)
         lufact.phash=lufact.A.phash
     else
         Pardiso.set_phase!(ps, Pardiso.NUM_FACT)
     end
+    Pardiso.fix_iparm!(ps, :N)
     Pardiso.pardiso(ps, Tv[], Acsc, Tv[])
     lufact
 end

src/umfpack_lu.jl

@@ -7,13 +7,13 @@ end
 
 """
 ```
-LUFactorization()
+LUFactorization(;valuetype=Float64, indextype=Int64)
 LUFactorization(matrix)
 ```
         
 Default Julia LU Factorization based on umfpack.
 """
-LUFactorization()=LUFactorization{Float64,Int64}(nothing,nothing,0)
+LUFactorization(;valuetype::Type=Float64,indextype::Type=Int64)=LUFactorization{valuetype,indextype}(nothing,nothing,0)
 
 
 function update!(lufact::LUFactorization)