Actual source code: fieldsplit.c

petsc-3.5.2 2014-09-08
Report Typos and Errors
  2: #include <petsc-private/pcimpl.h>     /*I "petscpc.h" I*/
  3: #include <petscdm.h>

  5: /*
  6:   There is a nice discussion of block preconditioners in

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

410:       if (stokes) {
411:         IS       zerodiags,rest;
412:         PetscInt nmin,nmax;

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

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

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


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

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

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

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

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

500:     jac->issetup = PETSC_TRUE;

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

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

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

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

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

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

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

629:   if (jac->type == PC_COMPOSITE_SCHUR) {
630:     IS          ccis;
631:     PetscInt    rstart,rend;
632:     char        lscname[256];
633:     PetscObject LSC_L;

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

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

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

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

680:       /* extract the A01 and A10 matrices */
681:       ilink = jac->head;
682:       ISComplement(ilink->is_col,rstart,rend,&ccis);
683:       if (jac->offdiag_use_amat) {
684:         MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
685:       } else {
686:         MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
687:       }
688:       ISDestroy(&ccis);
689:       ilink = ilink->next;
690:       ISComplement(ilink->is_col,rstart,rend,&ccis);
691:       if (jac->offdiag_use_amat) {
692:         MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
693:       } else {
694:         MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
695:       }
696:       ISDestroy(&ccis);

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

711:       PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);
712:       PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
713:       if (flg) {
714:         DM  dmInner;
715:         KSP kspInner;

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

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

741:       PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);
742:       PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
743:       if (flg) {
744:         DM dmInner;

746:         PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
747:         KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);
748:         KSPSetOptionsPrefix(jac->kspupper, schurprefix);
749:         KSPGetDM(jac->head->ksp, &dmInner);
750:         KSPSetDM(jac->kspupper, dmInner);
751:         KSPSetDMActive(jac->kspupper, PETSC_FALSE);
752:         KSPSetFromOptions(jac->kspupper);
753:         KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0]);
754:         VecDuplicate(jac->head->x, &jac->head->z);
755:       } else {
756:         jac->kspupper = jac->head->ksp;
757:         PetscObjectReference((PetscObject) jac->head->ksp);
758:       }

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

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

813:   jac->suboptionsset = PETSC_TRUE;
814:   return(0);
815: }

817: #define FieldSplitSplitSolveAdd(ilink,xx,yy) \
818:   (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
819:    VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
820:    KSPSolve(ilink->ksp,ilink->x,ilink->y) || \
821:    VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \
822:    VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE))

826: static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y)
827: {
828:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
829:   PetscErrorCode    ierr;
830:   PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
831:   KSP               kspA   = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;

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

889:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
890:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
891:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);

893:     if (kspUpper == kspA) {
894:       MatMult(jac->B,ilinkD->y,ilinkA->y);
895:       VecAXPY(ilinkA->x,-1.0,ilinkA->y);
896:       KSPSolve(kspA,ilinkA->x,ilinkA->y);
897:     } else {
898:       KSPSolve(kspA,ilinkA->x,ilinkA->y);
899:       MatMult(jac->B,ilinkD->y,ilinkA->x);
900:       KSPSolve(kspUpper,ilinkA->x,ilinkA->z);
901:       VecAXPY(ilinkA->y,-1.0,ilinkA->z);
902:     }
903:     VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
904:     VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
905:   }
906:   return(0);
907: }

911: static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y)
912: {
913:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
914:   PetscErrorCode    ierr;
915:   PC_FieldSplitLink ilink = jac->head;
916:   PetscInt          cnt,bs;

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

975: #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \
976:   (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
977:    VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
978:    KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \
979:    VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \
980:    VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE))

984: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y)
985: {
986:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
987:   PetscErrorCode    ierr;
988:   PC_FieldSplitLink ilink = jac->head;
989:   PetscInt          bs;

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

1049: static PetscErrorCode PCReset_FieldSplit(PC pc)
1050: {
1051:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1052:   PetscErrorCode    ierr;
1053:   PC_FieldSplitLink ilink = jac->head,next;

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

1090: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
1091: {
1092:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1093:   PetscErrorCode    ierr;
1094:   PC_FieldSplitLink ilink = jac->head,next;

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

1122: static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc)
1123: {
1124:   PetscErrorCode  ierr;
1125:   PetscInt        bs;
1126:   PetscBool       flg,stokes = PETSC_FALSE;
1127:   PC_FieldSplit   *jac = (PC_FieldSplit*)pc->data;
1128:   PCCompositeType ctype;

1131:   PetscOptionsHead("FieldSplit options");
1132:   PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);
1133:   PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);
1134:   if (flg) {
1135:     PCFieldSplitSetBlockSize(pc,bs);
1136:   }

1138:   PetscOptionsBool("-pc_fieldsplit_diag_use_amat","Use Amat (not Pmat) to extract diagonal fieldsplit blocks", "PCFieldSplitSetDiagUseAmat",pc->useAmat,&jac->diag_use_amat,NULL);
1139:   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);
1140:   /* FIXME: No programmatic equivalent to the following. */
1141:   PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);
1142:   if (stokes) {
1143:     PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);
1144:     jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF;
1145:   }

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

1168: /*------------------------------------------------------------------------------------*/

1172: static PetscErrorCode  PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1173: {
1174:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1175:   PetscErrorCode    ierr;
1176:   PC_FieldSplitLink ilink,next = jac->head;
1177:   char              prefix[128];
1178:   PetscInt          i;

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

1201:   ilink->nfields = n;
1202:   ilink->next    = NULL;
1203:   KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1204:   PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1205:   KSPSetType(ilink->ksp,KSPPREONLY);
1206:   PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);

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

1211:   if (!next) {
1212:     jac->head       = ilink;
1213:     ilink->previous = NULL;
1214:   } else {
1215:     while (next->next) {
1216:       next = next->next;
1217:     }
1218:     next->next      = ilink;
1219:     ilink->previous = next;
1220:   }
1221:   jac->nsplits++;
1222:   return(0);
1223: }

1227: static PetscErrorCode  PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp)
1228: {
1229:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

1233:   PetscMalloc1(jac->nsplits,subksp);
1234:   MatSchurComplementGetKSP(jac->schur,*subksp);

1236:   (*subksp)[1] = jac->kspschur;
1237:   if (n) *n = jac->nsplits;
1238:   return(0);
1239: }

1243: static PetscErrorCode  PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1244: {
1245:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1246:   PetscErrorCode    ierr;
1247:   PetscInt          cnt   = 0;
1248:   PC_FieldSplitLink ilink = jac->head;

1251:   PetscMalloc1(jac->nsplits,subksp);
1252:   while (ilink) {
1253:     (*subksp)[cnt++] = ilink->ksp;
1254:     ilink            = ilink->next;
1255:   }
1256:   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);
1257:   if (n) *n = jac->nsplits;
1258:   return(0);
1259: }

1263: static PetscErrorCode  PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is)
1264: {
1265:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1266:   PetscErrorCode    ierr;
1267:   PC_FieldSplitLink ilink, next = jac->head;
1268:   char              prefix[128];

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

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

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

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

1316:     Logically Collective on PC

1318:     Input Parameters:
1319: +   pc  - the preconditioner context
1320: .   splitname - name of this split, if NULL the number of the split is used
1321: .   n - the number of fields in this split
1322: -   fields - the fields in this split

1324:     Level: intermediate

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

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

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

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

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

1342: @*/
1343: PetscErrorCode  PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1344: {

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

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

1361:     Logically Collective on PC

1363:     Input Parameters:
1364: +   pc  - the preconditioner object
1365: -   flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from

1367:     Options Database:
1368: .     -pc_fieldsplit_diag_use_amat

1370:     Level: intermediate

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

1374: @*/
1375: PetscErrorCode  PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg)
1376: {
1377:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1378:   PetscBool      isfs;

1383:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1384:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1385:   jac->diag_use_amat = flg;
1386:   return(0);
1387: }

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

1394:     Logically Collective on PC

1396:     Input Parameters:
1397: .   pc  - the preconditioner object

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


1403:     Level: intermediate

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

1407: @*/
1408: PetscErrorCode  PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg)
1409: {
1410:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1411:   PetscBool      isfs;

1417:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1418:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1419:   *flg = jac->diag_use_amat;
1420:   return(0);
1421: }

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

1428:     Logically Collective on PC

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

1434:     Options Database:
1435: .     -pc_fieldsplit_off_diag_use_amat

1437:     Level: intermediate

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

1441: @*/
1442: PetscErrorCode  PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg)
1443: {
1444:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;
1445:   PetscBool      isfs;

1450:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1451:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1452:   jac->offdiag_use_amat = flg;
1453:   return(0);
1454: }

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

1461:     Logically Collective on PC

1463:     Input Parameters:
1464: .   pc  - the preconditioner object

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


1470:     Level: intermediate

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

1474: @*/
1475: PetscErrorCode  PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg)
1476: {
1477:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1478:   PetscBool      isfs;

1484:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1485:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT);
1486:   *flg = jac->offdiag_use_amat;
1487:   return(0);
1488: }



1494: /*@
1495:     PCFieldSplitSetIS - Sets the exact elements for field

1497:     Logically Collective on PC

1499:     Input Parameters:
1500: +   pc  - the preconditioner context
1501: .   splitname - name of this split, if NULL the number of the split is used
1502: -   is - the index set that defines the vector elements in this field


1505:     Notes:
1506:     Use PCFieldSplitSetFields(), for fields defined by strided types.

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

1511:     Level: intermediate

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

1515: @*/
1516: PetscErrorCode  PCFieldSplitSetIS(PC pc,const char splitname[],IS is)
1517: {

1524:   PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));
1525:   return(0);
1526: }

1530: /*@
1531:     PCFieldSplitGetIS - Retrieves the elements for a field as an IS

1533:     Logically Collective on PC

1535:     Input Parameters:
1536: +   pc  - the preconditioner context
1537: -   splitname - name of this split

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

1542:     Level: intermediate

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

1546: @*/
1547: PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is)
1548: {

1555:   {
1556:     PC_FieldSplit     *jac  = (PC_FieldSplit*) pc->data;
1557:     PC_FieldSplitLink ilink = jac->head;
1558:     PetscBool         found;

1560:     *is = NULL;
1561:     while (ilink) {
1562:       PetscStrcmp(ilink->splitname, splitname, &found);
1563:       if (found) {
1564:         *is = ilink->is;
1565:         break;
1566:       }
1567:       ilink = ilink->next;
1568:     }
1569:   }
1570:   return(0);
1571: }

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

1579:     Logically Collective on PC

1581:     Input Parameters:
1582: +   pc  - the preconditioner context
1583: -   bs - the block size

1585:     Level: intermediate

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

1589: @*/
1590: PetscErrorCode  PCFieldSplitSetBlockSize(PC pc,PetscInt bs)
1591: {

1597:   PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));
1598:   return(0);
1599: }

1603: /*@C
1604:    PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits

1606:    Collective on KSP

1608:    Input Parameter:
1609: .  pc - the preconditioner context

1611:    Output Parameters:
1612: +  n - the number of splits
1613: -  pc - the array of KSP contexts

1615:    Note:
1616:    After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user
1617:    (not the KSP just the array that contains them).

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

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


1626:    Level: advanced

1628: .seealso: PCFIELDSPLIT
1629: @*/
1630: PetscErrorCode  PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
1631: {

1637:   PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
1638:   return(0);
1639: }

1643: /*@
1644:     PCFieldSplitSetSchurPre -  Indicates if the Schur complement is preconditioned by a preconditioner constructed by the
1645:       A11 matrix. Otherwise no preconditioner is used.

1647:     Collective on PC

1649:     Input Parameters:
1650: +   pc      - the preconditioner context
1651: .   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
1652: -   userpre - matrix to use for preconditioning, or NULL

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

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

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

1674:     Level: intermediate

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

1678: @*/
1679: PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
1680: {

1685:   PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));
1686:   return(0);
1687: }
1688: PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */

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

1696:     Logically Collective on PC

1698:     Input Parameters:
1699: .   pc      - the preconditioner context

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

1705:     Level: intermediate

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

1709: @*/
1710: PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
1711: {

1716:   PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));
1717:   return(0);
1718: }

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

1725:     Not collective

1727:     Input Parameter:
1728: .   pc      - the preconditioner context

1730:     Output Parameter:
1731: .   S       - the Schur complement matrix

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

1736:     Level: advanced

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

1740: @*/
1741: PetscErrorCode  PCFieldSplitSchurGetS(PC pc,Mat *S)
1742: {
1744:   const char*    t;
1745:   PetscBool      isfs;
1746:   PC_FieldSplit  *jac;

1750:   PetscObjectGetType((PetscObject)pc,&t);
1751:   PetscStrcmp(t,PCFIELDSPLIT,&isfs);
1752:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
1753:   jac = (PC_FieldSplit*)pc->data;
1754:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
1755:   if (S) *S = jac->schur;
1756:   return(0);
1757: }

1761: /*@
1762:     PCFieldSplitSchurRestoreS -  restores the MatSchurComplement object used by this PC

1764:     Not collective

1766:     Input Parameters:
1767: +   pc      - the preconditioner context
1768: .   S       - the Schur complement matrix

1770:     Level: advanced

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

1774: @*/
1775: PetscErrorCode  PCFieldSplitSchurRestoreS(PC pc,Mat *S)
1776: {
1778:   const char*    t;
1779:   PetscBool      isfs;
1780:   PC_FieldSplit  *jac;

1784:   PetscObjectGetType((PetscObject)pc,&t);
1785:   PetscStrcmp(t,PCFIELDSPLIT,&isfs);
1786:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
1787:   jac = (PC_FieldSplit*)pc->data;
1788:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
1789:   if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten");
1790:   return(0);
1791: }


1796: static PetscErrorCode  PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
1797: {
1798:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

1802:   jac->schurpre = ptype;
1803:   if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
1804:     MatDestroy(&jac->schur_user);
1805:     jac->schur_user = pre;
1806:     PetscObjectReference((PetscObject)jac->schur_user);
1807:   }
1808:   return(0);
1809: }

1813: static PetscErrorCode  PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
1814: {
1815:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

1818:   *ptype = jac->schurpre;
1819:   *pre   = jac->schur_user;
1820:   return(0);
1821: }

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

1828:     Collective on PC

1830:     Input Parameters:
1831: +   pc  - the preconditioner context
1832: -   ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default

1834:     Options Database:
1835: .     -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full


1838:     Level: intermediate

1840:     Notes:
1841:     The FULL factorization is

1843: $   (A   B)  = (1       0) (A   0) (1  Ainv*B)
1844: $   (C   D)    (C*Ainv  1) (0   S) (0     1  )

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

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

1854:     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
1855:     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).

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

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

1862: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType
1863: @*/
1864: PetscErrorCode  PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype)
1865: {

1870:   PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));
1871:   return(0);
1872: }

1876: static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype)
1877: {
1878:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

1881:   jac->schurfactorization = ftype;
1882:   return(0);
1883: }

1887: /*@C
1888:    PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement

1890:    Collective on KSP

1892:    Input Parameter:
1893: .  pc - the preconditioner context

1895:    Output Parameters:
1896: +  A00 - the (0,0) block
1897: .  A01 - the (0,1) block
1898: .  A10 - the (1,0) block
1899: -  A11 - the (1,1) block

1901:    Level: advanced

1903: .seealso: PCFIELDSPLIT
1904: @*/
1905: PetscErrorCode  PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11)
1906: {
1907:   PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;

1911:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach.");
1912:   if (A00) *A00 = jac->pmat[0];
1913:   if (A01) *A01 = jac->B;
1914:   if (A10) *A10 = jac->C;
1915:   if (A11) *A11 = jac->pmat[1];
1916:   return(0);
1917: }

1921: static PetscErrorCode  PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type)
1922: {
1923:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

1927:   jac->type = type;
1928:   if (type == PC_COMPOSITE_SCHUR) {
1929:     pc->ops->apply = PCApply_FieldSplit_Schur;
1930:     pc->ops->view  = PCView_FieldSplit_Schur;

1932:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);
1933:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);
1934:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);
1935:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);

1937:   } else {
1938:     pc->ops->apply = PCApply_FieldSplit;
1939:     pc->ops->view  = PCView_FieldSplit;

1941:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
1942:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);
1943:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);
1944:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);
1945:   }
1946:   return(0);
1947: }

1951: static PetscErrorCode  PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs)
1952: {
1953:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

1956:   if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs);
1957:   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);
1958:   jac->bs = bs;
1959:   return(0);
1960: }

1964: /*@
1965:    PCFieldSplitSetType - Sets the type of fieldsplit preconditioner.

1967:    Collective on PC

1969:    Input Parameter:
1970: .  pc - the preconditioner context
1971: .  type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR

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

1976:    Level: Intermediate

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

1980: .seealso: PCCompositeSetType()

1982: @*/
1983: PetscErrorCode  PCFieldSplitSetType(PC pc,PCCompositeType type)
1984: {

1989:   PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));
1990:   return(0);
1991: }

1995: /*@
1996:   PCFieldSplitGetType - Gets the type of fieldsplit preconditioner.

1998:   Not collective

2000:   Input Parameter:
2001: . pc - the preconditioner context

2003:   Output Parameter:
2004: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR

2006:   Level: Intermediate

2008: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
2009: .seealso: PCCompositeSetType()
2010: @*/
2011: PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
2012: {
2013:   PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;

2018:   *type = jac->type;
2019:   return(0);
2020: }

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

2027:    Logically Collective

2029:    Input Parameters:
2030: +  pc   - the preconditioner context
2031: -  flg  - boolean indicating whether to use field splits defined by the DM

2033:    Options Database Key:
2034: .  -pc_fieldsplit_dm_splits

2036:    Level: Intermediate

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

2040: .seealso: PCFieldSplitGetDMSplits()

2042: @*/
2043: PetscErrorCode  PCFieldSplitSetDMSplits(PC pc,PetscBool flg)
2044: {
2045:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
2046:   PetscBool      isfs;

2052:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2053:   if (isfs) {
2054:     jac->dm_splits = flg;
2055:   }
2056:   return(0);
2057: }


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

2065:    Logically Collective

2067:    Input Parameter:
2068: .  pc   - the preconditioner context

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

2073:    Level: Intermediate

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

2077: .seealso: PCFieldSplitSetDMSplits()

2079: @*/
2080: PetscErrorCode  PCFieldSplitGetDMSplits(PC pc,PetscBool* flg)
2081: {
2082:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
2083:   PetscBool      isfs;

2089:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
2090:   if (isfs) {
2091:     if(flg) *flg = jac->dm_splits;
2092:   }
2093:   return(0);
2094: }

2096: /* -------------------------------------------------------------------------------------*/
2097: /*MC
2098:    PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
2099:                   fields or groups of fields. See the users manual section "Solving Block Matrices" for more details.

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

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

2107:    Level: intermediate

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

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

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

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

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

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

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

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

2167:    Concepts: physics based preconditioners, block preconditioners

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

2175: PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
2176: {
2178:   PC_FieldSplit  *jac;

2181:   PetscNewLog(pc,&jac);

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

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

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

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