8
$\begingroup$

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)?

$\endgroup$
13
$\begingroup$

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:

  • forming a dense matrix $K_{ba}$,
  • also create another dense matrix $T$ of the same size.
  • Then either factor the matrix $K_{aa}$ directly with 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.
  • Then call MatMatSolve(A,Kba,T).
  • Then call MatMatMult(Kab,T,MAT_INITIAL_MATRIX,1.0,&S).
  • Now call MatAXPY(S,-1.0,Kbb,MAT_SUBSET_NONZERO).
  • Followed by 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.

$\endgroup$
  • 2
    $\begingroup$ At the risk of reopening a settled question, I think this answer is too pessimistic. I think computing a schur complement can be a reasonable thing to do in some circumstances (substructuring/domain decomposition in particular) and there are other algorithms that can compute them with lower time complexity than what's described here. See scicomp.stackexchange.com/questions/28934 for inspiration. $\endgroup$ – rchilton1980 Sep 17 '18 at 16:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.