Actual source code: schurm.c

  1: #include <../src/ksp/ksp/utils/schurm/schurm.h>

  3: const char *const MatSchurComplementAinvTypes[] = {"DIAG", "LUMP", "BLOCKDIAG", "FULL", "MatSchurComplementAinvType", "MAT_SCHUR_COMPLEMENT_AINV_", NULL};

  5: PetscErrorCode MatCreateVecs_SchurComplement(Mat N, Vec *right, Vec *left)
  6: {
  7:   Mat_SchurComplement *Na = (Mat_SchurComplement *)N->data;

  9:   PetscFunctionBegin;
 10:   if (Na->D) {
 11:     PetscCall(MatCreateVecs(Na->D, right, left));
 12:     PetscFunctionReturn(PETSC_SUCCESS);
 13:   }
 14:   if (right) PetscCall(MatCreateVecs(Na->B, right, NULL));
 15:   if (left) PetscCall(MatCreateVecs(Na->C, NULL, left));
 16:   PetscFunctionReturn(PETSC_SUCCESS);
 17: }

 19: PetscErrorCode MatView_SchurComplement(Mat N, PetscViewer viewer)
 20: {
 21:   Mat_SchurComplement *Na = (Mat_SchurComplement *)N->data;

 23:   PetscFunctionBegin;
 24:   PetscCall(PetscViewerASCIIPrintf(viewer, "Schur complement A11 - A10 inv(A00) A01\n"));
 25:   if (Na->D) {
 26:     PetscCall(PetscViewerASCIIPrintf(viewer, "A11\n"));
 27:     PetscCall(PetscViewerASCIIPushTab(viewer));
 28:     PetscCall(MatView(Na->D, viewer));
 29:     PetscCall(PetscViewerASCIIPopTab(viewer));
 30:   } else {
 31:     PetscCall(PetscViewerASCIIPrintf(viewer, "A11 = 0\n"));
 32:   }
 33:   PetscCall(PetscViewerASCIIPrintf(viewer, "A10\n"));
 34:   PetscCall(PetscViewerASCIIPushTab(viewer));
 35:   PetscCall(MatView(Na->C, viewer));
 36:   PetscCall(PetscViewerASCIIPopTab(viewer));
 37:   PetscCall(PetscViewerASCIIPrintf(viewer, "KSP of A00\n"));
 38:   PetscCall(PetscViewerASCIIPushTab(viewer));
 39:   PetscCall(KSPView(Na->ksp, viewer));
 40:   PetscCall(PetscViewerASCIIPopTab(viewer));
 41:   PetscCall(PetscViewerASCIIPrintf(viewer, "A01\n"));
 42:   PetscCall(PetscViewerASCIIPushTab(viewer));
 43:   PetscCall(MatView(Na->B, viewer));
 44:   PetscCall(PetscViewerASCIIPopTab(viewer));
 45:   PetscFunctionReturn(PETSC_SUCCESS);
 46: }

 48: /*
 49:            A11^T - A01^T ksptrans(A00,Ap00) A10^T
 50: */
 51: PetscErrorCode MatMultTranspose_SchurComplement(Mat N, Vec x, Vec y)
 52: {
 53:   Mat_SchurComplement *Na = (Mat_SchurComplement *)N->data;

 55:   PetscFunctionBegin;
 56:   if (!Na->work1) PetscCall(MatCreateVecs(Na->A, &Na->work1, NULL));
 57:   if (!Na->work2) PetscCall(MatCreateVecs(Na->A, &Na->work2, NULL));
 58:   PetscCall(MatMultTranspose(Na->C, x, Na->work1));
 59:   PetscCall(KSPSolveTranspose(Na->ksp, Na->work1, Na->work2));
 60:   PetscCall(MatMultTranspose(Na->B, Na->work2, y));
 61:   PetscCall(VecScale(y, -1.0));
 62:   if (Na->D) PetscCall(MatMultTransposeAdd(Na->D, x, y, y));
 63:   PetscFunctionReturn(PETSC_SUCCESS);
 64: }

 66: /*
 67:            A11 - A10 ksp(A00,Ap00) A01
 68: */
 69: PetscErrorCode MatMult_SchurComplement(Mat N, Vec x, Vec y)
 70: {
 71:   Mat_SchurComplement *Na = (Mat_SchurComplement *)N->data;

 73:   PetscFunctionBegin;
 74:   if (!Na->work1) PetscCall(MatCreateVecs(Na->A, &Na->work1, NULL));
 75:   if (!Na->work2) PetscCall(MatCreateVecs(Na->A, &Na->work2, NULL));
 76:   PetscCall(MatMult(Na->B, x, Na->work1));
 77:   PetscCall(KSPSolve(Na->ksp, Na->work1, Na->work2));
 78:   PetscCall(MatMult(Na->C, Na->work2, y));
 79:   PetscCall(VecScale(y, -1.0));
 80:   if (Na->D) PetscCall(MatMultAdd(Na->D, x, y, y));
 81:   PetscFunctionReturn(PETSC_SUCCESS);
 82: }

 84: /*
 85:            A11 - A10 ksp(A00,Ap00) A01
 86: */
 87: PetscErrorCode MatMultAdd_SchurComplement(Mat N, Vec x, Vec y, Vec z)
 88: {
 89:   Mat_SchurComplement *Na = (Mat_SchurComplement *)N->data;

 91:   PetscFunctionBegin;
 92:   if (!Na->work1) PetscCall(MatCreateVecs(Na->A, &Na->work1, NULL));
 93:   if (!Na->work2) PetscCall(MatCreateVecs(Na->A, &Na->work2, NULL));
 94:   PetscCall(MatMult(Na->B, x, Na->work1));
 95:   PetscCall(KSPSolve(Na->ksp, Na->work1, Na->work2));
 96:   if (y == z) {
 97:     PetscCall(VecScale(Na->work2, -1.0));
 98:     PetscCall(MatMultAdd(Na->C, Na->work2, z, z));
 99:   } else {
100:     PetscCall(MatMult(Na->C, Na->work2, z));
101:     PetscCall(VecAYPX(z, -1.0, y));
102:   }
103:   if (Na->D) PetscCall(MatMultAdd(Na->D, x, z, z));
104:   PetscFunctionReturn(PETSC_SUCCESS);
105: }

107: PetscErrorCode MatSetFromOptions_SchurComplement(Mat N, PetscOptionItems *PetscOptionsObject)
108: {
109:   Mat_SchurComplement *Na = (Mat_SchurComplement *)N->data;

111:   PetscFunctionBegin;
112:   PetscOptionsHeadBegin(PetscOptionsObject, "MatSchurComplementOptions");
113:   Na->ainvtype = MAT_SCHUR_COMPLEMENT_AINV_DIAG;
114:   PetscCall(PetscOptionsEnum("-mat_schur_complement_ainv_type", "Type of approximation for DIAGFORM(A00) used when assembling Sp = A11 - A10 inv(DIAGFORM(A00)) A01", "MatSchurComplementSetAinvType", MatSchurComplementAinvTypes, (PetscEnum)Na->ainvtype,
115:                              (PetscEnum *)&Na->ainvtype, NULL));
116:   PetscOptionsHeadEnd();
117:   PetscCall(KSPSetFromOptions(Na->ksp));
118:   PetscFunctionReturn(PETSC_SUCCESS);
119: }

121: PetscErrorCode MatDestroy_SchurComplement(Mat N)
122: {
123:   Mat_SchurComplement *Na = (Mat_SchurComplement *)N->data;

125:   PetscFunctionBegin;
126:   PetscCall(MatDestroy(&Na->A));
127:   PetscCall(MatDestroy(&Na->Ap));
128:   PetscCall(MatDestroy(&Na->B));
129:   PetscCall(MatDestroy(&Na->C));
130:   PetscCall(MatDestroy(&Na->D));
131:   PetscCall(VecDestroy(&Na->work1));
132:   PetscCall(VecDestroy(&Na->work2));
133:   PetscCall(KSPDestroy(&Na->ksp));
134:   PetscCall(PetscFree(N->data));
135:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_schurcomplement_seqdense_C", NULL));
136:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_schurcomplement_mpidense_C", NULL));
137:   PetscFunctionReturn(PETSC_SUCCESS);
138: }

140: /*@
141:       MatCreateSchurComplement - Creates a new `Mat` that behaves like the Schur complement of a matrix

143:    Collective

145:    Input Parameters:
146: +   A00  - the upper-left block of the original matrix A = [A00 A01; A10 A11]
147: .   Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A00^{-1}
148: .   A01  - the upper-right block of the original matrix A = [A00 A01; A10 A11]
149: .   A10  - the lower-left block of the original matrix A = [A00 A01; A10 A11]
150: -   A11  - (optional) the lower-right block of the original matrix A = [A00 A01; A10 A11]

152:    Output Parameter:
153: .   S - the matrix that behaves as the Schur complement S = A11 - A10 ksp(A00,Ap00) A01

155:    Level: intermediate

157:    Notes:
158:     The Schur complement is NOT explicitly formed! Rather, this function returns a virtual Schur complement
159:     that can compute the matrix-vector product by using formula S = A11 - A10 A^{-1} A01
160:     for Schur complement S and a `KSP` solver to approximate the action of A^{-1}.

162:     All four matrices must have the same MPI communicator.

164:     `A00` and  `A11` must be square matrices.

166:     `MatGetSchurComplement()` takes as arguments the index sets for the submatrices and returns both the virtual Schur complement (what this returns) plus
167:     a sparse approximation to the Schur complement (useful for building a preconditioner for the Schur complement) which can be obtained from this
168:     matrix with `MatSchurComplementGetPmat()`

170:     Developer Note:
171:     The API that includes `MatGetSchurComplement()`, `MatCreateSchurComplement()`, `MatSchurComplementGetPmat()` should be refactored to
172:     remove redundancy and be clearer and simpler.

174: .seealso: [](ch_ksp), `MatCreateNormal()`, `MatMult()`, `MatCreate()`, `MatSchurComplementGetKSP()`, `MatSchurComplementUpdateSubMatrices()`, `MatCreateTranspose()`, `MatGetSchurComplement()`,
175:           `MatSchurComplementGetPmat()`, `MatSchurComplementSetSubMatrices()`
176: @*/
177: PetscErrorCode MatCreateSchurComplement(Mat A00, Mat Ap00, Mat A01, Mat A10, Mat A11, Mat *S)
178: {
179:   PetscFunctionBegin;
180:   PetscCall(KSPInitializePackage());
181:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A00), S));
182:   PetscCall(MatSetType(*S, MATSCHURCOMPLEMENT));
183:   PetscCall(MatSchurComplementSetSubMatrices(*S, A00, Ap00, A01, A10, A11));
184:   PetscFunctionReturn(PETSC_SUCCESS);
185: }

187: /*@
188:       MatSchurComplementSetSubMatrices - Sets the matrices that define the Schur complement

190:    Collective

192:    Input Parameters:
193: +   S                - matrix obtained with `MatSetType`(S,`MATSCHURCOMPLEMENT`)
194: +   A00  - the upper-left block of the original matrix A = [A00 A01; A10 A11]
195: .   Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A00^{-1}
196: .   A01  - the upper-right block of the original matrix A = [A00 A01; A10 A11]
197: .   A10  - the lower-left block of the original matrix A = [A00 A01; A10 A11]
198: -   A11  - (optional) the lower-right block of the original matrix A = [A00 A01; A10 A11]

200:    Level: intermediate

202:    Notes:
203:      The Schur complement is NOT explicitly formed! Rather, this
204:      object performs the matrix-vector product of the Schur complement by using formula S = A11 - A10 ksp(A00,Ap00) A01

206:      All four matrices must have the same MPI communicator.

208:      `A00` and `A11` must be square matrices.

210:      This is to be used in the context of code such as
211: .vb
212:      MatSetType(S,MATSCHURCOMPLEMENT);
213:      MatSchurComplementSetSubMatrices(S,...);
214: .ve
215:     while `MatSchurComplementUpdateSubMatrices()` should only be called after `MatCreateSchurComplement()` or `MatSchurComplementSetSubMatrices()`

217: .seealso: [](ch_ksp), `Mat`, `MatCreateNormal()`, `MatMult()`, `MatCreate()`, `MatSchurComplementGetKSP()`, `MatSchurComplementUpdateSubMatrices()`, `MatCreateTranspose()`, `MatCreateSchurComplement()`, `MatGetSchurComplement()`
218: @*/
219: PetscErrorCode MatSchurComplementSetSubMatrices(Mat S, Mat A00, Mat Ap00, Mat A01, Mat A10, Mat A11)
220: {
221:   Mat_SchurComplement *Na = (Mat_SchurComplement *)S->data;
222:   PetscBool            isschur;

224:   PetscFunctionBegin;
225:   PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
226:   if (!isschur) PetscFunctionReturn(PETSC_SUCCESS);
227:   PetscCheck(!S->assembled, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONGSTATE, "Use MatSchurComplementUpdateSubMatrices() for already used matrix");
232:   PetscCheckSameComm(A00, 2, Ap00, 3);
233:   PetscCheckSameComm(A00, 2, A01, 4);
234:   PetscCheckSameComm(A00, 2, A10, 5);
235:   PetscCheck(A00->rmap->n == A00->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local rows of A00 %" PetscInt_FMT " do not equal local columns %" PetscInt_FMT, A00->rmap->n, A00->cmap->n);
236:   PetscCheck(A00->rmap->n == Ap00->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local rows of A00 %" PetscInt_FMT " do not equal local rows of Ap00 %" PetscInt_FMT, A00->rmap->n, Ap00->rmap->n);
237:   PetscCheck(Ap00->rmap->n == Ap00->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local rows of Ap00 %" PetscInt_FMT " do not equal local columns %" PetscInt_FMT, Ap00->rmap->n, Ap00->cmap->n);
238:   PetscCheck(A00->cmap->n == A01->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local columns of A00 %" PetscInt_FMT " do not equal local rows of A01 %" PetscInt_FMT, A00->cmap->n, A01->rmap->n);
239:   PetscCheck(A10->cmap->n == A00->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local columns of A10 %" PetscInt_FMT " do not equal local rows of A00 %" PetscInt_FMT, A10->cmap->n, A00->rmap->n);
240:   if (A11) {
242:     PetscCheckSameComm(A00, 2, A11, 6);
243:     PetscCheck(A10->rmap->n == A11->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local rows of A10 %" PetscInt_FMT " do not equal local rows A11 %" PetscInt_FMT, A10->rmap->n, A11->rmap->n);
244:   }

246:   PetscCall(MatSetSizes(S, A10->rmap->n, A01->cmap->n, A10->rmap->N, A01->cmap->N));
247:   PetscCall(PetscObjectReference((PetscObject)A00));
248:   PetscCall(PetscObjectReference((PetscObject)Ap00));
249:   PetscCall(PetscObjectReference((PetscObject)A01));
250:   PetscCall(PetscObjectReference((PetscObject)A10));
251:   Na->A  = A00;
252:   Na->Ap = Ap00;
253:   Na->B  = A01;
254:   Na->C  = A10;
255:   Na->D  = A11;
256:   if (A11) PetscCall(PetscObjectReference((PetscObject)A11));
257:   PetscCall(MatSetUp(S));
258:   PetscCall(KSPSetOperators(Na->ksp, A00, Ap00));
259:   S->assembled = PETSC_TRUE;
260:   PetscFunctionReturn(PETSC_SUCCESS);
261: }

263: /*@
264:   MatSchurComplementGetKSP - Gets the `KSP` object that is used to solve with A00 in the Schur complement matrix S = A11 - A10 ksp(A00,Ap00) A01

266:   Not Collective

268:   Input Parameter:
269: . S - matrix obtained with `MatCreateSchurComplement()` (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01

271:   Output Parameter:
272: . ksp - the linear solver object

274:   Options Database Key:
275: . -fieldsplit_<splitname_0>_XXX sets KSP and PC options for the 0-split solver inside the Schur complement used in `PCFIELDSPLIT`; default <splitname_0> is 0.

277:   Level: intermediate

279: .seealso: [](ch_ksp), `Mat`, `MatSchurComplementSetKSP()`, `MatCreateSchurComplement()`, `MatCreateNormal()`, `MatMult()`, `MatCreate()`
280: @*/
281: PetscErrorCode MatSchurComplementGetKSP(Mat S, KSP *ksp)
282: {
283:   Mat_SchurComplement *Na;
284:   PetscBool            isschur;

286:   PetscFunctionBegin;
288:   PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
289:   PetscCheck(isschur, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)S)->type_name);
291:   Na   = (Mat_SchurComplement *)S->data;
292:   *ksp = Na->ksp;
293:   PetscFunctionReturn(PETSC_SUCCESS);
294: }

296: /*@
297:   MatSchurComplementSetKSP - Sets the `KSP` object that is used to solve with A00 in the Schur complement matrix S = A11 - A10 ksp(A00,Ap00) A01

299:   Not Collective

301:   Input Parameters:
302: + S   - matrix created with `MatCreateSchurComplement()`
303: - ksp - the linear solver object

305:   Level: developer

307:   Developer Note:
308:     This is used in `PCFIELDSPLIT` to reuse the 0-split `KSP` to implement ksp(A00,Ap00) in S.
309:     The `KSP` operators are overwritten with A00 and Ap00 currently set in S.

311: .seealso: [](ch_ksp), `Mat`, `MatSchurComplementGetKSP()`, `MatCreateSchurComplement()`, `MatCreateNormal()`, `MatMult()`, `MatCreate()`, `MATSCHURCOMPLEMENT`
312: @*/
313: PetscErrorCode MatSchurComplementSetKSP(Mat S, KSP ksp)
314: {
315:   Mat_SchurComplement *Na;
316:   PetscBool            isschur;

318:   PetscFunctionBegin;
320:   PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
321:   if (!isschur) PetscFunctionReturn(PETSC_SUCCESS);
323:   Na = (Mat_SchurComplement *)S->data;
324:   PetscCall(PetscObjectReference((PetscObject)ksp));
325:   PetscCall(KSPDestroy(&Na->ksp));
326:   Na->ksp = ksp;
327:   PetscCall(KSPSetOperators(Na->ksp, Na->A, Na->Ap));
328:   PetscFunctionReturn(PETSC_SUCCESS);
329: }

331: /*@
332:       MatSchurComplementUpdateSubMatrices - Updates the Schur complement matrix object with new submatrices

334:    Collective

336:    Input Parameters:
337: +   S                - matrix obtained with `MatCreateSchurComplement()` (or `MatSchurSetSubMatrices()`) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
338: .   A00  - the upper-left block of the original matrix A = [A00 A01; A10 A11]
339: .   Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A00^{-1}
340: .   A01  - the upper-right block of the original matrix A = [A00 A01; A10 A11]
341: .   A10  - the lower-left block of the original matrix A = [A00 A01; A10 A11]
342: -   A11  - (optional) the lower-right block of the original matrix A = [A00 A01; A10 A11]

344:    Level: intermediate

346:    Notes:
347:      All four matrices must have the same MPI communicator

349:      `A00` and  `A11` must be square matrices

351:      All of the matrices provided must have the same sizes as was used with `MatCreateSchurComplement()` or `MatSchurComplementSetSubMatrices()`
352:      though they need not be the same matrices.

354:      This can only be called after `MatCreateSchurComplement()` or `MatSchurComplementSetSubMatrices()`, it cannot be called immediately after `MatSetType`(S,`MATSCHURCOMPLEMENT`);

356:    Developer Note:
357:      This code is almost identical to `MatSchurComplementSetSubMatrices()`. The API should be refactored.

359: .seealso: [](ch_ksp), `Mat`, `MatCreateNormal()`, `MatMult()`, `MatCreate()`, `MatSchurComplementGetKSP()`, `MatCreateSchurComplement()`
360: @*/
361: PetscErrorCode MatSchurComplementUpdateSubMatrices(Mat S, Mat A00, Mat Ap00, Mat A01, Mat A10, Mat A11)
362: {
363:   Mat_SchurComplement *Na = (Mat_SchurComplement *)S->data;
364:   PetscBool            isschur;

366:   PetscFunctionBegin;
368:   PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
369:   if (!isschur) PetscFunctionReturn(PETSC_SUCCESS);
370:   PetscCheck(S->assembled, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONGSTATE, "Use MatSchurComplementSetSubMatrices() for a new matrix");
375:   PetscCheckSameComm(A00, 2, Ap00, 3);
376:   PetscCheckSameComm(A00, 2, A01, 4);
377:   PetscCheckSameComm(A00, 2, A10, 5);
378:   PetscCheck(A00->rmap->n == A00->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local rows of A00 %" PetscInt_FMT " do not equal local columns %" PetscInt_FMT, A00->rmap->n, A00->cmap->n);
379:   PetscCheck(A00->rmap->n == Ap00->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local rows of A00 %" PetscInt_FMT " do not equal local rows of Ap00 %" PetscInt_FMT, A00->rmap->n, Ap00->rmap->n);
380:   PetscCheck(Ap00->rmap->n == Ap00->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local rows of Ap00 %" PetscInt_FMT " do not equal local columns %" PetscInt_FMT, Ap00->rmap->n, Ap00->cmap->n);
381:   PetscCheck(A00->cmap->n == A01->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local columns of A00 %" PetscInt_FMT " do not equal local rows of A01 %" PetscInt_FMT, A00->cmap->n, A01->rmap->n);
382:   PetscCheck(A10->cmap->n == A00->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local columns of A10 %" PetscInt_FMT " do not equal local rows of A00 %" PetscInt_FMT, A10->cmap->n, A00->rmap->n);
383:   if (A11) {
385:     PetscCheckSameComm(A00, 2, A11, 6);
386:     PetscCheck(A10->rmap->n == A11->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local rows of A10 %" PetscInt_FMT " do not equal local rows A11 %" PetscInt_FMT, A10->rmap->n, A11->rmap->n);
387:   }

389:   PetscCall(PetscObjectReference((PetscObject)A00));
390:   PetscCall(PetscObjectReference((PetscObject)Ap00));
391:   PetscCall(PetscObjectReference((PetscObject)A01));
392:   PetscCall(PetscObjectReference((PetscObject)A10));
393:   if (A11) PetscCall(PetscObjectReference((PetscObject)A11));

395:   PetscCall(MatDestroy(&Na->A));
396:   PetscCall(MatDestroy(&Na->Ap));
397:   PetscCall(MatDestroy(&Na->B));
398:   PetscCall(MatDestroy(&Na->C));
399:   PetscCall(MatDestroy(&Na->D));

401:   Na->A  = A00;
402:   Na->Ap = Ap00;
403:   Na->B  = A01;
404:   Na->C  = A10;
405:   Na->D  = A11;

407:   PetscCall(KSPSetOperators(Na->ksp, A00, Ap00));
408:   PetscFunctionReturn(PETSC_SUCCESS);
409: }

411: /*@C
412:   MatSchurComplementGetSubMatrices - Get the individual submatrices in the Schur complement

414:   Collective

416:   Input Parameter:
417: . S    - matrix obtained with `MatCreateSchurComplement()` (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01

419:   Output Parameters:
420: + A00  - the upper-left block of the original matrix A = [A00 A01; A10 A11]
421: . Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}
422: . A01  - the upper-right block of the original matrix A = [A00 A01; A10 A11]
423: . A10  - the lower-left block of the original matrix A = [A00 A01; A10 A11]
424: - A11  - (optional) the lower-right block of the original matrix A = [A00 A01; A10 A11]

426:   Level: intermediate

428:   Note:
429:   `A11` is optional, and thus can be `NULL`.

431:   The reference counts of the submatrices are not increased before they are returned and the matrices should not be modified or destroyed.

433: .seealso: [](ch_ksp), `MatCreateNormal()`, `MatMult()`, `MatCreate()`, `MatSchurComplementGetKSP()`, `MatCreateSchurComplement()`, `MatSchurComplementUpdateSubMatrices()`
434: @*/
435: PetscErrorCode MatSchurComplementGetSubMatrices(Mat S, Mat *A00, Mat *Ap00, Mat *A01, Mat *A10, Mat *A11)
436: {
437:   Mat_SchurComplement *Na = (Mat_SchurComplement *)S->data;
438:   PetscBool            flg;

440:   PetscFunctionBegin;
442:   PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &flg));
443:   PetscCheck(flg, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)S)->type_name);
444:   if (A00) *A00 = Na->A;
445:   if (Ap00) *Ap00 = Na->Ap;
446:   if (A01) *A01 = Na->B;
447:   if (A10) *A10 = Na->C;
448:   if (A11) *A11 = Na->D;
449:   PetscFunctionReturn(PETSC_SUCCESS);
450: }

452: #include <petsc/private/kspimpl.h>

454: /*@
455:   MatSchurComplementComputeExplicitOperator - Compute the Schur complement matrix explicitly

457:   Collective

459:   Input Parameter:
460: . M - the matrix obtained with `MatCreateSchurComplement()`

462:   Output Parameter:
463: . S - the Schur complement matrix

465:   Notes:
466:     This can be expensive, so it is mainly for testing

468:     Use `MatSchurComplementGetPmat()` to get a sparse approximation for the Schur complement suitable for use in building a preconditioner

470:   Level: advanced

472: .seealso: [](ch_ksp), `MatCreateSchurComplement()`, `MatSchurComplementUpdate()`, `MatSchurComplementGetPmat()`
473: @*/
474: PetscErrorCode MatSchurComplementComputeExplicitOperator(Mat A, Mat *S)
475: {
476:   Mat       B, C, D, E = NULL, Bd, AinvBd;
477:   KSP       ksp;
478:   PetscInt  n, N, m, M;
479:   PetscBool flg = PETSC_FALSE, set, symm;

481:   PetscFunctionBegin;
482:   PetscCall(MatSchurComplementGetSubMatrices(A, NULL, NULL, &B, &C, &D));
483:   PetscCall(MatSchurComplementGetKSP(A, &ksp));
484:   PetscCall(KSPSetUp(ksp));
485:   PetscCall(MatConvert(B, MATDENSE, MAT_INITIAL_MATRIX, &Bd));
486:   PetscCall(MatDuplicate(Bd, MAT_DO_NOT_COPY_VALUES, &AinvBd));
487:   PetscCall(KSPMatSolve(ksp, Bd, AinvBd));
488:   PetscCall(MatDestroy(&Bd));
489:   PetscCall(MatChop(AinvBd, PETSC_SMALL));
490:   if (D) {
491:     PetscCall(MatGetLocalSize(D, &m, &n));
492:     PetscCall(MatGetSize(D, &M, &N));
493:     PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A), m, n, M, N, NULL, S));
494:   }
495:   PetscCall(MatMatMult(C, AinvBd, D ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX, PETSC_DEFAULT, S));
496:   PetscCall(MatDestroy(&AinvBd));
497:   if (D) {
498:     PetscCall(PetscObjectTypeCompareAny((PetscObject)D, &flg, MATSEQSBAIJ, MATMPISBAIJ, ""));
499:     if (flg) {
500:       PetscCall(MatIsSymmetricKnown(A, &set, &symm));
501:       if (!set || !symm) PetscCall(MatConvert(D, MATBAIJ, MAT_INITIAL_MATRIX, &E)); /* convert the (1,1) block to nonsymmetric storage for MatAXPY() */
502:     }
503:     PetscCall(MatAXPY(*S, -1.0, E ? E : D, DIFFERENT_NONZERO_PATTERN)); /* calls Mat[Get|Restore]RowUpperTriangular(), so only the upper triangular part is valid with symmetric storage */
504:   }
505:   PetscCall(MatConvert(*S, !E && flg ? MATSBAIJ : MATAIJ, MAT_INPLACE_MATRIX, S)); /* if A is symmetric and the (1,1) block is a MatSBAIJ, return S as a MatSBAIJ */
506:   PetscCall(MatScale(*S, -1.0));
507:   PetscCall(MatDestroy(&E));
508:   PetscFunctionReturn(PETSC_SUCCESS);
509: }

511: /* Developer Notes:
512:     This should be implemented with a MatCreate_SchurComplement() as that is the standard design for new Mat classes. */
513: PetscErrorCode MatGetSchurComplement_Basic(Mat mat, IS isrow0, IS iscol0, IS isrow1, IS iscol1, MatReuse mreuse, Mat *S, MatSchurComplementAinvType ainvtype, MatReuse preuse, Mat *Sp)
514: {
515:   Mat      A = NULL, Ap = NULL, B = NULL, C = NULL, D = NULL;
516:   MatReuse reuse;

518:   PetscFunctionBegin;
528:   if (mreuse == MAT_IGNORE_MATRIX && preuse == MAT_IGNORE_MATRIX) PetscFunctionReturn(PETSC_SUCCESS);

532:   PetscCheck(!mat->factortype, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");

534:   reuse = MAT_INITIAL_MATRIX;
535:   if (mreuse == MAT_REUSE_MATRIX) {
536:     PetscCall(MatSchurComplementGetSubMatrices(*S, &A, &Ap, &B, &C, &D));
537:     PetscCheck(A && Ap && B && C, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Attempting to reuse matrix but Schur complement matrices unset");
538:     PetscCheck(A == Ap, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Preconditioning matrix does not match operator");
539:     PetscCall(MatDestroy(&Ap)); /* get rid of extra reference */
540:     reuse = MAT_REUSE_MATRIX;
541:   }
542:   PetscCall(MatCreateSubMatrix(mat, isrow0, iscol0, reuse, &A));
543:   PetscCall(MatCreateSubMatrix(mat, isrow0, iscol1, reuse, &B));
544:   PetscCall(MatCreateSubMatrix(mat, isrow1, iscol0, reuse, &C));
545:   PetscCall(MatCreateSubMatrix(mat, isrow1, iscol1, reuse, &D));
546:   switch (mreuse) {
547:   case MAT_INITIAL_MATRIX:
548:     PetscCall(MatCreateSchurComplement(A, A, B, C, D, S));
549:     break;
550:   case MAT_REUSE_MATRIX:
551:     PetscCall(MatSchurComplementUpdateSubMatrices(*S, A, A, B, C, D));
552:     break;
553:   default:
554:     PetscCheck(mreuse == MAT_IGNORE_MATRIX, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Unrecognized value of mreuse %d", (int)mreuse);
555:   }
556:   if (preuse != MAT_IGNORE_MATRIX) PetscCall(MatCreateSchurComplementPmat(A, B, C, D, ainvtype, preuse, Sp));
557:   PetscCall(MatDestroy(&A));
558:   PetscCall(MatDestroy(&B));
559:   PetscCall(MatDestroy(&C));
560:   PetscCall(MatDestroy(&D));
561:   PetscFunctionReturn(PETSC_SUCCESS);
562: }

564: /*@
565:     MatGetSchurComplement - Obtain the Schur complement from eliminating part of the matrix in another part.

567:     Collective

569:     Input Parameters:
570: +   A      - matrix in which the complement is to be taken
571: .   isrow0 - rows to eliminate
572: .   iscol0 - columns to eliminate, (isrow0,iscol0) should be square and nonsingular
573: .   isrow1 - rows in which the Schur complement is formed
574: .   iscol1 - columns in which the Schur complement is formed
575: .   mreuse - `MAT_INITIAL_MATRIX` or `MAT_REUSE_MATRIX`, use `MAT_IGNORE_MATRIX` to put nothing in S
576: .   ainvtype - the type of approximation used for the inverse of the (0,0) block used in forming Sp:
577:                        `MAT_SCHUR_COMPLEMENT_AINV_DIAG`, `MAT_SCHUR_COMPLEMENT_AINV_LUMP`, `MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG`, or `MAT_SCHUR_COMPLEMENT_AINV_FULL`
578: -   preuse - `MAT_INITIAL_MATRIX` or `MAT_REUSE_MATRIX`, use `MAT_IGNORE_MATRIX` to put nothing in Sp

580:     Output Parameters:
581: +   S      - exact Schur complement, often of type `MATSCHURCOMPLEMENT` which is difficult to use for preconditioning
582: -   Sp     - approximate Schur complement from which a preconditioner can be built A11 - A10 inv(DIAGFORM(A00)) A01

584:     Level: advanced

586:     Notes:
587:     Since the real Schur complement is usually dense, providing a good approximation to Sp usually requires
588:     application-specific information.

590:     Sometimes users would like to provide problem-specific data in the Schur complement, usually only for special row
591:     and column index sets.  In that case, the user should call `PetscObjectComposeFunction()` on the *S matrix and pass mreuse of `MAT_REUSE_MATRIX` to set
592:     "MatGetSchurComplement_C" to their function.  If their function needs to fall back to the default implementation, it
593:     should call `MatGetSchurComplement_Basic()`.

595:     `MatCreateSchurComplement()` takes as arguments the four submatrices and returns the virtual Schur complement (what this function returns in S).

597:     `MatSchurComplementGetPmat()` takes the virtual Schur complement and returns an explicit approximate Schur complement (what this returns in Sp).

599:     In other words calling `MatCreateSchurComplement()` followed by `MatSchurComplementGetPmat()` produces the same output as this function but with slightly different
600:     inputs. The actually submatrices of the original block matrix instead of index sets to the submatrices.

602:     Developer Note:
603:     The API that includes `MatGetSchurComplement()`, `MatCreateSchurComplement()`, `MatSchurComplementGetPmat()` should be refactored to
604:     remove redundancy and be clearer and simpler.

606: .seealso: [](ch_ksp), `MatCreateSubMatrix()`, `PCFIELDSPLIT`, `MatCreateSchurComplement()`, `MatSchurComplementAinvType`
607: @*/
608: PetscErrorCode MatGetSchurComplement(Mat A, IS isrow0, IS iscol0, IS isrow1, IS iscol1, MatReuse mreuse, Mat *S, MatSchurComplementAinvType ainvtype, MatReuse preuse, Mat *Sp)
609: {
610:   PetscErrorCode (*f)(Mat, IS, IS, IS, IS, MatReuse, Mat *, MatReuse, Mat *) = NULL;

612:   PetscFunctionBegin;
624:   PetscCheck(!A->factortype, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
625:   if (mreuse == MAT_REUSE_MATRIX) { /* This is the only situation, in which we can demand that the user pass a non-NULL pointer to non-garbage in S. */
626:     PetscCall(PetscObjectQueryFunction((PetscObject)*S, "MatGetSchurComplement_C", &f));
627:   }
628:   if (f) PetscCall((*f)(A, isrow0, iscol0, isrow1, iscol1, mreuse, S, preuse, Sp));
629:   else PetscCall(MatGetSchurComplement_Basic(A, isrow0, iscol0, isrow1, iscol1, mreuse, S, ainvtype, preuse, Sp));
630:   PetscFunctionReturn(PETSC_SUCCESS);
631: }

633: /*@
634:     MatSchurComplementSetAinvType - set the type of approximation used for the inverse of the (0,0) block used in forming Sp in `MatSchurComplementGetPmat()`

636:     Not Collective

638:     Input Parameters:
639: +   S        - matrix obtained with `MatCreateSchurComplement()` (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
640: -   ainvtype - type of approximation to be used to form approximate Schur complement Sp = A11 - A10 inv(DIAGFORM(A00)) A01:
641:                       `MAT_SCHUR_COMPLEMENT_AINV_DIAG`, `MAT_SCHUR_COMPLEMENT_AINV_LUMP`, `MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG`, or `MAT_SCHUR_COMPLEMENT_AINV_FULL`

643:     Options Database Key:
644:     -mat_schur_complement_ainv_type diag | lump | blockdiag | full

646:     Level: advanced

648: .seealso: [](ch_ksp), `MatSchurComplementAinvType`, `MatCreateSchurComplement()`, `MatGetSchurComplement()`, `MatSchurComplementGetPmat()`, `MatSchurComplementGetAinvType()`
649: @*/
650: PetscErrorCode MatSchurComplementSetAinvType(Mat S, MatSchurComplementAinvType ainvtype)
651: {
652:   PetscBool            isschur;
653:   Mat_SchurComplement *schur;

655:   PetscFunctionBegin;
657:   PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
658:   if (!isschur) PetscFunctionReturn(PETSC_SUCCESS);
660:   schur = (Mat_SchurComplement *)S->data;
661:   PetscCheck(ainvtype == MAT_SCHUR_COMPLEMENT_AINV_DIAG || ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP || ainvtype == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG || ainvtype == MAT_SCHUR_COMPLEMENT_AINV_FULL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown MatSchurComplementAinvType: %d", (int)ainvtype);
662:   schur->ainvtype = ainvtype;
663:   PetscFunctionReturn(PETSC_SUCCESS);
664: }

666: /*@
667:     MatSchurComplementGetAinvType - get the type of approximation for the inverse of the (0,0) block used in forming Sp in `MatSchurComplementGetPmat()`

669:     Not Collective

671:     Input Parameter:
672: .   S      - matrix obtained with `MatCreateSchurComplement()` (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01

674:     Output Parameter:
675: .   ainvtype - type of approximation used to form approximate Schur complement Sp = A11 - A10 inv(DIAGFORM(A00)) A01:
676:                       `MAT_SCHUR_COMPLEMENT_AINV_DIAG`, `MAT_SCHUR_COMPLEMENT_AINV_LUMP`, `MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG`, or `MAT_SCHUR_COMPLEMENT_AINV_FULL`

678:     Level: advanced

680: .seealso: [](ch_ksp), `MatSchurComplementAinvType`, `MatCreateSchurComplement()`, `MatGetSchurComplement()`, `MatSchurComplementGetPmat()`, `MatSchurComplementSetAinvType()`
681: @*/
682: PetscErrorCode MatSchurComplementGetAinvType(Mat S, MatSchurComplementAinvType *ainvtype)
683: {
684:   PetscBool            isschur;
685:   Mat_SchurComplement *schur;

687:   PetscFunctionBegin;
689:   PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
690:   PetscCheck(isschur, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)S)->type_name);
691:   schur = (Mat_SchurComplement *)S->data;
692:   if (ainvtype) *ainvtype = schur->ainvtype;
693:   PetscFunctionReturn(PETSC_SUCCESS);
694: }

696: /*@
697:     MatCreateSchurComplementPmat - create a preconditioning matrix for the Schur complement by explicitly assembling the sparse matrix
698:         Sp = A11 - A10 inv(DIAGFORM(A00)) A01

700:     Collective

702:     Input Parameters:
703: +   A00      - the upper-left part of the original matrix A = [A00 A01; A10 A11]
704: .   A01      - (optional) the upper-right part of the original matrix A = [A00 A01; A10 A11]
705: .   A10      - (optional) the lower-left part of the original matrix A = [A00 A01; A10 A11]
706: .   A11      - (optional) the lower-right part of the original matrix A = [A00 A01; A10 A11]
707: .   ainvtype - type of approximation for DIAGFORM(A00) used when forming Sp = A11 - A10 inv(DIAGFORM(A00)) A01. See MatSchurComplementAinvType.
708: -   preuse   - `MAT_INITIAL_MATRIX` for a new Sp, or `MAT_REUSE_MATRIX` to reuse an existing Sp, or `MAT_IGNORE_MATRIX` to put nothing in Sp

710:     Output Parameter:
711: -   Sp    - approximate Schur complement suitable for preconditioning the true Schur complement S = A11 - A10 inv(A00) A01

713:     Level: advanced

715: .seealso: [](ch_ksp), `MatCreateSchurComplement()`, `MatGetSchurComplement()`, `MatSchurComplementGetPmat()`, `MatSchurComplementAinvType`
716: @*/
717: PetscErrorCode MatCreateSchurComplementPmat(Mat A00, Mat A01, Mat A10, Mat A11, MatSchurComplementAinvType ainvtype, MatReuse preuse, Mat *Sp)
718: {
719:   PetscInt N00;

721:   PetscFunctionBegin;
722:   /* Use an appropriate approximate inverse of A00 to form A11 - A10 inv(DIAGFORM(A00)) A01; a NULL A01, A10 or A11 indicates a zero matrix. */
723:   /* TODO: Perhaps should create an appropriately-sized zero matrix of the same type as A00? */
725:   if (preuse == MAT_IGNORE_MATRIX) PetscFunctionReturn(PETSC_SUCCESS);

727:   /* A zero size A00 or empty A01 or A10 imply S = A11. */
728:   PetscCall(MatGetSize(A00, &N00, NULL));
729:   if (!A01 || !A10 || !N00) {
730:     if (preuse == MAT_INITIAL_MATRIX) {
731:       PetscCall(MatDuplicate(A11, MAT_COPY_VALUES, Sp));
732:     } else { /* MAT_REUSE_MATRIX */
733:       /* TODO: when can we pass SAME_NONZERO_PATTERN? */
734:       PetscCall(MatCopy(A11, *Sp, DIFFERENT_NONZERO_PATTERN));
735:     }
736:   } else {
737:     Mat AdB;
738:     Vec diag;

740:     if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP || ainvtype == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
741:       PetscCall(MatDuplicate(A01, MAT_COPY_VALUES, &AdB));
742:       PetscCall(MatCreateVecs(A00, &diag, NULL));
743:       if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP) {
744:         PetscCall(MatGetRowSum(A00, diag));
745:       } else {
746:         PetscCall(MatGetDiagonal(A00, diag));
747:       }
748:       PetscCall(VecReciprocal(diag));
749:       PetscCall(MatDiagonalScale(AdB, diag, NULL));
750:       PetscCall(VecDestroy(&diag));
751:     } else if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG) {
752:       Mat      A00_inv;
753:       MatType  type;
754:       MPI_Comm comm;

756:       PetscCall(PetscObjectGetComm((PetscObject)A00, &comm));
757:       PetscCall(MatGetType(A00, &type));
758:       PetscCall(MatCreate(comm, &A00_inv));
759:       PetscCall(MatSetType(A00_inv, type));
760:       PetscCall(MatInvertBlockDiagonalMat(A00, A00_inv));
761:       PetscCall(MatMatMult(A00_inv, A01, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &AdB));
762:       PetscCall(MatDestroy(&A00_inv));
763:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown MatSchurComplementAinvType: %d", ainvtype);
764:     /* Cannot really reuse Sp in MatMatMult() because of MatAYPX() -->
765:          MatAXPY() --> MatHeaderReplace() --> MatDestroy_XXX_MatMatMult()  */
766:     if (preuse == MAT_REUSE_MATRIX) PetscCall(MatDestroy(Sp));
767:     PetscCall(MatMatMult(A10, AdB, MAT_INITIAL_MATRIX, PETSC_DEFAULT, Sp));
768:     if (!A11) {
769:       PetscCall(MatScale(*Sp, -1.0));
770:     } else {
771:       /* TODO: when can we pass SAME_NONZERO_PATTERN? */
772:       PetscCall(MatAYPX(*Sp, -1, A11, DIFFERENT_NONZERO_PATTERN));
773:     }
774:     PetscCall(MatDestroy(&AdB));
775:   }
776:   PetscFunctionReturn(PETSC_SUCCESS);
777: }

779: PetscErrorCode MatSchurComplementGetPmat_Basic(Mat S, MatReuse preuse, Mat *Sp)
780: {
781:   Mat                  A, B, C, D;
782:   Mat_SchurComplement *schur = (Mat_SchurComplement *)S->data;

784:   PetscFunctionBegin;
785:   if (preuse == MAT_IGNORE_MATRIX) PetscFunctionReturn(PETSC_SUCCESS);
786:   PetscCall(MatSchurComplementGetSubMatrices(S, &A, NULL, &B, &C, &D));
787:   PetscCheck(A, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONGSTATE, "Schur complement component matrices unset");
788:   if (schur->ainvtype != MAT_SCHUR_COMPLEMENT_AINV_FULL) PetscCall(MatCreateSchurComplementPmat(A, B, C, D, schur->ainvtype, preuse, Sp));
789:   else {
790:     if (preuse == MAT_REUSE_MATRIX) PetscCall(MatDestroy(Sp));
791:     PetscCall(MatSchurComplementComputeExplicitOperator(S, Sp));
792:   }
793:   PetscFunctionReturn(PETSC_SUCCESS);
794: }

796: /*@
797:     MatSchurComplementGetPmat - Obtain a preconditioning matrix for the Schur complement by assembling Sp = A11 - A10 inv(DIAGFORM(A00)) A01

799:     Collective

801:     Input Parameters:
802: +   S      - matrix obtained with MatCreateSchurComplement() (or equivalent) that implements the action of A11 - A10 ksp(A00,Ap00) A01
803: -   preuse - `MAT_INITIAL_MATRIX` for a new Sp, or `MAT_REUSE_MATRIX` to reuse an existing Sp, or `MAT_IGNORE_MATRIX` to put nothing in Sp

805:     Output Parameter:
806: -   Sp     - approximate Schur complement suitable for preconditioning the exact Schur complement S = A11 - A10 inv(A00) A01

808:     Level: advanced

810:     Notes:
811:     The approximation of Sp depends on the the argument passed to to `MatSchurComplementSetAinvType()`
812:     `MAT_SCHUR_COMPLEMENT_AINV_DIAG`, `MAT_SCHUR_COMPLEMENT_AINV_LUMP`, `MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG`, or `MAT_SCHUR_COMPLEMENT_AINV_FULL`
813:     -mat_schur_complement_ainv_type <diag,lump,blockdiag,full>

815:     Sometimes users would like to provide problem-specific data in the Schur complement, usually only
816:     for special row and column index sets.  In that case, the user should call `PetscObjectComposeFunction()` to set
817:     "MatSchurComplementGetPmat_C" to their function.  If their function needs to fall back to the default implementation,
818:     it should call `MatSchurComplementGetPmat_Basic()`.

820:     Developer Note:
821:     The API that includes `MatGetSchurComplement()`, `MatCreateSchurComplement()`, `MatSchurComplementGetPmat()` should be refactored to
822:     remove redundancy and be clearer and simpler.

824:     This routine should be called `MatSchurComplementCreatePmat()`

826: .seealso: [](ch_ksp), `MatCreateSubMatrix()`, `PCFIELDSPLIT`, `MatGetSchurComplement()`, `MatCreateSchurComplement()`, `MatSchurComplementSetAinvType()`
827: @*/
828: PetscErrorCode MatSchurComplementGetPmat(Mat S, MatReuse preuse, Mat *Sp)
829: {
830:   PetscErrorCode (*f)(Mat, MatReuse, Mat *);

832:   PetscFunctionBegin;
836:   if (preuse != MAT_IGNORE_MATRIX) {
838:     if (preuse == MAT_INITIAL_MATRIX) *Sp = NULL;
840:   }
841:   PetscCheck(!S->factortype, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");

843:   PetscCall(PetscObjectQueryFunction((PetscObject)S, "MatSchurComplementGetPmat_C", &f));
844:   if (f) PetscCall((*f)(S, preuse, Sp));
845:   else PetscCall(MatSchurComplementGetPmat_Basic(S, preuse, Sp));
846:   PetscFunctionReturn(PETSC_SUCCESS);
847: }

849: static PetscErrorCode MatProductNumeric_SchurComplement_Dense(Mat C)
850: {
851:   Mat_Product         *product = C->product;
852:   Mat_SchurComplement *Na      = (Mat_SchurComplement *)product->A->data;
853:   Mat                  work1, work2;
854:   PetscScalar         *v;
855:   PetscInt             lda;

857:   PetscFunctionBegin;
858:   PetscCall(MatMatMult(Na->B, product->B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &work1));
859:   PetscCall(MatDuplicate(work1, MAT_DO_NOT_COPY_VALUES, &work2));
860:   PetscCall(KSPMatSolve(Na->ksp, work1, work2));
861:   PetscCall(MatDestroy(&work1));
862:   PetscCall(MatDenseGetArrayWrite(C, &v));
863:   PetscCall(MatDenseGetLDA(C, &lda));
864:   PetscCall(MatCreateDense(PetscObjectComm((PetscObject)C), C->rmap->n, C->cmap->n, C->rmap->N, C->cmap->N, v, &work1));
865:   PetscCall(MatDenseSetLDA(work1, lda));
866:   PetscCall(MatMatMult(Na->C, work2, MAT_REUSE_MATRIX, PETSC_DEFAULT, &work1));
867:   PetscCall(MatDenseRestoreArrayWrite(C, &v));
868:   PetscCall(MatDestroy(&work2));
869:   PetscCall(MatDestroy(&work1));
870:   if (Na->D) {
871:     PetscCall(MatMatMult(Na->D, product->B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &work1));
872:     PetscCall(MatAYPX(C, -1.0, work1, SAME_NONZERO_PATTERN));
873:     PetscCall(MatDestroy(&work1));
874:   } else PetscCall(MatScale(C, -1.0));
875:   PetscFunctionReturn(PETSC_SUCCESS);
876: }

878: static PetscErrorCode MatProductSymbolic_SchurComplement_Dense(Mat C)
879: {
880:   Mat_Product *product = C->product;
881:   Mat          A = product->A, B = product->B;
882:   PetscInt     m = A->rmap->n, n = B->cmap->n, M = A->rmap->N, N = B->cmap->N;
883:   PetscBool    flg;

885:   PetscFunctionBegin;
886:   PetscCall(MatSetSizes(C, m, n, M, N));
887:   PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)C, &flg, MATSEQDENSE, MATMPIDENSE, ""));
888:   if (!flg) {
889:     PetscCall(MatSetType(C, ((PetscObject)B)->type_name));
890:     C->ops->productsymbolic = MatProductSymbolic_SchurComplement_Dense;
891:   }
892:   PetscCall(MatSetUp(C));
893:   C->ops->productnumeric = MatProductNumeric_SchurComplement_Dense;
894:   PetscFunctionReturn(PETSC_SUCCESS);
895: }

897: static PetscErrorCode MatProductSetFromOptions_Dense_AB(Mat C)
898: {
899:   PetscFunctionBegin;
900:   C->ops->productsymbolic = MatProductSymbolic_SchurComplement_Dense;
901:   PetscFunctionReturn(PETSC_SUCCESS);
902: }

904: static PetscErrorCode MatProductSetFromOptions_SchurComplement_Dense(Mat C)
905: {
906:   Mat_Product *product = C->product;

908:   PetscFunctionBegin;
909:   PetscCheck(product->type == MATPRODUCT_AB, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Not for product type %s", MatProductTypes[product->type]);
910:   PetscCall(MatProductSetFromOptions_Dense_AB(C));
911:   PetscFunctionReturn(PETSC_SUCCESS);
912: }

914: /*MC
915:   MATSCHURCOMPLEMENT -  "schurcomplement" - Matrix type that behaves like the Schur complement of a matrix.

917:   Level: intermediate

919: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatType`, `MatCreateSchurComplement()`, `MatSchurComplementComputeExplicitOperator()`,
920:           `MatSchurComplementGetSubMatrices()`, `MatSchurComplementGetKSP()`
921: M*/
922: PETSC_EXTERN PetscErrorCode MatCreate_SchurComplement(Mat N)
923: {
924:   Mat_SchurComplement *Na;

926:   PetscFunctionBegin;
927:   PetscCall(PetscNew(&Na));
928:   N->data = (void *)Na;

930:   N->ops->destroy        = MatDestroy_SchurComplement;
931:   N->ops->getvecs        = MatCreateVecs_SchurComplement;
932:   N->ops->view           = MatView_SchurComplement;
933:   N->ops->mult           = MatMult_SchurComplement;
934:   N->ops->multtranspose  = MatMultTranspose_SchurComplement;
935:   N->ops->multadd        = MatMultAdd_SchurComplement;
936:   N->ops->setfromoptions = MatSetFromOptions_SchurComplement;
937:   N->assembled           = PETSC_FALSE;
938:   N->preallocated        = PETSC_FALSE;

940:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)N), &Na->ksp));
941:   PetscCall(PetscObjectChangeTypeName((PetscObject)N, MATSCHURCOMPLEMENT));
942:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_schurcomplement_seqdense_C", MatProductSetFromOptions_SchurComplement_Dense));
943:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_schurcomplement_mpidense_C", MatProductSetFromOptions_SchurComplement_Dense));
944:   PetscFunctionReturn(PETSC_SUCCESS);
945: }