Actual source code: fieldsplit.c

petsc-master 2018-08-13
Report Typos and Errors


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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

403:       if (jac->detect) {
404:         IS       zerodiags,rest;
405:         PetscInt nmin,nmax;

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

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

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

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

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

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

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

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

480:     jac->issetup = PETSC_TRUE;

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

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

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

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

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

567:     for (i=0; i<nsplit; i++) {
568:       Mat pmat;

570:       /* Check for preconditioning matrix attached to IS */
571:       PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &pmat);
572:       if (!pmat) {
573:         MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,scall,&jac->pmat[i]);
574:       }
575:       ilink = ilink->next;
576:     }
577:   }
578:   if (jac->diag_use_amat) {
579:     ilink = jac->head;
580:     if (!jac->mat) {
581:       PetscMalloc1(nsplit,&jac->mat);
582:       for (i=0; i<nsplit; i++) {
583:         MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i]);
584:         ilink = ilink->next;
585:       }
586:     } else {
587:       MatReuse scall;
588:       if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
589:         for (i=0; i<nsplit; i++) {
590:           MatDestroy(&jac->mat[i]);
591:         }
592:         scall = MAT_INITIAL_MATRIX;
593:       } else scall = MAT_REUSE_MATRIX;

595:       for (i=0; i<nsplit; i++) {
596:         if (jac->mat[i]) {MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,scall,&jac->mat[i]);}
597:         ilink = ilink->next;
598:       }
599:     }
600:   } else {
601:     jac->mat = jac->pmat;
602:   }

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

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

616:   if (jac->type != PC_COMPOSITE_ADDITIVE  && jac->type != PC_COMPOSITE_SCHUR) {
617:     /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */
618:     /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */
619:     ilink = jac->head;
620:     if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) {
621:       /* 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 */
622:       if (!jac->Afield) {
623:         PetscCalloc1(nsplit,&jac->Afield);
624:         if (jac->offdiag_use_amat) {
625:           MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);
626:         } else {
627:           MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);
628:         }
629:       } else {
630:         MatReuse scall;
631:         if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
632:           for (i=0; i<nsplit; i++) {
633:             MatDestroy(&jac->Afield[1]);
634:           }
635:           scall = MAT_INITIAL_MATRIX;
636:         } else scall = MAT_REUSE_MATRIX;

638:         if (jac->offdiag_use_amat) {
639:           MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,scall,&jac->Afield[1]);
640:         } else {
641:           MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,scall,&jac->Afield[1]);
642:         }
643:       }
644:     } else {
645:       if (!jac->Afield) {
646:         PetscMalloc1(nsplit,&jac->Afield);
647:         for (i=0; i<nsplit; i++) {
648:           if (jac->offdiag_use_amat) {
649:             MatCreateSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
650:           } else {
651:             MatCreateSubMatrix(pc->pmat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
652:           }
653:           ilink = ilink->next;
654:         }
655:       } else {
656:         MatReuse scall;
657:         if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
658:           for (i=0; i<nsplit; i++) {
659:             MatDestroy(&jac->Afield[i]);
660:           }
661:           scall = MAT_INITIAL_MATRIX;
662:         } else scall = MAT_REUSE_MATRIX;

664:         for (i=0; i<nsplit; i++) {
665:           if (jac->offdiag_use_amat) {
666:             MatCreateSubMatrix(pc->mat,ilink->is,NULL,scall,&jac->Afield[i]);
667:           } else {
668:             MatCreateSubMatrix(pc->pmat,ilink->is,NULL,scall,&jac->Afield[i]);
669:           }
670:           ilink = ilink->next;
671:         }
672:       }
673:     }
674:   }

676:   if (jac->type == PC_COMPOSITE_SCHUR) {
677:     IS          ccis;
678:     PetscBool   isspd;
679:     PetscInt    rstart,rend;
680:     char        lscname[256];
681:     PetscObject LSC_L;

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

685:     /* If pc->mat is SPD, don't scale by -1 the Schur complement */
686:     if (jac->schurscale == (PetscScalar)-1.0) {
687:       MatGetOption(pc->pmat,MAT_SPD,&isspd);
688:       jac->schurscale = (isspd == PETSC_TRUE) ? 1.0 : -1.0;
689:     }

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

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

698:       MatSchurComplementGetKSP(jac->schur, &kspInner);
699:       ilink = jac->head;
700:       ISComplement(ilink->is_col,rstart,rend,&ccis);
701:       if (jac->offdiag_use_amat) {
702:         MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);
703:       } else {
704:         MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);
705:       }
706:       ISDestroy(&ccis);
707:       ilink = ilink->next;
708:       ISComplement(ilink->is_col,rstart,rend,&ccis);
709:       if (jac->offdiag_use_amat) {
710:         MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);
711:       } else {
712:         MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);
713:       }
714:       ISDestroy(&ccis);
715:       MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
716:       if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
717:         MatDestroy(&jac->schurp);
718:         MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);
719:       }
720:       if (kspA != kspInner) {
721:         KSPSetOperators(kspA,jac->mat[0],jac->pmat[0]);
722:       }
723:       if (kspUpper != kspA) {
724:         KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0]);
725:       }
726:       KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
727:     } else {
728:       const char   *Dprefix;
729:       char         schurprefix[256], schurmatprefix[256];
730:       char         schurtestoption[256];
731:       MatNullSpace sp;
732:       PetscBool    flg;
733:       KSP          kspt;

735:       /* extract the A01 and A10 matrices */
736:       ilink = jac->head;
737:       ISComplement(ilink->is_col,rstart,rend,&ccis);
738:       if (jac->offdiag_use_amat) {
739:         MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
740:       } else {
741:         MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
742:       }
743:       ISDestroy(&ccis);
744:       ilink = ilink->next;
745:       ISComplement(ilink->is_col,rstart,rend,&ccis);
746:       if (jac->offdiag_use_amat) {
747:         MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
748:       } else {
749:         MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
750:       }
751:       ISDestroy(&ccis);

753:       /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
754:       MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);
755:       MatSetType(jac->schur,MATSCHURCOMPLEMENT);
756:       MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
757:       PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
758:       MatSetOptionsPrefix(jac->schur,schurmatprefix);
759:       MatSchurComplementGetKSP(jac->schur,&kspt);
760:       KSPSetOptionsPrefix(kspt,schurmatprefix);

762:       /* Note: this is not true in general */
763:       MatGetNullSpace(jac->mat[1], &sp);
764:       if (sp) {
765:         MatSetNullSpace(jac->schur, sp);
766:       }

768:       PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);
769:       PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
770:       if (flg) {
771:         DM  dmInner;
772:         KSP kspInner;
773:         PC  pcInner;

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

781:         /* Set DM for new solver */
782:         KSPGetDM(jac->head->ksp, &dmInner);
783:         KSPSetDM(kspInner, dmInner);
784:         KSPSetDMActive(kspInner, PETSC_FALSE);

786:         /* Defaults to PCKSP as preconditioner */
787:         KSPGetPC(kspInner, &pcInner);
788:         PCSetType(pcInner, PCKSP);
789:         PCKSPSetKSP(pcInner, jac->head->ksp);
790:       } else {
791:          /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or
792:           * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact,
793:           * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for
794:           * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make
795:           * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used
796:           * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */
797:         KSPSetType(jac->head->ksp,KSPGMRES);
798:         MatSchurComplementSetKSP(jac->schur,jac->head->ksp);
799:       }
800:       KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0]);
801:       KSPSetFromOptions(jac->head->ksp);
802:       MatSetFromOptions(jac->schur);

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

809:         MatSchurComplementGetKSP(jac->schur, &kspInner);
810:         KSPGetPC(kspInner, &pcInner);
811:         PetscObjectTypeCompare((PetscObject)pcInner, PCKSP, &flg);
812:         if (flg) {
813:           KSP ksp;

815:           PCKSPGetKSP(pcInner, &ksp);
816:           if (ksp == jac->head->ksp) {
817:             PCSetUseAmat(pcInner, PETSC_TRUE);
818:           }
819:         }
820:       }
821:       PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);
822:       PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
823:       if (flg) {
824:         DM dmInner;

826:         PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
827:         KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);
828:         KSPSetErrorIfNotConverged(jac->kspupper,pc->erroriffailure);
829:         KSPSetOptionsPrefix(jac->kspupper, schurprefix);
830:         KSPGetDM(jac->head->ksp, &dmInner);
831:         KSPSetDM(jac->kspupper, dmInner);
832:         KSPSetDMActive(jac->kspupper, PETSC_FALSE);
833:         KSPSetFromOptions(jac->kspupper);
834:         KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0]);
835:         VecDuplicate(jac->head->x, &jac->head->z);
836:       } else {
837:         jac->kspupper = jac->head->ksp;
838:         PetscObjectReference((PetscObject) jac->head->ksp);
839:       }

841:       if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
842:         MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);
843:       }
844:       KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);
845:       KSPSetErrorIfNotConverged(jac->kspschur,pc->erroriffailure);
846:       PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);
847:       PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);
848:       if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
849:         PC pcschur;
850:         KSPGetPC(jac->kspschur,&pcschur);
851:         PCSetType(pcschur,PCNONE);
852:         /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
853:       } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) {
854:         MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user);
855:       }
856:       KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
857:       KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);
858:       KSPSetOptionsPrefix(jac->kspschur,         Dprefix);
859:       /* propagate DM */
860:       {
861:         DM sdm;
862:         KSPGetDM(jac->head->next->ksp, &sdm);
863:         if (sdm) {
864:           KSPSetDM(jac->kspschur, sdm);
865:           KSPSetDMActive(jac->kspschur, PETSC_FALSE);
866:         }
867:       }
868:       /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
869:       /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */
870:       KSPSetFromOptions(jac->kspschur);
871:     }
872:     MatAssemblyBegin(jac->schur,MAT_FINAL_ASSEMBLY);
873:     MatAssemblyEnd(jac->schur,MAT_FINAL_ASSEMBLY);

875:     /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */
876:     PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);
877:     PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);
878:     if (!LSC_L) {PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);}
879:     if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);}
880:     PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);
881:     PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);
882:     if (!LSC_L) {PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);}
883:     if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);}
884:   } else {
885:     /* set up the individual splits' PCs */
886:     i     = 0;
887:     ilink = jac->head;
888:     while (ilink) {
889:       KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i]);
890:       /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
891:       if (!jac->suboptionsset) {KSPSetFromOptions(ilink->ksp);}
892:       i++;
893:       ilink = ilink->next;
894:     }
895:   }

897:   jac->suboptionsset = PETSC_TRUE;
898:   return(0);
899: }

901: #define FieldSplitSplitSolveAdd(ilink,xx,yy) \
902:   (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
903:    VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
904:    PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\
905:    KSPSolve(ilink->ksp,ilink->x,ilink->y) ||                               \
906:    PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\
907:    VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) ||  \
908:    VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE))

910: static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y)
911: {
912:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
913:   PetscErrorCode    ierr;
914:   PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
915:   KSP               kspA   = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;

918:   switch (jac->schurfactorization) {
919:   case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
920:     /* [A00 0; 0 -S], positive definite, suitable for MINRES */
921:     VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
922:     VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
923:     VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
924:     PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
925:     KSPSolve(kspA,ilinkA->x,ilinkA->y);
926:     PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
927:     VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
928:     VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
929:     PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
930:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
931:     PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
932:     VecScale(ilinkD->y,jac->schurscale);
933:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
934:     VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
935:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
936:     break;
937:   case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
938:     /* [A00 0; A10 S], suitable for left preconditioning */
939:     VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
940:     VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
941:     PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
942:     KSPSolve(kspA,ilinkA->x,ilinkA->y);
943:     PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
944:     MatMult(jac->C,ilinkA->y,ilinkD->x);
945:     VecScale(ilinkD->x,-1.);
946:     VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
947:     VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
948:     VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
949:     PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
950:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
951:     PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
952:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
953:     VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
954:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
955:     break;
956:   case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
957:     /* [A00 A01; 0 S], suitable for right preconditioning */
958:     VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
959:     VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
960:     PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
961:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
962:     PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);    MatMult(jac->B,ilinkD->y,ilinkA->x);
963:     VecScale(ilinkA->x,-1.);
964:     VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
965:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
966:     VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
967:     PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
968:     KSPSolve(kspA,ilinkA->x,ilinkA->y);
969:     PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
970:     VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
971:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
972:     VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
973:     break;
974:   case PC_FIELDSPLIT_SCHUR_FACT_FULL:
975:     /* [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 */
976:     VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
977:     VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
978:     PetscLogEventBegin(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);
979:     KSPSolve(kspLower,ilinkA->x,ilinkA->y);
980:     PetscLogEventEnd(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);
981:     MatMult(jac->C,ilinkA->y,ilinkD->x);
982:     VecScale(ilinkD->x,-1.0);
983:     VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
984:     VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);

986:     PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
987:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
988:     PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);
989:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
990:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);

992:     if (kspUpper == kspA) {
993:       MatMult(jac->B,ilinkD->y,ilinkA->y);
994:       VecAXPY(ilinkA->x,-1.0,ilinkA->y);
995:       PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
996:       KSPSolve(kspA,ilinkA->x,ilinkA->y);
997:       PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
998:     } else {
999:       PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1000:       KSPSolve(kspA,ilinkA->x,ilinkA->y);
1001:       PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);
1002:       MatMult(jac->B,ilinkD->y,ilinkA->x);
1003:       PetscLogEventBegin(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);
1004:       KSPSolve(kspUpper,ilinkA->x,ilinkA->z);
1005:       PetscLogEventEnd(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);
1006:       VecAXPY(ilinkA->y,-1.0,ilinkA->z);
1007:     }
1008:     VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1009:     VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
1010:   }
1011:   return(0);
1012: }

1014: static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y)
1015: {
1016:   PC_FieldSplit      *jac = (PC_FieldSplit*)pc->data;
1017:   PetscErrorCode     ierr;
1018:   PC_FieldSplitLink  ilink = jac->head;
1019:   PetscInt           cnt,bs;
1020:   KSPConvergedReason reason;

1023:   if (jac->type == PC_COMPOSITE_ADDITIVE) {
1024:     if (jac->defaultsplit) {
1025:       VecGetBlockSize(x,&bs);
1026:       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);
1027:       VecGetBlockSize(y,&bs);
1028:       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);
1029:       VecStrideGatherAll(x,jac->x,INSERT_VALUES);
1030:       while (ilink) {
1031:         PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1032:         KSPSolve(ilink->ksp,ilink->x,ilink->y);
1033:         PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1034:         KSPGetConvergedReason(ilink->ksp,&reason);
1035:         if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1036:           pc->failedreason = PC_SUBPC_ERROR;
1037:         }
1038:         ilink = ilink->next;
1039:       }
1040:       VecStrideScatterAll(jac->y,y,INSERT_VALUES);
1041:     } else {
1042:       VecSet(y,0.0);
1043:       while (ilink) {
1044:         FieldSplitSplitSolveAdd(ilink,x,y);
1045:         KSPGetConvergedReason(ilink->ksp,&reason);
1046:         if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1047:           pc->failedreason = PC_SUBPC_ERROR;
1048:         }
1049:         ilink = ilink->next;
1050:       }
1051:     }
1052:   } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) {
1053:     VecSet(y,0.0);
1054:     /* solve on first block for first block variables */
1055:     VecScatterBegin(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
1056:     VecScatterEnd(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
1057:     PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1058:     KSPSolve(ilink->ksp,ilink->x,ilink->y);
1059:     PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1060:     KSPGetConvergedReason(ilink->ksp,&reason);
1061:     if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1062:       pc->failedreason = PC_SUBPC_ERROR;
1063:     }
1064:     VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1065:     VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);

1067:     /* compute the residual only onto second block variables using first block variables */
1068:     MatMult(jac->Afield[1],ilink->y,ilink->next->x);
1069:     ilink = ilink->next;
1070:     VecScale(ilink->x,-1.0);
1071:     VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1072:     VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);

1074:     /* solve on second block variables */
1075:     PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1076:     KSPSolve(ilink->ksp,ilink->x,ilink->y);
1077:     PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1078:     KSPGetConvergedReason(ilink->ksp,&reason);
1079:     if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1080:       pc->failedreason = PC_SUBPC_ERROR;
1081:     }
1082:     VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1083:     VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1084:   } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1085:     if (!jac->w1) {
1086:       VecDuplicate(x,&jac->w1);
1087:       VecDuplicate(x,&jac->w2);
1088:     }
1089:     VecSet(y,0.0);
1090:     FieldSplitSplitSolveAdd(ilink,x,y);
1091:     KSPGetConvergedReason(ilink->ksp,&reason);
1092:     if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1093:       pc->failedreason = PC_SUBPC_ERROR;
1094:     }
1095:     cnt  = 1;
1096:     while (ilink->next) {
1097:       ilink = ilink->next;
1098:       /* compute the residual only over the part of the vector needed */
1099:       MatMult(jac->Afield[cnt++],y,ilink->x);
1100:       VecScale(ilink->x,-1.0);
1101:       VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1102:       VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1103:       PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1104:       KSPSolve(ilink->ksp,ilink->x,ilink->y);
1105:       PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1106:       KSPGetConvergedReason(ilink->ksp,&reason);
1107:       if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1108:         pc->failedreason = PC_SUBPC_ERROR;
1109:       }
1110:       VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1111:       VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1112:     }
1113:     if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1114:       cnt -= 2;
1115:       while (ilink->previous) {
1116:         ilink = ilink->previous;
1117:         /* compute the residual only over the part of the vector needed */
1118:         MatMult(jac->Afield[cnt--],y,ilink->x);
1119:         VecScale(ilink->x,-1.0);
1120:         VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1121:         VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1122:         PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1123:         KSPSolve(ilink->ksp,ilink->x,ilink->y);
1124:         PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1125:         KSPGetConvergedReason(ilink->ksp,&reason);
1126:         if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1127:           pc->failedreason = PC_SUBPC_ERROR;
1128:         }
1129:         VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1130:         VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1131:       }
1132:     }
1133:   } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type);
1134:   return(0);
1135: }

1137: #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \
1138:   (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1139:    VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1140:    PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \
1141:    KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) ||                  \
1142:    PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \
1143:    VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \
1144:    VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE))

1146: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y)
1147: {
1148:   PC_FieldSplit      *jac = (PC_FieldSplit*)pc->data;
1149:   PetscErrorCode     ierr;
1150:   PC_FieldSplitLink  ilink = jac->head;
1151:   PetscInt           bs;
1152:   KSPConvergedReason reason;

1155:   if (jac->type == PC_COMPOSITE_ADDITIVE) {
1156:     if (jac->defaultsplit) {
1157:       VecGetBlockSize(x,&bs);
1158:       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);
1159:       VecGetBlockSize(y,&bs);
1160:       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);
1161:       VecStrideGatherAll(x,jac->x,INSERT_VALUES);
1162:       while (ilink) {
1163:         PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1164:         KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);
1165:         PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);
1166:         KSPGetConvergedReason(ilink->ksp,&reason);
1167:         if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1168:           pc->failedreason = PC_SUBPC_ERROR;
1169:         }
1170:         ilink = ilink->next;
1171:       }
1172:       VecStrideScatterAll(jac->y,y,INSERT_VALUES);
1173:     } else {
1174:       VecSet(y,0.0);
1175:       while (ilink) {
1176:         FieldSplitSplitSolveAddTranspose(ilink,x,y);
1177:         KSPGetConvergedReason(ilink->ksp,&reason);
1178:         if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1179:           pc->failedreason = PC_SUBPC_ERROR;
1180:         }
1181:         ilink = ilink->next;
1182:       }
1183:     }
1184:   } else {
1185:     if (!jac->w1) {
1186:       VecDuplicate(x,&jac->w1);
1187:       VecDuplicate(x,&jac->w2);
1188:     }
1189:     VecSet(y,0.0);
1190:     if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1191:       FieldSplitSplitSolveAddTranspose(ilink,x,y);
1192:       KSPGetConvergedReason(ilink->ksp,&reason);
1193:       if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1194:         pc->failedreason = PC_SUBPC_ERROR;
1195:       }
1196:       while (ilink->next) {
1197:         ilink = ilink->next;
1198:         MatMultTranspose(pc->mat,y,jac->w1);
1199:         VecWAXPY(jac->w2,-1.0,jac->w1,x);
1200:         FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1201:       }
1202:       while (ilink->previous) {
1203:         ilink = ilink->previous;
1204:         MatMultTranspose(pc->mat,y,jac->w1);
1205:         VecWAXPY(jac->w2,-1.0,jac->w1,x);
1206:         FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1207:       }
1208:     } else {
1209:       while (ilink->next) {   /* get to last entry in linked list */
1210:         ilink = ilink->next;
1211:       }
1212:       FieldSplitSplitSolveAddTranspose(ilink,x,y);
1213:       KSPGetConvergedReason(ilink->ksp,&reason);
1214:       if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1215:         pc->failedreason = PC_SUBPC_ERROR;
1216:       }
1217:       while (ilink->previous) {
1218:         ilink = ilink->previous;
1219:         MatMultTranspose(pc->mat,y,jac->w1);
1220:         VecWAXPY(jac->w2,-1.0,jac->w1,x);
1221:         FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1222:       }
1223:     }
1224:   }
1225:   return(0);
1226: }

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

1235:   while (ilink) {
1236:     KSPDestroy(&ilink->ksp);
1237:     VecDestroy(&ilink->x);
1238:     VecDestroy(&ilink->y);
1239:     VecDestroy(&ilink->z);
1240:     VecScatterDestroy(&ilink->sctx);
1241:     ISDestroy(&ilink->is);
1242:     ISDestroy(&ilink->is_col);
1243:     PetscFree(ilink->splitname);
1244:     PetscFree(ilink->fields);
1245:     PetscFree(ilink->fields_col);
1246:     next  = ilink->next;
1247:     PetscFree(ilink);
1248:     ilink = next;
1249:   }
1250:   jac->head = NULL;
1251:   PetscFree2(jac->x,jac->y);
1252:   if (jac->mat && jac->mat != jac->pmat) {
1253:     MatDestroyMatrices(jac->nsplits,&jac->mat);
1254:   } else if (jac->mat) {
1255:     jac->mat = NULL;
1256:   }
1257:   if (jac->pmat) {MatDestroyMatrices(jac->nsplits,&jac->pmat);}
1258:   if (jac->Afield) {MatDestroyMatrices(jac->nsplits,&jac->Afield);}
1259:   jac->nsplits = 0;
1260:   VecDestroy(&jac->w1);
1261:   VecDestroy(&jac->w2);
1262:   MatDestroy(&jac->schur);
1263:   MatDestroy(&jac->schurp);
1264:   MatDestroy(&jac->schur_user);
1265:   KSPDestroy(&jac->kspschur);
1266:   KSPDestroy(&jac->kspupper);
1267:   MatDestroy(&jac->B);
1268:   MatDestroy(&jac->C);
1269:   jac->isrestrict = PETSC_FALSE;
1270:   return(0);
1271: }

1273: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
1274: {
1275:   PetscErrorCode    ierr;

1278:   PCReset_FieldSplit(pc);
1279:   PetscFree(pc->data);
1280:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",NULL);
1281:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);
1282:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);
1283:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);
1284:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);
1285:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);
1286:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);
1287:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);
1288:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);
1289:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",NULL);
1290:   return(0);
1291: }

1293: static PetscErrorCode PCSetFromOptions_FieldSplit(PetscOptionItems *PetscOptionsObject,PC pc)
1294: {
1295:   PetscErrorCode  ierr;
1296:   PetscInt        bs;
1297:   PetscBool       flg;
1298:   PC_FieldSplit   *jac = (PC_FieldSplit*)pc->data;
1299:   PCCompositeType ctype;

1302:   PetscOptionsHead(PetscOptionsObject,"FieldSplit options");
1303:   PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);
1304:   PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);
1305:   if (flg) {
1306:     PCFieldSplitSetBlockSize(pc,bs);
1307:   }
1308:   jac->diag_use_amat = pc->useAmat;
1309:   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);
1310:   jac->offdiag_use_amat = pc->useAmat;
1311:   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);
1312:   PetscOptionsBool("-pc_fieldsplit_detect_saddle_point","Form 2-way split by detecting zero diagonal entries", "PCFieldSplitSetDetectSaddlePoint",jac->detect,&jac->detect,NULL);
1313:   PCFieldSplitSetDetectSaddlePoint(pc,jac->detect); /* Sets split type and Schur PC type */
1314:   PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);
1315:   if (flg) {
1316:     PCFieldSplitSetType(pc,ctype);
1317:   }
1318:   /* Only setup fields once */
1319:   if ((jac->bs > 0) && (jac->nsplits == 0)) {
1320:     /* only allow user to set fields from command line if bs is already known.
1321:        otherwise user can set them in PCFieldSplitSetDefaults() */
1322:     PCFieldSplitSetRuntimeSplits_Private(pc);
1323:     if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
1324:   }
1325:   if (jac->type == PC_COMPOSITE_SCHUR) {
1326:     PetscOptionsGetEnum(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);
1327:     if (flg) {PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");}
1328:     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);
1329:     PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);
1330:     PetscOptionsScalar("-pc_fieldsplit_schur_scale","Scale Schur complement","PCFieldSplitSetSchurScale",jac->schurscale,&jac->schurscale,NULL);
1331:   }
1332:   PetscOptionsTail();
1333:   return(0);
1334: }

1336: /*------------------------------------------------------------------------------------*/

1338: static PetscErrorCode  PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1339: {
1340:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1341:   PetscErrorCode    ierr;
1342:   PC_FieldSplitLink ilink,next = jac->head;
1343:   char              prefix[128];
1344:   PetscInt          i;

1347:   if (jac->splitdefined) {
1348:     PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1349:     return(0);
1350:   }
1351:   for (i=0; i<n; i++) {
1352:     if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs);
1353:     if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]);
1354:   }
1355:   PetscNew(&ilink);
1356:   if (splitname) {
1357:     PetscStrallocpy(splitname,&ilink->splitname);
1358:   } else {
1359:     PetscMalloc1(3,&ilink->splitname);
1360:     PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);
1361:   }
1362:   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 */
1363:   PetscMalloc1(n,&ilink->fields);
1364:   PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));
1365:   PetscMalloc1(n,&ilink->fields_col);
1366:   PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));

1368:   ilink->nfields = n;
1369:   ilink->next    = NULL;
1370:   KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1371:   KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1372:   PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1373:   KSPSetType(ilink->ksp,KSPPREONLY);
1374:   PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);

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

1379:   if (!next) {
1380:     jac->head       = ilink;
1381:     ilink->previous = NULL;
1382:   } else {
1383:     while (next->next) {
1384:       next = next->next;
1385:     }
1386:     next->next      = ilink;
1387:     ilink->previous = next;
1388:   }
1389:   jac->nsplits++;
1390:   return(0);
1391: }

1393: static PetscErrorCode  PCFieldSplitSchurGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1394: {
1395:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

1399:   *subksp = NULL;
1400:   if (n) *n = 0;
1401:   if (jac->type == PC_COMPOSITE_SCHUR) {
1402:     PetscInt nn;

1404:     if (!jac->schur) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitSchurGetSubKSP()");
1405:     if (jac->nsplits != 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unexpected number of splits %D != 2",jac->nsplits);
1406:     nn   = jac->nsplits + (jac->kspupper != jac->head->ksp ? 1 : 0);
1407:     PetscMalloc1(nn,subksp);
1408:     (*subksp)[0] = jac->head->ksp;
1409:     (*subksp)[1] = jac->kspschur;
1410:     if (jac->kspupper != jac->head->ksp) (*subksp)[2] = jac->kspupper;
1411:     if (n) *n = nn;
1412:   }
1413:   return(0);
1414: }

1416: static PetscErrorCode  PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp)
1417: {
1418:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

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

1426:   (*subksp)[1] = jac->kspschur;
1427:   if (n) *n = jac->nsplits;
1428:   return(0);
1429: }

1431: static PetscErrorCode  PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1432: {
1433:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1434:   PetscErrorCode    ierr;
1435:   PetscInt          cnt   = 0;
1436:   PC_FieldSplitLink ilink = jac->head;

1439:   PetscMalloc1(jac->nsplits,subksp);
1440:   while (ilink) {
1441:     (*subksp)[cnt++] = ilink->ksp;
1442:     ilink            = ilink->next;
1443:   }
1444:   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);
1445:   if (n) *n = jac->nsplits;
1446:   return(0);
1447: }

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

1452:     Input Parameters:
1453: +   pc  - the preconditioner context
1454: +   is - the index set that defines the indices to which the fieldsplit is to be restricted

1456:     Level: advanced

1458: @*/
1459: PetscErrorCode  PCFieldSplitRestrictIS(PC pc,IS isy)
1460: {

1466:   PetscTryMethod(pc,"PCFieldSplitRestrictIS_C",(PC,IS),(pc,isy));
1467:   return(0);
1468: }


1471: static PetscErrorCode  PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy)
1472: {
1473:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1474:   PetscErrorCode    ierr;
1475:   PC_FieldSplitLink ilink = jac->head, next;
1476:   PetscInt          localsize,size,sizez,i;
1477:   const PetscInt    *ind, *indz;
1478:   PetscInt          *indc, *indcz;
1479:   PetscBool         flg;

1482:   ISGetLocalSize(isy,&localsize);
1483:   MPI_Scan(&localsize,&size,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)isy));
1484:   size -= localsize;
1485:   while(ilink) {
1486:     IS isrl,isr;
1487:     PC subpc;
1488:     ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl);
1489:     ISGetLocalSize(isrl,&localsize);
1490:     PetscMalloc1(localsize,&indc);
1491:     ISGetIndices(isrl,&ind);
1492:     PetscMemcpy(indc,ind,localsize*sizeof(PetscInt));
1493:     ISRestoreIndices(isrl,&ind);
1494:     ISDestroy(&isrl);
1495:     for (i=0; i<localsize; i++) *(indc+i) += size;
1496:     ISCreateGeneral(PetscObjectComm((PetscObject)isy),localsize,indc,PETSC_OWN_POINTER,&isr);
1497:     PetscObjectReference((PetscObject)isr);
1498:     ISDestroy(&ilink->is);
1499:     ilink->is     = isr;
1500:     PetscObjectReference((PetscObject)isr);
1501:     ISDestroy(&ilink->is_col);
1502:     ilink->is_col = isr;
1503:     ISDestroy(&isr);
1504:     KSPGetPC(ilink->ksp, &subpc);
1505:     PetscObjectTypeCompare((PetscObject)subpc,PCFIELDSPLIT,&flg);
1506:     if(flg) {
1507:       IS iszl,isz;
1508:       MPI_Comm comm;
1509:       ISGetLocalSize(ilink->is,&localsize);
1510:       comm   = PetscObjectComm((PetscObject)ilink->is);
1511:       ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl);
1512:       MPI_Scan(&localsize,&sizez,1,MPIU_INT,MPI_SUM,comm);
1513:       sizez -= localsize;
1514:       ISGetLocalSize(iszl,&localsize);
1515:       PetscMalloc1(localsize,&indcz);
1516:       ISGetIndices(iszl,&indz);
1517:       PetscMemcpy(indcz,indz,localsize*sizeof(PetscInt));
1518:       ISRestoreIndices(iszl,&indz);
1519:       ISDestroy(&iszl);
1520:       for (i=0; i<localsize; i++) *(indcz+i) += sizez;
1521:       ISCreateGeneral(comm,localsize,indcz,PETSC_OWN_POINTER,&isz);
1522:       PCFieldSplitRestrictIS(subpc,isz);
1523:       ISDestroy(&isz);
1524:     }
1525:     next = ilink->next;
1526:     ilink = next;
1527:   }
1528:   jac->isrestrict = PETSC_TRUE;
1529:   return(0);
1530: }

1532: static PetscErrorCode  PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is)
1533: {
1534:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1535:   PetscErrorCode    ierr;
1536:   PC_FieldSplitLink ilink, next = jac->head;
1537:   char              prefix[128];

1540:   if (jac->splitdefined) {
1541:     PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1542:     return(0);
1543:   }
1544:   PetscNew(&ilink);
1545:   if (splitname) {
1546:     PetscStrallocpy(splitname,&ilink->splitname);
1547:   } else {
1548:     PetscMalloc1(8,&ilink->splitname);
1549:     PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);
1550:   }
1551:   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 */
1552:   PetscObjectReference((PetscObject)is);
1553:   ISDestroy(&ilink->is);
1554:   ilink->is     = is;
1555:   PetscObjectReference((PetscObject)is);
1556:   ISDestroy(&ilink->is_col);
1557:   ilink->is_col = is;
1558:   ilink->next   = NULL;
1559:   KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1560:   KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1561:   PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1562:   KSPSetType(ilink->ksp,KSPPREONLY);
1563:   PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);

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

1568:   if (!next) {
1569:     jac->head       = ilink;
1570:     ilink->previous = NULL;
1571:   } else {
1572:     while (next->next) {
1573:       next = next->next;
1574:     }
1575:     next->next      = ilink;
1576:     ilink->previous = next;
1577:   }
1578:   jac->nsplits++;
1579:   return(0);
1580: }

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

1585:     Logically Collective on PC

1587:     Input Parameters:
1588: +   pc  - the preconditioner context
1589: .   splitname - name of this split, if NULL the number of the split is used
1590: .   n - the number of fields in this split
1591: -   fields - the fields in this split

1593:     Level: intermediate

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

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

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

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

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

1612: @*/
1613: PetscErrorCode  PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1614: {

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

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

1629:     Logically Collective on PC

1631:     Input Parameters:
1632: +   pc  - the preconditioner object
1633: -   flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from

1635:     Options Database:
1636: .     -pc_fieldsplit_diag_use_amat

1638:     Level: intermediate

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

1642: @*/
1643: PetscErrorCode  PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg)
1644: {
1645:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1646:   PetscBool      isfs;

1651:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1652:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1653:   jac->diag_use_amat = flg;
1654:   return(0);
1655: }

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

1660:     Logically Collective on PC

1662:     Input Parameters:
1663: .   pc  - the preconditioner object

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


1669:     Level: intermediate

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

1673: @*/
1674: PetscErrorCode  PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg)
1675: {
1676:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1677:   PetscBool      isfs;

1683:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1684:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1685:   *flg = jac->diag_use_amat;
1686:   return(0);
1687: }

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

1692:     Logically Collective on PC

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

1698:     Options Database:
1699: .     -pc_fieldsplit_off_diag_use_amat

1701:     Level: intermediate

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

1705: @*/
1706: PetscErrorCode  PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg)
1707: {
1708:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1709:   PetscBool      isfs;

1714:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1715:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1716:   jac->offdiag_use_amat = flg;
1717:   return(0);
1718: }

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

1723:     Logically Collective on PC

1725:     Input Parameters:
1726: .   pc  - the preconditioner object

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


1732:     Level: intermediate

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

1736: @*/
1737: PetscErrorCode  PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg)
1738: {
1739:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1740:   PetscBool      isfs;

1746:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1747:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1748:   *flg = jac->offdiag_use_amat;
1749:   return(0);
1750: }



1754: /*@C
1755:     PCFieldSplitSetIS - Sets the exact elements for field

1757:     Logically Collective on PC

1759:     Input Parameters:
1760: +   pc  - the preconditioner context
1761: .   splitname - name of this split, if NULL the number of the split is used
1762: -   is - the index set that defines the vector elements in this field


1765:     Notes:
1766:     Use PCFieldSplitSetFields(), for fields defined by strided types.

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

1771:     Level: intermediate

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

1775: @*/
1776: PetscErrorCode  PCFieldSplitSetIS(PC pc,const char splitname[],IS is)
1777: {

1784:   PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));
1785:   return(0);
1786: }

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

1791:     Logically Collective on PC

1793:     Input Parameters:
1794: +   pc  - the preconditioner context
1795: -   splitname - name of this split

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

1800:     Level: intermediate

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

1804: @*/
1805: PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is)
1806: {

1813:   {
1814:     PC_FieldSplit     *jac  = (PC_FieldSplit*) pc->data;
1815:     PC_FieldSplitLink ilink = jac->head;
1816:     PetscBool         found;

1818:     *is = NULL;
1819:     while (ilink) {
1820:       PetscStrcmp(ilink->splitname, splitname, &found);
1821:       if (found) {
1822:         *is = ilink->is;
1823:         break;
1824:       }
1825:       ilink = ilink->next;
1826:     }
1827:   }
1828:   return(0);
1829: }

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

1835:     Logically Collective on PC

1837:     Input Parameters:
1838: +   pc  - the preconditioner context
1839: -   bs - the block size

1841:     Level: intermediate

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

1845: @*/
1846: PetscErrorCode  PCFieldSplitSetBlockSize(PC pc,PetscInt bs)
1847: {

1853:   PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));
1854:   return(0);
1855: }

1857: /*@C
1858:    PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits

1860:    Collective on KSP

1862:    Input Parameter:
1863: .  pc - the preconditioner context

1865:    Output Parameters:
1866: +  n - the number of splits
1867: -  subksp - the array of KSP contexts

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

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

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

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


1884:    Level: advanced

1886: .seealso: PCFIELDSPLIT
1887: @*/
1888: PetscErrorCode  PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
1889: {

1895:   PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
1896:   return(0);
1897: }

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

1902:    Collective on KSP

1904:    Input Parameter:
1905: .  pc - the preconditioner context

1907:    Output Parameters:
1908: +  n - the number of splits
1909: -  subksp - the array of KSP contexts

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

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

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

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

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

1928:    Level: advanced

1930: .seealso: PCFIELDSPLIT
1931: @*/
1932: PetscErrorCode  PCFieldSplitSchurGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
1933: {

1939:   PetscUseMethod(pc,"PCFieldSplitSchurGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
1940:   return(0);
1941: }

1943: /*@
1944:     PCFieldSplitSetSchurPre -  Indicates what operator is used to construct the preconditioner for the Schur complement.
1945:       A11 matrix. Otherwise no preconditioner is used.

1947:     Collective on PC

1949:     Input Parameters:
1950: +   pc      - the preconditioner context
1951: .   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
1952:               PC_FIELDSPLIT_SCHUR_PRE_SELFP, and PC_FIELDSPLIT_SCHUR_PRE_FULL
1953: -   userpre - matrix to use for preconditioning, or NULL

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

1958:     Notes:
1959: $    If ptype is
1960: $        a11 then the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner
1961: $             matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix
1962: $        self the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix:
1963: $             The only preconditioner that currently works with this symbolic respresentation matrix object is the PCLSC
1964: $             preconditioner
1965: $        user then the preconditioner for the Schur complement is generated from the user provided matrix (pre argument
1966: $             to this function).
1967: $        selfp then the preconditioning for the Schur complement is generated from an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01
1968: $             This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be
1969: $             lumped before extracting the diagonal using the additional option -fieldsplit_1_mat_schur_complement_ainv_type lump
1970: $        full then the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation computed internally by PCFIELDSPLIT (this is expensive)
1971: $             useful mostly as a test that the Schur complement approach can work for your problem

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

1977:     Level: intermediate

1979: .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType,
1980:           MatSchurComplementSetAinvType(), PCLSC

1982: @*/
1983: PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
1984: {

1989:   PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));
1990:   return(0);
1991: }

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

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

1999:     Logically Collective on PC

2001:     Input Parameters:
2002: .   pc      - the preconditioner context

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

2008:     Level: intermediate

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

2012: @*/
2013: PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
2014: {

2019:   PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));
2020:   return(0);
2021: }

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

2026:     Not collective

2028:     Input Parameter:
2029: .   pc      - the preconditioner context

2031:     Output Parameter:
2032: .   S       - the Schur complement matrix

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

2037:     Level: advanced

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

2041: @*/
2042: PetscErrorCode  PCFieldSplitSchurGetS(PC pc,Mat *S)
2043: {
2045:   const char*    t;
2046:   PetscBool      isfs;
2047:   PC_FieldSplit  *jac;

2051:   PetscObjectGetType((PetscObject)pc,&t);
2052:   PetscStrcmp(t,PCFIELDSPLIT,&isfs);
2053:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
2054:   jac = (PC_FieldSplit*)pc->data;
2055:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
2056:   if (S) *S = jac->schur;
2057:   return(0);
2058: }

2060: /*@
2061:     PCFieldSplitSchurRestoreS -  restores the MatSchurComplement object used by this PC

2063:     Not collective

2065:     Input Parameters:
2066: +   pc      - the preconditioner context
2067: .   S       - the Schur complement matrix

2069:     Level: advanced

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

2073: @*/
2074: PetscErrorCode  PCFieldSplitSchurRestoreS(PC pc,Mat *S)
2075: {
2077:   const char*    t;
2078:   PetscBool      isfs;
2079:   PC_FieldSplit  *jac;

2083:   PetscObjectGetType((PetscObject)pc,&t);
2084:   PetscStrcmp(t,PCFIELDSPLIT,&isfs);
2085:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
2086:   jac = (PC_FieldSplit*)pc->data;
2087:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
2088:   if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten");
2089:   return(0);
2090: }


2093: static PetscErrorCode  PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
2094: {
2095:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

2099:   jac->schurpre = ptype;
2100:   if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
2101:     MatDestroy(&jac->schur_user);
2102:     jac->schur_user = pre;
2103:     PetscObjectReference((PetscObject)jac->schur_user);
2104:   }
2105:   return(0);
2106: }

2108: static PetscErrorCode  PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
2109: {
2110:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

2113:   *ptype = jac->schurpre;
2114:   *pre   = jac->schur_user;
2115:   return(0);
2116: }

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

2121:     Collective on PC

2123:     Input Parameters:
2124: +   pc  - the preconditioner context
2125: -   ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default

2127:     Options Database:
2128: .     -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full


2131:     Level: intermediate

2133:     Notes:
2134:     The FULL factorization is

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

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

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

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

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

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

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

2158: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCFieldSplitSetSchurScale()
2159: @*/
2160: PetscErrorCode  PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype)
2161: {

2166:   PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));
2167:   return(0);
2168: }

2170: static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype)
2171: {
2172:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2175:   jac->schurfactorization = ftype;
2176:   return(0);
2177: }

2179: /*@
2180:     PCFieldSplitSetSchurScale -  Controls the sign flip of S for PC_FIELDSPLIT_SCHUR_FACT_DIAG.

2182:     Collective on PC

2184:     Input Parameters:
2185: +   pc    - the preconditioner context
2186: -   scale - scaling factor for the Schur complement

2188:     Options Database:
2189: .     -pc_fieldsplit_schur_scale - default is -1.0

2191:     Level: intermediate

2193: .seealso: PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurFactType, PCFieldSplitSetSchurScale()
2194: @*/
2195: PetscErrorCode PCFieldSplitSetSchurScale(PC pc,PetscScalar scale)
2196: {

2202:   PetscTryMethod(pc,"PCFieldSplitSetSchurScale_C",(PC,PetscScalar),(pc,scale));
2203:   return(0);
2204: }

2206: static PetscErrorCode PCFieldSplitSetSchurScale_FieldSplit(PC pc,PetscScalar scale)
2207: {
2208:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2211:   jac->schurscale = scale;
2212:   return(0);
2213: }

2215: /*@C
2216:    PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement

2218:    Collective on KSP

2220:    Input Parameter:
2221: .  pc - the preconditioner context

2223:    Output Parameters:
2224: +  A00 - the (0,0) block
2225: .  A01 - the (0,1) block
2226: .  A10 - the (1,0) block
2227: -  A11 - the (1,1) block

2229:    Level: advanced

2231: .seealso: PCFIELDSPLIT
2232: @*/
2233: PetscErrorCode  PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11)
2234: {
2235:   PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;

2239:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach.");
2240:   if (A00) *A00 = jac->pmat[0];
2241:   if (A01) *A01 = jac->B;
2242:   if (A10) *A10 = jac->C;
2243:   if (A11) *A11 = jac->pmat[1];
2244:   return(0);
2245: }

2247: static PetscErrorCode  PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type)
2248: {
2249:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

2253:   jac->type = type;
2254:   if (type == PC_COMPOSITE_SCHUR) {
2255:     pc->ops->apply = PCApply_FieldSplit_Schur;
2256:     pc->ops->view  = PCView_FieldSplit_Schur;

2258:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);
2259:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);
2260:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);
2261:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);
2262:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",PCFieldSplitSetSchurScale_FieldSplit);

2264:   } else {
2265:     pc->ops->apply = PCApply_FieldSplit;
2266:     pc->ops->view  = PCView_FieldSplit;

2268:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2269:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);
2270:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);
2271:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);
2272:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",0);
2273:   }
2274:   return(0);
2275: }

2277: static PetscErrorCode  PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs)
2278: {
2279:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2282:   if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs);
2283:   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);
2284:   jac->bs = bs;
2285:   return(0);
2286: }

2288: /*@
2289:    PCFieldSplitSetType - Sets the type of fieldsplit preconditioner.

2291:    Collective on PC

2293:    Input Parameter:
2294: .  pc - the preconditioner context
2295: .  type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR

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

2300:    Level: Intermediate

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

2304: .seealso: PCCompositeSetType()

2306: @*/
2307: PetscErrorCode  PCFieldSplitSetType(PC pc,PCCompositeType type)
2308: {

2313:   PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));
2314:   return(0);
2315: }

2317: /*@
2318:   PCFieldSplitGetType - Gets the type of fieldsplit preconditioner.

2320:   Not collective

2322:   Input Parameter:
2323: . pc - the preconditioner context

2325:   Output Parameter:
2326: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR

2328:   Level: Intermediate

2330: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
2331: .seealso: PCCompositeSetType()
2332: @*/
2333: PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
2334: {
2335:   PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;

2340:   *type = jac->type;
2341:   return(0);
2342: }

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

2347:    Logically Collective

2349:    Input Parameters:
2350: +  pc   - the preconditioner context
2351: -  flg  - boolean indicating whether to use field splits defined by the DM

2353:    Options Database Key:
2354: .  -pc_fieldsplit_dm_splits

2356:    Level: Intermediate

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

2360: .seealso: PCFieldSplitGetDMSplits()

2362: @*/
2363: PetscErrorCode  PCFieldSplitSetDMSplits(PC pc,PetscBool flg)
2364: {
2365:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
2366:   PetscBool      isfs;

2372:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2373:   if (isfs) {
2374:     jac->dm_splits = flg;
2375:   }
2376:   return(0);
2377: }


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

2383:    Logically Collective

2385:    Input Parameter:
2386: .  pc   - the preconditioner context

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

2391:    Level: Intermediate

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

2395: .seealso: PCFieldSplitSetDMSplits()

2397: @*/
2398: PetscErrorCode  PCFieldSplitGetDMSplits(PC pc,PetscBool* flg)
2399: {
2400:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
2401:   PetscBool      isfs;

2407:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2408:   if (isfs) {
2409:     if(flg) *flg = jac->dm_splits;
2410:   }
2411:   return(0);
2412: }

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

2417:    Logically Collective

2419:    Input Parameter:
2420: .  pc   - the preconditioner context

2422:    Output Parameter:
2423: .  flg  - boolean indicating whether to detect fields or not

2425:    Level: Intermediate

2427: .seealso: PCFIELDSPLIT, PCFieldSplitSetDetectSaddlePoint()

2429: @*/
2430: PetscErrorCode PCFieldSplitGetDetectSaddlePoint(PC pc,PetscBool *flg)
2431: {
2432:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2435:   *flg = jac->detect;
2436:   return(0);
2437: }

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

2442:    Logically Collective

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

2447:    Input Parameter:
2448: .  pc   - the preconditioner context

2450:    Output Parameter:
2451: .  flg  - boolean indicating whether to detect fields or not

2453:    Options Database Key:
2454: .  -pc_fieldsplit_detect_saddle_point

2456:    Level: Intermediate

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

2460: @*/
2461: PetscErrorCode PCFieldSplitSetDetectSaddlePoint(PC pc,PetscBool flg)
2462: {
2463:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

2467:   jac->detect = flg;
2468:   if (jac->detect) {
2469:     PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);
2470:     PCFieldSplitSetSchurPre(pc,PC_FIELDSPLIT_SCHUR_PRE_SELF,NULL);
2471:   }
2472:   return(0);
2473: }

2475: /* -------------------------------------------------------------------------------------*/
2476: /*MC
2477:    PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
2478:                   fields or groups of fields. See the users manual section "Solving Block Matrices" for more details.

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

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

2486:    Level: intermediate

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

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

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

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

2509: $     For the Schur complement preconditioner if J = ( A00 A01 )
2510: $                                                    ( A10 A11 )
2511: $     the preconditioner using full factorization is
2512: $              ( I   -ksp(A00) A01 ) ( inv(A00)     0  ) (     I          0  )
2513: $              ( 0         I       ) (   0      ksp(S) ) ( -A10 ksp(A00)  I  )
2514:      where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_.  S is the Schur complement
2515: $              S = A11 - A10 ksp(A00) A01
2516:      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
2517:      in providing the SECOND split or 1 if not give). For PCFieldSplitGetSubKSP() when field number is 0,
2518:      it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default
2519:      A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() for all the possible ways to construct the preconditioner for S.

2521:      The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above,
2522:      diag gives
2523: $              ( inv(A00)     0   )
2524: $              (   0      -ksp(S) )
2525:      note that slightly counter intuitively there is a negative in front of the ksp(S) so that the preconditioner is positive definite. For SPD matrices J, the sign flip
2526:      can be turned off with PCFieldSplitSetSchurScale() or by command line -pc_fieldsplit_schur_scale 1.0. The lower factorization is the inverse of
2527: $              (  A00   0 )
2528: $              (  A10   S )
2529:      where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of
2530: $              ( A00 A01 )
2531: $              (  0   S  )
2532:      where again the inverses of A00 and S are applied using KSPs.

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

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

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

2544:    Concepts: physics based preconditioners, block preconditioners

2546:    There is a nice discussion of block preconditioners in

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

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

2555: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC,
2556:            PCFieldSplitGetSubKSP(), PCFieldSplitSchurGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(),
2557:           MatSchurComplementSetAinvType(), PCFieldSplitSetSchurScale(),
2558:           PCFieldSplitSetDetectSaddlePoint()
2559: M*/

2561: PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
2562: {
2564:   PC_FieldSplit  *jac;

2567:   PetscNewLog(pc,&jac);

2569:   jac->bs                 = -1;
2570:   jac->nsplits            = 0;
2571:   jac->type               = PC_COMPOSITE_MULTIPLICATIVE;
2572:   jac->schurpre           = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
2573:   jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
2574:   jac->schurscale         = -1.0;
2575:   jac->dm_splits          = PETSC_TRUE;
2576:   jac->detect             = PETSC_FALSE;

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

2580:   pc->ops->apply           = PCApply_FieldSplit;
2581:   pc->ops->applytranspose  = PCApplyTranspose_FieldSplit;
2582:   pc->ops->setup           = PCSetUp_FieldSplit;
2583:   pc->ops->reset           = PCReset_FieldSplit;
2584:   pc->ops->destroy         = PCDestroy_FieldSplit;
2585:   pc->ops->setfromoptions  = PCSetFromOptions_FieldSplit;
2586:   pc->ops->view            = PCView_FieldSplit;
2587:   pc->ops->applyrichardson = 0;

2589:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",PCFieldSplitSchurGetSubKSP_FieldSplit);
2590:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2591:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);
2592:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);
2593:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);
2594:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);
2595:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",PCFieldSplitRestrictIS_FieldSplit);
2596:   return(0);
2597: }