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 `\`. 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/ExtendableSparse.jl b/src/ExtendableSparse.jl index 6857d96..cc55a8e 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..79dbb86 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,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 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 new file mode 100644 index 0000000..d4e9af7 --- /dev/null +++ b/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