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