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
@@ -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
 
 

src/preconbuilders.jl

@@ -294,3 +294,112 @@ 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.
+
+Creates a preconditioner of the form
+```
+[ ≈ A       0      ]
+[  0   ≈ -α BᵀA⁻¹B ]
+```
+
+where α is 1 by default (can be adjusted by `S_factor`).
+
+In the Schur factor, we use `Diagonal(A)` for efficiency.
+
+# 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 full `lu`)
+- `S_factor`: An additional scalar factor for the Schur complement block. Can be used to obtain positive definite preconditioners
+
+Note: All factorizations need to be able for operate on vector views, as
+
+Returns a function `prec(M, p)` that creates a `SchurComplementPreconditioner`.
+"""
+function SchurComplementPreconBuilder(dofs_first_block, A_factorization, S_factorization = lu; verbosity = 0, S_factor = 1.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 = -S_factor * 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
+    # do it in parallel
+    t1 = Threads.@spawn @views ldiv!(u[part1], p.A_fac, v[part1])
+    t2 = Threads.@spawn @views ldiv!(u[part2], p.S_fac, v[part2])
+    fetch.([t1, t2])
+
+    return u
+end
+
+function LinearAlgebra.ldiv!(p::SchurComplementPreconditioner, v)
+
+    (part1, part2) = p.partitions
+    # do it in parallel
+    t1 = Threads.@spawn @views ldiv!(p.A_fac, v[part1])
+    t2 = Threads.@spawn @views ldiv!(p.S_fac, v[part2])
+    fetch.([t1, t2])
+
+    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,20 @@ 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
+
+        # S_factor = -1.0 to make the precon SPD
+        sol = cg(M, c, Pl = SchurComplementPreconBuilder(n^2, lu, S_factor = -1.0)(M))
+        @test sol ≈ sol1
+    end
+
     return
 end