Add Schur Complement Preconditioner

Author: pjaap

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 `\`.

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 <juergen.fuhrmann@wias-berlin.de>", "Daniel Runge"]
 
 [deps]

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

src/preconbuilders.jl

@@ -294,3 +294,98 @@ 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

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