Block preconditioning

* Generic block preconditioner

* Removed element types from factorization wrapper

* Code reorganization

Author: j-fu

Project.toml

@@ -1,7 +1,7 @@
 name = "ExtendableSparse"
 uuid = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
 authors = ["Juergen Fuhrmann <[email protected]>"]
-version = "1.0.1"
+version = "1.1.0"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"

docs/src/changes.md

@@ -1,4 +1,10 @@
 # Changes
+## v1.1, May 3, 2023
+  - AbstractFactorization and subtypes are now without element and index type information.
+    They wrap more concretely typed info. This shall allow to construct a preconditioner
+    without knowing matrix and element type.
+  - First steps to block preconditioning
+  - src directory restructured
 
 ## v1.0.1, Jan 22, 2023
   - support of AbstractSparseMatrixCSC interface

docs/src/iter.md

@@ -5,8 +5,8 @@ This functionality probably will be reduced in favor of LinearSolve.jl.
 ## Factorizations
 
 In this package, preconditioners and LU factorizations are both seen
-as complete or approximate _factorizations_. Correspondingly we provide a common  API for
-their creation.
+as complete or approximate _factorizations_. Correspondingly we provide a common  API for them.
+
 
 ```@autodocs
 Modules = [ExtendableSparse]
@@ -122,7 +122,15 @@ nm1, nm2
 
 ### API
 
+#### Recommended 
+```@autodocs
+Modules = [ExtendableSparse]
+Pages = ["iluzero.jl","ilut.jl","amg.jl","blockpreconditioner.jl"]
+```
+
+
+#### Experimental
 ```@autodocs
 Modules = [ExtendableSparse]
-Pages = ["jacobi.jl","parallel_jacobi.jl","iluzero.jl","ilu0.jl","ilut.jl","amg.jl"]
+Pages = ["jacobi.jl","parallel_jacobi.jl","ilu0.jl",]
 ```

docs/src/linearsolve.md

@@ -38,11 +38,26 @@ Also, the iterative method interface works with ExtendableSparse.
 ```@example
 using ExtendableSparse # hide
 using LinearSolve # hide
+using SparseArrays # hide
 using ILUZero # hide
 A = fdrand(10, 10, 10; matrixtype = ExtendableSparseMatrix)
 x = ones(1000)
 b = A * x
 y = LinearSolve.solve(LinearProblem(A, b), KrylovJL_CG();
-                      Pl = ILUZero.ilu0(sparse(A))).u
+                      Pl = ILUZero.ilu0(SparseMatrixCSC(A))).u
+sum(y)
+```
+
+However, ExtendableSparse provides a number of wrappers around preconditioners
+from various Julia packages.
+```@example
+using ExtendableSparse # hide
+using LinearSolve # hide
+using ILUZero # hide
+A = fdrand(10, 10, 10; matrixtype = ExtendableSparseMatrix)
+x = ones(1000)
+b = A * x
+y = LinearSolve.solve(LinearProblem(A, b), KrylovJL_CG();
+                      Pl = ILUZeroPreconditioner(A)).u
 sum(y)
 ```

src/ExtendableSparse.jl

@@ -18,38 +18,41 @@ using DocStringExtensions
 
 import SparseArrays: AbstractSparseMatrixCSC, rowvals, getcolptr, nonzeros
 
-include("sparsematrixcsc.jl")
-include("sparsematrixlnk.jl")
-include("extendable.jl")
+include("matrix/sparsematrixcsc.jl")
+include("matrix/sparsematrixlnk.jl")
+include("matrix/extendable.jl")
 
 export SparseMatrixLNK,
        ExtendableSparseMatrix, flush!, nnz, updateindex!, rawupdateindex!, colptrs, sparse
 
-include("factorizations.jl")
+include("factorizations/factorizations.jl")
 
 export JacobiPreconditioner,
-       ILU0Preconditioner,
-       ILUZeroPreconditioner,
-       ParallelJacobiPreconditioner,
-       ParallelILU0Preconditioner,
-       reorderlinsys
+    ILU0Preconditioner,
+    ILUZeroPreconditioner,
+    ParallelJacobiPreconditioner,
+    ParallelILU0Preconditioner,
+    BlockPreconditioner,allow_views,
+    reorderlinsys
+
 export AbstractFactorization, LUFactorization, CholeskyFactorization
 export issolver
 export factorize!, update!
 
-include("simple_iteration.jl")
+include("factorizations/simple_iteration.jl")
 export simple, simple!
 
-include("sprand.jl")
+include("matrix/sprand.jl")
 export sprand!, sprand_sdd!, fdrand, fdrand!, fdrand_coo, solverbenchmark
 
 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")
+    @require Pardiso="46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2" include("factorizations/pardiso_lu.jl")
+    @require IncompleteLU="40713840-3770-5561-ab4c-a76e7d0d7895" include("factorizations/ilut.jl")
+    @require AlgebraicMultigrid="2169fc97-5a83-5252-b627-83903c6c433c" include("factorizations/amg.jl")
 end
 
 export ILUTPreconditioner, AMGPreconditioner
 export PardisoLU, MKLPardisoLU, SparspakLU
 
+
 end # module

src/amg.jl

@@ -0,0 +1,34 @@
+mutable struct AMGPreconditioner <: AbstractPreconditioner
+    A::ExtendableSparseMatrix
+    factorization::AlgebraicMultigrid.Preconditioner
+    max_levels::Int
+    max_coarse::Int
+    function AMGPreconditioner(; max_levels = 10, max_coarse = 10)
+        precon = new()
+        precon.max_levels = max_levels
+        precon.max_coarse = max_coarse
+        precon
+    end
+end
+
+"""
+```
+AMGPreconditioner(;max_levels=10, max_coarse=10)
+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)
+"""
+function AMGPreconditioner end
+
+@eval begin
+    @makefrommatrix AMGPreconditioner
+end
+
+function update!(precon::AMGPreconditioner)
+    @inbounds flush!(precon.A)
+    precon.factorization = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(precon.A.cscmatrix))
+end
+
+allow_views(::AMGPreconditioner)=true
+allow_views(::Type{AMGPreconditioner})=true

src/factorizations/amg.jl

@@ -0,0 +1,101 @@
+mutable struct BlockPreconditioner <: AbstractPreconditioner
+    A::ExtendableSparseMatrix
+    factorization
+    phash::UInt64
+    partitioning::Union{Nothing,Vector{AbstractVector}}
+    facts::Vector
+    function BlockPreconditioner(;partitioning=nothing, factorization=ExtendableSparse.LUFactorization) 
+        p = new()
+        p.phash = 0
+        p.partitioning=partitioning
+        p.factorization=factorization
+        p
+    end
+end
+
+
+
+"""
+     BlockPreconditioner(;partitioning, factorization=LUFactorization)
+    
+Create a block preconditioner from partition of unknowns given by `partitioning`, a vector of AbstractVectors describing the
+indices of the partitions of the matrix. For a matrix of size `n x n`, e.g. partitioning could be `[ 1:n÷2, (n÷2+1):n]`
+or [ 1:2:n, 2:2:n].
+Factorization is a callable (Function or struct) which allows to create a factorization (with `ldiv!` methods) from a submatrix of A.
+"""
+function BlockPreconditioner end
+
+"""
+    allow_views(::preconditioner_type)
+
+Factorizations on matrix partitions within a block preconditioner may or may not work with array views.
+E.g. the umfpack factorization cannot work with views, while ILUZeroPreconditioner can.
+Implementing a method for `allow_views` returning `false` resp. `true` allows to dispatch to the proper case.
+"""
+allow_views(::Any)=false
+
+
+function update!(precon::BlockPreconditioner)
+    flush!(precon.A)
+    nall=sum(length,precon.partitioning)
+    n=size(precon.A,1)
+    if nall!=n
+	@warn "sum(length,partitioning)=$(nall) but n=$(n)"
+    end
+
+    if isnothing(precon.partitioning)
+        partitioning=[1:n]
+    end
+
+    np=length(precon.partitioning)
+    precon.facts=Vector{Any}(undef,np)
+    Threads.@threads for ipart=1:np
+        factorization=deepcopy(precon.factorization)
+        AP=precon.A[precon.partitioning[ipart],precon.partitioning[ipart]]
+        FP=factorization(AP)
+        precon.facts[ipart]=FP
+    end
+end
+
+
+
+
+function LinearAlgebra.ldiv!(p::BlockPreconditioner,v) 
+    partitioning=p.partitioning
+    facts=p.facts
+    np=length(partitioning)
+
+    if allow_views(p.factorization)
+        Threads.@threads for ipart=1:np
+	    ldiv!(facts[ipart],view(v,partitioning[ipart]))
+        end
+    else
+        Threads.@threads for ipart=1:np
+            vv=v[partitioning[ipart]]
+	    ldiv!(facts[ipart],vv)
+            view(v,partitioning[ipart]).=vv
+        end
+    end
+    v
+end
+
+function LinearAlgebra.ldiv!(u,p::BlockPreconditioner,v)
+    partitioning=p.partitioning
+    facts=p.facts
+    np=length(partitioning)
+    
+    if allow_views(p.factorization)
+        Threads.@threads for ipart=1:np
+	    ldiv!(view(u,partitioning[ipart]),facts[ipart],view(v,partitioning[ipart]))
+        end
+    else
+        Threads.@threads for ipart=1:np
+            uu=u[partitioning[ipart]]
+	    ldiv!(uu,facts[ipart],v[partitioning[ipart]])
+            view(u,partitioning[ipart]).=uu
+        end
+    end
+    u
+end
+
+Base.eltype(p::BlockPreconditioner)=eltype(p.facts[1])

src/factorizations/blockpreconditioner.jl

@@ -1,6 +1,6 @@
-mutable struct CholeskyFactorization{Tv, Ti} <: AbstractLUFactorization{Tv, Ti}
-    A::Union{ExtendableSparseMatrix{Tv, Ti}, Nothing}
-    fact::Union{SuiteSparse.CHOLMOD.Factor{Tv}, Nothing}
+mutable struct CholeskyFactorization <: AbstractLUFactorization
+    A::Union{ExtendableSparseMatrix, Nothing}
+    fact::Union{SuiteSparse.CHOLMOD.Factor, Nothing}
     phash::UInt64
     A64::Any
 end
@@ -11,15 +11,15 @@ CholeskyFactorization(matrix)
 
 Default Cholesky factorization via cholmod.
 """
-function CholeskyFactorization(; valuetype::Type = Float64, indextype::Type = Int64)
-    CholeskyFactorization{valuetype, indextype}(nothing, nothing, 0, nothing)
+function CholeskyFactorization()
+    CholeskyFactorization(nothing, nothing, 0, nothing)
 end
 
 function update!(cholfact::CholeskyFactorization)
     A = cholfact.A
     flush!(A)
     if A.phash != cholfact.phash
-        cholfact.A64 = Symmetric(SparseMatrixCSC{Float64, Int64}(A.cscmatrix))
+        cholfact.A64 = Symmetric(A.cscmatrix)
         cholfact.fact = cholesky(cholfact.A64)
         cholfact.phash = A.phash
     else