Actual source code: fieldsplit.c

  1: #include <petsc/private/pcimpl.h>
  2: #include <petsc/private/kspimpl.h>
  3: #include <petscdm.h>

  5: const char *const PCFieldSplitSchurPreTypes[]  = {"SELF", "SELFP", "A11", "USER", "FULL", "PCFieldSplitSchurPreType", "PC_FIELDSPLIT_SCHUR_PRE_", NULL};
  6: const char *const PCFieldSplitSchurFactTypes[] = {"DIAG", "LOWER", "UPPER", "FULL", "PCFieldSplitSchurFactType", "PC_FIELDSPLIT_SCHUR_FACT_", NULL};

  8: PetscLogEvent KSP_Solve_FS_0, KSP_Solve_FS_1, KSP_Solve_FS_S, KSP_Solve_FS_U, KSP_Solve_FS_L, KSP_Solve_FS_2, KSP_Solve_FS_3, KSP_Solve_FS_4;

 10: typedef struct _PC_FieldSplitLink *PC_FieldSplitLink;
 11: struct _PC_FieldSplitLink {
 12:   KSP               ksp;
 13:   Vec               x, y, z;
 14:   char             *splitname;
 15:   PetscInt          nfields;
 16:   PetscInt         *fields, *fields_col;
 17:   VecScatter        sctx;
 18:   IS                is, is_col;
 19:   PC_FieldSplitLink next, previous;
 20:   PetscLogEvent     event;

 22:   /* Used only when setting coordinates with PCSetCoordinates */
 23:   PetscInt   dim;
 24:   PetscInt   ndofs;
 25:   PetscReal *coords;
 26: };

 28: typedef struct {
 29:   PCCompositeType type;
 30:   PetscBool       defaultsplit; /* Flag for a system with a set of 'k' scalar fields with the same layout (and bs = k) */
 31:   PetscBool       splitdefined; /* Flag is set after the splits have been defined, to prevent more splits from being added */
 32:   PetscInt        bs;           /* Block size for IS and Mat structures */
 33:   PetscInt        nsplits;      /* Number of field divisions defined */
 34:   Vec            *x, *y, w1, w2;
 35:   Mat            *mat;    /* The diagonal block for each split */
 36:   Mat            *pmat;   /* The preconditioning diagonal block for each split */
 37:   Mat            *Afield; /* The rows of the matrix associated with each split */
 38:   PetscBool       issetup;

 40:   /* Only used when Schur complement preconditioning is used */
 41:   Mat                       B;          /* The (0,1) block */
 42:   Mat                       C;          /* The (1,0) block */
 43:   Mat                       schur;      /* The Schur complement S = A11 - A10 A00^{-1} A01, the KSP here, kspinner, is H_1 in [El08] */
 44:   Mat                       schurp;     /* Assembled approximation to S built by MatSchurComplement to be used as a preconditioning matrix when solving with S */
 45:   Mat                       schur_user; /* User-provided preconditioning matrix for the Schur complement */
 46:   PCFieldSplitSchurPreType  schurpre;   /* Determines which preconditioning matrix is used for the Schur complement */
 47:   PCFieldSplitSchurFactType schurfactorization;
 48:   KSP                       kspschur;   /* The solver for S */
 49:   KSP                       kspupper;   /* The solver for A in the upper diagonal part of the factorization (H_2 in [El08]) */
 50:   PetscScalar               schurscale; /* Scaling factor for the Schur complement solution with DIAG factorization */

 52:   /* Only used when Golub-Kahan bidiagonalization preconditioning is used */
 53:   Mat          H;           /* The modified matrix H = A00 + nu*A01*A01'              */
 54:   PetscReal    gkbtol;      /* Stopping tolerance for lower bound estimate            */
 55:   PetscInt     gkbdelay;    /* The delay window for the stopping criterion            */
 56:   PetscReal    gkbnu;       /* Parameter for augmented Lagrangian H = A + nu*A01*A01' */
 57:   PetscInt     gkbmaxit;    /* Maximum number of iterations for outer loop            */
 58:   PetscBool    gkbmonitor;  /* Monitor for gkb iterations and the lower bound error   */
 59:   PetscViewer  gkbviewer;   /* Viewer context for gkbmonitor                          */
 60:   Vec          u, v, d, Hu; /* Work vectors for the GKB algorithm                     */
 61:   PetscScalar *vecz;        /* Contains intermediate values, eg for lower bound       */

 63:   PC_FieldSplitLink head;
 64:   PetscBool         isrestrict;       /* indicates PCFieldSplitRestrictIS() has been last called on this object, hack */
 65:   PetscBool         suboptionsset;    /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */
 66:   PetscBool         dm_splits;        /* Whether to use DMCreateFieldDecomposition() whenever possible */
 67:   PetscBool         diag_use_amat;    /* Whether to extract diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
 68:   PetscBool         offdiag_use_amat; /* Whether to extract off-diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
 69:   PetscBool         detect;           /* Whether to form 2-way split by finding zero diagonal entries */
 70:   PetscBool         coordinates_set;  /* Whether PCSetCoordinates has been called */
 71: } PC_FieldSplit;

 73: /*
 74:     Note:
 75:     there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of
 76:    inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the
 77:    PC you could change this.
 78: */

 80: /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it.  This way the
 81: * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */
 82: static Mat FieldSplitSchurPre(PC_FieldSplit *jac)
 83: {
 84:   switch (jac->schurpre) {
 85:   case PC_FIELDSPLIT_SCHUR_PRE_SELF:
 86:     return jac->schur;
 87:   case PC_FIELDSPLIT_SCHUR_PRE_SELFP:
 88:     return jac->schurp;
 89:   case PC_FIELDSPLIT_SCHUR_PRE_A11:
 90:     return jac->pmat[1];
 91:   case PC_FIELDSPLIT_SCHUR_PRE_FULL: /* We calculate this and store it in schur_user */
 92:   case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */
 93:   default:
 94:     return jac->schur_user ? jac->schur_user : jac->pmat[1];
 95:   }
 96: }

 98: #include <petscdraw.h>
 99: static PetscErrorCode PCView_FieldSplit(PC pc, PetscViewer viewer)
100: {
101:   PC_FieldSplit    *jac = (PC_FieldSplit *)pc->data;
102:   PetscBool         iascii, isdraw;
103:   PetscInt          i, j;
104:   PC_FieldSplitLink ilink = jac->head;

106:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
107:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
108:   if (iascii) {
109:     if (jac->bs > 0) {
110:       PetscViewerASCIIPrintf(viewer, "  FieldSplit with %s composition: total splits = %" PetscInt_FMT ", blocksize = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits, jac->bs);
111:     } else {
112:       PetscViewerASCIIPrintf(viewer, "  FieldSplit with %s composition: total splits = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits);
113:     }
114:     if (pc->useAmat) PetscViewerASCIIPrintf(viewer, "  using Amat (not Pmat) as operator for blocks\n");
115:     if (jac->diag_use_amat) PetscViewerASCIIPrintf(viewer, "  using Amat (not Pmat) as operator for diagonal blocks\n");
116:     if (jac->offdiag_use_amat) PetscViewerASCIIPrintf(viewer, "  using Amat (not Pmat) as operator for off-diagonal blocks\n");
117:     PetscViewerASCIIPrintf(viewer, "  Solver info for each split is in the following KSP objects:\n");
118:     for (i = 0; i < jac->nsplits; i++) {
119:       if (ilink->fields) {
120:         PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Fields ", i);
121:         PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
122:         for (j = 0; j < ilink->nfields; j++) {
123:           if (j > 0) PetscViewerASCIIPrintf(viewer, ",");
124:           PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, ilink->fields[j]);
125:         }
126:         PetscViewerASCIIPrintf(viewer, "\n");
127:         PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
128:       } else {
129:         PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Defined by IS\n", i);
130:       }
131:       KSPView(ilink->ksp, viewer);
132:       ilink = ilink->next;
133:     }
134:   }

136:   if (isdraw) {
137:     PetscDraw draw;
138:     PetscReal x, y, w, wd;

140:     PetscViewerDrawGetDraw(viewer, 0, &draw);
141:     PetscDrawGetCurrentPoint(draw, &x, &y);
142:     w  = 2 * PetscMin(1.0 - x, x);
143:     wd = w / (jac->nsplits + 1);
144:     x  = x - wd * (jac->nsplits - 1) / 2.0;
145:     for (i = 0; i < jac->nsplits; i++) {
146:       PetscDrawPushCurrentPoint(draw, x, y);
147:       KSPView(ilink->ksp, viewer);
148:       PetscDrawPopCurrentPoint(draw);
149:       x += wd;
150:       ilink = ilink->next;
151:     }
152:   }
153:   return 0;
154: }

156: static PetscErrorCode PCView_FieldSplit_Schur(PC pc, PetscViewer viewer)
157: {
158:   PC_FieldSplit             *jac = (PC_FieldSplit *)pc->data;
159:   PetscBool                  iascii, isdraw;
160:   PetscInt                   i, j;
161:   PC_FieldSplitLink          ilink = jac->head;
162:   MatSchurComplementAinvType atype;

164:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
165:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
166:   if (iascii) {
167:     if (jac->bs > 0) {
168:       PetscViewerASCIIPrintf(viewer, "  FieldSplit with Schur preconditioner, blocksize = %" PetscInt_FMT ", factorization %s\n", jac->bs, PCFieldSplitSchurFactTypes[jac->schurfactorization]);
169:     } else {
170:       PetscViewerASCIIPrintf(viewer, "  FieldSplit with Schur preconditioner, factorization %s\n", PCFieldSplitSchurFactTypes[jac->schurfactorization]);
171:     }
172:     if (pc->useAmat) PetscViewerASCIIPrintf(viewer, "  using Amat (not Pmat) as operator for blocks\n");
173:     switch (jac->schurpre) {
174:     case PC_FIELDSPLIT_SCHUR_PRE_SELF:
175:       PetscViewerASCIIPrintf(viewer, "  Preconditioner for the Schur complement formed from S itself\n");
176:       break;
177:     case PC_FIELDSPLIT_SCHUR_PRE_SELFP:
178:       MatSchurComplementGetAinvType(jac->schur, &atype);
179:       PetscViewerASCIIPrintf(viewer, "  Preconditioner for the Schur complement formed from Sp, an assembled approximation to S, which uses A00's %sinverse\n", atype == MAT_SCHUR_COMPLEMENT_AINV_DIAG ? "diagonal's " : (atype == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG ? "block diagonal's " : (atype == MAT_SCHUR_COMPLEMENT_AINV_FULL ? "full " : "lumped diagonal's ")));
180:       break;
181:     case PC_FIELDSPLIT_SCHUR_PRE_A11:
182:       PetscViewerASCIIPrintf(viewer, "  Preconditioner for the Schur complement formed from A11\n");
183:       break;
184:     case PC_FIELDSPLIT_SCHUR_PRE_FULL:
185:       PetscViewerASCIIPrintf(viewer, "  Preconditioner for the Schur complement formed from the exact Schur complement\n");
186:       break;
187:     case PC_FIELDSPLIT_SCHUR_PRE_USER:
188:       if (jac->schur_user) {
189:         PetscViewerASCIIPrintf(viewer, "  Preconditioner for the Schur complement formed from user provided matrix\n");
190:       } else {
191:         PetscViewerASCIIPrintf(viewer, "  Preconditioner for the Schur complement formed from A11\n");
192:       }
193:       break;
194:     default:
195:       SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre);
196:     }
197:     PetscViewerASCIIPrintf(viewer, "  Split info:\n");
198:     PetscViewerASCIIPushTab(viewer);
199:     for (i = 0; i < jac->nsplits; i++) {
200:       if (ilink->fields) {
201:         PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Fields ", i);
202:         PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
203:         for (j = 0; j < ilink->nfields; j++) {
204:           if (j > 0) PetscViewerASCIIPrintf(viewer, ",");
205:           PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, ilink->fields[j]);
206:         }
207:         PetscViewerASCIIPrintf(viewer, "\n");
208:         PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
209:       } else {
210:         PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Defined by IS\n", i);
211:       }
212:       ilink = ilink->next;
213:     }
214:     PetscViewerASCIIPrintf(viewer, "KSP solver for A00 block\n");
215:     PetscViewerASCIIPushTab(viewer);
216:     if (jac->head) {
217:       KSPView(jac->head->ksp, viewer);
218:     } else PetscViewerASCIIPrintf(viewer, "  not yet available\n");
219:     PetscViewerASCIIPopTab(viewer);
220:     if (jac->head && jac->kspupper != jac->head->ksp) {
221:       PetscViewerASCIIPrintf(viewer, "KSP solver for upper A00 in upper triangular factor \n");
222:       PetscViewerASCIIPushTab(viewer);
223:       if (jac->kspupper) KSPView(jac->kspupper, viewer);
224:       else PetscViewerASCIIPrintf(viewer, "  not yet available\n");
225:       PetscViewerASCIIPopTab(viewer);
226:     }
227:     PetscViewerASCIIPrintf(viewer, "KSP solver for S = A11 - A10 inv(A00) A01 \n");
228:     PetscViewerASCIIPushTab(viewer);
229:     if (jac->kspschur) {
230:       KSPView(jac->kspschur, viewer);
231:     } else {
232:       PetscViewerASCIIPrintf(viewer, "  not yet available\n");
233:     }
234:     PetscViewerASCIIPopTab(viewer);
235:     PetscViewerASCIIPopTab(viewer);
236:   } else if (isdraw && jac->head) {
237:     PetscDraw draw;
238:     PetscReal x, y, w, wd, h;
239:     PetscInt  cnt = 2;
240:     char      str[32];

242:     PetscViewerDrawGetDraw(viewer, 0, &draw);
243:     PetscDrawGetCurrentPoint(draw, &x, &y);
244:     if (jac->kspupper != jac->head->ksp) cnt++;
245:     w  = 2 * PetscMin(1.0 - x, x);
246:     wd = w / (cnt + 1);

248:     PetscSNPrintf(str, 32, "Schur fact. %s", PCFieldSplitSchurFactTypes[jac->schurfactorization]);
249:     PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h);
250:     y -= h;
251:     if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_USER && !jac->schur_user) {
252:       PetscSNPrintf(str, 32, "Prec. for Schur from %s", PCFieldSplitSchurPreTypes[PC_FIELDSPLIT_SCHUR_PRE_A11]);
253:     } else {
254:       PetscSNPrintf(str, 32, "Prec. for Schur from %s", PCFieldSplitSchurPreTypes[jac->schurpre]);
255:     }
256:     PetscDrawStringBoxed(draw, x + wd * (cnt - 1) / 2.0, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h);
257:     y -= h;
258:     x = x - wd * (cnt - 1) / 2.0;

260:     PetscDrawPushCurrentPoint(draw, x, y);
261:     KSPView(jac->head->ksp, viewer);
262:     PetscDrawPopCurrentPoint(draw);
263:     if (jac->kspupper != jac->head->ksp) {
264:       x += wd;
265:       PetscDrawPushCurrentPoint(draw, x, y);
266:       KSPView(jac->kspupper, viewer);
267:       PetscDrawPopCurrentPoint(draw);
268:     }
269:     x += wd;
270:     PetscDrawPushCurrentPoint(draw, x, y);
271:     KSPView(jac->kspschur, viewer);
272:     PetscDrawPopCurrentPoint(draw);
273:   }
274:   return 0;
275: }

277: static PetscErrorCode PCView_FieldSplit_GKB(PC pc, PetscViewer viewer)
278: {
279:   PC_FieldSplit    *jac = (PC_FieldSplit *)pc->data;
280:   PetscBool         iascii, isdraw;
281:   PetscInt          i, j;
282:   PC_FieldSplitLink ilink = jac->head;

284:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
285:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
286:   if (iascii) {
287:     if (jac->bs > 0) {
288:       PetscViewerASCIIPrintf(viewer, "  FieldSplit with %s composition: total splits = %" PetscInt_FMT ", blocksize = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits, jac->bs);
289:     } else {
290:       PetscViewerASCIIPrintf(viewer, "  FieldSplit with %s composition: total splits = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits);
291:     }
292:     if (pc->useAmat) PetscViewerASCIIPrintf(viewer, "  using Amat (not Pmat) as operator for blocks\n");
293:     if (jac->diag_use_amat) PetscViewerASCIIPrintf(viewer, "  using Amat (not Pmat) as operator for diagonal blocks\n");
294:     if (jac->offdiag_use_amat) PetscViewerASCIIPrintf(viewer, "  using Amat (not Pmat) as operator for off-diagonal blocks\n");

296:     PetscViewerASCIIPrintf(viewer, "  Stopping tolerance=%.1e, delay in error estimate=%" PetscInt_FMT ", maximum iterations=%" PetscInt_FMT "\n", (double)jac->gkbtol, jac->gkbdelay, jac->gkbmaxit);
297:     PetscViewerASCIIPrintf(viewer, "  Solver info for H = A00 + nu*A01*A01' matrix:\n");
298:     PetscViewerASCIIPushTab(viewer);

300:     if (ilink->fields) {
301:       PetscViewerASCIIPrintf(viewer, "Split number 0 Fields ");
302:       PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
303:       for (j = 0; j < ilink->nfields; j++) {
304:         if (j > 0) PetscViewerASCIIPrintf(viewer, ",");
305:         PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, ilink->fields[j]);
306:       }
307:       PetscViewerASCIIPrintf(viewer, "\n");
308:       PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
309:     } else {
310:       PetscViewerASCIIPrintf(viewer, "Split number 0 Defined by IS\n");
311:     }
312:     KSPView(ilink->ksp, viewer);

314:     PetscViewerASCIIPopTab(viewer);
315:   }

317:   if (isdraw) {
318:     PetscDraw draw;
319:     PetscReal x, y, w, wd;

321:     PetscViewerDrawGetDraw(viewer, 0, &draw);
322:     PetscDrawGetCurrentPoint(draw, &x, &y);
323:     w  = 2 * PetscMin(1.0 - x, x);
324:     wd = w / (jac->nsplits + 1);
325:     x  = x - wd * (jac->nsplits - 1) / 2.0;
326:     for (i = 0; i < jac->nsplits; i++) {
327:       PetscDrawPushCurrentPoint(draw, x, y);
328:       KSPView(ilink->ksp, viewer);
329:       PetscDrawPopCurrentPoint(draw);
330:       x += wd;
331:       ilink = ilink->next;
332:     }
333:   }
334:   return 0;
335: }

337: /* Precondition: jac->bs is set to a meaningful value */
338: static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc)
339: {
340:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
341:   PetscInt       i, nfields, *ifields, nfields_col, *ifields_col;
342:   PetscBool      flg, flg_col;
343:   char           optionname[128], splitname[8], optionname_col[128];

345:   PetscMalloc1(jac->bs, &ifields);
346:   PetscMalloc1(jac->bs, &ifields_col);
347:   for (i = 0, flg = PETSC_TRUE;; i++) {
348:     PetscSNPrintf(splitname, sizeof(splitname), "%" PetscInt_FMT, i);
349:     PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%" PetscInt_FMT "_fields", i);
350:     PetscSNPrintf(optionname_col, sizeof(optionname_col), "-pc_fieldsplit_%" PetscInt_FMT "_fields_col", i);
351:     nfields     = jac->bs;
352:     nfields_col = jac->bs;
353:     PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg);
354:     PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, optionname_col, ifields_col, &nfields_col, &flg_col);
355:     if (!flg) break;
356:     else if (flg && !flg_col) {
358:       PCFieldSplitSetFields(pc, splitname, nfields, ifields, ifields);
359:     } else {
362:       PCFieldSplitSetFields(pc, splitname, nfields, ifields, ifields_col);
363:     }
364:   }
365:   if (i > 0) {
366:     /* Makes command-line setting of splits take precedence over setting them in code.
367:        Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would
368:        create new splits, which would probably not be what the user wanted. */
369:     jac->splitdefined = PETSC_TRUE;
370:   }
371:   PetscFree(ifields);
372:   PetscFree(ifields_col);
373:   return 0;
374: }

376: static PetscErrorCode PCFieldSplitSetDefaults(PC pc)
377: {
378:   PC_FieldSplit    *jac                = (PC_FieldSplit *)pc->data;
379:   PC_FieldSplitLink ilink              = jac->head;
380:   PetscBool         fieldsplit_default = PETSC_FALSE, coupling = PETSC_FALSE;
381:   PetscInt          i;

383:   /*
384:    Kinda messy, but at least this now uses DMCreateFieldDecomposition().
385:    Should probably be rewritten.
386:    */
387:   if (!ilink) {
388:     PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_fieldsplit_detect_coupling", &coupling, NULL);
389:     if (pc->dm && jac->dm_splits && !jac->detect && !coupling) {
390:       PetscInt  numFields, f, i, j;
391:       char    **fieldNames;
392:       IS       *fields;
393:       DM       *dms;
394:       DM        subdm[128];
395:       PetscBool flg;

397:       DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms);
398:       /* Allow the user to prescribe the splits */
399:       for (i = 0, flg = PETSC_TRUE;; i++) {
400:         PetscInt ifields[128];
401:         IS       compField;
402:         char     optionname[128], splitname[8];
403:         PetscInt nfields = numFields;

405:         PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%" PetscInt_FMT "_fields", i);
406:         PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg);
407:         if (!flg) break;
409:         DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i]);
410:         if (nfields == 1) {
411:           PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField);
412:         } else {
413:           PetscSNPrintf(splitname, sizeof(splitname), "%" PetscInt_FMT, i);
414:           PCFieldSplitSetIS(pc, splitname, compField);
415:         }
416:         ISDestroy(&compField);
417:         for (j = 0; j < nfields; ++j) {
418:           f = ifields[j];
419:           PetscFree(fieldNames[f]);
420:           ISDestroy(&fields[f]);
421:         }
422:       }
423:       if (i == 0) {
424:         for (f = 0; f < numFields; ++f) {
425:           PCFieldSplitSetIS(pc, fieldNames[f], fields[f]);
426:           PetscFree(fieldNames[f]);
427:           ISDestroy(&fields[f]);
428:         }
429:       } else {
430:         for (j = 0; j < numFields; j++) DMDestroy(dms + j);
431:         PetscFree(dms);
432:         PetscMalloc1(i, &dms);
433:         for (j = 0; j < i; ++j) dms[j] = subdm[j];
434:       }
435:       PetscFree(fieldNames);
436:       PetscFree(fields);
437:       if (dms) {
438:         PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n");
439:         for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) {
440:           const char *prefix;
441:           PetscObjectGetOptionsPrefix((PetscObject)(ilink->ksp), &prefix);
442:           PetscObjectSetOptionsPrefix((PetscObject)(dms[i]), prefix);
443:           KSPSetDM(ilink->ksp, dms[i]);
444:           KSPSetDMActive(ilink->ksp, PETSC_FALSE);
445:           {
446:             PetscErrorCode (*func)(KSP, Mat, Mat, void *);
447:             void *ctx;

449:             DMKSPGetComputeOperators(pc->dm, &func, &ctx);
450:             DMKSPSetComputeOperators(dms[i], func, ctx);
451:           }
452:           PetscObjectIncrementTabLevel((PetscObject)dms[i], (PetscObject)ilink->ksp, 0);
453:           DMDestroy(&dms[i]);
454:         }
455:         PetscFree(dms);
456:       }
457:     } else {
458:       if (jac->bs <= 0) {
459:         if (pc->pmat) {
460:           MatGetBlockSize(pc->pmat, &jac->bs);
461:         } else jac->bs = 1;
462:       }

464:       if (jac->detect) {
465:         IS       zerodiags, rest;
466:         PetscInt nmin, nmax;

468:         MatGetOwnershipRange(pc->mat, &nmin, &nmax);
469:         if (jac->diag_use_amat) {
470:           MatFindZeroDiagonals(pc->mat, &zerodiags);
471:         } else {
472:           MatFindZeroDiagonals(pc->pmat, &zerodiags);
473:         }
474:         ISComplement(zerodiags, nmin, nmax, &rest);
475:         PCFieldSplitSetIS(pc, "0", rest);
476:         PCFieldSplitSetIS(pc, "1", zerodiags);
477:         ISDestroy(&zerodiags);
478:         ISDestroy(&rest);
479:       } else if (coupling) {
480:         IS       coupling, rest;
481:         PetscInt nmin, nmax;

483:         MatGetOwnershipRange(pc->mat, &nmin, &nmax);
484:         if (jac->offdiag_use_amat) {
485:           MatFindOffBlockDiagonalEntries(pc->mat, &coupling);
486:         } else {
487:           MatFindOffBlockDiagonalEntries(pc->pmat, &coupling);
488:         }
489:         ISCreateStride(PetscObjectComm((PetscObject)pc->mat), nmax - nmin, nmin, 1, &rest);
490:         ISSetIdentity(rest);
491:         PCFieldSplitSetIS(pc, "0", rest);
492:         PCFieldSplitSetIS(pc, "1", coupling);
493:         ISDestroy(&coupling);
494:         ISDestroy(&rest);
495:       } else {
496:         PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_fieldsplit_default", &fieldsplit_default, NULL);
497:         if (!fieldsplit_default) {
498:           /* Allow user to set fields from command line,  if bs was known at the time of PCSetFromOptions_FieldSplit()
499:            then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */
500:           PCFieldSplitSetRuntimeSplits_Private(pc);
501:           if (jac->splitdefined) PetscInfo(pc, "Splits defined using the options database\n");
502:         }
503:         if ((fieldsplit_default || !jac->splitdefined) && !jac->isrestrict) {
504:           Mat       M = pc->pmat;
505:           PetscBool isnest;

507:           PetscInfo(pc, "Using default splitting of fields\n");
508:           PetscObjectTypeCompare((PetscObject)pc->pmat, MATNEST, &isnest);
509:           if (!isnest) {
510:             M = pc->mat;
511:             PetscObjectTypeCompare((PetscObject)pc->mat, MATNEST, &isnest);
512:           }
513:           if (isnest) {
514:             IS      *fields;
515:             PetscInt nf;

517:             MatNestGetSize(M, &nf, NULL);
518:             PetscMalloc1(nf, &fields);
519:             MatNestGetISs(M, fields, NULL);
520:             for (i = 0; i < nf; i++) PCFieldSplitSetIS(pc, NULL, fields[i]);
521:             PetscFree(fields);
522:           } else {
523:             for (i = 0; i < jac->bs; i++) {
524:               char splitname[8];
525:               PetscSNPrintf(splitname, sizeof(splitname), "%" PetscInt_FMT, i);
526:               PCFieldSplitSetFields(pc, splitname, 1, &i, &i);
527:             }
528:             jac->defaultsplit = PETSC_TRUE;
529:           }
530:         }
531:       }
532:     }
533:   } else if (jac->nsplits == 1) {
534:     if (ilink->is) {
535:       IS       is2;
536:       PetscInt nmin, nmax;

538:       MatGetOwnershipRange(pc->mat, &nmin, &nmax);
539:       ISComplement(ilink->is, nmin, nmax, &is2);
540:       PCFieldSplitSetIS(pc, "1", is2);
541:       ISDestroy(&is2);
542:     } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Must provide at least two sets of fields to PCFieldSplit()");
543:   }

546:   return 0;
547: }

549: static PetscErrorCode MatGolubKahanComputeExplicitOperator(Mat A, Mat B, Mat C, Mat *H, PetscReal gkbnu)
550: {
551:   Mat       BT, T;
552:   PetscReal nrmT, nrmB;

554:   MatHermitianTranspose(C, MAT_INITIAL_MATRIX, &T); /* Test if augmented matrix is symmetric */
555:   MatAXPY(T, -1.0, B, DIFFERENT_NONZERO_PATTERN);
556:   MatNorm(T, NORM_1, &nrmT);
557:   MatNorm(B, NORM_1, &nrmB);

560:   /* Compute augmented Lagrangian matrix H = A00 + nu*A01*A01'. This corresponds to */
561:   /* setting N := 1/nu*I in [Ar13].                                                 */
562:   MatHermitianTranspose(B, MAT_INITIAL_MATRIX, &BT);
563:   MatMatMult(B, BT, MAT_INITIAL_MATRIX, PETSC_DEFAULT, H); /* H = A01*A01'          */
564:   MatAYPX(*H, gkbnu, A, DIFFERENT_NONZERO_PATTERN);        /* H = A00 + nu*A01*A01' */

566:   MatDestroy(&BT);
567:   MatDestroy(&T);
568:   return 0;
569: }

571: PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(PetscOptions, const char pre[], const char name[], const char *value[], PetscBool *flg);

573: static PetscErrorCode PCSetUp_FieldSplit(PC pc)
574: {
575:   PC_FieldSplit    *jac = (PC_FieldSplit *)pc->data;
576:   PC_FieldSplitLink ilink;
577:   PetscInt          i, nsplit;
578:   PetscBool         sorted, sorted_col;

580:   pc->failedreason = PC_NOERROR;
581:   PCFieldSplitSetDefaults(pc);
582:   nsplit = jac->nsplits;
583:   ilink  = jac->head;

585:   /* get the matrices for each split */
586:   if (!jac->issetup) {
587:     PetscInt rstart, rend, nslots, bs;

589:     jac->issetup = PETSC_TRUE;

591:     /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */
592:     if (jac->defaultsplit || !ilink->is) {
593:       if (jac->bs <= 0) jac->bs = nsplit;
594:     }

596:     /*  MatCreateSubMatrix() for [S]BAIJ matrices can only work if the indices include entire blocks of the matrix */
597:     MatGetBlockSize(pc->pmat, &bs);
598:     if (bs > 1 && (jac->bs <= bs || jac->bs % bs)) {
599:       PetscBool blk;

601:       PetscObjectTypeCompareAny((PetscObject)pc->pmat, &blk, MATBAIJ, MATSBAIJ, MATSEQBAIJ, MATSEQSBAIJ, MATMPIBAIJ, MATMPISBAIJ, NULL);
603:     }

605:     bs = jac->bs;
606:     MatGetOwnershipRange(pc->pmat, &rstart, &rend);
607:     nslots = (rend - rstart) / bs;
608:     for (i = 0; i < nsplit; i++) {
609:       if (jac->defaultsplit) {
610:         ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + i, nsplit, &ilink->is);
611:         ISDuplicate(ilink->is, &ilink->is_col);
612:       } else if (!ilink->is) {
613:         if (ilink->nfields > 1) {
614:           PetscInt *ii, *jj, j, k, nfields = ilink->nfields, *fields = ilink->fields, *fields_col = ilink->fields_col;
615:           PetscMalloc1(ilink->nfields * nslots, &ii);
616:           PetscMalloc1(ilink->nfields * nslots, &jj);
617:           for (j = 0; j < nslots; j++) {
618:             for (k = 0; k < nfields; k++) {
619:               ii[nfields * j + k] = rstart + bs * j + fields[k];
620:               jj[nfields * j + k] = rstart + bs * j + fields_col[k];
621:             }
622:           }
623:           ISCreateGeneral(PetscObjectComm((PetscObject)pc), nslots * nfields, ii, PETSC_OWN_POINTER, &ilink->is);
624:           ISCreateGeneral(PetscObjectComm((PetscObject)pc), nslots * nfields, jj, PETSC_OWN_POINTER, &ilink->is_col);
625:           ISSetBlockSize(ilink->is, nfields);
626:           ISSetBlockSize(ilink->is_col, nfields);
627:         } else {
628:           ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + ilink->fields[0], bs, &ilink->is);
629:           ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + ilink->fields_col[0], bs, &ilink->is_col);
630:         }
631:       }
632:       ISSorted(ilink->is, &sorted);
633:       if (ilink->is_col) ISSorted(ilink->is_col, &sorted_col);
635:       ilink = ilink->next;
636:     }
637:   }

639:   ilink = jac->head;
640:   if (!jac->pmat) {
641:     Vec xtmp;

643:     MatCreateVecs(pc->pmat, &xtmp, NULL);
644:     PetscMalloc1(nsplit, &jac->pmat);
645:     PetscMalloc2(nsplit, &jac->x, nsplit, &jac->y);
646:     for (i = 0; i < nsplit; i++) {
647:       MatNullSpace sp;

649:       /* Check for preconditioning matrix attached to IS */
650:       PetscObjectQuery((PetscObject)ilink->is, "pmat", (PetscObject *)&jac->pmat[i]);
651:       if (jac->pmat[i]) {
652:         PetscObjectReference((PetscObject)jac->pmat[i]);
653:         if (jac->type == PC_COMPOSITE_SCHUR) {
654:           jac->schur_user = jac->pmat[i];

656:           PetscObjectReference((PetscObject)jac->schur_user);
657:         }
658:       } else {
659:         const char *prefix;
660:         MatCreateSubMatrix(pc->pmat, ilink->is, ilink->is_col, MAT_INITIAL_MATRIX, &jac->pmat[i]);
661:         KSPGetOptionsPrefix(ilink->ksp, &prefix);
662:         MatSetOptionsPrefix(jac->pmat[i], prefix);
663:         MatViewFromOptions(jac->pmat[i], NULL, "-mat_view");
664:       }
665:       /* create work vectors for each split */
666:       MatCreateVecs(jac->pmat[i], &jac->x[i], &jac->y[i]);
667:       ilink->x = jac->x[i];
668:       ilink->y = jac->y[i];
669:       ilink->z = NULL;
670:       /* compute scatter contexts needed by multiplicative versions and non-default splits */
671:       VecScatterCreate(xtmp, ilink->is, jac->x[i], NULL, &ilink->sctx);
672:       PetscObjectQuery((PetscObject)ilink->is, "nearnullspace", (PetscObject *)&sp);
673:       if (sp) MatSetNearNullSpace(jac->pmat[i], sp);
674:       ilink = ilink->next;
675:     }
676:     VecDestroy(&xtmp);
677:   } else {
678:     MatReuse scall;
679:     if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
680:       for (i = 0; i < nsplit; i++) MatDestroy(&jac->pmat[i]);
681:       scall = MAT_INITIAL_MATRIX;
682:     } else scall = MAT_REUSE_MATRIX;

684:     for (i = 0; i < nsplit; i++) {
685:       Mat pmat;

687:       /* Check for preconditioning matrix attached to IS */
688:       PetscObjectQuery((PetscObject)ilink->is, "pmat", (PetscObject *)&pmat);
689:       if (!pmat) MatCreateSubMatrix(pc->pmat, ilink->is, ilink->is_col, scall, &jac->pmat[i]);
690:       ilink = ilink->next;
691:     }
692:   }
693:   if (jac->diag_use_amat) {
694:     ilink = jac->head;
695:     if (!jac->mat) {
696:       PetscMalloc1(nsplit, &jac->mat);
697:       for (i = 0; i < nsplit; i++) {
698:         MatCreateSubMatrix(pc->mat, ilink->is, ilink->is_col, MAT_INITIAL_MATRIX, &jac->mat[i]);
699:         ilink = ilink->next;
700:       }
701:     } else {
702:       MatReuse scall;
703:       if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
704:         for (i = 0; i < nsplit; i++) MatDestroy(&jac->mat[i]);
705:         scall = MAT_INITIAL_MATRIX;
706:       } else scall = MAT_REUSE_MATRIX;

708:       for (i = 0; i < nsplit; i++) {
709:         MatCreateSubMatrix(pc->mat, ilink->is, ilink->is_col, scall, &jac->mat[i]);
710:         ilink = ilink->next;
711:       }
712:     }
713:   } else {
714:     jac->mat = jac->pmat;
715:   }

717:   /* Check for null space attached to IS */
718:   ilink = jac->head;
719:   for (i = 0; i < nsplit; i++) {
720:     MatNullSpace sp;

722:     PetscObjectQuery((PetscObject)ilink->is, "nullspace", (PetscObject *)&sp);
723:     if (sp) MatSetNullSpace(jac->mat[i], sp);
724:     ilink = ilink->next;
725:   }

727:   if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR && jac->type != PC_COMPOSITE_GKB) {
728:     /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */
729:     /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */
730:     ilink = jac->head;
731:     if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) {
732:       /* special case need where Afield[0] is not needed and only certain columns of Afield[1] are needed since update is only on those rows of the solution */
733:       if (!jac->Afield) {
734:         PetscCalloc1(nsplit, &jac->Afield);
735:         if (jac->offdiag_use_amat) {
736:           MatCreateSubMatrix(pc->mat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->Afield[1]);
737:         } else {
738:           MatCreateSubMatrix(pc->pmat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->Afield[1]);
739:         }
740:       } else {
741:         MatReuse scall;

743:         if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
744:           MatDestroy(&jac->Afield[1]);
745:           scall = MAT_INITIAL_MATRIX;
746:         } else scall = MAT_REUSE_MATRIX;

748:         if (jac->offdiag_use_amat) {
749:           MatCreateSubMatrix(pc->mat, ilink->next->is, ilink->is, scall, &jac->Afield[1]);
750:         } else {
751:           MatCreateSubMatrix(pc->pmat, ilink->next->is, ilink->is, scall, &jac->Afield[1]);
752:         }
753:       }
754:     } else {
755:       if (!jac->Afield) {
756:         PetscMalloc1(nsplit, &jac->Afield);
757:         for (i = 0; i < nsplit; i++) {
758:           if (jac->offdiag_use_amat) {
759:             MatCreateSubMatrix(pc->mat, ilink->is, NULL, MAT_INITIAL_MATRIX, &jac->Afield[i]);
760:           } else {
761:             MatCreateSubMatrix(pc->pmat, ilink->is, NULL, MAT_INITIAL_MATRIX, &jac->Afield[i]);
762:           }
763:           ilink = ilink->next;
764:         }
765:       } else {
766:         MatReuse scall;
767:         if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
768:           for (i = 0; i < nsplit; i++) MatDestroy(&jac->Afield[i]);
769:           scall = MAT_INITIAL_MATRIX;
770:         } else scall = MAT_REUSE_MATRIX;

772:         for (i = 0; i < nsplit; i++) {
773:           if (jac->offdiag_use_amat) {
774:             MatCreateSubMatrix(pc->mat, ilink->is, NULL, scall, &jac->Afield[i]);
775:           } else {
776:             MatCreateSubMatrix(pc->pmat, ilink->is, NULL, scall, &jac->Afield[i]);
777:           }
778:           ilink = ilink->next;
779:         }
780:       }
781:     }
782:   }

784:   if (jac->type == PC_COMPOSITE_SCHUR) {
785:     IS          ccis;
786:     PetscBool   isset, isspd;
787:     PetscInt    rstart, rend;
788:     char        lscname[256];
789:     PetscObject LSC_L;


793:     /* If pc->mat is SPD, don't scale by -1 the Schur complement */
794:     if (jac->schurscale == (PetscScalar)-1.0) {
795:       MatIsSPDKnown(pc->pmat, &isset, &isspd);
796:       jac->schurscale = (isset && isspd) ? 1.0 : -1.0;
797:     }

799:     /* When extracting off-diagonal submatrices, we take complements from this range */
800:     MatGetOwnershipRangeColumn(pc->mat, &rstart, &rend);

802:     if (jac->schur) {
803:       KSP      kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper;
804:       MatReuse scall;

806:       if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
807:         scall = MAT_INITIAL_MATRIX;
808:         MatDestroy(&jac->B);
809:         MatDestroy(&jac->C);
810:       } else scall = MAT_REUSE_MATRIX;

812:       MatSchurComplementGetKSP(jac->schur, &kspInner);
813:       ilink = jac->head;
814:       ISComplement(ilink->is_col, rstart, rend, &ccis);
815:       if (jac->offdiag_use_amat) {
816:         MatCreateSubMatrix(pc->mat, ilink->is, ccis, scall, &jac->B);
817:       } else {
818:         MatCreateSubMatrix(pc->pmat, ilink->is, ccis, scall, &jac->B);
819:       }
820:       ISDestroy(&ccis);
821:       ilink = ilink->next;
822:       ISComplement(ilink->is_col, rstart, rend, &ccis);
823:       if (jac->offdiag_use_amat) {
824:         MatCreateSubMatrix(pc->mat, ilink->is, ccis, scall, &jac->C);
825:       } else {
826:         MatCreateSubMatrix(pc->pmat, ilink->is, ccis, scall, &jac->C);
827:       }
828:       ISDestroy(&ccis);
829:       MatSchurComplementUpdateSubMatrices(jac->schur, jac->mat[0], jac->pmat[0], jac->B, jac->C, jac->mat[1]);
830:       if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
831:         MatDestroy(&jac->schurp);
832:         MatSchurComplementGetPmat(jac->schur, MAT_INITIAL_MATRIX, &jac->schurp);
833:       }
834:       if (kspA != kspInner) KSPSetOperators(kspA, jac->mat[0], jac->pmat[0]);
835:       if (kspUpper != kspA) KSPSetOperators(kspUpper, jac->mat[0], jac->pmat[0]);
836:       KSPSetOperators(jac->kspschur, jac->schur, FieldSplitSchurPre(jac));
837:     } else {
838:       const char  *Dprefix;
839:       char         schurprefix[256], schurmatprefix[256];
840:       char         schurtestoption[256];
841:       MatNullSpace sp;
842:       PetscBool    flg;
843:       KSP          kspt;

845:       /* extract the A01 and A10 matrices */
846:       ilink = jac->head;
847:       ISComplement(ilink->is_col, rstart, rend, &ccis);
848:       if (jac->offdiag_use_amat) {
849:         MatCreateSubMatrix(pc->mat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->B);
850:       } else {
851:         MatCreateSubMatrix(pc->pmat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->B);
852:       }
853:       ISDestroy(&ccis);
854:       ilink = ilink->next;
855:       ISComplement(ilink->is_col, rstart, rend, &ccis);
856:       if (jac->offdiag_use_amat) {
857:         MatCreateSubMatrix(pc->mat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->C);
858:       } else {
859:         MatCreateSubMatrix(pc->pmat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->C);
860:       }
861:       ISDestroy(&ccis);

863:       /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
864:       MatCreate(((PetscObject)jac->mat[0])->comm, &jac->schur);
865:       MatSetType(jac->schur, MATSCHURCOMPLEMENT);
866:       MatSchurComplementSetSubMatrices(jac->schur, jac->mat[0], jac->pmat[0], jac->B, jac->C, jac->mat[1]);
867:       PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
868:       MatSetOptionsPrefix(jac->schur, schurmatprefix);
869:       MatSchurComplementGetKSP(jac->schur, &kspt);
870:       KSPSetOptionsPrefix(kspt, schurmatprefix);

872:       /* Note: this is not true in general */
873:       MatGetNullSpace(jac->mat[1], &sp);
874:       if (sp) MatSetNullSpace(jac->schur, sp);

876:       PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);
877:       PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
878:       if (flg) {
879:         DM  dmInner;
880:         KSP kspInner;
881:         PC  pcInner;

883:         MatSchurComplementGetKSP(jac->schur, &kspInner);
884:         KSPReset(kspInner);
885:         KSPSetOperators(kspInner, jac->mat[0], jac->pmat[0]);
886:         PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
887:         /* Indent this deeper to emphasize the "inner" nature of this solver. */
888:         PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject)pc, 2);
889:         PetscObjectIncrementTabLevel((PetscObject)kspInner->pc, (PetscObject)pc, 2);
890:         KSPSetOptionsPrefix(kspInner, schurprefix);

892:         /* Set DM for new solver */
893:         KSPGetDM(jac->head->ksp, &dmInner);
894:         KSPSetDM(kspInner, dmInner);
895:         KSPSetDMActive(kspInner, PETSC_FALSE);

897:         /* Defaults to PCKSP as preconditioner */
898:         KSPGetPC(kspInner, &pcInner);
899:         PCSetType(pcInner, PCKSP);
900:         PCKSPSetKSP(pcInner, jac->head->ksp);
901:       } else {
902:         /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or
903:           * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact,
904:           * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for
905:           * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make
906:           * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used
907:           * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */
908:         KSPSetType(jac->head->ksp, KSPGMRES);
909:         MatSchurComplementSetKSP(jac->schur, jac->head->ksp);
910:       }
911:       KSPSetOperators(jac->head->ksp, jac->mat[0], jac->pmat[0]);
912:       KSPSetFromOptions(jac->head->ksp);
913:       MatSetFromOptions(jac->schur);

915:       PetscObjectTypeCompare((PetscObject)jac->schur, MATSCHURCOMPLEMENT, &flg);
916:       if (flg) { /* Need to do this otherwise PCSetUp_KSP will overwrite the amat of jac->head->ksp */
917:         KSP kspInner;
918:         PC  pcInner;

920:         MatSchurComplementGetKSP(jac->schur, &kspInner);
921:         KSPGetPC(kspInner, &pcInner);
922:         PetscObjectTypeCompare((PetscObject)pcInner, PCKSP, &flg);
923:         if (flg) {
924:           KSP ksp;

926:           PCKSPGetKSP(pcInner, &ksp);
927:           if (ksp == jac->head->ksp) PCSetUseAmat(pcInner, PETSC_TRUE);
928:         }
929:       }
930:       PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);
931:       PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
932:       if (flg) {
933:         DM dmInner;

935:         PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
936:         KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);
937:         KSPSetErrorIfNotConverged(jac->kspupper, pc->erroriffailure);
938:         KSPSetOptionsPrefix(jac->kspupper, schurprefix);
939:         PetscObjectIncrementTabLevel((PetscObject)jac->kspupper, (PetscObject)pc, 1);
940:         PetscObjectIncrementTabLevel((PetscObject)jac->kspupper->pc, (PetscObject)pc, 1);
941:         KSPGetDM(jac->head->ksp, &dmInner);
942:         KSPSetDM(jac->kspupper, dmInner);
943:         KSPSetDMActive(jac->kspupper, PETSC_FALSE);
944:         KSPSetFromOptions(jac->kspupper);
945:         KSPSetOperators(jac->kspupper, jac->mat[0], jac->pmat[0]);
946:         VecDuplicate(jac->head->x, &jac->head->z);
947:       } else {
948:         jac->kspupper = jac->head->ksp;
949:         PetscObjectReference((PetscObject)jac->head->ksp);
950:       }

952:       if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) MatSchurComplementGetPmat(jac->schur, MAT_INITIAL_MATRIX, &jac->schurp);
953:       KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspschur);
954:       KSPSetErrorIfNotConverged(jac->kspschur, pc->erroriffailure);
955:       PetscObjectIncrementTabLevel((PetscObject)jac->kspschur, (PetscObject)pc, 1);
956:       if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
957:         PC pcschur;
958:         KSPGetPC(jac->kspschur, &pcschur);
959:         PCSetType(pcschur, PCNONE);
960:         /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
961:       } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) {
962:         MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user);
963:       }
964:       KSPSetOperators(jac->kspschur, jac->schur, FieldSplitSchurPre(jac));
965:       KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);
966:       KSPSetOptionsPrefix(jac->kspschur, Dprefix);
967:       /* propagate DM */
968:       {
969:         DM sdm;
970:         KSPGetDM(jac->head->next->ksp, &sdm);
971:         if (sdm) {
972:           KSPSetDM(jac->kspschur, sdm);
973:           KSPSetDMActive(jac->kspschur, PETSC_FALSE);
974:         }
975:       }
976:       /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
977:       /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */
978:       KSPSetFromOptions(jac->kspschur);
979:     }
980:     MatAssemblyBegin(jac->schur, MAT_FINAL_ASSEMBLY);
981:     MatAssemblyEnd(jac->schur, MAT_FINAL_ASSEMBLY);

983:     /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */
984:     PetscSNPrintf(lscname, sizeof(lscname), "%s_LSC_L", ilink->splitname);
985:     PetscObjectQuery((PetscObject)pc->mat, lscname, (PetscObject *)&LSC_L);
986:     if (!LSC_L) PetscObjectQuery((PetscObject)pc->pmat, lscname, (PetscObject *)&LSC_L);
987:     if (LSC_L) PetscObjectCompose((PetscObject)jac->schur, "LSC_L", (PetscObject)LSC_L);
988:     PetscSNPrintf(lscname, sizeof(lscname), "%s_LSC_Lp", ilink->splitname);
989:     PetscObjectQuery((PetscObject)pc->pmat, lscname, (PetscObject *)&LSC_L);
990:     if (!LSC_L) PetscObjectQuery((PetscObject)pc->mat, lscname, (PetscObject *)&LSC_L);
991:     if (LSC_L) PetscObjectCompose((PetscObject)jac->schur, "LSC_Lp", (PetscObject)LSC_L);
992:   } else if (jac->type == PC_COMPOSITE_GKB) {
993:     IS       ccis;
994:     PetscInt rstart, rend;


998:     ilink = jac->head;

1000:     /* When extracting off-diagonal submatrices, we take complements from this range */
1001:     MatGetOwnershipRangeColumn(pc->mat, &rstart, &rend);

1003:     ISComplement(ilink->is_col, rstart, rend, &ccis);
1004:     if (jac->offdiag_use_amat) {
1005:       MatCreateSubMatrix(pc->mat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->B);
1006:     } else {
1007:       MatCreateSubMatrix(pc->pmat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->B);
1008:     }
1009:     ISDestroy(&ccis);
1010:     /* Create work vectors for GKB algorithm */
1011:     VecDuplicate(ilink->x, &jac->u);
1012:     VecDuplicate(ilink->x, &jac->Hu);
1013:     VecDuplicate(ilink->x, &jac->w2);
1014:     ilink = ilink->next;
1015:     ISComplement(ilink->is_col, rstart, rend, &ccis);
1016:     if (jac->offdiag_use_amat) {
1017:       MatCreateSubMatrix(pc->mat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->C);
1018:     } else {
1019:       MatCreateSubMatrix(pc->pmat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->C);
1020:     }
1021:     ISDestroy(&ccis);
1022:     /* Create work vectors for GKB algorithm */
1023:     VecDuplicate(ilink->x, &jac->v);
1024:     VecDuplicate(ilink->x, &jac->d);
1025:     VecDuplicate(ilink->x, &jac->w1);
1026:     MatGolubKahanComputeExplicitOperator(jac->mat[0], jac->B, jac->C, &jac->H, jac->gkbnu);
1027:     PetscCalloc1(jac->gkbdelay, &jac->vecz);

1029:     ilink = jac->head;
1030:     KSPSetOperators(ilink->ksp, jac->H, jac->H);
1031:     if (!jac->suboptionsset) KSPSetFromOptions(ilink->ksp);
1032:     /* Create gkb_monitor context */
1033:     if (jac->gkbmonitor) {
1034:       PetscInt tablevel;
1035:       PetscViewerCreate(PETSC_COMM_WORLD, &jac->gkbviewer);
1036:       PetscViewerSetType(jac->gkbviewer, PETSCVIEWERASCII);
1037:       PetscObjectGetTabLevel((PetscObject)ilink->ksp, &tablevel);
1038:       PetscViewerASCIISetTab(jac->gkbviewer, tablevel);
1039:       PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)ilink->ksp, 1);
1040:     }
1041:   } else {
1042:     /* set up the individual splits' PCs */
1043:     i     = 0;
1044:     ilink = jac->head;
1045:     while (ilink) {
1046:       KSPSetOperators(ilink->ksp, jac->mat[i], jac->pmat[i]);
1047:       /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
1048:       if (!jac->suboptionsset) KSPSetFromOptions(ilink->ksp);
1049:       i++;
1050:       ilink = ilink->next;
1051:     }
1052:   }

1054:   /* Set coordinates to the sub PC objects whenever these are set */
1055:   if (jac->coordinates_set) {
1056:     PC pc_coords;
1057:     if (jac->type == PC_COMPOSITE_SCHUR) {
1058:       // Head is first block.
1059:       KSPGetPC(jac->head->ksp, &pc_coords);
1060:       PCSetCoordinates(pc_coords, jac->head->dim, jac->head->ndofs, jac->head->coords);
1061:       // Second one is Schur block, but its KSP object is in kspschur.
1062:       KSPGetPC(jac->kspschur, &pc_coords);
1063:       PCSetCoordinates(pc_coords, jac->head->next->dim, jac->head->next->ndofs, jac->head->next->coords);
1064:     } else if (jac->type == PC_COMPOSITE_GKB) {
1065:       PetscInfo(pc, "Warning: Setting coordinates does nothing for the GKB Fieldpslit preconditioner");
1066:     } else {
1067:       ilink = jac->head;
1068:       while (ilink) {
1069:         KSPGetPC(ilink->ksp, &pc_coords);
1070:         PCSetCoordinates(pc_coords, ilink->dim, ilink->ndofs, ilink->coords);
1071:         ilink = ilink->next;
1072:       }
1073:     }
1074:   }

1076:   jac->suboptionsset = PETSC_TRUE;
1077:   return 0;
1078: }

1080: #define FieldSplitSplitSolveAdd(ilink, xx, yy) \
1081:   (VecScatterBegin(ilink->sctx, xx, ilink->x, INSERT_VALUES, SCATTER_FORWARD) || VecScatterEnd(ilink->sctx, xx, ilink->x, INSERT_VALUES, SCATTER_FORWARD) || PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL) || \
1082:    KSPSolve(ilink->ksp, ilink->x, ilink->y) || KSPCheckSolve(ilink->ksp, pc, ilink->y) || PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL) || VecScatterBegin(ilink->sctx, ilink->y, yy, ADD_VALUES, SCATTER_REVERSE) || \
1083:    VecScatterEnd(ilink->sctx, ilink->y, yy, ADD_VALUES, SCATTER_REVERSE))

1085: static PetscErrorCode PCApply_FieldSplit_Schur(PC pc, Vec x, Vec y)
1086: {
1087:   PC_FieldSplit    *jac    = (PC_FieldSplit *)pc->data;
1088:   PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
1089:   KSP               kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;

1091:   switch (jac->schurfactorization) {
1092:   case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
1093:     /* [A00 0; 0 -S], positive definite, suitable for MINRES */
1094:     VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD);
1095:     VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD);
1096:     VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD);
1097:     PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL);
1098:     KSPSolve(kspA, ilinkA->x, ilinkA->y);
1099:     KSPCheckSolve(kspA, pc, ilinkA->y);
1100:     PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL);
1101:     VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE);
1102:     VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD);
1103:     PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL);
1104:     KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y);
1105:     KSPCheckSolve(jac->kspschur, pc, ilinkD->y);
1106:     PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL);
1107:     VecScale(ilinkD->y, jac->schurscale);
1108:     VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE);
1109:     VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE);
1110:     VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE);
1111:     break;
1112:   case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
1113:     /* [A00 0; A10 S], suitable for left preconditioning */
1114:     VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD);
1115:     VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD);
1116:     PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL);
1117:     KSPSolve(kspA, ilinkA->x, ilinkA->y);
1118:     KSPCheckSolve(kspA, pc, ilinkA->y);
1119:     PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL);
1120:     MatMult(jac->C, ilinkA->y, ilinkD->x);
1121:     VecScale(ilinkD->x, -1.);
1122:     VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD);
1123:     VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE);
1124:     VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD);
1125:     PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL);
1126:     KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y);
1127:     KSPCheckSolve(jac->kspschur, pc, ilinkD->y);
1128:     PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL);
1129:     VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE);
1130:     VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE);
1131:     VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE);
1132:     break;
1133:   case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
1134:     /* [A00 A01; 0 S], suitable for right preconditioning */
1135:     VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD);
1136:     VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD);
1137:     PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL);
1138:     KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y);
1139:     KSPCheckSolve(jac->kspschur, pc, ilinkD->y);
1140:     PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL);
1141:     MatMult(jac->B, ilinkD->y, ilinkA->x);
1142:     VecScale(ilinkA->x, -1.);
1143:     VecScatterBegin(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD);
1144:     VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE);
1145:     VecScatterEnd(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD);
1146:     PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL);
1147:     KSPSolve(kspA, ilinkA->x, ilinkA->y);
1148:     KSPCheckSolve(kspA, pc, ilinkA->y);
1149:     PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL);
1150:     VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE);
1151:     VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE);
1152:     VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE);
1153:     break;
1154:   case PC_FIELDSPLIT_SCHUR_FACT_FULL:
1155:     /* [1 0; A10 A00^{-1} 1] [A00 0; 0 S] [1 A00^{-1}A01; 0 1] */
1156:     VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD);
1157:     VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD);
1158:     PetscLogEventBegin(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->y, NULL);
1159:     KSPSolve(kspLower, ilinkA->x, ilinkA->y);
1160:     KSPCheckSolve(kspLower, pc, ilinkA->y);
1161:     PetscLogEventEnd(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->y, NULL);
1162:     MatMult(jac->C, ilinkA->y, ilinkD->x);
1163:     VecScale(ilinkD->x, -1.0);
1164:     VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD);
1165:     VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD);

1167:     PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL);
1168:     KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y);
1169:     KSPCheckSolve(jac->kspschur, pc, ilinkD->y);
1170:     PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL);
1171:     VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE);

1173:     if (kspUpper == kspA) {
1174:       MatMult(jac->B, ilinkD->y, ilinkA->y);
1175:       VecAXPY(ilinkA->x, -1.0, ilinkA->y);
1176:       PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL);
1177:       KSPSolve(kspA, ilinkA->x, ilinkA->y);
1178:       KSPCheckSolve(kspA, pc, ilinkA->y);
1179:       PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL);
1180:     } else {
1181:       PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL);
1182:       KSPSolve(kspA, ilinkA->x, ilinkA->y);
1183:       KSPCheckSolve(kspA, pc, ilinkA->y);
1184:       MatMult(jac->B, ilinkD->y, ilinkA->x);
1185:       PetscLogEventBegin(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->z, NULL);
1186:       KSPSolve(kspUpper, ilinkA->x, ilinkA->z);
1187:       KSPCheckSolve(kspUpper, pc, ilinkA->z);
1188:       PetscLogEventEnd(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->z, NULL);
1189:       VecAXPY(ilinkA->y, -1.0, ilinkA->z);
1190:     }
1191:     VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE);
1192:     VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE);
1193:     VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE);
1194:   }
1195:   return 0;
1196: }

1198: static PetscErrorCode PCApply_FieldSplit(PC pc, Vec x, Vec y)
1199: {
1200:   PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
1201:   PC_FieldSplitLink ilink = jac->head;
1202:   PetscInt          cnt, bs;

1204:   if (jac->type == PC_COMPOSITE_ADDITIVE) {
1205:     if (jac->defaultsplit) {
1206:       VecGetBlockSize(x, &bs);
1208:       VecGetBlockSize(y, &bs);
1210:       VecStrideGatherAll(x, jac->x, INSERT_VALUES);
1211:       while (ilink) {
1212:         PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1213:         KSPSolve(ilink->ksp, ilink->x, ilink->y);
1214:         KSPCheckSolve(ilink->ksp, pc, ilink->y);
1215:         PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1216:         ilink = ilink->next;
1217:       }
1218:       VecStrideScatterAll(jac->y, y, INSERT_VALUES);
1219:     } else {
1220:       VecSet(y, 0.0);
1221:       while (ilink) {
1222:         FieldSplitSplitSolveAdd(ilink, x, y);
1223:         ilink = ilink->next;
1224:       }
1225:     }
1226:   } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) {
1227:     VecSet(y, 0.0);
1228:     /* solve on first block for first block variables */
1229:     VecScatterBegin(ilink->sctx, x, ilink->x, INSERT_VALUES, SCATTER_FORWARD);
1230:     VecScatterEnd(ilink->sctx, x, ilink->x, INSERT_VALUES, SCATTER_FORWARD);
1231:     PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1232:     KSPSolve(ilink->ksp, ilink->x, ilink->y);
1233:     KSPCheckSolve(ilink->ksp, pc, ilink->y);
1234:     PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1235:     VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE);
1236:     VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE);

1238:     /* compute the residual only onto second block variables using first block variables */
1239:     MatMult(jac->Afield[1], ilink->y, ilink->next->x);
1240:     ilink = ilink->next;
1241:     VecScale(ilink->x, -1.0);
1242:     VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD);
1243:     VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD);

1245:     /* solve on second block variables */
1246:     PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1247:     KSPSolve(ilink->ksp, ilink->x, ilink->y);
1248:     KSPCheckSolve(ilink->ksp, pc, ilink->y);
1249:     PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1250:     VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE);
1251:     VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE);
1252:   } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1253:     if (!jac->w1) {
1254:       VecDuplicate(x, &jac->w1);
1255:       VecDuplicate(x, &jac->w2);
1256:     }
1257:     VecSet(y, 0.0);
1258:     FieldSplitSplitSolveAdd(ilink, x, y);
1259:     cnt = 1;
1260:     while (ilink->next) {
1261:       ilink = ilink->next;
1262:       /* compute the residual only over the part of the vector needed */
1263:       MatMult(jac->Afield[cnt++], y, ilink->x);
1264:       VecScale(ilink->x, -1.0);
1265:       VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD);
1266:       VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD);
1267:       PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1268:       KSPSolve(ilink->ksp, ilink->x, ilink->y);
1269:       KSPCheckSolve(ilink->ksp, pc, ilink->y);
1270:       PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1271:       VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE);
1272:       VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE);
1273:     }
1274:     if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1275:       cnt -= 2;
1276:       while (ilink->previous) {
1277:         ilink = ilink->previous;
1278:         /* compute the residual only over the part of the vector needed */
1279:         MatMult(jac->Afield[cnt--], y, ilink->x);
1280:         VecScale(ilink->x, -1.0);
1281:         VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD);
1282:         VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD);
1283:         PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1284:         KSPSolve(ilink->ksp, ilink->x, ilink->y);
1285:         KSPCheckSolve(ilink->ksp, pc, ilink->y);
1286:         PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1287:         VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE);
1288:         VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE);
1289:       }
1290:     }
1291:   } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Unsupported or unknown composition %d", (int)jac->type);
1292:   return 0;
1293: }

1295: static PetscErrorCode PCApply_FieldSplit_GKB(PC pc, Vec x, Vec y)
1296: {
1297:   PC_FieldSplit    *jac    = (PC_FieldSplit *)pc->data;
1298:   PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
1299:   KSP               ksp = ilinkA->ksp;
1300:   Vec               u, v, Hu, d, work1, work2;
1301:   PetscScalar       alpha, z, nrmz2, *vecz;
1302:   PetscReal         lowbnd, nu, beta;
1303:   PetscInt          j, iterGKB;

1305:   VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD);
1306:   VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD);
1307:   VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD);
1308:   VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD);

1310:   u     = jac->u;
1311:   v     = jac->v;
1312:   Hu    = jac->Hu;
1313:   d     = jac->d;
1314:   work1 = jac->w1;
1315:   work2 = jac->w2;
1316:   vecz  = jac->vecz;

1318:   /* Change RHS to comply with matrix regularization H = A + nu*B*B' */
1319:   /* Add q = q + nu*B*b */
1320:   if (jac->gkbnu) {
1321:     nu = jac->gkbnu;
1322:     VecScale(ilinkD->x, jac->gkbnu);
1323:     MatMultAdd(jac->B, ilinkD->x, ilinkA->x, ilinkA->x); /* q = q + nu*B*b */
1324:   } else {
1325:     /* Situation when no augmented Lagrangian is used. Then we set inner  */
1326:     /* matrix N = I in [Ar13], and thus nu = 1.                           */
1327:     nu = 1;
1328:   }

1330:   /* Transform rhs from [q,tilde{b}] to [0,b] */
1331:   PetscLogEventBegin(ilinkA->event, ksp, ilinkA->x, ilinkA->y, NULL);
1332:   KSPSolve(ksp, ilinkA->x, ilinkA->y);
1333:   KSPCheckSolve(ksp, pc, ilinkA->y);
1334:   PetscLogEventEnd(ilinkA->event, ksp, ilinkA->x, ilinkA->y, NULL);
1335:   MatMultHermitianTranspose(jac->B, ilinkA->y, work1);
1336:   VecAXPBY(work1, 1.0 / nu, -1.0, ilinkD->x); /* c = b - B'*x        */

1338:   /* First step of algorithm */
1339:   VecNorm(work1, NORM_2, &beta); /* beta = sqrt(nu*c'*c)*/
1340:   KSPCheckDot(ksp, beta);
1341:   beta = PetscSqrtReal(nu) * beta;
1342:   VecAXPBY(v, nu / beta, 0.0, work1); /* v = nu/beta *c      */
1343:   MatMult(jac->B, v, work2);          /* u = H^{-1}*B*v      */
1344:   PetscLogEventBegin(ilinkA->event, ksp, work2, u, NULL);
1345:   KSPSolve(ksp, work2, u);
1346:   KSPCheckSolve(ksp, pc, u);
1347:   PetscLogEventEnd(ilinkA->event, ksp, work2, u, NULL);
1348:   MatMult(jac->H, u, Hu); /* alpha = u'*H*u      */
1349:   VecDot(Hu, u, &alpha);
1350:   KSPCheckDot(ksp, alpha);
1352:   alpha = PetscSqrtReal(PetscAbsScalar(alpha));
1353:   VecScale(u, 1.0 / alpha);
1354:   VecAXPBY(d, 1.0 / alpha, 0.0, v); /* v = nu/beta *c      */

1356:   z       = beta / alpha;
1357:   vecz[1] = z;

1359:   /* Computation of first iterate x(1) and p(1) */
1360:   VecAXPY(ilinkA->y, z, u);
1361:   VecCopy(d, ilinkD->y);
1362:   VecScale(ilinkD->y, -z);

1364:   iterGKB = 1;
1365:   lowbnd  = 2 * jac->gkbtol;
1366:   if (jac->gkbmonitor) PetscViewerASCIIPrintf(jac->gkbviewer, "%3" PetscInt_FMT " GKB Lower bound estimate %14.12e\n", iterGKB, (double)lowbnd);

1368:   while (iterGKB < jac->gkbmaxit && lowbnd > jac->gkbtol) {
1369:     iterGKB += 1;
1370:     MatMultHermitianTranspose(jac->B, u, work1); /* v <- nu*(B'*u-alpha/nu*v) */
1371:     VecAXPBY(v, nu, -alpha, work1);
1372:     VecNorm(v, NORM_2, &beta); /* beta = sqrt(nu)*v'*v      */
1373:     beta = beta / PetscSqrtReal(nu);
1374:     VecScale(v, 1.0 / beta);
1375:     MatMult(jac->B, v, work2); /* u <- H^{-1}*(B*v-beta*H*u) */
1376:     MatMult(jac->H, u, Hu);
1377:     VecAXPY(work2, -beta, Hu);
1378:     PetscLogEventBegin(ilinkA->event, ksp, work2, u, NULL);
1379:     KSPSolve(ksp, work2, u);
1380:     KSPCheckSolve(ksp, pc, u);
1381:     PetscLogEventEnd(ilinkA->event, ksp, work2, u, NULL);
1382:     MatMult(jac->H, u, Hu); /* alpha = u'*H*u            */
1383:     VecDot(Hu, u, &alpha);
1384:     KSPCheckDot(ksp, alpha);
1386:     alpha = PetscSqrtReal(PetscAbsScalar(alpha));
1387:     VecScale(u, 1.0 / alpha);

1389:     z       = -beta / alpha * z; /* z <- beta/alpha*z     */
1390:     vecz[0] = z;

1392:     /* Computation of new iterate x(i+1) and p(i+1) */
1393:     VecAXPBY(d, 1.0 / alpha, -beta / alpha, v); /* d = (v-beta*d)/alpha */
1394:     VecAXPY(ilinkA->y, z, u);                   /* r = r + z*u          */
1395:     VecAXPY(ilinkD->y, -z, d);                  /* p = p - z*d          */
1396:     MatMult(jac->H, ilinkA->y, Hu);             /* ||u||_H = u'*H*u     */
1397:     VecDot(Hu, ilinkA->y, &nrmz2);

1399:     /* Compute Lower Bound estimate */
1400:     if (iterGKB > jac->gkbdelay) {
1401:       lowbnd = 0.0;
1402:       for (j = 0; j < jac->gkbdelay; j++) lowbnd += PetscAbsScalar(vecz[j] * vecz[j]);
1403:       lowbnd = PetscSqrtReal(lowbnd / PetscAbsScalar(nrmz2));
1404:     }

1406:     for (j = 0; j < jac->gkbdelay - 1; j++) vecz[jac->gkbdelay - j - 1] = vecz[jac->gkbdelay - j - 2];
1407:     if (jac->gkbmonitor) PetscViewerASCIIPrintf(jac->gkbviewer, "%3" PetscInt_FMT " GKB Lower bound estimate %14.12e\n", iterGKB, (double)lowbnd);
1408:   }

1410:   VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE);
1411:   VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE);
1412:   VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE);
1413:   VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE);

1415:   return 0;
1416: }

1418: #define FieldSplitSplitSolveAddTranspose(ilink, xx, yy) \
1419:   (VecScatterBegin(ilink->sctx, xx, ilink->y, INSERT_VALUES, SCATTER_FORWARD) || VecScatterEnd(ilink->sctx, xx, ilink->y, INSERT_VALUES, SCATTER_FORWARD) || PetscLogEventBegin(ilink->event, ilink->ksp, ilink->y, ilink->x, NULL) || \
1420:    KSPSolveTranspose(ilink->ksp, ilink->y, ilink->x) || KSPCheckSolve(ilink->ksp, pc, ilink->x) || PetscLogEventEnd(ilink->event, ilink->ksp, ilink->y, ilink->x, NULL) || VecScatterBegin(ilink->sctx, ilink->x, yy, ADD_VALUES, SCATTER_REVERSE) || \
1421:    VecScatterEnd(ilink->sctx, ilink->x, yy, ADD_VALUES, SCATTER_REVERSE))

1423: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc, Vec x, Vec y)
1424: {
1425:   PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
1426:   PC_FieldSplitLink ilink = jac->head;
1427:   PetscInt          bs;

1429:   if (jac->type == PC_COMPOSITE_ADDITIVE) {
1430:     if (jac->defaultsplit) {
1431:       VecGetBlockSize(x, &bs);
1433:       VecGetBlockSize(y, &bs);
1435:       VecStrideGatherAll(x, jac->x, INSERT_VALUES);
1436:       while (ilink) {
1437:         PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1438:         KSPSolveTranspose(ilink->ksp, ilink->x, ilink->y);
1439:         KSPCheckSolve(ilink->ksp, pc, ilink->y);
1440:         PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1441:         ilink = ilink->next;
1442:       }
1443:       VecStrideScatterAll(jac->y, y, INSERT_VALUES);
1444:     } else {
1445:       VecSet(y, 0.0);
1446:       while (ilink) {
1447:         FieldSplitSplitSolveAddTranspose(ilink, x, y);
1448:         ilink = ilink->next;
1449:       }
1450:     }
1451:   } else {
1452:     if (!jac->w1) {
1453:       VecDuplicate(x, &jac->w1);
1454:       VecDuplicate(x, &jac->w2);
1455:     }
1456:     VecSet(y, 0.0);
1457:     if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1458:       FieldSplitSplitSolveAddTranspose(ilink, x, y);
1459:       while (ilink->next) {
1460:         ilink = ilink->next;
1461:         MatMultTranspose(pc->mat, y, jac->w1);
1462:         VecWAXPY(jac->w2, -1.0, jac->w1, x);
1463:         FieldSplitSplitSolveAddTranspose(ilink, jac->w2, y);
1464:       }
1465:       while (ilink->previous) {
1466:         ilink = ilink->previous;
1467:         MatMultTranspose(pc->mat, y, jac->w1);
1468:         VecWAXPY(jac->w2, -1.0, jac->w1, x);
1469:         FieldSplitSplitSolveAddTranspose(ilink, jac->w2, y);
1470:       }
1471:     } else {
1472:       while (ilink->next) { /* get to last entry in linked list */
1473:         ilink = ilink->next;
1474:       }
1475:       FieldSplitSplitSolveAddTranspose(ilink, x, y);
1476:       while (ilink->previous) {
1477:         ilink = ilink->previous;
1478:         MatMultTranspose(pc->mat, y, jac->w1);
1479:         VecWAXPY(jac->w2, -1.0, jac->w1, x);
1480:         FieldSplitSplitSolveAddTranspose(ilink, jac->w2, y);
1481:       }
1482:     }
1483:   }
1484:   return 0;
1485: }

1487: static PetscErrorCode PCReset_FieldSplit(PC pc)
1488: {
1489:   PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
1490:   PC_FieldSplitLink ilink = jac->head, next;

1492:   while (ilink) {
1493:     KSPDestroy(&ilink->ksp);
1494:     VecDestroy(&ilink->x);
1495:     VecDestroy(&ilink->y);
1496:     VecDestroy(&ilink->z);
1497:     VecScatterDestroy(&ilink->sctx);
1498:     ISDestroy(&ilink->is);
1499:     ISDestroy(&ilink->is_col);
1500:     PetscFree(ilink->splitname);
1501:     PetscFree(ilink->fields);
1502:     PetscFree(ilink->fields_col);
1503:     next = ilink->next;
1504:     PetscFree(ilink);
1505:     ilink = next;
1506:   }
1507:   jac->head = NULL;
1508:   PetscFree2(jac->x, jac->y);
1509:   if (jac->mat && jac->mat != jac->pmat) {
1510:     MatDestroyMatrices(jac->nsplits, &jac->mat);
1511:   } else if (jac->mat) {
1512:     jac->mat = NULL;
1513:   }
1514:   if (jac->pmat) MatDestroyMatrices(jac->nsplits, &jac->pmat);
1515:   if (jac->Afield) MatDestroyMatrices(jac->nsplits, &jac->Afield);
1516:   jac->nsplits = 0;
1517:   VecDestroy(&jac->w1);
1518:   VecDestroy(&jac->w2);
1519:   MatDestroy(&jac->schur);
1520:   MatDestroy(&jac->schurp);
1521:   MatDestroy(&jac->schur_user);
1522:   KSPDestroy(&jac->kspschur);
1523:   KSPDestroy(&jac->kspupper);
1524:   MatDestroy(&jac->B);
1525:   MatDestroy(&jac->C);
1526:   MatDestroy(&jac->H);
1527:   VecDestroy(&jac->u);
1528:   VecDestroy(&jac->v);
1529:   VecDestroy(&jac->Hu);
1530:   VecDestroy(&jac->d);
1531:   PetscFree(jac->vecz);
1532:   PetscViewerDestroy(&jac->gkbviewer);
1533:   jac->isrestrict = PETSC_FALSE;
1534:   return 0;
1535: }

1537: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
1538: {
1539:   PCReset_FieldSplit(pc);
1540:   PetscFree(pc->data);
1541:   PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL);
1542:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetFields_C", NULL);
1543:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetIS_C", NULL);
1544:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetType_C", NULL);
1545:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetBlockSize_C", NULL);
1546:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitRestrictIS_C", NULL);
1547:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSchurGetSubKSP_C", NULL);
1548:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", NULL);

1550:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", NULL);
1551:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", NULL);
1552:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", NULL);
1553:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", NULL);
1554:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", NULL);
1555:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", NULL);
1556:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", NULL);
1557:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", NULL);
1558:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", NULL);
1559:   return 0;
1560: }

1562: static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc, PetscOptionItems *PetscOptionsObject)
1563: {
1564:   PetscInt        bs;
1565:   PetscBool       flg;
1566:   PC_FieldSplit  *jac = (PC_FieldSplit *)pc->data;
1567:   PCCompositeType ctype;

1569:   PetscOptionsHeadBegin(PetscOptionsObject, "FieldSplit options");
1570:   PetscOptionsBool("-pc_fieldsplit_dm_splits", "Whether to use DMCreateFieldDecomposition() for splits", "PCFieldSplitSetDMSplits", jac->dm_splits, &jac->dm_splits, NULL);
1571:   PetscOptionsInt("-pc_fieldsplit_block_size", "Blocksize that defines number of fields", "PCFieldSplitSetBlockSize", jac->bs, &bs, &flg);
1572:   if (flg) PCFieldSplitSetBlockSize(pc, bs);
1573:   jac->diag_use_amat = pc->useAmat;
1574:   PetscOptionsBool("-pc_fieldsplit_diag_use_amat", "Use Amat (not Pmat) to extract diagonal fieldsplit blocks", "PCFieldSplitSetDiagUseAmat", jac->diag_use_amat, &jac->diag_use_amat, NULL);
1575:   jac->offdiag_use_amat = pc->useAmat;
1576:   PetscOptionsBool("-pc_fieldsplit_off_diag_use_amat", "Use Amat (not Pmat) to extract off-diagonal fieldsplit blocks", "PCFieldSplitSetOffDiagUseAmat", jac->offdiag_use_amat, &jac->offdiag_use_amat, NULL);
1577:   PetscOptionsBool("-pc_fieldsplit_detect_saddle_point", "Form 2-way split by detecting zero diagonal entries", "PCFieldSplitSetDetectSaddlePoint", jac->detect, &jac->detect, NULL);
1578:   PCFieldSplitSetDetectSaddlePoint(pc, jac->detect); /* Sets split type and Schur PC type */
1579:   PetscOptionsEnum("-pc_fieldsplit_type", "Type of composition", "PCFieldSplitSetType", PCCompositeTypes, (PetscEnum)jac->type, (PetscEnum *)&ctype, &flg);
1580:   if (flg) PCFieldSplitSetType(pc, ctype);
1581:   /* Only setup fields once */
1582:   if ((jac->bs > 0) && (jac->nsplits == 0)) {
1583:     /* only allow user to set fields from command line if bs is already known.
1584:        otherwise user can set them in PCFieldSplitSetDefaults() */
1585:     PCFieldSplitSetRuntimeSplits_Private(pc);
1586:     if (jac->splitdefined) PetscInfo(pc, "Splits defined using the options database\n");
1587:   }
1588:   if (jac->type == PC_COMPOSITE_SCHUR) {
1589:     PetscOptionsGetEnum(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_fieldsplit_schur_factorization_type", PCFieldSplitSchurFactTypes, (PetscEnum *)&jac->schurfactorization, &flg);
1590:     if (flg) PetscInfo(pc, "Deprecated use of -pc_fieldsplit_schur_factorization_type\n");
1591:     PetscOptionsEnum("-pc_fieldsplit_schur_fact_type", "Which off-diagonal parts of the block factorization to use", "PCFieldSplitSetSchurFactType", PCFieldSplitSchurFactTypes, (PetscEnum)jac->schurfactorization, (PetscEnum *)&jac->schurfactorization, NULL);
1592:     PetscOptionsEnum("-pc_fieldsplit_schur_precondition", "How to build preconditioner for Schur complement", "PCFieldSplitSetSchurPre", PCFieldSplitSchurPreTypes, (PetscEnum)jac->schurpre, (PetscEnum *)&jac->schurpre, NULL);
1593:     PetscOptionsScalar("-pc_fieldsplit_schur_scale", "Scale Schur complement", "PCFieldSplitSetSchurScale", jac->schurscale, &jac->schurscale, NULL);
1594:   } else if (jac->type == PC_COMPOSITE_GKB) {
1595:     PetscOptionsReal("-pc_fieldsplit_gkb_tol", "The tolerance for the lower bound stopping criterion", "PCFieldSplitGKBTol", jac->gkbtol, &jac->gkbtol, NULL);
1596:     PetscOptionsInt("-pc_fieldsplit_gkb_delay", "The delay value for lower bound criterion", "PCFieldSplitGKBDelay", jac->gkbdelay, &jac->gkbdelay, NULL);
1597:     PetscOptionsReal("-pc_fieldsplit_gkb_nu", "Parameter in augmented Lagrangian approach", "PCFieldSplitGKBNu", jac->gkbnu, &jac->gkbnu, NULL);
1599:     PetscOptionsInt("-pc_fieldsplit_gkb_maxit", "Maximum allowed number of iterations", "PCFieldSplitGKBMaxit", jac->gkbmaxit, &jac->gkbmaxit, NULL);
1600:     PetscOptionsBool("-pc_fieldsplit_gkb_monitor", "Prints number of GKB iterations and error", "PCFieldSplitGKB", jac->gkbmonitor, &jac->gkbmonitor, NULL);
1601:   }
1602:   /*
1603:     In the initial call to this routine the sub-solver data structures do not exist so we cannot call KSPSetFromOptions() on them yet.
1604:     But after the initial setup of ALL the layers of sub-solvers is completed we do want to call KSPSetFromOptions() on the sub-solvers every time it
1605:     is called on the outer solver in case changes were made in the options database

1607:     But even after PCSetUp_FieldSplit() is called all the options inside the inner levels of sub-solvers may still not have been set thus we only call the KSPSetFromOptions()
1608:     if we know that the entire stack of sub-solvers below this have been complete instantiated, we check this by seeing if any solver iterations are complete.
1609:     Without this extra check test p2p1fetidp_olof_full and others fail with incorrect matrix types.

1611:     There could be a negative side effect of calling the KSPSetFromOptions() below.

1613:     If one captured the PetscObjectState of the options database one could skip these calls if the database has not changed from the previous call
1614:   */
1615:   if (jac->issetup) {
1616:     PC_FieldSplitLink ilink = jac->head;
1617:     if (jac->type == PC_COMPOSITE_SCHUR) {
1618:       if (jac->kspupper && jac->kspupper->totalits > 0) KSPSetFromOptions(jac->kspupper);
1619:       if (jac->kspschur && jac->kspschur->totalits > 0) KSPSetFromOptions(jac->kspschur);
1620:     }
1621:     while (ilink) {
1622:       if (ilink->ksp->totalits > 0) KSPSetFromOptions(ilink->ksp);
1623:       ilink = ilink->next;
1624:     }
1625:   }
1626:   PetscOptionsHeadEnd();
1627:   return 0;
1628: }

1630: static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc, const char splitname[], PetscInt n, const PetscInt *fields, const PetscInt *fields_col)
1631: {
1632:   PC_FieldSplit    *jac = (PC_FieldSplit *)pc->data;
1633:   PC_FieldSplitLink ilink, next = jac->head;
1634:   char              prefix[128];
1635:   PetscInt          i;

1637:   if (jac->splitdefined) {
1638:     PetscInfo(pc, "Ignoring new split \"%s\" because the splits have already been defined\n", splitname);
1639:     return 0;
1640:   }
1641:   for (i = 0; i < n; i++) {
1644:   }
1645:   PetscNew(&ilink);
1646:   if (splitname) {
1647:     PetscStrallocpy(splitname, &ilink->splitname);
1648:   } else {
1649:     PetscMalloc1(3, &ilink->splitname);
1650:     PetscSNPrintf(ilink->splitname, 2, "%" PetscInt_FMT, jac->nsplits);
1651:   }
1652:   ilink->event = jac->nsplits < 5 ? KSP_Solve_FS_0 + jac->nsplits : KSP_Solve_FS_0 + 4; /* Any split great than 4 gets logged in the 4th split */
1653:   PetscMalloc1(n, &ilink->fields);
1654:   PetscArraycpy(ilink->fields, fields, n);
1655:   PetscMalloc1(n, &ilink->fields_col);
1656:   PetscArraycpy(ilink->fields_col, fields_col, n);

1658:   ilink->nfields = n;
1659:   ilink->next    = NULL;
1660:   KSPCreate(PetscObjectComm((PetscObject)pc), &ilink->ksp);
1661:   KSPSetErrorIfNotConverged(ilink->ksp, pc->erroriffailure);
1662:   PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)pc, 1);
1663:   KSPSetType(ilink->ksp, KSPPREONLY);

1665:   PetscSNPrintf(prefix, sizeof(prefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
1666:   KSPSetOptionsPrefix(ilink->ksp, prefix);

1668:   if (!next) {
1669:     jac->head       = ilink;
1670:     ilink->previous = NULL;
1671:   } else {
1672:     while (next->next) next = next->next;
1673:     next->next      = ilink;
1674:     ilink->previous = next;
1675:   }
1676:   jac->nsplits++;
1677:   return 0;
1678: }

1680: static PetscErrorCode PCFieldSplitSchurGetSubKSP_FieldSplit(PC pc, PetscInt *n, KSP **subksp)
1681: {
1682:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

1684:   *subksp = NULL;
1685:   if (n) *n = 0;
1686:   if (jac->type == PC_COMPOSITE_SCHUR) {
1687:     PetscInt nn;

1691:     nn = jac->nsplits + (jac->kspupper != jac->head->ksp ? 1 : 0);
1692:     PetscMalloc1(nn, subksp);
1693:     (*subksp)[0] = jac->head->ksp;
1694:     (*subksp)[1] = jac->kspschur;
1695:     if (jac->kspupper != jac->head->ksp) (*subksp)[2] = jac->kspupper;
1696:     if (n) *n = nn;
1697:   }
1698:   return 0;
1699: }

1701: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc, PetscInt *n, KSP **subksp)
1702: {
1703:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

1706:   PetscMalloc1(jac->nsplits, subksp);
1707:   MatSchurComplementGetKSP(jac->schur, *subksp);

1709:   (*subksp)[1] = jac->kspschur;
1710:   if (n) *n = jac->nsplits;
1711:   return 0;
1712: }

1714: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc, PetscInt *n, KSP **subksp)
1715: {
1716:   PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
1717:   PetscInt          cnt   = 0;
1718:   PC_FieldSplitLink ilink = jac->head;

1720:   PetscMalloc1(jac->nsplits, subksp);
1721:   while (ilink) {
1722:     (*subksp)[cnt++] = ilink->ksp;
1723:     ilink            = ilink->next;
1724:   }
1726:   if (n) *n = jac->nsplits;
1727:   return 0;
1728: }

1730: /*@C
1731:     PCFieldSplitRestrictIS - Restricts the fieldsplit `IS`s to be within a given `IS`.

1733:     Input Parameters:
1734: +   pc  - the preconditioner context
1735: -   is - the index set that defines the indices to which the fieldsplit is to be restricted

1737:     Level: advanced

1739:     Developer Note:
1740:     It seems the resulting `IS`s will not cover the entire space, so
1741:     how can they define a convergent preconditioner? Needs explaining.

1743: .seealso: [](sec_block_matrices), `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`
1744: @*/
1745: PetscErrorCode PCFieldSplitRestrictIS(PC pc, IS isy)
1746: {
1749:   PetscTryMethod(pc, "PCFieldSplitRestrictIS_C", (PC, IS), (pc, isy));
1750:   return 0;
1751: }

1753: static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy)
1754: {
1755:   PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
1756:   PC_FieldSplitLink ilink = jac->head, next;
1757:   PetscInt          localsize, size, sizez, i;
1758:   const PetscInt   *ind, *indz;
1759:   PetscInt         *indc, *indcz;
1760:   PetscBool         flg;

1762:   ISGetLocalSize(isy, &localsize);
1763:   MPI_Scan(&localsize, &size, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)isy));
1764:   size -= localsize;
1765:   while (ilink) {
1766:     IS isrl, isr;
1767:     PC subpc;
1768:     ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl);
1769:     ISGetLocalSize(isrl, &localsize);
1770:     PetscMalloc1(localsize, &indc);
1771:     ISGetIndices(isrl, &ind);
1772:     PetscArraycpy(indc, ind, localsize);
1773:     ISRestoreIndices(isrl, &ind);
1774:     ISDestroy(&isrl);
1775:     for (i = 0; i < localsize; i++) *(indc + i) += size;
1776:     ISCreateGeneral(PetscObjectComm((PetscObject)isy), localsize, indc, PETSC_OWN_POINTER, &isr);
1777:     PetscObjectReference((PetscObject)isr);
1778:     ISDestroy(&ilink->is);
1779:     ilink->is = isr;
1780:     PetscObjectReference((PetscObject)isr);
1781:     ISDestroy(&ilink->is_col);
1782:     ilink->is_col = isr;
1783:     ISDestroy(&isr);
1784:     KSPGetPC(ilink->ksp, &subpc);
1785:     PetscObjectTypeCompare((PetscObject)subpc, PCFIELDSPLIT, &flg);
1786:     if (flg) {
1787:       IS       iszl, isz;
1788:       MPI_Comm comm;
1789:       ISGetLocalSize(ilink->is, &localsize);
1790:       comm = PetscObjectComm((PetscObject)ilink->is);
1791:       ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl);
1792:       MPI_Scan(&localsize, &sizez, 1, MPIU_INT, MPI_SUM, comm);
1793:       sizez -= localsize;
1794:       ISGetLocalSize(iszl, &localsize);
1795:       PetscMalloc1(localsize, &indcz);
1796:       ISGetIndices(iszl, &indz);
1797:       PetscArraycpy(indcz, indz, localsize);
1798:       ISRestoreIndices(iszl, &indz);
1799:       ISDestroy(&iszl);
1800:       for (i = 0; i < localsize; i++) *(indcz + i) += sizez;
1801:       ISCreateGeneral(comm, localsize, indcz, PETSC_OWN_POINTER, &isz);
1802:       PCFieldSplitRestrictIS(subpc, isz);
1803:       ISDestroy(&isz);
1804:     }
1805:     next  = ilink->next;
1806:     ilink = next;
1807:   }
1808:   jac->isrestrict = PETSC_TRUE;
1809:   return 0;
1810: }

1812: static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc, const char splitname[], IS is)
1813: {
1814:   PC_FieldSplit    *jac = (PC_FieldSplit *)pc->data;
1815:   PC_FieldSplitLink ilink, next = jac->head;
1816:   char              prefix[128];

1818:   if (jac->splitdefined) {
1819:     PetscInfo(pc, "Ignoring new split \"%s\" because the splits have already been defined\n", splitname);
1820:     return 0;
1821:   }
1822:   PetscNew(&ilink);
1823:   if (splitname) {
1824:     PetscStrallocpy(splitname, &ilink->splitname);
1825:   } else {
1826:     PetscMalloc1(8, &ilink->splitname);
1827:     PetscSNPrintf(ilink->splitname, 7, "%" PetscInt_FMT, jac->nsplits);
1828:   }
1829:   ilink->event = jac->nsplits < 5 ? KSP_Solve_FS_0 + jac->nsplits : KSP_Solve_FS_0 + 4; /* Any split great than 4 gets logged in the 4th split */
1830:   PetscObjectReference((PetscObject)is);
1831:   ISDestroy(&ilink->is);
1832:   ilink->is = is;
1833:   PetscObjectReference((PetscObject)is);
1834:   ISDestroy(&ilink->is_col);
1835:   ilink->is_col = is;
1836:   ilink->next   = NULL;
1837:   KSPCreate(PetscObjectComm((PetscObject)pc), &ilink->ksp);
1838:   KSPSetErrorIfNotConverged(ilink->ksp, pc->erroriffailure);
1839:   PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)pc, 1);
1840:   KSPSetType(ilink->ksp, KSPPREONLY);

1842:   PetscSNPrintf(prefix, sizeof(prefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
1843:   KSPSetOptionsPrefix(ilink->ksp, prefix);

1845:   if (!next) {
1846:     jac->head       = ilink;
1847:     ilink->previous = NULL;
1848:   } else {
1849:     while (next->next) next = next->next;
1850:     next->next      = ilink;
1851:     ilink->previous = next;
1852:   }
1853:   jac->nsplits++;
1854:   return 0;
1855: }

1857: /*@C
1858:     PCFieldSplitSetFields - Sets the fields that define one particular split in `PCFIELDSPLIT`

1860:     Logically Collective

1862:     Input Parameters:
1863: +   pc  - the preconditioner context
1864: .   splitname - name of this split, if `NULL` the number of the split is used
1865: .   n - the number of fields in this split
1866: .   fields - the fields in this split
1867: -   fields_col - generally the same as fields, if it does not match fields then the matrix block that is solved for this set of fields comes from an off-diagonal block
1868:                  of the matrix and fields_col provides the column indices for that block

1870:     Level: intermediate

1872:     Notes:
1873:     Use `PCFieldSplitSetIS()` to set a  general set of indices as a split.

1875:      `PCFieldSplitSetFields()` is for defining fields as strided blocks. For example, if the block
1876:      size is three then one can define a split as 0, or 1 or 2 or 0,1 or 0,2 or 1,2 which mean
1877:      0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x....
1878:      where the numbered entries indicate what is in the split.

1880:      This function is called once per split (it creates a new split each time).  Solve options
1881:      for this split will be available under the prefix `-fieldsplit_SPLITNAME_`.

1883:    `PCFieldSplitSetIS()` does not support having a fields_col different from fields

1885:    Developer Note:
1886:    This routine does not actually create the `IS` representing the split, that is delayed
1887:    until `PCSetUp_FieldSplit()`, because information about the vector/matrix layouts may not be
1888:    available when this routine is called.

1890: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetBlockSize()`, `PCFieldSplitSetIS()`, `PCFieldSplitRestrictIS()`
1891: @*/
1892: PetscErrorCode PCFieldSplitSetFields(PC pc, const char splitname[], PetscInt n, const PetscInt *fields, const PetscInt *fields_col)
1893: {
1898:   PetscTryMethod(pc, "PCFieldSplitSetFields_C", (PC, const char[], PetscInt, const PetscInt *, const PetscInt *), (pc, splitname, n, fields, fields_col));
1899:   return 0;
1900: }

1902: /*@
1903:     PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) to build
1904:     the sub-matrices associated with each split. Where `KSPSetOperators`(ksp,Amat,Pmat)) was used to supply the operators.

1906:     Logically Collective

1908:     Input Parameters:
1909: +   pc  - the preconditioner object
1910: -   flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from

1912:     Options Database Keys:
1913: .   -pc_fieldsplit_diag_use_amat - use the Amat to provide the diagonal blocks

1915:     Level: intermediate

1917: .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitGetDiagUseAmat()`, `PCFieldSplitSetOffDiagUseAmat()`, `PCFIELDSPLIT`
1918: @*/
1919: PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc, PetscBool flg)
1920: {
1921:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1922:   PetscBool      isfs;

1925:   PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs);
1927:   jac->diag_use_amat = flg;
1928:   return 0;
1929: }

1931: /*@
1932:     PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) to build
1933:     the sub-matrices associated with each split.  Where `KSPSetOperators`(ksp,Amat,Pmat)) was used to supply the operators.

1935:     Logically Collective

1937:     Input Parameters:
1938: .   pc  - the preconditioner object

1940:     Output Parameters:
1941: .   flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from

1943:     Level: intermediate

1945: .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitSetDiagUseAmat()`, `PCFieldSplitGetOffDiagUseAmat()`, `PCFIELDSPLIT`
1946: @*/
1947: PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc, PetscBool *flg)
1948: {
1949:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1950:   PetscBool      isfs;

1954:   PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs);
1956:   *flg = jac->diag_use_amat;
1957:   return 0;
1958: }

1960: /*@
1961:     PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) to build
1962:     the sub-matrices associated with each split.  Where `KSPSetOperators`(ksp,Amat,Pmat)) was used to supply the operators.

1964:     Logically Collective

1966:     Input Parameters:
1967: +   pc  - the preconditioner object
1968: -   flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from

1970:     Options Database Keys:
1971: .     -pc_fieldsplit_off_diag_use_amat <bool> - use the Amat to extract the off-diagonal blocks

1973:     Level: intermediate

1975: .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitGetOffDiagUseAmat()`, `PCFieldSplitSetDiagUseAmat()`, `PCFIELDSPLIT`
1976: @*/
1977: PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc, PetscBool flg)
1978: {
1979:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1980:   PetscBool      isfs;

1983:   PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs);
1985:   jac->offdiag_use_amat = flg;
1986:   return 0;
1987: }

1989: /*@
1990:     PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) to build
1991:     the sub-matrices associated with each split.  Where `KSPSetOperators`(ksp,Amat,Pmat)) was used to supply the operators.

1993:     Logically Collective

1995:     Input Parameters:
1996: .   pc  - the preconditioner object

1998:     Output Parameters:
1999: .   flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from

2001:     Level: intermediate

2003: .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitSetOffDiagUseAmat()`, `PCFieldSplitGetDiagUseAmat()`, `PCFIELDSPLIT`
2004: @*/
2005: PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc, PetscBool *flg)
2006: {
2007:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2008:   PetscBool      isfs;

2012:   PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs);
2014:   *flg = jac->offdiag_use_amat;
2015:   return 0;
2016: }

2018: /*@C
2019:     PCFieldSplitSetIS - Sets the exact elements for a split in a `PCFIELDSPLIT`

2021:     Logically Collective

2023:     Input Parameters:
2024: +   pc  - the preconditioner context
2025: .   splitname - name of this split, if `NULL` the number of the split is used
2026: -   is - the index set that defines the elements in this split

2028:     Level: intermediate

2030:     Notes:
2031:     Use `PCFieldSplitSetFields()`, for splits defined by strided types.

2033:     This function is called once per split (it creates a new split each time).  Solve options
2034:     for this split will be available under the prefix -fieldsplit_SPLITNAME_.

2036: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetBlockSize()`
2037: @*/
2038: PetscErrorCode PCFieldSplitSetIS(PC pc, const char splitname[], IS is)
2039: {
2043:   PetscTryMethod(pc, "PCFieldSplitSetIS_C", (PC, const char[], IS), (pc, splitname, is));
2044:   return 0;
2045: }

2047: /*@C
2048:     PCFieldSplitGetIS - Retrieves the elements for a split as an `IS`

2050:     Logically Collective

2052:     Input Parameters:
2053: +   pc  - the preconditioner context
2054: -   splitname - name of this split

2056:     Output Parameter:
2057: -   is - the index set that defines the elements in this split, or `NULL` if the split is not found

2059:     Level: intermediate

2061: .seealso:  [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetIS()`
2062: @*/
2063: PetscErrorCode PCFieldSplitGetIS(PC pc, const char splitname[], IS *is)
2064: {
2068:   {
2069:     PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
2070:     PC_FieldSplitLink ilink = jac->head;
2071:     PetscBool         found;

2073:     *is = NULL;
2074:     while (ilink) {
2075:       PetscStrcmp(ilink->splitname, splitname, &found);
2076:       if (found) {
2077:         *is = ilink->is;
2078:         break;
2079:       }
2080:       ilink = ilink->next;
2081:     }
2082:   }
2083:   return 0;
2084: }

2086: /*@C
2087:     PCFieldSplitGetISByIndex - Retrieves the elements for a given split as an `IS`

2089:     Logically Collective

2091:     Input Parameters:
2092: +   pc  - the preconditioner context
2093: -   index - index of this split

2095:     Output Parameter:
2096: -   is - the index set that defines the elements in this split

2098:     Level: intermediate

2100: .seealso:  [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitGetIS()`, `PCFieldSplitSetIS()`
2101: @*/
2102: PetscErrorCode PCFieldSplitGetISByIndex(PC pc, PetscInt index, IS *is)
2103: {
2107:   {
2108:     PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
2109:     PC_FieldSplitLink ilink = jac->head;
2110:     PetscInt          i     = 0;

2113:     while (i < index) {
2114:       ilink = ilink->next;
2115:       ++i;
2116:     }
2117:     PCFieldSplitGetIS(pc, ilink->splitname, is);
2118:   }
2119:   return 0;
2120: }

2122: /*@
2123:     PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the
2124:       fieldsplit preconditioner when calling `PCFieldSplitSetIS()`. If not set the matrix block size is used.

2126:     Logically Collective

2128:     Input Parameters:
2129: +   pc  - the preconditioner context
2130: -   bs - the block size

2132:     Level: intermediate

2134: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`
2135: @*/
2136: PetscErrorCode PCFieldSplitSetBlockSize(PC pc, PetscInt bs)
2137: {
2140:   PetscTryMethod(pc, "PCFieldSplitSetBlockSize_C", (PC, PetscInt), (pc, bs));
2141:   return 0;
2142: }

2144: /*@C
2145:    PCFieldSplitGetSubKSP - Gets the `KSP` contexts for all splits

2147:    Collective

2149:    Input Parameter:
2150: .  pc - the preconditioner context

2152:    Output Parameters:
2153: +  n - the number of splits
2154: -  subksp - the array of `KSP` contexts

2156:    Level: advanced

2158:    Notes:
2159:    After `PCFieldSplitGetSubKSP()` the array of `KSP`s is to be freed by the user with `PetscFree()`
2160:    (not the `KSP`, just the array that contains them).

2162:    You must call `PCSetUp()` before calling `PCFieldSplitGetSubKSP()`.

2164:    If the fieldsplit is of type `PC_COMPOSITE_SCHUR`, it returns the `KSP` object used inside the
2165:    Schur complement and the `KSP` object used to iterate over the Schur complement.
2166:    To access all the `KSP` objects used in `PC_COMPOSITE_SCHUR`, use `PCFieldSplitSchurGetSubKSP()`.

2168:    If the fieldsplit is of type `PC_COMPOSITE_GKB`, it returns the `KSP` object used to solve the
2169:    inner linear system defined by the matrix H in each loop.

2171:    Fortran Usage:
2172:    You must pass in a `KSP` array that is large enough to contain all the `KSP`s.
2173:    You can call `PCFieldSplitGetSubKSP`(pc,n,`PETSC_NULL_KSP`,ierr) to determine how large the
2174:    `KSP` array must be.

2176:    Developer Note:
2177:    There should be a `PCFieldSplitRestoreSubKSP()` instead of requiring the user to call `PetscFree()`

2179:    The Fortran interface should be modernized to return directly the array of values.

2181: .seealso:  [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`, `PCFieldSplitSchurGetSubKSP()`
2182: @*/
2183: PetscErrorCode PCFieldSplitGetSubKSP(PC pc, PetscInt *n, KSP *subksp[])
2184: {
2187:   PetscUseMethod(pc, "PCFieldSplitGetSubKSP_C", (PC, PetscInt *, KSP **), (pc, n, subksp));
2188:   return 0;
2189: }

2191: /*@C
2192:    PCFieldSplitSchurGetSubKSP - Gets the `KSP` contexts used inside the Schur complement based `PCFIELDSPLIT`

2194:    Collective

2196:    Input Parameter:
2197: .  pc - the preconditioner context

2199:    Output Parameters:
2200: +  n - the number of splits
2201: -  subksp - the array of `KSP` contexts

2203:    Level: advanced

2205:    Notes:
2206:    After `PCFieldSplitSchurGetSubKSP()` the array of `KSP`s is to be freed by the user with `PetscFree()`
2207:    (not the `KSP` just the array that contains them).

2209:    You must call `PCSetUp()` before calling `PCFieldSplitSchurGetSubKSP()`.

2211:    If the fieldsplit type is of type `PC_COMPOSITE_SCHUR`, it returns (in order)
2212: +  1  - the `KSP` used for the (1,1) block
2213: .  2  - the `KSP` used for the Schur complement (not the one used for the interior Schur solver)
2214: -  3  - the `KSP` used for the (1,1) block in the upper triangular factor (if different from that of the (1,1) block).

2216:    It returns a null array if the fieldsplit is not of type `PC_COMPOSITE_SCHUR`; in this case, you should use `PCFieldSplitGetSubKSP()`.

2218:    Fortran Note:
2219:    You must pass in a `KSP` array that is large enough to contain all the local `KSP`s.
2220:    You can call `PCFieldSplitSchurGetSubKSP`(pc,n,`PETSC_NULL_KSP`,ierr) to determine how large the
2221:    `KSP` array must be.

2223:    Developer Notes:
2224:    There should be a `PCFieldSplitRestoreSubKSP()` instead of requiring the user to call `PetscFree()`

2226:    Should the functionality of `PCFieldSplitSchurGetSubKSP()` and `PCFieldSplitGetSubKSP()` be merged?

2228:    The Fortran interface should be modernized to return directly the array of values.

2230: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`, `PCFieldSplitGetSubKSP()`
2231: @*/
2232: PetscErrorCode PCFieldSplitSchurGetSubKSP(PC pc, PetscInt *n, KSP *subksp[])
2233: {
2236:   PetscUseMethod(pc, "PCFieldSplitSchurGetSubKSP_C", (PC, PetscInt *, KSP **), (pc, n, subksp));
2237:   return 0;
2238: }

2240: /*@
2241:     PCFieldSplitSetSchurPre -  Indicates from what operator the preconditioner is constructucted for the Schur complement.
2242:       The default is the A11 matrix.

2244:     Collective

2246:     Input Parameters:
2247: +   pc      - the preconditioner context
2248: .   ptype   - which matrix to use for preconditioning the Schur complement: `PC_FIELDSPLIT_SCHUR_PRE_A11` (default),
2249:               `PC_FIELDSPLIT_SCHUR_PRE_SELF`, `PC_FIELDSPLIT_SCHUR_PRE_USER`,
2250:               `PC_FIELDSPLIT_SCHUR_PRE_SELFP`, and `PC_FIELDSPLIT_SCHUR_PRE_FULL`
2251: -   pre - matrix to use for preconditioning, or `NULL`

2253:     Options Database Keys:
2254: +    -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11. See notes for meaning of various arguments
2255: -    -fieldsplit_1_pc_type <pctype> - the preconditioner algorithm that is used to construct the preconditioner from the operator

2257:      Level: intermediate

2259:     Notes:
2260:     If ptype is
2261: +     a11 - the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner
2262:      matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix
2263: .     self - the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix:
2264:           The only preconditioner that currently works with this symbolic representation matrix object is the `PCLSC`
2265:           preconditioner
2266: .     user - the preconditioner for the Schur complement is generated from the user provided matrix (pre argument
2267:           to this function).
2268: .     selfp - the preconditioning for the Schur complement is generated from an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01
2269:           This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be
2270:           lumped before extracting the diagonal using the additional option -fieldsplit_1_mat_schur_complement_ainv_type lump
2271: -     full - the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation
2272:       computed internally by `PCFIELDSPLIT` (this is expensive)
2273:       useful mostly as a test that the Schur complement approach can work for your problem

2275:      When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense
2276:     with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and
2277:     -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement.

2279: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSchurPre()`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`,
2280:           `MatSchurComplementSetAinvType()`, `PCLSC`,
2281:           `PCFieldSplitSchurPreType`
2282: @*/
2283: PetscErrorCode PCFieldSplitSetSchurPre(PC pc, PCFieldSplitSchurPreType ptype, Mat pre)
2284: {
2286:   PetscTryMethod(pc, "PCFieldSplitSetSchurPre_C", (PC, PCFieldSplitSchurPreType, Mat), (pc, ptype, pre));
2287:   return 0;
2288: }

2290: PetscErrorCode PCFieldSplitSchurPrecondition(PC pc, PCFieldSplitSchurPreType ptype, Mat pre)
2291: {
2292:   return PCFieldSplitSetSchurPre(pc, ptype, pre);
2293: } /* Deprecated name */

2295: /*@
2296:     PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be
2297:     preconditioned.  See `PCFieldSplitSetSchurPre()` for details.

2299:     Logically Collective

2301:     Input Parameter:
2302: .   pc      - the preconditioner context

2304:     Output Parameters:
2305: +   ptype   - which matrix to use for preconditioning the Schur complement: `PC_FIELDSPLIT_SCHUR_PRE_A11`, `PC_FIELDSPLIT_SCHUR_PRE_SELF`, `PC_FIELDSPLIT_PRE_USER`
2306: -   pre - matrix to use for preconditioning (with `PC_FIELDSPLIT_PRE_USER`), or NULL

2308:     Level: intermediate

2310: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitSetSchurPre()`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, `PCLSC`,
2311:           `PCFieldSplitSchurPreType`
2312: @*/
2313: PetscErrorCode PCFieldSplitGetSchurPre(PC pc, PCFieldSplitSchurPreType *ptype, Mat *pre)
2314: {
2316:   PetscUseMethod(pc, "PCFieldSplitGetSchurPre_C", (PC, PCFieldSplitSchurPreType *, Mat *), (pc, ptype, pre));
2317:   return 0;
2318: }

2320: /*@
2321:     PCFieldSplitSchurGetS -  extract the `MATSCHURCOMPLEMENT` object used by this `PCFIELDSPLIT` in case it needs to be configured separately

2323:     Not collective

2325:     Input Parameter:
2326: .   pc      - the preconditioner context

2328:     Output Parameter:
2329: .   S       - the Schur complement matrix

2331:     Level: advanced

2333:     Note:
2334:     This matrix should not be destroyed using `MatDestroy()`; rather, use `PCFieldSplitSchurRestoreS()`.

2336: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurPre()`, `MATSCHURCOMPLEMENT`, `PCFieldSplitSchurRestoreS()`,
2337:           `MatCreateSchurComplement()`, `MatSchurComplementGetKSP()`, `MatSchurComplementComputeExplicitOperator()`, `MatGetSchurComplement()`
2338: @*/
2339: PetscErrorCode PCFieldSplitSchurGetS(PC pc, Mat *S)
2340: {
2341:   const char    *t;
2342:   PetscBool      isfs;
2343:   PC_FieldSplit *jac;

2346:   PetscObjectGetType((PetscObject)pc, &t);
2347:   PetscStrcmp(t, PCFIELDSPLIT, &isfs);
2349:   jac = (PC_FieldSplit *)pc->data;
2351:   if (S) *S = jac->schur;
2352:   return 0;
2353: }

2355: /*@
2356:     PCFieldSplitSchurRestoreS -  returns the `MATSCHURCOMPLEMENT` matrix used by this `PC`

2358:     Not collective

2360:     Input Parameters:
2361: +   pc      - the preconditioner context
2362: -   S       - the Schur complement matrix

2364:     Level: advanced

2366: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurPre()`, `MatSchurComplement`, `PCFieldSplitSchurGetS()`
2367: @*/
2368: PetscErrorCode PCFieldSplitSchurRestoreS(PC pc, Mat *S)
2369: {
2370:   const char    *t;
2371:   PetscBool      isfs;
2372:   PC_FieldSplit *jac;

2375:   PetscObjectGetType((PetscObject)pc, &t);
2376:   PetscStrcmp(t, PCFIELDSPLIT, &isfs);
2378:   jac = (PC_FieldSplit *)pc->data;
2381:   return 0;
2382: }

2384: static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc, PCFieldSplitSchurPreType ptype, Mat pre)
2385: {
2386:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

2388:   jac->schurpre = ptype;
2389:   if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
2390:     MatDestroy(&jac->schur_user);
2391:     jac->schur_user = pre;
2392:     PetscObjectReference((PetscObject)jac->schur_user);
2393:   }
2394:   return 0;
2395: }

2397: static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc, PCFieldSplitSchurPreType *ptype, Mat *pre)
2398: {
2399:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

2401:   *ptype = jac->schurpre;
2402:   *pre   = jac->schur_user;
2403:   return 0;
2404: }

2406: /*@
2407:     PCFieldSplitSetSchurFactType -  sets which blocks of the approximate block factorization to retain in the preconditioner

2409:     Collective

2411:     Input Parameters:
2412: +   pc  - the preconditioner context
2413: -   ftype - which blocks of factorization to retain, `PC_FIELDSPLIT_SCHUR_FACT_FULL` is default

2415:     Options Database Key:
2416: .     -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> - default is full

2418:     Level: intermediate

2420:     Notes:
2421:     The FULL factorization is

2423: .vb
2424:    (A   B)  = (1       0) (A   0) (1  Ainv*B)  = L D U
2425:    (C   E)    (C*Ainv  1) (0   S) (0     1)
2426: .vb
2427:     where S = E - C*Ainv*B. In practice, the full factorization is applied via block triangular solves with the grouping L*(D*U). UPPER uses D*U, LOWER uses L*D,
2428:     and DIAG is the diagonal part with the sign of S flipped (because this makes the preconditioner positive definite for many formulations, thus allowing the use of KSPMINRES). Sign flipping of S can be turned off with PCFieldSplitSetSchurScale().

2430:     If A and S are solved exactly
2431: .vb
2432:       *) FULL factorization is a direct solver.
2433:       *) The preconditioned operator with LOWER or UPPER has all eigenvalues equal to 1 and minimal polynomial of degree 2, so `KSPGMRES` converges in 2 iterations.
2434:       *) With DIAG, the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so `KSPGMRES` converges in at most 4 iterations.
2435: .ve

2437:     If the iteration count is very low, consider using `KSPFGMRES` or `KSPGCR` which can use one less preconditioner
2438:     application in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice.

2440:     For symmetric problems in which A is positive definite and S is negative definite, DIAG can be used with `KSPMINRES`.

2442:     A flexible method like `KSPFGMRES` or `KSPGCR` must be used if the fieldsplit preconditioner is nonlinear (e.g. a few iterations of a Krylov method is used to solve with A or S).

2444:     References:
2445: +   * - Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000).
2446: -   * - Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001).

2448: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurScale()`
2449: @*/
2450: PetscErrorCode PCFieldSplitSetSchurFactType(PC pc, PCFieldSplitSchurFactType ftype)
2451: {
2453:   PetscTryMethod(pc, "PCFieldSplitSetSchurFactType_C", (PC, PCFieldSplitSchurFactType), (pc, ftype));
2454:   return 0;
2455: }

2457: static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc, PCFieldSplitSchurFactType ftype)
2458: {
2459:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

2461:   jac->schurfactorization = ftype;
2462:   return 0;
2463: }

2465: /*@
2466:     PCFieldSplitSetSchurScale -  Controls the sign flip of S for `PC_FIELDSPLIT_SCHUR_FACT_DIAG`.

2468:     Collective

2470:     Input Parameters:
2471: +   pc    - the preconditioner context
2472: -   scale - scaling factor for the Schur complement

2474:     Options Database Key:
2475: .   -pc_fieldsplit_schur_scale - default is -1.0

2477:     Level: intermediate

2479: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurFactType`, `PCFieldSplitSetSchurScale()`, `PCFieldSplitSetSchurFactType()`
2480: @*/
2481: PetscErrorCode PCFieldSplitSetSchurScale(PC pc, PetscScalar scale)
2482: {
2485:   PetscTryMethod(pc, "PCFieldSplitSetSchurScale_C", (PC, PetscScalar), (pc, scale));
2486:   return 0;
2487: }

2489: static PetscErrorCode PCFieldSplitSetSchurScale_FieldSplit(PC pc, PetscScalar scale)
2490: {
2491:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

2493:   jac->schurscale = scale;
2494:   return 0;
2495: }

2497: /*@C
2498:    PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement

2500:    Collective

2502:    Input Parameter:
2503: .  pc - the preconditioner context

2505:    Output Parameters:
2506: +  A00 - the (0,0) block
2507: .  A01 - the (0,1) block
2508: .  A10 - the (1,0) block
2509: -  A11 - the (1,1) block

2511:    Level: advanced

2513: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `MatSchurComplementGetSubMatrices()`, `MatSchurComplementSetSubMatrices()`
2514: @*/
2515: PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc, Mat *A00, Mat *A01, Mat *A10, Mat *A11)
2516: {
2517:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

2521:   if (A00) *A00 = jac->pmat[0];
2522:   if (A01) *A01 = jac->B;
2523:   if (A10) *A10 = jac->C;
2524:   if (A11) *A11 = jac->pmat[1];
2525:   return 0;
2526: }

2528: /*@
2529:     PCFieldSplitSetGKBTol -  Sets the solver tolerance for the generalized Golub-Kahan bidiagonalization preconditioner in `PCFIELDSPLIT`

2531:     Collective

2533:     Input Parameters:
2534: +   pc        - the preconditioner context
2535: -   tolerance - the solver tolerance

2537:     Options Database Key:
2538: .   -pc_fieldsplit_gkb_tol - default is 1e-5

2540:     Level: intermediate

2542:     Note:
2543:     The generalized GKB algorithm uses a lower bound estimate of the error in energy norm as stopping criterion.
2544:     It stops once the lower bound estimate undershoots the required solver tolerance. Although the actual error might be bigger than
2545:     this estimate, the stopping criterion is satisfactory in practical cases [A13].

2547:     References:
2548:     [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013.

2550: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBNu()`, `PCFieldSplitSetGKBMaxit()`
2551: @*/
2552: PetscErrorCode PCFieldSplitSetGKBTol(PC pc, PetscReal tolerance)
2553: {
2556:   PetscTryMethod(pc, "PCFieldSplitSetGKBTol_C", (PC, PetscReal), (pc, tolerance));
2557:   return 0;
2558: }

2560: static PetscErrorCode PCFieldSplitSetGKBTol_FieldSplit(PC pc, PetscReal tolerance)
2561: {
2562:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

2564:   jac->gkbtol = tolerance;
2565:   return 0;
2566: }

2568: /*@
2569:     PCFieldSplitSetGKBMaxit -  Sets the maximum number of iterations for the generalized Golub-Kahan bidiagonalization preconditioner in `PCFIELDSPLIT`

2571:     Collective

2573:     Input Parameters:
2574: +   pc     - the preconditioner context
2575: -   maxit  - the maximum number of iterations

2577:     Options Database Key:
2578: .   -pc_fieldsplit_gkb_maxit - default is 100

2580:     Level: intermediate

2582: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBNu()`
2583: @*/
2584: PetscErrorCode PCFieldSplitSetGKBMaxit(PC pc, PetscInt maxit)
2585: {
2588:   PetscTryMethod(pc, "PCFieldSplitSetGKBMaxit_C", (PC, PetscInt), (pc, maxit));
2589:   return 0;
2590: }

2592: static PetscErrorCode PCFieldSplitSetGKBMaxit_FieldSplit(PC pc, PetscInt maxit)
2593: {
2594:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

2596:   jac->gkbmaxit = maxit;
2597:   return 0;
2598: }

2600: /*@
2601:     PCFieldSplitSetGKBDelay -  Sets the delay in the lower bound error estimate in the generalized Golub-Kahan bidiagonalization in `PCFIELDSPLIT`
2602:     preconditioner.

2604:     Collective

2606:     Input Parameters:
2607: +   pc     - the preconditioner context
2608: -   delay  - the delay window in the lower bound estimate

2610:     Options Database Key:
2611: .   -pc_fieldsplit_gkb_delay - default is 5

2613:     Level: intermediate

2615:     Note:
2616:     The algorithm uses a lower bound estimate of the error in energy norm as stopping criterion. The lower bound of the error ||u-u^k||_H
2617:     is expressed as a truncated sum. The error at iteration k can only be measured at iteration (k + delay), and thus the algorithm needs
2618:     at least (delay + 1) iterations to stop. For more details on the generalized Golub-Kahan bidiagonalization method and its lower bound stopping criterion, please refer to

2620:     References:
2621:     [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013.

2623: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBNu()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBMaxit()`
2624: @*/
2625: PetscErrorCode PCFieldSplitSetGKBDelay(PC pc, PetscInt delay)
2626: {
2629:   PetscTryMethod(pc, "PCFieldSplitSetGKBDelay_C", (PC, PetscInt), (pc, delay));
2630:   return 0;
2631: }

2633: static PetscErrorCode PCFieldSplitSetGKBDelay_FieldSplit(PC pc, PetscInt delay)
2634: {
2635:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

2637:   jac->gkbdelay = delay;
2638:   return 0;
2639: }

2641: /*@
2642:     PCFieldSplitSetGKBNu -  Sets the scalar value nu >= 0 in the transformation H = A00 + nu*A01*A01' of the (1,1) block in the Golub-Kahan bidiagonalization preconditioner
2643:     in `PCFIELDSPLIT`

2645:     Collective

2647:     Input Parameters:
2648: +   pc     - the preconditioner context
2649: -   nu     - the shift parameter

2651:     Options Database Keys:
2652: .   -pc_fieldsplit_gkb_nu - default is 1

2654:     Level: intermediate

2656:     Notes:
2657:     This shift is in general done to obtain better convergence properties for the outer loop of the algorithm. This is often achieved by choosing nu sufficiently big. However,
2658:     if nu is chosen too big, the matrix H might be badly conditioned and the solution of the linear system Hx = b in the inner loop becomes difficult. It is therefore
2659:     necessary to find a good balance in between the convergence of the inner and outer loop.

2661:     For nu = 0, no shift is done. In this case A00 has to be positive definite. The matrix N in [Ar13] is then chosen as identity.

2663:     References:
2664:     [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013.

2666: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBMaxit()`
2667: @*/
2668: PetscErrorCode PCFieldSplitSetGKBNu(PC pc, PetscReal nu)
2669: {
2672:   PetscTryMethod(pc, "PCFieldSplitSetGKBNu_C", (PC, PetscReal), (pc, nu));
2673:   return 0;
2674: }

2676: static PetscErrorCode PCFieldSplitSetGKBNu_FieldSplit(PC pc, PetscReal nu)
2677: {
2678:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

2680:   jac->gkbnu = nu;
2681:   return 0;
2682: }

2684: static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc, PCCompositeType type)
2685: {
2686:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

2688:   jac->type = type;
2689:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", NULL);
2690:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", NULL);
2691:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", NULL);
2692:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", NULL);
2693:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", NULL);
2694:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", NULL);
2695:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", NULL);
2696:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", NULL);
2697:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", NULL);

2699:   if (type == PC_COMPOSITE_SCHUR) {
2700:     pc->ops->apply = PCApply_FieldSplit_Schur;
2701:     pc->ops->view  = PCView_FieldSplit_Schur;

2703:     PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit_Schur);
2704:     PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", PCFieldSplitSetSchurPre_FieldSplit);
2705:     PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", PCFieldSplitGetSchurPre_FieldSplit);
2706:     PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", PCFieldSplitSetSchurFactType_FieldSplit);
2707:     PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", PCFieldSplitSetSchurScale_FieldSplit);
2708:   } else if (type == PC_COMPOSITE_GKB) {
2709:     pc->ops->apply = PCApply_FieldSplit_GKB;
2710:     pc->ops->view  = PCView_FieldSplit_GKB;

2712:     PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit);
2713:     PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", PCFieldSplitSetGKBTol_FieldSplit);
2714:     PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", PCFieldSplitSetGKBMaxit_FieldSplit);
2715:     PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", PCFieldSplitSetGKBNu_FieldSplit);
2716:     PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", PCFieldSplitSetGKBDelay_FieldSplit);
2717:   } else {
2718:     pc->ops->apply = PCApply_FieldSplit;
2719:     pc->ops->view  = PCView_FieldSplit;

2721:     PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit);
2722:   }
2723:   return 0;
2724: }

2726: static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc, PetscInt bs)
2727: {
2728:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

2732:   jac->bs = bs;
2733:   return 0;
2734: }

2736: static PetscErrorCode PCSetCoordinates_FieldSplit(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[])
2737: {
2738:   PC_FieldSplit    *jac           = (PC_FieldSplit *)pc->data;
2739:   PC_FieldSplitLink ilink_current = jac->head;
2740:   IS                is_owned;

2742:   jac->coordinates_set = PETSC_TRUE; // Internal flag
2743:   MatGetOwnershipIS(pc->mat, &is_owned, PETSC_NULL);

2745:   while (ilink_current) {
2746:     // For each IS, embed it to get local coords indces
2747:     IS              is_coords;
2748:     PetscInt        ndofs_block;
2749:     const PetscInt *block_dofs_enumeration; // Numbering of the dofs relevant to the current block

2751:     // Setting drop to true for safety. It should make no difference.
2752:     ISEmbed(ilink_current->is, is_owned, PETSC_TRUE, &is_coords);
2753:     ISGetLocalSize(is_coords, &ndofs_block);
2754:     ISGetIndices(is_coords, &block_dofs_enumeration);

2756:     // Allocate coordinates vector and set it directly
2757:     PetscMalloc1(ndofs_block * dim, &(ilink_current->coords));
2758:     for (PetscInt dof = 0; dof < ndofs_block; ++dof) {
2759:       for (PetscInt d = 0; d < dim; ++d) (ilink_current->coords)[dim * dof + d] = coords[dim * block_dofs_enumeration[dof] + d];
2760:     }
2761:     ilink_current->dim   = dim;
2762:     ilink_current->ndofs = ndofs_block;
2763:     ISRestoreIndices(is_coords, &block_dofs_enumeration);
2764:     ISDestroy(&is_coords);
2765:     ilink_current = ilink_current->next;
2766:   }
2767:   ISDestroy(&is_owned);
2768:   return 0;
2769: }

2771: /*@
2772:    PCFieldSplitSetType - Sets the type, `PCCompositeType`, of a `PCFIELDSPLIT`

2774:    Collective

2776:    Input Parameters:
2777: +  pc - the preconditioner context
2778: -  type - `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE` (default), `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`

2780:    Options Database Key:
2781: .  -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type

2783:    Level: Intermediate

2785: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCCompositeType`, `PCCompositeGetType()`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`,
2786:           `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`
2787: @*/
2788: PetscErrorCode PCFieldSplitSetType(PC pc, PCCompositeType type)
2789: {
2791:   PetscTryMethod(pc, "PCFieldSplitSetType_C", (PC, PCCompositeType), (pc, type));
2792:   return 0;
2793: }

2795: /*@
2796:   PCFieldSplitGetType - Gets the type, `PCCompositeType`, of a `PCFIELDSPLIT`

2798:   Not collective

2800:   Input Parameter:
2801: . pc - the preconditioner context

2803:   Output Parameter:
2804: . type - `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE` (default), `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`

2806:   Level: Intermediate

2808: .seealso: [](sec_block_matrices), `PC`, `PCCompositeSetType()`, `PCFIELDSPLIT`, `PCCompositeType`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`,
2809:           `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`
2810: @*/
2811: PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
2812: {
2813:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

2817:   *type = jac->type;
2818:   return 0;
2819: }

2821: /*@
2822:    PCFieldSplitSetDMSplits - Flags whether `DMCreateFieldDecomposition()` should be used to define the splits in a `PCFIELDSPLIT`, whenever possible.

2824:    Logically Collective

2826:    Input Parameters:
2827: +  pc   - the preconditioner context
2828: -  flg  - boolean indicating whether to use field splits defined by the `DM`

2830:    Options Database Key:
2831: .  -pc_fieldsplit_dm_splits <bool> - use the field splits defined by the `DM`

2833:    Level: Intermediate

2835: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitGetDMSplits()`, `PCFieldSplitSetFields()`, `PCFieldsplitSetIS()`
2836: @*/
2837: PetscErrorCode PCFieldSplitSetDMSplits(PC pc, PetscBool flg)
2838: {
2839:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2840:   PetscBool      isfs;

2844:   PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs);
2845:   if (isfs) jac->dm_splits = flg;
2846:   return 0;
2847: }

2849: /*@
2850:    PCFieldSplitGetDMSplits - Returns flag indicating whether `DMCreateFieldDecomposition()` should be used to define the splits in a `PCFIELDSPLIT`, whenever possible.

2852:    Logically Collective

2854:    Input Parameter:
2855: .  pc   - the preconditioner context

2857:    Output Parameter:
2858: .  flg  - boolean indicating whether to use field splits defined by the `DM`

2860:    Level: Intermediate

2862: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetDMSplits()`, `PCFieldSplitSetFields()`, `PCFieldsplitSetIS()`
2863: @*/
2864: PetscErrorCode PCFieldSplitGetDMSplits(PC pc, PetscBool *flg)
2865: {
2866:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2867:   PetscBool      isfs;

2871:   PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs);
2872:   if (isfs) {
2873:     if (flg) *flg = jac->dm_splits;
2874:   }
2875:   return 0;
2876: }

2878: /*@
2879:    PCFieldSplitGetDetectSaddlePoint - Returns flag indicating whether `PCFIELDSPLIT` will attempt to automatically determine fields based on zero diagonal entries.

2881:    Logically Collective

2883:    Input Parameter:
2884: .  pc   - the preconditioner context

2886:    Output Parameter:
2887: .  flg  - boolean indicating whether to detect fields or not

2889:    Level: Intermediate

2891: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetDetectSaddlePoint()`
2892: @*/
2893: PetscErrorCode PCFieldSplitGetDetectSaddlePoint(PC pc, PetscBool *flg)
2894: {
2895:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

2897:   *flg = jac->detect;
2898:   return 0;
2899: }

2901: /*@
2902:    PCFieldSplitSetDetectSaddlePoint - Sets flag indicating whether `PCFIELDSPLIT` will attempt to automatically determine fields based on zero diagonal entries.

2904:    Logically Collective

2906:    Input Parameter:
2907: .  pc   - the preconditioner context

2909:    Output Parameter:
2910: .  flg  - boolean indicating whether to detect fields or not

2912:    Options Database Key:
2913: .  -pc_fieldsplit_detect_saddle_point <bool> - detect and use the saddle point

2915:    Level: Intermediate

2917:  Note:
2918:    Also sets the split type to `PC_COMPOSITE_SCHUR` (see `PCFieldSplitSetType()`) and the Schur preconditioner type to `PC_FIELDSPLIT_SCHUR_PRE_SELF` (see `PCFieldSplitSetSchurPre()`).

2920: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitGetDetectSaddlePoint()`, `PCFieldSplitSetType()`, `PCFieldSplitSetSchurPre()`, `PC_FIELDSPLIT_SCHUR_PRE_SELF`
2921: @*/
2922: PetscErrorCode PCFieldSplitSetDetectSaddlePoint(PC pc, PetscBool flg)
2923: {
2924:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

2926:   jac->detect = flg;
2927:   if (jac->detect) {
2928:     PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
2929:     PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_SELF, NULL);
2930:   }
2931:   return 0;
2932: }

2934: /*MC
2935:    PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
2936:    collections of variables (that may overlap) called splits. See [the users manual section on "Solving Block Matrices"](sec_block_matrices) for more details.

2938:    Options Database Keys:
2939: +   -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the `%d`'th split
2940: .   -pc_fieldsplit_default - automatically add any fields to additional splits that have not
2941:                               been supplied explicitly by `-pc_fieldsplit_%d_fields`
2942: .   -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields)
2943: .   -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur,gkb> - type of relaxation or factorization splitting
2944: .   -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11; see `PCFieldSplitSetSchurPre()`
2945: .   -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> - set factorization type when using `-pc_fieldsplit_type schur`; see `PCFieldSplitSetSchurFactType()`
2946: -   -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero diagonal and uses Schur complement with no preconditioner as the solver

2948:      Options prefixes for inner solvers when using the Schur complement preconditioner are `-fieldsplit_0_` and `-fieldsplit_1_` .
2949:      The options prefix for the inner solver when using the Golub-Kahan biadiagonalization preconditioner is `-fieldsplit_0_`
2950:      For all other solvers they are `-fieldsplit_%d_` for the `%d`'th field; use `-fieldsplit_` for all fields.

2952:      To set options on the solvers for each block append `-fieldsplit_` to all the `PC`
2953:      options database keys. For example, `-fieldsplit_pc_type ilu` `-fieldsplit_pc_factor_levels 1`

2955:      To set the options on the solvers separate for each block call `PCFieldSplitGetSubKSP()`
2956:       and set the options directly on the resulting `KSP` object

2958:     Level: intermediate

2960:    Notes:
2961:     Use `PCFieldSplitSetFields()` to set splits defined by "strided" entries and `PCFieldSplitSetIS()`
2962:      to define a split by an arbitrary collection of entries.

2964:       If no splits are set the default is used. The splits are defined by entries strided by bs,
2965:       beginning at 0 then 1, etc to bs-1. The block size can be set with `PCFieldSplitSetBlockSize()`,
2966:       if this is not called the block size defaults to the blocksize of the second matrix passed
2967:       to `KSPSetOperators()`/`PCSetOperators()`.

2969:       For the Schur complement preconditioner if

2971:       ```{math}
2972:       J = \left[\begin{array}{cc} A_{00} & A_{01} \\ A_{10} & A_{11} \end{array}\right]
2973:       ```

2975:       the preconditioner using `full` factorization is logically
2976:       ```{math}
2977:       \left[\begin{array}{cc} I & -\text{ksp}(A_{00}) \\ 0 & I \end{array}\right] \left[\begin{array}{cc} \text{inv}(A_{00}) & 0 \\ 0 & \text{ksp}(S) \end{array}\right] \left[\begin{array}{cc} I & 0 \\ -A_{10} \text{ksp}(A_{00}) & I \end{array}\right]
2978:       ```
2979:      where the action of $\text{inv}(A_{00})$ is applied using the KSP solver with prefix `-fieldsplit_0_`.  $S$ is the Schur complement
2980:      ```{math}
2981:      S = A_{11} - A_{10} \text{ksp}(A_{00}) A_{01}
2982:      ```
2983:      which is usually dense and not stored explicitly.  The action of $\text{ksp}(S)$ is computed using the KSP solver with prefix `-fieldsplit_splitname_` (where `splitname` was given
2984:      in providing the SECOND split or 1 if not given). For `PCFieldSplitGetSubKSP()` when field number is 0,
2985:      it returns the KSP associated with `-fieldsplit_0_` while field number 1 gives `-fieldsplit_1_` KSP. By default
2986:      $A_{11}$ is used to construct a preconditioner for $S$, use `PCFieldSplitSetSchurPre()` for all the possible ways to construct the preconditioner for $S$.

2988:      The factorization type is set using `-pc_fieldsplit_schur_fact_type <diag, lower, upper, full>`. `full` is shown above,
2989:      `diag` gives
2990:       ```{math}
2991:       \left[\begin{array}{cc} \text{inv}(A_{00}) & 0 \\  0 & -\text{ksp}(S) \end{array}\right]
2992:       ```
2993:      Note that, slightly counter intuitively, there is a negative in front of the $\text{ksp}(S)$  so that the preconditioner is positive definite. For SPD matrices $J$, the sign flip
2994:      can be turned off with `PCFieldSplitSetSchurScale()` or by command line `-pc_fieldsplit_schur_scale 1.0`. The `lower` factorization is the inverse of
2995:       ```{math}
2996:       \left[\begin{array}{cc} A_{00} & 0 \\  A_{10} & S \end{array}\right]
2997:       ```
2998:      where the inverses of A_{00} and S are applied using KSPs. The upper factorization is the inverse of
2999:       ```{math}
3000:       \left[\begin{array}{cc} A_{00} & A_{01} \\  0 & S \end{array}\right]
3001:       ```
3002:      where again the inverses of $A_{00}$ and $S$ are applied using `KSP`s.

3004:      If only one set of indices (one `IS`) is provided with `PCFieldSplitSetIS()` then the complement of that `IS`
3005:      is used automatically for a second block.

3007:      The fieldsplit preconditioner cannot currently be used with the `MATBAIJ` or `MATSBAIJ` data formats if the blocksize is larger than 1.
3008:      Generally it should be used with the `MATAIJ` format.

3010:      The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see,
3011:      for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling {cite}`Wesseling2009`.
3012:      One can also use `PCFIELDSPLIT`
3013:      inside a smoother resulting in "Distributive Smoothers".

3015:      See "A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations" {cite}`elman2008tcp`.

3017:      The Constrained Pressure Preconditioner (CPR) can be implemented using `PCCOMPOSITE` with `PCGALERKIN`. CPR first solves an $R A P$ subsystem, updates the
3018:      residual on all variables (`PCCompositeSetType(pc,PC_COMPOSITE_MULTIPLICATIVE)`), and then applies a simple ILU like preconditioner on all the variables.

3020:      The generalized Golub-Kahan bidiagonalization preconditioner (GKB) can be applied to symmetric $2 \times 2$ block matrices of the shape
3021:      ```{math}
3022:      \left[\begin{array}{cc} A_{00} & A_{01} \\ A_{01}' & 0 \end{array}\right]
3023:      ```
3024:      with $A_{00}$ positive semi-definite. The implementation follows {cite}`Arioli2013`. Therein, we choose $N := 1/\nu * I$ and the $(1,1)$-block of the matrix is modified to $H = _{A00} + \nu*A_{01}*A_{01}'$.
3025:      A linear system $Hx = b$ has to be solved in each iteration of the GKB algorithm. This solver is chosen with the option prefix `-fieldsplit_0_`.

3027:    Developer Note:
3028:    The Schur complement functionality of `PCFIELDSPLIT` should likely be factored into its own `PC` thus simplifying the implementation of the preconditioners and their
3029:    user API.

3031:      References:
3032:      ```{bibliography}
3033:      :filter: docname in docnames
3034:      ```

3036: .seealso: [](sec_block_matrices), `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCLSC`,
3037:           `PCFieldSplitGetSubKSP()`, `PCFieldSplitSchurGetSubKSP()`, `PCFieldSplitSetFields()`,
3038:           `PCFieldSplitSetType()`, `PCFieldSplitSetIS()`, `PCFieldSplitSetSchurPre()`, `PCFieldSplitSetSchurFactType()`,
3039:           `MatSchurComplementSetAinvType()`, `PCFieldSplitSetSchurScale()`, `PCFieldSplitSetDetectSaddlePoint()`
3040: M*/

3042: PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
3043: {
3044:   PC_FieldSplit *jac;

3046:   PetscNew(&jac);

3048:   jac->bs                 = -1;
3049:   jac->nsplits            = 0;
3050:   jac->type               = PC_COMPOSITE_MULTIPLICATIVE;
3051:   jac->schurpre           = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
3052:   jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
3053:   jac->schurscale         = -1.0;
3054:   jac->dm_splits          = PETSC_TRUE;
3055:   jac->detect             = PETSC_FALSE;
3056:   jac->gkbtol             = 1e-5;
3057:   jac->gkbdelay           = 5;
3058:   jac->gkbnu              = 1;
3059:   jac->gkbmaxit           = 100;
3060:   jac->gkbmonitor         = PETSC_FALSE;
3061:   jac->coordinates_set    = PETSC_FALSE;

3063:   pc->data = (void *)jac;

3065:   pc->ops->apply           = PCApply_FieldSplit;
3066:   pc->ops->applytranspose  = PCApplyTranspose_FieldSplit;
3067:   pc->ops->setup           = PCSetUp_FieldSplit;
3068:   pc->ops->reset           = PCReset_FieldSplit;
3069:   pc->ops->destroy         = PCDestroy_FieldSplit;
3070:   pc->ops->setfromoptions  = PCSetFromOptions_FieldSplit;
3071:   pc->ops->view            = PCView_FieldSplit;
3072:   pc->ops->applyrichardson = NULL;

3074:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSchurGetSubKSP_C", PCFieldSplitSchurGetSubKSP_FieldSplit);
3075:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit);
3076:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetFields_C", PCFieldSplitSetFields_FieldSplit);
3077:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetIS_C", PCFieldSplitSetIS_FieldSplit);
3078:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetType_C", PCFieldSplitSetType_FieldSplit);
3079:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetBlockSize_C", PCFieldSplitSetBlockSize_FieldSplit);
3080:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitRestrictIS_C", PCFieldSplitRestrictIS_FieldSplit);
3081:   PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_FieldSplit);
3082:   return 0;
3083: }