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 solver for A00 block viewable with the additional option -%sksp_view\n", ((PetscObject)Na->ksp)->prefix ? ((PetscObject)Na->ksp)->prefix : NULL));
38: PetscCall(PetscViewerASCIIPrintf(viewer, "A01\n"));
39: PetscCall(PetscViewerASCIIPushTab(viewer));
40: PetscCall(MatView(Na->B, viewer));
41: PetscCall(PetscViewerASCIIPopTab(viewer));
42: PetscFunctionReturn(PETSC_SUCCESS);
43: }
45: /*
46: A11^T - A01^T ksptrans(A00,Ap00) A10^T
47: */
48: PetscErrorCode MatMultTranspose_SchurComplement(Mat N, Vec x, Vec y)
49: {
50: Mat_SchurComplement *Na = (Mat_SchurComplement *)N->data;
52: PetscFunctionBegin;
53: if (!Na->work1) PetscCall(MatCreateVecs(Na->A, &Na->work1, NULL));
54: if (!Na->work2) PetscCall(MatCreateVecs(Na->A, &Na->work2, NULL));
55: PetscCall(MatMultTranspose(Na->C, x, Na->work1));
56: PetscCall(KSPSolveTranspose(Na->ksp, Na->work1, Na->work2));
57: PetscCall(MatMultTranspose(Na->B, Na->work2, y));
58: PetscCall(VecScale(y, -1.0));
59: if (Na->D) PetscCall(MatMultTransposeAdd(Na->D, x, y, y));
60: PetscFunctionReturn(PETSC_SUCCESS);
61: }
63: /*
64: A11 - A10 ksp(A00,Ap00) A01
65: */
66: PetscErrorCode MatMult_SchurComplement(Mat N, Vec x, Vec y)
67: {
68: Mat_SchurComplement *Na = (Mat_SchurComplement *)N->data;
70: PetscFunctionBegin;
71: if (!Na->work1) PetscCall(MatCreateVecs(Na->A, &Na->work1, NULL));
72: if (!Na->work2) PetscCall(MatCreateVecs(Na->A, &Na->work2, NULL));
73: PetscCall(MatMult(Na->B, x, Na->work1));
74: PetscCall(KSPSolve(Na->ksp, Na->work1, Na->work2));
75: PetscCall(MatMult(Na->C, Na->work2, y));
76: PetscCall(VecScale(y, -1.0));
77: if (Na->D) PetscCall(MatMultAdd(Na->D, x, y, y));
78: PetscFunctionReturn(PETSC_SUCCESS);
79: }
81: /*
82: A11 - A10 ksp(A00,Ap00) A01
83: */
84: PetscErrorCode MatMultAdd_SchurComplement(Mat N, Vec x, Vec y, Vec z)
85: {
86: Mat_SchurComplement *Na = (Mat_SchurComplement *)N->data;
88: PetscFunctionBegin;
89: if (!Na->work1) PetscCall(MatCreateVecs(Na->A, &Na->work1, NULL));
90: if (!Na->work2) PetscCall(MatCreateVecs(Na->A, &Na->work2, NULL));
91: PetscCall(MatMult(Na->B, x, Na->work1));
92: PetscCall(KSPSolve(Na->ksp, Na->work1, Na->work2));
93: if (y == z) {
94: PetscCall(VecScale(Na->work2, -1.0));
95: PetscCall(MatMultAdd(Na->C, Na->work2, z, z));
96: } else {
97: PetscCall(MatMult(Na->C, Na->work2, z));
98: PetscCall(VecAYPX(z, -1.0, y));
99: }
100: if (Na->D) PetscCall(MatMultAdd(Na->D, x, z, z));
101: PetscFunctionReturn(PETSC_SUCCESS);
102: }
104: PetscErrorCode MatSetFromOptions_SchurComplement(Mat N, PetscOptionItems *PetscOptionsObject)
105: {
106: Mat_SchurComplement *Na = (Mat_SchurComplement *)N->data;
108: PetscFunctionBegin;
109: PetscOptionsHeadBegin(PetscOptionsObject, "MatSchurComplementOptions");
110: Na->ainvtype = MAT_SCHUR_COMPLEMENT_AINV_DIAG;
111: 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,
112: (PetscEnum *)&Na->ainvtype, NULL));
113: PetscOptionsHeadEnd();
114: PetscCall(KSPSetFromOptions(Na->ksp));
115: PetscFunctionReturn(PETSC_SUCCESS);
116: }
118: PetscErrorCode MatDestroy_SchurComplement(Mat N)
119: {
120: Mat_SchurComplement *Na = (Mat_SchurComplement *)N->data;
122: PetscFunctionBegin;
123: PetscCall(MatDestroy(&Na->A));
124: PetscCall(MatDestroy(&Na->Ap));
125: PetscCall(MatDestroy(&Na->B));
126: PetscCall(MatDestroy(&Na->C));
127: PetscCall(MatDestroy(&Na->D));
128: PetscCall(VecDestroy(&Na->work1));
129: PetscCall(VecDestroy(&Na->work2));
130: PetscCall(KSPDestroy(&Na->ksp));
131: PetscCall(PetscFree(N->data));
132: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_schurcomplement_seqdense_C", NULL));
133: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_schurcomplement_mpidense_C", NULL));
134: PetscFunctionReturn(PETSC_SUCCESS);
135: }
137: /*@
138: MatCreateSchurComplement - Creates a new `Mat` that behaves like the Schur complement of a matrix
140: Collective
142: Input Parameters:
143: + A00 - the upper-left block of the original matrix A = [A00 A01; A10 A11]
144: . Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A00^{-1}
145: . A01 - the upper-right block of the original matrix A = [A00 A01; A10 A11]
146: . A10 - the lower-left block of the original matrix A = [A00 A01; A10 A11]
147: - A11 - (optional) the lower-right block of the original matrix A = [A00 A01; A10 A11]
149: Output Parameter:
150: . S - the matrix that behaves as the Schur complement S = A11 - A10 ksp(A00,Ap00) A01
152: Level: intermediate
154: Notes:
155: The Schur complement is NOT explicitly formed! Rather, this function returns a virtual Schur complement
156: that can compute the matrix-vector product by using formula S = A11 - A10 A^{-1} A01
157: for Schur complement S and a `KSP` solver to approximate the action of A^{-1}.
159: All four matrices must have the same MPI communicator.
161: `A00` and `A11` must be square matrices.
163: `MatGetSchurComplement()` takes as arguments the index sets for the submatrices and returns both the virtual Schur complement (what this returns) plus
164: a sparse approximation to the Schur complement (useful for building a preconditioner for the Schur complement) which can be obtained from this
165: matrix with `MatSchurComplementGetPmat()`
167: Developer Notes:
168: The API that includes `MatGetSchurComplement()`, `MatCreateSchurComplement()`, `MatSchurComplementGetPmat()` should be refactored to
169: remove redundancy and be clearer and simpler.
171: .seealso: [](ch_ksp), `MatCreateNormal()`, `MatMult()`, `MatCreate()`, `MatSchurComplementGetKSP()`, `MatSchurComplementUpdateSubMatrices()`, `MatCreateTranspose()`, `MatGetSchurComplement()`,
172: `MatSchurComplementGetPmat()`, `MatSchurComplementSetSubMatrices()`
173: @*/
174: PetscErrorCode MatCreateSchurComplement(Mat A00, Mat Ap00, Mat A01, Mat A10, Mat A11, Mat *S)
175: {
176: PetscFunctionBegin;
177: PetscCall(KSPInitializePackage());
178: PetscCall(MatCreate(PetscObjectComm((PetscObject)A00), S));
179: PetscCall(MatSetType(*S, MATSCHURCOMPLEMENT));
180: PetscCall(MatSchurComplementSetSubMatrices(*S, A00, Ap00, A01, A10, A11));
181: PetscFunctionReturn(PETSC_SUCCESS);
182: }
184: /*@
185: MatSchurComplementSetSubMatrices - Sets the matrices that define the Schur complement
187: Collective
189: Input Parameters:
190: + S - matrix obtained with `MatSetType`(S,`MATSCHURCOMPLEMENT`)
191: . A00 - the upper-left block of the original matrix A = [A00 A01; A10 A11]
192: . Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A00^{-1}
193: . A01 - the upper-right block of the original matrix A = [A00 A01; A10 A11]
194: . A10 - the lower-left block of the original matrix A = [A00 A01; A10 A11]
195: - A11 - (optional) the lower-right block of the original matrix A = [A00 A01; A10 A11]
197: Level: intermediate
199: Notes:
200: The Schur complement is NOT explicitly formed! Rather, this
201: object performs the matrix-vector product of the Schur complement by using formula S = A11 - A10 ksp(A00,Ap00) A01
203: All four matrices must have the same MPI communicator.
205: `A00` and `A11` must be square matrices.
207: This is to be used in the context of code such as
208: .vb
209: MatSetType(S,MATSCHURCOMPLEMENT);
210: MatSchurComplementSetSubMatrices(S,...);
211: .ve
212: while `MatSchurComplementUpdateSubMatrices()` should only be called after `MatCreateSchurComplement()` or `MatSchurComplementSetSubMatrices()`
214: .seealso: [](ch_ksp), `Mat`, `MatCreateNormal()`, `MatMult()`, `MatCreate()`, `MatSchurComplementGetKSP()`, `MatSchurComplementUpdateSubMatrices()`, `MatCreateTranspose()`, `MatCreateSchurComplement()`, `MatGetSchurComplement()`
215: @*/
216: PetscErrorCode MatSchurComplementSetSubMatrices(Mat S, Mat A00, Mat Ap00, Mat A01, Mat A10, Mat A11)
217: {
218: Mat_SchurComplement *Na = (Mat_SchurComplement *)S->data;
219: PetscBool isschur;
221: PetscFunctionBegin;
222: PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
223: if (!isschur) PetscFunctionReturn(PETSC_SUCCESS);
224: PetscCheck(!S->assembled, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONGSTATE, "Use MatSchurComplementUpdateSubMatrices() for already used matrix");
229: PetscCheckSameComm(A00, 2, Ap00, 3);
230: PetscCheckSameComm(A00, 2, A01, 4);
231: PetscCheckSameComm(A00, 2, A10, 5);
232: 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);
233: 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);
234: 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);
235: 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);
236: 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);
237: if (A11) {
239: PetscCheckSameComm(A00, 2, A11, 6);
240: 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);
241: }
243: PetscCall(MatSetSizes(S, A10->rmap->n, A01->cmap->n, A10->rmap->N, A01->cmap->N));
244: PetscCall(PetscObjectReference((PetscObject)A00));
245: PetscCall(PetscObjectReference((PetscObject)Ap00));
246: PetscCall(PetscObjectReference((PetscObject)A01));
247: PetscCall(PetscObjectReference((PetscObject)A10));
248: Na->A = A00;
249: Na->Ap = Ap00;
250: Na->B = A01;
251: Na->C = A10;
252: Na->D = A11;
253: if (A11) PetscCall(PetscObjectReference((PetscObject)A11));
254: PetscCall(MatSetUp(S));
255: PetscCall(KSPSetOperators(Na->ksp, A00, Ap00));
256: S->assembled = PETSC_TRUE;
257: PetscFunctionReturn(PETSC_SUCCESS);
258: }
260: /*@
261: MatSchurComplementGetKSP - Gets the `KSP` object that is used to solve with A00 in the Schur complement matrix S = A11 - A10 ksp(A00,Ap00) A01
263: Not Collective
265: Input Parameter:
266: . S - matrix obtained with `MatCreateSchurComplement()` (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
268: Output Parameter:
269: . ksp - the linear solver object
271: Options Database Key:
272: . -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.
274: Level: intermediate
276: .seealso: [](ch_ksp), `Mat`, `MatSchurComplementSetKSP()`, `MatCreateSchurComplement()`, `MatCreateNormal()`, `MatMult()`, `MatCreate()`
277: @*/
278: PetscErrorCode MatSchurComplementGetKSP(Mat S, KSP *ksp)
279: {
280: Mat_SchurComplement *Na;
281: PetscBool isschur;
283: PetscFunctionBegin;
285: PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
286: PetscCheck(isschur, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)S)->type_name);
287: PetscAssertPointer(ksp, 2);
288: Na = (Mat_SchurComplement *)S->data;
289: *ksp = Na->ksp;
290: PetscFunctionReturn(PETSC_SUCCESS);
291: }
293: /*@
294: MatSchurComplementSetKSP - Sets the `KSP` object that is used to solve with A00 in the Schur complement matrix S = A11 - A10 ksp(A00,Ap00) A01
296: Not Collective
298: Input Parameters:
299: + S - matrix created with `MatCreateSchurComplement()`
300: - ksp - the linear solver object
302: Level: developer
304: Developer Notes:
305: This is used in `PCFIELDSPLIT` to reuse the 0-split `KSP` to implement ksp(A00,Ap00) in S.
306: The `KSP` operators are overwritten with A00 and Ap00 currently set in S.
308: .seealso: [](ch_ksp), `Mat`, `MatSchurComplementGetKSP()`, `MatCreateSchurComplement()`, `MatCreateNormal()`, `MatMult()`, `MatCreate()`, `MATSCHURCOMPLEMENT`
309: @*/
310: PetscErrorCode MatSchurComplementSetKSP(Mat S, KSP ksp)
311: {
312: Mat_SchurComplement *Na;
313: PetscBool isschur;
315: PetscFunctionBegin;
317: PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
318: if (!isschur) PetscFunctionReturn(PETSC_SUCCESS);
320: Na = (Mat_SchurComplement *)S->data;
321: PetscCall(PetscObjectReference((PetscObject)ksp));
322: PetscCall(KSPDestroy(&Na->ksp));
323: Na->ksp = ksp;
324: PetscCall(KSPSetOperators(Na->ksp, Na->A, Na->Ap));
325: PetscFunctionReturn(PETSC_SUCCESS);
326: }
328: /*@
329: MatSchurComplementUpdateSubMatrices - Updates the Schur complement matrix object with new submatrices
331: Collective
333: Input Parameters:
334: + S - matrix obtained with `MatCreateSchurComplement()` (or `MatSchurSetSubMatrices()`) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
335: . A00 - the upper-left block of the original matrix A = [A00 A01; A10 A11]
336: . Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A00^{-1}
337: . A01 - the upper-right block of the original matrix A = [A00 A01; A10 A11]
338: . A10 - the lower-left block of the original matrix A = [A00 A01; A10 A11]
339: - A11 - (optional) the lower-right block of the original matrix A = [A00 A01; A10 A11]
341: Level: intermediate
343: Notes:
344: All four matrices must have the same MPI communicator
346: `A00` and `A11` must be square matrices
348: All of the matrices provided must have the same sizes as was used with `MatCreateSchurComplement()` or `MatSchurComplementSetSubMatrices()`
349: though they need not be the same matrices.
351: This can only be called after `MatCreateSchurComplement()` or `MatSchurComplementSetSubMatrices()`, it cannot be called immediately after `MatSetType`(S,`MATSCHURCOMPLEMENT`);
353: Developer Notes:
354: This code is almost identical to `MatSchurComplementSetSubMatrices()`. The API should be refactored.
356: .seealso: [](ch_ksp), `Mat`, `MatCreateNormal()`, `MatMult()`, `MatCreate()`, `MatSchurComplementGetKSP()`, `MatCreateSchurComplement()`
357: @*/
358: PetscErrorCode MatSchurComplementUpdateSubMatrices(Mat S, Mat A00, Mat Ap00, Mat A01, Mat A10, Mat A11)
359: {
360: Mat_SchurComplement *Na = (Mat_SchurComplement *)S->data;
361: PetscBool isschur;
363: PetscFunctionBegin;
365: PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
366: if (!isschur) PetscFunctionReturn(PETSC_SUCCESS);
367: PetscCheck(S->assembled, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONGSTATE, "Use MatSchurComplementSetSubMatrices() for a new matrix");
372: PetscCheckSameComm(A00, 2, Ap00, 3);
373: PetscCheckSameComm(A00, 2, A01, 4);
374: PetscCheckSameComm(A00, 2, A10, 5);
375: 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);
376: 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);
377: 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);
378: 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);
379: 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);
380: if (A11) {
382: PetscCheckSameComm(A00, 2, A11, 6);
383: 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);
384: }
386: PetscCall(PetscObjectReference((PetscObject)A00));
387: PetscCall(PetscObjectReference((PetscObject)Ap00));
388: PetscCall(PetscObjectReference((PetscObject)A01));
389: PetscCall(PetscObjectReference((PetscObject)A10));
390: if (A11) PetscCall(PetscObjectReference((PetscObject)A11));
392: PetscCall(MatDestroy(&Na->A));
393: PetscCall(MatDestroy(&Na->Ap));
394: PetscCall(MatDestroy(&Na->B));
395: PetscCall(MatDestroy(&Na->C));
396: PetscCall(MatDestroy(&Na->D));
398: Na->A = A00;
399: Na->Ap = Ap00;
400: Na->B = A01;
401: Na->C = A10;
402: Na->D = A11;
404: PetscCall(KSPSetOperators(Na->ksp, A00, Ap00));
405: PetscFunctionReturn(PETSC_SUCCESS);
406: }
408: /*@C
409: MatSchurComplementGetSubMatrices - Get the individual submatrices in the Schur complement
411: Collective
413: Input Parameter:
414: . S - matrix obtained with `MatCreateSchurComplement()` (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
416: Output Parameters:
417: + A00 - the upper-left block of the original matrix A = [A00 A01; A10 A11]
418: . Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}
419: . A01 - the upper-right block of the original matrix A = [A00 A01; A10 A11]
420: . A10 - the lower-left block of the original matrix A = [A00 A01; A10 A11]
421: - A11 - (optional) the lower-right block of the original matrix A = [A00 A01; A10 A11]
423: Level: intermediate
425: Note:
426: Use `NULL` for any unneeded output argument.
428: The reference counts of the submatrices are not increased before they are returned and the matrices should not be modified or destroyed.
430: .seealso: [](ch_ksp), `MatCreateNormal()`, `MatMult()`, `MatCreate()`, `MatSchurComplementGetKSP()`, `MatCreateSchurComplement()`, `MatSchurComplementUpdateSubMatrices()`
431: @*/
432: PetscErrorCode MatSchurComplementGetSubMatrices(Mat S, Mat *A00, Mat *Ap00, Mat *A01, Mat *A10, Mat *A11)
433: {
434: Mat_SchurComplement *Na = (Mat_SchurComplement *)S->data;
435: PetscBool flg;
437: PetscFunctionBegin;
439: PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &flg));
440: PetscCheck(flg, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)S)->type_name);
441: if (A00) *A00 = Na->A;
442: if (Ap00) *Ap00 = Na->Ap;
443: if (A01) *A01 = Na->B;
444: if (A10) *A10 = Na->C;
445: if (A11) *A11 = Na->D;
446: PetscFunctionReturn(PETSC_SUCCESS);
447: }
449: #include <petsc/private/kspimpl.h>
451: /*@
452: MatSchurComplementComputeExplicitOperator - Compute the Schur complement matrix explicitly
454: Collective
456: Input Parameter:
457: . A - the matrix obtained with `MatCreateSchurComplement()`
459: Output Parameter:
460: . S - the Schur complement matrix
462: Notes:
463: This can be expensive, so it is mainly for testing
465: Use `MatSchurComplementGetPmat()` to get a sparse approximation for the Schur complement suitable for use in building a preconditioner
467: Level: advanced
469: .seealso: [](ch_ksp), `MatCreateSchurComplement()`, `MatSchurComplementUpdateSubMatrices()`, `MatSchurComplementGetPmat()`
470: @*/
471: PetscErrorCode MatSchurComplementComputeExplicitOperator(Mat A, Mat *S)
472: {
473: Mat B, C, D, E = NULL, Bd, AinvBd;
474: MatType mtype;
475: VecType vtype;
476: KSP ksp;
477: PetscInt n, N, m, M;
478: PetscBool flg = PETSC_FALSE, set, symm;
480: PetscFunctionBegin;
481: PetscCall(MatSchurComplementGetSubMatrices(A, NULL, NULL, &B, &C, &D));
482: PetscCall(MatSchurComplementGetKSP(A, &ksp));
483: PetscCall(KSPSetUp(ksp));
484: PetscCall(MatGetVecType(B, &vtype));
485: PetscCall(MatGetLocalSize(B, &m, &n));
486: PetscCall(MatGetSize(B, &M, &N));
487: PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)A), vtype, m, n, M, N, -1, NULL, &AinvBd));
488: PetscCall(MatGetType(AinvBd, &mtype));
489: PetscCall(MatConvert(B, mtype, MAT_INITIAL_MATRIX, &Bd));
490: PetscCall(KSPMatSolve(ksp, Bd, AinvBd));
491: PetscCall(MatDestroy(&Bd));
492: PetscCall(MatFilter(AinvBd, PETSC_SMALL, PETSC_FALSE, PETSC_FALSE));
493: if (D) {
494: PetscCall(MatGetLocalSize(D, &m, &n));
495: PetscCall(MatGetSize(D, &M, &N));
496: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A), m, n, M, N, NULL, S));
497: }
498: PetscCall(MatMatMult(C, AinvBd, D ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX, PETSC_DEFAULT, S));
499: PetscCall(MatDestroy(&AinvBd));
500: if (D) {
501: PetscCall(PetscObjectTypeCompareAny((PetscObject)D, &flg, MATSEQSBAIJ, MATMPISBAIJ, ""));
502: if (flg) {
503: PetscCall(MatIsSymmetricKnown(A, &set, &symm));
504: if (!set || !symm) PetscCall(MatConvert(D, MATBAIJ, MAT_INITIAL_MATRIX, &E)); /* convert the (1,1) block to nonsymmetric storage for MatAXPY() */
505: }
506: 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 */
507: }
508: 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 */
509: PetscCall(MatScale(*S, -1.0));
510: PetscCall(MatDestroy(&E));
511: PetscFunctionReturn(PETSC_SUCCESS);
512: }
514: /* Developer Notes:
515: This should be implemented with a MatCreate_SchurComplement() as that is the standard design for new Mat classes. */
516: PetscErrorCode MatGetSchurComplement_Basic(Mat mat, IS isrow0, IS iscol0, IS isrow1, IS iscol1, MatReuse mreuse, Mat *S, MatSchurComplementAinvType ainvtype, MatReuse preuse, Mat *Sp)
517: {
518: Mat A = NULL, Ap = NULL, B = NULL, C = NULL, D = NULL;
519: MatReuse reuse;
521: PetscFunctionBegin;
531: if (mreuse == MAT_IGNORE_MATRIX && preuse == MAT_IGNORE_MATRIX) PetscFunctionReturn(PETSC_SUCCESS);
535: PetscCheck(!mat->factortype, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
537: reuse = MAT_INITIAL_MATRIX;
538: if (mreuse == MAT_REUSE_MATRIX) {
539: PetscCall(MatSchurComplementGetSubMatrices(*S, &A, &Ap, &B, &C, &D));
540: PetscCheck(A && Ap && B && C, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Attempting to reuse matrix but Schur complement matrices unset");
541: PetscCheck(A == Ap, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Preconditioning matrix does not match operator");
542: PetscCall(MatDestroy(&Ap)); /* get rid of extra reference */
543: reuse = MAT_REUSE_MATRIX;
544: }
545: PetscCall(MatCreateSubMatrix(mat, isrow0, iscol0, reuse, &A));
546: PetscCall(MatCreateSubMatrix(mat, isrow0, iscol1, reuse, &B));
547: PetscCall(MatCreateSubMatrix(mat, isrow1, iscol0, reuse, &C));
548: PetscCall(MatCreateSubMatrix(mat, isrow1, iscol1, reuse, &D));
549: switch (mreuse) {
550: case MAT_INITIAL_MATRIX:
551: PetscCall(MatCreateSchurComplement(A, A, B, C, D, S));
552: break;
553: case MAT_REUSE_MATRIX:
554: PetscCall(MatSchurComplementUpdateSubMatrices(*S, A, A, B, C, D));
555: break;
556: default:
557: PetscCheck(mreuse == MAT_IGNORE_MATRIX, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Unrecognized value of mreuse %d", (int)mreuse);
558: }
559: if (preuse != MAT_IGNORE_MATRIX) PetscCall(MatCreateSchurComplementPmat(A, B, C, D, ainvtype, preuse, Sp));
560: PetscCall(MatDestroy(&A));
561: PetscCall(MatDestroy(&B));
562: PetscCall(MatDestroy(&C));
563: PetscCall(MatDestroy(&D));
564: PetscFunctionReturn(PETSC_SUCCESS);
565: }
567: /*@
568: MatGetSchurComplement - Obtain the Schur complement from eliminating part of the matrix in another part.
570: Collective
572: Input Parameters:
573: + A - matrix in which the complement is to be taken
574: . isrow0 - rows to eliminate
575: . iscol0 - columns to eliminate, (isrow0,iscol0) should be square and nonsingular
576: . isrow1 - rows in which the Schur complement is formed
577: . iscol1 - columns in which the Schur complement is formed
578: . mreuse - `MAT_INITIAL_MATRIX` or `MAT_REUSE_MATRIX`, use `MAT_IGNORE_MATRIX` to put nothing in S
579: . ainvtype - the type of approximation used for the inverse of the (0,0) block used in forming Sp:
580: `MAT_SCHUR_COMPLEMENT_AINV_DIAG`, `MAT_SCHUR_COMPLEMENT_AINV_LUMP`, `MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG`, or `MAT_SCHUR_COMPLEMENT_AINV_FULL`
581: - preuse - `MAT_INITIAL_MATRIX` or `MAT_REUSE_MATRIX`, use `MAT_IGNORE_MATRIX` to put nothing in Sp
583: Output Parameters:
584: + S - exact Schur complement, often of type `MATSCHURCOMPLEMENT` which is difficult to use for preconditioning
585: - Sp - approximate Schur complement from which a preconditioner can be built A11 - A10 inv(DIAGFORM(A00)) A01
587: Level: advanced
589: Notes:
590: Since the real Schur complement is usually dense, providing a good approximation to Sp usually requires
591: application-specific information.
593: Sometimes users would like to provide problem-specific data in the Schur complement, usually only for special row
594: 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
595: "MatGetSchurComplement_C" to their function. If their function needs to fall back to the default implementation, it
596: should call `MatGetSchurComplement_Basic()`.
598: `MatCreateSchurComplement()` takes as arguments the four submatrices and returns the virtual Schur complement (what this function returns in S).
600: `MatSchurComplementGetPmat()` takes the virtual Schur complement and returns an explicit approximate Schur complement (what this returns in Sp).
602: In other words calling `MatCreateSchurComplement()` followed by `MatSchurComplementGetPmat()` produces the same output as this function but with slightly different
603: inputs. The actually submatrices of the original block matrix instead of index sets to the submatrices.
605: Developer Notes:
606: The API that includes `MatGetSchurComplement()`, `MatCreateSchurComplement()`, `MatSchurComplementGetPmat()` should be refactored to
607: remove redundancy and be clearer and simpler.
609: .seealso: [](ch_ksp), `MatCreateSubMatrix()`, `PCFIELDSPLIT`, `MatCreateSchurComplement()`, `MatSchurComplementAinvType`
610: @*/
611: PetscErrorCode MatGetSchurComplement(Mat A, IS isrow0, IS iscol0, IS isrow1, IS iscol1, MatReuse mreuse, Mat *S, MatSchurComplementAinvType ainvtype, MatReuse preuse, Mat *Sp)
612: {
613: PetscErrorCode (*f)(Mat, IS, IS, IS, IS, MatReuse, Mat *, MatReuse, Mat *) = NULL;
615: PetscFunctionBegin;
627: PetscCheck(!A->factortype, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
628: 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. */
629: PetscCall(PetscObjectQueryFunction((PetscObject)*S, "MatGetSchurComplement_C", &f));
630: }
631: if (f) PetscCall((*f)(A, isrow0, iscol0, isrow1, iscol1, mreuse, S, preuse, Sp));
632: else PetscCall(MatGetSchurComplement_Basic(A, isrow0, iscol0, isrow1, iscol1, mreuse, S, ainvtype, preuse, Sp));
633: PetscFunctionReturn(PETSC_SUCCESS);
634: }
636: /*@
637: MatSchurComplementSetAinvType - set the type of approximation used for the inverse of the (0,0) block used in forming Sp in `MatSchurComplementGetPmat()`
639: Not Collective
641: Input Parameters:
642: + S - matrix obtained with `MatCreateSchurComplement()` (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
643: - ainvtype - type of approximation to be used to form approximate Schur complement Sp = A11 - A10 inv(DIAGFORM(A00)) A01:
644: `MAT_SCHUR_COMPLEMENT_AINV_DIAG`, `MAT_SCHUR_COMPLEMENT_AINV_LUMP`, `MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG`, or `MAT_SCHUR_COMPLEMENT_AINV_FULL`
646: Options Database Key:
647: . -mat_schur_complement_ainv_type diag | lump | blockdiag | full - set schur complement type
649: Level: advanced
651: .seealso: [](ch_ksp), `MatSchurComplementAinvType`, `MatCreateSchurComplement()`, `MatGetSchurComplement()`, `MatSchurComplementGetPmat()`, `MatSchurComplementGetAinvType()`
652: @*/
653: PetscErrorCode MatSchurComplementSetAinvType(Mat S, MatSchurComplementAinvType ainvtype)
654: {
655: PetscBool isschur;
656: Mat_SchurComplement *schur;
658: PetscFunctionBegin;
660: PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
661: if (!isschur) PetscFunctionReturn(PETSC_SUCCESS);
663: schur = (Mat_SchurComplement *)S->data;
664: 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);
665: schur->ainvtype = ainvtype;
666: PetscFunctionReturn(PETSC_SUCCESS);
667: }
669: /*@
670: MatSchurComplementGetAinvType - get the type of approximation for the inverse of the (0,0) block used in forming Sp in `MatSchurComplementGetPmat()`
672: Not Collective
674: Input Parameter:
675: . S - matrix obtained with `MatCreateSchurComplement()` (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
677: Output Parameter:
678: . ainvtype - type of approximation used to form approximate Schur complement Sp = A11 - A10 inv(DIAGFORM(A00)) A01:
679: `MAT_SCHUR_COMPLEMENT_AINV_DIAG`, `MAT_SCHUR_COMPLEMENT_AINV_LUMP`, `MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG`, or `MAT_SCHUR_COMPLEMENT_AINV_FULL`
681: Level: advanced
683: .seealso: [](ch_ksp), `MatSchurComplementAinvType`, `MatCreateSchurComplement()`, `MatGetSchurComplement()`, `MatSchurComplementGetPmat()`, `MatSchurComplementSetAinvType()`
684: @*/
685: PetscErrorCode MatSchurComplementGetAinvType(Mat S, MatSchurComplementAinvType *ainvtype)
686: {
687: PetscBool isschur;
688: Mat_SchurComplement *schur;
690: PetscFunctionBegin;
692: PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
693: PetscCheck(isschur, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)S)->type_name);
694: schur = (Mat_SchurComplement *)S->data;
695: if (ainvtype) *ainvtype = schur->ainvtype;
696: PetscFunctionReturn(PETSC_SUCCESS);
697: }
699: /*@
700: MatCreateSchurComplementPmat - create a preconditioning matrix for the Schur complement by explicitly assembling the sparse matrix
701: Sp = A11 - A10 inv(DIAGFORM(A00)) A01
703: Collective
705: Input Parameters:
706: + A00 - the upper-left part of the original matrix A = [A00 A01; A10 A11]
707: . A01 - (optional) the upper-right part of the original matrix A = [A00 A01; A10 A11]
708: . A10 - (optional) the lower-left part of the original matrix A = [A00 A01; A10 A11]
709: . A11 - (optional) the lower-right part of the original matrix A = [A00 A01; A10 A11]
710: . ainvtype - type of approximation for DIAGFORM(A00) used when forming Sp = A11 - A10 inv(DIAGFORM(A00)) A01. See MatSchurComplementAinvType.
711: - 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
713: Output Parameter:
714: . Sp - approximate Schur complement suitable for preconditioning the true Schur complement S = A11 - A10 inv(A00) A01
716: Level: advanced
718: .seealso: [](ch_ksp), `MatCreateSchurComplement()`, `MatGetSchurComplement()`, `MatSchurComplementGetPmat()`, `MatSchurComplementAinvType`
719: @*/
720: PetscErrorCode MatCreateSchurComplementPmat(Mat A00, Mat A01, Mat A10, Mat A11, MatSchurComplementAinvType ainvtype, MatReuse preuse, Mat *Sp)
721: {
722: PetscInt N00;
724: PetscFunctionBegin;
725: /* 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. */
726: /* TODO: Perhaps should create an appropriately-sized zero matrix of the same type as A00? */
728: if (preuse == MAT_IGNORE_MATRIX) PetscFunctionReturn(PETSC_SUCCESS);
730: /* A zero size A00 or empty A01 or A10 imply S = A11. */
731: PetscCall(MatGetSize(A00, &N00, NULL));
732: if (!A01 || !A10 || !N00) {
733: if (preuse == MAT_INITIAL_MATRIX) {
734: PetscCall(MatDuplicate(A11, MAT_COPY_VALUES, Sp));
735: } else { /* MAT_REUSE_MATRIX */
736: /* TODO: when can we pass SAME_NONZERO_PATTERN? */
737: PetscCall(MatCopy(A11, *Sp, DIFFERENT_NONZERO_PATTERN));
738: }
739: } else {
740: Mat AdB, T;
741: Vec diag;
742: PetscBool flg;
744: if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP || ainvtype == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
745: PetscCall(PetscObjectTypeCompare((PetscObject)A01, MATTRANSPOSEVIRTUAL, &flg));
746: if (flg) {
747: PetscCall(MatTransposeGetMat(A01, &T));
748: PetscCall(MatTranspose(T, MAT_INITIAL_MATRIX, &AdB));
749: } else {
750: PetscCall(PetscObjectTypeCompare((PetscObject)A01, MATHERMITIANTRANSPOSEVIRTUAL, &flg));
751: if (flg) {
752: PetscCall(MatHermitianTransposeGetMat(A01, &T));
753: PetscCall(MatHermitianTranspose(T, MAT_INITIAL_MATRIX, &AdB));
754: }
755: }
756: if (!flg) PetscCall(MatDuplicate(A01, MAT_COPY_VALUES, &AdB));
757: else {
758: PetscScalar shift, scale;
760: PetscCall(MatShellGetScalingShifts(A01, &shift, &scale, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
761: PetscCall(MatShift(AdB, shift));
762: PetscCall(MatScale(AdB, scale));
763: }
764: PetscCall(MatCreateVecs(A00, &diag, NULL));
765: if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP) {
766: PetscCall(MatGetRowSum(A00, diag));
767: } else {
768: PetscCall(MatGetDiagonal(A00, diag));
769: }
770: PetscCall(VecReciprocal(diag));
771: PetscCall(MatDiagonalScale(AdB, diag, NULL));
772: PetscCall(VecDestroy(&diag));
773: } else if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG) {
774: Mat A00_inv;
775: MatType type;
776: MPI_Comm comm;
778: PetscCall(PetscObjectGetComm((PetscObject)A00, &comm));
779: PetscCall(MatGetType(A00, &type));
780: PetscCall(MatCreate(comm, &A00_inv));
781: PetscCall(MatSetType(A00_inv, type));
782: PetscCall(MatInvertBlockDiagonalMat(A00, A00_inv));
783: PetscCall(MatMatMult(A00_inv, A01, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &AdB));
784: PetscCall(MatDestroy(&A00_inv));
785: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown MatSchurComplementAinvType: %d", ainvtype);
786: /* Cannot really reuse Sp in MatMatMult() because of MatAYPX() -->
787: MatAXPY() --> MatHeaderReplace() --> MatDestroy_XXX_MatMatMult() */
788: if (preuse == MAT_REUSE_MATRIX) PetscCall(MatDestroy(Sp));
789: PetscCall(MatMatMult(A10, AdB, MAT_INITIAL_MATRIX, PETSC_DEFAULT, Sp));
790: if (!A11) {
791: PetscCall(MatScale(*Sp, -1.0));
792: } else {
793: /* TODO: when can we pass SAME_NONZERO_PATTERN? */
794: PetscCall(MatAYPX(*Sp, -1, A11, DIFFERENT_NONZERO_PATTERN));
795: }
796: PetscCall(MatDestroy(&AdB));
797: }
798: PetscFunctionReturn(PETSC_SUCCESS);
799: }
801: static PetscErrorCode MatSchurComplementGetPmat_Basic(Mat S, MatReuse preuse, Mat *Sp)
802: {
803: Mat A, B, C, D;
804: Mat_SchurComplement *schur = (Mat_SchurComplement *)S->data;
806: PetscFunctionBegin;
807: if (preuse == MAT_IGNORE_MATRIX) PetscFunctionReturn(PETSC_SUCCESS);
808: PetscCall(MatSchurComplementGetSubMatrices(S, &A, NULL, &B, &C, &D));
809: PetscCheck(A, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONGSTATE, "Schur complement component matrices unset");
810: if (schur->ainvtype != MAT_SCHUR_COMPLEMENT_AINV_FULL) PetscCall(MatCreateSchurComplementPmat(A, B, C, D, schur->ainvtype, preuse, Sp));
811: else {
812: if (preuse == MAT_REUSE_MATRIX) PetscCall(MatDestroy(Sp));
813: PetscCall(MatSchurComplementComputeExplicitOperator(S, Sp));
814: }
815: PetscFunctionReturn(PETSC_SUCCESS);
816: }
818: /*@
819: MatSchurComplementGetPmat - Obtain a preconditioning matrix for the Schur complement by assembling Sp = A11 - A10 inv(DIAGFORM(A00)) A01
821: Collective
823: Input Parameters:
824: + S - matrix obtained with MatCreateSchurComplement() (or equivalent) that implements the action of A11 - A10 ksp(A00,Ap00) A01
825: - 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
827: Output Parameter:
828: . Sp - approximate Schur complement suitable for preconditioning the exact Schur complement S = A11 - A10 inv(A00) A01
830: Level: advanced
832: Notes:
833: The approximation of Sp depends on the argument passed to `MatSchurComplementSetAinvType()`
834: `MAT_SCHUR_COMPLEMENT_AINV_DIAG`, `MAT_SCHUR_COMPLEMENT_AINV_LUMP`, `MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG`, or `MAT_SCHUR_COMPLEMENT_AINV_FULL`
835: -mat_schur_complement_ainv_type <diag,lump,blockdiag,full>
837: Sometimes users would like to provide problem-specific data in the Schur complement, usually only
838: for special row and column index sets. In that case, the user should call `PetscObjectComposeFunction()` to set
839: "MatSchurComplementGetPmat_C" to their function. If their function needs to fall back to the default implementation,
840: it should call `MatSchurComplementGetPmat_Basic()`.
842: Developer Notes:
843: The API that includes `MatGetSchurComplement()`, `MatCreateSchurComplement()`, `MatSchurComplementGetPmat()` should be refactored to
844: remove redundancy and be clearer and simpler.
846: This routine should be called `MatSchurComplementCreatePmat()`
848: .seealso: [](ch_ksp), `MatCreateSubMatrix()`, `PCFIELDSPLIT`, `MatGetSchurComplement()`, `MatCreateSchurComplement()`, `MatSchurComplementSetAinvType()`
849: @*/
850: PetscErrorCode MatSchurComplementGetPmat(Mat S, MatReuse preuse, Mat *Sp)
851: {
852: PetscErrorCode (*f)(Mat, MatReuse, Mat *);
854: PetscFunctionBegin;
858: if (preuse != MAT_IGNORE_MATRIX) {
859: PetscAssertPointer(Sp, 3);
860: if (preuse == MAT_INITIAL_MATRIX) *Sp = NULL;
862: }
863: PetscCheck(!S->factortype, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
865: PetscCall(PetscObjectQueryFunction((PetscObject)S, "MatSchurComplementGetPmat_C", &f));
866: if (f) PetscCall((*f)(S, preuse, Sp));
867: else PetscCall(MatSchurComplementGetPmat_Basic(S, preuse, Sp));
868: PetscFunctionReturn(PETSC_SUCCESS);
869: }
871: static PetscErrorCode MatProductNumeric_SchurComplement_Dense(Mat C)
872: {
873: Mat_Product *product = C->product;
874: Mat_SchurComplement *Na = (Mat_SchurComplement *)product->A->data;
875: Mat work1, work2;
876: PetscScalar *v;
877: PetscInt lda;
879: PetscFunctionBegin;
880: PetscCall(MatMatMult(Na->B, product->B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &work1));
881: PetscCall(MatDuplicate(work1, MAT_DO_NOT_COPY_VALUES, &work2));
882: PetscCall(KSPMatSolve(Na->ksp, work1, work2));
883: PetscCall(MatDestroy(&work1));
884: PetscCall(MatDenseGetArrayWrite(C, &v));
885: PetscCall(MatDenseGetLDA(C, &lda));
886: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)C), C->rmap->n, C->cmap->n, C->rmap->N, C->cmap->N, v, &work1));
887: PetscCall(MatDenseSetLDA(work1, lda));
888: PetscCall(MatMatMult(Na->C, work2, MAT_REUSE_MATRIX, PETSC_DEFAULT, &work1));
889: PetscCall(MatDenseRestoreArrayWrite(C, &v));
890: PetscCall(MatDestroy(&work2));
891: PetscCall(MatDestroy(&work1));
892: if (Na->D) {
893: PetscCall(MatMatMult(Na->D, product->B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &work1));
894: PetscCall(MatAYPX(C, -1.0, work1, SAME_NONZERO_PATTERN));
895: PetscCall(MatDestroy(&work1));
896: } else PetscCall(MatScale(C, -1.0));
897: PetscFunctionReturn(PETSC_SUCCESS);
898: }
900: static PetscErrorCode MatProductSymbolic_SchurComplement_Dense(Mat C)
901: {
902: Mat_Product *product = C->product;
903: Mat A = product->A, B = product->B;
904: PetscInt m = A->rmap->n, n = B->cmap->n, M = A->rmap->N, N = B->cmap->N;
905: PetscBool flg;
907: PetscFunctionBegin;
908: PetscCall(MatSetSizes(C, m, n, M, N));
909: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)C, &flg, MATSEQDENSE, MATMPIDENSE, ""));
910: if (!flg) {
911: PetscCall(MatSetType(C, ((PetscObject)B)->type_name));
912: C->ops->productsymbolic = MatProductSymbolic_SchurComplement_Dense;
913: }
914: PetscCall(MatSetUp(C));
915: C->ops->productnumeric = MatProductNumeric_SchurComplement_Dense;
916: PetscFunctionReturn(PETSC_SUCCESS);
917: }
919: static PetscErrorCode MatProductSetFromOptions_Dense_AB(Mat C)
920: {
921: PetscFunctionBegin;
922: C->ops->productsymbolic = MatProductSymbolic_SchurComplement_Dense;
923: PetscFunctionReturn(PETSC_SUCCESS);
924: }
926: static PetscErrorCode MatProductSetFromOptions_SchurComplement_Dense(Mat C)
927: {
928: Mat_Product *product = C->product;
930: PetscFunctionBegin;
931: PetscCheck(product->type == MATPRODUCT_AB, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Not for product type %s", MatProductTypes[product->type]);
932: PetscCall(MatProductSetFromOptions_Dense_AB(C));
933: PetscFunctionReturn(PETSC_SUCCESS);
934: }
936: /*MC
937: MATSCHURCOMPLEMENT - "schurcomplement" - Matrix type that behaves like the Schur complement of a matrix.
939: Level: intermediate
941: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatType`, `MatCreateSchurComplement()`, `MatSchurComplementComputeExplicitOperator()`,
942: `MatSchurComplementGetSubMatrices()`, `MatSchurComplementGetKSP()`
943: M*/
944: PETSC_EXTERN PetscErrorCode MatCreate_SchurComplement(Mat N)
945: {
946: Mat_SchurComplement *Na;
948: PetscFunctionBegin;
949: PetscCall(PetscNew(&Na));
950: N->data = (void *)Na;
952: N->ops->destroy = MatDestroy_SchurComplement;
953: N->ops->getvecs = MatCreateVecs_SchurComplement;
954: N->ops->view = MatView_SchurComplement;
955: N->ops->mult = MatMult_SchurComplement;
956: N->ops->multtranspose = MatMultTranspose_SchurComplement;
957: N->ops->multadd = MatMultAdd_SchurComplement;
958: N->ops->setfromoptions = MatSetFromOptions_SchurComplement;
959: N->assembled = PETSC_FALSE;
960: N->preallocated = PETSC_FALSE;
962: PetscCall(KSPCreate(PetscObjectComm((PetscObject)N), &Na->ksp));
963: PetscCall(PetscObjectChangeTypeName((PetscObject)N, MATSCHURCOMPLEMENT));
964: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_schurcomplement_seqdense_C", MatProductSetFromOptions_SchurComplement_Dense));
965: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_schurcomplement_mpidense_C", MatProductSetFromOptions_SchurComplement_Dense));
966: PetscFunctionReturn(PETSC_SUCCESS);
967: }