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:     }
595:     bs = jac->bs;
596:     MatGetOwnershipRange(pc->pmat, &rstart, &rend);
597:     nslots = (rend - rstart) / bs;
598:     for (i = 0; i < nsplit; i++) {
599:       if (jac->defaultsplit) {
600:         ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + i, nsplit, &ilink->is);
601:         ISDuplicate(ilink->is, &ilink->is_col);
602:       } else if (!ilink->is) {
603:         if (ilink->nfields > 1) {
604:           PetscInt *ii, *jj, j, k, nfields = ilink->nfields, *fields = ilink->fields, *fields_col = ilink->fields_col;
605:           PetscMalloc1(ilink->nfields * nslots, &ii);
606:           PetscMalloc1(ilink->nfields * nslots, &jj);
607:           for (j = 0; j < nslots; j++) {
608:             for (k = 0; k < nfields; k++) {
609:               ii[nfields * j + k] = rstart + bs * j + fields[k];
610:               jj[nfields * j + k] = rstart + bs * j + fields_col[k];
611:             }
612:           }
613:           ISCreateGeneral(PetscObjectComm((PetscObject)pc), nslots * nfields, ii, PETSC_OWN_POINTER, &ilink->is);
614:           ISCreateGeneral(PetscObjectComm((PetscObject)pc), nslots * nfields, jj, PETSC_OWN_POINTER, &ilink->is_col);
615:           ISSetBlockSize(ilink->is, nfields);
616:           ISSetBlockSize(ilink->is_col, nfields);
617:         } else {
618:           ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + ilink->fields[0], bs, &ilink->is);
619:           ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + ilink->fields_col[0], bs, &ilink->is_col);
620:         }
621:       }
622:       ISSorted(ilink->is, &sorted);
623:       if (ilink->is_col) ISSorted(ilink->is_col, &sorted_col);
625:       ilink = ilink->next;
626:     }
627:   }

629:   ilink = jac->head;
630:   if (!jac->pmat) {
631:     Vec xtmp;

633:     MatCreateVecs(pc->pmat, &xtmp, NULL);
634:     PetscMalloc1(nsplit, &jac->pmat);
635:     PetscMalloc2(nsplit, &jac->x, nsplit, &jac->y);
636:     for (i = 0; i < nsplit; i++) {
637:       MatNullSpace sp;

639:       /* Check for preconditioning matrix attached to IS */
640:       PetscObjectQuery((PetscObject)ilink->is, "pmat", (PetscObject *)&jac->pmat[i]);
641:       if (jac->pmat[i]) {
642:         PetscObjectReference((PetscObject)jac->pmat[i]);
643:         if (jac->type == PC_COMPOSITE_SCHUR) {
644:           jac->schur_user = jac->pmat[i];

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

674:     for (i = 0; i < nsplit; i++) {
675:       Mat pmat;

677:       /* Check for preconditioning matrix attached to IS */
678:       PetscObjectQuery((PetscObject)ilink->is, "pmat", (PetscObject *)&pmat);
679:       if (!pmat) MatCreateSubMatrix(pc->pmat, ilink->is, ilink->is_col, scall, &jac->pmat[i]);
680:       ilink = ilink->next;
681:     }
682:   }
683:   if (jac->diag_use_amat) {
684:     ilink = jac->head;
685:     if (!jac->mat) {
686:       PetscMalloc1(nsplit, &jac->mat);
687:       for (i = 0; i < nsplit; i++) {
688:         MatCreateSubMatrix(pc->mat, ilink->is, ilink->is_col, MAT_INITIAL_MATRIX, &jac->mat[i]);
689:         ilink = ilink->next;
690:       }
691:     } else {
692:       MatReuse scall;
693:       if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
694:         for (i = 0; i < nsplit; i++) MatDestroy(&jac->mat[i]);
695:         scall = MAT_INITIAL_MATRIX;
696:       } else scall = MAT_REUSE_MATRIX;

698:       for (i = 0; i < nsplit; i++) {
699:         MatCreateSubMatrix(pc->mat, ilink->is, ilink->is_col, scall, &jac->mat[i]);
700:         ilink = ilink->next;
701:       }
702:     }
703:   } else {
704:     jac->mat = jac->pmat;
705:   }

707:   /* Check for null space attached to IS */
708:   ilink = jac->head;
709:   for (i = 0; i < nsplit; i++) {
710:     MatNullSpace sp;

712:     PetscObjectQuery((PetscObject)ilink->is, "nullspace", (PetscObject *)&sp);
713:     if (sp) MatSetNullSpace(jac->mat[i], sp);
714:     ilink = ilink->next;
715:   }

717:   if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR && jac->type != PC_COMPOSITE_GKB) {
718:     /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */
719:     /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */
720:     ilink = jac->head;
721:     if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) {
722:       /* 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 */
723:       if (!jac->Afield) {
724:         PetscCalloc1(nsplit, &jac->Afield);
725:         if (jac->offdiag_use_amat) {
726:           MatCreateSubMatrix(pc->mat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->Afield[1]);
727:         } else {
728:           MatCreateSubMatrix(pc->pmat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->Afield[1]);
729:         }
730:       } else {
731:         MatReuse scall;

733:         if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
734:           MatDestroy(&jac->Afield[1]);
735:           scall = MAT_INITIAL_MATRIX;
736:         } else scall = MAT_REUSE_MATRIX;

738:         if (jac->offdiag_use_amat) {
739:           MatCreateSubMatrix(pc->mat, ilink->next->is, ilink->is, scall, &jac->Afield[1]);
740:         } else {
741:           MatCreateSubMatrix(pc->pmat, ilink->next->is, ilink->is, scall, &jac->Afield[1]);
742:         }
743:       }
744:     } else {
745:       if (!jac->Afield) {
746:         PetscMalloc1(nsplit, &jac->Afield);
747:         for (i = 0; i < nsplit; i++) {
748:           if (jac->offdiag_use_amat) {
749:             MatCreateSubMatrix(pc->mat, ilink->is, NULL, MAT_INITIAL_MATRIX, &jac->Afield[i]);
750:           } else {
751:             MatCreateSubMatrix(pc->pmat, ilink->is, NULL, MAT_INITIAL_MATRIX, &jac->Afield[i]);
752:           }
753:           ilink = ilink->next;
754:         }
755:       } else {
756:         MatReuse scall;
757:         if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
758:           for (i = 0; i < nsplit; i++) MatDestroy(&jac->Afield[i]);
759:           scall = MAT_INITIAL_MATRIX;
760:         } else scall = MAT_REUSE_MATRIX;

762:         for (i = 0; i < nsplit; i++) {
763:           if (jac->offdiag_use_amat) {
764:             MatCreateSubMatrix(pc->mat, ilink->is, NULL, scall, &jac->Afield[i]);
765:           } else {
766:             MatCreateSubMatrix(pc->pmat, ilink->is, NULL, scall, &jac->Afield[i]);
767:           }
768:           ilink = ilink->next;
769:         }
770:       }
771:     }
772:   }

774:   if (jac->type == PC_COMPOSITE_SCHUR) {
775:     IS          ccis;
776:     PetscBool   isset, isspd;
777:     PetscInt    rstart, rend;
778:     char        lscname[256];
779:     PetscObject LSC_L;


783:     /* If pc->mat is SPD, don't scale by -1 the Schur complement */
784:     if (jac->schurscale == (PetscScalar)-1.0) {
785:       MatIsSPDKnown(pc->pmat, &isset, &isspd);
786:       jac->schurscale = (isset && isspd) ? 1.0 : -1.0;
787:     }

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

792:     if (jac->schur) {
793:       KSP      kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper;
794:       MatReuse scall;

796:       if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
797:         scall = MAT_INITIAL_MATRIX;
798:         MatDestroy(&jac->B);
799:         MatDestroy(&jac->C);
800:       } else scall = MAT_REUSE_MATRIX;

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

835:       /* extract the A01 and A10 matrices */
836:       ilink = jac->head;
837:       ISComplement(ilink->is_col, rstart, rend, &ccis);
838:       if (jac->offdiag_use_amat) {
839:         MatCreateSubMatrix(pc->mat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->B);
840:       } else {
841:         MatCreateSubMatrix(pc->pmat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->B);
842:       }
843:       ISDestroy(&ccis);
844:       ilink = ilink->next;
845:       ISComplement(ilink->is_col, rstart, rend, &ccis);
846:       if (jac->offdiag_use_amat) {
847:         MatCreateSubMatrix(pc->mat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->C);
848:       } else {
849:         MatCreateSubMatrix(pc->pmat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->C);
850:       }
851:       ISDestroy(&ccis);

853:       /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
854:       MatCreate(((PetscObject)jac->mat[0])->comm, &jac->schur);
855:       MatSetType(jac->schur, MATSCHURCOMPLEMENT);
856:       MatSchurComplementSetSubMatrices(jac->schur, jac->mat[0], jac->pmat[0], jac->B, jac->C, jac->mat[1]);
857:       PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
858:       MatSetOptionsPrefix(jac->schur, schurmatprefix);
859:       MatSchurComplementGetKSP(jac->schur, &kspt);
860:       KSPSetOptionsPrefix(kspt, schurmatprefix);

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

866:       PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);
867:       PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
868:       if (flg) {
869:         DM  dmInner;
870:         KSP kspInner;
871:         PC  pcInner;

873:         MatSchurComplementGetKSP(jac->schur, &kspInner);
874:         KSPReset(kspInner);
875:         KSPSetOperators(kspInner, jac->mat[0], jac->pmat[0]);
876:         PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
877:         /* Indent this deeper to emphasize the "inner" nature of this solver. */
878:         PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject)pc, 2);
879:         PetscObjectIncrementTabLevel((PetscObject)kspInner->pc, (PetscObject)pc, 2);
880:         KSPSetOptionsPrefix(kspInner, schurprefix);

882:         /* Set DM for new solver */
883:         KSPGetDM(jac->head->ksp, &dmInner);
884:         KSPSetDM(kspInner, dmInner);
885:         KSPSetDMActive(kspInner, PETSC_FALSE);

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

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

910:         MatSchurComplementGetKSP(jac->schur, &kspInner);
911:         KSPGetPC(kspInner, &pcInner);
912:         PetscObjectTypeCompare((PetscObject)pcInner, PCKSP, &flg);
913:         if (flg) {
914:           KSP ksp;

916:           PCKSPGetKSP(pcInner, &ksp);
917:           if (ksp == jac->head->ksp) PCSetUseAmat(pcInner, PETSC_TRUE);
918:         }
919:       }
920:       PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);
921:       PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
922:       if (flg) {
923:         DM dmInner;

925:         PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
926:         KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);
927:         KSPSetErrorIfNotConverged(jac->kspupper, pc->erroriffailure);
928:         KSPSetOptionsPrefix(jac->kspupper, schurprefix);
929:         PetscObjectIncrementTabLevel((PetscObject)jac->kspupper, (PetscObject)pc, 1);
930:         PetscObjectIncrementTabLevel((PetscObject)jac->kspupper->pc, (PetscObject)pc, 1);
931:         KSPGetDM(jac->head->ksp, &dmInner);
932:         KSPSetDM(jac->kspupper, dmInner);
933:         KSPSetDMActive(jac->kspupper, PETSC_FALSE);
934:         KSPSetFromOptions(jac->kspupper);
935:         KSPSetOperators(jac->kspupper, jac->mat[0], jac->pmat[0]);
936:         VecDuplicate(jac->head->x, &jac->head->z);
937:       } else {
938:         jac->kspupper = jac->head->ksp;
939:         PetscObjectReference((PetscObject)jac->head->ksp);
940:       }

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

973:     /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */
974:     PetscSNPrintf(lscname, sizeof(lscname), "%s_LSC_L", ilink->splitname);
975:     PetscObjectQuery((PetscObject)pc->mat, lscname, (PetscObject *)&LSC_L);
976:     if (!LSC_L) PetscObjectQuery((PetscObject)pc->pmat, lscname, (PetscObject *)&LSC_L);
977:     if (LSC_L) PetscObjectCompose((PetscObject)jac->schur, "LSC_L", (PetscObject)LSC_L);
978:     PetscSNPrintf(lscname, sizeof(lscname), "%s_LSC_Lp", ilink->splitname);
979:     PetscObjectQuery((PetscObject)pc->pmat, lscname, (PetscObject *)&LSC_L);
980:     if (!LSC_L) PetscObjectQuery((PetscObject)pc->mat, lscname, (PetscObject *)&LSC_L);
981:     if (LSC_L) PetscObjectCompose((PetscObject)jac->schur, "LSC_Lp", (PetscObject)LSC_L);
982:   } else if (jac->type == PC_COMPOSITE_GKB) {
983:     IS       ccis;
984:     PetscInt rstart, rend;


988:     ilink = jac->head;

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

993:     ISComplement(ilink->is_col, rstart, rend, &ccis);
994:     if (jac->offdiag_use_amat) {
995:       MatCreateSubMatrix(pc->mat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->B);
996:     } else {
997:       MatCreateSubMatrix(pc->pmat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->B);
998:     }
999:     ISDestroy(&ccis);
1000:     /* Create work vectors for GKB algorithm */
1001:     VecDuplicate(ilink->x, &jac->u);
1002:     VecDuplicate(ilink->x, &jac->Hu);
1003:     VecDuplicate(ilink->x, &jac->w2);
1004:     ilink = ilink->next;
1005:     ISComplement(ilink->is_col, rstart, rend, &ccis);
1006:     if (jac->offdiag_use_amat) {
1007:       MatCreateSubMatrix(pc->mat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->C);
1008:     } else {
1009:       MatCreateSubMatrix(pc->pmat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->C);
1010:     }
1011:     ISDestroy(&ccis);
1012:     /* Create work vectors for GKB algorithm */
1013:     VecDuplicate(ilink->x, &jac->v);
1014:     VecDuplicate(ilink->x, &jac->d);
1015:     VecDuplicate(ilink->x, &jac->w1);
1016:     MatGolubKahanComputeExplicitOperator(jac->mat[0], jac->B, jac->C, &jac->H, jac->gkbnu);
1017:     PetscCalloc1(jac->gkbdelay, &jac->vecz);

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

1044:   /* Set coordinates to the sub PC objects whenever these are set */
1045:   if (jac->coordinates_set) {
1046:     PC pc_coords;
1047:     if (jac->type == PC_COMPOSITE_SCHUR) {
1048:       // Head is first block.
1049:       KSPGetPC(jac->head->ksp, &pc_coords);
1050:       PCSetCoordinates(pc_coords, jac->head->dim, jac->head->ndofs, jac->head->coords);
1051:       // Second one is Schur block, but its KSP object is in kspschur.
1052:       KSPGetPC(jac->kspschur, &pc_coords);
1053:       PCSetCoordinates(pc_coords, jac->head->next->dim, jac->head->next->ndofs, jac->head->next->coords);
1054:     } else if (jac->type == PC_COMPOSITE_GKB) {
1055:       PetscInfo(pc, "Warning: Setting coordinates does nothing for the GKB Fieldpslit preconditioner");
1056:     } else {
1057:       ilink = jac->head;
1058:       while (ilink) {
1059:         KSPGetPC(ilink->ksp, &pc_coords);
1060:         PCSetCoordinates(pc_coords, ilink->dim, ilink->ndofs, ilink->coords);
1061:         ilink = ilink->next;
1062:       }
1063:     }
1064:   }

1066:   jac->suboptionsset = PETSC_TRUE;
1067:   return 0;
1068: }

1070: #define FieldSplitSplitSolveAdd(ilink, xx, yy) \
1071:   (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) || \
1072:    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) || \
1073:    VecScatterEnd(ilink->sctx, ilink->y, yy, ADD_VALUES, SCATTER_REVERSE))

1075: static PetscErrorCode PCApply_FieldSplit_Schur(PC pc, Vec x, Vec y)
1076: {
1077:   PC_FieldSplit    *jac    = (PC_FieldSplit *)pc->data;
1078:   PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
1079:   KSP               kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;

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

1157:     PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL);
1158:     KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y);
1159:     KSPCheckSolve(jac->kspschur, pc, ilinkD->y);
1160:     PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL);
1161:     VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE);

1163:     if (kspUpper == kspA) {
1164:       MatMult(jac->B, ilinkD->y, ilinkA->y);
1165:       VecAXPY(ilinkA->x, -1.0, ilinkA->y);
1166:       PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL);
1167:       KSPSolve(kspA, ilinkA->x, ilinkA->y);
1168:       KSPCheckSolve(kspA, pc, ilinkA->y);
1169:       PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL);
1170:     } else {
1171:       PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL);
1172:       KSPSolve(kspA, ilinkA->x, ilinkA->y);
1173:       KSPCheckSolve(kspA, pc, ilinkA->y);
1174:       MatMult(jac->B, ilinkD->y, ilinkA->x);
1175:       PetscLogEventBegin(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->z, NULL);
1176:       KSPSolve(kspUpper, ilinkA->x, ilinkA->z);
1177:       KSPCheckSolve(kspUpper, pc, ilinkA->z);
1178:       PetscLogEventEnd(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->z, NULL);
1179:       VecAXPY(ilinkA->y, -1.0, ilinkA->z);
1180:     }
1181:     VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE);
1182:     VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE);
1183:     VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE);
1184:   }
1185:   return 0;
1186: }

1188: static PetscErrorCode PCApply_FieldSplit(PC pc, Vec x, Vec y)
1189: {
1190:   PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
1191:   PC_FieldSplitLink ilink = jac->head;
1192:   PetscInt          cnt, bs;

1194:   if (jac->type == PC_COMPOSITE_ADDITIVE) {
1195:     if (jac->defaultsplit) {
1196:       VecGetBlockSize(x, &bs);
1198:       VecGetBlockSize(y, &bs);
1200:       VecStrideGatherAll(x, jac->x, INSERT_VALUES);
1201:       while (ilink) {
1202:         PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1203:         KSPSolve(ilink->ksp, ilink->x, ilink->y);
1204:         KSPCheckSolve(ilink->ksp, pc, ilink->y);
1205:         PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1206:         ilink = ilink->next;
1207:       }
1208:       VecStrideScatterAll(jac->y, y, INSERT_VALUES);
1209:     } else {
1210:       VecSet(y, 0.0);
1211:       while (ilink) {
1212:         FieldSplitSplitSolveAdd(ilink, x, y);
1213:         ilink = ilink->next;
1214:       }
1215:     }
1216:   } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) {
1217:     VecSet(y, 0.0);
1218:     /* solve on first block for first block variables */
1219:     VecScatterBegin(ilink->sctx, x, ilink->x, INSERT_VALUES, SCATTER_FORWARD);
1220:     VecScatterEnd(ilink->sctx, x, ilink->x, INSERT_VALUES, SCATTER_FORWARD);
1221:     PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1222:     KSPSolve(ilink->ksp, ilink->x, ilink->y);
1223:     KSPCheckSolve(ilink->ksp, pc, ilink->y);
1224:     PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL);
1225:     VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE);
1226:     VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE);

1228:     /* compute the residual only onto second block variables using first block variables */
1229:     MatMult(jac->Afield[1], ilink->y, ilink->next->x);
1230:     ilink = ilink->next;
1231:     VecScale(ilink->x, -1.0);
1232:     VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD);
1233:     VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD);

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

1285: static PetscErrorCode PCApply_FieldSplit_GKB(PC pc, Vec x, Vec y)
1286: {
1287:   PC_FieldSplit    *jac    = (PC_FieldSplit *)pc->data;
1288:   PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
1289:   KSP               ksp = ilinkA->ksp;
1290:   Vec               u, v, Hu, d, work1, work2;
1291:   PetscScalar       alpha, z, nrmz2, *vecz;
1292:   PetscReal         lowbnd, nu, beta;
1293:   PetscInt          j, iterGKB;

1295:   VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD);
1296:   VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD);
1297:   VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD);
1298:   VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD);

1300:   u     = jac->u;
1301:   v     = jac->v;
1302:   Hu    = jac->Hu;
1303:   d     = jac->d;
1304:   work1 = jac->w1;
1305:   work2 = jac->w2;
1306:   vecz  = jac->vecz;

1308:   /* Change RHS to comply with matrix regularization H = A + nu*B*B' */
1309:   /* Add q = q + nu*B*b */
1310:   if (jac->gkbnu) {
1311:     nu = jac->gkbnu;
1312:     VecScale(ilinkD->x, jac->gkbnu);
1313:     MatMultAdd(jac->B, ilinkD->x, ilinkA->x, ilinkA->x); /* q = q + nu*B*b */
1314:   } else {
1315:     /* Situation when no augmented Lagrangian is used. Then we set inner  */
1316:     /* matrix N = I in [Ar13], and thus nu = 1.                           */
1317:     nu = 1;
1318:   }

1320:   /* Transform rhs from [q,tilde{b}] to [0,b] */
1321:   PetscLogEventBegin(ilinkA->event, ksp, ilinkA->x, ilinkA->y, NULL);
1322:   KSPSolve(ksp, ilinkA->x, ilinkA->y);
1323:   KSPCheckSolve(ksp, pc, ilinkA->y);
1324:   PetscLogEventEnd(ilinkA->event, ksp, ilinkA->x, ilinkA->y, NULL);
1325:   MatMultHermitianTranspose(jac->B, ilinkA->y, work1);
1326:   VecAXPBY(work1, 1.0 / nu, -1.0, ilinkD->x); /* c = b - B'*x        */

1328:   /* First step of algorithm */
1329:   VecNorm(work1, NORM_2, &beta); /* beta = sqrt(nu*c'*c)*/
1330:   KSPCheckDot(ksp, beta);
1331:   beta = PetscSqrtReal(nu) * beta;
1332:   VecAXPBY(v, nu / beta, 0.0, work1); /* v = nu/beta *c      */
1333:   MatMult(jac->B, v, work2);          /* u = H^{-1}*B*v      */
1334:   PetscLogEventBegin(ilinkA->event, ksp, work2, u, NULL);
1335:   KSPSolve(ksp, work2, u);
1336:   KSPCheckSolve(ksp, pc, u);
1337:   PetscLogEventEnd(ilinkA->event, ksp, work2, u, NULL);
1338:   MatMult(jac->H, u, Hu); /* alpha = u'*H*u      */
1339:   VecDot(Hu, u, &alpha);
1340:   KSPCheckDot(ksp, alpha);
1342:   alpha = PetscSqrtReal(PetscAbsScalar(alpha));
1343:   VecScale(u, 1.0 / alpha);
1344:   VecAXPBY(d, 1.0 / alpha, 0.0, v); /* v = nu/beta *c      */

1346:   z       = beta / alpha;
1347:   vecz[1] = z;

1349:   /* Computation of first iterate x(1) and p(1) */
1350:   VecAXPY(ilinkA->y, z, u);
1351:   VecCopy(d, ilinkD->y);
1352:   VecScale(ilinkD->y, -z);

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

1358:   while (iterGKB < jac->gkbmaxit && lowbnd > jac->gkbtol) {
1359:     iterGKB += 1;
1360:     MatMultHermitianTranspose(jac->B, u, work1); /* v <- nu*(B'*u-alpha/nu*v) */
1361:     VecAXPBY(v, nu, -alpha, work1);
1362:     VecNorm(v, NORM_2, &beta); /* beta = sqrt(nu)*v'*v      */
1363:     beta = beta / PetscSqrtReal(nu);
1364:     VecScale(v, 1.0 / beta);
1365:     MatMult(jac->B, v, work2); /* u <- H^{-1}*(B*v-beta*H*u) */
1366:     MatMult(jac->H, u, Hu);
1367:     VecAXPY(work2, -beta, Hu);
1368:     PetscLogEventBegin(ilinkA->event, ksp, work2, u, NULL);
1369:     KSPSolve(ksp, work2, u);
1370:     KSPCheckSolve(ksp, pc, u);
1371:     PetscLogEventEnd(ilinkA->event, ksp, work2, u, NULL);
1372:     MatMult(jac->H, u, Hu); /* alpha = u'*H*u            */
1373:     VecDot(Hu, u, &alpha);
1374:     KSPCheckDot(ksp, alpha);
1376:     alpha = PetscSqrtReal(PetscAbsScalar(alpha));
1377:     VecScale(u, 1.0 / alpha);

1379:     z       = -beta / alpha * z; /* z <- beta/alpha*z     */
1380:     vecz[0] = z;

1382:     /* Computation of new iterate x(i+1) and p(i+1) */
1383:     VecAXPBY(d, 1.0 / alpha, -beta / alpha, v); /* d = (v-beta*d)/alpha */
1384:     VecAXPY(ilinkA->y, z, u);                   /* r = r + z*u          */
1385:     VecAXPY(ilinkD->y, -z, d);                  /* p = p - z*d          */
1386:     MatMult(jac->H, ilinkA->y, Hu);             /* ||u||_H = u'*H*u     */
1387:     VecDot(Hu, ilinkA->y, &nrmz2);

1389:     /* Compute Lower Bound estimate */
1390:     if (iterGKB > jac->gkbdelay) {
1391:       lowbnd = 0.0;
1392:       for (j = 0; j < jac->gkbdelay; j++) lowbnd += PetscAbsScalar(vecz[j] * vecz[j]);
1393:       lowbnd = PetscSqrtReal(lowbnd / PetscAbsScalar(nrmz2));
1394:     }

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

1400:   VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE);
1401:   VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE);
1402:   VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE);
1403:   VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE);

1405:   return 0;
1406: }

1408: #define FieldSplitSplitSolveAddTranspose(ilink, xx, yy) \
1409:   (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) || \
1410:    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) || \
1411:    VecScatterEnd(ilink->sctx, ilink->x, yy, ADD_VALUES, SCATTER_REVERSE))

1413: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc, Vec x, Vec y)
1414: {
1415:   PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
1416:   PC_FieldSplitLink ilink = jac->head;
1417:   PetscInt          bs;

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

1477: static PetscErrorCode PCReset_FieldSplit(PC pc)
1478: {
1479:   PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
1480:   PC_FieldSplitLink ilink = jac->head, next;

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

1527: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
1528: {
1529:   PCReset_FieldSplit(pc);
1530:   PetscFree(pc->data);
1531:   PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL);
1532:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetFields_C", NULL);
1533:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetIS_C", NULL);
1534:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetType_C", NULL);
1535:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetBlockSize_C", NULL);
1536:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitRestrictIS_C", NULL);
1537:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSchurGetSubKSP_C", NULL);
1538:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", NULL);

1540:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", NULL);
1541:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", NULL);
1542:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", NULL);
1543:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", NULL);
1544:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", NULL);
1545:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", NULL);
1546:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", NULL);
1547:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", NULL);
1548:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", NULL);
1549:   return 0;
1550: }

1552: static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc, PetscOptionItems *PetscOptionsObject)
1553: {
1554:   PetscInt        bs;
1555:   PetscBool       flg;
1556:   PC_FieldSplit  *jac = (PC_FieldSplit *)pc->data;
1557:   PCCompositeType ctype;

1559:   PetscOptionsHeadBegin(PetscOptionsObject, "FieldSplit options");
1560:   PetscOptionsBool("-pc_fieldsplit_dm_splits", "Whether to use DMCreateFieldDecomposition() for splits", "PCFieldSplitSetDMSplits", jac->dm_splits, &jac->dm_splits, NULL);
1561:   PetscOptionsInt("-pc_fieldsplit_block_size", "Blocksize that defines number of fields", "PCFieldSplitSetBlockSize", jac->bs, &bs, &flg);
1562:   if (flg) PCFieldSplitSetBlockSize(pc, bs);
1563:   jac->diag_use_amat = pc->useAmat;
1564:   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);
1565:   jac->offdiag_use_amat = pc->useAmat;
1566:   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);
1567:   PetscOptionsBool("-pc_fieldsplit_detect_saddle_point", "Form 2-way split by detecting zero diagonal entries", "PCFieldSplitSetDetectSaddlePoint", jac->detect, &jac->detect, NULL);
1568:   PCFieldSplitSetDetectSaddlePoint(pc, jac->detect); /* Sets split type and Schur PC type */
1569:   PetscOptionsEnum("-pc_fieldsplit_type", "Type of composition", "PCFieldSplitSetType", PCCompositeTypes, (PetscEnum)jac->type, (PetscEnum *)&ctype, &flg);
1570:   if (flg) PCFieldSplitSetType(pc, ctype);
1571:   /* Only setup fields once */
1572:   if ((jac->bs > 0) && (jac->nsplits == 0)) {
1573:     /* only allow user to set fields from command line if bs is already known.
1574:        otherwise user can set them in PCFieldSplitSetDefaults() */
1575:     PCFieldSplitSetRuntimeSplits_Private(pc);
1576:     if (jac->splitdefined) PetscInfo(pc, "Splits defined using the options database\n");
1577:   }
1578:   if (jac->type == PC_COMPOSITE_SCHUR) {
1579:     PetscOptionsGetEnum(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_fieldsplit_schur_factorization_type", PCFieldSplitSchurFactTypes, (PetscEnum *)&jac->schurfactorization, &flg);
1580:     if (flg) PetscInfo(pc, "Deprecated use of -pc_fieldsplit_schur_factorization_type\n");
1581:     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);
1582:     PetscOptionsEnum("-pc_fieldsplit_schur_precondition", "How to build preconditioner for Schur complement", "PCFieldSplitSetSchurPre", PCFieldSplitSchurPreTypes, (PetscEnum)jac->schurpre, (PetscEnum *)&jac->schurpre, NULL);
1583:     PetscOptionsScalar("-pc_fieldsplit_schur_scale", "Scale Schur complement", "PCFieldSplitSetSchurScale", jac->schurscale, &jac->schurscale, NULL);
1584:   } else if (jac->type == PC_COMPOSITE_GKB) {
1585:     PetscOptionsReal("-pc_fieldsplit_gkb_tol", "The tolerance for the lower bound stopping criterion", "PCFieldSplitGKBTol", jac->gkbtol, &jac->gkbtol, NULL);
1586:     PetscOptionsInt("-pc_fieldsplit_gkb_delay", "The delay value for lower bound criterion", "PCFieldSplitGKBDelay", jac->gkbdelay, &jac->gkbdelay, NULL);
1587:     PetscOptionsReal("-pc_fieldsplit_gkb_nu", "Parameter in augmented Lagrangian approach", "PCFieldSplitGKBNu", jac->gkbnu, &jac->gkbnu, NULL);
1589:     PetscOptionsInt("-pc_fieldsplit_gkb_maxit", "Maximum allowed number of iterations", "PCFieldSplitGKBMaxit", jac->gkbmaxit, &jac->gkbmaxit, NULL);
1590:     PetscOptionsBool("-pc_fieldsplit_gkb_monitor", "Prints number of GKB iterations and error", "PCFieldSplitGKB", jac->gkbmonitor, &jac->gkbmonitor, NULL);
1591:   }
1592:   /*
1593:     In the initial call to this routine the sub-solver data structures do not exist so we cannot call KSPSetFromOptions() on them yet.
1594:     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
1595:     is called on the outer solver in case changes were made in the options database

1597:     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()
1598:     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.
1599:     Without this extra check test p2p1fetidp_olof_full and others fail with incorrect matrix types.

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

1603:     If one captured the PetscObjectState of the options database one could skip these calls if the database has not changed from the previous call
1604:   */
1605:   if (jac->issetup) {
1606:     PC_FieldSplitLink ilink = jac->head;
1607:     if (jac->type == PC_COMPOSITE_SCHUR) {
1608:       if (jac->kspupper && jac->kspupper->totalits > 0) KSPSetFromOptions(jac->kspupper);
1609:       if (jac->kspschur && jac->kspschur->totalits > 0) KSPSetFromOptions(jac->kspschur);
1610:     }
1611:     while (ilink) {
1612:       if (ilink->ksp->totalits > 0) KSPSetFromOptions(ilink->ksp);
1613:       ilink = ilink->next;
1614:     }
1615:   }
1616:   PetscOptionsHeadEnd();
1617:   return 0;
1618: }

1620: static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc, const char splitname[], PetscInt n, const PetscInt *fields, const PetscInt *fields_col)
1621: {
1622:   PC_FieldSplit    *jac = (PC_FieldSplit *)pc->data;
1623:   PC_FieldSplitLink ilink, next = jac->head;
1624:   char              prefix[128];
1625:   PetscInt          i;

1627:   if (jac->splitdefined) {
1628:     PetscInfo(pc, "Ignoring new split \"%s\" because the splits have already been defined\n", splitname);
1629:     return 0;
1630:   }
1631:   for (i = 0; i < n; i++) {
1634:   }
1635:   PetscNew(&ilink);
1636:   if (splitname) {
1637:     PetscStrallocpy(splitname, &ilink->splitname);
1638:   } else {
1639:     PetscMalloc1(3, &ilink->splitname);
1640:     PetscSNPrintf(ilink->splitname, 2, "%" PetscInt_FMT, jac->nsplits);
1641:   }
1642:   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 */
1643:   PetscMalloc1(n, &ilink->fields);
1644:   PetscArraycpy(ilink->fields, fields, n);
1645:   PetscMalloc1(n, &ilink->fields_col);
1646:   PetscArraycpy(ilink->fields_col, fields_col, n);

1648:   ilink->nfields = n;
1649:   ilink->next    = NULL;
1650:   KSPCreate(PetscObjectComm((PetscObject)pc), &ilink->ksp);
1651:   KSPSetErrorIfNotConverged(ilink->ksp, pc->erroriffailure);
1652:   PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)pc, 1);
1653:   KSPSetType(ilink->ksp, KSPPREONLY);

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

1658:   if (!next) {
1659:     jac->head       = ilink;
1660:     ilink->previous = NULL;
1661:   } else {
1662:     while (next->next) next = next->next;
1663:     next->next      = ilink;
1664:     ilink->previous = next;
1665:   }
1666:   jac->nsplits++;
1667:   return 0;
1668: }

1670: static PetscErrorCode PCFieldSplitSchurGetSubKSP_FieldSplit(PC pc, PetscInt *n, KSP **subksp)
1671: {
1672:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

1674:   *subksp = NULL;
1675:   if (n) *n = 0;
1676:   if (jac->type == PC_COMPOSITE_SCHUR) {
1677:     PetscInt nn;

1681:     nn = jac->nsplits + (jac->kspupper != jac->head->ksp ? 1 : 0);
1682:     PetscMalloc1(nn, subksp);
1683:     (*subksp)[0] = jac->head->ksp;
1684:     (*subksp)[1] = jac->kspschur;
1685:     if (jac->kspupper != jac->head->ksp) (*subksp)[2] = jac->kspupper;
1686:     if (n) *n = nn;
1687:   }
1688:   return 0;
1689: }

1691: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc, PetscInt *n, KSP **subksp)
1692: {
1693:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

1696:   PetscMalloc1(jac->nsplits, subksp);
1697:   MatSchurComplementGetKSP(jac->schur, *subksp);

1699:   (*subksp)[1] = jac->kspschur;
1700:   if (n) *n = jac->nsplits;
1701:   return 0;
1702: }

1704: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc, PetscInt *n, KSP **subksp)
1705: {
1706:   PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
1707:   PetscInt          cnt   = 0;
1708:   PC_FieldSplitLink ilink = jac->head;

1710:   PetscMalloc1(jac->nsplits, subksp);
1711:   while (ilink) {
1712:     (*subksp)[cnt++] = ilink->ksp;
1713:     ilink            = ilink->next;
1714:   }
1716:   if (n) *n = jac->nsplits;
1717:   return 0;
1718: }

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

1723:     Input Parameters:
1724: +   pc  - the preconditioner context
1725: -   is - the index set that defines the indices to which the fieldsplit is to be restricted

1727:     Level: advanced

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

1733: .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`
1734: @*/
1735: PetscErrorCode PCFieldSplitRestrictIS(PC pc, IS isy)
1736: {
1739:   PetscTryMethod(pc, "PCFieldSplitRestrictIS_C", (PC, IS), (pc, isy));
1740:   return 0;
1741: }

1743: static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy)
1744: {
1745:   PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
1746:   PC_FieldSplitLink ilink = jac->head, next;
1747:   PetscInt          localsize, size, sizez, i;
1748:   const PetscInt   *ind, *indz;
1749:   PetscInt         *indc, *indcz;
1750:   PetscBool         flg;

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

1802: static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc, const char splitname[], IS is)
1803: {
1804:   PC_FieldSplit    *jac = (PC_FieldSplit *)pc->data;
1805:   PC_FieldSplitLink ilink, next = jac->head;
1806:   char              prefix[128];

1808:   if (jac->splitdefined) {
1809:     PetscInfo(pc, "Ignoring new split \"%s\" because the splits have already been defined\n", splitname);
1810:     return 0;
1811:   }
1812:   PetscNew(&ilink);
1813:   if (splitname) {
1814:     PetscStrallocpy(splitname, &ilink->splitname);
1815:   } else {
1816:     PetscMalloc1(8, &ilink->splitname);
1817:     PetscSNPrintf(ilink->splitname, 7, "%" PetscInt_FMT, jac->nsplits);
1818:   }
1819:   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 */
1820:   PetscObjectReference((PetscObject)is);
1821:   ISDestroy(&ilink->is);
1822:   ilink->is = is;
1823:   PetscObjectReference((PetscObject)is);
1824:   ISDestroy(&ilink->is_col);
1825:   ilink->is_col = is;
1826:   ilink->next   = NULL;
1827:   KSPCreate(PetscObjectComm((PetscObject)pc), &ilink->ksp);
1828:   KSPSetErrorIfNotConverged(ilink->ksp, pc->erroriffailure);
1829:   PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)pc, 1);
1830:   KSPSetType(ilink->ksp, KSPPREONLY);

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

1835:   if (!next) {
1836:     jac->head       = ilink;
1837:     ilink->previous = NULL;
1838:   } else {
1839:     while (next->next) next = next->next;
1840:     next->next      = ilink;
1841:     ilink->previous = next;
1842:   }
1843:   jac->nsplits++;
1844:   return 0;
1845: }

1847: /*@C
1848:     PCFieldSplitSetFields - Sets the fields that define one particular split in the field split preconditioner

1850:     Logically Collective on pc

1852:     Input Parameters:
1853: +   pc  - the preconditioner context
1854: .   splitname - name of this split, if NULL the number of the split is used
1855: .   n - the number of fields in this split
1856: -   fields - the fields in this split

1858:     Level: intermediate

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

1863:      `PCFieldSplitSetFields()` is for defining fields as strided blocks. For example, if the block
1864:      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
1865:      0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x....
1866:      where the numbered entries indicate what is in the split.

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

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

1876: .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetBlockSize()`, `PCFieldSplitSetIS()`, `PCFieldSplitRestrictIS()`
1877: @*/
1878: PetscErrorCode PCFieldSplitSetFields(PC pc, const char splitname[], PetscInt n, const PetscInt *fields, const PetscInt *fields_col)
1879: {
1884:   PetscTryMethod(pc, "PCFieldSplitSetFields_C", (PC, const char[], PetscInt, const PetscInt *, const PetscInt *), (pc, splitname, n, fields, fields_col));
1885:   return 0;
1886: }

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

1892:     Logically Collective on pc

1894:     Input Parameters:
1895: +   pc  - the preconditioner object
1896: -   flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from

1898:     Options Database Keys:
1899: .   -pc_fieldsplit_diag_use_amat - use the Amat to provide the diagonal blocks

1901:     Level: intermediate

1903: .seealso: `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitGetDiagUseAmat()`, `PCFieldSplitSetOffDiagUseAmat()`, `PCFIELDSPLIT`
1904: @*/
1905: PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc, PetscBool flg)
1906: {
1907:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1908:   PetscBool      isfs;

1911:   PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs);
1913:   jac->diag_use_amat = flg;
1914:   return 0;
1915: }

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

1921:     Logically Collective on pc

1923:     Input Parameters:
1924: .   pc  - the preconditioner object

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

1929:     Level: intermediate

1931: .seealso: `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitSetDiagUseAmat()`, `PCFieldSplitGetOffDiagUseAmat()`, `PCFIELDSPLIT`
1932: @*/
1933: PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc, PetscBool *flg)
1934: {
1935:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1936:   PetscBool      isfs;

1940:   PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs);
1942:   *flg = jac->diag_use_amat;
1943:   return 0;
1944: }

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

1950:     Logically Collective on pc

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

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

1959:     Level: intermediate

1961: .seealso: `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitGetOffDiagUseAmat()`, `PCFieldSplitSetDiagUseAmat()`, `PCFIELDSPLIT`
1962: @*/
1963: PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc, PetscBool flg)
1964: {
1965:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1966:   PetscBool      isfs;

1969:   PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs);
1971:   jac->offdiag_use_amat = flg;
1972:   return 0;
1973: }

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

1979:     Logically Collective on pc

1981:     Input Parameters:
1982: .   pc  - the preconditioner object

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

1987:     Level: intermediate

1989: .seealso: `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitSetOffDiagUseAmat()`, `PCFieldSplitGetDiagUseAmat()`, `PCFIELDSPLIT`
1990: @*/
1991: PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc, PetscBool *flg)
1992: {
1993:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
1994:   PetscBool      isfs;

1998:   PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs);
2000:   *flg = jac->offdiag_use_amat;
2001:   return 0;
2002: }

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

2007:     Logically Collective on pc

2009:     Input Parameters:
2010: +   pc  - the preconditioner context
2011: .   splitname - name of this split, if NULL the number of the split is used
2012: -   is - the index set that defines the elements in this split

2014:     Notes:
2015:     Use `PCFieldSplitSetFields()`, for splits defined by strided types.

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

2020:     Level: intermediate

2022: .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetBlockSize()`
2023: @*/
2024: PetscErrorCode PCFieldSplitSetIS(PC pc, const char splitname[], IS is)
2025: {
2029:   PetscTryMethod(pc, "PCFieldSplitSetIS_C", (PC, const char[], IS), (pc, splitname, is));
2030:   return 0;
2031: }

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

2036:     Logically Collective on pc

2038:     Input Parameters:
2039: +   pc  - the preconditioner context
2040: -   splitname - name of this split

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

2045:     Level: intermediate

2047: .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetIS()`
2048: @*/
2049: PetscErrorCode PCFieldSplitGetIS(PC pc, const char splitname[], IS *is)
2050: {
2054:   {
2055:     PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
2056:     PC_FieldSplitLink ilink = jac->head;
2057:     PetscBool         found;

2059:     *is = NULL;
2060:     while (ilink) {
2061:       PetscStrcmp(ilink->splitname, splitname, &found);
2062:       if (found) {
2063:         *is = ilink->is;
2064:         break;
2065:       }
2066:       ilink = ilink->next;
2067:     }
2068:   }
2069:   return 0;
2070: }

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

2075:     Logically Collective on pc

2077:     Input Parameters:
2078: +   pc  - the preconditioner context
2079: -   index - index of this split

2081:     Output Parameter:
2082: -   is - the index set that defines the elements in this split

2084:     Level: intermediate

2086: .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitGetIS()`, `PCFieldSplitSetIS()`
2087: @*/
2088: PetscErrorCode PCFieldSplitGetISByIndex(PC pc, PetscInt index, IS *is)
2089: {
2093:   {
2094:     PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
2095:     PC_FieldSplitLink ilink = jac->head;
2096:     PetscInt          i     = 0;

2099:     while (i < index) {
2100:       ilink = ilink->next;
2101:       ++i;
2102:     }
2103:     PCFieldSplitGetIS(pc, ilink->splitname, is);
2104:   }
2105:   return 0;
2106: }

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

2112:     Logically Collective on pc

2114:     Input Parameters:
2115: +   pc  - the preconditioner context
2116: -   bs - the block size

2118:     Level: intermediate

2120: .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`
2121: @*/
2122: PetscErrorCode PCFieldSplitSetBlockSize(PC pc, PetscInt bs)
2123: {
2126:   PetscTryMethod(pc, "PCFieldSplitSetBlockSize_C", (PC, PetscInt), (pc, bs));
2127:   return 0;
2128: }

2130: /*@C
2131:    PCFieldSplitGetSubKSP - Gets the `KSP` contexts for all splits

2133:    Collective on pc

2135:    Input Parameter:
2136: .  pc - the preconditioner context

2138:    Output Parameters:
2139: +  n - the number of splits
2140: -  subksp - the array of `KSP` contexts

2142:    Level: advanced

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

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

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

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

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

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

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

2167: .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`, `PCFieldSplitSchurGetSubKSP()`
2168: @*/
2169: PetscErrorCode PCFieldSplitGetSubKSP(PC pc, PetscInt *n, KSP *subksp[])
2170: {
2173:   PetscUseMethod(pc, "PCFieldSplitGetSubKSP_C", (PC, PetscInt *, KSP **), (pc, n, subksp));
2174:   return 0;
2175: }

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

2180:    Collective on `KSP`

2182:    Input Parameter:
2183: .  pc - the preconditioner context

2185:    Output Parameters:
2186: +  n - the number of splits
2187: -  subksp - the array of `KSP` contexts

2189:    Level: advanced

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

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

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

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

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

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

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

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

2216: .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`, `PCFieldSplitGetSubKSP()`
2217: @*/
2218: PetscErrorCode PCFieldSplitSchurGetSubKSP(PC pc, PetscInt *n, KSP *subksp[])
2219: {
2222:   PetscUseMethod(pc, "PCFieldSplitSchurGetSubKSP_C", (PC, PetscInt *, KSP **), (pc, n, subksp));
2223:   return 0;
2224: }

2226: /*@
2227:     PCFieldSplitSetSchurPre -  Indicates from what operator the preconditioner is constructucted for the Schur complement.
2228:       The default is the A11 matrix.

2230:     Collective on pc

2232:     Input Parameters:
2233: +   pc      - the preconditioner context
2234: .   ptype   - which matrix to use for preconditioning the Schur complement: `PC_FIELDSPLIT_SCHUR_PRE_A11` (default),
2235:               `PC_FIELDSPLIT_SCHUR_PRE_SELF`, `PC_FIELDSPLIT_SCHUR_PRE_USER`,
2236:               `PC_FIELDSPLIT_SCHUR_PRE_SELFP`, and `PC_FIELDSPLIT_SCHUR_PRE_FULL`
2237: -   pre - matrix to use for preconditioning, or NULL

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

2243:     Notes:
2244:     If ptype is
2245: +     a11 - the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner
2246:      matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix
2247: .     self - the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix:
2248:           The only preconditioner that currently works with this symbolic respresentation matrix object is the PCLSC
2249:           preconditioner
2250: .     user - the preconditioner for the Schur complement is generated from the user provided matrix (pre argument
2251:           to this function).
2252: .     selfp - the preconditioning for the Schur complement is generated from an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01
2253:           This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be
2254:           lumped before extracting the diagonal using the additional option -fieldsplit_1_mat_schur_complement_ainv_type lump
2255: -     full - the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation
2256:       computed internally by `PCFIELDSPLIT` (this is expensive)
2257:       useful mostly as a test that the Schur complement approach can work for your problem

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

2263:     Level: intermediate

2265: .seealso: `PCFieldSplitGetSchurPre()`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`,
2266:           `MatSchurComplementSetAinvType()`, `PCLSC`,
2267:           `PCFieldSplitSchurPreType`
2268: @*/
2269: PetscErrorCode PCFieldSplitSetSchurPre(PC pc, PCFieldSplitSchurPreType ptype, Mat pre)
2270: {
2272:   PetscTryMethod(pc, "PCFieldSplitSetSchurPre_C", (PC, PCFieldSplitSchurPreType, Mat), (pc, ptype, pre));
2273:   return 0;
2274: }

2276: PetscErrorCode PCFieldSplitSchurPrecondition(PC pc, PCFieldSplitSchurPreType ptype, Mat pre)
2277: {
2278:   return PCFieldSplitSetSchurPre(pc, ptype, pre);
2279: } /* Deprecated name */

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

2285:     Logically Collective on pc

2287:     Input Parameter:
2288: .   pc      - the preconditioner context

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

2294:     Level: intermediate

2296: .seealso: `PCFieldSplitSetSchurPre()`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, `PCLSC`,
2297:           `PCFieldSplitSchurPreType`
2298: @*/
2299: PetscErrorCode PCFieldSplitGetSchurPre(PC pc, PCFieldSplitSchurPreType *ptype, Mat *pre)
2300: {
2302:   PetscUseMethod(pc, "PCFieldSplitGetSchurPre_C", (PC, PCFieldSplitSchurPreType *, Mat *), (pc, ptype, pre));
2303:   return 0;
2304: }

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

2309:     Not collective

2311:     Input Parameter:
2312: .   pc      - the preconditioner context

2314:     Output Parameter:
2315: .   S       - the Schur complement matrix

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

2320:     Level: advanced

2322: .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurPre()`, `MATSCHURCOMPLEMENT`, `PCFieldSplitSchurRestoreS()`,
2323:           `MatCreateSchurComplement()`, `MatSchurComplementGetKSP()`, `MatSchurComplementComputeExplicitOperator()`, `MatGetSchurComplement()`
2324: @*/
2325: PetscErrorCode PCFieldSplitSchurGetS(PC pc, Mat *S)
2326: {
2327:   const char    *t;
2328:   PetscBool      isfs;
2329:   PC_FieldSplit *jac;

2332:   PetscObjectGetType((PetscObject)pc, &t);
2333:   PetscStrcmp(t, PCFIELDSPLIT, &isfs);
2335:   jac = (PC_FieldSplit *)pc->data;
2337:   if (S) *S = jac->schur;
2338:   return 0;
2339: }

2341: /*@
2342:     PCFieldSplitSchurRestoreS -  returns the `MATSCHURCOMPLEMENT` matrix used by this `PC`

2344:     Not collective

2346:     Input Parameters:
2347: +   pc      - the preconditioner context
2348: -   S       - the Schur complement matrix

2350:     Level: advanced

2352: .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurPre()`, `MatSchurComplement`, `PCFieldSplitSchurGetS()`
2353: @*/
2354: PetscErrorCode PCFieldSplitSchurRestoreS(PC pc, Mat *S)
2355: {
2356:   const char    *t;
2357:   PetscBool      isfs;
2358:   PC_FieldSplit *jac;

2361:   PetscObjectGetType((PetscObject)pc, &t);
2362:   PetscStrcmp(t, PCFIELDSPLIT, &isfs);
2364:   jac = (PC_FieldSplit *)pc->data;
2367:   return 0;
2368: }

2370: static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc, PCFieldSplitSchurPreType ptype, Mat pre)
2371: {
2372:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

2374:   jac->schurpre = ptype;
2375:   if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
2376:     MatDestroy(&jac->schur_user);
2377:     jac->schur_user = pre;
2378:     PetscObjectReference((PetscObject)jac->schur_user);
2379:   }
2380:   return 0;
2381: }

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

2387:   *ptype = jac->schurpre;
2388:   *pre   = jac->schur_user;
2389:   return 0;
2390: }

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

2395:     Collective on pc

2397:     Input Parameters:
2398: +   pc  - the preconditioner context
2399: -   ftype - which blocks of factorization to retain, `PC_FIELDSPLIT_SCHUR_FACT_FULL` is default

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

2404:     Level: intermediate

2406:     Notes:
2407:     The FULL factorization is

2409: .vb
2410:    (A   B)  = (1       0) (A   0) (1  Ainv*B)  = L D U
2411:    (C   E)    (C*Ainv  1) (0   S) (0     1)
2412: .vb
2413:     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,
2414:     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().

2416:     If A and S are solved exactly
2417: .vb
2418:       *) FULL factorization is a direct solver.
2419:       *) 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.
2420:       *) 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.
2421: .ve

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

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

2428:     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).

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

2434: .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurScale()`
2435: @*/
2436: PetscErrorCode PCFieldSplitSetSchurFactType(PC pc, PCFieldSplitSchurFactType ftype)
2437: {
2439:   PetscTryMethod(pc, "PCFieldSplitSetSchurFactType_C", (PC, PCFieldSplitSchurFactType), (pc, ftype));
2440:   return 0;
2441: }

2443: static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc, PCFieldSplitSchurFactType ftype)
2444: {
2445:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

2447:   jac->schurfactorization = ftype;
2448:   return 0;
2449: }

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

2454:     Collective on pc

2456:     Input Parameters:
2457: +   pc    - the preconditioner context
2458: -   scale - scaling factor for the Schur complement

2460:     Options Database Key:
2461: .     -pc_fieldsplit_schur_scale - default is -1.0

2463:     Level: intermediate

2465: .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurFactType`, `PCFieldSplitSetSchurScale()`, `PCFieldSplitSetSchurFactType()`
2466: @*/
2467: PetscErrorCode PCFieldSplitSetSchurScale(PC pc, PetscScalar scale)
2468: {
2471:   PetscTryMethod(pc, "PCFieldSplitSetSchurScale_C", (PC, PetscScalar), (pc, scale));
2472:   return 0;
2473: }

2475: static PetscErrorCode PCFieldSplitSetSchurScale_FieldSplit(PC pc, PetscScalar scale)
2476: {
2477:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

2479:   jac->schurscale = scale;
2480:   return 0;
2481: }

2483: /*@C
2484:    PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement

2486:    Collective on pc

2488:    Input Parameter:
2489: .  pc - the preconditioner context

2491:    Output Parameters:
2492: +  A00 - the (0,0) block
2493: .  A01 - the (0,1) block
2494: .  A10 - the (1,0) block
2495: -  A11 - the (1,1) block

2497:    Level: advanced

2499: .seealso: `PCFIELDSPLIT`, `MatSchurComplementGetSubMatrices()`, `MatSchurComplementSetSubMatrices()`
2500: @*/
2501: PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc, Mat *A00, Mat *A01, Mat *A10, Mat *A11)
2502: {
2503:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

2507:   if (A00) *A00 = jac->pmat[0];
2508:   if (A01) *A01 = jac->B;
2509:   if (A10) *A10 = jac->C;
2510:   if (A11) *A11 = jac->pmat[1];
2511:   return 0;
2512: }

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

2517:     Collective on pc

2519:     Input Parameters:
2520: +   pc        - the preconditioner context
2521: -   tolerance - the solver tolerance

2523:     Options Database Key:
2524: .   -pc_fieldsplit_gkb_tol - default is 1e-5

2526:     Level: intermediate

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

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

2535: .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBNu()`, `PCFieldSplitSetGKBMaxit()`
2536: @*/
2537: PetscErrorCode PCFieldSplitSetGKBTol(PC pc, PetscReal tolerance)
2538: {
2541:   PetscTryMethod(pc, "PCFieldSplitSetGKBTol_C", (PC, PetscReal), (pc, tolerance));
2542:   return 0;
2543: }

2545: static PetscErrorCode PCFieldSplitSetGKBTol_FieldSplit(PC pc, PetscReal tolerance)
2546: {
2547:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

2549:   jac->gkbtol = tolerance;
2550:   return 0;
2551: }

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

2556:     Collective on pc

2558:     Input Parameters:
2559: +   pc     - the preconditioner context
2560: -   maxit  - the maximum number of iterations

2562:     Options Database Key:
2563: .     -pc_fieldsplit_gkb_maxit - default is 100

2565:     Level: intermediate

2567: .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBNu()`
2568: @*/
2569: PetscErrorCode PCFieldSplitSetGKBMaxit(PC pc, PetscInt maxit)
2570: {
2573:   PetscTryMethod(pc, "PCFieldSplitSetGKBMaxit_C", (PC, PetscInt), (pc, maxit));
2574:   return 0;
2575: }

2577: static PetscErrorCode PCFieldSplitSetGKBMaxit_FieldSplit(PC pc, PetscInt maxit)
2578: {
2579:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

2581:   jac->gkbmaxit = maxit;
2582:   return 0;
2583: }

2585: /*@
2586:     PCFieldSplitSetGKBDelay -  Sets the delay in the lower bound error estimate in the generalized Golub-Kahan bidiagonalization in `PCFIELDSPLIT`
2587:     preconditioner.

2589:     Collective on pc

2591:     Input Parameters:
2592: +   pc     - the preconditioner context
2593: -   delay  - the delay window in the lower bound estimate

2595:     Options Database Key:
2596: .     -pc_fieldsplit_gkb_delay - default is 5

2598:     Level: intermediate

2600:     Note:
2601:     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
2602:     is expressed as a truncated sum. The error at iteration k can only be measured at iteration (k + delay), and thus the algorithm needs
2603:     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

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

2607: .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetGKBNu()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBMaxit()`
2608: @*/
2609: PetscErrorCode PCFieldSplitSetGKBDelay(PC pc, PetscInt delay)
2610: {
2613:   PetscTryMethod(pc, "PCFieldSplitSetGKBDelay_C", (PC, PetscInt), (pc, delay));
2614:   return 0;
2615: }

2617: static PetscErrorCode PCFieldSplitSetGKBDelay_FieldSplit(PC pc, PetscInt delay)
2618: {
2619:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

2621:   jac->gkbdelay = delay;
2622:   return 0;
2623: }

2625: /*@
2626:     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
2627:     in `PCFIELDSPLIT`

2629:     Collective on pc

2631:     Input Parameters:
2632: +   pc     - the preconditioner context
2633: -   nu     - the shift parameter

2635:     Options Database Keys:
2636: .     -pc_fieldsplit_gkb_nu - default is 1

2638:     Level: intermediate

2640:     Notes:
2641:     This shift is in general done to obtain better convergence properties for the outer loop of the algorithm. This is often achieved by chosing nu sufficiently big. However,
2642:     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 gets difficult. It is therefore
2643:     necessary to find a good balance in between the convergence of the inner and outer loop.

2645:     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.

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

2649: .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBMaxit()`
2650: @*/
2651: PetscErrorCode PCFieldSplitSetGKBNu(PC pc, PetscReal nu)
2652: {
2655:   PetscTryMethod(pc, "PCFieldSplitSetGKBNu_C", (PC, PetscReal), (pc, nu));
2656:   return 0;
2657: }

2659: static PetscErrorCode PCFieldSplitSetGKBNu_FieldSplit(PC pc, PetscReal nu)
2660: {
2661:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

2663:   jac->gkbnu = nu;
2664:   return 0;
2665: }

2667: static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc, PCCompositeType type)
2668: {
2669:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

2671:   jac->type = type;
2672:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", NULL);
2673:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", NULL);
2674:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", NULL);
2675:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", NULL);
2676:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", NULL);
2677:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", NULL);
2678:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", NULL);
2679:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", NULL);
2680:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", NULL);

2682:   if (type == PC_COMPOSITE_SCHUR) {
2683:     pc->ops->apply = PCApply_FieldSplit_Schur;
2684:     pc->ops->view  = PCView_FieldSplit_Schur;

2686:     PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit_Schur);
2687:     PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", PCFieldSplitSetSchurPre_FieldSplit);
2688:     PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", PCFieldSplitGetSchurPre_FieldSplit);
2689:     PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", PCFieldSplitSetSchurFactType_FieldSplit);
2690:     PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", PCFieldSplitSetSchurScale_FieldSplit);
2691:   } else if (type == PC_COMPOSITE_GKB) {
2692:     pc->ops->apply = PCApply_FieldSplit_GKB;
2693:     pc->ops->view  = PCView_FieldSplit_GKB;

2695:     PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit);
2696:     PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", PCFieldSplitSetGKBTol_FieldSplit);
2697:     PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", PCFieldSplitSetGKBMaxit_FieldSplit);
2698:     PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", PCFieldSplitSetGKBNu_FieldSplit);
2699:     PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", PCFieldSplitSetGKBDelay_FieldSplit);
2700:   } else {
2701:     pc->ops->apply = PCApply_FieldSplit;
2702:     pc->ops->view  = PCView_FieldSplit;

2704:     PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit);
2705:   }
2706:   return 0;
2707: }

2709: static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc, PetscInt bs)
2710: {
2711:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

2715:   jac->bs = bs;
2716:   return 0;
2717: }

2719: static PetscErrorCode PCSetCoordinates_FieldSplit(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[])
2720: {
2721:   PC_FieldSplit    *jac           = (PC_FieldSplit *)pc->data;
2722:   PC_FieldSplitLink ilink_current = jac->head;
2723:   IS                is_owned;

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

2728:   while (ilink_current) {
2729:     // For each IS, embed it to get local coords indces
2730:     IS              is_coords;
2731:     PetscInt        ndofs_block;
2732:     const PetscInt *block_dofs_enumeration; // Numbering of the dofs relevant to the current block

2734:     // Setting drop to true for safety. It should make no difference.
2735:     ISEmbed(ilink_current->is, is_owned, PETSC_TRUE, &is_coords);
2736:     ISGetLocalSize(is_coords, &ndofs_block);
2737:     ISGetIndices(is_coords, &block_dofs_enumeration);

2739:     // Allocate coordinates vector and set it directly
2740:     PetscMalloc1(ndofs_block * dim, &(ilink_current->coords));
2741:     for (PetscInt dof = 0; dof < ndofs_block; ++dof) {
2742:       for (PetscInt d = 0; d < dim; ++d) (ilink_current->coords)[dim * dof + d] = coords[dim * block_dofs_enumeration[dof] + d];
2743:     }
2744:     ilink_current->dim   = dim;
2745:     ilink_current->ndofs = ndofs_block;
2746:     ISRestoreIndices(is_coords, &block_dofs_enumeration);
2747:     ISDestroy(&is_coords);
2748:     ilink_current = ilink_current->next;
2749:   }
2750:   ISDestroy(&is_owned);
2751:   return 0;
2752: }

2754: /*@
2755:    PCFieldSplitSetType - Sets the type, `PCCompositeType`, of a `PCFIELDSPLIT`

2757:    Collective on pc

2759:    Input Parameters:
2760: +  pc - the preconditioner context
2761: -  type - `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE` (default), `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`

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

2766:    Level: Intermediate

2768: .seealso: `PCFIELDSPLIT`, `PCCompositeType`, `PCCompositeGetType()`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`,
2769:           `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`
2770: @*/
2771: PetscErrorCode PCFieldSplitSetType(PC pc, PCCompositeType type)
2772: {
2774:   PetscTryMethod(pc, "PCFieldSplitSetType_C", (PC, PCCompositeType), (pc, type));
2775:   return 0;
2776: }

2778: /*@
2779:   PCFieldSplitGetType - Gets the type, `PCCompositeType`, of a `PCFIELDSPLIT`

2781:   Not collective

2783:   Input Parameter:
2784: . pc - the preconditioner context

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

2789:   Level: Intermediate

2791: .seealso: `PCCompositeSetType()`, `PCFIELDSPLIT`, `PCCompositeType`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`,
2792:           `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`
2793: @*/
2794: PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
2795: {
2796:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

2800:   *type = jac->type;
2801:   return 0;
2802: }

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

2807:    Logically Collective on pc

2809:    Input Parameters:
2810: +  pc   - the preconditioner context
2811: -  flg  - boolean indicating whether to use field splits defined by the `DM`

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

2816:    Level: Intermediate

2818: .seealso: `PCFIELDSPLIT`, `PCFieldSplitGetDMSplits()`, `PCFieldSplitSetFields()`, `PCFieldsplitSetIS()`
2819: @*/
2820: PetscErrorCode PCFieldSplitSetDMSplits(PC pc, PetscBool flg)
2821: {
2822:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2823:   PetscBool      isfs;

2827:   PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs);
2828:   if (isfs) jac->dm_splits = flg;
2829:   return 0;
2830: }

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

2835:    Logically Collective

2837:    Input Parameter:
2838: .  pc   - the preconditioner context

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

2843:    Level: Intermediate

2845: .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetDMSplits()`, `PCFieldSplitSetFields()`, `PCFieldsplitSetIS()`
2846: @*/
2847: PetscErrorCode PCFieldSplitGetDMSplits(PC pc, PetscBool *flg)
2848: {
2849:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
2850:   PetscBool      isfs;

2854:   PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs);
2855:   if (isfs) {
2856:     if (flg) *flg = jac->dm_splits;
2857:   }
2858:   return 0;
2859: }

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

2864:    Logically Collective

2866:    Input Parameter:
2867: .  pc   - the preconditioner context

2869:    Output Parameter:
2870: .  flg  - boolean indicating whether to detect fields or not

2872:    Level: Intermediate

2874: .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetDetectSaddlePoint()`
2875: @*/
2876: PetscErrorCode PCFieldSplitGetDetectSaddlePoint(PC pc, PetscBool *flg)
2877: {
2878:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

2880:   *flg = jac->detect;
2881:   return 0;
2882: }

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

2887:    Logically Collective

2889:    Input Parameter:
2890: .  pc   - the preconditioner context

2892:    Output Parameter:
2893: .  flg  - boolean indicating whether to detect fields or not

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

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

2901:    Level: Intermediate

2903: .seealso: `PCFIELDSPLIT`, `PCFieldSplitGetDetectSaddlePoint()`, `PCFieldSplitSetType()`, `PCFieldSplitSetSchurPre()`, `PC_FIELDSPLIT_SCHUR_PRE_SELF`
2904: @*/
2905: PetscErrorCode PCFieldSplitSetDetectSaddlePoint(PC pc, PetscBool flg)
2906: {
2907:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

2909:   jac->detect = flg;
2910:   if (jac->detect) {
2911:     PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
2912:     PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_SELF, NULL);
2913:   }
2914:   return 0;
2915: }

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

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

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

2927:    Level: intermediate

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

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

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

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

2952:       For the Schur complement preconditioner if

2954:       ```{math}
2955:       J = \left[\begin{array}{cc} A_{00} & A_{01} \\ A_{10} & A_{11} \end{array}\right]
2956:       ```

2958:       the preconditioner using `full` factorization is logically
2959:       ```{math}
2960:       \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]
2961:       ```
2962:      where the action of $\text{inv}(A_{00})$ is applied using the KSP solver with prefix `-fieldsplit_0_`.  $S$ is the Schur complement
2963:      ```{math}
2964:      S = A_{11} - A_{10} \text{ksp}(A_{00}) A_{01}
2965:      ```
2966:      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
2967:      in providing the SECOND split or 1 if not given). For `PCFieldSplitGetSub\text{ksp}()` when field number is 0,
2968:      it returns the KSP associated with `-fieldsplit_0_` while field number 1 gives `-fieldsplit_1_` KSP. By default
2969:      $A_{11}$ is used to construct a preconditioner for $S$, use `PCFieldSplitSetSchurPre()` for all the possible ways to construct the preconditioner for $S$.

2971:      The factorization type is set using `-pc_fieldsplit_schur_fact_type <diag, lower, upper, full>`. `full` is shown above,
2972:      `diag` gives
2973:       ```{math}
2974:       \left[\begin{array}{cc} \text{inv}(A_{00}) & 0 \\  0 & -\text{ksp}(S) \end{array}\right]
2975:       ```
2976:      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
2977:      can be turned off with `PCFieldSplitSetSchurScale()` or by command line `-pc_fieldsplit_schur_scale 1.0`. The `lower` factorization is the inverse of
2978:       ```{math}
2979:       \left[\begin{array}{cc} A_{00} & 0 \\  A_{10} & S \end{array}\right]
2980:       ```
2981:      where the inverses of A_{00} and S are applied using KSPs. The upper factorization is the inverse of
2982:       ```{math}
2983:       \left[\begin{array}{cc} A_{00} & A_{01} \\  0 & S \end{array}\right]
2984:       ```
2985:      where again the inverses of $A_{00}$ and $S$ are applied using `KSP`s.

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

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

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

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

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

3003:      The generalized Golub-Kahan bidiagonalization preconditioner (GKB) can be applied to symmetric $2 \times 2$ block matrices of the shape
3004:      ```{math}
3005:      \left[\begin{array}{cc} A_{00} & A_{01} \\ A_{01}' & 0 \end{array}\right]
3006:      ```
3007:      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}'$.
3008:      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_`.

3010:      References:

3012:      ```{bibliography}
3013:      :filter: docname in docnames
3014:      ```

3016:    Developer Note:
3017:    The Schur complement functionality of `PCFIELDSPLIT` should likely be factored into its own `PC` thus simplify the implementation of the preconditioners and their
3018:    user API.

3020: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCLSC`,
3021:           `PCFieldSplitGetSubKSP()`, `PCFieldSplitSchurGetSubKSP()`, `PCFieldSplitSetFields()`,
3022:           `PCFieldSplitSetType()`, `PCFieldSplitSetIS()`, `PCFieldSplitSetSchurPre()`, `PCFieldSplitSetSchurFactType()`,
3023:           `MatSchurComplementSetAinvType()`, `PCFieldSplitSetSchurScale()`, `PCFieldSplitSetDetectSaddlePoint()`
3024: M*/

3026: PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
3027: {
3028:   PC_FieldSplit *jac;

3030:   PetscNew(&jac);

3032:   jac->bs                 = -1;
3033:   jac->nsplits            = 0;
3034:   jac->type               = PC_COMPOSITE_MULTIPLICATIVE;
3035:   jac->schurpre           = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
3036:   jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
3037:   jac->schurscale         = -1.0;
3038:   jac->dm_splits          = PETSC_TRUE;
3039:   jac->detect             = PETSC_FALSE;
3040:   jac->gkbtol             = 1e-5;
3041:   jac->gkbdelay           = 5;
3042:   jac->gkbnu              = 1;
3043:   jac->gkbmaxit           = 100;
3044:   jac->gkbmonitor         = PETSC_FALSE;
3045:   jac->coordinates_set    = PETSC_FALSE;

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

3049:   pc->ops->apply           = PCApply_FieldSplit;
3050:   pc->ops->applytranspose  = PCApplyTranspose_FieldSplit;
3051:   pc->ops->setup           = PCSetUp_FieldSplit;
3052:   pc->ops->reset           = PCReset_FieldSplit;
3053:   pc->ops->destroy         = PCDestroy_FieldSplit;
3054:   pc->ops->setfromoptions  = PCSetFromOptions_FieldSplit;
3055:   pc->ops->view            = PCView_FieldSplit;
3056:   pc->ops->applyrichardson = NULL;

3058:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSchurGetSubKSP_C", PCFieldSplitSchurGetSubKSP_FieldSplit);
3059:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit);
3060:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetFields_C", PCFieldSplitSetFields_FieldSplit);
3061:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetIS_C", PCFieldSplitSetIS_FieldSplit);
3062:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetType_C", PCFieldSplitSetType_FieldSplit);
3063:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetBlockSize_C", PCFieldSplitSetBlockSize_FieldSplit);
3064:   PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitRestrictIS_C", PCFieldSplitRestrictIS_FieldSplit);
3065:   PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_FieldSplit);
3066:   return 0;
3067: }