Actual source code: fieldsplit.c
1: #include <petsc/private/pcimpl.h>
2: #include <petsc/private/kspimpl.h>
3: #include <petsc/private/matimpl.h>
4: #include <petscdm.h>
5: #include <petscdevice.h>
6: #if PetscDefined(HAVE_CUDA)
7: #include <petscdevice_cuda.h>
8: #endif
9: #if PetscDefined(HAVE_HIP)
10: #include <petscdevice_hip.h>
11: #endif
13: const char *const PCFieldSplitSchurPreTypes[] = {"SELF", "SELFP", "A11", "USER", "FULL", "PCFieldSplitSchurPreType", "PC_FIELDSPLIT_SCHUR_PRE_", NULL};
14: const char *const PCFieldSplitSchurFactTypes[] = {"DIAG", "LOWER", "UPPER", "FULL", "PCFieldSplitSchurFactType", "PC_FIELDSPLIT_SCHUR_FACT_", NULL};
16: PetscLogEvent KSP_Solve_FS_0, KSP_Solve_FS_1, KSP_Solve_FS_S, KSP_Solve_FS_U, KSP_Solve_FS_L, KSP_Solve_FS_2, KSP_Solve_FS_3, KSP_Solve_FS_4;
18: typedef struct _PC_FieldSplitLink *PC_FieldSplitLink;
19: struct _PC_FieldSplitLink {
20: KSP ksp;
21: Vec x, y, z;
22: Mat X, Y, Z;
23: char *splitname;
24: PetscInt nfields;
25: PetscInt *fields, *fields_col;
26: VecScatter sctx;
27: IS is, is_col;
28: PC_FieldSplitLink next, previous;
29: PetscLogEvent event;
31: /* Used only when setting coordinates with PCSetCoordinates */
32: PetscInt dim;
33: PetscInt ndofs;
34: PetscReal *coords;
35: };
37: typedef struct {
38: PCCompositeType type;
39: PetscBool defaultsplit; /* Flag for a system with a set of 'k' scalar fields with the same layout (and bs = k) */
40: PetscBool splitdefined; /* Flag is set after the splits have been defined, to prevent more splits from being added */
41: PetscInt bs; /* Block size for IS and Mat structures */
42: PetscInt nsplits; /* Number of field divisions defined */
43: Vec *x, *y, w1, w2;
44: Mat *mat; /* The diagonal block for each split */
45: Mat *pmat; /* The preconditioning diagonal block for each split */
46: Mat *Afield; /* The rows of the matrix associated with each split */
47: PetscBool issetup;
49: /* Only used when Schur complement preconditioning is used */
50: Mat B; /* The (0,1) block */
51: Mat C; /* The (1,0) block */
52: Mat schur; /* The Schur complement S = A11 - A10 A00^{-1} A01, the KSP here, kspinner, is H_1 in [El08] */
53: Mat schurp; /* Assembled approximation to S built by MatSchurComplement to be used as a matrix for constructing the preconditioner when solving with S */
54: Mat schur_user; /* User-provided matrix for constructing the preconditioner for the Schur complement */
55: PCFieldSplitSchurPreType schurpre; /* Determines which matrix is used for the Schur complement */
56: PCFieldSplitSchurFactType schurfactorization;
57: KSP kspschur; /* The solver for S */
58: KSP kspupper; /* The solver for A in the upper diagonal part of the factorization (H_2 in [El08]) */
59: PetscScalar schurscale; /* Scaling factor for the Schur complement solution with DIAG factorization */
61: /* Only used when Golub-Kahan bidiagonalization preconditioning is used */
62: Mat H; /* The modified matrix H = A00 + nu*A01*A01' */
63: PetscReal gkbtol; /* Stopping tolerance for lower bound estimate */
64: PetscInt gkbdelay; /* The delay window for the stopping criterion */
65: PetscReal gkbnu; /* Parameter for augmented Lagrangian H = A + nu*A01*A01' */
66: PetscInt gkbmaxit; /* Maximum number of iterations for outer loop */
67: PetscBool gkbmonitor; /* Monitor for gkb iterations and the lower bound error */
68: PetscViewer gkbviewer; /* Viewer context for gkbmonitor */
69: Vec u, v, d, Hu; /* Work vectors for the GKB algorithm */
70: PetscScalar *vecz; /* Contains intermediate values, eg for lower bound */
72: PC_FieldSplitLink head;
73: PetscBool isrestrict; /* indicates PCFieldSplitRestrictIS() has been last called on this object, hack */
74: PetscBool suboptionsset; /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */
75: PetscBool dm_splits; /* Whether to use DMCreateFieldDecomposition() whenever possible */
76: PetscBool diag_use_amat; /* Whether to extract diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
77: PetscBool offdiag_use_amat; /* Whether to extract off-diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
78: PetscBool detect; /* Whether to form 2-way split by finding zero diagonal entries */
79: PetscBool coordinates_set; /* Whether PCSetCoordinates has been called */
80: } PC_FieldSplit;
82: /*
83: Note:
84: there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of
85: inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the
86: PC you could change this.
87: */
89: /* This helper is so that setting a user-provided matrix is orthogonal to choosing to use it. This way the
90: * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */
91: static Mat FieldSplitSchurPre(PC_FieldSplit *jac)
92: {
93: switch (jac->schurpre) {
94: case PC_FIELDSPLIT_SCHUR_PRE_SELF:
95: return jac->schur;
96: case PC_FIELDSPLIT_SCHUR_PRE_SELFP:
97: return jac->schurp;
98: case PC_FIELDSPLIT_SCHUR_PRE_A11:
99: return jac->pmat[1];
100: case PC_FIELDSPLIT_SCHUR_PRE_FULL: /* We calculate this and store it in schur_user */
101: case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */
102: default:
103: return jac->schur_user ? jac->schur_user : jac->pmat[1];
104: }
105: }
107: #include <petscdraw.h>
108: static PetscErrorCode PCView_FieldSplit(PC pc, PetscViewer viewer)
109: {
110: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
111: PetscBool isascii, isdraw;
112: PetscInt i, j;
113: PC_FieldSplitLink ilink = jac->head;
115: PetscFunctionBegin;
116: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
117: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
118: if (isascii) {
119: if (jac->bs > 0) {
120: PetscCall(PetscViewerASCIIPrintf(viewer, " FieldSplit with %s composition: total splits = %" PetscInt_FMT ", blocksize = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits, jac->bs));
121: } else {
122: PetscCall(PetscViewerASCIIPrintf(viewer, " FieldSplit with %s composition: total splits = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits));
123: }
124: if (pc->useAmat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for blocks\n"));
125: if (jac->diag_use_amat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for diagonal blocks\n"));
126: if (jac->offdiag_use_amat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for off-diagonal blocks\n"));
127: PetscCall(PetscViewerASCIIPrintf(viewer, " Solver info for each split is in the following KSP objects:\n"));
128: for (i = 0; i < jac->nsplits; i++) {
129: if (ilink->fields) {
130: PetscCall(PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Fields ", i));
131: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
132: for (j = 0; j < ilink->nfields; j++) {
133: if (j > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ","));
134: PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, ilink->fields[j]));
135: }
136: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
137: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
138: } else {
139: PetscCall(PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Defined by IS\n", i));
140: }
141: PetscCall(KSPView(ilink->ksp, viewer));
142: ilink = ilink->next;
143: }
144: }
146: if (isdraw) {
147: PetscDraw draw;
148: PetscReal x, y, w, wd;
150: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
151: PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
152: w = 2 * PetscMin(1.0 - x, x);
153: wd = w / (jac->nsplits + 1);
154: x = x - wd * (jac->nsplits - 1) / 2.0;
155: for (i = 0; i < jac->nsplits; i++) {
156: PetscCall(PetscDrawPushCurrentPoint(draw, x, y));
157: PetscCall(KSPView(ilink->ksp, viewer));
158: PetscCall(PetscDrawPopCurrentPoint(draw));
159: x += wd;
160: ilink = ilink->next;
161: }
162: }
163: PetscFunctionReturn(PETSC_SUCCESS);
164: }
166: static PetscErrorCode PCView_FieldSplit_Schur(PC pc, PetscViewer viewer)
167: {
168: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
169: PetscBool isascii, isdraw;
170: PetscInt i, j;
171: PC_FieldSplitLink ilink = jac->head;
172: MatSchurComplementAinvType atype;
174: PetscFunctionBegin;
175: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
176: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
177: if (isascii) {
178: if (jac->bs > 0) {
179: PetscCall(PetscViewerASCIIPrintf(viewer, " FieldSplit with Schur preconditioner, blocksize = %" PetscInt_FMT ", factorization %s\n", jac->bs, PCFieldSplitSchurFactTypes[jac->schurfactorization]));
180: } else {
181: PetscCall(PetscViewerASCIIPrintf(viewer, " FieldSplit with Schur preconditioner, factorization %s\n", PCFieldSplitSchurFactTypes[jac->schurfactorization]));
182: }
183: if (pc->useAmat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for blocks\n"));
184: switch (jac->schurpre) {
185: case PC_FIELDSPLIT_SCHUR_PRE_SELF:
186: PetscCall(PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from S itself\n"));
187: break;
188: case PC_FIELDSPLIT_SCHUR_PRE_SELFP:
189: if (jac->schur) {
190: PetscCall(MatSchurComplementGetAinvType(jac->schur, &atype));
191: PetscCall(PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from Sp, an assembled approximation to S, which uses A00's %sinverse\n", atype == MAT_SCHUR_COMPLEMENT_AINV_DIAG ? "diagonal's " : (atype == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG ? "block diagonal's " : (atype == MAT_SCHUR_COMPLEMENT_AINV_FULL ? "full " : "lumped diagonal's "))));
192: }
193: break;
194: case PC_FIELDSPLIT_SCHUR_PRE_A11:
195: PetscCall(PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from A11\n"));
196: break;
197: case PC_FIELDSPLIT_SCHUR_PRE_FULL:
198: PetscCall(PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from the exact Schur complement\n"));
199: break;
200: case PC_FIELDSPLIT_SCHUR_PRE_USER:
201: if (jac->schur_user) {
202: PetscCall(PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from user provided matrix\n"));
203: } else {
204: PetscCall(PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from A11\n"));
205: }
206: break;
207: default:
208: SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre);
209: }
210: PetscCall(PetscViewerASCIIPrintf(viewer, " Split info:\n"));
211: PetscCall(PetscViewerASCIIPushTab(viewer));
212: for (i = 0; i < jac->nsplits; i++) {
213: if (ilink->fields) {
214: PetscCall(PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Fields ", i));
215: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
216: for (j = 0; j < ilink->nfields; j++) {
217: if (j > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ","));
218: PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, ilink->fields[j]));
219: }
220: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
221: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
222: } else {
223: PetscCall(PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Defined by IS\n", i));
224: }
225: ilink = ilink->next;
226: }
227: PetscCall(PetscViewerASCIIPrintf(viewer, "KSP solver for A00 block\n"));
228: PetscCall(PetscViewerASCIIPushTab(viewer));
229: if (jac->head) PetscCall(KSPView(jac->head->ksp, viewer));
230: else PetscCall(PetscViewerASCIIPrintf(viewer, " not yet available\n"));
231: PetscCall(PetscViewerASCIIPopTab(viewer));
232: if (jac->head && jac->kspupper != jac->head->ksp) {
233: PetscCall(PetscViewerASCIIPrintf(viewer, "KSP solver for upper A00 in upper triangular factor\n"));
234: PetscCall(PetscViewerASCIIPushTab(viewer));
235: if (jac->kspupper) PetscCall(KSPView(jac->kspupper, viewer));
236: else PetscCall(PetscViewerASCIIPrintf(viewer, " not yet available\n"));
237: PetscCall(PetscViewerASCIIPopTab(viewer));
238: }
239: PetscCall(PetscViewerASCIIPrintf(viewer, "KSP solver for S = A11 - A10 inv(A00) A01\n"));
240: PetscCall(PetscViewerASCIIPushTab(viewer));
241: if (jac->kspschur) {
242: PetscCall(KSPView(jac->kspschur, viewer));
243: } else {
244: PetscCall(PetscViewerASCIIPrintf(viewer, " not yet available\n"));
245: }
246: PetscCall(PetscViewerASCIIPopTab(viewer));
247: PetscCall(PetscViewerASCIIPopTab(viewer));
248: } else if (isdraw && jac->head) {
249: PetscDraw draw;
250: PetscReal x, y, w, wd, h;
251: PetscInt cnt = 2;
252: char str[32];
254: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
255: PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
256: if (jac->kspupper != jac->head->ksp) cnt++;
257: w = 2 * PetscMin(1.0 - x, x);
258: wd = w / (cnt + 1);
260: PetscCall(PetscSNPrintf(str, 32, "Schur fact. %s", PCFieldSplitSchurFactTypes[jac->schurfactorization]));
261: PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h));
262: y -= h;
263: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_USER && !jac->schur_user) {
264: PetscCall(PetscSNPrintf(str, 32, "Prec. for Schur from %s", PCFieldSplitSchurPreTypes[PC_FIELDSPLIT_SCHUR_PRE_A11]));
265: } else {
266: PetscCall(PetscSNPrintf(str, 32, "Prec. for Schur from %s", PCFieldSplitSchurPreTypes[jac->schurpre]));
267: }
268: PetscCall(PetscDrawStringBoxed(draw, x + wd * (cnt - 1) / 2.0, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h));
269: y -= h;
270: x = x - wd * (cnt - 1) / 2.0;
272: PetscCall(PetscDrawPushCurrentPoint(draw, x, y));
273: PetscCall(KSPView(jac->head->ksp, viewer));
274: PetscCall(PetscDrawPopCurrentPoint(draw));
275: if (jac->kspupper != jac->head->ksp) {
276: x += wd;
277: PetscCall(PetscDrawPushCurrentPoint(draw, x, y));
278: PetscCall(KSPView(jac->kspupper, viewer));
279: PetscCall(PetscDrawPopCurrentPoint(draw));
280: }
281: x += wd;
282: PetscCall(PetscDrawPushCurrentPoint(draw, x, y));
283: PetscCall(KSPView(jac->kspschur, viewer));
284: PetscCall(PetscDrawPopCurrentPoint(draw));
285: }
286: PetscFunctionReturn(PETSC_SUCCESS);
287: }
289: static PetscErrorCode PCView_FieldSplit_GKB(PC pc, PetscViewer viewer)
290: {
291: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
292: PetscBool isascii, isdraw;
293: PetscInt i, j;
294: PC_FieldSplitLink ilink = jac->head;
296: PetscFunctionBegin;
297: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
298: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
299: if (isascii) {
300: if (jac->bs > 0) {
301: PetscCall(PetscViewerASCIIPrintf(viewer, " FieldSplit with %s composition: total splits = %" PetscInt_FMT ", blocksize = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits, jac->bs));
302: } else {
303: PetscCall(PetscViewerASCIIPrintf(viewer, " FieldSplit with %s composition: total splits = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits));
304: }
305: if (pc->useAmat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for blocks\n"));
306: if (jac->diag_use_amat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for diagonal blocks\n"));
307: if (jac->offdiag_use_amat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for off-diagonal blocks\n"));
309: PetscCall(PetscViewerASCIIPrintf(viewer, " Stopping tolerance=%.1e, delay in error estimate=%" PetscInt_FMT ", maximum iterations=%" PetscInt_FMT "\n", (double)jac->gkbtol, jac->gkbdelay, jac->gkbmaxit));
310: PetscCall(PetscViewerASCIIPrintf(viewer, " Solver info for H = A00 + nu*A01*A01' matrix:\n"));
311: PetscCall(PetscViewerASCIIPushTab(viewer));
313: if (ilink->fields) {
314: PetscCall(PetscViewerASCIIPrintf(viewer, "Split number 0 Fields "));
315: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
316: for (j = 0; j < ilink->nfields; j++) {
317: if (j > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ","));
318: PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, ilink->fields[j]));
319: }
320: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
321: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
322: } else {
323: PetscCall(PetscViewerASCIIPrintf(viewer, "Split number 0 Defined by IS\n"));
324: }
325: PetscCall(KSPView(ilink->ksp, viewer));
327: PetscCall(PetscViewerASCIIPopTab(viewer));
328: }
330: if (isdraw) {
331: PetscDraw draw;
332: PetscReal x, y, w, wd;
334: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
335: PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
336: w = 2 * PetscMin(1.0 - x, x);
337: wd = w / (jac->nsplits + 1);
338: x = x - wd * (jac->nsplits - 1) / 2.0;
339: for (i = 0; i < jac->nsplits; i++) {
340: PetscCall(PetscDrawPushCurrentPoint(draw, x, y));
341: PetscCall(KSPView(ilink->ksp, viewer));
342: PetscCall(PetscDrawPopCurrentPoint(draw));
343: x += wd;
344: ilink = ilink->next;
345: }
346: }
347: PetscFunctionReturn(PETSC_SUCCESS);
348: }
350: /* Precondition: jac->bs is set to a meaningful value or MATNEST */
351: static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc)
352: {
353: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
354: PetscInt bs, i, nfields, *ifields, nfields_col, *ifields_col;
355: PetscBool flg, flg_col, mnest;
356: char optionname[128], splitname[8], optionname_col[128];
358: PetscFunctionBegin;
359: PetscCall(PetscObjectTypeCompare((PetscObject)pc->mat, MATNEST, &mnest));
360: if (mnest) {
361: PetscCall(MatNestGetSize(pc->pmat, &bs, NULL));
362: } else {
363: bs = jac->bs;
364: }
365: PetscCall(PetscMalloc2(bs, &ifields, bs, &ifields_col));
366: for (i = 0, flg = PETSC_TRUE;; i++) {
367: PetscCall(PetscSNPrintf(splitname, sizeof(splitname), "%" PetscInt_FMT, i));
368: PetscCall(PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%" PetscInt_FMT "_fields", i));
369: PetscCall(PetscSNPrintf(optionname_col, sizeof(optionname_col), "-pc_fieldsplit_%" PetscInt_FMT "_fields_col", i));
370: nfields = bs;
371: nfields_col = bs;
372: PetscCall(PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg));
373: PetscCall(PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, optionname_col, ifields_col, &nfields_col, &flg_col));
374: if (!flg) break;
375: else if (flg && !flg_col) {
376: PetscCheck(nfields, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot list zero fields");
377: PetscCall(PCFieldSplitSetFields(pc, splitname, nfields, ifields, ifields));
378: } else {
379: PetscCheck(nfields && nfields_col, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot list zero fields");
380: PetscCheck(nfields == nfields_col, PETSC_COMM_SELF, PETSC_ERR_USER, "Number of row and column fields must match");
381: PetscCall(PCFieldSplitSetFields(pc, splitname, nfields, ifields, ifields_col));
382: }
383: }
384: if (i > 0) {
385: /* Makes command-line setting of splits take precedence over setting them in code.
386: Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would
387: create new splits, which would probably not be what the user wanted. */
388: jac->splitdefined = PETSC_TRUE;
389: }
390: PetscCall(PetscFree2(ifields, ifields_col));
391: PetscFunctionReturn(PETSC_SUCCESS);
392: }
394: static PetscErrorCode PCFieldSplitSetDefaults(PC pc)
395: {
396: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
397: PC_FieldSplitLink ilink = jac->head;
398: PetscBool fieldsplit_default = PETSC_FALSE, coupling = PETSC_FALSE;
399: PetscInt i;
401: PetscFunctionBegin;
402: /*
403: Kinda messy, but at least this now uses DMCreateFieldDecomposition().
404: Should probably be rewritten.
405: */
406: if (!ilink) {
407: PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_fieldsplit_detect_coupling", &coupling, NULL));
408: if (pc->dm && jac->dm_splits && !jac->detect && !coupling) {
409: PetscInt numFields, f, i, j;
410: char **fieldNames;
411: IS *fields;
412: DM *dms;
413: DM subdm[128];
414: PetscBool flg;
416: PetscCall(DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms));
417: /* Allow the user to prescribe the splits */
418: for (i = 0, flg = PETSC_TRUE;; i++) {
419: PetscInt ifields[128];
420: IS compField;
421: char optionname[128], splitname[8];
422: PetscInt nfields = numFields;
424: PetscCall(PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%" PetscInt_FMT "_fields", i));
425: PetscCall(PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg));
426: if (!flg) break;
427: PetscCheck(numFields <= 128, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot currently support %" PetscInt_FMT " > 128 fields", numFields);
428: PetscCall(DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i]));
429: if (nfields == 1) {
430: PetscCall(PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField));
431: } else {
432: PetscCall(PetscSNPrintf(splitname, sizeof(splitname), "%" PetscInt_FMT, i));
433: PetscCall(PCFieldSplitSetIS(pc, splitname, compField));
434: }
435: PetscCall(ISDestroy(&compField));
436: for (j = 0; j < nfields; ++j) {
437: f = ifields[j];
438: PetscCall(PetscFree(fieldNames[f]));
439: PetscCall(ISDestroy(&fields[f]));
440: }
441: }
442: if (i == 0) {
443: for (f = 0; f < numFields; ++f) {
444: PetscCall(PCFieldSplitSetIS(pc, fieldNames[f], fields[f]));
445: PetscCall(PetscFree(fieldNames[f]));
446: PetscCall(ISDestroy(&fields[f]));
447: }
448: } else {
449: for (j = 0; j < numFields; j++) PetscCall(DMDestroy(dms + j));
450: PetscCall(PetscFree(dms));
451: PetscCall(PetscMalloc1(i, &dms));
452: for (j = 0; j < i; ++j) dms[j] = subdm[j];
453: }
454: PetscCall(PetscFree(fieldNames));
455: PetscCall(PetscFree(fields));
456: if (dms) {
457: PetscCall(PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n"));
458: for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) {
459: const char *prefix;
460: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ilink->ksp, &prefix));
461: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)dms[i], prefix));
462: PetscCall(KSPSetDM(ilink->ksp, dms[i]));
463: PetscCall(KSPSetDMActive(ilink->ksp, KSP_DMACTIVE_ALL, PETSC_FALSE));
464: PetscCall(PetscObjectIncrementTabLevel((PetscObject)dms[i], (PetscObject)ilink->ksp, 0));
465: PetscCall(DMDestroy(&dms[i]));
466: }
467: PetscCall(PetscFree(dms));
468: }
469: } else {
470: if (jac->bs <= 0) {
471: if (pc->pmat) PetscCall(MatGetBlockSize(pc->pmat, &jac->bs));
472: else jac->bs = 1;
473: }
475: if (jac->detect) {
476: IS zerodiags, rest;
477: PetscInt nmin, nmax;
479: PetscCall(MatGetOwnershipRange(pc->mat, &nmin, &nmax));
480: if (jac->diag_use_amat) {
481: PetscCall(MatFindZeroDiagonals(pc->mat, &zerodiags));
482: } else {
483: PetscCall(MatFindZeroDiagonals(pc->pmat, &zerodiags));
484: }
485: PetscCall(ISComplement(zerodiags, nmin, nmax, &rest));
486: PetscCall(PCFieldSplitSetIS(pc, "0", rest));
487: PetscCall(PCFieldSplitSetIS(pc, "1", zerodiags));
488: PetscCall(ISDestroy(&zerodiags));
489: PetscCall(ISDestroy(&rest));
490: } else if (coupling) {
491: IS coupling, rest;
492: PetscInt nmin, nmax;
494: PetscCall(MatGetOwnershipRange(pc->mat, &nmin, &nmax));
495: if (jac->offdiag_use_amat) {
496: PetscCall(MatFindOffBlockDiagonalEntries(pc->mat, &coupling));
497: } else {
498: PetscCall(MatFindOffBlockDiagonalEntries(pc->pmat, &coupling));
499: }
500: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc->mat), nmax - nmin, nmin, 1, &rest));
501: PetscCall(ISSetIdentity(rest));
502: PetscCall(PCFieldSplitSetIS(pc, "0", rest));
503: PetscCall(PCFieldSplitSetIS(pc, "1", coupling));
504: PetscCall(ISDestroy(&coupling));
505: PetscCall(ISDestroy(&rest));
506: } else {
507: PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_fieldsplit_default", &fieldsplit_default, NULL));
508: if (!fieldsplit_default) {
509: /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit()
510: then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */
511: PetscCall(PCFieldSplitSetRuntimeSplits_Private(pc));
512: if (jac->splitdefined) PetscCall(PetscInfo(pc, "Splits defined using the options database\n"));
513: }
514: if ((fieldsplit_default || !jac->splitdefined) && !jac->isrestrict) {
515: Mat M = pc->pmat;
516: PetscBool isnest;
517: PetscInt nf;
519: PetscCall(PetscInfo(pc, "Using default splitting of fields\n"));
520: PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATNEST, &isnest));
521: if (!isnest) {
522: M = pc->mat;
523: PetscCall(PetscObjectTypeCompare((PetscObject)pc->mat, MATNEST, &isnest));
524: }
525: if (!isnest) nf = jac->bs;
526: else PetscCall(MatNestGetSize(M, &nf, NULL));
527: for (i = 0; i < nf; i++) {
528: char splitname[8];
530: PetscCall(PetscSNPrintf(splitname, sizeof(splitname), "%" PetscInt_FMT, i));
531: PetscCall(PCFieldSplitSetFields(pc, splitname, 1, &i, &i));
532: }
533: jac->defaultsplit = PETSC_TRUE;
534: }
535: }
536: }
537: } else if (jac->nsplits == 1) {
538: IS is2;
539: PetscInt nmin, nmax;
541: PetscCheck(ilink->is, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Must provide at least two sets of fields to PCFieldSplit()");
542: PetscCall(MatGetOwnershipRange(pc->mat, &nmin, &nmax));
543: PetscCall(ISComplement(ilink->is, nmin, nmax, &is2));
544: PetscCall(PCFieldSplitSetIS(pc, "1", is2));
545: PetscCall(ISDestroy(&is2));
546: }
548: PetscCheck(jac->nsplits >= 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Unhandled case, must have at least two fields, not %" PetscInt_FMT, jac->nsplits);
549: PetscFunctionReturn(PETSC_SUCCESS);
550: }
552: static PetscErrorCode MatGolubKahanComputeExplicitOperator(Mat A, Mat B, Mat C, Mat *H, PetscReal gkbnu)
553: {
554: Mat BT, T;
555: PetscReal nrmT, nrmB;
557: PetscFunctionBegin;
558: PetscCall(MatHermitianTranspose(C, MAT_INITIAL_MATRIX, &T)); /* Test if augmented matrix is symmetric */
559: PetscCall(MatAXPY(T, -1.0, B, DIFFERENT_NONZERO_PATTERN));
560: PetscCall(MatNorm(T, NORM_1, &nrmT));
561: PetscCall(MatNorm(B, NORM_1, &nrmB));
562: PetscCheck(nrmB <= 0 || nrmT / nrmB < PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix is not symmetric/Hermitian, GKB is not applicable.");
564: /* Compute augmented Lagrangian matrix H = A00 + nu*A01*A01'. This corresponds to */
565: /* setting N := 1/nu*I in [Ar13]. */
566: PetscCall(MatHermitianTranspose(B, MAT_INITIAL_MATRIX, &BT));
567: PetscCall(MatMatMult(B, BT, MAT_INITIAL_MATRIX, PETSC_CURRENT, H)); /* H = A01*A01' */
568: PetscCall(MatAYPX(*H, gkbnu, A, DIFFERENT_NONZERO_PATTERN)); /* H = A00 + nu*A01*A01' */
570: PetscCall(MatDestroy(&BT));
571: PetscCall(MatDestroy(&T));
572: PetscFunctionReturn(PETSC_SUCCESS);
573: }
575: PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(PetscOptions, const char pre[], const char name[], const char *option[], const char *value[], PetscBool *flg);
577: static PetscErrorCode PCSetUp_FieldSplit(PC pc)
578: {
579: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
580: PC_FieldSplitLink ilink;
581: PetscInt i, nsplit;
582: PetscBool matnest = PETSC_FALSE;
584: PetscFunctionBegin;
585: pc->failedreason = PC_NOERROR;
586: PetscCall(PCFieldSplitSetDefaults(pc));
587: nsplit = jac->nsplits;
588: ilink = jac->head;
589: if (pc->pmat) PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATNEST, &matnest));
591: /* get the matrices for each split */
592: if (!jac->issetup) {
593: PetscInt rstart, rend, nslots, bs;
595: jac->issetup = PETSC_TRUE;
597: /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */
598: if (jac->defaultsplit || !ilink->is) {
599: if (jac->bs <= 0) jac->bs = nsplit;
600: }
602: /* MatCreateSubMatrix() for [S]BAIJ matrices can only work if the indices include entire blocks of the matrix */
603: PetscCall(MatGetBlockSize(pc->pmat, &bs));
604: if (bs > 1 && (jac->bs <= bs || jac->bs % bs)) {
605: PetscBool blk;
607: PetscCall(PetscObjectTypeCompareAny((PetscObject)pc->pmat, &blk, MATBAIJ, MATSBAIJ, MATSEQBAIJ, MATSEQSBAIJ, MATMPIBAIJ, MATMPISBAIJ, NULL));
608: PetscCheck(!blk, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Cannot use MATBAIJ with PCFIELDSPLIT and currently set matrix and PC blocksizes");
609: }
611: if (!matnest) { /* use the matrix blocksize and stride IS to determine the index sets that define the submatrices */
612: bs = jac->bs;
613: PetscCall(MatGetOwnershipRange(pc->pmat, &rstart, &rend));
614: nslots = (rend - rstart) / bs;
615: for (i = 0; i < nsplit; i++) {
616: if (jac->defaultsplit) {
617: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + i, nsplit, &ilink->is));
618: PetscCall(PetscObjectReference((PetscObject)ilink->is));
619: ilink->is_col = ilink->is;
620: } else if (!ilink->is) {
621: PetscBool same_fields = PETSC_TRUE;
623: for (PetscInt k = 0; k < ilink->nfields; k++) {
624: if (ilink->fields[k] != ilink->fields_col[k]) same_fields = PETSC_FALSE;
625: }
627: if (ilink->nfields > 1) {
628: PetscInt *ii, *jj, j, k, nfields = ilink->nfields, *fields = ilink->fields, *fields_col = ilink->fields_col;
630: PetscCall(PetscMalloc1(ilink->nfields * nslots, &ii));
631: if (!same_fields) PetscCall(PetscMalloc1(ilink->nfields * nslots, &jj));
632: for (j = 0; j < nslots; j++) {
633: for (k = 0; k < nfields; k++) {
634: ii[nfields * j + k] = rstart + bs * j + fields[k];
635: if (!same_fields) jj[nfields * j + k] = rstart + bs * j + fields_col[k];
636: }
637: }
638: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), nslots * nfields, ii, PETSC_OWN_POINTER, &ilink->is));
639: if (!same_fields) PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), nslots * nfields, jj, PETSC_OWN_POINTER, &ilink->is_col));
640: else {
641: PetscCall(PetscObjectReference((PetscObject)ilink->is));
642: ilink->is_col = ilink->is;
643: }
644: PetscCall(ISSetBlockSize(ilink->is, nfields));
645: PetscCall(ISSetBlockSize(ilink->is_col, nfields));
646: } else {
647: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + ilink->fields[0], bs, &ilink->is));
648: if (!same_fields) PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + ilink->fields_col[0], bs, &ilink->is_col));
649: else {
650: PetscCall(PetscObjectReference((PetscObject)ilink->is));
651: ilink->is_col = ilink->is;
652: }
653: }
654: }
655: ilink = ilink->next;
656: }
657: } else { /* use the IS that define the MATNEST to determine the index sets that define the submatrices */
658: IS *rowis, *colis, *ises = NULL;
659: PetscInt mis, nis;
661: PetscCall(MatNestGetSize(pc->pmat, &mis, &nis));
662: PetscCall(PetscMalloc2(mis, &rowis, nis, &colis));
663: PetscCall(MatNestGetISs(pc->pmat, rowis, colis));
664: if (!jac->defaultsplit) PetscCall(PetscMalloc1(mis, &ises));
666: for (i = 0; i < nsplit; i++) {
667: if (jac->defaultsplit) {
668: PetscCall(ISDuplicate(rowis[i], &ilink->is));
669: PetscCall(PetscObjectReference((PetscObject)ilink->is));
670: ilink->is_col = ilink->is;
671: } else if (!ilink->is) {
672: if (ilink->nfields > 1) {
673: for (PetscInt j = 0; j < ilink->nfields; j++) ises[j] = rowis[ilink->fields[j]];
674: PetscCall(ISConcatenate(PetscObjectComm((PetscObject)pc), ilink->nfields, ises, &ilink->is));
675: } else {
676: PetscCall(ISDuplicate(rowis[ilink->fields[0]], &ilink->is));
677: }
678: PetscCall(PetscObjectReference((PetscObject)ilink->is));
679: ilink->is_col = ilink->is;
680: }
681: ilink = ilink->next;
682: }
683: PetscCall(PetscFree2(rowis, colis));
684: PetscCall(PetscFree(ises));
685: }
686: }
688: ilink = jac->head;
689: if (!jac->pmat) {
690: Vec xtmp;
692: PetscCall(MatCreateVecs(pc->pmat, &xtmp, NULL));
693: PetscCall(PetscMalloc1(nsplit, &jac->pmat));
694: PetscCall(PetscMalloc2(nsplit, &jac->x, nsplit, &jac->y));
695: for (i = 0; i < nsplit; i++) {
696: MatNullSpace sp;
698: /* Check for matrix attached to IS */
699: PetscCall(PetscObjectQuery((PetscObject)ilink->is, "pmat", (PetscObject *)&jac->pmat[i]));
700: if (jac->pmat[i]) {
701: PetscCall(PetscObjectReference((PetscObject)jac->pmat[i]));
702: if (jac->type == PC_COMPOSITE_SCHUR) {
703: jac->schur_user = jac->pmat[i];
705: PetscCall(PetscObjectReference((PetscObject)jac->schur_user));
706: }
707: } else {
708: const char *prefix;
709: PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ilink->is_col, MAT_INITIAL_MATRIX, &jac->pmat[i]));
710: PetscCall(MatGetOptionsPrefix(jac->pmat[i], &prefix));
711: if (!prefix) {
712: PetscCall(KSPGetOptionsPrefix(ilink->ksp, &prefix));
713: PetscCall(MatSetOptionsPrefix(jac->pmat[i], prefix));
714: }
715: PetscCall(MatSetFromOptions(jac->pmat[i]));
716: PetscCall(MatViewFromOptions(jac->pmat[i], NULL, "-mat_view"));
717: }
718: /* create work vectors for each split */
719: PetscCall(MatCreateVecs(jac->pmat[i], &jac->x[i], &jac->y[i]));
720: ilink->x = jac->x[i];
721: ilink->y = jac->y[i];
722: ilink->z = NULL;
723: /* compute scatter contexts needed by multiplicative versions and non-default splits */
724: PetscCall(VecScatterCreate(xtmp, ilink->is, jac->x[i], NULL, &ilink->sctx));
725: PetscCall(PetscObjectQuery((PetscObject)ilink->is, "nearnullspace", (PetscObject *)&sp));
726: if (sp) PetscCall(MatSetNearNullSpace(jac->pmat[i], sp));
727: ilink = ilink->next;
728: }
729: PetscCall(VecDestroy(&xtmp));
730: } else {
731: MatReuse scall;
732: MatNullSpace *nullsp = NULL;
734: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
735: PetscCall(MatGetNullSpaces(nsplit, jac->pmat, &nullsp));
736: for (i = 0; i < nsplit; i++) PetscCall(MatDestroy(&jac->pmat[i]));
737: scall = MAT_INITIAL_MATRIX;
738: } else scall = MAT_REUSE_MATRIX;
740: for (i = 0; i < nsplit; i++) {
741: Mat pmat;
743: /* Check for matrix attached to IS */
744: PetscCall(PetscObjectQuery((PetscObject)ilink->is, "pmat", (PetscObject *)&pmat));
745: if (!pmat) PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ilink->is_col, scall, &jac->pmat[i]));
746: ilink = ilink->next;
747: }
748: if (nullsp) PetscCall(MatRestoreNullSpaces(nsplit, jac->pmat, &nullsp));
749: }
750: if (jac->diag_use_amat) {
751: ilink = jac->head;
752: if (!jac->mat) {
753: PetscCall(PetscMalloc1(nsplit, &jac->mat));
754: for (i = 0; i < nsplit; i++) {
755: PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ilink->is_col, MAT_INITIAL_MATRIX, &jac->mat[i]));
756: ilink = ilink->next;
757: }
758: } else {
759: MatReuse scall;
760: MatNullSpace *nullsp = NULL;
762: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
763: PetscCall(MatGetNullSpaces(nsplit, jac->mat, &nullsp));
764: for (i = 0; i < nsplit; i++) PetscCall(MatDestroy(&jac->mat[i]));
765: scall = MAT_INITIAL_MATRIX;
766: } else scall = MAT_REUSE_MATRIX;
768: for (i = 0; i < nsplit; i++) {
769: PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ilink->is_col, scall, &jac->mat[i]));
770: ilink = ilink->next;
771: }
772: if (nullsp) PetscCall(MatRestoreNullSpaces(nsplit, jac->mat, &nullsp));
773: }
774: } else {
775: jac->mat = jac->pmat;
776: }
778: /* Check for null space attached to IS */
779: ilink = jac->head;
780: for (i = 0; i < nsplit; i++) {
781: MatNullSpace sp;
783: PetscCall(PetscObjectQuery((PetscObject)ilink->is, "nullspace", (PetscObject *)&sp));
784: if (sp) PetscCall(MatSetNullSpace(jac->mat[i], sp));
785: ilink = ilink->next;
786: }
788: if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR && jac->type != PC_COMPOSITE_GKB) {
789: /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */
790: /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */
791: ilink = jac->head;
792: if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) {
793: /* special case need where Afield[0] is not needed and only certain columns of Afield[1] are needed since update is only on those rows of the solution */
794: if (!jac->Afield) {
795: PetscCall(PetscCalloc1(nsplit, &jac->Afield));
796: if (jac->offdiag_use_amat) {
797: PetscCall(MatCreateSubMatrix(pc->mat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->Afield[1]));
798: } else {
799: PetscCall(MatCreateSubMatrix(pc->pmat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->Afield[1]));
800: }
801: } else {
802: MatReuse scall;
804: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
805: PetscCall(MatDestroy(&jac->Afield[1]));
806: scall = MAT_INITIAL_MATRIX;
807: } else scall = MAT_REUSE_MATRIX;
809: if (jac->offdiag_use_amat) {
810: PetscCall(MatCreateSubMatrix(pc->mat, ilink->next->is, ilink->is, scall, &jac->Afield[1]));
811: } else {
812: PetscCall(MatCreateSubMatrix(pc->pmat, ilink->next->is, ilink->is, scall, &jac->Afield[1]));
813: }
814: }
815: } else {
816: if (!jac->Afield) {
817: PetscCall(PetscMalloc1(nsplit, &jac->Afield));
818: for (i = 0; i < nsplit; i++) {
819: if (jac->offdiag_use_amat) {
820: PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, NULL, MAT_INITIAL_MATRIX, &jac->Afield[i]));
821: } else {
822: PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, NULL, MAT_INITIAL_MATRIX, &jac->Afield[i]));
823: }
824: ilink = ilink->next;
825: }
826: } else {
827: MatReuse scall;
828: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
829: for (i = 0; i < nsplit; i++) PetscCall(MatDestroy(&jac->Afield[i]));
830: scall = MAT_INITIAL_MATRIX;
831: } else scall = MAT_REUSE_MATRIX;
833: for (i = 0; i < nsplit; i++) {
834: if (jac->offdiag_use_amat) {
835: PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, NULL, scall, &jac->Afield[i]));
836: } else {
837: PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, NULL, scall, &jac->Afield[i]));
838: }
839: ilink = ilink->next;
840: }
841: }
842: }
843: }
845: if (jac->type == PC_COMPOSITE_SCHUR) {
846: PetscBool isset, isspd = PETSC_FALSE, issym = PETSC_FALSE, flg;
847: char lscname[256];
848: PetscObject LSC_L;
850: PetscCheck(nsplit == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "To use Schur complement preconditioner you must have exactly 2 fields");
852: /* If pc->mat is SPD, don't scale by -1 the Schur complement */
853: PetscCall(MatIsSPDKnown(pc->pmat, &isset, &isspd));
854: if (jac->schurscale == (PetscScalar)-1.0) jac->schurscale = (isset && isspd) ? 1.0 : -1.0;
855: PetscCall(MatIsSymmetricKnown(pc->pmat, &isset, &issym));
857: PetscCall(PetscObjectTypeCompareAny(jac->offdiag_use_amat ? (PetscObject)pc->mat : (PetscObject)pc->pmat, &flg, MATSEQSBAIJ, MATMPISBAIJ, ""));
859: if (jac->schur) {
860: KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper;
861: MatReuse scall;
863: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
864: scall = MAT_INITIAL_MATRIX;
865: PetscCall(MatDestroy(&jac->B));
866: PetscCall(MatDestroy(&jac->C));
867: } else scall = MAT_REUSE_MATRIX;
869: PetscCall(MatSchurComplementGetKSP(jac->schur, &kspInner));
870: ilink = jac->head;
871: PetscCall(MatCreateSubMatrix(jac->offdiag_use_amat ? pc->mat : pc->pmat, ilink->is, ilink->next->is, scall, &jac->B));
872: if (!flg) PetscCall(MatCreateSubMatrix(jac->offdiag_use_amat ? pc->mat : pc->pmat, ilink->next->is, ilink->is, scall, &jac->C));
873: else {
874: PetscCall(MatIsHermitianKnown(jac->offdiag_use_amat ? pc->mat : pc->pmat, &isset, &flg));
875: if (isset && flg) PetscCall(MatCreateHermitianTranspose(jac->B, &jac->C));
876: else PetscCall(MatCreateTranspose(jac->B, &jac->C));
877: }
878: ilink = ilink->next;
879: PetscCall(MatSchurComplementUpdateSubMatrices(jac->schur, jac->mat[0], jac->pmat[0], jac->B, jac->C, jac->mat[1]));
880: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
881: PetscCall(MatDestroy(&jac->schurp));
882: PetscCall(MatSchurComplementGetPmat(jac->schur, MAT_INITIAL_MATRIX, &jac->schurp));
883: } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL && jac->kspupper != jac->head->ksp) {
884: PetscCall(MatDestroy(&jac->schur_user));
885: PetscCall(MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user));
886: }
887: if (kspA != kspInner) PetscCall(KSPSetOperators(kspA, jac->mat[0], jac->pmat[0]));
888: if (kspUpper != kspA) PetscCall(KSPSetOperators(kspUpper, jac->mat[0], jac->pmat[0]));
889: PetscCall(KSPSetOperators(jac->kspschur, jac->schur, FieldSplitSchurPre(jac)));
890: } else {
891: const char *Dprefix;
892: char schurprefix[256], schurmatprefix[256];
893: char schurtestoption[256];
894: MatNullSpace sp;
895: KSP kspt;
897: /* extract the A01 and A10 matrices */
898: ilink = jac->head;
899: PetscCall(MatCreateSubMatrix(jac->offdiag_use_amat ? pc->mat : pc->pmat, ilink->is, ilink->next->is, MAT_INITIAL_MATRIX, &jac->B));
900: if (!flg) PetscCall(MatCreateSubMatrix(jac->offdiag_use_amat ? pc->mat : pc->pmat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->C));
901: else {
902: PetscCall(MatIsHermitianKnown(jac->offdiag_use_amat ? pc->mat : pc->pmat, &isset, &flg));
903: if (isset && flg) PetscCall(MatCreateHermitianTranspose(jac->B, &jac->C));
904: else PetscCall(MatCreateTranspose(jac->B, &jac->C));
905: }
906: ilink = ilink->next;
907: /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
908: PetscCall(MatCreate(((PetscObject)jac->mat[0])->comm, &jac->schur));
909: PetscCall(MatSetType(jac->schur, MATSCHURCOMPLEMENT));
910: PetscCall(MatSchurComplementSetSubMatrices(jac->schur, jac->mat[0], jac->pmat[0], jac->B, jac->C, jac->mat[1]));
911: PetscCall(PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname));
912: PetscCall(MatSetOptionsPrefix(jac->schur, schurmatprefix));
913: PetscCall(MatSchurComplementGetKSP(jac->schur, &kspt));
914: PetscCall(KSPSetOptionsPrefix(kspt, schurmatprefix));
916: /* Note: this is not true in general */
917: PetscCall(MatGetNullSpace(jac->mat[1], &sp));
918: if (sp) PetscCall(MatSetNullSpace(jac->schur, sp));
920: PetscCall(PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname));
921: PetscCall(PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, schurtestoption, NULL, NULL, &flg));
922: if (flg) {
923: DM dmInner;
924: KSP kspInner;
925: PC pcInner;
927: PetscCall(MatSchurComplementGetKSP(jac->schur, &kspInner));
928: PetscCall(KSPReset(kspInner));
929: PetscCall(KSPSetOperators(kspInner, jac->mat[0], jac->pmat[0]));
930: PetscCall(PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname));
931: /* Indent this deeper to emphasize the "inner" nature of this solver. */
932: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject)pc, 2));
933: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspInner->pc, (PetscObject)pc, 2));
934: PetscCall(KSPSetOptionsPrefix(kspInner, schurprefix));
936: /* Set DM for new solver */
937: PetscCall(KSPGetDM(jac->head->ksp, &dmInner));
938: PetscCall(KSPSetDM(kspInner, dmInner));
939: PetscCall(KSPSetDMActive(kspInner, KSP_DMACTIVE_ALL, PETSC_FALSE));
941: /* Defaults to PCKSP as preconditioner */
942: PetscCall(KSPGetPC(kspInner, &pcInner));
943: PetscCall(PCSetType(pcInner, PCKSP));
944: PetscCall(PCKSPSetKSP(pcInner, jac->head->ksp));
945: } else {
946: /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or
947: * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact,
948: * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for
949: * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make
950: * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used
951: * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */
952: PetscCall(KSPSetType(jac->head->ksp, KSPGMRES));
953: PetscCall(MatSchurComplementSetKSP(jac->schur, jac->head->ksp));
954: }
955: PetscCall(KSPSetOperators(jac->head->ksp, jac->mat[0], jac->pmat[0]));
956: PetscCall(KSPSetFromOptions(jac->head->ksp));
957: PetscCall(MatSetFromOptions(jac->schur));
959: PetscCall(PetscObjectTypeCompare((PetscObject)jac->schur, MATSCHURCOMPLEMENT, &flg));
960: if (flg) { /* Need to do this otherwise PCSetUp_KSP will overwrite the amat of jac->head->ksp */
961: KSP kspInner;
962: PC pcInner;
964: PetscCall(MatSchurComplementGetKSP(jac->schur, &kspInner));
965: PetscCall(KSPGetPC(kspInner, &pcInner));
966: PetscCall(PetscObjectTypeCompare((PetscObject)pcInner, PCKSP, &flg));
967: if (flg) {
968: KSP ksp;
970: PetscCall(PCKSPGetKSP(pcInner, &ksp));
971: if (ksp == jac->head->ksp) PetscCall(PCSetUseAmat(pcInner, PETSC_TRUE));
972: }
973: }
974: PetscCall(PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname));
975: PetscCall(PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, schurtestoption, NULL, NULL, &flg));
976: if (flg) {
977: DM dmInner;
979: PetscCall(PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname));
980: PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper));
981: PetscCall(KSPSetNestLevel(jac->kspupper, pc->kspnestlevel));
982: PetscCall(KSPSetErrorIfNotConverged(jac->kspupper, pc->erroriffailure));
983: PetscCall(KSPSetOptionsPrefix(jac->kspupper, schurprefix));
984: PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->kspupper, (PetscObject)pc, 1));
985: PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->kspupper->pc, (PetscObject)pc, 1));
986: PetscCall(KSPGetDM(jac->head->ksp, &dmInner));
987: PetscCall(KSPSetDM(jac->kspupper, dmInner));
988: PetscCall(KSPSetDMActive(jac->kspupper, KSP_DMACTIVE_ALL, PETSC_FALSE));
989: PetscCall(KSPSetFromOptions(jac->kspupper));
990: PetscCall(KSPSetOperators(jac->kspupper, jac->mat[0], jac->pmat[0]));
991: PetscCall(VecDuplicate(jac->head->x, &jac->head->z));
992: } else {
993: jac->kspupper = jac->head->ksp;
994: PetscCall(PetscObjectReference((PetscObject)jac->head->ksp));
995: }
997: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) PetscCall(MatSchurComplementGetPmat(jac->schur, MAT_INITIAL_MATRIX, &jac->schurp));
998: PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspschur));
999: PetscCall(KSPSetNestLevel(jac->kspschur, pc->kspnestlevel));
1000: PetscCall(KSPSetErrorIfNotConverged(jac->kspschur, pc->erroriffailure));
1001: PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->kspschur, (PetscObject)pc, 1));
1002: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
1003: PC pcschur;
1004: PetscCall(KSPGetPC(jac->kspschur, &pcschur));
1005: PetscCall(PCSetType(pcschur, PCNONE));
1006: /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
1007: } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) {
1008: if (jac->schurfactorization != PC_FIELDSPLIT_SCHUR_FACT_FULL || jac->kspupper != jac->head->ksp) PetscCall(MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user));
1009: }
1010: PetscCall(KSPSetOperators(jac->kspschur, jac->schur, FieldSplitSchurPre(jac)));
1011: PetscCall(KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix));
1012: PetscCall(KSPSetOptionsPrefix(jac->kspschur, Dprefix));
1013: /* propagate DM */
1014: {
1015: DM sdm;
1016: PetscCall(KSPGetDM(jac->head->next->ksp, &sdm));
1017: if (sdm) {
1018: PetscCall(KSPSetDM(jac->kspschur, sdm));
1019: PetscCall(KSPSetDMActive(jac->kspschur, KSP_DMACTIVE_ALL, PETSC_FALSE));
1020: }
1021: }
1022: /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
1023: /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */
1024: PetscCall(KSPSetFromOptions(jac->kspschur));
1025: }
1026: PetscCall(MatAssemblyBegin(jac->schur, MAT_FINAL_ASSEMBLY));
1027: PetscCall(MatAssemblyEnd(jac->schur, MAT_FINAL_ASSEMBLY));
1028: if (issym) PetscCall(MatSetOption(jac->schur, MAT_SYMMETRIC, PETSC_TRUE));
1029: if (isspd) PetscCall(MatSetOption(jac->schur, MAT_SPD, PETSC_TRUE));
1031: /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */
1032: PetscCall(PetscSNPrintf(lscname, sizeof(lscname), "%s_LSC_L", ilink->splitname));
1033: PetscCall(PetscObjectQuery((PetscObject)pc->mat, lscname, &LSC_L));
1034: if (!LSC_L) PetscCall(PetscObjectQuery((PetscObject)pc->pmat, lscname, &LSC_L));
1035: if (LSC_L) PetscCall(PetscObjectCompose((PetscObject)jac->schur, "LSC_L", LSC_L));
1036: PetscCall(PetscSNPrintf(lscname, sizeof(lscname), "%s_LSC_Lp", ilink->splitname));
1037: PetscCall(PetscObjectQuery((PetscObject)pc->pmat, lscname, &LSC_L));
1038: if (!LSC_L) PetscCall(PetscObjectQuery((PetscObject)pc->mat, lscname, &LSC_L));
1039: if (LSC_L) PetscCall(PetscObjectCompose((PetscObject)jac->schur, "LSC_Lp", LSC_L));
1040: } else if (jac->type == PC_COMPOSITE_GKB) {
1041: PetscCheck(nsplit == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "To use GKB preconditioner you must have exactly 2 fields");
1042: ilink = jac->head;
1043: PetscCall(MatCreateSubMatrix(jac->offdiag_use_amat ? pc->mat : pc->pmat, ilink->is, ilink->next->is, MAT_INITIAL_MATRIX, &jac->B));
1044: /* Create work vectors for GKB algorithm */
1045: PetscCall(VecDuplicate(ilink->x, &jac->u));
1046: PetscCall(VecDuplicate(ilink->x, &jac->Hu));
1047: PetscCall(VecDuplicate(ilink->x, &jac->w2));
1048: PetscCall(MatCreateSubMatrix(jac->offdiag_use_amat ? pc->mat : pc->pmat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->C));
1049: ilink = ilink->next;
1050: /* Create work vectors for GKB algorithm */
1051: PetscCall(VecDuplicate(ilink->x, &jac->v));
1052: PetscCall(VecDuplicate(ilink->x, &jac->d));
1053: PetscCall(VecDuplicate(ilink->x, &jac->w1));
1054: PetscCall(MatGolubKahanComputeExplicitOperator(jac->mat[0], jac->B, jac->C, &jac->H, jac->gkbnu));
1055: PetscCall(PetscCalloc1(jac->gkbdelay, &jac->vecz));
1057: ilink = jac->head;
1058: PetscCall(KSPSetOperators(ilink->ksp, jac->H, jac->H));
1059: if (!jac->suboptionsset) PetscCall(KSPSetFromOptions(ilink->ksp));
1060: /* Create gkb_monitor context */
1061: if (jac->gkbmonitor) {
1062: PetscInt tablevel;
1063: PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &jac->gkbviewer));
1064: PetscCall(PetscViewerSetType(jac->gkbviewer, PETSCVIEWERASCII));
1065: PetscCall(PetscObjectGetTabLevel((PetscObject)ilink->ksp, &tablevel));
1066: PetscCall(PetscViewerASCIISetTab(jac->gkbviewer, tablevel));
1067: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)ilink->ksp, 1));
1068: }
1069: } else {
1070: /* set up the individual splits' PCs */
1071: i = 0;
1072: ilink = jac->head;
1073: while (ilink) {
1074: PetscCall(KSPSetOperators(ilink->ksp, jac->mat[i], jac->pmat[i]));
1075: /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
1076: if (!jac->suboptionsset) PetscCall(KSPSetFromOptions(ilink->ksp));
1077: i++;
1078: ilink = ilink->next;
1079: }
1080: }
1082: /* Set coordinates to the sub PC objects whenever these are set */
1083: if (jac->coordinates_set) {
1084: PC pc_coords;
1085: if (jac->type == PC_COMPOSITE_SCHUR) {
1086: // Head is first block.
1087: PetscCall(KSPGetPC(jac->head->ksp, &pc_coords));
1088: PetscCall(PCSetCoordinates(pc_coords, jac->head->dim, jac->head->ndofs, jac->head->coords));
1089: // Second one is Schur block, but its KSP object is in kspschur.
1090: PetscCall(KSPGetPC(jac->kspschur, &pc_coords));
1091: PetscCall(PCSetCoordinates(pc_coords, jac->head->next->dim, jac->head->next->ndofs, jac->head->next->coords));
1092: } else if (jac->type == PC_COMPOSITE_GKB) {
1093: PetscCall(PetscInfo(pc, "Warning: Setting coordinates does nothing for the GKB Fieldpslit preconditioner\n"));
1094: } else {
1095: ilink = jac->head;
1096: while (ilink) {
1097: PetscCall(KSPGetPC(ilink->ksp, &pc_coords));
1098: PetscCall(PCSetCoordinates(pc_coords, ilink->dim, ilink->ndofs, ilink->coords));
1099: ilink = ilink->next;
1100: }
1101: }
1102: }
1104: jac->suboptionsset = PETSC_TRUE;
1105: PetscFunctionReturn(PETSC_SUCCESS);
1106: }
1108: static PetscErrorCode PCSetUpOnBlocks_FieldSplit_Schur(PC pc)
1109: {
1110: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1111: PC_FieldSplitLink ilinkA = jac->head;
1112: KSP kspA = ilinkA->ksp, kspUpper = jac->kspupper;
1114: PetscFunctionBegin;
1115: if (jac->schurfactorization == PC_FIELDSPLIT_SCHUR_FACT_FULL && kspUpper != kspA) {
1116: PetscCall(KSPSetUp(kspUpper));
1117: PetscCall(KSPSetUpOnBlocks(kspUpper));
1118: }
1119: PetscCall(KSPSetUp(kspA));
1120: PetscCall(KSPSetUpOnBlocks(kspA));
1121: if (jac->schurpre != PC_FIELDSPLIT_SCHUR_PRE_FULL) {
1122: PetscCall(KSPSetUp(jac->kspschur));
1123: PetscCall(KSPSetUpOnBlocks(jac->kspschur));
1124: } else if (kspUpper == kspA && jac->schurfactorization == PC_FIELDSPLIT_SCHUR_FACT_FULL) {
1125: Mat A;
1126: PetscInt m, M, N;
1127: VecType vtype;
1128: PetscMemType mtype;
1129: PetscScalar *array;
1131: PetscCall(MatGetSize(jac->B, &M, &N));
1132: PetscCall(MatGetLocalSize(jac->B, &m, NULL));
1133: PetscCall(MatGetVecType(jac->B, &vtype));
1134: PetscCall(VecGetArrayAndMemType(ilinkA->x, &array, &mtype));
1135: PetscCall(VecRestoreArrayAndMemType(ilinkA->x, &array));
1136: PetscCall(PetscObjectQuery((PetscObject)jac->schur, "AinvB", (PetscObject *)&A));
1137: if (A) {
1138: PetscInt P;
1140: PetscCall(MatGetSize(A, NULL, &P));
1141: if (P < N + 1) { // need to recreate AinvB, otherwise, the Schur complement won't be updated
1142: PetscCall(PetscObjectCompose((PetscObject)jac->schur, "AinvB", NULL));
1143: A = NULL;
1144: }
1145: }
1146: if (!A) {
1147: if (PetscMemTypeHost(mtype) || (!PetscDefined(HAVE_CUDA) && !PetscDefined(HAVE_HIP))) PetscCall(PetscMalloc1(m * (N + 1), &array));
1148: #if PetscDefined(HAVE_CUDA)
1149: else if (PetscMemTypeCUDA(mtype)) PetscCallCUDA(cudaMalloc((void **)&array, sizeof(PetscScalar) * m * (N + 1)));
1150: #endif
1151: #if PetscDefined(HAVE_HIP)
1152: else if (PetscMemTypeHIP(mtype)) PetscCallHIP(hipMalloc((void **)&array, sizeof(PetscScalar) * m * (N + 1)));
1153: #endif
1154: PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)jac->schur), vtype, m, PETSC_DECIDE, M, N + 1, -1, array, &A)); // number of columns of the Schur complement plus one
1155: PetscCall(PetscObjectCompose((PetscObject)jac->schur, "AinvB", (PetscObject)A));
1156: PetscCall(MatDestroy(&A));
1157: }
1158: }
1159: PetscFunctionReturn(PETSC_SUCCESS);
1160: }
1162: static PetscErrorCode PCSetUpOnBlocks_FieldSplit(PC pc)
1163: {
1164: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1165: PC_FieldSplitLink ilink = jac->head;
1167: PetscFunctionBegin;
1168: while (ilink) {
1169: PetscCall(KSPSetUp(ilink->ksp));
1170: PetscCall(KSPSetUpOnBlocks(ilink->ksp));
1171: ilink = ilink->next;
1172: }
1173: PetscFunctionReturn(PETSC_SUCCESS);
1174: }
1176: static PetscErrorCode PCSetUpOnBlocks_FieldSplit_GKB(PC pc)
1177: {
1178: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1179: PC_FieldSplitLink ilinkA = jac->head;
1180: KSP ksp = ilinkA->ksp;
1182: PetscFunctionBegin;
1183: PetscCall(KSPSetUp(ksp));
1184: PetscCall(KSPSetUpOnBlocks(ksp));
1185: PetscFunctionReturn(PETSC_SUCCESS);
1186: }
1188: static PetscErrorCode PCApply_FieldSplit_Schur(PC pc, Vec x, Vec y)
1189: {
1190: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1191: PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
1192: KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;
1193: Mat AinvB = NULL;
1194: PetscInt N, P;
1196: PetscFunctionBegin;
1197: switch (jac->schurfactorization) {
1198: case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
1199: /* [A00 0; 0 -S], positive definite, suitable for MINRES */
1200: PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1201: PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
1202: PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1203: PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1204: PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y));
1205: PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
1206: PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1207: PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1208: PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
1209: PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1210: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
1211: PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y));
1212: PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
1213: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
1214: PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1215: PetscCall(VecScale(ilinkD->y, jac->schurscale));
1216: PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1217: PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1218: PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1219: break;
1220: case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
1221: /* [A00 0; A10 S], suitable for left preconditioning */
1222: PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1223: PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1224: PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1225: PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y));
1226: PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
1227: PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1228: PetscCall(MatMult(jac->C, ilinkA->y, ilinkD->x));
1229: PetscCall(VecScale(ilinkD->x, -1.));
1230: PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
1231: PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1232: PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
1233: PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1234: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
1235: PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y));
1236: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
1237: PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
1238: PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1239: PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1240: PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1241: PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1242: break;
1243: case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
1244: /* [A00 A01; 0 S], suitable for right preconditioning */
1245: PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
1246: PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
1247: PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1248: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
1249: PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y));
1250: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
1251: PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
1252: PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1253: PetscCall(MatMult(jac->B, ilinkD->y, ilinkA->x));
1254: PetscCall(VecScale(ilinkA->x, -1.));
1255: PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD));
1256: PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1257: PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD));
1258: PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1259: PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y));
1260: PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
1261: PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1262: PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1263: PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1264: PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1265: break;
1266: case PC_FIELDSPLIT_SCHUR_FACT_FULL:
1267: /* [1 0; A10 A00^{-1} 1] [A00 0; 0 S] [1 A00^{-1}A01; 0 1] */
1268: PetscCall(MatGetSize(jac->B, NULL, &P));
1269: N = P;
1270: PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1271: PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1272: PetscCall(PetscLogEventBegin(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->y, NULL));
1273: if (kspUpper == kspA) {
1274: PetscCall(PetscObjectQuery((PetscObject)jac->schur, "AinvB", (PetscObject *)&AinvB));
1275: if (AinvB) {
1276: PetscCall(MatGetSize(AinvB, NULL, &N));
1277: if (N > P) { // first time PCApply_FieldSplit_Schur() is called
1278: PetscMemType mtype;
1279: Vec c = NULL;
1280: PetscScalar *array;
1281: PetscInt m, M;
1283: PetscCall(MatGetSize(jac->B, &M, NULL));
1284: PetscCall(MatGetLocalSize(jac->B, &m, NULL));
1285: PetscCall(MatDenseGetArrayAndMemType(AinvB, &array, &mtype));
1286: if (PetscMemTypeHost(mtype) || (!PetscDefined(HAVE_CUDA) && !PetscDefined(HAVE_HIP))) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)jac->schur), 1, m, M, array + m * P, &c));
1287: #if PetscDefined(HAVE_CUDA)
1288: else if (PetscMemTypeCUDA(mtype)) PetscCall(VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)jac->schur), 1, m, M, array + m * P, &c));
1289: #endif
1290: #if PetscDefined(HAVE_HIP)
1291: else if (PetscMemTypeHIP(mtype)) PetscCall(VecCreateMPIHIPWithArray(PetscObjectComm((PetscObject)jac->schur), 1, m, M, array + m * P, &c));
1292: #endif
1293: PetscCall(MatDenseRestoreArrayAndMemType(AinvB, &array));
1294: PetscCall(VecCopy(ilinkA->x, c));
1295: PetscCall(MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user));
1296: PetscCall(KSPSetOperators(jac->kspschur, jac->schur, jac->schur_user));
1297: PetscCall(VecCopy(c, ilinkA->y)); // retrieve the solution as the last column of the composed Mat
1298: PetscCall(VecDestroy(&c));
1299: }
1300: }
1301: }
1302: if (N == P) PetscCall(KSPSolve(kspLower, ilinkA->x, ilinkA->y));
1303: PetscCall(KSPCheckSolve(kspLower, pc, ilinkA->y));
1304: PetscCall(PetscLogEventEnd(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->y, NULL));
1305: PetscCall(MatMult(jac->C, ilinkA->y, ilinkD->x));
1306: PetscCall(VecScale(ilinkD->x, -1.0));
1307: PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
1308: PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
1310: PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1311: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
1312: PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y));
1313: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
1314: PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
1315: PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1316: PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1318: if (kspUpper == kspA) {
1319: if (!AinvB) {
1320: PetscCall(MatMult(jac->B, ilinkD->y, ilinkA->y));
1321: PetscCall(VecAXPY(ilinkA->x, -1.0, ilinkA->y));
1322: PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1323: PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y));
1324: PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
1325: PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1326: } else PetscCall(MatMultAdd(AinvB, ilinkD->y, ilinkA->y, ilinkA->y));
1327: } else {
1328: PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1329: PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y));
1330: PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
1331: PetscCall(MatMult(jac->B, ilinkD->y, ilinkA->x));
1332: PetscCall(PetscLogEventBegin(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->z, NULL));
1333: PetscCall(KSPSolve(kspUpper, ilinkA->x, ilinkA->z));
1334: PetscCall(KSPCheckSolve(kspUpper, pc, ilinkA->z));
1335: PetscCall(PetscLogEventEnd(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->z, NULL));
1336: PetscCall(VecAXPY(ilinkA->y, -1.0, ilinkA->z));
1337: }
1338: PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1339: PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1340: PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1341: }
1342: PetscFunctionReturn(PETSC_SUCCESS);
1343: }
1345: /*
1346: PCFieldSplitCreateWorkMats_Private - Allocate per-field dense work matrices for multi-RHS
1348: Input Parameters:
1349: + pc - the PC context
1350: - X - matrix to copy column-layout from
1352: Notes:
1353: If matrices already exist with correct column count, they are reused.
1354: If column count changed, old matrices are destroyed and new ones created.
1355: */
1356: static PetscErrorCode PCFieldSplitCreateWorkMats_Private(PC pc, Mat X)
1357: {
1358: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1359: PC_FieldSplitLink ilink = jac->head;
1360: PetscInt mx, Mx, my, My, N;
1362: PetscFunctionBegin;
1363: while (ilink) {
1364: /* check if reallocation needed (previous allocation with wrong column count) */
1365: if (ilink->X) {
1366: PetscCall(MatGetSize(ilink->X, NULL, &N));
1367: if (N != X->cmap->N) {
1368: PetscCall(MatDestroy(&ilink->X));
1369: PetscCall(MatDestroy(&ilink->Y));
1370: PetscCall(MatDestroy(&ilink->Z));
1371: }
1372: }
1373: /* create if needed */
1374: if (!ilink->X) {
1375: VecType xtype, ytype;
1377: PetscCall(VecGetType(ilink->x, &xtype));
1378: PetscCall(VecGetType(ilink->y, &ytype));
1379: PetscCall(VecGetLocalSize(ilink->x, &mx));
1380: PetscCall(VecGetSize(ilink->x, &Mx));
1381: PetscCall(VecGetLocalSize(ilink->y, &my));
1382: PetscCall(VecGetSize(ilink->y, &My));
1383: /* use default lda */
1384: PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)pc), xtype, mx, X->cmap->n, Mx, X->cmap->N, -1, NULL, &ilink->X));
1385: PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)pc), ytype, my, X->cmap->n, My, X->cmap->N, -1, NULL, &ilink->Y));
1386: }
1387: ilink = ilink->next;
1388: }
1389: PetscFunctionReturn(PETSC_SUCCESS);
1390: }
1392: static PetscErrorCode PCMatApply_FieldSplit_Schur(PC pc, Mat X, Mat Y)
1393: {
1394: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1395: PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
1396: KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;
1397: Mat AinvB = NULL;
1398: PetscInt N, P;
1400: PetscFunctionBegin;
1401: /* create working matrices with the correct number of columns */
1402: PetscCall(PCFieldSplitCreateWorkMats_Private(pc, X));
1403: switch (jac->schurfactorization) {
1404: case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
1405: /* [A00 0; 0 -S], positive definite, suitable for MINRES */
1406: PetscCall(MatDenseScatter_Private(ilinkA->sctx, X, ilinkA->X, INSERT_VALUES, SCATTER_FORWARD));
1407: PetscCall(MatDenseScatter_Private(ilinkD->sctx, X, ilinkD->X, INSERT_VALUES, SCATTER_FORWARD));
1408: PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->X, ilinkA->Y, NULL));
1409: PetscCall(KSPMatSolve(kspA, ilinkA->X, ilinkA->Y));
1410: PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->X, ilinkA->Y, NULL));
1411: PetscCall(MatDenseScatter_Private(ilinkA->sctx, ilinkA->Y, Y, INSERT_VALUES, SCATTER_REVERSE));
1412: PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->X, ilinkD->Y, NULL));
1413: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
1414: PetscCall(KSPMatSolve(jac->kspschur, ilinkD->X, ilinkD->Y));
1415: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
1416: PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->X, ilinkD->Y, NULL));
1417: PetscCall(MatScale(ilinkD->Y, jac->schurscale));
1418: PetscCall(MatDenseScatter_Private(ilinkD->sctx, ilinkD->Y, Y, INSERT_VALUES, SCATTER_REVERSE));
1419: break;
1420: case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
1421: /* [A00 0; A10 S], suitable for left preconditioning */
1422: PetscCall(MatDenseScatter_Private(ilinkA->sctx, X, ilinkA->X, INSERT_VALUES, SCATTER_FORWARD));
1423: PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->X, ilinkA->Y, NULL));
1424: PetscCall(KSPMatSolve(kspA, ilinkA->X, ilinkA->Y));
1425: PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->X, ilinkA->Y, NULL));
1426: PetscCall(MatMatMult(jac->C, ilinkA->Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilinkD->X));
1427: PetscCall(MatScale(ilinkD->X, -1.0));
1428: PetscCall(MatDenseScatter_Private(ilinkD->sctx, X, ilinkD->X, ADD_VALUES, SCATTER_FORWARD));
1429: PetscCall(MatDenseScatter_Private(ilinkA->sctx, ilinkA->Y, Y, INSERT_VALUES, SCATTER_REVERSE));
1430: PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->X, ilinkD->Y, NULL));
1431: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
1432: PetscCall(KSPMatSolve(jac->kspschur, ilinkD->X, ilinkD->Y));
1433: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
1434: PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->X, ilinkD->Y, NULL));
1435: PetscCall(MatDenseScatter_Private(ilinkD->sctx, ilinkD->Y, Y, INSERT_VALUES, SCATTER_REVERSE));
1436: break;
1437: case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
1438: /* [A00 A01; 0 S], suitable for right preconditioning */
1439: PetscCall(MatDenseScatter_Private(ilinkD->sctx, X, ilinkD->X, INSERT_VALUES, SCATTER_FORWARD));
1440: PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->X, ilinkD->Y, NULL));
1441: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
1442: PetscCall(KSPMatSolve(jac->kspschur, ilinkD->X, ilinkD->Y));
1443: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
1444: PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->X, ilinkD->Y, NULL));
1445: PetscCall(MatMatMult(jac->B, ilinkD->Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilinkA->X));
1446: PetscCall(MatScale(ilinkA->X, -1.0));
1447: PetscCall(MatDenseScatter_Private(ilinkA->sctx, X, ilinkA->X, ADD_VALUES, SCATTER_FORWARD));
1448: PetscCall(MatDenseScatter_Private(ilinkD->sctx, ilinkD->Y, Y, INSERT_VALUES, SCATTER_REVERSE));
1449: PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->X, ilinkA->Y, NULL));
1450: PetscCall(KSPMatSolve(kspA, ilinkA->X, ilinkA->Y));
1451: PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->X, ilinkA->Y, NULL));
1452: PetscCall(MatDenseScatter_Private(ilinkA->sctx, ilinkA->Y, Y, INSERT_VALUES, SCATTER_REVERSE));
1453: break;
1454: case PC_FIELDSPLIT_SCHUR_FACT_FULL:
1455: /* [1 0; A10 A00^{-1} 1] [A00 0; 0 S] [1 A00^{-1}A01; 0 1] */
1456: PetscCall(MatGetSize(jac->B, NULL, &P));
1457: N = P;
1458: PetscCall(MatDenseScatter_Private(ilinkA->sctx, X, ilinkA->X, INSERT_VALUES, SCATTER_FORWARD));
1459: PetscCall(PetscLogEventBegin(KSP_Solve_FS_L, kspLower, ilinkA->X, ilinkA->Y, NULL));
1460: if (kspUpper == kspA) {
1461: PetscCall(PetscObjectQuery((PetscObject)jac->schur, "AinvB", (PetscObject *)&AinvB));
1462: if (AinvB) {
1463: PetscCall(MatGetSize(AinvB, NULL, &N));
1464: if (N > P) { // first time PCApply_FieldSplit_Schur() is called
1465: PetscMemType mtype;
1466: Mat C = NULL;
1467: PetscScalar *array;
1468: PetscInt m, M, q, Q, p;
1470: PetscCall(MatGetSize(jac->B, &M, NULL));
1471: PetscCall(MatGetLocalSize(jac->B, &m, NULL));
1472: PetscCall(MatGetSize(X, NULL, &Q));
1473: PetscCall(MatGetLocalSize(X, NULL, &q));
1474: PetscCall(MatDenseGetArrayAndMemType(AinvB, &array, &mtype));
1475: if (N != P + Q) {
1476: Mat replace;
1478: PetscCall(MatGetLocalSize(jac->B, NULL, &p));
1479: if (PetscMemTypeHost(mtype) || (!PetscDefined(HAVE_CUDA) && !PetscDefined(HAVE_HIP))) {
1480: PetscCall(PetscFree(array));
1481: PetscCall(PetscMalloc1(m * (P + Q), &array));
1482: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)jac->schur), m, PETSC_DECIDE, M, P + Q, array, &replace));
1483: }
1484: #if PetscDefined(HAVE_CUDA)
1485: else if (PetscMemTypeCUDA(mtype)) {
1486: PetscCallCUDA(cudaFree(array));
1487: PetscCallCUDA(cudaMalloc((void **)&array, sizeof(PetscScalar) * m * (P + Q)));
1488: PetscCall(MatCreateDenseCUDA(PetscObjectComm((PetscObject)jac->schur), m, PETSC_DECIDE, M, P + Q, array, &replace));
1489: }
1490: #endif
1491: #if PetscDefined(HAVE_HIP)
1492: else if (PetscMemTypeHIP(mtype)) {
1493: PetscCallHIP(hipFree(array));
1494: PetscCallHIP(hipMalloc((void **)&array, sizeof(PetscScalar) * m * (P + Q)));
1495: PetscCall(MatCreateDenseHIP(PetscObjectComm((PetscObject)jac->schur), m, PETSC_DECIDE, M, P + Q, array, &replace));
1496: }
1497: #endif
1498: PetscCall(MatHeaderReplace(AinvB, &replace));
1499: }
1500: if (PetscMemTypeHost(mtype) || (!PetscDefined(HAVE_CUDA) && !PetscDefined(HAVE_HIP))) PetscCall(MatCreateDense(PetscObjectComm((PetscObject)jac->schur), m, q, M, Q, array + m * P, &C));
1501: #if PetscDefined(HAVE_CUDA)
1502: else if (PetscMemTypeCUDA(mtype)) PetscCall(MatCreateDenseCUDA(PetscObjectComm((PetscObject)jac->schur), m, q, M, Q, array + m * P, &C));
1503: #endif
1504: #if PetscDefined(HAVE_HIP)
1505: else if (PetscMemTypeHIP(mtype)) PetscCall(MatCreateDenseHIP(PetscObjectComm((PetscObject)jac->schur), m, q, M, Q, array + m * P, &C));
1506: #endif
1507: PetscCall(MatDenseRestoreArrayAndMemType(AinvB, &array));
1508: PetscCall(MatCopy(ilinkA->X, C, SAME_NONZERO_PATTERN));
1509: PetscCall(MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user));
1510: PetscCall(KSPSetOperators(jac->kspschur, jac->schur, jac->schur_user));
1511: PetscCall(MatCopy(C, ilinkA->Y, SAME_NONZERO_PATTERN)); // retrieve solutions as last columns of the composed Mat
1512: PetscCall(MatDestroy(&C));
1513: }
1514: }
1515: }
1516: if (N == P) PetscCall(KSPMatSolve(kspLower, ilinkA->X, ilinkA->Y));
1517: PetscCall(PetscLogEventEnd(KSP_Solve_FS_L, kspLower, ilinkA->X, ilinkA->Y, NULL));
1518: PetscCall(MatMatMult(jac->C, ilinkA->Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilinkD->X));
1519: PetscCall(MatScale(ilinkD->X, -1.0));
1520: PetscCall(MatDenseScatter_Private(ilinkD->sctx, X, ilinkD->X, ADD_VALUES, SCATTER_FORWARD));
1522: PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->X, ilinkD->Y, NULL));
1523: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
1524: PetscCall(KSPMatSolve(jac->kspschur, ilinkD->X, ilinkD->Y));
1525: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
1526: PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->X, ilinkD->Y, NULL));
1527: PetscCall(MatDenseScatter_Private(ilinkD->sctx, ilinkD->Y, Y, INSERT_VALUES, SCATTER_REVERSE));
1529: if (kspUpper == kspA) {
1530: if (!AinvB) {
1531: PetscCall(MatMatMult(jac->B, ilinkD->Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilinkA->Y));
1532: PetscCall(MatAXPY(ilinkA->X, -1.0, ilinkA->Y, SAME_NONZERO_PATTERN));
1533: PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->X, ilinkA->Y, NULL));
1534: PetscCall(KSPMatSolve(kspA, ilinkA->X, ilinkA->Y));
1535: PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->X, ilinkA->Y, NULL));
1536: } else {
1537: PetscCall(MatMatMult(AinvB, ilinkD->Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilinkA->X));
1538: PetscCall(MatAXPY(ilinkA->Y, 1.0, ilinkA->X, SAME_NONZERO_PATTERN));
1539: }
1540: } else {
1541: PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->X, ilinkA->Y, NULL));
1542: PetscCall(KSPMatSolve(kspA, ilinkA->X, ilinkA->Y));
1543: PetscCall(MatMatMult(jac->B, ilinkD->Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilinkA->X));
1544: if (!ilinkA->Z) PetscCall(MatDuplicate(ilinkA->X, MAT_DO_NOT_COPY_VALUES, &ilinkA->Z));
1545: PetscCall(PetscLogEventBegin(KSP_Solve_FS_U, kspUpper, ilinkA->X, ilinkA->Z, NULL));
1546: PetscCall(KSPMatSolve(kspUpper, ilinkA->X, ilinkA->Z));
1547: PetscCall(PetscLogEventEnd(KSP_Solve_FS_U, kspUpper, ilinkA->X, ilinkA->Z, NULL));
1548: PetscCall(MatAXPY(ilinkA->Y, -1.0, ilinkA->Z, SAME_NONZERO_PATTERN));
1549: }
1550: PetscCall(MatDenseScatter_Private(ilinkA->sctx, ilinkA->Y, Y, INSERT_VALUES, SCATTER_REVERSE));
1551: }
1552: PetscFunctionReturn(PETSC_SUCCESS);
1553: }
1555: static PetscErrorCode PCApplyTranspose_FieldSplit_Schur(PC pc, Vec x, Vec y)
1556: {
1557: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1558: PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
1559: KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;
1561: PetscFunctionBegin;
1562: switch (jac->schurfactorization) {
1563: case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
1564: /* [A00 0; 0 -S], positive definite, suitable for MINRES */
1565: PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1566: PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
1567: PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1568: PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1569: PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
1570: PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
1571: PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1572: PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1573: PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
1574: PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1575: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
1576: PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y));
1577: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
1578: PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
1579: PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1580: PetscCall(VecScale(ilinkD->y, jac->schurscale));
1581: PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1582: PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1583: PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1584: break;
1585: case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
1586: PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1587: PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1588: PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1589: PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
1590: PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
1591: PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1592: PetscCall(MatMultTranspose(jac->B, ilinkA->y, ilinkD->x));
1593: PetscCall(VecScale(ilinkD->x, -1.));
1594: PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
1595: PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1596: PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
1597: PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1598: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
1599: PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y));
1600: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
1601: PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
1602: PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1603: PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1604: PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1605: PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1606: break;
1607: case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
1608: PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
1609: PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
1610: PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1611: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
1612: PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y));
1613: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
1614: PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
1615: PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1616: PetscCall(MatMultTranspose(jac->C, ilinkD->y, ilinkA->x));
1617: PetscCall(VecScale(ilinkA->x, -1.));
1618: PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD));
1619: PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1620: PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD));
1621: PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1622: PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
1623: PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
1624: PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1625: PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1626: PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1627: PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1628: break;
1629: case PC_FIELDSPLIT_SCHUR_FACT_FULL:
1630: PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1631: PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1632: PetscCall(PetscLogEventBegin(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->y, NULL));
1633: PetscCall(KSPSolveTranspose(kspUpper, ilinkA->x, ilinkA->y));
1634: PetscCall(KSPCheckSolve(kspUpper, pc, ilinkA->y));
1635: PetscCall(PetscLogEventEnd(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->y, NULL));
1636: PetscCall(MatMultTranspose(jac->B, ilinkA->y, ilinkD->x));
1637: PetscCall(VecScale(ilinkD->x, -1.0));
1638: PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
1639: PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
1641: PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1642: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
1643: PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y));
1644: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
1645: PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
1646: PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1647: PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1649: if (kspLower == kspA) {
1650: PetscCall(MatMultTranspose(jac->C, ilinkD->y, ilinkA->y));
1651: PetscCall(VecAXPY(ilinkA->x, -1.0, ilinkA->y));
1652: PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1653: PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
1654: PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
1655: PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1656: } else {
1657: PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1658: PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
1659: PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
1660: PetscCall(MatMultTranspose(jac->C, ilinkD->y, ilinkA->x));
1661: PetscCall(PetscLogEventBegin(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->z, NULL));
1662: PetscCall(KSPSolveTranspose(kspLower, ilinkA->x, ilinkA->z));
1663: PetscCall(KSPCheckSolve(kspLower, pc, ilinkA->z));
1664: PetscCall(PetscLogEventEnd(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->z, NULL));
1665: PetscCall(VecAXPY(ilinkA->y, -1.0, ilinkA->z));
1666: }
1667: PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1668: PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1669: PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1670: }
1671: PetscFunctionReturn(PETSC_SUCCESS);
1672: }
1674: #define FieldSplitSplitSolveAdd(ilink, xx, yy) \
1675: ((PetscErrorCode)(VecScatterBegin(ilink->sctx, xx, ilink->x, INSERT_VALUES, SCATTER_FORWARD) || VecScatterEnd(ilink->sctx, xx, ilink->x, INSERT_VALUES, SCATTER_FORWARD) || PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL) || \
1676: KSPSolve(ilink->ksp, ilink->x, ilink->y) || KSPCheckSolve(ilink->ksp, pc, ilink->y) || PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL) || VecScatterBegin(ilink->sctx, ilink->y, yy, ADD_VALUES, SCATTER_REVERSE) || \
1677: VecScatterEnd(ilink->sctx, ilink->y, yy, ADD_VALUES, SCATTER_REVERSE)))
1679: static PetscErrorCode PCApply_FieldSplit(PC pc, Vec x, Vec y)
1680: {
1681: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1682: PC_FieldSplitLink ilink = jac->head;
1683: PetscInt cnt, bs;
1685: PetscFunctionBegin;
1686: if (jac->type == PC_COMPOSITE_ADDITIVE) {
1687: PetscBool matnest;
1689: PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATNEST, &matnest));
1690: if (jac->defaultsplit && !matnest) {
1691: PetscCall(VecGetBlockSize(x, &bs));
1692: PetscCheck(jac->bs <= 0 || bs == jac->bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Blocksize of x vector %" PetscInt_FMT " does not match fieldsplit blocksize %" PetscInt_FMT, bs, jac->bs);
1693: PetscCall(VecGetBlockSize(y, &bs));
1694: PetscCheck(jac->bs <= 0 || bs == jac->bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Blocksize of y vector %" PetscInt_FMT " does not match fieldsplit blocksize %" PetscInt_FMT, bs, jac->bs);
1695: PetscCall(VecStrideGatherAll(x, jac->x, INSERT_VALUES));
1696: while (ilink) {
1697: PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
1698: PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y));
1699: PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
1700: PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
1701: ilink = ilink->next;
1702: }
1703: PetscCall(VecStrideScatterAll(jac->y, y, INSERT_VALUES));
1704: } else {
1705: PetscCall(VecSet(y, 0.0));
1706: while (ilink) {
1707: PetscCall(FieldSplitSplitSolveAdd(ilink, x, y));
1708: ilink = ilink->next;
1709: }
1710: }
1711: } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) {
1712: PetscCall(VecSet(y, 0.0));
1713: /* solve on first block for first block variables */
1714: PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, INSERT_VALUES, SCATTER_FORWARD));
1715: PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, INSERT_VALUES, SCATTER_FORWARD));
1716: PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
1717: PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y));
1718: PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
1719: PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
1720: PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
1721: PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
1723: /* compute the residual only onto second block variables using first block variables */
1724: PetscCall(MatMult(jac->Afield[1], ilink->y, ilink->next->x));
1725: ilink = ilink->next;
1726: PetscCall(VecScale(ilink->x, -1.0));
1727: PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));
1728: PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));
1730: /* solve on second block variables */
1731: PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
1732: PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y));
1733: PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
1734: PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
1735: PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
1736: PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
1737: } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1738: if (!jac->w1) {
1739: PetscCall(VecDuplicate(x, &jac->w1));
1740: PetscCall(VecDuplicate(x, &jac->w2));
1741: }
1742: PetscCall(VecSet(y, 0.0));
1743: PetscCall(FieldSplitSplitSolveAdd(ilink, x, y));
1744: cnt = 1;
1745: while (ilink->next) {
1746: ilink = ilink->next;
1747: /* compute the residual only over the part of the vector needed */
1748: PetscCall(MatMult(jac->Afield[cnt++], y, ilink->x));
1749: PetscCall(VecScale(ilink->x, -1.0));
1750: PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));
1751: PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));
1752: PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
1753: PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y));
1754: PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
1755: PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
1756: PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
1757: PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
1758: }
1759: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1760: cnt -= 2;
1761: while (ilink->previous) {
1762: ilink = ilink->previous;
1763: /* compute the residual only over the part of the vector needed */
1764: PetscCall(MatMult(jac->Afield[cnt--], y, ilink->x));
1765: PetscCall(VecScale(ilink->x, -1.0));
1766: PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));
1767: PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));
1768: PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
1769: PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y));
1770: PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
1771: PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
1772: PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
1773: PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
1774: }
1775: }
1776: } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Unsupported or unknown composition %d", (int)jac->type);
1777: PetscFunctionReturn(PETSC_SUCCESS);
1778: }
1780: static PetscErrorCode PCMatApply_FieldSplit(PC pc, Mat X, Mat Y)
1781: {
1782: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1783: PC_FieldSplitLink ilink = jac->head;
1784: PetscInt cnt;
1786: PetscFunctionBegin;
1787: /* create working matrices with the correct number of columns */
1788: PetscCall(PCFieldSplitCreateWorkMats_Private(pc, X));
1789: if (jac->type == PC_COMPOSITE_ADDITIVE) {
1790: PetscCall(MatZeroEntries(Y));
1791: while (ilink) {
1792: PetscCall(MatDenseScatter_Private(ilink->sctx, X, ilink->X, INSERT_VALUES, SCATTER_FORWARD));
1793: PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
1794: PetscCall(KSPMatSolve(ilink->ksp, ilink->X, ilink->Y));
1795: PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
1796: PetscCall(MatDenseScatter_Private(ilink->sctx, ilink->Y, Y, ADD_VALUES, SCATTER_REVERSE));
1797: ilink = ilink->next;
1798: }
1799: } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) {
1800: PetscCall(MatZeroEntries(Y));
1801: PetscCall(MatDenseScatter_Private(ilink->sctx, X, ilink->X, INSERT_VALUES, SCATTER_FORWARD));
1802: PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
1803: PetscCall(KSPMatSolve(ilink->ksp, ilink->X, ilink->Y));
1804: PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
1805: PetscCall(MatDenseScatter_Private(ilink->sctx, ilink->Y, Y, ADD_VALUES, SCATTER_REVERSE));
1807: /* compute the residual only onto second block variables using first block variables */
1808: PetscCall(MatMatMult(jac->Afield[1], ilink->Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilink->next->X));
1809: ilink = ilink->next;
1810: PetscCall(MatScale(ilink->X, -1.0));
1811: PetscCall(MatDenseScatter_Private(ilink->sctx, X, ilink->X, ADD_VALUES, SCATTER_FORWARD));
1813: /* solve on second block variables */
1814: PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
1815: PetscCall(KSPMatSolve(ilink->ksp, ilink->X, ilink->Y));
1816: PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
1817: PetscCall(MatDenseScatter_Private(ilink->sctx, ilink->Y, Y, ADD_VALUES, SCATTER_REVERSE));
1818: } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1819: /* general multiplicative with any number of splits */
1820: PetscCall(MatZeroEntries(Y));
1821: /* first split */
1822: PetscCall(MatDenseScatter_Private(ilink->sctx, X, ilink->X, INSERT_VALUES, SCATTER_FORWARD));
1823: PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
1824: PetscCall(KSPMatSolve(ilink->ksp, ilink->X, ilink->Y));
1825: PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
1826: PetscCall(MatDenseScatter_Private(ilink->sctx, ilink->Y, Y, ADD_VALUES, SCATTER_REVERSE));
1827: cnt = 1;
1828: /* forward sweep */
1829: while (ilink->next) {
1830: ilink = ilink->next;
1831: /* compute the residual only over the part of the vector needed */
1832: PetscCall(MatMatMult(jac->Afield[cnt++], Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilink->X));
1833: PetscCall(MatScale(ilink->X, -1.0));
1834: PetscCall(MatDenseScatter_Private(ilink->sctx, X, ilink->X, ADD_VALUES, SCATTER_FORWARD));
1835: PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
1836: PetscCall(KSPMatSolve(ilink->ksp, ilink->X, ilink->Y));
1837: PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
1838: PetscCall(MatDenseScatter_Private(ilink->sctx, ilink->Y, Y, ADD_VALUES, SCATTER_REVERSE));
1839: }
1840: /* backward sweep for symmetric multiplicative */
1841: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1842: cnt -= 2;
1843: while (ilink->previous) {
1844: ilink = ilink->previous;
1845: /* compute the residual only over the part of the vector needed */
1846: PetscCall(MatMatMult(jac->Afield[cnt--], Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilink->X));
1847: PetscCall(MatScale(ilink->X, -1.0));
1848: PetscCall(MatDenseScatter_Private(ilink->sctx, X, ilink->X, ADD_VALUES, SCATTER_FORWARD));
1849: PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
1850: PetscCall(KSPMatSolve(ilink->ksp, ilink->X, ilink->Y));
1851: PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
1852: PetscCall(MatDenseScatter_Private(ilink->sctx, ilink->Y, Y, ADD_VALUES, SCATTER_REVERSE));
1853: }
1854: }
1855: } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "PCMatApply() not implemented for this fieldsplit type");
1856: PetscFunctionReturn(PETSC_SUCCESS);
1857: }
1859: static PetscErrorCode PCApply_FieldSplit_GKB(PC pc, Vec x, Vec y)
1860: {
1861: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1862: PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
1863: KSP ksp = ilinkA->ksp;
1864: Vec u, v, Hu, d, work1, work2;
1865: PetscScalar alpha, z, nrmz2, *vecz;
1866: PetscReal lowbnd, nu, beta;
1867: PetscInt j, iterGKB;
1869: PetscFunctionBegin;
1870: PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1871: PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
1872: PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1873: PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
1875: u = jac->u;
1876: v = jac->v;
1877: Hu = jac->Hu;
1878: d = jac->d;
1879: work1 = jac->w1;
1880: work2 = jac->w2;
1881: vecz = jac->vecz;
1883: /* Change RHS to comply with matrix regularization H = A + nu*B*B' */
1884: /* Add q = q + nu*B*b */
1885: if (jac->gkbnu) {
1886: nu = jac->gkbnu;
1887: PetscCall(VecScale(ilinkD->x, jac->gkbnu));
1888: PetscCall(MatMultAdd(jac->B, ilinkD->x, ilinkA->x, ilinkA->x)); /* q = q + nu*B*b */
1889: } else {
1890: /* Situation when no augmented Lagrangian is used. Then we set inner */
1891: /* matrix N = I in [Ar13], and thus nu = 1. */
1892: nu = 1;
1893: }
1895: /* Transform rhs from [q,tilde{b}] to [0,b] */
1896: PetscCall(PetscLogEventBegin(ilinkA->event, ksp, ilinkA->x, ilinkA->y, NULL));
1897: PetscCall(KSPSolve(ksp, ilinkA->x, ilinkA->y));
1898: PetscCall(KSPCheckSolve(ksp, pc, ilinkA->y));
1899: PetscCall(PetscLogEventEnd(ilinkA->event, ksp, ilinkA->x, ilinkA->y, NULL));
1900: PetscCall(MatMultHermitianTranspose(jac->B, ilinkA->y, work1));
1901: PetscCall(VecAXPBY(work1, 1.0 / nu, -1.0, ilinkD->x)); /* c = b - B'*x */
1903: /* First step of algorithm */
1904: PetscCall(VecNorm(work1, NORM_2, &beta)); /* beta = sqrt(nu*c'*c)*/
1905: KSPCheckDot(ksp, beta);
1906: beta = PetscSqrtReal(nu) * beta;
1907: PetscCall(VecAXPBY(v, nu / beta, 0.0, work1)); /* v = nu/beta *c */
1908: PetscCall(MatMult(jac->B, v, work2)); /* u = H^{-1}*B*v */
1909: PetscCall(PetscLogEventBegin(ilinkA->event, ksp, work2, u, NULL));
1910: PetscCall(KSPSolve(ksp, work2, u));
1911: PetscCall(KSPCheckSolve(ksp, pc, u));
1912: PetscCall(PetscLogEventEnd(ilinkA->event, ksp, work2, u, NULL));
1913: PetscCall(MatMult(jac->H, u, Hu)); /* alpha = u'*H*u */
1914: PetscCall(VecDot(Hu, u, &alpha));
1915: KSPCheckDot(ksp, alpha);
1916: PetscCheck(PetscRealPart(alpha) > 0.0, PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "GKB preconditioner diverged, H is not positive definite");
1917: alpha = PetscSqrtReal(PetscAbsScalar(alpha));
1918: PetscCall(VecScale(u, 1.0 / alpha));
1919: PetscCall(VecAXPBY(d, 1.0 / alpha, 0.0, v)); /* v = nu/beta *c */
1921: z = beta / alpha;
1922: vecz[1] = z;
1924: /* Computation of first iterate x(1) and p(1) */
1925: PetscCall(VecAXPY(ilinkA->y, z, u));
1926: PetscCall(VecCopy(d, ilinkD->y));
1927: PetscCall(VecScale(ilinkD->y, -z));
1929: iterGKB = 1;
1930: lowbnd = 2 * jac->gkbtol;
1931: if (jac->gkbmonitor) PetscCall(PetscViewerASCIIPrintf(jac->gkbviewer, "%3" PetscInt_FMT " GKB Lower bound estimate %14.12e\n", iterGKB, (double)lowbnd));
1933: while (iterGKB < jac->gkbmaxit && lowbnd > jac->gkbtol) {
1934: iterGKB += 1;
1935: PetscCall(MatMultHermitianTranspose(jac->B, u, work1)); /* v <- nu*(B'*u-alpha/nu*v) */
1936: PetscCall(VecAXPBY(v, nu, -alpha, work1));
1937: PetscCall(VecNorm(v, NORM_2, &beta)); /* beta = sqrt(nu)*v'*v */
1938: beta = beta / PetscSqrtReal(nu);
1939: PetscCall(VecScale(v, 1.0 / beta));
1940: PetscCall(MatMult(jac->B, v, work2)); /* u <- H^{-1}*(B*v-beta*H*u) */
1941: PetscCall(MatMult(jac->H, u, Hu));
1942: PetscCall(VecAXPY(work2, -beta, Hu));
1943: PetscCall(PetscLogEventBegin(ilinkA->event, ksp, work2, u, NULL));
1944: PetscCall(KSPSolve(ksp, work2, u));
1945: PetscCall(KSPCheckSolve(ksp, pc, u));
1946: PetscCall(PetscLogEventEnd(ilinkA->event, ksp, work2, u, NULL));
1947: PetscCall(MatMult(jac->H, u, Hu)); /* alpha = u'*H*u */
1948: PetscCall(VecDot(Hu, u, &alpha));
1949: KSPCheckDot(ksp, alpha);
1950: PetscCheck(PetscRealPart(alpha) > 0.0, PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "GKB preconditioner diverged, H is not positive definite");
1951: alpha = PetscSqrtReal(PetscAbsScalar(alpha));
1952: PetscCall(VecScale(u, 1.0 / alpha));
1954: z = -beta / alpha * z; /* z <- beta/alpha*z */
1955: vecz[0] = z;
1957: /* Computation of new iterate x(i+1) and p(i+1) */
1958: PetscCall(VecAXPBY(d, 1.0 / alpha, -beta / alpha, v)); /* d = (v-beta*d)/alpha */
1959: PetscCall(VecAXPY(ilinkA->y, z, u)); /* r = r + z*u */
1960: PetscCall(VecAXPY(ilinkD->y, -z, d)); /* p = p - z*d */
1961: PetscCall(MatMult(jac->H, ilinkA->y, Hu)); /* ||u||_H = u'*H*u */
1962: PetscCall(VecDot(Hu, ilinkA->y, &nrmz2));
1964: /* Compute Lower Bound estimate */
1965: if (iterGKB > jac->gkbdelay) {
1966: lowbnd = 0.0;
1967: for (j = 0; j < jac->gkbdelay; j++) lowbnd += PetscAbsScalar(vecz[j] * vecz[j]);
1968: lowbnd = PetscSqrtReal(lowbnd / PetscAbsScalar(nrmz2));
1969: }
1971: for (j = 0; j < jac->gkbdelay - 1; j++) vecz[jac->gkbdelay - j - 1] = vecz[jac->gkbdelay - j - 2];
1972: if (jac->gkbmonitor) PetscCall(PetscViewerASCIIPrintf(jac->gkbviewer, "%3" PetscInt_FMT " GKB Lower bound estimate %14.12e\n", iterGKB, (double)lowbnd));
1973: }
1975: PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1976: PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1977: PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1978: PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1979: PetscFunctionReturn(PETSC_SUCCESS);
1980: }
1982: #define FieldSplitSplitSolveAddTranspose(ilink, xx, yy) \
1983: ((PetscErrorCode)(VecScatterBegin(ilink->sctx, xx, ilink->y, INSERT_VALUES, SCATTER_FORWARD) || VecScatterEnd(ilink->sctx, xx, ilink->y, INSERT_VALUES, SCATTER_FORWARD) || PetscLogEventBegin(ilink->event, ilink->ksp, ilink->y, ilink->x, NULL) || \
1984: KSPSolveTranspose(ilink->ksp, ilink->y, ilink->x) || KSPCheckSolve(ilink->ksp, pc, ilink->x) || PetscLogEventEnd(ilink->event, ilink->ksp, ilink->y, ilink->x, NULL) || VecScatterBegin(ilink->sctx, ilink->x, yy, ADD_VALUES, SCATTER_REVERSE) || \
1985: VecScatterEnd(ilink->sctx, ilink->x, yy, ADD_VALUES, SCATTER_REVERSE)))
1987: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc, Vec x, Vec y)
1988: {
1989: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1990: PC_FieldSplitLink ilink = jac->head;
1991: PetscInt bs;
1993: PetscFunctionBegin;
1994: if (jac->type == PC_COMPOSITE_ADDITIVE) {
1995: PetscBool matnest;
1997: PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATNEST, &matnest));
1998: if (jac->defaultsplit && !matnest) {
1999: PetscCall(VecGetBlockSize(x, &bs));
2000: PetscCheck(jac->bs <= 0 || bs == jac->bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Blocksize of x vector %" PetscInt_FMT " does not match fieldsplit blocksize %" PetscInt_FMT, bs, jac->bs);
2001: PetscCall(VecGetBlockSize(y, &bs));
2002: PetscCheck(jac->bs <= 0 || bs == jac->bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Blocksize of y vector %" PetscInt_FMT " does not match fieldsplit blocksize %" PetscInt_FMT, bs, jac->bs);
2003: PetscCall(VecStrideGatherAll(x, jac->x, INSERT_VALUES));
2004: while (ilink) {
2005: PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
2006: PetscCall(KSPSolveTranspose(ilink->ksp, ilink->x, ilink->y));
2007: PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
2008: PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
2009: ilink = ilink->next;
2010: }
2011: PetscCall(VecStrideScatterAll(jac->y, y, INSERT_VALUES));
2012: } else {
2013: PetscCall(VecSet(y, 0.0));
2014: while (ilink) {
2015: PetscCall(FieldSplitSplitSolveAddTranspose(ilink, x, y));
2016: ilink = ilink->next;
2017: }
2018: }
2019: } else {
2020: if (!jac->w1) {
2021: PetscCall(VecDuplicate(x, &jac->w1));
2022: PetscCall(VecDuplicate(x, &jac->w2));
2023: }
2024: PetscCall(VecSet(y, 0.0));
2025: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
2026: PetscCall(FieldSplitSplitSolveAddTranspose(ilink, x, y));
2027: while (ilink->next) {
2028: ilink = ilink->next;
2029: PetscCall(MatMultTranspose(pc->mat, y, jac->w1));
2030: PetscCall(VecWAXPY(jac->w2, -1.0, jac->w1, x));
2031: PetscCall(FieldSplitSplitSolveAddTranspose(ilink, jac->w2, y));
2032: }
2033: while (ilink->previous) {
2034: ilink = ilink->previous;
2035: PetscCall(MatMultTranspose(pc->mat, y, jac->w1));
2036: PetscCall(VecWAXPY(jac->w2, -1.0, jac->w1, x));
2037: PetscCall(FieldSplitSplitSolveAddTranspose(ilink, jac->w2, y));
2038: }
2039: } else {
2040: while (ilink->next) { /* get to last entry in linked list */
2041: ilink = ilink->next;
2042: }
2043: PetscCall(FieldSplitSplitSolveAddTranspose(ilink, x, y));
2044: while (ilink->previous) {
2045: ilink = ilink->previous;
2046: PetscCall(MatMultTranspose(pc->mat, y, jac->w1));
2047: PetscCall(VecWAXPY(jac->w2, -1.0, jac->w1, x));
2048: PetscCall(FieldSplitSplitSolveAddTranspose(ilink, jac->w2, y));
2049: }
2050: }
2051: }
2052: PetscFunctionReturn(PETSC_SUCCESS);
2053: }
2055: static PetscErrorCode PCReset_FieldSplit(PC pc)
2056: {
2057: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2058: PC_FieldSplitLink ilink = jac->head, next;
2060: PetscFunctionBegin;
2061: while (ilink) {
2062: PetscCall(KSPDestroy(&ilink->ksp));
2063: PetscCall(VecDestroy(&ilink->x));
2064: PetscCall(VecDestroy(&ilink->y));
2065: PetscCall(VecDestroy(&ilink->z));
2066: PetscCall(MatDestroy(&ilink->X));
2067: PetscCall(MatDestroy(&ilink->Y));
2068: PetscCall(MatDestroy(&ilink->Z));
2069: PetscCall(VecScatterDestroy(&ilink->sctx));
2070: PetscCall(ISDestroy(&ilink->is));
2071: PetscCall(ISDestroy(&ilink->is_col));
2072: PetscCall(PetscFree(ilink->splitname));
2073: PetscCall(PetscFree(ilink->fields));
2074: PetscCall(PetscFree(ilink->fields_col));
2075: next = ilink->next;
2076: PetscCall(PetscFree(ilink));
2077: ilink = next;
2078: }
2079: jac->head = NULL;
2080: PetscCall(PetscFree2(jac->x, jac->y));
2081: if (jac->mat && jac->mat != jac->pmat) {
2082: PetscCall(MatDestroyMatrices(jac->nsplits, &jac->mat));
2083: } else if (jac->mat) {
2084: jac->mat = NULL;
2085: }
2086: if (jac->pmat) PetscCall(MatDestroyMatrices(jac->nsplits, &jac->pmat));
2087: if (jac->Afield) PetscCall(MatDestroyMatrices(jac->nsplits, &jac->Afield));
2088: jac->nsplits = 0;
2089: PetscCall(VecDestroy(&jac->w1));
2090: PetscCall(VecDestroy(&jac->w2));
2091: if (jac->schur) PetscCall(PetscObjectCompose((PetscObject)jac->schur, "AinvB", NULL));
2092: PetscCall(MatDestroy(&jac->schur));
2093: PetscCall(MatDestroy(&jac->schurp));
2094: PetscCall(MatDestroy(&jac->schur_user));
2095: PetscCall(KSPDestroy(&jac->kspschur));
2096: PetscCall(KSPDestroy(&jac->kspupper));
2097: PetscCall(MatDestroy(&jac->B));
2098: PetscCall(MatDestroy(&jac->C));
2099: PetscCall(MatDestroy(&jac->H));
2100: PetscCall(VecDestroy(&jac->u));
2101: PetscCall(VecDestroy(&jac->v));
2102: PetscCall(VecDestroy(&jac->Hu));
2103: PetscCall(VecDestroy(&jac->d));
2104: PetscCall(PetscFree(jac->vecz));
2105: PetscCall(PetscViewerDestroy(&jac->gkbviewer));
2106: jac->isrestrict = PETSC_FALSE;
2107: PetscFunctionReturn(PETSC_SUCCESS);
2108: }
2110: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
2111: {
2112: PetscFunctionBegin;
2113: PetscCall(PCReset_FieldSplit(pc));
2114: PetscCall(PetscFree(pc->data));
2115: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
2116: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetFields_C", NULL));
2117: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetIS_C", NULL));
2118: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetType_C", NULL));
2119: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetBlockSize_C", NULL));
2120: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitRestrictIS_C", NULL));
2121: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSchurGetSubKSP_C", NULL));
2122: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", NULL));
2123: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", NULL));
2124: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", NULL));
2125: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", NULL));
2126: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", NULL));
2127: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", NULL));
2128: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", NULL));
2129: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", NULL));
2130: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", NULL));
2131: PetscFunctionReturn(PETSC_SUCCESS);
2132: }
2134: static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc, PetscOptionItems PetscOptionsObject)
2135: {
2136: PetscInt bs;
2137: PetscBool flg;
2138: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2139: PCCompositeType ctype;
2141: PetscFunctionBegin;
2142: PetscOptionsHeadBegin(PetscOptionsObject, "FieldSplit options");
2143: PetscCall(PetscOptionsBool("-pc_fieldsplit_dm_splits", "Whether to use DMCreateFieldDecomposition() for splits", "PCFieldSplitSetDMSplits", jac->dm_splits, &jac->dm_splits, NULL));
2144: PetscCall(PetscOptionsInt("-pc_fieldsplit_block_size", "Blocksize that defines number of fields", "PCFieldSplitSetBlockSize", jac->bs, &bs, &flg));
2145: if (flg) PetscCall(PCFieldSplitSetBlockSize(pc, bs));
2146: jac->diag_use_amat = pc->useAmat;
2147: PetscCall(PetscOptionsBool("-pc_fieldsplit_diag_use_amat", "Use Amat (not Pmat) to extract diagonal fieldsplit blocks", "PCFieldSplitSetDiagUseAmat", jac->diag_use_amat, &jac->diag_use_amat, NULL));
2148: jac->offdiag_use_amat = pc->useAmat;
2149: PetscCall(PetscOptionsBool("-pc_fieldsplit_off_diag_use_amat", "Use Amat (not Pmat) to extract off-diagonal fieldsplit blocks", "PCFieldSplitSetOffDiagUseAmat", jac->offdiag_use_amat, &jac->offdiag_use_amat, NULL));
2150: PetscCall(PetscOptionsBool("-pc_fieldsplit_detect_saddle_point", "Form 2-way split by detecting zero diagonal entries", "PCFieldSplitSetDetectSaddlePoint", jac->detect, &jac->detect, NULL));
2151: PetscCall(PCFieldSplitSetDetectSaddlePoint(pc, jac->detect)); /* Sets split type and Schur PC type */
2152: PetscCall(PetscOptionsEnum("-pc_fieldsplit_type", "Type of composition", "PCFieldSplitSetType", PCCompositeTypes, (PetscEnum)jac->type, (PetscEnum *)&ctype, &flg));
2153: if (flg) PetscCall(PCFieldSplitSetType(pc, ctype));
2154: /* Only setup fields once */
2155: if (jac->bs > 0 && jac->nsplits == 0) {
2156: /* only allow user to set fields from command line.
2157: otherwise user can set them in PCFieldSplitSetDefaults() */
2158: PetscCall(PCFieldSplitSetRuntimeSplits_Private(pc));
2159: if (jac->splitdefined) PetscCall(PetscInfo(pc, "Splits defined using the options database\n"));
2160: }
2161: if (jac->type == PC_COMPOSITE_SCHUR) {
2162: PetscCall(PetscOptionsGetEnum(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_fieldsplit_schur_factorization_type", PCFieldSplitSchurFactTypes, (PetscEnum *)&jac->schurfactorization, &flg));
2163: if (flg) PetscCall(PetscInfo(pc, "Deprecated use of -pc_fieldsplit_schur_factorization_type\n"));
2164: PetscCall(PetscOptionsEnum("-pc_fieldsplit_schur_fact_type", "Which off-diagonal parts of the block factorization to use", "PCFieldSplitSetSchurFactType", PCFieldSplitSchurFactTypes, (PetscEnum)jac->schurfactorization, (PetscEnum *)&jac->schurfactorization, NULL));
2165: PetscCall(PetscOptionsEnum("-pc_fieldsplit_schur_precondition", "How to build preconditioner for Schur complement", "PCFieldSplitSetSchurPre", PCFieldSplitSchurPreTypes, (PetscEnum)jac->schurpre, (PetscEnum *)&jac->schurpre, NULL));
2166: PetscCall(PetscOptionsScalar("-pc_fieldsplit_schur_scale", "Scale Schur complement", "PCFieldSplitSetSchurScale", jac->schurscale, &jac->schurscale, NULL));
2167: } else if (jac->type == PC_COMPOSITE_GKB) {
2168: PetscCall(PetscOptionsReal("-pc_fieldsplit_gkb_tol", "The tolerance for the lower bound stopping criterion", "PCFieldSplitSetGKBTol", jac->gkbtol, &jac->gkbtol, NULL));
2169: PetscCall(PetscOptionsInt("-pc_fieldsplit_gkb_delay", "The delay value for lower bound criterion", "PCFieldSplitSetGKBDelay", jac->gkbdelay, &jac->gkbdelay, NULL));
2170: PetscCall(PetscOptionsBoundedReal("-pc_fieldsplit_gkb_nu", "Parameter in augmented Lagrangian approach", "PCFieldSplitSetGKBNu", jac->gkbnu, &jac->gkbnu, NULL, 0.0));
2171: PetscCall(PetscOptionsInt("-pc_fieldsplit_gkb_maxit", "Maximum allowed number of iterations", "PCFieldSplitSetGKBMaxit", jac->gkbmaxit, &jac->gkbmaxit, NULL));
2172: PetscCall(PetscOptionsBool("-pc_fieldsplit_gkb_monitor", "Prints number of GKB iterations and error", "PCFieldSplitGKB", jac->gkbmonitor, &jac->gkbmonitor, NULL));
2173: }
2174: /*
2175: In the initial call to this routine the sub-solver data structures do not exist so we cannot call KSPSetFromOptions() on them yet.
2176: But after the initial setup of ALL the layers of sub-solvers is completed we do want to call KSPSetFromOptions() on the sub-solvers every time it
2177: is called on the outer solver in case changes were made in the options database
2179: But even after PCSetUp_FieldSplit() is called all the options inside the inner levels of sub-solvers may still not have been set thus we only call the KSPSetFromOptions()
2180: if we know that the entire stack of sub-solvers below this have been complete instantiated, we check this by seeing if any solver iterations are complete.
2181: Without this extra check test p2p1fetidp_olof_full and others fail with incorrect matrix types.
2183: There could be a negative side effect of calling the KSPSetFromOptions() below.
2185: If one captured the PetscObjectState of the options database one could skip these calls if the database has not changed from the previous call
2186: */
2187: if (jac->issetup) {
2188: PC_FieldSplitLink ilink = jac->head;
2189: if (jac->type == PC_COMPOSITE_SCHUR) {
2190: if (jac->kspupper && jac->kspupper->totalits > 0) PetscCall(KSPSetFromOptions(jac->kspupper));
2191: if (jac->kspschur && jac->kspschur->totalits > 0) PetscCall(KSPSetFromOptions(jac->kspschur));
2192: }
2193: while (ilink) {
2194: if (ilink->ksp->totalits > 0) PetscCall(KSPSetFromOptions(ilink->ksp));
2195: ilink = ilink->next;
2196: }
2197: }
2198: PetscOptionsHeadEnd();
2199: PetscFunctionReturn(PETSC_SUCCESS);
2200: }
2202: static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc, const char splitname[], PetscInt n, const PetscInt *fields, const PetscInt *fields_col)
2203: {
2204: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2205: PC_FieldSplitLink ilink, next = jac->head;
2206: char prefix[128];
2207: PetscInt i;
2208: PetscLogEvent nse;
2210: PetscFunctionBegin;
2211: if (jac->splitdefined) {
2212: PetscCall(PetscInfo(pc, "Ignoring new split \"%s\" because the splits have already been defined\n", splitname));
2213: PetscFunctionReturn(PETSC_SUCCESS);
2214: }
2215: for (i = 0; i < n; i++) PetscCheck(fields[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative field %" PetscInt_FMT " requested", fields[i]);
2216: PetscCall(PetscNew(&ilink));
2217: if (splitname) {
2218: PetscCall(PetscStrallocpy(splitname, &ilink->splitname));
2219: } else {
2220: PetscCall(PetscMalloc1(3, &ilink->splitname));
2221: PetscCall(PetscSNPrintf(ilink->splitname, 2, "%" PetscInt_FMT, jac->nsplits));
2222: }
2223: PetscCall(PetscMPIIntCast(jac->nsplits, &nse));
2224: ilink->event = jac->nsplits < 5 ? KSP_Solve_FS_0 + nse : KSP_Solve_FS_0 + 4; /* Splits greater than 4 logged in 4th split */
2225: PetscCall(PetscMalloc1(n, &ilink->fields));
2226: PetscCall(PetscArraycpy(ilink->fields, fields, n));
2227: PetscCall(PetscMalloc1(n, &ilink->fields_col));
2228: PetscCall(PetscArraycpy(ilink->fields_col, fields_col, n));
2230: ilink->nfields = n;
2231: ilink->next = NULL;
2232: PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &ilink->ksp));
2233: PetscCall(KSPSetNestLevel(ilink->ksp, pc->kspnestlevel));
2234: PetscCall(KSPSetErrorIfNotConverged(ilink->ksp, pc->erroriffailure));
2235: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)pc, 1));
2236: PetscCall(KSPSetType(ilink->ksp, KSPPREONLY));
2238: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname));
2239: PetscCall(KSPSetOptionsPrefix(ilink->ksp, prefix));
2241: if (!next) {
2242: jac->head = ilink;
2243: ilink->previous = NULL;
2244: } else {
2245: while (next->next) next = next->next;
2246: next->next = ilink;
2247: ilink->previous = next;
2248: }
2249: jac->nsplits++;
2250: PetscFunctionReturn(PETSC_SUCCESS);
2251: }
2253: static PetscErrorCode PCFieldSplitSchurGetSubKSP_FieldSplit(PC pc, PetscInt *n, KSP **subksp)
2254: {
2255: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2257: PetscFunctionBegin;
2258: *subksp = NULL;
2259: if (n) *n = 0;
2260: if (jac->type == PC_COMPOSITE_SCHUR) {
2261: PetscInt nn;
2263: PetscCheck(jac->schur, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitSchurGetSubKSP()");
2264: PetscCheck(jac->nsplits == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Unexpected number of splits %" PetscInt_FMT " != 2", jac->nsplits);
2265: nn = jac->nsplits + (jac->kspupper != jac->head->ksp ? 1 : 0);
2266: PetscCall(PetscMalloc1(nn, subksp));
2267: (*subksp)[0] = jac->head->ksp;
2268: (*subksp)[1] = jac->kspschur;
2269: if (jac->kspupper != jac->head->ksp) (*subksp)[2] = jac->kspupper;
2270: if (n) *n = nn;
2271: }
2272: PetscFunctionReturn(PETSC_SUCCESS);
2273: }
2275: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc, PetscInt *n, KSP **subksp)
2276: {
2277: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2279: PetscFunctionBegin;
2280: PetscCheck(jac->schur, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitGetSubKSP()");
2281: PetscCall(PetscMalloc1(jac->nsplits, subksp));
2282: PetscCall(MatSchurComplementGetKSP(jac->schur, *subksp));
2284: (*subksp)[1] = jac->kspschur;
2285: if (n) *n = jac->nsplits;
2286: PetscFunctionReturn(PETSC_SUCCESS);
2287: }
2289: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc, PetscInt *n, KSP **subksp)
2290: {
2291: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2292: PetscInt cnt = 0;
2293: PC_FieldSplitLink ilink = jac->head;
2295: PetscFunctionBegin;
2296: PetscCall(PetscMalloc1(jac->nsplits, subksp));
2297: while (ilink) {
2298: (*subksp)[cnt++] = ilink->ksp;
2299: ilink = ilink->next;
2300: }
2301: PetscCheck(cnt == jac->nsplits, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Corrupt PCFIELDSPLIT object: number of splits in linked list %" PetscInt_FMT " does not match number in object %" PetscInt_FMT, cnt, jac->nsplits);
2302: if (n) *n = jac->nsplits;
2303: PetscFunctionReturn(PETSC_SUCCESS);
2304: }
2306: /*@
2307: PCFieldSplitRestrictIS - Restricts the fieldsplit `IS`s to be within a given `IS`.
2309: Input Parameters:
2310: + pc - the preconditioner context
2311: - isy - the index set that defines the indices to which the fieldsplit is to be restricted
2313: Level: advanced
2315: Developer Notes:
2316: It seems the resulting `IS`s will not cover the entire space, so
2317: how can they define a convergent preconditioner? Needs explaining.
2319: .seealso: [](sec_block_matrices), `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`
2320: @*/
2321: PetscErrorCode PCFieldSplitRestrictIS(PC pc, IS isy)
2322: {
2323: PetscFunctionBegin;
2326: PetscTryMethod(pc, "PCFieldSplitRestrictIS_C", (PC, IS), (pc, isy));
2327: PetscFunctionReturn(PETSC_SUCCESS);
2328: }
2330: static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy)
2331: {
2332: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2333: PC_FieldSplitLink ilink = jac->head, next;
2334: PetscInt localsize, size, sizez, i;
2335: const PetscInt *ind, *indz;
2336: PetscInt *indc, *indcz;
2337: PetscBool flg;
2339: PetscFunctionBegin;
2340: PetscCall(ISGetLocalSize(isy, &localsize));
2341: PetscCallMPI(MPI_Scan(&localsize, &size, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)isy)));
2342: size -= localsize;
2343: while (ilink) {
2344: IS isrl, isr;
2345: PC subpc;
2346: PetscCall(ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl));
2347: PetscCall(ISGetLocalSize(isrl, &localsize));
2348: PetscCall(PetscMalloc1(localsize, &indc));
2349: PetscCall(ISGetIndices(isrl, &ind));
2350: PetscCall(PetscArraycpy(indc, ind, localsize));
2351: PetscCall(ISRestoreIndices(isrl, &ind));
2352: PetscCall(ISDestroy(&isrl));
2353: for (i = 0; i < localsize; i++) *(indc + i) += size;
2354: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)isy), localsize, indc, PETSC_OWN_POINTER, &isr));
2355: PetscCall(PetscObjectReference((PetscObject)isr));
2356: PetscCall(ISDestroy(&ilink->is));
2357: ilink->is = isr;
2358: PetscCall(PetscObjectReference((PetscObject)isr));
2359: PetscCall(ISDestroy(&ilink->is_col));
2360: ilink->is_col = isr;
2361: PetscCall(ISDestroy(&isr));
2362: PetscCall(KSPGetPC(ilink->ksp, &subpc));
2363: PetscCall(PetscObjectTypeCompare((PetscObject)subpc, PCFIELDSPLIT, &flg));
2364: if (flg) {
2365: IS iszl, isz;
2366: MPI_Comm comm;
2367: PetscCall(ISGetLocalSize(ilink->is, &localsize));
2368: comm = PetscObjectComm((PetscObject)ilink->is);
2369: PetscCall(ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl));
2370: PetscCallMPI(MPI_Scan(&localsize, &sizez, 1, MPIU_INT, MPI_SUM, comm));
2371: sizez -= localsize;
2372: PetscCall(ISGetLocalSize(iszl, &localsize));
2373: PetscCall(PetscMalloc1(localsize, &indcz));
2374: PetscCall(ISGetIndices(iszl, &indz));
2375: PetscCall(PetscArraycpy(indcz, indz, localsize));
2376: PetscCall(ISRestoreIndices(iszl, &indz));
2377: PetscCall(ISDestroy(&iszl));
2378: for (i = 0; i < localsize; i++) *(indcz + i) += sizez;
2379: PetscCall(ISCreateGeneral(comm, localsize, indcz, PETSC_OWN_POINTER, &isz));
2380: PetscCall(PCFieldSplitRestrictIS(subpc, isz));
2381: PetscCall(ISDestroy(&isz));
2382: }
2383: next = ilink->next;
2384: ilink = next;
2385: }
2386: jac->isrestrict = PETSC_TRUE;
2387: PetscFunctionReturn(PETSC_SUCCESS);
2388: }
2390: static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc, const char splitname[], IS is)
2391: {
2392: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2393: PC_FieldSplitLink ilink, next = jac->head;
2394: char prefix[128];
2395: PetscLogEvent nse;
2397: PetscFunctionBegin;
2398: if (jac->splitdefined) {
2399: PetscCall(PetscInfo(pc, "Ignoring new split \"%s\" because the splits have already been defined\n", splitname));
2400: PetscFunctionReturn(PETSC_SUCCESS);
2401: }
2402: PetscCall(PetscNew(&ilink));
2403: if (splitname) {
2404: PetscCall(PetscStrallocpy(splitname, &ilink->splitname));
2405: } else {
2406: PetscCall(PetscMalloc1(8, &ilink->splitname));
2407: PetscCall(PetscSNPrintf(ilink->splitname, 7, "%" PetscInt_FMT, jac->nsplits));
2408: }
2409: PetscCall(PetscMPIIntCast(jac->nsplits, &nse));
2410: ilink->event = jac->nsplits < 5 ? KSP_Solve_FS_0 + nse : KSP_Solve_FS_0 + 4; /* Splits greater than 4 logged in 4th split */
2411: PetscCall(PetscObjectReference((PetscObject)is));
2412: PetscCall(ISDestroy(&ilink->is));
2413: ilink->is = is;
2414: PetscCall(PetscObjectReference((PetscObject)is));
2415: PetscCall(ISDestroy(&ilink->is_col));
2416: ilink->is_col = is;
2417: ilink->next = NULL;
2418: PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &ilink->ksp));
2419: PetscCall(KSPSetNestLevel(ilink->ksp, pc->kspnestlevel));
2420: PetscCall(KSPSetErrorIfNotConverged(ilink->ksp, pc->erroriffailure));
2421: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)pc, 1));
2422: PetscCall(KSPSetType(ilink->ksp, KSPPREONLY));
2424: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname));
2425: PetscCall(KSPSetOptionsPrefix(ilink->ksp, prefix));
2427: if (!next) {
2428: jac->head = ilink;
2429: ilink->previous = NULL;
2430: } else {
2431: while (next->next) next = next->next;
2432: next->next = ilink;
2433: ilink->previous = next;
2434: }
2435: jac->nsplits++;
2436: PetscFunctionReturn(PETSC_SUCCESS);
2437: }
2439: /*@
2440: PCFieldSplitSetFields - Sets the fields that define one particular split in `PCFIELDSPLIT`
2442: Logically Collective
2444: Input Parameters:
2445: + pc - the preconditioner context
2446: . splitname - name of this split, if `NULL` the number of the split is used
2447: . n - the number of fields in this split
2448: . fields - the fields in this split
2449: - fields_col - generally the same as `fields`, if it does not match `fields` then the submatrix that is solved for this set of fields comes from an off-diagonal block
2450: of the matrix and `fields_col` provides the column indices for that block
2452: Options Database Key:
2453: . -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the `%d`'th split
2455: Level: intermediate
2457: Notes:
2458: Use `PCFieldSplitSetIS()` to set a general set of indices as a split.
2460: If the matrix used to construct the preconditioner is `MATNEST` then field i refers to the `is_row[i]` `IS` passed to `MatCreateNest()`.
2462: If the matrix used to construct the preconditioner is not `MATNEST` then
2463: `PCFieldSplitSetFields()` is for defining fields as strided blocks (based on the block size provided to the matrix with `MatSetBlockSize()` or
2464: to the `PC` with `PCFieldSplitSetBlockSize()`). For example, if the block
2465: size is three then one can define a split as 0, or 1 or 2 or 0,1 or 0,2 or 1,2 which mean
2466: 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x23x56x8.. x12x45x78x....
2467: where the numbered entries indicate what is in the split.
2469: This function is called once per split (it creates a new split each time). Solve options
2470: for this split will be available under the prefix `-fieldsplit_SPLITNAME_`.
2472: `PCFieldSplitSetIS()` does not support having a `fields_col` different from `fields`
2474: Developer Notes:
2475: This routine does not actually create the `IS` representing the split, that is delayed
2476: until `PCSetUp_FieldSplit()`, because information about the vector/matrix layouts may not be
2477: available when this routine is called.
2479: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetBlockSize()`, `PCFieldSplitSetIS()`, `PCFieldSplitRestrictIS()`,
2480: `MatSetBlockSize()`, `MatCreateNest()`
2481: @*/
2482: PetscErrorCode PCFieldSplitSetFields(PC pc, const char splitname[], PetscInt n, const PetscInt fields[], const PetscInt fields_col[])
2483: {
2484: PetscFunctionBegin;
2486: PetscAssertPointer(splitname, 2);
2487: PetscCheck(n >= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Provided number of fields %" PetscInt_FMT " in split \"%s\" not positive", n, splitname);
2488: PetscAssertPointer(fields, 4);
2489: PetscTryMethod(pc, "PCFieldSplitSetFields_C", (PC, const char[], PetscInt, const PetscInt *, const PetscInt *), (pc, splitname, n, fields, fields_col));
2490: PetscFunctionReturn(PETSC_SUCCESS);
2491: }
2493: /*@
2494: PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) to build
2495: the sub-matrices associated with each split. Where `KSPSetOperators`(ksp,Amat,Pmat) was used to supply the operators.
2497: Logically Collective
2499: Input Parameters:
2500: + pc - the preconditioner object
2501: - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from
2503: Options Database Key:
2504: . -pc_fieldsplit_diag_use_amat - use the Amat to provide the diagonal blocks
2506: Level: intermediate
2508: .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitGetDiagUseAmat()`, `PCFieldSplitSetOffDiagUseAmat()`, `PCFIELDSPLIT`
2509: @*/
2510: PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc, PetscBool flg)
2511: {
2512: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2513: PetscBool isfs;
2515: PetscFunctionBegin;
2517: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
2518: PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT);
2519: jac->diag_use_amat = flg;
2520: PetscFunctionReturn(PETSC_SUCCESS);
2521: }
2523: /*@
2524: PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) to build
2525: the sub-matrices associated with each split. Where `KSPSetOperators`(ksp,Amat,Pmat) was used to supply the operators.
2527: Logically Collective
2529: Input Parameter:
2530: . pc - the preconditioner object
2532: Output Parameter:
2533: . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from
2535: Level: intermediate
2537: .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitSetDiagUseAmat()`, `PCFieldSplitGetOffDiagUseAmat()`, `PCFIELDSPLIT`
2538: @*/
2539: PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc, PetscBool *flg)
2540: {
2541: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2542: PetscBool isfs;
2544: PetscFunctionBegin;
2546: PetscAssertPointer(flg, 2);
2547: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
2548: PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT);
2549: *flg = jac->diag_use_amat;
2550: PetscFunctionReturn(PETSC_SUCCESS);
2551: }
2553: /*@
2554: PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) to build
2555: the sub-matrices associated with each split. Where `KSPSetOperators`(ksp,Amat,Pmat) was used to supply the operators.
2557: Logically Collective
2559: Input Parameters:
2560: + pc - the preconditioner object
2561: - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from
2563: Options Database Key:
2564: . -pc_fieldsplit_off_diag_use_amat <bool> - use the Amat to extract the off-diagonal blocks
2566: Level: intermediate
2568: .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitGetOffDiagUseAmat()`, `PCFieldSplitSetDiagUseAmat()`, `PCFIELDSPLIT`
2569: @*/
2570: PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc, PetscBool flg)
2571: {
2572: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2573: PetscBool isfs;
2575: PetscFunctionBegin;
2577: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
2578: PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT);
2579: jac->offdiag_use_amat = flg;
2580: PetscFunctionReturn(PETSC_SUCCESS);
2581: }
2583: /*@
2584: PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) to build
2585: the sub-matrices associated with each split. Where `KSPSetOperators`(ksp,Amat,Pmat) was used to supply the operators.
2587: Logically Collective
2589: Input Parameter:
2590: . pc - the preconditioner object
2592: Output Parameter:
2593: . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from
2595: Level: intermediate
2597: .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitSetOffDiagUseAmat()`, `PCFieldSplitGetDiagUseAmat()`, `PCFIELDSPLIT`
2598: @*/
2599: PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc, PetscBool *flg)
2600: {
2601: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2602: PetscBool isfs;
2604: PetscFunctionBegin;
2606: PetscAssertPointer(flg, 2);
2607: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
2608: PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT);
2609: *flg = jac->offdiag_use_amat;
2610: PetscFunctionReturn(PETSC_SUCCESS);
2611: }
2613: /*@
2614: PCFieldSplitSetIS - Sets the exact elements for a split in a `PCFIELDSPLIT`
2616: Logically Collective
2618: Input Parameters:
2619: + pc - the preconditioner context
2620: . splitname - name of this split, if `NULL` the number of the split is used
2621: - is - the index set that defines the elements in this split
2623: Level: intermediate
2625: Notes:
2626: Use `PCFieldSplitSetFields()`, for splits defined by strided `IS` based on the matrix block size or the `is_rows[]` passed into `MATNEST`
2628: This function is called once per split (it creates a new split each time). Solve options
2629: for this split will be available under the prefix -fieldsplit_SPLITNAME_.
2631: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetBlockSize()`, `PCFieldSplitSetFields()`
2632: @*/
2633: PetscErrorCode PCFieldSplitSetIS(PC pc, const char splitname[], IS is)
2634: {
2635: PetscFunctionBegin;
2637: if (splitname) PetscAssertPointer(splitname, 2);
2639: PetscTryMethod(pc, "PCFieldSplitSetIS_C", (PC, const char[], IS), (pc, splitname, is));
2640: PetscFunctionReturn(PETSC_SUCCESS);
2641: }
2643: /*@
2644: PCFieldSplitGetIS - Retrieves the elements for a split as an `IS`
2646: Logically Collective
2648: Input Parameters:
2649: + pc - the preconditioner context
2650: - splitname - name of this split
2652: Output Parameter:
2653: . is - the index set that defines the elements in this split, or `NULL` if the split is not found
2655: Level: intermediate
2657: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetIS()`, `PCFieldSplitGetISByIndex()`
2658: @*/
2659: PetscErrorCode PCFieldSplitGetIS(PC pc, const char splitname[], IS *is)
2660: {
2661: PetscFunctionBegin;
2663: PetscAssertPointer(splitname, 2);
2664: PetscAssertPointer(is, 3);
2665: {
2666: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2667: PC_FieldSplitLink ilink = jac->head;
2668: PetscBool found;
2670: *is = NULL;
2671: while (ilink) {
2672: PetscCall(PetscStrcmp(ilink->splitname, splitname, &found));
2673: if (found) {
2674: *is = ilink->is;
2675: break;
2676: }
2677: ilink = ilink->next;
2678: }
2679: }
2680: PetscFunctionReturn(PETSC_SUCCESS);
2681: }
2683: /*@
2684: PCFieldSplitGetISByIndex - Retrieves the elements for a given split as an `IS`
2686: Logically Collective
2688: Input Parameters:
2689: + pc - the preconditioner context
2690: - index - index of this split
2692: Output Parameter:
2693: . is - the index set that defines the elements in this split
2695: Level: intermediate
2697: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitGetIS()`, `PCFieldSplitSetIS()`,
2699: @*/
2700: PetscErrorCode PCFieldSplitGetISByIndex(PC pc, PetscInt index, IS *is)
2701: {
2702: PetscFunctionBegin;
2703: PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative field %" PetscInt_FMT " requested", index);
2705: PetscAssertPointer(is, 3);
2706: {
2707: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2708: PC_FieldSplitLink ilink = jac->head;
2709: PetscInt i = 0;
2710: PetscCheck(index < jac->nsplits, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field %" PetscInt_FMT " requested but only %" PetscInt_FMT " exist", index, jac->nsplits);
2712: while (i < index) {
2713: ilink = ilink->next;
2714: ++i;
2715: }
2716: PetscCall(PCFieldSplitGetIS(pc, ilink->splitname, is));
2717: }
2718: PetscFunctionReturn(PETSC_SUCCESS);
2719: }
2721: /*@
2722: PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the
2723: fieldsplit preconditioner when calling `PCFieldSplitSetFields()`. If not set the matrix block size is used.
2725: Logically Collective
2727: Input Parameters:
2728: + pc - the preconditioner context
2729: - bs - the block size
2731: Level: intermediate
2733: Note:
2734: If the matrix is a `MATNEST` then the `is_rows[]` passed to `MatCreateNest()` determines the fields.
2736: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`
2737: @*/
2738: PetscErrorCode PCFieldSplitSetBlockSize(PC pc, PetscInt bs)
2739: {
2740: PetscFunctionBegin;
2743: PetscTryMethod(pc, "PCFieldSplitSetBlockSize_C", (PC, PetscInt), (pc, bs));
2744: PetscFunctionReturn(PETSC_SUCCESS);
2745: }
2747: /*@C
2748: PCFieldSplitGetSubKSP - Gets the `KSP` contexts for all splits
2750: Collective
2752: Input Parameter:
2753: . pc - the preconditioner context
2755: Output Parameters:
2756: + n - the number of splits
2757: - subksp - the array of `KSP` contexts
2759: Level: advanced
2761: Notes:
2762: After `PCFieldSplitGetSubKSP()` the array of `KSP`s is to be freed by the user with `PetscFree()`
2763: (not the `KSP`, just the array that contains them).
2765: You must call `PCSetUp()` before calling `PCFieldSplitGetSubKSP()`.
2767: If the fieldsplit is of type `PC_COMPOSITE_SCHUR`, it returns the `KSP` object used inside the
2768: Schur complement and the `KSP` object used to iterate over the Schur complement.
2769: To access all the `KSP` objects used in `PC_COMPOSITE_SCHUR`, use `PCFieldSplitSchurGetSubKSP()`.
2771: If the fieldsplit is of type `PC_COMPOSITE_GKB`, it returns the `KSP` object used to solve the
2772: inner linear system defined by the matrix H in each loop.
2774: Fortran Note:
2775: Call `PCFieldSplitRestoreSubKSP()` when the array of `KSP` is no longer needed
2777: Developer Notes:
2778: There should be a `PCFieldSplitRestoreSubKSP()` instead of requiring the user to call `PetscFree()`
2780: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`, `PCFieldSplitSchurGetSubKSP()`
2781: @*/
2782: PetscErrorCode PCFieldSplitGetSubKSP(PC pc, PetscInt *n, KSP *subksp[])
2783: {
2784: PetscFunctionBegin;
2786: if (n) PetscAssertPointer(n, 2);
2787: PetscUseMethod(pc, "PCFieldSplitGetSubKSP_C", (PC, PetscInt *, KSP **), (pc, n, subksp));
2788: PetscFunctionReturn(PETSC_SUCCESS);
2789: }
2791: /*@C
2792: PCFieldSplitSchurGetSubKSP - Gets the `KSP` contexts used inside the Schur complement based `PCFIELDSPLIT`
2794: Collective
2796: Input Parameter:
2797: . pc - the preconditioner context
2799: Output Parameters:
2800: + n - the number of splits
2801: - subksp - the array of `KSP` contexts
2803: Level: advanced
2805: Notes:
2806: After `PCFieldSplitSchurGetSubKSP()` the array of `KSP`s is to be freed by the user with `PetscFree()`
2807: (not the `KSP` just the array that contains them).
2809: You must call `PCSetUp()` before calling `PCFieldSplitSchurGetSubKSP()`.
2811: If the fieldsplit type is of type `PC_COMPOSITE_SCHUR`, it returns (in order)
2812: + 1 - the `KSP` used for the (1,1) block
2813: . 2 - the `KSP` used for the Schur complement (not the one used for the interior Schur solver)
2814: - 3 - the `KSP` used for the (1,1) block in the upper triangular factor (if different from that of the (1,1) block).
2816: It returns a null array if the fieldsplit is not of type `PC_COMPOSITE_SCHUR`; in this case, you should use `PCFieldSplitGetSubKSP()`.
2818: Fortran Note:
2819: Call `PCFieldSplitSchurRestoreSubKSP()` when the array of `KSP` is no longer needed
2821: Developer Notes:
2822: There should be a `PCFieldSplitRestoreSubKSP()` instead of requiring the user to call `PetscFree()`
2824: Should the functionality of `PCFieldSplitSchurGetSubKSP()` and `PCFieldSplitGetSubKSP()` be merged?
2826: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`, `PCFieldSplitGetSubKSP()`
2827: @*/
2828: PetscErrorCode PCFieldSplitSchurGetSubKSP(PC pc, PetscInt *n, KSP *subksp[])
2829: {
2830: PetscFunctionBegin;
2832: if (n) PetscAssertPointer(n, 2);
2833: PetscUseMethod(pc, "PCFieldSplitSchurGetSubKSP_C", (PC, PetscInt *, KSP **), (pc, n, subksp));
2834: PetscFunctionReturn(PETSC_SUCCESS);
2835: }
2837: /*@
2838: PCFieldSplitSetSchurPre - Indicates from what operator the preconditioner is constructed for the Schur complement.
2839: The default is the A11 matrix.
2841: Collective
2843: Input Parameters:
2844: + pc - the preconditioner context
2845: . ptype - which matrix to use for preconditioning the Schur complement: `PC_FIELDSPLIT_SCHUR_PRE_A11` (default),
2846: `PC_FIELDSPLIT_SCHUR_PRE_SELF`, `PC_FIELDSPLIT_SCHUR_PRE_USER`,
2847: `PC_FIELDSPLIT_SCHUR_PRE_SELFP`, and `PC_FIELDSPLIT_SCHUR_PRE_FULL`
2848: - pre - matrix to use for preconditioning, or `NULL`
2850: Options Database Keys:
2851: + -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is `a11`. See notes for meaning of various arguments
2852: - -fieldsplit_1_pc_type <pctype> - the preconditioner algorithm that is used to construct the preconditioner from the operator
2854: Level: intermediate
2856: Notes:
2857: If ptype is
2858: + a11 - the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner
2859: matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix
2860: . self - the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix:
2861: The only preconditioners that currently work with this symbolic representation matrix object are `PCLSC` and `PCHPDDM`
2862: . user - the preconditioner for the Schur complement is generated from the user provided matrix (pre argument
2863: to this function).
2864: . selfp - the preconditioning for the Schur complement is generated from an explicitly-assembled approximation $ Sp = A11 - A10 inv(diag(A00)) A01 $
2865: This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be
2866: lumped before extracting the diagonal using the additional option `-fieldsplit_1_mat_schur_complement_ainv_type lump`
2867: - full - the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation
2868: computed internally by `PCFIELDSPLIT` (this is expensive)
2869: useful mostly as a test that the Schur complement approach can work for your problem
2871: When solving a saddle point problem, where the A11 block is identically zero, using `a11` as the ptype only makes sense
2872: with the additional option `-fieldsplit_1_pc_type none`. Usually for saddle point problems one would use a `ptype` of `self` and
2873: `-fieldsplit_1_pc_type lsc` which uses the least squares commutator to compute a preconditioner for the Schur complement.
2875: Developer Note:
2876: The name of this function and the option `-pc_fieldsplit_schur_precondition` are inconsistent; precondition should be used everywhere.
2878: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSchurPre()`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`,
2879: `MatSchurComplementSetAinvType()`, `PCLSC`, `PCFieldSplitSetSchurFactType()`
2880: @*/
2881: PetscErrorCode PCFieldSplitSetSchurPre(PC pc, PCFieldSplitSchurPreType ptype, Mat pre)
2882: {
2883: PetscFunctionBegin;
2885: PetscTryMethod(pc, "PCFieldSplitSetSchurPre_C", (PC, PCFieldSplitSchurPreType, Mat), (pc, ptype, pre));
2886: PetscFunctionReturn(PETSC_SUCCESS);
2887: }
2889: PetscErrorCode PCFieldSplitSchurPrecondition(PC pc, PCFieldSplitSchurPreType ptype, Mat pre)
2890: {
2891: return PCFieldSplitSetSchurPre(pc, ptype, pre);
2892: } /* Deprecated name */
2894: /*@
2895: PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be
2896: preconditioned. See `PCFieldSplitSetSchurPre()` for details.
2898: Logically Collective
2900: Input Parameter:
2901: . pc - the preconditioner context
2903: Output Parameters:
2904: + ptype - which matrix to use for preconditioning the Schur complement: `PC_FIELDSPLIT_SCHUR_PRE_A11`, `PC_FIELDSPLIT_SCHUR_PRE_SELF`, `PC_FIELDSPLIT_SCHUR_PRE_USER`
2905: - pre - matrix to use for preconditioning (with `PC_FIELDSPLIT_SCHUR_PRE_USER`), or `NULL`
2907: Level: intermediate
2909: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitSetSchurPre()`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, `PCLSC`
2910: @*/
2911: PetscErrorCode PCFieldSplitGetSchurPre(PC pc, PCFieldSplitSchurPreType *ptype, Mat *pre)
2912: {
2913: PetscFunctionBegin;
2915: PetscUseMethod(pc, "PCFieldSplitGetSchurPre_C", (PC, PCFieldSplitSchurPreType *, Mat *), (pc, ptype, pre));
2916: PetscFunctionReturn(PETSC_SUCCESS);
2917: }
2919: /*@
2920: PCFieldSplitSchurGetS - extract the `MATSCHURCOMPLEMENT` object used by this `PCFIELDSPLIT` in case it needs to be configured separately
2922: Not Collective
2924: Input Parameter:
2925: . pc - the preconditioner context
2927: Output Parameter:
2928: . S - the Schur complement matrix
2930: Level: advanced
2932: Note:
2933: This matrix should not be destroyed using `MatDestroy()`; rather, use `PCFieldSplitSchurRestoreS()`.
2935: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurPre()`, `MATSCHURCOMPLEMENT`, `PCFieldSplitSchurRestoreS()`,
2936: `MatCreateSchurComplement()`, `MatSchurComplementGetKSP()`, `MatSchurComplementComputeExplicitOperator()`, `MatGetSchurComplement()`
2937: @*/
2938: PetscErrorCode PCFieldSplitSchurGetS(PC pc, Mat *S)
2939: {
2940: const char *t;
2941: PetscBool isfs;
2942: PC_FieldSplit *jac;
2944: PetscFunctionBegin;
2946: PetscCall(PetscObjectGetType((PetscObject)pc, &t));
2947: PetscCall(PetscStrcmp(t, PCFIELDSPLIT, &isfs));
2948: PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PC of type PCFIELDSPLIT, got %s instead", t);
2949: jac = (PC_FieldSplit *)pc->data;
2950: PetscCheck(jac->type == PC_COMPOSITE_SCHUR, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PCFIELDSPLIT of type SCHUR, got %d instead", jac->type);
2951: if (S) *S = jac->schur;
2952: PetscFunctionReturn(PETSC_SUCCESS);
2953: }
2955: /*@
2956: PCFieldSplitSchurRestoreS - returns the `MATSCHURCOMPLEMENT` matrix used by this `PC`
2958: Not Collective
2960: Input Parameters:
2961: + pc - the preconditioner context
2962: - S - the Schur complement matrix
2964: Level: advanced
2966: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurPre()`, `MatSchurComplement`, `PCFieldSplitSchurGetS()`
2967: @*/
2968: PetscErrorCode PCFieldSplitSchurRestoreS(PC pc, Mat *S)
2969: {
2970: const char *t;
2971: PetscBool isfs;
2972: PC_FieldSplit *jac;
2974: PetscFunctionBegin;
2976: PetscCall(PetscObjectGetType((PetscObject)pc, &t));
2977: PetscCall(PetscStrcmp(t, PCFIELDSPLIT, &isfs));
2978: PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PC of type PCFIELDSPLIT, got %s instead", t);
2979: jac = (PC_FieldSplit *)pc->data;
2980: PetscCheck(jac->type == PC_COMPOSITE_SCHUR, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PCFIELDSPLIT of type SCHUR, got %d instead", jac->type);
2981: PetscCheck(S && (*S == jac->schur), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "MatSchurComplement restored is not the same as gotten");
2982: PetscFunctionReturn(PETSC_SUCCESS);
2983: }
2985: static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc, PCFieldSplitSchurPreType ptype, Mat pre)
2986: {
2987: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2989: PetscFunctionBegin;
2990: jac->schurpre = ptype;
2991: if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
2992: PetscCall(MatDestroy(&jac->schur_user));
2993: jac->schur_user = pre;
2994: PetscCall(PetscObjectReference((PetscObject)jac->schur_user));
2995: }
2996: PetscFunctionReturn(PETSC_SUCCESS);
2997: }
2999: static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc, PCFieldSplitSchurPreType *ptype, Mat *pre)
3000: {
3001: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
3003: PetscFunctionBegin;
3004: if (ptype) *ptype = jac->schurpre;
3005: if (pre) *pre = jac->schur_user;
3006: PetscFunctionReturn(PETSC_SUCCESS);
3007: }
3009: /*@
3010: PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain in the preconditioner {cite}`murphy2000note` and {cite}`ipsen2001note`
3012: Collective
3014: Input Parameters:
3015: + pc - the preconditioner context
3016: - ftype - which blocks of factorization to retain, `PC_FIELDSPLIT_SCHUR_FACT_FULL` is default
3018: Options Database Key:
3019: . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> - default is `full`
3021: Level: intermediate
3023: Notes:
3024: The `full` factorization is
3026: ```{math}
3027: \left(\begin{array}{cc} A & B \\
3028: C & E \\
3029: \end{array}\right) =
3030: \left(\begin{array}{cc} I & 0 \\
3031: C A^{-1} & I \\
3032: \end{array}\right)
3033: \left(\begin{array}{cc} A & 0 \\
3034: 0 & S \\
3035: \end{array}\right)
3036: \left(\begin{array}{cc} I & A^{-1}B \\
3037: 0 & I \\
3038: \end{array}\right) = L D U,
3039: ```
3041: where $ S = E - C A^{-1} B $. In practice, the full factorization is applied via block triangular solves with the grouping $L(DU)$. `upper` uses $DU$, `lower` uses $LD$,
3042: and `diag` is the diagonal part with the sign of $S$ flipped (because this makes the preconditioner positive definite for many formulations,
3043: thus allowing the use of `KSPMINRES)`. Sign flipping of $S$ can be turned off with `PCFieldSplitSetSchurScale()`.
3045: If $A$ and $S$ are solved exactly
3046: + 1 - `full` factorization is a direct solver.
3047: . 2 - The preconditioned operator with `lower` or `upper` has all eigenvalues equal to 1 and minimal polynomial of degree 2, so `KSPGMRES` converges in 2 iterations.
3048: - 3 - With `diag`, the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so `KSPGMRES` converges in at most 4 iterations.
3050: If the iteration count is very low, consider using `KSPFGMRES` or `KSPGCR` which can use one less preconditioner
3051: application in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice.
3053: For symmetric problems in which $A$ is positive definite and $S$ is negative definite, `diag` can be used with `KSPMINRES`.
3055: A flexible method like `KSPFGMRES` or `KSPGCR`, [](sec_flexibleksp), must be used if the fieldsplit preconditioner is nonlinear (e.g., a few iterations of a Krylov method is used to solve with $A$ or $S$).
3057: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurScale()`,
3058: [](sec_flexibleksp), `PCFieldSplitSetSchurPre()`
3059: @*/
3060: PetscErrorCode PCFieldSplitSetSchurFactType(PC pc, PCFieldSplitSchurFactType ftype)
3061: {
3062: PetscFunctionBegin;
3064: PetscTryMethod(pc, "PCFieldSplitSetSchurFactType_C", (PC, PCFieldSplitSchurFactType), (pc, ftype));
3065: PetscFunctionReturn(PETSC_SUCCESS);
3066: }
3068: static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc, PCFieldSplitSchurFactType ftype)
3069: {
3070: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
3072: PetscFunctionBegin;
3073: jac->schurfactorization = ftype;
3074: PetscFunctionReturn(PETSC_SUCCESS);
3075: }
3077: /*@
3078: PCFieldSplitSetSchurScale - Controls the sign flip of S for `PC_FIELDSPLIT_SCHUR_FACT_DIAG`.
3080: Collective
3082: Input Parameters:
3083: + pc - the preconditioner context
3084: - scale - scaling factor for the Schur complement
3086: Options Database Key:
3087: . -pc_fieldsplit_schur_scale <scale> - default is -1.0
3089: Level: intermediate
3091: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurFactType`, `PCFieldSplitSetSchurFactType()`
3092: @*/
3093: PetscErrorCode PCFieldSplitSetSchurScale(PC pc, PetscScalar scale)
3094: {
3095: PetscFunctionBegin;
3098: PetscTryMethod(pc, "PCFieldSplitSetSchurScale_C", (PC, PetscScalar), (pc, scale));
3099: PetscFunctionReturn(PETSC_SUCCESS);
3100: }
3102: static PetscErrorCode PCFieldSplitSetSchurScale_FieldSplit(PC pc, PetscScalar scale)
3103: {
3104: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
3106: PetscFunctionBegin;
3107: jac->schurscale = scale;
3108: PetscFunctionReturn(PETSC_SUCCESS);
3109: }
3111: /*@C
3112: PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement
3114: Collective
3116: Input Parameter:
3117: . pc - the preconditioner context
3119: Output Parameters:
3120: + A00 - the (0,0) block
3121: . A01 - the (0,1) block
3122: . A10 - the (1,0) block
3123: - A11 - the (1,1) block
3125: Level: advanced
3127: Note:
3128: Use `NULL` for any unneeded output arguments
3130: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `MatSchurComplementGetSubMatrices()`, `MatSchurComplementSetSubMatrices()`
3131: @*/
3132: PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc, Mat *A00, Mat *A01, Mat *A10, Mat *A11)
3133: {
3134: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
3136: PetscFunctionBegin;
3138: PetscCheck(jac->type == PC_COMPOSITE_SCHUR, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach.");
3139: if (A00) *A00 = jac->pmat[0];
3140: if (A01) *A01 = jac->B;
3141: if (A10) *A10 = jac->C;
3142: if (A11) *A11 = jac->pmat[1];
3143: PetscFunctionReturn(PETSC_SUCCESS);
3144: }
3146: /*@
3147: PCFieldSplitSetGKBTol - Sets the solver tolerance for the generalized Golub-Kahan bidiagonalization preconditioner {cite}`arioli2013` in `PCFIELDSPLIT`
3149: Collective
3151: Input Parameters:
3152: + pc - the preconditioner context
3153: - tolerance - the solver tolerance
3155: Options Database Key:
3156: . -pc_fieldsplit_gkb_tol <tolerance> - default is 1e-5
3158: Level: intermediate
3160: Note:
3161: The generalized GKB algorithm {cite}`arioli2013` uses a lower bound estimate of the error in energy norm as stopping criterion.
3162: It stops once the lower bound estimate undershoots the required solver tolerance. Although the actual error might be bigger than
3163: this estimate, the stopping criterion is satisfactory in practical cases.
3165: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBNu()`, `PCFieldSplitSetGKBMaxit()`
3166: @*/
3167: PetscErrorCode PCFieldSplitSetGKBTol(PC pc, PetscReal tolerance)
3168: {
3169: PetscFunctionBegin;
3172: PetscTryMethod(pc, "PCFieldSplitSetGKBTol_C", (PC, PetscReal), (pc, tolerance));
3173: PetscFunctionReturn(PETSC_SUCCESS);
3174: }
3176: static PetscErrorCode PCFieldSplitSetGKBTol_FieldSplit(PC pc, PetscReal tolerance)
3177: {
3178: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
3180: PetscFunctionBegin;
3181: jac->gkbtol = tolerance;
3182: PetscFunctionReturn(PETSC_SUCCESS);
3183: }
3185: /*@
3186: PCFieldSplitSetGKBMaxit - Sets the maximum number of iterations for the generalized Golub-Kahan bidiagonalization preconditioner {cite}`arioli2013` in `PCFIELDSPLIT`
3188: Collective
3190: Input Parameters:
3191: + pc - the preconditioner context
3192: - maxit - the maximum number of iterations
3194: Options Database Key:
3195: . -pc_fieldsplit_gkb_maxit <maxit> - default is 100
3197: Level: intermediate
3199: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBNu()`
3200: @*/
3201: PetscErrorCode PCFieldSplitSetGKBMaxit(PC pc, PetscInt maxit)
3202: {
3203: PetscFunctionBegin;
3206: PetscTryMethod(pc, "PCFieldSplitSetGKBMaxit_C", (PC, PetscInt), (pc, maxit));
3207: PetscFunctionReturn(PETSC_SUCCESS);
3208: }
3210: static PetscErrorCode PCFieldSplitSetGKBMaxit_FieldSplit(PC pc, PetscInt maxit)
3211: {
3212: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
3214: PetscFunctionBegin;
3215: jac->gkbmaxit = maxit;
3216: PetscFunctionReturn(PETSC_SUCCESS);
3217: }
3219: /*@
3220: PCFieldSplitSetGKBDelay - Sets the delay in the lower bound error estimate in the generalized Golub-Kahan bidiagonalization {cite}`arioli2013` in `PCFIELDSPLIT`
3221: preconditioner.
3223: Collective
3225: Input Parameters:
3226: + pc - the preconditioner context
3227: - delay - the delay window in the lower bound estimate
3229: Options Database Key:
3230: . -pc_fieldsplit_gkb_delay <delay> - default is 5
3232: Level: intermediate
3234: Notes:
3235: The algorithm uses a lower bound estimate of the error in energy norm as stopping criterion. The lower bound of the error $ ||u-u^k||_H $
3236: is expressed as a truncated sum. The error at iteration k can only be measured at iteration (k + `delay`), and thus the algorithm needs
3237: at least (`delay` + 1) iterations to stop.
3239: For more details on the generalized Golub-Kahan bidiagonalization method and its lower bound stopping criterion, please refer to {cite}`arioli2013`
3241: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBNu()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBMaxit()`
3242: @*/
3243: PetscErrorCode PCFieldSplitSetGKBDelay(PC pc, PetscInt delay)
3244: {
3245: PetscFunctionBegin;
3248: PetscTryMethod(pc, "PCFieldSplitSetGKBDelay_C", (PC, PetscInt), (pc, delay));
3249: PetscFunctionReturn(PETSC_SUCCESS);
3250: }
3252: static PetscErrorCode PCFieldSplitSetGKBDelay_FieldSplit(PC pc, PetscInt delay)
3253: {
3254: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
3256: PetscFunctionBegin;
3257: jac->gkbdelay = delay;
3258: PetscFunctionReturn(PETSC_SUCCESS);
3259: }
3261: /*@
3262: PCFieldSplitSetGKBNu - Sets the scalar value nu >= 0 in the transformation H = A00 + nu*A01*A01' of the (1,1) block in the
3263: Golub-Kahan bidiagonalization preconditioner {cite}`arioli2013` in `PCFIELDSPLIT`
3265: Collective
3267: Input Parameters:
3268: + pc - the preconditioner context
3269: - nu - the shift parameter
3271: Options Database Key:
3272: . -pc_fieldsplit_gkb_nu <nu> - default is 1
3274: Level: intermediate
3276: Notes:
3277: This shift is in general done to obtain better convergence properties for the outer loop of the algorithm. This is often achieved by choosing `nu` sufficiently large. However,
3278: if `nu` is chosen too large, the matrix H might be badly conditioned and the solution of the linear system $Hx = b$ in the inner loop becomes difficult. It is therefore
3279: necessary to find a good balance in between the convergence of the inner and outer loop.
3281: For `nu` = 0, no shift is done. In this case A00 has to be positive definite. The matrix N in {cite}`arioli2013` is then chosen as identity.
3283: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBMaxit()`
3284: @*/
3285: PetscErrorCode PCFieldSplitSetGKBNu(PC pc, PetscReal nu)
3286: {
3287: PetscFunctionBegin;
3290: PetscTryMethod(pc, "PCFieldSplitSetGKBNu_C", (PC, PetscReal), (pc, nu));
3291: PetscFunctionReturn(PETSC_SUCCESS);
3292: }
3294: static PetscErrorCode PCFieldSplitSetGKBNu_FieldSplit(PC pc, PetscReal nu)
3295: {
3296: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
3298: PetscFunctionBegin;
3299: jac->gkbnu = nu;
3300: PetscFunctionReturn(PETSC_SUCCESS);
3301: }
3303: static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc, PCCompositeType type)
3304: {
3305: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
3307: PetscFunctionBegin;
3308: jac->type = type;
3309: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", NULL));
3310: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", NULL));
3311: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", NULL));
3312: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", NULL));
3313: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", NULL));
3314: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", NULL));
3315: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", NULL));
3316: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", NULL));
3317: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", NULL));
3319: if (type == PC_COMPOSITE_SCHUR) {
3320: pc->ops->apply = PCApply_FieldSplit_Schur;
3321: pc->ops->applytranspose = PCApplyTranspose_FieldSplit_Schur;
3322: pc->ops->matapply = PCMatApply_FieldSplit_Schur;
3323: pc->ops->view = PCView_FieldSplit_Schur;
3324: pc->ops->setuponblocks = PCSetUpOnBlocks_FieldSplit_Schur;
3326: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit_Schur));
3327: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", PCFieldSplitSetSchurPre_FieldSplit));
3328: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", PCFieldSplitGetSchurPre_FieldSplit));
3329: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", PCFieldSplitSetSchurFactType_FieldSplit));
3330: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", PCFieldSplitSetSchurScale_FieldSplit));
3331: } else if (type == PC_COMPOSITE_GKB) {
3332: pc->ops->apply = PCApply_FieldSplit_GKB;
3333: pc->ops->applytranspose = NULL;
3334: pc->ops->matapply = NULL;
3335: pc->ops->view = PCView_FieldSplit_GKB;
3336: pc->ops->setuponblocks = PCSetUpOnBlocks_FieldSplit_GKB;
3338: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit));
3339: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", PCFieldSplitSetGKBTol_FieldSplit));
3340: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", PCFieldSplitSetGKBMaxit_FieldSplit));
3341: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", PCFieldSplitSetGKBNu_FieldSplit));
3342: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", PCFieldSplitSetGKBDelay_FieldSplit));
3343: } else {
3344: pc->ops->apply = PCApply_FieldSplit;
3345: pc->ops->applytranspose = PCApplyTranspose_FieldSplit;
3346: pc->ops->matapply = PCMatApply_FieldSplit;
3347: pc->ops->view = PCView_FieldSplit;
3348: pc->ops->setuponblocks = PCSetUpOnBlocks_FieldSplit;
3350: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit));
3351: }
3352: PetscFunctionReturn(PETSC_SUCCESS);
3353: }
3355: static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc, PetscInt bs)
3356: {
3357: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
3359: PetscFunctionBegin;
3360: PetscCheck(bs >= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Blocksize must be positive, you gave %" PetscInt_FMT, bs);
3361: PetscCheck(jac->bs <= 0 || jac->bs == bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Cannot change fieldsplit blocksize from %" PetscInt_FMT " to %" PetscInt_FMT " after it has been set", jac->bs, bs);
3362: jac->bs = bs;
3363: PetscFunctionReturn(PETSC_SUCCESS);
3364: }
3366: static PetscErrorCode PCSetCoordinates_FieldSplit(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[])
3367: {
3368: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
3369: PC_FieldSplitLink ilink_current = jac->head;
3370: IS is_owned;
3372: PetscFunctionBegin;
3373: jac->coordinates_set = PETSC_TRUE; // Internal flag
3374: PetscCall(MatGetOwnershipIS(pc->mat, &is_owned, NULL));
3376: while (ilink_current) {
3377: // For each IS, embed it to get local coords indces
3378: IS is_coords;
3379: PetscInt ndofs_block;
3380: const PetscInt *block_dofs_enumeration; // Numbering of the dofs relevant to the current block
3382: // Setting drop to true for safety. It should make no difference.
3383: PetscCall(ISEmbed(ilink_current->is, is_owned, PETSC_TRUE, &is_coords));
3384: PetscCall(ISGetLocalSize(is_coords, &ndofs_block));
3385: PetscCall(ISGetIndices(is_coords, &block_dofs_enumeration));
3387: // Allocate coordinates vector and set it directly
3388: PetscCall(PetscMalloc1(ndofs_block * dim, &ilink_current->coords));
3389: for (PetscInt dof = 0; dof < ndofs_block; ++dof) {
3390: for (PetscInt d = 0; d < dim; ++d) (ilink_current->coords)[dim * dof + d] = coords[dim * block_dofs_enumeration[dof] + d];
3391: }
3392: ilink_current->dim = dim;
3393: ilink_current->ndofs = ndofs_block;
3394: PetscCall(ISRestoreIndices(is_coords, &block_dofs_enumeration));
3395: PetscCall(ISDestroy(&is_coords));
3396: ilink_current = ilink_current->next;
3397: }
3398: PetscCall(ISDestroy(&is_owned));
3399: PetscFunctionReturn(PETSC_SUCCESS);
3400: }
3402: /*@
3403: PCFieldSplitSetType - Sets the type, `PCCompositeType`, of a `PCFIELDSPLIT`
3405: Collective
3407: Input Parameters:
3408: + pc - the preconditioner context
3409: - type - `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE` (default), `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`,
3410: `PC_COMPOSITE_GKB`
3412: Options Database Key:
3413: . -pc_fieldsplit_type <one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type
3415: Level: intermediate
3417: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCCompositeType`, `PCCompositeGetType()`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`,
3418: `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`, `PCFieldSplitSetSchurFactType()`
3419: @*/
3420: PetscErrorCode PCFieldSplitSetType(PC pc, PCCompositeType type)
3421: {
3422: PetscFunctionBegin;
3424: PetscTryMethod(pc, "PCFieldSplitSetType_C", (PC, PCCompositeType), (pc, type));
3425: PetscFunctionReturn(PETSC_SUCCESS);
3426: }
3428: /*@
3429: PCFieldSplitGetType - Gets the type, `PCCompositeType`, of a `PCFIELDSPLIT`
3431: Not collective
3433: Input Parameter:
3434: . pc - the preconditioner context
3436: Output Parameter:
3437: . type - `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE` (default), `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`
3439: Level: intermediate
3441: .seealso: [](sec_block_matrices), `PC`, `PCCompositeSetType()`, `PCFIELDSPLIT`, `PCCompositeType`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`,
3442: `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`
3443: @*/
3444: PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
3445: {
3446: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
3448: PetscFunctionBegin;
3450: PetscAssertPointer(type, 2);
3451: *type = jac->type;
3452: PetscFunctionReturn(PETSC_SUCCESS);
3453: }
3455: /*@
3456: PCFieldSplitSetDMSplits - Flags whether `DMCreateFieldDecomposition()` should be used to define the splits in a `PCFIELDSPLIT`, whenever possible.
3458: Logically Collective
3460: Input Parameters:
3461: + pc - the preconditioner context
3462: - flg - boolean indicating whether to use field splits defined by the `DM`
3464: Options Database Key:
3465: . -pc_fieldsplit_dm_splits <bool> - use the field splits defined by the `DM`
3467: Level: intermediate
3469: Developer Note:
3470: The name should be `PCFieldSplitSetUseDMSplits()`, similar change to options database
3472: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitGetDMSplits()`, `DMCreateFieldDecomposition()`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`
3473: @*/
3474: PetscErrorCode PCFieldSplitSetDMSplits(PC pc, PetscBool flg)
3475: {
3476: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
3477: PetscBool isfs;
3479: PetscFunctionBegin;
3482: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
3483: if (isfs) jac->dm_splits = flg;
3484: PetscFunctionReturn(PETSC_SUCCESS);
3485: }
3487: /*@
3488: PCFieldSplitGetDMSplits - Returns flag indicating whether `DMCreateFieldDecomposition()` should be used to define the splits in a `PCFIELDSPLIT`, whenever possible.
3490: Logically Collective
3492: Input Parameter:
3493: . pc - the preconditioner context
3495: Output Parameter:
3496: . flg - boolean indicating whether to use field splits defined by the `DM`
3498: Level: intermediate
3500: Developer Note:
3501: The name should be `PCFieldSplitGetUseDMSplits()`
3503: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetDMSplits()`, `DMCreateFieldDecomposition()`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`
3504: @*/
3505: PetscErrorCode PCFieldSplitGetDMSplits(PC pc, PetscBool *flg)
3506: {
3507: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
3508: PetscBool isfs;
3510: PetscFunctionBegin;
3512: PetscAssertPointer(flg, 2);
3513: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
3514: if (isfs) {
3515: if (flg) *flg = jac->dm_splits;
3516: }
3517: PetscFunctionReturn(PETSC_SUCCESS);
3518: }
3520: /*@
3521: PCFieldSplitGetDetectSaddlePoint - Returns flag indicating whether `PCFIELDSPLIT` will attempt to automatically determine fields based on zero diagonal entries.
3523: Logically Collective
3525: Input Parameter:
3526: . pc - the preconditioner context
3528: Output Parameter:
3529: . flg - boolean indicating whether to detect fields or not
3531: Level: intermediate
3533: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetDetectSaddlePoint()`
3534: @*/
3535: PetscErrorCode PCFieldSplitGetDetectSaddlePoint(PC pc, PetscBool *flg)
3536: {
3537: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
3539: PetscFunctionBegin;
3540: *flg = jac->detect;
3541: PetscFunctionReturn(PETSC_SUCCESS);
3542: }
3544: /*@
3545: PCFieldSplitSetDetectSaddlePoint - Sets flag indicating whether `PCFIELDSPLIT` will attempt to automatically determine fields based on zero diagonal entries.
3547: Logically Collective
3549: Input Parameter:
3550: . pc - the preconditioner context
3552: Output Parameter:
3553: . flg - boolean indicating whether to detect fields or not
3555: Options Database Key:
3556: . -pc_fieldsplit_detect_saddle_point <bool> - detect and use the saddle point
3558: Level: intermediate
3560: Note:
3561: Also sets the split type to `PC_COMPOSITE_SCHUR` (see `PCFieldSplitSetType()`) and the Schur preconditioner type to `PC_FIELDSPLIT_SCHUR_PRE_SELF` (see `PCFieldSplitSetSchurPre()`).
3563: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitGetDetectSaddlePoint()`, `PCFieldSplitSetType()`, `PCFieldSplitSetSchurPre()`, `PC_FIELDSPLIT_SCHUR_PRE_SELF`
3564: @*/
3565: PetscErrorCode PCFieldSplitSetDetectSaddlePoint(PC pc, PetscBool flg)
3566: {
3567: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
3569: PetscFunctionBegin;
3570: jac->detect = flg;
3571: if (jac->detect) {
3572: PetscCall(PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR));
3573: PetscCall(PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_SELF, NULL));
3574: }
3575: PetscFunctionReturn(PETSC_SUCCESS);
3576: }
3578: /*MC
3579: PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
3580: collections of variables (that may overlap) called fields or splits. Each field often represents a different continuum variable
3581: represented on a grid, such as velocity, pressure, or temperature.
3582: In the literature these are sometimes called block preconditioners; but should not be confused with `PCBJACOBI`.
3583: See [the users manual section on "Solving Block Matrices"](sec_block_matrices) for more details.
3585: Options Database Keys:
3586: + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the `%d`'th split
3587: . -pc_fieldsplit_default - automatically add any fields to additional splits that have not
3588: been supplied explicitly by `-pc_fieldsplit_%d_fields`
3589: . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields)
3590: when the matrix is not of `MatType` `MATNEST`
3591: . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur,gkb> - type of relaxation or factorization splitting
3592: . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is `a11`; see `PCFieldSplitSetSchurPre()`
3593: . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> - set factorization type when using `-pc_fieldsplit_type schur`;
3594: see `PCFieldSplitSetSchurFactType()`
3595: . -pc_fieldsplit_dm_splits <true,false> (default is true) - Whether to use `DMCreateFieldDecomposition()` for splits
3596: - -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero diagonal and uses Schur complement with no preconditioner as the solver
3598: Options prefixes for inner solvers when using the Schur complement preconditioner are `-fieldsplit_0_` and `-fieldsplit_1_` .
3599: The options prefix for the inner solver when using the Golub-Kahan biadiagonalization preconditioner is `-fieldsplit_0_`
3600: For all other solvers they are `-fieldsplit_%d_` for the `%d`'th field; use `-fieldsplit_` for all fields.
3602: To set options on the solvers for all blocks, prepend `-fieldsplit_` to all the `PC`
3603: options database keys. For example, `-fieldsplit_pc_type ilu` `-fieldsplit_pc_factor_levels 1`.
3605: To set the options on the solvers separate for each block call `PCFieldSplitGetSubKSP()`
3606: and set the options directly on the resulting `KSP` object
3608: Level: intermediate
3610: Notes:
3611: Use `PCFieldSplitSetFields()` to set splits defined by "strided" entries or with a `MATNEST` and `PCFieldSplitSetIS()`
3612: to define a split by an arbitrary collection of entries.
3614: If no splits are set, the default is used. If a `DM` is associated with the `PC` and it supports
3615: `DMCreateFieldDecomposition()`, then that is used for the default. Otherwise if the matrix is not `MATNEST`, the splits are defined by entries strided by bs,
3616: beginning at 0 then 1, etc to bs-1. The block size can be set with `PCFieldSplitSetBlockSize()`,
3617: if this is not called the block size defaults to the blocksize of the second matrix passed
3618: to `KSPSetOperators()`/`PCSetOperators()`.
3620: For the Schur complement preconditioner if
3621: ```{math}
3622: J = \left[\begin{array}{cc} A_{00} & A_{01} \\ A_{10} & A_{11} \end{array}\right]
3623: ```
3625: the preconditioner using `full` factorization is logically
3626: ```{math}
3627: \left[\begin{array}{cc} I & -\text{ksp}(A_{00}) A_{01} \\ 0 & I \end{array}\right] \left[\begin{array}{cc} \text{ksp}(A_{00}) & 0 \\ 0 & \text{ksp}(S) \end{array}\right] \left[\begin{array}{cc} I & 0 \\ -A_{10} \text{ksp}(A_{00}) & I \end{array}\right]
3628: ```
3629: where the action of $\text{ksp}(A_{00})$ is applied using the `KSP` solver with prefix `-fieldsplit_0_`. $S$ is the Schur complement
3630: ```{math}
3631: S = A_{11} - A_{10} \text{ksp}(A_{00}) A_{01}
3632: ```
3633: which is usually dense and not stored explicitly. The action of $\text{ksp}(S)$ is computed using the `KSP` solver with prefix `-fieldsplit_splitname_` (where `splitname`
3634: was given in providing the SECOND split or 1 if not given). Accordingly, if using `PCFieldSplitGetSubKSP()`, the array of sub-`KSP` contexts will hold two `KSP`s: at its
3635: 0th index, the `KSP` associated with `-fieldsplit_0_`, and at its 1st index, the `KSP` corresponding to `-fieldsplit_1_`.
3636: By default, $A_{11}$ is used to construct a preconditioner for $S$, use `PCFieldSplitSetSchurPre()` for all the possible ways to construct the preconditioner for $S$.
3638: The factorization type is set using `-pc_fieldsplit_schur_fact_type <diag, lower, upper, full>`. `full` is shown above,
3639: `diag` gives
3640: ```{math}
3641: \left[\begin{array}{cc} \text{ksp}(A_{00}) & 0 \\ 0 & -\text{ksp}(S) \end{array}\right]
3642: ```
3643: Note that, slightly counter intuitively, there is a negative in front of the $\text{ksp}(S)$ so that the preconditioner is positive definite. For SPD matrices $J$, the sign flip
3644: can be turned off with `PCFieldSplitSetSchurScale()` or by command line `-pc_fieldsplit_schur_scale 1.0`. The `lower` factorization is the inverse of
3645: ```{math}
3646: \left[\begin{array}{cc} A_{00} & 0 \\ A_{10} & S \end{array}\right]
3647: ```
3648: where the inverses of $A_{00}$ and $S$ are applied using `KSP`s. The upper factorization is the inverse of
3649: ```{math}
3650: \left[\begin{array}{cc} A_{00} & A_{01} \\ 0 & S \end{array}\right]
3651: ```
3652: where again the inverses of $A_{00}$ and $S$ are applied using `KSP`s.
3654: If only one set of indices (one `IS`) is provided with `PCFieldSplitSetIS()` then the complement of that `IS`
3655: is used automatically for a second submatrix.
3657: The fieldsplit preconditioner cannot currently be used with the `MATBAIJ` or `MATSBAIJ` data formats if the blocksize is larger than 1.
3658: Generally it should be used with the `MATAIJ` or `MATNEST` `MatType`
3660: The forms of these preconditioners are closely related, if not identical, to forms derived as "Distributive Iterations", see,
3661: for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling {cite}`wesseling2009`.
3662: One can also use `PCFIELDSPLIT` inside a smoother resulting in "Distributive Smoothers".
3664: See "A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations" {cite}`elman2008tcp`.
3666: The Constrained Pressure Preconditioner (CPR) can be implemented using `PCCOMPOSITE` with `PCGALERKIN`. CPR first solves an $R A P$ subsystem, updates the
3667: residual on all variables (`PCCompositeSetType(pc,PC_COMPOSITE_MULTIPLICATIVE)`), and then applies a simple ILU like preconditioner on all the variables.
3669: The generalized Golub-Kahan bidiagonalization preconditioner (GKB) can be applied to symmetric $2 \times 2$ block matrices of the shape
3670: ```{math}
3671: \left[\begin{array}{cc} A_{00} & A_{01} \\ A_{01}' & 0 \end{array}\right]
3672: ```
3673: with $A_{00}$ positive semi-definite. The implementation follows {cite}`arioli2013`. Therein, we choose $N := 1/\nu * I$ and the $(1,1)$-block of the matrix is modified to $H = _{A00} + \nu*A_{01}*A_{01}'$.
3674: A linear system $Hx = b$ has to be solved in each iteration of the GKB algorithm. This solver is chosen with the option prefix `-fieldsplit_0_`.
3676: Some `PCFIELDSPLIT` variants are called physics-based preconditioners, since the preconditioner takes into account the underlying physics of the
3677: problem. But this nomenclature is not well-defined.
3679: Developer Note:
3680: The Schur complement functionality of `PCFIELDSPLIT` should likely be factored into its own `PC` thus simplifying the implementation of the preconditioners and their
3681: user API.
3683: .seealso: [](sec_block_matrices), `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCLSC`,
3684: `PCFieldSplitGetSubKSP()`, `PCFieldSplitSchurGetSubKSP()`, `PCFieldSplitSetFields()`,
3685: `PCFieldSplitSetType()`, `PCFieldSplitSetIS()`, `PCFieldSplitSetSchurPre()`, `PCFieldSplitSetSchurFactType()`,
3686: `MatSchurComplementSetAinvType()`, `PCFieldSplitSetSchurScale()`, `PCFieldSplitSetDetectSaddlePoint()`
3687: M*/
3689: PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
3690: {
3691: PC_FieldSplit *jac;
3693: PetscFunctionBegin;
3694: PetscCall(PetscNew(&jac));
3696: jac->bs = -1;
3697: jac->type = PC_COMPOSITE_MULTIPLICATIVE;
3698: jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
3699: jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
3700: jac->schurscale = -1.0;
3701: jac->dm_splits = PETSC_TRUE;
3702: jac->gkbtol = 1e-5;
3703: jac->gkbdelay = 5;
3704: jac->gkbnu = 1;
3705: jac->gkbmaxit = 100;
3707: pc->data = (void *)jac;
3709: pc->ops->setup = PCSetUp_FieldSplit;
3710: pc->ops->reset = PCReset_FieldSplit;
3711: pc->ops->destroy = PCDestroy_FieldSplit;
3712: pc->ops->setfromoptions = PCSetFromOptions_FieldSplit;
3713: pc->ops->applyrichardson = NULL;
3715: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSchurGetSubKSP_C", PCFieldSplitSchurGetSubKSP_FieldSplit));
3716: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetFields_C", PCFieldSplitSetFields_FieldSplit));
3717: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetIS_C", PCFieldSplitSetIS_FieldSplit));
3718: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetType_C", PCFieldSplitSetType_FieldSplit));
3719: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetBlockSize_C", PCFieldSplitSetBlockSize_FieldSplit));
3720: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitRestrictIS_C", PCFieldSplitRestrictIS_FieldSplit));
3721: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_FieldSplit));
3723: /* Initialize function pointers */
3724: PetscCall(PCFieldSplitSetType(pc, jac->type));
3725: PetscFunctionReturn(PETSC_SUCCESS);
3726: }