Actual source code: fieldsplit.c

petsc-master 2016-02-06
Report Typos and Errors
  2: #include <petsc/private/pcimpl.h>     /*I "petscpc.h" I*/
  3: #include <petscdm.h>


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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

494:     jac->issetup = PETSC_TRUE;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

711:       /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
712:       MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);
713:       MatSetType(jac->schur,MATSCHURCOMPLEMENT);
714:       MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
715:       PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
716:       /* Note that the inner KSP is NOT going to inherit this prefix, and if it did, it would be reset just below.  Is that what we want? */
717:       MatSetOptionsPrefix(jac->schur,schurmatprefix);
718:       MatSetFromOptions(jac->schur);
719:       MatGetNullSpace(jac->mat[1], &sp);
720:       if (sp) {
721:         MatSetNullSpace(jac->schur, sp);
722:       }

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

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

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

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

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

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

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

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

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

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

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

904:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
905:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
906:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);

908:     if (kspUpper == kspA) {
909:       MatMult(jac->B,ilinkD->y,ilinkA->y);
910:       VecAXPY(ilinkA->x,-1.0,ilinkA->y);
911:       KSPSolve(kspA,ilinkA->x,ilinkA->y);
912:     } else {
913:       KSPSolve(kspA,ilinkA->x,ilinkA->y);
914:       MatMult(jac->B,ilinkD->y,ilinkA->x);
915:       KSPSolve(kspUpper,ilinkA->x,ilinkA->z);
916:       VecAXPY(ilinkA->y,-1.0,ilinkA->z);
917:     }
918:     VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
919:     VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
920:   }
921:   return(0);
922: }

926: static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y)
927: {
928:   PC_FieldSplit      *jac = (PC_FieldSplit*)pc->data;
929:   PetscErrorCode     ierr;
930:   PC_FieldSplitLink  ilink = jac->head;
931:   PetscInt           cnt,bs;
932:   KSPConvergedReason reason;

935:   if (jac->type == PC_COMPOSITE_ADDITIVE) {
936:     if (jac->defaultsplit) {
937:       VecGetBlockSize(x,&bs);
938:       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);
939:       VecGetBlockSize(y,&bs);
940:       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);
941:       VecStrideGatherAll(x,jac->x,INSERT_VALUES);
942:       while (ilink) {
943:         KSPSolve(ilink->ksp,ilink->x,ilink->y);
944:         KSPGetConvergedReason(ilink->ksp,&reason);
945:         if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
946:           pc->failedreason = PC_SUBPC_ERROR;
947:         }
948:         ilink = ilink->next;
949:       }
950:       VecStrideScatterAll(jac->y,y,INSERT_VALUES);
951:     } else {
952:       VecSet(y,0.0);
953:       while (ilink) {
954:         FieldSplitSplitSolveAdd(ilink,x,y);
955:         KSPGetConvergedReason(ilink->ksp,&reason);
956:         if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
957:           pc->failedreason = PC_SUBPC_ERROR;
958:         }
959:         ilink = ilink->next;
960:       }
961:     }
962:   } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) {
963:     VecSet(y,0.0);
964:     /* solve on first block for first block variables */
965:     VecScatterBegin(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
966:     VecScatterEnd(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);
967:     KSPSolve(ilink->ksp,ilink->x,ilink->y);
968:     KSPGetConvergedReason(ilink->ksp,&reason);
969:     if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
970:       pc->failedreason = PC_SUBPC_ERROR;
971:     }
972:     VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
973:     VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);

975:     /* compute the residual only onto second block variables using first block variables */
976:     MatMult(jac->Afield[1],ilink->y,ilink->next->x);
977:     ilink = ilink->next;
978:     VecScale(ilink->x,-1.0);
979:     VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
980:     VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);

982:     /* solve on second block variables */
983:     KSPSolve(ilink->ksp,ilink->x,ilink->y);
984:     KSPGetConvergedReason(ilink->ksp,&reason);
985:     if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
986:       pc->failedreason = PC_SUBPC_ERROR;
987:     }
988:     VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
989:     VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
990:   } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
991:     if (!jac->w1) {
992:       VecDuplicate(x,&jac->w1);
993:       VecDuplicate(x,&jac->w2);
994:     }
995:     VecSet(y,0.0);
996:     FieldSplitSplitSolveAdd(ilink,x,y);
997:     KSPGetConvergedReason(ilink->ksp,&reason);
998:     if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
999:       pc->failedreason = PC_SUBPC_ERROR;
1000:     }
1001:     cnt  = 1;
1002:     while (ilink->next) {
1003:       ilink = ilink->next;
1004:       /* compute the residual only over the part of the vector needed */
1005:       MatMult(jac->Afield[cnt++],y,ilink->x);
1006:       VecScale(ilink->x,-1.0);
1007:       VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1008:       VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1009:       KSPSolve(ilink->ksp,ilink->x,ilink->y);
1010:       KSPGetConvergedReason(ilink->ksp,&reason);
1011:       if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1012:         pc->failedreason = PC_SUBPC_ERROR;
1013:       }
1014:       VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1015:       VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1016:     }
1017:     if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1018:       cnt -= 2;
1019:       while (ilink->previous) {
1020:         ilink = ilink->previous;
1021:         /* compute the residual only over the part of the vector needed */
1022:         MatMult(jac->Afield[cnt--],y,ilink->x);
1023:         VecScale(ilink->x,-1.0);
1024:         VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1025:         VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
1026:         KSPSolve(ilink->ksp,ilink->x,ilink->y);
1027:         KSPGetConvergedReason(ilink->ksp,&reason);
1028:         if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1029:           pc->failedreason = PC_SUBPC_ERROR;
1030:         }
1031:         VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1032:         VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
1033:       }
1034:     }
1035:   } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type);
1036:   return(0);
1037: }

1039: #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \
1040:   (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1041:    VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
1042:    KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \
1043:    VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \
1044:    VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE))

1048: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y)
1049: {
1050:   PC_FieldSplit      *jac = (PC_FieldSplit*)pc->data;
1051:   PetscErrorCode     ierr;
1052:   PC_FieldSplitLink  ilink = jac->head;
1053:   PetscInt           bs;
1054:   KSPConvergedReason reason;

1057:   if (jac->type == PC_COMPOSITE_ADDITIVE) {
1058:     if (jac->defaultsplit) {
1059:       VecGetBlockSize(x,&bs);
1060:       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);
1061:       VecGetBlockSize(y,&bs);
1062:       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);
1063:       VecStrideGatherAll(x,jac->x,INSERT_VALUES);
1064:       while (ilink) {
1065:         KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);
1066:         KSPGetConvergedReason(ilink->ksp,&reason);
1067:         if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1068:           pc->failedreason = PC_SUBPC_ERROR;
1069:         }
1070:         ilink = ilink->next;
1071:       }
1072:       VecStrideScatterAll(jac->y,y,INSERT_VALUES);
1073:     } else {
1074:       VecSet(y,0.0);
1075:       while (ilink) {
1076:         FieldSplitSplitSolveAddTranspose(ilink,x,y);
1077:         KSPGetConvergedReason(ilink->ksp,&reason);
1078:         if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1079:           pc->failedreason = PC_SUBPC_ERROR;
1080:         }
1081:         ilink = ilink->next;
1082:       }
1083:     }
1084:   } else {
1085:     if (!jac->w1) {
1086:       VecDuplicate(x,&jac->w1);
1087:       VecDuplicate(x,&jac->w2);
1088:     }
1089:     VecSet(y,0.0);
1090:     if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1091:       FieldSplitSplitSolveAddTranspose(ilink,x,y);
1092:       KSPGetConvergedReason(ilink->ksp,&reason);
1093:       if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1094:         pc->failedreason = PC_SUBPC_ERROR;
1095:       }
1096:       while (ilink->next) {
1097:         ilink = ilink->next;
1098:         MatMultTranspose(pc->mat,y,jac->w1);
1099:         VecWAXPY(jac->w2,-1.0,jac->w1,x);
1100:         FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1101:       }
1102:       while (ilink->previous) {
1103:         ilink = ilink->previous;
1104:         MatMultTranspose(pc->mat,y,jac->w1);
1105:         VecWAXPY(jac->w2,-1.0,jac->w1,x);
1106:         FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1107:       }
1108:     } else {
1109:       while (ilink->next) {   /* get to last entry in linked list */
1110:         ilink = ilink->next;
1111:       }
1112:       FieldSplitSplitSolveAddTranspose(ilink,x,y);
1113:       KSPGetConvergedReason(ilink->ksp,&reason);
1114:       if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
1115:         pc->failedreason = PC_SUBPC_ERROR;
1116:       }
1117:       while (ilink->previous) {
1118:         ilink = ilink->previous;
1119:         MatMultTranspose(pc->mat,y,jac->w1);
1120:         VecWAXPY(jac->w2,-1.0,jac->w1,x);
1121:         FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1122:       }
1123:     }
1124:   }
1125:   return(0);
1126: }

1130: static PetscErrorCode PCReset_FieldSplit(PC pc)
1131: {
1132:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1133:   PetscErrorCode    ierr;
1134:   PC_FieldSplitLink ilink = jac->head,next;

1137:   while (ilink) {
1138:     KSPReset(ilink->ksp);
1139:     VecDestroy(&ilink->x);
1140:     VecDestroy(&ilink->y);
1141:     VecDestroy(&ilink->z);
1142:     VecScatterDestroy(&ilink->sctx);
1143:     ISDestroy(&ilink->is);
1144:     ISDestroy(&ilink->is_col);
1145:     next  = ilink->next;
1146:     ilink = next;
1147:   }
1148:   PetscFree2(jac->x,jac->y);
1149:   if (jac->mat && jac->mat != jac->pmat) {
1150:     MatDestroyMatrices(jac->nsplits,&jac->mat);
1151:   } else if (jac->mat) {
1152:     jac->mat = NULL;
1153:   }
1154:   if (jac->pmat) {MatDestroyMatrices(jac->nsplits,&jac->pmat);}
1155:   if (jac->Afield) {MatDestroyMatrices(jac->nsplits,&jac->Afield);}
1156:   VecDestroy(&jac->w1);
1157:   VecDestroy(&jac->w2);
1158:   MatDestroy(&jac->schur);
1159:   MatDestroy(&jac->schurp);
1160:   MatDestroy(&jac->schur_user);
1161:   KSPDestroy(&jac->kspschur);
1162:   KSPDestroy(&jac->kspupper);
1163:   MatDestroy(&jac->B);
1164:   MatDestroy(&jac->C);
1165:   jac->reset = PETSC_TRUE;
1166:   return(0);
1167: }

1171: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
1172: {
1173:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1174:   PetscErrorCode    ierr;
1175:   PC_FieldSplitLink ilink = jac->head,next;

1178:   PCReset_FieldSplit(pc);
1179:   while (ilink) {
1180:     KSPDestroy(&ilink->ksp);
1181:     next  = ilink->next;
1182:     PetscFree(ilink->splitname);
1183:     PetscFree(ilink->fields);
1184:     PetscFree(ilink->fields_col);
1185:     PetscFree(ilink);
1186:     ilink = next;
1187:   }
1188:   PetscFree2(jac->x,jac->y);
1189:   PetscFree(pc->data);
1190:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);
1191:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);
1192:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);
1193:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);
1194:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);
1195:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);
1196:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);
1197:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);
1198:   return(0);
1199: }

1203: static PetscErrorCode PCSetFromOptions_FieldSplit(PetscOptionItems *PetscOptionsObject,PC pc)
1204: {
1205:   PetscErrorCode  ierr;
1206:   PetscInt        bs;
1207:   PetscBool       flg,stokes = PETSC_FALSE;
1208:   PC_FieldSplit   *jac = (PC_FieldSplit*)pc->data;
1209:   PCCompositeType ctype;

1212:   PetscOptionsHead(PetscOptionsObject,"FieldSplit options");
1213:   PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);
1214:   PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);
1215:   if (flg) {
1216:     PCFieldSplitSetBlockSize(pc,bs);
1217:   }
1218:   jac->diag_use_amat = pc->useAmat;
1219:   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);
1220:   jac->offdiag_use_amat = pc->useAmat;
1221:   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);
1222:   /* FIXME: No programmatic equivalent to the following. */
1223:   PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);
1224:   if (stokes) {
1225:     PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);
1226:     jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF;
1227:   }

1229:   PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);
1230:   if (flg) {
1231:     PCFieldSplitSetType(pc,ctype);
1232:   }
1233:   /* Only setup fields once */
1234:   if ((jac->bs > 0) && (jac->nsplits == 0)) {
1235:     /* only allow user to set fields from command line if bs is already known.
1236:        otherwise user can set them in PCFieldSplitSetDefaults() */
1237:     PCFieldSplitSetRuntimeSplits_Private(pc);
1238:     if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
1239:   }
1240:   if (jac->type == PC_COMPOSITE_SCHUR) {
1241:     PetscOptionsGetEnum(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);
1242:     if (flg) {PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");}
1243:     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);
1244:     PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);
1245:   }
1246:   PetscOptionsTail();
1247:   return(0);
1248: }

1250: /*------------------------------------------------------------------------------------*/

1254: static PetscErrorCode  PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1255: {
1256:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1257:   PetscErrorCode    ierr;
1258:   PC_FieldSplitLink ilink,next = jac->head;
1259:   char              prefix[128];
1260:   PetscInt          i;

1263:   if (jac->splitdefined) {
1264:     PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1265:     return(0);
1266:   }
1267:   for (i=0; i<n; i++) {
1268:     if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs);
1269:     if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]);
1270:   }
1271:   PetscNew(&ilink);
1272:   if (splitname) {
1273:     PetscStrallocpy(splitname,&ilink->splitname);
1274:   } else {
1275:     PetscMalloc1(3,&ilink->splitname);
1276:     PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);
1277:   }
1278:   PetscMalloc1(n,&ilink->fields);
1279:   PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));
1280:   PetscMalloc1(n,&ilink->fields_col);
1281:   PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));

1283:   ilink->nfields = n;
1284:   ilink->next    = NULL;
1285:   KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1286:   KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1287:   PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1288:   KSPSetType(ilink->ksp,KSPPREONLY);
1289:   PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);

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

1294:   if (!next) {
1295:     jac->head       = ilink;
1296:     ilink->previous = NULL;
1297:   } else {
1298:     while (next->next) {
1299:       next = next->next;
1300:     }
1301:     next->next      = ilink;
1302:     ilink->previous = next;
1303:   }
1304:   jac->nsplits++;
1305:   return(0);
1306: }

1310: static PetscErrorCode  PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp)
1311: {
1312:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

1316:   PetscMalloc1(jac->nsplits,subksp);
1317:   MatSchurComplementGetKSP(jac->schur,*subksp);

1319:   (*subksp)[1] = jac->kspschur;
1320:   if (n) *n = jac->nsplits;
1321:   return(0);
1322: }

1326: static PetscErrorCode  PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1327: {
1328:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1329:   PetscErrorCode    ierr;
1330:   PetscInt          cnt   = 0;
1331:   PC_FieldSplitLink ilink = jac->head;

1334:   PetscMalloc1(jac->nsplits,subksp);
1335:   while (ilink) {
1336:     (*subksp)[cnt++] = ilink->ksp;
1337:     ilink            = ilink->next;
1338:   }
1339:   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);
1340:   if (n) *n = jac->nsplits;
1341:   return(0);
1342: }

1346: static PetscErrorCode  PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is)
1347: {
1348:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1349:   PetscErrorCode    ierr;
1350:   PC_FieldSplitLink ilink, next = jac->head;
1351:   char              prefix[128];

1354:   if (jac->splitdefined) {
1355:     PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1356:     return(0);
1357:   }
1358:   PetscNew(&ilink);
1359:   if (splitname) {
1360:     PetscStrallocpy(splitname,&ilink->splitname);
1361:   } else {
1362:     PetscMalloc1(8,&ilink->splitname);
1363:     PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);
1364:   }
1365:   PetscObjectReference((PetscObject)is);
1366:   ISDestroy(&ilink->is);
1367:   ilink->is     = is;
1368:   PetscObjectReference((PetscObject)is);
1369:   ISDestroy(&ilink->is_col);
1370:   ilink->is_col = is;
1371:   ilink->next   = NULL;
1372:   KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1373:   KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);
1374:   PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1375:   KSPSetType(ilink->ksp,KSPPREONLY);
1376:   PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);

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

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

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

1400:     Logically Collective on PC

1402:     Input Parameters:
1403: +   pc  - the preconditioner context
1404: .   splitname - name of this split, if NULL the number of the split is used
1405: .   n - the number of fields in this split
1406: -   fields - the fields in this split

1408:     Level: intermediate

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

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

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

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

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

1426: @*/
1427: PetscErrorCode  PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1428: {

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

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

1445:     Logically Collective on PC

1447:     Input Parameters:
1448: +   pc  - the preconditioner object
1449: -   flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from

1451:     Options Database:
1452: .     -pc_fieldsplit_diag_use_amat

1454:     Level: intermediate

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

1458: @*/
1459: PetscErrorCode  PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg)
1460: {
1461:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1462:   PetscBool      isfs;

1467:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1468:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1469:   jac->diag_use_amat = flg;
1470:   return(0);
1471: }

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

1478:     Logically Collective on PC

1480:     Input Parameters:
1481: .   pc  - the preconditioner object

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


1487:     Level: intermediate

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

1491: @*/
1492: PetscErrorCode  PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg)
1493: {
1494:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1495:   PetscBool      isfs;

1501:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1502:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1503:   *flg = jac->diag_use_amat;
1504:   return(0);
1505: }

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

1512:     Logically Collective on PC

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

1518:     Options Database:
1519: .     -pc_fieldsplit_off_diag_use_amat

1521:     Level: intermediate

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

1525: @*/
1526: PetscErrorCode  PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg)
1527: {
1528:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1529:   PetscBool      isfs;

1534:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1535:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1536:   jac->offdiag_use_amat = flg;
1537:   return(0);
1538: }

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

1545:     Logically Collective on PC

1547:     Input Parameters:
1548: .   pc  - the preconditioner object

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


1554:     Level: intermediate

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

1558: @*/
1559: PetscErrorCode  PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg)
1560: {
1561:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1562:   PetscBool      isfs;

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



1578: /*@C
1579:     PCFieldSplitSetIS - Sets the exact elements for field

1581:     Logically Collective on PC

1583:     Input Parameters:
1584: +   pc  - the preconditioner context
1585: .   splitname - name of this split, if NULL the number of the split is used
1586: -   is - the index set that defines the vector elements in this field


1589:     Notes:
1590:     Use PCFieldSplitSetFields(), for fields defined by strided types.

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

1595:     Level: intermediate

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

1599: @*/
1600: PetscErrorCode  PCFieldSplitSetIS(PC pc,const char splitname[],IS is)
1601: {

1608:   PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));
1609:   return(0);
1610: }

1614: /*@
1615:     PCFieldSplitGetIS - Retrieves the elements for a field as an IS

1617:     Logically Collective on PC

1619:     Input Parameters:
1620: +   pc  - the preconditioner context
1621: -   splitname - name of this split

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

1626:     Level: intermediate

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

1630: @*/
1631: PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is)
1632: {

1639:   {
1640:     PC_FieldSplit     *jac  = (PC_FieldSplit*) pc->data;
1641:     PC_FieldSplitLink ilink = jac->head;
1642:     PetscBool         found;

1644:     *is = NULL;
1645:     while (ilink) {
1646:       PetscStrcmp(ilink->splitname, splitname, &found);
1647:       if (found) {
1648:         *is = ilink->is;
1649:         break;
1650:       }
1651:       ilink = ilink->next;
1652:     }
1653:   }
1654:   return(0);
1655: }

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

1663:     Logically Collective on PC

1665:     Input Parameters:
1666: +   pc  - the preconditioner context
1667: -   bs - the block size

1669:     Level: intermediate

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

1673: @*/
1674: PetscErrorCode  PCFieldSplitSetBlockSize(PC pc,PetscInt bs)
1675: {

1681:   PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));
1682:   return(0);
1683: }

1687: /*@C
1688:    PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits

1690:    Collective on KSP

1692:    Input Parameter:
1693: .  pc - the preconditioner context

1695:    Output Parameters:
1696: +  n - the number of splits
1697: -  pc - the array of KSP contexts

1699:    Note:
1700:    After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user
1701:    (not the KSP just the array that contains them).

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

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


1710:    Level: advanced

1712: .seealso: PCFIELDSPLIT
1713: @*/
1714: PetscErrorCode  PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
1715: {

1721:   PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
1722:   return(0);
1723: }

1727: /*@
1728:     PCFieldSplitSetSchurPre -  Indicates if the Schur complement is preconditioned by a preconditioner constructed by the
1729:       A11 matrix. Otherwise no preconditioner is used.

1731:     Collective on PC

1733:     Input Parameters:
1734: +   pc      - the preconditioner context
1735: .   ptype   - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11 (default), PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER
1736: -   userpre - matrix to use for preconditioning, or NULL

1738:     Options Database:
1739: .     -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> default is a11

1741:     Notes:
1742: $    If ptype is
1743: $        user then the preconditioner for the Schur complement is generated from the user provided matrix (pre argument
1744: $             to this function).
1745: $        a11 then the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner
1746: $             matrix associated with the Schur complement (i.e. A11), not he Schur complement matrix
1747: $        full then the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation computed internally by PFIELDSPLIT (this is expensive)
1748: $        self the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix:
1749: $             The only preconditioner that currently works with this symbolic respresentation matrix object is the PCLSC
1750: $             preconditioner
1751: $        selfp then the preconditioning for the Schur complement is generated from an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01
1752: $             This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be
1753: $             lumped before extracting the diagonal: -fieldsplit_1_mat_schur_complement_ainv_type lump

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

1759:     Level: intermediate

1761: .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType,
1762:           MatSchurComplementSetAinvType(), PCLSC

1764: @*/
1765: PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
1766: {

1771:   PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));
1772:   return(0);
1773: }
1774: PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */

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

1782:     Logically Collective on PC

1784:     Input Parameters:
1785: .   pc      - the preconditioner context

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

1791:     Level: intermediate

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

1795: @*/
1796: PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
1797: {

1802:   PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));
1803:   return(0);
1804: }

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

1811:     Not collective

1813:     Input Parameter:
1814: .   pc      - the preconditioner context

1816:     Output Parameter:
1817: .   S       - the Schur complement matrix

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

1822:     Level: advanced

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

1826: @*/
1827: PetscErrorCode  PCFieldSplitSchurGetS(PC pc,Mat *S)
1828: {
1830:   const char*    t;
1831:   PetscBool      isfs;
1832:   PC_FieldSplit  *jac;

1836:   PetscObjectGetType((PetscObject)pc,&t);
1837:   PetscStrcmp(t,PCFIELDSPLIT,&isfs);
1838:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
1839:   jac = (PC_FieldSplit*)pc->data;
1840:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
1841:   if (S) *S = jac->schur;
1842:   return(0);
1843: }

1847: /*@
1848:     PCFieldSplitSchurRestoreS -  restores the MatSchurComplement object used by this PC

1850:     Not collective

1852:     Input Parameters:
1853: +   pc      - the preconditioner context
1854: .   S       - the Schur complement matrix

1856:     Level: advanced

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

1860: @*/
1861: PetscErrorCode  PCFieldSplitSchurRestoreS(PC pc,Mat *S)
1862: {
1864:   const char*    t;
1865:   PetscBool      isfs;
1866:   PC_FieldSplit  *jac;

1870:   PetscObjectGetType((PetscObject)pc,&t);
1871:   PetscStrcmp(t,PCFIELDSPLIT,&isfs);
1872:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
1873:   jac = (PC_FieldSplit*)pc->data;
1874:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
1875:   if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten");
1876:   return(0);
1877: }


1882: static PetscErrorCode  PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
1883: {
1884:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

1888:   jac->schurpre = ptype;
1889:   if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
1890:     MatDestroy(&jac->schur_user);
1891:     jac->schur_user = pre;
1892:     PetscObjectReference((PetscObject)jac->schur_user);
1893:   }
1894:   return(0);
1895: }

1899: static PetscErrorCode  PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
1900: {
1901:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

1904:   *ptype = jac->schurpre;
1905:   *pre   = jac->schur_user;
1906:   return(0);
1907: }

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

1914:     Collective on PC

1916:     Input Parameters:
1917: +   pc  - the preconditioner context
1918: -   ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default

1920:     Options Database:
1921: .     -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full


1924:     Level: intermediate

1926:     Notes:
1927:     The FULL factorization is

1929: $   (A   B)  = (1       0) (A   0) (1  Ainv*B)
1930: $   (C   D)    (C*Ainv  1) (0   S) (0     1  )

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

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

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

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

1947: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType
1948: @*/
1949: PetscErrorCode  PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype)
1950: {

1955:   PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));
1956:   return(0);
1957: }

1961: static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype)
1962: {
1963:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

1966:   jac->schurfactorization = ftype;
1967:   return(0);
1968: }

1972: /*@C
1973:    PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement

1975:    Collective on KSP

1977:    Input Parameter:
1978: .  pc - the preconditioner context

1980:    Output Parameters:
1981: +  A00 - the (0,0) block
1982: .  A01 - the (0,1) block
1983: .  A10 - the (1,0) block
1984: -  A11 - the (1,1) block

1986:    Level: advanced

1988: .seealso: PCFIELDSPLIT
1989: @*/
1990: PetscErrorCode  PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11)
1991: {
1992:   PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;

1996:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach.");
1997:   if (A00) *A00 = jac->pmat[0];
1998:   if (A01) *A01 = jac->B;
1999:   if (A10) *A10 = jac->C;
2000:   if (A11) *A11 = jac->pmat[1];
2001:   return(0);
2002: }

2006: static PetscErrorCode  PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type)
2007: {
2008:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

2012:   jac->type = type;
2013:   if (type == PC_COMPOSITE_SCHUR) {
2014:     pc->ops->apply = PCApply_FieldSplit_Schur;
2015:     pc->ops->view  = PCView_FieldSplit_Schur;

2017:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);
2018:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);
2019:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);
2020:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);

2022:   } else {
2023:     pc->ops->apply = PCApply_FieldSplit;
2024:     pc->ops->view  = PCView_FieldSplit;

2026:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2027:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);
2028:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);
2029:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);
2030:   }
2031:   return(0);
2032: }

2036: static PetscErrorCode  PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs)
2037: {
2038:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

2041:   if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs);
2042:   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);
2043:   jac->bs = bs;
2044:   return(0);
2045: }

2049: /*@
2050:    PCFieldSplitSetType - Sets the type of fieldsplit preconditioner.

2052:    Collective on PC

2054:    Input Parameter:
2055: .  pc - the preconditioner context
2056: .  type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR

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

2061:    Level: Intermediate

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

2065: .seealso: PCCompositeSetType()

2067: @*/
2068: PetscErrorCode  PCFieldSplitSetType(PC pc,PCCompositeType type)
2069: {

2074:   PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));
2075:   return(0);
2076: }

2080: /*@
2081:   PCFieldSplitGetType - Gets the type of fieldsplit preconditioner.

2083:   Not collective

2085:   Input Parameter:
2086: . pc - the preconditioner context

2088:   Output Parameter:
2089: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR

2091:   Level: Intermediate

2093: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
2094: .seealso: PCCompositeSetType()
2095: @*/
2096: PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
2097: {
2098:   PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;

2103:   *type = jac->type;
2104:   return(0);
2105: }

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

2112:    Logically Collective

2114:    Input Parameters:
2115: +  pc   - the preconditioner context
2116: -  flg  - boolean indicating whether to use field splits defined by the DM

2118:    Options Database Key:
2119: .  -pc_fieldsplit_dm_splits

2121:    Level: Intermediate

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

2125: .seealso: PCFieldSplitGetDMSplits()

2127: @*/
2128: PetscErrorCode  PCFieldSplitSetDMSplits(PC pc,PetscBool flg)
2129: {
2130:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
2131:   PetscBool      isfs;

2137:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2138:   if (isfs) {
2139:     jac->dm_splits = flg;
2140:   }
2141:   return(0);
2142: }


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

2150:    Logically Collective

2152:    Input Parameter:
2153: .  pc   - the preconditioner context

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

2158:    Level: Intermediate

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

2162: .seealso: PCFieldSplitSetDMSplits()

2164: @*/
2165: PetscErrorCode  PCFieldSplitGetDMSplits(PC pc,PetscBool* flg)
2166: {
2167:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
2168:   PetscBool      isfs;

2174:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2175:   if (isfs) {
2176:     if(flg) *flg = jac->dm_splits;
2177:   }
2178:   return(0);
2179: }

2181: /* -------------------------------------------------------------------------------------*/
2182: /*MC
2183:    PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
2184:                   fields or groups of fields. See the users manual section "Solving Block Matrices" for more details.

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

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

2192:    Level: intermediate

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

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

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

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

2215: $     For the Schur complement preconditioner if J = ( A00 A01 )
2216: $                                                    ( A10 A11 )
2217: $     the preconditioner using full factorization is
2218: $              ( I   -ksp(A00) A01 ) ( inv(A00)     0  ) (     I          0  )
2219: $              ( 0         I       ) (   0      ksp(S) ) ( -A10 ksp(A00)  I  )
2220:      where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_.  S is the Schur complement
2221: $              S = A11 - A10 ksp(A00) A01
2222:      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
2223:      in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0,
2224:      it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default
2225:      A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() for all the possible ways to construct the preconditioner for S.

2227:      The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above,
2228:      diag gives
2229: $              ( inv(A00)     0   )
2230: $              (   0      -ksp(S) )
2231:      note that slightly counter intuitively there is a negative in front of the ksp(S) so that the preconditioner is positive definite. The lower factorization is the inverse of
2232: $              (  A00   0 )
2233: $              (  A10   S )
2234:      where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of
2235: $              ( A00 A01 )
2236: $              (  0   S  )
2237:      where again the inverses of A00 and S are applied using KSPs.

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

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

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

2249:    Concepts: physics based preconditioners, block preconditioners

2251:    There is a nice discussion of block preconditioners in

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

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

2261: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC,
2262:            PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(),
2263:            MatSchurComplementSetAinvType()
2264: M*/

2268: PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
2269: {
2271:   PC_FieldSplit  *jac;

2274:   PetscNewLog(pc,&jac);

2276:   jac->bs                 = -1;
2277:   jac->nsplits            = 0;
2278:   jac->type               = PC_COMPOSITE_MULTIPLICATIVE;
2279:   jac->schurpre           = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
2280:   jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
2281:   jac->dm_splits          = PETSC_TRUE;

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

2285:   pc->ops->apply           = PCApply_FieldSplit;
2286:   pc->ops->applytranspose  = PCApplyTranspose_FieldSplit;
2287:   pc->ops->setup           = PCSetUp_FieldSplit;
2288:   pc->ops->reset           = PCReset_FieldSplit;
2289:   pc->ops->destroy         = PCDestroy_FieldSplit;
2290:   pc->ops->setfromoptions  = PCSetFromOptions_FieldSplit;
2291:   pc->ops->view            = PCView_FieldSplit;
2292:   pc->ops->applyrichardson = 0;

2294:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2295:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);
2296:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);
2297:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);
2298:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);
2299:   return(0);
2300: }