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: }