Actual source code: fieldsplit.c
1: #include <petsc/private/pcimpl.h>
2: #include <petsc/private/kspimpl.h>
3: #include <petscdm.h>
5: const char *const PCFieldSplitSchurPreTypes[] = {"SELF", "SELFP", "A11", "USER", "FULL", "PCFieldSplitSchurPreType", "PC_FIELDSPLIT_SCHUR_PRE_", NULL};
6: const char *const PCFieldSplitSchurFactTypes[] = {"DIAG", "LOWER", "UPPER", "FULL", "PCFieldSplitSchurFactType", "PC_FIELDSPLIT_SCHUR_FACT_", NULL};
8: 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;
10: typedef struct _PC_FieldSplitLink *PC_FieldSplitLink;
11: struct _PC_FieldSplitLink {
12: KSP ksp;
13: Vec x, y, z;
14: char *splitname;
15: PetscInt nfields;
16: PetscInt *fields, *fields_col;
17: VecScatter sctx;
18: IS is, is_col;
19: PC_FieldSplitLink next, previous;
20: PetscLogEvent event;
22: /* Used only when setting coordinates with PCSetCoordinates */
23: PetscInt dim;
24: PetscInt ndofs;
25: PetscReal *coords;
26: };
28: typedef struct {
29: PCCompositeType type;
30: PetscBool defaultsplit; /* Flag for a system with a set of 'k' scalar fields with the same layout (and bs = k) */
31: PetscBool splitdefined; /* Flag is set after the splits have been defined, to prevent more splits from being added */
32: PetscInt bs; /* Block size for IS and Mat structures */
33: PetscInt nsplits; /* Number of field divisions defined */
34: Vec *x, *y, w1, w2;
35: Mat *mat; /* The diagonal block for each split */
36: Mat *pmat; /* The preconditioning diagonal block for each split */
37: Mat *Afield; /* The rows of the matrix associated with each split */
38: PetscBool issetup;
40: /* Only used when Schur complement preconditioning is used */
41: Mat B; /* The (0,1) block */
42: Mat C; /* The (1,0) block */
43: Mat schur; /* The Schur complement S = A11 - A10 A00^{-1} A01, the KSP here, kspinner, is H_1 in [El08] */
44: Mat schurp; /* Assembled approximation to S built by MatSchurComplement to be used as a preconditioning matrix when solving with S */
45: Mat schur_user; /* User-provided preconditioning matrix for the Schur complement */
46: PCFieldSplitSchurPreType schurpre; /* Determines which preconditioning matrix is used for the Schur complement */
47: PCFieldSplitSchurFactType schurfactorization;
48: KSP kspschur; /* The solver for S */
49: KSP kspupper; /* The solver for A in the upper diagonal part of the factorization (H_2 in [El08]) */
50: PetscScalar schurscale; /* Scaling factor for the Schur complement solution with DIAG factorization */
52: /* Only used when Golub-Kahan bidiagonalization preconditioning is used */
53: Mat H; /* The modified matrix H = A00 + nu*A01*A01' */
54: PetscReal gkbtol; /* Stopping tolerance for lower bound estimate */
55: PetscInt gkbdelay; /* The delay window for the stopping criterion */
56: PetscReal gkbnu; /* Parameter for augmented Lagrangian H = A + nu*A01*A01' */
57: PetscInt gkbmaxit; /* Maximum number of iterations for outer loop */
58: PetscBool gkbmonitor; /* Monitor for gkb iterations and the lower bound error */
59: PetscViewer gkbviewer; /* Viewer context for gkbmonitor */
60: Vec u, v, d, Hu; /* Work vectors for the GKB algorithm */
61: PetscScalar *vecz; /* Contains intermediate values, eg for lower bound */
63: PC_FieldSplitLink head;
64: PetscBool isrestrict; /* indicates PCFieldSplitRestrictIS() has been last called on this object, hack */
65: PetscBool suboptionsset; /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */
66: PetscBool dm_splits; /* Whether to use DMCreateFieldDecomposition() whenever possible */
67: PetscBool diag_use_amat; /* Whether to extract diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
68: PetscBool offdiag_use_amat; /* Whether to extract off-diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
69: PetscBool detect; /* Whether to form 2-way split by finding zero diagonal entries */
70: PetscBool coordinates_set; /* Whether PCSetCoordinates has been called */
71: } PC_FieldSplit;
73: /*
74: Note:
75: there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of
76: inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the
77: PC you could change this.
78: */
80: /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it. This way the
81: * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */
82: static Mat FieldSplitSchurPre(PC_FieldSplit *jac)
83: {
84: switch (jac->schurpre) {
85: case PC_FIELDSPLIT_SCHUR_PRE_SELF:
86: return jac->schur;
87: case PC_FIELDSPLIT_SCHUR_PRE_SELFP:
88: return jac->schurp;
89: case PC_FIELDSPLIT_SCHUR_PRE_A11:
90: return jac->pmat[1];
91: case PC_FIELDSPLIT_SCHUR_PRE_FULL: /* We calculate this and store it in schur_user */
92: case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */
93: default:
94: return jac->schur_user ? jac->schur_user : jac->pmat[1];
95: }
96: }
98: #include <petscdraw.h>
99: static PetscErrorCode PCView_FieldSplit(PC pc, PetscViewer viewer)
100: {
101: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
102: PetscBool iascii, isdraw;
103: PetscInt i, j;
104: PC_FieldSplitLink ilink = jac->head;
106: PetscFunctionBegin;
107: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
108: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
109: if (iascii) {
110: if (jac->bs > 0) {
111: PetscCall(PetscViewerASCIIPrintf(viewer, " FieldSplit with %s composition: total splits = %" PetscInt_FMT ", blocksize = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits, jac->bs));
112: } else {
113: PetscCall(PetscViewerASCIIPrintf(viewer, " FieldSplit with %s composition: total splits = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits));
114: }
115: if (pc->useAmat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for blocks\n"));
116: if (jac->diag_use_amat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for diagonal blocks\n"));
117: if (jac->offdiag_use_amat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for off-diagonal blocks\n"));
118: PetscCall(PetscViewerASCIIPrintf(viewer, " Solver info for each split is in the following KSP objects:\n"));
119: for (i = 0; i < jac->nsplits; i++) {
120: if (ilink->fields) {
121: PetscCall(PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Fields ", i));
122: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
123: for (j = 0; j < ilink->nfields; j++) {
124: if (j > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ","));
125: PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, ilink->fields[j]));
126: }
127: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
128: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
129: } else {
130: PetscCall(PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Defined by IS\n", i));
131: }
132: PetscCall(KSPView(ilink->ksp, viewer));
133: ilink = ilink->next;
134: }
135: }
137: if (isdraw) {
138: PetscDraw draw;
139: PetscReal x, y, w, wd;
141: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
142: PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
143: w = 2 * PetscMin(1.0 - x, x);
144: wd = w / (jac->nsplits + 1);
145: x = x - wd * (jac->nsplits - 1) / 2.0;
146: for (i = 0; i < jac->nsplits; i++) {
147: PetscCall(PetscDrawPushCurrentPoint(draw, x, y));
148: PetscCall(KSPView(ilink->ksp, viewer));
149: PetscCall(PetscDrawPopCurrentPoint(draw));
150: x += wd;
151: ilink = ilink->next;
152: }
153: }
154: PetscFunctionReturn(PETSC_SUCCESS);
155: }
157: static PetscErrorCode PCView_FieldSplit_Schur(PC pc, PetscViewer viewer)
158: {
159: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
160: PetscBool iascii, isdraw;
161: PetscInt i, j;
162: PC_FieldSplitLink ilink = jac->head;
163: MatSchurComplementAinvType atype;
165: PetscFunctionBegin;
166: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
167: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
168: if (iascii) {
169: if (jac->bs > 0) {
170: PetscCall(PetscViewerASCIIPrintf(viewer, " FieldSplit with Schur preconditioner, blocksize = %" PetscInt_FMT ", factorization %s\n", jac->bs, PCFieldSplitSchurFactTypes[jac->schurfactorization]));
171: } else {
172: PetscCall(PetscViewerASCIIPrintf(viewer, " FieldSplit with Schur preconditioner, factorization %s\n", PCFieldSplitSchurFactTypes[jac->schurfactorization]));
173: }
174: if (pc->useAmat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for blocks\n"));
175: switch (jac->schurpre) {
176: case PC_FIELDSPLIT_SCHUR_PRE_SELF:
177: PetscCall(PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from S itself\n"));
178: break;
179: case PC_FIELDSPLIT_SCHUR_PRE_SELFP:
180: if (jac->schur) {
181: PetscCall(MatSchurComplementGetAinvType(jac->schur, &atype));
182: 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 "))));
183: }
184: break;
185: case PC_FIELDSPLIT_SCHUR_PRE_A11:
186: PetscCall(PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from A11\n"));
187: break;
188: case PC_FIELDSPLIT_SCHUR_PRE_FULL:
189: PetscCall(PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from the exact Schur complement\n"));
190: break;
191: case PC_FIELDSPLIT_SCHUR_PRE_USER:
192: if (jac->schur_user) {
193: PetscCall(PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from user provided matrix\n"));
194: } else {
195: PetscCall(PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from A11\n"));
196: }
197: break;
198: default:
199: SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre);
200: }
201: PetscCall(PetscViewerASCIIPrintf(viewer, " Split info:\n"));
202: PetscCall(PetscViewerASCIIPushTab(viewer));
203: for (i = 0; i < jac->nsplits; i++) {
204: if (ilink->fields) {
205: PetscCall(PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Fields ", i));
206: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
207: for (j = 0; j < ilink->nfields; j++) {
208: if (j > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ","));
209: PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, ilink->fields[j]));
210: }
211: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
212: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
213: } else {
214: PetscCall(PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Defined by IS\n", i));
215: }
216: ilink = ilink->next;
217: }
218: PetscCall(PetscViewerASCIIPrintf(viewer, "KSP solver for A00 block\n"));
219: PetscCall(PetscViewerASCIIPushTab(viewer));
220: if (jac->head) {
221: PetscCall(KSPView(jac->head->ksp, viewer));
222: } else PetscCall(PetscViewerASCIIPrintf(viewer, " not yet available\n"));
223: PetscCall(PetscViewerASCIIPopTab(viewer));
224: if (jac->head && jac->kspupper != jac->head->ksp) {
225: PetscCall(PetscViewerASCIIPrintf(viewer, "KSP solver for upper A00 in upper triangular factor\n"));
226: PetscCall(PetscViewerASCIIPushTab(viewer));
227: if (jac->kspupper) PetscCall(KSPView(jac->kspupper, viewer));
228: else PetscCall(PetscViewerASCIIPrintf(viewer, " not yet available\n"));
229: PetscCall(PetscViewerASCIIPopTab(viewer));
230: }
231: PetscCall(PetscViewerASCIIPrintf(viewer, "KSP solver for S = A11 - A10 inv(A00) A01\n"));
232: PetscCall(PetscViewerASCIIPushTab(viewer));
233: if (jac->kspschur) {
234: PetscCall(KSPView(jac->kspschur, viewer));
235: } else {
236: PetscCall(PetscViewerASCIIPrintf(viewer, " not yet available\n"));
237: }
238: PetscCall(PetscViewerASCIIPopTab(viewer));
239: PetscCall(PetscViewerASCIIPopTab(viewer));
240: } else if (isdraw && jac->head) {
241: PetscDraw draw;
242: PetscReal x, y, w, wd, h;
243: PetscInt cnt = 2;
244: char str[32];
246: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
247: PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
248: if (jac->kspupper != jac->head->ksp) cnt++;
249: w = 2 * PetscMin(1.0 - x, x);
250: wd = w / (cnt + 1);
252: PetscCall(PetscSNPrintf(str, 32, "Schur fact. %s", PCFieldSplitSchurFactTypes[jac->schurfactorization]));
253: PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h));
254: y -= h;
255: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_USER && !jac->schur_user) {
256: PetscCall(PetscSNPrintf(str, 32, "Prec. for Schur from %s", PCFieldSplitSchurPreTypes[PC_FIELDSPLIT_SCHUR_PRE_A11]));
257: } else {
258: PetscCall(PetscSNPrintf(str, 32, "Prec. for Schur from %s", PCFieldSplitSchurPreTypes[jac->schurpre]));
259: }
260: PetscCall(PetscDrawStringBoxed(draw, x + wd * (cnt - 1) / 2.0, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h));
261: y -= h;
262: x = x - wd * (cnt - 1) / 2.0;
264: PetscCall(PetscDrawPushCurrentPoint(draw, x, y));
265: PetscCall(KSPView(jac->head->ksp, viewer));
266: PetscCall(PetscDrawPopCurrentPoint(draw));
267: if (jac->kspupper != jac->head->ksp) {
268: x += wd;
269: PetscCall(PetscDrawPushCurrentPoint(draw, x, y));
270: PetscCall(KSPView(jac->kspupper, viewer));
271: PetscCall(PetscDrawPopCurrentPoint(draw));
272: }
273: x += wd;
274: PetscCall(PetscDrawPushCurrentPoint(draw, x, y));
275: PetscCall(KSPView(jac->kspschur, viewer));
276: PetscCall(PetscDrawPopCurrentPoint(draw));
277: }
278: PetscFunctionReturn(PETSC_SUCCESS);
279: }
281: static PetscErrorCode PCView_FieldSplit_GKB(PC pc, PetscViewer viewer)
282: {
283: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
284: PetscBool iascii, isdraw;
285: PetscInt i, j;
286: PC_FieldSplitLink ilink = jac->head;
288: PetscFunctionBegin;
289: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
290: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
291: if (iascii) {
292: if (jac->bs > 0) {
293: PetscCall(PetscViewerASCIIPrintf(viewer, " FieldSplit with %s composition: total splits = %" PetscInt_FMT ", blocksize = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits, jac->bs));
294: } else {
295: PetscCall(PetscViewerASCIIPrintf(viewer, " FieldSplit with %s composition: total splits = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits));
296: }
297: if (pc->useAmat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for blocks\n"));
298: if (jac->diag_use_amat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for diagonal blocks\n"));
299: if (jac->offdiag_use_amat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for off-diagonal blocks\n"));
301: PetscCall(PetscViewerASCIIPrintf(viewer, " Stopping tolerance=%.1e, delay in error estimate=%" PetscInt_FMT ", maximum iterations=%" PetscInt_FMT "\n", (double)jac->gkbtol, jac->gkbdelay, jac->gkbmaxit));
302: PetscCall(PetscViewerASCIIPrintf(viewer, " Solver info for H = A00 + nu*A01*A01' matrix:\n"));
303: PetscCall(PetscViewerASCIIPushTab(viewer));
305: if (ilink->fields) {
306: PetscCall(PetscViewerASCIIPrintf(viewer, "Split number 0 Fields "));
307: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
308: for (j = 0; j < ilink->nfields; j++) {
309: if (j > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ","));
310: PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, ilink->fields[j]));
311: }
312: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
313: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
314: } else {
315: PetscCall(PetscViewerASCIIPrintf(viewer, "Split number 0 Defined by IS\n"));
316: }
317: PetscCall(KSPView(ilink->ksp, viewer));
319: PetscCall(PetscViewerASCIIPopTab(viewer));
320: }
322: if (isdraw) {
323: PetscDraw draw;
324: PetscReal x, y, w, wd;
326: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
327: PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
328: w = 2 * PetscMin(1.0 - x, x);
329: wd = w / (jac->nsplits + 1);
330: x = x - wd * (jac->nsplits - 1) / 2.0;
331: for (i = 0; i < jac->nsplits; i++) {
332: PetscCall(PetscDrawPushCurrentPoint(draw, x, y));
333: PetscCall(KSPView(ilink->ksp, viewer));
334: PetscCall(PetscDrawPopCurrentPoint(draw));
335: x += wd;
336: ilink = ilink->next;
337: }
338: }
339: PetscFunctionReturn(PETSC_SUCCESS);
340: }
342: /* Precondition: jac->bs is set to a meaningful value */
343: static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc)
344: {
345: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
346: PetscInt i, nfields, *ifields, nfields_col, *ifields_col;
347: PetscBool flg, flg_col;
348: char optionname[128], splitname[8], optionname_col[128];
350: PetscFunctionBegin;
351: PetscCall(PetscMalloc1(jac->bs, &ifields));
352: PetscCall(PetscMalloc1(jac->bs, &ifields_col));
353: for (i = 0, flg = PETSC_TRUE;; i++) {
354: PetscCall(PetscSNPrintf(splitname, sizeof(splitname), "%" PetscInt_FMT, i));
355: PetscCall(PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%" PetscInt_FMT "_fields", i));
356: PetscCall(PetscSNPrintf(optionname_col, sizeof(optionname_col), "-pc_fieldsplit_%" PetscInt_FMT "_fields_col", i));
357: nfields = jac->bs;
358: nfields_col = jac->bs;
359: PetscCall(PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg));
360: PetscCall(PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, optionname_col, ifields_col, &nfields_col, &flg_col));
361: if (!flg) break;
362: else if (flg && !flg_col) {
363: PetscCheck(nfields, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot list zero fields");
364: PetscCall(PCFieldSplitSetFields(pc, splitname, nfields, ifields, ifields));
365: } else {
366: PetscCheck(nfields && nfields_col, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot list zero fields");
367: PetscCheck(nfields == nfields_col, PETSC_COMM_SELF, PETSC_ERR_USER, "Number of row and column fields must match");
368: PetscCall(PCFieldSplitSetFields(pc, splitname, nfields, ifields, ifields_col));
369: }
370: }
371: if (i > 0) {
372: /* Makes command-line setting of splits take precedence over setting them in code.
373: Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would
374: create new splits, which would probably not be what the user wanted. */
375: jac->splitdefined = PETSC_TRUE;
376: }
377: PetscCall(PetscFree(ifields));
378: PetscCall(PetscFree(ifields_col));
379: PetscFunctionReturn(PETSC_SUCCESS);
380: }
382: static PetscErrorCode PCFieldSplitSetDefaults(PC pc)
383: {
384: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
385: PC_FieldSplitLink ilink = jac->head;
386: PetscBool fieldsplit_default = PETSC_FALSE, coupling = PETSC_FALSE;
387: PetscInt i;
389: PetscFunctionBegin;
390: /*
391: Kinda messy, but at least this now uses DMCreateFieldDecomposition().
392: Should probably be rewritten.
393: */
394: if (!ilink) {
395: PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_fieldsplit_detect_coupling", &coupling, NULL));
396: if (pc->dm && jac->dm_splits && !jac->detect && !coupling) {
397: PetscInt numFields, f, i, j;
398: char **fieldNames;
399: IS *fields;
400: DM *dms;
401: DM subdm[128];
402: PetscBool flg;
404: PetscCall(DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms));
405: /* Allow the user to prescribe the splits */
406: for (i = 0, flg = PETSC_TRUE;; i++) {
407: PetscInt ifields[128];
408: IS compField;
409: char optionname[128], splitname[8];
410: PetscInt nfields = numFields;
412: PetscCall(PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%" PetscInt_FMT "_fields", i));
413: PetscCall(PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg));
414: if (!flg) break;
415: PetscCheck(numFields <= 128, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot currently support %" PetscInt_FMT " > 128 fields", numFields);
416: PetscCall(DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i]));
417: if (nfields == 1) {
418: PetscCall(PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField));
419: } else {
420: PetscCall(PetscSNPrintf(splitname, sizeof(splitname), "%" PetscInt_FMT, i));
421: PetscCall(PCFieldSplitSetIS(pc, splitname, compField));
422: }
423: PetscCall(ISDestroy(&compField));
424: for (j = 0; j < nfields; ++j) {
425: f = ifields[j];
426: PetscCall(PetscFree(fieldNames[f]));
427: PetscCall(ISDestroy(&fields[f]));
428: }
429: }
430: if (i == 0) {
431: for (f = 0; f < numFields; ++f) {
432: PetscCall(PCFieldSplitSetIS(pc, fieldNames[f], fields[f]));
433: PetscCall(PetscFree(fieldNames[f]));
434: PetscCall(ISDestroy(&fields[f]));
435: }
436: } else {
437: for (j = 0; j < numFields; j++) PetscCall(DMDestroy(dms + j));
438: PetscCall(PetscFree(dms));
439: PetscCall(PetscMalloc1(i, &dms));
440: for (j = 0; j < i; ++j) dms[j] = subdm[j];
441: }
442: PetscCall(PetscFree(fieldNames));
443: PetscCall(PetscFree(fields));
444: if (dms) {
445: PetscCall(PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n"));
446: for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) {
447: const char *prefix;
448: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ilink->ksp, &prefix));
449: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)dms[i], prefix));
450: PetscCall(KSPSetDM(ilink->ksp, dms[i]));
451: PetscCall(KSPSetDMActive(ilink->ksp, PETSC_FALSE));
452: PetscCall(PetscObjectIncrementTabLevel((PetscObject)dms[i], (PetscObject)ilink->ksp, 0));
453: PetscCall(DMDestroy(&dms[i]));
454: }
455: PetscCall(PetscFree(dms));
456: }
457: } else {
458: if (jac->bs <= 0) {
459: if (pc->pmat) {
460: PetscCall(MatGetBlockSize(pc->pmat, &jac->bs));
461: } else jac->bs = 1;
462: }
464: if (jac->detect) {
465: IS zerodiags, rest;
466: PetscInt nmin, nmax;
468: PetscCall(MatGetOwnershipRange(pc->mat, &nmin, &nmax));
469: if (jac->diag_use_amat) {
470: PetscCall(MatFindZeroDiagonals(pc->mat, &zerodiags));
471: } else {
472: PetscCall(MatFindZeroDiagonals(pc->pmat, &zerodiags));
473: }
474: PetscCall(ISComplement(zerodiags, nmin, nmax, &rest));
475: PetscCall(PCFieldSplitSetIS(pc, "0", rest));
476: PetscCall(PCFieldSplitSetIS(pc, "1", zerodiags));
477: PetscCall(ISDestroy(&zerodiags));
478: PetscCall(ISDestroy(&rest));
479: } else if (coupling) {
480: IS coupling, rest;
481: PetscInt nmin, nmax;
483: PetscCall(MatGetOwnershipRange(pc->mat, &nmin, &nmax));
484: if (jac->offdiag_use_amat) {
485: PetscCall(MatFindOffBlockDiagonalEntries(pc->mat, &coupling));
486: } else {
487: PetscCall(MatFindOffBlockDiagonalEntries(pc->pmat, &coupling));
488: }
489: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc->mat), nmax - nmin, nmin, 1, &rest));
490: PetscCall(ISSetIdentity(rest));
491: PetscCall(PCFieldSplitSetIS(pc, "0", rest));
492: PetscCall(PCFieldSplitSetIS(pc, "1", coupling));
493: PetscCall(ISDestroy(&coupling));
494: PetscCall(ISDestroy(&rest));
495: } else {
496: PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_fieldsplit_default", &fieldsplit_default, NULL));
497: if (!fieldsplit_default) {
498: /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit()
499: then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */
500: PetscCall(PCFieldSplitSetRuntimeSplits_Private(pc));
501: if (jac->splitdefined) PetscCall(PetscInfo(pc, "Splits defined using the options database\n"));
502: }
503: if ((fieldsplit_default || !jac->splitdefined) && !jac->isrestrict) {
504: Mat M = pc->pmat;
505: PetscBool isnest;
507: PetscCall(PetscInfo(pc, "Using default splitting of fields\n"));
508: PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATNEST, &isnest));
509: if (!isnest) {
510: M = pc->mat;
511: PetscCall(PetscObjectTypeCompare((PetscObject)pc->mat, MATNEST, &isnest));
512: }
513: if (isnest) {
514: IS *fields;
515: PetscInt nf;
517: PetscCall(MatNestGetSize(M, &nf, NULL));
518: PetscCall(PetscMalloc1(nf, &fields));
519: PetscCall(MatNestGetISs(M, fields, NULL));
520: for (i = 0; i < nf; i++) PetscCall(PCFieldSplitSetIS(pc, NULL, fields[i]));
521: PetscCall(PetscFree(fields));
522: } else {
523: for (i = 0; i < jac->bs; i++) {
524: char splitname[8];
525: PetscCall(PetscSNPrintf(splitname, sizeof(splitname), "%" PetscInt_FMT, i));
526: PetscCall(PCFieldSplitSetFields(pc, splitname, 1, &i, &i));
527: }
528: jac->defaultsplit = PETSC_TRUE;
529: }
530: }
531: }
532: }
533: } else if (jac->nsplits == 1) {
534: IS is2;
535: PetscInt nmin, nmax;
537: PetscCheck(ilink->is, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Must provide at least two sets of fields to PCFieldSplit()");
538: PetscCall(MatGetOwnershipRange(pc->mat, &nmin, &nmax));
539: PetscCall(ISComplement(ilink->is, nmin, nmax, &is2));
540: PetscCall(PCFieldSplitSetIS(pc, "1", is2));
541: PetscCall(ISDestroy(&is2));
542: }
544: PetscCheck(jac->nsplits >= 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Unhandled case, must have at least two fields, not %" PetscInt_FMT, jac->nsplits);
545: PetscFunctionReturn(PETSC_SUCCESS);
546: }
548: static PetscErrorCode MatGolubKahanComputeExplicitOperator(Mat A, Mat B, Mat C, Mat *H, PetscReal gkbnu)
549: {
550: Mat BT, T;
551: PetscReal nrmT, nrmB;
553: PetscFunctionBegin;
554: PetscCall(MatHermitianTranspose(C, MAT_INITIAL_MATRIX, &T)); /* Test if augmented matrix is symmetric */
555: PetscCall(MatAXPY(T, -1.0, B, DIFFERENT_NONZERO_PATTERN));
556: PetscCall(MatNorm(T, NORM_1, &nrmT));
557: PetscCall(MatNorm(B, NORM_1, &nrmB));
558: PetscCheck(nrmB <= 0 || nrmT / nrmB < PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix is not symmetric/hermitian, GKB is not applicable.");
560: /* Compute augmented Lagrangian matrix H = A00 + nu*A01*A01'. This corresponds to */
561: /* setting N := 1/nu*I in [Ar13]. */
562: PetscCall(MatHermitianTranspose(B, MAT_INITIAL_MATRIX, &BT));
563: PetscCall(MatMatMult(B, BT, MAT_INITIAL_MATRIX, PETSC_DEFAULT, H)); /* H = A01*A01' */
564: PetscCall(MatAYPX(*H, gkbnu, A, DIFFERENT_NONZERO_PATTERN)); /* H = A00 + nu*A01*A01' */
566: PetscCall(MatDestroy(&BT));
567: PetscCall(MatDestroy(&T));
568: PetscFunctionReturn(PETSC_SUCCESS);
569: }
571: PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(PetscOptions, const char pre[], const char name[], const char *option[], const char *value[], PetscBool *flg);
573: static PetscErrorCode PCSetUp_FieldSplit(PC pc)
574: {
575: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
576: PC_FieldSplitLink ilink;
577: PetscInt i, nsplit;
578: PetscBool sorted, sorted_col;
580: PetscFunctionBegin;
581: pc->failedreason = PC_NOERROR;
582: PetscCall(PCFieldSplitSetDefaults(pc));
583: nsplit = jac->nsplits;
584: ilink = jac->head;
586: /* get the matrices for each split */
587: if (!jac->issetup) {
588: PetscInt rstart, rend, nslots, bs;
590: jac->issetup = PETSC_TRUE;
592: /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */
593: if (jac->defaultsplit || !ilink->is) {
594: if (jac->bs <= 0) jac->bs = nsplit;
595: }
597: /* MatCreateSubMatrix() for [S]BAIJ matrices can only work if the indices include entire blocks of the matrix */
598: PetscCall(MatGetBlockSize(pc->pmat, &bs));
599: if (bs > 1 && (jac->bs <= bs || jac->bs % bs)) {
600: PetscBool blk;
602: PetscCall(PetscObjectTypeCompareAny((PetscObject)pc->pmat, &blk, MATBAIJ, MATSBAIJ, MATSEQBAIJ, MATSEQSBAIJ, MATMPIBAIJ, MATMPISBAIJ, NULL));
603: PetscCheck(!blk, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Cannot use MATBAIJ with PCFIELDSPLIT and currently set matrix and PC blocksizes");
604: }
606: bs = jac->bs;
607: PetscCall(MatGetOwnershipRange(pc->pmat, &rstart, &rend));
608: nslots = (rend - rstart) / bs;
609: for (i = 0; i < nsplit; i++) {
610: if (jac->defaultsplit) {
611: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + i, nsplit, &ilink->is));
612: PetscCall(ISDuplicate(ilink->is, &ilink->is_col));
613: } else if (!ilink->is) {
614: if (ilink->nfields > 1) {
615: PetscInt *ii, *jj, j, k, nfields = ilink->nfields, *fields = ilink->fields, *fields_col = ilink->fields_col;
616: PetscCall(PetscMalloc1(ilink->nfields * nslots, &ii));
617: PetscCall(PetscMalloc1(ilink->nfields * nslots, &jj));
618: for (j = 0; j < nslots; j++) {
619: for (k = 0; k < nfields; k++) {
620: ii[nfields * j + k] = rstart + bs * j + fields[k];
621: jj[nfields * j + k] = rstart + bs * j + fields_col[k];
622: }
623: }
624: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), nslots * nfields, ii, PETSC_OWN_POINTER, &ilink->is));
625: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), nslots * nfields, jj, PETSC_OWN_POINTER, &ilink->is_col));
626: PetscCall(ISSetBlockSize(ilink->is, nfields));
627: PetscCall(ISSetBlockSize(ilink->is_col, nfields));
628: } else {
629: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + ilink->fields[0], bs, &ilink->is));
630: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + ilink->fields_col[0], bs, &ilink->is_col));
631: }
632: }
633: PetscCall(ISSorted(ilink->is, &sorted));
634: if (ilink->is_col) PetscCall(ISSorted(ilink->is_col, &sorted_col));
635: PetscCheck(sorted && sorted_col, PETSC_COMM_SELF, PETSC_ERR_USER, "Fields must be sorted when creating split");
636: ilink = ilink->next;
637: }
638: }
640: ilink = jac->head;
641: if (!jac->pmat) {
642: Vec xtmp;
644: PetscCall(MatCreateVecs(pc->pmat, &xtmp, NULL));
645: PetscCall(PetscMalloc1(nsplit, &jac->pmat));
646: PetscCall(PetscMalloc2(nsplit, &jac->x, nsplit, &jac->y));
647: for (i = 0; i < nsplit; i++) {
648: MatNullSpace sp;
650: /* Check for preconditioning matrix attached to IS */
651: PetscCall(PetscObjectQuery((PetscObject)ilink->is, "pmat", (PetscObject *)&jac->pmat[i]));
652: if (jac->pmat[i]) {
653: PetscCall(PetscObjectReference((PetscObject)jac->pmat[i]));
654: if (jac->type == PC_COMPOSITE_SCHUR) {
655: jac->schur_user = jac->pmat[i];
657: PetscCall(PetscObjectReference((PetscObject)jac->schur_user));
658: }
659: } else {
660: const char *prefix;
661: PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ilink->is_col, MAT_INITIAL_MATRIX, &jac->pmat[i]));
662: PetscCall(MatGetOptionsPrefix(jac->pmat[i], &prefix));
663: if (!prefix) {
664: PetscCall(KSPGetOptionsPrefix(ilink->ksp, &prefix));
665: PetscCall(MatSetOptionsPrefix(jac->pmat[i], prefix));
666: }
667: PetscCall(MatSetFromOptions(jac->pmat[i]));
668: PetscCall(MatViewFromOptions(jac->pmat[i], NULL, "-mat_view"));
669: }
670: /* create work vectors for each split */
671: PetscCall(MatCreateVecs(jac->pmat[i], &jac->x[i], &jac->y[i]));
672: ilink->x = jac->x[i];
673: ilink->y = jac->y[i];
674: ilink->z = NULL;
675: /* compute scatter contexts needed by multiplicative versions and non-default splits */
676: PetscCall(VecScatterCreate(xtmp, ilink->is, jac->x[i], NULL, &ilink->sctx));
677: PetscCall(PetscObjectQuery((PetscObject)ilink->is, "nearnullspace", (PetscObject *)&sp));
678: if (sp) PetscCall(MatSetNearNullSpace(jac->pmat[i], sp));
679: ilink = ilink->next;
680: }
681: PetscCall(VecDestroy(&xtmp));
682: } else {
683: MatReuse scall;
684: MatNullSpace *nullsp = NULL;
686: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
687: PetscCall(MatGetNullSpaces(nsplit, jac->pmat, &nullsp));
688: for (i = 0; i < nsplit; i++) PetscCall(MatDestroy(&jac->pmat[i]));
689: scall = MAT_INITIAL_MATRIX;
690: } else scall = MAT_REUSE_MATRIX;
692: for (i = 0; i < nsplit; i++) {
693: Mat pmat;
695: /* Check for preconditioning matrix attached to IS */
696: PetscCall(PetscObjectQuery((PetscObject)ilink->is, "pmat", (PetscObject *)&pmat));
697: if (!pmat) PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ilink->is_col, scall, &jac->pmat[i]));
698: ilink = ilink->next;
699: }
700: if (nullsp) PetscCall(MatRestoreNullSpaces(nsplit, jac->pmat, &nullsp));
701: }
702: if (jac->diag_use_amat) {
703: ilink = jac->head;
704: if (!jac->mat) {
705: PetscCall(PetscMalloc1(nsplit, &jac->mat));
706: for (i = 0; i < nsplit; i++) {
707: PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ilink->is_col, MAT_INITIAL_MATRIX, &jac->mat[i]));
708: ilink = ilink->next;
709: }
710: } else {
711: MatReuse scall;
712: MatNullSpace *nullsp = NULL;
714: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
715: PetscCall(MatGetNullSpaces(nsplit, jac->mat, &nullsp));
716: for (i = 0; i < nsplit; i++) PetscCall(MatDestroy(&jac->mat[i]));
717: scall = MAT_INITIAL_MATRIX;
718: } else scall = MAT_REUSE_MATRIX;
720: for (i = 0; i < nsplit; i++) {
721: PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ilink->is_col, scall, &jac->mat[i]));
722: ilink = ilink->next;
723: }
724: if (nullsp) PetscCall(MatRestoreNullSpaces(nsplit, jac->mat, &nullsp));
725: }
726: } else {
727: jac->mat = jac->pmat;
728: }
730: /* Check for null space attached to IS */
731: ilink = jac->head;
732: for (i = 0; i < nsplit; i++) {
733: MatNullSpace sp;
735: PetscCall(PetscObjectQuery((PetscObject)ilink->is, "nullspace", (PetscObject *)&sp));
736: if (sp) PetscCall(MatSetNullSpace(jac->mat[i], sp));
737: ilink = ilink->next;
738: }
740: if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR && jac->type != PC_COMPOSITE_GKB) {
741: /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */
742: /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */
743: ilink = jac->head;
744: if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) {
745: /* 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 */
746: if (!jac->Afield) {
747: PetscCall(PetscCalloc1(nsplit, &jac->Afield));
748: if (jac->offdiag_use_amat) {
749: PetscCall(MatCreateSubMatrix(pc->mat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->Afield[1]));
750: } else {
751: PetscCall(MatCreateSubMatrix(pc->pmat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->Afield[1]));
752: }
753: } else {
754: MatReuse scall;
756: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
757: PetscCall(MatDestroy(&jac->Afield[1]));
758: scall = MAT_INITIAL_MATRIX;
759: } else scall = MAT_REUSE_MATRIX;
761: if (jac->offdiag_use_amat) {
762: PetscCall(MatCreateSubMatrix(pc->mat, ilink->next->is, ilink->is, scall, &jac->Afield[1]));
763: } else {
764: PetscCall(MatCreateSubMatrix(pc->pmat, ilink->next->is, ilink->is, scall, &jac->Afield[1]));
765: }
766: }
767: } else {
768: if (!jac->Afield) {
769: PetscCall(PetscMalloc1(nsplit, &jac->Afield));
770: for (i = 0; i < nsplit; i++) {
771: if (jac->offdiag_use_amat) {
772: PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, NULL, MAT_INITIAL_MATRIX, &jac->Afield[i]));
773: } else {
774: PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, NULL, MAT_INITIAL_MATRIX, &jac->Afield[i]));
775: }
776: ilink = ilink->next;
777: }
778: } else {
779: MatReuse scall;
780: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
781: for (i = 0; i < nsplit; i++) PetscCall(MatDestroy(&jac->Afield[i]));
782: scall = MAT_INITIAL_MATRIX;
783: } else scall = MAT_REUSE_MATRIX;
785: for (i = 0; i < nsplit; i++) {
786: if (jac->offdiag_use_amat) {
787: PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, NULL, scall, &jac->Afield[i]));
788: } else {
789: PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, NULL, scall, &jac->Afield[i]));
790: }
791: ilink = ilink->next;
792: }
793: }
794: }
795: }
797: if (jac->type == PC_COMPOSITE_SCHUR) {
798: IS ccis;
799: PetscBool isset, isspd;
800: PetscInt rstart, rend;
801: char lscname[256];
802: PetscObject LSC_L;
803: PetscBool set, flg;
805: PetscCheck(nsplit == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "To use Schur complement preconditioner you must have exactly 2 fields");
807: /* If pc->mat is SPD, don't scale by -1 the Schur complement */
808: if (jac->schurscale == (PetscScalar)-1.0) {
809: PetscCall(MatIsSPDKnown(pc->pmat, &isset, &isspd));
810: jac->schurscale = (isset && isspd) ? 1.0 : -1.0;
811: }
813: /* When extracting off-diagonal submatrices, we take complements from this range */
814: PetscCall(MatGetOwnershipRangeColumn(pc->mat, &rstart, &rend));
815: PetscCall(PetscObjectTypeCompareAny(jac->offdiag_use_amat ? (PetscObject)pc->mat : (PetscObject)pc->pmat, &flg, MATSEQSBAIJ, MATMPISBAIJ, ""));
817: if (jac->schur) {
818: KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper;
819: MatReuse scall;
821: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
822: scall = MAT_INITIAL_MATRIX;
823: PetscCall(MatDestroy(&jac->B));
824: PetscCall(MatDestroy(&jac->C));
825: } else scall = MAT_REUSE_MATRIX;
827: PetscCall(MatSchurComplementGetKSP(jac->schur, &kspInner));
828: ilink = jac->head;
829: PetscCall(ISComplement(ilink->is_col, rstart, rend, &ccis));
830: if (jac->offdiag_use_amat) {
831: PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ccis, scall, &jac->B));
832: } else {
833: PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ccis, scall, &jac->B));
834: }
835: PetscCall(ISDestroy(&ccis));
836: if (!flg) {
837: ilink = ilink->next;
838: PetscCall(ISComplement(ilink->is_col, rstart, rend, &ccis));
839: if (jac->offdiag_use_amat) {
840: PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ccis, scall, &jac->C));
841: } else {
842: PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ccis, scall, &jac->C));
843: }
844: PetscCall(ISDestroy(&ccis));
845: } else {
846: PetscCall(MatIsHermitianKnown(jac->offdiag_use_amat ? pc->mat : pc->pmat, &set, &flg));
847: if (set && flg) PetscCall(MatCreateHermitianTranspose(jac->B, &jac->C));
848: else PetscCall(MatCreateTranspose(jac->B, &jac->C));
849: }
850: PetscCall(MatSchurComplementUpdateSubMatrices(jac->schur, jac->mat[0], jac->pmat[0], jac->B, jac->C, jac->mat[1]));
851: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
852: PetscCall(MatDestroy(&jac->schurp));
853: PetscCall(MatSchurComplementGetPmat(jac->schur, MAT_INITIAL_MATRIX, &jac->schurp));
854: } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) {
855: PetscCall(MatDestroy(&jac->schur_user));
856: PetscCall(MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user));
857: }
858: if (kspA != kspInner) PetscCall(KSPSetOperators(kspA, jac->mat[0], jac->pmat[0]));
859: if (kspUpper != kspA) PetscCall(KSPSetOperators(kspUpper, jac->mat[0], jac->pmat[0]));
860: PetscCall(KSPSetOperators(jac->kspschur, jac->schur, FieldSplitSchurPre(jac)));
861: } else {
862: const char *Dprefix;
863: char schurprefix[256], schurmatprefix[256];
864: char schurtestoption[256];
865: MatNullSpace sp;
866: KSP kspt;
868: /* extract the A01 and A10 matrices */
869: ilink = jac->head;
870: PetscCall(ISComplement(ilink->is_col, rstart, rend, &ccis));
871: if (jac->offdiag_use_amat) {
872: PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->B));
873: } else {
874: PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->B));
875: }
876: PetscCall(ISDestroy(&ccis));
877: ilink = ilink->next;
878: if (!flg) {
879: PetscCall(ISComplement(ilink->is_col, rstart, rend, &ccis));
880: if (jac->offdiag_use_amat) {
881: PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->C));
882: } else {
883: PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->C));
884: }
885: PetscCall(ISDestroy(&ccis));
886: } else {
887: PetscCall(MatIsHermitianKnown(jac->offdiag_use_amat ? pc->mat : pc->pmat, &set, &flg));
888: if (set && flg) PetscCall(MatCreateHermitianTranspose(jac->B, &jac->C));
889: else PetscCall(MatCreateTranspose(jac->B, &jac->C));
890: }
891: /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
892: PetscCall(MatCreate(((PetscObject)jac->mat[0])->comm, &jac->schur));
893: PetscCall(MatSetType(jac->schur, MATSCHURCOMPLEMENT));
894: PetscCall(MatSchurComplementSetSubMatrices(jac->schur, jac->mat[0], jac->pmat[0], jac->B, jac->C, jac->mat[1]));
895: PetscCall(PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname));
896: PetscCall(MatSetOptionsPrefix(jac->schur, schurmatprefix));
897: PetscCall(MatSchurComplementGetKSP(jac->schur, &kspt));
898: PetscCall(KSPSetOptionsPrefix(kspt, schurmatprefix));
900: /* Note: this is not true in general */
901: PetscCall(MatGetNullSpace(jac->mat[1], &sp));
902: if (sp) PetscCall(MatSetNullSpace(jac->schur, sp));
904: PetscCall(PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname));
905: PetscCall(PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, schurtestoption, NULL, NULL, &flg));
906: if (flg) {
907: DM dmInner;
908: KSP kspInner;
909: PC pcInner;
911: PetscCall(MatSchurComplementGetKSP(jac->schur, &kspInner));
912: PetscCall(KSPReset(kspInner));
913: PetscCall(KSPSetOperators(kspInner, jac->mat[0], jac->pmat[0]));
914: PetscCall(PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname));
915: /* Indent this deeper to emphasize the "inner" nature of this solver. */
916: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject)pc, 2));
917: PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspInner->pc, (PetscObject)pc, 2));
918: PetscCall(KSPSetOptionsPrefix(kspInner, schurprefix));
920: /* Set DM for new solver */
921: PetscCall(KSPGetDM(jac->head->ksp, &dmInner));
922: PetscCall(KSPSetDM(kspInner, dmInner));
923: PetscCall(KSPSetDMActive(kspInner, PETSC_FALSE));
925: /* Defaults to PCKSP as preconditioner */
926: PetscCall(KSPGetPC(kspInner, &pcInner));
927: PetscCall(PCSetType(pcInner, PCKSP));
928: PetscCall(PCKSPSetKSP(pcInner, jac->head->ksp));
929: } else {
930: /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or
931: * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact,
932: * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for
933: * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make
934: * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used
935: * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */
936: PetscCall(KSPSetType(jac->head->ksp, KSPGMRES));
937: PetscCall(MatSchurComplementSetKSP(jac->schur, jac->head->ksp));
938: }
939: PetscCall(KSPSetOperators(jac->head->ksp, jac->mat[0], jac->pmat[0]));
940: PetscCall(KSPSetFromOptions(jac->head->ksp));
941: PetscCall(MatSetFromOptions(jac->schur));
943: PetscCall(PetscObjectTypeCompare((PetscObject)jac->schur, MATSCHURCOMPLEMENT, &flg));
944: if (flg) { /* Need to do this otherwise PCSetUp_KSP will overwrite the amat of jac->head->ksp */
945: KSP kspInner;
946: PC pcInner;
948: PetscCall(MatSchurComplementGetKSP(jac->schur, &kspInner));
949: PetscCall(KSPGetPC(kspInner, &pcInner));
950: PetscCall(PetscObjectTypeCompare((PetscObject)pcInner, PCKSP, &flg));
951: if (flg) {
952: KSP ksp;
954: PetscCall(PCKSPGetKSP(pcInner, &ksp));
955: if (ksp == jac->head->ksp) PetscCall(PCSetUseAmat(pcInner, PETSC_TRUE));
956: }
957: }
958: PetscCall(PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname));
959: PetscCall(PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, schurtestoption, NULL, NULL, &flg));
960: if (flg) {
961: DM dmInner;
963: PetscCall(PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname));
964: PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper));
965: PetscCall(KSPSetNestLevel(jac->kspupper, pc->kspnestlevel));
966: PetscCall(KSPSetErrorIfNotConverged(jac->kspupper, pc->erroriffailure));
967: PetscCall(KSPSetOptionsPrefix(jac->kspupper, schurprefix));
968: PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->kspupper, (PetscObject)pc, 1));
969: PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->kspupper->pc, (PetscObject)pc, 1));
970: PetscCall(KSPGetDM(jac->head->ksp, &dmInner));
971: PetscCall(KSPSetDM(jac->kspupper, dmInner));
972: PetscCall(KSPSetDMActive(jac->kspupper, PETSC_FALSE));
973: PetscCall(KSPSetFromOptions(jac->kspupper));
974: PetscCall(KSPSetOperators(jac->kspupper, jac->mat[0], jac->pmat[0]));
975: PetscCall(VecDuplicate(jac->head->x, &jac->head->z));
976: } else {
977: jac->kspupper = jac->head->ksp;
978: PetscCall(PetscObjectReference((PetscObject)jac->head->ksp));
979: }
981: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) PetscCall(MatSchurComplementGetPmat(jac->schur, MAT_INITIAL_MATRIX, &jac->schurp));
982: PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspschur));
983: PetscCall(KSPSetNestLevel(jac->kspschur, pc->kspnestlevel));
984: PetscCall(KSPSetErrorIfNotConverged(jac->kspschur, pc->erroriffailure));
985: PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->kspschur, (PetscObject)pc, 1));
986: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
987: PC pcschur;
988: PetscCall(KSPGetPC(jac->kspschur, &pcschur));
989: PetscCall(PCSetType(pcschur, PCNONE));
990: /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
991: } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) {
992: PetscCall(MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user));
993: }
994: PetscCall(KSPSetOperators(jac->kspschur, jac->schur, FieldSplitSchurPre(jac)));
995: PetscCall(KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix));
996: PetscCall(KSPSetOptionsPrefix(jac->kspschur, Dprefix));
997: /* propagate DM */
998: {
999: DM sdm;
1000: PetscCall(KSPGetDM(jac->head->next->ksp, &sdm));
1001: if (sdm) {
1002: PetscCall(KSPSetDM(jac->kspschur, sdm));
1003: PetscCall(KSPSetDMActive(jac->kspschur, PETSC_FALSE));
1004: }
1005: }
1006: /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
1007: /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */
1008: PetscCall(KSPSetFromOptions(jac->kspschur));
1009: }
1010: PetscCall(MatAssemblyBegin(jac->schur, MAT_FINAL_ASSEMBLY));
1011: PetscCall(MatAssemblyEnd(jac->schur, MAT_FINAL_ASSEMBLY));
1013: /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */
1014: PetscCall(PetscSNPrintf(lscname, sizeof(lscname), "%s_LSC_L", ilink->splitname));
1015: PetscCall(PetscObjectQuery((PetscObject)pc->mat, lscname, (PetscObject *)&LSC_L));
1016: if (!LSC_L) PetscCall(PetscObjectQuery((PetscObject)pc->pmat, lscname, (PetscObject *)&LSC_L));
1017: if (LSC_L) PetscCall(PetscObjectCompose((PetscObject)jac->schur, "LSC_L", (PetscObject)LSC_L));
1018: PetscCall(PetscSNPrintf(lscname, sizeof(lscname), "%s_LSC_Lp", ilink->splitname));
1019: PetscCall(PetscObjectQuery((PetscObject)pc->pmat, lscname, (PetscObject *)&LSC_L));
1020: if (!LSC_L) PetscCall(PetscObjectQuery((PetscObject)pc->mat, lscname, (PetscObject *)&LSC_L));
1021: if (LSC_L) PetscCall(PetscObjectCompose((PetscObject)jac->schur, "LSC_Lp", (PetscObject)LSC_L));
1022: } else if (jac->type == PC_COMPOSITE_GKB) {
1023: IS ccis;
1024: PetscInt rstart, rend;
1026: PetscCheck(nsplit == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "To use GKB preconditioner you must have exactly 2 fields");
1028: ilink = jac->head;
1030: /* When extracting off-diagonal submatrices, we take complements from this range */
1031: PetscCall(MatGetOwnershipRangeColumn(pc->mat, &rstart, &rend));
1033: PetscCall(ISComplement(ilink->is_col, rstart, rend, &ccis));
1034: if (jac->offdiag_use_amat) {
1035: PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->B));
1036: } else {
1037: PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->B));
1038: }
1039: PetscCall(ISDestroy(&ccis));
1040: /* Create work vectors for GKB algorithm */
1041: PetscCall(VecDuplicate(ilink->x, &jac->u));
1042: PetscCall(VecDuplicate(ilink->x, &jac->Hu));
1043: PetscCall(VecDuplicate(ilink->x, &jac->w2));
1044: ilink = ilink->next;
1045: PetscCall(ISComplement(ilink->is_col, rstart, rend, &ccis));
1046: if (jac->offdiag_use_amat) {
1047: PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->C));
1048: } else {
1049: PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->C));
1050: }
1051: PetscCall(ISDestroy(&ccis));
1052: /* Create work vectors for GKB algorithm */
1053: PetscCall(VecDuplicate(ilink->x, &jac->v));
1054: PetscCall(VecDuplicate(ilink->x, &jac->d));
1055: PetscCall(VecDuplicate(ilink->x, &jac->w1));
1056: PetscCall(MatGolubKahanComputeExplicitOperator(jac->mat[0], jac->B, jac->C, &jac->H, jac->gkbnu));
1057: PetscCall(PetscCalloc1(jac->gkbdelay, &jac->vecz));
1059: ilink = jac->head;
1060: PetscCall(KSPSetOperators(ilink->ksp, jac->H, jac->H));
1061: if (!jac->suboptionsset) PetscCall(KSPSetFromOptions(ilink->ksp));
1062: /* Create gkb_monitor context */
1063: if (jac->gkbmonitor) {
1064: PetscInt tablevel;
1065: PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &jac->gkbviewer));
1066: PetscCall(PetscViewerSetType(jac->gkbviewer, PETSCVIEWERASCII));
1067: PetscCall(PetscObjectGetTabLevel((PetscObject)ilink->ksp, &tablevel));
1068: PetscCall(PetscViewerASCIISetTab(jac->gkbviewer, tablevel));
1069: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)ilink->ksp, 1));
1070: }
1071: } else {
1072: /* set up the individual splits' PCs */
1073: i = 0;
1074: ilink = jac->head;
1075: while (ilink) {
1076: PetscCall(KSPSetOperators(ilink->ksp, jac->mat[i], jac->pmat[i]));
1077: /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
1078: if (!jac->suboptionsset) PetscCall(KSPSetFromOptions(ilink->ksp));
1079: i++;
1080: ilink = ilink->next;
1081: }
1082: }
1084: /* Set coordinates to the sub PC objects whenever these are set */
1085: if (jac->coordinates_set) {
1086: PC pc_coords;
1087: if (jac->type == PC_COMPOSITE_SCHUR) {
1088: // Head is first block.
1089: PetscCall(KSPGetPC(jac->head->ksp, &pc_coords));
1090: PetscCall(PCSetCoordinates(pc_coords, jac->head->dim, jac->head->ndofs, jac->head->coords));
1091: // Second one is Schur block, but its KSP object is in kspschur.
1092: PetscCall(KSPGetPC(jac->kspschur, &pc_coords));
1093: PetscCall(PCSetCoordinates(pc_coords, jac->head->next->dim, jac->head->next->ndofs, jac->head->next->coords));
1094: } else if (jac->type == PC_COMPOSITE_GKB) {
1095: PetscCall(PetscInfo(pc, "Warning: Setting coordinates does nothing for the GKB Fieldpslit preconditioner\n"));
1096: } else {
1097: ilink = jac->head;
1098: while (ilink) {
1099: PetscCall(KSPGetPC(ilink->ksp, &pc_coords));
1100: PetscCall(PCSetCoordinates(pc_coords, ilink->dim, ilink->ndofs, ilink->coords));
1101: ilink = ilink->next;
1102: }
1103: }
1104: }
1106: jac->suboptionsset = PETSC_TRUE;
1107: PetscFunctionReturn(PETSC_SUCCESS);
1108: }
1110: #define FieldSplitSplitSolveAdd(ilink, xx, yy) \
1111: ((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) || \
1112: 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) || \
1113: VecScatterEnd(ilink->sctx, ilink->y, yy, ADD_VALUES, SCATTER_REVERSE)))
1115: static PetscErrorCode PCApply_FieldSplit_Schur(PC pc, Vec x, Vec y)
1116: {
1117: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1118: PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
1119: KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;
1121: PetscFunctionBegin;
1122: switch (jac->schurfactorization) {
1123: case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
1124: /* [A00 0; 0 -S], positive definite, suitable for MINRES */
1125: PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1126: PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
1127: PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1128: PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1129: PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y));
1130: PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
1131: PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1132: PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1133: PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
1134: PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1135: PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y));
1136: PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
1137: PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1138: PetscCall(VecScale(ilinkD->y, jac->schurscale));
1139: PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1140: PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1141: PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1142: break;
1143: case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
1144: /* [A00 0; A10 S], suitable for left preconditioning */
1145: PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1146: PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1147: PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1148: PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y));
1149: PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
1150: PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1151: PetscCall(MatMult(jac->C, ilinkA->y, ilinkD->x));
1152: PetscCall(VecScale(ilinkD->x, -1.));
1153: PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
1154: PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1155: PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
1156: PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1157: PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y));
1158: PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
1159: PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1160: PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1161: PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1162: PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1163: break;
1164: case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
1165: /* [A00 A01; 0 S], suitable for right preconditioning */
1166: PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
1167: PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
1168: PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1169: PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y));
1170: PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
1171: PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1172: PetscCall(MatMult(jac->B, ilinkD->y, ilinkA->x));
1173: PetscCall(VecScale(ilinkA->x, -1.));
1174: PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD));
1175: PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1176: PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD));
1177: PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1178: PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y));
1179: PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
1180: PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1181: PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1182: PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1183: PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1184: break;
1185: case PC_FIELDSPLIT_SCHUR_FACT_FULL:
1186: /* [1 0; A10 A00^{-1} 1] [A00 0; 0 S] [1 A00^{-1}A01; 0 1] */
1187: PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1188: PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1189: PetscCall(PetscLogEventBegin(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->y, NULL));
1190: PetscCall(KSPSolve(kspLower, ilinkA->x, ilinkA->y));
1191: PetscCall(KSPCheckSolve(kspLower, pc, ilinkA->y));
1192: PetscCall(PetscLogEventEnd(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->y, NULL));
1193: PetscCall(MatMult(jac->C, ilinkA->y, ilinkD->x));
1194: PetscCall(VecScale(ilinkD->x, -1.0));
1195: PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
1196: PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
1198: PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1199: PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y));
1200: PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
1201: PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1202: PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1204: if (kspUpper == kspA) {
1205: PetscCall(MatMult(jac->B, ilinkD->y, ilinkA->y));
1206: PetscCall(VecAXPY(ilinkA->x, -1.0, ilinkA->y));
1207: PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1208: PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y));
1209: PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
1210: PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1211: } else {
1212: PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1213: PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y));
1214: PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
1215: PetscCall(MatMult(jac->B, ilinkD->y, ilinkA->x));
1216: PetscCall(PetscLogEventBegin(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->z, NULL));
1217: PetscCall(KSPSolve(kspUpper, ilinkA->x, ilinkA->z));
1218: PetscCall(KSPCheckSolve(kspUpper, pc, ilinkA->z));
1219: PetscCall(PetscLogEventEnd(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->z, NULL));
1220: PetscCall(VecAXPY(ilinkA->y, -1.0, ilinkA->z));
1221: }
1222: PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1223: PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1224: PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1225: }
1226: PetscFunctionReturn(PETSC_SUCCESS);
1227: }
1229: static PetscErrorCode PCApplyTranspose_FieldSplit_Schur(PC pc, Vec x, Vec y)
1230: {
1231: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1232: PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
1233: KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;
1235: PetscFunctionBegin;
1236: switch (jac->schurfactorization) {
1237: case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
1238: /* [A00 0; 0 -S], positive definite, suitable for MINRES */
1239: PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1240: PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
1241: PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1242: PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1243: PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
1244: PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
1245: PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1246: PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1247: PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
1248: PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1249: PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y));
1250: PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
1251: PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1252: PetscCall(VecScale(ilinkD->y, jac->schurscale));
1253: PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1254: PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1255: PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1256: break;
1257: case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
1258: PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1259: PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1260: PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1261: PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
1262: PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
1263: PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1264: PetscCall(MatMultTranspose(jac->B, ilinkA->y, ilinkD->x));
1265: PetscCall(VecScale(ilinkD->x, -1.));
1266: PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
1267: PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1268: PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
1269: PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1270: PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y));
1271: PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
1272: PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1273: PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1274: PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1275: PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1276: break;
1277: case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
1278: PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
1279: PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
1280: PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1281: PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y));
1282: PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
1283: PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1284: PetscCall(MatMultTranspose(jac->C, ilinkD->y, ilinkA->x));
1285: PetscCall(VecScale(ilinkA->x, -1.));
1286: PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD));
1287: PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1288: PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD));
1289: PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1290: PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
1291: PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
1292: PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1293: PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1294: PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1295: PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1296: break;
1297: case PC_FIELDSPLIT_SCHUR_FACT_FULL:
1298: PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1299: PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1300: PetscCall(PetscLogEventBegin(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->y, NULL));
1301: PetscCall(KSPSolveTranspose(kspUpper, ilinkA->x, ilinkA->y));
1302: PetscCall(KSPCheckSolve(kspUpper, pc, ilinkA->y));
1303: PetscCall(PetscLogEventEnd(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->y, NULL));
1304: PetscCall(MatMultTranspose(jac->B, ilinkA->y, ilinkD->x));
1305: PetscCall(VecScale(ilinkD->x, -1.0));
1306: PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
1307: PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
1309: PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1310: PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y));
1311: PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
1312: PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
1313: PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1315: if (kspLower == kspA) {
1316: PetscCall(MatMultTranspose(jac->C, ilinkD->y, ilinkA->y));
1317: PetscCall(VecAXPY(ilinkA->x, -1.0, ilinkA->y));
1318: PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1319: PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
1320: PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
1321: PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1322: } else {
1323: PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
1324: PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
1325: PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
1326: PetscCall(MatMultTranspose(jac->C, ilinkD->y, ilinkA->x));
1327: PetscCall(PetscLogEventBegin(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->z, NULL));
1328: PetscCall(KSPSolveTranspose(kspLower, ilinkA->x, ilinkA->z));
1329: PetscCall(KSPCheckSolve(kspLower, pc, ilinkA->z));
1330: PetscCall(PetscLogEventEnd(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->z, NULL));
1331: PetscCall(VecAXPY(ilinkA->y, -1.0, ilinkA->z));
1332: }
1333: PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1334: PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1335: PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1336: }
1337: PetscFunctionReturn(PETSC_SUCCESS);
1338: }
1340: static PetscErrorCode PCApply_FieldSplit(PC pc, Vec x, Vec y)
1341: {
1342: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1343: PC_FieldSplitLink ilink = jac->head;
1344: PetscInt cnt, bs;
1346: PetscFunctionBegin;
1347: if (jac->type == PC_COMPOSITE_ADDITIVE) {
1348: if (jac->defaultsplit) {
1349: PetscCall(VecGetBlockSize(x, &bs));
1350: 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);
1351: PetscCall(VecGetBlockSize(y, &bs));
1352: 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);
1353: PetscCall(VecStrideGatherAll(x, jac->x, INSERT_VALUES));
1354: while (ilink) {
1355: PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
1356: PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y));
1357: PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
1358: PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
1359: ilink = ilink->next;
1360: }
1361: PetscCall(VecStrideScatterAll(jac->y, y, INSERT_VALUES));
1362: } else {
1363: PetscCall(VecSet(y, 0.0));
1364: while (ilink) {
1365: PetscCall(FieldSplitSplitSolveAdd(ilink, x, y));
1366: ilink = ilink->next;
1367: }
1368: }
1369: } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) {
1370: PetscCall(VecSet(y, 0.0));
1371: /* solve on first block for first block variables */
1372: PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, INSERT_VALUES, SCATTER_FORWARD));
1373: PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, INSERT_VALUES, SCATTER_FORWARD));
1374: PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
1375: PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y));
1376: PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
1377: PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
1378: PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
1379: PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
1381: /* compute the residual only onto second block variables using first block variables */
1382: PetscCall(MatMult(jac->Afield[1], ilink->y, ilink->next->x));
1383: ilink = ilink->next;
1384: PetscCall(VecScale(ilink->x, -1.0));
1385: PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));
1386: PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));
1388: /* solve on second block variables */
1389: PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
1390: PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y));
1391: PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
1392: PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
1393: PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
1394: PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
1395: } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1396: if (!jac->w1) {
1397: PetscCall(VecDuplicate(x, &jac->w1));
1398: PetscCall(VecDuplicate(x, &jac->w2));
1399: }
1400: PetscCall(VecSet(y, 0.0));
1401: PetscCall(FieldSplitSplitSolveAdd(ilink, x, y));
1402: cnt = 1;
1403: while (ilink->next) {
1404: ilink = ilink->next;
1405: /* compute the residual only over the part of the vector needed */
1406: PetscCall(MatMult(jac->Afield[cnt++], y, ilink->x));
1407: PetscCall(VecScale(ilink->x, -1.0));
1408: PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));
1409: PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));
1410: PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
1411: PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y));
1412: PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
1413: PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
1414: PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
1415: PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
1416: }
1417: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1418: cnt -= 2;
1419: while (ilink->previous) {
1420: ilink = ilink->previous;
1421: /* compute the residual only over the part of the vector needed */
1422: PetscCall(MatMult(jac->Afield[cnt--], y, ilink->x));
1423: PetscCall(VecScale(ilink->x, -1.0));
1424: PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));
1425: PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));
1426: PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
1427: PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y));
1428: PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
1429: PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
1430: PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
1431: PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
1432: }
1433: }
1434: } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Unsupported or unknown composition %d", (int)jac->type);
1435: PetscFunctionReturn(PETSC_SUCCESS);
1436: }
1438: static PetscErrorCode PCApply_FieldSplit_GKB(PC pc, Vec x, Vec y)
1439: {
1440: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1441: PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
1442: KSP ksp = ilinkA->ksp;
1443: Vec u, v, Hu, d, work1, work2;
1444: PetscScalar alpha, z, nrmz2, *vecz;
1445: PetscReal lowbnd, nu, beta;
1446: PetscInt j, iterGKB;
1448: PetscFunctionBegin;
1449: PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1450: PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
1451: PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
1452: PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
1454: u = jac->u;
1455: v = jac->v;
1456: Hu = jac->Hu;
1457: d = jac->d;
1458: work1 = jac->w1;
1459: work2 = jac->w2;
1460: vecz = jac->vecz;
1462: /* Change RHS to comply with matrix regularization H = A + nu*B*B' */
1463: /* Add q = q + nu*B*b */
1464: if (jac->gkbnu) {
1465: nu = jac->gkbnu;
1466: PetscCall(VecScale(ilinkD->x, jac->gkbnu));
1467: PetscCall(MatMultAdd(jac->B, ilinkD->x, ilinkA->x, ilinkA->x)); /* q = q + nu*B*b */
1468: } else {
1469: /* Situation when no augmented Lagrangian is used. Then we set inner */
1470: /* matrix N = I in [Ar13], and thus nu = 1. */
1471: nu = 1;
1472: }
1474: /* Transform rhs from [q,tilde{b}] to [0,b] */
1475: PetscCall(PetscLogEventBegin(ilinkA->event, ksp, ilinkA->x, ilinkA->y, NULL));
1476: PetscCall(KSPSolve(ksp, ilinkA->x, ilinkA->y));
1477: PetscCall(KSPCheckSolve(ksp, pc, ilinkA->y));
1478: PetscCall(PetscLogEventEnd(ilinkA->event, ksp, ilinkA->x, ilinkA->y, NULL));
1479: PetscCall(MatMultHermitianTranspose(jac->B, ilinkA->y, work1));
1480: PetscCall(VecAXPBY(work1, 1.0 / nu, -1.0, ilinkD->x)); /* c = b - B'*x */
1482: /* First step of algorithm */
1483: PetscCall(VecNorm(work1, NORM_2, &beta)); /* beta = sqrt(nu*c'*c)*/
1484: KSPCheckDot(ksp, beta);
1485: beta = PetscSqrtReal(nu) * beta;
1486: PetscCall(VecAXPBY(v, nu / beta, 0.0, work1)); /* v = nu/beta *c */
1487: PetscCall(MatMult(jac->B, v, work2)); /* u = H^{-1}*B*v */
1488: PetscCall(PetscLogEventBegin(ilinkA->event, ksp, work2, u, NULL));
1489: PetscCall(KSPSolve(ksp, work2, u));
1490: PetscCall(KSPCheckSolve(ksp, pc, u));
1491: PetscCall(PetscLogEventEnd(ilinkA->event, ksp, work2, u, NULL));
1492: PetscCall(MatMult(jac->H, u, Hu)); /* alpha = u'*H*u */
1493: PetscCall(VecDot(Hu, u, &alpha));
1494: KSPCheckDot(ksp, alpha);
1495: PetscCheck(PetscRealPart(alpha) > 0.0, PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "GKB preconditioner diverged, H is not positive definite");
1496: alpha = PetscSqrtReal(PetscAbsScalar(alpha));
1497: PetscCall(VecScale(u, 1.0 / alpha));
1498: PetscCall(VecAXPBY(d, 1.0 / alpha, 0.0, v)); /* v = nu/beta *c */
1500: z = beta / alpha;
1501: vecz[1] = z;
1503: /* Computation of first iterate x(1) and p(1) */
1504: PetscCall(VecAXPY(ilinkA->y, z, u));
1505: PetscCall(VecCopy(d, ilinkD->y));
1506: PetscCall(VecScale(ilinkD->y, -z));
1508: iterGKB = 1;
1509: lowbnd = 2 * jac->gkbtol;
1510: if (jac->gkbmonitor) PetscCall(PetscViewerASCIIPrintf(jac->gkbviewer, "%3" PetscInt_FMT " GKB Lower bound estimate %14.12e\n", iterGKB, (double)lowbnd));
1512: while (iterGKB < jac->gkbmaxit && lowbnd > jac->gkbtol) {
1513: iterGKB += 1;
1514: PetscCall(MatMultHermitianTranspose(jac->B, u, work1)); /* v <- nu*(B'*u-alpha/nu*v) */
1515: PetscCall(VecAXPBY(v, nu, -alpha, work1));
1516: PetscCall(VecNorm(v, NORM_2, &beta)); /* beta = sqrt(nu)*v'*v */
1517: beta = beta / PetscSqrtReal(nu);
1518: PetscCall(VecScale(v, 1.0 / beta));
1519: PetscCall(MatMult(jac->B, v, work2)); /* u <- H^{-1}*(B*v-beta*H*u) */
1520: PetscCall(MatMult(jac->H, u, Hu));
1521: PetscCall(VecAXPY(work2, -beta, Hu));
1522: PetscCall(PetscLogEventBegin(ilinkA->event, ksp, work2, u, NULL));
1523: PetscCall(KSPSolve(ksp, work2, u));
1524: PetscCall(KSPCheckSolve(ksp, pc, u));
1525: PetscCall(PetscLogEventEnd(ilinkA->event, ksp, work2, u, NULL));
1526: PetscCall(MatMult(jac->H, u, Hu)); /* alpha = u'*H*u */
1527: PetscCall(VecDot(Hu, u, &alpha));
1528: KSPCheckDot(ksp, alpha);
1529: PetscCheck(PetscRealPart(alpha) > 0.0, PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "GKB preconditioner diverged, H is not positive definite");
1530: alpha = PetscSqrtReal(PetscAbsScalar(alpha));
1531: PetscCall(VecScale(u, 1.0 / alpha));
1533: z = -beta / alpha * z; /* z <- beta/alpha*z */
1534: vecz[0] = z;
1536: /* Computation of new iterate x(i+1) and p(i+1) */
1537: PetscCall(VecAXPBY(d, 1.0 / alpha, -beta / alpha, v)); /* d = (v-beta*d)/alpha */
1538: PetscCall(VecAXPY(ilinkA->y, z, u)); /* r = r + z*u */
1539: PetscCall(VecAXPY(ilinkD->y, -z, d)); /* p = p - z*d */
1540: PetscCall(MatMult(jac->H, ilinkA->y, Hu)); /* ||u||_H = u'*H*u */
1541: PetscCall(VecDot(Hu, ilinkA->y, &nrmz2));
1543: /* Compute Lower Bound estimate */
1544: if (iterGKB > jac->gkbdelay) {
1545: lowbnd = 0.0;
1546: for (j = 0; j < jac->gkbdelay; j++) lowbnd += PetscAbsScalar(vecz[j] * vecz[j]);
1547: lowbnd = PetscSqrtReal(lowbnd / PetscAbsScalar(nrmz2));
1548: }
1550: for (j = 0; j < jac->gkbdelay - 1; j++) vecz[jac->gkbdelay - j - 1] = vecz[jac->gkbdelay - j - 2];
1551: if (jac->gkbmonitor) PetscCall(PetscViewerASCIIPrintf(jac->gkbviewer, "%3" PetscInt_FMT " GKB Lower bound estimate %14.12e\n", iterGKB, (double)lowbnd));
1552: }
1554: PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1555: PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
1556: PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1557: PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
1558: PetscFunctionReturn(PETSC_SUCCESS);
1559: }
1561: #define FieldSplitSplitSolveAddTranspose(ilink, xx, yy) \
1562: ((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) || \
1563: 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) || \
1564: VecScatterEnd(ilink->sctx, ilink->x, yy, ADD_VALUES, SCATTER_REVERSE)))
1566: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc, Vec x, Vec y)
1567: {
1568: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1569: PC_FieldSplitLink ilink = jac->head;
1570: PetscInt bs;
1572: PetscFunctionBegin;
1573: if (jac->type == PC_COMPOSITE_ADDITIVE) {
1574: if (jac->defaultsplit) {
1575: PetscCall(VecGetBlockSize(x, &bs));
1576: 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);
1577: PetscCall(VecGetBlockSize(y, &bs));
1578: 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);
1579: PetscCall(VecStrideGatherAll(x, jac->x, INSERT_VALUES));
1580: while (ilink) {
1581: PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
1582: PetscCall(KSPSolveTranspose(ilink->ksp, ilink->x, ilink->y));
1583: PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
1584: PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
1585: ilink = ilink->next;
1586: }
1587: PetscCall(VecStrideScatterAll(jac->y, y, INSERT_VALUES));
1588: } else {
1589: PetscCall(VecSet(y, 0.0));
1590: while (ilink) {
1591: PetscCall(FieldSplitSplitSolveAddTranspose(ilink, x, y));
1592: ilink = ilink->next;
1593: }
1594: }
1595: } else {
1596: if (!jac->w1) {
1597: PetscCall(VecDuplicate(x, &jac->w1));
1598: PetscCall(VecDuplicate(x, &jac->w2));
1599: }
1600: PetscCall(VecSet(y, 0.0));
1601: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1602: PetscCall(FieldSplitSplitSolveAddTranspose(ilink, x, y));
1603: while (ilink->next) {
1604: ilink = ilink->next;
1605: PetscCall(MatMultTranspose(pc->mat, y, jac->w1));
1606: PetscCall(VecWAXPY(jac->w2, -1.0, jac->w1, x));
1607: PetscCall(FieldSplitSplitSolveAddTranspose(ilink, jac->w2, y));
1608: }
1609: while (ilink->previous) {
1610: ilink = ilink->previous;
1611: PetscCall(MatMultTranspose(pc->mat, y, jac->w1));
1612: PetscCall(VecWAXPY(jac->w2, -1.0, jac->w1, x));
1613: PetscCall(FieldSplitSplitSolveAddTranspose(ilink, jac->w2, y));
1614: }
1615: } else {
1616: while (ilink->next) { /* get to last entry in linked list */
1617: ilink = ilink->next;
1618: }
1619: PetscCall(FieldSplitSplitSolveAddTranspose(ilink, x, y));
1620: while (ilink->previous) {
1621: ilink = ilink->previous;
1622: PetscCall(MatMultTranspose(pc->mat, y, jac->w1));
1623: PetscCall(VecWAXPY(jac->w2, -1.0, jac->w1, x));
1624: PetscCall(FieldSplitSplitSolveAddTranspose(ilink, jac->w2, y));
1625: }
1626: }
1627: }
1628: PetscFunctionReturn(PETSC_SUCCESS);
1629: }
1631: static PetscErrorCode PCReset_FieldSplit(PC pc)
1632: {
1633: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1634: PC_FieldSplitLink ilink = jac->head, next;
1636: PetscFunctionBegin;
1637: while (ilink) {
1638: PetscCall(KSPDestroy(&ilink->ksp));
1639: PetscCall(VecDestroy(&ilink->x));
1640: PetscCall(VecDestroy(&ilink->y));
1641: PetscCall(VecDestroy(&ilink->z));
1642: PetscCall(VecScatterDestroy(&ilink->sctx));
1643: PetscCall(ISDestroy(&ilink->is));
1644: PetscCall(ISDestroy(&ilink->is_col));
1645: PetscCall(PetscFree(ilink->splitname));
1646: PetscCall(PetscFree(ilink->fields));
1647: PetscCall(PetscFree(ilink->fields_col));
1648: next = ilink->next;
1649: PetscCall(PetscFree(ilink));
1650: ilink = next;
1651: }
1652: jac->head = NULL;
1653: PetscCall(PetscFree2(jac->x, jac->y));
1654: if (jac->mat && jac->mat != jac->pmat) {
1655: PetscCall(MatDestroyMatrices(jac->nsplits, &jac->mat));
1656: } else if (jac->mat) {
1657: jac->mat = NULL;
1658: }
1659: if (jac->pmat) PetscCall(MatDestroyMatrices(jac->nsplits, &jac->pmat));
1660: if (jac->Afield) PetscCall(MatDestroyMatrices(jac->nsplits, &jac->Afield));
1661: jac->nsplits = 0;
1662: PetscCall(VecDestroy(&jac->w1));
1663: PetscCall(VecDestroy(&jac->w2));
1664: PetscCall(MatDestroy(&jac->schur));
1665: PetscCall(MatDestroy(&jac->schurp));
1666: PetscCall(MatDestroy(&jac->schur_user));
1667: PetscCall(KSPDestroy(&jac->kspschur));
1668: PetscCall(KSPDestroy(&jac->kspupper));
1669: PetscCall(MatDestroy(&jac->B));
1670: PetscCall(MatDestroy(&jac->C));
1671: PetscCall(MatDestroy(&jac->H));
1672: PetscCall(VecDestroy(&jac->u));
1673: PetscCall(VecDestroy(&jac->v));
1674: PetscCall(VecDestroy(&jac->Hu));
1675: PetscCall(VecDestroy(&jac->d));
1676: PetscCall(PetscFree(jac->vecz));
1677: PetscCall(PetscViewerDestroy(&jac->gkbviewer));
1678: jac->isrestrict = PETSC_FALSE;
1679: PetscFunctionReturn(PETSC_SUCCESS);
1680: }
1682: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
1683: {
1684: PetscFunctionBegin;
1685: PetscCall(PCReset_FieldSplit(pc));
1686: PetscCall(PetscFree(pc->data));
1687: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
1688: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetFields_C", NULL));
1689: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetIS_C", NULL));
1690: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetType_C", NULL));
1691: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetBlockSize_C", NULL));
1692: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitRestrictIS_C", NULL));
1693: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSchurGetSubKSP_C", NULL));
1694: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", NULL));
1696: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", NULL));
1697: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", NULL));
1698: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", NULL));
1699: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", NULL));
1700: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", NULL));
1701: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", NULL));
1702: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", NULL));
1703: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", NULL));
1704: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", NULL));
1705: PetscFunctionReturn(PETSC_SUCCESS);
1706: }
1708: static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc, PetscOptionItems *PetscOptionsObject)
1709: {
1710: PetscInt bs;
1711: PetscBool flg;
1712: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1713: PCCompositeType ctype;
1715: PetscFunctionBegin;
1716: PetscOptionsHeadBegin(PetscOptionsObject, "FieldSplit options");
1717: PetscCall(PetscOptionsBool("-pc_fieldsplit_dm_splits", "Whether to use DMCreateFieldDecomposition() for splits", "PCFieldSplitSetDMSplits", jac->dm_splits, &jac->dm_splits, NULL));
1718: PetscCall(PetscOptionsInt("-pc_fieldsplit_block_size", "Blocksize that defines number of fields", "PCFieldSplitSetBlockSize", jac->bs, &bs, &flg));
1719: if (flg) PetscCall(PCFieldSplitSetBlockSize(pc, bs));
1720: jac->diag_use_amat = pc->useAmat;
1721: 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));
1722: jac->offdiag_use_amat = pc->useAmat;
1723: 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));
1724: PetscCall(PetscOptionsBool("-pc_fieldsplit_detect_saddle_point", "Form 2-way split by detecting zero diagonal entries", "PCFieldSplitSetDetectSaddlePoint", jac->detect, &jac->detect, NULL));
1725: PetscCall(PCFieldSplitSetDetectSaddlePoint(pc, jac->detect)); /* Sets split type and Schur PC type */
1726: PetscCall(PetscOptionsEnum("-pc_fieldsplit_type", "Type of composition", "PCFieldSplitSetType", PCCompositeTypes, (PetscEnum)jac->type, (PetscEnum *)&ctype, &flg));
1727: if (flg) PetscCall(PCFieldSplitSetType(pc, ctype));
1728: /* Only setup fields once */
1729: if ((jac->bs > 0) && (jac->nsplits == 0)) {
1730: /* only allow user to set fields from command line if bs is already known.
1731: otherwise user can set them in PCFieldSplitSetDefaults() */
1732: PetscCall(PCFieldSplitSetRuntimeSplits_Private(pc));
1733: if (jac->splitdefined) PetscCall(PetscInfo(pc, "Splits defined using the options database\n"));
1734: }
1735: if (jac->type == PC_COMPOSITE_SCHUR) {
1736: PetscCall(PetscOptionsGetEnum(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_fieldsplit_schur_factorization_type", PCFieldSplitSchurFactTypes, (PetscEnum *)&jac->schurfactorization, &flg));
1737: if (flg) PetscCall(PetscInfo(pc, "Deprecated use of -pc_fieldsplit_schur_factorization_type\n"));
1738: 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));
1739: PetscCall(PetscOptionsEnum("-pc_fieldsplit_schur_precondition", "How to build preconditioner for Schur complement", "PCFieldSplitSetSchurPre", PCFieldSplitSchurPreTypes, (PetscEnum)jac->schurpre, (PetscEnum *)&jac->schurpre, NULL));
1740: PetscCall(PetscOptionsScalar("-pc_fieldsplit_schur_scale", "Scale Schur complement", "PCFieldSplitSetSchurScale", jac->schurscale, &jac->schurscale, NULL));
1741: } else if (jac->type == PC_COMPOSITE_GKB) {
1742: PetscCall(PetscOptionsReal("-pc_fieldsplit_gkb_tol", "The tolerance for the lower bound stopping criterion", "PCFieldSplitSetGKBTol", jac->gkbtol, &jac->gkbtol, NULL));
1743: PetscCall(PetscOptionsInt("-pc_fieldsplit_gkb_delay", "The delay value for lower bound criterion", "PCFieldSplitSetGKBDelay", jac->gkbdelay, &jac->gkbdelay, NULL));
1744: PetscCall(PetscOptionsBoundedReal("-pc_fieldsplit_gkb_nu", "Parameter in augmented Lagrangian approach", "PCFieldSplitSetGKBNu", jac->gkbnu, &jac->gkbnu, NULL, 0.0));
1745: PetscCall(PetscOptionsInt("-pc_fieldsplit_gkb_maxit", "Maximum allowed number of iterations", "PCFieldSplitSetGKBMaxit", jac->gkbmaxit, &jac->gkbmaxit, NULL));
1746: PetscCall(PetscOptionsBool("-pc_fieldsplit_gkb_monitor", "Prints number of GKB iterations and error", "PCFieldSplitGKB", jac->gkbmonitor, &jac->gkbmonitor, NULL));
1747: }
1748: /*
1749: In the initial call to this routine the sub-solver data structures do not exist so we cannot call KSPSetFromOptions() on them yet.
1750: 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
1751: is called on the outer solver in case changes were made in the options database
1753: 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()
1754: 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.
1755: Without this extra check test p2p1fetidp_olof_full and others fail with incorrect matrix types.
1757: There could be a negative side effect of calling the KSPSetFromOptions() below.
1759: If one captured the PetscObjectState of the options database one could skip these calls if the database has not changed from the previous call
1760: */
1761: if (jac->issetup) {
1762: PC_FieldSplitLink ilink = jac->head;
1763: if (jac->type == PC_COMPOSITE_SCHUR) {
1764: if (jac->kspupper && jac->kspupper->totalits > 0) PetscCall(KSPSetFromOptions(jac->kspupper));
1765: if (jac->kspschur && jac->kspschur->totalits > 0) PetscCall(KSPSetFromOptions(jac->kspschur));
1766: }
1767: while (ilink) {
1768: if (ilink->ksp->totalits > 0) PetscCall(KSPSetFromOptions(ilink->ksp));
1769: ilink = ilink->next;
1770: }
1771: }
1772: PetscOptionsHeadEnd();
1773: PetscFunctionReturn(PETSC_SUCCESS);
1774: }
1776: static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc, const char splitname[], PetscInt n, const PetscInt *fields, const PetscInt *fields_col)
1777: {
1778: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1779: PC_FieldSplitLink ilink, next = jac->head;
1780: char prefix[128];
1781: PetscInt i;
1783: PetscFunctionBegin;
1784: if (jac->splitdefined) {
1785: PetscCall(PetscInfo(pc, "Ignoring new split \"%s\" because the splits have already been defined\n", splitname));
1786: PetscFunctionReturn(PETSC_SUCCESS);
1787: }
1788: for (i = 0; i < n; i++) {
1789: PetscCheck(fields[i] < jac->bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field %" PetscInt_FMT " requested but only %" PetscInt_FMT " exist", fields[i], jac->bs);
1790: PetscCheck(fields[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative field %" PetscInt_FMT " requested", fields[i]);
1791: }
1792: PetscCall(PetscNew(&ilink));
1793: if (splitname) {
1794: PetscCall(PetscStrallocpy(splitname, &ilink->splitname));
1795: } else {
1796: PetscCall(PetscMalloc1(3, &ilink->splitname));
1797: PetscCall(PetscSNPrintf(ilink->splitname, 2, "%" PetscInt_FMT, jac->nsplits));
1798: }
1799: ilink->event = jac->nsplits < 5 ? KSP_Solve_FS_0 + jac->nsplits : KSP_Solve_FS_0 + 4; /* Any split great than 4 gets logged in the 4th split */
1800: PetscCall(PetscMalloc1(n, &ilink->fields));
1801: PetscCall(PetscArraycpy(ilink->fields, fields, n));
1802: PetscCall(PetscMalloc1(n, &ilink->fields_col));
1803: PetscCall(PetscArraycpy(ilink->fields_col, fields_col, n));
1805: ilink->nfields = n;
1806: ilink->next = NULL;
1807: PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &ilink->ksp));
1808: PetscCall(KSPSetNestLevel(ilink->ksp, pc->kspnestlevel));
1809: PetscCall(KSPSetErrorIfNotConverged(ilink->ksp, pc->erroriffailure));
1810: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)pc, 1));
1811: PetscCall(KSPSetType(ilink->ksp, KSPPREONLY));
1813: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname));
1814: PetscCall(KSPSetOptionsPrefix(ilink->ksp, prefix));
1816: if (!next) {
1817: jac->head = ilink;
1818: ilink->previous = NULL;
1819: } else {
1820: while (next->next) next = next->next;
1821: next->next = ilink;
1822: ilink->previous = next;
1823: }
1824: jac->nsplits++;
1825: PetscFunctionReturn(PETSC_SUCCESS);
1826: }
1828: static PetscErrorCode PCFieldSplitSchurGetSubKSP_FieldSplit(PC pc, PetscInt *n, KSP **subksp)
1829: {
1830: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1832: PetscFunctionBegin;
1833: *subksp = NULL;
1834: if (n) *n = 0;
1835: if (jac->type == PC_COMPOSITE_SCHUR) {
1836: PetscInt nn;
1838: PetscCheck(jac->schur, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitSchurGetSubKSP()");
1839: PetscCheck(jac->nsplits == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Unexpected number of splits %" PetscInt_FMT " != 2", jac->nsplits);
1840: nn = jac->nsplits + (jac->kspupper != jac->head->ksp ? 1 : 0);
1841: PetscCall(PetscMalloc1(nn, subksp));
1842: (*subksp)[0] = jac->head->ksp;
1843: (*subksp)[1] = jac->kspschur;
1844: if (jac->kspupper != jac->head->ksp) (*subksp)[2] = jac->kspupper;
1845: if (n) *n = nn;
1846: }
1847: PetscFunctionReturn(PETSC_SUCCESS);
1848: }
1850: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc, PetscInt *n, KSP **subksp)
1851: {
1852: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1854: PetscFunctionBegin;
1855: PetscCheck(jac->schur, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitGetSubKSP()");
1856: PetscCall(PetscMalloc1(jac->nsplits, subksp));
1857: PetscCall(MatSchurComplementGetKSP(jac->schur, *subksp));
1859: (*subksp)[1] = jac->kspschur;
1860: if (n) *n = jac->nsplits;
1861: PetscFunctionReturn(PETSC_SUCCESS);
1862: }
1864: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc, PetscInt *n, KSP **subksp)
1865: {
1866: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1867: PetscInt cnt = 0;
1868: PC_FieldSplitLink ilink = jac->head;
1870: PetscFunctionBegin;
1871: PetscCall(PetscMalloc1(jac->nsplits, subksp));
1872: while (ilink) {
1873: (*subksp)[cnt++] = ilink->ksp;
1874: ilink = ilink->next;
1875: }
1876: 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);
1877: if (n) *n = jac->nsplits;
1878: PetscFunctionReturn(PETSC_SUCCESS);
1879: }
1881: /*@C
1882: PCFieldSplitRestrictIS - Restricts the fieldsplit `IS`s to be within a given `IS`.
1884: Input Parameters:
1885: + pc - the preconditioner context
1886: - isy - the index set that defines the indices to which the fieldsplit is to be restricted
1888: Level: advanced
1890: Developer Notes:
1891: It seems the resulting `IS`s will not cover the entire space, so
1892: how can they define a convergent preconditioner? Needs explaining.
1894: .seealso: [](sec_block_matrices), `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`
1895: @*/
1896: PetscErrorCode PCFieldSplitRestrictIS(PC pc, IS isy)
1897: {
1898: PetscFunctionBegin;
1901: PetscTryMethod(pc, "PCFieldSplitRestrictIS_C", (PC, IS), (pc, isy));
1902: PetscFunctionReturn(PETSC_SUCCESS);
1903: }
1905: static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy)
1906: {
1907: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1908: PC_FieldSplitLink ilink = jac->head, next;
1909: PetscInt localsize, size, sizez, i;
1910: const PetscInt *ind, *indz;
1911: PetscInt *indc, *indcz;
1912: PetscBool flg;
1914: PetscFunctionBegin;
1915: PetscCall(ISGetLocalSize(isy, &localsize));
1916: PetscCallMPI(MPI_Scan(&localsize, &size, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)isy)));
1917: size -= localsize;
1918: while (ilink) {
1919: IS isrl, isr;
1920: PC subpc;
1921: PetscCall(ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl));
1922: PetscCall(ISGetLocalSize(isrl, &localsize));
1923: PetscCall(PetscMalloc1(localsize, &indc));
1924: PetscCall(ISGetIndices(isrl, &ind));
1925: PetscCall(PetscArraycpy(indc, ind, localsize));
1926: PetscCall(ISRestoreIndices(isrl, &ind));
1927: PetscCall(ISDestroy(&isrl));
1928: for (i = 0; i < localsize; i++) *(indc + i) += size;
1929: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)isy), localsize, indc, PETSC_OWN_POINTER, &isr));
1930: PetscCall(PetscObjectReference((PetscObject)isr));
1931: PetscCall(ISDestroy(&ilink->is));
1932: ilink->is = isr;
1933: PetscCall(PetscObjectReference((PetscObject)isr));
1934: PetscCall(ISDestroy(&ilink->is_col));
1935: ilink->is_col = isr;
1936: PetscCall(ISDestroy(&isr));
1937: PetscCall(KSPGetPC(ilink->ksp, &subpc));
1938: PetscCall(PetscObjectTypeCompare((PetscObject)subpc, PCFIELDSPLIT, &flg));
1939: if (flg) {
1940: IS iszl, isz;
1941: MPI_Comm comm;
1942: PetscCall(ISGetLocalSize(ilink->is, &localsize));
1943: comm = PetscObjectComm((PetscObject)ilink->is);
1944: PetscCall(ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl));
1945: PetscCallMPI(MPI_Scan(&localsize, &sizez, 1, MPIU_INT, MPI_SUM, comm));
1946: sizez -= localsize;
1947: PetscCall(ISGetLocalSize(iszl, &localsize));
1948: PetscCall(PetscMalloc1(localsize, &indcz));
1949: PetscCall(ISGetIndices(iszl, &indz));
1950: PetscCall(PetscArraycpy(indcz, indz, localsize));
1951: PetscCall(ISRestoreIndices(iszl, &indz));
1952: PetscCall(ISDestroy(&iszl));
1953: for (i = 0; i < localsize; i++) *(indcz + i) += sizez;
1954: PetscCall(ISCreateGeneral(comm, localsize, indcz, PETSC_OWN_POINTER, &isz));
1955: PetscCall(PCFieldSplitRestrictIS(subpc, isz));
1956: PetscCall(ISDestroy(&isz));
1957: }
1958: next = ilink->next;
1959: ilink = next;
1960: }
1961: jac->isrestrict = PETSC_TRUE;
1962: PetscFunctionReturn(PETSC_SUCCESS);
1963: }
1965: static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc, const char splitname[], IS is)
1966: {
1967: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1968: PC_FieldSplitLink ilink, next = jac->head;
1969: char prefix[128];
1971: PetscFunctionBegin;
1972: if (jac->splitdefined) {
1973: PetscCall(PetscInfo(pc, "Ignoring new split \"%s\" because the splits have already been defined\n", splitname));
1974: PetscFunctionReturn(PETSC_SUCCESS);
1975: }
1976: PetscCall(PetscNew(&ilink));
1977: if (splitname) {
1978: PetscCall(PetscStrallocpy(splitname, &ilink->splitname));
1979: } else {
1980: PetscCall(PetscMalloc1(8, &ilink->splitname));
1981: PetscCall(PetscSNPrintf(ilink->splitname, 7, "%" PetscInt_FMT, jac->nsplits));
1982: }
1983: ilink->event = jac->nsplits < 5 ? KSP_Solve_FS_0 + jac->nsplits : KSP_Solve_FS_0 + 4; /* Any split great than 4 gets logged in the 4th split */
1984: PetscCall(PetscObjectReference((PetscObject)is));
1985: PetscCall(ISDestroy(&ilink->is));
1986: ilink->is = is;
1987: PetscCall(PetscObjectReference((PetscObject)is));
1988: PetscCall(ISDestroy(&ilink->is_col));
1989: ilink->is_col = is;
1990: ilink->next = NULL;
1991: PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &ilink->ksp));
1992: PetscCall(KSPSetNestLevel(ilink->ksp, pc->kspnestlevel));
1993: PetscCall(KSPSetErrorIfNotConverged(ilink->ksp, pc->erroriffailure));
1994: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)pc, 1));
1995: PetscCall(KSPSetType(ilink->ksp, KSPPREONLY));
1997: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname));
1998: PetscCall(KSPSetOptionsPrefix(ilink->ksp, prefix));
2000: if (!next) {
2001: jac->head = ilink;
2002: ilink->previous = NULL;
2003: } else {
2004: while (next->next) next = next->next;
2005: next->next = ilink;
2006: ilink->previous = next;
2007: }
2008: jac->nsplits++;
2009: PetscFunctionReturn(PETSC_SUCCESS);
2010: }
2012: /*@C
2013: PCFieldSplitSetFields - Sets the fields that define one particular split in `PCFIELDSPLIT`
2015: Logically Collective
2017: Input Parameters:
2018: + pc - the preconditioner context
2019: . splitname - name of this split, if `NULL` the number of the split is used
2020: . n - the number of fields in this split
2021: . fields - the fields in this split
2022: - fields_col - generally the same as fields, if it does not match fields then the matrix block that is solved for this set of fields comes from an off-diagonal block
2023: of the matrix and fields_col provides the column indices for that block
2025: Level: intermediate
2027: Notes:
2028: Use `PCFieldSplitSetIS()` to set a general set of indices as a split.
2030: `PCFieldSplitSetFields()` is for defining fields as strided blocks. For example, if the block
2031: 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
2032: 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x....
2033: where the numbered entries indicate what is in the split.
2035: This function is called once per split (it creates a new split each time). Solve options
2036: for this split will be available under the prefix `-fieldsplit_SPLITNAME_`.
2038: `PCFieldSplitSetIS()` does not support having a fields_col different from fields
2040: Developer Notes:
2041: This routine does not actually create the `IS` representing the split, that is delayed
2042: until `PCSetUp_FieldSplit()`, because information about the vector/matrix layouts may not be
2043: available when this routine is called.
2045: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetBlockSize()`, `PCFieldSplitSetIS()`, `PCFieldSplitRestrictIS()`
2046: @*/
2047: PetscErrorCode PCFieldSplitSetFields(PC pc, const char splitname[], PetscInt n, const PetscInt *fields, const PetscInt *fields_col)
2048: {
2049: PetscFunctionBegin;
2051: PetscAssertPointer(splitname, 2);
2052: PetscCheck(n >= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Provided number of fields %" PetscInt_FMT " in split \"%s\" not positive", n, splitname);
2053: PetscAssertPointer(fields, 4);
2054: PetscTryMethod(pc, "PCFieldSplitSetFields_C", (PC, const char[], PetscInt, const PetscInt *, const PetscInt *), (pc, splitname, n, fields, fields_col));
2055: PetscFunctionReturn(PETSC_SUCCESS);
2056: }
2058: /*@
2059: PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) to build
2060: the sub-matrices associated with each split. Where `KSPSetOperators`(ksp,Amat,Pmat) was used to supply the operators.
2062: Logically Collective
2064: Input Parameters:
2065: + pc - the preconditioner object
2066: - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from
2068: Options Database Key:
2069: . -pc_fieldsplit_diag_use_amat - use the Amat to provide the diagonal blocks
2071: Level: intermediate
2073: .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitGetDiagUseAmat()`, `PCFieldSplitSetOffDiagUseAmat()`, `PCFIELDSPLIT`
2074: @*/
2075: PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc, PetscBool flg)
2076: {
2077: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2078: PetscBool isfs;
2080: PetscFunctionBegin;
2082: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
2083: PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT);
2084: jac->diag_use_amat = flg;
2085: PetscFunctionReturn(PETSC_SUCCESS);
2086: }
2088: /*@
2089: PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) to build
2090: the sub-matrices associated with each split. Where `KSPSetOperators`(ksp,Amat,Pmat) was used to supply the operators.
2092: Logically Collective
2094: Input Parameter:
2095: . pc - the preconditioner object
2097: Output Parameter:
2098: . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from
2100: Level: intermediate
2102: .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitSetDiagUseAmat()`, `PCFieldSplitGetOffDiagUseAmat()`, `PCFIELDSPLIT`
2103: @*/
2104: PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc, PetscBool *flg)
2105: {
2106: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2107: PetscBool isfs;
2109: PetscFunctionBegin;
2111: PetscAssertPointer(flg, 2);
2112: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
2113: PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT);
2114: *flg = jac->diag_use_amat;
2115: PetscFunctionReturn(PETSC_SUCCESS);
2116: }
2118: /*@
2119: PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) to build
2120: the sub-matrices associated with each split. Where `KSPSetOperators`(ksp,Amat,Pmat) was used to supply the operators.
2122: Logically Collective
2124: Input Parameters:
2125: + pc - the preconditioner object
2126: - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from
2128: Options Database Key:
2129: . -pc_fieldsplit_off_diag_use_amat <bool> - use the Amat to extract the off-diagonal blocks
2131: Level: intermediate
2133: .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitGetOffDiagUseAmat()`, `PCFieldSplitSetDiagUseAmat()`, `PCFIELDSPLIT`
2134: @*/
2135: PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc, PetscBool flg)
2136: {
2137: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2138: PetscBool isfs;
2140: PetscFunctionBegin;
2142: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
2143: PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT);
2144: jac->offdiag_use_amat = flg;
2145: PetscFunctionReturn(PETSC_SUCCESS);
2146: }
2148: /*@
2149: PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) to build
2150: the sub-matrices associated with each split. Where `KSPSetOperators`(ksp,Amat,Pmat) was used to supply the operators.
2152: Logically Collective
2154: Input Parameter:
2155: . pc - the preconditioner object
2157: Output Parameter:
2158: . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from
2160: Level: intermediate
2162: .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitSetOffDiagUseAmat()`, `PCFieldSplitGetDiagUseAmat()`, `PCFIELDSPLIT`
2163: @*/
2164: PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc, PetscBool *flg)
2165: {
2166: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2167: PetscBool isfs;
2169: PetscFunctionBegin;
2171: PetscAssertPointer(flg, 2);
2172: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
2173: PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT);
2174: *flg = jac->offdiag_use_amat;
2175: PetscFunctionReturn(PETSC_SUCCESS);
2176: }
2178: /*@C
2179: PCFieldSplitSetIS - Sets the exact elements for a split in a `PCFIELDSPLIT`
2181: Logically Collective
2183: Input Parameters:
2184: + pc - the preconditioner context
2185: . splitname - name of this split, if `NULL` the number of the split is used
2186: - is - the index set that defines the elements in this split
2188: Level: intermediate
2190: Notes:
2191: Use `PCFieldSplitSetFields()`, for splits defined by strided types.
2193: This function is called once per split (it creates a new split each time). Solve options
2194: for this split will be available under the prefix -fieldsplit_SPLITNAME_.
2196: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetBlockSize()`
2197: @*/
2198: PetscErrorCode PCFieldSplitSetIS(PC pc, const char splitname[], IS is)
2199: {
2200: PetscFunctionBegin;
2202: if (splitname) PetscAssertPointer(splitname, 2);
2204: PetscTryMethod(pc, "PCFieldSplitSetIS_C", (PC, const char[], IS), (pc, splitname, is));
2205: PetscFunctionReturn(PETSC_SUCCESS);
2206: }
2208: /*@C
2209: PCFieldSplitGetIS - Retrieves the elements for a split as an `IS`
2211: Logically Collective
2213: Input Parameters:
2214: + pc - the preconditioner context
2215: - splitname - name of this split
2217: Output Parameter:
2218: . is - the index set that defines the elements in this split, or `NULL` if the split is not found
2220: Level: intermediate
2222: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetIS()`
2223: @*/
2224: PetscErrorCode PCFieldSplitGetIS(PC pc, const char splitname[], IS *is)
2225: {
2226: PetscFunctionBegin;
2228: PetscAssertPointer(splitname, 2);
2229: PetscAssertPointer(is, 3);
2230: {
2231: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2232: PC_FieldSplitLink ilink = jac->head;
2233: PetscBool found;
2235: *is = NULL;
2236: while (ilink) {
2237: PetscCall(PetscStrcmp(ilink->splitname, splitname, &found));
2238: if (found) {
2239: *is = ilink->is;
2240: break;
2241: }
2242: ilink = ilink->next;
2243: }
2244: }
2245: PetscFunctionReturn(PETSC_SUCCESS);
2246: }
2248: /*@C
2249: PCFieldSplitGetISByIndex - Retrieves the elements for a given split as an `IS`
2251: Logically Collective
2253: Input Parameters:
2254: + pc - the preconditioner context
2255: - index - index of this split
2257: Output Parameter:
2258: . is - the index set that defines the elements in this split
2260: Level: intermediate
2262: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitGetIS()`, `PCFieldSplitSetIS()`
2263: @*/
2264: PetscErrorCode PCFieldSplitGetISByIndex(PC pc, PetscInt index, IS *is)
2265: {
2266: PetscFunctionBegin;
2267: PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative field %" PetscInt_FMT " requested", index);
2269: PetscAssertPointer(is, 3);
2270: {
2271: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2272: PC_FieldSplitLink ilink = jac->head;
2273: PetscInt i = 0;
2274: PetscCheck(index < jac->nsplits, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field %" PetscInt_FMT " requested but only %" PetscInt_FMT " exist", index, jac->nsplits);
2276: while (i < index) {
2277: ilink = ilink->next;
2278: ++i;
2279: }
2280: PetscCall(PCFieldSplitGetIS(pc, ilink->splitname, is));
2281: }
2282: PetscFunctionReturn(PETSC_SUCCESS);
2283: }
2285: /*@
2286: PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the
2287: fieldsplit preconditioner when calling `PCFieldSplitSetIS()`. If not set the matrix block size is used.
2289: Logically Collective
2291: Input Parameters:
2292: + pc - the preconditioner context
2293: - bs - the block size
2295: Level: intermediate
2297: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`
2298: @*/
2299: PetscErrorCode PCFieldSplitSetBlockSize(PC pc, PetscInt bs)
2300: {
2301: PetscFunctionBegin;
2304: PetscTryMethod(pc, "PCFieldSplitSetBlockSize_C", (PC, PetscInt), (pc, bs));
2305: PetscFunctionReturn(PETSC_SUCCESS);
2306: }
2308: /*@C
2309: PCFieldSplitGetSubKSP - Gets the `KSP` contexts for all splits
2311: Collective
2313: Input Parameter:
2314: . pc - the preconditioner context
2316: Output Parameters:
2317: + n - the number of splits
2318: - subksp - the array of `KSP` contexts
2320: Level: advanced
2322: Notes:
2323: After `PCFieldSplitGetSubKSP()` the array of `KSP`s is to be freed by the user with `PetscFree()`
2324: (not the `KSP`, just the array that contains them).
2326: You must call `PCSetUp()` before calling `PCFieldSplitGetSubKSP()`.
2328: If the fieldsplit is of type `PC_COMPOSITE_SCHUR`, it returns the `KSP` object used inside the
2329: Schur complement and the `KSP` object used to iterate over the Schur complement.
2330: To access all the `KSP` objects used in `PC_COMPOSITE_SCHUR`, use `PCFieldSplitSchurGetSubKSP()`.
2332: If the fieldsplit is of type `PC_COMPOSITE_GKB`, it returns the `KSP` object used to solve the
2333: inner linear system defined by the matrix H in each loop.
2335: Fortran Notes:
2336: You must pass in a `KSP` array that is large enough to contain all the `KSP`s.
2337: You can call `PCFieldSplitGetSubKSP`(pc,n,`PETSC_NULL_KSP`,ierr) to determine how large the
2338: `KSP` array must be.
2340: Developer Notes:
2341: There should be a `PCFieldSplitRestoreSubKSP()` instead of requiring the user to call `PetscFree()`
2343: The Fortran interface should be modernized to return directly the array of values.
2345: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`, `PCFieldSplitSchurGetSubKSP()`
2346: @*/
2347: PetscErrorCode PCFieldSplitGetSubKSP(PC pc, PetscInt *n, KSP *subksp[])
2348: {
2349: PetscFunctionBegin;
2351: if (n) PetscAssertPointer(n, 2);
2352: PetscUseMethod(pc, "PCFieldSplitGetSubKSP_C", (PC, PetscInt *, KSP **), (pc, n, subksp));
2353: PetscFunctionReturn(PETSC_SUCCESS);
2354: }
2356: /*@C
2357: PCFieldSplitSchurGetSubKSP - Gets the `KSP` contexts used inside the Schur complement based `PCFIELDSPLIT`
2359: Collective
2361: Input Parameter:
2362: . pc - the preconditioner context
2364: Output Parameters:
2365: + n - the number of splits
2366: - subksp - the array of `KSP` contexts
2368: Level: advanced
2370: Notes:
2371: After `PCFieldSplitSchurGetSubKSP()` the array of `KSP`s is to be freed by the user with `PetscFree()`
2372: (not the `KSP` just the array that contains them).
2374: You must call `PCSetUp()` before calling `PCFieldSplitSchurGetSubKSP()`.
2376: If the fieldsplit type is of type `PC_COMPOSITE_SCHUR`, it returns (in order)
2377: + 1 - the `KSP` used for the (1,1) block
2378: . 2 - the `KSP` used for the Schur complement (not the one used for the interior Schur solver)
2379: - 3 - the `KSP` used for the (1,1) block in the upper triangular factor (if different from that of the (1,1) block).
2381: It returns a null array if the fieldsplit is not of type `PC_COMPOSITE_SCHUR`; in this case, you should use `PCFieldSplitGetSubKSP()`.
2383: Fortran Notes:
2384: You must pass in a `KSP` array that is large enough to contain all the local `KSP`s.
2385: You can call `PCFieldSplitSchurGetSubKSP`(pc,n,`PETSC_NULL_KSP`,ierr) to determine how large the
2386: `KSP` array must be.
2388: Developer Notes:
2389: There should be a `PCFieldSplitRestoreSubKSP()` instead of requiring the user to call `PetscFree()`
2391: Should the functionality of `PCFieldSplitSchurGetSubKSP()` and `PCFieldSplitGetSubKSP()` be merged?
2393: The Fortran interface should be modernized to return directly the array of values.
2395: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`, `PCFieldSplitGetSubKSP()`
2396: @*/
2397: PetscErrorCode PCFieldSplitSchurGetSubKSP(PC pc, PetscInt *n, KSP *subksp[])
2398: {
2399: PetscFunctionBegin;
2401: if (n) PetscAssertPointer(n, 2);
2402: PetscUseMethod(pc, "PCFieldSplitSchurGetSubKSP_C", (PC, PetscInt *, KSP **), (pc, n, subksp));
2403: PetscFunctionReturn(PETSC_SUCCESS);
2404: }
2406: /*@
2407: PCFieldSplitSetSchurPre - Indicates from what operator the preconditioner is constructed for the Schur complement.
2408: The default is the A11 matrix.
2410: Collective
2412: Input Parameters:
2413: + pc - the preconditioner context
2414: . ptype - which matrix to use for preconditioning the Schur complement: `PC_FIELDSPLIT_SCHUR_PRE_A11` (default),
2415: `PC_FIELDSPLIT_SCHUR_PRE_SELF`, `PC_FIELDSPLIT_SCHUR_PRE_USER`,
2416: `PC_FIELDSPLIT_SCHUR_PRE_SELFP`, and `PC_FIELDSPLIT_SCHUR_PRE_FULL`
2417: - pre - matrix to use for preconditioning, or `NULL`
2419: Options Database Keys:
2420: + -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is `a11`. See notes for meaning of various arguments
2421: - -fieldsplit_1_pc_type <pctype> - the preconditioner algorithm that is used to construct the preconditioner from the operator
2423: Level: intermediate
2425: Notes:
2426: If ptype is
2427: + a11 - the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner
2428: matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix
2429: . self - the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix:
2430: The only preconditioners that currently work with this symbolic representation matrix object are `PCLSC` and `PCHPDDM`
2431: . user - the preconditioner for the Schur complement is generated from the user provided matrix (pre argument
2432: to this function).
2433: . selfp - the preconditioning for the Schur complement is generated from an explicitly-assembled approximation $ Sp = A11 - A10 inv(diag(A00)) A01 $
2434: This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be
2435: lumped before extracting the diagonal using the additional option `-fieldsplit_1_mat_schur_complement_ainv_type lump`
2436: - full - the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation
2437: computed internally by `PCFIELDSPLIT` (this is expensive)
2438: useful mostly as a test that the Schur complement approach can work for your problem
2440: When solving a saddle point problem, where the A11 block is identically zero, using `a11` as the ptype only makes sense
2441: with the additional option `-fieldsplit_1_pc_type none`. Usually for saddle point problems one would use a `ptype` of `self` and
2442: `-fieldsplit_1_pc_type lsc` which uses the least squares commutator to compute a preconditioner for the Schur complement.
2444: Developer Note:
2445: The name of this function and the option `-pc_fieldsplit_schur_precondition` are inconsistent; precondition should be used everywhere.
2447: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSchurPre()`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`,
2448: `MatSchurComplementSetAinvType()`, `PCLSC`, `PCFieldSplitSetSchurFactType()`
2449: @*/
2450: PetscErrorCode PCFieldSplitSetSchurPre(PC pc, PCFieldSplitSchurPreType ptype, Mat pre)
2451: {
2452: PetscFunctionBegin;
2454: PetscTryMethod(pc, "PCFieldSplitSetSchurPre_C", (PC, PCFieldSplitSchurPreType, Mat), (pc, ptype, pre));
2455: PetscFunctionReturn(PETSC_SUCCESS);
2456: }
2458: PetscErrorCode PCFieldSplitSchurPrecondition(PC pc, PCFieldSplitSchurPreType ptype, Mat pre)
2459: {
2460: return PCFieldSplitSetSchurPre(pc, ptype, pre);
2461: } /* Deprecated name */
2463: /*@
2464: PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be
2465: preconditioned. See `PCFieldSplitSetSchurPre()` for details.
2467: Logically Collective
2469: Input Parameter:
2470: . pc - the preconditioner context
2472: Output Parameters:
2473: + 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`
2474: - pre - matrix to use for preconditioning (with `PC_FIELDSPLIT_SCHUR_PRE_USER`), or `NULL`
2476: Level: intermediate
2478: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitSetSchurPre()`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, `PCLSC`
2479: @*/
2480: PetscErrorCode PCFieldSplitGetSchurPre(PC pc, PCFieldSplitSchurPreType *ptype, Mat *pre)
2481: {
2482: PetscFunctionBegin;
2484: PetscUseMethod(pc, "PCFieldSplitGetSchurPre_C", (PC, PCFieldSplitSchurPreType *, Mat *), (pc, ptype, pre));
2485: PetscFunctionReturn(PETSC_SUCCESS);
2486: }
2488: /*@
2489: PCFieldSplitSchurGetS - extract the `MATSCHURCOMPLEMENT` object used by this `PCFIELDSPLIT` in case it needs to be configured separately
2491: Not Collective
2493: Input Parameter:
2494: . pc - the preconditioner context
2496: Output Parameter:
2497: . S - the Schur complement matrix
2499: Level: advanced
2501: Note:
2502: This matrix should not be destroyed using `MatDestroy()`; rather, use `PCFieldSplitSchurRestoreS()`.
2504: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurPre()`, `MATSCHURCOMPLEMENT`, `PCFieldSplitSchurRestoreS()`,
2505: `MatCreateSchurComplement()`, `MatSchurComplementGetKSP()`, `MatSchurComplementComputeExplicitOperator()`, `MatGetSchurComplement()`
2506: @*/
2507: PetscErrorCode PCFieldSplitSchurGetS(PC pc, Mat *S)
2508: {
2509: const char *t;
2510: PetscBool isfs;
2511: PC_FieldSplit *jac;
2513: PetscFunctionBegin;
2515: PetscCall(PetscObjectGetType((PetscObject)pc, &t));
2516: PetscCall(PetscStrcmp(t, PCFIELDSPLIT, &isfs));
2517: PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PC of type PCFIELDSPLIT, got %s instead", t);
2518: jac = (PC_FieldSplit *)pc->data;
2519: PetscCheck(jac->type == PC_COMPOSITE_SCHUR, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PCFIELDSPLIT of type SCHUR, got %d instead", jac->type);
2520: if (S) *S = jac->schur;
2521: PetscFunctionReturn(PETSC_SUCCESS);
2522: }
2524: /*@
2525: PCFieldSplitSchurRestoreS - returns the `MATSCHURCOMPLEMENT` matrix used by this `PC`
2527: Not Collective
2529: Input Parameters:
2530: + pc - the preconditioner context
2531: - S - the Schur complement matrix
2533: Level: advanced
2535: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurPre()`, `MatSchurComplement`, `PCFieldSplitSchurGetS()`
2536: @*/
2537: PetscErrorCode PCFieldSplitSchurRestoreS(PC pc, Mat *S)
2538: {
2539: const char *t;
2540: PetscBool isfs;
2541: PC_FieldSplit *jac;
2543: PetscFunctionBegin;
2545: PetscCall(PetscObjectGetType((PetscObject)pc, &t));
2546: PetscCall(PetscStrcmp(t, PCFIELDSPLIT, &isfs));
2547: PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PC of type PCFIELDSPLIT, got %s instead", t);
2548: jac = (PC_FieldSplit *)pc->data;
2549: PetscCheck(jac->type == PC_COMPOSITE_SCHUR, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PCFIELDSPLIT of type SCHUR, got %d instead", jac->type);
2550: PetscCheck(S && (*S == jac->schur), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "MatSchurComplement restored is not the same as gotten");
2551: PetscFunctionReturn(PETSC_SUCCESS);
2552: }
2554: static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc, PCFieldSplitSchurPreType ptype, Mat pre)
2555: {
2556: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2558: PetscFunctionBegin;
2559: jac->schurpre = ptype;
2560: if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
2561: PetscCall(MatDestroy(&jac->schur_user));
2562: jac->schur_user = pre;
2563: PetscCall(PetscObjectReference((PetscObject)jac->schur_user));
2564: }
2565: PetscFunctionReturn(PETSC_SUCCESS);
2566: }
2568: static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc, PCFieldSplitSchurPreType *ptype, Mat *pre)
2569: {
2570: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2572: PetscFunctionBegin;
2573: if (ptype) *ptype = jac->schurpre;
2574: if (pre) *pre = jac->schur_user;
2575: PetscFunctionReturn(PETSC_SUCCESS);
2576: }
2578: /*@
2579: PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain in the preconditioner {cite}`murphy2000note` and {cite}`ipsen2001note`
2581: Collective
2583: Input Parameters:
2584: + pc - the preconditioner context
2585: - ftype - which blocks of factorization to retain, `PC_FIELDSPLIT_SCHUR_FACT_FULL` is default
2587: Options Database Key:
2588: . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> - default is `full`
2590: Level: intermediate
2592: Notes:
2593: The FULL factorization is
2595: ```{math}
2596: \left(\begin{array}{cc} A & B \\
2597: C & E \\
2598: \end{array}\right) =
2599: \left(\begin{array}{cc} 1 & 0 \\
2600: C*A^{-1} & I \\
2601: \end{array}\right)
2602: \left(\begin{array}{cc} A & 0 \\
2603: 0 & S \\
2604: \end{array}\right)
2605: \left(\begin{array}{cc} I & A^{-1}B \\
2606: 0 & I \\
2607: \end{array}\right) = L D U.
2608: ```
2610: where $ S = E - C*A^{-1}*B $. In practice, the full factorization is applied via block triangular solves with the grouping $L*(D*U)$. UPPER uses $D*U$, LOWER uses $L*D$,
2611: and DIAG is the diagonal part with the sign of $ S $ flipped (because this makes the preconditioner positive definite for many formulations,
2612: thus allowing the use of `KSPMINRES)`. Sign flipping of $ S $ can be turned off with `PCFieldSplitSetSchurScale()`.
2614: If $A$ and $S$ are solved exactly
2615: + 1 - FULL factorization is a direct solver.
2616: . 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.
2617: - 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.
2619: If the iteration count is very low, consider using `KSPFGMRES` or `KSPGCR` which can use one less preconditioner
2620: application in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice.
2622: For symmetric problems in which $A$ is positive definite and $S$ is negative definite, DIAG can be used with `KSPMINRES`.
2624: 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).
2626: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurScale()`,
2627: [](sec_flexibleksp), `PCFieldSplitSetSchurPre()`
2628: @*/
2629: PetscErrorCode PCFieldSplitSetSchurFactType(PC pc, PCFieldSplitSchurFactType ftype)
2630: {
2631: PetscFunctionBegin;
2633: PetscTryMethod(pc, "PCFieldSplitSetSchurFactType_C", (PC, PCFieldSplitSchurFactType), (pc, ftype));
2634: PetscFunctionReturn(PETSC_SUCCESS);
2635: }
2637: static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc, PCFieldSplitSchurFactType ftype)
2638: {
2639: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2641: PetscFunctionBegin;
2642: jac->schurfactorization = ftype;
2643: PetscFunctionReturn(PETSC_SUCCESS);
2644: }
2646: /*@
2647: PCFieldSplitSetSchurScale - Controls the sign flip of S for `PC_FIELDSPLIT_SCHUR_FACT_DIAG`.
2649: Collective
2651: Input Parameters:
2652: + pc - the preconditioner context
2653: - scale - scaling factor for the Schur complement
2655: Options Database Key:
2656: . -pc_fieldsplit_schur_scale <scale> - default is -1.0
2658: Level: intermediate
2660: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurFactType`, `PCFieldSplitSetSchurFactType()`
2661: @*/
2662: PetscErrorCode PCFieldSplitSetSchurScale(PC pc, PetscScalar scale)
2663: {
2664: PetscFunctionBegin;
2667: PetscTryMethod(pc, "PCFieldSplitSetSchurScale_C", (PC, PetscScalar), (pc, scale));
2668: PetscFunctionReturn(PETSC_SUCCESS);
2669: }
2671: static PetscErrorCode PCFieldSplitSetSchurScale_FieldSplit(PC pc, PetscScalar scale)
2672: {
2673: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2675: PetscFunctionBegin;
2676: jac->schurscale = scale;
2677: PetscFunctionReturn(PETSC_SUCCESS);
2678: }
2680: /*@C
2681: PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement
2683: Collective
2685: Input Parameter:
2686: . pc - the preconditioner context
2688: Output Parameters:
2689: + A00 - the (0,0) block
2690: . A01 - the (0,1) block
2691: . A10 - the (1,0) block
2692: - A11 - the (1,1) block
2694: Level: advanced
2696: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `MatSchurComplementGetSubMatrices()`, `MatSchurComplementSetSubMatrices()`
2697: @*/
2698: PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc, Mat *A00, Mat *A01, Mat *A10, Mat *A11)
2699: {
2700: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2702: PetscFunctionBegin;
2704: PetscCheck(jac->type == PC_COMPOSITE_SCHUR, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach.");
2705: if (A00) *A00 = jac->pmat[0];
2706: if (A01) *A01 = jac->B;
2707: if (A10) *A10 = jac->C;
2708: if (A11) *A11 = jac->pmat[1];
2709: PetscFunctionReturn(PETSC_SUCCESS);
2710: }
2712: /*@
2713: PCFieldSplitSetGKBTol - Sets the solver tolerance for the generalized Golub-Kahan bidiagonalization preconditioner {cite}`arioli2013` in `PCFIELDSPLIT`
2715: Collective
2717: Input Parameters:
2718: + pc - the preconditioner context
2719: - tolerance - the solver tolerance
2721: Options Database Key:
2722: . -pc_fieldsplit_gkb_tol <tolerance> - default is 1e-5
2724: Level: intermediate
2726: Note:
2727: The generalized GKB algorithm {cite}`arioli2013` uses a lower bound estimate of the error in energy norm as stopping criterion.
2728: It stops once the lower bound estimate undershoots the required solver tolerance. Although the actual error might be bigger than
2729: this estimate, the stopping criterion is satisfactory in practical cases.
2731: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBNu()`, `PCFieldSplitSetGKBMaxit()`
2732: @*/
2733: PetscErrorCode PCFieldSplitSetGKBTol(PC pc, PetscReal tolerance)
2734: {
2735: PetscFunctionBegin;
2738: PetscTryMethod(pc, "PCFieldSplitSetGKBTol_C", (PC, PetscReal), (pc, tolerance));
2739: PetscFunctionReturn(PETSC_SUCCESS);
2740: }
2742: static PetscErrorCode PCFieldSplitSetGKBTol_FieldSplit(PC pc, PetscReal tolerance)
2743: {
2744: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2746: PetscFunctionBegin;
2747: jac->gkbtol = tolerance;
2748: PetscFunctionReturn(PETSC_SUCCESS);
2749: }
2751: /*@
2752: PCFieldSplitSetGKBMaxit - Sets the maximum number of iterations for the generalized Golub-Kahan bidiagonalization preconditioner in `PCFIELDSPLIT`
2754: Collective
2756: Input Parameters:
2757: + pc - the preconditioner context
2758: - maxit - the maximum number of iterations
2760: Options Database Key:
2761: . -pc_fieldsplit_gkb_maxit <maxit> - default is 100
2763: Level: intermediate
2765: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBNu()`
2766: @*/
2767: PetscErrorCode PCFieldSplitSetGKBMaxit(PC pc, PetscInt maxit)
2768: {
2769: PetscFunctionBegin;
2772: PetscTryMethod(pc, "PCFieldSplitSetGKBMaxit_C", (PC, PetscInt), (pc, maxit));
2773: PetscFunctionReturn(PETSC_SUCCESS);
2774: }
2776: static PetscErrorCode PCFieldSplitSetGKBMaxit_FieldSplit(PC pc, PetscInt maxit)
2777: {
2778: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2780: PetscFunctionBegin;
2781: jac->gkbmaxit = maxit;
2782: PetscFunctionReturn(PETSC_SUCCESS);
2783: }
2785: /*@
2786: PCFieldSplitSetGKBDelay - Sets the delay in the lower bound error estimate in the generalized Golub-Kahan bidiagonalization {cite}`arioli2013` in `PCFIELDSPLIT`
2787: preconditioner.
2789: Collective
2791: Input Parameters:
2792: + pc - the preconditioner context
2793: - delay - the delay window in the lower bound estimate
2795: Options Database Key:
2796: . -pc_fieldsplit_gkb_delay <delay> - default is 5
2798: Level: intermediate
2800: Notes:
2801: 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 $
2802: is expressed as a truncated sum. The error at iteration k can only be measured at iteration (k + `delay`), and thus the algorithm needs
2803: at least (`delay` + 1) iterations to stop.
2805: For more details on the generalized Golub-Kahan bidiagonalization method and its lower bound stopping criterion, please refer to {cite}`arioli2013`
2807: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBNu()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBMaxit()`
2808: @*/
2809: PetscErrorCode PCFieldSplitSetGKBDelay(PC pc, PetscInt delay)
2810: {
2811: PetscFunctionBegin;
2814: PetscTryMethod(pc, "PCFieldSplitSetGKBDelay_C", (PC, PetscInt), (pc, delay));
2815: PetscFunctionReturn(PETSC_SUCCESS);
2816: }
2818: static PetscErrorCode PCFieldSplitSetGKBDelay_FieldSplit(PC pc, PetscInt delay)
2819: {
2820: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2822: PetscFunctionBegin;
2823: jac->gkbdelay = delay;
2824: PetscFunctionReturn(PETSC_SUCCESS);
2825: }
2827: /*@
2828: PCFieldSplitSetGKBNu - Sets the scalar value nu >= 0 in the transformation H = A00 + nu*A01*A01' of the (1,1) block in the
2829: Golub-Kahan bidiagonalization preconditioner {cite}`arioli2013` in `PCFIELDSPLIT`
2831: Collective
2833: Input Parameters:
2834: + pc - the preconditioner context
2835: - nu - the shift parameter
2837: Options Database Key:
2838: . -pc_fieldsplit_gkb_nu <nu> - default is 1
2840: Level: intermediate
2842: Notes:
2843: 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,
2844: 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
2845: necessary to find a good balance in between the convergence of the inner and outer loop.
2847: 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.
2849: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBMaxit()`
2850: @*/
2851: PetscErrorCode PCFieldSplitSetGKBNu(PC pc, PetscReal nu)
2852: {
2853: PetscFunctionBegin;
2856: PetscTryMethod(pc, "PCFieldSplitSetGKBNu_C", (PC, PetscReal), (pc, nu));
2857: PetscFunctionReturn(PETSC_SUCCESS);
2858: }
2860: static PetscErrorCode PCFieldSplitSetGKBNu_FieldSplit(PC pc, PetscReal nu)
2861: {
2862: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2864: PetscFunctionBegin;
2865: jac->gkbnu = nu;
2866: PetscFunctionReturn(PETSC_SUCCESS);
2867: }
2869: static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc, PCCompositeType type)
2870: {
2871: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2873: PetscFunctionBegin;
2874: jac->type = type;
2875: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", NULL));
2876: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", NULL));
2877: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", NULL));
2878: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", NULL));
2879: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", NULL));
2880: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", NULL));
2881: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", NULL));
2882: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", NULL));
2883: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", NULL));
2885: if (type == PC_COMPOSITE_SCHUR) {
2886: pc->ops->apply = PCApply_FieldSplit_Schur;
2887: pc->ops->applytranspose = PCApplyTranspose_FieldSplit_Schur;
2888: pc->ops->view = PCView_FieldSplit_Schur;
2890: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit_Schur));
2891: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", PCFieldSplitSetSchurPre_FieldSplit));
2892: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", PCFieldSplitGetSchurPre_FieldSplit));
2893: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", PCFieldSplitSetSchurFactType_FieldSplit));
2894: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", PCFieldSplitSetSchurScale_FieldSplit));
2895: } else if (type == PC_COMPOSITE_GKB) {
2896: pc->ops->apply = PCApply_FieldSplit_GKB;
2897: pc->ops->view = PCView_FieldSplit_GKB;
2899: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit));
2900: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", PCFieldSplitSetGKBTol_FieldSplit));
2901: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", PCFieldSplitSetGKBMaxit_FieldSplit));
2902: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", PCFieldSplitSetGKBNu_FieldSplit));
2903: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", PCFieldSplitSetGKBDelay_FieldSplit));
2904: } else {
2905: pc->ops->apply = PCApply_FieldSplit;
2906: pc->ops->view = PCView_FieldSplit;
2908: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit));
2909: }
2910: PetscFunctionReturn(PETSC_SUCCESS);
2911: }
2913: static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc, PetscInt bs)
2914: {
2915: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2917: PetscFunctionBegin;
2918: PetscCheck(bs >= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Blocksize must be positive, you gave %" PetscInt_FMT, bs);
2919: 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);
2920: jac->bs = bs;
2921: PetscFunctionReturn(PETSC_SUCCESS);
2922: }
2924: static PetscErrorCode PCSetCoordinates_FieldSplit(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[])
2925: {
2926: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2927: PC_FieldSplitLink ilink_current = jac->head;
2928: IS is_owned;
2930: PetscFunctionBegin;
2931: jac->coordinates_set = PETSC_TRUE; // Internal flag
2932: PetscCall(MatGetOwnershipIS(pc->mat, &is_owned, NULL));
2934: while (ilink_current) {
2935: // For each IS, embed it to get local coords indces
2936: IS is_coords;
2937: PetscInt ndofs_block;
2938: const PetscInt *block_dofs_enumeration; // Numbering of the dofs relevant to the current block
2940: // Setting drop to true for safety. It should make no difference.
2941: PetscCall(ISEmbed(ilink_current->is, is_owned, PETSC_TRUE, &is_coords));
2942: PetscCall(ISGetLocalSize(is_coords, &ndofs_block));
2943: PetscCall(ISGetIndices(is_coords, &block_dofs_enumeration));
2945: // Allocate coordinates vector and set it directly
2946: PetscCall(PetscMalloc1(ndofs_block * dim, &ilink_current->coords));
2947: for (PetscInt dof = 0; dof < ndofs_block; ++dof) {
2948: for (PetscInt d = 0; d < dim; ++d) (ilink_current->coords)[dim * dof + d] = coords[dim * block_dofs_enumeration[dof] + d];
2949: }
2950: ilink_current->dim = dim;
2951: ilink_current->ndofs = ndofs_block;
2952: PetscCall(ISRestoreIndices(is_coords, &block_dofs_enumeration));
2953: PetscCall(ISDestroy(&is_coords));
2954: ilink_current = ilink_current->next;
2955: }
2956: PetscCall(ISDestroy(&is_owned));
2957: PetscFunctionReturn(PETSC_SUCCESS);
2958: }
2960: /*@
2961: PCFieldSplitSetType - Sets the type, `PCCompositeType`, of a `PCFIELDSPLIT`
2963: Collective
2965: Input Parameters:
2966: + pc - the preconditioner context
2967: - type - `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE` (default), `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`,
2968: `PC_COMPOSITE_GKB`
2970: Options Database Key:
2971: . -pc_fieldsplit_type <one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type
2973: Level: intermediate
2975: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCCompositeType`, `PCCompositeGetType()`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`,
2976: `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`, `PCFieldSplitSetSchurFactType()`
2977: @*/
2978: PetscErrorCode PCFieldSplitSetType(PC pc, PCCompositeType type)
2979: {
2980: PetscFunctionBegin;
2982: PetscTryMethod(pc, "PCFieldSplitSetType_C", (PC, PCCompositeType), (pc, type));
2983: PetscFunctionReturn(PETSC_SUCCESS);
2984: }
2986: /*@
2987: PCFieldSplitGetType - Gets the type, `PCCompositeType`, of a `PCFIELDSPLIT`
2989: Not collective
2991: Input Parameter:
2992: . pc - the preconditioner context
2994: Output Parameter:
2995: . type - `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE` (default), `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`
2997: Level: intermediate
2999: .seealso: [](sec_block_matrices), `PC`, `PCCompositeSetType()`, `PCFIELDSPLIT`, `PCCompositeType`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`,
3000: `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`
3001: @*/
3002: PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
3003: {
3004: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
3006: PetscFunctionBegin;
3008: PetscAssertPointer(type, 2);
3009: *type = jac->type;
3010: PetscFunctionReturn(PETSC_SUCCESS);
3011: }
3013: /*@
3014: PCFieldSplitSetDMSplits - Flags whether `DMCreateFieldDecomposition()` should be used to define the splits in a `PCFIELDSPLIT`, whenever possible.
3016: Logically Collective
3018: Input Parameters:
3019: + pc - the preconditioner context
3020: - flg - boolean indicating whether to use field splits defined by the `DM`
3022: Options Database Key:
3023: . -pc_fieldsplit_dm_splits <bool> - use the field splits defined by the `DM`
3025: Level: intermediate
3027: Developer Note:
3028: The name should be `PCFieldSplitSetUseDMSplits()`, similar change to options database
3030: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitGetDMSplits()`, `DMCreateFieldDecomposition()`, `PCFieldSplitSetFields()`, `PCFieldsplitSetIS()`
3031: @*/
3032: PetscErrorCode PCFieldSplitSetDMSplits(PC pc, PetscBool flg)
3033: {
3034: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
3035: PetscBool isfs;
3037: PetscFunctionBegin;
3040: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
3041: if (isfs) jac->dm_splits = flg;
3042: PetscFunctionReturn(PETSC_SUCCESS);
3043: }
3045: /*@
3046: PCFieldSplitGetDMSplits - Returns flag indicating whether `DMCreateFieldDecomposition()` should be used to define the splits in a `PCFIELDSPLIT`, whenever possible.
3048: Logically Collective
3050: Input Parameter:
3051: . pc - the preconditioner context
3053: Output Parameter:
3054: . flg - boolean indicating whether to use field splits defined by the `DM`
3056: Level: intermediate
3058: Developer Note:
3059: The name should be `PCFieldSplitGetUseDMSplits()`
3061: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetDMSplits()`, `DMCreateFieldDecomposition()`, `PCFieldSplitSetFields()`, `PCFieldsplitSetIS()`
3062: @*/
3063: PetscErrorCode PCFieldSplitGetDMSplits(PC pc, PetscBool *flg)
3064: {
3065: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
3066: PetscBool isfs;
3068: PetscFunctionBegin;
3070: PetscAssertPointer(flg, 2);
3071: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
3072: if (isfs) {
3073: if (flg) *flg = jac->dm_splits;
3074: }
3075: PetscFunctionReturn(PETSC_SUCCESS);
3076: }
3078: /*@
3079: PCFieldSplitGetDetectSaddlePoint - Returns flag indicating whether `PCFIELDSPLIT` will attempt to automatically determine fields based on zero diagonal entries.
3081: Logically Collective
3083: Input Parameter:
3084: . pc - the preconditioner context
3086: Output Parameter:
3087: . flg - boolean indicating whether to detect fields or not
3089: Level: intermediate
3091: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetDetectSaddlePoint()`
3092: @*/
3093: PetscErrorCode PCFieldSplitGetDetectSaddlePoint(PC pc, PetscBool *flg)
3094: {
3095: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
3097: PetscFunctionBegin;
3098: *flg = jac->detect;
3099: PetscFunctionReturn(PETSC_SUCCESS);
3100: }
3102: /*@
3103: PCFieldSplitSetDetectSaddlePoint - Sets flag indicating whether `PCFIELDSPLIT` will attempt to automatically determine fields based on zero diagonal entries.
3105: Logically Collective
3107: Input Parameter:
3108: . pc - the preconditioner context
3110: Output Parameter:
3111: . flg - boolean indicating whether to detect fields or not
3113: Options Database Key:
3114: . -pc_fieldsplit_detect_saddle_point <bool> - detect and use the saddle point
3116: Level: intermediate
3118: Note:
3119: Also sets the split type to `PC_COMPOSITE_SCHUR` (see `PCFieldSplitSetType()`) and the Schur preconditioner type to `PC_FIELDSPLIT_SCHUR_PRE_SELF` (see `PCFieldSplitSetSchurPre()`).
3121: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitGetDetectSaddlePoint()`, `PCFieldSplitSetType()`, `PCFieldSplitSetSchurPre()`, `PC_FIELDSPLIT_SCHUR_PRE_SELF`
3122: @*/
3123: PetscErrorCode PCFieldSplitSetDetectSaddlePoint(PC pc, PetscBool flg)
3124: {
3125: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
3127: PetscFunctionBegin;
3128: jac->detect = flg;
3129: if (jac->detect) {
3130: PetscCall(PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR));
3131: PetscCall(PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_SELF, NULL));
3132: }
3133: PetscFunctionReturn(PETSC_SUCCESS);
3134: }
3136: /*MC
3137: PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
3138: collections of variables (that may overlap) called splits. See [the users manual section on "Solving Block Matrices"](sec_block_matrices) for more details.
3140: Options Database Keys:
3141: + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the `%d`'th split
3142: . -pc_fieldsplit_default - automatically add any fields to additional splits that have not
3143: been supplied explicitly by `-pc_fieldsplit_%d_fields`
3144: . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields)
3145: . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur,gkb> - type of relaxation or factorization splitting
3146: . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is `a11`; see `PCFieldSplitSetSchurPre()`
3147: . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> - set factorization type when using `-pc_fieldsplit_type schur`;
3148: see `PCFieldSplitSetSchurFactType()`
3149: . -pc_fieldsplit_dm_splits <true,false> (default is true) - Whether to use `DMCreateFieldDecomposition()` for splits
3150: - -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero diagonal and uses Schur complement with no preconditioner as the solver
3152: Options prefixes for inner solvers when using the Schur complement preconditioner are `-fieldsplit_0_` and `-fieldsplit_1_` .
3153: The options prefix for the inner solver when using the Golub-Kahan biadiagonalization preconditioner is `-fieldsplit_0_`
3154: For all other solvers they are `-fieldsplit_%d_` for the `%d`'th field; use `-fieldsplit_` for all fields.
3156: To set options on the solvers for each block append `-fieldsplit_` to all the `PC`
3157: options database keys. For example, `-fieldsplit_pc_type ilu` `-fieldsplit_pc_factor_levels 1`
3159: To set the options on the solvers separate for each block call `PCFieldSplitGetSubKSP()`
3160: and set the options directly on the resulting `KSP` object
3162: Level: intermediate
3164: Notes:
3165: Use `PCFieldSplitSetFields()` to set splits defined by "strided" entries and `PCFieldSplitSetIS()`
3166: to define a split by an arbitrary collection of entries.
3168: If no splits are set, the default is used. If a `DM` is associated with the `PC` and it supports
3169: `DMCreateFieldDecomposition()`, then that is used for the default. Otherwise, the splits are defined by entries strided by bs,
3170: beginning at 0 then 1, etc to bs-1. The block size can be set with `PCFieldSplitSetBlockSize()`,
3171: if this is not called the block size defaults to the blocksize of the second matrix passed
3172: to `KSPSetOperators()`/`PCSetOperators()`.
3174: For the Schur complement preconditioner if
3176: ```{math}
3177: J = \left[\begin{array}{cc} A_{00} & A_{01} \\ A_{10} & A_{11} \end{array}\right]
3178: ```
3180: the preconditioner using `full` factorization is logically
3181: ```{math}
3182: \left[\begin{array}{cc} I & -\text{ksp}(A_{00}) A_{01} \\ 0 & I \end{array}\right] \left[\begin{array}{cc} \text{inv}(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]
3183: ```
3184: where the action of $\text{inv}(A_{00})$ is applied using the KSP solver with prefix `-fieldsplit_0_`. $S$ is the Schur complement
3185: ```{math}
3186: S = A_{11} - A_{10} \text{ksp}(A_{00}) A_{01}
3187: ```
3188: 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` was given
3189: in providing the SECOND split or 1 if not given). For `PCFieldSplitGetSubKSP()` when field number is 0,
3190: it returns the KSP associated with `-fieldsplit_0_` while field number 1 gives `-fieldsplit_1_` KSP. By default
3191: $A_{11}$ is used to construct a preconditioner for $S$, use `PCFieldSplitSetSchurPre()` for all the possible ways to construct the preconditioner for $S$.
3193: The factorization type is set using `-pc_fieldsplit_schur_fact_type <diag, lower, upper, full>`. `full` is shown above,
3194: `diag` gives
3195: ```{math}
3196: \left[\begin{array}{cc} \text{inv}(A_{00}) & 0 \\ 0 & -\text{ksp}(S) \end{array}\right]
3197: ```
3198: 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
3199: can be turned off with `PCFieldSplitSetSchurScale()` or by command line `-pc_fieldsplit_schur_scale 1.0`. The `lower` factorization is the inverse of
3200: ```{math}
3201: \left[\begin{array}{cc} A_{00} & 0 \\ A_{10} & S \end{array}\right]
3202: ```
3203: where the inverses of A_{00} and S are applied using KSPs. The upper factorization is the inverse of
3204: ```{math}
3205: \left[\begin{array}{cc} A_{00} & A_{01} \\ 0 & S \end{array}\right]
3206: ```
3207: where again the inverses of $A_{00}$ and $S$ are applied using `KSP`s.
3209: If only one set of indices (one `IS`) is provided with `PCFieldSplitSetIS()` then the complement of that `IS`
3210: is used automatically for a second block.
3212: The fieldsplit preconditioner cannot currently be used with the `MATBAIJ` or `MATSBAIJ` data formats if the blocksize is larger than 1.
3213: Generally it should be used with the `MATAIJ` format.
3215: The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see,
3216: for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling {cite}`wesseling2009`.
3217: One can also use `PCFIELDSPLIT`
3218: inside a smoother resulting in "Distributive Smoothers".
3220: See "A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations" {cite}`elman2008tcp`.
3222: The Constrained Pressure Preconditioner (CPR) can be implemented using `PCCOMPOSITE` with `PCGALERKIN`. CPR first solves an $R A P$ subsystem, updates the
3223: residual on all variables (`PCCompositeSetType(pc,PC_COMPOSITE_MULTIPLICATIVE)`), and then applies a simple ILU like preconditioner on all the variables.
3225: The generalized Golub-Kahan bidiagonalization preconditioner (GKB) can be applied to symmetric $2 \times 2$ block matrices of the shape
3226: ```{math}
3227: \left[\begin{array}{cc} A_{00} & A_{01} \\ A_{01}' & 0 \end{array}\right]
3228: ```
3229: 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}'$.
3230: 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_`.
3232: Developer Note:
3233: The Schur complement functionality of `PCFIELDSPLIT` should likely be factored into its own `PC` thus simplifying the implementation of the preconditioners and their
3234: user API.
3236: .seealso: [](sec_block_matrices), `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCLSC`,
3237: `PCFieldSplitGetSubKSP()`, `PCFieldSplitSchurGetSubKSP()`, `PCFieldSplitSetFields()`,
3238: `PCFieldSplitSetType()`, `PCFieldSplitSetIS()`, `PCFieldSplitSetSchurPre()`, `PCFieldSplitSetSchurFactType()`,
3239: `MatSchurComplementSetAinvType()`, `PCFieldSplitSetSchurScale()`, `PCFieldSplitSetDetectSaddlePoint()`
3240: M*/
3242: PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
3243: {
3244: PC_FieldSplit *jac;
3246: PetscFunctionBegin;
3247: PetscCall(PetscNew(&jac));
3249: jac->bs = -1;
3250: jac->nsplits = 0;
3251: jac->type = PC_COMPOSITE_MULTIPLICATIVE;
3252: jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
3253: jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
3254: jac->schurscale = -1.0;
3255: jac->dm_splits = PETSC_TRUE;
3256: jac->detect = PETSC_FALSE;
3257: jac->gkbtol = 1e-5;
3258: jac->gkbdelay = 5;
3259: jac->gkbnu = 1;
3260: jac->gkbmaxit = 100;
3261: jac->gkbmonitor = PETSC_FALSE;
3262: jac->coordinates_set = PETSC_FALSE;
3264: pc->data = (void *)jac;
3266: pc->ops->apply = PCApply_FieldSplit;
3267: pc->ops->applytranspose = PCApplyTranspose_FieldSplit;
3268: pc->ops->setup = PCSetUp_FieldSplit;
3269: pc->ops->reset = PCReset_FieldSplit;
3270: pc->ops->destroy = PCDestroy_FieldSplit;
3271: pc->ops->setfromoptions = PCSetFromOptions_FieldSplit;
3272: pc->ops->view = PCView_FieldSplit;
3273: pc->ops->applyrichardson = NULL;
3275: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSchurGetSubKSP_C", PCFieldSplitSchurGetSubKSP_FieldSplit));
3276: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit));
3277: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetFields_C", PCFieldSplitSetFields_FieldSplit));
3278: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetIS_C", PCFieldSplitSetIS_FieldSplit));
3279: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetType_C", PCFieldSplitSetType_FieldSplit));
3280: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetBlockSize_C", PCFieldSplitSetBlockSize_FieldSplit));
3281: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitRestrictIS_C", PCFieldSplitRestrictIS_FieldSplit));
3282: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_FieldSplit));
3283: PetscFunctionReturn(PETSC_SUCCESS);
3284: }