Add product preconditioning

CHANGELOG.md

@@ -1,5 +1,17 @@
 # Changelog
 
+## [2.3.0] - 2026-06-16
+- Product and Identity preconditioners and PreconBuilders
+- Block preconditioner and PreconBuilder now allow for vector 
+  of preconditioners, possibly different for each block
+
+## [2.2.0] - 2026-06-11
+- Base.similar method for GenericSparseMatrix
+
+## [2.1.0] - 2026-06-06
+- Base.sum for ExtendableSparse
+- 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 `\`.

Project.toml

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

src/ExtendableSparse.jl

@@ -43,6 +43,7 @@ export eliminate_dirichlet, eliminate_dirichlet!, mark_dirichlet
 
 include("preconbuilders.jl")
 export LinearSolvePreconBuilder, BlockPreconBuilder, JacobiPreconBuilder
+export ProductPreconBuilder, IdentityPreconBuilder
 @public ILUZeroPreconBuilder, ILUTPreconBuilder
 
 

src/preconbuilders.jl

@@ -102,14 +102,18 @@ end
 
 mutable struct BlockPreconditioner
     A::AbstractMatrix
-    factorization
+    factorizations
     partitioning::Union{Nothing, Vector{AbstractVector}}
     facts::Vector
-    function BlockPreconditioner(A; partitioning = nothing, factorization = nothing)
+    function BlockPreconditioner(A; partitioning = nothing, factorization = nothing, factorizations = nothing)
         p = new()
         p.A = A
         p.partitioning = partitioning
-        p.factorization = factorization
+        if !isnothing(factorization)
+            p.factorizations = [factorization for i in 1:length(partitioning)]
+        else
+            p.factorizations = factorizations
+        end
         update!(p)
         return p
     end
@@ -151,10 +155,8 @@ function update!(precon::BlockPreconditioner)
     np = length(precon.partitioning)
     precon.facts = Vector{Any}(undef, np)
     return Threads.@threads for ipart in 1:np
-        factorization = deepcopy(precon.factorization)
         AP = precon.A[precon.partitioning[ipart], precon.partitioning[ipart]]
-        FP = factorization(AP)
-        precon.facts[ipart] = FP
+        precon.facts[ipart] = precon.factorizations[ipart](AP)
     end
 end
 
@@ -164,12 +166,10 @@ function LinearAlgebra.ldiv!(p::BlockPreconditioner, v)
     facts = p.facts
     np = length(partitioning)
 
-    if allow_views(p.factorization)
-        Threads.@threads for ipart in 1:np
+    Threads.@threads for ipart in 1:np
+        if allow_views(p.factorizations[ipart])
             ldiv!(facts[ipart], view(v, partitioning[ipart]))
-        end
-    else
-        Threads.@threads for ipart in 1:np
+        else
             vv = v[partitioning[ipart]]
             ldiv!(facts[ipart], vv)
             view(v, partitioning[ipart]) .= vv
@@ -182,12 +182,10 @@ 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 in 1:np
+    Threads.@threads for ipart in 1:np
+        if allow_views(p.factorizations[ipart])
             ldiv!(view(u, partitioning[ipart]), facts[ipart], view(v, partitioning[ipart]))
-        end
-    else
-        Threads.@threads for ipart in 1:np
+        else
             uu = u[partitioning[ipart]]
             ldiv!(uu, facts[ipart], v[partitioning[ipart]])
             view(u, partitioning[ipart]) .= uu
@@ -209,7 +207,8 @@ from partition of unknowns.
 - `partitioning(A)`  shall return 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].
 
-- `precs(A,p)` shall return a left precondioner for a matrix block.
+- `precs(A,p)` shall return a left precondioner for a matrix block. It may be one function used for each partition
+   or a vector of functions - one for each partition.
 """
 Base.@kwdef mutable struct BlockPreconBuilder
     precs = UMFPACKPreconBuilder()
@@ -218,8 +217,14 @@ end
 
 function (blockprecs::BlockPreconBuilder)(A, p)
     (; precs, partitioning) = blockprecs
-    factorization = A -> precs(A, p)[1]
-    bp = BlockPreconditioner(A; partitioning = partitioning(A), factorization)
+    Apart = partitioning(A)
+    npart = length(Apart)
+    if isa(precs, Vector)
+        factorizations = [A -> precs[i](A, p)[1] for i in 1:npart]
+    else
+        factorizations = [A -> precs(A, p)[1] for i in 1:npart]
+    end
+    bp = BlockPreconditioner(A; partitioning = Apart, factorizations)
     return (bp, LinearAlgebra.I)
 end
 
@@ -294,3 +299,69 @@ LinearAlgebra.ldiv!(fact::JacobiPreconditioner, v) = ldiv!(fact.factorization, v
 
 allow_views(::JacobiPreconditioner) = true
 allow_views(::Type{JacobiPreconditioner}) = true
+
+"""
+    struct ProductPreconditioner
+
+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*}
+```
+"""
+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
+end
+
+function LinearAlgebra.ldiv!(p::ProductPreconditioner, v)
+    u = ldiv!(copy(v), p, v)
+    v .= u
+    return v
+end
+
+
+"""
+    struct ProductPreconBuilder
+
+LinearSolve `precs` compatible preonditioner constructor for [`ProductPreconditioner`](@ref)
+"""
+Base.@kwdef mutable struct ProductPreconBuilder
+    precs1 = JacobiPreconBuilder()
+    precs2 = JacobiPreconBuilder()
+end
+
+function (prodprecs::ProductPreconBuilder)(A, p)
+    M1 = prodprecs.precs1(A, p)[1]
+    M2 = prodprecs.precs2(A, p)[1]
+    return ProductPreconditioner(A, M1, M2), LinearAlgebra.I
+end
+
+"""
+    struct IdentityPreconBuilder
+
+LinearSolve `precs` compatible preonditioner constructor for trivial preonditioner.
+"""
+struct IdentityPreconBuilder
+end
+
+function (::IdentityPreconBuilder)(A, p)
+    return LinearAlgebra.I, LinearAlgebra.I
+end

test/alltests.jl

@@ -108,3 +108,7 @@ end
 @testset "Block" begin
     include("test_block.jl")
 end
+
+@testset "Product" begin
+    include("test_prod.jl")
+end

test/test_prod.jl

@@ -0,0 +1,71 @@
+module test_prod
+using Test
+using ExtendableSparse
+using ExtendableSparse: ILUZeroPreconBuilder, IdentityPreconBuilder, ProductPreconBuilder
+using IterativeSolvers
+using LinearAlgebra
+using Random
+
+function main(; n = 100)
+    seed = 1340
+    Random.seed!(seed)
+    A = fdrand(n, n)
+    sol0 = ones(n^2)
+    b = A * sol0
+
+    # need to reinit random so the shadow residual is the same for all calls
+    Random.seed!(seed)
+    sol1, hist1 = bicgstabl(A, b, log = true)
+    @test sol1 ≈ sol0 rtol = 1.0e-5
+    @show hist1
+
+
+    IdB = IdentityPreconBuilder()
+    IluB = ILUZeroPreconBuilder()
+    ProdB1 = ProductPreconBuilder(IdB, IluB)
+    ProdB2 = ProductPreconBuilder(IluB, IdB)
+    ProdB3 = ProductPreconBuilder(IluB, IluB)
+
+    IdP, _ = IdB(A, 0)
+    Random.seed!(seed)
+    sol2, hist2 = bicgstabl(A, b, Pl = IdP, log = true)
+    @test sol2 ≈ sol0  rtol = 1.0e-5
+    @show hist2
+    @test hist2.iters ≈ hist1.iters atol = 5
+
+    IluP, _ = IluB(A, 0)
+    Random.seed!(seed)
+    sol3, hist3 = bicgstabl(A, b, Pl = IluP, log = true)
+    @test sol3 ≈ sol0  rtol = 1.0e-5
+    @show hist3
+    @test hist3.iters < hist2.iters
+
+    ProdP1, _ = ProdB1(A, 0)
+    Random.seed!(seed)
+    sol4, hist4 = bicgstabl(A, b, Pl = ProdP1, log = true)
+    @show hist4
+    @test sol4 ≈ sol0  rtol = 1.0e-5
+    @test hist4.iters < hist3.iters
+
+    ProdP2, _ = ProdB2(A, 0)
+    Random.seed!(seed)
+    sol5, hist5 = bicgstabl(A, b, Pl = ProdP2, log = true)
+    @show hist5
+    @test sol5 ≈ sol0  rtol = 1.0e-5
+    @test hist5.iters ≈ hist4.iters atol = 7
+
+    ProdP3, _ = ProdB3(A, 0)
+    Random.seed!(seed)
+    sol6, hist6 = bicgstabl(A, b, Pl = ProdP3, log = true)
+    @show hist6
+    @test sol6 ≈ sol0  rtol = 1.0e-5
+    @test hist6.iters < hist5.iters
+
+
+    return
+end
+
+main()
+
+
+end