Actual source code: schurm.c
1: #include <../src/ksp/ksp/utils/schurm/schurm.h>
2: #include <../src/mat/impls/shell/shell.h>
4: const char *const MatSchurComplementAinvTypes[] = {"DIAG", "LUMP", "BLOCKDIAG", "FULL", "MatSchurComplementAinvType", "MAT_SCHUR_COMPLEMENT_AINV_", NULL};
6: PetscErrorCode MatCreateVecs_SchurComplement(Mat N, Vec *right, Vec *left)
7: {
8: Mat_SchurComplement *Na = (Mat_SchurComplement *)N->data;
10: PetscFunctionBegin;
11: if (Na->D) {
12: PetscCall(MatCreateVecs(Na->D, right, left));
13: PetscFunctionReturn(PETSC_SUCCESS);
14: }
15: if (right) PetscCall(MatCreateVecs(Na->B, right, NULL));
16: if (left) PetscCall(MatCreateVecs(Na->C, NULL, left));
17: PetscFunctionReturn(PETSC_SUCCESS);
18: }
20: PetscErrorCode MatView_SchurComplement(Mat N, PetscViewer viewer)
21: {
22: Mat_SchurComplement *Na = (Mat_SchurComplement *)N->data;
24: PetscFunctionBegin;
25: PetscCall(PetscViewerASCIIPrintf(viewer, "Schur complement A11 - A10 inv(A00) A01\n"));
26: if (Na->D) {
27: PetscCall(PetscViewerASCIIPrintf(viewer, "A11\n"));
28: PetscCall(PetscViewerASCIIPushTab(viewer));
29: PetscCall(MatView(Na->D, viewer));
30: PetscCall(PetscViewerASCIIPopTab(viewer));
31: } else {
32: PetscCall(PetscViewerASCIIPrintf(viewer, "A11 = 0\n"));
33: }
34: PetscCall(PetscViewerASCIIPrintf(viewer, "A10\n"));
35: PetscCall(PetscViewerASCIIPushTab(viewer));
36: PetscCall(MatView(Na->C, viewer));
37: PetscCall(PetscViewerASCIIPopTab(viewer));
38: PetscCall(PetscViewerASCIIPrintf(viewer, "KSP solver for A00 block viewable with the additional option -%sksp_view\n", ((PetscObject)Na->ksp)->prefix ? ((PetscObject)Na->ksp)->prefix : NULL));
39: PetscCall(PetscViewerASCIIPrintf(viewer, "A01\n"));
40: PetscCall(PetscViewerASCIIPushTab(viewer));
41: PetscCall(MatView(Na->B, viewer));
42: PetscCall(PetscViewerASCIIPopTab(viewer));
43: PetscFunctionReturn(PETSC_SUCCESS);
44: }
46: /*
47: A11^T - A01^T ksptrans(A00,Ap00) A10^T
48: */
49: PetscErrorCode MatMultTranspose_SchurComplement(Mat N, Vec x, Vec y)
50: {
51: Mat_SchurComplement *Na = (Mat_SchurComplement *)N->data;
53: PetscFunctionBegin;
54: if (!Na->work1) PetscCall(MatCreateVecs(Na->A, &Na->work1, NULL));
55: if (!Na->work2) PetscCall(MatCreateVecs(Na->A, &Na->work2, NULL));
56: PetscCall(MatMultTranspose(Na->C, x, Na->work1));
57: PetscCall(KSPSolveTranspose(Na->ksp, Na->work1, Na->work2));
58: PetscCall(MatMultTranspose(Na->B, Na->work2, y));
59: PetscCall(VecScale(y, -1.0));
60: if (Na->D) PetscCall(MatMultTransposeAdd(Na->D, x, y, y));
61: PetscFunctionReturn(PETSC_SUCCESS);
62: }
64: /*
65: A11 - A10 ksp(A00,Ap00) A01
66: */
67: PetscErrorCode MatMult_SchurComplement(Mat N, Vec x, Vec y)
68: {
69: Mat_SchurComplement *Na = (Mat_SchurComplement *)N->data;
71: PetscFunctionBegin;
72: if (!Na->work1) PetscCall(MatCreateVecs(Na->A, &Na->work1, NULL));
73: if (!Na->work2) PetscCall(MatCreateVecs(Na->A, &Na->work2, NULL));
74: PetscCall(MatMult(Na->B, x, Na->work1));
75: PetscCall(KSPSolve(Na->ksp, Na->work1, Na->work2));
76: PetscCall(MatMult(Na->C, Na->work2, y));
77: PetscCall(VecScale(y, -1.0));
78: if (Na->D) PetscCall(MatMultAdd(Na->D, x, y, y));
79: PetscFunctionReturn(PETSC_SUCCESS);
80: }
82: /*
83: A11 - A10 ksp(A00,Ap00) A01
84: */
85: PetscErrorCode MatMultAdd_SchurComplement(Mat N, Vec x, Vec y, Vec z)
86: {
87: Mat_SchurComplement *Na = (Mat_SchurComplement *)N->data;
89: PetscFunctionBegin;
90: if (!Na->work1) PetscCall(MatCreateVecs(Na->A, &Na->work1, NULL));
91: if (!Na->work2) PetscCall(MatCreateVecs(Na->A, &Na->work2, NULL));
92: PetscCall(MatMult(Na->B, x, Na->work1));
93: PetscCall(KSPSolve(Na->ksp, Na->work1, Na->work2));
94: if (y == z) {
95: PetscCall(VecScale(Na->work2, -1.0));
96: PetscCall(MatMultAdd(Na->C, Na->work2, z, z));
97: } else {
98: PetscCall(MatMult(Na->C, Na->work2, z));
99: PetscCall(VecAYPX(z, -1.0, y));
100: }
101: if (Na->D) PetscCall(MatMultAdd(Na->D, x, z, z));
102: PetscFunctionReturn(PETSC_SUCCESS);
103: }
105: PetscErrorCode MatSetFromOptions_SchurComplement(Mat N, PetscOptionItems *PetscOptionsObject)
106: {
107: Mat_SchurComplement *Na = (Mat_SchurComplement *)N->data;
109: PetscFunctionBegin;
110: PetscOptionsHeadBegin(PetscOptionsObject, "MatSchurComplementOptions");
111: Na->ainvtype = MAT_SCHUR_COMPLEMENT_AINV_DIAG;
112: 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,
113: (PetscEnum *)&Na->ainvtype, NULL));
114: PetscOptionsHeadEnd();
115: PetscCall(KSPSetFromOptions(Na->ksp));
116: PetscFunctionReturn(PETSC_SUCCESS);
117: }
119: PetscErrorCode MatDestroy_SchurComplement(Mat N)
120: {
121: Mat_SchurComplement *Na = (Mat_SchurComplement *)N->data;
123: PetscFunctionBegin;
124: PetscCall(MatDestroy(&Na->A));
125: PetscCall(MatDestroy(&Na->Ap));
126: PetscCall(MatDestroy(&Na->B));
127: PetscCall(MatDestroy(&Na->C));
128: PetscCall(MatDestroy(&Na->D));
129: PetscCall(VecDestroy(&Na->work1));
130: PetscCall(VecDestroy(&Na->work2));
131: PetscCall(KSPDestroy(&Na->ksp));
132: PetscCall(PetscFree(N->data));
133: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_schurcomplement_seqdense_C", NULL));
134: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_schurcomplement_mpidense_C", NULL));
135: PetscFunctionReturn(PETSC_SUCCESS);
136: }
138: /*@
139: MatCreateSchurComplement - Creates a new `Mat` that behaves like the Schur complement of a matrix
141: Collective
143: Input Parameters:
144: + A00 - the upper-left block of the original matrix A = [A00 A01; A10 A11]
145: . Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A00^{-1}
146: . A01 - the upper-right block of the original matrix A = [A00 A01; A10 A11]
147: . A10 - the lower-left block of the original matrix A = [A00 A01; A10 A11]
148: - A11 - (optional) the lower-right block of the original matrix A = [A00 A01; A10 A11]
150: Output Parameter:
151: . S - the matrix that behaves as the Schur complement S = A11 - A10 ksp(A00,Ap00) A01
153: Level: intermediate
155: Notes:
156: The Schur complement is NOT explicitly formed! Rather, this function returns a virtual Schur complement
157: that can compute the matrix-vector product by using formula S = A11 - A10 A^{-1} A01
158: for Schur complement S and a `KSP` solver to approximate the action of A^{-1}.
160: All four matrices must have the same MPI communicator.
162: `A00` and `A11` must be square matrices.
164: `MatGetSchurComplement()` takes as arguments the index sets for the submatrices and returns both the virtual Schur complement (what this returns) plus
165: a sparse approximation to the Schur complement (useful for building a preconditioner for the Schur complement) which can be obtained from this
166: matrix with `MatSchurComplementGetPmat()`
168: Developer Notes:
169: The API that includes `MatGetSchurComplement()`, `MatCreateSchurComplement()`, `MatSchurComplementGetPmat()` should be refactored to
170: remove redundancy and be clearer and simpler.
172: .seealso: [](ch_ksp), `MatCreateNormal()`, `MatMult()`, `MatCreate()`, `MatSchurComplementGetKSP()`, `MatSchurComplementUpdateSubMatrices()`, `MatCreateTranspose()`, `MatGetSchurComplement()`,
173: `MatSchurComplementGetPmat()`, `MatSchurComplementSetSubMatrices()`
174: @*/
175: PetscErrorCode MatCreateSchurComplement(Mat A00, Mat Ap00, Mat A01, Mat A10, Mat A11, Mat *S)
176: {
177: PetscFunctionBegin;
178: PetscCall(KSPInitializePackage());
179: PetscCall(MatCreate(PetscObjectComm((PetscObject)A00), S));
180: PetscCall(MatSetType(*S, MATSCHURCOMPLEMENT));
181: PetscCall(MatSchurComplementSetSubMatrices(*S, A00, Ap00, A01, A10, A11));
182: PetscFunctionReturn(PETSC_SUCCESS);
183: }
185: /*@
186: MatSchurComplementSetSubMatrices - Sets the matrices that define the Schur complement
188: Collective
190: Input Parameters:
191: + S - matrix obtained with `MatSetType`(S,`MATSCHURCOMPLEMENT`)
192: . A00 - the upper-left block of the original matrix A = [A00 A01; A10 A11]
193: . Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A00^{-1}
194: . A01 - the upper-right block of the original matrix A = [A00 A01; A10 A11]
195: . A10 - the lower-left block of the original matrix A = [A00 A01; A10 A11]
196: - A11 - (optional) the lower-right block of the original matrix A = [A00 A01; A10 A11]
198: Level: intermediate
200: Notes:
201: The Schur complement is NOT explicitly formed! Rather, this
202: object performs the matrix-vector product of the Schur complement by using formula S = A11 - A10 ksp(A00,Ap00) A01
204: All four matrices must have the same MPI communicator.
206: `A00` and `A11` must be square matrices.
208: This is to be used in the context of code such as
209: .vb
210: MatSetType(S,MATSCHURCOMPLEMENT);
211: MatSchurComplementSetSubMatrices(S,...);
212: .ve
213: while `MatSchurComplementUpdateSubMatrices()` should only be called after `MatCreateSchurComplement()` or `MatSchurComplementSetSubMatrices()`
215: .seealso: [](ch_ksp), `Mat`, `MatCreateNormal()`, `MatMult()`, `MatCreate()`, `MatSchurComplementGetKSP()`, `MatSchurComplementUpdateSubMatrices()`, `MatCreateTranspose()`, `MatCreateSchurComplement()`, `MatGetSchurComplement()`
216: @*/
217: PetscErrorCode MatSchurComplementSetSubMatrices(Mat S, Mat A00, Mat Ap00, Mat A01, Mat A10, Mat A11)
218: {
219: Mat_SchurComplement *Na = (Mat_SchurComplement *)S->data;
220: PetscBool isschur;
222: PetscFunctionBegin;
223: PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
224: if (!isschur) PetscFunctionReturn(PETSC_SUCCESS);
225: PetscCheck(!S->assembled, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONGSTATE, "Use MatSchurComplementUpdateSubMatrices() for already used matrix");
230: PetscCheckSameComm(A00, 2, Ap00, 3);
231: PetscCheckSameComm(A00, 2, A01, 4);
232: PetscCheckSameComm(A00, 2, A10, 5);
233: 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);
234: 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);
235: 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);
236: 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);
237: 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);
238: if (A11) {
240: PetscCheckSameComm(A00, 2, A11, 6);
241: 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);
242: }
244: PetscCall(MatSetSizes(S, A10->rmap->n, A01->cmap->n, A10->rmap->N, A01->cmap->N));
245: PetscCall(PetscObjectReference((PetscObject)A00));
246: PetscCall(PetscObjectReference((PetscObject)Ap00));
247: PetscCall(PetscObjectReference((PetscObject)A01));
248: PetscCall(PetscObjectReference((PetscObject)A10));
249: Na->A = A00;
250: Na->Ap = Ap00;
251: Na->B = A01;
252: Na->C = A10;
253: Na->D = A11;
254: if (A11) PetscCall(PetscObjectReference((PetscObject)A11));
255: PetscCall(MatSetUp(S));
256: PetscCall(KSPSetOperators(Na->ksp, A00, Ap00));
257: S->assembled = PETSC_TRUE;
258: PetscFunctionReturn(PETSC_SUCCESS);
259: }
261: /*@
262: MatSchurComplementGetKSP - Gets the `KSP` object that is used to solve with A00 in the Schur complement matrix S = A11 - A10 ksp(A00,Ap00) A01
264: Not Collective
266: Input Parameter:
267: . S - matrix obtained with `MatCreateSchurComplement()` (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
269: Output Parameter:
270: . ksp - the linear solver object
272: Options Database Key:
273: . -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.
275: Level: intermediate
277: .seealso: [](ch_ksp), `Mat`, `MatSchurComplementSetKSP()`, `MatCreateSchurComplement()`, `MatCreateNormal()`, `MatMult()`, `MatCreate()`
278: @*/
279: PetscErrorCode MatSchurComplementGetKSP(Mat S, KSP *ksp)
280: {
281: Mat_SchurComplement *Na;
282: PetscBool isschur;
284: PetscFunctionBegin;
286: PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
287: PetscCheck(isschur, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)S)->type_name);
288: PetscAssertPointer(ksp, 2);
289: Na = (Mat_SchurComplement *)S->data;
290: *ksp = Na->ksp;
291: PetscFunctionReturn(PETSC_SUCCESS);
292: }
294: /*@
295: MatSchurComplementSetKSP - Sets the `KSP` object that is used to solve with A00 in the Schur complement matrix S = A11 - A10 ksp(A00,Ap00) A01
297: Not Collective
299: Input Parameters:
300: + S - matrix created with `MatCreateSchurComplement()`
301: - ksp - the linear solver object
303: Level: developer
305: Developer Notes:
306: This is used in `PCFIELDSPLIT` to reuse the 0-split `KSP` to implement ksp(A00,Ap00) in S.
307: The `KSP` operators are overwritten with A00 and Ap00 currently set in S.
309: .seealso: [](ch_ksp), `Mat`, `MatSchurComplementGetKSP()`, `MatCreateSchurComplement()`, `MatCreateNormal()`, `MatMult()`, `MatCreate()`, `MATSCHURCOMPLEMENT`
310: @*/
311: PetscErrorCode MatSchurComplementSetKSP(Mat S, KSP ksp)
312: {
313: Mat_SchurComplement *Na;
314: PetscBool isschur;
316: PetscFunctionBegin;
318: PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
319: if (!isschur) PetscFunctionReturn(PETSC_SUCCESS);
321: Na = (Mat_SchurComplement *)S->data;
322: PetscCall(PetscObjectReference((PetscObject)ksp));
323: PetscCall(KSPDestroy(&Na->ksp));
324: Na->ksp = ksp;
325: PetscCall(KSPSetOperators(Na->ksp, Na->A, Na->Ap));
326: PetscFunctionReturn(PETSC_SUCCESS);
327: }
329: /*@
330: MatSchurComplementUpdateSubMatrices - Updates the Schur complement matrix object with new submatrices
332: Collective
334: Input Parameters:
335: + S - matrix obtained with `MatCreateSchurComplement()` (or `MatSchurSetSubMatrices()`) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
336: . A00 - the upper-left block of the original matrix A = [A00 A01; A10 A11]
337: . Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A00^{-1}
338: . A01 - the upper-right block of the original matrix A = [A00 A01; A10 A11]
339: . A10 - the lower-left block of the original matrix A = [A00 A01; A10 A11]
340: - A11 - (optional) the lower-right block of the original matrix A = [A00 A01; A10 A11]
342: Level: intermediate
344: Notes:
345: All four matrices must have the same MPI communicator
347: `A00` and `A11` must be square matrices
349: All of the matrices provided must have the same sizes as was used with `MatCreateSchurComplement()` or `MatSchurComplementSetSubMatrices()`
350: though they need not be the same matrices.
352: This can only be called after `MatCreateSchurComplement()` or `MatSchurComplementSetSubMatrices()`, it cannot be called immediately after `MatSetType`(S,`MATSCHURCOMPLEMENT`);
354: Developer Notes:
355: This code is almost identical to `MatSchurComplementSetSubMatrices()`. The API should be refactored.
357: .seealso: [](ch_ksp), `Mat`, `MatCreateNormal()`, `MatMult()`, `MatCreate()`, `MatSchurComplementGetKSP()`, `MatCreateSchurComplement()`
358: @*/
359: PetscErrorCode MatSchurComplementUpdateSubMatrices(Mat S, Mat A00, Mat Ap00, Mat A01, Mat A10, Mat A11)
360: {
361: Mat_SchurComplement *Na = (Mat_SchurComplement *)S->data;
362: PetscBool isschur;
364: PetscFunctionBegin;
366: PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
367: if (!isschur) PetscFunctionReturn(PETSC_SUCCESS);
368: PetscCheck(S->assembled, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONGSTATE, "Use MatSchurComplementSetSubMatrices() for a new matrix");
373: PetscCheckSameComm(A00, 2, Ap00, 3);
374: PetscCheckSameComm(A00, 2, A01, 4);
375: PetscCheckSameComm(A00, 2, A10, 5);
376: 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);
377: 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);
378: 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);
379: 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);
380: 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);
381: if (A11) {
383: PetscCheckSameComm(A00, 2, A11, 6);
384: 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);
385: }
387: PetscCall(PetscObjectReference((PetscObject)A00));
388: PetscCall(PetscObjectReference((PetscObject)Ap00));
389: PetscCall(PetscObjectReference((PetscObject)A01));
390: PetscCall(PetscObjectReference((PetscObject)A10));
391: if (A11) PetscCall(PetscObjectReference((PetscObject)A11));
393: PetscCall(MatDestroy(&Na->A));
394: PetscCall(MatDestroy(&Na->Ap));
395: PetscCall(MatDestroy(&Na->B));
396: PetscCall(MatDestroy(&Na->C));
397: PetscCall(MatDestroy(&Na->D));
399: Na->A = A00;
400: Na->Ap = Ap00;
401: Na->B = A01;
402: Na->C = A10;
403: Na->D = A11;
405: PetscCall(KSPSetOperators(Na->ksp, A00, Ap00));
406: PetscFunctionReturn(PETSC_SUCCESS);
407: }
409: /*@C
410: MatSchurComplementGetSubMatrices - Get the individual submatrices in the Schur complement
412: Collective
414: Input Parameter:
415: . S - matrix obtained with `MatCreateSchurComplement()` (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
417: Output Parameters:
418: + A00 - the upper-left block of the original matrix A = [A00 A01; A10 A11]
419: . Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}
420: . A01 - the upper-right block of the original matrix A = [A00 A01; A10 A11]
421: . A10 - the lower-left block of the original matrix A = [A00 A01; A10 A11]
422: - A11 - (optional) the lower-right block of the original matrix A = [A00 A01; A10 A11]
424: Level: intermediate
426: Note:
427: Use `NULL` for any unneeded output argument.
429: The reference counts of the submatrices are not increased before they are returned and the matrices should not be modified or destroyed.
431: .seealso: [](ch_ksp), `MatCreateNormal()`, `MatMult()`, `MatCreate()`, `MatSchurComplementGetKSP()`, `MatCreateSchurComplement()`, `MatSchurComplementUpdateSubMatrices()`
432: @*/
433: PetscErrorCode MatSchurComplementGetSubMatrices(Mat S, Mat *A00, Mat *Ap00, Mat *A01, Mat *A10, Mat *A11)
434: {
435: Mat_SchurComplement *Na = (Mat_SchurComplement *)S->data;
436: PetscBool flg;
438: PetscFunctionBegin;
440: PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &flg));
441: PetscCheck(flg, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)S)->type_name);
442: if (A00) *A00 = Na->A;
443: if (Ap00) *Ap00 = Na->Ap;
444: if (A01) *A01 = Na->B;
445: if (A10) *A10 = Na->C;
446: if (A11) *A11 = Na->D;
447: PetscFunctionReturn(PETSC_SUCCESS);
448: }
450: #include <petsc/private/kspimpl.h>
452: /*@
453: MatSchurComplementComputeExplicitOperator - Compute the Schur complement matrix explicitly
455: Collective
457: Input Parameter:
458: . A - the matrix obtained with `MatCreateSchurComplement()`
460: Output Parameter:
461: . S - the Schur complement matrix
463: Notes:
464: This can be expensive, so it is mainly for testing
466: Use `MatSchurComplementGetPmat()` to get a sparse approximation for the Schur complement suitable for use in building a preconditioner
468: Level: advanced
470: .seealso: [](ch_ksp), `MatCreateSchurComplement()`, `MatSchurComplementUpdateSubMatrices()`, `MatSchurComplementGetPmat()`
471: @*/
472: PetscErrorCode MatSchurComplementComputeExplicitOperator(Mat A, Mat *S)
473: {
474: Mat B, C, D, E = NULL, Bd, AinvBd;
475: KSP ksp;
476: PetscInt n, N, m, M;
477: PetscBool flg = PETSC_FALSE, set, symm;
479: PetscFunctionBegin;
480: PetscCall(MatSchurComplementGetSubMatrices(A, NULL, NULL, &B, &C, &D));
481: PetscCall(MatSchurComplementGetKSP(A, &ksp));
482: PetscCall(KSPSetUp(ksp));
483: PetscCall(MatConvert(B, MATDENSE, MAT_INITIAL_MATRIX, &Bd));
484: PetscCall(MatDuplicate(Bd, MAT_DO_NOT_COPY_VALUES, &AinvBd));
485: PetscCall(KSPMatSolve(ksp, Bd, AinvBd));
486: PetscCall(MatDestroy(&Bd));
487: PetscCall(MatFilter(AinvBd, PETSC_SMALL, PETSC_FALSE, PETSC_FALSE));
488: if (D) {
489: PetscCall(MatGetLocalSize(D, &m, &n));
490: PetscCall(MatGetSize(D, &M, &N));
491: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A), m, n, M, N, NULL, S));
492: }
493: PetscCall(MatMatMult(C, AinvBd, D ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX, PETSC_DEFAULT, S));
494: PetscCall(MatDestroy(&AinvBd));
495: if (D) {
496: PetscCall(PetscObjectTypeCompareAny((PetscObject)D, &flg, MATSEQSBAIJ, MATMPISBAIJ, ""));
497: if (flg) {
498: PetscCall(MatIsSymmetricKnown(A, &set, &symm));
499: if (!set || !symm) PetscCall(MatConvert(D, MATBAIJ, MAT_INITIAL_MATRIX, &E)); /* convert the (1,1) block to nonsymmetric storage for MatAXPY() */
500: }
501: 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 */
502: }
503: 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 */
504: PetscCall(MatScale(*S, -1.0));
505: PetscCall(MatDestroy(&E));
506: PetscFunctionReturn(PETSC_SUCCESS);
507: }
509: /* Developer Notes:
510: This should be implemented with a MatCreate_SchurComplement() as that is the standard design for new Mat classes. */
511: PetscErrorCode MatGetSchurComplement_Basic(Mat mat, IS isrow0, IS iscol0, IS isrow1, IS iscol1, MatReuse mreuse, Mat *S, MatSchurComplementAinvType ainvtype, MatReuse preuse, Mat *Sp)
512: {
513: Mat A = NULL, Ap = NULL, B = NULL, C = NULL, D = NULL;
514: MatReuse reuse;
516: PetscFunctionBegin;
526: if (mreuse == MAT_IGNORE_MATRIX && preuse == MAT_IGNORE_MATRIX) PetscFunctionReturn(PETSC_SUCCESS);
530: PetscCheck(!mat->factortype, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
532: reuse = MAT_INITIAL_MATRIX;
533: if (mreuse == MAT_REUSE_MATRIX) {
534: PetscCall(MatSchurComplementGetSubMatrices(*S, &A, &Ap, &B, &C, &D));
535: PetscCheck(A && Ap && B && C, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Attempting to reuse matrix but Schur complement matrices unset");
536: PetscCheck(A == Ap, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Preconditioning matrix does not match operator");
537: PetscCall(MatDestroy(&Ap)); /* get rid of extra reference */
538: reuse = MAT_REUSE_MATRIX;
539: }
540: PetscCall(MatCreateSubMatrix(mat, isrow0, iscol0, reuse, &A));
541: PetscCall(MatCreateSubMatrix(mat, isrow0, iscol1, reuse, &B));
542: PetscCall(MatCreateSubMatrix(mat, isrow1, iscol0, reuse, &C));
543: PetscCall(MatCreateSubMatrix(mat, isrow1, iscol1, reuse, &D));
544: switch (mreuse) {
545: case MAT_INITIAL_MATRIX:
546: PetscCall(MatCreateSchurComplement(A, A, B, C, D, S));
547: break;
548: case MAT_REUSE_MATRIX:
549: PetscCall(MatSchurComplementUpdateSubMatrices(*S, A, A, B, C, D));
550: break;
551: default:
552: PetscCheck(mreuse == MAT_IGNORE_MATRIX, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Unrecognized value of mreuse %d", (int)mreuse);
553: }
554: if (preuse != MAT_IGNORE_MATRIX) PetscCall(MatCreateSchurComplementPmat(A, B, C, D, ainvtype, preuse, Sp));
555: PetscCall(MatDestroy(&A));
556: PetscCall(MatDestroy(&B));
557: PetscCall(MatDestroy(&C));
558: PetscCall(MatDestroy(&D));
559: PetscFunctionReturn(PETSC_SUCCESS);
560: }
562: /*@
563: MatGetSchurComplement - Obtain the Schur complement from eliminating part of the matrix in another part.
565: Collective
567: Input Parameters:
568: + A - matrix in which the complement is to be taken
569: . isrow0 - rows to eliminate
570: . iscol0 - columns to eliminate, (isrow0,iscol0) should be square and nonsingular
571: . isrow1 - rows in which the Schur complement is formed
572: . iscol1 - columns in which the Schur complement is formed
573: . mreuse - `MAT_INITIAL_MATRIX` or `MAT_REUSE_MATRIX`, use `MAT_IGNORE_MATRIX` to put nothing in S
574: . ainvtype - the type of approximation used for the inverse of the (0,0) block used in forming Sp:
575: `MAT_SCHUR_COMPLEMENT_AINV_DIAG`, `MAT_SCHUR_COMPLEMENT_AINV_LUMP`, `MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG`, or `MAT_SCHUR_COMPLEMENT_AINV_FULL`
576: - preuse - `MAT_INITIAL_MATRIX` or `MAT_REUSE_MATRIX`, use `MAT_IGNORE_MATRIX` to put nothing in Sp
578: Output Parameters:
579: + S - exact Schur complement, often of type `MATSCHURCOMPLEMENT` which is difficult to use for preconditioning
580: - Sp - approximate Schur complement from which a preconditioner can be built A11 - A10 inv(DIAGFORM(A00)) A01
582: Level: advanced
584: Notes:
585: Since the real Schur complement is usually dense, providing a good approximation to Sp usually requires
586: application-specific information.
588: Sometimes users would like to provide problem-specific data in the Schur complement, usually only for special row
589: 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
590: "MatGetSchurComplement_C" to their function. If their function needs to fall back to the default implementation, it
591: should call `MatGetSchurComplement_Basic()`.
593: `MatCreateSchurComplement()` takes as arguments the four submatrices and returns the virtual Schur complement (what this function returns in S).
595: `MatSchurComplementGetPmat()` takes the virtual Schur complement and returns an explicit approximate Schur complement (what this returns in Sp).
597: In other words calling `MatCreateSchurComplement()` followed by `MatSchurComplementGetPmat()` produces the same output as this function but with slightly different
598: inputs. The actually submatrices of the original block matrix instead of index sets to the submatrices.
600: Developer Notes:
601: The API that includes `MatGetSchurComplement()`, `MatCreateSchurComplement()`, `MatSchurComplementGetPmat()` should be refactored to
602: remove redundancy and be clearer and simpler.
604: .seealso: [](ch_ksp), `MatCreateSubMatrix()`, `PCFIELDSPLIT`, `MatCreateSchurComplement()`, `MatSchurComplementAinvType`
605: @*/
606: PetscErrorCode MatGetSchurComplement(Mat A, IS isrow0, IS iscol0, IS isrow1, IS iscol1, MatReuse mreuse, Mat *S, MatSchurComplementAinvType ainvtype, MatReuse preuse, Mat *Sp)
607: {
608: PetscErrorCode (*f)(Mat, IS, IS, IS, IS, MatReuse, Mat *, MatReuse, Mat *) = NULL;
610: PetscFunctionBegin;
622: PetscCheck(!A->factortype, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
623: 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. */
624: PetscCall(PetscObjectQueryFunction((PetscObject)*S, "MatGetSchurComplement_C", &f));
625: }
626: if (f) PetscCall((*f)(A, isrow0, iscol0, isrow1, iscol1, mreuse, S, preuse, Sp));
627: else PetscCall(MatGetSchurComplement_Basic(A, isrow0, iscol0, isrow1, iscol1, mreuse, S, ainvtype, preuse, Sp));
628: PetscFunctionReturn(PETSC_SUCCESS);
629: }
631: /*@
632: MatSchurComplementSetAinvType - set the type of approximation used for the inverse of the (0,0) block used in forming Sp in `MatSchurComplementGetPmat()`
634: Not Collective
636: Input Parameters:
637: + S - matrix obtained with `MatCreateSchurComplement()` (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
638: - ainvtype - type of approximation to be used to form approximate Schur complement Sp = A11 - A10 inv(DIAGFORM(A00)) A01:
639: `MAT_SCHUR_COMPLEMENT_AINV_DIAG`, `MAT_SCHUR_COMPLEMENT_AINV_LUMP`, `MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG`, or `MAT_SCHUR_COMPLEMENT_AINV_FULL`
641: Options Database Key:
642: . -mat_schur_complement_ainv_type diag | lump | blockdiag | full - set schur complement type
644: Level: advanced
646: .seealso: [](ch_ksp), `MatSchurComplementAinvType`, `MatCreateSchurComplement()`, `MatGetSchurComplement()`, `MatSchurComplementGetPmat()`, `MatSchurComplementGetAinvType()`
647: @*/
648: PetscErrorCode MatSchurComplementSetAinvType(Mat S, MatSchurComplementAinvType ainvtype)
649: {
650: PetscBool isschur;
651: Mat_SchurComplement *schur;
653: PetscFunctionBegin;
655: PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
656: if (!isschur) PetscFunctionReturn(PETSC_SUCCESS);
658: schur = (Mat_SchurComplement *)S->data;
659: 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);
660: schur->ainvtype = ainvtype;
661: PetscFunctionReturn(PETSC_SUCCESS);
662: }
664: /*@
665: MatSchurComplementGetAinvType - get the type of approximation for the inverse of the (0,0) block used in forming Sp in `MatSchurComplementGetPmat()`
667: Not Collective
669: Input Parameter:
670: . S - matrix obtained with `MatCreateSchurComplement()` (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
672: Output Parameter:
673: . ainvtype - type of approximation used to form approximate Schur complement Sp = A11 - A10 inv(DIAGFORM(A00)) A01:
674: `MAT_SCHUR_COMPLEMENT_AINV_DIAG`, `MAT_SCHUR_COMPLEMENT_AINV_LUMP`, `MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG`, or `MAT_SCHUR_COMPLEMENT_AINV_FULL`
676: Level: advanced
678: .seealso: [](ch_ksp), `MatSchurComplementAinvType`, `MatCreateSchurComplement()`, `MatGetSchurComplement()`, `MatSchurComplementGetPmat()`, `MatSchurComplementSetAinvType()`
679: @*/
680: PetscErrorCode MatSchurComplementGetAinvType(Mat S, MatSchurComplementAinvType *ainvtype)
681: {
682: PetscBool isschur;
683: Mat_SchurComplement *schur;
685: PetscFunctionBegin;
687: PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
688: PetscCheck(isschur, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)S)->type_name);
689: schur = (Mat_SchurComplement *)S->data;
690: if (ainvtype) *ainvtype = schur->ainvtype;
691: PetscFunctionReturn(PETSC_SUCCESS);
692: }
694: /*@
695: MatCreateSchurComplementPmat - create a preconditioning matrix for the Schur complement by explicitly assembling the sparse matrix
696: Sp = A11 - A10 inv(DIAGFORM(A00)) A01
698: Collective
700: Input Parameters:
701: + A00 - the upper-left part of the original matrix A = [A00 A01; A10 A11]
702: . A01 - (optional) the upper-right part of the original matrix A = [A00 A01; A10 A11]
703: . A10 - (optional) the lower-left part of the original matrix A = [A00 A01; A10 A11]
704: . A11 - (optional) the lower-right part of the original matrix A = [A00 A01; A10 A11]
705: . ainvtype - type of approximation for DIAGFORM(A00) used when forming Sp = A11 - A10 inv(DIAGFORM(A00)) A01. See MatSchurComplementAinvType.
706: - 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
708: Output Parameter:
709: . Sp - approximate Schur complement suitable for preconditioning the true Schur complement S = A11 - A10 inv(A00) A01
711: Level: advanced
713: .seealso: [](ch_ksp), `MatCreateSchurComplement()`, `MatGetSchurComplement()`, `MatSchurComplementGetPmat()`, `MatSchurComplementAinvType`
714: @*/
715: PetscErrorCode MatCreateSchurComplementPmat(Mat A00, Mat A01, Mat A10, Mat A11, MatSchurComplementAinvType ainvtype, MatReuse preuse, Mat *Sp)
716: {
717: PetscInt N00;
719: PetscFunctionBegin;
720: /* 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. */
721: /* TODO: Perhaps should create an appropriately-sized zero matrix of the same type as A00? */
723: if (preuse == MAT_IGNORE_MATRIX) PetscFunctionReturn(PETSC_SUCCESS);
725: /* A zero size A00 or empty A01 or A10 imply S = A11. */
726: PetscCall(MatGetSize(A00, &N00, NULL));
727: if (!A01 || !A10 || !N00) {
728: if (preuse == MAT_INITIAL_MATRIX) {
729: PetscCall(MatDuplicate(A11, MAT_COPY_VALUES, Sp));
730: } else { /* MAT_REUSE_MATRIX */
731: /* TODO: when can we pass SAME_NONZERO_PATTERN? */
732: PetscCall(MatCopy(A11, *Sp, DIFFERENT_NONZERO_PATTERN));
733: }
734: } else {
735: Mat AdB, T;
736: Vec diag;
737: PetscBool flg;
739: if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP || ainvtype == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
740: PetscCall(PetscObjectTypeCompare((PetscObject)A01, MATTRANSPOSEVIRTUAL, &flg));
741: if (flg) {
742: PetscCall(MatTransposeGetMat(A01, &T));
743: PetscCall(MatTranspose(T, MAT_INITIAL_MATRIX, &AdB));
744: } else {
745: PetscCall(PetscObjectTypeCompare((PetscObject)A01, MATHERMITIANTRANSPOSEVIRTUAL, &flg));
746: if (flg) {
747: PetscCall(MatHermitianTransposeGetMat(A01, &T));
748: PetscCall(MatHermitianTranspose(T, MAT_INITIAL_MATRIX, &AdB));
749: }
750: }
751: if (!flg) PetscCall(MatDuplicate(A01, MAT_COPY_VALUES, &AdB));
752: else {
753: PetscCheck(!((Mat_Shell *)A01->data)->zrows && !((Mat_Shell *)A01->data)->zcols, PetscObjectComm((PetscObject)A01), PETSC_ERR_SUP, "Cannot call MatCreateSchurComplementPmat() if MatZeroRows() or MatZeroRowsColumns() has been called on the input Mat");
754: PetscCheck(!((Mat_Shell *)A01->data)->axpy, PetscObjectComm((PetscObject)A01), PETSC_ERR_SUP, "Cannot call MatCreateSchurComplementPmat() if MatAXPY() has been called on the input Mat");
755: PetscCheck(!((Mat_Shell *)A01->data)->left && !((Mat_Shell *)A01->data)->right, PetscObjectComm((PetscObject)A01), PETSC_ERR_SUP, "Cannot call MatCreateSchurComplementPmat() if MatDiagonalScale() has been called on the input Mat");
756: PetscCheck(!((Mat_Shell *)A01->data)->dshift, PetscObjectComm((PetscObject)A01), PETSC_ERR_SUP, "Cannot call MatCreateSchurComplementPmat() if MatDiagonalSet() has been called on the input Mat");
757: PetscCall(MatScale(AdB, ((Mat_Shell *)A01->data)->vscale));
758: PetscCall(MatShift(AdB, ((Mat_Shell *)A01->data)->vshift));
759: }
760: PetscCall(MatCreateVecs(A00, &diag, NULL));
761: if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP) {
762: PetscCall(MatGetRowSum(A00, diag));
763: } else {
764: PetscCall(MatGetDiagonal(A00, diag));
765: }
766: PetscCall(VecReciprocal(diag));
767: PetscCall(MatDiagonalScale(AdB, diag, NULL));
768: PetscCall(VecDestroy(&diag));
769: } else if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG) {
770: Mat A00_inv;
771: MatType type;
772: MPI_Comm comm;
774: PetscCall(PetscObjectGetComm((PetscObject)A00, &comm));
775: PetscCall(MatGetType(A00, &type));
776: PetscCall(MatCreate(comm, &A00_inv));
777: PetscCall(MatSetType(A00_inv, type));
778: PetscCall(MatInvertBlockDiagonalMat(A00, A00_inv));
779: PetscCall(MatMatMult(A00_inv, A01, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &AdB));
780: PetscCall(MatDestroy(&A00_inv));
781: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown MatSchurComplementAinvType: %d", ainvtype);
782: /* Cannot really reuse Sp in MatMatMult() because of MatAYPX() -->
783: MatAXPY() --> MatHeaderReplace() --> MatDestroy_XXX_MatMatMult() */
784: if (preuse == MAT_REUSE_MATRIX) PetscCall(MatDestroy(Sp));
785: PetscCall(MatMatMult(A10, AdB, MAT_INITIAL_MATRIX, PETSC_DEFAULT, Sp));
786: if (!A11) {
787: PetscCall(MatScale(*Sp, -1.0));
788: } else {
789: /* TODO: when can we pass SAME_NONZERO_PATTERN? */
790: PetscCall(MatAYPX(*Sp, -1, A11, DIFFERENT_NONZERO_PATTERN));
791: }
792: PetscCall(MatDestroy(&AdB));
793: }
794: PetscFunctionReturn(PETSC_SUCCESS);
795: }
797: static PetscErrorCode MatSchurComplementGetPmat_Basic(Mat S, MatReuse preuse, Mat *Sp)
798: {
799: Mat A, B, C, D;
800: Mat_SchurComplement *schur = (Mat_SchurComplement *)S->data;
802: PetscFunctionBegin;
803: if (preuse == MAT_IGNORE_MATRIX) PetscFunctionReturn(PETSC_SUCCESS);
804: PetscCall(MatSchurComplementGetSubMatrices(S, &A, NULL, &B, &C, &D));
805: PetscCheck(A, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONGSTATE, "Schur complement component matrices unset");
806: if (schur->ainvtype != MAT_SCHUR_COMPLEMENT_AINV_FULL) PetscCall(MatCreateSchurComplementPmat(A, B, C, D, schur->ainvtype, preuse, Sp));
807: else {
808: if (preuse == MAT_REUSE_MATRIX) PetscCall(MatDestroy(Sp));
809: PetscCall(MatSchurComplementComputeExplicitOperator(S, Sp));
810: }
811: PetscFunctionReturn(PETSC_SUCCESS);
812: }
814: /*@
815: MatSchurComplementGetPmat - Obtain a preconditioning matrix for the Schur complement by assembling Sp = A11 - A10 inv(DIAGFORM(A00)) A01
817: Collective
819: Input Parameters:
820: + S - matrix obtained with MatCreateSchurComplement() (or equivalent) that implements the action of A11 - A10 ksp(A00,Ap00) A01
821: - 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
823: Output Parameter:
824: . Sp - approximate Schur complement suitable for preconditioning the exact Schur complement S = A11 - A10 inv(A00) A01
826: Level: advanced
828: Notes:
829: The approximation of Sp depends on the argument passed to `MatSchurComplementSetAinvType()`
830: `MAT_SCHUR_COMPLEMENT_AINV_DIAG`, `MAT_SCHUR_COMPLEMENT_AINV_LUMP`, `MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG`, or `MAT_SCHUR_COMPLEMENT_AINV_FULL`
831: -mat_schur_complement_ainv_type <diag,lump,blockdiag,full>
833: Sometimes users would like to provide problem-specific data in the Schur complement, usually only
834: for special row and column index sets. In that case, the user should call `PetscObjectComposeFunction()` to set
835: "MatSchurComplementGetPmat_C" to their function. If their function needs to fall back to the default implementation,
836: it should call `MatSchurComplementGetPmat_Basic()`.
838: Developer Notes:
839: The API that includes `MatGetSchurComplement()`, `MatCreateSchurComplement()`, `MatSchurComplementGetPmat()` should be refactored to
840: remove redundancy and be clearer and simpler.
842: This routine should be called `MatSchurComplementCreatePmat()`
844: .seealso: [](ch_ksp), `MatCreateSubMatrix()`, `PCFIELDSPLIT`, `MatGetSchurComplement()`, `MatCreateSchurComplement()`, `MatSchurComplementSetAinvType()`
845: @*/
846: PetscErrorCode MatSchurComplementGetPmat(Mat S, MatReuse preuse, Mat *Sp)
847: {
848: PetscErrorCode (*f)(Mat, MatReuse, Mat *);
850: PetscFunctionBegin;
854: if (preuse != MAT_IGNORE_MATRIX) {
855: PetscAssertPointer(Sp, 3);
856: if (preuse == MAT_INITIAL_MATRIX) *Sp = NULL;
858: }
859: PetscCheck(!S->factortype, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
861: PetscCall(PetscObjectQueryFunction((PetscObject)S, "MatSchurComplementGetPmat_C", &f));
862: if (f) PetscCall((*f)(S, preuse, Sp));
863: else PetscCall(MatSchurComplementGetPmat_Basic(S, preuse, Sp));
864: PetscFunctionReturn(PETSC_SUCCESS);
865: }
867: static PetscErrorCode MatProductNumeric_SchurComplement_Dense(Mat C)
868: {
869: Mat_Product *product = C->product;
870: Mat_SchurComplement *Na = (Mat_SchurComplement *)product->A->data;
871: Mat work1, work2;
872: PetscScalar *v;
873: PetscInt lda;
875: PetscFunctionBegin;
876: PetscCall(MatMatMult(Na->B, product->B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &work1));
877: PetscCall(MatDuplicate(work1, MAT_DO_NOT_COPY_VALUES, &work2));
878: PetscCall(KSPMatSolve(Na->ksp, work1, work2));
879: PetscCall(MatDestroy(&work1));
880: PetscCall(MatDenseGetArrayWrite(C, &v));
881: PetscCall(MatDenseGetLDA(C, &lda));
882: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)C), C->rmap->n, C->cmap->n, C->rmap->N, C->cmap->N, v, &work1));
883: PetscCall(MatDenseSetLDA(work1, lda));
884: PetscCall(MatMatMult(Na->C, work2, MAT_REUSE_MATRIX, PETSC_DEFAULT, &work1));
885: PetscCall(MatDenseRestoreArrayWrite(C, &v));
886: PetscCall(MatDestroy(&work2));
887: PetscCall(MatDestroy(&work1));
888: if (Na->D) {
889: PetscCall(MatMatMult(Na->D, product->B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &work1));
890: PetscCall(MatAYPX(C, -1.0, work1, SAME_NONZERO_PATTERN));
891: PetscCall(MatDestroy(&work1));
892: } else PetscCall(MatScale(C, -1.0));
893: PetscFunctionReturn(PETSC_SUCCESS);
894: }
896: static PetscErrorCode MatProductSymbolic_SchurComplement_Dense(Mat C)
897: {
898: Mat_Product *product = C->product;
899: Mat A = product->A, B = product->B;
900: PetscInt m = A->rmap->n, n = B->cmap->n, M = A->rmap->N, N = B->cmap->N;
901: PetscBool flg;
903: PetscFunctionBegin;
904: PetscCall(MatSetSizes(C, m, n, M, N));
905: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)C, &flg, MATSEQDENSE, MATMPIDENSE, ""));
906: if (!flg) {
907: PetscCall(MatSetType(C, ((PetscObject)B)->type_name));
908: C->ops->productsymbolic = MatProductSymbolic_SchurComplement_Dense;
909: }
910: PetscCall(MatSetUp(C));
911: C->ops->productnumeric = MatProductNumeric_SchurComplement_Dense;
912: PetscFunctionReturn(PETSC_SUCCESS);
913: }
915: static PetscErrorCode MatProductSetFromOptions_Dense_AB(Mat C)
916: {
917: PetscFunctionBegin;
918: C->ops->productsymbolic = MatProductSymbolic_SchurComplement_Dense;
919: PetscFunctionReturn(PETSC_SUCCESS);
920: }
922: static PetscErrorCode MatProductSetFromOptions_SchurComplement_Dense(Mat C)
923: {
924: Mat_Product *product = C->product;
926: PetscFunctionBegin;
927: PetscCheck(product->type == MATPRODUCT_AB, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Not for product type %s", MatProductTypes[product->type]);
928: PetscCall(MatProductSetFromOptions_Dense_AB(C));
929: PetscFunctionReturn(PETSC_SUCCESS);
930: }
932: /*MC
933: MATSCHURCOMPLEMENT - "schurcomplement" - Matrix type that behaves like the Schur complement of a matrix.
935: Level: intermediate
937: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatType`, `MatCreateSchurComplement()`, `MatSchurComplementComputeExplicitOperator()`,
938: `MatSchurComplementGetSubMatrices()`, `MatSchurComplementGetKSP()`
939: M*/
940: PETSC_EXTERN PetscErrorCode MatCreate_SchurComplement(Mat N)
941: {
942: Mat_SchurComplement *Na;
944: PetscFunctionBegin;
945: PetscCall(PetscNew(&Na));
946: N->data = (void *)Na;
948: N->ops->destroy = MatDestroy_SchurComplement;
949: N->ops->getvecs = MatCreateVecs_SchurComplement;
950: N->ops->view = MatView_SchurComplement;
951: N->ops->mult = MatMult_SchurComplement;
952: N->ops->multtranspose = MatMultTranspose_SchurComplement;
953: N->ops->multadd = MatMultAdd_SchurComplement;
954: N->ops->setfromoptions = MatSetFromOptions_SchurComplement;
955: N->assembled = PETSC_FALSE;
956: N->preallocated = PETSC_FALSE;
958: PetscCall(KSPCreate(PetscObjectComm((PetscObject)N), &Na->ksp));
959: PetscCall(PetscObjectChangeTypeName((PetscObject)N, MATSCHURCOMPLEMENT));
960: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_schurcomplement_seqdense_C", MatProductSetFromOptions_SchurComplement_Dense));
961: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_schurcomplement_mpidense_C", MatProductSetFromOptions_SchurComplement_Dense));
962: PetscFunctionReturn(PETSC_SUCCESS);
963: }