Actual source code: fieldsplit.c

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

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


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

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

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

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

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

493:     jac->issetup = PETSC_TRUE;

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

531:   ilink = jac->head;
532:   if (!jac->pmat) {
533:     Vec xtmp;

535:     MatGetVecs(pc->pmat,&xtmp,NULL);
536:     PetscMalloc1(nsplit,&jac->pmat);
537:     PetscMalloc2(nsplit,&jac->x,nsplit,&jac->y);
538:     for (i=0; i<nsplit; i++) {
539:       MatNullSpace sp;

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

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

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

604:   if (jac->type != PC_COMPOSITE_ADDITIVE  && jac->type != PC_COMPOSITE_SCHUR) {
605:     /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */
606:     /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */
607:     ilink = jac->head;
608:     if (!jac->Afield) {
609:       PetscMalloc1(nsplit,&jac->Afield);
610:       for (i=0; i<nsplit; i++) {
611:         MatGetSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
612:         ilink = ilink->next;
613:       }
614:     } else {
615:       for (i=0; i<nsplit; i++) {
616:         MatGetSubMatrix(pc->mat,ilink->is,NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);
617:         ilink = ilink->next;
618:       }
619:     }
620:   }

622:   if (jac->type == PC_COMPOSITE_SCHUR) {
623:     IS          ccis;
624:     PetscInt    rstart,rend;
625:     char        lscname[256];
626:     PetscObject LSC_L;

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

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

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

637:       MatSchurComplementGetKSP(jac->schur, &kspInner);
638:       ilink = jac->head;
639:       ISComplement(ilink->is_col,rstart,rend,&ccis);
640:       if (jac->offdiag_use_amat) {
641:         MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);
642:       } else {
643:         MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);
644:       }
645:       ISDestroy(&ccis);
646:       ilink = ilink->next;
647:       ISComplement(ilink->is_col,rstart,rend,&ccis);
648:       if (jac->offdiag_use_amat) {
649:         MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);
650:       } else {
651:         MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);
652:       }
653:       ISDestroy(&ccis);
654:       MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
655:       if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
656:         MatDestroy(&jac->schurp);
657:         MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);
658:       }
659:       if (kspA != kspInner) {
660:         KSPSetOperators(kspA,jac->mat[0],jac->pmat[0]);
661:       }
662:       if (kspUpper != kspA) {
663:         KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0]);
664:       }
665:       KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
666:     } else {
667:       const char   *Dprefix;
668:       char         schurprefix[256], schurmatprefix[256];
669:       char         schurtestoption[256];
670:       MatNullSpace sp;
671:       PetscBool    flg;

673:       /* extract the A01 and A10 matrices */
674:       ilink = jac->head;
675:       ISComplement(ilink->is_col,rstart,rend,&ccis);
676:       if (jac->offdiag_use_amat) {
677:         MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
678:       } else {
679:         MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
680:       }
681:       ISDestroy(&ccis);
682:       ilink = ilink->next;
683:       ISComplement(ilink->is_col,rstart,rend,&ccis);
684:       if (jac->offdiag_use_amat) {
685:         MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
686:       } else {
687:         MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
688:       }
689:       ISDestroy(&ccis);

691:       /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
692:       MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);
693:       MatSetType(jac->schur,MATSCHURCOMPLEMENT);
694:       MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
695:       PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
696:       /* 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? */
697:       MatSetOptionsPrefix(jac->schur,schurmatprefix);
698:       MatSetFromOptions(jac->schur);
699:       MatGetNullSpace(jac->pmat[1], &sp);
700:       if (sp) {
701:         MatSetNullSpace(jac->schur, sp);
702:       }

704:       PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);
705:       PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
706:       if (flg) {
707:         DM  dmInner;
708:         KSP kspInner;

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

716:         /* Set DM for new solver */
717:         KSPGetDM(jac->head->ksp, &dmInner);
718:         KSPSetDM(kspInner, dmInner);
719:         KSPSetDMActive(kspInner, PETSC_FALSE);
720:       } else {
721:          /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or
722:           * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact,
723:           * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for
724:           * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make
725:           * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used
726:           * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */
727:         KSPSetType(jac->head->ksp,KSPGMRES);
728:         MatSchurComplementSetKSP(jac->schur,jac->head->ksp);
729:       }
730:       KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0]);
731:       KSPSetFromOptions(jac->head->ksp);
732:       MatSetFromOptions(jac->schur);

734:       PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);
735:       PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
736:       if (flg) {
737:         DM dmInner;

739:         PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
740:         KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);
741:         KSPSetOptionsPrefix(jac->kspupper, schurprefix);
742:         KSPGetDM(jac->head->ksp, &dmInner);
743:         KSPSetDM(jac->kspupper, dmInner);
744:         KSPSetDMActive(jac->kspupper, PETSC_FALSE);
745:         KSPSetFromOptions(jac->kspupper);
746:         KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0]);
747:         VecDuplicate(jac->head->x, &jac->head->z);
748:       } else {
749:         jac->kspupper = jac->head->ksp;
750:         PetscObjectReference((PetscObject) jac->head->ksp);
751:       }

753:       if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
754:         MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);
755:       }
756:       KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);
757:       PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);
758:       PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);
759:       if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
760:         PC pcschur;
761:         KSPGetPC(jac->kspschur,&pcschur);
762:         PCSetType(pcschur,PCNONE);
763:         /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
764:       } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) {
765:         MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user);
766:       }
767:       KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
768:       KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);
769:       KSPSetOptionsPrefix(jac->kspschur,         Dprefix);
770:       /* propogate DM */
771:       {
772:         DM sdm;
773:         KSPGetDM(jac->head->next->ksp, &sdm);
774:         if (sdm) {
775:           KSPSetDM(jac->kspschur, sdm);
776:           KSPSetDMActive(jac->kspschur, PETSC_FALSE);
777:         }
778:       }
779:       /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
780:       /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */
781:       KSPSetFromOptions(jac->kspschur);
782:     }

784:     /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */
785:     PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);
786:     PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);
787:     if (!LSC_L) {PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);}
788:     if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);}
789:     PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);
790:     PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);
791:     if (!LSC_L) {PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);}
792:     if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);}
793:   } else {
794:     /* set up the individual splits' PCs */
795:     i     = 0;
796:     ilink = jac->head;
797:     while (ilink) {
798:       KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i]);
799:       /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
800:       if (!jac->suboptionsset) {KSPSetFromOptions(ilink->ksp);}
801:       i++;
802:       ilink = ilink->next;
803:     }
804:   }

806:   jac->suboptionsset = PETSC_TRUE;
807:   return(0);
808: }

810: #define FieldSplitSplitSolveAdd(ilink,xx,yy) \
811:   (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
812:    VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
813:    KSPSolve(ilink->ksp,ilink->x,ilink->y) || \
814:    VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \
815:    VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE))

819: static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y)
820: {
821:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
822:   PetscErrorCode    ierr;
823:   PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
824:   KSP               kspA   = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;

827:   switch (jac->schurfactorization) {
828:   case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
829:     /* [A00 0; 0 -S], positive definite, suitable for MINRES */
830:     VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
831:     VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
832:     VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
833:     KSPSolve(kspA,ilinkA->x,ilinkA->y);
834:     VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
835:     VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
836:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
837:     VecScale(ilinkD->y,-1.);
838:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
839:     VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
840:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
841:     break;
842:   case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
843:     /* [A00 0; A10 S], suitable for left preconditioning */
844:     VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
845:     VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
846:     KSPSolve(kspA,ilinkA->x,ilinkA->y);
847:     MatMult(jac->C,ilinkA->y,ilinkD->x);
848:     VecScale(ilinkD->x,-1.);
849:     VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
850:     VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
851:     VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
852:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
853:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
854:     VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
855:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
856:     break;
857:   case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
858:     /* [A00 A01; 0 S], suitable for right preconditioning */
859:     VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
860:     VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
861:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
862:     MatMult(jac->B,ilinkD->y,ilinkA->x);
863:     VecScale(ilinkA->x,-1.);
864:     VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
865:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
866:     VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
867:     KSPSolve(kspA,ilinkA->x,ilinkA->y);
868:     VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
869:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
870:     VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
871:     break;
872:   case PC_FIELDSPLIT_SCHUR_FACT_FULL:
873:     /* [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 */
874:     VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
875:     VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
876:     KSPSolve(kspLower,ilinkA->x,ilinkA->y);
877:     MatMult(jac->C,ilinkA->y,ilinkD->x);
878:     VecScale(ilinkD->x,-1.0);
879:     VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
880:     VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);

882:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
883:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
884:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);

886:     if (kspUpper == kspA) {
887:       MatMult(jac->B,ilinkD->y,ilinkA->y);
888:       VecAXPY(ilinkA->x,-1.0,ilinkA->y);
889:       KSPSolve(kspA,ilinkA->x,ilinkA->y);
890:     } else {
891:       KSPSolve(kspA,ilinkA->x,ilinkA->y);
892:       MatMult(jac->B,ilinkD->y,ilinkA->x);
893:       KSPSolve(kspUpper,ilinkA->x,ilinkA->z);
894:       VecAXPY(ilinkA->y,-1.0,ilinkA->z);
895:     }
896:     VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
897:     VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
898:   }
899:   return(0);
900: }

904: static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y)
905: {
906:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
907:   PetscErrorCode    ierr;
908:   PC_FieldSplitLink ilink = jac->head;
909:   PetscInt          cnt,bs;

912:   if (jac->type == PC_COMPOSITE_ADDITIVE) {
913:     if (jac->defaultsplit) {
914:       VecGetBlockSize(x,&bs);
915:       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);
916:       VecGetBlockSize(y,&bs);
917:       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);
918:       VecStrideGatherAll(x,jac->x,INSERT_VALUES);
919:       while (ilink) {
920:         KSPSolve(ilink->ksp,ilink->x,ilink->y);
921:         ilink = ilink->next;
922:       }
923:       VecStrideScatterAll(jac->y,y,INSERT_VALUES);
924:     } else {
925:       VecSet(y,0.0);
926:       while (ilink) {
927:         FieldSplitSplitSolveAdd(ilink,x,y);
928:         ilink = ilink->next;
929:       }
930:     }
931:   } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
932:     if (!jac->w1) {
933:       VecDuplicate(x,&jac->w1);
934:       VecDuplicate(x,&jac->w2);
935:     }
936:     VecSet(y,0.0);
937:     FieldSplitSplitSolveAdd(ilink,x,y);
938:     cnt  = 1;
939:     while (ilink->next) {
940:       ilink = ilink->next;
941:       /* compute the residual only over the part of the vector needed */
942:       MatMult(jac->Afield[cnt++],y,ilink->x);
943:       VecScale(ilink->x,-1.0);
944:       VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
945:       VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
946:       KSPSolve(ilink->ksp,ilink->x,ilink->y);
947:       VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
948:       VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
949:     }
950:     if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
951:       cnt -= 2;
952:       while (ilink->previous) {
953:         ilink = ilink->previous;
954:         /* compute the residual only over the part of the vector needed */
955:         MatMult(jac->Afield[cnt--],y,ilink->x);
956:         VecScale(ilink->x,-1.0);
957:         VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
958:         VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
959:         KSPSolve(ilink->ksp,ilink->x,ilink->y);
960:         VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
961:         VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
962:       }
963:     }
964:   } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type);
965:   return(0);
966: }

968: #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \
969:   (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
970:    VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
971:    KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \
972:    VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \
973:    VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE))

977: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y)
978: {
979:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
980:   PetscErrorCode    ierr;
981:   PC_FieldSplitLink ilink = jac->head;
982:   PetscInt          bs;

985:   if (jac->type == PC_COMPOSITE_ADDITIVE) {
986:     if (jac->defaultsplit) {
987:       VecGetBlockSize(x,&bs);
988:       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);
989:       VecGetBlockSize(y,&bs);
990:       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);
991:       VecStrideGatherAll(x,jac->x,INSERT_VALUES);
992:       while (ilink) {
993:         KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);
994:         ilink = ilink->next;
995:       }
996:       VecStrideScatterAll(jac->y,y,INSERT_VALUES);
997:     } else {
998:       VecSet(y,0.0);
999:       while (ilink) {
1000:         FieldSplitSplitSolveAddTranspose(ilink,x,y);
1001:         ilink = ilink->next;
1002:       }
1003:     }
1004:   } else {
1005:     if (!jac->w1) {
1006:       VecDuplicate(x,&jac->w1);
1007:       VecDuplicate(x,&jac->w2);
1008:     }
1009:     VecSet(y,0.0);
1010:     if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
1011:       FieldSplitSplitSolveAddTranspose(ilink,x,y);
1012:       while (ilink->next) {
1013:         ilink = ilink->next;
1014:         MatMultTranspose(pc->mat,y,jac->w1);
1015:         VecWAXPY(jac->w2,-1.0,jac->w1,x);
1016:         FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1017:       }
1018:       while (ilink->previous) {
1019:         ilink = ilink->previous;
1020:         MatMultTranspose(pc->mat,y,jac->w1);
1021:         VecWAXPY(jac->w2,-1.0,jac->w1,x);
1022:         FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1023:       }
1024:     } else {
1025:       while (ilink->next) {   /* get to last entry in linked list */
1026:         ilink = ilink->next;
1027:       }
1028:       FieldSplitSplitSolveAddTranspose(ilink,x,y);
1029:       while (ilink->previous) {
1030:         ilink = ilink->previous;
1031:         MatMultTranspose(pc->mat,y,jac->w1);
1032:         VecWAXPY(jac->w2,-1.0,jac->w1,x);
1033:         FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1034:       }
1035:     }
1036:   }
1037:   return(0);
1038: }

1042: static PetscErrorCode PCReset_FieldSplit(PC pc)
1043: {
1044:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1045:   PetscErrorCode    ierr;
1046:   PC_FieldSplitLink ilink = jac->head,next;

1049:   while (ilink) {
1050:     KSPReset(ilink->ksp);
1051:     VecDestroy(&ilink->x);
1052:     VecDestroy(&ilink->y);
1053:     VecDestroy(&ilink->z);
1054:     VecScatterDestroy(&ilink->sctx);
1055:     ISDestroy(&ilink->is);
1056:     ISDestroy(&ilink->is_col);
1057:     next  = ilink->next;
1058:     ilink = next;
1059:   }
1060:   PetscFree2(jac->x,jac->y);
1061:   if (jac->mat && jac->mat != jac->pmat) {
1062:     MatDestroyMatrices(jac->nsplits,&jac->mat);
1063:   } else if (jac->mat) {
1064:     jac->mat = NULL;
1065:   }
1066:   if (jac->pmat) {MatDestroyMatrices(jac->nsplits,&jac->pmat);}
1067:   if (jac->Afield) {MatDestroyMatrices(jac->nsplits,&jac->Afield);}
1068:   VecDestroy(&jac->w1);
1069:   VecDestroy(&jac->w2);
1070:   MatDestroy(&jac->schur);
1071:   MatDestroy(&jac->schurp);
1072:   MatDestroy(&jac->schur_user);
1073:   KSPDestroy(&jac->kspschur);
1074:   KSPDestroy(&jac->kspupper);
1075:   MatDestroy(&jac->B);
1076:   MatDestroy(&jac->C);
1077:   jac->reset = PETSC_TRUE;
1078:   return(0);
1079: }

1083: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
1084: {
1085:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1086:   PetscErrorCode    ierr;
1087:   PC_FieldSplitLink ilink = jac->head,next;

1090:   PCReset_FieldSplit(pc);
1091:   while (ilink) {
1092:     KSPDestroy(&ilink->ksp);
1093:     next  = ilink->next;
1094:     PetscFree(ilink->splitname);
1095:     PetscFree(ilink->fields);
1096:     PetscFree(ilink->fields_col);
1097:     PetscFree(ilink);
1098:     ilink = next;
1099:   }
1100:   PetscFree2(jac->x,jac->y);
1101:   PetscFree(pc->data);
1102:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);
1103:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);
1104:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);
1105:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);
1106:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);
1107:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);
1108:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);
1109:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);
1110:   return(0);
1111: }

1115: static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc)
1116: {
1117:   PetscErrorCode  ierr;
1118:   PetscInt        bs;
1119:   PetscBool       flg,stokes = PETSC_FALSE;
1120:   PC_FieldSplit   *jac = (PC_FieldSplit*)pc->data;
1121:   PCCompositeType ctype;

1124:   PetscOptionsHead("FieldSplit options");
1125:   PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);
1126:   PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);
1127:   if (flg) {
1128:     PCFieldSplitSetBlockSize(pc,bs);
1129:   }

1131:   PetscOptionsBool("-pc_fieldsplit_diag_use_amat","Use Amat (not Pmat) to extract diagonal fieldsplit blocks", "PCFieldSplitSetDiagUseAmat",pc->useAmat,&jac->diag_use_amat,NULL);
1132:   PetscOptionsBool("-pc_fieldsplit_off_diag_use_amat","Use Amat (not Pmat) to extract off-diagonal fieldsplit blocks", "PCFieldSplitSetOffDiagUseAmat",pc->useAmat,&jac->offdiag_use_amat,NULL);
1133:   /* FIXME: No programmatic equivalent to the following. */
1134:   PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);
1135:   if (stokes) {
1136:     PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);
1137:     jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF;
1138:   }

1140:   PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);
1141:   if (flg) {
1142:     PCFieldSplitSetType(pc,ctype);
1143:   }
1144:   /* Only setup fields once */
1145:   if ((jac->bs > 0) && (jac->nsplits == 0)) {
1146:     /* only allow user to set fields from command line if bs is already known.
1147:        otherwise user can set them in PCFieldSplitSetDefaults() */
1148:     PCFieldSplitSetRuntimeSplits_Private(pc);
1149:     if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
1150:   }
1151:   if (jac->type == PC_COMPOSITE_SCHUR) {
1152:     PetscOptionsGetEnum(((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);
1153:     if (flg) {PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");}
1154:     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);
1155:     PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);
1156:   }
1157:   PetscOptionsTail();
1158:   return(0);
1159: }

1161: /*------------------------------------------------------------------------------------*/

1165: static PetscErrorCode  PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1166: {
1167:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1168:   PetscErrorCode    ierr;
1169:   PC_FieldSplitLink ilink,next = jac->head;
1170:   char              prefix[128];
1171:   PetscInt          i;

1174:   if (jac->splitdefined) {
1175:     PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1176:     return(0);
1177:   }
1178:   for (i=0; i<n; i++) {
1179:     if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs);
1180:     if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]);
1181:   }
1182:   PetscNew(&ilink);
1183:   if (splitname) {
1184:     PetscStrallocpy(splitname,&ilink->splitname);
1185:   } else {
1186:     PetscMalloc1(3,&ilink->splitname);
1187:     PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);
1188:   }
1189:   PetscMalloc1(n,&ilink->fields);
1190:   PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));
1191:   PetscMalloc1(n,&ilink->fields_col);
1192:   PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));

1194:   ilink->nfields = n;
1195:   ilink->next    = NULL;
1196:   KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1197:   PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1198:   KSPSetType(ilink->ksp,KSPPREONLY);
1199:   PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);

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

1204:   if (!next) {
1205:     jac->head       = ilink;
1206:     ilink->previous = NULL;
1207:   } else {
1208:     while (next->next) {
1209:       next = next->next;
1210:     }
1211:     next->next      = ilink;
1212:     ilink->previous = next;
1213:   }
1214:   jac->nsplits++;
1215:   return(0);
1216: }

1220: static PetscErrorCode  PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp)
1221: {
1222:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

1226:   PetscMalloc1(jac->nsplits,subksp);
1227:   MatSchurComplementGetKSP(jac->schur,*subksp);

1229:   (*subksp)[1] = jac->kspschur;
1230:   if (n) *n = jac->nsplits;
1231:   return(0);
1232: }

1236: static PetscErrorCode  PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1237: {
1238:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1239:   PetscErrorCode    ierr;
1240:   PetscInt          cnt   = 0;
1241:   PC_FieldSplitLink ilink = jac->head;

1244:   PetscMalloc1(jac->nsplits,subksp);
1245:   while (ilink) {
1246:     (*subksp)[cnt++] = ilink->ksp;
1247:     ilink            = ilink->next;
1248:   }
1249:   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);
1250:   if (n) *n = jac->nsplits;
1251:   return(0);
1252: }

1256: static PetscErrorCode  PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is)
1257: {
1258:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1259:   PetscErrorCode    ierr;
1260:   PC_FieldSplitLink ilink, next = jac->head;
1261:   char              prefix[128];

1264:   if (jac->splitdefined) {
1265:     PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1266:     return(0);
1267:   }
1268:   PetscNew(&ilink);
1269:   if (splitname) {
1270:     PetscStrallocpy(splitname,&ilink->splitname);
1271:   } else {
1272:     PetscMalloc1(8,&ilink->splitname);
1273:     PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);
1274:   }
1275:   PetscObjectReference((PetscObject)is);
1276:   ISDestroy(&ilink->is);
1277:   ilink->is     = is;
1278:   PetscObjectReference((PetscObject)is);
1279:   ISDestroy(&ilink->is_col);
1280:   ilink->is_col = is;
1281:   ilink->next   = NULL;
1282:   KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1283:   PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1284:   KSPSetType(ilink->ksp,KSPPREONLY);
1285:   PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);

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

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

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

1309:     Logically Collective on PC

1311:     Input Parameters:
1312: +   pc  - the preconditioner context
1313: .   splitname - name of this split, if NULL the number of the split is used
1314: .   n - the number of fields in this split
1315: -   fields - the fields in this split

1317:     Level: intermediate

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

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

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

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

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

1335: @*/
1336: PetscErrorCode  PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1337: {

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

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

1354:     Logically Collective on PC

1356:     Input Parameters:
1357: +   pc  - the preconditioner object
1358: -   flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from

1360:     Options Database:
1361: .     -pc_fieldsplit_diag_use_amat

1363:     Level: intermediate

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

1367: @*/
1368: PetscErrorCode  PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg)
1369: {
1370:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1371:   PetscBool      isfs;

1376:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1377:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1378:   jac->diag_use_amat = flg;
1379:   return(0);
1380: }

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

1387:     Logically Collective on PC

1389:     Input Parameters:
1390: .   pc  - the preconditioner object

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


1396:     Level: intermediate

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

1400: @*/
1401: PetscErrorCode  PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg)
1402: {
1403:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1404:   PetscBool      isfs;

1410:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1411:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1412:   *flg = jac->diag_use_amat;
1413:   return(0);
1414: }

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

1421:     Logically Collective on PC

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

1427:     Options Database:
1428: .     -pc_fieldsplit_off_diag_use_amat

1430:     Level: intermediate

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

1434: @*/
1435: PetscErrorCode  PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg)
1436: {
1437:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1438:   PetscBool      isfs;

1443:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1444:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1445:   jac->offdiag_use_amat = flg;
1446:   return(0);
1447: }

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

1454:     Logically Collective on PC

1456:     Input Parameters:
1457: .   pc  - the preconditioner object

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


1463:     Level: intermediate

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

1467: @*/
1468: PetscErrorCode  PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg)
1469: {
1470:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1471:   PetscBool      isfs;

1477:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1478:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1479:   *flg = jac->offdiag_use_amat;
1480:   return(0);
1481: }



1487: /*@
1488:     PCFieldSplitSetIS - Sets the exact elements for field

1490:     Logically Collective on PC

1492:     Input Parameters:
1493: +   pc  - the preconditioner context
1494: .   splitname - name of this split, if NULL the number of the split is used
1495: -   is - the index set that defines the vector elements in this field


1498:     Notes:
1499:     Use PCFieldSplitSetFields(), for fields defined by strided types.

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

1504:     Level: intermediate

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

1508: @*/
1509: PetscErrorCode  PCFieldSplitSetIS(PC pc,const char splitname[],IS is)
1510: {

1517:   PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));
1518:   return(0);
1519: }

1523: /*@
1524:     PCFieldSplitGetIS - Retrieves the elements for a field as an IS

1526:     Logically Collective on PC

1528:     Input Parameters:
1529: +   pc  - the preconditioner context
1530: -   splitname - name of this split

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

1535:     Level: intermediate

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

1539: @*/
1540: PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is)
1541: {

1548:   {
1549:     PC_FieldSplit     *jac  = (PC_FieldSplit*) pc->data;
1550:     PC_FieldSplitLink ilink = jac->head;
1551:     PetscBool         found;

1553:     *is = NULL;
1554:     while (ilink) {
1555:       PetscStrcmp(ilink->splitname, splitname, &found);
1556:       if (found) {
1557:         *is = ilink->is;
1558:         break;
1559:       }
1560:       ilink = ilink->next;
1561:     }
1562:   }
1563:   return(0);
1564: }

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

1572:     Logically Collective on PC

1574:     Input Parameters:
1575: +   pc  - the preconditioner context
1576: -   bs - the block size

1578:     Level: intermediate

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

1582: @*/
1583: PetscErrorCode  PCFieldSplitSetBlockSize(PC pc,PetscInt bs)
1584: {

1590:   PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));
1591:   return(0);
1592: }

1596: /*@C
1597:    PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits

1599:    Collective on KSP

1601:    Input Parameter:
1602: .  pc - the preconditioner context

1604:    Output Parameters:
1605: +  n - the number of splits
1606: -  pc - the array of KSP contexts

1608:    Note:
1609:    After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user
1610:    (not the KSP just the array that contains them).

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

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


1619:    Level: advanced

1621: .seealso: PCFIELDSPLIT
1622: @*/
1623: PetscErrorCode  PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
1624: {

1630:   PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
1631:   return(0);
1632: }

1636: /*@
1637:     PCFieldSplitSetSchurPre -  Indicates if the Schur complement is preconditioned by a preconditioner constructed by the
1638:       A11 matrix. Otherwise no preconditioner is used.

1640:     Collective on PC

1642:     Input Parameters:
1643: +   pc      - the preconditioner context
1644: .   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
1645: -   userpre - matrix to use for preconditioning, or NULL

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

1650:     Notes:
1651: $    If ptype is
1652: $        user then the preconditioner for the Schur complement is generated by the provided matrix (pre argument
1653: $             to this function).
1654: $        a11 then the preconditioner for the Schur complement is generated by the block diagonal part of the preconditioner
1655: $             matrix associated with the Schur complement (i.e. A11), not he Schur complement matrix
1656: $        full then the preconditioner uses the exact Schur complement (this is expensive)
1657: $        self the preconditioner for the Schur complement is generated from the Schur complement matrix itself:
1658: $             The only preconditioner that currently works directly with the Schur complement matrix object is the PCLSC
1659: $             preconditioner
1660: $        selfp then the preconditioning matrix is an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01
1661: $             This is only a good preconditioner when diag(A00) is a good preconditioner for A00.

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

1667:     Level: intermediate

1669: .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC

1671: @*/
1672: PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
1673: {

1678:   PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));
1679:   return(0);
1680: }
1681: PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */

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

1689:     Logically Collective on PC

1691:     Input Parameters:
1692: .   pc      - the preconditioner context

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

1698:     Level: intermediate

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

1702: @*/
1703: PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
1704: {

1709:   PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));
1710:   return(0);
1711: }

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

1718:     Not collective

1720:     Input Parameter:
1721: .   pc      - the preconditioner context

1723:     Output Parameter:
1724: .   S       - the Schur complement matrix

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

1729:     Level: advanced

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

1733: @*/
1734: PetscErrorCode  PCFieldSplitSchurGetS(PC pc,Mat *S)
1735: {
1737:   const char*    t;
1738:   PetscBool      isfs;
1739:   PC_FieldSplit  *jac;

1743:   PetscObjectGetType((PetscObject)pc,&t);
1744:   PetscStrcmp(t,PCFIELDSPLIT,&isfs);
1745:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
1746:   jac = (PC_FieldSplit*)pc->data;
1747:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
1748:   if (S) *S = jac->schur;
1749:   return(0);
1750: }

1754: /*@
1755:     PCFieldSplitSchurRestoreS -  restores the MatSchurComplement object used by this PC

1757:     Not collective

1759:     Input Parameters:
1760: +   pc      - the preconditioner context
1761: .   S       - the Schur complement matrix

1763:     Level: advanced

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

1767: @*/
1768: PetscErrorCode  PCFieldSplitSchurRestoreS(PC pc,Mat *S)
1769: {
1771:   const char*    t;
1772:   PetscBool      isfs;
1773:   PC_FieldSplit  *jac;

1777:   PetscObjectGetType((PetscObject)pc,&t);
1778:   PetscStrcmp(t,PCFIELDSPLIT,&isfs);
1779:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
1780:   jac = (PC_FieldSplit*)pc->data;
1781:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
1782:   if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten");
1783:   return(0);
1784: }


1789: static PetscErrorCode  PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
1790: {
1791:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

1795:   jac->schurpre = ptype;
1796:   if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
1797:     MatDestroy(&jac->schur_user);
1798:     jac->schur_user = pre;
1799:     PetscObjectReference((PetscObject)jac->schur_user);
1800:   }
1801:   return(0);
1802: }

1806: static PetscErrorCode  PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
1807: {
1808:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

1811:   *ptype = jac->schurpre;
1812:   *pre   = jac->schur_user;
1813:   return(0);
1814: }

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

1821:     Collective on PC

1823:     Input Parameters:
1824: +   pc  - the preconditioner context
1825: -   ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default

1827:     Options Database:
1828: .     -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full


1831:     Level: intermediate

1833:     Notes:
1834:     The FULL factorization is

1836: $   (A   B)  = (1       0) (A   0) (1  Ainv*B)
1837: $   (C   D)    (C*Ainv  1) (0   S) (0     1  )

1839:     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,
1840:     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).

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

1847:     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
1848:     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).

1850:     References:
1851:     Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000) pp. 1969-1972.

1853:     Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001), pp. 1050-1051.

1855: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType
1856: @*/
1857: PetscErrorCode  PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype)
1858: {

1863:   PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));
1864:   return(0);
1865: }

1869: static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype)
1870: {
1871:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

1874:   jac->schurfactorization = ftype;
1875:   return(0);
1876: }

1880: /*@C
1881:    PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement

1883:    Collective on KSP

1885:    Input Parameter:
1886: .  pc - the preconditioner context

1888:    Output Parameters:
1889: +  A00 - the (0,0) block
1890: .  A01 - the (0,1) block
1891: .  A10 - the (1,0) block
1892: -  A11 - the (1,1) block

1894:    Level: advanced

1896: .seealso: PCFIELDSPLIT
1897: @*/
1898: PetscErrorCode  PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11)
1899: {
1900:   PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;

1904:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach.");
1905:   if (A00) *A00 = jac->pmat[0];
1906:   if (A01) *A01 = jac->B;
1907:   if (A10) *A10 = jac->C;
1908:   if (A11) *A11 = jac->pmat[1];
1909:   return(0);
1910: }

1914: static PetscErrorCode  PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type)
1915: {
1916:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

1920:   jac->type = type;
1921:   if (type == PC_COMPOSITE_SCHUR) {
1922:     pc->ops->apply = PCApply_FieldSplit_Schur;
1923:     pc->ops->view  = PCView_FieldSplit_Schur;

1925:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);
1926:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);
1927:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);
1928:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);

1930:   } else {
1931:     pc->ops->apply = PCApply_FieldSplit;
1932:     pc->ops->view  = PCView_FieldSplit;

1934:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
1935:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);
1936:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);
1937:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);
1938:   }
1939:   return(0);
1940: }

1944: static PetscErrorCode  PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs)
1945: {
1946:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

1949:   if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs);
1950:   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);
1951:   jac->bs = bs;
1952:   return(0);
1953: }

1957: /*@
1958:    PCFieldSplitSetType - Sets the type of fieldsplit preconditioner.

1960:    Collective on PC

1962:    Input Parameter:
1963: .  pc - the preconditioner context
1964: .  type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR

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

1969:    Level: Intermediate

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

1973: .seealso: PCCompositeSetType()

1975: @*/
1976: PetscErrorCode  PCFieldSplitSetType(PC pc,PCCompositeType type)
1977: {

1982:   PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));
1983:   return(0);
1984: }

1988: /*@
1989:   PCFieldSplitGetType - Gets the type of fieldsplit preconditioner.

1991:   Not collective

1993:   Input Parameter:
1994: . pc - the preconditioner context

1996:   Output Parameter:
1997: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR

1999:   Level: Intermediate

2001: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
2002: .seealso: PCCompositeSetType()
2003: @*/
2004: PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
2005: {
2006:   PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;

2011:   *type = jac->type;
2012:   return(0);
2013: }

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

2020:    Logically Collective

2022:    Input Parameters:
2023: +  pc   - the preconditioner context
2024: -  flg  - boolean indicating whether to use field splits defined by the DM

2026:    Options Database Key:
2027: .  -pc_fieldsplit_dm_splits

2029:    Level: Intermediate

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

2033: .seealso: PCFieldSplitGetDMSplits()

2035: @*/
2036: PetscErrorCode  PCFieldSplitSetDMSplits(PC pc,PetscBool flg)
2037: {
2038:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
2039:   PetscBool      isfs;

2045:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2046:   if (isfs) {
2047:     jac->dm_splits = flg;
2048:   }
2049:   return(0);
2050: }


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

2058:    Logically Collective

2060:    Input Parameter:
2061: .  pc   - the preconditioner context

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

2066:    Level: Intermediate

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

2070: .seealso: PCFieldSplitSetDMSplits()

2072: @*/
2073: PetscErrorCode  PCFieldSplitGetDMSplits(PC pc,PetscBool* flg)
2074: {
2075:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
2076:   PetscBool      isfs;

2082:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2083:   if (isfs) {
2084:     if(flg) *flg = jac->dm_splits;
2085:   }
2086:   return(0);
2087: }

2089: /* -------------------------------------------------------------------------------------*/
2090: /*MC
2091:    PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
2092:                   fields or groups of fields. See the users manual section "Solving Block Matrices" for more details.

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

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

2100:    Level: intermediate

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

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

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

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

2123: $     For the Schur complement preconditioner if J = ( A00 A01 )
2124: $                                                    ( A10 A11 )
2125: $     the preconditioner using full factorization is
2126: $              ( I   -ksp(A00) A01 ) ( inv(A00)     0  ) (     I          0  )
2127: $              ( 0         I       ) (   0      ksp(S) ) ( -A10 ksp(A00)  I  )
2128:      where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_.  S is the Schur complement
2129: $              S = A11 - A10 ksp(A00) A01
2130:      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
2131:      in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0,
2132:      it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default
2133:      A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() to turn on or off this
2134:      option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc.
2135:      When option -fieldsplit_schur_precondition selfp is given, an approximation to S is assembled --
2136:      Sp = A11 - A10 inv(diag(A00)) A01, which has type AIJ and can be used with a variety of preconditioners
2137:      (e.g., -fieldsplit_1_pc_type asm).
2138:      The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above,
2139:      diag gives
2140: $              ( inv(A00)     0   )
2141: $              (   0      -ksp(S) )
2142:      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
2143: $              (  A00   0 )
2144: $              (  A10   S )
2145:      where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of
2146: $              ( A00 A01 )
2147: $              (  0   S  )
2148:      where again the inverses of A00 and S are applied using KSPs.

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

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

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

2160:    Concepts: physics based preconditioners, block preconditioners

2162:    There is a nice discussion of block preconditioners in

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

2168: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC,
2169:            PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre()
2170: M*/

2174: PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
2175: {
2177:   PC_FieldSplit  *jac;

2180:   PetscNewLog(pc,&jac);

2182:   jac->bs                 = -1;
2183:   jac->nsplits            = 0;
2184:   jac->type               = PC_COMPOSITE_MULTIPLICATIVE;
2185:   jac->schurpre           = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
2186:   jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
2187:   jac->dm_splits          = PETSC_TRUE;

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

2191:   pc->ops->apply           = PCApply_FieldSplit;
2192:   pc->ops->applytranspose  = PCApplyTranspose_FieldSplit;
2193:   pc->ops->setup           = PCSetUp_FieldSplit;
2194:   pc->ops->reset           = PCReset_FieldSplit;
2195:   pc->ops->destroy         = PCDestroy_FieldSplit;
2196:   pc->ops->setfromoptions  = PCSetFromOptions_FieldSplit;
2197:   pc->ops->view            = PCView_FieldSplit;
2198:   pc->ops->applyrichardson = 0;

2200:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2201:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);
2202:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);
2203:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);
2204:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);
2205:   return(0);
2206: }