Infrastructure for factorizations including Pardiso

Infrastructure for factorizations including Pardiso
* Reusable LU factorization
* infrastructure for Pardiso, AlgebraicMultigrid, IncompleteLU etc. based on Require.jl
* minmal julia version -> 1.5 (because of lu!)
* Unittests for packages loaded via require
* Documentation update
* Define generic default methods for factorizations, specialize
only if necessary.

.github/workflows/ci.yml

@@ -15,7 +15,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.3' # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'.
+          - '1.5' # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'.
           - '1'   # Leave this line unchanged. '1' will automatically expand to the latest stable 1.x release of Julia.
 #          - 'nightly'
         os:

Project.toml

@@ -1,18 +1,21 @@
 name = "ExtendableSparse"
 uuid = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
 authors = ["Juergen Fuhrmann <[email protected]>"]
-version = "0.4.2"
+version = "0.5.0"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 DocStringExtensions = "^0.8.0"
-julia = "1"
+Requires = "^1.1.3"
+julia = "^1.5"
 
 [extras]
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"

docs/Project.toml

@@ -1,7 +1,9 @@
 [deps]
+AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8"
 ExtendableSparse = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
-
+IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895"
+Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"

docs/make.jl

@@ -1,6 +1,5 @@
 push!(LOAD_PATH,"../src/")
-using Documenter, ExtendableSparse
-
+using Documenter, ExtendableSparse,Pardiso,AlgebraicMultigrid,IncompleteLU
 
 makedocs(sitename="ExtendableSparse.jl",
          modules = [ExtendableSparse],
@@ -13,6 +12,7 @@ makedocs(sitename="ExtendableSparse.jl",
              "example.md",
              "extsparse.md",
              "iter.md",
+             "internal.md",
              "changes.md",
          ])
 

docs/src/changes.md

@@ -1,5 +1,12 @@
 # Changes
-
+## v0.5, April 10, 2021
+- Introduce lu/lu! , factorize/factorize!, unifying LU factorizations and preconditioners
+- Interface packages: Pardiso, AlgebraicMultigrid, IncompleteLU via Requires.jl
+## v0.4, March 2021
+- Fix handling of Symmetrix matrices
+- `rawupdateindex` does not check for entering zeros
+- Compare with `COO` method
+- Benchmarks in documentation
 ## v0.3.7, March 20, 2021
 - Added parallel jacobi preconditioner (thanks, @jkr)
 - Fixes ldiv

docs/src/extsparse.md

@@ -1,7 +1,14 @@
 # Sparse matrix handling
 
+## Matrix creation and update API
 ```@autodocs
 Modules = [ExtendableSparse]
-Pages = ["sparsematrixlnk.jl","sparsematrixcsc.jl","extendable.jl", "sprand.jl"]
+Pages = ["extendable.jl"]
+```
+
+## Test matrix creation
+```@autodocs
+Modules = [ExtendableSparse]
+Pages = ["sprand.jl"]
 ```
 

docs/src/internal.md

@@ -0,0 +1,15 @@
+# Internal API
+
+
+## Linked List Sparse Matrix format
+```@autodocs
+Modules = [ExtendableSparse]
+Pages = ["sparsematrixlnk.jl"]
+```
+
+## Some methods for SparseMatrixCSC
+```@autodocs
+Modules = [ExtendableSparse]
+Pages = ["sparsematrixcsc.jl"]
+```
+

docs/src/iter.md

@@ -1,5 +1,27 @@
-# Preconditioners
+# Factorizations & Preconditioners 
+
+## Factorizations
+
+In this package, preconditioners and LU factorizations are subcategories are both seen
+as complete or approximate _factorizations_. Correspondingly there is a common API for
+their creation.
+
+Factorizations from these package know the matrices which have been factorized.
+
 ```@autodocs
 Modules = [ExtendableSparse]
-Pages = ["preconditioners.jl", "ilu0.jl", "jacobi.jl", "parallel_jacobi.jl","simple_iteration.jl"]
+Pages = ["factorizations.jl"]
 ```
+
+## LU Factorizations
+```@autodocs
+Modules = [ExtendableSparse]
+Pages = ["umfpack_lu.jl", "pardiso_lu.jl"]
+```
+
+## Preconditioners
+```@autodocs
+Modules = [ExtendableSparse]
+Pages = ["jacobi.jl","ilu0.jl","parallel_jacobi.jl","ilut.jl","amg.jl"]
+```
+

src/ExtendableSparse.jl

@@ -1,19 +1,41 @@
 module ExtendableSparse
-using DocStringExtensions
 using SparseArrays
 using LinearAlgebra
+using SuiteSparse
+using Requires
+
+using DocStringExtensions
+
+
 
 include("sparsematrixcsc.jl")
 include("sparsematrixlnk.jl")
 include("extendable.jl")
+
 export SparseMatrixLNK,ExtendableSparseMatrix,flush!,nnz, updateindex!, rawupdateindex!, colptrs
 
-include("preconditioners.jl")
+include("factorizations.jl")
 export JacobiPreconditioner, ILU0Preconditioner, ParallelJacobiPreconditioner
+export issolver
+export factorize,factorize!, update!
+export ILUTPreconditioner, AMGPreconditioner
 
 include("simple_iteration.jl")
 export simple,simple!
 
+
+
+
 include("sprand.jl")
 export sprand!,sprand_sdd!, fdrand,fdrand!,fdrand_coo
+
+
+function __init__()
+    @require Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2" include("pardiso_lu.jl")
+    @require IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895" include("ilut.jl")
+    @require AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c" include("amg.jl")
+end
+
+
+
 end # module

src/amg.jl

@@ -0,0 +1,26 @@
+"""
+$(TYPEDEF)
+
+"""
+mutable struct AMGPreconditioner{Tv, Ti} <: AbstractExtendableSparsePreconditioner{Tv,Ti}
+    A::ExtendableSparseMatrix{Tv,Ti}
+    fact
+end
+
+"""
+```
+AMGPreconditioner(extmatrix)
+```
+"""
+function AMGPreconditioner(A::ExtendableSparseMatrix{Tv,Ti}) where {Tv,Ti}
+    @inbounds flush!(A)
+    p=AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(A.cscmatrix))
+    AMGPreconditioner{Tv,Ti}(A,p)
+end
+
+function update!(precon::AMGPreconditioner{Tv,Ti}) where {Tv,Ti}
+    @inbounds flush!(precon.A)
+    precon.fact=AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(precon.A.cscmatrix))
+end
+
+

src/extendable.jl

@@ -19,60 +19,46 @@ mutable struct ExtendableSparseMatrix{Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv
     lnkmatrix::Union{SparseMatrixLNK{Tv,Ti},Nothing}
 
     """
-    Time stamp of last pattern update
+    Pattern hash
     """
-    pattern_timestamp::Float64
+    phash::UInt64
 end
 
 
 """
-$(SIGNATURES)
+```
+ExtendableSparseMatrix(Tv,Ti,m,n)
+ExtendableSparseMatrix(Tv,m,n)
+ExtendableSparseMatrix(m,n)
+```
+Create empty ExtendableSparseMatrix. This is equivalent to `spzeros(m,n)` for
+`SparseMartrixCSC`.
 
-Create empty ExtendableSparseMatrix.
 """
 function ExtendableSparseMatrix{Tv,Ti}(m, n) where{Tv,Ti<:Integer}
-    ExtendableSparseMatrix{Tv,Ti}(spzeros(Tv,Ti,m,n),nothing,time())
+    ExtendableSparseMatrix{Tv,Ti}(spzeros(Tv,Ti,m,n),nothing,0)
 end
 
-"""
-$(SIGNATURES)
-
-Create empty ExtendableSparseMatrix.
-"""
 function ExtendableSparseMatrix(valuetype::Type{Tv},indextype::Type{Ti},m, n) where{Tv,Ti<:Integer}
     ExtendableSparseMatrix{Tv,Ti}(m,n)
 end
 
-"""
-$(SIGNATURES)
-
-Create empty ExtendablSparseMatrix.
-This is a pendant to spzeros.
-"""
 ExtendableSparseMatrix(valuetype::Type{Tv},m, n) where{Tv}=ExtendableSparseMatrix{Tv,Int}(m,n)
-
-
-"""
-$(SIGNATURES)
-
-Create empty ExtendableSparseMatrix.
-This is a pendant to spzeros.
-"""
 ExtendableSparseMatrix(m, n)=ExtendableSparseMatrix{Float64,Int}(m,n)
 
-
 """
 $(SIGNATURES)
 
-  Create ExtendableSparseMatrix from SparseMatrixCSC
+ Create ExtendableSparseMatrix from SparseMatrixCSC
 """
 function ExtendableSparseMatrix(csc::SparseMatrixCSC{Tv,Ti}) where{Tv,Ti<:Integer}
-    return ExtendableSparseMatrix{Tv,Ti}(csc, nothing, time())
+    return ExtendableSparseMatrix{Tv,Ti}(csc, nothing, phash(csc))
 end
 
 """
 $(SIGNATURES)
-  Create similar extendableSparseMatrix
+
+Create similar but emtpy extendableSparseMatrix
 """
 Base.similar(m::ExtendableSparseMatrix{Tv,Ti}) where {Tv,Ti}=ExtendableSparseMatrix{Tv,Ti}(size(m)...)
 
@@ -119,7 +105,6 @@ end
 $(SIGNATURES)
 Like [`updateindex!`](@ref) but without 
 checking if v is zero.
-
 """
 function rawupdateindex!(ext::ExtendableSparseMatrix{Tv,Ti}, op,v, i,j) where {Tv,Ti<:Integer}
     k=findindex(ext.cscmatrix,i,j)
@@ -141,7 +126,6 @@ $(SIGNATURES)
 Find index in CSC matrix and set value if it exists. Otherwise,
 set index in extension if `v` is nonzero.
 """
-
 function Base.setindex!(ext::ExtendableSparseMatrix{Tv,Ti}, v, i,j) where {Tv,Ti}
     k=findindex(ext.cscmatrix,i,j)
     if k>0
@@ -199,8 +183,6 @@ function Base.show(io::IO,::MIME"text/plain",ext::ExtendableSparseMatrix)
     end
 end
 
-
-
 """
 $(SIGNATURES)
 
@@ -211,7 +193,7 @@ function flush!(ext::ExtendableSparseMatrix)
     if ext.lnkmatrix!=nothing && nnz(ext.lnkmatrix)>0
         ext.cscmatrix=ext.lnkmatrix+ext.cscmatrix
         ext.lnkmatrix=nothing
-        ext.pattern_timestamp=time()
+        ext.phash=phash(ext.cscmatrix)
     end
     return ext
 end
@@ -273,17 +255,6 @@ function SparseArrays.findnz(ext::ExtendableSparseMatrix)
 end
 
 
-"""
-$(SIGNATURES)
-
-[`flush!`](@ref) and return LU factorization of ext.cscmatrix
-"""
-function LinearAlgebra.lu(ext::ExtendableSparseMatrix)
-    @inbounds flush!(ext)
-    return LinearAlgebra.lu(ext.cscmatrix)
-end
-
-
 
 """
 $(SIGNATURES)
@@ -347,7 +318,7 @@ $(SIGNATURES)
 [`flush!`](@ref) and calculate norm from cscmatrix
 """
 function LinearAlgebra.norm(A::ExtendableSparseMatrix, p::Real=2)
-    @time @inbounds flush!(A)
+    @inbounds flush!(A)
     return LinearAlgebra.norm(A.cscmatrix,p)
 end