How can I compute the Schur complement:
$$ S = K_{bb} - K_{ba} K_{aa}^{-1} K_{ab} $$
where
$$ K=\begin{pmatrix} K_{aa} & K_{ab} \\ K_{ba} & K_{bb} \end{pmatrix} $$
(in some ordering) is a PETSc matrix (Mat
)?
It is very expensive to compute the Schur complement of a matrix and very rarely needed in practice. PETSc highly recommends avoiding algorithms that need it. The Schur complement of a matrix (dense or sparse) is essentially always dense, so begin by:
MatLUFactor()
or MatCholeskyFactor()
, or use MatGetFactor()
followed by MatLUFactorSymbolic()
followed by MatLUFactorNumeric()
if you wish to use and external solver package like SuperLU_Dist. Call the result A
.MatMatSolve(A,Kba,T)
.MatMatMult(Kab,T,MAT_INITIAL_MATRIX,1.0,&S)
.MatAXPY(S,-1.0,Kbb,MAT_SUBSET_NONZERO)
.MatScale(S,-1.0)
For computing Schur complements like this it does not make sense to use the KSP
iterative solvers since for solving many moderate size problems using a direct factorization is much faster than iterative solvers. As you can see, this requires a great deal of work space and computation so is best avoided. However, it is not necessary to assemble the Schur complement $S$ in order to solve systems with it. Use MatCreateSchurComplement(Kaa,Kaa_pre,Kab,Kba,Kbb,&S)
to create a matrix that applies the action of S
(using Kaa_pre
to solve with Kaa
), but does not assemble.
Alternatively, if you already have a block matrix $$ K = \begin{pmatrix} K_{aa} & K_{ab} \\ K_{ba} & K_{bb} \end{pmatrix} $$ (in some ordering), then you can create index sets (IS
) isa and isb to address each block, then use MatGetSchurComplement()
to create the Schur complement and/or an approximation suitable for preconditioning. Since $S$ is generally dense, standard preconditioning methods cannot typically be applied directly to Schur complements. There are many approaches to preconditioning Schur complements including using the SIMPLE approximation $ K_{bb} - K_{ba} \operatorname{diag}(K_{aa})^{-1} K_{ab} $ to create a sparse matrix that approximates the Schur complement (this is returned by default for the optional "preconditioning" matrix in MatGetSchurComplement()
).
Another alternative is to interpret the matrices as differential operators and apply approximate commutator arguments to find a spectrally equivalent operation that can be applied efficiently (see the "PCD" preconditioners from Elman, Silvester, and Wathen). A variant of this is the least squares commutator, which is closely related to the Moore-Penrose pseudoinverse, and is available in PCLSC
which operates on matrices of type MATSCHURCOMPLEMENT
.