Actual source code: fieldsplit.c

petsc-master 2019-08-18
Report Typos and Errors


  3:  #include <petsc/private/pcimpl.h>
  4: #include <petsc/private/kspimpl.h>    /*  This is needed to provide the appropriate PETSC_EXTERN for KSP_Solve_FS ....*/
  5:  #include <petscdm.h>

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

 10: 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;

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

333:     PetscViewerASCIIPopTab(viewer);
334:   }

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

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


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

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

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

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

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

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

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

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

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

505:         MatGetOwnershipRange(pc->mat,&nmin,&nmax);
506:         MatFindOffBlockDiagonalEntries(pc->mat,&coupling);
507:         ISCreateStride(PetscObjectComm((PetscObject)pc->mat),nmax-nmin,nmin,1,&rest);
508:         ISSetIdentity(rest);
509:         PCFieldSplitSetIS(pc,"0",rest);
510:         PCFieldSplitSetIS(pc,"1",coupling);
511:         ISDestroy(&coupling);
512:         ISDestroy(&rest);
513:       } else {
514:         PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&fieldsplit_default,NULL);
515:         if (!fieldsplit_default) {
516:           /* Allow user to set fields from command line,  if bs was known at the time of PCSetFromOptions_FieldSplit()
517:            then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */
518:           PCFieldSplitSetRuntimeSplits_Private(pc);
519:           if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
520:         }
521:         if ((fieldsplit_default || !jac->splitdefined) && !jac->isrestrict) {
522:           PetscInfo(pc,"Using default splitting of fields\n");
523:           for (i=0; i<jac->bs; i++) {
524:             char splitname[8];
525:             PetscSNPrintf(splitname,sizeof(splitname),"%D",i);
526:             PCFieldSplitSetFields(pc,splitname,1,&i,&i);
527:           }
528:           jac->defaultsplit = PETSC_TRUE;
529:         }
530:       }
531:     }
532:   } else if (jac->nsplits == 1) {
533:     if (ilink->is) {
534:       IS       is2;
535:       PetscInt nmin,nmax;

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

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

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

555:   MatHermitianTranspose(C,MAT_INITIAL_MATRIX,&T);            /* Test if augmented matrix is symmetric */
556:   MatAXPY(T,-1.0,B,DIFFERENT_NONZERO_PATTERN);
557:   MatNorm(T,NORM_1,&nrmT);
558:   MatNorm(B,NORM_1,&nrmB);
559:   if (nrmB > 0) {
560:     if (nrmT/nrmB >= PETSC_SMALL) {
561:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix is not symmetric/hermitian, GKB is not applicable.");
562:     }
563:   }
564:   /* Compute augmented Lagrangian matrix H = A00 + nu*A01*A01'. This corresponds to */
565:   /* setting N := 1/nu*I in [Ar13].                                                 */
566:   MatHermitianTranspose(B,MAT_INITIAL_MATRIX,&BT);
567:   MatMatMult(B,BT,MAT_INITIAL_MATRIX,PETSC_DEFAULT,H);       /* H = A01*A01'          */
568:   MatAYPX(*H,gkbnu,A,DIFFERENT_NONZERO_PATTERN);             /* H = A00 + nu*A01*A01' */

570:   MatDestroy(&BT);
571:   MatDestroy(&T);
572:   return(0);
573: }

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

577: static PetscErrorCode PCSetUp_FieldSplit(PC pc)
578: {
579:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
580:   PetscErrorCode    ierr;
581:   PC_FieldSplitLink ilink;
582:   PetscInt          i,nsplit;
583:   PetscBool         sorted, sorted_col;

586:   pc->failedreason = PC_NOERROR;
587:   PCFieldSplitSetDefaults(pc);
588:   nsplit = jac->nsplits;
589:   ilink  = jac->head;

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

595:     jac->issetup = PETSC_TRUE;

597:     /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */
598:     if (jac->defaultsplit || !ilink->is) {
599:       if (jac->bs <= 0) jac->bs = nsplit;
600:     }
601:     bs     = jac->bs;
602:     MatGetOwnershipRange(pc->pmat,&rstart,&rend);
603:     nslots = (rend - rstart)/bs;
604:     for (i=0; i<nsplit; i++) {
605:       if (jac->defaultsplit) {
606:         ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+i,nsplit,&ilink->is);
607:         ISDuplicate(ilink->is,&ilink->is_col);
608:       } else if (!ilink->is) {
609:         if (ilink->nfields > 1) {
610:           PetscInt *ii,*jj,j,k,nfields = ilink->nfields,*fields = ilink->fields,*fields_col = ilink->fields_col;
611:           PetscMalloc1(ilink->nfields*nslots,&ii);
612:           PetscMalloc1(ilink->nfields*nslots,&jj);
613:           for (j=0; j<nslots; j++) {
614:             for (k=0; k<nfields; k++) {
615:               ii[nfields*j + k] = rstart + bs*j + fields[k];
616:               jj[nfields*j + k] = rstart + bs*j + fields_col[k];
617:             }
618:           }
619:           ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);
620:           ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,jj,PETSC_OWN_POINTER,&ilink->is_col);
621:           ISSetBlockSize(ilink->is, nfields);
622:           ISSetBlockSize(ilink->is_col, nfields);
623:         } else {
624:           ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields[0],bs,&ilink->is);
625:           ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields_col[0],bs,&ilink->is_col);
626:         }
627:       }
628:       ISSorted(ilink->is,&sorted);
629:       if (ilink->is_col) { ISSorted(ilink->is_col,&sorted_col); }
630:       if (!sorted || !sorted_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split");
631:       ilink = ilink->next;
632:     }
633:   }

635:   ilink = jac->head;
636:   if (!jac->pmat) {
637:     Vec xtmp;

639:     MatCreateVecs(pc->pmat,&xtmp,NULL);
640:     PetscMalloc1(nsplit,&jac->pmat);
641:     PetscMalloc2(nsplit,&jac->x,nsplit,&jac->y);
642:     for (i=0; i<nsplit; i++) {
643:       MatNullSpace sp;

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

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

682:     for (i=0; i<nsplit; i++) {
683:       Mat pmat;

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

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

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

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

731:   if (jac->type != PC_COMPOSITE_ADDITIVE  && jac->type != PC_COMPOSITE_SCHUR && jac->type != PC_COMPOSITE_GKB) {
732:     /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */
733:     /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */
734:     ilink = jac->head;
735:     if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) {
736:       /* 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 */
737:       if (!jac->Afield) {
738:         PetscCalloc1(nsplit,&jac->Afield);
739:         if (jac->offdiag_use_amat) {
740:           MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);
741:         } else {
742:           MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);
743:         }
744:       } else {
745:         MatReuse scall;
746:         if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
747:           for (i=0; i<nsplit; i++) {
748:             MatDestroy(&jac->Afield[1]);
749:           }
750:           scall = MAT_INITIAL_MATRIX;
751:         } else scall = MAT_REUSE_MATRIX;

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

779:         for (i=0; i<nsplit; i++) {
780:           if (jac->offdiag_use_amat) {
781:             MatCreateSubMatrix(pc->mat,ilink->is,NULL,scall,&jac->Afield[i]);
782:           } else {
783:             MatCreateSubMatrix(pc->pmat,ilink->is,NULL,scall,&jac->Afield[i]);
784:           }
785:           ilink = ilink->next;
786:         }
787:       }
788:     }
789:   }

791:   if (jac->type == PC_COMPOSITE_SCHUR) {
792:     IS          ccis;
793:     PetscBool   isspd;
794:     PetscInt    rstart,rend;
795:     char        lscname[256];
796:     PetscObject LSC_L;

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

800:     /* If pc->mat is SPD, don't scale by -1 the Schur complement */
801:     if (jac->schurscale == (PetscScalar)-1.0) {
802:       MatGetOption(pc->pmat,MAT_SPD,&isspd);
803:       jac->schurscale = (isspd == PETSC_TRUE) ? 1.0 : -1.0;
804:     }

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

809:     /* need to handle case when one is resetting up the preconditioner */
810:     if (jac->schur) {
811:       KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper;

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

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

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

877:       /* Note: this is not true in general */
878:       MatGetNullSpace(jac->mat[1], &sp);
879:       if (sp) {
880:         MatSetNullSpace(jac->schur, sp);
881:       }

883:       PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);
884:       PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
885:       if (flg) {
886:         DM  dmInner;
887:         KSP kspInner;
888:         PC  pcInner;

890:         MatSchurComplementGetKSP(jac->schur, &kspInner);
891:         KSPReset(kspInner);
892:         KSPSetOperators(kspInner,jac->mat[0],jac->pmat[0]);
893:         PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
894:         /* Indent this deeper to emphasize the "inner" nature of this solver. */
895:         PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject) pc, 2);
896:         PetscObjectIncrementTabLevel((PetscObject)kspInner->pc, (PetscObject) pc, 2);
897:         KSPSetOptionsPrefix(kspInner, schurprefix);

899:         /* Set DM for new solver */
900:         KSPGetDM(jac->head->ksp, &dmInner);
901:         KSPSetDM(kspInner, dmInner);
902:         KSPSetDMActive(kspInner, PETSC_FALSE);

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

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

927:         MatSchurComplementGetKSP(jac->schur, &kspInner);
928:         KSPGetPC(kspInner, &pcInner);
929:         PetscObjectTypeCompare((PetscObject)pcInner, PCKSP, &flg);
930:         if (flg) {
931:           KSP ksp;

933:           PCKSPGetKSP(pcInner, &ksp);
934:           if (ksp == jac->head->ksp) {
935:             PCSetUseAmat(pcInner, PETSC_TRUE);
936:           }
937:         }
938:       }
939:       PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);
940:       PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
941:       if (flg) {
942:         DM dmInner;

944:         PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
945:         KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);
946:         KSPSetErrorIfNotConverged(jac->kspupper,pc->erroriffailure);
947:         KSPSetOptionsPrefix(jac->kspupper, schurprefix);
948:         PetscObjectIncrementTabLevel((PetscObject)jac->kspupper, (PetscObject) pc, 1);
949:         PetscObjectIncrementTabLevel((PetscObject)jac->kspupper->pc, (PetscObject) pc, 1);
950:         KSPGetDM(jac->head->ksp, &dmInner);
951:         KSPSetDM(jac->kspupper, dmInner);
952:         KSPSetDMActive(jac->kspupper, PETSC_FALSE);
953:         KSPSetFromOptions(jac->kspupper);
954:         KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0]);
955:         VecDuplicate(jac->head->x, &jac->head->z);
956:       } else {
957:         jac->kspupper = jac->head->ksp;
958:         PetscObjectReference((PetscObject) jac->head->ksp);
959:       }

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

995:     /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */
996:     PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);
997:     PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);
998:     if (!LSC_L) {PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);}
999:     if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);}
1000:     PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);
1001:     PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);
1002:     if (!LSC_L) {PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);}
1003:     if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);}
1004:   } else if (jac->type == PC_COMPOSITE_GKB) {
1005:     IS          ccis;
1006:     PetscInt    rstart,rend;

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

1010:     ilink = jac->head;

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

1015:     ISComplement(ilink->is_col,rstart,rend,&ccis);
1016:     if (jac->offdiag_use_amat) {
1017:      MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
1018:     } else {
1019:       MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
1020:     }
1021:     ISDestroy(&ccis);
1022:     /* Create work vectors for GKB algorithm */
1023:     VecDuplicate(ilink->x,&jac->u);
1024:     VecDuplicate(ilink->x,&jac->Hu);
1025:     VecDuplicate(ilink->x,&jac->w2);
1026:     ilink = ilink->next;
1027:     ISComplement(ilink->is_col,rstart,rend,&ccis);
1028:     if (jac->offdiag_use_amat) {
1029:       MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
1030:     } else {
1031:       MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
1032:     }
1033:     ISDestroy(&ccis);
1034:     /* Create work vectors for GKB algorithm */
1035:     VecDuplicate(ilink->x,&jac->v);
1036:     VecDuplicate(ilink->x,&jac->d);
1037:     VecDuplicate(ilink->x,&jac->w1);
1038:     MatGolubKahanComputeExplicitOperator(jac->mat[0],jac->B,jac->C,&jac->H,jac->gkbnu);
1039:     PetscCalloc1(jac->gkbdelay,&jac->vecz);

1041:     ilink = jac->head;
1042:     KSPSetOperators(ilink->ksp,jac->H,jac->H);
1043:     if (!jac->suboptionsset) {KSPSetFromOptions(ilink->ksp);}
1044:     /* Create gkb_monitor context */
1045:     if (jac->gkbmonitor) {
1046:       PetscInt  tablevel;
1047:       PetscViewerCreate(PETSC_COMM_WORLD,&jac->gkbviewer);
1048:       PetscViewerSetType(jac->gkbviewer,PETSCVIEWERASCII);
1049:       PetscObjectGetTabLevel((PetscObject)ilink->ksp,&tablevel);
1050:       PetscViewerASCIISetTab(jac->gkbviewer,tablevel);
1051:       PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)ilink->ksp,1);
1052:     }
1053:   } else {
1054:     /* set up the individual splits' PCs */
1055:     i     = 0;
1056:     ilink = jac->head;
1057:     while (ilink) {
1058:       KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i]);
1059:       /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
1060:       if (!jac->suboptionsset) {KSPSetFromOptions(ilink->ksp);}
1061:       i++;
1062:       ilink = ilink->next;
1063:     }
1064:   }

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

1070: #define FieldSplitSplitSolveAdd(ilink,xx,yy) \
1071:   (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
1072:    VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
1073:    PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\
1074:    KSPSolve(ilink->ksp,ilink->x,ilink->y) ||                               \
1075:    KSPCheckSolve(ilink->ksp,pc,ilink->y)  || \
1076:    PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\
1077:    VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) ||  \
1078:    VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE))

1080: static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y)
1081: {
1082:   PC_FieldSplit      *jac = (PC_FieldSplit*)pc->data;
1083:   PetscErrorCode     ierr;
1084:   PC_FieldSplitLink  ilinkA = jac->head, ilinkD = ilinkA->next;
1085:   KSP                kspA   = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;

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

1163:     PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1164:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
1165:     KSPCheckSolve(jac->kspschur,pc,ilinkD->y);
1166:     PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
1167:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);

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

1194: static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y)
1195: {
1196:   PC_FieldSplit      *jac = (PC_FieldSplit*)pc->data;
1197:   PetscErrorCode     ierr;
1198:   PC_FieldSplitLink  ilink = jac->head;
1199:   PetscInt           cnt,bs;

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

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

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


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

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

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

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

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

1339:   /* First step of algorithm */
1340:   VecNorm(work1,NORM_2,&beta);                   /* beta = sqrt(nu*c'*c)*/
1341:   KSPCheckDot(ksp,beta);
1342:   beta  = PetscSqrtScalar(nu)*beta;
1343:   VecAXPBY(v,nu/beta,0.0,work1);                   /* v = nu/beta *c      */
1344:   MatMult(jac->B,v,work2);                       /* u = H^{-1}*B*v      */
1345:   PetscLogEventBegin(ilinkA->event,ksp,work2,u,NULL);
1346:   KSPSolve(ksp,work2,u);
1347:   KSPCheckSolve(ksp,pc,u);
1348:   PetscLogEventEnd(ilinkA->event,ksp,work2,u,NULL);
1349:   MatMult(jac->H,u,Hu);                          /* alpha = u'*H*u      */
1350:   VecDot(Hu,u,&alpha);
1351:   KSPCheckDot(ksp,alpha);
1352:   if (PetscRealPart(alpha) <= 0.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_NOT_CONVERGED,"GKB preconditioner diverged, H is not positive definite");
1353:   alpha = PetscSqrtScalar(PetscAbsScalar(alpha));
1354:   VecScale(u,1.0/alpha);
1355:   VecAXPBY(d,1.0/alpha,0.0,v);                       /* v = nu/beta *c      */

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

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

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

1370:   while (iterGKB < jac->gkbmaxit && lowbnd > jac->gkbtol) {
1371:     iterGKB += 1;
1372:     MatMultHermitianTranspose(jac->B,u,work1); /* v <- nu*(B'*u-alpha/nu*v) */
1373:     VecAXPBY(v,nu,-alpha,work1);
1374:     VecNorm(v,NORM_2,&beta);                   /* beta = sqrt(nu)*v'*v      */
1375:     beta  = beta/PetscSqrtScalar(nu);
1376:     VecScale(v,1.0/beta);
1377:     MatMult(jac->B,v,work2);                  /* u <- H^{-1}*(B*v-beta*H*u) */
1378:     MatMult(jac->H,u,Hu);
1379:     VecAXPY(work2,-beta,Hu);
1380:     PetscLogEventBegin(ilinkA->event,ksp,work2,u,NULL);
1381:     KSPSolve(ksp,work2,u);
1382:     KSPCheckSolve(ksp,pc,u);
1383:     PetscLogEventEnd(ilinkA->event,ksp,work2,u,NULL);
1384:     MatMult(jac->H,u,Hu);                      /* alpha = u'*H*u            */
1385:     VecDot(Hu,u,&alpha);
1386:     KSPCheckDot(ksp,alpha);
1387:     if (PetscRealPart(alpha) <= 0.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_NOT_CONVERGED,"GKB preconditioner diverged, H is not positive definite");
1388:     alpha = PetscSqrtScalar(PetscAbsScalar(alpha));
1389:     VecScale(u,1.0/alpha);

1391:     z = -beta/alpha*z;                                            /* z <- beta/alpha*z     */
1392:     vecz[0] = z;

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

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

1410:     for (j=0; j<jac->gkbdelay-1; j++) {
1411:       vecz[jac->gkbdelay-j-1] = vecz[jac->gkbdelay-j-2];
1412:     }
1413:     if (jac->gkbmonitor) {
1414:       PetscViewerASCIIPrintf(jac->gkbviewer,"%3D GKB Lower bound estimate %14.12e\n",iterGKB,lowbnd);
1415:     }
1416:   }

1418:   VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1419:   VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1420:   VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
1421:   VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);

1423:   return(0);
1424: }


1427: #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \
1428:   (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1429:    VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1430:    PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \
1431:    KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) ||                  \
1432:    KSPCheckSolve(ilink->ksp,pc,ilink->x) || \
1433:    PetscLogEventEnd(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) ||   \
1434:    VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \
1435:    VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE))

1437: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y)
1438: {
1439:   PC_FieldSplit      *jac = (PC_FieldSplit*)pc->data;
1440:   PetscErrorCode     ierr;
1441:   PC_FieldSplitLink  ilink = jac->head;
1442:   PetscInt           bs;

1445:   if (jac->type == PC_COMPOSITE_ADDITIVE) {
1446:     if (jac->defaultsplit) {
1447:       VecGetBlockSize(x,&bs);
1448:       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);
1449:       VecGetBlockSize(y,&bs);
1450:       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);
1451:       VecStrideGatherAll(x,jac->x,INSERT_VALUES);
1452:       while (ilink) {
1453:         PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1454:         KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);
1455:         KSPCheckSolve(ilink->ksp,pc,ilink->y);
1456:         PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1457:         ilink = ilink->next;
1458:       }
1459:       VecStrideScatterAll(jac->y,y,INSERT_VALUES);
1460:     } else {
1461:       VecSet(y,0.0);
1462:       while (ilink) {
1463:         FieldSplitSplitSolveAddTranspose(ilink,x,y);
1464:         ilink = ilink->next;
1465:       }
1466:     }
1467:   } else {
1468:     if (!jac->w1) {
1469:       VecDuplicate(x,&jac->w1);
1470:       VecDuplicate(x,&jac->w2);
1471:     }
1472:     VecSet(y,0.0);
1473:     if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1474:       FieldSplitSplitSolveAddTranspose(ilink,x,y);
1475:       while (ilink->next) {
1476:         ilink = ilink->next;
1477:         MatMultTranspose(pc->mat,y,jac->w1);
1478:         VecWAXPY(jac->w2,-1.0,jac->w1,x);
1479:         FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1480:       }
1481:       while (ilink->previous) {
1482:         ilink = ilink->previous;
1483:         MatMultTranspose(pc->mat,y,jac->w1);
1484:         VecWAXPY(jac->w2,-1.0,jac->w1,x);
1485:         FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1486:       }
1487:     } else {
1488:       while (ilink->next) {   /* get to last entry in linked list */
1489:         ilink = ilink->next;
1490:       }
1491:       FieldSplitSplitSolveAddTranspose(ilink,x,y);
1492:       while (ilink->previous) {
1493:         ilink = ilink->previous;
1494:         MatMultTranspose(pc->mat,y,jac->w1);
1495:         VecWAXPY(jac->w2,-1.0,jac->w1,x);
1496:         FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1497:       }
1498:     }
1499:   }
1500:   return(0);
1501: }

1503: static PetscErrorCode PCReset_FieldSplit(PC pc)
1504: {
1505:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1506:   PetscErrorCode    ierr;
1507:   PC_FieldSplitLink ilink = jac->head,next;

1510:   while (ilink) {
1511:     KSPDestroy(&ilink->ksp);
1512:     VecDestroy(&ilink->x);
1513:     VecDestroy(&ilink->y);
1514:     VecDestroy(&ilink->z);
1515:     VecScatterDestroy(&ilink->sctx);
1516:     ISDestroy(&ilink->is);
1517:     ISDestroy(&ilink->is_col);
1518:     PetscFree(ilink->splitname);
1519:     PetscFree(ilink->fields);
1520:     PetscFree(ilink->fields_col);
1521:     next  = ilink->next;
1522:     PetscFree(ilink);
1523:     ilink = next;
1524:   }
1525:   jac->head = NULL;
1526:   PetscFree2(jac->x,jac->y);
1527:   if (jac->mat && jac->mat != jac->pmat) {
1528:     MatDestroyMatrices(jac->nsplits,&jac->mat);
1529:   } else if (jac->mat) {
1530:     jac->mat = NULL;
1531:   }
1532:   if (jac->pmat) {MatDestroyMatrices(jac->nsplits,&jac->pmat);}
1533:   if (jac->Afield) {MatDestroyMatrices(jac->nsplits,&jac->Afield);}
1534:   jac->nsplits = 0;
1535:   VecDestroy(&jac->w1);
1536:   VecDestroy(&jac->w2);
1537:   MatDestroy(&jac->schur);
1538:   MatDestroy(&jac->schurp);
1539:   MatDestroy(&jac->schur_user);
1540:   KSPDestroy(&jac->kspschur);
1541:   KSPDestroy(&jac->kspupper);
1542:   MatDestroy(&jac->B);
1543:   MatDestroy(&jac->C);
1544:   MatDestroy(&jac->H);
1545:   VecDestroy(&jac->u);
1546:   VecDestroy(&jac->v);
1547:   VecDestroy(&jac->Hu);
1548:   VecDestroy(&jac->d);
1549:   PetscFree(jac->vecz);
1550:   PetscViewerDestroy(&jac->gkbviewer);
1551:   jac->isrestrict = PETSC_FALSE;
1552:   return(0);
1553: }

1555: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
1556: {
1557:   PetscErrorCode    ierr;

1560:   PCReset_FieldSplit(pc);
1561:   PetscFree(pc->data);
1562:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",NULL);
1563:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);
1564:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);
1565:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);
1566:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);
1567:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);
1568:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);
1569:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);
1570:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);
1571:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",NULL);
1572:   return(0);
1573: }

1575: static PetscErrorCode PCSetFromOptions_FieldSplit(PetscOptionItems *PetscOptionsObject,PC pc)
1576: {
1577:   PetscErrorCode  ierr;
1578:   PetscInt        bs;
1579:   PetscBool       flg;
1580:   PC_FieldSplit   *jac = (PC_FieldSplit*)pc->data;
1581:   PCCompositeType ctype;

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

1625: /*------------------------------------------------------------------------------------*/

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

1636:   if (jac->splitdefined) {
1637:     PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1638:     return(0);
1639:   }
1640:   for (i=0; i<n; i++) {
1641:     if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs);
1642:     if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]);
1643:   }
1644:   PetscNew(&ilink);
1645:   if (splitname) {
1646:     PetscStrallocpy(splitname,&ilink->splitname);
1647:   } else {
1648:     PetscMalloc1(3,&ilink->splitname);
1649:     PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);
1650:   }
1651:   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 */
1652:   PetscMalloc1(n,&ilink->fields);
1653:   PetscArraycpy(ilink->fields,fields,n);
1654:   PetscMalloc1(n,&ilink->fields_col);
1655:   PetscArraycpy(ilink->fields_col,fields_col,n);

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

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

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

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

1688:   *subksp = NULL;
1689:   if (n) *n = 0;
1690:   if (jac->type == PC_COMPOSITE_SCHUR) {
1691:     PetscInt nn;

1693:     if (!jac->schur) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitSchurGetSubKSP()");
1694:     if (jac->nsplits != 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unexpected number of splits %D != 2",jac->nsplits);
1695:     nn   = jac->nsplits + (jac->kspupper != jac->head->ksp ? 1 : 0);
1696:     PetscMalloc1(nn,subksp);
1697:     (*subksp)[0] = jac->head->ksp;
1698:     (*subksp)[1] = jac->kspschur;
1699:     if (jac->kspupper != jac->head->ksp) (*subksp)[2] = jac->kspupper;
1700:     if (n) *n = nn;
1701:   }
1702:   return(0);
1703: }

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

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

1715:   (*subksp)[1] = jac->kspschur;
1716:   if (n) *n = jac->nsplits;
1717:   return(0);
1718: }

1720: static PetscErrorCode  PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1721: {
1722:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1723:   PetscErrorCode    ierr;
1724:   PetscInt          cnt   = 0;
1725:   PC_FieldSplitLink ilink = jac->head;

1728:   PetscMalloc1(jac->nsplits,subksp);
1729:   while (ilink) {
1730:     (*subksp)[cnt++] = ilink->ksp;
1731:     ilink            = ilink->next;
1732:   }
1733:   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);
1734:   if (n) *n = jac->nsplits;
1735:   return(0);
1736: }

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

1741:     Input Parameters:
1742: +   pc  - the preconditioner context
1743: -   is - the index set that defines the indices to which the fieldsplit is to be restricted

1745:     Level: advanced

1747: @*/
1748: PetscErrorCode  PCFieldSplitRestrictIS(PC pc,IS isy)
1749: {

1755:   PetscTryMethod(pc,"PCFieldSplitRestrictIS_C",(PC,IS),(pc,isy));
1756:   return(0);
1757: }


1760: static PetscErrorCode  PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy)
1761: {
1762:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1763:   PetscErrorCode    ierr;
1764:   PC_FieldSplitLink ilink = jac->head, next;
1765:   PetscInt          localsize,size,sizez,i;
1766:   const PetscInt    *ind, *indz;
1767:   PetscInt          *indc, *indcz;
1768:   PetscBool         flg;

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

1821: static PetscErrorCode  PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is)
1822: {
1823:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1824:   PetscErrorCode    ierr;
1825:   PC_FieldSplitLink ilink, next = jac->head;
1826:   char              prefix[128];

1829:   if (jac->splitdefined) {
1830:     PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1831:     return(0);
1832:   }
1833:   PetscNew(&ilink);
1834:   if (splitname) {
1835:     PetscStrallocpy(splitname,&ilink->splitname);
1836:   } else {
1837:     PetscMalloc1(8,&ilink->splitname);
1838:     PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);
1839:   }
1840:   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 */
1841:   PetscObjectReference((PetscObject)is);
1842:   ISDestroy(&ilink->is);
1843:   ilink->is     = is;
1844:   PetscObjectReference((PetscObject)is);
1845:   ISDestroy(&ilink->is_col);
1846:   ilink->is_col = is;
1847:   ilink->next   = NULL;
1848:   KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1849:   KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1850:   PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1851:   KSPSetType(ilink->ksp,KSPPREONLY);
1852:   PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);

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

1857:   if (!next) {
1858:     jac->head       = ilink;
1859:     ilink->previous = NULL;
1860:   } else {
1861:     while (next->next) {
1862:       next = next->next;
1863:     }
1864:     next->next      = ilink;
1865:     ilink->previous = next;
1866:   }
1867:   jac->nsplits++;
1868:   return(0);
1869: }

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

1874:     Logically Collective on PC

1876:     Input Parameters:
1877: +   pc  - the preconditioner context
1878: .   splitname - name of this split, if NULL the number of the split is used
1879: .   n - the number of fields in this split
1880: -   fields - the fields in this split

1882:     Level: intermediate

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

1887:      The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block
1888:      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
1889:      0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x....
1890:      where the numbered entries indicate what is in the field.

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

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

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

1901: @*/
1902: PetscErrorCode  PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1903: {

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

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

1918:     Logically Collective on PC

1920:     Input Parameters:
1921: +   pc  - the preconditioner object
1922: -   flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from

1924:     Options Database:
1925: .     -pc_fieldsplit_diag_use_amat

1927:     Level: intermediate

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

1931: @*/
1932: PetscErrorCode  PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg)
1933: {
1934:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1935:   PetscBool      isfs;

1940:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1941:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1942:   jac->diag_use_amat = flg;
1943:   return(0);
1944: }

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

1949:     Logically Collective on PC

1951:     Input Parameters:
1952: .   pc  - the preconditioner object

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


1958:     Level: intermediate

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

1962: @*/
1963: PetscErrorCode  PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg)
1964: {
1965:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1966:   PetscBool      isfs;

1972:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1973:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1974:   *flg = jac->diag_use_amat;
1975:   return(0);
1976: }

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

1981:     Logically Collective on PC

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

1987:     Options Database:
1988: .     -pc_fieldsplit_off_diag_use_amat

1990:     Level: intermediate

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

1994: @*/
1995: PetscErrorCode  PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg)
1996: {
1997:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1998:   PetscBool      isfs;

2003:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2004:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
2005:   jac->offdiag_use_amat = flg;
2006:   return(0);
2007: }

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

2012:     Logically Collective on PC

2014:     Input Parameters:
2015: .   pc  - the preconditioner object

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


2021:     Level: intermediate

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

2025: @*/
2026: PetscErrorCode  PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg)
2027: {
2028:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
2029:   PetscBool      isfs;

2035:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2036:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
2037:   *flg = jac->offdiag_use_amat;
2038:   return(0);
2039: }



2043: /*@C
2044:     PCFieldSplitSetIS - Sets the exact elements for field

2046:     Logically Collective on PC

2048:     Input Parameters:
2049: +   pc  - the preconditioner context
2050: .   splitname - name of this split, if NULL the number of the split is used
2051: -   is - the index set that defines the vector elements in this field


2054:     Notes:
2055:     Use PCFieldSplitSetFields(), for fields defined by strided types.

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

2060:     Level: intermediate

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

2064: @*/
2065: PetscErrorCode  PCFieldSplitSetIS(PC pc,const char splitname[],IS is)
2066: {

2073:   PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));
2074:   return(0);
2075: }

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

2080:     Logically Collective on PC

2082:     Input Parameters:
2083: +   pc  - the preconditioner context
2084: -   splitname - name of this split

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

2089:     Level: intermediate

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

2093: @*/
2094: PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is)
2095: {

2102:   {
2103:     PC_FieldSplit     *jac  = (PC_FieldSplit*) pc->data;
2104:     PC_FieldSplitLink ilink = jac->head;
2105:     PetscBool         found;

2107:     *is = NULL;
2108:     while (ilink) {
2109:       PetscStrcmp(ilink->splitname, splitname, &found);
2110:       if (found) {
2111:         *is = ilink->is;
2112:         break;
2113:       }
2114:       ilink = ilink->next;
2115:     }
2116:   }
2117:   return(0);
2118: }

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

2123:     Logically Collective on PC

2125:     Input Parameters:
2126: +   pc  - the preconditioner context
2127: -   index - index of this split

2129:     Output Parameter:
2130: -   is - the index set that defines the vector elements in this field

2132:     Level: intermediate

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

2136: @*/
2137: PetscErrorCode PCFieldSplitGetISByIndex(PC pc,PetscInt index,IS *is)
2138: {

2142:   if (index < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",index);
2145:   {
2146:     PC_FieldSplit     *jac  = (PC_FieldSplit*) pc->data;
2147:     PC_FieldSplitLink ilink = jac->head;
2148:     PetscInt          i     = 0;
2149:     if (index >= jac->nsplits) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",index,jac->nsplits);

2151:     while (i < index) {
2152:       ilink = ilink->next;
2153:       ++i;
2154:     }
2155:     PCFieldSplitGetIS(pc,ilink->splitname,is);
2156:   }
2157:   return(0);
2158: }

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

2164:     Logically Collective on PC

2166:     Input Parameters:
2167: +   pc  - the preconditioner context
2168: -   bs - the block size

2170:     Level: intermediate

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

2174: @*/
2175: PetscErrorCode  PCFieldSplitSetBlockSize(PC pc,PetscInt bs)
2176: {

2182:   PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));
2183:   return(0);
2184: }

2186: /*@C
2187:    PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits

2189:    Collective on KSP

2191:    Input Parameter:
2192: .  pc - the preconditioner context

2194:    Output Parameters:
2195: +  n - the number of splits
2196: -  subksp - the array of KSP contexts

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

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

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

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

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


2216:    Level: advanced

2218: .seealso: PCFIELDSPLIT
2219: @*/
2220: PetscErrorCode  PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
2221: {

2227:   PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
2228:   return(0);
2229: }

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

2234:    Collective on KSP

2236:    Input Parameter:
2237: .  pc - the preconditioner context

2239:    Output Parameters:
2240: +  n - the number of splits
2241: -  subksp - the array of KSP contexts

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

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

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

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

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

2260:    Level: advanced

2262: .seealso: PCFIELDSPLIT
2263: @*/
2264: PetscErrorCode  PCFieldSplitSchurGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
2265: {

2271:   PetscUseMethod(pc,"PCFieldSplitSchurGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
2272:   return(0);
2273: }

2275: /*@
2276:     PCFieldSplitSetSchurPre -  Indicates what operator is used to construct the preconditioner for the Schur complement.
2277:       A11 matrix. Otherwise no preconditioner is used.

2279:     Collective on PC

2281:     Input Parameters:
2282: +   pc      - the preconditioner context
2283: .   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
2284:               PC_FIELDSPLIT_SCHUR_PRE_SELFP, and PC_FIELDSPLIT_SCHUR_PRE_FULL
2285: -   userpre - matrix to use for preconditioning, or NULL

2287:     Options Database:
2288: .     -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11. See notes for meaning of various arguments

2290:     Notes:
2291: $    If ptype is
2292: $        a11 then the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner
2293: $             matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix
2294: $        self the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix:
2295: $             The only preconditioner that currently works with this symbolic respresentation matrix object is the PCLSC
2296: $             preconditioner
2297: $        user then the preconditioner for the Schur complement is generated from the user provided matrix (pre argument
2298: $             to this function).
2299: $        selfp then the preconditioning for the Schur complement is generated from an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01
2300: $             This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be
2301: $             lumped before extracting the diagonal using the additional option -fieldsplit_1_mat_schur_complement_ainv_type lump
2302: $        full then the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation computed internally by PCFIELDSPLIT (this is expensive)
2303: $             useful mostly as a test that the Schur complement approach can work for your problem

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

2309:     Level: intermediate

2311: .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType,
2312:           MatSchurComplementSetAinvType(), PCLSC

2314: @*/
2315: PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
2316: {

2321:   PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));
2322:   return(0);
2323: }

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

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

2331:     Logically Collective on PC

2333:     Input Parameters:
2334: .   pc      - the preconditioner context

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

2340:     Level: intermediate

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

2344: @*/
2345: PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
2346: {

2351:   PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));
2352:   return(0);
2353: }

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

2358:     Not collective

2360:     Input Parameter:
2361: .   pc      - the preconditioner context

2363:     Output Parameter:
2364: .   S       - the Schur complement matrix

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

2369:     Level: advanced

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

2373: @*/
2374: PetscErrorCode  PCFieldSplitSchurGetS(PC pc,Mat *S)
2375: {
2377:   const char*    t;
2378:   PetscBool      isfs;
2379:   PC_FieldSplit  *jac;

2383:   PetscObjectGetType((PetscObject)pc,&t);
2384:   PetscStrcmp(t,PCFIELDSPLIT,&isfs);
2385:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
2386:   jac = (PC_FieldSplit*)pc->data;
2387:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
2388:   if (S) *S = jac->schur;
2389:   return(0);
2390: }

2392: /*@
2393:     PCFieldSplitSchurRestoreS -  restores the MatSchurComplement object used by this PC

2395:     Not collective

2397:     Input Parameters:
2398: +   pc      - the preconditioner context
2399: -   S       - the Schur complement matrix

2401:     Level: advanced

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

2405: @*/
2406: PetscErrorCode  PCFieldSplitSchurRestoreS(PC pc,Mat *S)
2407: {
2409:   const char*    t;
2410:   PetscBool      isfs;
2411:   PC_FieldSplit  *jac;

2415:   PetscObjectGetType((PetscObject)pc,&t);
2416:   PetscStrcmp(t,PCFIELDSPLIT,&isfs);
2417:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
2418:   jac = (PC_FieldSplit*)pc->data;
2419:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
2420:   if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten");
2421:   return(0);
2422: }


2425: static PetscErrorCode  PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
2426: {
2427:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

2431:   jac->schurpre = ptype;
2432:   if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
2433:     MatDestroy(&jac->schur_user);
2434:     jac->schur_user = pre;
2435:     PetscObjectReference((PetscObject)jac->schur_user);
2436:   }
2437:   return(0);
2438: }

2440: static PetscErrorCode  PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
2441: {
2442:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

2445:   *ptype = jac->schurpre;
2446:   *pre   = jac->schur_user;
2447:   return(0);
2448: }

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

2453:     Collective on PC

2455:     Input Parameters:
2456: +   pc  - the preconditioner context
2457: -   ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default

2459:     Options Database:
2460: .     -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full


2463:     Level: intermediate

2465:     Notes:
2466:     The FULL factorization is

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

2471:     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,
2472:     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().

2474: $    If A and S are solved exactly
2475: $      *) FULL factorization is a direct solver.
2476: $      *) 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.
2477: $      *) 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.

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

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

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

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

2490: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCFieldSplitSetSchurScale()
2491: @*/
2492: PetscErrorCode  PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype)
2493: {

2498:   PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));
2499:   return(0);
2500: }

2502: static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype)
2503: {
2504:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2507:   jac->schurfactorization = ftype;
2508:   return(0);
2509: }

2511: /*@
2512:     PCFieldSplitSetSchurScale -  Controls the sign flip of S for PC_FIELDSPLIT_SCHUR_FACT_DIAG.

2514:     Collective on PC

2516:     Input Parameters:
2517: +   pc    - the preconditioner context
2518: -   scale - scaling factor for the Schur complement

2520:     Options Database:
2521: .     -pc_fieldsplit_schur_scale - default is -1.0

2523:     Level: intermediate

2525: .seealso: PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurFactType, PCFieldSplitSetSchurScale()
2526: @*/
2527: PetscErrorCode PCFieldSplitSetSchurScale(PC pc,PetscScalar scale)
2528: {

2534:   PetscTryMethod(pc,"PCFieldSplitSetSchurScale_C",(PC,PetscScalar),(pc,scale));
2535:   return(0);
2536: }

2538: static PetscErrorCode PCFieldSplitSetSchurScale_FieldSplit(PC pc,PetscScalar scale)
2539: {
2540:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2543:   jac->schurscale = scale;
2544:   return(0);
2545: }

2547: /*@C
2548:    PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement

2550:    Collective on KSP

2552:    Input Parameter:
2553: .  pc - the preconditioner context

2555:    Output Parameters:
2556: +  A00 - the (0,0) block
2557: .  A01 - the (0,1) block
2558: .  A10 - the (1,0) block
2559: -  A11 - the (1,1) block

2561:    Level: advanced

2563: .seealso: PCFIELDSPLIT
2564: @*/
2565: PetscErrorCode  PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11)
2566: {
2567:   PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;

2571:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach.");
2572:   if (A00) *A00 = jac->pmat[0];
2573:   if (A01) *A01 = jac->B;
2574:   if (A10) *A10 = jac->C;
2575:   if (A11) *A11 = jac->pmat[1];
2576:   return(0);
2577: }

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

2582:     Collective on PC

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

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

2591:     Input Parameters:
2592: +   pc        - the preconditioner context
2593: -   tolerance - the solver tolerance

2595:     Options Database:
2596: .     -pc_fieldsplit_gkb_tol - default is 1e-5

2598:     Level: intermediate

2600: .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBDelay(), PCFieldSplitSetGKBNu(), PCFieldSplitSetGKBMaxit()
2601: @*/
2602: PetscErrorCode PCFieldSplitSetGKBTol(PC pc,PetscReal tolerance)
2603: {

2609:   PetscTryMethod(pc,"PCFieldSplitSetGKBTol_C",(PC,PetscReal),(pc,tolerance));
2610:   return(0);
2611: }

2613: static PetscErrorCode PCFieldSplitSetGKBTol_FieldSplit(PC pc,PetscReal tolerance)
2614: {
2615:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2618:   jac->gkbtol = tolerance;
2619:   return(0);
2620: }


2623: /*@
2624:     PCFieldSplitSetGKBMaxit -  Sets the maximum number of iterations for the generalized Golub-Kahan bidiagonalization
2625:     preconditioner.

2627:     Collective on PC

2629:     Input Parameters:
2630: +   pc     - the preconditioner context
2631: -   maxit  - the maximum number of iterations

2633:     Options Database:
2634: .     -pc_fieldsplit_gkb_maxit - default is 100

2636:     Level: intermediate

2638: .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBDelay(), PCFieldSplitSetGKBTol(), PCFieldSplitSetGKBNu()
2639: @*/
2640: PetscErrorCode PCFieldSplitSetGKBMaxit(PC pc,PetscInt maxit)
2641: {

2647:   PetscTryMethod(pc,"PCFieldSplitSetGKBMaxit_C",(PC,PetscInt),(pc,maxit));
2648:   return(0);
2649: }

2651: static PetscErrorCode PCFieldSplitSetGKBMaxit_FieldSplit(PC pc,PetscInt maxit)
2652: {
2653:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2656:   jac->gkbmaxit = maxit;
2657:   return(0);
2658: }

2660: /*@
2661:     PCFieldSplitSetGKBDelay -  Sets the delay in the lower bound error estimate in the generalized Golub-Kahan bidiagonalization
2662:     preconditioner.

2664:     Collective on PC

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

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

2673:     Input Parameters:
2674: +   pc     - the preconditioner context
2675: -   delay  - the delay window in the lower bound estimate

2677:     Options Database:
2678: .     -pc_fieldsplit_gkb_delay - default is 5

2680:     Level: intermediate

2682: .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBNu(), PCFieldSplitSetGKBTol(), PCFieldSplitSetGKBMaxit()
2683: @*/
2684: PetscErrorCode PCFieldSplitSetGKBDelay(PC pc,PetscInt delay)
2685: {

2691:   PetscTryMethod(pc,"PCFieldSplitSetGKBDelay_C",(PC,PetscInt),(pc,delay));
2692:   return(0);
2693: }

2695: static PetscErrorCode PCFieldSplitSetGKBDelay_FieldSplit(PC pc,PetscInt delay)
2696: {
2697:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2700:   jac->gkbdelay = delay;
2701:   return(0);
2702: }

2704: /*@
2705:     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.

2707:     Collective on PC

2709:     Notes:
2710:     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,
2711:     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
2712:     necessary to find a good balance in between the convergence of the inner and outer loop.

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

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

2718:     Input Parameters:
2719: +   pc     - the preconditioner context
2720: -   nu     - the shift parameter

2722:     Options Database:
2723: .     -pc_fieldsplit_gkb_nu - default is 1

2725:     Level: intermediate

2727: .seealso: PCFIELDSPLIT, PCFieldSplitSetGKBDelay(), PCFieldSplitSetGKBTol(), PCFieldSplitSetGKBMaxit()
2728: @*/
2729: PetscErrorCode PCFieldSplitSetGKBNu(PC pc,PetscReal nu)
2730: {

2736:   PetscTryMethod(pc,"PCFieldSplitSetGKBNu_C",(PC,PetscReal),(pc,nu));
2737:   return(0);
2738: }

2740: static PetscErrorCode PCFieldSplitSetGKBNu_FieldSplit(PC pc,PetscReal nu)
2741: {
2742:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2745:   jac->gkbnu = nu;
2746:   return(0);
2747: }


2750: static PetscErrorCode  PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type)
2751: {
2752:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

2756:   jac->type = type;

2758:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",0);
2759:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);
2760:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);
2761:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);
2762:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",0);
2763:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBTol_C",0);
2764:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBMaxit_C",0);
2765:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBNu_C",0);
2766:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBDelay_C",0);

2768:   if (type == PC_COMPOSITE_SCHUR) {
2769:     pc->ops->apply = PCApply_FieldSplit_Schur;
2770:     pc->ops->view  = PCView_FieldSplit_Schur;

2772:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);
2773:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);
2774:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);
2775:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);
2776:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",PCFieldSplitSetSchurScale_FieldSplit);
2777:   } else if (type == PC_COMPOSITE_GKB){
2778:     pc->ops->apply = PCApply_FieldSplit_GKB;
2779:     pc->ops->view  = PCView_FieldSplit_GKB;

2781:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2782:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBTol_C",PCFieldSplitSetGKBTol_FieldSplit);
2783:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBMaxit_C",PCFieldSplitSetGKBMaxit_FieldSplit);
2784:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBNu_C",PCFieldSplitSetGKBNu_FieldSplit);
2785:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBDelay_C",PCFieldSplitSetGKBDelay_FieldSplit);
2786:   } else {
2787:     pc->ops->apply = PCApply_FieldSplit;
2788:     pc->ops->view  = PCView_FieldSplit;

2790:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2791:   }
2792:   return(0);
2793: }

2795: static PetscErrorCode  PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs)
2796: {
2797:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2800:   if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs);
2801:   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);
2802:   jac->bs = bs;
2803:   return(0);
2804: }

2806: /*@
2807:    PCFieldSplitSetType - Sets the type of fieldsplit preconditioner.

2809:    Collective on PC

2811:    Input Parameter:
2812: +  pc - the preconditioner context
2813: -  type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR

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

2818:    Level: Intermediate

2820: .seealso: PCCompositeSetType()

2822: @*/
2823: PetscErrorCode  PCFieldSplitSetType(PC pc,PCCompositeType type)
2824: {

2829:   PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));
2830:   return(0);
2831: }

2833: /*@
2834:   PCFieldSplitGetType - Gets the type of fieldsplit preconditioner.

2836:   Not collective

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

2841:   Output Parameter:
2842: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR

2844:   Level: Intermediate

2846: .seealso: PCCompositeSetType()
2847: @*/
2848: PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
2849: {
2850:   PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;

2855:   *type = jac->type;
2856:   return(0);
2857: }

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

2862:    Logically Collective

2864:    Input Parameters:
2865: +  pc   - the preconditioner context
2866: -  flg  - boolean indicating whether to use field splits defined by the DM

2868:    Options Database Key:
2869: .  -pc_fieldsplit_dm_splits

2871:    Level: Intermediate

2873: .seealso: PCFieldSplitGetDMSplits()

2875: @*/
2876: PetscErrorCode  PCFieldSplitSetDMSplits(PC pc,PetscBool flg)
2877: {
2878:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
2879:   PetscBool      isfs;

2885:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2886:   if (isfs) {
2887:     jac->dm_splits = flg;
2888:   }
2889:   return(0);
2890: }


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

2896:    Logically Collective

2898:    Input Parameter:
2899: .  pc   - the preconditioner context

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

2904:    Level: Intermediate

2906: .seealso: PCFieldSplitSetDMSplits()

2908: @*/
2909: PetscErrorCode  PCFieldSplitGetDMSplits(PC pc,PetscBool* flg)
2910: {
2911:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
2912:   PetscBool      isfs;

2918:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2919:   if (isfs) {
2920:     if(flg) *flg = jac->dm_splits;
2921:   }
2922:   return(0);
2923: }

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

2928:    Logically Collective

2930:    Input Parameter:
2931: .  pc   - the preconditioner context

2933:    Output Parameter:
2934: .  flg  - boolean indicating whether to detect fields or not

2936:    Level: Intermediate

2938: .seealso: PCFIELDSPLIT, PCFieldSplitSetDetectSaddlePoint()

2940: @*/
2941: PetscErrorCode PCFieldSplitGetDetectSaddlePoint(PC pc,PetscBool *flg)
2942: {
2943:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2946:   *flg = jac->detect;
2947:   return(0);
2948: }

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

2953:    Logically Collective

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

2958:    Input Parameter:
2959: .  pc   - the preconditioner context

2961:    Output Parameter:
2962: .  flg  - boolean indicating whether to detect fields or not

2964:    Options Database Key:
2965: .  -pc_fieldsplit_detect_saddle_point

2967:    Level: Intermediate

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

2971: @*/
2972: PetscErrorCode PCFieldSplitSetDetectSaddlePoint(PC pc,PetscBool flg)
2973: {
2974:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

2978:   jac->detect = flg;
2979:   if (jac->detect) {
2980:     PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);
2981:     PCFieldSplitSetSchurPre(pc,PC_FIELDSPLIT_SCHUR_PRE_SELF,NULL);
2982:   }
2983:   return(0);
2984: }

2986: /* -------------------------------------------------------------------------------------*/
2987: /*MC
2988:    PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
2989:                   fields or groups of fields. See the users manual section "Solving Block Matrices" for more details.

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

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

2997:    Level: intermediate

2999:    Options Database Keys:
3000: +   -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split
3001: .   -pc_fieldsplit_default - automatically add any fields to additional splits that have not
3002:                               been supplied explicitly by -pc_fieldsplit_%d_fields
3003: .   -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields)
3004: .   -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur,gkb> - type of relaxation or factorization splitting
3005: .   -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11; see PCFieldSplitSetSchurPre()
3006: .   -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero diagonal and uses Schur complement with no preconditioner as the solver

3008: .    Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_
3009:      for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields
3010: -    Options prefix for inner solver when using Golub Kahan biadiagonalization preconditioner is -fieldsplit_0_

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

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

3021: $     For the Schur complement preconditioner if J = ( A00 A01 )
3022: $                                                    ( A10 A11 )
3023: $     the preconditioner using full factorization is
3024: $              ( I   -ksp(A00) A01 ) ( inv(A00)     0  ) (     I          0  )
3025: $              ( 0         I       ) (   0      ksp(S) ) ( -A10 ksp(A00)  I  )
3026:      where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_.  S is the Schur complement
3027: $              S = A11 - A10 ksp(A00) A01
3028:      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
3029:      in providing the SECOND split or 1 if not give). For PCFieldSplitGetSubKSP() when field number is 0,
3030:      it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default
3031:      A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() for all the possible ways to construct the preconditioner for S.

3033:      The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above,
3034:      diag gives
3035: $              ( inv(A00)     0   )
3036: $              (   0      -ksp(S) )
3037:      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
3038:      can be turned off with PCFieldSplitSetSchurScale() or by command line -pc_fieldsplit_schur_scale 1.0. The lower factorization is the inverse of
3039: $              (  A00   0 )
3040: $              (  A10   S )
3041:      where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of
3042: $              ( A00 A01 )
3043: $              (  0   S  )
3044:      where again the inverses of A00 and S are applied using KSPs.

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

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

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

3056:    There is a nice discussion of block preconditioners in

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

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

3065:    The generalized Golub-Kahan bidiagonalization preconditioner (gkb) can be applied to symmetric 2x2 block matrices of the shape
3066: $        ( A00  A01 )
3067: $        ( A01' 0   )
3068:    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'.
3069:    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_.

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

3073: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC,
3074:            PCFieldSplitGetSubKSP(), PCFieldSplitSchurGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(),
3075:           MatSchurComplementSetAinvType(), PCFieldSplitSetSchurScale(),
3076:           PCFieldSplitSetDetectSaddlePoint()
3077: M*/

3079: PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
3080: {
3082:   PC_FieldSplit  *jac;

3085:   PetscNewLog(pc,&jac);

3087:   jac->bs                 = -1;
3088:   jac->nsplits            = 0;
3089:   jac->type               = PC_COMPOSITE_MULTIPLICATIVE;
3090:   jac->schurpre           = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
3091:   jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
3092:   jac->schurscale         = -1.0;
3093:   jac->dm_splits          = PETSC_TRUE;
3094:   jac->detect             = PETSC_FALSE;
3095:   jac->gkbtol             = 1e-5;
3096:   jac->gkbdelay           = 5;
3097:   jac->gkbnu              = 1;
3098:   jac->gkbmaxit           = 100;
3099:   jac->gkbmonitor         = PETSC_FALSE;

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

3103:   pc->ops->apply           = PCApply_FieldSplit;
3104:   pc->ops->applytranspose  = PCApplyTranspose_FieldSplit;
3105:   pc->ops->setup           = PCSetUp_FieldSplit;
3106:   pc->ops->reset           = PCReset_FieldSplit;
3107:   pc->ops->destroy         = PCDestroy_FieldSplit;
3108:   pc->ops->setfromoptions  = PCSetFromOptions_FieldSplit;
3109:   pc->ops->view            = PCView_FieldSplit;
3110:   pc->ops->applyrichardson = 0;

3112:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",PCFieldSplitSchurGetSubKSP_FieldSplit);
3113:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
3114:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);
3115:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);
3116:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);
3117:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);
3118:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",PCFieldSplitRestrictIS_FieldSplit);
3119:   return(0);
3120: }