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