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: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_anytype_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 - matrix from which the preconditioner is constructed 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 - matrix from which the preconditioner is constructed 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: PetscCall(PetscObjectReference((PetscObject)A11));
250: Na->A = A00;
251: Na->Ap = Ap00;
252: Na->B = A01;
253: Na->C = A10;
254: Na->D = 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 - matrix from which the preconditioner is constructed 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: /*@
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 - matrix from which the preconditioner is constructed 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: Level: advanced
465: Notes:
466: This can be expensive when `S` is large, 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: `S` will automatically have the same prefix as `A` appended by `explicit_operator_`,
471: there are three options available: `-fieldsplit_1_explicit_operator_mat_type`,
472: `-fieldsplit_1_explicit_operator_mat_symmetric`, and `-fieldsplit_1_explicit_operator_mat_hermitian`
474: Developer Note:
475: The three aforementioned should not be parsed and used in this routine, but rather in `MatSetFromOptions()`
477: .seealso: [](ch_ksp), `MatCreateSchurComplement()`, `MatSchurComplementUpdateSubMatrices()`, `MatSchurComplementGetPmat()`
478: @*/
479: PetscErrorCode MatSchurComplementComputeExplicitOperator(Mat A, Mat *S)
480: {
481: Mat P, B, C, D, E = NULL, Bd, AinvBd, sub = NULL;
482: MatType mtype;
483: VecType vtype;
484: KSP ksp;
485: PetscInt n, N, m, M;
486: PetscBool flg = PETSC_FALSE, set, symm;
487: char prefix[256], type[256];
489: PetscFunctionBegin;
490: PetscAssertPointer(S, 2);
492: PetscCall(PetscObjectQuery((PetscObject)A, "AinvB", (PetscObject *)&AinvBd));
493: set = (PetscBool)(AinvBd != NULL);
494: if (set && AinvBd->cmap->N == -1) PetscFunctionReturn(PETSC_SUCCESS); // early bail out if composed Mat is uninitialized
495: PetscCall(MatSchurComplementGetSubMatrices(A, &P, NULL, &B, &C, &D));
496: PetscCall(MatGetVecType(B, &vtype));
497: PetscCall(MatGetLocalSize(B, &m, &n));
498: PetscCall(MatSchurComplementGetKSP(A, &ksp));
499: PetscCall(KSPSetUp(ksp));
500: if (set) {
501: PetscCheck(AinvBd->cmap->N >= A->cmap->N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_SIZ, "Composed Mat should have at least as many columns as the Schur complement (%" PetscInt_FMT " >= %" PetscInt_FMT ")", AinvBd->cmap->N, A->cmap->N);
502: PetscCall(MatGetType(AinvBd, &mtype));
503: if (AinvBd->cmap->N > A->cmap->N) {
504: Mat s[2];
506: PetscCall(MatDuplicate(AinvBd, MAT_DO_NOT_COPY_VALUES, &Bd));
507: PetscCall(MatDenseGetSubMatrix(Bd, PETSC_DECIDE, PETSC_DECIDE, A->cmap->N, AinvBd->cmap->N, s));
508: PetscCall(MatDenseGetSubMatrix(AinvBd, PETSC_DECIDE, PETSC_DECIDE, A->cmap->N, AinvBd->cmap->N, s + 1));
509: PetscCall(MatCopy(s[1], s[0], SAME_NONZERO_PATTERN)); // copy the last columns of the composed Mat, which are likely the input columns of PCApply_FieldSplit_Schur()
510: PetscCall(MatDenseRestoreSubMatrix(AinvBd, s + 1));
511: PetscCall(MatDenseRestoreSubMatrix(Bd, s));
512: PetscCall(MatDenseGetSubMatrix(Bd, PETSC_DECIDE, PETSC_DECIDE, 0, A->cmap->N, &sub));
513: PetscCall(MatConvert(B, mtype, MAT_REUSE_MATRIX, &sub)); // copy A01 into the first columns of the block of RHS of KSPMatSolve()
514: PetscCall(MatDenseRestoreSubMatrix(Bd, &sub));
515: } else PetscCall(MatConvert(B, mtype, MAT_INITIAL_MATRIX, &Bd));
516: } else {
517: PetscCall(MatGetSize(B, &M, &N));
518: PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)A), vtype, m, n, M, N, PETSC_DECIDE, NULL, &AinvBd));
519: PetscCall(MatGetType(AinvBd, &mtype));
520: PetscCall(MatConvert(B, mtype, MAT_INITIAL_MATRIX, &Bd));
521: }
522: PetscCall(KSPMatSolve(ksp, Bd, AinvBd));
523: if (set && AinvBd->cmap->N > A->cmap->N) {
524: Mat AinvB;
525: PetscScalar *v;
526: PetscMemType type;
528: PetscCall(MatDenseGetArrayWriteAndMemType(AinvBd, &v, &type)); // no easy way to resize a Mat, so create a new one with the same data pointer
529: PetscCall(MatCreateDenseWithMemType(PetscObjectComm((PetscObject)A), type, AinvBd->rmap->n, A->cmap->n, AinvBd->rmap->N, A->cmap->N, PETSC_DECIDE, v, &AinvB));
530: PetscCall(MatDenseReplaceArrayWithMemType(AinvB, type, v)); // let MatDestroy() free the data pointer
531: PetscCall(MatDenseRestoreArrayWriteAndMemType(AinvBd, &v));
532: PetscCall(MatHeaderReplace(AinvBd, &AinvB)); // replace the input composed Mat with just A00^-1 A01 (trailing columns are removed)
533: }
534: PetscCall(MatDestroy(&Bd));
535: if (!set) PetscCall(MatFilter(AinvBd, PETSC_SMALL, PETSC_FALSE, PETSC_FALSE));
536: if (D && !*S) {
537: PetscCall(MatGetLocalSize(D, &m, &n));
538: PetscCall(MatGetSize(D, &M, &N));
539: PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)A), vtype, m, n, M, N, PETSC_DECIDE, NULL, S));
540: } else if (*S) {
541: PetscCall(MatGetType(AinvBd, &mtype));
542: PetscCall(MatSetType(*S, mtype));
543: }
544: PetscCall(MatMatMult(C, AinvBd, *S ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX, PETSC_DETERMINE, S));
545: if (!set) PetscCall(MatDestroy(&AinvBd));
546: else {
547: PetscCall(MatScale(AinvBd, -1.0));
548: PetscCall(MatFilter(AinvBd, PETSC_MACHINE_EPSILON, PETSC_FALSE, PETSC_FALSE));
549: PetscCall(MatFilter(*S, PETSC_MACHINE_EPSILON, PETSC_FALSE, PETSC_FALSE));
550: }
551: if (D) {
552: PetscCall(PetscObjectTypeCompareAny((PetscObject)D, &flg, MATSEQSBAIJ, MATMPISBAIJ, ""));
553: if (flg) {
554: PetscCall(MatIsSymmetricKnown(A, &set, &symm));
555: if (!set || !symm) PetscCall(MatConvert(D, MATBAIJ, MAT_INITIAL_MATRIX, &E)); /* convert the (1,1) block to nonsymmetric storage for MatAXPY() */
556: }
557: 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 */
558: if (!E && flg) PetscCall(MatConvert(*S, MATSBAIJ, MAT_INPLACE_MATRIX, S)); /* if A is symmetric and the (1,1) block is a MatSBAIJ, return S as a MatSBAIJ since the lower triangular part is invalid */
559: }
560: PetscCall(MatDestroy(&E));
561: PetscCall(MatScale(*S, -1.0));
562: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%sexplicit_operator_", ((PetscObject)A)->prefix ? ((PetscObject)A)->prefix : ""));
563: PetscCall(MatSetOptionsPrefix(*S, prefix));
564: PetscObjectOptionsBegin((PetscObject)*S);
565: PetscCall(PetscOptionsFList("-mat_type", "Matrix type", "MatSetType", MatList, !E && flg ? MATSBAIJ : mtype, type, 256, &set));
566: if (set) PetscCall(MatConvert(*S, type, MAT_INPLACE_MATRIX, S));
567: flg = PETSC_FALSE;
568: PetscCall(PetscOptionsBool("-mat_symmetric", "Sets the MAT_SYMMETRIC option", "MatSetOption", flg, &flg, &set));
569: if (set) PetscCall(MatSetOption(*S, MAT_SYMMETRIC, flg));
570: if (PetscDefined(USE_COMPLEX)) {
571: flg = PETSC_FALSE;
572: PetscCall(PetscOptionsBool("-mat_hermitian", "Sets the MAT_HERMITIAN option", "MatSetOption", flg, &flg, &set));
573: if (set) PetscCall(MatSetOption(*S, MAT_HERMITIAN, flg));
574: }
575: PetscOptionsEnd();
576: PetscFunctionReturn(PETSC_SUCCESS);
577: }
579: /* Developer Notes:
580: This should be implemented with a MatCreate_SchurComplement() as that is the standard design for new Mat classes. */
581: PetscErrorCode MatGetSchurComplement_Basic(Mat mat, IS isrow0, IS iscol0, IS isrow1, IS iscol1, MatReuse mreuse, Mat *S, MatSchurComplementAinvType ainvtype, MatReuse preuse, Mat *Sp)
582: {
583: Mat A = NULL, Ap = NULL, B = NULL, C = NULL, D = NULL;
584: MatReuse reuse;
586: PetscFunctionBegin;
596: if (mreuse == MAT_IGNORE_MATRIX && preuse == MAT_IGNORE_MATRIX) PetscFunctionReturn(PETSC_SUCCESS);
600: PetscCheck(!mat->factortype, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
602: reuse = MAT_INITIAL_MATRIX;
603: if (mreuse == MAT_REUSE_MATRIX) {
604: PetscCall(MatSchurComplementGetSubMatrices(*S, &A, &Ap, &B, &C, &D));
605: PetscCheck(A && Ap && B && C, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Attempting to reuse matrix but Schur complement matrices unset");
606: PetscCheck(A == Ap, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Matrix for constructing the preconditioner does not match operator");
607: PetscCall(MatDestroy(&Ap)); /* get rid of extra reference */
608: reuse = MAT_REUSE_MATRIX;
609: }
610: PetscCall(MatCreateSubMatrix(mat, isrow0, iscol0, reuse, &A));
611: PetscCall(MatCreateSubMatrix(mat, isrow0, iscol1, reuse, &B));
612: PetscCall(MatCreateSubMatrix(mat, isrow1, iscol0, reuse, &C));
613: PetscCall(MatCreateSubMatrix(mat, isrow1, iscol1, reuse, &D));
614: switch (mreuse) {
615: case MAT_INITIAL_MATRIX:
616: PetscCall(MatCreateSchurComplement(A, A, B, C, D, S));
617: break;
618: case MAT_REUSE_MATRIX:
619: PetscCall(MatSchurComplementUpdateSubMatrices(*S, A, A, B, C, D));
620: break;
621: default:
622: PetscCheck(mreuse == MAT_IGNORE_MATRIX, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Unrecognized value of mreuse %d", (int)mreuse);
623: }
624: if (preuse != MAT_IGNORE_MATRIX) PetscCall(MatCreateSchurComplementPmat(A, B, C, D, ainvtype, preuse, Sp));
625: PetscCall(MatDestroy(&A));
626: PetscCall(MatDestroy(&B));
627: PetscCall(MatDestroy(&C));
628: PetscCall(MatDestroy(&D));
629: PetscFunctionReturn(PETSC_SUCCESS);
630: }
632: /*@
633: MatGetSchurComplement - Obtain the Schur complement from eliminating part of the matrix in another part.
635: Collective
637: Input Parameters:
638: + A - matrix in which the complement is to be taken
639: . isrow0 - rows to eliminate
640: . iscol0 - columns to eliminate, (isrow0,iscol0) should be square and nonsingular
641: . isrow1 - rows in which the Schur complement is formed
642: . iscol1 - columns in which the Schur complement is formed
643: . mreuse - `MAT_INITIAL_MATRIX` or `MAT_REUSE_MATRIX`, use `MAT_IGNORE_MATRIX` to put nothing in `S`
644: . ainvtype - the type of approximation used for the inverse of the (0,0) block used in forming `Sp`:
645: `MAT_SCHUR_COMPLEMENT_AINV_DIAG`, `MAT_SCHUR_COMPLEMENT_AINV_LUMP`, `MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG`, or `MAT_SCHUR_COMPLEMENT_AINV_FULL`
646: - preuse - `MAT_INITIAL_MATRIX` or `MAT_REUSE_MATRIX`, use `MAT_IGNORE_MATRIX` to put nothing in `Sp`
648: Output Parameters:
649: + S - exact Schur complement, often of type `MATSCHURCOMPLEMENT` which is difficult to use for preconditioning
650: - Sp - approximate Schur complement from which a preconditioner can be built $A11 - A10 inv(DIAGFORM(A00)) A01$
652: Level: advanced
654: Notes:
655: Since the real Schur complement is usually dense, providing a good approximation to `Sp` usually requires
656: application-specific information.
658: Sometimes users would like to provide problem-specific data in the Schur complement, usually only for special row
659: 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
660: "MatGetSchurComplement_C" to their function. If their function needs to fall back to the default implementation, it
661: should call `MatGetSchurComplement_Basic()`.
663: `MatCreateSchurComplement()` takes as arguments the four submatrices and returns the virtual Schur complement (what this function returns in S).
665: `MatSchurComplementGetPmat()` takes the virtual Schur complement and returns an explicit approximate Schur complement (what this returns in Sp).
667: In other words calling `MatCreateSchurComplement()` followed by `MatSchurComplementGetPmat()` produces the same output as this function but with slightly different
668: inputs. The actually submatrices of the original block matrix instead of index sets to the submatrices.
670: Developer Notes:
671: The API that includes `MatGetSchurComplement()`, `MatCreateSchurComplement()`, `MatSchurComplementGetPmat()` should be refactored to
672: remove redundancy and be clearer and simpler.
674: .seealso: [](ch_ksp), `MatCreateSubMatrix()`, `PCFIELDSPLIT`, `MatCreateSchurComplement()`, `MatSchurComplementAinvType`
675: @*/
676: PetscErrorCode MatGetSchurComplement(Mat A, IS isrow0, IS iscol0, IS isrow1, IS iscol1, MatReuse mreuse, Mat *S, MatSchurComplementAinvType ainvtype, MatReuse preuse, Mat *Sp)
677: {
678: PetscErrorCode (*f)(Mat, IS, IS, IS, IS, MatReuse, Mat *, MatReuse, Mat *) = NULL;
680: PetscFunctionBegin;
692: PetscCheck(!A->factortype, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
693: 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. */
694: PetscCall(PetscObjectQueryFunction((PetscObject)*S, "MatGetSchurComplement_C", &f));
695: }
696: if (f) PetscCall((*f)(A, isrow0, iscol0, isrow1, iscol1, mreuse, S, preuse, Sp));
697: else PetscCall(MatGetSchurComplement_Basic(A, isrow0, iscol0, isrow1, iscol1, mreuse, S, ainvtype, preuse, Sp));
698: PetscFunctionReturn(PETSC_SUCCESS);
699: }
701: /*@
702: MatSchurComplementSetAinvType - set the type of approximation used for the inverse of the (0,0) block used in forming `Sp` in `MatSchurComplementGetPmat()`
704: Not Collective
706: Input Parameters:
707: + S - matrix obtained with `MatCreateSchurComplement()` (or equivalent) and implementing the action of $A11 - A10 ksp(A00,Ap00) A01$
708: - ainvtype - type of approximation to be used to form approximate Schur complement $Sp = A11 - A10 inv(DIAGFORM(A00)) A01$:
709: `MAT_SCHUR_COMPLEMENT_AINV_DIAG`, `MAT_SCHUR_COMPLEMENT_AINV_LUMP`, `MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG`, or `MAT_SCHUR_COMPLEMENT_AINV_FULL`
711: Options Database Key:
712: . -mat_schur_complement_ainv_type diag | lump | blockdiag | full - set schur complement type
714: Level: advanced
716: .seealso: [](ch_ksp), `MatSchurComplementAinvType`, `MatCreateSchurComplement()`, `MatGetSchurComplement()`, `MatSchurComplementGetPmat()`, `MatSchurComplementGetAinvType()`
717: @*/
718: PetscErrorCode MatSchurComplementSetAinvType(Mat S, MatSchurComplementAinvType ainvtype)
719: {
720: PetscBool isschur;
721: Mat_SchurComplement *schur;
723: PetscFunctionBegin;
725: PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
726: if (!isschur) PetscFunctionReturn(PETSC_SUCCESS);
728: schur = (Mat_SchurComplement *)S->data;
729: 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);
730: schur->ainvtype = ainvtype;
731: PetscFunctionReturn(PETSC_SUCCESS);
732: }
734: /*@
735: MatSchurComplementGetAinvType - get the type of approximation for the inverse of the (0,0) block used in forming `Sp` in `MatSchurComplementGetPmat()`
737: Not Collective
739: Input Parameter:
740: . S - matrix obtained with `MatCreateSchurComplement()` (or equivalent) and implementing the action of $A11 - A10 ksp(A00,Ap00) A01$
742: Output Parameter:
743: . ainvtype - type of approximation used to form approximate Schur complement Sp = A11 - A10 inv(DIAGFORM(A00)) A01:
744: `MAT_SCHUR_COMPLEMENT_AINV_DIAG`, `MAT_SCHUR_COMPLEMENT_AINV_LUMP`, `MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG`, or `MAT_SCHUR_COMPLEMENT_AINV_FULL`
746: Level: advanced
748: .seealso: [](ch_ksp), `MatSchurComplementAinvType`, `MatCreateSchurComplement()`, `MatGetSchurComplement()`, `MatSchurComplementGetPmat()`, `MatSchurComplementSetAinvType()`
749: @*/
750: PetscErrorCode MatSchurComplementGetAinvType(Mat S, MatSchurComplementAinvType *ainvtype)
751: {
752: PetscBool isschur;
753: Mat_SchurComplement *schur;
755: PetscFunctionBegin;
757: PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
758: PetscCheck(isschur, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)S)->type_name);
759: schur = (Mat_SchurComplement *)S->data;
760: if (ainvtype) *ainvtype = schur->ainvtype;
761: PetscFunctionReturn(PETSC_SUCCESS);
762: }
764: /*@
765: MatCreateSchurComplementPmat - create a matrix for preconditioning the Schur complement by explicitly assembling the sparse matrix
766: $Sp = A11 - A10 inv(DIAGFORM(A00)) A01$
768: Collective
770: Input Parameters:
771: + A00 - the upper-left part of the original matrix $A = [A00 A01; A10 A11]$
772: . A01 - (optional) the upper-right part of the original matrix $A = [A00 A01; A10 A11]$
773: . A10 - (optional) the lower-left part of the original matrix $A = [A00 A01; A10 A11]$
774: . A11 - (optional) the lower-right part of the original matrix $A = [A00 A01; A10 A11]$
775: . ainvtype - type of approximation for DIAGFORM(A00) used when forming $Sp = A11 - A10 inv(DIAGFORM(A00)) A01$. See `MatSchurComplementAinvType`.
776: - 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`
778: Output Parameter:
779: . Sp - approximate Schur complement suitable for constructing a preconditioner for the true Schur complement $S = A11 - A10 inv(A00) A01$
781: Level: advanced
783: .seealso: [](ch_ksp), `MatCreateSchurComplement()`, `MatGetSchurComplement()`, `MatSchurComplementGetPmat()`, `MatSchurComplementAinvType`
784: @*/
785: PetscErrorCode MatCreateSchurComplementPmat(Mat A00, Mat A01, Mat A10, Mat A11, MatSchurComplementAinvType ainvtype, MatReuse preuse, Mat *Sp)
786: {
787: PetscInt N00;
789: PetscFunctionBegin;
790: /* 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. */
791: /* TODO: Perhaps should create an appropriately-sized zero matrix of the same type as A00? */
793: if (preuse == MAT_IGNORE_MATRIX) PetscFunctionReturn(PETSC_SUCCESS);
795: /* A zero size A00 or empty A01 or A10 imply S = A11. */
796: PetscCall(MatGetSize(A00, &N00, NULL));
797: if (!A01 || !A10 || !N00) {
798: if (preuse == MAT_INITIAL_MATRIX) {
799: PetscCall(MatDuplicate(A11, MAT_COPY_VALUES, Sp));
800: } else { /* MAT_REUSE_MATRIX */
801: /* TODO: when can we pass SAME_NONZERO_PATTERN? */
802: PetscCall(MatCopy(A11, *Sp, DIFFERENT_NONZERO_PATTERN));
803: }
804: } else {
805: Mat AdB, T;
806: Vec diag;
807: PetscBool flg;
809: if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP || ainvtype == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
810: PetscCall(PetscObjectTypeCompare((PetscObject)A01, MATTRANSPOSEVIRTUAL, &flg));
811: if (flg) {
812: PetscCall(MatTransposeGetMat(A01, &T));
813: PetscCall(MatTranspose(T, MAT_INITIAL_MATRIX, &AdB));
814: } else {
815: PetscCall(PetscObjectTypeCompare((PetscObject)A01, MATHERMITIANTRANSPOSEVIRTUAL, &flg));
816: if (flg) {
817: PetscCall(MatHermitianTransposeGetMat(A01, &T));
818: PetscCall(MatHermitianTranspose(T, MAT_INITIAL_MATRIX, &AdB));
819: }
820: }
821: if (!flg) PetscCall(MatDuplicate(A01, MAT_COPY_VALUES, &AdB));
822: else {
823: PetscScalar shift, scale;
825: 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));
826: PetscCall(MatShift(AdB, shift));
827: PetscCall(MatScale(AdB, scale));
828: }
829: PetscCall(MatCreateVecs(A00, &diag, NULL));
830: if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP) {
831: PetscCall(MatGetRowSum(A00, diag));
832: } else {
833: PetscCall(MatGetDiagonal(A00, diag));
834: }
835: PetscCall(VecReciprocal(diag));
836: PetscCall(MatDiagonalScale(AdB, diag, NULL));
837: PetscCall(VecDestroy(&diag));
838: } else if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG) {
839: Mat A00_inv;
840: MatType type;
841: MPI_Comm comm;
843: PetscCall(PetscObjectGetComm((PetscObject)A00, &comm));
844: PetscCall(MatGetType(A00, &type));
845: PetscCall(MatCreate(comm, &A00_inv));
846: PetscCall(MatSetType(A00_inv, type));
847: PetscCall(MatInvertBlockDiagonalMat(A00, A00_inv));
848: PetscCall(MatMatMult(A00_inv, A01, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &AdB));
849: PetscCall(MatDestroy(&A00_inv));
850: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown MatSchurComplementAinvType: %d", ainvtype);
851: /* Cannot really reuse Sp in MatMatMult() because of MatAYPX() -->
852: MatAXPY() --> MatHeaderReplace() --> MatDestroy_XXX_MatMatMult() */
853: if (preuse == MAT_REUSE_MATRIX) PetscCall(MatDestroy(Sp));
854: PetscCall(MatMatMult(A10, AdB, MAT_INITIAL_MATRIX, PETSC_DETERMINE, Sp));
855: PetscCall(MatScale(*Sp, -1.0));
856: if (A11) { /* TODO: when can we pass SAME_NONZERO_PATTERN? */
857: PetscCall(MatAXPY(*Sp, 1.0, A11, DIFFERENT_NONZERO_PATTERN));
858: }
859: PetscCall(MatDestroy(&AdB));
860: }
861: PetscFunctionReturn(PETSC_SUCCESS);
862: }
864: static PetscErrorCode MatSchurComplementGetPmat_Basic(Mat S, MatReuse preuse, Mat *Sp)
865: {
866: Mat A, B, C, D;
867: Mat_SchurComplement *schur = (Mat_SchurComplement *)S->data;
868: MatNullSpace sp;
870: PetscFunctionBegin;
871: if (preuse == MAT_IGNORE_MATRIX) PetscFunctionReturn(PETSC_SUCCESS);
872: PetscCall(MatSchurComplementGetSubMatrices(S, &A, NULL, &B, &C, &D));
873: PetscCheck(A, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONGSTATE, "Schur complement component matrices unset");
874: if (schur->ainvtype != MAT_SCHUR_COMPLEMENT_AINV_FULL) PetscCall(MatCreateSchurComplementPmat(A, B, C, D, schur->ainvtype, preuse, Sp));
875: else {
876: if (preuse == MAT_REUSE_MATRIX) PetscCall(MatDestroy(Sp));
877: PetscCall(MatSchurComplementComputeExplicitOperator(S, Sp));
878: }
879: /* If the Schur complement has a nullspace, then Sp nullspace contains it, independently of the ainv type */
880: PetscCall(MatGetNullSpace(S, &sp));
881: if (sp) PetscCall(MatSetNullSpace(*Sp, sp));
882: PetscCall(MatGetTransposeNullSpace(S, &sp));
883: if (sp) PetscCall(MatSetTransposeNullSpace(*Sp, sp));
884: PetscFunctionReturn(PETSC_SUCCESS);
885: }
887: /*@
888: MatSchurComplementGetPmat - Obtain a matrix for preconditioning the Schur complement by assembling $Sp = A11 - A10 inv(DIAGFORM(A00)) A01$
890: Collective
892: Input Parameters:
893: + S - matrix obtained with MatCreateSchurComplement() (or equivalent) that implements the action of $A11 - A10 ksp(A00,Ap00) A01$
894: - 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`
896: Output Parameter:
897: . Sp - approximate Schur complement suitable for preconditioning the exact Schur complement $S = A11 - A10 inv(A00) A01$
899: Level: advanced
901: Notes:
902: The approximation of `Sp` depends on the argument passed to `MatSchurComplementSetAinvType()`
903: `MAT_SCHUR_COMPLEMENT_AINV_DIAG`, `MAT_SCHUR_COMPLEMENT_AINV_LUMP`, `MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG`, or `MAT_SCHUR_COMPLEMENT_AINV_FULL`
904: -mat_schur_complement_ainv_type <diag,lump,blockdiag,full>
906: Sometimes users would like to provide problem-specific data in the Schur complement, usually only
907: for special row and column index sets. In that case, the user should call `PetscObjectComposeFunction()` to set
908: "MatSchurComplementGetPmat_C" to their function. If their function needs to fall back to the default implementation,
909: it should call `MatSchurComplementGetPmat_Basic()`.
911: Developer Notes:
912: The API that includes `MatGetSchurComplement()`, `MatCreateSchurComplement()`, `MatSchurComplementGetPmat()` should be refactored to
913: remove redundancy and be clearer and simpler.
915: This routine should be called `MatSchurComplementCreatePmat()`
917: .seealso: [](ch_ksp), `MatCreateSubMatrix()`, `PCFIELDSPLIT`, `MatGetSchurComplement()`, `MatCreateSchurComplement()`, `MatSchurComplementSetAinvType()`
918: @*/
919: PetscErrorCode MatSchurComplementGetPmat(Mat S, MatReuse preuse, Mat *Sp)
920: {
921: PetscErrorCode (*f)(Mat, MatReuse, Mat *);
923: PetscFunctionBegin;
927: if (preuse != MAT_IGNORE_MATRIX) {
928: PetscAssertPointer(Sp, 3);
929: if (preuse == MAT_INITIAL_MATRIX) *Sp = NULL;
931: }
932: PetscCheck(!S->factortype, PetscObjectComm((PetscObject)S), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
934: PetscCall(PetscObjectQueryFunction((PetscObject)S, "MatSchurComplementGetPmat_C", &f));
935: if (f) PetscCall((*f)(S, preuse, Sp));
936: else PetscCall(MatSchurComplementGetPmat_Basic(S, preuse, Sp));
937: PetscFunctionReturn(PETSC_SUCCESS);
938: }
940: static PetscErrorCode MatProductNumeric_SchurComplement_Dense(Mat C)
941: {
942: Mat_Product *product = C->product;
943: Mat_SchurComplement *Na = (Mat_SchurComplement *)product->A->data;
944: Mat work1, work2;
945: PetscScalar *v;
946: PetscInt lda;
948: PetscFunctionBegin;
949: PetscCall(MatMatMult(Na->B, product->B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &work1));
950: PetscCall(MatDuplicate(work1, MAT_DO_NOT_COPY_VALUES, &work2));
951: PetscCall(KSPMatSolve(Na->ksp, work1, work2));
952: PetscCall(MatDestroy(&work1));
953: PetscCall(MatDenseGetArrayWrite(C, &v));
954: PetscCall(MatDenseGetLDA(C, &lda));
955: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)C), C->rmap->n, C->cmap->n, C->rmap->N, C->cmap->N, v, &work1));
956: PetscCall(MatDenseSetLDA(work1, lda));
957: PetscCall(MatMatMult(Na->C, work2, MAT_REUSE_MATRIX, PETSC_DETERMINE, &work1));
958: PetscCall(MatDenseRestoreArrayWrite(C, &v));
959: PetscCall(MatDestroy(&work2));
960: PetscCall(MatDestroy(&work1));
961: if (Na->D) {
962: PetscCall(MatMatMult(Na->D, product->B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &work1));
963: PetscCall(MatAYPX(C, -1.0, work1, SAME_NONZERO_PATTERN));
964: PetscCall(MatDestroy(&work1));
965: } else PetscCall(MatScale(C, -1.0));
966: PetscFunctionReturn(PETSC_SUCCESS);
967: }
969: static PetscErrorCode MatProductSymbolic_SchurComplement_Dense(Mat C)
970: {
971: Mat_Product *product = C->product;
972: Mat A = product->A, B = product->B;
973: PetscInt m = A->rmap->n, n = B->cmap->n, M = A->rmap->N, N = B->cmap->N;
974: PetscBool flg;
976: PetscFunctionBegin;
977: PetscCall(MatSetSizes(C, m, n, M, N));
978: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)C, &flg, MATSEQDENSE, MATMPIDENSE, ""));
979: if (!flg) {
980: PetscCall(MatSetType(C, ((PetscObject)B)->type_name));
981: C->ops->productsymbolic = MatProductSymbolic_SchurComplement_Dense;
982: }
983: PetscCall(MatSetUp(C));
984: C->ops->productnumeric = MatProductNumeric_SchurComplement_Dense;
985: PetscFunctionReturn(PETSC_SUCCESS);
986: }
988: static PetscErrorCode MatProductSetFromOptions_SchurComplement_Dense(Mat C)
989: {
990: Mat_Product *product = C->product;
992: PetscFunctionBegin;
993: if (product->type != MATPRODUCT_AB) PetscFunctionReturn(PETSC_SUCCESS);
994: C->ops->productsymbolic = MatProductSymbolic_SchurComplement_Dense;
995: PetscFunctionReturn(PETSC_SUCCESS);
996: }
998: static PetscErrorCode MatProductNumeric_SchurComplement_Any(Mat C)
999: {
1000: Mat_SchurComplement *Na = (Mat_SchurComplement *)C->data;
1002: PetscFunctionBegin;
1003: if (Na->D && Na->D->product) PetscCall(MatProductNumeric(Na->D));
1004: if (Na->B->product) PetscCall(MatProductNumeric(Na->B));
1005: if (Na->C->product) PetscCall(MatProductNumeric(Na->C));
1006: C->assembled = PETSC_TRUE;
1007: PetscFunctionReturn(PETSC_SUCCESS);
1008: }
1010: static PetscErrorCode MatProductSymbolic_SchurComplement_Any(Mat C)
1011: {
1012: Mat_SchurComplement *Na = (Mat_SchurComplement *)C->data;
1014: PetscFunctionBegin;
1015: if (Na->D && Na->D->product) PetscCall(MatProductSymbolic(Na->D));
1016: if (Na->B->product) PetscCall(MatProductSymbolic(Na->B));
1017: if (Na->C->product) PetscCall(MatProductSymbolic(Na->C));
1018: C->ops->productnumeric = MatProductNumeric_SchurComplement_Any;
1019: C->preallocated = PETSC_TRUE;
1020: C->assembled = PETSC_FALSE;
1021: PetscFunctionReturn(PETSC_SUCCESS);
1022: }
1024: static PetscErrorCode MatProductSetFromOptions_SchurComplement_Any(Mat C)
1025: {
1026: Mat_Product *product = C->product;
1027: Mat_SchurComplement *Na, *Ca;
1028: Mat B = product->B, S = product->A, pB = NULL, pC = NULL, pD = NULL;
1029: KSP ksp;
1030: PetscInt m = PETSC_DECIDE, n = PETSC_DECIDE, M = PETSC_DECIDE, N = PETSC_DECIDE;
1031: MatProductType pbtype = MATPRODUCT_UNSPECIFIED, pctype = MATPRODUCT_UNSPECIFIED;
1032: PetscBool isschur;
1034: PetscFunctionBegin;
1035: if (product->type == MATPRODUCT_ABC || product->type == MATPRODUCT_AtB) PetscFunctionReturn(PETSC_SUCCESS);
1036: /* A * S not yet supported (should be easy though) */
1037: PetscCall(PetscObjectTypeCompare((PetscObject)S, MATSCHURCOMPLEMENT, &isschur));
1038: if (!isschur) PetscFunctionReturn(PETSC_SUCCESS);
1040: Na = (Mat_SchurComplement *)S->data;
1041: if (Na->D) {
1042: PetscCall(MatProductCreate(Na->D, B, NULL, &pD));
1043: PetscCall(MatProductSetType(pD, product->type));
1044: PetscCall(MatProductSetFromOptions(pD));
1045: }
1046: if (pD && !pD->ops->productsymbolic) {
1047: PetscCall(MatDestroy(&pD));
1048: PetscFunctionReturn(PETSC_SUCCESS);
1049: }
1051: /* S = A11 - A10 M A01 */
1052: switch (product->type) {
1053: case MATPRODUCT_AB: /* A11 B - A10 * M * A01 * B */
1054: pbtype = product->type;
1055: PetscCall(PetscObjectReference((PetscObject)Na->C));
1056: pC = Na->C;
1057: m = S->rmap->n;
1058: M = S->rmap->N;
1059: n = B->cmap->n;
1060: N = B->cmap->N;
1061: break;
1062: case MATPRODUCT_ABt: /* A11 B^t - A10 * M * A01 * B^t */
1063: pbtype = product->type;
1064: PetscCall(PetscObjectReference((PetscObject)Na->C));
1065: pC = Na->C;
1066: m = S->rmap->n;
1067: M = S->rmap->N;
1068: n = B->rmap->n;
1069: N = B->rmap->N;
1070: break;
1071: case MATPRODUCT_PtAP: /* Pt A11 P - Pt * A10 * M * A01 * P */
1072: pbtype = MATPRODUCT_AB;
1073: pctype = MATPRODUCT_AtB;
1074: m = B->cmap->n;
1075: M = B->cmap->N;
1076: n = B->cmap->n;
1077: N = B->cmap->N;
1078: break;
1079: case MATPRODUCT_RARt: /* R A11 Rt - R * A10 * M * A01 * Rt */
1080: pbtype = MATPRODUCT_ABt;
1081: pctype = MATPRODUCT_AB;
1082: m = B->rmap->n;
1083: M = B->rmap->N;
1084: n = B->rmap->n;
1085: N = B->rmap->N;
1086: break;
1087: default:
1088: break;
1089: }
1090: PetscCall(MatProductCreate(Na->B, B, NULL, &pB));
1091: PetscCall(MatProductSetType(pB, pbtype));
1092: PetscCall(MatProductSetFromOptions(pB));
1093: if (!pB->ops->productsymbolic) {
1094: PetscCall(MatDestroy(&pB));
1095: PetscFunctionReturn(PETSC_SUCCESS);
1096: }
1097: if (pC == NULL) { /* Some work can in principle be saved here if we recognize symmetry */
1098: PetscCall(MatProductCreate(B, Na->C, NULL, &pC));
1099: PetscCall(MatProductSetType(pC, pctype));
1100: PetscCall(MatProductSetFromOptions(pC));
1101: if (!pC->ops->productsymbolic) {
1102: PetscCall(MatDestroy(&pC));
1103: PetscFunctionReturn(PETSC_SUCCESS);
1104: }
1105: }
1106: PetscCall(MatSetType(C, MATSCHURCOMPLEMENT));
1107: PetscCall(MatSetSizes(C, m, n, M, N));
1108: PetscCall(PetscLayoutSetUp(C->rmap));
1109: PetscCall(PetscLayoutSetUp(C->cmap));
1110: PetscCall(PetscObjectReference((PetscObject)Na->A));
1111: PetscCall(PetscObjectReference((PetscObject)Na->Ap));
1112: Ca = (Mat_SchurComplement *)C->data;
1113: Ca->A = Na->A;
1114: Ca->Ap = Na->Ap;
1115: Ca->B = pB;
1116: Ca->C = pC;
1117: Ca->D = pD;
1118: C->ops->productsymbolic = MatProductSymbolic_SchurComplement_Any;
1119: PetscCall(MatSchurComplementGetKSP(S, &ksp));
1120: PetscCall(MatSchurComplementSetKSP(C, ksp));
1121: PetscFunctionReturn(PETSC_SUCCESS);
1122: }
1124: /*MC
1125: MATSCHURCOMPLEMENT - "schurcomplement" - Matrix type that behaves like the Schur complement of a matrix.
1127: Level: intermediate
1129: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatType`, `MatCreateSchurComplement()`, `MatSchurComplementComputeExplicitOperator()`,
1130: `MatSchurComplementGetSubMatrices()`, `MatSchurComplementGetKSP()`
1131: M*/
1132: PETSC_EXTERN PetscErrorCode MatCreate_SchurComplement(Mat N)
1133: {
1134: Mat_SchurComplement *Na;
1136: PetscFunctionBegin;
1137: PetscCall(PetscNew(&Na));
1138: N->data = (void *)Na;
1140: N->ops->destroy = MatDestroy_SchurComplement;
1141: N->ops->getvecs = MatCreateVecs_SchurComplement;
1142: N->ops->view = MatView_SchurComplement;
1143: N->ops->mult = MatMult_SchurComplement;
1144: N->ops->multtranspose = MatMultTranspose_SchurComplement;
1145: N->ops->multadd = MatMultAdd_SchurComplement;
1146: N->ops->setfromoptions = MatSetFromOptions_SchurComplement;
1147: N->assembled = PETSC_FALSE;
1148: N->preallocated = PETSC_FALSE;
1150: PetscCall(KSPCreate(PetscObjectComm((PetscObject)N), &Na->ksp));
1151: PetscCall(PetscObjectChangeTypeName((PetscObject)N, MATSCHURCOMPLEMENT));
1152: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_schurcomplement_seqdense_C", MatProductSetFromOptions_SchurComplement_Dense));
1153: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_schurcomplement_mpidense_C", MatProductSetFromOptions_SchurComplement_Dense));
1154: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_SchurComplement_Any));
1155: PetscFunctionReturn(PETSC_SUCCESS);
1156: }