Update block preconditioning

.codespellrc

@@ -1,2 +1,2 @@
 [codespell]
-ignore-words-list = missings,rcall,linke,fo,coo,alledges
+ignore-words-list = missings,rcall,linke,fo,coo,alledges,bu

CHANGELOG.md

@@ -1,5 +1,10 @@
 # Changelog
 
+## [2.4.0] - 2026-07-12
+- JacobiPreconBuilder now allows for blocksize
+- More tests
+- Internal reorganization (simplified BlockPreconditioner etc.)
+
 ## [2.3.0] - 2026-06-16
 - Product and Identity preconditioners and PreconBuilders
 - Block preconditioner and PreconBuilder now allow for vector 
@@ -13,7 +18,6 @@
 - fix allocations in lnk+csc matrix addition
 
 ## [2.0.0] - 2026-01-06
-
 - Remove solver + precon API which is not based on precs or directly overloading `\`.
   Fully rely on LinearSolve (besides `\`)
 - Move AMGBuilder etc to corresponding packages (depending on the PRs)

Project.toml

@@ -1,6 +1,6 @@
 name = "ExtendableSparse"
 uuid = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
-version = "2.3.1"
+version = "2.4.0"
 authors = ["Juergen Fuhrmann <juergen.fuhrmann@wias-berlin.de>", "Daniel Runge"]
 
 [deps]

src/ExtendableSparse.jl

@@ -41,6 +41,7 @@ export flush!, updateindex!, rawupdateindex!, reset!, nnznew
 export eliminate_dirichlet, eliminate_dirichlet!, mark_dirichlet
 
 
+include("preconditioners.jl")
 include("preconbuilders.jl")
 export LinearSolvePreconBuilder, BlockPreconBuilder, JacobiPreconBuilder
 export ProductPreconBuilder, IdentityPreconBuilder

src/preconbuilders.jl

@@ -1,3 +1,9 @@
+#
+# This file defines `preconbuilders` which allow to specify `precs` preconditioner constructors
+# for iterative solvers as they are handled in LinearSolve.jl.
+# See https://docs.sciml.ai/LinearSolve/stable/basics/Preconditioners/#Specifying-Preconditioners
+#
+
 """
         LinearSolvePreconBuilder(; method=UMFPACKFactorization())
 
@@ -11,17 +17,7 @@ end
 
 
 """
-    JacobiPreconBuilder()
-
-Return callable object constructing a left Jacobi preconditioner
-to be passed as the `precs` parameter to iterative methods wrapped by LinearSolve.jl.
-"""
-struct JacobiPreconBuilder end
-(::JacobiPreconBuilder)(A::AbstractSparseMatrixCSC, p) = (JacobiPreconditioner(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A))), LinearAlgebra.I)
-
-
-"""
-    ILUZeroPreconBuilder(;blocksize=1)
+    ILUZeroPreconBuilder(; blocksize = 1)
 
 Return callable object constructing a left zero fill-in ILU preconditioner 
 using [ILUZero.jl](https://github.com/mcovalt/ILUZero.jl)
@@ -30,66 +26,17 @@ Base.@kwdef mutable struct ILUZeroPreconBuilder
     blocksize::Int = 1
 end
 
-struct ILUBlockPrecon{N, NN, Tv, Ti}
-    ilu0::ILUZero.ILU0Precon{SMatrix{N, N, Tv, NN}, Ti, SVector{N, Tv}}
-end
-
-function LinearAlgebra.ldiv!(
-        Y::Vector{Tv},
-        A::ILUBlockPrecon{N, NN, Tv, Ti},
-        B::Vector{Tv}
-    ) where {N, NN, Tv, Ti}
-    BY = reinterpret(SVector{N, Tv}, Y)
-    BB = reinterpret(SVector{N, Tv}, B)
-    ldiv!(BY, A.ilu0, BB)
-    return Y
-end
-
-"""
-    pointblock(matrix,blocksize)
-
-Create a pointblock matrix.
-"""
-function pointblock(A0::ExtendableSparseMatrixCSC{Tv, Ti}, blocksize) where {Tv, Ti}
-    A = SparseMatrixCSC(A0)
-    colptr = A.colptr
-    rowval = A.rowval
-    nzval = A.nzval
-    n = A.n
-    block = zeros(Tv, blocksize, blocksize)
-    nblock = n ÷ blocksize
-    b = SMatrix{blocksize, blocksize}(block)
-    Tb = typeof(b)
-    Ab = ExtendableSparseMatrixCSC{Tb, Ti}(nblock, nblock)
-
-
-    for i in 1:n
-        for k in colptr[i]:(colptr[i + 1] - 1)
-            j = rowval[k]
-            iblock = (i - 1) ÷ blocksize + 1
-            jblock = (j - 1) ÷ blocksize + 1
-            ii = (i - 1) % blocksize + 1
-            jj = (j - 1) % blocksize + 1
-            block[ii, jj] = nzval[k]
-            rawupdateindex!(Ab, +, SMatrix{blocksize, blocksize}(block), iblock, jblock)
-            block[ii, jj] = zero(Tv)
-        end
-    end
-    return flush!(Ab)
-end
-
 function (b::ILUZeroPreconBuilder)(A0, p)
     A = SparseMatrixCSC(size(A0)..., getcolptr(A0), rowvals(A0), nonzeros(A0))
     return if b.blocksize == 1
         (ILUZero.ilu0(A), LinearAlgebra.I)
     else
-        (ILUBlockPrecon(ILUZero.ilu0(pointblock(A, b.blocksize), SVector{b.blocksize, eltype(A)})), LinearAlgebra.I)
+        (ILU0BlockPrecon(ILUZero.ilu0(pointblock(A, b.blocksize), SVector{b.blocksize, eltype(A)})), LinearAlgebra.I)
     end
 end
 
-
 """
-    ILUTPreconBuilder(; droptol=0.1)
+    ILUTPreconBuilder(; droptol = 0.1)
 
 Return callable object constructing a left ILUT preconditioner 
 using [IncompleteLU.jl](https://github.com/haampie/IncompleteLU.jl)
@@ -100,103 +47,6 @@ end
 (::ILUTPreconBuilder)(A, p) = error("import IncompleteLU.jl in order to use ILUTBuilder")
 
 
-mutable struct BlockPreconditioner
-    A::AbstractMatrix
-    factorizations
-    partitioning::Union{Nothing, Vector{AbstractVector}}
-    facts::Vector
-    function BlockPreconditioner(A; partitioning = nothing, factorization = nothing, factorizations = nothing)
-        p = new()
-        p.A = A
-        p.partitioning = partitioning
-        if !isnothing(factorization)
-            p.factorizations = [factorization for i in 1:length(partitioning)]
-        else
-            p.factorizations = factorizations
-        end
-        update!(p)
-        return p
-    end
-end
-
-
-"""
-     BlockPreconditioner(;partitioning, factorization)
-
-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)
-    return Threads.@threads for ipart in 1:np
-        AP = precon.A[precon.partitioning[ipart], precon.partitioning[ipart]]
-        precon.facts[ipart] = precon.factorizations[ipart](AP)
-    end
-end
-
-
-function LinearAlgebra.ldiv!(p::BlockPreconditioner, v)
-    partitioning = p.partitioning
-    facts = p.facts
-    np = length(partitioning)
-
-    Threads.@threads for ipart in 1:np
-        if allow_views(p.factorizations[ipart])
-            ldiv!(facts[ipart], view(v, partitioning[ipart]))
-        else
-            vv = v[partitioning[ipart]]
-            ldiv!(facts[ipart], vv)
-            view(v, partitioning[ipart]) .= vv
-        end
-    end
-    return v
-end
-
-function LinearAlgebra.ldiv!(u, p::BlockPreconditioner, v)
-    partitioning = p.partitioning
-    facts = p.facts
-    np = length(partitioning)
-    Threads.@threads for ipart in 1:np
-        if allow_views(p.factorizations[ipart])
-            ldiv!(view(u, partitioning[ipart]), facts[ipart], view(v, partitioning[ipart]))
-        else
-            uu = u[partitioning[ipart]]
-            ldiv!(uu, facts[ipart], v[partitioning[ipart]])
-            view(u, partitioning[ipart]) .= uu
-        end
-    end
-    return u
-end
-
-Base.eltype(p::BlockPreconditioner) = eltype(p.facts[1])
-
-
 """
      BlockPreconBuilder(;precs=UMFPACKPreconBuilder(),  
                          partitioning = A -> [1:size(A,1)]
@@ -228,120 +78,25 @@ function (blockprecs::BlockPreconBuilder)(A, p)
     return (bp, LinearAlgebra.I)
 end
 
-"""
-    Allow array for precs => different precoms
-"""
-
-
-mutable struct _JacobiPreconditioner{Tv}
-    invdiag::Vector{Tv}
-end
-
-function jacobi(A::SparseMatrixCSC{Tv, Ti}) where {Tv, Ti}
-    invdiag = Array{Tv, 1}(undef, A.n)
-    n = A.n
-    @inbounds for i in 1:n
-        invdiag[i] = one(Tv) / A[i, i]
-    end
-    return _JacobiPreconditioner(invdiag)
-end
-
-function jacobi!(p::_JacobiPreconditioner{Tv}, A::SparseMatrixCSC{Tv, Ti}) where {Tv, Ti}
-    n = A.n
-    @inbounds for i in 1:n
-        p.invdiag[i] = one(Tv) / A[i, i]
-    end
-    return p
-end
-
-mutable struct JacobiPreconditioner
-    A::AbstractMatrix
-    factorization::Union{_JacobiPreconditioner, Nothing}
-    function JacobiPreconditioner(A)
-        p = new()
-        p.A = A
-        p.factorization = nothing
-        update!(p)
-        return p
-    end
-end
-
-function LinearAlgebra.ldiv!(u, p::_JacobiPreconditioner, v)
-    n = length(p.invdiag)
-    for i in 1:n
-        @inbounds u[i] = p.invdiag[i] * v[i]
-    end
-    return u
-end
-
-LinearAlgebra.ldiv!(p::_JacobiPreconditioner, v) = ldiv!(v, p, v)
-
-
-"""
-```
-JacobiPreconditioner()
-JacobiPreconditioner(matrix)
-```
-
-Jacobi preconditioner.
-"""
-function JacobiPreconditioner end
-
-function update!(p::JacobiPreconditioner)
-    flush!(p.A)
-    Tv = eltype(p.A)
-    p.factorization = jacobi(SparseMatrixCSC(p.A))
-    return p
-end
-
-LinearAlgebra.ldiv!(u, fact::JacobiPreconditioner, v) = ldiv!(u, fact.factorization, v)
-LinearAlgebra.ldiv!(fact::JacobiPreconditioner, v) = ldiv!(fact.factorization, v)
-
-allow_views(::JacobiPreconditioner) = true
-allow_views(::Type{JacobiPreconditioner}) = true
 
 """
-    struct ProductPreconditioner
+    JacobiPreconBuilder(; blocksize = 1)
 
-Product of two left preconditioning steps.
-The operation ``u=M^{-1}v`` is defined by two simple iteration steps
-with two different preconditioners ``M_1`` and ``M_2``:
-Let ``u_0=0``. Then calculate
-```math    
- \\begin{align*}
-   u_1&= u_0 - M_1^{-1}(Au_0 - v) = M_1^{-1}v\\\\
-   u  &= u_1 - M_2^{-1}(Au_1 - v)
-  \\end{align*}
-```
+Return callable object constructing a left Jacobi preconditioner
+to be passed as the `precs` parameter to iterative methods wrapped by LinearSolve.jl.
 """
-Base.@kwdef struct ProductPreconditioner{TA, TM1, TM2}
-    A::TA
-    M1::TM1
-    M2::TM2
-end
-
-function LinearAlgebra.ldiv!(u, p::ProductPreconditioner, v)
-    (; A, M1, M2) = p
-    u1 = similar(u)
-    u2 = similar(u)
-    ldiv!(u1, M1, v)
-    mul!(u2, A, u1)
-    ldiv!(u, M2, v - u2)
-    u .+= u1
-    return u
+Base.@kwdef struct JacobiPreconBuilder
+    blocksize::Int = 1
 end
 
-function LinearAlgebra.ldiv!(p::ProductPreconditioner, v)
-    u = ldiv!(copy(v), p, v)
-    v .= u
-    return v
-end
+(b::JacobiPreconBuilder)(A::AbstractSparseMatrixCSC, p) = (JacobiPreconditioner(A; blocksize = b.blocksize), LinearAlgebra.I)
 
 
 """
-    struct ProductPreconBuilder
+    ProductPreconBuilder(precs1, precs2)
 
-LinearSolve `precs` compatible preonditioner constructor for [`ProductPreconditioner`](@ref)
+Return LinearSolve `precs` compatible callable object which constructs
+a  [`ProductPreconditioner`](@ref) from the  `precs1` and `precs2` preconbuilders.
 """
 Base.@kwdef mutable struct ProductPreconBuilder
     precs1 = JacobiPreconBuilder()