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