Actual source code: fieldsplit.c

petsc-master 2016-09-27
Report Typos and Errors
 2:  #include <petsc/private/pcimpl.h>
  3: #include <petsc/private/kspimpl.h>    /*  This is needed to provide the appropriate PETSC_EXTERN for KSP_Solve_FS ....*/
 4:  #include <petscdm.h>

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

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

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

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

 36:   /* Only used when Schur complement preconditioning is used */
 37:   Mat                       B;                     /* The (0,1) block */
 38:   Mat                       C;                     /* The (1,0) block */
 39:   Mat                       schur;                 /* The Schur complement S = A11 - A10 A00^{-1} A01, the KSP here, kspinner, is H_1 in [El08] */
 40:   Mat                       schurp;                /* Assembled approximation to S built by MatSchurComplement to be used as a preconditioning matrix when solving with S */
 41:   Mat                       schur_user;            /* User-provided preconditioning matrix for the Schur complement */
 42:   PCFieldSplitSchurPreType  schurpre;              /* Determines which preconditioning matrix is used for the Schur complement */
 43:   PCFieldSplitSchurFactType schurfactorization;
 44:   KSP                       kspschur;              /* The solver for S */
 45:   KSP                       kspupper;              /* The solver for A in the upper diagonal part of the factorization (H_2 in [El08]) */
 46:   PC_FieldSplitLink         head;
 47:   PetscBool                 reset;                  /* indicates PCReset() has been last called on this object, hack */
 48:   PetscBool                 isrestrict;             /* indicates PCFieldSplitRestrictIS() has been last called on this object, hack */
 49:   PetscBool                 suboptionsset;          /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */
 50:   PetscBool                 dm_splits;              /* Whether to use DMCreateFieldDecomposition() whenever possible */
 51:   PetscBool                 diag_use_amat;          /* Whether to extract diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
 52:   PetscBool                 offdiag_use_amat;       /* Whether to extract off-diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
 53: } PC_FieldSplit;

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

 61: /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it.  This way the
 62: * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */
 63: static Mat FieldSplitSchurPre(PC_FieldSplit *jac)
 64: {
 65:   switch (jac->schurpre) {
 66:   case PC_FIELDSPLIT_SCHUR_PRE_SELF: return jac->schur;
 67:   case PC_FIELDSPLIT_SCHUR_PRE_SELFP: return jac->schurp;
 68:   case PC_FIELDSPLIT_SCHUR_PRE_A11: return jac->pmat[1];
 69:   case PC_FIELDSPLIT_SCHUR_PRE_FULL: /* We calculate this and store it in schur_user */
 70:   case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */
 71:   default:
 72:     return jac->schur_user ? jac->schur_user : jac->pmat[1];
 73:   }
 74: }


 77:  #include <petscdraw.h>
 80: static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer)
 81: {
 82:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
 83:   PetscErrorCode    ierr;
 84:   PetscBool         iascii,isdraw;
 85:   PetscInt          i,j;
 86:   PC_FieldSplitLink ilink = jac->head;

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

129:  if (isdraw) {
130:     PetscDraw draw;
131:     PetscReal x,y,w,wd;

133:     PetscViewerDrawGetDraw(viewer,0,&draw);
134:     PetscDrawGetCurrentPoint(draw,&x,&y);
135:     w    = 2*PetscMin(1.0 - x,x);
136:     wd   = w/(jac->nsplits + 1);
137:     x    = x - wd*(jac->nsplits-1)/2.0;
138:     for (i=0; i<jac->nsplits; i++) {
139:       PetscDrawPushCurrentPoint(draw,x,y);
140:       KSPView(ilink->ksp,viewer);
141:       PetscDrawPopCurrentPoint(draw);
142:       x    += wd;
143:       ilink = ilink->next;
144:     }
145:   }
146:   return(0);
147: }

151: static PetscErrorCode PCView_FieldSplit_Schur(PC pc,PetscViewer viewer)
152: {
153:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
154:   PetscErrorCode    ierr;
155:   PetscBool         iascii,isdraw;
156:   PetscInt          i,j;
157:   PC_FieldSplitLink ilink = jac->head;

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

237:     PetscViewerDrawGetDraw(viewer,0,&draw);
238:     PetscDrawGetCurrentPoint(draw,&x,&y);
239:     if (jac->kspupper != jac->head->ksp) cnt++;
240:     w  = 2*PetscMin(1.0 - x,x);
241:     wd = w/(cnt + 1);

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

255:     PetscDrawPushCurrentPoint(draw,x,y);
256:     KSPView(jac->head->ksp,viewer);
257:     PetscDrawPopCurrentPoint(draw);
258:     if (jac->kspupper != jac->head->ksp) {
259:       x   += wd;
260:       PetscDrawPushCurrentPoint(draw,x,y);
261:       KSPView(jac->kspupper,viewer);
262:       PetscDrawPopCurrentPoint(draw);
263:     }
264:     x   += wd;
265:     PetscDrawPushCurrentPoint(draw,x,y);
266:     KSPView(jac->kspschur,viewer);
267:     PetscDrawPopCurrentPoint(draw);
268:   }
269:   return(0);
270: }

274: /* Precondition: jac->bs is set to a meaningful value */
275: static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc)
276: {
278:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
279:   PetscInt       i,nfields,*ifields,nfields_col,*ifields_col;
280:   PetscBool      flg,flg_col;
281:   char           optionname[128],splitname[8],optionname_col[128];

284:   PetscMalloc1(jac->bs,&ifields);
285:   PetscMalloc1(jac->bs,&ifields_col);
286:   for (i=0,flg=PETSC_TRUE;; i++) {
287:     PetscSNPrintf(splitname,sizeof(splitname),"%D",i);
288:     PetscSNPrintf(optionname,sizeof(optionname),"-pc_fieldsplit_%D_fields",i);
289:     PetscSNPrintf(optionname_col,sizeof(optionname_col),"-pc_fieldsplit_%D_fields_col",i);
290:     nfields     = jac->bs;
291:     nfields_col = jac->bs;
292:     PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname,ifields,&nfields,&flg);
293:     PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname_col,ifields_col,&nfields_col,&flg_col);
294:     if (!flg) break;
295:     else if (flg && !flg_col) {
296:       if (!nfields) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields");
297:       PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields);
298:     } else {
299:       if (!nfields || !nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields");
300:       if (nfields != nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Number of row and column fields must match");
301:       PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields_col);
302:     }
303:   }
304:   if (i > 0) {
305:     /* Makes command-line setting of splits take precedence over setting them in code.
306:        Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would
307:        create new splits, which would probably not be what the user wanted. */
308:     jac->splitdefined = PETSC_TRUE;
309:   }
310:   PetscFree(ifields);
311:   PetscFree(ifields_col);
312:   return(0);
313: }

317: static PetscErrorCode PCFieldSplitSetDefaults(PC pc)
318: {
319:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
320:   PetscErrorCode    ierr;
321:   PC_FieldSplitLink ilink = jac->head;
322:   PetscBool         fieldsplit_default = PETSC_FALSE,stokes = PETSC_FALSE,coupling = PETSC_FALSE;
323:   PetscInt          i;

326:   /*
327:    Kinda messy, but at least this now uses DMCreateFieldDecomposition() even with jac->reset.
328:    Should probably be rewritten.
329:    */
330:   if (!ilink || jac->reset) {
331:     PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);
332:     PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_coupling",&coupling,NULL);
333:     if (pc->dm && jac->dm_splits && !stokes && !coupling) {
334:       PetscInt  numFields, f, i, j;
335:       char      **fieldNames;
336:       IS        *fields;
337:       DM        *dms;
338:       DM        subdm[128];
339:       PetscBool flg;

341:       DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms);
342:       /* Allow the user to prescribe the splits */
343:       for (i = 0, flg = PETSC_TRUE;; i++) {
344:         PetscInt ifields[128];
345:         IS       compField;
346:         char     optionname[128], splitname[8];
347:         PetscInt nfields = numFields;

349:         PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%D_fields", i);
350:         PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg);
351:         if (!flg) break;
352:         if (numFields > 128) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot currently support %d > 128 fields", numFields);
353:         DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i]);
354:         if (nfields == 1) {
355:           PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField);
356:         } else {
357:           PetscSNPrintf(splitname, sizeof(splitname), "%D", i);
358:           PCFieldSplitSetIS(pc, splitname, compField);
359:         }
360:         ISDestroy(&compField);
361:         for (j = 0; j < nfields; ++j) {
362:           f    = ifields[j];
363:           PetscFree(fieldNames[f]);
364:           ISDestroy(&fields[f]);
365:         }
366:       }
367:       if (i == 0) {
368:         for (f = 0; f < numFields; ++f) {
369:           PCFieldSplitSetIS(pc, fieldNames[f], fields[f]);
370:           PetscFree(fieldNames[f]);
371:           ISDestroy(&fields[f]);
372:         }
373:       } else {
374:         for (j=0; j<numFields; j++) {
375:           DMDestroy(dms+j);
376:         }
377:         PetscFree(dms);
378:         PetscMalloc1(i, &dms);
379:         for (j = 0; j < i; ++j) dms[j] = subdm[j];
380:       }
381:       PetscFree(fieldNames);
382:       PetscFree(fields);
383:       if (dms) {
384:         PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n");
385:         for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) {
386:           const char *prefix;
387:           PetscObjectGetOptionsPrefix((PetscObject)(ilink->ksp),&prefix);
388:           PetscObjectSetOptionsPrefix((PetscObject)(dms[i]), prefix);
389:           KSPSetDM(ilink->ksp, dms[i]);
390:           KSPSetDMActive(ilink->ksp, PETSC_FALSE);
391:           PetscObjectIncrementTabLevel((PetscObject)dms[i],(PetscObject)ilink->ksp,0);
392:           DMDestroy(&dms[i]);
393:         }
394:         PetscFree(dms);
395:       }
396:     } else {
397:       if (jac->bs <= 0) {
398:         if (pc->pmat) {
399:           MatGetBlockSize(pc->pmat,&jac->bs);
400:         } else jac->bs = 1;
401:       }

403:       if (stokes) {
404:         IS       zerodiags,rest;
405:         PetscInt nmin,nmax;

407:         MatGetOwnershipRange(pc->mat,&nmin,&nmax);
408:         MatFindZeroDiagonals(pc->mat,&zerodiags);
409:         ISComplement(zerodiags,nmin,nmax,&rest);
410:         if (jac->reset) {
411:           if (!jac->head) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"At reset jac must have head from previous setup");
412:           jac->head->is       = rest;
413:           jac->head->next->is = zerodiags;
414:         } else {
415:           PCFieldSplitSetIS(pc,"0",rest);
416:           PCFieldSplitSetIS(pc,"1",zerodiags);
417:         }
418:         ISDestroy(&zerodiags);
419:         ISDestroy(&rest);
420:       } else if (coupling) {
421:         IS       coupling,rest;
422:         PetscInt nmin,nmax;

424:         MatGetOwnershipRange(pc->mat,&nmin,&nmax);
425:         MatFindOffBlockDiagonalEntries(pc->mat,&coupling);
426:         ISCreateStride(PetscObjectComm((PetscObject)pc->mat),nmax-nmin,nmin,1,&rest);
427:         ISSetIdentity(rest);
428:         if (jac->reset) {
429:           if (!jac->head) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"At reset jac must have head from previous setup");
430:           jac->head->is       = rest;
431:           jac->head->next->is = coupling;
432:         } else {
433:           PCFieldSplitSetIS(pc,"0",rest);
434:           PCFieldSplitSetIS(pc,"1",coupling);
435:         }
436:         ISDestroy(&coupling);
437:         ISDestroy(&rest);
438:       } else {
439:         if (jac->reset && !jac->isrestrict) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cases not yet handled when PCReset() was used");
440:         PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&fieldsplit_default,NULL);
441:         if (!fieldsplit_default) {
442:           /* Allow user to set fields from command line,  if bs was known at the time of PCSetFromOptions_FieldSplit()
443:            then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */
444:           PCFieldSplitSetRuntimeSplits_Private(pc);
445:           if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
446:         }
447:         if ((fieldsplit_default || !jac->splitdefined) && !jac->isrestrict) {
448:           PetscInfo(pc,"Using default splitting of fields\n");
449:           for (i=0; i<jac->bs; i++) {
450:             char splitname[8];
451:             PetscSNPrintf(splitname,sizeof(splitname),"%D",i);
452:             PCFieldSplitSetFields(pc,splitname,1,&i,&i);
453:           }
454:           jac->defaultsplit = PETSC_TRUE;
455:         }
456:       }
457:     }
458:   } else if (jac->nsplits == 1) {
459:     if (ilink->is) {
460:       IS       is2;
461:       PetscInt nmin,nmax;

463:       MatGetOwnershipRange(pc->mat,&nmin,&nmax);
464:       ISComplement(ilink->is,nmin,nmax,&is2);
465:       PCFieldSplitSetIS(pc,"1",is2);
466:       ISDestroy(&is2);
467:     } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()");
468:   }

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

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

478: static PetscErrorCode PCSetUp_FieldSplit(PC pc)
479: {
480:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
481:   PetscErrorCode    ierr;
482:   PC_FieldSplitLink ilink;
483:   PetscInt          i,nsplit;
484:   PetscBool         sorted, sorted_col;

487:   PCFieldSplitSetDefaults(pc);
488:   nsplit = jac->nsplits;
489:   ilink  = jac->head;

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

495:     jac->issetup = PETSC_TRUE;

497:     /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */
498:     if (jac->defaultsplit || !ilink->is) {
499:       if (jac->bs <= 0) jac->bs = nsplit;
500:     }
501:     bs     = jac->bs;
502:     MatGetOwnershipRange(pc->pmat,&rstart,&rend);
503:     nslots = (rend - rstart)/bs;
504:     for (i=0; i<nsplit; i++) {
505:       if (jac->defaultsplit) {
506:         ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+i,nsplit,&ilink->is);
507:         ISDuplicate(ilink->is,&ilink->is_col);
508:       } else if (!ilink->is) {
509:         if (ilink->nfields > 1) {
510:           PetscInt *ii,*jj,j,k,nfields = ilink->nfields,*fields = ilink->fields,*fields_col = ilink->fields_col;
511:           PetscMalloc1(ilink->nfields*nslots,&ii);
512:           PetscMalloc1(ilink->nfields*nslots,&jj);
513:           for (j=0; j<nslots; j++) {
514:             for (k=0; k<nfields; k++) {
515:               ii[nfields*j + k] = rstart + bs*j + fields[k];
516:               jj[nfields*j + k] = rstart + bs*j + fields_col[k];
517:             }
518:           }
519:           ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);
520:           ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,jj,PETSC_OWN_POINTER,&ilink->is_col);
521:           ISSetBlockSize(ilink->is, nfields);
522:           ISSetBlockSize(ilink->is_col, nfields);
523:         } else {
524:           ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields[0],bs,&ilink->is);
525:           ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields_col[0],bs,&ilink->is_col);
526:        }
527:       }
528:       ISSorted(ilink->is,&sorted);
529:       if (ilink->is_col) { ISSorted(ilink->is_col,&sorted_col); }
530:       if (!sorted || !sorted_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split");
531:       ilink = ilink->next;
532:     }
533:   }

535:   ilink = jac->head;
536:   if (!jac->pmat) {
537:     Vec xtmp;

539:     MatCreateVecs(pc->pmat,&xtmp,NULL);
540:     PetscMalloc1(nsplit,&jac->pmat);
541:     PetscMalloc2(nsplit,&jac->x,nsplit,&jac->y);
542:     for (i=0; i<nsplit; i++) {
543:       MatNullSpace sp;

545:       /* Check for preconditioning matrix attached to IS */
546:       PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &jac->pmat[i]);
547:       if (jac->pmat[i]) {
548:         PetscObjectReference((PetscObject) jac->pmat[i]);
549:         if (jac->type == PC_COMPOSITE_SCHUR) {
550:           jac->schur_user = jac->pmat[i];

552:           PetscObjectReference((PetscObject) jac->schur_user);
553:         }
554:       } else {
555:         const char *prefix;
556:         MatGetSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->pmat[i]);
557:         KSPGetOptionsPrefix(ilink->ksp,&prefix);
558:         MatSetOptionsPrefix(jac->pmat[i],prefix);
559:         MatViewFromOptions(jac->pmat[i],NULL,"-mat_view");
560:       }
561:       /* create work vectors for each split */
562:       MatCreateVecs(jac->pmat[i],&jac->x[i],&jac->y[i]);
563:       ilink->x = jac->x[i]; ilink->y = jac->y[i]; ilink->z = NULL;
564:       /* compute scatter contexts needed by multiplicative versions and non-default splits */
565:       VecScatterCreate(xtmp,ilink->is,jac->x[i],NULL,&ilink->sctx);
566:       PetscObjectQuery((PetscObject) ilink->is, "nearnullspace", (PetscObject*) &sp);
567:       if (sp) {
568:         MatSetNearNullSpace(jac->pmat[i], sp);
569:       }
570:       ilink = ilink->next;
571:     }
572:     VecDestroy(&xtmp);
573:   } else {
574:     for (i=0; i<nsplit; i++) {
575:       Mat pmat;

577:       /* Check for preconditioning matrix attached to IS */
578:       PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &pmat);
579:       if (!pmat) {
580:         MatGetSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->pmat[i]);
581:       }
582:       ilink = ilink->next;
583:     }
584:   }
585:   if (jac->diag_use_amat) {
586:     ilink = jac->head;
587:     if (!jac->mat) {
588:       PetscMalloc1(nsplit,&jac->mat);
589:       for (i=0; i<nsplit; i++) {
590:         MatGetSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i]);
591:         ilink = ilink->next;
592:       }
593:     } else {
594:       for (i=0; i<nsplit; i++) {
595:         if (jac->mat[i]) {MatGetSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->mat[i]);}
596:         ilink = ilink->next;
597:       }
598:     }
599:   } else {
600:     jac->mat = jac->pmat;
601:   }

603:   /* Check for null space attached to IS */
604:   ilink = jac->head;
605:   for (i=0; i<nsplit; i++) {
606:     MatNullSpace sp;

608:     PetscObjectQuery((PetscObject) ilink->is, "nullspace", (PetscObject*) &sp);
609:     if (sp) {
610:       MatSetNullSpace(jac->mat[i], sp);
611:     }
612:     ilink = ilink->next;
613:   }

615:   if (jac->type != PC_COMPOSITE_ADDITIVE  && jac->type != PC_COMPOSITE_SCHUR) {
616:     /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */
617:     /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */
618:     ilink = jac->head;
619:     if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) {
620:       /* 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 */
621:       if (!jac->Afield) {
622:         PetscCalloc1(nsplit,&jac->Afield);
623:         if (jac->offdiag_use_amat) {
624:           MatGetSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);
625:         } else {
626:           MatGetSubMatrix(pc->pmat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);
627:         }
628:       } else {
629:         if (jac->offdiag_use_amat) {
630:           MatGetSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_REUSE_MATRIX,&jac->Afield[1]);
631:         } else {
632:           MatGetSubMatrix(pc->pmat,ilink->next->is,ilink->is,MAT_REUSE_MATRIX,&jac->Afield[1]);
633:         }
634:       }
635:     } else {
636:       if (!jac->Afield) {
637:         PetscMalloc1(nsplit,&jac->Afield);
638:         for (i=0; i<nsplit; i++) {
639:           if (jac->offdiag_use_amat) {
640:             MatGetSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
641:           } else {
642:             MatGetSubMatrix(pc->pmat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
643:           }
644:           ilink = ilink->next;
645:         }
646:       } else {
647:         for (i=0; i<nsplit; i++) {
648:           if (jac->offdiag_use_amat) {
649:             MatGetSubMatrix(pc->mat,ilink->is,NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);
650:           } else {
651:             MatGetSubMatrix(pc->pmat,ilink->is,NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);
652:           }
653:           ilink = ilink->next;
654:         }
655:       }
656:     }
657:   }

659:   if (jac->type == PC_COMPOSITE_SCHUR) {
660:     IS          ccis;
661:     PetscInt    rstart,rend;
662:     char        lscname[256];
663:     PetscObject LSC_L;

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

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

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

674:       MatSchurComplementGetKSP(jac->schur, &kspInner);
675:       ilink = jac->head;
676:       ISComplement(ilink->is_col,rstart,rend,&ccis);
677:       if (jac->offdiag_use_amat) {
678:         MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);
679:       } else {
680:         MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);
681:       }
682:       ISDestroy(&ccis);
683:       ilink = ilink->next;
684:       ISComplement(ilink->is_col,rstart,rend,&ccis);
685:       if (jac->offdiag_use_amat) {
686:         MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);
687:       } else {
688:         MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);
689:       }
690:       ISDestroy(&ccis);
691:       MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
692:       if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
693:         MatDestroy(&jac->schurp);
694:         MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);
695:       }
696:       if (kspA != kspInner) {
697:         KSPSetOperators(kspA,jac->mat[0],jac->pmat[0]);
698:       }
699:       if (kspUpper != kspA) {
700:         KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0]);
701:       }
702:       KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
703:     } else {
704:       const char   *Dprefix;
705:       char         schurprefix[256], schurmatprefix[256];
706:       char         schurtestoption[256];
707:       MatNullSpace sp;
708:       PetscBool    flg;

710:       /* extract the A01 and A10 matrices */
711:       ilink = jac->head;
712:       ISComplement(ilink->is_col,rstart,rend,&ccis);
713:       if (jac->offdiag_use_amat) {
714:         MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
715:       } else {
716:         MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
717:       }
718:       ISDestroy(&ccis);
719:       ilink = ilink->next;
720:       ISComplement(ilink->is_col,rstart,rend,&ccis);
721:       if (jac->offdiag_use_amat) {
722:         MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
723:       } else {
724:         MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
725:       }
726:       ISDestroy(&ccis);

728:       /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
729:       MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);
730:       MatSetType(jac->schur,MATSCHURCOMPLEMENT);
731:       MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
732:       PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
733:       /* Note that the inner KSP is NOT going to inherit this prefix, and if it did, it would be reset just below.  Is that what we want? */
734:       MatSetOptionsPrefix(jac->schur,schurmatprefix);
735:       MatSetFromOptions(jac->schur);
736:       MatGetNullSpace(jac->mat[1], &sp);
737:       if (sp) {
738:         MatSetNullSpace(jac->schur, sp);
739:       }

741:       PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);
742:       PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
743:       if (flg) {
744:         DM  dmInner;
745:         KSP kspInner;

747:         MatSchurComplementGetKSP(jac->schur, &kspInner);
748:         PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
749:         /* Indent this deeper to emphasize the "inner" nature of this solver. */
750:         PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject) pc, 2);
751:         KSPSetOptionsPrefix(kspInner, schurprefix);

753:         /* Set DM for new solver */
754:         KSPGetDM(jac->head->ksp, &dmInner);
755:         KSPSetDM(kspInner, dmInner);
756:         KSPSetDMActive(kspInner, PETSC_FALSE);
757:       } else {
758:          /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or
759:           * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact,
760:           * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for
761:           * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make
762:           * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used
763:           * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */
764:         KSPSetType(jac->head->ksp,KSPGMRES);
765:         MatSchurComplementSetKSP(jac->schur,jac->head->ksp);
766:       }
767:       KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0]);
768:       KSPSetFromOptions(jac->head->ksp);
769:       MatSetFromOptions(jac->schur);

771:       PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);
772:       PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
773:       if (flg) {
774:         DM dmInner;

776:         PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
777:         KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);
778:         KSPSetErrorIfNotConverged(jac->kspupper,pc->erroriffailure);
779:         KSPSetOptionsPrefix(jac->kspupper, schurprefix);
780:         KSPGetDM(jac->head->ksp, &dmInner);
781:         KSPSetDM(jac->kspupper, dmInner);
782:         KSPSetDMActive(jac->kspupper, PETSC_FALSE);
783:         KSPSetFromOptions(jac->kspupper);
784:         KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0]);
785:         VecDuplicate(jac->head->x, &jac->head->z);
786:       } else {
787:         jac->kspupper = jac->head->ksp;
788:         PetscObjectReference((PetscObject) jac->head->ksp);
789:       }

791:       if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
792:         MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);
793:       }
794:       KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);
795:       KSPSetErrorIfNotConverged(jac->kspschur,pc->erroriffailure);
796:       PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);
797:       PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);
798:       if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
799:         PC pcschur;
800:         KSPGetPC(jac->kspschur,&pcschur);
801:         PCSetType(pcschur,PCNONE);
802:         /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
803:       } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) {
804:         MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user);
805:       }
806:       KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
807:       KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);
808:       KSPSetOptionsPrefix(jac->kspschur,         Dprefix);
809:       /* propogate DM */
810:       {
811:         DM sdm;
812:         KSPGetDM(jac->head->next->ksp, &sdm);
813:         if (sdm) {
814:           KSPSetDM(jac->kspschur, sdm);
815:           KSPSetDMActive(jac->kspschur, PETSC_FALSE);
816:         }
817:       }
818:       /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
819:       /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */
820:       KSPSetFromOptions(jac->kspschur);
821:     }

823:     /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */
824:     PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);
825:     PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);
826:     if (!LSC_L) {PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);}
827:     if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);}
828:     PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);
829:     PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);
830:     if (!LSC_L) {PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);}
831:     if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);}
832:   } else {
833:     /* set up the individual splits' PCs */
834:     i     = 0;
835:     ilink = jac->head;
836:     while (ilink) {
837:       KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i]);
838:       /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
839:       if (!jac->suboptionsset) {KSPSetFromOptions(ilink->ksp);}
840:       i++;
841:       ilink = ilink->next;
842:     }
843:   }

845:   jac->suboptionsset = PETSC_TRUE;
846:   return(0);
847: }

849: #define FieldSplitSplitSolveAdd(ilink,xx,yy) \
850:   (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
851:    VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
852:    PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\
853:    KSPSolve(ilink->ksp,ilink->x,ilink->y) ||                               \
854:    PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\
855:    VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) ||  \
856:    VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE))

860: static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y)
861: {
862:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
863:   PetscErrorCode    ierr;
864:   PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
865:   KSP               kspA   = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;

868:   switch (jac->schurfactorization) {
869:   case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
870:     /* [A00 0; 0 -S], positive definite, suitable for MINRES */
871:     VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
872:     VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
873:     VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
874:     PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
875:     KSPSolve(kspA,ilinkA->x,ilinkA->y);
876:     PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
877:     VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
878:     VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
879:     PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
880:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
881:     PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
882:     VecScale(ilinkD->y,-1.);
883:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
884:     VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
885:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
886:     break;
887:   case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
888:     /* [A00 0; A10 S], suitable for left preconditioning */
889:     VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
890:     VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
891:     PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
892:     KSPSolve(kspA,ilinkA->x,ilinkA->y);
893:     PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
894:     MatMult(jac->C,ilinkA->y,ilinkD->x);
895:     VecScale(ilinkD->x,-1.);
896:     VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
897:     VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
898:     VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
899:     PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
900:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
901:     PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
902:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
903:     VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
904:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
905:     break;
906:   case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
907:     /* [A00 A01; 0 S], suitable for right preconditioning */
908:     VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
909:     VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
910:     PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
911:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
912:     PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);    MatMult(jac->B,ilinkD->y,ilinkA->x);
913:     VecScale(ilinkA->x,-1.);
914:     VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
915:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
916:     VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
917:     PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
918:     KSPSolve(kspA,ilinkA->x,ilinkA->y);
919:     PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
920:     VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
921:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
922:     VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
923:     break;
924:   case PC_FIELDSPLIT_SCHUR_FACT_FULL:
925:     /* [1 0; A10 A00^{-1} 1] [A00 0; 0 S] [1 A00^{-1}A01; 0 1], an exact solve if applied exactly, needs one extra solve with A */
926:     VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
927:     VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
928:     PetscLogEventBegin(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);
929:     KSPSolve(kspLower,ilinkA->x,ilinkA->y);
930:     PetscLogEventEnd(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);
931:     MatMult(jac->C,ilinkA->y,ilinkD->x);
932:     VecScale(ilinkD->x,-1.0);
933:     VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
934:     VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);

936:     PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
937:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
938:     PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
939:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
940:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);

942:     if (kspUpper == kspA) {
943:       MatMult(jac->B,ilinkD->y,ilinkA->y);
944:       VecAXPY(ilinkA->x,-1.0,ilinkA->y);
945:       PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
946:       KSPSolve(kspA,ilinkA->x,ilinkA->y);
947:       PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
948:     } else {
949:       PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
950:       KSPSolve(kspA,ilinkA->x,ilinkA->y);
951:       PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
952:       MatMult(jac->B,ilinkD->y,ilinkA->x);
953:       PetscLogEventBegin(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);
954:       KSPSolve(kspUpper,ilinkA->x,ilinkA->z);
955:       PetscLogEventEnd(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);
956:       VecAXPY(ilinkA->y,-1.0,ilinkA->z);
957:     }
958:     VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
959:     VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
960:   }
961:   return(0);
962: }

966: static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y)
967: {
968:   PC_FieldSplit      *jac = (PC_FieldSplit*)pc->data;
969:   PetscErrorCode     ierr;
970:   PC_FieldSplitLink  ilink = jac->head;
971:   PetscInt           cnt,bs;
972:   KSPConvergedReason reason;

975:   if (jac->type == PC_COMPOSITE_ADDITIVE) {
976:     if (jac->defaultsplit) {
977:       VecGetBlockSize(x,&bs);
978:       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);
979:       VecGetBlockSize(y,&bs);
980:       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);
981:       VecStrideGatherAll(x,jac->x,INSERT_VALUES);
982:       while (ilink) {
983:         PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
984:         KSPSolve(ilink->ksp,ilink->x,ilink->y);
985:         PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
986:         KSPGetConvergedReason(ilink->ksp,&reason);
987:         if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
988:           pc->failedreason = PC_SUBPC_ERROR;
989:         }
990:         ilink = ilink->next;
991:       }
992:       VecStrideScatterAll(jac->y,y,INSERT_VALUES);
993:     } else {
994:       VecSet(y,0.0);
995:       while (ilink) {
996:         FieldSplitSplitSolveAdd(ilink,x,y);
997:         KSPGetConvergedReason(ilink->ksp,&reason);
998:         if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
999:           pc->failedreason = PC_SUBPC_ERROR;
1000:         }
1001:         ilink = ilink->next;
1002:       }
1003:     }
1004:   } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) {
1005:     VecSet(y,0.0);
1006:     /* solve on first block for first block variables */
1007:     VecScatterBegin(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
1008:     VecScatterEnd(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
1009:     PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1010:     KSPSolve(ilink->ksp,ilink->x,ilink->y);
1011:     PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1012:     KSPGetConvergedReason(ilink->ksp,&reason);
1013:     if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1014:       pc->failedreason = PC_SUBPC_ERROR;
1015:     }
1016:     VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1017:     VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);

1019:     /* compute the residual only onto second block variables using first block variables */
1020:     MatMult(jac->Afield[1],ilink->y,ilink->next->x);
1021:     ilink = ilink->next;
1022:     VecScale(ilink->x,-1.0);
1023:     VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1024:     VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);

1026:     /* solve on second block variables */
1027:     PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1028:     KSPSolve(ilink->ksp,ilink->x,ilink->y);
1029:     PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1030:     KSPGetConvergedReason(ilink->ksp,&reason);
1031:     if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1032:       pc->failedreason = PC_SUBPC_ERROR;
1033:     }
1034:     VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1035:     VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1036:   } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1037:     if (!jac->w1) {
1038:       VecDuplicate(x,&jac->w1);
1039:       VecDuplicate(x,&jac->w2);
1040:     }
1041:     VecSet(y,0.0);
1042:     FieldSplitSplitSolveAdd(ilink,x,y);
1043:     KSPGetConvergedReason(ilink->ksp,&reason);
1044:     if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1045:       pc->failedreason = PC_SUBPC_ERROR;
1046:     }
1047:     cnt  = 1;
1048:     while (ilink->next) {
1049:       ilink = ilink->next;
1050:       /* compute the residual only over the part of the vector needed */
1051:       MatMult(jac->Afield[cnt++],y,ilink->x);
1052:       VecScale(ilink->x,-1.0);
1053:       VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1054:       VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1055:       PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1056:       KSPSolve(ilink->ksp,ilink->x,ilink->y);
1057:       PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1058:       KSPGetConvergedReason(ilink->ksp,&reason);
1059:       if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1060:         pc->failedreason = PC_SUBPC_ERROR;
1061:       }
1062:       VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1063:       VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1064:     }
1065:     if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1066:       cnt -= 2;
1067:       while (ilink->previous) {
1068:         ilink = ilink->previous;
1069:         /* compute the residual only over the part of the vector needed */
1070:         MatMult(jac->Afield[cnt--],y,ilink->x);
1071:         VecScale(ilink->x,-1.0);
1072:         VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1073:         VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1074:         PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1075:         KSPSolve(ilink->ksp,ilink->x,ilink->y);
1076:         PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1077:         KSPGetConvergedReason(ilink->ksp,&reason);
1078:         if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1079:           pc->failedreason = PC_SUBPC_ERROR;
1080:         }
1081:         VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1082:         VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1083:       }
1084:     }
1085:   } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type);
1086:   return(0);
1087: }

1089: #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \
1090:   (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1091:    VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1092:    PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \
1093:    KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) ||                  \
1094:    PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \
1095:    VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \
1096:    VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE))

1100: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y)
1101: {
1102:   PC_FieldSplit      *jac = (PC_FieldSplit*)pc->data;
1103:   PetscErrorCode     ierr;
1104:   PC_FieldSplitLink  ilink = jac->head;
1105:   PetscInt           bs;
1106:   KSPConvergedReason reason;

1109:   if (jac->type == PC_COMPOSITE_ADDITIVE) {
1110:     if (jac->defaultsplit) {
1111:       VecGetBlockSize(x,&bs);
1112:       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);
1113:       VecGetBlockSize(y,&bs);
1114:       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);
1115:       VecStrideGatherAll(x,jac->x,INSERT_VALUES);
1116:       while (ilink) {
1117:         PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1118:         KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);
1119:         PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1120:         KSPGetConvergedReason(ilink->ksp,&reason);
1121:         if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1122:           pc->failedreason = PC_SUBPC_ERROR;
1123:         }
1124:         ilink = ilink->next;
1125:       }
1126:       VecStrideScatterAll(jac->y,y,INSERT_VALUES);
1127:     } else {
1128:       VecSet(y,0.0);
1129:       while (ilink) {
1130:         FieldSplitSplitSolveAddTranspose(ilink,x,y);
1131:         KSPGetConvergedReason(ilink->ksp,&reason);
1132:         if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1133:           pc->failedreason = PC_SUBPC_ERROR;
1134:         }
1135:         ilink = ilink->next;
1136:       }
1137:     }
1138:   } else {
1139:     if (!jac->w1) {
1140:       VecDuplicate(x,&jac->w1);
1141:       VecDuplicate(x,&jac->w2);
1142:     }
1143:     VecSet(y,0.0);
1144:     if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1145:       FieldSplitSplitSolveAddTranspose(ilink,x,y);
1146:       KSPGetConvergedReason(ilink->ksp,&reason);
1147:       if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1148:         pc->failedreason = PC_SUBPC_ERROR;
1149:       }
1150:       while (ilink->next) {
1151:         ilink = ilink->next;
1152:         MatMultTranspose(pc->mat,y,jac->w1);
1153:         VecWAXPY(jac->w2,-1.0,jac->w1,x);
1154:         FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1155:       }
1156:       while (ilink->previous) {
1157:         ilink = ilink->previous;
1158:         MatMultTranspose(pc->mat,y,jac->w1);
1159:         VecWAXPY(jac->w2,-1.0,jac->w1,x);
1160:         FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1161:       }
1162:     } else {
1163:       while (ilink->next) {   /* get to last entry in linked list */
1164:         ilink = ilink->next;
1165:       }
1166:       FieldSplitSplitSolveAddTranspose(ilink,x,y);
1167:       KSPGetConvergedReason(ilink->ksp,&reason);
1168:       if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1169:         pc->failedreason = PC_SUBPC_ERROR;
1170:       }
1171:       while (ilink->previous) {
1172:         ilink = ilink->previous;
1173:         MatMultTranspose(pc->mat,y,jac->w1);
1174:         VecWAXPY(jac->w2,-1.0,jac->w1,x);
1175:         FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1176:       }
1177:     }
1178:   }
1179:   return(0);
1180: }

1184: static PetscErrorCode PCReset_FieldSplit(PC pc)
1185: {
1186:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1187:   PetscErrorCode    ierr;
1188:   PC_FieldSplitLink ilink = jac->head,next;

1191:   while (ilink) {
1192:     KSPReset(ilink->ksp);
1193:     VecDestroy(&ilink->x);
1194:     VecDestroy(&ilink->y);
1195:     VecDestroy(&ilink->z);
1196:     VecScatterDestroy(&ilink->sctx);
1197:     if (!ilink->is_orig) {             /* save the original IS */
1198:       PetscObjectReference((PetscObject)ilink->is);
1199:       ilink->is_orig = ilink->is;
1200:     }
1201:     ISDestroy(&ilink->is);
1202:     ISDestroy(&ilink->is_col);
1203:     next  = ilink->next;
1204:     ilink = next;
1205:   }
1206:   PetscFree2(jac->x,jac->y);
1207:   if (jac->mat && jac->mat != jac->pmat) {
1208:     MatDestroyMatrices(jac->nsplits,&jac->mat);
1209:   } else if (jac->mat) {
1210:     jac->mat = NULL;
1211:   }
1212:   if (jac->pmat) {MatDestroyMatrices(jac->nsplits,&jac->pmat);}
1213:   if (jac->Afield) {MatDestroyMatrices(jac->nsplits,&jac->Afield);}
1214:   VecDestroy(&jac->w1);
1215:   VecDestroy(&jac->w2);
1216:   MatDestroy(&jac->schur);
1217:   MatDestroy(&jac->schurp);
1218:   MatDestroy(&jac->schur_user);
1219:   KSPDestroy(&jac->kspschur);
1220:   KSPDestroy(&jac->kspupper);
1221:   MatDestroy(&jac->B);
1222:   MatDestroy(&jac->C);
1223:   jac->reset = PETSC_TRUE;
1224:   jac->isrestrict = PETSC_FALSE;
1225:   return(0);
1226: }

1230: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
1231: {
1232:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1233:   PetscErrorCode    ierr;
1234:   PC_FieldSplitLink ilink = jac->head,next;

1237:   PCReset_FieldSplit(pc);
1238:   while (ilink) {
1239:     KSPDestroy(&ilink->ksp);
1240:     ISDestroy(&ilink->is_orig);
1241:     next  = ilink->next;
1242:     PetscFree(ilink->splitname);
1243:     PetscFree(ilink->fields);
1244:     PetscFree(ilink->fields_col);
1245:     PetscFree(ilink);
1246:     ilink = next;
1247:   }
1248:   PetscFree2(jac->x,jac->y);
1249:   PetscFree(pc->data);
1250:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);
1251:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);
1252:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);
1253:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);
1254:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);
1255:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);
1256:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);
1257:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);
1258:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",NULL);
1259:   return(0);
1260: }

1264: static PetscErrorCode PCSetFromOptions_FieldSplit(PetscOptionItems *PetscOptionsObject,PC pc)
1265: {
1266:   PetscErrorCode  ierr;
1267:   PetscInt        bs;
1268:   PetscBool       flg,stokes = PETSC_FALSE;
1269:   PC_FieldSplit   *jac = (PC_FieldSplit*)pc->data;
1270:   PCCompositeType ctype;

1273:   PetscOptionsHead(PetscOptionsObject,"FieldSplit options");
1274:   PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);
1275:   PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);
1276:   if (flg) {
1277:     PCFieldSplitSetBlockSize(pc,bs);
1278:   }
1279:   jac->diag_use_amat = pc->useAmat;
1280:   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);
1281:   jac->offdiag_use_amat = pc->useAmat;
1282:   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);
1283:   /* FIXME: No programmatic equivalent to the following. */
1284:   PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);
1285:   if (stokes) {
1286:     PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);
1287:     jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF;
1288:   }

1290:   PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);
1291:   if (flg) {
1292:     PCFieldSplitSetType(pc,ctype);
1293:   }
1294:   /* Only setup fields once */
1295:   if ((jac->bs > 0) && (jac->nsplits == 0)) {
1296:     /* only allow user to set fields from command line if bs is already known.
1297:        otherwise user can set them in PCFieldSplitSetDefaults() */
1298:     PCFieldSplitSetRuntimeSplits_Private(pc);
1299:     if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
1300:   }
1301:   if (jac->type == PC_COMPOSITE_SCHUR) {
1302:     PetscOptionsGetEnum(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);
1303:     if (flg) {PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");}
1304:     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);
1305:     PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);
1306:   }
1307:   PetscOptionsTail();
1308:   return(0);
1309: }

1311: /*------------------------------------------------------------------------------------*/

1315: static PetscErrorCode  PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1316: {
1317:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1318:   PetscErrorCode    ierr;
1319:   PC_FieldSplitLink ilink,next = jac->head;
1320:   char              prefix[128];
1321:   PetscInt          i;

1324:   if (jac->splitdefined) {
1325:     PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1326:     return(0);
1327:   }
1328:   for (i=0; i<n; i++) {
1329:     if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs);
1330:     if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]);
1331:   }
1332:   PetscNew(&ilink);
1333:   if (splitname) {
1334:     PetscStrallocpy(splitname,&ilink->splitname);
1335:   } else {
1336:     PetscMalloc1(3,&ilink->splitname);
1337:     PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);
1338:   }
1339:   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 */
1340:   PetscMalloc1(n,&ilink->fields);
1341:   PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));
1342:   PetscMalloc1(n,&ilink->fields_col);
1343:   PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));

1345:   ilink->nfields = n;
1346:   ilink->next    = NULL;
1347:   KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1348:   KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1349:   PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1350:   KSPSetType(ilink->ksp,KSPPREONLY);
1351:   PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);

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

1356:   if (!next) {
1357:     jac->head       = ilink;
1358:     ilink->previous = NULL;
1359:   } else {
1360:     while (next->next) {
1361:       next = next->next;
1362:     }
1363:     next->next      = ilink;
1364:     ilink->previous = next;
1365:   }
1366:   jac->nsplits++;
1367:   return(0);
1368: }

1372: static PetscErrorCode  PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp)
1373: {
1374:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

1378:   PetscMalloc1(jac->nsplits,subksp);
1379:   MatSchurComplementGetKSP(jac->schur,*subksp);

1381:   (*subksp)[1] = jac->kspschur;
1382:   if (n) *n = jac->nsplits;
1383:   return(0);
1384: }

1388: static PetscErrorCode  PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1389: {
1390:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1391:   PetscErrorCode    ierr;
1392:   PetscInt          cnt   = 0;
1393:   PC_FieldSplitLink ilink = jac->head;

1396:   PetscMalloc1(jac->nsplits,subksp);
1397:   while (ilink) {
1398:     (*subksp)[cnt++] = ilink->ksp;
1399:     ilink            = ilink->next;
1400:   }
1401:   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);
1402:   if (n) *n = jac->nsplits;
1403:   return(0);
1404: }

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

1411:     Input Parameters:
1412: +   pc  - the preconditioner context
1413: +   is - the index set that defines the indices to which the fieldsplit is to be restricted

1415:     Level: advanced

1417: @*/
1418: PetscErrorCode  PCFieldSplitRestrictIS(PC pc,IS isy)
1419: {

1425:   PetscTryMethod(pc,"PCFieldSplitRestrictIS_C",(PC,IS),(pc,isy));
1426:   return(0);
1427: }


1432: static PetscErrorCode  PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy)
1433: {
1434:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1435:   PetscErrorCode    ierr;
1436:   PC_FieldSplitLink ilink = jac->head, next;
1437:   PetscInt          localsize,size,sizez,i;
1438:   const PetscInt    *ind, *indz;
1439:   PetscInt          *indc, *indcz;
1440:   PetscBool         flg;

1443:   ISGetLocalSize(isy,&localsize);
1444:   MPI_Scan(&localsize,&size,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)isy));
1445:   size -= localsize;
1446:   while(ilink) {
1447:     IS isrl,isr;
1448:     PC subpc;
1449:     if (jac->reset) {
1450:       ISEmbed(ilink->is_orig, isy, PETSC_TRUE, &isrl);
1451:     } else {
1452:       ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl);
1453:     }
1454:     ISGetLocalSize(isrl,&localsize);
1455:     PetscMalloc1(localsize,&indc);
1456:     ISGetIndices(isrl,&ind);
1457:     PetscMemcpy(indc,ind,localsize*sizeof(PetscInt));
1458:     ISRestoreIndices(isrl,&ind);
1459:     ISDestroy(&isrl);
1460:     for (i=0; i<localsize; i++) *(indc+i) += size;
1461:     ISCreateGeneral(PetscObjectComm((PetscObject)isy),localsize,indc,PETSC_OWN_POINTER,&isr);
1462:     PetscObjectReference((PetscObject)isr);
1463:     ISDestroy(&ilink->is);
1464:     ilink->is     = isr;
1465:     PetscObjectReference((PetscObject)isr);
1466:     ISDestroy(&ilink->is_col);
1467:     ilink->is_col = isr;
1468:     ISDestroy(&isr);
1469:     KSPGetPC(ilink->ksp, &subpc);
1470:     PetscObjectTypeCompare((PetscObject)subpc,PCFIELDSPLIT,&flg);
1471:     if(flg) {
1472:       IS iszl,isz;
1473:       MPI_Comm comm;
1474:       if (jac->reset) {
1475:         ISGetLocalSize(ilink->is_orig,&localsize);
1476:         comm = PetscObjectComm((PetscObject)ilink->is_orig);
1477:         ISEmbed(isy, ilink->is_orig, PETSC_TRUE, &iszl);
1478:       } else {
1479:         ISGetLocalSize(ilink->is,&localsize);
1480:         comm = PetscObjectComm((PetscObject)ilink->is);
1481:         ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl);
1482:       }
1483:       MPI_Scan(&localsize,&sizez,1,MPIU_INT,MPI_SUM,comm);
1484:       sizez -= localsize;
1485:       ISGetLocalSize(iszl,&localsize);
1486:       PetscMalloc1(localsize,&indcz);
1487:       ISGetIndices(iszl,&indz);
1488:       PetscMemcpy(indcz,indz,localsize*sizeof(PetscInt));
1489:       ISRestoreIndices(iszl,&indz);
1490:       ISDestroy(&iszl);
1491:       for (i=0; i<localsize; i++) *(indcz+i) += sizez;
1492:       ISCreateGeneral(comm,localsize,indcz,PETSC_OWN_POINTER,&isz);
1493:       PCFieldSplitRestrictIS(subpc,isz);
1494:       ISDestroy(&isz);
1495:     }
1496:     next = ilink->next;
1497:     ilink = next;
1498:   }
1499:   jac->isrestrict = PETSC_TRUE;
1500:   return(0);
1501: }

1505: static PetscErrorCode  PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is)
1506: {
1507:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1508:   PetscErrorCode    ierr;
1509:   PC_FieldSplitLink ilink, next = jac->head;
1510:   char              prefix[128];

1513:   if (jac->splitdefined) {
1514:     PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1515:     return(0);
1516:   }
1517:   PetscNew(&ilink);
1518:   if (splitname) {
1519:     PetscStrallocpy(splitname,&ilink->splitname);
1520:   } else {
1521:     PetscMalloc1(8,&ilink->splitname);
1522:     PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);
1523:   }
1524:   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 */
1525:   PetscObjectReference((PetscObject)is);
1526:   ISDestroy(&ilink->is);
1527:   ilink->is     = is;
1528:   PetscObjectReference((PetscObject)is);
1529:   ISDestroy(&ilink->is_col);
1530:   ilink->is_col = is;
1531:   ilink->next   = NULL;
1532:   KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1533:   KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1534:   PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1535:   KSPSetType(ilink->ksp,KSPPREONLY);
1536:   PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);

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

1541:   if (!next) {
1542:     jac->head       = ilink;
1543:     ilink->previous = NULL;
1544:   } else {
1545:     while (next->next) {
1546:       next = next->next;
1547:     }
1548:     next->next      = ilink;
1549:     ilink->previous = next;
1550:   }
1551:   jac->nsplits++;
1552:   return(0);
1553: }

1557: /*@
1558:     PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner

1560:     Logically Collective on PC

1562:     Input Parameters:
1563: +   pc  - the preconditioner context
1564: .   splitname - name of this split, if NULL the number of the split is used
1565: .   n - the number of fields in this split
1566: -   fields - the fields in this split

1568:     Level: intermediate

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

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

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

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

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

1586: @*/
1587: PetscErrorCode  PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1588: {

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

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

1605:     Logically Collective on PC

1607:     Input Parameters:
1608: +   pc  - the preconditioner object
1609: -   flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from

1611:     Options Database:
1612: .     -pc_fieldsplit_diag_use_amat

1614:     Level: intermediate

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

1618: @*/
1619: PetscErrorCode  PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg)
1620: {
1621:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1622:   PetscBool      isfs;

1627:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1628:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1629:   jac->diag_use_amat = flg;
1630:   return(0);
1631: }

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

1638:     Logically Collective on PC

1640:     Input Parameters:
1641: .   pc  - the preconditioner object

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


1647:     Level: intermediate

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

1651: @*/
1652: PetscErrorCode  PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg)
1653: {
1654:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1655:   PetscBool      isfs;

1661:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1662:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1663:   *flg = jac->diag_use_amat;
1664:   return(0);
1665: }

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

1672:     Logically Collective on PC

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

1678:     Options Database:
1679: .     -pc_fieldsplit_off_diag_use_amat

1681:     Level: intermediate

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

1685: @*/
1686: PetscErrorCode  PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg)
1687: {
1688:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1689:   PetscBool      isfs;

1694:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1695:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1696:   jac->offdiag_use_amat = flg;
1697:   return(0);
1698: }

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

1705:     Logically Collective on PC

1707:     Input Parameters:
1708: .   pc  - the preconditioner object

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


1714:     Level: intermediate

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

1718: @*/
1719: PetscErrorCode  PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg)
1720: {
1721:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1722:   PetscBool      isfs;

1728:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1729:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1730:   *flg = jac->offdiag_use_amat;
1731:   return(0);
1732: }



1738: /*@C
1739:     PCFieldSplitSetIS - Sets the exact elements for field

1741:     Logically Collective on PC

1743:     Input Parameters:
1744: +   pc  - the preconditioner context
1745: .   splitname - name of this split, if NULL the number of the split is used
1746: -   is - the index set that defines the vector elements in this field


1749:     Notes:
1750:     Use PCFieldSplitSetFields(), for fields defined by strided types.

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

1755:     Level: intermediate

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

1759: @*/
1760: PetscErrorCode  PCFieldSplitSetIS(PC pc,const char splitname[],IS is)
1761: {

1768:   PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));
1769:   return(0);
1770: }

1774: /*@
1775:     PCFieldSplitGetIS - Retrieves the elements for a field as an IS

1777:     Logically Collective on PC

1779:     Input Parameters:
1780: +   pc  - the preconditioner context
1781: -   splitname - name of this split

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

1786:     Level: intermediate

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

1790: @*/
1791: PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is)
1792: {

1799:   {
1800:     PC_FieldSplit     *jac  = (PC_FieldSplit*) pc->data;
1801:     PC_FieldSplitLink ilink = jac->head;
1802:     PetscBool         found;

1804:     *is = NULL;
1805:     while (ilink) {
1806:       PetscStrcmp(ilink->splitname, splitname, &found);
1807:       if (found) {
1808:         *is = ilink->is;
1809:         break;
1810:       }
1811:       ilink = ilink->next;
1812:     }
1813:   }
1814:   return(0);
1815: }

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

1823:     Logically Collective on PC

1825:     Input Parameters:
1826: +   pc  - the preconditioner context
1827: -   bs - the block size

1829:     Level: intermediate

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

1833: @*/
1834: PetscErrorCode  PCFieldSplitSetBlockSize(PC pc,PetscInt bs)
1835: {

1841:   PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));
1842:   return(0);
1843: }

1847: /*@C
1848:    PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits

1850:    Collective on KSP

1852:    Input Parameter:
1853: .  pc - the preconditioner context

1855:    Output Parameters:
1856: +  n - the number of splits
1857: -  pc - the array of KSP contexts

1859:    Note:
1860:    After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user
1861:    (not the KSP just the array that contains them).

1863:    You must call KSPSetUp() before calling PCFieldSplitGetSubKSP().

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


1870:    Level: advanced

1872: .seealso: PCFIELDSPLIT
1873: @*/
1874: PetscErrorCode  PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
1875: {

1881:   PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
1882:   return(0);
1883: }

1887: /*@
1888:     PCFieldSplitSetSchurPre -  Indicates what operator is used to construct the preconditioner for the Schur complement.
1889:       A11 matrix. Otherwise no preconditioner is used.

1891:     Collective on PC

1893:     Input Parameters:
1894: +   pc      - the preconditioner context
1895: .   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 
1896:               PC_FIELDSPLIT_SCHUR_PRE_SELFP, and PC_FIELDSPLIT_SCHUR_PRE_FULL
1897: -   userpre - matrix to use for preconditioning, or NULL

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

1902:     Notes:
1903: $    If ptype is
1904: $        a11 then the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner
1905: $             matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix
1906: $        self the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix:
1907: $             The only preconditioner that currently works with this symbolic respresentation matrix object is the PCLSC
1908: $             preconditioner
1909: $        user then the preconditioner for the Schur complement is generated from the user provided matrix (pre argument
1910: $             to this function).
1911: $        selfp then the preconditioning for the Schur complement is generated from an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01
1912: $             This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be
1913: $             lumped before extracting the diagonal using the additional option -fieldsplit_1_mat_schur_complement_ainv_type lump
1914: $        full then the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation computed internally by PFIELDSPLIT (this is expensive)
1915: $             useful mostly as a test that the Schur complement approach can work for your problem

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

1921:     Level: intermediate

1923: .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType,
1924:           MatSchurComplementSetAinvType(), PCLSC

1926: @*/
1927: PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
1928: {

1933:   PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));
1934:   return(0);
1935: }
1936: PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */

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

1944:     Logically Collective on PC

1946:     Input Parameters:
1947: .   pc      - the preconditioner context

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

1953:     Level: intermediate

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

1957: @*/
1958: PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
1959: {

1964:   PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));
1965:   return(0);
1966: }

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

1973:     Not collective

1975:     Input Parameter:
1976: .   pc      - the preconditioner context

1978:     Output Parameter:
1979: .   S       - the Schur complement matrix

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

1984:     Level: advanced

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

1988: @*/
1989: PetscErrorCode  PCFieldSplitSchurGetS(PC pc,Mat *S)
1990: {
1992:   const char*    t;
1993:   PetscBool      isfs;
1994:   PC_FieldSplit  *jac;

1998:   PetscObjectGetType((PetscObject)pc,&t);
1999:   PetscStrcmp(t,PCFIELDSPLIT,&isfs);
2000:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
2001:   jac = (PC_FieldSplit*)pc->data;
2002:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
2003:   if (S) *S = jac->schur;
2004:   return(0);
2005: }

2009: /*@
2010:     PCFieldSplitSchurRestoreS -  restores the MatSchurComplement object used by this PC

2012:     Not collective

2014:     Input Parameters:
2015: +   pc      - the preconditioner context
2016: .   S       - the Schur complement matrix

2018:     Level: advanced

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

2022: @*/
2023: PetscErrorCode  PCFieldSplitSchurRestoreS(PC pc,Mat *S)
2024: {
2026:   const char*    t;
2027:   PetscBool      isfs;
2028:   PC_FieldSplit  *jac;

2032:   PetscObjectGetType((PetscObject)pc,&t);
2033:   PetscStrcmp(t,PCFIELDSPLIT,&isfs);
2034:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
2035:   jac = (PC_FieldSplit*)pc->data;
2036:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
2037:   if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten");
2038:   return(0);
2039: }


2044: static PetscErrorCode  PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
2045: {
2046:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

2050:   jac->schurpre = ptype;
2051:   if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
2052:     MatDestroy(&jac->schur_user);
2053:     jac->schur_user = pre;
2054:     PetscObjectReference((PetscObject)jac->schur_user);
2055:   }
2056:   return(0);
2057: }

2061: static PetscErrorCode  PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
2062: {
2063:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

2066:   *ptype = jac->schurpre;
2067:   *pre   = jac->schur_user;
2068:   return(0);
2069: }

2073: /*@
2074:     PCFieldSplitSetSchurFactType -  sets which blocks of the approximate block factorization to retain

2076:     Collective on PC

2078:     Input Parameters:
2079: +   pc  - the preconditioner context
2080: -   ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default

2082:     Options Database:
2083: .     -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full


2086:     Level: intermediate

2088:     Notes:
2089:     The FULL factorization is

2091: $   (A   B)  = (1       0) (A   0) (1  Ainv*B)
2092: $   (C   D)    (C*Ainv  1) (0   S) (0     1  )

2094:     where S = D - 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,
2095:     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).

2097:     If applied exactly, FULL factorization is a direct solver. The preconditioned operator with LOWER or UPPER has all eigenvalues equal to 1 and minimal polynomial
2098:     of degree 2, so KSPGMRES converges in 2 iterations. If the iteration count is very low, consider using KSPFGMRES or KSPGCR which can use one less preconditioner
2099:     application in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice. With DIAG,
2100:     the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations.

2102:     For symmetric problems in which A is positive definite and S is negative definite, DIAG can be used with KSPMINRES. Note that a flexible method like KSPFGMRES
2103:     or KSPGCR must be used if the fieldsplit preconditioner is nonlinear (e.g. a few iterations of a Krylov method is used inside a split).

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

2109: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType
2110: @*/
2111: PetscErrorCode  PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype)
2112: {

2117:   PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));
2118:   return(0);
2119: }

2123: static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype)
2124: {
2125:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2128:   jac->schurfactorization = ftype;
2129:   return(0);
2130: }

2134: /*@C
2135:    PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement

2137:    Collective on KSP

2139:    Input Parameter:
2140: .  pc - the preconditioner context

2142:    Output Parameters:
2143: +  A00 - the (0,0) block
2144: .  A01 - the (0,1) block
2145: .  A10 - the (1,0) block
2146: -  A11 - the (1,1) block

2148:    Level: advanced

2150: .seealso: PCFIELDSPLIT
2151: @*/
2152: PetscErrorCode  PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11)
2153: {
2154:   PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;

2158:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach.");
2159:   if (A00) *A00 = jac->pmat[0];
2160:   if (A01) *A01 = jac->B;
2161:   if (A10) *A10 = jac->C;
2162:   if (A11) *A11 = jac->pmat[1];
2163:   return(0);
2164: }

2168: static PetscErrorCode  PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type)
2169: {
2170:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

2174:   jac->type = type;
2175:   if (type == PC_COMPOSITE_SCHUR) {
2176:     pc->ops->apply = PCApply_FieldSplit_Schur;
2177:     pc->ops->view  = PCView_FieldSplit_Schur;

2179:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);
2180:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);
2181:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);
2182:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);

2184:   } else {
2185:     pc->ops->apply = PCApply_FieldSplit;
2186:     pc->ops->view  = PCView_FieldSplit;

2188:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2189:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);
2190:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);
2191:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);
2192:   }
2193:   return(0);
2194: }

2198: static PetscErrorCode  PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs)
2199: {
2200:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2203:   if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs);
2204:   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);
2205:   jac->bs = bs;
2206:   return(0);
2207: }

2211: /*@
2212:    PCFieldSplitSetType - Sets the type of fieldsplit preconditioner.

2214:    Collective on PC

2216:    Input Parameter:
2217: .  pc - the preconditioner context
2218: .  type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR

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

2223:    Level: Intermediate

2225: .keywords: PC, set, type, composite preconditioner, additive, multiplicative

2227: .seealso: PCCompositeSetType()

2229: @*/
2230: PetscErrorCode  PCFieldSplitSetType(PC pc,PCCompositeType type)
2231: {

2236:   PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));
2237:   return(0);
2238: }

2242: /*@
2243:   PCFieldSplitGetType - Gets the type of fieldsplit preconditioner.

2245:   Not collective

2247:   Input Parameter:
2248: . pc - the preconditioner context

2250:   Output Parameter:
2251: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR

2253:   Level: Intermediate

2255: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
2256: .seealso: PCCompositeSetType()
2257: @*/
2258: PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
2259: {
2260:   PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;

2265:   *type = jac->type;
2266:   return(0);
2267: }

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

2274:    Logically Collective

2276:    Input Parameters:
2277: +  pc   - the preconditioner context
2278: -  flg  - boolean indicating whether to use field splits defined by the DM

2280:    Options Database Key:
2281: .  -pc_fieldsplit_dm_splits

2283:    Level: Intermediate

2285: .keywords: PC, DM, composite preconditioner, additive, multiplicative

2287: .seealso: PCFieldSplitGetDMSplits()

2289: @*/
2290: PetscErrorCode  PCFieldSplitSetDMSplits(PC pc,PetscBool flg)
2291: {
2292:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
2293:   PetscBool      isfs;

2299:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2300:   if (isfs) {
2301:     jac->dm_splits = flg;
2302:   }
2303:   return(0);
2304: }


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

2312:    Logically Collective

2314:    Input Parameter:
2315: .  pc   - the preconditioner context

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

2320:    Level: Intermediate

2322: .keywords: PC, DM, composite preconditioner, additive, multiplicative

2324: .seealso: PCFieldSplitSetDMSplits()

2326: @*/
2327: PetscErrorCode  PCFieldSplitGetDMSplits(PC pc,PetscBool* flg)
2328: {
2329:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
2330:   PetscBool      isfs;

2336:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2337:   if (isfs) {
2338:     if(flg) *flg = jac->dm_splits;
2339:   }
2340:   return(0);
2341: }

2343: /* -------------------------------------------------------------------------------------*/
2344: /*MC
2345:    PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
2346:                   fields or groups of fields. See the users manual section "Solving Block Matrices" for more details.

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

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

2354:    Level: intermediate

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

2365: -    Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_
2366:      for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields

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

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

2377: $     For the Schur complement preconditioner if J = ( A00 A01 )
2378: $                                                    ( A10 A11 )
2379: $     the preconditioner using full factorization is
2380: $              ( I   -ksp(A00) A01 ) ( inv(A00)     0  ) (     I          0  )
2381: $              ( 0         I       ) (   0      ksp(S) ) ( -A10 ksp(A00)  I  )
2382:      where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_.  S is the Schur complement
2383: $              S = A11 - A10 ksp(A00) A01
2384:      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
2385:      in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0,
2386:      it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default
2387:      A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() for all the possible ways to construct the preconditioner for S.

2389:      The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above,
2390:      diag gives
2391: $              ( inv(A00)     0   )
2392: $              (   0      -ksp(S) )
2393:      note that slightly counter intuitively there is a negative in front of the ksp(S) so that the preconditioner is positive definite. The lower factorization is the inverse of
2394: $              (  A00   0 )
2395: $              (  A10   S )
2396:      where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of
2397: $              ( A00 A01 )
2398: $              (  0   S  )
2399:      where again the inverses of A00 and S are applied using KSPs.

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

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

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

2411:    Concepts: physics based preconditioners, block preconditioners

2413:    There is a nice discussion of block preconditioners in

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

2419:    The Constrained Pressure Preconditioner (CPR) does not appear to be currently implementable directly with PCFIELDSPLIT. CPR solves first the Schur complemented pressure equation, updates the
2420:    residual on all variables and then applies a simple ILU like preconditioner on all the variables. So it is very much like the full Schur complement with selfp representing the Schur complement but instead
2421:    of backsolving for the saturations in the last step it solves a full coupled (ILU) system for updates to all the variables.

2423: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC,
2424:            PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(),
2425:            MatSchurComplementSetAinvType()
2426: M*/

2430: PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
2431: {
2433:   PC_FieldSplit  *jac;

2436:   PetscNewLog(pc,&jac);

2438:   jac->bs                 = -1;
2439:   jac->nsplits            = 0;
2440:   jac->type               = PC_COMPOSITE_MULTIPLICATIVE;
2441:   jac->schurpre           = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
2442:   jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
2443:   jac->dm_splits          = PETSC_TRUE;

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

2447:   pc->ops->apply           = PCApply_FieldSplit;
2448:   pc->ops->applytranspose  = PCApplyTranspose_FieldSplit;
2449:   pc->ops->setup           = PCSetUp_FieldSplit;
2450:   pc->ops->reset           = PCReset_FieldSplit;
2451:   pc->ops->destroy         = PCDestroy_FieldSplit;
2452:   pc->ops->setfromoptions  = PCSetFromOptions_FieldSplit;
2453:   pc->ops->view            = PCView_FieldSplit;
2454:   pc->ops->applyrichardson = 0;

2456:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2457:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);
2458:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);
2459:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);
2460:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);
2461:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",PCFieldSplitRestrictIS_FieldSplit);
2462:   return(0);
2463: }