Add more functionality (#22)

- Addition, multiplication etc. for ExtendableSparse
- ILUZero.jl as dependency for ILUZeroPreconditioner
- fix ILU0Preconditioner backsubstitution for unsymmetric matrices

Author: j-fu

Project.toml

@@ -1,10 +1,11 @@
 name = "ExtendableSparse"
 uuid = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
 authors = ["Juergen Fuhrmann <juergen.fuhrmann@wias-berlin.de>"]
-version = "0.9.2"
+version = "0.9.3"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+ILUZero = "88f59080-6952-5380-9ea5-54057fb9a43f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
@@ -15,9 +16,10 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 DocStringExtensions = "0.8.0,0.9"
+ILUZero= "0.2"
 Requires = "1.1.3"
+Sparspak = "0.3.0"
 julia = "1.6"
-Sparspak= "0.3.0"
 
 [extras]
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"

docs/src/iter.md

@@ -124,6 +124,6 @@ nm1,nm2
 
 ```@autodocs
 Modules = [ExtendableSparse]
-Pages = ["jacobi.jl","ilu0.jl","parallel_jacobi.jl","ilut.jl","amg.jl"]
+Pages = ["jacobi.jl","parallel_jacobi.jl","iluzero.jl","ilu0.jl","ilut.jl","amg.jl"]
 ```
 

src/ExtendableSparse.jl

@@ -2,6 +2,7 @@ module ExtendableSparse
 using SparseArrays
 using LinearAlgebra
 using Sparspak
+using ILUZero
 
 
 # Define our own constant here in order to be able to
@@ -28,19 +29,16 @@ export SparseMatrixLNK,ExtendableSparseMatrix,flush!,nnz, updateindex!, rawupdat
 
 
 include("factorizations.jl")
-export JacobiPreconditioner, ILU0Preconditioner, ParallelJacobiPreconditioner, ParallelILU0Preconditioner, reorderlinsys
+
+export JacobiPreconditioner, ILU0Preconditioner,ILUZeroPreconditioner, ParallelJacobiPreconditioner, ParallelILU0Preconditioner, reorderlinsys
 export AbstractFactorization,LUFactorization, CholeskyFactorization
 export issolver
 export factorize!, update!
-export ILUTPreconditioner, AMGPreconditioner
-export PardisoLU, MKLPardisoLU,SparspakLU
 
 include("simple_iteration.jl")
 export simple,simple!
 
 
-
-
 include("sprand.jl")
 export sprand!,sprand_sdd!, fdrand,fdrand!,fdrand_coo, solverbenchmark
 
@@ -51,6 +49,8 @@ function __init__()
     @require AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c" include("amg.jl")
 end
 
+export ILUTPreconditioner, AMGPreconditioner
+export PardisoLU, MKLPardisoLU,SparspakLU
 
 
 end # module

src/extendable.jl

@@ -44,6 +44,7 @@ function ExtendableSparseMatrix(valuetype::Type{Tv},indextype::Type{Ti},m, n) wh
 end
 
 ExtendableSparseMatrix(valuetype::Type{Tv},m, n) where{Tv}=ExtendableSparseMatrix{Tv,Int}(m,n)
+
 ExtendableSparseMatrix(m, n)=ExtendableSparseMatrix{Float64,Int}(m,n)
 
 """
@@ -55,6 +56,59 @@ function ExtendableSparseMatrix(csc::SparseMatrixCSC{Tv,Ti}) where{Tv,Ti<:Intege
     return ExtendableSparseMatrix{Tv,Ti}(csc, nothing, phash(csc))
 end
 
+
+"""
+$(SIGNATURES)
+
+ Create ExtendableSparseMatrix from Diagonal
+"""
+ExtendableSparseMatrix(D::Diagonal)=ExtendableSparseMatrix(sparse(D))
+
+"""
+$(SIGNATURES)
+
+ Create ExtendableSparseMatrix from AbstractMatrix, dropping all zero entries.
+ This is the equivalent to `sparse(A)`.
+"""
+ExtendableSparseMatrix(A::AbstractMatrix)=ExtendableSparseMatrix(sparse(A))
+
+"""
+    ExtendableSparseMatrix(I,J,V)
+    ExtendableSparseMatrix(I,J,V,m,n)
+    ExtendableSparseMatrix(I,J,V,combine::Function)
+    ExtendableSparseMatrix(I,J,V,m,n,combine::Function)
+
+Create ExtendableSparseMatrix from triplet (COO) data.
+"""
+ExtendableSparseMatrix(I,J,V::AbstractVector)=ExtendableSparseMatrix(sparse(I,J,V))
+ExtendableSparseMatrix(I,J,V::AbstractVector,m,n)=ExtendableSparseMatrix(sparse(I,J,V,m,n))
+ExtendableSparseMatrix(I,J,V::AbstractVector,combine::Function)=ExtendableSparseMatrix(sparse(I,J,V,combine))
+ExtendableSparseMatrix(I,J,V::AbstractVector,m,n,combine::Function)=ExtendableSparseMatrix(sparse(I,J,V,m,n,combine))
+
+
+# THese are probably too much...
+# function Base.transpose(A::ExtendableSparseMatrix)
+#     flush!(A)
+#     ExtendableSparseMatrix(Base.transpose(sparse(A)))
+# end
+# function Base.adjoint(A::ExtendableSparseMatrix)
+#     flush!(A)
+#     ExtendableSparseMatrix(Base.adjoint(sparse(A)))
+# end
+# function SparseArrays.sparse(text::LinearAlgebra.Transpose{Tv,ExtendableSparseMatrix{Tv,Ti}}) where {Tv,Ti}
+#     transpose(sparse(parent(text)))
+# end
+
+"""
+$(SIGNATURES)
+
+ Create SparseMatrixCSC from ExtendableSparseMatrix
+"""
+function SparseArrays.sparse(A::ExtendableSparseMatrix)
+    flush!(A)
+    A.cscmatrix
+end
+
 """
 $(SIGNATURES)
 
@@ -211,7 +265,7 @@ $(SIGNATURES)
 [`flush!`](@ref) and return number of nonzeros in ext.cscmatrix.
 """
 function SparseArrays.nnz(ext::ExtendableSparseMatrix)
-    @inbounds flush!(ext)
+     flush!(ext)
     return nnz(ext.cscmatrix)
 end
 
@@ -222,7 +276,7 @@ $(SIGNATURES)
 [`flush!`](@ref) and return nonzeros in ext.cscmatrix.
 """
 function SparseArrays.nonzeros(ext::ExtendableSparseMatrix)
-    @inbounds flush!(ext)
+     flush!(ext)
     return nonzeros(ext.cscmatrix)
 end
 
@@ -234,7 +288,7 @@ $(SIGNATURES)
 [`flush!`](@ref) and return rowvals in ext.cscmatrix.
 """
 function SparseArrays.rowvals(ext::ExtendableSparseMatrix)
-    @inbounds flush!(ext)
+     flush!(ext)
     rowvals(ext.cscmatrix)
 end
 
@@ -244,13 +298,16 @@ function SparseArrays.getrowval(S::ExtendableSparseMatrix)
 end
 
 
+
+
+
 """
 $(SIGNATURES)
 
 [`flush!`](@ref) and return colptr of  in ext.cscmatrix.
 """
 function SparseArrays.getcolptr(ext::ExtendableSparseMatrix)
-    @inbounds flush!(ext)
+     flush!(ext)
     return ext.cscmatrix.colptr
 end
 
@@ -261,7 +318,7 @@ $(SIGNATURES)
 [`flush!`](@ref) and return findnz(ext.cscmatrix).
 """
 function SparseArrays.findnz(ext::ExtendableSparseMatrix)
-    @inbounds flush!(ext)
+     flush!(ext)
     return findnz(ext.cscmatrix)
 end
 
@@ -334,7 +391,7 @@ $(SIGNATURES)
 [`flush!`](@ref) and ldiv with ext.cscmatrix
 """
 function  LinearAlgebra.ldiv!(r, ext::ExtendableSparse.ExtendableSparseMatrix, x)
-    @inbounds flush!(ext)
+     flush!(ext)
     return LinearAlgebra.ldiv!(r,ext.cscmatrix,x)
 end
 
@@ -344,7 +401,7 @@ $(SIGNATURES)
 [`flush!`](@ref) and multiply with ext.cscmatrix
 """
 function  LinearAlgebra.mul!(r,ext::ExtendableSparse.ExtendableSparseMatrix, x)
-    @inbounds flush!(ext)
+     flush!(ext)
     return LinearAlgebra.mul!(r,ext.cscmatrix,x)
 end
 
@@ -355,7 +412,7 @@ $(SIGNATURES)
 [`flush!`](@ref) and calculate norm from cscmatrix
 """
 function LinearAlgebra.norm(A::ExtendableSparseMatrix, p::Real=2)
-    @inbounds flush!(A)
+     flush!(A)
     return LinearAlgebra.norm(A.cscmatrix,p)
 end
 
@@ -365,7 +422,7 @@ $(SIGNATURES)
 [`flush!`](@ref) and calculate opnorm from cscmatrix
 """
 function LinearAlgebra.opnorm(A::ExtendableSparseMatrix, p::Real=2)
-    @inbounds flush!(A)
+     flush!(A)
     return LinearAlgebra.opnorm(A.cscmatrix,p)
 end
 
@@ -375,7 +432,7 @@ $(SIGNATURES)
 [`flush!`](@ref) and calculate cond from cscmatrix
 """
 function LinearAlgebra.cond(A::ExtendableSparseMatrix, p::Real=2)
-    @inbounds flush!(A)
+     flush!(A)
     return LinearAlgebra.cond(A.cscmatrix,p)
 end
 
@@ -385,7 +442,7 @@ $(SIGNATURES)
 [`flush!`](@ref) and check for symmetry of cscmatrix
 """
 function LinearAlgebra.issymmetric(A::ExtendableSparseMatrix)
-    @inbounds flush!(A)
+     flush!(A)
     return LinearAlgebra.issymmetric(A.cscmatrix)
 end
 
@@ -397,17 +454,53 @@ $(SIGNATURES)
 Add SparseMatrixCSC matrix and [`ExtendableSparseMatrix`](@ref)  ext.
 """
 function Base.:+(ext::ExtendableSparseMatrix,csc::SparseMatrixCSC)
-    @inbounds flush!(ext)
+     flush!(ext)
     return ext.cscmatrix+csc
 end
 
+function Base.:+(A::ExtendableSparseMatrix,B::ExtendableSparseMatrix)
+    flush!(A)
+    flush!(B)
+    return ExtendableSparseMatrix(A.cscmatrix+B.cscmatrix)
+end
+
+function Base.:-(A::ExtendableSparseMatrix,B::ExtendableSparseMatrix)
+    flush!(A)
+    flush!(B)
+    return ExtendableSparseMatrix(A.cscmatrix-B.cscmatrix)
+end
+
+
+function Base.:*(A::ExtendableSparseMatrix,B::ExtendableSparseMatrix)
+     flush!(A)
+     flush!(B)
+    return ExtendableSparseMatrix(A.cscmatrix*B.cscmatrix)
+end
+
+"""
+$(SIGNATURES)
+"""
+function Base.:*(d::Diagonal,ext::ExtendableSparseMatrix)
+    flush!(ext)
+    return ExtendableSparseMatrix(d*ext.cscmatrix)
+end
+
+"""
+$(SIGNATURES)
+"""
+function Base.:*(ext::ExtendableSparseMatrix,d::Diagonal)
+    flush!(ext)
+    return ExtendableSparseMatrix(ext.cscmatrix*d)
+end
+
+
 """
 $(SIGNATURES)
 
 Subtract  SparseMatrixCSC matrix from  [`ExtendableSparseMatrix`](@ref)  ext.
 """
 function Base.:-(ext::ExtendableSparseMatrix,csc::SparseMatrixCSC)
-    @inbounds flush!(ext)
+     flush!(ext)
     return ext.cscmatrix-csc
 end
 
@@ -417,7 +510,7 @@ $(SIGNATURES)
 Subtract  [`ExtendableSparseMatrix`](@ref)  ext from  SparseMatrixCSC.
 """
 function Base.:-(csc::SparseMatrixCSC,ext::ExtendableSparseMatrix)
-    @inbounds flush!(ext)
+     flush!(ext)
     return csc - ext.cscmatrix
 end
 
@@ -426,17 +519,10 @@ end
 $(SIGNATURES)
 """
 function SparseArrays.dropzeros!(ext::ExtendableSparseMatrix)
-    @inbounds flush!(ext)
+     flush!(ext)
     dropzeros!(ext.cscmatrix)
 end
 
-"""
-$(SIGNATURES)
-"""
-function Base.:*(d::Diagonal,ext::ExtendableSparseMatrix)
-    @inbounds flush!(ext)
-    return ExtendableSparseMatrix(d*ext.cscmatrix)
-end
 
 """
 $(SIGNATURES)

src/factorizations.jl

@@ -16,7 +16,7 @@ and  provide a method
   The idea is that, depending if the matrix pattern has changed, 
   different steps are needed to update the preconditioner.
 
-  Moreover, they have the ExtendableSparseMatrix as a field, ensuring 
+  Moreover, they have the ExtendableSparseMatrix `A` as a field and an `update!` method, ensuring 
   consistency after construction.
 """
 abstract type AbstractFactorization{Tv, Ti} end
@@ -146,15 +146,16 @@ end
 
 include("jacobi.jl")
 include("ilu0.jl")
+include("iluzero.jl")
 include("parallel_jacobi.jl")
 include("parallel_ilu0.jl")
 
 include("SparspakCSCInterface.jl")
 include("sparspak.jl")
 
 @eval begin
-
     @makefrommatrix ILU0Preconditioner
+    @makefrommatrix ILUZeroPreconditioner
     @makefrommatrix JacobiPreconditioner
     @makefrommatrix ParallelJacobiPreconditioner
     @makefrommatrix ParallelILU0Preconditioner