From e1dc36608dff6aaf2818a4ad5f49fcd61d11cce0 Mon Sep 17 00:00:00 2001 From: Patrick Jaap Date: Fri, 13 Feb 2026 14:08:54 +0100 Subject: [PATCH 1/6] Add Schur Complement Preconditioner --- CHANGELOG.md | 4 ++ Project.toml | 2 +- src/ExtendableSparse.jl | 2 +- src/preconbuilders.jl | 94 +++++++++++++++++++++++++++++++++++++++++ test/test_block.jl | 20 +++++++-- 5 files changed, 117 insertions(+), 5 deletions(-) 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..99fd287 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 diff --git a/src/preconbuilders.jl b/src/preconbuilders.jl index 93d4994..18515cb 100644 --- a/src/preconbuilders.jl +++ b/src/preconbuilders.jl @@ -294,3 +294,97 @@ 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) + + # this is the resulting preconditioner + function prec(M, p) + + # 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) + + # compute the (dense!) Schur Matrix + # S ≈ - Bᵀ A⁻¹ B + S = zeros(n2, n2) + + @info "Computing the Schur complement..." + ldiv_buffer = zeros(n1) + for col in 1:n2 + # TODO: the following is not thread parallel? + @views ldiv!(ldiv_buffer, A_fac, B[:, col]) + @views mul!(S[:, col], B', ldiv_buffer) + end + S .*= -1.0 + @info "...done" + + S_fac = S_factorization(sparse(S)) + + 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..0532863 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, nothing)) + @test sol ≈ sol1 + end + return end From d83841db566a0adde3a7e67362038af718c9ece4 Mon Sep 17 00:00:00 2001 From: Patrick Jaap Date: Fri, 13 Feb 2026 15:35:44 +0100 Subject: [PATCH 2/6] export SchurComplementPreconBuilder --- src/ExtendableSparse.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ExtendableSparse.jl b/src/ExtendableSparse.jl index 99fd287..49ca0ad 100644 --- a/src/ExtendableSparse.jl +++ b/src/ExtendableSparse.jl @@ -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 From f4e4fe7f1ec15285dd4c8b7b658357ac5c413817 Mon Sep 17 00:00:00 2001 From: Patrick Jaap Date: Mon, 16 Feb 2026 16:44:27 +0100 Subject: [PATCH 3/6] SchurComplement: compute S in parallel --- Project.toml | 1 + src/ExtendableSparse.jl | 1 + src/preconbuilders.jl | 41 +++++++++++++++++++++++++++++++---------- 3 files changed, 33 insertions(+), 10 deletions(-) diff --git a/Project.toml b/Project.toml index 592c6c4..664808c 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ version = "2.1.0" authors = ["Juergen Fuhrmann ", "Daniel Runge"] [deps] +ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" ILUZero = "88f59080-6952-5380-9ea5-54057fb9a43f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/src/ExtendableSparse.jl b/src/ExtendableSparse.jl index 49ca0ad..ae86a03 100644 --- a/src/ExtendableSparse.jl +++ b/src/ExtendableSparse.jl @@ -5,6 +5,7 @@ $(read(joinpath(@__DIR__, "..", "README.md"), String)) """ module ExtendableSparse +using ChunkSplitters: chunks using DocStringExtensions: DocStringExtensions, SIGNATURES, TYPEDEF, TYPEDFIELDS, TYPEDSIGNATURES using ILUZero: ILUZero using LinearAlgebra: LinearAlgebra, Diagonal, Hermitian, Symmetric, Tridiagonal, convert, mul!, ldiv!, lu diff --git a/src/preconbuilders.jl b/src/preconbuilders.jl index 18515cb..5fecd8c 100644 --- a/src/preconbuilders.jl +++ b/src/preconbuilders.jl @@ -348,21 +348,42 @@ function SchurComplementPreconBuilder(dofs_first_block, A_factorization, S_facto # first factorization A_fac = A_factorization(A) - # compute the (dense!) Schur Matrix + # compute the (dense!) Schur Matrix in parallel # S ≈ - Bᵀ A⁻¹ B + + function col_loop(col_chunk) + + # sigh, some factorizations (like ilu0) store internal buffers. Hence, every thread needs a copy. + local A_fac_copy = deepcopy(A_fac) + + # local buffers and result + local S_chunk = zeros(n2, length(col_chunk)) + local ldiv_buffer = zeros(n1) + + for (icol, col) in enumerate(col_chunk) + @views ldiv!(ldiv_buffer, A_fac_copy, B[:, col]) + @views mul!(S_chunk[:, icol], B', ldiv_buffer) + end + S_chunk .*= -1.0 + + return S_chunk + end + + # split columns into chunks and assemble S + col_chunks = chunks(1:n2, n = Threads.nthreads()) + tasks = map(col_chunks) do col_chunk + Threads.@spawn col_loop(col_chunk) + end + S_chunks = fetch.(tasks) + S = zeros(n2, n2) - @info "Computing the Schur complement..." - ldiv_buffer = zeros(n1) - for col in 1:n2 - # TODO: the following is not thread parallel? - @views ldiv!(ldiv_buffer, A_fac, B[:, col]) - @views mul!(S[:, col], B', ldiv_buffer) + for (col_chunk, S_chunk) in zip(col_chunks, S_chunks) + S[:, col_chunk] .= S_chunk end - S .*= -1.0 - @info "...done" - S_fac = S_factorization(sparse(S)) + # factorize S + S_fac = S_factorization(S) return SchurComplementPreconditioner([1:n1, (n1 + 1):(n1 + n2)], A_fac, S_fac) end From f6957982ce18534faae19551989f37225a5ad26b Mon Sep 17 00:00:00 2001 From: Patrick Jaap Date: Tue, 17 Feb 2026 17:16:54 +0100 Subject: [PATCH 4/6] SchurComplementPrecon: Use Diagonal(A) for the Schur block --- src/preconbuilders.jl | 39 +++++---------------------------------- test/test_block.jl | 2 +- 2 files changed, 6 insertions(+), 35 deletions(-) diff --git a/src/preconbuilders.jl b/src/preconbuilders.jl index 5fecd8c..7b6dd0e 100644 --- a/src/preconbuilders.jl +++ b/src/preconbuilders.jl @@ -332,7 +332,7 @@ Returns a function `prec(M, p)` that creates a `SchurComplementPreconditioner`. function SchurComplementPreconBuilder(dofs_first_block, A_factorization, S_factorization = lu) # this is the resulting preconditioner - function prec(M, p) + function prec(M) # We have # M = [ A B ] → n1 dofs @@ -348,39 +348,10 @@ function SchurComplementPreconBuilder(dofs_first_block, A_factorization, S_facto # first factorization A_fac = A_factorization(A) - # compute the (dense!) Schur Matrix in parallel - # S ≈ - Bᵀ A⁻¹ B - - function col_loop(col_chunk) - - # sigh, some factorizations (like ilu0) store internal buffers. Hence, every thread needs a copy. - local A_fac_copy = deepcopy(A_fac) - - # local buffers and result - local S_chunk = zeros(n2, length(col_chunk)) - local ldiv_buffer = zeros(n1) - - for (icol, col) in enumerate(col_chunk) - @views ldiv!(ldiv_buffer, A_fac_copy, B[:, col]) - @views mul!(S_chunk[:, icol], B', ldiv_buffer) - end - S_chunk .*= -1.0 - - return S_chunk - end - - # split columns into chunks and assemble S - col_chunks = chunks(1:n2, n = Threads.nthreads()) - tasks = map(col_chunks) do col_chunk - Threads.@spawn col_loop(col_chunk) - end - S_chunks = fetch.(tasks) - - S = zeros(n2, n2) - - for (col_chunk, S_chunk) in zip(col_chunks, S_chunks) - S[:, col_chunk] .= S_chunk - end + # 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) # factorize S S_fac = S_factorization(S) diff --git a/test/test_block.jl b/test/test_block.jl index 0532863..4b53240 100644 --- a/test/test_block.jl +++ b/test/test_block.jl @@ -39,7 +39,7 @@ function main(; n = 100) sol1 = ones(n^2 + m^2) c = M * sol1 - sol = cg(M, c, Pl = SchurComplementPreconBuilder(n^2, lu)(M, nothing)) + sol = cg(M, c, Pl = SchurComplementPreconBuilder(n^2, lu)(M)) @test sol ≈ sol1 end From 1049927b10b17ed488ea76798c771b395d46956e Mon Sep 17 00:00:00 2001 From: Patrick Jaap Date: Wed, 18 Feb 2026 10:02:49 +0100 Subject: [PATCH 5/6] deps: remove Chunksplitters --- Project.toml | 1 - src/ExtendableSparse.jl | 1 - 2 files changed, 2 deletions(-) diff --git a/Project.toml b/Project.toml index 664808c..592c6c4 100644 --- a/Project.toml +++ b/Project.toml @@ -4,7 +4,6 @@ version = "2.1.0" authors = ["Juergen Fuhrmann ", "Daniel Runge"] [deps] -ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" ILUZero = "88f59080-6952-5380-9ea5-54057fb9a43f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/src/ExtendableSparse.jl b/src/ExtendableSparse.jl index ae86a03..49ca0ad 100644 --- a/src/ExtendableSparse.jl +++ b/src/ExtendableSparse.jl @@ -5,7 +5,6 @@ $(read(joinpath(@__DIR__, "..", "README.md"), String)) """ module ExtendableSparse -using ChunkSplitters: chunks using DocStringExtensions: DocStringExtensions, SIGNATURES, TYPEDEF, TYPEDFIELDS, TYPEDSIGNATURES using ILUZero: ILUZero using LinearAlgebra: LinearAlgebra, Diagonal, Hermitian, Symmetric, Tridiagonal, convert, mul!, ldiv!, lu From e3a8e219e907cb31e44b379453448d2073803862 Mon Sep 17 00:00:00 2001 From: Patrick Jaap Date: Wed, 18 Feb 2026 14:42:04 +0100 Subject: [PATCH 6/6] SchurComplementPrecon: add verbosity kwarg --- src/preconbuilders.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/preconbuilders.jl b/src/preconbuilders.jl index 7b6dd0e..1b1e76e 100644 --- a/src/preconbuilders.jl +++ b/src/preconbuilders.jl @@ -329,7 +329,7 @@ Factory function creating preconditioners for saddle-point problems. Returns a function `prec(M, p)` that creates a `SchurComplementPreconditioner`. """ -function SchurComplementPreconBuilder(dofs_first_block, A_factorization, S_factorization = lu) +function SchurComplementPreconBuilder(dofs_first_block, A_factorization, S_factorization = lu; verbosity = 0) # this is the resulting preconditioner function prec(M) @@ -348,14 +348,20 @@ function SchurComplementPreconBuilder(dofs_first_block, A_factorization, S_facto # 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