From b8468acfdca4dca088ee5cd6e86d1cfa836b0089 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Fri, 17 Apr 2026 21:04:20 +0200 Subject: [PATCH 1/3] draft for Product preconditioner --- src/ExtendableSparse.jl | 1 + src/preconbuilders.jl | 95 ++++++++++++++++++++++++++++++++--------- test/test_prod.jl | 59 +++++++++++++++++++++++++ 3 files changed, 136 insertions(+), 19 deletions(-) create mode 100644 test/test_prod.jl diff --git a/src/ExtendableSparse.jl b/src/ExtendableSparse.jl index 00ddbd7..a210cf4 100644 --- a/src/ExtendableSparse.jl +++ b/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 diff --git a/src/preconbuilders.jl b/src/preconbuilders.jl index 93d4994..081e832 100644 --- a/src/preconbuilders.jl +++ b/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,55 @@ LinearAlgebra.ldiv!(fact::JacobiPreconditioner, v) = ldiv!(fact.factorization, v allow_views(::JacobiPreconditioner) = true allow_views(::Type{JacobiPreconditioner}) = true + + +Base.@kwdef struct ProductPreconditioner + A::AbstractMatrix + precon1 + precon2 +end + +Base.@kwdef mutable struct ProductPreconBuilder + precs1 = JacobiPreconBuilder() + precs2 = JacobiPreconBuilder() +end + +function (prodprecs::ProductPreconBuilder)(A, p) + precon1 = prodprecs.precs1(A, p)[1] + precon2 = prodprecs.precs2(A, p)[1] + return ProductPreconditioner(A, precon1, precon2), LinearAlgebra.I +end + +function LinearAlgebra.ldiv!(u, p::ProductPreconditioner, v) + u1 = similar(u) + u2 = similar(u) + u3 = similar(u) + ldiv!(u1, p.precon1, v) + ldiv!(u2, p.precon2, v) + mul!(u3, p.A, -u1) + ldiv!(u, p.precon2, u3) + u .+= u1 + u .+= u2 + return u +end + +function LinearAlgebra.ldiv!(p::ProductPreconditioner, v) + u1 = similar(v) + u2 = similar(v) + u3 = similar(v) + ldiv!(u1, p.precon1, v) + ldiv!(u2, p.precon2, v) + mul!(u3, p.A, -u1) + ldiv!(v, p.precon2, u3) + v .+= u1 + v .+= u2 + return v +end + + +struct IdentityPreconBuilder +end + +function (::IdentityPreconBuilder)(A, p) + return LinearAlgebra.I, LinearAlgebra.I +end diff --git a/test/test_prod.jl b/test/test_prod.jl new file mode 100644 index 0000000..27c5af4 --- /dev/null +++ b/test/test_prod.jl @@ -0,0 +1,59 @@ +module test_prod +using Test +using ExtendableSparse +using ExtendableSparse: ILUZeroPreconBuilder, IdentityPreconBuilder, ProductPreconBuilder +using IterativeSolvers +using LinearAlgebra +using Random + +function main(; n = 100) + Random.seed!(2026) + A = fdrand(n, n) + sol0 = ones(n^2) + b = A * sol0 + + + 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) + sol2, hist2 = bicgstabl(A, b, Pl = IdP, log = true) + @test sol2 ≈ sol0 rtol = 1.0e-5 + @show hist2 + + IluP, _ = IluB(A, 0) + sol3, hist3 = bicgstabl(A, b, Pl = IluP, log = true) + @test sol3 ≈ sol0 + @show hist3 + + ProdP1, _ = ProdB1(A, 0) + sol4, hist4 = bicgstabl(A, b, Pl = ProdP1, log = true) + @show hist4 + @test sol4 ≈ sol0 + + ProdP2, _ = ProdB2(A, 0) + sol5, hist5 = bicgstabl(A, b, Pl = ProdP2, log = true) + @show hist5 + @test sol5 ≈ sol0 + + ProdP3, _ = ProdB3(A, 0) + sol6, hist6 = bicgstabl(A, b, Pl = ProdP3, log = true) + @show hist6 + @test sol6 ≈ sol0 + + + return +end + +main() + + +end From be5af941ba386eee0f18180e50331aeeaa23a952 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Tue, 16 Jun 2026 21:17:52 +0200 Subject: [PATCH 2/3] Fix & document ProductPreconditioner --- Project.toml | 2 +- src/preconbuilders.jl | 76 +++++++++++++++++++++++++------------------ test/alltests.jl | 4 +++ test/test_prod.jl | 24 ++++++++++---- 4 files changed, 68 insertions(+), 38 deletions(-) diff --git a/Project.toml b/Project.toml index 5693471..b7cd762 100644 --- a/Project.toml +++ b/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 ", "Daniel Runge"] [deps] diff --git a/src/preconbuilders.jl b/src/preconbuilders.jl index 081e832..79dbb86 100644 --- a/src/preconbuilders.jl +++ b/src/preconbuilders.jl @@ -300,51 +300,65 @@ LinearAlgebra.ldiv!(fact::JacobiPreconditioner, v) = ldiv!(fact.factorization, v allow_views(::JacobiPreconditioner) = true allow_views(::Type{JacobiPreconditioner}) = true - -Base.@kwdef struct ProductPreconditioner - A::AbstractMatrix - precon1 - precon2 -end - -Base.@kwdef mutable struct ProductPreconBuilder - precs1 = JacobiPreconBuilder() - precs2 = JacobiPreconBuilder() -end - -function (prodprecs::ProductPreconBuilder)(A, p) - precon1 = prodprecs.precs1(A, p)[1] - precon2 = prodprecs.precs2(A, p)[1] - return ProductPreconditioner(A, precon1, precon2), LinearAlgebra.I +""" + 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) - u3 = similar(u) - ldiv!(u1, p.precon1, v) - ldiv!(u2, p.precon2, v) - mul!(u3, p.A, -u1) - ldiv!(u, p.precon2, u3) + ldiv!(u1, M1, v) + mul!(u2, A, u1) + ldiv!(u, M2, v - u2) u .+= u1 - u .+= u2 return u end function LinearAlgebra.ldiv!(p::ProductPreconditioner, v) - u1 = similar(v) - u2 = similar(v) - u3 = similar(v) - ldiv!(u1, p.precon1, v) - ldiv!(u2, p.precon2, v) - mul!(u3, p.A, -u1) - ldiv!(v, p.precon2, u3) - v .+= u1 - v .+= u2 + 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 diff --git a/test/alltests.jl b/test/alltests.jl index 0204826..99b4e84 100644 --- a/test/alltests.jl +++ b/test/alltests.jl @@ -108,3 +108,7 @@ end @testset "Block" begin include("test_block.jl") end + +@testset "Product" begin + include("test_prod.jl") +end diff --git a/test/test_prod.jl b/test/test_prod.jl index 27c5af4..d4e9af7 100644 --- a/test/test_prod.jl +++ b/test/test_prod.jl @@ -7,12 +7,14 @@ using LinearAlgebra using Random function main(; n = 100) - Random.seed!(2026) + 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 @@ -25,29 +27,39 @@ function main(; n = 100) 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 + @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 + @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 + @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 + @test sol6 ≈ sol0 rtol = 1.0e-5 + @test hist6.iters < hist5.iters return From 07c6726a3700d78e5835e8893a9ad89cfc9a2a4f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Tue, 16 Jun 2026 21:24:22 +0200 Subject: [PATCH 3/3] update changelog --- CHANGELOG.md | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index a150f63..5ca57a8 100644 --- a/CHANGELOG.md +++ b/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 `\`.