Actual source code: fieldsplit.c

petsc-master 2020-08-09
Report Typos and Errors
  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;
 21: };

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

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

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

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

 67: /*
 68:     Notes:
 69:     there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of
 70:    inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the
 71:    PC you could change this.
 72: */

 74: /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it.  This way the
 75: * Section 1.5 Writing Application Codes with PETSc-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */
 76: static Mat FieldSplitSchurPre(PC_FieldSplit *jac)
 77: {
 78:   switch (jac->schurpre) {
 79:   case PC_FIELDSPLIT_SCHUR_PRE_SELF: return jac->schur;
 80:   case PC_FIELDSPLIT_SCHUR_PRE_SELFP: return jac->schurp;
 81:   case PC_FIELDSPLIT_SCHUR_PRE_A11: return jac->pmat[1];
 82:   case PC_FIELDSPLIT_SCHUR_PRE_FULL: /* We calculate this and store it in schur_user */
 83:   case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */
 84:   default:
 85:     return jac->schur_user ? jac->schur_user : jac->pmat[1];
 86:   }
 87: }


 90:  #include <petscdraw.h>
 91: static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer)
 92: {
 93:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
 94:   PetscErrorCode    ierr;
 95:   PetscBool         iascii,isdraw;
 96:   PetscInt          i,j;
 97:   PC_FieldSplitLink ilink = jac->head;

100:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
101:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
102:   if (iascii) {
103:     if (jac->bs > 0) {
104:       PetscViewerASCIIPrintf(viewer,"  FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);
105:     } else {
106:       PetscViewerASCIIPrintf(viewer,"  FieldSplit with %s composition: total splits = %D\n",PCCompositeTypes[jac->type],jac->nsplits);
107:     }
108:     if (pc->useAmat) {
109:       PetscViewerASCIIPrintf(viewer,"  using Amat (not Pmat) as operator for blocks\n");
110:     }
111:     if (jac->diag_use_amat) {
112:       PetscViewerASCIIPrintf(viewer,"  using Amat (not Pmat) as operator for diagonal blocks\n");
113:     }
114:     if (jac->offdiag_use_amat) {
115:       PetscViewerASCIIPrintf(viewer,"  using Amat (not Pmat) as operator for off-diagonal blocks\n");
116:     }
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 %D Fields ",i);
121:         PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
122:         for (j=0; j<ilink->nfields; j++) {
123:           if (j > 0) {
124:             PetscViewerASCIIPrintf(viewer,",");
125:           }
126:           PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);
127:         }
128:         PetscViewerASCIIPrintf(viewer,"\n");
129:         PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
130:       } else {
131:         PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);
132:       }
133:       KSPView(ilink->ksp,viewer);
134:       ilink = ilink->next;
135:     }
136:   }

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

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

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

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

249:     PetscViewerDrawGetDraw(viewer,0,&draw);
250:     PetscDrawGetCurrentPoint(draw,&x,&y);
251:     if (jac->kspupper != jac->head->ksp) cnt++;
252:     w  = 2*PetscMin(1.0 - x,x);
253:     wd = w/(cnt + 1);

255:     PetscSNPrintf(str,32,"Schur fact. %s",PCFieldSplitSchurFactTypes[jac->schurfactorization]);
256:     PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
257:     y   -= h;
258:     if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_USER &&  !jac->schur_user) {
259:       PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[PC_FIELDSPLIT_SCHUR_PRE_A11]);
260:     } else {
261:       PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[jac->schurpre]);
262:     }
263:     PetscDrawStringBoxed(draw,x+wd*(cnt-1)/2.0,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
264:     y   -= h;
265:     x    = x - wd*(cnt-1)/2.0;

267:     PetscDrawPushCurrentPoint(draw,x,y);
268:     KSPView(jac->head->ksp,viewer);
269:     PetscDrawPopCurrentPoint(draw);
270:     if (jac->kspupper != jac->head->ksp) {
271:       x   += wd;
272:       PetscDrawPushCurrentPoint(draw,x,y);
273:       KSPView(jac->kspupper,viewer);
274:       PetscDrawPopCurrentPoint(draw);
275:     }
276:     x   += wd;
277:     PetscDrawPushCurrentPoint(draw,x,y);
278:     KSPView(jac->kspschur,viewer);
279:     PetscDrawPopCurrentPoint(draw);
280:   }
281:   return(0);
282: }

284: static PetscErrorCode PCView_FieldSplit_GKB(PC pc,PetscViewer viewer)
285: {
286:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
287:   PetscErrorCode    ierr;
288:   PetscBool         iascii,isdraw;
289:   PetscInt          i,j;
290:   PC_FieldSplitLink ilink = jac->head;

293:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
294:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
295:   if (iascii) {
296:     if (jac->bs > 0) {
297:       PetscViewerASCIIPrintf(viewer,"  FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);
298:     } else {
299:       PetscViewerASCIIPrintf(viewer,"  FieldSplit with %s composition: total splits = %D\n",PCCompositeTypes[jac->type],jac->nsplits);
300:     }
301:     if (pc->useAmat) {
302:       PetscViewerASCIIPrintf(viewer,"  using Amat (not Pmat) as operator for blocks\n");
303:     }
304:     if (jac->diag_use_amat) {
305:       PetscViewerASCIIPrintf(viewer,"  using Amat (not Pmat) as operator for diagonal blocks\n");
306:     }
307:     if (jac->offdiag_use_amat) {
308:       PetscViewerASCIIPrintf(viewer,"  using Amat (not Pmat) as operator for off-diagonal blocks\n");
309:     }

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

315:     if (ilink->fields) {
316:       PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",0);
317:       PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
318:       for (j=0; j<ilink->nfields; j++) {
319:         if (j > 0) {
320:           PetscViewerASCIIPrintf(viewer,",");
321:         }
322:         PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);
323:       }
324:       PetscViewerASCIIPrintf(viewer,"\n");
325:       PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
326:     } else {
327:         PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",0);
328:     }
329:     KSPView(ilink->ksp,viewer);

331:     PetscViewerASCIIPopTab(viewer);
332:   }

334:  if (isdraw) {
335:     PetscDraw draw;
336:     PetscReal x,y,w,wd;

338:     PetscViewerDrawGetDraw(viewer,0,&draw);
339:     PetscDrawGetCurrentPoint(draw,&x,&y);
340:     w    = 2*PetscMin(1.0 - x,x);
341:     wd   = w/(jac->nsplits + 1);
342:     x    = x - wd*(jac->nsplits-1)/2.0;
343:     for (i=0; i<jac->nsplits; i++) {
344:       PetscDrawPushCurrentPoint(draw,x,y);
345:       KSPView(ilink->ksp,viewer);
346:       PetscDrawPopCurrentPoint(draw);
347:       x    += wd;
348:       ilink = ilink->next;
349:     }
350:   }
351:   return(0);
352: }


355: /* Precondition: jac->bs is set to a meaningful value */
356: static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc)
357: {
359:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
360:   PetscInt       i,nfields,*ifields,nfields_col,*ifields_col;
361:   PetscBool      flg,flg_col;
362:   char           optionname[128],splitname[8],optionname_col[128];

365:   PetscMalloc1(jac->bs,&ifields);
366:   PetscMalloc1(jac->bs,&ifields_col);
367:   for (i=0,flg=PETSC_TRUE;; i++) {
368:     PetscSNPrintf(splitname,sizeof(splitname),"%D",i);
369:     PetscSNPrintf(optionname,sizeof(optionname),"-pc_fieldsplit_%D_fields",i);
370:     PetscSNPrintf(optionname_col,sizeof(optionname_col),"-pc_fieldsplit_%D_fields_col",i);
371:     nfields     = jac->bs;
372:     nfields_col = jac->bs;
373:     PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname,ifields,&nfields,&flg);
374:     PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname_col,ifields_col,&nfields_col,&flg_col);
375:     if (!flg) break;
376:     else if (flg && !flg_col) {
377:       if (!nfields) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields");
378:       PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields);
379:     } else {
380:       if (!nfields || !nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields");
381:       if (nfields != nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Number of row and column fields must match");
382:       PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields_col);
383:     }
384:   }
385:   if (i > 0) {
386:     /* Makes command-line setting of splits take precedence over setting them in code.
387:        Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would
388:        create new splits, which would probably not be what the user wanted. */
389:     jac->splitdefined = PETSC_TRUE;
390:   }
391:   PetscFree(ifields);
392:   PetscFree(ifields_col);
393:   return(0);
394: }

396: static PetscErrorCode PCFieldSplitSetDefaults(PC pc)
397: {
398:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
399:   PetscErrorCode    ierr;
400:   PC_FieldSplitLink ilink = jac->head;
401:   PetscBool         fieldsplit_default = PETSC_FALSE,coupling = PETSC_FALSE;
402:   PetscInt          i;

405:   /*
406:    Kinda messy, but at least this now uses DMCreateFieldDecomposition().
407:    Should probably be rewritten.
408:    */
409:   if (!ilink) {
410:     PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_coupling",&coupling,NULL);
411:     if (pc->dm && jac->dm_splits && !jac->detect && !coupling) {
412:       PetscInt  numFields, f, i, j;
413:       char      **fieldNames;
414:       IS        *fields;
415:       DM        *dms;
416:       DM        subdm[128];
417:       PetscBool flg;

419:       DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms);
420:       /* Allow the user to prescribe the splits */
421:       for (i = 0, flg = PETSC_TRUE;; i++) {
422:         PetscInt ifields[128];
423:         IS       compField;
424:         char     optionname[128], splitname[8];
425:         PetscInt nfields = numFields;

427:         PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%D_fields", i);
428:         PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg);
429:         if (!flg) break;
430:         if (numFields > 128) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot currently support %d > 128 fields", numFields);
431:         DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i]);
432:         if (nfields == 1) {
433:           PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField);
434:         } else {
435:           PetscSNPrintf(splitname, sizeof(splitname), "%D", i);
436:           PCFieldSplitSetIS(pc, splitname, compField);
437:         }
438:         ISDestroy(&compField);
439:         for (j = 0; j < nfields; ++j) {
440:           f    = ifields[j];
441:           PetscFree(fieldNames[f]);
442:           ISDestroy(&fields[f]);
443:         }
444:       }
445:       if (i == 0) {
446:         for (f = 0; f < numFields; ++f) {
447:           PCFieldSplitSetIS(pc, fieldNames[f], fields[f]);
448:           PetscFree(fieldNames[f]);
449:           ISDestroy(&fields[f]);
450:         }
451:       } else {
452:         for (j=0; j<numFields; j++) {
453:           DMDestroy(dms+j);
454:         }
455:         PetscFree(dms);
456:         PetscMalloc1(i, &dms);
457:         for (j = 0; j < i; ++j) dms[j] = subdm[j];
458:       }
459:       PetscFree(fieldNames);
460:       PetscFree(fields);
461:       if (dms) {
462:         PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n");
463:         for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) {
464:           const char *prefix;
465:           PetscObjectGetOptionsPrefix((PetscObject)(ilink->ksp),&prefix);
466:           PetscObjectSetOptionsPrefix((PetscObject)(dms[i]), prefix);
467:           KSPSetDM(ilink->ksp, dms[i]);
468:           KSPSetDMActive(ilink->ksp, PETSC_FALSE);
469:           {
470:             PetscErrorCode (*func)(KSP,Mat,Mat,void*);
471:             void            *ctx;

473:             DMKSPGetComputeOperators(pc->dm, &func, &ctx);
474:             DMKSPSetComputeOperators(dms[i],  func,  ctx);
475:           }
476:           PetscObjectIncrementTabLevel((PetscObject)dms[i],(PetscObject)ilink->ksp,0);
477:           DMDestroy(&dms[i]);
478:         }
479:         PetscFree(dms);
480:       }
481:     } else {
482:       if (jac->bs <= 0) {
483:         if (pc->pmat) {
484:           MatGetBlockSize(pc->pmat,&jac->bs);
485:         } else jac->bs = 1;
486:       }

488:       if (jac->detect) {
489:         IS       zerodiags,rest;
490:         PetscInt nmin,nmax;

492:         MatGetOwnershipRange(pc->mat,&nmin,&nmax);
493:         MatFindZeroDiagonals(pc->mat,&zerodiags);
494:         ISComplement(zerodiags,nmin,nmax,&rest);
495:         PCFieldSplitSetIS(pc,"0",rest);
496:         PCFieldSplitSetIS(pc,"1",zerodiags);
497:         ISDestroy(&zerodiags);
498:         ISDestroy(&rest);
499:       } else if (coupling) {
500:         IS       coupling,rest;
501:         PetscInt nmin,nmax;

503:         MatGetOwnershipRange(pc->mat,&nmin,&nmax);
504:         MatFindOffBlockDiagonalEntries(pc->mat,&coupling);
505:         ISCreateStride(PetscObjectComm((PetscObject)pc->mat),nmax-nmin,nmin,1,&rest);
506:         ISSetIdentity(rest);
507:         PCFieldSplitSetIS(pc,"0",rest);
508:         PCFieldSplitSetIS(pc,"1",coupling);
509:         ISDestroy(&coupling);
510:         ISDestroy(&rest);
511:       } else {
512:         PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&fieldsplit_default,NULL);
513:         if (!fieldsplit_default) {
514:           /* Allow user to set fields from command line,  if bs was known at the time of PCSetFromOptions_FieldSplit()
515:            then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */
516:           PCFieldSplitSetRuntimeSplits_Private(pc);
517:           if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
518:         }
519:         if ((fieldsplit_default || !jac->splitdefined) && !jac->isrestrict) {
520:           Mat       M = pc->pmat;
521:           PetscBool isnest;

523:           PetscInfo(pc,"Using default splitting of fields\n");
524:           PetscObjectTypeCompare((PetscObject)pc->pmat,MATNEST,&isnest);
525:           if (!isnest) {
526:             M    = pc->mat;
527:             PetscObjectTypeCompare((PetscObject)pc->mat,MATNEST,&isnest);
528:           }
529:           if (isnest) {
530:             IS       *fields;
531:             PetscInt nf;

533:             MatNestGetSize(M,&nf,NULL);
534:             PetscMalloc1(nf,&fields);
535:             MatNestGetISs(M,fields,NULL);
536:             for (i=0;i<nf;i++) {
537:               PCFieldSplitSetIS(pc,NULL,fields[i]);
538:             }
539:             PetscFree(fields);
540:           } else {
541:             for (i=0; i<jac->bs; i++) {
542:               char splitname[8];
543:               PetscSNPrintf(splitname,sizeof(splitname),"%D",i);
544:               PCFieldSplitSetFields(pc,splitname,1,&i,&i);
545:             }
546:             jac->defaultsplit = PETSC_TRUE;
547:           }
548:         }
549:       }
550:     }
551:   } else if (jac->nsplits == 1) {
552:     if (ilink->is) {
553:       IS       is2;
554:       PetscInt nmin,nmax;

556:       MatGetOwnershipRange(pc->mat,&nmin,&nmax);
557:       ISComplement(ilink->is,nmin,nmax,&is2);
558:       PCFieldSplitSetIS(pc,"1",is2);
559:       ISDestroy(&is2);
560:     } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()");
561:   }

563:   if (jac->nsplits < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unhandled case, must have at least two fields, not %d", jac->nsplits);
564:   return(0);
565: }

567: static PetscErrorCode MatGolubKahanComputeExplicitOperator(Mat A,Mat B,Mat C,Mat *H,PetscReal gkbnu)
568: {
569:   PetscErrorCode    ierr;
570:   Mat               BT,T;
571:   PetscReal         nrmT,nrmB;

574:   MatHermitianTranspose(C,MAT_INITIAL_MATRIX,&T);            /* Test if augmented matrix is symmetric */
575:   MatAXPY(T,-1.0,B,DIFFERENT_NONZERO_PATTERN);
576:   MatNorm(T,NORM_1,&nrmT);
577:   MatNorm(B,NORM_1,&nrmB);
578:   if (nrmB > 0) {
579:     if (nrmT/nrmB >= PETSC_SMALL) {
580:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix is not symmetric/hermitian, GKB is not applicable.");
581:     }
582:   }
583:   /* Compute augmented Lagrangian matrix H = A00 + nu*A01*A01'. This corresponds to */
584:   /* setting N := 1/nu*I in [Ar13].                                                 */
585:   MatHermitianTranspose(B,MAT_INITIAL_MATRIX,&BT);
586:   MatMatMult(B,BT,MAT_INITIAL_MATRIX,PETSC_DEFAULT,H);       /* H = A01*A01'          */
587:   MatAYPX(*H,gkbnu,A,DIFFERENT_NONZERO_PATTERN);             /* H = A00 + nu*A01*A01' */

589:   MatDestroy(&BT);
590:   MatDestroy(&T);
591:   return(0);
592: }

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

596: static PetscErrorCode PCSetUp_FieldSplit(PC pc)
597: {
598:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
599:   PetscErrorCode    ierr;
600:   PC_FieldSplitLink ilink;
601:   PetscInt          i,nsplit;
602:   PetscBool         sorted, sorted_col;

605:   pc->failedreason = PC_NOERROR;
606:   PCFieldSplitSetDefaults(pc);
607:   nsplit = jac->nsplits;
608:   ilink  = jac->head;

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

614:     jac->issetup = PETSC_TRUE;

616:     /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */
617:     if (jac->defaultsplit || !ilink->is) {
618:       if (jac->bs <= 0) jac->bs = nsplit;
619:     }
620:     bs     = jac->bs;
621:     MatGetOwnershipRange(pc->pmat,&rstart,&rend);
622:     nslots = (rend - rstart)/bs;
623:     for (i=0; i<nsplit; i++) {
624:       if (jac->defaultsplit) {
625:         ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+i,nsplit,&ilink->is);
626:         ISDuplicate(ilink->is,&ilink->is_col);
627:       } else if (!ilink->is) {
628:         if (ilink->nfields > 1) {
629:           PetscInt *ii,*jj,j,k,nfields = ilink->nfields,*fields = ilink->fields,*fields_col = ilink->fields_col;
630:           PetscMalloc1(ilink->nfields*nslots,&ii);
631:           PetscMalloc1(ilink->nfields*nslots,&jj);
632:           for (j=0; j<nslots; j++) {
633:             for (k=0; k<nfields; k++) {
634:               ii[nfields*j + k] = rstart + bs*j + fields[k];
635:               jj[nfields*j + k] = rstart + bs*j + fields_col[k];
636:             }
637:           }
638:           ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);
639:           ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,jj,PETSC_OWN_POINTER,&ilink->is_col);
640:           ISSetBlockSize(ilink->is, nfields);
641:           ISSetBlockSize(ilink->is_col, nfields);
642:         } else {
643:           ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields[0],bs,&ilink->is);
644:           ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields_col[0],bs,&ilink->is_col);
645:         }
646:       }
647:       ISSorted(ilink->is,&sorted);
648:       if (ilink->is_col) { ISSorted(ilink->is_col,&sorted_col); }
649:       if (!sorted || !sorted_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split");
650:       ilink = ilink->next;
651:     }
652:   }

654:   ilink = jac->head;
655:   if (!jac->pmat) {
656:     Vec xtmp;

658:     MatCreateVecs(pc->pmat,&xtmp,NULL);
659:     PetscMalloc1(nsplit,&jac->pmat);
660:     PetscMalloc2(nsplit,&jac->x,nsplit,&jac->y);
661:     for (i=0; i<nsplit; i++) {
662:       MatNullSpace sp;

664:       /* Check for preconditioning matrix attached to IS */
665:       PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &jac->pmat[i]);
666:       if (jac->pmat[i]) {
667:         PetscObjectReference((PetscObject) jac->pmat[i]);
668:         if (jac->type == PC_COMPOSITE_SCHUR) {
669:           jac->schur_user = jac->pmat[i];

671:           PetscObjectReference((PetscObject) jac->schur_user);
672:         }
673:       } else {
674:         const char *prefix;
675:         MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->pmat[i]);
676:         KSPGetOptionsPrefix(ilink->ksp,&prefix);
677:         MatSetOptionsPrefix(jac->pmat[i],prefix);
678:         MatViewFromOptions(jac->pmat[i],NULL,"-mat_view");
679:       }
680:       /* create work vectors for each split */
681:       MatCreateVecs(jac->pmat[i],&jac->x[i],&jac->y[i]);
682:       ilink->x = jac->x[i]; ilink->y = jac->y[i]; ilink->z = NULL;
683:       /* compute scatter contexts needed by multiplicative versions and non-default splits */
684:       VecScatterCreate(xtmp,ilink->is,jac->x[i],NULL,&ilink->sctx);
685:       PetscObjectQuery((PetscObject) ilink->is, "nearnullspace", (PetscObject*) &sp);
686:       if (sp) {
687:         MatSetNearNullSpace(jac->pmat[i], sp);
688:       }
689:       ilink = ilink->next;
690:     }
691:     VecDestroy(&xtmp);
692:   } else {
693:     MatReuse scall;
694:     if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
695:       for (i=0; i<nsplit; i++) {
696:         MatDestroy(&jac->pmat[i]);
697:       }
698:       scall = MAT_INITIAL_MATRIX;
699:     } else scall = MAT_REUSE_MATRIX;

701:     for (i=0; i<nsplit; i++) {
702:       Mat pmat;

704:       /* Check for preconditioning matrix attached to IS */
705:       PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &pmat);
706:       if (!pmat) {
707:         MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,scall,&jac->pmat[i]);
708:       }
709:       ilink = ilink->next;
710:     }
711:   }
712:   if (jac->diag_use_amat) {
713:     ilink = jac->head;
714:     if (!jac->mat) {
715:       PetscMalloc1(nsplit,&jac->mat);
716:       for (i=0; i<nsplit; i++) {
717:         MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i]);
718:         ilink = ilink->next;
719:       }
720:     } else {
721:       MatReuse scall;
722:       if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
723:         for (i=0; i<nsplit; i++) {
724:           MatDestroy(&jac->mat[i]);
725:         }
726:         scall = MAT_INITIAL_MATRIX;
727:       } else scall = MAT_REUSE_MATRIX;

729:       for (i=0; i<nsplit; i++) {
730:         MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,scall,&jac->mat[i]);
731:         ilink = ilink->next;
732:       }
733:     }
734:   } else {
735:     jac->mat = jac->pmat;
736:   }

738:   /* Check for null space attached to IS */
739:   ilink = jac->head;
740:   for (i=0; i<nsplit; i++) {
741:     MatNullSpace sp;

743:     PetscObjectQuery((PetscObject) ilink->is, "nullspace", (PetscObject*) &sp);
744:     if (sp) {
745:       MatSetNullSpace(jac->mat[i], sp);
746:     }
747:     ilink = ilink->next;
748:   }

750:   if (jac->type != PC_COMPOSITE_ADDITIVE  && jac->type != PC_COMPOSITE_SCHUR && jac->type != PC_COMPOSITE_GKB) {
751:     /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */
752:     /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */
753:     ilink = jac->head;
754:     if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) {
755:       /* 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 */
756:       if (!jac->Afield) {
757:         PetscCalloc1(nsplit,&jac->Afield);
758:         if (jac->offdiag_use_amat) {
759:           MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);
760:         } else {
761:           MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);
762:         }
763:       } else {
764:         MatReuse scall;

766:         if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
767:           MatDestroy(&jac->Afield[1]);
768:           scall = MAT_INITIAL_MATRIX;
769:         } else scall = MAT_REUSE_MATRIX;

771:         if (jac->offdiag_use_amat) {
772:           MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,scall,&jac->Afield[1]);
773:         } else {
774:           MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,scall,&jac->Afield[1]);
775:         }
776:       }
777:     } else {
778:       if (!jac->Afield) {
779:         PetscMalloc1(nsplit,&jac->Afield);
780:         for (i=0; i<nsplit; i++) {
781:           if (jac->offdiag_use_amat) {
782:             MatCreateSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
783:           } else {
784:             MatCreateSubMatrix(pc->pmat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
785:           }
786:           ilink = ilink->next;
787:         }
788:       } else {
789:         MatReuse scall;
790:         if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
791:           for (i=0; i<nsplit; i++) {
792:             MatDestroy(&jac->Afield[i]);
793:           }
794:           scall = MAT_INITIAL_MATRIX;
795:         } else scall = MAT_REUSE_MATRIX;

797:         for (i=0; i<nsplit; i++) {
798:           if (jac->offdiag_use_amat) {
799:             MatCreateSubMatrix(pc->mat,ilink->is,NULL,scall,&jac->Afield[i]);
800:           } else {
801:             MatCreateSubMatrix(pc->pmat,ilink->is,NULL,scall,&jac->Afield[i]);
802:           }
803:           ilink = ilink->next;
804:         }
805:       }
806:     }
807:   }

809:   if (jac->type == PC_COMPOSITE_SCHUR) {
810:     IS          ccis;
811:     PetscBool   isspd;
812:     PetscInt    rstart,rend;
813:     char        lscname[256];
814:     PetscObject LSC_L;

816:     if (nsplit != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields");

818:     /* If pc->mat is SPD, don't scale by -1 the Schur complement */
819:     if (jac->schurscale == (PetscScalar)-1.0) {
820:       MatGetOption(pc->pmat,MAT_SPD,&isspd);
821:       jac->schurscale = (isspd == PETSC_TRUE) ? 1.0 : -1.0;
822:     }

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

827:     if (jac->schur) {
828:       KSP      kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper;
829:       MatReuse scall;

831:       if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
832:         scall = MAT_INITIAL_MATRIX;
833:         MatDestroy(&jac->B);
834:         MatDestroy(&jac->C);
835:       } else scall = MAT_REUSE_MATRIX;

837:       MatSchurComplementGetKSP(jac->schur, &kspInner);
838:       ilink = jac->head;
839:       ISComplement(ilink->is_col,rstart,rend,&ccis);
840:       if (jac->offdiag_use_amat) {
841:         MatCreateSubMatrix(pc->mat,ilink->is,ccis,scall,&jac->B);
842:       } else {
843:         MatCreateSubMatrix(pc->pmat,ilink->is,ccis,scall,&jac->B);
844:       }
845:       ISDestroy(&ccis);
846:       ilink = ilink->next;
847:       ISComplement(ilink->is_col,rstart,rend,&ccis);
848:       if (jac->offdiag_use_amat) {
849:         MatCreateSubMatrix(pc->mat,ilink->is,ccis,scall,&jac->C);
850:       } else {
851:         MatCreateSubMatrix(pc->pmat,ilink->is,ccis,scall,&jac->C);
852:       }
853:       ISDestroy(&ccis);
854:       MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
855:       if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
856:         MatDestroy(&jac->schurp);
857:         MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);
858:       }
859:       if (kspA != kspInner) {
860:         KSPSetOperators(kspA,jac->mat[0],jac->pmat[0]);
861:       }
862:       if (kspUpper != kspA) {
863:         KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0]);
864:       }
865:       KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
866:     } else {
867:       const char   *Dprefix;
868:       char         schurprefix[256], schurmatprefix[256];
869:       char         schurtestoption[256];
870:       MatNullSpace sp;
871:       PetscBool    flg;
872:       KSP          kspt;

874:       /* extract the A01 and A10 matrices */
875:       ilink = jac->head;
876:       ISComplement(ilink->is_col,rstart,rend,&ccis);
877:       if (jac->offdiag_use_amat) {
878:         MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
879:       } else {
880:         MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
881:       }
882:       ISDestroy(&ccis);
883:       ilink = ilink->next;
884:       ISComplement(ilink->is_col,rstart,rend,&ccis);
885:       if (jac->offdiag_use_amat) {
886:         MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
887:       } else {
888:         MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
889:       }
890:       ISDestroy(&ccis);

892:       /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
893:       MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);
894:       MatSetType(jac->schur,MATSCHURCOMPLEMENT);
895:       MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
896:       PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
897:       MatSetOptionsPrefix(jac->schur,schurmatprefix);
898:       MatSchurComplementGetKSP(jac->schur,&kspt);
899:       KSPSetOptionsPrefix(kspt,schurmatprefix);

901:       /* Note: this is not true in general */
902:       MatGetNullSpace(jac->mat[1], &sp);
903:       if (sp) {
904:         MatSetNullSpace(jac->schur, sp);
905:       }

907:       PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);
908:       PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
909:       if (flg) {
910:         DM  dmInner;
911:         KSP kspInner;
912:         PC  pcInner;

914:         MatSchurComplementGetKSP(jac->schur, &kspInner);
915:         KSPReset(kspInner);
916:         KSPSetOperators(kspInner,jac->mat[0],jac->pmat[0]);
917:         PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
918:         /* Indent this deeper to emphasize the "inner" nature of this solver. */
919:         PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject) pc, 2);
920:         PetscObjectIncrementTabLevel((PetscObject)kspInner->pc, (PetscObject) pc, 2);
921:         KSPSetOptionsPrefix(kspInner, schurprefix);

923:         /* Set DM for new solver */
924:         KSPGetDM(jac->head->ksp, &dmInner);
925:         KSPSetDM(kspInner, dmInner);
926:         KSPSetDMActive(kspInner, PETSC_FALSE);

928:         /* Defaults to PCKSP as preconditioner */
929:         KSPGetPC(kspInner, &pcInner);
930:         PCSetType(pcInner, PCKSP);
931:         PCKSPSetKSP(pcInner, jac->head->ksp);
932:       } else {
933:          /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or
934:           * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact,
935:           * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for
936:           * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make
937:           * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used
938:           * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */
939:         KSPSetType(jac->head->ksp,KSPGMRES);
940:         MatSchurComplementSetKSP(jac->schur,jac->head->ksp);
941:       }
942:       KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0]);
943:       KSPSetFromOptions(jac->head->ksp);
944:       MatSetFromOptions(jac->schur);

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

951:         MatSchurComplementGetKSP(jac->schur, &kspInner);
952:         KSPGetPC(kspInner, &pcInner);
953:         PetscObjectTypeCompare((PetscObject)pcInner, PCKSP, &flg);
954:         if (flg) {
955:           KSP ksp;

957:           PCKSPGetKSP(pcInner, &ksp);
958:           if (ksp == jac->head->ksp) {
959:             PCSetUseAmat(pcInner, PETSC_TRUE);
960:           }
961:         }
962:       }
963:       PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);
964:       PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
965:       if (flg) {
966:         DM dmInner;

968:         PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
969:         KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);
970:         KSPSetErrorIfNotConverged(jac->kspupper,pc->erroriffailure);
971:         KSPSetOptionsPrefix(jac->kspupper, schurprefix);
972:         PetscObjectIncrementTabLevel((PetscObject)jac->kspupper, (PetscObject) pc, 1);
973:         PetscObjectIncrementTabLevel((PetscObject)jac->kspupper->pc, (PetscObject) pc, 1);
974:         KSPGetDM(jac->head->ksp, &dmInner);
975:         KSPSetDM(jac->kspupper, dmInner);
976:         KSPSetDMActive(jac->kspupper, PETSC_FALSE);
977:         KSPSetFromOptions(jac->kspupper);
978:         KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0]);
979:         VecDuplicate(jac->head->x, &jac->head->z);
980:       } else {
981:         jac->kspupper = jac->head->ksp;
982:         PetscObjectReference((PetscObject) jac->head->ksp);
983:       }

985:       if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
986:         MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);
987:       }
988:       KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);
989:       KSPSetErrorIfNotConverged(jac->kspschur,pc->erroriffailure);
990:       PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);
991:       PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);
992:       if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
993:         PC pcschur;
994:         KSPGetPC(jac->kspschur,&pcschur);
995:         PCSetType(pcschur,PCNONE);
996:         /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
997:       } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) {
998:         MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user);
999:       }
1000:       KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
1001:       KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);
1002:       KSPSetOptionsPrefix(jac->kspschur,         Dprefix);
1003:       /* propagate DM */
1004:       {
1005:         DM sdm;
1006:         KSPGetDM(jac->head->next->ksp, &sdm);
1007:         if (sdm) {
1008:           KSPSetDM(jac->kspschur, sdm);
1009:           KSPSetDMActive(jac->kspschur, PETSC_FALSE);
1010:         }
1011:       }
1012:       /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
1013:       /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */
1014:       KSPSetFromOptions(jac->kspschur);
1015:     }
1016:     MatAssemblyBegin(jac->schur,MAT_FINAL_ASSEMBLY);
1017:     MatAssemblyEnd(jac->schur,MAT_FINAL_ASSEMBLY);

1019:     /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */
1020:     PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);
1021:     PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);
1022:     if (!LSC_L) {PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);}
1023:     if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);}
1024:     PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);
1025:     PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);
1026:     if (!LSC_L) {PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);}
1027:     if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);}
1028:   } else if (jac->type == PC_COMPOSITE_GKB) {
1029:     IS          ccis;
1030:     PetscInt    rstart,rend;

1032:     if (nsplit != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use GKB preconditioner you must have exactly 2 fields");

1034:     ilink = jac->head;

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

1039:     ISComplement(ilink->is_col,rstart,rend,&ccis);
1040:     if (jac->offdiag_use_amat) {
1041:       MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
1042:     } else {
1043:       MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
1044:     }
1045:     ISDestroy(&ccis);
1046:     /* Create work vectors for GKB algorithm */
1047:     VecDuplicate(ilink->x,&jac->u);
1048:     VecDuplicate(ilink->x,&jac->Hu);
1049:     VecDuplicate(ilink->x,&jac->w2);
1050:     ilink = ilink->next;
1051:     ISComplement(ilink->is_col,rstart,rend,&ccis);
1052:     if (jac->offdiag_use_amat) {
1053:       MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
1054:     } else {
1055:       MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
1056:     }
1057:     ISDestroy(&ccis);
1058:     /* Create work vectors for GKB algorithm */
1059:     VecDuplicate(ilink->x,&jac->v);
1060:     VecDuplicate(ilink->x,&jac->d);
1061:     VecDuplicate(ilink->x,&jac->w1);
1062:     MatGolubKahanComputeExplicitOperator(jac->mat[0],jac->B,jac->C,&jac->H,jac->gkbnu);
1063:     PetscCalloc1(jac->gkbdelay,&jac->vecz);

1065:     ilink = jac->head;
1066:     KSPSetOperators(ilink->ksp,jac->H,jac->H);
1067:     if (!jac->suboptionsset) {KSPSetFromOptions(ilink->ksp);}
1068:     /* Create gkb_monitor context */
1069:     if (jac->gkbmonitor) {
1070:       PetscInt  tablevel;
1071:       PetscViewerCreate(PETSC_COMM_WORLD,&jac->gkbviewer);
1072:       PetscViewerSetType(jac->gkbviewer,PETSCVIEWERASCII);
1073:       PetscObjectGetTabLevel((PetscObject)ilink->ksp,&tablevel);
1074:       PetscViewerASCIISetTab(jac->gkbviewer,tablevel);
1075:       PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)ilink->ksp,1);
1076:     }
1077:   } else {
1078:     /* set up the individual splits' PCs */
1079:     i     = 0;
1080:     ilink = jac->head;
1081:     while (ilink) {
1082:       KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i]);
1083:       /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
1084:       if (!jac->suboptionsset) {KSPSetFromOptions(ilink->ksp);}
1085:       i++;
1086:       ilink = ilink->next;
1087:     }
1088:   }

1090:   jac->suboptionsset = PETSC_TRUE;
1091:   return(0);
1092: }

1094: #define FieldSplitSplitSolveAdd(ilink,xx,yy) \
1095:   (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
1096:    VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
1097:    PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\
1098:    KSPSolve(ilink->ksp,ilink->x,ilink->y) ||                               \
1099:    KSPCheckSolve(ilink->ksp,pc,ilink->y)  || \
1100:    PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\
1101:    VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) ||  \
1102:    VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE))

1104: static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y)
1105: {
1106:   PC_FieldSplit      *jac = (PC_FieldSplit*)pc->data;
1107:   PetscErrorCode     ierr;
1108:   PC_FieldSplitLink  ilinkA = jac->head, ilinkD = ilinkA->next;
1109:   KSP                kspA   = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;

1112:   switch (jac->schurfactorization) {
1113:   case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
1114:     /* [A00 0; 0 -S], positive definite, suitable for MINRES */
1115:     VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1116:     VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1117:     VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1118:     PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1119:     KSPSolve(kspA,ilinkA->x,ilinkA->y);
1120:     KSPCheckSolve(kspA,pc,ilinkA->y);
1121:     PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1122:     VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1123:     VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1124:     PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1125:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
1126:     KSPCheckSolve(jac->kspschur,pc,ilinkD->y);
1127:     PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1128:     VecScale(ilinkD->y,jac->schurscale);
1129:     VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1130:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1131:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1132:     break;
1133:   case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
1134:     /* [A00 0; A10 S], suitable for left preconditioning */
1135:     VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1136:     VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1137:     PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1138:     KSPSolve(kspA,ilinkA->x,ilinkA->y);
1139:     KSPCheckSolve(kspA,pc,ilinkA->y);
1140:     PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1141:     MatMult(jac->C,ilinkA->y,ilinkD->x);
1142:     VecScale(ilinkD->x,-1.);
1143:     VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
1144:     VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1145:     VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
1146:     PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1147:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
1148:     KSPCheckSolve(jac->kspschur,pc,ilinkD->y);
1149:     PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1150:     VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1151:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1152:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1153:     break;
1154:   case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
1155:     /* [A00 A01; 0 S], suitable for right preconditioning */
1156:     VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1157:     VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1158:     PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1159:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
1160:     KSPCheckSolve(jac->kspschur,pc,ilinkD->y);
1161:     PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);    MatMult(jac->B,ilinkD->y,ilinkA->x);
1162:     VecScale(ilinkA->x,-1.);
1163:     VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
1164:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1165:     VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
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:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1171:     VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1172:     VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1173:     break;
1174:   case PC_FIELDSPLIT_SCHUR_FACT_FULL:
1175:     /* [1 0; A10 A00^{-1} 1] [A00 0; 0 S] [1 A00^{-1}A01; 0 1] */
1176:     VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1177:     VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1178:     PetscLogEventBegin(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);
1179:     KSPSolve(kspLower,ilinkA->x,ilinkA->y);
1180:     KSPCheckSolve(kspLower,pc,ilinkA->y);
1181:     PetscLogEventEnd(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);
1182:     MatMult(jac->C,ilinkA->y,ilinkD->x);
1183:     VecScale(ilinkD->x,-1.0);
1184:     VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
1185:     VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);

1187:     PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1188:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
1189:     KSPCheckSolve(jac->kspschur,pc,ilinkD->y);
1190:     PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1191:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);

1193:     if (kspUpper == kspA) {
1194:       MatMult(jac->B,ilinkD->y,ilinkA->y);
1195:       VecAXPY(ilinkA->x,-1.0,ilinkA->y);
1196:       PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1197:       KSPSolve(kspA,ilinkA->x,ilinkA->y);
1198:       KSPCheckSolve(kspA,pc,ilinkA->y);
1199:       PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1200:     } else {
1201:       PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1202:       KSPSolve(kspA,ilinkA->x,ilinkA->y);
1203:       KSPCheckSolve(kspA,pc,ilinkA->y);
1204:       MatMult(jac->B,ilinkD->y,ilinkA->x);
1205:       PetscLogEventBegin(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);
1206:       KSPSolve(kspUpper,ilinkA->x,ilinkA->z);
1207:       KSPCheckSolve(kspUpper,pc,ilinkA->z);
1208:       PetscLogEventEnd(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);
1209:       VecAXPY(ilinkA->y,-1.0,ilinkA->z);
1210:     }
1211:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1212:     VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1213:     VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1214:   }
1215:   return(0);
1216: }

1218: static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y)
1219: {
1220:   PC_FieldSplit      *jac = (PC_FieldSplit*)pc->data;
1221:   PetscErrorCode     ierr;
1222:   PC_FieldSplitLink  ilink = jac->head;
1223:   PetscInt           cnt,bs;

1226:   if (jac->type == PC_COMPOSITE_ADDITIVE) {
1227:     if (jac->defaultsplit) {
1228:       VecGetBlockSize(x,&bs);
1229:       if (jac->bs > 0 && bs != jac->bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of x vector %D does not match fieldsplit blocksize %D",bs,jac->bs);
1230:       VecGetBlockSize(y,&bs);
1231:       if (jac->bs > 0 && bs != jac->bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of y vector %D does not match fieldsplit blocksize %D",bs,jac->bs);
1232:       VecStrideGatherAll(x,jac->x,INSERT_VALUES);
1233:       while (ilink) {
1234:         PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1235:         KSPSolve(ilink->ksp,ilink->x,ilink->y);
1236:         KSPCheckSolve(ilink->ksp,pc,ilink->y);
1237:         PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1238:         ilink = ilink->next;
1239:       }
1240:       VecStrideScatterAll(jac->y,y,INSERT_VALUES);
1241:     } else {
1242:       VecSet(y,0.0);
1243:       while (ilink) {
1244:         FieldSplitSplitSolveAdd(ilink,x,y);
1245:         ilink = ilink->next;
1246:       }
1247:     }
1248:   } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) {
1249:     VecSet(y,0.0);
1250:     /* solve on first block for first block variables */
1251:     VecScatterBegin(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
1252:     VecScatterEnd(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
1253:     PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1254:     KSPSolve(ilink->ksp,ilink->x,ilink->y);
1255:     KSPCheckSolve(ilink->ksp,pc,ilink->y);
1256:     PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1257:     VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1258:     VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);

1260:     /* compute the residual only onto second block variables using first block variables */
1261:     MatMult(jac->Afield[1],ilink->y,ilink->next->x);
1262:     ilink = ilink->next;
1263:     VecScale(ilink->x,-1.0);
1264:     VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1265:     VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);

1267:     /* solve on second block variables */
1268:     PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1269:     KSPSolve(ilink->ksp,ilink->x,ilink->y);
1270:     KSPCheckSolve(ilink->ksp,pc,ilink->y);
1271:     PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1272:     VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1273:     VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1274:   } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1275:     if (!jac->w1) {
1276:       VecDuplicate(x,&jac->w1);
1277:       VecDuplicate(x,&jac->w2);
1278:     }
1279:     VecSet(y,0.0);
1280:     FieldSplitSplitSolveAdd(ilink,x,y);
1281:     cnt  = 1;
1282:     while (ilink->next) {
1283:       ilink = ilink->next;
1284:       /* compute the residual only over the part of the vector needed */
1285:       MatMult(jac->Afield[cnt++],y,ilink->x);
1286:       VecScale(ilink->x,-1.0);
1287:       VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1288:       VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1289:       PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1290:       KSPSolve(ilink->ksp,ilink->x,ilink->y);
1291:       KSPCheckSolve(ilink->ksp,pc,ilink->y);
1292:       PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1293:       VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1294:       VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1295:     }
1296:     if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1297:       cnt -= 2;
1298:       while (ilink->previous) {
1299:         ilink = ilink->previous;
1300:         /* compute the residual only over the part of the vector needed */
1301:         MatMult(jac->Afield[cnt--],y,ilink->x);
1302:         VecScale(ilink->x,-1.0);
1303:         VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1304:         VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1305:         PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1306:         KSPSolve(ilink->ksp,ilink->x,ilink->y);
1307:         KSPCheckSolve(ilink->ksp,pc,ilink->y);
1308:         PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1309:         VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1310:         VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1311:       }
1312:     }
1313:   } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type);
1314:   return(0);
1315: }


1318: static PetscErrorCode PCApply_FieldSplit_GKB(PC pc,Vec x,Vec y)
1319: {
1320:   PC_FieldSplit      *jac = (PC_FieldSplit*)pc->data;
1321:   PetscErrorCode     ierr;
1322:   PC_FieldSplitLink  ilinkA = jac->head,ilinkD = ilinkA->next;
1323:   KSP                ksp = ilinkA->ksp;
1324:   Vec                u,v,Hu,d,work1,work2;
1325:   PetscScalar        alpha,z,nrmz2,*vecz;
1326:   PetscReal          lowbnd,nu,beta;
1327:   PetscInt           j,iterGKB;

1330:   VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1331:   VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
1332:   VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
1333:   VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);

1335:   u     = jac->u;
1336:   v     = jac->v;
1337:   Hu    = jac->Hu;
1338:   d     = jac->d;
1339:   work1 = jac->w1;
1340:   work2 = jac->w2;
1341:   vecz  = jac->vecz;

1343:   /* Change RHS to comply with matrix regularization H = A + nu*B*B' */
1344:   /* Add q = q + nu*B*b */
1345:   if (jac->gkbnu) {
1346:     nu = jac->gkbnu;
1347:     VecScale(ilinkD->x,jac->gkbnu);
1348:     MatMultAdd(jac->B,ilinkD->x,ilinkA->x,ilinkA->x);            /* q = q + nu*B*b */
1349:   } else {
1350:     /* Situation when no augmented Lagrangian is used. Then we set inner  */
1351:     /* matrix N = I in [Ar13], and thus nu = 1.                           */
1352:     nu = 1;
1353:   }

1355:   /* Transform rhs from [q,tilde{b}] to [0,b] */
1356:   PetscLogEventBegin(ilinkA->event,ksp,ilinkA->x,ilinkA->y,NULL);
1357:   KSPSolve(ksp,ilinkA->x,ilinkA->y);
1358:   KSPCheckSolve(ksp,pc,ilinkA->y);
1359:   PetscLogEventEnd(ilinkA->event,ksp,ilinkA->x,ilinkA->y,NULL);
1360:   MatMultHermitianTranspose(jac->B,ilinkA->y,work1);
1361:   VecAXPBY(work1,1.0/nu,-1.0,ilinkD->x);            /* c = b - B'*x        */

1363:   /* First step of algorithm */
1364:   VecNorm(work1,NORM_2,&beta);                   /* beta = sqrt(nu*c'*c)*/
1365:   KSPCheckDot(ksp,beta);
1366:   beta  = PetscSqrtScalar(nu)*beta;
1367:   VecAXPBY(v,nu/beta,0.0,work1);                   /* v = nu/beta *c      */
1368:   MatMult(jac->B,v,work2);                       /* u = H^{-1}*B*v      */
1369:   PetscLogEventBegin(ilinkA->event,ksp,work2,u,NULL);
1370:   KSPSolve(ksp,work2,u);
1371:   KSPCheckSolve(ksp,pc,u);
1372:   PetscLogEventEnd(ilinkA->event,ksp,work2,u,NULL);
1373:   MatMult(jac->H,u,Hu);                          /* alpha = u'*H*u      */
1374:   VecDot(Hu,u,&alpha);
1375:   KSPCheckDot(ksp,alpha);
1376:   if (PetscRealPart(alpha) <= 0.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_NOT_CONVERGED,"GKB preconditioner diverged, H is not positive definite");
1377:   alpha = PetscSqrtScalar(PetscAbsScalar(alpha));
1378:   VecScale(u,1.0/alpha);
1379:   VecAXPBY(d,1.0/alpha,0.0,v);                       /* v = nu/beta *c      */

1381:   z = beta/alpha;
1382:   vecz[1] = z;

1384:   /* Computation of first iterate x(1) and p(1) */
1385:   VecAXPY(ilinkA->y,z,u);
1386:   VecCopy(d,ilinkD->y);
1387:   VecScale(ilinkD->y,-z);

1389:   iterGKB = 1; lowbnd = 2*jac->gkbtol;
1390:   if (jac->gkbmonitor) {
1391:       PetscViewerASCIIPrintf(jac->gkbviewer,"%3D GKB Lower bound estimate %14.12e\n",iterGKB,lowbnd);
1392:   }

1394:   while (iterGKB < jac->gkbmaxit && lowbnd > jac->gkbtol) {
1395:     iterGKB += 1;
1396:     MatMultHermitianTranspose(jac->B,u,work1); /* v <- nu*(B'*u-alpha/nu*v) */
1397:     VecAXPBY(v,nu,-alpha,work1);
1398:     VecNorm(v,NORM_2,&beta);                   /* beta = sqrt(nu)*v'*v      */
1399:     beta  = beta/PetscSqrtScalar(nu);
1400:     VecScale(v,1.0/beta);
1401:     MatMult(jac->B,v,work2);                  /* u <- H^{-1}*(B*v-beta*H*u) */
1402:     MatMult(jac->H,u,Hu);
1403:     VecAXPY(work2,-beta,Hu);
1404:     PetscLogEventBegin(ilinkA->event,ksp,work2,u,NULL);
1405:     KSPSolve(ksp,work2,u);
1406:     KSPCheckSolve(ksp,pc,u);
1407:     PetscLogEventEnd(ilinkA->event,ksp,work2,u,NULL);
1408:     MatMult(jac->H,u,Hu);                      /* alpha = u'*H*u            */
1409:     VecDot(Hu,u,&alpha);
1410:     KSPCheckDot(ksp,alpha);
1411:     if (PetscRealPart(alpha) <= 0.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_NOT_CONVERGED,"GKB preconditioner diverged, H is not positive definite");
1412:     alpha = PetscSqrtScalar(PetscAbsScalar(alpha));
1413:     VecScale(u,1.0/alpha);

1415:     z = -beta/alpha*z;                                            /* z <- beta/alpha*z     */
1416:     vecz[0] = z;

1418:     /* Computation of new iterate x(i+1) and p(i+1) */
1419:     VecAXPBY(d,1.0/alpha,-beta/alpha,v);       /* d = (v-beta*d)/alpha */
1420:     VecAXPY(ilinkA->y,z,u);                  /* r = r + z*u          */
1421:     VecAXPY(ilinkD->y,-z,d);                 /* p = p - z*d          */
1422:     MatMult(jac->H,ilinkA->y,Hu);            /* ||u||_H = u'*H*u     */
1423:     VecDot(Hu,ilinkA->y,&nrmz2);

1425:     /* Compute Lower Bound estimate */
1426:     if (iterGKB > jac->gkbdelay) {
1427:       lowbnd = 0.0;
1428:       for (j=0; j<jac->gkbdelay; j++) {
1429:         lowbnd += PetscAbsScalar(vecz[j]*vecz[j]);
1430:       }
1431:       lowbnd = PetscSqrtScalar(lowbnd/PetscAbsScalar(nrmz2));
1432:     }

1434:     for (j=0; j<jac->gkbdelay-1; j++) {
1435:       vecz[jac->gkbdelay-j-1] = vecz[jac->gkbdelay-j-2];
1436:     }
1437:     if (jac->gkbmonitor) {
1438:       PetscViewerASCIIPrintf(jac->gkbviewer,"%3D GKB Lower bound estimate %14.12e\n",iterGKB,lowbnd);
1439:     }
1440:   }

1442:   VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1443:   VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1444:   VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1445:   VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);

1447:   return(0);
1448: }


1451: #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \
1452:   (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1453:    VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1454:    PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \
1455:    KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) ||                  \
1456:    KSPCheckSolve(ilink->ksp,pc,ilink->x) || \
1457:    PetscLogEventEnd(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) ||   \
1458:    VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \
1459:    VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE))

1461: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y)
1462: {
1463:   PC_FieldSplit      *jac = (PC_FieldSplit*)pc->data;
1464:   PetscErrorCode     ierr;
1465:   PC_FieldSplitLink  ilink = jac->head;
1466:   PetscInt           bs;

1469:   if (jac->type == PC_COMPOSITE_ADDITIVE) {
1470:     if (jac->defaultsplit) {
1471:       VecGetBlockSize(x,&bs);
1472:       if (jac->bs > 0 && bs != jac->bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of x vector %D does not match fieldsplit blocksize %D",bs,jac->bs);
1473:       VecGetBlockSize(y,&bs);
1474:       if (jac->bs > 0 && bs != jac->bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of y vector %D does not match fieldsplit blocksize %D",bs,jac->bs);
1475:       VecStrideGatherAll(x,jac->x,INSERT_VALUES);
1476:       while (ilink) {
1477:         PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1478:         KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);
1479:         KSPCheckSolve(ilink->ksp,pc,ilink->y);
1480:         PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1481:         ilink = ilink->next;
1482:       }
1483:       VecStrideScatterAll(jac->y,y,INSERT_VALUES);
1484:     } else {
1485:       VecSet(y,0.0);
1486:       while (ilink) {
1487:         FieldSplitSplitSolveAddTranspose(ilink,x,y);
1488:         ilink = ilink->next;
1489:       }
1490:     }
1491:   } else {
1492:     if (!jac->w1) {
1493:       VecDuplicate(x,&jac->w1);
1494:       VecDuplicate(x,&jac->w2);
1495:     }
1496:     VecSet(y,0.0);
1497:     if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1498:       FieldSplitSplitSolveAddTranspose(ilink,x,y);
1499:       while (ilink->next) {
1500:         ilink = ilink->next;
1501:         MatMultTranspose(pc->mat,y,jac->w1);
1502:         VecWAXPY(jac->w2,-1.0,jac->w1,x);
1503:         FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1504:       }
1505:       while (ilink->previous) {
1506:         ilink = ilink->previous;
1507:         MatMultTranspose(pc->mat,y,jac->w1);
1508:         VecWAXPY(jac->w2,-1.0,jac->w1,x);
1509:         FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1510:       }
1511:     } else {
1512:       while (ilink->next) {   /* get to last entry in linked list */
1513:         ilink = ilink->next;
1514:       }
1515:       FieldSplitSplitSolveAddTranspose(ilink,x,y);
1516:       while (ilink->previous) {
1517:         ilink = ilink->previous;
1518:         MatMultTranspose(pc->mat,y,jac->w1);
1519:         VecWAXPY(jac->w2,-1.0,jac->w1,x);
1520:         FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1521:       }
1522:     }
1523:   }
1524:   return(0);
1525: }

1527: static PetscErrorCode PCReset_FieldSplit(PC pc)
1528: {
1529:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1530:   PetscErrorCode    ierr;
1531:   PC_FieldSplitLink ilink = jac->head,next;

1534:   while (ilink) {
1535:     KSPDestroy(&ilink->ksp);
1536:     VecDestroy(&ilink->x);
1537:     VecDestroy(&ilink->y);
1538:     VecDestroy(&ilink->z);
1539:     VecScatterDestroy(&ilink->sctx);
1540:     ISDestroy(&ilink->is);
1541:     ISDestroy(&ilink->is_col);
1542:     PetscFree(ilink->splitname);
1543:     PetscFree(ilink->fields);
1544:     PetscFree(ilink->fields_col);
1545:     next  = ilink->next;
1546:     PetscFree(ilink);
1547:     ilink = next;
1548:   }
1549:   jac->head = NULL;
1550:   PetscFree2(jac->x,jac->y);
1551:   if (jac->mat && jac->mat != jac->pmat) {
1552:     MatDestroyMatrices(jac->nsplits,&jac->mat);
1553:   } else if (jac->mat) {
1554:     jac->mat = NULL;
1555:   }
1556:   if (jac->pmat) {MatDestroyMatrices(jac->nsplits,&jac->pmat);}
1557:   if (jac->Afield) {MatDestroyMatrices(jac->nsplits,&jac->Afield);}
1558:   jac->nsplits = 0;
1559:   VecDestroy(&jac->w1);
1560:   VecDestroy(&jac->w2);
1561:   MatDestroy(&jac->schur);
1562:   MatDestroy(&jac->schurp);
1563:   MatDestroy(&jac->schur_user);
1564:   KSPDestroy(&jac->kspschur);
1565:   KSPDestroy(&jac->kspupper);
1566:   MatDestroy(&jac->B);
1567:   MatDestroy(&jac->C);
1568:   MatDestroy(&jac->H);
1569:   VecDestroy(&jac->u);
1570:   VecDestroy(&jac->v);
1571:   VecDestroy(&jac->Hu);
1572:   VecDestroy(&jac->d);
1573:   PetscFree(jac->vecz);
1574:   PetscViewerDestroy(&jac->gkbviewer);
1575:   jac->isrestrict = PETSC_FALSE;
1576:   return(0);
1577: }

1579: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
1580: {
1581:   PetscErrorCode    ierr;

1584:   PCReset_FieldSplit(pc);
1585:   PetscFree(pc->data);
1586:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",NULL);
1587:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);
1588:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);
1589:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);
1590:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);
1591:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);
1592:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);
1593:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);
1594:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);
1595:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",NULL);
1596:   return(0);
1597: }

1599: static PetscErrorCode PCSetFromOptions_FieldSplit(PetscOptionItems *PetscOptionsObject,PC pc)
1600: {
1601:   PetscErrorCode  ierr;
1602:   PetscInt        bs;
1603:   PetscBool       flg;
1604:   PC_FieldSplit   *jac = (PC_FieldSplit*)pc->data;
1605:   PCCompositeType ctype;

1608:   PetscOptionsHead(PetscOptionsObject,"FieldSplit options");
1609:   PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);
1610:   PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);
1611:   if (flg) {
1612:     PCFieldSplitSetBlockSize(pc,bs);
1613:   }
1614:   jac->diag_use_amat = pc->useAmat;
1615:   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);
1616:   jac->offdiag_use_amat = pc->useAmat;
1617:   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);
1618:   PetscOptionsBool("-pc_fieldsplit_detect_saddle_point","Form 2-way split by detecting zero diagonal entries", "PCFieldSplitSetDetectSaddlePoint",jac->detect,&jac->detect,NULL);
1619:   PCFieldSplitSetDetectSaddlePoint(pc,jac->detect); /* Sets split type and Schur PC type */
1620:   PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);
1621:   if (flg) {
1622:     PCFieldSplitSetType(pc,ctype);
1623:   }
1624:   /* Only setup fields once */
1625:   if ((jac->bs > 0) && (jac->nsplits == 0)) {
1626:     /* only allow user to set fields from command line if bs is already known.
1627:        otherwise user can set them in PCFieldSplitSetDefaults() */
1628:     PCFieldSplitSetRuntimeSplits_Private(pc);
1629:     if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
1630:   }
1631:   if (jac->type == PC_COMPOSITE_SCHUR) {
1632:     PetscOptionsGetEnum(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);
1633:     if (flg) {PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");}
1634:     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);
1635:     PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);
1636:     PetscOptionsScalar("-pc_fieldsplit_schur_scale","Scale Schur complement","PCFieldSplitSetSchurScale",jac->schurscale,&jac->schurscale,NULL);
1637:   } else if (jac->type == PC_COMPOSITE_GKB) {
1638:     PetscOptionsReal("-pc_fieldsplit_gkb_tol","The tolerance for the lower bound stopping criterion","PCFieldSplitGKBTol",jac->gkbtol,&jac->gkbtol,NULL);
1639:     PetscOptionsInt("-pc_fieldsplit_gkb_delay","The delay value for lower bound criterion","PCFieldSplitGKBDelay",jac->gkbdelay,&jac->gkbdelay,NULL);
1640:     PetscOptionsReal("-pc_fieldsplit_gkb_nu","Parameter in augmented Lagrangian approach","PCFieldSplitGKBNu",jac->gkbnu,&jac->gkbnu,NULL);
1641:     if (jac->gkbnu < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nu cannot be less than 0: value %f",jac->gkbnu);
1642:     PetscOptionsInt("-pc_fieldsplit_gkb_maxit","Maximum allowed number of iterations","PCFieldSplitGKBMaxit",jac->gkbmaxit,&jac->gkbmaxit,NULL);
1643:     PetscOptionsBool("-pc_fieldsplit_gkb_monitor","Prints number of GKB iterations and error","PCFieldSplitGKB",jac->gkbmonitor,&jac->gkbmonitor,NULL);
1644:   }
1645:   PetscOptionsTail();
1646:   return(0);
1647: }

1649: /*------------------------------------------------------------------------------------*/

1651: static PetscErrorCode  PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1652: {
1653:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1654:   PetscErrorCode    ierr;
1655:   PC_FieldSplitLink ilink,next = jac->head;
1656:   char              prefix[128];
1657:   PetscInt          i;

1660:   if (jac->splitdefined) {
1661:     PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1662:     return(0);
1663:   }
1664:   for (i=0; i<n; i++) {
1665:     if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs);
1666:     if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]);
1667:   }
1668:   PetscNew(&ilink);
1669:   if (splitname) {
1670:     PetscStrallocpy(splitname,&ilink->splitname);
1671:   } else {
1672:     PetscMalloc1(3,&ilink->splitname);
1673:     PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);
1674:   }
1675:   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 */
1676:   PetscMalloc1(n,&ilink->fields);
1677:   PetscArraycpy(ilink->fields,fields,n);
1678:   PetscMalloc1(n,&ilink->fields_col);
1679:   PetscArraycpy(ilink->fields_col,fields_col,n);

1681:   ilink->nfields = n;
1682:   ilink->next    = NULL;
1683:   KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1684:   KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1685:   PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1686:   KSPSetType(ilink->ksp,KSPPREONLY);
1687:   PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);

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

1692:   if (!next) {
1693:     jac->head       = ilink;
1694:     ilink->previous = NULL;
1695:   } else {
1696:     while (next->next) {
1697:       next = next->next;
1698:     }
1699:     next->next      = ilink;
1700:     ilink->previous = next;
1701:   }
1702:   jac->nsplits++;
1703:   return(0);
1704: }

1706: static PetscErrorCode  PCFieldSplitSchurGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1707: {
1708:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

1712:   *subksp = NULL;
1713:   if (n) *n = 0;
1714:   if (jac->type == PC_COMPOSITE_SCHUR) {
1715:     PetscInt nn;

1717:     if (!jac->schur) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitSchurGetSubKSP()");
1718:     if (jac->nsplits != 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unexpected number of splits %D != 2",jac->nsplits);
1719:     nn   = jac->nsplits + (jac->kspupper != jac->head->ksp ? 1 : 0);
1720:     PetscMalloc1(nn,subksp);
1721:     (*subksp)[0] = jac->head->ksp;
1722:     (*subksp)[1] = jac->kspschur;
1723:     if (jac->kspupper != jac->head->ksp) (*subksp)[2] = jac->kspupper;
1724:     if (n) *n = nn;
1725:   }
1726:   return(0);
1727: }

1729: static PetscErrorCode  PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp)
1730: {
1731:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

1735:   if (!jac->schur) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitGetSubKSP()");
1736:   PetscMalloc1(jac->nsplits,subksp);
1737:   MatSchurComplementGetKSP(jac->schur,*subksp);

1739:   (*subksp)[1] = jac->kspschur;
1740:   if (n) *n = jac->nsplits;
1741:   return(0);
1742: }

1744: static PetscErrorCode  PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1745: {
1746:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1747:   PetscErrorCode    ierr;
1748:   PetscInt          cnt   = 0;
1749:   PC_FieldSplitLink ilink = jac->head;

1752:   PetscMalloc1(jac->nsplits,subksp);
1753:   while (ilink) {
1754:     (*subksp)[cnt++] = ilink->ksp;
1755:     ilink            = ilink->next;
1756:   }
1757:   if (cnt != jac->nsplits) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt PCFIELDSPLIT object: number of splits in linked list %D does not match number in object %D",cnt,jac->nsplits);
1758:   if (n) *n = jac->nsplits;
1759:   return(0);
1760: }

1762: /*@C
1763:     PCFieldSplitRestrictIS - Restricts the fieldsplit ISs to be within a given IS.

1765:     Input Parameters:
1766: +   pc  - the preconditioner context
1767: -   is - the index set that defines the indices to which the fieldsplit is to be restricted

1769:     Level: advanced

1771: @*/
1772: PetscErrorCode  PCFieldSplitRestrictIS(PC pc,IS isy)
1773: {

1779:   PetscTryMethod(pc,"PCFieldSplitRestrictIS_C",(PC,IS),(pc,isy));
1780:   return(0);
1781: }


1784: static PetscErrorCode  PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy)
1785: {
1786:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1787:   PetscErrorCode    ierr;
1788:   PC_FieldSplitLink ilink = jac->head, next;
1789:   PetscInt          localsize,size,sizez,i;
1790:   const PetscInt    *ind, *indz;
1791:   PetscInt          *indc, *indcz;
1792:   PetscBool         flg;

1795:   ISGetLocalSize(isy,&localsize);
1796:   MPI_Scan(&localsize,&size,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)isy));
1797:   size -= localsize;
1798:   while(ilink) {
1799:     IS isrl,isr;
1800:     PC subpc;
1801:     ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl);
1802:     ISGetLocalSize(isrl,&localsize);
1803:     PetscMalloc1(localsize,&indc);
1804:     ISGetIndices(isrl,&ind);
1805:     PetscArraycpy(indc,ind,localsize);
1806:     ISRestoreIndices(isrl,&ind);
1807:     ISDestroy(&isrl);
1808:     for (i=0; i<localsize; i++) *(indc+i) += size;
1809:     ISCreateGeneral(PetscObjectComm((PetscObject)isy),localsize,indc,PETSC_OWN_POINTER,&isr);
1810:     PetscObjectReference((PetscObject)isr);
1811:     ISDestroy(&ilink->is);
1812:     ilink->is     = isr;
1813:     PetscObjectReference((PetscObject)isr);
1814:     ISDestroy(&ilink->is_col);
1815:     ilink->is_col = isr;
1816:     ISDestroy(&isr);
1817:     KSPGetPC(ilink->ksp, &subpc);
1818:     PetscObjectTypeCompare((PetscObject)subpc,PCFIELDSPLIT,&flg);
1819:     if(flg) {
1820:       IS iszl,isz;
1821:       MPI_Comm comm;
1822:       ISGetLocalSize(ilink->is,&localsize);
1823:       comm   = PetscObjectComm((PetscObject)ilink->is);
1824:       ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl);
1825:       MPI_Scan(&localsize,&sizez,1,MPIU_INT,MPI_SUM,comm);
1826:       sizez -= localsize;
1827:       ISGetLocalSize(iszl,&localsize);
1828:       PetscMalloc1(localsize,&indcz);
1829:       ISGetIndices(iszl,&indz);
1830:       PetscArraycpy(indcz,indz,localsize);
1831:       ISRestoreIndices(iszl,&indz);
1832:       ISDestroy(&iszl);
1833:       for (i=0; i<localsize; i++) *(indcz+i) += sizez;
1834:       ISCreateGeneral(comm,localsize,indcz,PETSC_OWN_POINTER,&isz);
1835:       PCFieldSplitRestrictIS(subpc,isz);
1836:       ISDestroy(&isz);
1837:     }
1838:     next = ilink->next;
1839:     ilink = next;
1840:   }
1841:   jac->isrestrict = PETSC_TRUE;
1842:   return(0);
1843: }

1845: static PetscErrorCode  PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is)
1846: {
1847:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1848:   PetscErrorCode    ierr;
1849:   PC_FieldSplitLink ilink, next = jac->head;
1850:   char              prefix[128];

1853:   if (jac->splitdefined) {
1854:     PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1855:     return(0);
1856:   }
1857:   PetscNew(&ilink);
1858:   if (splitname) {
1859:     PetscStrallocpy(splitname,&ilink->splitname);
1860:   } else {
1861:     PetscMalloc1(8,&ilink->splitname);
1862:     PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);
1863:   }
1864:   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 */
1865:   PetscObjectReference((PetscObject)is);
1866:   ISDestroy(&ilink->is);
1867:   ilink->is     = is;
1868:   PetscObjectReference((PetscObject)is);
1869:   ISDestroy(&ilink->is_col);
1870:   ilink->is_col = is;
1871:   ilink->next   = NULL;
1872:   KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1873:   KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1874:   PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1875:   KSPSetType(ilink->ksp,KSPPREONLY);
1876:   PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);

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

1881:   if (!next) {
1882:     jac->head       = ilink;
1883:     ilink->previous = NULL;
1884:   } else {
1885:     while (next->next) {
1886:       next = next->next;
1887:     }
1888:     next->next      = ilink;
1889:     ilink->previous = next;
1890:   }
1891:   jac->nsplits++;
1892:   return(0);
1893: }

1895: /*@C
1896:     PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner

1898:     Logically Collective on PC

1900:     Input Parameters:
1901: +   pc  - the preconditioner context
1902: .   splitname - name of this split, if NULL the number of the split is used
1903: .   n - the number of fields in this split
1904: -   fields - the fields in this split

1906:     Level: intermediate

1908:     Notes:
1909:     Use PCFieldSplitSetIS() to set a completely general set of indices as a field.

1911:      The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block
1912:      size is three then one can define a field as 0, or 1 or 2 or 0,1 or 0,2 or 1,2 which mean
1913:      0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x....
1914:      where the numbered entries indicate what is in the field.

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

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

1923: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS()

1925: @*/
1926: PetscErrorCode  PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1927: {

1933:   if (n < 1) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname);
1935:   PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col));
1936:   return(0);
1937: }

1939: /*@
1940:     PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat)

1942:     Logically Collective on PC

1944:     Input Parameters:
1945: +   pc  - the preconditioner object
1946: -   flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from

1948:     Options Database:
1949: .     -pc_fieldsplit_diag_use_amat

1951:     Level: intermediate

1953: .seealso: PCFieldSplitGetDiagUseAmat(), PCFieldSplitSetOffDiagUseAmat(), PCFIELDSPLIT

1955: @*/
1956: PetscErrorCode  PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg)
1957: {
1958:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1959:   PetscBool      isfs;

1964:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1965:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1966:   jac->diag_use_amat = flg;
1967:   return(0);
1968: }

1970: /*@
1971:     PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat)

1973:     Logically Collective on PC

1975:     Input Parameters:
1976: .   pc  - the preconditioner object

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


1982:     Level: intermediate

1984: .seealso: PCFieldSplitSetDiagUseAmat(), PCFieldSplitGetOffDiagUseAmat(), PCFIELDSPLIT

1986: @*/
1987: PetscErrorCode  PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg)
1988: {
1989:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1990:   PetscBool      isfs;

1996:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1997:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1998:   *flg = jac->diag_use_amat;
1999:   return(0);
2000: }

2002: /*@
2003:     PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat)

2005:     Logically Collective on PC

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

2011:     Options Database:
2012: .     -pc_fieldsplit_off_diag_use_amat

2014:     Level: intermediate

2016: .seealso: PCFieldSplitGetOffDiagUseAmat(), PCFieldSplitSetDiagUseAmat(), PCFIELDSPLIT

2018: @*/
2019: PetscErrorCode  PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg)
2020: {
2021:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
2022:   PetscBool      isfs;

2027:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2028:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
2029:   jac->offdiag_use_amat = flg;
2030:   return(0);
2031: }

2033: /*@
2034:     PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat)

2036:     Logically Collective on PC

2038:     Input Parameters:
2039: .   pc  - the preconditioner object

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


2045:     Level: intermediate

2047: .seealso: PCFieldSplitSetOffDiagUseAmat(), PCFieldSplitGetDiagUseAmat(), PCFIELDSPLIT

2049: @*/
2050: PetscErrorCode  PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg)
2051: {
2052:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
2053:   PetscBool      isfs;

2059:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2060:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
2061:   *flg = jac->offdiag_use_amat;
2062:   return(0);
2063: }



2067: /*@C
2068:     PCFieldSplitSetIS - Sets the exact elements for field

2070:     Logically Collective on PC

2072:     Input Parameters:
2073: +   pc  - the preconditioner context
2074: .   splitname - name of this split, if NULL the number of the split is used
2075: -   is - the index set that defines the vector elements in this field


2078:     Notes:
2079:     Use PCFieldSplitSetFields(), for fields defined by strided types.

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

2084:     Level: intermediate

2086: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize()

2088: @*/
2089: PetscErrorCode  PCFieldSplitSetIS(PC pc,const char splitname[],IS is)
2090: {

2097:   PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));
2098:   return(0);
2099: }

2101: /*@C
2102:     PCFieldSplitGetIS - Retrieves the elements for a field as an IS

2104:     Logically Collective on PC

2106:     Input Parameters:
2107: +   pc  - the preconditioner context
2108: -   splitname - name of this split

2110:     Output Parameter:
2111: -   is - the index set that defines the vector elements in this field, or NULL if the field is not found

2113:     Level: intermediate

2115: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS()

2117: @*/
2118: PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is)
2119: {

2126:   {
2127:     PC_FieldSplit     *jac  = (PC_FieldSplit*) pc->data;
2128:     PC_FieldSplitLink ilink = jac->head;
2129:     PetscBool         found;

2131:     *is = NULL;
2132:     while (ilink) {
2133:       PetscStrcmp(ilink->splitname, splitname, &found);
2134:       if (found) {
2135:         *is = ilink->is;
2136:         break;
2137:       }
2138:       ilink = ilink->next;
2139:     }
2140:   }
2141:   return(0);
2142: }

2144: /*@C
2145:     PCFieldSplitGetISByIndex - Retrieves the elements for a given index field as an IS

2147:     Logically Collective on PC

2149:     Input Parameters:
2150: +   pc  - the preconditioner context
2151: -   index - index of this split

2153:     Output Parameter:
2154: -   is - the index set that defines the vector elements in this field

2156:     Level: intermediate

2158: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitGetIS(), PCFieldSplitSetIS()

2160: @*/
2161: PetscErrorCode PCFieldSplitGetISByIndex(PC pc,PetscInt index,IS *is)
2162: {

2166:   if (index < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",index);
2169:   {
2170:     PC_FieldSplit     *jac  = (PC_FieldSplit*) pc->data;
2171:     PC_FieldSplitLink ilink = jac->head;
2172:     PetscInt          i     = 0;
2173:     if (index >= jac->nsplits) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",index,jac->nsplits);

2175:     while (i < index) {
2176:       ilink = ilink->next;
2177:       ++i;
2178:     }
2179:     PCFieldSplitGetIS(pc,ilink->splitname,is);
2180:   }
2181:   return(0);
2182: }

2184: /*@
2185:     PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the
2186:       fieldsplit preconditioner. If not set the matrix block size is used.

2188:     Logically Collective on PC

2190:     Input Parameters:
2191: +   pc  - the preconditioner context
2192: -   bs - the block size

2194:     Level: intermediate

2196: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields()

2198: @*/
2199: PetscErrorCode  PCFieldSplitSetBlockSize(PC pc,PetscInt bs)
2200: {

2206:   PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));
2207:   return(0);
2208: }

2210: /*@C
2211:    PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits

2213:    Collective on KSP

2215:    Input Parameter:
2216: .  pc - the preconditioner context

2218:    Output Parameters:
2219: +  n - the number of splits
2220: -  subksp - the array of KSP contexts

2222:    Note:
2223:    After PCFieldSplitGetSubKSP() the array of KSPs is to be freed by the user with PetscFree()
2224:    (not the KSP just the array that contains them).

2226:    You must call PCSetUp() before calling PCFieldSplitGetSubKSP().

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

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

2235:    Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs.
2236:       You can call PCFieldSplitGetSubKSP(pc,n,PETSC_NULL_KSP,ierr) to determine how large the
2237:       KSP array must be.


2240:    Level: advanced

2242: .seealso: PCFIELDSPLIT
2243: @*/
2244: PetscErrorCode  PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
2245: {

2251:   PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
2252:   return(0);
2253: }

2255: /*@C
2256:    PCFieldSplitSchurGetSubKSP - Gets the KSP contexts used inside the Schur complement based PCFIELDSPLIT

2258:    Collective on KSP

2260:    Input Parameter:
2261: .  pc - the preconditioner context

2263:    Output Parameters:
2264: +  n - the number of splits
2265: -  subksp - the array of KSP contexts

2267:    Note:
2268:    After PCFieldSplitSchurGetSubKSP() the array of KSPs is to be freed by the user with PetscFree()
2269:    (not the KSP just the array that contains them).

2271:    You must call PCSetUp() before calling PCFieldSplitSchurGetSubKSP().

2273:    If the fieldsplit type is of type PC_COMPOSITE_SCHUR, it returns (in order)
2274:    - the KSP used for the (1,1) block
2275:    - the KSP used for the Schur complement (not the one used for the interior Schur solver)
2276:    - the KSP used for the (1,1) block in the upper triangular factor (if different from that of the (1,1) block).

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

2280:    Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs.
2281:       You can call PCFieldSplitSchurGetSubKSP(pc,n,PETSC_NULL_KSP,ierr) to determine how large the
2282:       KSP array must be.

2284:    Level: advanced

2286: .seealso: PCFIELDSPLIT
2287: @*/
2288: PetscErrorCode  PCFieldSplitSchurGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
2289: {

2295:   PetscUseMethod(pc,"PCFieldSplitSchurGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
2296:   return(0);
2297: }

2299: /*@
2300:     PCFieldSplitSetSchurPre -  Indicates from what operator the preconditioner is constructucted for the Schur complement.
2301:       The default is the A11 matrix. 

2303:     Collective on PC

2305:     Input Parameters:
2306: +   pc      - the preconditioner context
2307: .   ptype   - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11 (default), PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_SCHUR_PRE_USER
2308:               PC_FIELDSPLIT_SCHUR_PRE_SELFP, and PC_FIELDSPLIT_SCHUR_PRE_FULL
2309: -   userpre - matrix to use for preconditioning, or NULL

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

2315:     Notes:
2316: $    If ptype is
2317: $        a11 - the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner
2318: $        matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix
2319: $        self - the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix:
2320: $             The only preconditioner that currently works with this symbolic respresentation matrix object is the PCLSC
2321: $             preconditioner
2322: $        user - the preconditioner for the Schur complement is generated from the user provided matrix (pre argument
2323: $             to this function).
2324: $        selfp - the preconditioning for the Schur complement is generated from an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01
2325: $             This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be
2326: $             lumped before extracting the diagonal using the additional option -fieldsplit_1_mat_schur_complement_ainv_type lump
2327: $        full - the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation computed internally by PCFIELDSPLIT (this is expensive)
2328: $             useful mostly as a test that the Schur complement approach can work for your problem

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

2334:     Level: intermediate

2336: .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType,
2337:           MatSchurComplementSetAinvType(), PCLSC

2339: @*/
2340: PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
2341: {

2346:   PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));
2347:   return(0);
2348: }

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

2352: /*@
2353:     PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be
2354:     preconditioned.  See PCFieldSplitSetSchurPre() for details.

2356:     Logically Collective on PC

2358:     Input Parameters:
2359: .   pc      - the preconditioner context

2361:     Output Parameters:
2362: +   ptype   - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11, PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER
2363: -   userpre - matrix to use for preconditioning (with PC_FIELDSPLIT_PRE_USER), or NULL

2365:     Level: intermediate

2367: .seealso: PCFieldSplitSetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC

2369: @*/
2370: PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
2371: {

2376:   PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));
2377:   return(0);
2378: }

2380: /*@
2381:     PCFieldSplitSchurGetS -  extract the MatSchurComplement object used by this PC in case it needs to be configured separately

2383:     Not collective

2385:     Input Parameter:
2386: .   pc      - the preconditioner context

2388:     Output Parameter:
2389: .   S       - the Schur complement matrix

2391:     Notes:
2392:     This matrix should not be destroyed using MatDestroy(); rather, use PCFieldSplitSchurRestoreS().

2394:     Level: advanced

2396: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurRestoreS()

2398: @*/
2399: PetscErrorCode  PCFieldSplitSchurGetS(PC pc,Mat *S)
2400: {
2402:   const char*    t;
2403:   PetscBool      isfs;
2404:   PC_FieldSplit  *jac;

2408:   PetscObjectGetType((PetscObject)pc,&t);
2409:   PetscStrcmp(t,PCFIELDSPLIT,&isfs);
2410:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
2411:   jac = (PC_FieldSplit*)pc->data;
2412:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
2413:   if (S) *S = jac->schur;
2414:   return(0);
2415: }

2417: /*@
2418:     PCFieldSplitSchurRestoreS -  restores the MatSchurComplement object used by this PC

2420:     Not collective

2422:     Input Parameters:
2423: +   pc      - the preconditioner context
2424: -   S       - the Schur complement matrix

2426:     Level: advanced

2428: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurGetS()

2430: @*/
2431: PetscErrorCode  PCFieldSplitSchurRestoreS(PC pc,Mat *S)
2432: {
2434:   const char*    t;
2435:   PetscBool      isfs;
2436:   PC_FieldSplit  *jac;

2440:   PetscObjectGetType((PetscObject)pc,&t);
2441:   PetscStrcmp(t,PCFIELDSPLIT,&isfs);
2442:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
2443:   jac = (PC_FieldSplit*)pc->data;
2444:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
2445:   if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten");
2446:   return(0);
2447: }


2450: static PetscErrorCode  PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
2451: {
2452:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

2456:   jac->schurpre = ptype;
2457:   if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
2458:     MatDestroy(&jac->schur_user);
2459:     jac->schur_user = pre;
2460:     PetscObjectReference((PetscObject)jac->schur_user);
2461:   }
2462:   return(0);
2463: }

2465: static PetscErrorCode  PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
2466: {
2467:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

2470:   *ptype = jac->schurpre;
2471:   *pre   = jac->schur_user;
2472:   return(0);
2473: }

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

2478:     Collective on PC

2480:     Input Parameters:
2481: +   pc  - the preconditioner context
2482: -   ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default

2484:     Options Database:
2485: .     -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full


2488:     Level: intermediate

2490:     Notes:
2491:     The FULL factorization is

2493: $   (A   B)  = (1       0) (A   0) (1  Ainv*B)  = L D U
2494: $   (C   E)    (C*Ainv  1) (0   S) (0     1  )

2496:     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,
2497:     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().

2499: $    If A and S are solved exactly
2500: $      *) FULL factorization is a direct solver.
2501: $      *) 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.
2502: $      *) 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.

2504:     If the iteration count is very low, consider using KSPFGMRES or KSPGCR which can use one less preconditioner
2505:     Section 1.5 Writing Application Codes with PETSc in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice.

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

2509:     Note that 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).

2511:     References:
2512: +   1. - Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000).
2513: -   2. - Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001).

2515: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCFieldSplitSetSchurScale()
2516: @*/
2517: PetscErrorCode  PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype)
2518: {

2523:   PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));
2524:   return(0);
2525: }

2527: static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype)
2528: {
2529:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2532:   jac->schurfactorization = ftype;
2533:   return(0);
2534: }

2536: /*@
2537:     PCFieldSplitSetSchurScale -  Controls the sign flip of S for PC_FIELDSPLIT_SCHUR_FACT_DIAG.

2539:     Collective on PC

2541:     Input Parameters:
2542: +   pc    - the preconditioner context
2543: -   scale - scaling factor for the Schur complement

2545:     Options Database:
2546: .     -pc_fieldsplit_schur_scale - default is -1.0

2548:     Level: intermediate

2550: .seealso: PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurFactType, PCFieldSplitSetSchurScale()
2551: @*/
2552: PetscErrorCode PCFieldSplitSetSchurScale(PC pc,PetscScalar scale)
2553: {

2559:   PetscTryMethod(pc,"PCFieldSplitSetSchurScale_C",(PC,PetscScalar),(pc,scale));
2560:   return(0);
2561: }

2563: static PetscErrorCode PCFieldSplitSetSchurScale_FieldSplit(PC pc,PetscScalar scale)
2564: {
2565:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2568:   jac->schurscale = scale;
2569:   return(0);
2570: }

2572: /*@C
2573:    PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement

2575:    Collective on KSP

2577:    Input Parameter:
2578: .  pc - the preconditioner context

2580:    Output Parameters:
2581: +  A00 - the (0,0) block
2582: .  A01 - the (0,1) block
2583: .  A10 - the (1,0) block
2584: -  A11 - the (1,1) block

2586:    Level: advanced

2588: .seealso: PCFIELDSPLIT
2589: @*/
2590: PetscErrorCode  PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11)
2591: {
2592:   PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;

2596:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach.");
2597:   if (A00) *A00 = jac->pmat[0];
2598:   if (A01) *A01 = jac->B;
2599:   if (A10) *A10 = jac->C;
2600:   if (A11) *A11 = jac->pmat[1];
2601:   return(0);
2602: }

2604: /*@
2605:     PCFieldSplitSetGKBTol -  Sets the solver tolerance for the generalized Golub-Kahan bidiagonalization preconditioner.

2607:     Collective on PC

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

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

2616:     Input Parameters:
2617: +   pc        - the preconditioner context
2618: -   tolerance - the solver tolerance

2620:     Options Database:
2621: .     -pc_fieldsplit_gkb_tol - default is 1e-5

2623:     Level: intermediate

2625: .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBDelay(), PCFieldSplitSetGKBNu(), PCFieldSplitSetGKBMaxit()
2626: @*/
2627: PetscErrorCode PCFieldSplitSetGKBTol(PC pc,PetscReal tolerance)
2628: {

2634:   PetscTryMethod(pc,"PCFieldSplitSetGKBTol_C",(PC,PetscReal),(pc,tolerance));
2635:   return(0);
2636: }

2638: static PetscErrorCode PCFieldSplitSetGKBTol_FieldSplit(PC pc,PetscReal tolerance)
2639: {
2640:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2643:   jac->gkbtol = tolerance;
2644:   return(0);
2645: }


2648: /*@
2649:     PCFieldSplitSetGKBMaxit -  Sets the maximum number of iterations for the generalized Golub-Kahan bidiagonalization
2650:     preconditioner.

2652:     Collective on PC

2654:     Input Parameters:
2655: +   pc     - the preconditioner context
2656: -   maxit  - the maximum number of iterations

2658:     Options Database:
2659: .     -pc_fieldsplit_gkb_maxit - default is 100

2661:     Level: intermediate

2663: .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBDelay(), PCFieldSplitSetGKBTol(), PCFieldSplitSetGKBNu()
2664: @*/
2665: PetscErrorCode PCFieldSplitSetGKBMaxit(PC pc,PetscInt maxit)
2666: {

2672:   PetscTryMethod(pc,"PCFieldSplitSetGKBMaxit_C",(PC,PetscInt),(pc,maxit));
2673:   return(0);
2674: }

2676: static PetscErrorCode PCFieldSplitSetGKBMaxit_FieldSplit(PC pc,PetscInt maxit)
2677: {
2678:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2681:   jac->gkbmaxit = maxit;
2682:   return(0);
2683: }

2685: /*@
2686:     PCFieldSplitSetGKBDelay -  Sets the delay in the lower bound error estimate in the generalized Golub-Kahan bidiagonalization
2687:     preconditioner.

2689:     Collective on PC

2691:     Notes:
2692:     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
2693:     is expressed as a truncated sum. The error at iteration k can only be measured at iteration (k + delay), and thus the algorithm needs
2694:     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

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

2698:     Input Parameters:
2699: +   pc     - the preconditioner context
2700: -   delay  - the delay window in the lower bound estimate

2702:     Options Database:
2703: .     -pc_fieldsplit_gkb_delay - default is 5

2705:     Level: intermediate

2707: .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBNu(), PCFieldSplitSetGKBTol(), PCFieldSplitSetGKBMaxit()
2708: @*/
2709: PetscErrorCode PCFieldSplitSetGKBDelay(PC pc,PetscInt delay)
2710: {

2716:   PetscTryMethod(pc,"PCFieldSplitSetGKBDelay_C",(PC,PetscInt),(pc,delay));
2717:   return(0);
2718: }

2720: static PetscErrorCode PCFieldSplitSetGKBDelay_FieldSplit(PC pc,PetscInt delay)
2721: {
2722:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2725:   jac->gkbdelay = delay;
2726:   return(0);
2727: }

2729: /*@
2730:     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.

2732:     Collective on PC

2734:     Notes:
2735:     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,
2736:     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
2737:     necessary to find a good balance in between the convergence of the inner and outer loop.

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

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

2743:     Input Parameters:
2744: +   pc     - the preconditioner context
2745: -   nu     - the shift parameter

2747:     Options Database:
2748: .     -pc_fieldsplit_gkb_nu - default is 1

2750:     Level: intermediate

2752: .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBDelay(), PCFieldSplitSetGKBTol(), PCFieldSplitSetGKBMaxit()
2753: @*/
2754: PetscErrorCode PCFieldSplitSetGKBNu(PC pc,PetscReal nu)
2755: {

2761:   PetscTryMethod(pc,"PCFieldSplitSetGKBNu_C",(PC,PetscReal),(pc,nu));
2762:   return(0);
2763: }

2765: static PetscErrorCode PCFieldSplitSetGKBNu_FieldSplit(PC pc,PetscReal nu)
2766: {
2767:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2770:   jac->gkbnu = nu;
2771:   return(0);
2772: }


2775: static PetscErrorCode  PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type)
2776: {
2777:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

2781:   jac->type = type;

2783:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",0);
2784:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);
2785:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);
2786:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);
2787:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",0);
2788:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBTol_C",0);
2789:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBMaxit_C",0);
2790:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBNu_C",0);
2791:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBDelay_C",0);

2793:   if (type == PC_COMPOSITE_SCHUR) {
2794:     pc->ops->apply = PCApply_FieldSplit_Schur;
2795:     pc->ops->view  = PCView_FieldSplit_Schur;

2797:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);
2798:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);
2799:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);
2800:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);
2801:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",PCFieldSplitSetSchurScale_FieldSplit);
2802:   } else if (type == PC_COMPOSITE_GKB){
2803:     pc->ops->apply = PCApply_FieldSplit_GKB;
2804:     pc->ops->view  = PCView_FieldSplit_GKB;

2806:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2807:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBTol_C",PCFieldSplitSetGKBTol_FieldSplit);
2808:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBMaxit_C",PCFieldSplitSetGKBMaxit_FieldSplit);
2809:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBNu_C",PCFieldSplitSetGKBNu_FieldSplit);
2810:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBDelay_C",PCFieldSplitSetGKBDelay_FieldSplit);
2811:   } else {
2812:     pc->ops->apply = PCApply_FieldSplit;
2813:     pc->ops->view  = PCView_FieldSplit;

2815:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2816:   }
2817:   return(0);
2818: }

2820: static PetscErrorCode  PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs)
2821: {
2822:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2825:   if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs);
2826:   if (jac->bs > 0 && jac->bs != bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot change fieldsplit blocksize from %D to %D after it has been set",jac->bs,bs);
2827:   jac->bs = bs;
2828:   return(0);
2829: }

2831: /*@
2832:    PCFieldSplitSetType - Sets the type of fieldsplit preconditioner.

2834:    Collective on PC

2836:    Input Parameter:
2837: +  pc - the preconditioner context
2838: -  type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR

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

2843:    Level: Intermediate

2845: .seealso: PCCompositeSetType()

2847: @*/
2848: PetscErrorCode  PCFieldSplitSetType(PC pc,PCCompositeType type)
2849: {

2854:   PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));
2855:   return(0);
2856: }

2858: /*@
2859:   PCFieldSplitGetType - Gets the type of fieldsplit preconditioner.

2861:   Not collective

2863:   Input Parameter:
2864: . pc - the preconditioner context

2866:   Output Parameter:
2867: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR

2869:   Level: Intermediate

2871: .seealso: PCCompositeSetType()
2872: @*/
2873: PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
2874: {
2875:   PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;

2880:   *type = jac->type;
2881:   return(0);
2882: }

2884: /*@
2885:    PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible.

2887:    Logically Collective

2889:    Input Parameters:
2890: +  pc   - the preconditioner context
2891: -  flg  - boolean indicating whether to use field splits defined by the DM

2893:    Options Database Key:
2894: .  -pc_fieldsplit_dm_splits

2896:    Level: Intermediate

2898: .seealso: PCFieldSplitGetDMSplits()

2900: @*/
2901: PetscErrorCode  PCFieldSplitSetDMSplits(PC pc,PetscBool flg)
2902: {
2903:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
2904:   PetscBool      isfs;

2910:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2911:   if (isfs) {
2912:     jac->dm_splits = flg;
2913:   }
2914:   return(0);
2915: }


2918: /*@
2919:    PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible.

2921:    Logically Collective

2923:    Input Parameter:
2924: .  pc   - the preconditioner context

2926:    Output Parameter:
2927: .  flg  - boolean indicating whether to use field splits defined by the DM

2929:    Level: Intermediate

2931: .seealso: PCFieldSplitSetDMSplits()

2933: @*/
2934: PetscErrorCode  PCFieldSplitGetDMSplits(PC pc,PetscBool* flg)
2935: {
2936:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
2937:   PetscBool      isfs;

2943:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2944:   if (isfs) {
2945:     if(flg) *flg = jac->dm_splits;
2946:   }
2947:   return(0);
2948: }

2950: /*@
2951:    PCFieldSplitGetDetectSaddlePoint - Returns flag indicating whether PCFieldSplit will attempt to automatically determine fields based on zero diagonal entries.

2953:    Logically Collective

2955:    Input Parameter:
2956: .  pc   - the preconditioner context

2958:    Output Parameter:
2959: .  flg  - boolean indicating whether to detect fields or not

2961:    Level: Intermediate

2963: .seealso: PCFIELDSPLIT, PCFieldSplitSetDetectSaddlePoint()

2965: @*/
2966: PetscErrorCode PCFieldSplitGetDetectSaddlePoint(PC pc,PetscBool *flg)
2967: {
2968:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2971:   *flg = jac->detect;
2972:   return(0);
2973: }

2975: /*@
2976:    PCFieldSplitSetDetectSaddlePoint - Sets flag indicating whether PCFieldSplit will attempt to automatically determine fields based on zero diagonal entries.

2978:    Logically Collective

2980:    Notes:
2981:    Also sets the split type to PC_COMPOSITE_SCHUR (see PCFieldSplitSetType()) and the Schur preconditioner type to PC_FIELDSPLIT_SCHUR_PRE_SELF (see PCFieldSplitSetSchurPre()).

2983:    Input Parameter:
2984: .  pc   - the preconditioner context

2986:    Output Parameter:
2987: .  flg  - boolean indicating whether to detect fields or not

2989:    Options Database Key:
2990: .  -pc_fieldsplit_detect_saddle_point

2992:    Level: Intermediate

2994: .seealso: PCFIELDSPLIT, PCFieldSplitSetDetectSaddlePoint(), PCFieldSplitSetType(), PCFieldSplitSetSchurPre()

2996: @*/
2997: PetscErrorCode PCFieldSplitSetDetectSaddlePoint(PC pc,PetscBool flg)
2998: {
2999:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

3003:   jac->detect = flg;
3004:   if (jac->detect) {
3005:     PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);
3006:     PCFieldSplitSetSchurPre(pc,PC_FIELDSPLIT_SCHUR_PRE_SELF,NULL);
3007:   }
3008:   return(0);
3009: }

3011: /* -------------------------------------------------------------------------------------*/
3012: /*MC
3013:    PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
3014:                   fields or groups of fields. See the users manual section "Solving Block Matrices" for more details.

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

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

3022:    Level: intermediate

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

3034:      Options prefixes for inner solvers when using the Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ .
3035:      For all other solvers they are -fieldsplit_%d_ for the dth field; use -fieldsplit_ for all fields.
3036:      The options prefix for the inner solver when using the Golub-Kahan biadiagonalization preconditioner is -fieldsplit_0_

3038:    Notes:
3039:     Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS()
3040:      to define a field by an arbitrary collection of entries.

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

3047: $     For the Schur complement preconditioner if J = ( A00 A01 )
3048: $                                                    ( A10 A11 )
3049: $     the preconditioner using full factorization is
3050: $              ( I   -ksp(A00) A01 ) ( inv(A00)     0  ) (     I          0  )
3051: $              ( 0         I       ) (   0      ksp(S) ) ( -A10 ksp(A00)  I  )
3052:      where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_.  S is the Schur complement
3053: $              S = A11 - A10 ksp(A00) A01
3054:      which is usually dense and not stored explicitly.  The action of ksp(S) is computed using the KSP solver with prefix -fieldsplit_splitname_ (where splitname was given
3055:      in providing the SECOND split or 1 if not give). For PCFieldSplitGetSubKSP() when field number is 0,
3056:      it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default
3057:      A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() for all the possible ways to construct the preconditioner for S.

3059:      The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above,
3060:      diag gives
3061: $              ( inv(A00)     0   )
3062: $              (   0      -ksp(S) )
3063:      note that slightly counter intuitively there is a negative in front of the ksp(S) so that the preconditioner is positive definite. For SPD matrices J, the sign flip
3064:      can be turned off with PCFieldSplitSetSchurScale() or by command line -pc_fieldsplit_schur_scale 1.0. The lower factorization is the inverse of
3065: $              (  A00   0 )
3066: $              (  A10   S )
3067:      where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of
3068: $              ( A00 A01 )
3069: $              (  0   S  )
3070:      where again the inverses of A00 and S are applied using KSPs.

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

3075:      The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1.
3076:      Generally it should be used with the AIJ format.

3078:      The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see,
3079:      for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT
3080:      inside a smoother resulting in "Distributive Smoothers".

3082:    There is a nice discussion of block preconditioners in

3084: [El08] A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations
3085:        Howard Elman, V.E. Howle, John Shadid, Robert Shuttleworth, Ray Tuminaro, Journal of Computational Physics 227 (2008) 1790--1808
3086:        http://chess.cs.umd.edu/~elman/papers/tax.pdf

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

3091:    The generalized Golub-Kahan bidiagonalization preconditioner (gkb) can be applied to symmetric 2x2 block matrices of the shape
3092: $        ( A00  A01 )
3093: $        ( A01' 0   )
3094:    with A00 positive semi-definite. The implementation follows [Ar13]. Therein, we choose N := 1/nu * I and the (1,1)-block of the matrix is modified to H = A00 + nu*A01*A01'.
3095:    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_.

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

3099: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, PCLSC,
3100:            PCFieldSplitGetSubKSP(), PCFieldSplitSchurGetSubKSP(), PCFieldSplitSetFields(),
3101:            PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(), PCFieldSplitSetSchurFactType(),
3102:            MatSchurComplementSetAinvType(), PCFieldSplitSetSchurScale(), PCFieldSplitSetDetectSaddlePoint()
3103: M*/

3105: PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
3106: {
3108:   PC_FieldSplit  *jac;

3111:   PetscNewLog(pc,&jac);

3113:   jac->bs                 = -1;
3114:   jac->nsplits            = 0;
3115:   jac->type               = PC_COMPOSITE_MULTIPLICATIVE;
3116:   jac->schurpre           = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
3117:   jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
3118:   jac->schurscale         = -1.0;
3119:   jac->dm_splits          = PETSC_TRUE;
3120:   jac->detect             = PETSC_FALSE;
3121:   jac->gkbtol             = 1e-5;
3122:   jac->gkbdelay           = 5;
3123:   jac->gkbnu              = 1;
3124:   jac->gkbmaxit           = 100;
3125:   jac->gkbmonitor         = PETSC_FALSE;

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

3129:   pc->ops->apply           = PCApply_FieldSplit;
3130:   pc->ops->applytranspose  = PCApplyTranspose_FieldSplit;
3131:   pc->ops->setup           = PCSetUp_FieldSplit;
3132:   pc->ops->reset           = PCReset_FieldSplit;
3133:   pc->ops->destroy         = PCDestroy_FieldSplit;
3134:   pc->ops->setfromoptions  = PCSetFromOptions_FieldSplit;
3135:   pc->ops->view            = PCView_FieldSplit;
3136:   pc->ops->applyrichardson = NULL;

3138:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",PCFieldSplitSchurGetSubKSP_FieldSplit);
3139:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
3140:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);
3141:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);
3142:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);
3143:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);
3144:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",PCFieldSplitRestrictIS_FieldSplit);
3145:   return(0);
3146: }