diff --git a/CHANGELOG.md b/CHANGELOG.md index a150f63..c672e5b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,9 @@ # Changelog +## [2.1.0] + +- Add `SchurComplementPreconditioner` and `SchurComplementPreconBuilder` + ## [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 6edec0d..592c6c4 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "ExtendableSparse" uuid = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3" -version = "2.0.1" +version = "2.1.0" authors = ["Juergen Fuhrmann ", "Daniel Runge"] [deps] diff --git a/src/ExtendableSparse.jl b/src/ExtendableSparse.jl index 00ddbd7..49ca0ad 100644 --- a/src/ExtendableSparse.jl +++ b/src/ExtendableSparse.jl @@ -7,7 +7,7 @@ module ExtendableSparse using DocStringExtensions: DocStringExtensions, SIGNATURES, TYPEDEF, TYPEDFIELDS, TYPEDSIGNATURES using ILUZero: ILUZero -using LinearAlgebra: LinearAlgebra, Diagonal, Hermitian, Symmetric, Tridiagonal, convert, mul!, ldiv! +using LinearAlgebra: LinearAlgebra, Diagonal, Hermitian, Symmetric, Tridiagonal, convert, mul!, ldiv!, lu using SparseArrays: SparseArrays, AbstractSparseMatrix, AbstractSparseMatrixCSC, SparseMatrixCSC using SparseArrays: dropzeros!, findnz, nzrange, sparse, spzeros, rowvals, getcolptr, nonzeros, nnz using Sparspak: sparspaklu @@ -42,7 +42,7 @@ export eliminate_dirichlet, eliminate_dirichlet!, mark_dirichlet include("preconbuilders.jl") -export LinearSolvePreconBuilder, BlockPreconBuilder, JacobiPreconBuilder +export LinearSolvePreconBuilder, BlockPreconBuilder, JacobiPreconBuilder, SchurComplementPreconBuilder @public ILUZeroPreconBuilder, ILUTPreconBuilder diff --git a/src/preconbuilders.jl b/src/preconbuilders.jl index 93d4994..1b1e76e 100644 --- a/src/preconbuilders.jl +++ b/src/preconbuilders.jl @@ -294,3 +294,95 @@ LinearAlgebra.ldiv!(fact::JacobiPreconditioner, v) = ldiv!(fact.factorization, v allow_views(::JacobiPreconditioner) = true allow_views(::Type{JacobiPreconditioner}) = true + + +""" + SchurComplementPreconditioner{AT, ST} + +Preconditioner for saddle-point problems of the form: +``` +[ A B ] +[ Bᵀ 0 ] +``` + +# Fields +- `partitions`: Vector of two ranges defining matrix partitioning +- `A_fac::AT`: Factorization of the A-block (top-left) +- `S_fac::ST`: Factorization of the Schur complement (approximation of -BᵀA⁻¹B) +""" +struct SchurComplementPreconditioner{AT, ST} + partitions + A_fac::AT + S_fac::ST +end + + +""" + SchurComplementPreconBuilder(dofs_first_block, A_factorization, S_factorization = lu) + +Factory function creating preconditioners for saddle-point problems. + +# Arguments +- `dofs_first_block`: Number of degrees of freedom in the first block (size of A matrix) +- `A_factorization`: Factorization method for the A-block (e.g., `ilu0, Diagonal, lu`) +- `S_factorization`: Factorization method for the Schur complement (defaults to `lu`) + +Returns a function `prec(M, p)` that creates a `SchurComplementPreconditioner`. +""" +function SchurComplementPreconBuilder(dofs_first_block, A_factorization, S_factorization = lu; verbosity = 0) + + # this is the resulting preconditioner + function prec(M) + + # We have + # M = [ A B ] → n1 dofs + # [ Bᵀ 0 ] → n2 dofs + + n1 = dofs_first_block + n2 = size(M, 1) - n1 + + # I see no other way than creating this expensive copy + A = M[1:n1, 1:n1] + B = M[1:n1, (n1 + 1):end] + + # first factorization + A_fac = A_factorization(A) + + verbosity > 0 && @info "SchurComplementPreconBuilder: A ($n1×$n1) is factorized" + + # compute the Schur Matrix + # S ≈ -BᵀA⁻¹B + # we use the diagonal of A: this is _very_ performant and creates a sparse result + S = - B' * (Diagonal(A) \ B) + + verbosity > 0 && @info "SchurComplementPreconBuilder: S ($n2×$n2) is computed" + + # factorize S + S_fac = S_factorization(S) + + verbosity > 0 && @info "SchurComplementPreconBuilder: S ($n2×$n2) is factorized" + + return SchurComplementPreconditioner([1:n1, (n1 + 1):(n1 + n2)], A_fac, S_fac) + end + + return prec +end + + +function LinearAlgebra.ldiv!(u, p::SchurComplementPreconditioner, v) + + (part1, part2) = p.partitions + @views ldiv!(u[part1], p.A_fac, v[part1]) + @views ldiv!(u[part2], p.S_fac, v[part2]) + + return u +end + +function LinearAlgebra.ldiv!(p::SchurComplementPreconditioner, v) + + (part1, part2) = p.partitions + @views ldiv!(p.A_fac, v[part1]) + @views ldiv!(p.S_fac, v[part2]) + + return v +end diff --git a/test/test_block.jl b/test/test_block.jl index f147cb3..4b53240 100644 --- a/test/test_block.jl +++ b/test/test_block.jl @@ -1,12 +1,13 @@ module test_block -using Test +using AMGCLWrap using ExtendableSparse -using ExtendableSparse: BlockPreconditioner, jacobi +using ExtendableSparse: BlockPreconditioner, jacobi, SchurComplementPreconBuilder using ILUZero, AlgebraicMultigrid using IterativeSolvers using LinearAlgebra +using SparseArrays using Sparspak -using AMGCLWrap +using Test ExtendableSparse.allow_views(::typeof(ilu0)) = true @@ -29,6 +30,19 @@ function main(; n = 100) sol = cg(A, b, Pl = BlockPreconditioner(A; partitioning, factorization = sparspaklu)) @test sol ≈ sol0 + # Schur complement: create a saddle point system + let + m = n ÷ 10 + B = I[1:(n^2), 1:(m^2)] + M = [ A B; B' spzeros(m^2, m^2)] + + sol1 = ones(n^2 + m^2) + c = M * sol1 + + sol = cg(M, c, Pl = SchurComplementPreconBuilder(n^2, lu)(M)) + @test sol ≈ sol1 + end + return end