Actual source code: fieldsplit.c

petsc-master 2017-04-26
Report Typos and Errors


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

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

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

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

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

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

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

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


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

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

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

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

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

157:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
158:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
159:   if (iascii) {
160:     if (jac->bs > 0) {
161:       PetscViewerASCIIPrintf(viewer,"  FieldSplit with Schur preconditioner, blocksize = %D, factorization %s\n",jac->bs,PCFieldSplitSchurFactTypes[jac->schurfactorization]);
162:     } else {
163:       PetscViewerASCIIPrintf(viewer,"  FieldSplit with Schur preconditioner, factorization %s\n",PCFieldSplitSchurFactTypes[jac->schurfactorization]);
164:     }
165:     if (pc->useAmat) {
166:       PetscViewerASCIIPrintf(viewer,"  using Amat (not Pmat) as operator for blocks\n");
167:     }
168:     switch (jac->schurpre) {
169:     case PC_FIELDSPLIT_SCHUR_PRE_SELF:
170:       PetscViewerASCIIPrintf(viewer,"  Preconditioner for the Schur complement formed from S itself\n");
171:       break;
172:     case PC_FIELDSPLIT_SCHUR_PRE_SELFP:
173:       MatSchurComplementGetAinvType(jac->schur,&atype);
174:       PetscViewerASCIIPrintf(viewer,"  Preconditioner for the Schur complement formed from Sp, an assembled approximation to S, which uses A00's %sdiagonal's inverse\n",atype == MAT_SCHUR_COMPLEMENT_AINV_DIAG ? "" : "lumped ");break;
175:     case PC_FIELDSPLIT_SCHUR_PRE_A11:
176:       PetscViewerASCIIPrintf(viewer,"  Preconditioner for the Schur complement formed from A11\n");
177:       break;
178:     case PC_FIELDSPLIT_SCHUR_PRE_FULL:
179:       PetscViewerASCIIPrintf(viewer,"  Preconditioner for the Schur complement formed from the exact Schur complement\n");
180:       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: }

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

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

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

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

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

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

400:       if (stokes) {
401:         IS       zerodiags,rest;
402:         PetscInt nmin,nmax;

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

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

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

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

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

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

469:   PCFieldSplitSetDefaults(pc);
470:   nsplit = jac->nsplits;
471:   ilink  = jac->head;

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

477:     jac->issetup = PETSC_TRUE;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

692:       /* extract the A01 and A10 matrices */
693:       ilink = jac->head;
694:       ISComplement(ilink->is_col,rstart,rend,&ccis);
695:       if (jac->offdiag_use_amat) {
696:         MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
697:       } else {
698:         MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
699:       }
700:       ISDestroy(&ccis);
701:       ilink = ilink->next;
702:       ISComplement(ilink->is_col,rstart,rend,&ccis);
703:       if (jac->offdiag_use_amat) {
704:         MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
705:       } else {
706:         MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
707:       }
708:       ISDestroy(&ccis);

710:       /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
711:       MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);
712:       MatSetType(jac->schur,MATSCHURCOMPLEMENT);
713:       MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
714:       PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
715:       /* 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? */
716:       MatSetOptionsPrefix(jac->schur,schurmatprefix);
717:       MatSetFromOptions(jac->schur);
718:       MatGetNullSpace(jac->mat[1], &sp);
719:       if (sp) {
720:         MatSetNullSpace(jac->schur, sp);
721:       }

723:       PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);
724:       PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
725:       if (flg) {
726:         DM  dmInner;
727:         KSP kspInner;

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

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

753:       PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);
754:       PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
755:       if (flg) {
756:         DM dmInner;

758:         PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
759:         KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);
760:         KSPSetErrorIfNotConverged(jac->kspupper,pc->erroriffailure);
761:         KSPSetOptionsPrefix(jac->kspupper, schurprefix);
762:         KSPGetDM(jac->head->ksp, &dmInner);
763:         KSPSetDM(jac->kspupper, dmInner);
764:         KSPSetDMActive(jac->kspupper, PETSC_FALSE);
765:         KSPSetFromOptions(jac->kspupper);
766:         KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0]);
767:         VecDuplicate(jac->head->x, &jac->head->z);
768:       } else {
769:         jac->kspupper = jac->head->ksp;
770:         PetscObjectReference((PetscObject) jac->head->ksp);
771:       }

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

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

827:   jac->suboptionsset = PETSC_TRUE;
828:   return(0);
829: }

831: #define FieldSplitSplitSolveAdd(ilink,xx,yy) \
832:   (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
833:    VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
834:    PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\
835:    KSPSolve(ilink->ksp,ilink->x,ilink->y) ||                               \
836:    PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\
837:    VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) ||  \
838:    VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE))

840: static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y)
841: {
842:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
843:   PetscErrorCode    ierr;
844:   PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
845:   KSP               kspA   = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;

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

916:     PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
917:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
918:     PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
919:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
920:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);

922:     if (kspUpper == kspA) {
923:       MatMult(jac->B,ilinkD->y,ilinkA->y);
924:       VecAXPY(ilinkA->x,-1.0,ilinkA->y);
925:       PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
926:       KSPSolve(kspA,ilinkA->x,ilinkA->y);
927:       PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
928:     } else {
929:       PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
930:       KSPSolve(kspA,ilinkA->x,ilinkA->y);
931:       PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
932:       MatMult(jac->B,ilinkD->y,ilinkA->x);
933:       PetscLogEventBegin(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);
934:       KSPSolve(kspUpper,ilinkA->x,ilinkA->z);
935:       PetscLogEventEnd(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);
936:       VecAXPY(ilinkA->y,-1.0,ilinkA->z);
937:     }
938:     VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
939:     VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
940:   }
941:   return(0);
942: }

944: static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y)
945: {
946:   PC_FieldSplit      *jac = (PC_FieldSplit*)pc->data;
947:   PetscErrorCode     ierr;
948:   PC_FieldSplitLink  ilink = jac->head;
949:   PetscInt           cnt,bs;
950:   KSPConvergedReason reason;

953:   if (jac->type == PC_COMPOSITE_ADDITIVE) {
954:     if (jac->defaultsplit) {
955:       VecGetBlockSize(x,&bs);
956:       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);
957:       VecGetBlockSize(y,&bs);
958:       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);
959:       VecStrideGatherAll(x,jac->x,INSERT_VALUES);
960:       while (ilink) {
961:         PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
962:         KSPSolve(ilink->ksp,ilink->x,ilink->y);
963:         PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
964:         KSPGetConvergedReason(ilink->ksp,&reason);
965:         if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
966:           pc->failedreason = PC_SUBPC_ERROR;
967:         }
968:         ilink = ilink->next;
969:       }
970:       VecStrideScatterAll(jac->y,y,INSERT_VALUES);
971:     } else {
972:       VecSet(y,0.0);
973:       while (ilink) {
974:         FieldSplitSplitSolveAdd(ilink,x,y);
975:         KSPGetConvergedReason(ilink->ksp,&reason);
976:         if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
977:           pc->failedreason = PC_SUBPC_ERROR;
978:         }
979:         ilink = ilink->next;
980:       }
981:     }
982:   } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) {
983:     VecSet(y,0.0);
984:     /* solve on first block for first block variables */
985:     VecScatterBegin(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
986:     VecScatterEnd(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
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:     VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
995:     VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);

997:     /* compute the residual only onto second block variables using first block variables */
998:     MatMult(jac->Afield[1],ilink->y,ilink->next->x);
999:     ilink = ilink->next;
1000:     VecScale(ilink->x,-1.0);
1001:     VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1002:     VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);

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

1067: #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \
1068:   (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1069:    VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1070:    PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \
1071:    KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) ||                  \
1072:    PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \
1073:    VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \
1074:    VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE))

1076: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y)
1077: {
1078:   PC_FieldSplit      *jac = (PC_FieldSplit*)pc->data;
1079:   PetscErrorCode     ierr;
1080:   PC_FieldSplitLink  ilink = jac->head;
1081:   PetscInt           bs;
1082:   KSPConvergedReason reason;

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

1158: static PetscErrorCode PCReset_FieldSplit(PC pc)
1159: {
1160:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1161:   PetscErrorCode    ierr;
1162:   PC_FieldSplitLink ilink = jac->head,next;

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

1203: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
1204: {
1205:   PetscErrorCode    ierr;

1208:   PCReset_FieldSplit(pc);
1209:   PetscFree(pc->data);
1210:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);
1211:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);
1212:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);
1213:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);
1214:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);
1215:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);
1216:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);
1217:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);
1218:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",NULL);
1219:   return(0);
1220: }

1222: static PetscErrorCode PCSetFromOptions_FieldSplit(PetscOptionItems *PetscOptionsObject,PC pc)
1223: {
1224:   PetscErrorCode  ierr;
1225:   PetscInt        bs;
1226:   PetscBool       flg,stokes = PETSC_FALSE;
1227:   PC_FieldSplit   *jac = (PC_FieldSplit*)pc->data;
1228:   PCCompositeType ctype;

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

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

1269: /*------------------------------------------------------------------------------------*/

1271: static PetscErrorCode  PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1272: {
1273:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1274:   PetscErrorCode    ierr;
1275:   PC_FieldSplitLink ilink,next = jac->head;
1276:   char              prefix[128];
1277:   PetscInt          i;

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

1301:   ilink->nfields = n;
1302:   ilink->next    = NULL;
1303:   KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1304:   KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1305:   PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1306:   KSPSetType(ilink->ksp,KSPPREONLY);
1307:   PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);

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

1312:   if (!next) {
1313:     jac->head       = ilink;
1314:     ilink->previous = NULL;
1315:   } else {
1316:     while (next->next) {
1317:       next = next->next;
1318:     }
1319:     next->next      = ilink;
1320:     ilink->previous = next;
1321:   }
1322:   jac->nsplits++;
1323:   return(0);
1324: }

1326: static PetscErrorCode  PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp)
1327: {
1328:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

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

1336:   (*subksp)[1] = jac->kspschur;
1337:   if (n) *n = jac->nsplits;
1338:   return(0);
1339: }

1341: static PetscErrorCode  PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1342: {
1343:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1344:   PetscErrorCode    ierr;
1345:   PetscInt          cnt   = 0;
1346:   PC_FieldSplitLink ilink = jac->head;

1349:   PetscMalloc1(jac->nsplits,subksp);
1350:   while (ilink) {
1351:     (*subksp)[cnt++] = ilink->ksp;
1352:     ilink            = ilink->next;
1353:   }
1354:   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);
1355:   if (n) *n = jac->nsplits;
1356:   return(0);
1357: }

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

1362:     Input Parameters:
1363: +   pc  - the preconditioner context
1364: +   is - the index set that defines the indices to which the fieldsplit is to be restricted

1366:     Level: advanced

1368: @*/
1369: PetscErrorCode  PCFieldSplitRestrictIS(PC pc,IS isy)
1370: {

1376:   PetscTryMethod(pc,"PCFieldSplitRestrictIS_C",(PC,IS),(pc,isy));
1377:   return(0);
1378: }


1381: static PetscErrorCode  PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy)
1382: {
1383:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1384:   PetscErrorCode    ierr;
1385:   PC_FieldSplitLink ilink = jac->head, next;
1386:   PetscInt          localsize,size,sizez,i;
1387:   const PetscInt    *ind, *indz;
1388:   PetscInt          *indc, *indcz;
1389:   PetscBool         flg;

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

1442: static PetscErrorCode  PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is)
1443: {
1444:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1445:   PetscErrorCode    ierr;
1446:   PC_FieldSplitLink ilink, next = jac->head;
1447:   char              prefix[128];

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

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

1478:   if (!next) {
1479:     jac->head       = ilink;
1480:     ilink->previous = NULL;
1481:   } else {
1482:     while (next->next) {
1483:       next = next->next;
1484:     }
1485:     next->next      = ilink;
1486:     ilink->previous = next;
1487:   }
1488:   jac->nsplits++;
1489:   return(0);
1490: }

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

1495:     Logically Collective on PC

1497:     Input Parameters:
1498: +   pc  - the preconditioner context
1499: .   splitname - name of this split, if NULL the number of the split is used
1500: .   n - the number of fields in this split
1501: -   fields - the fields in this split

1503:     Level: intermediate

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

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

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

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

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

1521: @*/
1522: PetscErrorCode  PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1523: {

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

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

1538:     Logically Collective on PC

1540:     Input Parameters:
1541: +   pc  - the preconditioner object
1542: -   flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from

1544:     Options Database:
1545: .     -pc_fieldsplit_diag_use_amat

1547:     Level: intermediate

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

1551: @*/
1552: PetscErrorCode  PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg)
1553: {
1554:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1555:   PetscBool      isfs;

1560:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1561:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1562:   jac->diag_use_amat = flg;
1563:   return(0);
1564: }

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

1569:     Logically Collective on PC

1571:     Input Parameters:
1572: .   pc  - the preconditioner object

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


1578:     Level: intermediate

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

1582: @*/
1583: PetscErrorCode  PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg)
1584: {
1585:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1586:   PetscBool      isfs;

1592:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1593:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1594:   *flg = jac->diag_use_amat;
1595:   return(0);
1596: }

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

1601:     Logically Collective on PC

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

1607:     Options Database:
1608: .     -pc_fieldsplit_off_diag_use_amat

1610:     Level: intermediate

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

1614: @*/
1615: PetscErrorCode  PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg)
1616: {
1617:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1618:   PetscBool      isfs;

1623:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1624:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1625:   jac->offdiag_use_amat = flg;
1626:   return(0);
1627: }

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

1632:     Logically Collective on PC

1634:     Input Parameters:
1635: .   pc  - the preconditioner object

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


1641:     Level: intermediate

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

1645: @*/
1646: PetscErrorCode  PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg)
1647: {
1648:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1649:   PetscBool      isfs;

1655:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1656:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1657:   *flg = jac->offdiag_use_amat;
1658:   return(0);
1659: }



1663: /*@C
1664:     PCFieldSplitSetIS - Sets the exact elements for field

1666:     Logically Collective on PC

1668:     Input Parameters:
1669: +   pc  - the preconditioner context
1670: .   splitname - name of this split, if NULL the number of the split is used
1671: -   is - the index set that defines the vector elements in this field


1674:     Notes:
1675:     Use PCFieldSplitSetFields(), for fields defined by strided types.

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

1680:     Level: intermediate

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

1684: @*/
1685: PetscErrorCode  PCFieldSplitSetIS(PC pc,const char splitname[],IS is)
1686: {

1693:   PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));
1694:   return(0);
1695: }

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

1700:     Logically Collective on PC

1702:     Input Parameters:
1703: +   pc  - the preconditioner context
1704: -   splitname - name of this split

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

1709:     Level: intermediate

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

1713: @*/
1714: PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is)
1715: {

1722:   {
1723:     PC_FieldSplit     *jac  = (PC_FieldSplit*) pc->data;
1724:     PC_FieldSplitLink ilink = jac->head;
1725:     PetscBool         found;

1727:     *is = NULL;
1728:     while (ilink) {
1729:       PetscStrcmp(ilink->splitname, splitname, &found);
1730:       if (found) {
1731:         *is = ilink->is;
1732:         break;
1733:       }
1734:       ilink = ilink->next;
1735:     }
1736:   }
1737:   return(0);
1738: }

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

1744:     Logically Collective on PC

1746:     Input Parameters:
1747: +   pc  - the preconditioner context
1748: -   bs - the block size

1750:     Level: intermediate

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

1754: @*/
1755: PetscErrorCode  PCFieldSplitSetBlockSize(PC pc,PetscInt bs)
1756: {

1762:   PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));
1763:   return(0);
1764: }

1766: /*@C
1767:    PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits

1769:    Collective on KSP

1771:    Input Parameter:
1772: .  pc - the preconditioner context

1774:    Output Parameters:
1775: +  n - the number of splits
1776: -  subksp - the array of KSP contexts

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

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

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


1789:    Level: advanced

1791: .seealso: PCFIELDSPLIT
1792: @*/
1793: PetscErrorCode  PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
1794: {

1800:   PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
1801:   return(0);
1802: }

1804: /*@
1805:     PCFieldSplitSetSchurPre -  Indicates what operator is used to construct the preconditioner for the Schur complement.
1806:       A11 matrix. Otherwise no preconditioner is used.

1808:     Collective on PC

1810:     Input Parameters:
1811: +   pc      - the preconditioner context
1812: .   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 
1813:               PC_FIELDSPLIT_SCHUR_PRE_SELFP, and PC_FIELDSPLIT_SCHUR_PRE_FULL
1814: -   userpre - matrix to use for preconditioning, or NULL

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

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

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

1838:     Level: intermediate

1840: .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType,
1841:           MatSchurComplementSetAinvType(), PCLSC

1843: @*/
1844: PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
1845: {

1850:   PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));
1851:   return(0);
1852: }
1853: PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */

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

1859:     Logically Collective on PC

1861:     Input Parameters:
1862: .   pc      - the preconditioner context

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

1868:     Level: intermediate

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

1872: @*/
1873: PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
1874: {

1879:   PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));
1880:   return(0);
1881: }

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

1886:     Not collective

1888:     Input Parameter:
1889: .   pc      - the preconditioner context

1891:     Output Parameter:
1892: .   S       - the Schur complement matrix

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

1897:     Level: advanced

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

1901: @*/
1902: PetscErrorCode  PCFieldSplitSchurGetS(PC pc,Mat *S)
1903: {
1905:   const char*    t;
1906:   PetscBool      isfs;
1907:   PC_FieldSplit  *jac;

1911:   PetscObjectGetType((PetscObject)pc,&t);
1912:   PetscStrcmp(t,PCFIELDSPLIT,&isfs);
1913:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
1914:   jac = (PC_FieldSplit*)pc->data;
1915:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
1916:   if (S) *S = jac->schur;
1917:   return(0);
1918: }

1920: /*@
1921:     PCFieldSplitSchurRestoreS -  restores the MatSchurComplement object used by this PC

1923:     Not collective

1925:     Input Parameters:
1926: +   pc      - the preconditioner context
1927: .   S       - the Schur complement matrix

1929:     Level: advanced

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

1933: @*/
1934: PetscErrorCode  PCFieldSplitSchurRestoreS(PC pc,Mat *S)
1935: {
1937:   const char*    t;
1938:   PetscBool      isfs;
1939:   PC_FieldSplit  *jac;

1943:   PetscObjectGetType((PetscObject)pc,&t);
1944:   PetscStrcmp(t,PCFIELDSPLIT,&isfs);
1945:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
1946:   jac = (PC_FieldSplit*)pc->data;
1947:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
1948:   if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten");
1949:   return(0);
1950: }


1953: static PetscErrorCode  PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
1954: {
1955:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

1959:   jac->schurpre = ptype;
1960:   if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
1961:     MatDestroy(&jac->schur_user);
1962:     jac->schur_user = pre;
1963:     PetscObjectReference((PetscObject)jac->schur_user);
1964:   }
1965:   return(0);
1966: }

1968: static PetscErrorCode  PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
1969: {
1970:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

1973:   *ptype = jac->schurpre;
1974:   *pre   = jac->schur_user;
1975:   return(0);
1976: }

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

1981:     Collective on PC

1983:     Input Parameters:
1984: +   pc  - the preconditioner context
1985: -   ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default

1987:     Options Database:
1988: .     -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full


1991:     Level: intermediate

1993:     Notes:
1994:     The FULL factorization is

1996: $   (A   B)  = (1       0) (A   0) (1  Ainv*B)
1997: $   (C   D)    (C*Ainv  1) (0   S) (0     1  )

1999:     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,
2000:     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).

2002:     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
2003:     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
2004:     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,
2005:     the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations.

2007:     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
2008:     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).

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

2014: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType
2015: @*/
2016: PetscErrorCode  PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype)
2017: {

2022:   PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));
2023:   return(0);
2024: }

2026: static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype)
2027: {
2028:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2031:   jac->schurfactorization = ftype;
2032:   return(0);
2033: }

2035: /*@C
2036:    PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement

2038:    Collective on KSP

2040:    Input Parameter:
2041: .  pc - the preconditioner context

2043:    Output Parameters:
2044: +  A00 - the (0,0) block
2045: .  A01 - the (0,1) block
2046: .  A10 - the (1,0) block
2047: -  A11 - the (1,1) block

2049:    Level: advanced

2051: .seealso: PCFIELDSPLIT
2052: @*/
2053: PetscErrorCode  PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11)
2054: {
2055:   PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;

2059:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach.");
2060:   if (A00) *A00 = jac->pmat[0];
2061:   if (A01) *A01 = jac->B;
2062:   if (A10) *A10 = jac->C;
2063:   if (A11) *A11 = jac->pmat[1];
2064:   return(0);
2065: }

2067: static PetscErrorCode  PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type)
2068: {
2069:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

2073:   jac->type = type;
2074:   if (type == PC_COMPOSITE_SCHUR) {
2075:     pc->ops->apply = PCApply_FieldSplit_Schur;
2076:     pc->ops->view  = PCView_FieldSplit_Schur;

2078:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);
2079:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);
2080:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);
2081:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);

2083:   } else {
2084:     pc->ops->apply = PCApply_FieldSplit;
2085:     pc->ops->view  = PCView_FieldSplit;

2087:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2088:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);
2089:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);
2090:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);
2091:   }
2092:   return(0);
2093: }

2095: static PetscErrorCode  PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs)
2096: {
2097:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2100:   if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs);
2101:   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);
2102:   jac->bs = bs;
2103:   return(0);
2104: }

2106: /*@
2107:    PCFieldSplitSetType - Sets the type of fieldsplit preconditioner.

2109:    Collective on PC

2111:    Input Parameter:
2112: .  pc - the preconditioner context
2113: .  type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR

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

2118:    Level: Intermediate

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

2122: .seealso: PCCompositeSetType()

2124: @*/
2125: PetscErrorCode  PCFieldSplitSetType(PC pc,PCCompositeType type)
2126: {

2131:   PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));
2132:   return(0);
2133: }

2135: /*@
2136:   PCFieldSplitGetType - Gets the type of fieldsplit preconditioner.

2138:   Not collective

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

2143:   Output Parameter:
2144: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR

2146:   Level: Intermediate

2148: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
2149: .seealso: PCCompositeSetType()
2150: @*/
2151: PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
2152: {
2153:   PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;

2158:   *type = jac->type;
2159:   return(0);
2160: }

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

2165:    Logically Collective

2167:    Input Parameters:
2168: +  pc   - the preconditioner context
2169: -  flg  - boolean indicating whether to use field splits defined by the DM

2171:    Options Database Key:
2172: .  -pc_fieldsplit_dm_splits

2174:    Level: Intermediate

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

2178: .seealso: PCFieldSplitGetDMSplits()

2180: @*/
2181: PetscErrorCode  PCFieldSplitSetDMSplits(PC pc,PetscBool flg)
2182: {
2183:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
2184:   PetscBool      isfs;

2190:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2191:   if (isfs) {
2192:     jac->dm_splits = flg;
2193:   }
2194:   return(0);
2195: }


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

2201:    Logically Collective

2203:    Input Parameter:
2204: .  pc   - the preconditioner context

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

2209:    Level: Intermediate

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

2213: .seealso: PCFieldSplitSetDMSplits()

2215: @*/
2216: PetscErrorCode  PCFieldSplitGetDMSplits(PC pc,PetscBool* flg)
2217: {
2218:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
2219:   PetscBool      isfs;

2225:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2226:   if (isfs) {
2227:     if(flg) *flg = jac->dm_splits;
2228:   }
2229:   return(0);
2230: }

2232: /* -------------------------------------------------------------------------------------*/
2233: /*MC
2234:    PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
2235:                   fields or groups of fields. See the users manual section "Solving Block Matrices" for more details.

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

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

2243:    Level: intermediate

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

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

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

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

2266: $     For the Schur complement preconditioner if J = ( A00 A01 )
2267: $                                                    ( A10 A11 )
2268: $     the preconditioner using full factorization is
2269: $              ( I   -ksp(A00) A01 ) ( inv(A00)     0  ) (     I          0  )
2270: $              ( 0         I       ) (   0      ksp(S) ) ( -A10 ksp(A00)  I  )
2271:      where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_.  S is the Schur complement
2272: $              S = A11 - A10 ksp(A00) A01
2273:      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
2274:      in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0,
2275:      it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default
2276:      A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() for all the possible ways to construct the preconditioner for S.

2278:      The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above,
2279:      diag gives
2280: $              ( inv(A00)     0   )
2281: $              (   0      -ksp(S) )
2282:      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
2283: $              (  A00   0 )
2284: $              (  A10   S )
2285:      where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of
2286: $              ( A00 A01 )
2287: $              (  0   S  )
2288:      where again the inverses of A00 and S are applied using KSPs.

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

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

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

2300:    Concepts: physics based preconditioners, block preconditioners

2302:    There is a nice discussion of block preconditioners in

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

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

2311: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC,
2312:            PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(),
2313:            MatSchurComplementSetAinvType()
2314: M*/

2316: PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
2317: {
2319:   PC_FieldSplit  *jac;

2322:   PetscNewLog(pc,&jac);

2324:   jac->bs                 = -1;
2325:   jac->nsplits            = 0;
2326:   jac->type               = PC_COMPOSITE_MULTIPLICATIVE;
2327:   jac->schurpre           = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
2328:   jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
2329:   jac->dm_splits          = PETSC_TRUE;

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

2333:   pc->ops->apply           = PCApply_FieldSplit;
2334:   pc->ops->applytranspose  = PCApplyTranspose_FieldSplit;
2335:   pc->ops->setup           = PCSetUp_FieldSplit;
2336:   pc->ops->reset           = PCReset_FieldSplit;
2337:   pc->ops->destroy         = PCDestroy_FieldSplit;
2338:   pc->ops->setfromoptions  = PCSetFromOptions_FieldSplit;
2339:   pc->ops->view            = PCView_FieldSplit;
2340:   pc->ops->applyrichardson = 0;

2342:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2343:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);
2344:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);
2345:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);
2346:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);
2347:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",PCFieldSplitRestrictIS_FieldSplit);
2348:   return(0);
2349: }