First steps towards cholesky

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.5.1"
+version = "0.5.2"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"

src/cholmod_cholesky.jl

@@ -0,0 +1,40 @@
+"""
+$(TYPEDEF)
+
+Default Julia LU Factorization based on umfpack.
+"""
+mutable struct ExtendableSparseCholmodCholesky{Tv, Ti} <: AbstractExtendableSparseLU{Tv,Ti}
+    A::ExtendableSparseMatrix{Tv,Ti}
+    A64
+    fact::SuiteSparse.CHOLMOD.Factor{Tv}
+    phash::UInt64
+end
+
+"""
+```
+ExtendableSparseCholmodCholesky(A)
+```
+"""
+function ExtendableSparseCholmodCholesky(A::ExtendableSparseMatrix{Tv,Ti}) where {Tv,Ti}
+    flush!(A)
+    A64=Symmetric(SparseMatrixCSC{Float64,Int64}(A.cscmatrix))
+    ExtendableSparseCholmodCholesky(A,A64,cholesky(A64),A.phash)
+end
+
+function update!(cholfact::ExtendableSparseCholmodCholesky)
+    flush!(cholfact.A)
+    if cholfact.A.phash!=cholfact.phash
+        cholfact.A64=Symmetric{SparseMatrixCSC{Float64,Int64}(A.cscmatrix)}
+        cholfact.fact=cholesky(cholfact.fact,cholfact.A64)
+        cholfact.phash=cholfact.A.phash
+    else
+        cholfact.A64.data.nzval.=cholfact.A.cscmatrix.nzval
+        cholfact.fact=cholesky!(cholfact.fact,cholfact.A64)
+    end
+    cholfact
+end
+
+
+LinearAlgebra.ldiv!(fact::ExtendableSparseCholmodCholesky, v)=fact.fact\v
+LinearAlgebra.ldiv!(u,fact::ExtendableSparseCholmodCholesky, v)=u.=fact.fact\v
+

src/extendable.jl

@@ -342,6 +342,16 @@ function LinearAlgebra.cond(A::ExtendableSparseMatrix, p::Real=2)
     return LinearAlgebra.cond(A.cscmatrix,p)
 end
 
+"""
+$(SIGNATURES)
+
+[`flush!`](@ref) and check for symmetry of cscmatrix
+"""
+function LinearAlgebra.issymmetric(A::ExtendableSparseMatrix)
+    @inbounds flush!(A)
+    return LinearAlgebra.issymmetric(A.cscmatrix)
+end
+
 
 
 """
@@ -382,3 +392,5 @@ function SparseArrays.dropzeros!(ext::ExtendableSparseMatrix)
     @inbounds flush!(ext)
     dropzeros!(ext.cscmatrix)
 end
+
+

src/factorizations.jl

@@ -102,8 +102,15 @@ end
 ```
 factorize(matrix)
 factorize(matrix; kind=:umfpack)
+factorize(matrix; kind=:default)
 ```
-Default Julia LU factorization based on UMFPACK.
+Default Julia LU factorization based on SuiteSparse.UMFPACK.
+
+```
+factorize(matrix; kind=:cholmod)
+factorize(matrix; kind=:cholesky)
+```
+Default Cholesky factorization via SuiteSparse.CHOLMOD
 
 ```
 factorize(matrix; kind=:pardiso)
@@ -145,6 +152,9 @@ Create the  [`AMGPreconditioner`](@ref) wrapping the Ruge-Stüben AMG preconditi
 function factorize(A::ExtendableSparseMatrix; kwargs...)
     opt=options(;kwargs...)
     opt[:kind]==:umfpack && return  ExtendableSparseUmfpackLU(A)
+    opt[:kind]==:default && return  ExtendableSparseUmfpackLU(A)
+    opt[:kind]==:cholmod && return  ExtendableSparseCholmodCholesky(A)
+    opt[:kind]==:cholesky && return  ExtendableSparseCholmodCholesky(A)
     opt[:kind]==:pardiso && return  PardisoLU(A,ps=Pardiso.PardisoSolver())
     opt[:kind]==:mklpardiso && return  PardisoLU(A,ps=Pardiso.MKLPardisoSolver())
     if opt[:ensurelu]
@@ -237,3 +247,4 @@ include("jacobi.jl")
 include("ilu0.jl")
 include("parallel_jacobi.jl")
 include("umfpack_lu.jl")
+include("cholmod_cholesky.jl")

src/pardiso_lu.jl

@@ -9,6 +9,23 @@ mutable struct PardisoLU{Tv, Ti} <: AbstractExtendableSparseLU{Tv,Ti}
     phash::UInt64
 end
 
+
+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
+
 """
 ```
 PardisoLU(A; ps=Pardiso.MKLPardisoSolver)
@@ -17,8 +34,8 @@ PardisoLU(A; ps=Pardiso.MKLPardisoSolver)
 function PardisoLU(A::ExtendableSparseMatrix{Tv,Ti};ps::Pardiso.AbstractPardisoSolver=Pardiso.MKLPardisoSolver()) where {Tv,Ti}
     @inbounds flush!(A)
     Acsc=A.cscmatrix
-    eltype(Acsc) == Float64 ? Pardiso.set_matrixtype!(ps, Pardiso.REAL_NONSYM) : Pardiso.set_matrixtype!(ps,  Pardiso.COMPLEX_NONSYM)
     Pardiso.pardisoinit(ps)
+    Pardiso.set_matrixtype!(ps,A)
     Pardiso.fix_iparm!(ps, :N)
     Pardiso.set_phase!(ps, Pardiso.ANALYSIS_NUM_FACT)
     Pardiso.pardiso(ps, Tv[], Acsc, Tv[])
@@ -32,6 +49,7 @@ function update!(lufact::PardisoLU{Tv,Ti}) where {Tv, Ti}
     if lufact.phash!=lufact.A.phash
         Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
         Pardiso.pardiso(ps, Tv[], Acsc, Tv[])
+        Pardiso.set_matrixtype!(ps,lufact.A)
         Pardiso.set_phase!(ps, Pardiso.ANALYSIS_NUM_FACT)
         lufact.phash=lufact.A.phash
     else