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: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
107: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
108: if (iascii) {
109: if (jac->bs > 0) {
110: PetscViewerASCIIPrintf(viewer, " FieldSplit with %s composition: total splits = %" PetscInt_FMT ", blocksize = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits, jac->bs);
111: } else {
112: PetscViewerASCIIPrintf(viewer, " FieldSplit with %s composition: total splits = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits);
113: }
114: if (pc->useAmat) PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for blocks\n");
115: if (jac->diag_use_amat) PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for diagonal blocks\n");
116: if (jac->offdiag_use_amat) PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for off-diagonal blocks\n");
117: PetscViewerASCIIPrintf(viewer, " Solver info for each split is in the following KSP objects:\n");
118: for (i = 0; i < jac->nsplits; i++) {
119: if (ilink->fields) {
120: PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Fields ", i);
121: PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
122: for (j = 0; j < ilink->nfields; j++) {
123: if (j > 0) PetscViewerASCIIPrintf(viewer, ",");
124: PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, ilink->fields[j]);
125: }
126: PetscViewerASCIIPrintf(viewer, "\n");
127: PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
128: } else {
129: PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Defined by IS\n", i);
130: }
131: KSPView(ilink->ksp, viewer);
132: ilink = ilink->next;
133: }
134: }
136: if (isdraw) {
137: PetscDraw draw;
138: PetscReal x, y, w, wd;
140: PetscViewerDrawGetDraw(viewer, 0, &draw);
141: PetscDrawGetCurrentPoint(draw, &x, &y);
142: w = 2 * PetscMin(1.0 - x, x);
143: wd = w / (jac->nsplits + 1);
144: x = x - wd * (jac->nsplits - 1) / 2.0;
145: for (i = 0; i < jac->nsplits; i++) {
146: PetscDrawPushCurrentPoint(draw, x, y);
147: KSPView(ilink->ksp, viewer);
148: PetscDrawPopCurrentPoint(draw);
149: x += wd;
150: ilink = ilink->next;
151: }
152: }
153: return 0;
154: }
156: static PetscErrorCode PCView_FieldSplit_Schur(PC pc, PetscViewer viewer)
157: {
158: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
159: PetscBool iascii, isdraw;
160: PetscInt i, j;
161: PC_FieldSplitLink ilink = jac->head;
162: MatSchurComplementAinvType atype;
164: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
165: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
166: if (iascii) {
167: if (jac->bs > 0) {
168: PetscViewerASCIIPrintf(viewer, " FieldSplit with Schur preconditioner, blocksize = %" PetscInt_FMT ", factorization %s\n", jac->bs, PCFieldSplitSchurFactTypes[jac->schurfactorization]);
169: } else {
170: PetscViewerASCIIPrintf(viewer, " FieldSplit with Schur preconditioner, factorization %s\n", PCFieldSplitSchurFactTypes[jac->schurfactorization]);
171: }
172: if (pc->useAmat) PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for blocks\n");
173: switch (jac->schurpre) {
174: case PC_FIELDSPLIT_SCHUR_PRE_SELF:
175: PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from S itself\n");
176: break;
177: case PC_FIELDSPLIT_SCHUR_PRE_SELFP:
178: MatSchurComplementGetAinvType(jac->schur, &atype);
179: 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 ")));
180: break;
181: case PC_FIELDSPLIT_SCHUR_PRE_A11:
182: PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from A11\n");
183: break;
184: case PC_FIELDSPLIT_SCHUR_PRE_FULL:
185: PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from the exact Schur complement\n");
186: break;
187: case PC_FIELDSPLIT_SCHUR_PRE_USER:
188: if (jac->schur_user) {
189: PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from user provided matrix\n");
190: } else {
191: PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from A11\n");
192: }
193: break;
194: default:
195: SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre);
196: }
197: PetscViewerASCIIPrintf(viewer, " Split info:\n");
198: PetscViewerASCIIPushTab(viewer);
199: for (i = 0; i < jac->nsplits; i++) {
200: if (ilink->fields) {
201: PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Fields ", i);
202: PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
203: for (j = 0; j < ilink->nfields; j++) {
204: if (j > 0) PetscViewerASCIIPrintf(viewer, ",");
205: PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, ilink->fields[j]);
206: }
207: PetscViewerASCIIPrintf(viewer, "\n");
208: PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
209: } else {
210: PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Defined by IS\n", i);
211: }
212: ilink = ilink->next;
213: }
214: PetscViewerASCIIPrintf(viewer, "KSP solver for A00 block\n");
215: PetscViewerASCIIPushTab(viewer);
216: if (jac->head) {
217: KSPView(jac->head->ksp, viewer);
218: } else PetscViewerASCIIPrintf(viewer, " not yet available\n");
219: PetscViewerASCIIPopTab(viewer);
220: if (jac->head && jac->kspupper != jac->head->ksp) {
221: PetscViewerASCIIPrintf(viewer, "KSP solver for upper A00 in upper triangular factor \n");
222: PetscViewerASCIIPushTab(viewer);
223: if (jac->kspupper) KSPView(jac->kspupper, viewer);
224: else PetscViewerASCIIPrintf(viewer, " not yet available\n");
225: PetscViewerASCIIPopTab(viewer);
226: }
227: PetscViewerASCIIPrintf(viewer, "KSP solver for S = A11 - A10 inv(A00) A01 \n");
228: PetscViewerASCIIPushTab(viewer);
229: if (jac->kspschur) {
230: KSPView(jac->kspschur, viewer);
231: } else {
232: PetscViewerASCIIPrintf(viewer, " not yet available\n");
233: }
234: PetscViewerASCIIPopTab(viewer);
235: PetscViewerASCIIPopTab(viewer);
236: } else if (isdraw && jac->head) {
237: PetscDraw draw;
238: PetscReal x, y, w, wd, h;
239: PetscInt cnt = 2;
240: char str[32];
242: PetscViewerDrawGetDraw(viewer, 0, &draw);
243: PetscDrawGetCurrentPoint(draw, &x, &y);
244: if (jac->kspupper != jac->head->ksp) cnt++;
245: w = 2 * PetscMin(1.0 - x, x);
246: wd = w / (cnt + 1);
248: PetscSNPrintf(str, 32, "Schur fact. %s", PCFieldSplitSchurFactTypes[jac->schurfactorization]);
249: PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h);
250: y -= h;
251: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_USER && !jac->schur_user) {
252: PetscSNPrintf(str, 32, "Prec. for Schur from %s", PCFieldSplitSchurPreTypes[PC_FIELDSPLIT_SCHUR_PRE_A11]);
253: } else {
254: PetscSNPrintf(str, 32, "Prec. for Schur from %s", PCFieldSplitSchurPreTypes[jac->schurpre]);
255: }
256: PetscDrawStringBoxed(draw, x + wd * (cnt - 1) / 2.0, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h);
257: y -= h;
258: x = x - wd * (cnt - 1) / 2.0;
260: PetscDrawPushCurrentPoint(draw, x, y);
261: KSPView(jac->head->ksp, viewer);
262: PetscDrawPopCurrentPoint(draw);
263: if (jac->kspupper != jac->head->ksp) {
264: x += wd;
265: PetscDrawPushCurrentPoint(draw, x, y);
266: KSPView(jac->kspupper, viewer);
267: PetscDrawPopCurrentPoint(draw);
268: }
269: x += wd;
270: PetscDrawPushCurrentPoint(draw, x, y);
271: KSPView(jac->kspschur, viewer);
272: PetscDrawPopCurrentPoint(draw);
273: }
274: return 0;
275: }
277: static PetscErrorCode PCView_FieldSplit_GKB(PC pc, PetscViewer viewer)
278: {
279: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
280: PetscBool iascii, isdraw;
281: PetscInt i, j;
282: PC_FieldSplitLink ilink = jac->head;
284: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
285: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
286: if (iascii) {
287: if (jac->bs > 0) {
288: PetscViewerASCIIPrintf(viewer, " FieldSplit with %s composition: total splits = %" PetscInt_FMT ", blocksize = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits, jac->bs);
289: } else {
290: PetscViewerASCIIPrintf(viewer, " FieldSplit with %s composition: total splits = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits);
291: }
292: if (pc->useAmat) PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for blocks\n");
293: if (jac->diag_use_amat) PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for diagonal blocks\n");
294: if (jac->offdiag_use_amat) PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for off-diagonal blocks\n");
296: PetscViewerASCIIPrintf(viewer, " Stopping tolerance=%.1e, delay in error estimate=%" PetscInt_FMT ", maximum iterations=%" PetscInt_FMT "\n", (double)jac->gkbtol, jac->gkbdelay, jac->gkbmaxit);
297: PetscViewerASCIIPrintf(viewer, " Solver info for H = A00 + nu*A01*A01' matrix:\n");
298: PetscViewerASCIIPushTab(viewer);
300: if (ilink->fields) {
301: PetscViewerASCIIPrintf(viewer, "Split number 0 Fields ");
302: PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
303: for (j = 0; j < ilink->nfields; j++) {
304: if (j > 0) PetscViewerASCIIPrintf(viewer, ",");
305: PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, ilink->fields[j]);
306: }
307: PetscViewerASCIIPrintf(viewer, "\n");
308: PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
309: } else {
310: PetscViewerASCIIPrintf(viewer, "Split number 0 Defined by IS\n");
311: }
312: KSPView(ilink->ksp, viewer);
314: PetscViewerASCIIPopTab(viewer);
315: }
317: if (isdraw) {
318: PetscDraw draw;
319: PetscReal x, y, w, wd;
321: PetscViewerDrawGetDraw(viewer, 0, &draw);
322: PetscDrawGetCurrentPoint(draw, &x, &y);
323: w = 2 * PetscMin(1.0 - x, x);
324: wd = w / (jac->nsplits + 1);
325: x = x - wd * (jac->nsplits - 1) / 2.0;
326: for (i = 0; i < jac->nsplits; i++) {
327: PetscDrawPushCurrentPoint(draw, x, y);
328: KSPView(ilink->ksp, viewer);
329: PetscDrawPopCurrentPoint(draw);
330: x += wd;
331: ilink = ilink->next;
332: }
333: }
334: return 0;
335: }
337: /* Precondition: jac->bs is set to a meaningful value */
338: static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc)
339: {
340: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
341: PetscInt i, nfields, *ifields, nfields_col, *ifields_col;
342: PetscBool flg, flg_col;
343: char optionname[128], splitname[8], optionname_col[128];
345: PetscMalloc1(jac->bs, &ifields);
346: PetscMalloc1(jac->bs, &ifields_col);
347: for (i = 0, flg = PETSC_TRUE;; i++) {
348: PetscSNPrintf(splitname, sizeof(splitname), "%" PetscInt_FMT, i);
349: PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%" PetscInt_FMT "_fields", i);
350: PetscSNPrintf(optionname_col, sizeof(optionname_col), "-pc_fieldsplit_%" PetscInt_FMT "_fields_col", i);
351: nfields = jac->bs;
352: nfields_col = jac->bs;
353: PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg);
354: PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, optionname_col, ifields_col, &nfields_col, &flg_col);
355: if (!flg) break;
356: else if (flg && !flg_col) {
358: PCFieldSplitSetFields(pc, splitname, nfields, ifields, ifields);
359: } else {
362: PCFieldSplitSetFields(pc, splitname, nfields, ifields, ifields_col);
363: }
364: }
365: if (i > 0) {
366: /* Makes command-line setting of splits take precedence over setting them in code.
367: Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would
368: create new splits, which would probably not be what the user wanted. */
369: jac->splitdefined = PETSC_TRUE;
370: }
371: PetscFree(ifields);
372: PetscFree(ifields_col);
373: return 0;
374: }
376: static PetscErrorCode PCFieldSplitSetDefaults(PC pc)
377: {
378: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
379: PC_FieldSplitLink ilink = jac->head;
380: PetscBool fieldsplit_default = PETSC_FALSE, coupling = PETSC_FALSE;
381: PetscInt i;
383: /*
384: Kinda messy, but at least this now uses DMCreateFieldDecomposition().
385: Should probably be rewritten.
386: */
387: if (!ilink) {
388: PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_fieldsplit_detect_coupling", &coupling, NULL);
389: if (pc->dm && jac->dm_splits && !jac->detect && !coupling) {
390: PetscInt numFields, f, i, j;
391: char **fieldNames;
392: IS *fields;
393: DM *dms;
394: DM subdm[128];
395: PetscBool flg;
397: DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms);
398: /* Allow the user to prescribe the splits */
399: for (i = 0, flg = PETSC_TRUE;; i++) {
400: PetscInt ifields[128];
401: IS compField;
402: char optionname[128], splitname[8];
403: PetscInt nfields = numFields;
405: PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%" PetscInt_FMT "_fields", i);
406: PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg);
407: if (!flg) break;
409: DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i]);
410: if (nfields == 1) {
411: PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField);
412: } else {
413: PetscSNPrintf(splitname, sizeof(splitname), "%" PetscInt_FMT, i);
414: PCFieldSplitSetIS(pc, splitname, compField);
415: }
416: ISDestroy(&compField);
417: for (j = 0; j < nfields; ++j) {
418: f = ifields[j];
419: PetscFree(fieldNames[f]);
420: ISDestroy(&fields[f]);
421: }
422: }
423: if (i == 0) {
424: for (f = 0; f < numFields; ++f) {
425: PCFieldSplitSetIS(pc, fieldNames[f], fields[f]);
426: PetscFree(fieldNames[f]);
427: ISDestroy(&fields[f]);
428: }
429: } else {
430: for (j = 0; j < numFields; j++) DMDestroy(dms + j);
431: PetscFree(dms);
432: PetscMalloc1(i, &dms);
433: for (j = 0; j < i; ++j) dms[j] = subdm[j];
434: }
435: PetscFree(fieldNames);
436: PetscFree(fields);
437: if (dms) {
438: PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n");
439: for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) {
440: const char *prefix;
441: PetscObjectGetOptionsPrefix((PetscObject)(ilink->ksp), &prefix);
442: PetscObjectSetOptionsPrefix((PetscObject)(dms[i]), prefix);
443: KSPSetDM(ilink->ksp, dms[i]);
444: KSPSetDMActive(ilink->ksp, PETSC_FALSE);
445: {
446: PetscErrorCode (*func)(KSP, Mat, Mat, void *);
447: void *ctx;
449: DMKSPGetComputeOperators(pc->dm, &func, &ctx);
450: DMKSPSetComputeOperators(dms[i], func, ctx);
451: }
452: PetscObjectIncrementTabLevel((PetscObject)dms[i], (PetscObject)ilink->ksp, 0);
453: DMDestroy(&dms[i]);
454: }
455: PetscFree(dms);
456: }
457: } else {
458: if (jac->bs <= 0) {
459: if (pc->pmat) {
460: MatGetBlockSize(pc->pmat, &jac->bs);
461: } else jac->bs = 1;
462: }
464: if (jac->detect) {
465: IS zerodiags, rest;
466: PetscInt nmin, nmax;
468: MatGetOwnershipRange(pc->mat, &nmin, &nmax);
469: if (jac->diag_use_amat) {
470: MatFindZeroDiagonals(pc->mat, &zerodiags);
471: } else {
472: MatFindZeroDiagonals(pc->pmat, &zerodiags);
473: }
474: ISComplement(zerodiags, nmin, nmax, &rest);
475: PCFieldSplitSetIS(pc, "0", rest);
476: PCFieldSplitSetIS(pc, "1", zerodiags);
477: ISDestroy(&zerodiags);
478: ISDestroy(&rest);
479: } else if (coupling) {
480: IS coupling, rest;
481: PetscInt nmin, nmax;
483: MatGetOwnershipRange(pc->mat, &nmin, &nmax);
484: if (jac->offdiag_use_amat) {
485: MatFindOffBlockDiagonalEntries(pc->mat, &coupling);
486: } else {
487: MatFindOffBlockDiagonalEntries(pc->pmat, &coupling);
488: }
489: ISCreateStride(PetscObjectComm((PetscObject)pc->mat), nmax - nmin, nmin, 1, &rest);
490: ISSetIdentity(rest);
491: PCFieldSplitSetIS(pc, "0", rest);
492: PCFieldSplitSetIS(pc, "1", coupling);
493: ISDestroy(&coupling);
494: ISDestroy(&rest);
495: } else {
496: 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: PCFieldSplitSetRuntimeSplits_Private(pc);
501: if (jac->splitdefined) 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: PetscInfo(pc, "Using default splitting of fields\n");
508: PetscObjectTypeCompare((PetscObject)pc->pmat, MATNEST, &isnest);
509: if (!isnest) {
510: M = pc->mat;
511: PetscObjectTypeCompare((PetscObject)pc->mat, MATNEST, &isnest);
512: }
513: if (isnest) {
514: IS *fields;
515: PetscInt nf;
517: MatNestGetSize(M, &nf, NULL);
518: PetscMalloc1(nf, &fields);
519: MatNestGetISs(M, fields, NULL);
520: for (i = 0; i < nf; i++) PCFieldSplitSetIS(pc, NULL, fields[i]);
521: PetscFree(fields);
522: } else {
523: for (i = 0; i < jac->bs; i++) {
524: char splitname[8];
525: PetscSNPrintf(splitname, sizeof(splitname), "%" PetscInt_FMT, i);
526: PCFieldSplitSetFields(pc, splitname, 1, &i, &i);
527: }
528: jac->defaultsplit = PETSC_TRUE;
529: }
530: }
531: }
532: }
533: } else if (jac->nsplits == 1) {
534: if (ilink->is) {
535: IS is2;
536: PetscInt nmin, nmax;
538: MatGetOwnershipRange(pc->mat, &nmin, &nmax);
539: ISComplement(ilink->is, nmin, nmax, &is2);
540: PCFieldSplitSetIS(pc, "1", is2);
541: ISDestroy(&is2);
542: } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Must provide at least two sets of fields to PCFieldSplit()");
543: }
546: return 0;
547: }
549: static PetscErrorCode MatGolubKahanComputeExplicitOperator(Mat A, Mat B, Mat C, Mat *H, PetscReal gkbnu)
550: {
551: Mat BT, T;
552: PetscReal nrmT, nrmB;
554: MatHermitianTranspose(C, MAT_INITIAL_MATRIX, &T); /* Test if augmented matrix is symmetric */
555: MatAXPY(T, -1.0, B, DIFFERENT_NONZERO_PATTERN);
556: MatNorm(T, NORM_1, &nrmT);
557: MatNorm(B, NORM_1, &nrmB);
560: /* Compute augmented Lagrangian matrix H = A00 + nu*A01*A01'. This corresponds to */
561: /* setting N := 1/nu*I in [Ar13]. */
562: MatHermitianTranspose(B, MAT_INITIAL_MATRIX, &BT);
563: MatMatMult(B, BT, MAT_INITIAL_MATRIX, PETSC_DEFAULT, H); /* H = A01*A01' */
564: MatAYPX(*H, gkbnu, A, DIFFERENT_NONZERO_PATTERN); /* H = A00 + nu*A01*A01' */
566: MatDestroy(&BT);
567: MatDestroy(&T);
568: return 0;
569: }
571: PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(PetscOptions, const char pre[], const char name[], 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: pc->failedreason = PC_NOERROR;
581: PCFieldSplitSetDefaults(pc);
582: nsplit = jac->nsplits;
583: ilink = jac->head;
585: /* get the matrices for each split */
586: if (!jac->issetup) {
587: PetscInt rstart, rend, nslots, bs;
589: jac->issetup = PETSC_TRUE;
591: /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */
592: if (jac->defaultsplit || !ilink->is) {
593: if (jac->bs <= 0) jac->bs = nsplit;
594: }
596: /* MatCreateSubMatrix() for [S]BAIJ matrices can only work if the indices include entire blocks of the matrix */
597: MatGetBlockSize(pc->pmat, &bs);
598: if (bs > 1 && (jac->bs <= bs || jac->bs % bs)) {
599: PetscBool blk;
601: PetscObjectTypeCompareAny((PetscObject)pc->pmat, &blk, MATBAIJ, MATSBAIJ, MATSEQBAIJ, MATSEQSBAIJ, MATMPIBAIJ, MATMPISBAIJ, NULL);
603: }
605: bs = jac->bs;
606: MatGetOwnershipRange(pc->pmat, &rstart, &rend);
607: nslots = (rend - rstart) / bs;
608: for (i = 0; i < nsplit; i++) {
609: if (jac->defaultsplit) {
610: ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + i, nsplit, &ilink->is);
611: ISDuplicate(ilink->is, &ilink->is_col);
612: } else if (!ilink->is) {
613: if (ilink->nfields > 1) {
614: PetscInt *ii, *jj, j, k, nfields = ilink->nfields, *fields = ilink->fields, *fields_col = ilink->fields_col;
615: PetscMalloc1(ilink->nfields * nslots, &ii);
616: PetscMalloc1(ilink->nfields * nslots, &jj);
617: for (j = 0; j < nslots; j++) {
618: for (k = 0; k < nfields; k++) {
619: ii[nfields * j + k] = rstart + bs * j + fields[k];
620: jj[nfields * j + k] = rstart + bs * j + fields_col[k];
621: }
622: }
623: ISCreateGeneral(PetscObjectComm((PetscObject)pc), nslots * nfields, ii, PETSC_OWN_POINTER, &ilink->is);
624: ISCreateGeneral(PetscObjectComm((PetscObject)pc), nslots * nfields, jj, PETSC_OWN_POINTER, &ilink->is_col);
625: ISSetBlockSize(ilink->is, nfields);
626: ISSetBlockSize(ilink->is_col, nfields);
627: } else {
628: ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + ilink->fields[0], bs, &ilink->is);
629: ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + ilink->fields_col[0], bs, &ilink->is_col);
630: }
631: }
632: ISSorted(ilink->is, &sorted);
633: if (ilink->is_col) ISSorted(ilink->is_col, &sorted_col);
635: ilink = ilink->next;
636: }
637: }
639: ilink = jac->head;
640: if (!jac->pmat) {
641: Vec xtmp;
643: MatCreateVecs(pc->pmat, &xtmp, NULL);
644: PetscMalloc1(nsplit, &jac->pmat);
645: PetscMalloc2(nsplit, &jac->x, nsplit, &jac->y);
646: for (i = 0; i < nsplit; i++) {
647: MatNullSpace sp;
649: /* Check for preconditioning matrix attached to IS */
650: PetscObjectQuery((PetscObject)ilink->is, "pmat", (PetscObject *)&jac->pmat[i]);
651: if (jac->pmat[i]) {
652: PetscObjectReference((PetscObject)jac->pmat[i]);
653: if (jac->type == PC_COMPOSITE_SCHUR) {
654: jac->schur_user = jac->pmat[i];
656: PetscObjectReference((PetscObject)jac->schur_user);
657: }
658: } else {
659: const char *prefix;
660: MatCreateSubMatrix(pc->pmat, ilink->is, ilink->is_col, MAT_INITIAL_MATRIX, &jac->pmat[i]);
661: KSPGetOptionsPrefix(ilink->ksp, &prefix);
662: MatSetOptionsPrefix(jac->pmat[i], prefix);
663: MatViewFromOptions(jac->pmat[i], NULL, "-mat_view");
664: }
665: /* create work vectors for each split */
666: MatCreateVecs(jac->pmat[i], &jac->x[i], &jac->y[i]);
667: ilink->x = jac->x[i];
668: ilink->y = jac->y[i];
669: ilink->z = NULL;
670: /* compute scatter contexts needed by multiplicative versions and non-default splits */
671: VecScatterCreate(xtmp, ilink->is, jac->x[i], NULL, &ilink->sctx);
672: PetscObjectQuery((PetscObject)ilink->is, "nearnullspace", (PetscObject *)&sp);
673: if (sp) MatSetNearNullSpace(jac->pmat[i], sp);
674: ilink = ilink->next;
675: }
676: VecDestroy(&xtmp);
677: } else {
678: MatReuse scall;
679: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
680: for (i = 0; i < nsplit; i++) MatDestroy(&jac->pmat[i]);
681: scall = MAT_INITIAL_MATRIX;
682: } else scall = MAT_REUSE_MATRIX;
684: for (i = 0; i < nsplit; i++) {
685: Mat pmat;
687: /* Check for preconditioning matrix attached to IS */
688: PetscObjectQuery((PetscObject)ilink->is, "pmat", (PetscObject *)&pmat);
689: if (!pmat) MatCreateSubMatrix(pc->pmat, ilink->is, ilink->is_col, scall, &jac->pmat[i]);
690: ilink = ilink->next;
691: }
692: }
693: if (jac->diag_use_amat) {
694: ilink = jac->head;
695: if (!jac->mat) {
696: PetscMalloc1(nsplit, &jac->mat);
697: for (i = 0; i < nsplit; i++) {
698: MatCreateSubMatrix(pc->mat, ilink->is, ilink->is_col, MAT_INITIAL_MATRIX, &jac->mat[i]);
699: ilink = ilink->next;
700: }
701: } else {
702: MatReuse scall;
703: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
704: for (i = 0; i < nsplit; i++) MatDestroy(&jac->mat[i]);
705: scall = MAT_INITIAL_MATRIX;
706: } else scall = MAT_REUSE_MATRIX;
708: for (i = 0; i < nsplit; i++) {
709: MatCreateSubMatrix(pc->mat, ilink->is, ilink->is_col, scall, &jac->mat[i]);
710: ilink = ilink->next;
711: }
712: }
713: } else {
714: jac->mat = jac->pmat;
715: }
717: /* Check for null space attached to IS */
718: ilink = jac->head;
719: for (i = 0; i < nsplit; i++) {
720: MatNullSpace sp;
722: PetscObjectQuery((PetscObject)ilink->is, "nullspace", (PetscObject *)&sp);
723: if (sp) MatSetNullSpace(jac->mat[i], sp);
724: ilink = ilink->next;
725: }
727: if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR && jac->type != PC_COMPOSITE_GKB) {
728: /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */
729: /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */
730: ilink = jac->head;
731: if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) {
732: /* 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 */
733: if (!jac->Afield) {
734: PetscCalloc1(nsplit, &jac->Afield);
735: if (jac->offdiag_use_amat) {
736: MatCreateSubMatrix(pc->mat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->Afield[1]);
737: } else {
738: MatCreateSubMatrix(pc->pmat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->Afield[1]);
739: }
740: } else {
741: MatReuse scall;
743: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
744: MatDestroy(&jac->Afield[1]);
745: scall = MAT_INITIAL_MATRIX;
746: } else scall = MAT_REUSE_MATRIX;
748: if (jac->offdiag_use_amat) {
749: MatCreateSubMatrix(pc->mat, ilink->next->is, ilink->is, scall, &jac->Afield[1]);
750: } else {
751: MatCreateSubMatrix(pc->pmat, ilink->next->is, ilink->is, scall, &jac->Afield[1]);
752: }
753: }
754: } else {
755: if (!jac->Afield) {
756: PetscMalloc1(nsplit, &jac->Afield);
757: for (i = 0; i < nsplit; i++) {
758: if (jac->offdiag_use_amat) {
759: MatCreateSubMatrix(pc->mat, ilink->is, NULL, MAT_INITIAL_MATRIX, &jac->Afield[i]);
760: } else {
761: MatCreateSubMatrix(pc->pmat, ilink->is, NULL, MAT_INITIAL_MATRIX, &jac->Afield[i]);
762: }
763: ilink = ilink->next;
764: }
765: } else {
766: MatReuse scall;
767: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
768: for (i = 0; i < nsplit; i++) MatDestroy(&jac->Afield[i]);
769: scall = MAT_INITIAL_MATRIX;
770: } else scall = MAT_REUSE_MATRIX;
772: for (i = 0; i < nsplit; i++) {
773: if (jac->offdiag_use_amat) {
774: MatCreateSubMatrix(pc->mat, ilink->is, NULL, scall, &jac->Afield[i]);
775: } else {
776: MatCreateSubMatrix(pc->pmat, ilink->is, NULL, scall, &jac->Afield[i]);
777: }
778: ilink = ilink->next;
779: }
780: }
781: }
782: }
784: if (jac->type == PC_COMPOSITE_SCHUR) {
785: IS ccis;
786: PetscBool isset, isspd;
787: PetscInt rstart, rend;
788: char lscname[256];
789: PetscObject LSC_L;
793: /* If pc->mat is SPD, don't scale by -1 the Schur complement */
794: if (jac->schurscale == (PetscScalar)-1.0) {
795: MatIsSPDKnown(pc->pmat, &isset, &isspd);
796: jac->schurscale = (isset && isspd) ? 1.0 : -1.0;
797: }
799: /* When extracting off-diagonal submatrices, we take complements from this range */
800: MatGetOwnershipRangeColumn(pc->mat, &rstart, &rend);
802: if (jac->schur) {
803: KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper;
804: MatReuse scall;
806: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
807: scall = MAT_INITIAL_MATRIX;
808: MatDestroy(&jac->B);
809: MatDestroy(&jac->C);
810: } else scall = MAT_REUSE_MATRIX;
812: MatSchurComplementGetKSP(jac->schur, &kspInner);
813: ilink = jac->head;
814: ISComplement(ilink->is_col, rstart, rend, &ccis);
815: if (jac->offdiag_use_amat) {
816: MatCreateSubMatrix(pc->mat, ilink->is, ccis, scall, &jac->B);
817: } else {
818: MatCreateSubMatrix(pc->pmat, ilink->is, ccis, scall, &jac->B);
819: }
820: ISDestroy(&ccis);
821: ilink = ilink->next;
822: ISComplement(ilink->is_col, rstart, rend, &ccis);
823: if (jac->offdiag_use_amat) {
824: MatCreateSubMatrix(pc->mat, ilink->is, ccis, scall, &jac->C);
825: } else {
826: MatCreateSubMatrix(pc->pmat, ilink->is, ccis, scall, &jac->C);
827: }
828: ISDestroy(&ccis);
829: MatSchurComplementUpdateSubMatrices(jac->schur, jac->mat[0], jac->pmat[0], jac->B, jac->C, jac->mat[1]);
830: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
831: MatDestroy(&jac->schurp);
832: MatSchurComplementGetPmat(jac->schur, MAT_INITIAL_MATRIX, &jac->schurp);
833: }
834: if (kspA != kspInner) KSPSetOperators(kspA, jac->mat[0], jac->pmat[0]);
835: if (kspUpper != kspA) KSPSetOperators(kspUpper, jac->mat[0], jac->pmat[0]);
836: KSPSetOperators(jac->kspschur, jac->schur, FieldSplitSchurPre(jac));
837: } else {
838: const char *Dprefix;
839: char schurprefix[256], schurmatprefix[256];
840: char schurtestoption[256];
841: MatNullSpace sp;
842: PetscBool flg;
843: KSP kspt;
845: /* extract the A01 and A10 matrices */
846: ilink = jac->head;
847: ISComplement(ilink->is_col, rstart, rend, &ccis);
848: if (jac->offdiag_use_amat) {
849: MatCreateSubMatrix(pc->mat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->B);
850: } else {
851: MatCreateSubMatrix(pc->pmat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->B);
852: }
853: ISDestroy(&ccis);
854: ilink = ilink->next;
855: ISComplement(ilink->is_col, rstart, rend, &ccis);
856: if (jac->offdiag_use_amat) {
857: MatCreateSubMatrix(pc->mat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->C);
858: } else {
859: MatCreateSubMatrix(pc->pmat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->C);
860: }
861: ISDestroy(&ccis);
863: /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
864: MatCreate(((PetscObject)jac->mat[0])->comm, &jac->schur);
865: MatSetType(jac->schur, MATSCHURCOMPLEMENT);
866: MatSchurComplementSetSubMatrices(jac->schur, jac->mat[0], jac->pmat[0], jac->B, jac->C, jac->mat[1]);
867: PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
868: MatSetOptionsPrefix(jac->schur, schurmatprefix);
869: MatSchurComplementGetKSP(jac->schur, &kspt);
870: KSPSetOptionsPrefix(kspt, schurmatprefix);
872: /* Note: this is not true in general */
873: MatGetNullSpace(jac->mat[1], &sp);
874: if (sp) MatSetNullSpace(jac->schur, sp);
876: PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);
877: PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
878: if (flg) {
879: DM dmInner;
880: KSP kspInner;
881: PC pcInner;
883: MatSchurComplementGetKSP(jac->schur, &kspInner);
884: KSPReset(kspInner);
885: KSPSetOperators(kspInner, jac->mat[0], jac->pmat[0]);
886: PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
887: /* Indent this deeper to emphasize the "inner" nature of this solver. */
888: PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject)pc, 2);
889: PetscObjectIncrementTabLevel((PetscObject)kspInner->pc, (PetscObject)pc, 2);
890: KSPSetOptionsPrefix(kspInner, schurprefix);
892: /* Set DM for new solver */
893: KSPGetDM(jac->head->ksp, &dmInner);
894: KSPSetDM(kspInner, dmInner);
895: KSPSetDMActive(kspInner, PETSC_FALSE);
897: /* Defaults to PCKSP as preconditioner */
898: KSPGetPC(kspInner, &pcInner);
899: PCSetType(pcInner, PCKSP);
900: PCKSPSetKSP(pcInner, jac->head->ksp);
901: } else {
902: /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or
903: * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact,
904: * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for
905: * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make
906: * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used
907: * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */
908: KSPSetType(jac->head->ksp, KSPGMRES);
909: MatSchurComplementSetKSP(jac->schur, jac->head->ksp);
910: }
911: KSPSetOperators(jac->head->ksp, jac->mat[0], jac->pmat[0]);
912: KSPSetFromOptions(jac->head->ksp);
913: MatSetFromOptions(jac->schur);
915: PetscObjectTypeCompare((PetscObject)jac->schur, MATSCHURCOMPLEMENT, &flg);
916: if (flg) { /* Need to do this otherwise PCSetUp_KSP will overwrite the amat of jac->head->ksp */
917: KSP kspInner;
918: PC pcInner;
920: MatSchurComplementGetKSP(jac->schur, &kspInner);
921: KSPGetPC(kspInner, &pcInner);
922: PetscObjectTypeCompare((PetscObject)pcInner, PCKSP, &flg);
923: if (flg) {
924: KSP ksp;
926: PCKSPGetKSP(pcInner, &ksp);
927: if (ksp == jac->head->ksp) PCSetUseAmat(pcInner, PETSC_TRUE);
928: }
929: }
930: PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);
931: PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
932: if (flg) {
933: DM dmInner;
935: PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
936: KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);
937: KSPSetErrorIfNotConverged(jac->kspupper, pc->erroriffailure);
938: KSPSetOptionsPrefix(jac->kspupper, schurprefix);
939: PetscObjectIncrementTabLevel((PetscObject)jac->kspupper, (PetscObject)pc, 1);
940: PetscObjectIncrementTabLevel((PetscObject)jac->kspupper->pc, (PetscObject)pc, 1);
941: KSPGetDM(jac->head->ksp, &dmInner);
942: KSPSetDM(jac->kspupper, dmInner);
943: KSPSetDMActive(jac->kspupper, PETSC_FALSE);
944: KSPSetFromOptions(jac->kspupper);
945: KSPSetOperators(jac->kspupper, jac->mat[0], jac->pmat[0]);
946: VecDuplicate(jac->head->x, &jac->head->z);
947: } else {
948: jac->kspupper = jac->head->ksp;
949: PetscObjectReference((PetscObject)jac->head->ksp);
950: }
952: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) MatSchurComplementGetPmat(jac->schur, MAT_INITIAL_MATRIX, &jac->schurp);
953: KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspschur);
954: KSPSetErrorIfNotConverged(jac->kspschur, pc->erroriffailure);
955: PetscObjectIncrementTabLevel((PetscObject)jac->kspschur, (PetscObject)pc, 1);
956: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
957: PC pcschur;
958: KSPGetPC(jac->kspschur, &pcschur);
959: PCSetType(pcschur, PCNONE);
960: /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
961: } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) {
962: MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user);
963: }
964: KSPSetOperators(jac->kspschur, jac->schur, FieldSplitSchurPre(jac));
965: KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);
966: KSPSetOptionsPrefix(jac->kspschur, Dprefix);
967: /* propagate DM */
968: {
969: DM sdm;
970: KSPGetDM(jac->head->next->ksp, &sdm);
971: if (sdm) {
972: KSPSetDM(jac->kspschur, sdm);
973: KSPSetDMActive(jac->kspschur, PETSC_FALSE);
974: }
975: }
976: /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
977: /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */
978: KSPSetFromOptions(jac->kspschur);
979: }
980: MatAssemblyBegin(jac->schur, MAT_FINAL_ASSEMBLY);
981: MatAssemblyEnd(jac->schur, MAT_FINAL_ASSEMBLY);
983: /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */
984: PetscSNPrintf(lscname, sizeof(lscname), "%s_LSC_L", ilink->splitname);
985: PetscObjectQuery((PetscObject)pc->mat, lscname, (PetscObject *)&LSC_L);
986: if (!LSC_L) PetscObjectQuery((PetscObject)pc->pmat, lscname, (PetscObject *)&LSC_L);
987: if (LSC_L) PetscObjectCompose((PetscObject)jac->schur, "LSC_L", (PetscObject)LSC_L);
988: PetscSNPrintf(lscname, sizeof(lscname), "%s_LSC_Lp", ilink->splitname);
989: PetscObjectQuery((PetscObject)pc->pmat, lscname, (PetscObject *)&LSC_L);
990: if (!LSC_L) PetscObjectQuery((PetscObject)pc->mat, lscname, (PetscObject *)&LSC_L);
991: if (LSC_L) PetscObjectCompose((PetscObject)jac->schur, "LSC_Lp", (PetscObject)LSC_L);
992: } else if (jac->type == PC_COMPOSITE_GKB) {
993: IS ccis;
994: PetscInt rstart, rend;
998: ilink = jac->head;
1000: /* When extracting off-diagonal submatrices, we take complements from this range */
1001: MatGetOwnershipRangeColumn(pc->mat, &rstart, &rend);
1003: ISComplement(ilink->is_col, rstart, rend, &ccis);
1004: if (jac->offdiag_use_amat) {
1005: MatCreateSubMatrix(pc->mat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->B);
1006: } else {
1007: MatCreateSubMatrix(pc->pmat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->B);
1008: }
1009: ISDestroy(&ccis);
1010: /* Create work vectors for GKB algorithm */
1011: VecDuplicate(ilink->x, &jac->u);
1012: VecDuplicate(ilink->x, &jac->Hu);
1013: VecDuplicate(ilink->x, &jac->w2);
1014: ilink = ilink->next;
1015: ISComplement(ilink->is_col, rstart, rend, &ccis);
1016: if (jac->offdiag_use_amat) {
1017: MatCreateSubMatrix(pc->mat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->C);
1018: } else {
1019: MatCreateSubMatrix(pc->pmat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->C);
1020: }
1021: ISDestroy(&ccis);
1022: /* Create work vectors for GKB algorithm */
1023: VecDuplicate(ilink->x, &jac->v);
1024: VecDuplicate(ilink->x, &jac->d);
1025: VecDuplicate(ilink->x, &jac->w1);
1026: MatGolubKahanComputeExplicitOperator(jac->mat[0], jac->B, jac->C, &jac->H, jac->gkbnu);
1027: PetscCalloc1(jac->gkbdelay, &jac->vecz);
1029: ilink = jac->head;
1030: KSPSetOperators(ilink->ksp, jac->H, jac->H);
1031: if (!jac->suboptionsset) KSPSetFromOptions(ilink->ksp);
1032: /* Create gkb_monitor context */
1033: if (jac->gkbmonitor) {
1034: PetscInt tablevel;
1035: PetscViewerCreate(PETSC_COMM_WORLD, &jac->gkbviewer);
1036: PetscViewerSetType(jac->gkbviewer, PETSCVIEWERASCII);
1037: PetscObjectGetTabLevel((PetscObject)ilink->ksp, &tablevel);
1038: PetscViewerASCIISetTab(jac->gkbviewer, tablevel);
1039: PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)ilink->ksp, 1);
1040: }
1041: } else {
1042: /* set up the individual splits' PCs */
1043: i = 0;
1044: ilink = jac->head;
1045: while (ilink) {
1046: KSPSetOperators(ilink->ksp, jac->mat[i], jac->pmat[i]);
1047: /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
1048: if (!jac->suboptionsset) KSPSetFromOptions(ilink->ksp);
1049: i++;
1050: ilink = ilink->next;
1051: }
1052: }
1054: /* Set coordinates to the sub PC objects whenever these are set */
1055: if (jac->coordinates_set) {
1056: PC pc_coords;
1057: if (jac->type == PC_COMPOSITE_SCHUR) {
1058: // Head is first block.
1059: KSPGetPC(jac->head->ksp, &pc_coords);
1060: PCSetCoordinates(pc_coords, jac->head->dim, jac->head->ndofs, jac->head->coords);
1061: // Second one is Schur block, but its KSP object is in kspschur.
1062: KSPGetPC(jac->kspschur, &pc_coords);
1063: PCSetCoordinates(pc_coords, jac->head->next->dim, jac->head->next->ndofs, jac->head->next->coords);
1064: } else if (jac->type == PC_COMPOSITE_GKB) {
1065: PetscInfo(pc, "Warning: Setting coordinates does nothing for the GKB Fieldpslit preconditioner");
1066: } else {
1067: ilink = jac->head;
1068: while (ilink) {
1069: KSPGetPC(ilink->ksp, &pc_coords);
1070: PCSetCoordinates(pc_coords, ilink->dim, ilink->ndofs, ilink->coords);
1071: ilink = ilink->next;
1072: }
1073: }
1074: }
1076: jac->suboptionsset = PETSC_TRUE;
1077: return 0;
1078: }
1080: #define FieldSplitSplitSolveAdd(ilink, xx, yy) \
1081: (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) || \
1082: 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) || \
1083: VecScatterEnd(ilink->sctx, ilink->y, yy, ADD_VALUES, SCATTER_REVERSE))
1085: static PetscErrorCode PCApply_FieldSplit_Schur(PC pc, Vec x, Vec y)
1086: {
1087: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1088: PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
1089: KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;
1091: switch (jac->schurfactorization) {
1092: case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
1093: /* [A00 0; 0 -S], positive definite, suitable for MINRES */
1094: VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD);
1095: VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD);
1096: VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD);
1097: PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL);
1098: KSPSolve(kspA, ilinkA->x, ilinkA->y);
1099: KSPCheckSolve(kspA, pc, ilinkA->y);
1100: PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL);
1101: VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE);
1102: VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD);
1103: PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL);
1104: KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y);
1105: KSPCheckSolve(jac->kspschur, pc, ilinkD->y);
1106: PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL);
1107: VecScale(ilinkD->y, jac->schurscale);
1108: VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE);
1109: VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE);
1110: VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE);
1111: break;
1112: case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
1113: /* [A00 0; A10 S], suitable for left preconditioning */
1114: VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD);
1115: VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD);
1116: PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL);
1117: KSPSolve(kspA, ilinkA->x, ilinkA->y);
1118: KSPCheckSolve(kspA, pc, ilinkA->y);
1119: PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL);
1120: MatMult(jac->C, ilinkA->y, ilinkD->x);
1121: VecScale(ilinkD->x, -1.);
1122: VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD);
1123: VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE);
1124: VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD);
1125: PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL);
1126: KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y);
1127: KSPCheckSolve(jac->kspschur, pc, ilinkD->y);
1128: PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL);
1129: VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE);
1130: VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE);
1131: VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE);
1132: break;
1133: case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
1134: /* [A00 A01; 0 S], suitable for right preconditioning */
1135: VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD);
1136: VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD);
1137: PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL);
1138: KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y);
1139: KSPCheckSolve(jac->kspschur, pc, ilinkD->y);
1140: PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL);
1141: MatMult(jac->B, ilinkD->y, ilinkA->x);
1142: VecScale(ilinkA->x, -1.);
1143: VecScatterBegin(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD);
1144: VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE);
1145: VecScatterEnd(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD);
1146: PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL);
1147: KSPSolve(kspA, ilinkA->x, ilinkA->y);
1148: KSPCheckSolve(kspA, pc, ilinkA->y);
1149: PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL);
1150: VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE);
1151: VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE);
1152: VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE);
1153: break;
1154: case PC_FIELDSPLIT_SCHUR_FACT_FULL:
1155: /* [1 0; A10 A00^{-1} 1] [A00 0; 0 S] [1 A00^{-1}A01; 0 1] */
1156: VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD);
1157: VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD);
1158: PetscLogEventBegin(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->y, NULL);
1159: KSPSolve(kspLower, ilinkA->x, ilinkA->y);
1160: KSPCheckSolve(kspLower, pc, ilinkA->y);
1161: PetscLogEventEnd(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->y, NULL);
1162: MatMult(jac->C, ilinkA->y, ilinkD->x);
1163: VecScale(ilinkD->x, -1.0);
1164: VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD);
1165: VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD);
1167: PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL);
1168: KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y);
1169: KSPCheckSolve(jac->kspschur, pc, ilinkD->y);
1170: PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL);
1171: VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE);
1173: if (kspUpper == kspA) {
1174: MatMult(jac->B, ilinkD->y, ilinkA->y);
1175: VecAXPY(ilinkA->x, -1.0, ilinkA->y);
1176: PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL);
1177: KSPSolve(kspA, ilinkA->x, ilinkA->y);
1178: KSPCheckSolve(kspA, pc, ilinkA->y);
1179: PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL);
1180: } else {
1181: PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL);
1182: KSPSolve(kspA, ilinkA->x, ilinkA->y);
1183: KSPCheckSolve(kspA, pc, ilinkA->y);
1184: MatMult(jac->B, ilinkD->y, ilinkA->x);
1185: PetscLogEventBegin(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->z, NULL);
1186: KSPSolve(kspUpper, ilinkA->x, ilinkA->z);
1187: KSPCheckSolve(kspUpper, pc, ilinkA->z);
1188: PetscLogEventEnd(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->z, NULL);
1189: VecAXPY(ilinkA->y, -1.0, ilinkA->z);
1190: }
1191: VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE);
1192: VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE);
1193: VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE);
1194: }
1195: return 0;
1196: }
1198: static PetscErrorCode PCApply_FieldSplit(PC pc, Vec x, Vec y)
1199: {
1200: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1201: PC_FieldSplitLink ilink = jac->head;
1202: PetscInt cnt, bs;
1204: if (jac->type == PC_COMPOSITE_ADDITIVE) {
1205: if (jac->defaultsplit) {
1206: VecGetBlockSize(x, &bs);
1208: VecGetBlockSize(y, &bs);
1210: VecStrideGatherAll(x, jac->x, INSERT_VALUES);
1211: while (ilink) {
1212: PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1213: KSPSolve(ilink->ksp, ilink->x, ilink->y);
1214: KSPCheckSolve(ilink->ksp, pc, ilink->y);
1215: PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1216: ilink = ilink->next;
1217: }
1218: VecStrideScatterAll(jac->y, y, INSERT_VALUES);
1219: } else {
1220: VecSet(y, 0.0);
1221: while (ilink) {
1222: FieldSplitSplitSolveAdd(ilink, x, y);
1223: ilink = ilink->next;
1224: }
1225: }
1226: } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) {
1227: VecSet(y, 0.0);
1228: /* solve on first block for first block variables */
1229: VecScatterBegin(ilink->sctx, x, ilink->x, INSERT_VALUES, SCATTER_FORWARD);
1230: VecScatterEnd(ilink->sctx, x, ilink->x, INSERT_VALUES, SCATTER_FORWARD);
1231: PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1232: KSPSolve(ilink->ksp, ilink->x, ilink->y);
1233: KSPCheckSolve(ilink->ksp, pc, ilink->y);
1234: PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1235: VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE);
1236: VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE);
1238: /* compute the residual only onto second block variables using first block variables */
1239: MatMult(jac->Afield[1], ilink->y, ilink->next->x);
1240: ilink = ilink->next;
1241: VecScale(ilink->x, -1.0);
1242: VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD);
1243: VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD);
1245: /* solve on second block variables */
1246: PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1247: KSPSolve(ilink->ksp, ilink->x, ilink->y);
1248: KSPCheckSolve(ilink->ksp, pc, ilink->y);
1249: PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1250: VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE);
1251: VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE);
1252: } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1253: if (!jac->w1) {
1254: VecDuplicate(x, &jac->w1);
1255: VecDuplicate(x, &jac->w2);
1256: }
1257: VecSet(y, 0.0);
1258: FieldSplitSplitSolveAdd(ilink, x, y);
1259: cnt = 1;
1260: while (ilink->next) {
1261: ilink = ilink->next;
1262: /* compute the residual only over the part of the vector needed */
1263: MatMult(jac->Afield[cnt++], y, ilink->x);
1264: VecScale(ilink->x, -1.0);
1265: VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD);
1266: VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD);
1267: PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1268: KSPSolve(ilink->ksp, ilink->x, ilink->y);
1269: KSPCheckSolve(ilink->ksp, pc, ilink->y);
1270: PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1271: VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE);
1272: VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE);
1273: }
1274: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1275: cnt -= 2;
1276: while (ilink->previous) {
1277: ilink = ilink->previous;
1278: /* compute the residual only over the part of the vector needed */
1279: MatMult(jac->Afield[cnt--], y, ilink->x);
1280: VecScale(ilink->x, -1.0);
1281: VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD);
1282: VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD);
1283: PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1284: KSPSolve(ilink->ksp, ilink->x, ilink->y);
1285: KSPCheckSolve(ilink->ksp, pc, ilink->y);
1286: PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1287: VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE);
1288: VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE);
1289: }
1290: }
1291: } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Unsupported or unknown composition %d", (int)jac->type);
1292: return 0;
1293: }
1295: static PetscErrorCode PCApply_FieldSplit_GKB(PC pc, Vec x, Vec y)
1296: {
1297: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1298: PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
1299: KSP ksp = ilinkA->ksp;
1300: Vec u, v, Hu, d, work1, work2;
1301: PetscScalar alpha, z, nrmz2, *vecz;
1302: PetscReal lowbnd, nu, beta;
1303: PetscInt j, iterGKB;
1305: VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD);
1306: VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD);
1307: VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD);
1308: VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD);
1310: u = jac->u;
1311: v = jac->v;
1312: Hu = jac->Hu;
1313: d = jac->d;
1314: work1 = jac->w1;
1315: work2 = jac->w2;
1316: vecz = jac->vecz;
1318: /* Change RHS to comply with matrix regularization H = A + nu*B*B' */
1319: /* Add q = q + nu*B*b */
1320: if (jac->gkbnu) {
1321: nu = jac->gkbnu;
1322: VecScale(ilinkD->x, jac->gkbnu);
1323: MatMultAdd(jac->B, ilinkD->x, ilinkA->x, ilinkA->x); /* q = q + nu*B*b */
1324: } else {
1325: /* Situation when no augmented Lagrangian is used. Then we set inner */
1326: /* matrix N = I in [Ar13], and thus nu = 1. */
1327: nu = 1;
1328: }
1330: /* Transform rhs from [q,tilde{b}] to [0,b] */
1331: PetscLogEventBegin(ilinkA->event, ksp, ilinkA->x, ilinkA->y, NULL);
1332: KSPSolve(ksp, ilinkA->x, ilinkA->y);
1333: KSPCheckSolve(ksp, pc, ilinkA->y);
1334: PetscLogEventEnd(ilinkA->event, ksp, ilinkA->x, ilinkA->y, NULL);
1335: MatMultHermitianTranspose(jac->B, ilinkA->y, work1);
1336: VecAXPBY(work1, 1.0 / nu, -1.0, ilinkD->x); /* c = b - B'*x */
1338: /* First step of algorithm */
1339: VecNorm(work1, NORM_2, &beta); /* beta = sqrt(nu*c'*c)*/
1340: KSPCheckDot(ksp, beta);
1341: beta = PetscSqrtReal(nu) * beta;
1342: VecAXPBY(v, nu / beta, 0.0, work1); /* v = nu/beta *c */
1343: MatMult(jac->B, v, work2); /* u = H^{-1}*B*v */
1344: PetscLogEventBegin(ilinkA->event, ksp, work2, u, NULL);
1345: KSPSolve(ksp, work2, u);
1346: KSPCheckSolve(ksp, pc, u);
1347: PetscLogEventEnd(ilinkA->event, ksp, work2, u, NULL);
1348: MatMult(jac->H, u, Hu); /* alpha = u'*H*u */
1349: VecDot(Hu, u, &alpha);
1350: KSPCheckDot(ksp, alpha);
1352: alpha = PetscSqrtReal(PetscAbsScalar(alpha));
1353: VecScale(u, 1.0 / alpha);
1354: VecAXPBY(d, 1.0 / alpha, 0.0, v); /* v = nu/beta *c */
1356: z = beta / alpha;
1357: vecz[1] = z;
1359: /* Computation of first iterate x(1) and p(1) */
1360: VecAXPY(ilinkA->y, z, u);
1361: VecCopy(d, ilinkD->y);
1362: VecScale(ilinkD->y, -z);
1364: iterGKB = 1;
1365: lowbnd = 2 * jac->gkbtol;
1366: if (jac->gkbmonitor) PetscViewerASCIIPrintf(jac->gkbviewer, "%3" PetscInt_FMT " GKB Lower bound estimate %14.12e\n", iterGKB, (double)lowbnd);
1368: while (iterGKB < jac->gkbmaxit && lowbnd > jac->gkbtol) {
1369: iterGKB += 1;
1370: MatMultHermitianTranspose(jac->B, u, work1); /* v <- nu*(B'*u-alpha/nu*v) */
1371: VecAXPBY(v, nu, -alpha, work1);
1372: VecNorm(v, NORM_2, &beta); /* beta = sqrt(nu)*v'*v */
1373: beta = beta / PetscSqrtReal(nu);
1374: VecScale(v, 1.0 / beta);
1375: MatMult(jac->B, v, work2); /* u <- H^{-1}*(B*v-beta*H*u) */
1376: MatMult(jac->H, u, Hu);
1377: VecAXPY(work2, -beta, Hu);
1378: PetscLogEventBegin(ilinkA->event, ksp, work2, u, NULL);
1379: KSPSolve(ksp, work2, u);
1380: KSPCheckSolve(ksp, pc, u);
1381: PetscLogEventEnd(ilinkA->event, ksp, work2, u, NULL);
1382: MatMult(jac->H, u, Hu); /* alpha = u'*H*u */
1383: VecDot(Hu, u, &alpha);
1384: KSPCheckDot(ksp, alpha);
1386: alpha = PetscSqrtReal(PetscAbsScalar(alpha));
1387: VecScale(u, 1.0 / alpha);
1389: z = -beta / alpha * z; /* z <- beta/alpha*z */
1390: vecz[0] = z;
1392: /* Computation of new iterate x(i+1) and p(i+1) */
1393: VecAXPBY(d, 1.0 / alpha, -beta / alpha, v); /* d = (v-beta*d)/alpha */
1394: VecAXPY(ilinkA->y, z, u); /* r = r + z*u */
1395: VecAXPY(ilinkD->y, -z, d); /* p = p - z*d */
1396: MatMult(jac->H, ilinkA->y, Hu); /* ||u||_H = u'*H*u */
1397: VecDot(Hu, ilinkA->y, &nrmz2);
1399: /* Compute Lower Bound estimate */
1400: if (iterGKB > jac->gkbdelay) {
1401: lowbnd = 0.0;
1402: for (j = 0; j < jac->gkbdelay; j++) lowbnd += PetscAbsScalar(vecz[j] * vecz[j]);
1403: lowbnd = PetscSqrtReal(lowbnd / PetscAbsScalar(nrmz2));
1404: }
1406: for (j = 0; j < jac->gkbdelay - 1; j++) vecz[jac->gkbdelay - j - 1] = vecz[jac->gkbdelay - j - 2];
1407: if (jac->gkbmonitor) PetscViewerASCIIPrintf(jac->gkbviewer, "%3" PetscInt_FMT " GKB Lower bound estimate %14.12e\n", iterGKB, (double)lowbnd);
1408: }
1410: VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE);
1411: VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE);
1412: VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE);
1413: VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE);
1415: return 0;
1416: }
1418: #define FieldSplitSplitSolveAddTranspose(ilink, xx, yy) \
1419: (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) || \
1420: 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) || \
1421: VecScatterEnd(ilink->sctx, ilink->x, yy, ADD_VALUES, SCATTER_REVERSE))
1423: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc, Vec x, Vec y)
1424: {
1425: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1426: PC_FieldSplitLink ilink = jac->head;
1427: PetscInt bs;
1429: if (jac->type == PC_COMPOSITE_ADDITIVE) {
1430: if (jac->defaultsplit) {
1431: VecGetBlockSize(x, &bs);
1433: VecGetBlockSize(y, &bs);
1435: VecStrideGatherAll(x, jac->x, INSERT_VALUES);
1436: while (ilink) {
1437: PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1438: KSPSolveTranspose(ilink->ksp, ilink->x, ilink->y);
1439: KSPCheckSolve(ilink->ksp, pc, ilink->y);
1440: PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1441: ilink = ilink->next;
1442: }
1443: VecStrideScatterAll(jac->y, y, INSERT_VALUES);
1444: } else {
1445: VecSet(y, 0.0);
1446: while (ilink) {
1447: FieldSplitSplitSolveAddTranspose(ilink, x, y);
1448: ilink = ilink->next;
1449: }
1450: }
1451: } else {
1452: if (!jac->w1) {
1453: VecDuplicate(x, &jac->w1);
1454: VecDuplicate(x, &jac->w2);
1455: }
1456: VecSet(y, 0.0);
1457: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1458: FieldSplitSplitSolveAddTranspose(ilink, x, y);
1459: while (ilink->next) {
1460: ilink = ilink->next;
1461: MatMultTranspose(pc->mat, y, jac->w1);
1462: VecWAXPY(jac->w2, -1.0, jac->w1, x);
1463: FieldSplitSplitSolveAddTranspose(ilink, jac->w2, y);
1464: }
1465: while (ilink->previous) {
1466: ilink = ilink->previous;
1467: MatMultTranspose(pc->mat, y, jac->w1);
1468: VecWAXPY(jac->w2, -1.0, jac->w1, x);
1469: FieldSplitSplitSolveAddTranspose(ilink, jac->w2, y);
1470: }
1471: } else {
1472: while (ilink->next) { /* get to last entry in linked list */
1473: ilink = ilink->next;
1474: }
1475: FieldSplitSplitSolveAddTranspose(ilink, x, y);
1476: while (ilink->previous) {
1477: ilink = ilink->previous;
1478: MatMultTranspose(pc->mat, y, jac->w1);
1479: VecWAXPY(jac->w2, -1.0, jac->w1, x);
1480: FieldSplitSplitSolveAddTranspose(ilink, jac->w2, y);
1481: }
1482: }
1483: }
1484: return 0;
1485: }
1487: static PetscErrorCode PCReset_FieldSplit(PC pc)
1488: {
1489: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1490: PC_FieldSplitLink ilink = jac->head, next;
1492: while (ilink) {
1493: KSPDestroy(&ilink->ksp);
1494: VecDestroy(&ilink->x);
1495: VecDestroy(&ilink->y);
1496: VecDestroy(&ilink->z);
1497: VecScatterDestroy(&ilink->sctx);
1498: ISDestroy(&ilink->is);
1499: ISDestroy(&ilink->is_col);
1500: PetscFree(ilink->splitname);
1501: PetscFree(ilink->fields);
1502: PetscFree(ilink->fields_col);
1503: next = ilink->next;
1504: PetscFree(ilink);
1505: ilink = next;
1506: }
1507: jac->head = NULL;
1508: PetscFree2(jac->x, jac->y);
1509: if (jac->mat && jac->mat != jac->pmat) {
1510: MatDestroyMatrices(jac->nsplits, &jac->mat);
1511: } else if (jac->mat) {
1512: jac->mat = NULL;
1513: }
1514: if (jac->pmat) MatDestroyMatrices(jac->nsplits, &jac->pmat);
1515: if (jac->Afield) MatDestroyMatrices(jac->nsplits, &jac->Afield);
1516: jac->nsplits = 0;
1517: VecDestroy(&jac->w1);
1518: VecDestroy(&jac->w2);
1519: MatDestroy(&jac->schur);
1520: MatDestroy(&jac->schurp);
1521: MatDestroy(&jac->schur_user);
1522: KSPDestroy(&jac->kspschur);
1523: KSPDestroy(&jac->kspupper);
1524: MatDestroy(&jac->B);
1525: MatDestroy(&jac->C);
1526: MatDestroy(&jac->H);
1527: VecDestroy(&jac->u);
1528: VecDestroy(&jac->v);
1529: VecDestroy(&jac->Hu);
1530: VecDestroy(&jac->d);
1531: PetscFree(jac->vecz);
1532: PetscViewerDestroy(&jac->gkbviewer);
1533: jac->isrestrict = PETSC_FALSE;
1534: return 0;
1535: }
1537: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
1538: {
1539: PCReset_FieldSplit(pc);
1540: PetscFree(pc->data);
1541: PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL);
1542: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetFields_C", NULL);
1543: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetIS_C", NULL);
1544: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetType_C", NULL);
1545: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetBlockSize_C", NULL);
1546: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitRestrictIS_C", NULL);
1547: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSchurGetSubKSP_C", NULL);
1548: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", NULL);
1550: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", NULL);
1551: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", NULL);
1552: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", NULL);
1553: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", NULL);
1554: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", NULL);
1555: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", NULL);
1556: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", NULL);
1557: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", NULL);
1558: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", NULL);
1559: return 0;
1560: }
1562: static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc, PetscOptionItems *PetscOptionsObject)
1563: {
1564: PetscInt bs;
1565: PetscBool flg;
1566: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1567: PCCompositeType ctype;
1569: PetscOptionsHeadBegin(PetscOptionsObject, "FieldSplit options");
1570: PetscOptionsBool("-pc_fieldsplit_dm_splits", "Whether to use DMCreateFieldDecomposition() for splits", "PCFieldSplitSetDMSplits", jac->dm_splits, &jac->dm_splits, NULL);
1571: PetscOptionsInt("-pc_fieldsplit_block_size", "Blocksize that defines number of fields", "PCFieldSplitSetBlockSize", jac->bs, &bs, &flg);
1572: if (flg) PCFieldSplitSetBlockSize(pc, bs);
1573: jac->diag_use_amat = pc->useAmat;
1574: 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);
1575: jac->offdiag_use_amat = pc->useAmat;
1576: 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);
1577: PetscOptionsBool("-pc_fieldsplit_detect_saddle_point", "Form 2-way split by detecting zero diagonal entries", "PCFieldSplitSetDetectSaddlePoint", jac->detect, &jac->detect, NULL);
1578: PCFieldSplitSetDetectSaddlePoint(pc, jac->detect); /* Sets split type and Schur PC type */
1579: PetscOptionsEnum("-pc_fieldsplit_type", "Type of composition", "PCFieldSplitSetType", PCCompositeTypes, (PetscEnum)jac->type, (PetscEnum *)&ctype, &flg);
1580: if (flg) PCFieldSplitSetType(pc, ctype);
1581: /* Only setup fields once */
1582: if ((jac->bs > 0) && (jac->nsplits == 0)) {
1583: /* only allow user to set fields from command line if bs is already known.
1584: otherwise user can set them in PCFieldSplitSetDefaults() */
1585: PCFieldSplitSetRuntimeSplits_Private(pc);
1586: if (jac->splitdefined) PetscInfo(pc, "Splits defined using the options database\n");
1587: }
1588: if (jac->type == PC_COMPOSITE_SCHUR) {
1589: PetscOptionsGetEnum(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_fieldsplit_schur_factorization_type", PCFieldSplitSchurFactTypes, (PetscEnum *)&jac->schurfactorization, &flg);
1590: if (flg) PetscInfo(pc, "Deprecated use of -pc_fieldsplit_schur_factorization_type\n");
1591: 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);
1592: PetscOptionsEnum("-pc_fieldsplit_schur_precondition", "How to build preconditioner for Schur complement", "PCFieldSplitSetSchurPre", PCFieldSplitSchurPreTypes, (PetscEnum)jac->schurpre, (PetscEnum *)&jac->schurpre, NULL);
1593: PetscOptionsScalar("-pc_fieldsplit_schur_scale", "Scale Schur complement", "PCFieldSplitSetSchurScale", jac->schurscale, &jac->schurscale, NULL);
1594: } else if (jac->type == PC_COMPOSITE_GKB) {
1595: PetscOptionsReal("-pc_fieldsplit_gkb_tol", "The tolerance for the lower bound stopping criterion", "PCFieldSplitGKBTol", jac->gkbtol, &jac->gkbtol, NULL);
1596: PetscOptionsInt("-pc_fieldsplit_gkb_delay", "The delay value for lower bound criterion", "PCFieldSplitGKBDelay", jac->gkbdelay, &jac->gkbdelay, NULL);
1597: PetscOptionsReal("-pc_fieldsplit_gkb_nu", "Parameter in augmented Lagrangian approach", "PCFieldSplitGKBNu", jac->gkbnu, &jac->gkbnu, NULL);
1599: PetscOptionsInt("-pc_fieldsplit_gkb_maxit", "Maximum allowed number of iterations", "PCFieldSplitGKBMaxit", jac->gkbmaxit, &jac->gkbmaxit, NULL);
1600: PetscOptionsBool("-pc_fieldsplit_gkb_monitor", "Prints number of GKB iterations and error", "PCFieldSplitGKB", jac->gkbmonitor, &jac->gkbmonitor, NULL);
1601: }
1602: /*
1603: In the initial call to this routine the sub-solver data structures do not exist so we cannot call KSPSetFromOptions() on them yet.
1604: 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
1605: is called on the outer solver in case changes were made in the options database
1607: 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()
1608: 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.
1609: Without this extra check test p2p1fetidp_olof_full and others fail with incorrect matrix types.
1611: There could be a negative side effect of calling the KSPSetFromOptions() below.
1613: If one captured the PetscObjectState of the options database one could skip these calls if the database has not changed from the previous call
1614: */
1615: if (jac->issetup) {
1616: PC_FieldSplitLink ilink = jac->head;
1617: if (jac->type == PC_COMPOSITE_SCHUR) {
1618: if (jac->kspupper && jac->kspupper->totalits > 0) KSPSetFromOptions(jac->kspupper);
1619: if (jac->kspschur && jac->kspschur->totalits > 0) KSPSetFromOptions(jac->kspschur);
1620: }
1621: while (ilink) {
1622: if (ilink->ksp->totalits > 0) KSPSetFromOptions(ilink->ksp);
1623: ilink = ilink->next;
1624: }
1625: }
1626: PetscOptionsHeadEnd();
1627: return 0;
1628: }
1630: static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc, const char splitname[], PetscInt n, const PetscInt *fields, const PetscInt *fields_col)
1631: {
1632: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1633: PC_FieldSplitLink ilink, next = jac->head;
1634: char prefix[128];
1635: PetscInt i;
1637: if (jac->splitdefined) {
1638: PetscInfo(pc, "Ignoring new split \"%s\" because the splits have already been defined\n", splitname);
1639: return 0;
1640: }
1641: for (i = 0; i < n; i++) {
1644: }
1645: PetscNew(&ilink);
1646: if (splitname) {
1647: PetscStrallocpy(splitname, &ilink->splitname);
1648: } else {
1649: PetscMalloc1(3, &ilink->splitname);
1650: PetscSNPrintf(ilink->splitname, 2, "%" PetscInt_FMT, jac->nsplits);
1651: }
1652: 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 */
1653: PetscMalloc1(n, &ilink->fields);
1654: PetscArraycpy(ilink->fields, fields, n);
1655: PetscMalloc1(n, &ilink->fields_col);
1656: PetscArraycpy(ilink->fields_col, fields_col, n);
1658: ilink->nfields = n;
1659: ilink->next = NULL;
1660: KSPCreate(PetscObjectComm((PetscObject)pc), &ilink->ksp);
1661: KSPSetErrorIfNotConverged(ilink->ksp, pc->erroriffailure);
1662: PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)pc, 1);
1663: KSPSetType(ilink->ksp, KSPPREONLY);
1665: PetscSNPrintf(prefix, sizeof(prefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
1666: KSPSetOptionsPrefix(ilink->ksp, prefix);
1668: if (!next) {
1669: jac->head = ilink;
1670: ilink->previous = NULL;
1671: } else {
1672: while (next->next) next = next->next;
1673: next->next = ilink;
1674: ilink->previous = next;
1675: }
1676: jac->nsplits++;
1677: return 0;
1678: }
1680: static PetscErrorCode PCFieldSplitSchurGetSubKSP_FieldSplit(PC pc, PetscInt *n, KSP **subksp)
1681: {
1682: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1684: *subksp = NULL;
1685: if (n) *n = 0;
1686: if (jac->type == PC_COMPOSITE_SCHUR) {
1687: PetscInt nn;
1691: nn = jac->nsplits + (jac->kspupper != jac->head->ksp ? 1 : 0);
1692: PetscMalloc1(nn, subksp);
1693: (*subksp)[0] = jac->head->ksp;
1694: (*subksp)[1] = jac->kspschur;
1695: if (jac->kspupper != jac->head->ksp) (*subksp)[2] = jac->kspupper;
1696: if (n) *n = nn;
1697: }
1698: return 0;
1699: }
1701: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc, PetscInt *n, KSP **subksp)
1702: {
1703: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1706: PetscMalloc1(jac->nsplits, subksp);
1707: MatSchurComplementGetKSP(jac->schur, *subksp);
1709: (*subksp)[1] = jac->kspschur;
1710: if (n) *n = jac->nsplits;
1711: return 0;
1712: }
1714: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc, PetscInt *n, KSP **subksp)
1715: {
1716: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1717: PetscInt cnt = 0;
1718: PC_FieldSplitLink ilink = jac->head;
1720: PetscMalloc1(jac->nsplits, subksp);
1721: while (ilink) {
1722: (*subksp)[cnt++] = ilink->ksp;
1723: ilink = ilink->next;
1724: }
1726: if (n) *n = jac->nsplits;
1727: return 0;
1728: }
1730: /*@C
1731: PCFieldSplitRestrictIS - Restricts the fieldsplit `IS`s to be within a given `IS`.
1733: Input Parameters:
1734: + pc - the preconditioner context
1735: - is - the index set that defines the indices to which the fieldsplit is to be restricted
1737: Level: advanced
1739: Developer Note:
1740: It seems the resulting `IS`s will not cover the entire space, so
1741: how can they define a convergent preconditioner? Needs explaining.
1743: .seealso: [](sec_block_matrices), `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`
1744: @*/
1745: PetscErrorCode PCFieldSplitRestrictIS(PC pc, IS isy)
1746: {
1749: PetscTryMethod(pc, "PCFieldSplitRestrictIS_C", (PC, IS), (pc, isy));
1750: return 0;
1751: }
1753: static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy)
1754: {
1755: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1756: PC_FieldSplitLink ilink = jac->head, next;
1757: PetscInt localsize, size, sizez, i;
1758: const PetscInt *ind, *indz;
1759: PetscInt *indc, *indcz;
1760: PetscBool flg;
1762: ISGetLocalSize(isy, &localsize);
1763: MPI_Scan(&localsize, &size, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)isy));
1764: size -= localsize;
1765: while (ilink) {
1766: IS isrl, isr;
1767: PC subpc;
1768: ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl);
1769: ISGetLocalSize(isrl, &localsize);
1770: PetscMalloc1(localsize, &indc);
1771: ISGetIndices(isrl, &ind);
1772: PetscArraycpy(indc, ind, localsize);
1773: ISRestoreIndices(isrl, &ind);
1774: ISDestroy(&isrl);
1775: for (i = 0; i < localsize; i++) *(indc + i) += size;
1776: ISCreateGeneral(PetscObjectComm((PetscObject)isy), localsize, indc, PETSC_OWN_POINTER, &isr);
1777: PetscObjectReference((PetscObject)isr);
1778: ISDestroy(&ilink->is);
1779: ilink->is = isr;
1780: PetscObjectReference((PetscObject)isr);
1781: ISDestroy(&ilink->is_col);
1782: ilink->is_col = isr;
1783: ISDestroy(&isr);
1784: KSPGetPC(ilink->ksp, &subpc);
1785: PetscObjectTypeCompare((PetscObject)subpc, PCFIELDSPLIT, &flg);
1786: if (flg) {
1787: IS iszl, isz;
1788: MPI_Comm comm;
1789: ISGetLocalSize(ilink->is, &localsize);
1790: comm = PetscObjectComm((PetscObject)ilink->is);
1791: ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl);
1792: MPI_Scan(&localsize, &sizez, 1, MPIU_INT, MPI_SUM, comm);
1793: sizez -= localsize;
1794: ISGetLocalSize(iszl, &localsize);
1795: PetscMalloc1(localsize, &indcz);
1796: ISGetIndices(iszl, &indz);
1797: PetscArraycpy(indcz, indz, localsize);
1798: ISRestoreIndices(iszl, &indz);
1799: ISDestroy(&iszl);
1800: for (i = 0; i < localsize; i++) *(indcz + i) += sizez;
1801: ISCreateGeneral(comm, localsize, indcz, PETSC_OWN_POINTER, &isz);
1802: PCFieldSplitRestrictIS(subpc, isz);
1803: ISDestroy(&isz);
1804: }
1805: next = ilink->next;
1806: ilink = next;
1807: }
1808: jac->isrestrict = PETSC_TRUE;
1809: return 0;
1810: }
1812: static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc, const char splitname[], IS is)
1813: {
1814: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1815: PC_FieldSplitLink ilink, next = jac->head;
1816: char prefix[128];
1818: if (jac->splitdefined) {
1819: PetscInfo(pc, "Ignoring new split \"%s\" because the splits have already been defined\n", splitname);
1820: return 0;
1821: }
1822: PetscNew(&ilink);
1823: if (splitname) {
1824: PetscStrallocpy(splitname, &ilink->splitname);
1825: } else {
1826: PetscMalloc1(8, &ilink->splitname);
1827: PetscSNPrintf(ilink->splitname, 7, "%" PetscInt_FMT, jac->nsplits);
1828: }
1829: 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 */
1830: PetscObjectReference((PetscObject)is);
1831: ISDestroy(&ilink->is);
1832: ilink->is = is;
1833: PetscObjectReference((PetscObject)is);
1834: ISDestroy(&ilink->is_col);
1835: ilink->is_col = is;
1836: ilink->next = NULL;
1837: KSPCreate(PetscObjectComm((PetscObject)pc), &ilink->ksp);
1838: KSPSetErrorIfNotConverged(ilink->ksp, pc->erroriffailure);
1839: PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)pc, 1);
1840: KSPSetType(ilink->ksp, KSPPREONLY);
1842: PetscSNPrintf(prefix, sizeof(prefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
1843: KSPSetOptionsPrefix(ilink->ksp, prefix);
1845: if (!next) {
1846: jac->head = ilink;
1847: ilink->previous = NULL;
1848: } else {
1849: while (next->next) next = next->next;
1850: next->next = ilink;
1851: ilink->previous = next;
1852: }
1853: jac->nsplits++;
1854: return 0;
1855: }
1857: /*@C
1858: PCFieldSplitSetFields - Sets the fields that define one particular split in `PCFIELDSPLIT`
1860: Logically Collective
1862: Input Parameters:
1863: + pc - the preconditioner context
1864: . splitname - name of this split, if `NULL` the number of the split is used
1865: . n - the number of fields in this split
1866: . fields - the fields in this split
1867: - 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
1868: of the matrix and fields_col provides the column indices for that block
1870: Level: intermediate
1872: Notes:
1873: Use `PCFieldSplitSetIS()` to set a general set of indices as a split.
1875: `PCFieldSplitSetFields()` is for defining fields as strided blocks. For example, if the block
1876: 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
1877: 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x....
1878: where the numbered entries indicate what is in the split.
1880: This function is called once per split (it creates a new split each time). Solve options
1881: for this split will be available under the prefix `-fieldsplit_SPLITNAME_`.
1883: `PCFieldSplitSetIS()` does not support having a fields_col different from fields
1885: Developer Note:
1886: This routine does not actually create the `IS` representing the split, that is delayed
1887: until `PCSetUp_FieldSplit()`, because information about the vector/matrix layouts may not be
1888: available when this routine is called.
1890: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetBlockSize()`, `PCFieldSplitSetIS()`, `PCFieldSplitRestrictIS()`
1891: @*/
1892: PetscErrorCode PCFieldSplitSetFields(PC pc, const char splitname[], PetscInt n, const PetscInt *fields, const PetscInt *fields_col)
1893: {
1898: PetscTryMethod(pc, "PCFieldSplitSetFields_C", (PC, const char[], PetscInt, const PetscInt *, const PetscInt *), (pc, splitname, n, fields, fields_col));
1899: return 0;
1900: }
1902: /*@
1903: PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) to build
1904: the sub-matrices associated with each split. Where `KSPSetOperators`(ksp,Amat,Pmat)) was used to supply the operators.
1906: Logically Collective
1908: Input Parameters:
1909: + pc - the preconditioner object
1910: - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from
1912: Options Database Keys:
1913: . -pc_fieldsplit_diag_use_amat - use the Amat to provide the diagonal blocks
1915: Level: intermediate
1917: .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitGetDiagUseAmat()`, `PCFieldSplitSetOffDiagUseAmat()`, `PCFIELDSPLIT`
1918: @*/
1919: PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc, PetscBool flg)
1920: {
1921: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1922: PetscBool isfs;
1925: PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs);
1927: jac->diag_use_amat = flg;
1928: return 0;
1929: }
1931: /*@
1932: PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) to build
1933: the sub-matrices associated with each split. Where `KSPSetOperators`(ksp,Amat,Pmat)) was used to supply the operators.
1935: Logically Collective
1937: Input Parameters:
1938: . pc - the preconditioner object
1940: Output Parameters:
1941: . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from
1943: Level: intermediate
1945: .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitSetDiagUseAmat()`, `PCFieldSplitGetOffDiagUseAmat()`, `PCFIELDSPLIT`
1946: @*/
1947: PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc, PetscBool *flg)
1948: {
1949: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1950: PetscBool isfs;
1954: PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs);
1956: *flg = jac->diag_use_amat;
1957: return 0;
1958: }
1960: /*@
1961: PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) to build
1962: the sub-matrices associated with each split. Where `KSPSetOperators`(ksp,Amat,Pmat)) was used to supply the operators.
1964: Logically Collective
1966: Input Parameters:
1967: + pc - the preconditioner object
1968: - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from
1970: Options Database Keys:
1971: . -pc_fieldsplit_off_diag_use_amat <bool> - use the Amat to extract the off-diagonal blocks
1973: Level: intermediate
1975: .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitGetOffDiagUseAmat()`, `PCFieldSplitSetDiagUseAmat()`, `PCFIELDSPLIT`
1976: @*/
1977: PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc, PetscBool flg)
1978: {
1979: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1980: PetscBool isfs;
1983: PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs);
1985: jac->offdiag_use_amat = flg;
1986: return 0;
1987: }
1989: /*@
1990: PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) to build
1991: the sub-matrices associated with each split. Where `KSPSetOperators`(ksp,Amat,Pmat)) was used to supply the operators.
1993: Logically Collective
1995: Input Parameters:
1996: . pc - the preconditioner object
1998: Output Parameters:
1999: . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from
2001: Level: intermediate
2003: .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitSetOffDiagUseAmat()`, `PCFieldSplitGetDiagUseAmat()`, `PCFIELDSPLIT`
2004: @*/
2005: PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc, PetscBool *flg)
2006: {
2007: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2008: PetscBool isfs;
2012: PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs);
2014: *flg = jac->offdiag_use_amat;
2015: return 0;
2016: }
2018: /*@C
2019: PCFieldSplitSetIS - Sets the exact elements for a split in a `PCFIELDSPLIT`
2021: Logically Collective
2023: Input Parameters:
2024: + pc - the preconditioner context
2025: . splitname - name of this split, if `NULL` the number of the split is used
2026: - is - the index set that defines the elements in this split
2028: Level: intermediate
2030: Notes:
2031: Use `PCFieldSplitSetFields()`, for splits defined by strided types.
2033: This function is called once per split (it creates a new split each time). Solve options
2034: for this split will be available under the prefix -fieldsplit_SPLITNAME_.
2036: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetBlockSize()`
2037: @*/
2038: PetscErrorCode PCFieldSplitSetIS(PC pc, const char splitname[], IS is)
2039: {
2043: PetscTryMethod(pc, "PCFieldSplitSetIS_C", (PC, const char[], IS), (pc, splitname, is));
2044: return 0;
2045: }
2047: /*@C
2048: PCFieldSplitGetIS - Retrieves the elements for a split as an `IS`
2050: Logically Collective
2052: Input Parameters:
2053: + pc - the preconditioner context
2054: - splitname - name of this split
2056: Output Parameter:
2057: - is - the index set that defines the elements in this split, or `NULL` if the split is not found
2059: Level: intermediate
2061: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetIS()`
2062: @*/
2063: PetscErrorCode PCFieldSplitGetIS(PC pc, const char splitname[], IS *is)
2064: {
2068: {
2069: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2070: PC_FieldSplitLink ilink = jac->head;
2071: PetscBool found;
2073: *is = NULL;
2074: while (ilink) {
2075: PetscStrcmp(ilink->splitname, splitname, &found);
2076: if (found) {
2077: *is = ilink->is;
2078: break;
2079: }
2080: ilink = ilink->next;
2081: }
2082: }
2083: return 0;
2084: }
2086: /*@C
2087: PCFieldSplitGetISByIndex - Retrieves the elements for a given split as an `IS`
2089: Logically Collective
2091: Input Parameters:
2092: + pc - the preconditioner context
2093: - index - index of this split
2095: Output Parameter:
2096: - is - the index set that defines the elements in this split
2098: Level: intermediate
2100: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitGetIS()`, `PCFieldSplitSetIS()`
2101: @*/
2102: PetscErrorCode PCFieldSplitGetISByIndex(PC pc, PetscInt index, IS *is)
2103: {
2107: {
2108: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2109: PC_FieldSplitLink ilink = jac->head;
2110: PetscInt i = 0;
2113: while (i < index) {
2114: ilink = ilink->next;
2115: ++i;
2116: }
2117: PCFieldSplitGetIS(pc, ilink->splitname, is);
2118: }
2119: return 0;
2120: }
2122: /*@
2123: PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the
2124: fieldsplit preconditioner when calling `PCFieldSplitSetIS()`. If not set the matrix block size is used.
2126: Logically Collective
2128: Input Parameters:
2129: + pc - the preconditioner context
2130: - bs - the block size
2132: Level: intermediate
2134: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`
2135: @*/
2136: PetscErrorCode PCFieldSplitSetBlockSize(PC pc, PetscInt bs)
2137: {
2140: PetscTryMethod(pc, "PCFieldSplitSetBlockSize_C", (PC, PetscInt), (pc, bs));
2141: return 0;
2142: }
2144: /*@C
2145: PCFieldSplitGetSubKSP - Gets the `KSP` contexts for all splits
2147: Collective
2149: Input Parameter:
2150: . pc - the preconditioner context
2152: Output Parameters:
2153: + n - the number of splits
2154: - subksp - the array of `KSP` contexts
2156: Level: advanced
2158: Notes:
2159: After `PCFieldSplitGetSubKSP()` the array of `KSP`s is to be freed by the user with `PetscFree()`
2160: (not the `KSP`, just the array that contains them).
2162: You must call `PCSetUp()` before calling `PCFieldSplitGetSubKSP()`.
2164: If the fieldsplit is of type `PC_COMPOSITE_SCHUR`, it returns the `KSP` object used inside the
2165: Schur complement and the `KSP` object used to iterate over the Schur complement.
2166: To access all the `KSP` objects used in `PC_COMPOSITE_SCHUR`, use `PCFieldSplitSchurGetSubKSP()`.
2168: If the fieldsplit is of type `PC_COMPOSITE_GKB`, it returns the `KSP` object used to solve the
2169: inner linear system defined by the matrix H in each loop.
2171: Fortran Usage:
2172: You must pass in a `KSP` array that is large enough to contain all the `KSP`s.
2173: You can call `PCFieldSplitGetSubKSP`(pc,n,`PETSC_NULL_KSP`,ierr) to determine how large the
2174: `KSP` array must be.
2176: Developer Note:
2177: There should be a `PCFieldSplitRestoreSubKSP()` instead of requiring the user to call `PetscFree()`
2179: The Fortran interface should be modernized to return directly the array of values.
2181: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`, `PCFieldSplitSchurGetSubKSP()`
2182: @*/
2183: PetscErrorCode PCFieldSplitGetSubKSP(PC pc, PetscInt *n, KSP *subksp[])
2184: {
2187: PetscUseMethod(pc, "PCFieldSplitGetSubKSP_C", (PC, PetscInt *, KSP **), (pc, n, subksp));
2188: return 0;
2189: }
2191: /*@C
2192: PCFieldSplitSchurGetSubKSP - Gets the `KSP` contexts used inside the Schur complement based `PCFIELDSPLIT`
2194: Collective
2196: Input Parameter:
2197: . pc - the preconditioner context
2199: Output Parameters:
2200: + n - the number of splits
2201: - subksp - the array of `KSP` contexts
2203: Level: advanced
2205: Notes:
2206: After `PCFieldSplitSchurGetSubKSP()` the array of `KSP`s is to be freed by the user with `PetscFree()`
2207: (not the `KSP` just the array that contains them).
2209: You must call `PCSetUp()` before calling `PCFieldSplitSchurGetSubKSP()`.
2211: If the fieldsplit type is of type `PC_COMPOSITE_SCHUR`, it returns (in order)
2212: + 1 - the `KSP` used for the (1,1) block
2213: . 2 - the `KSP` used for the Schur complement (not the one used for the interior Schur solver)
2214: - 3 - the `KSP` used for the (1,1) block in the upper triangular factor (if different from that of the (1,1) block).
2216: It returns a null array if the fieldsplit is not of type `PC_COMPOSITE_SCHUR`; in this case, you should use `PCFieldSplitGetSubKSP()`.
2218: Fortran Note:
2219: You must pass in a `KSP` array that is large enough to contain all the local `KSP`s.
2220: You can call `PCFieldSplitSchurGetSubKSP`(pc,n,`PETSC_NULL_KSP`,ierr) to determine how large the
2221: `KSP` array must be.
2223: Developer Notes:
2224: There should be a `PCFieldSplitRestoreSubKSP()` instead of requiring the user to call `PetscFree()`
2226: Should the functionality of `PCFieldSplitSchurGetSubKSP()` and `PCFieldSplitGetSubKSP()` be merged?
2228: The Fortran interface should be modernized to return directly the array of values.
2230: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`, `PCFieldSplitGetSubKSP()`
2231: @*/
2232: PetscErrorCode PCFieldSplitSchurGetSubKSP(PC pc, PetscInt *n, KSP *subksp[])
2233: {
2236: PetscUseMethod(pc, "PCFieldSplitSchurGetSubKSP_C", (PC, PetscInt *, KSP **), (pc, n, subksp));
2237: return 0;
2238: }
2240: /*@
2241: PCFieldSplitSetSchurPre - Indicates from what operator the preconditioner is constructucted for the Schur complement.
2242: The default is the A11 matrix.
2244: Collective
2246: Input Parameters:
2247: + pc - the preconditioner context
2248: . ptype - which matrix to use for preconditioning the Schur complement: `PC_FIELDSPLIT_SCHUR_PRE_A11` (default),
2249: `PC_FIELDSPLIT_SCHUR_PRE_SELF`, `PC_FIELDSPLIT_SCHUR_PRE_USER`,
2250: `PC_FIELDSPLIT_SCHUR_PRE_SELFP`, and `PC_FIELDSPLIT_SCHUR_PRE_FULL`
2251: - pre - matrix to use for preconditioning, or `NULL`
2253: Options Database Keys:
2254: + -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11. See notes for meaning of various arguments
2255: - -fieldsplit_1_pc_type <pctype> - the preconditioner algorithm that is used to construct the preconditioner from the operator
2257: Level: intermediate
2259: Notes:
2260: If ptype is
2261: + a11 - the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner
2262: matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix
2263: . self - the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix:
2264: The only preconditioner that currently works with this symbolic representation matrix object is the `PCLSC`
2265: preconditioner
2266: . user - the preconditioner for the Schur complement is generated from the user provided matrix (pre argument
2267: to this function).
2268: . selfp - the preconditioning for the Schur complement is generated from an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01
2269: This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be
2270: lumped before extracting the diagonal using the additional option -fieldsplit_1_mat_schur_complement_ainv_type lump
2271: - full - the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation
2272: computed internally by `PCFIELDSPLIT` (this is expensive)
2273: useful mostly as a test that the Schur complement approach can work for your problem
2275: When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense
2276: with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and
2277: -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement.
2279: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSchurPre()`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`,
2280: `MatSchurComplementSetAinvType()`, `PCLSC`,
2281: `PCFieldSplitSchurPreType`
2282: @*/
2283: PetscErrorCode PCFieldSplitSetSchurPre(PC pc, PCFieldSplitSchurPreType ptype, Mat pre)
2284: {
2286: PetscTryMethod(pc, "PCFieldSplitSetSchurPre_C", (PC, PCFieldSplitSchurPreType, Mat), (pc, ptype, pre));
2287: return 0;
2288: }
2290: PetscErrorCode PCFieldSplitSchurPrecondition(PC pc, PCFieldSplitSchurPreType ptype, Mat pre)
2291: {
2292: return PCFieldSplitSetSchurPre(pc, ptype, pre);
2293: } /* Deprecated name */
2295: /*@
2296: PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be
2297: preconditioned. See `PCFieldSplitSetSchurPre()` for details.
2299: Logically Collective
2301: Input Parameter:
2302: . pc - the preconditioner context
2304: Output Parameters:
2305: + ptype - which matrix to use for preconditioning the Schur complement: `PC_FIELDSPLIT_SCHUR_PRE_A11`, `PC_FIELDSPLIT_SCHUR_PRE_SELF`, `PC_FIELDSPLIT_PRE_USER`
2306: - pre - matrix to use for preconditioning (with `PC_FIELDSPLIT_PRE_USER`), or NULL
2308: Level: intermediate
2310: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitSetSchurPre()`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, `PCLSC`,
2311: `PCFieldSplitSchurPreType`
2312: @*/
2313: PetscErrorCode PCFieldSplitGetSchurPre(PC pc, PCFieldSplitSchurPreType *ptype, Mat *pre)
2314: {
2316: PetscUseMethod(pc, "PCFieldSplitGetSchurPre_C", (PC, PCFieldSplitSchurPreType *, Mat *), (pc, ptype, pre));
2317: return 0;
2318: }
2320: /*@
2321: PCFieldSplitSchurGetS - extract the `MATSCHURCOMPLEMENT` object used by this `PCFIELDSPLIT` in case it needs to be configured separately
2323: Not collective
2325: Input Parameter:
2326: . pc - the preconditioner context
2328: Output Parameter:
2329: . S - the Schur complement matrix
2331: Level: advanced
2333: Note:
2334: This matrix should not be destroyed using `MatDestroy()`; rather, use `PCFieldSplitSchurRestoreS()`.
2336: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurPre()`, `MATSCHURCOMPLEMENT`, `PCFieldSplitSchurRestoreS()`,
2337: `MatCreateSchurComplement()`, `MatSchurComplementGetKSP()`, `MatSchurComplementComputeExplicitOperator()`, `MatGetSchurComplement()`
2338: @*/
2339: PetscErrorCode PCFieldSplitSchurGetS(PC pc, Mat *S)
2340: {
2341: const char *t;
2342: PetscBool isfs;
2343: PC_FieldSplit *jac;
2346: PetscObjectGetType((PetscObject)pc, &t);
2347: PetscStrcmp(t, PCFIELDSPLIT, &isfs);
2349: jac = (PC_FieldSplit *)pc->data;
2351: if (S) *S = jac->schur;
2352: return 0;
2353: }
2355: /*@
2356: PCFieldSplitSchurRestoreS - returns the `MATSCHURCOMPLEMENT` matrix used by this `PC`
2358: Not collective
2360: Input Parameters:
2361: + pc - the preconditioner context
2362: - S - the Schur complement matrix
2364: Level: advanced
2366: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurPre()`, `MatSchurComplement`, `PCFieldSplitSchurGetS()`
2367: @*/
2368: PetscErrorCode PCFieldSplitSchurRestoreS(PC pc, Mat *S)
2369: {
2370: const char *t;
2371: PetscBool isfs;
2372: PC_FieldSplit *jac;
2375: PetscObjectGetType((PetscObject)pc, &t);
2376: PetscStrcmp(t, PCFIELDSPLIT, &isfs);
2378: jac = (PC_FieldSplit *)pc->data;
2381: return 0;
2382: }
2384: static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc, PCFieldSplitSchurPreType ptype, Mat pre)
2385: {
2386: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2388: jac->schurpre = ptype;
2389: if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
2390: MatDestroy(&jac->schur_user);
2391: jac->schur_user = pre;
2392: PetscObjectReference((PetscObject)jac->schur_user);
2393: }
2394: return 0;
2395: }
2397: static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc, PCFieldSplitSchurPreType *ptype, Mat *pre)
2398: {
2399: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2401: *ptype = jac->schurpre;
2402: *pre = jac->schur_user;
2403: return 0;
2404: }
2406: /*@
2407: PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain in the preconditioner
2409: Collective
2411: Input Parameters:
2412: + pc - the preconditioner context
2413: - ftype - which blocks of factorization to retain, `PC_FIELDSPLIT_SCHUR_FACT_FULL` is default
2415: Options Database Key:
2416: . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> - default is full
2418: Level: intermediate
2420: Notes:
2421: The FULL factorization is
2423: .vb
2424: (A B) = (1 0) (A 0) (1 Ainv*B) = L D U
2425: (C E) (C*Ainv 1) (0 S) (0 1)
2426: .vb
2427: where S = E - C*Ainv*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,
2428: and DIAG is the diagonal part with the sign of S flipped (because this makes the preconditioner positive definite for many formulations, thus allowing the use of KSPMINRES). Sign flipping of S can be turned off with PCFieldSplitSetSchurScale().
2430: If A and S are solved exactly
2431: .vb
2432: *) FULL factorization is a direct solver.
2433: *) 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.
2434: *) 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.
2435: .ve
2437: If the iteration count is very low, consider using `KSPFGMRES` or `KSPGCR` which can use one less preconditioner
2438: application in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice.
2440: For symmetric problems in which A is positive definite and S is negative definite, DIAG can be used with `KSPMINRES`.
2442: A flexible method like `KSPFGMRES` or `KSPGCR` 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).
2444: References:
2445: + * - Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000).
2446: - * - Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001).
2448: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurScale()`
2449: @*/
2450: PetscErrorCode PCFieldSplitSetSchurFactType(PC pc, PCFieldSplitSchurFactType ftype)
2451: {
2453: PetscTryMethod(pc, "PCFieldSplitSetSchurFactType_C", (PC, PCFieldSplitSchurFactType), (pc, ftype));
2454: return 0;
2455: }
2457: static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc, PCFieldSplitSchurFactType ftype)
2458: {
2459: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2461: jac->schurfactorization = ftype;
2462: return 0;
2463: }
2465: /*@
2466: PCFieldSplitSetSchurScale - Controls the sign flip of S for `PC_FIELDSPLIT_SCHUR_FACT_DIAG`.
2468: Collective
2470: Input Parameters:
2471: + pc - the preconditioner context
2472: - scale - scaling factor for the Schur complement
2474: Options Database Key:
2475: . -pc_fieldsplit_schur_scale - default is -1.0
2477: Level: intermediate
2479: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurFactType`, `PCFieldSplitSetSchurScale()`, `PCFieldSplitSetSchurFactType()`
2480: @*/
2481: PetscErrorCode PCFieldSplitSetSchurScale(PC pc, PetscScalar scale)
2482: {
2485: PetscTryMethod(pc, "PCFieldSplitSetSchurScale_C", (PC, PetscScalar), (pc, scale));
2486: return 0;
2487: }
2489: static PetscErrorCode PCFieldSplitSetSchurScale_FieldSplit(PC pc, PetscScalar scale)
2490: {
2491: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2493: jac->schurscale = scale;
2494: return 0;
2495: }
2497: /*@C
2498: PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement
2500: Collective
2502: Input Parameter:
2503: . pc - the preconditioner context
2505: Output Parameters:
2506: + A00 - the (0,0) block
2507: . A01 - the (0,1) block
2508: . A10 - the (1,0) block
2509: - A11 - the (1,1) block
2511: Level: advanced
2513: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `MatSchurComplementGetSubMatrices()`, `MatSchurComplementSetSubMatrices()`
2514: @*/
2515: PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc, Mat *A00, Mat *A01, Mat *A10, Mat *A11)
2516: {
2517: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2521: if (A00) *A00 = jac->pmat[0];
2522: if (A01) *A01 = jac->B;
2523: if (A10) *A10 = jac->C;
2524: if (A11) *A11 = jac->pmat[1];
2525: return 0;
2526: }
2528: /*@
2529: PCFieldSplitSetGKBTol - Sets the solver tolerance for the generalized Golub-Kahan bidiagonalization preconditioner in `PCFIELDSPLIT`
2531: Collective
2533: Input Parameters:
2534: + pc - the preconditioner context
2535: - tolerance - the solver tolerance
2537: Options Database Key:
2538: . -pc_fieldsplit_gkb_tol - default is 1e-5
2540: Level: intermediate
2542: Note:
2543: The generalized GKB algorithm uses a lower bound estimate of the error in energy norm as stopping criterion.
2544: It stops once the lower bound estimate undershoots the required solver tolerance. Although the actual error might be bigger than
2545: this estimate, the stopping criterion is satisfactory in practical cases [A13].
2547: References:
2548: [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013.
2550: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBNu()`, `PCFieldSplitSetGKBMaxit()`
2551: @*/
2552: PetscErrorCode PCFieldSplitSetGKBTol(PC pc, PetscReal tolerance)
2553: {
2556: PetscTryMethod(pc, "PCFieldSplitSetGKBTol_C", (PC, PetscReal), (pc, tolerance));
2557: return 0;
2558: }
2560: static PetscErrorCode PCFieldSplitSetGKBTol_FieldSplit(PC pc, PetscReal tolerance)
2561: {
2562: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2564: jac->gkbtol = tolerance;
2565: return 0;
2566: }
2568: /*@
2569: PCFieldSplitSetGKBMaxit - Sets the maximum number of iterations for the generalized Golub-Kahan bidiagonalization preconditioner in `PCFIELDSPLIT`
2571: Collective
2573: Input Parameters:
2574: + pc - the preconditioner context
2575: - maxit - the maximum number of iterations
2577: Options Database Key:
2578: . -pc_fieldsplit_gkb_maxit - default is 100
2580: Level: intermediate
2582: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBNu()`
2583: @*/
2584: PetscErrorCode PCFieldSplitSetGKBMaxit(PC pc, PetscInt maxit)
2585: {
2588: PetscTryMethod(pc, "PCFieldSplitSetGKBMaxit_C", (PC, PetscInt), (pc, maxit));
2589: return 0;
2590: }
2592: static PetscErrorCode PCFieldSplitSetGKBMaxit_FieldSplit(PC pc, PetscInt maxit)
2593: {
2594: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2596: jac->gkbmaxit = maxit;
2597: return 0;
2598: }
2600: /*@
2601: PCFieldSplitSetGKBDelay - Sets the delay in the lower bound error estimate in the generalized Golub-Kahan bidiagonalization in `PCFIELDSPLIT`
2602: preconditioner.
2604: Collective
2606: Input Parameters:
2607: + pc - the preconditioner context
2608: - delay - the delay window in the lower bound estimate
2610: Options Database Key:
2611: . -pc_fieldsplit_gkb_delay - default is 5
2613: Level: intermediate
2615: Note:
2616: 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
2617: is expressed as a truncated sum. The error at iteration k can only be measured at iteration (k + delay), and thus the algorithm needs
2618: at least (delay + 1) iterations to stop. For more details on the generalized Golub-Kahan bidiagonalization method and its lower bound stopping criterion, please refer to
2620: References:
2621: [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013.
2623: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBNu()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBMaxit()`
2624: @*/
2625: PetscErrorCode PCFieldSplitSetGKBDelay(PC pc, PetscInt delay)
2626: {
2629: PetscTryMethod(pc, "PCFieldSplitSetGKBDelay_C", (PC, PetscInt), (pc, delay));
2630: return 0;
2631: }
2633: static PetscErrorCode PCFieldSplitSetGKBDelay_FieldSplit(PC pc, PetscInt delay)
2634: {
2635: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2637: jac->gkbdelay = delay;
2638: return 0;
2639: }
2641: /*@
2642: PCFieldSplitSetGKBNu - Sets the scalar value nu >= 0 in the transformation H = A00 + nu*A01*A01' of the (1,1) block in the Golub-Kahan bidiagonalization preconditioner
2643: in `PCFIELDSPLIT`
2645: Collective
2647: Input Parameters:
2648: + pc - the preconditioner context
2649: - nu - the shift parameter
2651: Options Database Keys:
2652: . -pc_fieldsplit_gkb_nu - default is 1
2654: Level: intermediate
2656: Notes:
2657: 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 big. However,
2658: if nu is chosen too big, 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
2659: necessary to find a good balance in between the convergence of the inner and outer loop.
2661: For nu = 0, no shift is done. In this case A00 has to be positive definite. The matrix N in [Ar13] is then chosen as identity.
2663: References:
2664: [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013.
2666: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBMaxit()`
2667: @*/
2668: PetscErrorCode PCFieldSplitSetGKBNu(PC pc, PetscReal nu)
2669: {
2672: PetscTryMethod(pc, "PCFieldSplitSetGKBNu_C", (PC, PetscReal), (pc, nu));
2673: return 0;
2674: }
2676: static PetscErrorCode PCFieldSplitSetGKBNu_FieldSplit(PC pc, PetscReal nu)
2677: {
2678: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2680: jac->gkbnu = nu;
2681: return 0;
2682: }
2684: static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc, PCCompositeType type)
2685: {
2686: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2688: jac->type = type;
2689: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", NULL);
2690: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", NULL);
2691: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", NULL);
2692: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", NULL);
2693: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", NULL);
2694: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", NULL);
2695: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", NULL);
2696: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", NULL);
2697: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", NULL);
2699: if (type == PC_COMPOSITE_SCHUR) {
2700: pc->ops->apply = PCApply_FieldSplit_Schur;
2701: pc->ops->view = PCView_FieldSplit_Schur;
2703: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit_Schur);
2704: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", PCFieldSplitSetSchurPre_FieldSplit);
2705: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", PCFieldSplitGetSchurPre_FieldSplit);
2706: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", PCFieldSplitSetSchurFactType_FieldSplit);
2707: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", PCFieldSplitSetSchurScale_FieldSplit);
2708: } else if (type == PC_COMPOSITE_GKB) {
2709: pc->ops->apply = PCApply_FieldSplit_GKB;
2710: pc->ops->view = PCView_FieldSplit_GKB;
2712: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit);
2713: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", PCFieldSplitSetGKBTol_FieldSplit);
2714: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", PCFieldSplitSetGKBMaxit_FieldSplit);
2715: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", PCFieldSplitSetGKBNu_FieldSplit);
2716: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", PCFieldSplitSetGKBDelay_FieldSplit);
2717: } else {
2718: pc->ops->apply = PCApply_FieldSplit;
2719: pc->ops->view = PCView_FieldSplit;
2721: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit);
2722: }
2723: return 0;
2724: }
2726: static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc, PetscInt bs)
2727: {
2728: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2732: jac->bs = bs;
2733: return 0;
2734: }
2736: static PetscErrorCode PCSetCoordinates_FieldSplit(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[])
2737: {
2738: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2739: PC_FieldSplitLink ilink_current = jac->head;
2740: IS is_owned;
2742: jac->coordinates_set = PETSC_TRUE; // Internal flag
2743: MatGetOwnershipIS(pc->mat, &is_owned, PETSC_NULL);
2745: while (ilink_current) {
2746: // For each IS, embed it to get local coords indces
2747: IS is_coords;
2748: PetscInt ndofs_block;
2749: const PetscInt *block_dofs_enumeration; // Numbering of the dofs relevant to the current block
2751: // Setting drop to true for safety. It should make no difference.
2752: ISEmbed(ilink_current->is, is_owned, PETSC_TRUE, &is_coords);
2753: ISGetLocalSize(is_coords, &ndofs_block);
2754: ISGetIndices(is_coords, &block_dofs_enumeration);
2756: // Allocate coordinates vector and set it directly
2757: PetscMalloc1(ndofs_block * dim, &(ilink_current->coords));
2758: for (PetscInt dof = 0; dof < ndofs_block; ++dof) {
2759: for (PetscInt d = 0; d < dim; ++d) (ilink_current->coords)[dim * dof + d] = coords[dim * block_dofs_enumeration[dof] + d];
2760: }
2761: ilink_current->dim = dim;
2762: ilink_current->ndofs = ndofs_block;
2763: ISRestoreIndices(is_coords, &block_dofs_enumeration);
2764: ISDestroy(&is_coords);
2765: ilink_current = ilink_current->next;
2766: }
2767: ISDestroy(&is_owned);
2768: return 0;
2769: }
2771: /*@
2772: PCFieldSplitSetType - Sets the type, `PCCompositeType`, of a `PCFIELDSPLIT`
2774: Collective
2776: Input Parameters:
2777: + pc - the preconditioner context
2778: - type - `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE` (default), `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`
2780: Options Database Key:
2781: . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type
2783: Level: Intermediate
2785: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCCompositeType`, `PCCompositeGetType()`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`,
2786: `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`
2787: @*/
2788: PetscErrorCode PCFieldSplitSetType(PC pc, PCCompositeType type)
2789: {
2791: PetscTryMethod(pc, "PCFieldSplitSetType_C", (PC, PCCompositeType), (pc, type));
2792: return 0;
2793: }
2795: /*@
2796: PCFieldSplitGetType - Gets the type, `PCCompositeType`, of a `PCFIELDSPLIT`
2798: Not collective
2800: Input Parameter:
2801: . pc - the preconditioner context
2803: Output Parameter:
2804: . type - `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE` (default), `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`
2806: Level: Intermediate
2808: .seealso: [](sec_block_matrices), `PC`, `PCCompositeSetType()`, `PCFIELDSPLIT`, `PCCompositeType`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`,
2809: `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`
2810: @*/
2811: PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
2812: {
2813: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2817: *type = jac->type;
2818: return 0;
2819: }
2821: /*@
2822: PCFieldSplitSetDMSplits - Flags whether `DMCreateFieldDecomposition()` should be used to define the splits in a `PCFIELDSPLIT`, whenever possible.
2824: Logically Collective
2826: Input Parameters:
2827: + pc - the preconditioner context
2828: - flg - boolean indicating whether to use field splits defined by the `DM`
2830: Options Database Key:
2831: . -pc_fieldsplit_dm_splits <bool> - use the field splits defined by the `DM`
2833: Level: Intermediate
2835: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitGetDMSplits()`, `PCFieldSplitSetFields()`, `PCFieldsplitSetIS()`
2836: @*/
2837: PetscErrorCode PCFieldSplitSetDMSplits(PC pc, PetscBool flg)
2838: {
2839: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2840: PetscBool isfs;
2844: PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs);
2845: if (isfs) jac->dm_splits = flg;
2846: return 0;
2847: }
2849: /*@
2850: PCFieldSplitGetDMSplits - Returns flag indicating whether `DMCreateFieldDecomposition()` should be used to define the splits in a `PCFIELDSPLIT`, whenever possible.
2852: Logically Collective
2854: Input Parameter:
2855: . pc - the preconditioner context
2857: Output Parameter:
2858: . flg - boolean indicating whether to use field splits defined by the `DM`
2860: Level: Intermediate
2862: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetDMSplits()`, `PCFieldSplitSetFields()`, `PCFieldsplitSetIS()`
2863: @*/
2864: PetscErrorCode PCFieldSplitGetDMSplits(PC pc, PetscBool *flg)
2865: {
2866: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2867: PetscBool isfs;
2871: PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs);
2872: if (isfs) {
2873: if (flg) *flg = jac->dm_splits;
2874: }
2875: return 0;
2876: }
2878: /*@
2879: PCFieldSplitGetDetectSaddlePoint - Returns flag indicating whether `PCFIELDSPLIT` will attempt to automatically determine fields based on zero diagonal entries.
2881: Logically Collective
2883: Input Parameter:
2884: . pc - the preconditioner context
2886: Output Parameter:
2887: . flg - boolean indicating whether to detect fields or not
2889: Level: Intermediate
2891: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetDetectSaddlePoint()`
2892: @*/
2893: PetscErrorCode PCFieldSplitGetDetectSaddlePoint(PC pc, PetscBool *flg)
2894: {
2895: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2897: *flg = jac->detect;
2898: return 0;
2899: }
2901: /*@
2902: PCFieldSplitSetDetectSaddlePoint - Sets flag indicating whether `PCFIELDSPLIT` will attempt to automatically determine fields based on zero diagonal entries.
2904: Logically Collective
2906: Input Parameter:
2907: . pc - the preconditioner context
2909: Output Parameter:
2910: . flg - boolean indicating whether to detect fields or not
2912: Options Database Key:
2913: . -pc_fieldsplit_detect_saddle_point <bool> - detect and use the saddle point
2915: Level: Intermediate
2917: Note:
2918: Also sets the split type to `PC_COMPOSITE_SCHUR` (see `PCFieldSplitSetType()`) and the Schur preconditioner type to `PC_FIELDSPLIT_SCHUR_PRE_SELF` (see `PCFieldSplitSetSchurPre()`).
2920: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitGetDetectSaddlePoint()`, `PCFieldSplitSetType()`, `PCFieldSplitSetSchurPre()`, `PC_FIELDSPLIT_SCHUR_PRE_SELF`
2921: @*/
2922: PetscErrorCode PCFieldSplitSetDetectSaddlePoint(PC pc, PetscBool flg)
2923: {
2924: PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2926: jac->detect = flg;
2927: if (jac->detect) {
2928: PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
2929: PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_SELF, NULL);
2930: }
2931: return 0;
2932: }
2934: /*MC
2935: PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
2936: collections of variables (that may overlap) called splits. See [the users manual section on "Solving Block Matrices"](sec_block_matrices) for more details.
2938: Options Database Keys:
2939: + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the `%d`'th split
2940: . -pc_fieldsplit_default - automatically add any fields to additional splits that have not
2941: been supplied explicitly by `-pc_fieldsplit_%d_fields`
2942: . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields)
2943: . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur,gkb> - type of relaxation or factorization splitting
2944: . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11; see `PCFieldSplitSetSchurPre()`
2945: . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> - set factorization type when using `-pc_fieldsplit_type schur`; see `PCFieldSplitSetSchurFactType()`
2946: - -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero diagonal and uses Schur complement with no preconditioner as the solver
2948: Options prefixes for inner solvers when using the Schur complement preconditioner are `-fieldsplit_0_` and `-fieldsplit_1_` .
2949: The options prefix for the inner solver when using the Golub-Kahan biadiagonalization preconditioner is `-fieldsplit_0_`
2950: For all other solvers they are `-fieldsplit_%d_` for the `%d`'th field; use `-fieldsplit_` for all fields.
2952: To set options on the solvers for each block append `-fieldsplit_` to all the `PC`
2953: options database keys. For example, `-fieldsplit_pc_type ilu` `-fieldsplit_pc_factor_levels 1`
2955: To set the options on the solvers separate for each block call `PCFieldSplitGetSubKSP()`
2956: and set the options directly on the resulting `KSP` object
2958: Level: intermediate
2960: Notes:
2961: Use `PCFieldSplitSetFields()` to set splits defined by "strided" entries and `PCFieldSplitSetIS()`
2962: to define a split by an arbitrary collection of entries.
2964: If no splits are set the default is used. The splits are defined by entries strided by bs,
2965: beginning at 0 then 1, etc to bs-1. The block size can be set with `PCFieldSplitSetBlockSize()`,
2966: if this is not called the block size defaults to the blocksize of the second matrix passed
2967: to `KSPSetOperators()`/`PCSetOperators()`.
2969: For the Schur complement preconditioner if
2971: ```{math}
2972: J = \left[\begin{array}{cc} A_{00} & A_{01} \\ A_{10} & A_{11} \end{array}\right]
2973: ```
2975: the preconditioner using `full` factorization is logically
2976: ```{math}
2977: \left[\begin{array}{cc} I & -\text{ksp}(A_{00}) \\ 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]
2978: ```
2979: where the action of $\text{inv}(A_{00})$ is applied using the KSP solver with prefix `-fieldsplit_0_`. $S$ is the Schur complement
2980: ```{math}
2981: S = A_{11} - A_{10} \text{ksp}(A_{00}) A_{01}
2982: ```
2983: 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
2984: in providing the SECOND split or 1 if not given). For `PCFieldSplitGetSubKSP()` when field number is 0,
2985: it returns the KSP associated with `-fieldsplit_0_` while field number 1 gives `-fieldsplit_1_` KSP. By default
2986: $A_{11}$ is used to construct a preconditioner for $S$, use `PCFieldSplitSetSchurPre()` for all the possible ways to construct the preconditioner for $S$.
2988: The factorization type is set using `-pc_fieldsplit_schur_fact_type <diag, lower, upper, full>`. `full` is shown above,
2989: `diag` gives
2990: ```{math}
2991: \left[\begin{array}{cc} \text{inv}(A_{00}) & 0 \\ 0 & -\text{ksp}(S) \end{array}\right]
2992: ```
2993: 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
2994: can be turned off with `PCFieldSplitSetSchurScale()` or by command line `-pc_fieldsplit_schur_scale 1.0`. The `lower` factorization is the inverse of
2995: ```{math}
2996: \left[\begin{array}{cc} A_{00} & 0 \\ A_{10} & S \end{array}\right]
2997: ```
2998: where the inverses of A_{00} and S are applied using KSPs. The upper factorization is the inverse of
2999: ```{math}
3000: \left[\begin{array}{cc} A_{00} & A_{01} \\ 0 & S \end{array}\right]
3001: ```
3002: where again the inverses of $A_{00}$ and $S$ are applied using `KSP`s.
3004: If only one set of indices (one `IS`) is provided with `PCFieldSplitSetIS()` then the complement of that `IS`
3005: is used automatically for a second block.
3007: The fieldsplit preconditioner cannot currently be used with the `MATBAIJ` or `MATSBAIJ` data formats if the blocksize is larger than 1.
3008: Generally it should be used with the `MATAIJ` format.
3010: The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see,
3011: for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling {cite}`Wesseling2009`.
3012: One can also use `PCFIELDSPLIT`
3013: inside a smoother resulting in "Distributive Smoothers".
3015: See "A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations" {cite}`elman2008tcp`.
3017: The Constrained Pressure Preconditioner (CPR) can be implemented using `PCCOMPOSITE` with `PCGALERKIN`. CPR first solves an $R A P$ subsystem, updates the
3018: residual on all variables (`PCCompositeSetType(pc,PC_COMPOSITE_MULTIPLICATIVE)`), and then applies a simple ILU like preconditioner on all the variables.
3020: The generalized Golub-Kahan bidiagonalization preconditioner (GKB) can be applied to symmetric $2 \times 2$ block matrices of the shape
3021: ```{math}
3022: \left[\begin{array}{cc} A_{00} & A_{01} \\ A_{01}' & 0 \end{array}\right]
3023: ```
3024: 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}'$.
3025: 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_`.
3027: Developer Note:
3028: The Schur complement functionality of `PCFIELDSPLIT` should likely be factored into its own `PC` thus simplifying the implementation of the preconditioners and their
3029: user API.
3031: References:
3032: ```{bibliography}
3033: :filter: docname in docnames
3034: ```
3036: .seealso: [](sec_block_matrices), `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCLSC`,
3037: `PCFieldSplitGetSubKSP()`, `PCFieldSplitSchurGetSubKSP()`, `PCFieldSplitSetFields()`,
3038: `PCFieldSplitSetType()`, `PCFieldSplitSetIS()`, `PCFieldSplitSetSchurPre()`, `PCFieldSplitSetSchurFactType()`,
3039: `MatSchurComplementSetAinvType()`, `PCFieldSplitSetSchurScale()`, `PCFieldSplitSetDetectSaddlePoint()`
3040: M*/
3042: PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
3043: {
3044: PC_FieldSplit *jac;
3046: PetscNew(&jac);
3048: jac->bs = -1;
3049: jac->nsplits = 0;
3050: jac->type = PC_COMPOSITE_MULTIPLICATIVE;
3051: jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
3052: jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
3053: jac->schurscale = -1.0;
3054: jac->dm_splits = PETSC_TRUE;
3055: jac->detect = PETSC_FALSE;
3056: jac->gkbtol = 1e-5;
3057: jac->gkbdelay = 5;
3058: jac->gkbnu = 1;
3059: jac->gkbmaxit = 100;
3060: jac->gkbmonitor = PETSC_FALSE;
3061: jac->coordinates_set = PETSC_FALSE;
3063: pc->data = (void *)jac;
3065: pc->ops->apply = PCApply_FieldSplit;
3066: pc->ops->applytranspose = PCApplyTranspose_FieldSplit;
3067: pc->ops->setup = PCSetUp_FieldSplit;
3068: pc->ops->reset = PCReset_FieldSplit;
3069: pc->ops->destroy = PCDestroy_FieldSplit;
3070: pc->ops->setfromoptions = PCSetFromOptions_FieldSplit;
3071: pc->ops->view = PCView_FieldSplit;
3072: pc->ops->applyrichardson = NULL;
3074: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSchurGetSubKSP_C", PCFieldSplitSchurGetSubKSP_FieldSplit);
3075: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit);
3076: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetFields_C", PCFieldSplitSetFields_FieldSplit);
3077: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetIS_C", PCFieldSplitSetIS_FieldSplit);
3078: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetType_C", PCFieldSplitSetType_FieldSplit);
3079: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetBlockSize_C", PCFieldSplitSetBlockSize_FieldSplit);
3080: PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitRestrictIS_C", PCFieldSplitRestrictIS_FieldSplit);
3081: PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_FieldSplit);
3082: return 0;
3083: }