Actual source code: fieldsplit.c

petsc-3.7.6 2017-04-24
Report Typos and Errors
  2: #include <petsc/private/pcimpl.h>     /*I "petscpc.h" I*/
  3: #include <petsc/private/kspimpl.h>
  4: #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,is_orig;
 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:   PC_FieldSplitLink         head;
 48:   PetscBool                 reset;                  /* indicates PCReset() has been last called on this object, hack */
 49:   PetscBool                 isrestrict;             /* indicates PCFieldSplitRestrictIS() has been last called on this object, hack */
 50:   PetscBool                 suboptionsset;          /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */
 51:   PetscBool                 dm_splits;              /* Whether to use DMCreateFieldDecomposition() whenever possible */
 52:   PetscBool                 diag_use_amat;          /* Whether to extract diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
 53:   PetscBool                 offdiag_use_amat;       /* Whether to extract off-diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
 54: } PC_FieldSplit;

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

408:       if (stokes) {
409:         IS       zerodiags,rest;
410:         PetscInt nmin,nmax;

412:         MatGetOwnershipRange(pc->mat,&nmin,&nmax);
413:         MatFindZeroDiagonals(pc->mat,&zerodiags);
414:         ISComplement(zerodiags,nmin,nmax,&rest);
415:         if (jac->reset) {
416:           jac->head->is       = rest;
417:           jac->head->next->is = zerodiags;
418:         } else {
419:           PCFieldSplitSetIS(pc,"0",rest);
420:           PCFieldSplitSetIS(pc,"1",zerodiags);
421:         }
422:         ISDestroy(&zerodiags);
423:         ISDestroy(&rest);
424:       } else if (coupling) {
425:         IS       coupling,rest;
426:         PetscInt nmin,nmax;

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

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


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

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

482: static PetscErrorCode PCSetUp_FieldSplit(PC pc)
483: {
484:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
485:   PetscErrorCode    ierr;
486:   PC_FieldSplitLink ilink;
487:   PetscInt          i,nsplit;
488:   PetscBool         sorted, sorted_col;

491:   PCFieldSplitSetDefaults(pc);
492:   nsplit = jac->nsplits;
493:   ilink  = jac->head;

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

499:     jac->issetup = PETSC_TRUE;

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

539:   ilink = jac->head;
540:   if (!jac->pmat) {
541:     Vec xtmp;

543:     MatCreateVecs(pc->pmat,&xtmp,NULL);
544:     PetscMalloc1(nsplit,&jac->pmat);
545:     PetscMalloc2(nsplit,&jac->x,nsplit,&jac->y);
546:     for (i=0; i<nsplit; i++) {
547:       MatNullSpace sp;

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

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

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

607:   /* Check for null space attached to IS */
608:   ilink = jac->head;
609:   for (i=0; i<nsplit; i++) {
610:     MatNullSpace sp;

612:     PetscObjectQuery((PetscObject) ilink->is, "nullspace", (PetscObject*) &sp);
613:     if (sp) {
614:       MatSetNullSpace(jac->mat[i], sp);
615:     }
616:     ilink = ilink->next;
617:   }

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

663:   if (jac->type == PC_COMPOSITE_SCHUR) {
664:     IS          ccis;
665:     PetscInt    rstart,rend;
666:     char        lscname[256];
667:     PetscObject LSC_L;

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

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

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

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

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

732:       /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
733:       MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);
734:       MatSetType(jac->schur,MATSCHURCOMPLEMENT);
735:       MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
736:       PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
737:       /* 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? */
738:       MatSetOptionsPrefix(jac->schur,schurmatprefix);
739:       MatSetFromOptions(jac->schur);
740:       MatGetNullSpace(jac->mat[1], &sp);
741:       if (sp) {
742:         MatSetNullSpace(jac->schur, sp);
743:       }

745:       PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);
746:       PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
747:       if (flg) {
748:         DM  dmInner;
749:         KSP kspInner;

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

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

775:       PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);
776:       PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
777:       if (flg) {
778:         DM dmInner;

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

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

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

849:   jac->suboptionsset = PETSC_TRUE;
850:   return(0);
851: }

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

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

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

940:     PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
941:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
942:     PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
943:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
944:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);

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

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

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

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

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

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

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

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

1188: static PetscErrorCode PCReset_FieldSplit(PC pc)
1189: {
1190:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1191:   PetscErrorCode    ierr;
1192:   PC_FieldSplitLink ilink = jac->head,next;

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

1234: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
1235: {
1236:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1237:   PetscErrorCode    ierr;
1238:   PC_FieldSplitLink ilink = jac->head,next;

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

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

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

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

1315: /*------------------------------------------------------------------------------------*/

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

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

1349:   ilink->nfields = n;
1350:   ilink->next    = NULL;
1351:   KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1352:   KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1353:   PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1354:   KSPSetType(ilink->ksp,KSPPREONLY);
1355:   PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);

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

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

1376: static PetscErrorCode  PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp)
1377: {
1378:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

1382:   PetscMalloc1(jac->nsplits,subksp);
1383:   MatSchurComplementGetKSP(jac->schur,*subksp);

1385:   (*subksp)[1] = jac->kspschur;
1386:   if (n) *n = jac->nsplits;
1387:   return(0);
1388: }

1392: static PetscErrorCode  PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1393: {
1394:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1395:   PetscErrorCode    ierr;
1396:   PetscInt          cnt   = 0;
1397:   PC_FieldSplitLink ilink = jac->head;

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

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

1415:     Input Parameters:
1416: +   pc  - the preconditioner context
1417: +   is - the index set that defines the indices to which the fieldsplit is to be restricted

1419:     Level: advanced

1421: @*/
1422: PetscErrorCode  PCFieldSplitRestrictIS(PC pc,IS isy)
1423: {

1429:   PetscTryMethod(pc,"PCFieldSplitRestrictIS_C",(PC,IS),(pc,isy));
1430:   return(0);
1431: }


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

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

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

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

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

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

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

1564:     Logically Collective on PC

1566:     Input Parameters:
1567: +   pc  - the preconditioner context
1568: .   splitname - name of this split, if NULL the number of the split is used
1569: .   n - the number of fields in this split
1570: -   fields - the fields in this split

1572:     Level: intermediate

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

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

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

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

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

1590: @*/
1591: PetscErrorCode  PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1592: {

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

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

1609:     Logically Collective on PC

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

1615:     Options Database:
1616: .     -pc_fieldsplit_diag_use_amat

1618:     Level: intermediate

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

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

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

1639: /*@
1640:     PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract 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 diagonal blocks from


1651:     Level: intermediate

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

1655: @*/
1656: PetscErrorCode  PCFieldSplitGetDiagUseAmat(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->diag_use_amat;
1668:   return(0);
1669: }

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

1676:     Logically Collective on PC

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

1682:     Options Database:
1683: .     -pc_fieldsplit_off_diag_use_amat

1685:     Level: intermediate

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

1689: @*/
1690: PetscErrorCode  PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg)
1691: {
1692:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1693:   PetscBool      isfs;

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

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

1709:     Logically Collective on PC

1711:     Input Parameters:
1712: .   pc  - the preconditioner object

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


1718:     Level: intermediate

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

1722: @*/
1723: PetscErrorCode  PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg)
1724: {
1725:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1726:   PetscBool      isfs;

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



1742: /*@C
1743:     PCFieldSplitSetIS - Sets the exact elements for field

1745:     Logically Collective on PC

1747:     Input Parameters:
1748: +   pc  - the preconditioner context
1749: .   splitname - name of this split, if NULL the number of the split is used
1750: -   is - the index set that defines the vector elements in this field


1753:     Notes:
1754:     Use PCFieldSplitSetFields(), for fields defined by strided types.

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

1759:     Level: intermediate

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

1763: @*/
1764: PetscErrorCode  PCFieldSplitSetIS(PC pc,const char splitname[],IS is)
1765: {

1772:   PetscUseMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));
1773:   return(0);
1774: }

1778: /*@
1779:     PCFieldSplitGetIS - Retrieves the elements for a field as an IS

1781:     Logically Collective on PC

1783:     Input Parameters:
1784: +   pc  - the preconditioner context
1785: -   splitname - name of this split

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

1790:     Level: intermediate

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

1794: @*/
1795: PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is)
1796: {

1803:   {
1804:     PC_FieldSplit     *jac  = (PC_FieldSplit*) pc->data;
1805:     PC_FieldSplitLink ilink = jac->head;
1806:     PetscBool         found;

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

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

1827:     Logically Collective on PC

1829:     Input Parameters:
1830: +   pc  - the preconditioner context
1831: -   bs - the block size

1833:     Level: intermediate

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

1837: @*/
1838: PetscErrorCode  PCFieldSplitSetBlockSize(PC pc,PetscInt bs)
1839: {

1845:   PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));
1846:   return(0);
1847: }

1851: /*@C
1852:    PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits

1854:    Collective on KSP

1856:    Input Parameter:
1857: .  pc - the preconditioner context

1859:    Output Parameters:
1860: +  n - the number of splits
1861: -  pc - the array of KSP contexts

1863:    Note:
1864:    After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user
1865:    (not the KSP just the array that contains them).

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

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


1874:    Level: advanced

1876: .seealso: PCFIELDSPLIT
1877: @*/
1878: PetscErrorCode  PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
1879: {

1885:   PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
1886:   return(0);
1887: }

1891: /*@
1892:     PCFieldSplitSetSchurPre -  Indicates if the Schur complement is preconditioned by a preconditioner constructed by the
1893:       A11 matrix. Otherwise no preconditioner is used.

1895:     Collective on PC

1897:     Input Parameters:
1898: +   pc      - the preconditioner context
1899: .   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 
1900:               PC_FIELDSPLIT_SCHUR_PRE_SELFP, and PC_FIELDSPLIT_SCHUR_PRE_FULL
1901: -   userpre - matrix to use for preconditioning, or NULL

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

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

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

1925:     Level: intermediate

1927: .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType,
1928:           MatSchurComplementSetAinvType(), PCLSC

1930: @*/
1931: PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
1932: {

1937:   PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));
1938:   return(0);
1939: }
1940: PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */

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

1948:     Logically Collective on PC

1950:     Input Parameters:
1951: .   pc      - the preconditioner context

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

1957:     Level: intermediate

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

1961: @*/
1962: PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
1963: {

1968:   PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));
1969:   return(0);
1970: }

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

1977:     Not collective

1979:     Input Parameter:
1980: .   pc      - the preconditioner context

1982:     Output Parameter:
1983: .   S       - the Schur complement matrix

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

1988:     Level: advanced

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

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

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

2013: /*@
2014:     PCFieldSplitSchurRestoreS -  restores the MatSchurComplement object used by this PC

2016:     Not collective

2018:     Input Parameters:
2019: +   pc      - the preconditioner context
2020: .   S       - the Schur complement matrix

2022:     Level: advanced

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

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

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


2048: static PetscErrorCode  PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
2049: {
2050:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

2054:   jac->schurpre = ptype;
2055:   if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
2056:     MatDestroy(&jac->schur_user);
2057:     jac->schur_user = pre;
2058:     PetscObjectReference((PetscObject)jac->schur_user);
2059:   }
2060:   return(0);
2061: }

2065: static PetscErrorCode  PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
2066: {
2067:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

2070:   *ptype = jac->schurpre;
2071:   *pre   = jac->schur_user;
2072:   return(0);
2073: }

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

2080:     Collective on PC

2082:     Input Parameters:
2083: +   pc  - the preconditioner context
2084: -   ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default

2086:     Options Database:
2087: .     -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full


2090:     Level: intermediate

2092:     Notes:
2093:     The FULL factorization is

2095: $   (A   B)  = (1       0) (A   0) (1  Ainv*B)
2096: $   (C   D)    (C*Ainv  1) (0   S) (0     1  )

2098:     where S = D - C*Ainv*B. In practice, the full factorization is applied via block triangular solves with the grouping L*(D*U). UPPER uses D*U, LOWER uses L*D,
2099:     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).

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

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

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

2113: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType
2114: @*/
2115: PetscErrorCode  PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype)
2116: {

2121:   PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));
2122:   return(0);
2123: }

2127: static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype)
2128: {
2129:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2132:   jac->schurfactorization = ftype;
2133:   return(0);
2134: }

2138: /*@C
2139:    PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement

2141:    Collective on KSP

2143:    Input Parameter:
2144: .  pc - the preconditioner context

2146:    Output Parameters:
2147: +  A00 - the (0,0) block
2148: .  A01 - the (0,1) block
2149: .  A10 - the (1,0) block
2150: -  A11 - the (1,1) block

2152:    Level: advanced

2154: .seealso: PCFIELDSPLIT
2155: @*/
2156: PetscErrorCode  PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11)
2157: {
2158:   PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;

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

2172: static PetscErrorCode  PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type)
2173: {
2174:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

2178:   jac->type = type;
2179:   if (type == PC_COMPOSITE_SCHUR) {
2180:     pc->ops->apply = PCApply_FieldSplit_Schur;
2181:     pc->ops->view  = PCView_FieldSplit_Schur;

2183:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);
2184:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);
2185:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);
2186:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);

2188:   } else {
2189:     pc->ops->apply = PCApply_FieldSplit;
2190:     pc->ops->view  = PCView_FieldSplit;

2192:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2193:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);
2194:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);
2195:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);
2196:   }
2197:   return(0);
2198: }

2202: static PetscErrorCode  PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs)
2203: {
2204:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

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

2215: /*@
2216:    PCFieldSplitSetType - Sets the type of fieldsplit preconditioner.

2218:    Collective on PC

2220:    Input Parameter:
2221: .  pc - the preconditioner context
2222: .  type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR

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

2227:    Level: Intermediate

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

2231: .seealso: PCCompositeSetType()

2233: @*/
2234: PetscErrorCode  PCFieldSplitSetType(PC pc,PCCompositeType type)
2235: {

2240:   PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));
2241:   return(0);
2242: }

2246: /*@
2247:   PCFieldSplitGetType - Gets the type of fieldsplit preconditioner.

2249:   Not collective

2251:   Input Parameter:
2252: . pc - the preconditioner context

2254:   Output Parameter:
2255: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR

2257:   Level: Intermediate

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

2269:   *type = jac->type;
2270:   return(0);
2271: }

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

2278:    Logically Collective

2280:    Input Parameters:
2281: +  pc   - the preconditioner context
2282: -  flg  - boolean indicating whether to use field splits defined by the DM

2284:    Options Database Key:
2285: .  -pc_fieldsplit_dm_splits

2287:    Level: Intermediate

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

2291: .seealso: PCFieldSplitGetDMSplits()

2293: @*/
2294: PetscErrorCode  PCFieldSplitSetDMSplits(PC pc,PetscBool flg)
2295: {
2296:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
2297:   PetscBool      isfs;

2303:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2304:   if (isfs) {
2305:     jac->dm_splits = flg;
2306:   }
2307:   return(0);
2308: }


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

2316:    Logically Collective

2318:    Input Parameter:
2319: .  pc   - the preconditioner context

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

2324:    Level: Intermediate

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

2328: .seealso: PCFieldSplitSetDMSplits()

2330: @*/
2331: PetscErrorCode  PCFieldSplitGetDMSplits(PC pc,PetscBool* flg)
2332: {
2333:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
2334:   PetscBool      isfs;

2340:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2341:   if (isfs) {
2342:     if(flg) *flg = jac->dm_splits;
2343:   }
2344:   return(0);
2345: }

2347: /* -------------------------------------------------------------------------------------*/
2348: /*MC
2349:    PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
2350:                   fields or groups of fields. See the users manual section "Solving Block Matrices" for more details.

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

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

2358:    Level: intermediate

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

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

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

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

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

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

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

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

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

2415:    Concepts: physics based preconditioners, block preconditioners

2417:    There is a nice discussion of block preconditioners in

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

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

2427: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC,
2428:            PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(),
2429:            MatSchurComplementSetAinvType()
2430: M*/

2434: PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
2435: {
2437:   PC_FieldSplit  *jac;

2440:   PetscNewLog(pc,&jac);

2442:   jac->bs                 = -1;
2443:   jac->nsplits            = 0;
2444:   jac->type               = PC_COMPOSITE_MULTIPLICATIVE;
2445:   jac->schurpre           = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
2446:   jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
2447:   jac->dm_splits          = PETSC_TRUE;

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

2451:   pc->ops->apply           = PCApply_FieldSplit;
2452:   pc->ops->applytranspose  = PCApplyTranspose_FieldSplit;
2453:   pc->ops->setup           = PCSetUp_FieldSplit;
2454:   pc->ops->reset           = PCReset_FieldSplit;
2455:   pc->ops->destroy         = PCDestroy_FieldSplit;
2456:   pc->ops->setfromoptions  = PCSetFromOptions_FieldSplit;
2457:   pc->ops->view            = PCView_FieldSplit;
2458:   pc->ops->applyrichardson = 0;

2460:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2461:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);
2462:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);
2463:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);
2464:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);
2465:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",PCFieldSplitRestrictIS_FieldSplit);
2466:   return(0);
2467: }