Actual source code: fieldsplit.c

petsc-master 2017-10-20
Report Typos and Errors


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

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

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

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

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

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

 49:   PC_FieldSplitLink         head;
 50:   PetscBool                 isrestrict;             /* indicates PCFieldSplitRestrictIS() has been last called on this object, hack */
 51:   PetscBool                 suboptionsset;          /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */
 52:   PetscBool                 dm_splits;              /* Whether to use DMCreateFieldDecomposition() whenever possible */
 53:   PetscBool                 diag_use_amat;          /* Whether to extract diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
 54:   PetscBool                 offdiag_use_amat;       /* Whether to extract off-diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
 55: } PC_FieldSplit;

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

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


 79:  #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: }

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

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

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

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

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

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

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

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

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

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

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

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

406:         MatGetOwnershipRange(pc->mat,&nmin,&nmax);
407:         MatFindZeroDiagonals(pc->mat,&zerodiags);
408:         ISComplement(zerodiags,nmin,nmax,&rest);
409:         PCFieldSplitSetIS(pc,"0",rest);
410:         PCFieldSplitSetIS(pc,"1",zerodiags);
411:         ISDestroy(&zerodiags);
412:         ISDestroy(&rest);
413:       } else if (coupling) {
414:         IS       coupling,rest;
415:         PetscInt nmin,nmax;

417:         MatGetOwnershipRange(pc->mat,&nmin,&nmax);
418:         MatFindOffBlockDiagonalEntries(pc->mat,&coupling);
419:         ISCreateStride(PetscObjectComm((PetscObject)pc->mat),nmax-nmin,nmin,1,&rest);
420:         ISSetIdentity(rest);
421:         PCFieldSplitSetIS(pc,"0",rest);
422:         PCFieldSplitSetIS(pc,"1",coupling);
423:         ISDestroy(&coupling);
424:         ISDestroy(&rest);
425:       } else {
426:         PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&fieldsplit_default,NULL);
427:         if (!fieldsplit_default) {
428:           /* Allow user to set fields from command line,  if bs was known at the time of PCSetFromOptions_FieldSplit()
429:            then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */
430:           PCFieldSplitSetRuntimeSplits_Private(pc);
431:           if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
432:         }
433:         if ((fieldsplit_default || !jac->splitdefined) && !jac->isrestrict) {
434:           PetscInfo(pc,"Using default splitting of fields\n");
435:           for (i=0; i<jac->bs; i++) {
436:             char splitname[8];
437:             PetscSNPrintf(splitname,sizeof(splitname),"%D",i);
438:             PCFieldSplitSetFields(pc,splitname,1,&i,&i);
439:           }
440:           jac->defaultsplit = PETSC_TRUE;
441:         }
442:       }
443:     }
444:   } else if (jac->nsplits == 1) {
445:     if (ilink->is) {
446:       IS       is2;
447:       PetscInt nmin,nmax;

449:       MatGetOwnershipRange(pc->mat,&nmin,&nmax);
450:       ISComplement(ilink->is,nmin,nmax,&is2);
451:       PCFieldSplitSetIS(pc,"1",is2);
452:       ISDestroy(&is2);
453:     } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()");
454:   }

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

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

462: static PetscErrorCode PCSetUp_FieldSplit(PC pc)
463: {
464:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
465:   PetscErrorCode    ierr;
466:   PC_FieldSplitLink ilink;
467:   PetscInt          i,nsplit;
468:   PetscBool         sorted, sorted_col;

471:   PCFieldSplitSetDefaults(pc);
472:   nsplit = jac->nsplits;
473:   ilink  = jac->head;

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

479:     jac->issetup = PETSC_TRUE;

481:     /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */
482:     if (jac->defaultsplit || !ilink->is) {
483:       if (jac->bs <= 0) jac->bs = nsplit;
484:     }
485:     bs     = jac->bs;
486:     MatGetOwnershipRange(pc->pmat,&rstart,&rend);
487:     nslots = (rend - rstart)/bs;
488:     for (i=0; i<nsplit; i++) {
489:       if (jac->defaultsplit) {
490:         ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+i,nsplit,&ilink->is);
491:         ISDuplicate(ilink->is,&ilink->is_col);
492:       } else if (!ilink->is) {
493:         if (ilink->nfields > 1) {
494:           PetscInt *ii,*jj,j,k,nfields = ilink->nfields,*fields = ilink->fields,*fields_col = ilink->fields_col;
495:           PetscMalloc1(ilink->nfields*nslots,&ii);
496:           PetscMalloc1(ilink->nfields*nslots,&jj);
497:           for (j=0; j<nslots; j++) {
498:             for (k=0; k<nfields; k++) {
499:               ii[nfields*j + k] = rstart + bs*j + fields[k];
500:               jj[nfields*j + k] = rstart + bs*j + fields_col[k];
501:             }
502:           }
503:           ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);
504:           ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,jj,PETSC_OWN_POINTER,&ilink->is_col);
505:           ISSetBlockSize(ilink->is, nfields);
506:           ISSetBlockSize(ilink->is_col, nfields);
507:         } else {
508:           ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields[0],bs,&ilink->is);
509:           ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields_col[0],bs,&ilink->is_col);
510:        }
511:       }
512:       ISSorted(ilink->is,&sorted);
513:       if (ilink->is_col) { ISSorted(ilink->is_col,&sorted_col); }
514:       if (!sorted || !sorted_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split");
515:       ilink = ilink->next;
516:     }
517:   }

519:   ilink = jac->head;
520:   if (!jac->pmat) {
521:     Vec xtmp;

523:     MatCreateVecs(pc->pmat,&xtmp,NULL);
524:     PetscMalloc1(nsplit,&jac->pmat);
525:     PetscMalloc2(nsplit,&jac->x,nsplit,&jac->y);
526:     for (i=0; i<nsplit; i++) {
527:       MatNullSpace sp;

529:       /* Check for preconditioning matrix attached to IS */
530:       PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &jac->pmat[i]);
531:       if (jac->pmat[i]) {
532:         PetscObjectReference((PetscObject) jac->pmat[i]);
533:         if (jac->type == PC_COMPOSITE_SCHUR) {
534:           jac->schur_user = jac->pmat[i];

536:           PetscObjectReference((PetscObject) jac->schur_user);
537:         }
538:       } else {
539:         const char *prefix;
540:         MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->pmat[i]);
541:         KSPGetOptionsPrefix(ilink->ksp,&prefix);
542:         MatSetOptionsPrefix(jac->pmat[i],prefix);
543:         MatViewFromOptions(jac->pmat[i],NULL,"-mat_view");
544:       }
545:       /* create work vectors for each split */
546:       MatCreateVecs(jac->pmat[i],&jac->x[i],&jac->y[i]);
547:       ilink->x = jac->x[i]; ilink->y = jac->y[i]; ilink->z = NULL;
548:       /* compute scatter contexts needed by multiplicative versions and non-default splits */
549:       VecScatterCreate(xtmp,ilink->is,jac->x[i],NULL,&ilink->sctx);
550:       PetscObjectQuery((PetscObject) ilink->is, "nearnullspace", (PetscObject*) &sp);
551:       if (sp) {
552:         MatSetNearNullSpace(jac->pmat[i], sp);
553:       }
554:       ilink = ilink->next;
555:     }
556:     VecDestroy(&xtmp);
557:   } else {
558:     for (i=0; i<nsplit; i++) {
559:       Mat pmat;

561:       /* Check for preconditioning matrix attached to IS */
562:       PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &pmat);
563:       if (!pmat) {
564:         MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->pmat[i]);
565:       }
566:       ilink = ilink->next;
567:     }
568:   }
569:   if (jac->diag_use_amat) {
570:     ilink = jac->head;
571:     if (!jac->mat) {
572:       PetscMalloc1(nsplit,&jac->mat);
573:       for (i=0; i<nsplit; i++) {
574:         MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i]);
575:         ilink = ilink->next;
576:       }
577:     } else {
578:       for (i=0; i<nsplit; i++) {
579:         if (jac->mat[i]) {MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->mat[i]);}
580:         ilink = ilink->next;
581:       }
582:     }
583:   } else {
584:     jac->mat = jac->pmat;
585:   }

587:   /* Check for null space attached to IS */
588:   ilink = jac->head;
589:   for (i=0; i<nsplit; i++) {
590:     MatNullSpace sp;

592:     PetscObjectQuery((PetscObject) ilink->is, "nullspace", (PetscObject*) &sp);
593:     if (sp) {
594:       MatSetNullSpace(jac->mat[i], sp);
595:     }
596:     ilink = ilink->next;
597:   }

599:   if (jac->type != PC_COMPOSITE_ADDITIVE  && jac->type != PC_COMPOSITE_SCHUR) {
600:     /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */
601:     /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */
602:     ilink = jac->head;
603:     if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) {
604:       /* 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 */
605:       if (!jac->Afield) {
606:         PetscCalloc1(nsplit,&jac->Afield);
607:         if (jac->offdiag_use_amat) {
608:           MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);
609:         } else {
610:           MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);
611:         }
612:       } else {
613:         if (jac->offdiag_use_amat) {
614:           MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_REUSE_MATRIX,&jac->Afield[1]);
615:         } else {
616:           MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,MAT_REUSE_MATRIX,&jac->Afield[1]);
617:         }
618:       }
619:     } else {
620:       if (!jac->Afield) {
621:         PetscMalloc1(nsplit,&jac->Afield);
622:         for (i=0; i<nsplit; i++) {
623:           if (jac->offdiag_use_amat) {
624:             MatCreateSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
625:           } else {
626:             MatCreateSubMatrix(pc->pmat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
627:           }
628:           ilink = ilink->next;
629:         }
630:       } else {
631:         for (i=0; i<nsplit; i++) {
632:           if (jac->offdiag_use_amat) {
633:             MatCreateSubMatrix(pc->mat,ilink->is,NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);
634:           } else {
635:             MatCreateSubMatrix(pc->pmat,ilink->is,NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);
636:           }
637:           ilink = ilink->next;
638:         }
639:       }
640:     }
641:   }

643:   if (jac->type == PC_COMPOSITE_SCHUR) {
644:     IS          ccis;
645:     PetscBool   isspd;
646:     PetscInt    rstart,rend;
647:     char        lscname[256];
648:     PetscObject LSC_L;

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

652:     /* If pc->mat is SPD, don't scale by -1 the Schur complement */
653:     if (jac->schurscale == (PetscScalar)-1.0) {
654:       MatGetOption(pc->pmat,MAT_SPD,&isspd);
655:       jac->schurscale = (isspd == PETSC_TRUE) ? 1.0 : -1.0;
656:     }

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

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

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

701:       /* extract the A01 and A10 matrices */
702:       ilink = jac->head;
703:       ISComplement(ilink->is_col,rstart,rend,&ccis);
704:       if (jac->offdiag_use_amat) {
705:         MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
706:       } else {
707:         MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
708:       }
709:       ISDestroy(&ccis);
710:       ilink = ilink->next;
711:       ISComplement(ilink->is_col,rstart,rend,&ccis);
712:       if (jac->offdiag_use_amat) {
713:         MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
714:       } else {
715:         MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
716:       }
717:       ISDestroy(&ccis);

719:       /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
720:       MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);
721:       MatSetType(jac->schur,MATSCHURCOMPLEMENT);
722:       MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
723:       PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
724:       /* 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? */
725:       MatSetOptionsPrefix(jac->schur,schurmatprefix);
726:       MatSetFromOptions(jac->schur);
727:       MatGetNullSpace(jac->mat[1], &sp);
728:       if (sp) {
729:         MatSetNullSpace(jac->schur, sp);
730:       }

732:       PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);
733:       PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
734:       if (flg) {
735:         DM  dmInner;
736:         KSP kspInner;

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

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

762:       PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);
763:       PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
764:       if (flg) {
765:         DM dmInner;

767:         PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
768:         KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);
769:         KSPSetErrorIfNotConverged(jac->kspupper,pc->erroriffailure);
770:         KSPSetOptionsPrefix(jac->kspupper, schurprefix);
771:         KSPGetDM(jac->head->ksp, &dmInner);
772:         KSPSetDM(jac->kspupper, dmInner);
773:         KSPSetDMActive(jac->kspupper, PETSC_FALSE);
774:         KSPSetFromOptions(jac->kspupper);
775:         KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0]);
776:         VecDuplicate(jac->head->x, &jac->head->z);
777:       } else {
778:         jac->kspupper = jac->head->ksp;
779:         PetscObjectReference((PetscObject) jac->head->ksp);
780:       }

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

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

836:   jac->suboptionsset = PETSC_TRUE;
837:   return(0);
838: }

840: #define FieldSplitSplitSolveAdd(ilink,xx,yy) \
841:   (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
842:    VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
843:    PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\
844:    KSPSolve(ilink->ksp,ilink->x,ilink->y) ||                               \
845:    PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\
846:    VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) ||  \
847:    VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE))

849: static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y)
850: {
851:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
852:   PetscErrorCode    ierr;
853:   PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
854:   KSP               kspA   = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;

857:   switch (jac->schurfactorization) {
858:   case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
859:     /* [A00 0; 0 -S], positive definite, suitable for MINRES */
860:     VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
861:     VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
862:     VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
863:     PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
864:     KSPSolve(kspA,ilinkA->x,ilinkA->y);
865:     PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
866:     VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
867:     VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
868:     PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
869:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
870:     PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
871:     VecScale(ilinkD->y,jac->schurscale);
872:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
873:     VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
874:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
875:     break;
876:   case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
877:     /* [A00 0; A10 S], suitable for left preconditioning */
878:     VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
879:     VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
880:     PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
881:     KSPSolve(kspA,ilinkA->x,ilinkA->y);
882:     PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
883:     MatMult(jac->C,ilinkA->y,ilinkD->x);
884:     VecScale(ilinkD->x,-1.);
885:     VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
886:     VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
887:     VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
888:     PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
889:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
890:     PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
891:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
892:     VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
893:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
894:     break;
895:   case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
896:     /* [A00 A01; 0 S], suitable for right preconditioning */
897:     VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
898:     VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_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);    MatMult(jac->B,ilinkD->y,ilinkA->x);
902:     VecScale(ilinkA->x,-1.);
903:     VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
904:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
905:     VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
906:     PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
907:     KSPSolve(kspA,ilinkA->x,ilinkA->y);
908:     PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
909:     VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
910:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
911:     VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
912:     break;
913:   case PC_FIELDSPLIT_SCHUR_FACT_FULL:
914:     /* [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 */
915:     VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
916:     VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
917:     PetscLogEventBegin(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);
918:     KSPSolve(kspLower,ilinkA->x,ilinkA->y);
919:     PetscLogEventEnd(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);
920:     MatMult(jac->C,ilinkA->y,ilinkD->x);
921:     VecScale(ilinkD->x,-1.0);
922:     VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
923:     VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);

925:     PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
926:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
927:     PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
928:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
929:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);

931:     if (kspUpper == kspA) {
932:       MatMult(jac->B,ilinkD->y,ilinkA->y);
933:       VecAXPY(ilinkA->x,-1.0,ilinkA->y);
934:       PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
935:       KSPSolve(kspA,ilinkA->x,ilinkA->y);
936:       PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
937:     } else {
938:       PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
939:       KSPSolve(kspA,ilinkA->x,ilinkA->y);
940:       PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
941:       MatMult(jac->B,ilinkD->y,ilinkA->x);
942:       PetscLogEventBegin(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);
943:       KSPSolve(kspUpper,ilinkA->x,ilinkA->z);
944:       PetscLogEventEnd(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);
945:       VecAXPY(ilinkA->y,-1.0,ilinkA->z);
946:     }
947:     VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
948:     VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
949:   }
950:   return(0);
951: }

953: static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y)
954: {
955:   PC_FieldSplit      *jac = (PC_FieldSplit*)pc->data;
956:   PetscErrorCode     ierr;
957:   PC_FieldSplitLink  ilink = jac->head;
958:   PetscInt           cnt,bs;
959:   KSPConvergedReason reason;

962:   if (jac->type == PC_COMPOSITE_ADDITIVE) {
963:     if (jac->defaultsplit) {
964:       VecGetBlockSize(x,&bs);
965:       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);
966:       VecGetBlockSize(y,&bs);
967:       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);
968:       VecStrideGatherAll(x,jac->x,INSERT_VALUES);
969:       while (ilink) {
970:         PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
971:         KSPSolve(ilink->ksp,ilink->x,ilink->y);
972:         PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
973:         KSPGetConvergedReason(ilink->ksp,&reason);
974:         if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
975:           pc->failedreason = PC_SUBPC_ERROR;
976:         }
977:         ilink = ilink->next;
978:       }
979:       VecStrideScatterAll(jac->y,y,INSERT_VALUES);
980:     } else {
981:       VecSet(y,0.0);
982:       while (ilink) {
983:         FieldSplitSplitSolveAdd(ilink,x,y);
984:         KSPGetConvergedReason(ilink->ksp,&reason);
985:         if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
986:           pc->failedreason = PC_SUBPC_ERROR;
987:         }
988:         ilink = ilink->next;
989:       }
990:     }
991:   } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) {
992:     VecSet(y,0.0);
993:     /* solve on first block for first block variables */
994:     VecScatterBegin(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
995:     VecScatterEnd(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
996:     PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
997:     KSPSolve(ilink->ksp,ilink->x,ilink->y);
998:     PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
999:     KSPGetConvergedReason(ilink->ksp,&reason);
1000:     if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1001:       pc->failedreason = PC_SUBPC_ERROR;
1002:     }
1003:     VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1004:     VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);

1006:     /* compute the residual only onto second block variables using first block variables */
1007:     MatMult(jac->Afield[1],ilink->y,ilink->next->x);
1008:     ilink = ilink->next;
1009:     VecScale(ilink->x,-1.0);
1010:     VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1011:     VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);

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

1076: #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \
1077:   (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1078:    VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1079:    PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \
1080:    KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) ||                  \
1081:    PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \
1082:    VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \
1083:    VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE))

1085: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y)
1086: {
1087:   PC_FieldSplit      *jac = (PC_FieldSplit*)pc->data;
1088:   PetscErrorCode     ierr;
1089:   PC_FieldSplitLink  ilink = jac->head;
1090:   PetscInt           bs;
1091:   KSPConvergedReason reason;

1094:   if (jac->type == PC_COMPOSITE_ADDITIVE) {
1095:     if (jac->defaultsplit) {
1096:       VecGetBlockSize(x,&bs);
1097:       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);
1098:       VecGetBlockSize(y,&bs);
1099:       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);
1100:       VecStrideGatherAll(x,jac->x,INSERT_VALUES);
1101:       while (ilink) {
1102:         PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1103:         KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);
1104:         PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1105:         KSPGetConvergedReason(ilink->ksp,&reason);
1106:         if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1107:           pc->failedreason = PC_SUBPC_ERROR;
1108:         }
1109:         ilink = ilink->next;
1110:       }
1111:       VecStrideScatterAll(jac->y,y,INSERT_VALUES);
1112:     } else {
1113:       VecSet(y,0.0);
1114:       while (ilink) {
1115:         FieldSplitSplitSolveAddTranspose(ilink,x,y);
1116:         KSPGetConvergedReason(ilink->ksp,&reason);
1117:         if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1118:           pc->failedreason = PC_SUBPC_ERROR;
1119:         }
1120:         ilink = ilink->next;
1121:       }
1122:     }
1123:   } else {
1124:     if (!jac->w1) {
1125:       VecDuplicate(x,&jac->w1);
1126:       VecDuplicate(x,&jac->w2);
1127:     }
1128:     VecSet(y,0.0);
1129:     if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
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:       while (ilink->next) {
1136:         ilink = ilink->next;
1137:         MatMultTranspose(pc->mat,y,jac->w1);
1138:         VecWAXPY(jac->w2,-1.0,jac->w1,x);
1139:         FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1140:       }
1141:       while (ilink->previous) {
1142:         ilink = ilink->previous;
1143:         MatMultTranspose(pc->mat,y,jac->w1);
1144:         VecWAXPY(jac->w2,-1.0,jac->w1,x);
1145:         FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1146:       }
1147:     } else {
1148:       while (ilink->next) {   /* get to last entry in linked list */
1149:         ilink = ilink->next;
1150:       }
1151:       FieldSplitSplitSolveAddTranspose(ilink,x,y);
1152:       KSPGetConvergedReason(ilink->ksp,&reason);
1153:       if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1154:         pc->failedreason = PC_SUBPC_ERROR;
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:     }
1163:   }
1164:   return(0);
1165: }

1167: static PetscErrorCode PCReset_FieldSplit(PC pc)
1168: {
1169:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1170:   PetscErrorCode    ierr;
1171:   PC_FieldSplitLink ilink = jac->head,next;

1174:   while (ilink) {
1175:     KSPDestroy(&ilink->ksp);
1176:     VecDestroy(&ilink->x);
1177:     VecDestroy(&ilink->y);
1178:     VecDestroy(&ilink->z);
1179:     VecScatterDestroy(&ilink->sctx);
1180:     ISDestroy(&ilink->is);
1181:     ISDestroy(&ilink->is_col);
1182:     PetscFree(ilink->splitname);
1183:     PetscFree(ilink->fields);
1184:     PetscFree(ilink->fields_col);
1185:     next  = ilink->next;
1186:     PetscFree(ilink);
1187:     ilink = next;
1188:   }
1189:   jac->head = NULL;
1190:   PetscFree2(jac->x,jac->y);
1191:   if (jac->mat && jac->mat != jac->pmat) {
1192:     MatDestroyMatrices(jac->nsplits,&jac->mat);
1193:   } else if (jac->mat) {
1194:     jac->mat = NULL;
1195:   }
1196:   if (jac->pmat) {MatDestroyMatrices(jac->nsplits,&jac->pmat);}
1197:   if (jac->Afield) {MatDestroyMatrices(jac->nsplits,&jac->Afield);}
1198:   jac->nsplits = 0;
1199:   VecDestroy(&jac->w1);
1200:   VecDestroy(&jac->w2);
1201:   MatDestroy(&jac->schur);
1202:   MatDestroy(&jac->schurp);
1203:   MatDestroy(&jac->schur_user);
1204:   KSPDestroy(&jac->kspschur);
1205:   KSPDestroy(&jac->kspupper);
1206:   MatDestroy(&jac->B);
1207:   MatDestroy(&jac->C);
1208:   jac->isrestrict = PETSC_FALSE;
1209:   return(0);
1210: }

1212: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
1213: {
1214:   PetscErrorCode    ierr;

1217:   PCReset_FieldSplit(pc);
1218:   PetscFree(pc->data);
1219:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);
1220:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);
1221:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);
1222:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);
1223:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);
1224:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);
1225:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);
1226:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);
1227:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",NULL);
1228:   return(0);
1229: }

1231: static PetscErrorCode PCSetFromOptions_FieldSplit(PetscOptionItems *PetscOptionsObject,PC pc)
1232: {
1233:   PetscErrorCode  ierr;
1234:   PetscInt        bs;
1235:   PetscBool       flg,stokes = PETSC_FALSE;
1236:   PC_FieldSplit   *jac = (PC_FieldSplit*)pc->data;
1237:   PCCompositeType ctype;

1240:   PetscOptionsHead(PetscOptionsObject,"FieldSplit options");
1241:   PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);
1242:   PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);
1243:   if (flg) {
1244:     PCFieldSplitSetBlockSize(pc,bs);
1245:   }
1246:   jac->diag_use_amat = pc->useAmat;
1247:   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);
1248:   jac->offdiag_use_amat = pc->useAmat;
1249:   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);
1250:   /* FIXME: No programmatic equivalent to the following. */
1251:   PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);
1252:   if (stokes) {
1253:     PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);
1254:     jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF;
1255:   }

1257:   PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);
1258:   if (flg) {
1259:     PCFieldSplitSetType(pc,ctype);
1260:   }
1261:   /* Only setup fields once */
1262:   if ((jac->bs > 0) && (jac->nsplits == 0)) {
1263:     /* only allow user to set fields from command line if bs is already known.
1264:        otherwise user can set them in PCFieldSplitSetDefaults() */
1265:     PCFieldSplitSetRuntimeSplits_Private(pc);
1266:     if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
1267:   }
1268:   if (jac->type == PC_COMPOSITE_SCHUR) {
1269:     PetscOptionsGetEnum(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);
1270:     if (flg) {PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");}
1271:     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);
1272:     PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);
1273:     PetscOptionsScalar("-pc_fieldsplit_schur_scale","Scale Schur complement","PCFieldSplitSetSchurScale",jac->schurscale,&jac->schurscale,NULL);
1274:   }
1275:   PetscOptionsTail();
1276:   return(0);
1277: }

1279: /*------------------------------------------------------------------------------------*/

1281: static PetscErrorCode  PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1282: {
1283:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1284:   PetscErrorCode    ierr;
1285:   PC_FieldSplitLink ilink,next = jac->head;
1286:   char              prefix[128];
1287:   PetscInt          i;

1290:   if (jac->splitdefined) {
1291:     PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1292:     return(0);
1293:   }
1294:   for (i=0; i<n; i++) {
1295:     if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs);
1296:     if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]);
1297:   }
1298:   PetscNew(&ilink);
1299:   if (splitname) {
1300:     PetscStrallocpy(splitname,&ilink->splitname);
1301:   } else {
1302:     PetscMalloc1(3,&ilink->splitname);
1303:     PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);
1304:   }
1305:   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 */
1306:   PetscMalloc1(n,&ilink->fields);
1307:   PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));
1308:   PetscMalloc1(n,&ilink->fields_col);
1309:   PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));

1311:   ilink->nfields = n;
1312:   ilink->next    = NULL;
1313:   KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1314:   KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1315:   PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1316:   KSPSetType(ilink->ksp,KSPPREONLY);
1317:   PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);

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

1322:   if (!next) {
1323:     jac->head       = ilink;
1324:     ilink->previous = NULL;
1325:   } else {
1326:     while (next->next) {
1327:       next = next->next;
1328:     }
1329:     next->next      = ilink;
1330:     ilink->previous = next;
1331:   }
1332:   jac->nsplits++;
1333:   return(0);
1334: }

1336: static PetscErrorCode  PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp)
1337: {
1338:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

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

1346:   (*subksp)[1] = jac->kspschur;
1347:   if (n) *n = jac->nsplits;
1348:   return(0);
1349: }

1351: static PetscErrorCode  PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1352: {
1353:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1354:   PetscErrorCode    ierr;
1355:   PetscInt          cnt   = 0;
1356:   PC_FieldSplitLink ilink = jac->head;

1359:   PetscMalloc1(jac->nsplits,subksp);
1360:   while (ilink) {
1361:     (*subksp)[cnt++] = ilink->ksp;
1362:     ilink            = ilink->next;
1363:   }
1364:   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);
1365:   if (n) *n = jac->nsplits;
1366:   return(0);
1367: }

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

1372:     Input Parameters:
1373: +   pc  - the preconditioner context
1374: +   is - the index set that defines the indices to which the fieldsplit is to be restricted

1376:     Level: advanced

1378: @*/
1379: PetscErrorCode  PCFieldSplitRestrictIS(PC pc,IS isy)
1380: {

1386:   PetscTryMethod(pc,"PCFieldSplitRestrictIS_C",(PC,IS),(pc,isy));
1387:   return(0);
1388: }


1391: static PetscErrorCode  PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy)
1392: {
1393:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1394:   PetscErrorCode    ierr;
1395:   PC_FieldSplitLink ilink = jac->head, next;
1396:   PetscInt          localsize,size,sizez,i;
1397:   const PetscInt    *ind, *indz;
1398:   PetscInt          *indc, *indcz;
1399:   PetscBool         flg;

1402:   ISGetLocalSize(isy,&localsize);
1403:   MPI_Scan(&localsize,&size,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)isy));
1404:   size -= localsize;
1405:   while(ilink) {
1406:     IS isrl,isr;
1407:     PC subpc;
1408:     ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl);
1409:     ISGetLocalSize(isrl,&localsize);
1410:     PetscMalloc1(localsize,&indc);
1411:     ISGetIndices(isrl,&ind);
1412:     PetscMemcpy(indc,ind,localsize*sizeof(PetscInt));
1413:     ISRestoreIndices(isrl,&ind);
1414:     ISDestroy(&isrl);
1415:     for (i=0; i<localsize; i++) *(indc+i) += size;
1416:     ISCreateGeneral(PetscObjectComm((PetscObject)isy),localsize,indc,PETSC_OWN_POINTER,&isr);
1417:     PetscObjectReference((PetscObject)isr);
1418:     ISDestroy(&ilink->is);
1419:     ilink->is     = isr;
1420:     PetscObjectReference((PetscObject)isr);
1421:     ISDestroy(&ilink->is_col);
1422:     ilink->is_col = isr;
1423:     ISDestroy(&isr);
1424:     KSPGetPC(ilink->ksp, &subpc);
1425:     PetscObjectTypeCompare((PetscObject)subpc,PCFIELDSPLIT,&flg);
1426:     if(flg) {
1427:       IS iszl,isz;
1428:       MPI_Comm comm;
1429:       ISGetLocalSize(ilink->is,&localsize);
1430:       comm   = PetscObjectComm((PetscObject)ilink->is);
1431:       ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl);
1432:       MPI_Scan(&localsize,&sizez,1,MPIU_INT,MPI_SUM,comm);
1433:       sizez -= localsize;
1434:       ISGetLocalSize(iszl,&localsize);
1435:       PetscMalloc1(localsize,&indcz);
1436:       ISGetIndices(iszl,&indz);
1437:       PetscMemcpy(indcz,indz,localsize*sizeof(PetscInt));
1438:       ISRestoreIndices(iszl,&indz);
1439:       ISDestroy(&iszl);
1440:       for (i=0; i<localsize; i++) *(indcz+i) += sizez;
1441:       ISCreateGeneral(comm,localsize,indcz,PETSC_OWN_POINTER,&isz);
1442:       PCFieldSplitRestrictIS(subpc,isz);
1443:       ISDestroy(&isz);
1444:     }
1445:     next = ilink->next;
1446:     ilink = next;
1447:   }
1448:   jac->isrestrict = PETSC_TRUE;
1449:   return(0);
1450: }

1452: static PetscErrorCode  PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is)
1453: {
1454:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1455:   PetscErrorCode    ierr;
1456:   PC_FieldSplitLink ilink, next = jac->head;
1457:   char              prefix[128];

1460:   if (jac->splitdefined) {
1461:     PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1462:     return(0);
1463:   }
1464:   PetscNew(&ilink);
1465:   if (splitname) {
1466:     PetscStrallocpy(splitname,&ilink->splitname);
1467:   } else {
1468:     PetscMalloc1(8,&ilink->splitname);
1469:     PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);
1470:   }
1471:   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 */
1472:   PetscObjectReference((PetscObject)is);
1473:   ISDestroy(&ilink->is);
1474:   ilink->is     = is;
1475:   PetscObjectReference((PetscObject)is);
1476:   ISDestroy(&ilink->is_col);
1477:   ilink->is_col = is;
1478:   ilink->next   = NULL;
1479:   KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1480:   KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1481:   PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1482:   KSPSetType(ilink->ksp,KSPPREONLY);
1483:   PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);

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

1488:   if (!next) {
1489:     jac->head       = ilink;
1490:     ilink->previous = NULL;
1491:   } else {
1492:     while (next->next) {
1493:       next = next->next;
1494:     }
1495:     next->next      = ilink;
1496:     ilink->previous = next;
1497:   }
1498:   jac->nsplits++;
1499:   return(0);
1500: }

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

1505:     Logically Collective on PC

1507:     Input Parameters:
1508: +   pc  - the preconditioner context
1509: .   splitname - name of this split, if NULL the number of the split is used
1510: .   n - the number of fields in this split
1511: -   fields - the fields in this split

1513:     Level: intermediate

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

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

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

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

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

1531: @*/
1532: PetscErrorCode  PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1533: {

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

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

1548:     Logically Collective on PC

1550:     Input Parameters:
1551: +   pc  - the preconditioner object
1552: -   flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from

1554:     Options Database:
1555: .     -pc_fieldsplit_diag_use_amat

1557:     Level: intermediate

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

1561: @*/
1562: PetscErrorCode  PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg)
1563: {
1564:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1565:   PetscBool      isfs;

1570:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1571:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1572:   jac->diag_use_amat = flg;
1573:   return(0);
1574: }

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

1579:     Logically Collective on PC

1581:     Input Parameters:
1582: .   pc  - the preconditioner object

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


1588:     Level: intermediate

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

1592: @*/
1593: PetscErrorCode  PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg)
1594: {
1595:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1596:   PetscBool      isfs;

1602:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1603:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1604:   *flg = jac->diag_use_amat;
1605:   return(0);
1606: }

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

1611:     Logically Collective on PC

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

1617:     Options Database:
1618: .     -pc_fieldsplit_off_diag_use_amat

1620:     Level: intermediate

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

1624: @*/
1625: PetscErrorCode  PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg)
1626: {
1627:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1628:   PetscBool      isfs;

1633:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1634:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1635:   jac->offdiag_use_amat = flg;
1636:   return(0);
1637: }

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

1642:     Logically Collective on PC

1644:     Input Parameters:
1645: .   pc  - the preconditioner object

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


1651:     Level: intermediate

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

1655: @*/
1656: PetscErrorCode  PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg)
1657: {
1658:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1659:   PetscBool      isfs;

1665:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1666:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1667:   *flg = jac->offdiag_use_amat;
1668:   return(0);
1669: }



1673: /*@C
1674:     PCFieldSplitSetIS - Sets the exact elements for field

1676:     Logically Collective on PC

1678:     Input Parameters:
1679: +   pc  - the preconditioner context
1680: .   splitname - name of this split, if NULL the number of the split is used
1681: -   is - the index set that defines the vector elements in this field


1684:     Notes:
1685:     Use PCFieldSplitSetFields(), for fields defined by strided types.

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

1690:     Level: intermediate

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

1694: @*/
1695: PetscErrorCode  PCFieldSplitSetIS(PC pc,const char splitname[],IS is)
1696: {

1703:   PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));
1704:   return(0);
1705: }

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

1710:     Logically Collective on PC

1712:     Input Parameters:
1713: +   pc  - the preconditioner context
1714: -   splitname - name of this split

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

1719:     Level: intermediate

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

1723: @*/
1724: PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is)
1725: {

1732:   {
1733:     PC_FieldSplit     *jac  = (PC_FieldSplit*) pc->data;
1734:     PC_FieldSplitLink ilink = jac->head;
1735:     PetscBool         found;

1737:     *is = NULL;
1738:     while (ilink) {
1739:       PetscStrcmp(ilink->splitname, splitname, &found);
1740:       if (found) {
1741:         *is = ilink->is;
1742:         break;
1743:       }
1744:       ilink = ilink->next;
1745:     }
1746:   }
1747:   return(0);
1748: }

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

1754:     Logically Collective on PC

1756:     Input Parameters:
1757: +   pc  - the preconditioner context
1758: -   bs - the block size

1760:     Level: intermediate

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

1764: @*/
1765: PetscErrorCode  PCFieldSplitSetBlockSize(PC pc,PetscInt bs)
1766: {

1772:   PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));
1773:   return(0);
1774: }

1776: /*@C
1777:    PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits

1779:    Collective on KSP

1781:    Input Parameter:
1782: .  pc - the preconditioner context

1784:    Output Parameters:
1785: +  n - the number of splits
1786: -  subksp - the array of KSP contexts

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

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

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


1799:    Level: advanced

1801: .seealso: PCFIELDSPLIT
1802: @*/
1803: PetscErrorCode  PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
1804: {

1810:   PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
1811:   return(0);
1812: }

1814: /*@
1815:     PCFieldSplitSetSchurPre -  Indicates what operator is used to construct the preconditioner for the Schur complement.
1816:       A11 matrix. Otherwise no preconditioner is used.

1818:     Collective on PC

1820:     Input Parameters:
1821: +   pc      - the preconditioner context
1822: .   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 
1823:               PC_FIELDSPLIT_SCHUR_PRE_SELFP, and PC_FIELDSPLIT_SCHUR_PRE_FULL
1824: -   userpre - matrix to use for preconditioning, or NULL

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

1829:     Notes:
1830: $    If ptype is
1831: $        a11 then the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner
1832: $             matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix
1833: $        self the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix:
1834: $             The only preconditioner that currently works with this symbolic respresentation matrix object is the PCLSC
1835: $             preconditioner
1836: $        user then the preconditioner for the Schur complement is generated from the user provided matrix (pre argument
1837: $             to this function).
1838: $        selfp then the preconditioning for the Schur complement is generated from an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01
1839: $             This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be
1840: $             lumped before extracting the diagonal using the additional option -fieldsplit_1_mat_schur_complement_ainv_type lump
1841: $        full then the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation computed internally by PFIELDSPLIT (this is expensive)
1842: $             useful mostly as a test that the Schur complement approach can work for your problem

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

1848:     Level: intermediate

1850: .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType,
1851:           MatSchurComplementSetAinvType(), PCLSC

1853: @*/
1854: PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
1855: {

1860:   PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));
1861:   return(0);
1862: }
1863: PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */

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

1869:     Logically Collective on PC

1871:     Input Parameters:
1872: .   pc      - the preconditioner context

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

1878:     Level: intermediate

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

1882: @*/
1883: PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
1884: {

1889:   PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));
1890:   return(0);
1891: }

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

1896:     Not collective

1898:     Input Parameter:
1899: .   pc      - the preconditioner context

1901:     Output Parameter:
1902: .   S       - the Schur complement matrix

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

1907:     Level: advanced

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

1911: @*/
1912: PetscErrorCode  PCFieldSplitSchurGetS(PC pc,Mat *S)
1913: {
1915:   const char*    t;
1916:   PetscBool      isfs;
1917:   PC_FieldSplit  *jac;

1921:   PetscObjectGetType((PetscObject)pc,&t);
1922:   PetscStrcmp(t,PCFIELDSPLIT,&isfs);
1923:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
1924:   jac = (PC_FieldSplit*)pc->data;
1925:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
1926:   if (S) *S = jac->schur;
1927:   return(0);
1928: }

1930: /*@
1931:     PCFieldSplitSchurRestoreS -  restores the MatSchurComplement object used by this PC

1933:     Not collective

1935:     Input Parameters:
1936: +   pc      - the preconditioner context
1937: .   S       - the Schur complement matrix

1939:     Level: advanced

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

1943: @*/
1944: PetscErrorCode  PCFieldSplitSchurRestoreS(PC pc,Mat *S)
1945: {
1947:   const char*    t;
1948:   PetscBool      isfs;
1949:   PC_FieldSplit  *jac;

1953:   PetscObjectGetType((PetscObject)pc,&t);
1954:   PetscStrcmp(t,PCFIELDSPLIT,&isfs);
1955:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
1956:   jac = (PC_FieldSplit*)pc->data;
1957:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
1958:   if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten");
1959:   return(0);
1960: }


1963: static PetscErrorCode  PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
1964: {
1965:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

1969:   jac->schurpre = ptype;
1970:   if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
1971:     MatDestroy(&jac->schur_user);
1972:     jac->schur_user = pre;
1973:     PetscObjectReference((PetscObject)jac->schur_user);
1974:   }
1975:   return(0);
1976: }

1978: static PetscErrorCode  PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
1979: {
1980:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

1983:   *ptype = jac->schurpre;
1984:   *pre   = jac->schur_user;
1985:   return(0);
1986: }

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

1991:     Collective on PC

1993:     Input Parameters:
1994: +   pc  - the preconditioner context
1995: -   ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default

1997:     Options Database:
1998: .     -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full


2001:     Level: intermediate

2003:     Notes:
2004:     The FULL factorization is

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

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

2012: $    If A and S are solved exactly
2013: $      *) FULL factorization is a direct solver.
2014: $      *) The preconditioned operator with LOWER or UPPER has all eigenvalues equal to 1 and minimal polynomial of degree 2, so KSPGMRES converges in 2 iterations.
2015: $      *) With DIAG, the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations.

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

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

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

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

2028: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCFieldSplitSetSchurScale()
2029: @*/
2030: PetscErrorCode  PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype)
2031: {

2036:   PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));
2037:   return(0);
2038: }

2040: static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype)
2041: {
2042:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2045:   jac->schurfactorization = ftype;
2046:   return(0);
2047: }

2049: /*@
2050:     PCFieldSplitSetSchurScale -  Controls the sign flip of S for PC_FIELDSPLIT_SCHUR_FACT_DIAG.

2052:     Collective on PC

2054:     Input Parameters:
2055: +   pc    - the preconditioner context
2056: -   scale - scaling factor for the Schur complement

2058:     Options Database:
2059: .     -pc_fieldsplit_schur_scale - default is -1.0

2061:     Level: intermediate

2063: .seealso: PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurFactType, PCFieldSplitSetSchurScale()
2064: @*/
2065: PetscErrorCode PCFieldSplitSetSchurScale(PC pc,PetscScalar scale)
2066: {

2072:   PetscTryMethod(pc,"PCFieldSplitSetSchurScale_C",(PC,PetscScalar),(pc,scale));
2073:   return(0);
2074: }

2076: static PetscErrorCode PCFieldSplitSetSchurScale_FieldSplit(PC pc,PetscScalar scale)
2077: {
2078:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2081:   jac->schurscale = scale;
2082:   return(0);
2083: }

2085: /*@C
2086:    PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement

2088:    Collective on KSP

2090:    Input Parameter:
2091: .  pc - the preconditioner context

2093:    Output Parameters:
2094: +  A00 - the (0,0) block
2095: .  A01 - the (0,1) block
2096: .  A10 - the (1,0) block
2097: -  A11 - the (1,1) block

2099:    Level: advanced

2101: .seealso: PCFIELDSPLIT
2102: @*/
2103: PetscErrorCode  PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11)
2104: {
2105:   PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;

2109:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach.");
2110:   if (A00) *A00 = jac->pmat[0];
2111:   if (A01) *A01 = jac->B;
2112:   if (A10) *A10 = jac->C;
2113:   if (A11) *A11 = jac->pmat[1];
2114:   return(0);
2115: }

2117: static PetscErrorCode  PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type)
2118: {
2119:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

2123:   jac->type = type;
2124:   if (type == PC_COMPOSITE_SCHUR) {
2125:     pc->ops->apply = PCApply_FieldSplit_Schur;
2126:     pc->ops->view  = PCView_FieldSplit_Schur;

2128:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);
2129:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);
2130:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);
2131:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);
2132:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",PCFieldSplitSetSchurScale_FieldSplit);

2134:   } else {
2135:     pc->ops->apply = PCApply_FieldSplit;
2136:     pc->ops->view  = PCView_FieldSplit;

2138:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2139:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);
2140:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);
2141:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);
2142:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",0);
2143:   }
2144:   return(0);
2145: }

2147: static PetscErrorCode  PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs)
2148: {
2149:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2152:   if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs);
2153:   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);
2154:   jac->bs = bs;
2155:   return(0);
2156: }

2158: /*@
2159:    PCFieldSplitSetType - Sets the type of fieldsplit preconditioner.

2161:    Collective on PC

2163:    Input Parameter:
2164: .  pc - the preconditioner context
2165: .  type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR

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

2170:    Level: Intermediate

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

2174: .seealso: PCCompositeSetType()

2176: @*/
2177: PetscErrorCode  PCFieldSplitSetType(PC pc,PCCompositeType type)
2178: {

2183:   PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));
2184:   return(0);
2185: }

2187: /*@
2188:   PCFieldSplitGetType - Gets the type of fieldsplit preconditioner.

2190:   Not collective

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

2195:   Output Parameter:
2196: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR

2198:   Level: Intermediate

2200: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
2201: .seealso: PCCompositeSetType()
2202: @*/
2203: PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
2204: {
2205:   PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;

2210:   *type = jac->type;
2211:   return(0);
2212: }

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

2217:    Logically Collective

2219:    Input Parameters:
2220: +  pc   - the preconditioner context
2221: -  flg  - boolean indicating whether to use field splits defined by the DM

2223:    Options Database Key:
2224: .  -pc_fieldsplit_dm_splits

2226:    Level: Intermediate

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

2230: .seealso: PCFieldSplitGetDMSplits()

2232: @*/
2233: PetscErrorCode  PCFieldSplitSetDMSplits(PC pc,PetscBool flg)
2234: {
2235:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
2236:   PetscBool      isfs;

2242:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2243:   if (isfs) {
2244:     jac->dm_splits = flg;
2245:   }
2246:   return(0);
2247: }


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

2253:    Logically Collective

2255:    Input Parameter:
2256: .  pc   - the preconditioner context

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

2261:    Level: Intermediate

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

2265: .seealso: PCFieldSplitSetDMSplits()

2267: @*/
2268: PetscErrorCode  PCFieldSplitGetDMSplits(PC pc,PetscBool* flg)
2269: {
2270:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
2271:   PetscBool      isfs;

2277:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2278:   if (isfs) {
2279:     if(flg) *flg = jac->dm_splits;
2280:   }
2281:   return(0);
2282: }

2284: /* -------------------------------------------------------------------------------------*/
2285: /*MC
2286:    PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
2287:                   fields or groups of fields. See the users manual section "Solving Block Matrices" for more details.

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

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

2295:    Level: intermediate

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

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

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

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

2318: $     For the Schur complement preconditioner if J = ( A00 A01 )
2319: $                                                    ( A10 A11 )
2320: $     the preconditioner using full factorization is
2321: $              ( I   -ksp(A00) A01 ) ( inv(A00)     0  ) (     I          0  )
2322: $              ( 0         I       ) (   0      ksp(S) ) ( -A10 ksp(A00)  I  )
2323:      where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_.  S is the Schur complement
2324: $              S = A11 - A10 ksp(A00) A01
2325:      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
2326:      in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0,
2327:      it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default
2328:      A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() for all the possible ways to construct the preconditioner for S.

2330:      The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above,
2331:      diag gives
2332: $              ( inv(A00)     0   )
2333: $              (   0      -ksp(S) )
2334:      note that slightly counter intuitively there is a negative in front of the ksp(S) so that the preconditioner is positive definite. For SPD matrices J, the sign flip
2335:      can be turned off with PCFieldSplitSetSchurScale() or by command line -pc_fieldsplit_schur_scale 1.0. The lower factorization is the inverse of
2336: $              (  A00   0 )
2337: $              (  A10   S )
2338:      where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of
2339: $              ( A00 A01 )
2340: $              (  0   S  )
2341:      where again the inverses of A00 and S are applied using KSPs.

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

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

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

2353:    Concepts: physics based preconditioners, block preconditioners

2355:    There is a nice discussion of block preconditioners in

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

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

2364: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC,
2365:            PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(),
2366:           MatSchurComplementSetAinvType(), PCFieldSplitSetSchurScale()
2367: M*/

2369: PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
2370: {
2372:   PC_FieldSplit  *jac;

2375:   PetscNewLog(pc,&jac);

2377:   jac->bs                 = -1;
2378:   jac->nsplits            = 0;
2379:   jac->type               = PC_COMPOSITE_MULTIPLICATIVE;
2380:   jac->schurpre           = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
2381:   jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
2382:   jac->schurscale         = -1.0;
2383:   jac->dm_splits          = PETSC_TRUE;

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

2387:   pc->ops->apply           = PCApply_FieldSplit;
2388:   pc->ops->applytranspose  = PCApplyTranspose_FieldSplit;
2389:   pc->ops->setup           = PCSetUp_FieldSplit;
2390:   pc->ops->reset           = PCReset_FieldSplit;
2391:   pc->ops->destroy         = PCDestroy_FieldSplit;
2392:   pc->ops->setfromoptions  = PCSetFromOptions_FieldSplit;
2393:   pc->ops->view            = PCView_FieldSplit;
2394:   pc->ops->applyrichardson = 0;

2396:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2397:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);
2398:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);
2399:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);
2400:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);
2401:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",PCFieldSplitRestrictIS_FieldSplit);
2402:   return(0);
2403: }