Actual source code: fieldsplit.c

petsc-dev 2014-04-16
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: } PC_FieldSplit;

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

492:     jac->issetup = PETSC_TRUE;

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

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

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

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

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

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

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

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

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

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

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

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

663:       /* extract the A01 and A10 matrices */
664:       ilink = jac->head;
665:       ISComplement(ilink->is_col,rstart,rend,&ccis);
666:       MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
667:       ISDestroy(&ccis);
668:       ilink = ilink->next;
669:       ISComplement(ilink->is_col,rstart,rend,&ccis);
670:       MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
671:       ISDestroy(&ccis);

673:       /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
674:       MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);
675:       MatSetType(jac->schur,MATSCHURCOMPLEMENT);
676:       MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);
677:       PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
678:       /* 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? */
679:       MatSetOptionsPrefix(jac->schur,schurmatprefix);
680:       MatSetFromOptions(jac->schur);
681:       MatGetNullSpace(jac->pmat[1], &sp);
682:       if (sp) {
683:         MatSetNullSpace(jac->schur, sp);
684:       }

686:       PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);
687:       PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
688:       if (flg) {
689:         DM  dmInner;
690:         KSP kspInner;

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

698:         /* Set DM for new solver */
699:         KSPGetDM(jac->head->ksp, &dmInner);
700:         KSPSetDM(kspInner, dmInner);
701:         KSPSetDMActive(kspInner, PETSC_FALSE);
702:       } else {
703:          /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or
704:           * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact,
705:           * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for
706:           * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make
707:           * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used
708:           * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */
709:         KSPSetType(jac->head->ksp,KSPGMRES);
710:         MatSchurComplementSetKSP(jac->schur,jac->head->ksp);
711:       }
712:       KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0]);
713:       KSPSetFromOptions(jac->head->ksp);
714:       MatSetFromOptions(jac->schur);

716:       PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);
717:       PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);
718:       if (flg) {
719:         DM dmInner;

721:         PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);
722:         KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);
723:         KSPSetOptionsPrefix(jac->kspupper, schurprefix);
724:         KSPGetDM(jac->head->ksp, &dmInner);
725:         KSPSetDM(jac->kspupper, dmInner);
726:         KSPSetDMActive(jac->kspupper, PETSC_FALSE);
727:         KSPSetFromOptions(jac->kspupper);
728:         KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0]);
729:         VecDuplicate(jac->head->x, &jac->head->z);
730:       } else {
731:         jac->kspupper = jac->head->ksp;
732:         PetscObjectReference((PetscObject) jac->head->ksp);
733:       }

735:       if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
736:         MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);
737:       }
738:       KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);
739:       PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);
740:       PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);
741:       if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
742:         PC pcschur;
743:         KSPGetPC(jac->kspschur,&pcschur);
744:         PCSetType(pcschur,PCNONE);
745:         /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
746:       } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) {
747:         MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user);
748:       }
749:       KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
750:       KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);
751:       KSPSetOptionsPrefix(jac->kspschur,         Dprefix);
752:       /* propogate DM */
753:       {
754:         DM sdm;
755:         KSPGetDM(jac->head->next->ksp, &sdm);
756:         if (sdm) {
757:           KSPSetDM(jac->kspschur, sdm);
758:           KSPSetDMActive(jac->kspschur, PETSC_FALSE);
759:         }
760:       }
761:       /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
762:       /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */
763:       KSPSetFromOptions(jac->kspschur);
764:     }

766:     /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */
767:     PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);
768:     PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);
769:     if (!LSC_L) {PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);}
770:     if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);}
771:     PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);
772:     PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);
773:     if (!LSC_L) {PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);}
774:     if (LSC_L) {PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);}
775:   } else {
776:     /* set up the individual splits' PCs */
777:     i     = 0;
778:     ilink = jac->head;
779:     while (ilink) {
780:       KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i]);
781:       /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
782:       if (!jac->suboptionsset) {KSPSetFromOptions(ilink->ksp);}
783:       i++;
784:       ilink = ilink->next;
785:     }
786:   }

788:   jac->suboptionsset = PETSC_TRUE;
789:   return(0);
790: }

792: #define FieldSplitSplitSolveAdd(ilink,xx,yy) \
793:   (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
794:    VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \
795:    KSPSolve(ilink->ksp,ilink->x,ilink->y) || \
796:    VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \
797:    VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE))

801: static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y)
802: {
803:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
804:   PetscErrorCode    ierr;
805:   PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
806:   KSP               kspA   = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;

809:   switch (jac->schurfactorization) {
810:   case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
811:     /* [A00 0; 0 -S], positive definite, suitable for MINRES */
812:     VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
813:     VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
814:     VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
815:     KSPSolve(kspA,ilinkA->x,ilinkA->y);
816:     VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
817:     VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
818:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
819:     VecScale(ilinkD->y,-1.);
820:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
821:     VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
822:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
823:     break;
824:   case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
825:     /* [A00 0; A10 S], suitable for left preconditioning */
826:     VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
827:     VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
828:     KSPSolve(kspA,ilinkA->x,ilinkA->y);
829:     MatMult(jac->C,ilinkA->y,ilinkD->x);
830:     VecScale(ilinkD->x,-1.);
831:     VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
832:     VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
833:     VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
834:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
835:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
836:     VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
837:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
838:     break;
839:   case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
840:     /* [A00 A01; 0 S], suitable for right preconditioning */
841:     VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
842:     VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);
843:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
844:     MatMult(jac->B,ilinkD->y,ilinkA->x);
845:     VecScale(ilinkA->x,-1.);
846:     VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
847:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
848:     VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);
849:     KSPSolve(kspA,ilinkA->x,ilinkA->y);
850:     VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
851:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
852:     VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
853:     break;
854:   case PC_FIELDSPLIT_SCHUR_FACT_FULL:
855:     /* [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 */
856:     VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
857:     VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);
858:     KSPSolve(kspLower,ilinkA->x,ilinkA->y);
859:     MatMult(jac->C,ilinkA->y,ilinkD->x);
860:     VecScale(ilinkD->x,-1.0);
861:     VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);
862:     VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);

864:     KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);
865:     VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);
866:     VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);

868:     if (kspUpper == kspA) {
869:       MatMult(jac->B,ilinkD->y,ilinkA->y);
870:       VecAXPY(ilinkA->x,-1.0,ilinkA->y);
871:       KSPSolve(kspA,ilinkA->x,ilinkA->y);
872:     } else {
873:       KSPSolve(kspA,ilinkA->x,ilinkA->y);
874:       MatMult(jac->B,ilinkD->y,ilinkA->x);
875:       KSPSolve(kspUpper,ilinkA->x,ilinkA->z);
876:       VecAXPY(ilinkA->y,-1.0,ilinkA->z);
877:     }
878:     VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
879:     VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);
880:   }
881:   return(0);
882: }

886: static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y)
887: {
888:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
889:   PetscErrorCode    ierr;
890:   PC_FieldSplitLink ilink = jac->head;
891:   PetscInt          cnt,bs;

894:   if (jac->type == PC_COMPOSITE_ADDITIVE) {
895:     if (jac->defaultsplit) {
896:       VecGetBlockSize(x,&bs);
897:       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);
898:       VecGetBlockSize(y,&bs);
899:       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);
900:       VecStrideGatherAll(x,jac->x,INSERT_VALUES);
901:       while (ilink) {
902:         KSPSolve(ilink->ksp,ilink->x,ilink->y);
903:         ilink = ilink->next;
904:       }
905:       VecStrideScatterAll(jac->y,y,INSERT_VALUES);
906:     } else {
907:       VecSet(y,0.0);
908:       while (ilink) {
909:         FieldSplitSplitSolveAdd(ilink,x,y);
910:         ilink = ilink->next;
911:       }
912:     }
913:   } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
914:     if (!jac->w1) {
915:       VecDuplicate(x,&jac->w1);
916:       VecDuplicate(x,&jac->w2);
917:     }
918:     VecSet(y,0.0);
919:     FieldSplitSplitSolveAdd(ilink,x,y);
920:     cnt  = 1;
921:     while (ilink->next) {
922:       ilink = ilink->next;
923:       /* compute the residual only over the part of the vector needed */
924:       MatMult(jac->Afield[cnt++],y,ilink->x);
925:       VecScale(ilink->x,-1.0);
926:       VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
927:       VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
928:       KSPSolve(ilink->ksp,ilink->x,ilink->y);
929:       VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
930:       VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
931:     }
932:     if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
933:       cnt -= 2;
934:       while (ilink->previous) {
935:         ilink = ilink->previous;
936:         /* compute the residual only over the part of the vector needed */
937:         MatMult(jac->Afield[cnt--],y,ilink->x);
938:         VecScale(ilink->x,-1.0);
939:         VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
940:         VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);
941:         KSPSolve(ilink->ksp,ilink->x,ilink->y);
942:         VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
943:         VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);
944:       }
945:     }
946:   } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type);
947:   return(0);
948: }

950: #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \
951:   (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
952:    VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \
953:    KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \
954:    VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \
955:    VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE))

959: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y)
960: {
961:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
962:   PetscErrorCode    ierr;
963:   PC_FieldSplitLink ilink = jac->head;
964:   PetscInt          bs;

967:   if (jac->type == PC_COMPOSITE_ADDITIVE) {
968:     if (jac->defaultsplit) {
969:       VecGetBlockSize(x,&bs);
970:       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);
971:       VecGetBlockSize(y,&bs);
972:       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);
973:       VecStrideGatherAll(x,jac->x,INSERT_VALUES);
974:       while (ilink) {
975:         KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);
976:         ilink = ilink->next;
977:       }
978:       VecStrideScatterAll(jac->y,y,INSERT_VALUES);
979:     } else {
980:       VecSet(y,0.0);
981:       while (ilink) {
982:         FieldSplitSplitSolveAddTranspose(ilink,x,y);
983:         ilink = ilink->next;
984:       }
985:     }
986:   } else {
987:     if (!jac->w1) {
988:       VecDuplicate(x,&jac->w1);
989:       VecDuplicate(x,&jac->w2);
990:     }
991:     VecSet(y,0.0);
992:     if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
993:       FieldSplitSplitSolveAddTranspose(ilink,x,y);
994:       while (ilink->next) {
995:         ilink = ilink->next;
996:         MatMultTranspose(pc->mat,y,jac->w1);
997:         VecWAXPY(jac->w2,-1.0,jac->w1,x);
998:         FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
999:       }
1000:       while (ilink->previous) {
1001:         ilink = ilink->previous;
1002:         MatMultTranspose(pc->mat,y,jac->w1);
1003:         VecWAXPY(jac->w2,-1.0,jac->w1,x);
1004:         FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1005:       }
1006:     } else {
1007:       while (ilink->next) {   /* get to last entry in linked list */
1008:         ilink = ilink->next;
1009:       }
1010:       FieldSplitSplitSolveAddTranspose(ilink,x,y);
1011:       while (ilink->previous) {
1012:         ilink = ilink->previous;
1013:         MatMultTranspose(pc->mat,y,jac->w1);
1014:         VecWAXPY(jac->w2,-1.0,jac->w1,x);
1015:         FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
1016:       }
1017:     }
1018:   }
1019:   return(0);
1020: }

1024: static PetscErrorCode PCReset_FieldSplit(PC pc)
1025: {
1026:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1027:   PetscErrorCode    ierr;
1028:   PC_FieldSplitLink ilink = jac->head,next;

1031:   while (ilink) {
1032:     KSPReset(ilink->ksp);
1033:     VecDestroy(&ilink->x);
1034:     VecDestroy(&ilink->y);
1035:     VecDestroy(&ilink->z);
1036:     VecScatterDestroy(&ilink->sctx);
1037:     ISDestroy(&ilink->is);
1038:     ISDestroy(&ilink->is_col);
1039:     next  = ilink->next;
1040:     ilink = next;
1041:   }
1042:   PetscFree2(jac->x,jac->y);
1043:   if (jac->mat && jac->mat != jac->pmat) {
1044:     MatDestroyMatrices(jac->nsplits,&jac->mat);
1045:   } else if (jac->mat) {
1046:     jac->mat = NULL;
1047:   }
1048:   if (jac->pmat) {MatDestroyMatrices(jac->nsplits,&jac->pmat);}
1049:   if (jac->Afield) {MatDestroyMatrices(jac->nsplits,&jac->Afield);}
1050:   VecDestroy(&jac->w1);
1051:   VecDestroy(&jac->w2);
1052:   MatDestroy(&jac->schur);
1053:   MatDestroy(&jac->schurp);
1054:   MatDestroy(&jac->schur_user);
1055:   KSPDestroy(&jac->kspschur);
1056:   KSPDestroy(&jac->kspupper);
1057:   MatDestroy(&jac->B);
1058:   MatDestroy(&jac->C);
1059:   jac->reset = PETSC_TRUE;
1060:   return(0);
1061: }

1065: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
1066: {
1067:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1068:   PetscErrorCode    ierr;
1069:   PC_FieldSplitLink ilink = jac->head,next;

1072:   PCReset_FieldSplit(pc);
1073:   while (ilink) {
1074:     KSPDestroy(&ilink->ksp);
1075:     next  = ilink->next;
1076:     PetscFree(ilink->splitname);
1077:     PetscFree(ilink->fields);
1078:     PetscFree(ilink->fields_col);
1079:     PetscFree(ilink);
1080:     ilink = next;
1081:   }
1082:   PetscFree2(jac->x,jac->y);
1083:   PetscFree(pc->data);
1084:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);
1085:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);
1086:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);
1087:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);
1088:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);
1089:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);
1090:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);
1091:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);
1092:   return(0);
1093: }

1097: static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc)
1098: {
1099:   PetscErrorCode  ierr;
1100:   PetscInt        bs;
1101:   PetscBool       flg,stokes = PETSC_FALSE;
1102:   PC_FieldSplit   *jac = (PC_FieldSplit*)pc->data;
1103:   PCCompositeType ctype;

1106:   PetscOptionsHead("FieldSplit options");
1107:   PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);
1108:   PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);
1109:   if (flg) {
1110:     PCFieldSplitSetBlockSize(pc,bs);
1111:   }

1113:   PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);
1114:   if (stokes) {
1115:     PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);
1116:     jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF;
1117:   }

1119:   PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);
1120:   if (flg) {
1121:     PCFieldSplitSetType(pc,ctype);
1122:   }
1123:   /* Only setup fields once */
1124:   if ((jac->bs > 0) && (jac->nsplits == 0)) {
1125:     /* only allow user to set fields from command line if bs is already known.
1126:        otherwise user can set them in PCFieldSplitSetDefaults() */
1127:     PCFieldSplitSetRuntimeSplits_Private(pc);
1128:     if (jac->splitdefined) {PetscInfo(pc,"Splits defined using the options database\n");}
1129:   }
1130:   if (jac->type == PC_COMPOSITE_SCHUR) {
1131:     PetscOptionsGetEnum(((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);
1132:     if (flg) {PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");}
1133:     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);
1134:     PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);
1135:   }
1136:   PetscOptionsTail();
1137:   return(0);
1138: }

1140: /*------------------------------------------------------------------------------------*/

1144: static PetscErrorCode  PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1145: {
1146:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1147:   PetscErrorCode    ierr;
1148:   PC_FieldSplitLink ilink,next = jac->head;
1149:   char              prefix[128];
1150:   PetscInt          i;

1153:   if (jac->splitdefined) {
1154:     PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1155:     return(0);
1156:   }
1157:   for (i=0; i<n; i++) {
1158:     if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs);
1159:     if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]);
1160:   }
1161:   PetscNew(&ilink);
1162:   if (splitname) {
1163:     PetscStrallocpy(splitname,&ilink->splitname);
1164:   } else {
1165:     PetscMalloc1(3,&ilink->splitname);
1166:     PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);
1167:   }
1168:   PetscMalloc1(n,&ilink->fields);
1169:   PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));
1170:   PetscMalloc1(n,&ilink->fields_col);
1171:   PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));

1173:   ilink->nfields = n;
1174:   ilink->next    = NULL;
1175:   KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1176:   PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1177:   KSPSetType(ilink->ksp,KSPPREONLY);
1178:   PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);

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

1183:   if (!next) {
1184:     jac->head       = ilink;
1185:     ilink->previous = NULL;
1186:   } else {
1187:     while (next->next) {
1188:       next = next->next;
1189:     }
1190:     next->next      = ilink;
1191:     ilink->previous = next;
1192:   }
1193:   jac->nsplits++;
1194:   return(0);
1195: }

1199: static PetscErrorCode  PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp)
1200: {
1201:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

1205:   PetscMalloc1(jac->nsplits,subksp);
1206:   MatSchurComplementGetKSP(jac->schur,*subksp);

1208:   (*subksp)[1] = jac->kspschur;
1209:   if (n) *n = jac->nsplits;
1210:   return(0);
1211: }

1215: static PetscErrorCode  PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
1216: {
1217:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1218:   PetscErrorCode    ierr;
1219:   PetscInt          cnt   = 0;
1220:   PC_FieldSplitLink ilink = jac->head;

1223:   PetscMalloc1(jac->nsplits,subksp);
1224:   while (ilink) {
1225:     (*subksp)[cnt++] = ilink->ksp;
1226:     ilink            = ilink->next;
1227:   }
1228:   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);
1229:   if (n) *n = jac->nsplits;
1230:   return(0);
1231: }

1235: static PetscErrorCode  PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is)
1236: {
1237:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
1238:   PetscErrorCode    ierr;
1239:   PC_FieldSplitLink ilink, next = jac->head;
1240:   char              prefix[128];

1243:   if (jac->splitdefined) {
1244:     PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);
1245:     return(0);
1246:   }
1247:   PetscNew(&ilink);
1248:   if (splitname) {
1249:     PetscStrallocpy(splitname,&ilink->splitname);
1250:   } else {
1251:     PetscMalloc1(8,&ilink->splitname);
1252:     PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);
1253:   }
1254:   PetscObjectReference((PetscObject)is);
1255:   ISDestroy(&ilink->is);
1256:   ilink->is     = is;
1257:   PetscObjectReference((PetscObject)is);
1258:   ISDestroy(&ilink->is_col);
1259:   ilink->is_col = is;
1260:   ilink->next   = NULL;
1261:   KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);
1262:   PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);
1263:   KSPSetType(ilink->ksp,KSPPREONLY);
1264:   PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);

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

1269:   if (!next) {
1270:     jac->head       = ilink;
1271:     ilink->previous = NULL;
1272:   } else {
1273:     while (next->next) {
1274:       next = next->next;
1275:     }
1276:     next->next      = ilink;
1277:     ilink->previous = next;
1278:   }
1279:   jac->nsplits++;
1280:   return(0);
1281: }

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

1288:     Logically Collective on PC

1290:     Input Parameters:
1291: +   pc  - the preconditioner context
1292: .   splitname - name of this split, if NULL the number of the split is used
1293: .   n - the number of fields in this split
1294: -   fields - the fields in this split

1296:     Level: intermediate

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

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

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

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

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

1314: @*/
1315: PetscErrorCode  PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col)
1316: {

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

1330: /*@
1331:     PCFieldSplitSetIS - Sets the exact elements for field

1333:     Logically Collective on PC

1335:     Input Parameters:
1336: +   pc  - the preconditioner context
1337: .   splitname - name of this split, if NULL the number of the split is used
1338: -   is - the index set that defines the vector elements in this field


1341:     Notes:
1342:     Use PCFieldSplitSetFields(), for fields defined by strided types.

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

1347:     Level: intermediate

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

1351: @*/
1352: PetscErrorCode  PCFieldSplitSetIS(PC pc,const char splitname[],IS is)
1353: {

1360:   PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));
1361:   return(0);
1362: }

1366: /*@
1367:     PCFieldSplitGetIS - Retrieves the elements for a field as an IS

1369:     Logically Collective on PC

1371:     Input Parameters:
1372: +   pc  - the preconditioner context
1373: -   splitname - name of this split

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

1378:     Level: intermediate

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

1382: @*/
1383: PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is)
1384: {

1391:   {
1392:     PC_FieldSplit     *jac  = (PC_FieldSplit*) pc->data;
1393:     PC_FieldSplitLink ilink = jac->head;
1394:     PetscBool         found;

1396:     *is = NULL;
1397:     while (ilink) {
1398:       PetscStrcmp(ilink->splitname, splitname, &found);
1399:       if (found) {
1400:         *is = ilink->is;
1401:         break;
1402:       }
1403:       ilink = ilink->next;
1404:     }
1405:   }
1406:   return(0);
1407: }

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

1415:     Logically Collective on PC

1417:     Input Parameters:
1418: +   pc  - the preconditioner context
1419: -   bs - the block size

1421:     Level: intermediate

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

1425: @*/
1426: PetscErrorCode  PCFieldSplitSetBlockSize(PC pc,PetscInt bs)
1427: {

1433:   PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));
1434:   return(0);
1435: }

1439: /*@C
1440:    PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits

1442:    Collective on KSP

1444:    Input Parameter:
1445: .  pc - the preconditioner context

1447:    Output Parameters:
1448: +  n - the number of splits
1449: -  pc - the array of KSP contexts

1451:    Note:
1452:    After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user
1453:    (not the KSP just the array that contains them).

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

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


1462:    Level: advanced

1464: .seealso: PCFIELDSPLIT
1465: @*/
1466: PetscErrorCode  PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
1467: {

1473:   PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));
1474:   return(0);
1475: }

1479: /*@
1480:     PCFieldSplitSetSchurPre -  Indicates if the Schur complement is preconditioned by a preconditioner constructed by the
1481:       A11 matrix. Otherwise no preconditioner is used.

1483:     Collective on PC

1485:     Input Parameters:
1486: +   pc      - the preconditioner context
1487: .   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
1488: -   userpre - matrix to use for preconditioning, or NULL

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

1493:     Notes:
1494: $    If ptype is
1495: $        user then the preconditioner for the Schur complement is generated by the provided matrix (pre argument
1496: $             to this function).
1497: $        a11 then the preconditioner for the Schur complement is generated by the block diagonal part of the original
1498: $             matrix associated with the Schur complement (i.e. A11)
1499: $        full then the preconditioner uses the exact Schur complement (this is expensive)
1500: $        self the preconditioner for the Schur complement is generated from the Schur complement matrix itself:
1501: $             The only preconditioner that currently works directly with the Schur complement matrix object is the PCLSC
1502: $             preconditioner
1503: $        selfp then the preconditioning matrix is an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01
1504: $             This is only a good preconditioner when diag(A00) is a good preconditioner for A00.

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

1510:     Level: intermediate

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

1514: @*/
1515: PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
1516: {

1521:   PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));
1522:   return(0);
1523: }
1524: PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */

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

1532:     Logically Collective on PC

1534:     Input Parameters:
1535: .   pc      - the preconditioner context

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

1541:     Level: intermediate

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

1545: @*/
1546: PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
1547: {

1552:   PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));
1553:   return(0);
1554: }

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

1561:     Not collective

1563:     Input Parameter:
1564: .   pc      - the preconditioner context

1566:     Output Parameter:
1567: .   S       - the Schur complement matrix

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

1572:     Level: advanced

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

1576: @*/
1577: PetscErrorCode  PCFieldSplitSchurGetS(PC pc,Mat *S)
1578: {
1580:   const char*    t;
1581:   PetscBool      isfs;
1582:   PC_FieldSplit  *jac;

1586:   PetscObjectGetType((PetscObject)pc,&t);
1587:   PetscStrcmp(t,PCFIELDSPLIT,&isfs);
1588:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
1589:   jac = (PC_FieldSplit*)pc->data;
1590:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
1591:   if (S) *S = jac->schur;
1592:   return(0);
1593: }

1597: /*@
1598:     PCFieldSplitSchurRestoreS -  restores the MatSchurComplement object used by this PC

1600:     Not collective

1602:     Input Parameters:
1603: +   pc      - the preconditioner context
1604: .   S       - the Schur complement matrix

1606:     Level: advanced

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

1610: @*/
1611: PetscErrorCode  PCFieldSplitSchurRestoreS(PC pc,Mat *S)
1612: {
1614:   const char*    t;
1615:   PetscBool      isfs;
1616:   PC_FieldSplit  *jac;

1620:   PetscObjectGetType((PetscObject)pc,&t);
1621:   PetscStrcmp(t,PCFIELDSPLIT,&isfs);
1622:   if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t);
1623:   jac = (PC_FieldSplit*)pc->data;
1624:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type);
1625:   if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten");
1626:   return(0);
1627: }


1632: static PetscErrorCode  PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre)
1633: {
1634:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

1638:   jac->schurpre = ptype;
1639:   if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
1640:     MatDestroy(&jac->schur_user);
1641:     jac->schur_user = pre;
1642:     PetscObjectReference((PetscObject)jac->schur_user);
1643:   }
1644:   return(0);
1645: }

1649: static PetscErrorCode  PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre)
1650: {
1651:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

1654:   *ptype = jac->schurpre;
1655:   *pre   = jac->schur_user;
1656:   return(0);
1657: }

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

1664:     Collective on PC

1666:     Input Parameters:
1667: +   pc  - the preconditioner context
1668: -   ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default

1670:     Options Database:
1671: .     -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full


1674:     Level: intermediate

1676:     Notes:
1677:     The FULL factorization is

1679: $   (A   B)  = (1       0) (A   0) (1  Ainv*B)
1680: $   (C   D)    (C*Ainv  1) (0   S) (0     1  )

1682:     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,
1683:     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).

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

1690:     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
1691:     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).

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

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

1698: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType
1699: @*/
1700: PetscErrorCode  PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype)
1701: {

1706:   PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));
1707:   return(0);
1708: }

1712: static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype)
1713: {
1714:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

1717:   jac->schurfactorization = ftype;
1718:   return(0);
1719: }

1723: /*@C
1724:    PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement

1726:    Collective on KSP

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

1731:    Output Parameters:
1732: +  A00 - the (0,0) block
1733: .  A01 - the (0,1) block
1734: .  A10 - the (1,0) block
1735: -  A11 - the (1,1) block

1737:    Level: advanced

1739: .seealso: PCFIELDSPLIT
1740: @*/
1741: PetscErrorCode  PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11)
1742: {
1743:   PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;

1747:   if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach.");
1748:   if (A00) *A00 = jac->pmat[0];
1749:   if (A01) *A01 = jac->B;
1750:   if (A10) *A10 = jac->C;
1751:   if (A11) *A11 = jac->pmat[1];
1752:   return(0);
1753: }

1757: static PetscErrorCode  PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type)
1758: {
1759:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

1763:   jac->type = type;
1764:   if (type == PC_COMPOSITE_SCHUR) {
1765:     pc->ops->apply = PCApply_FieldSplit_Schur;
1766:     pc->ops->view  = PCView_FieldSplit_Schur;

1768:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);
1769:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);
1770:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);
1771:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);

1773:   } else {
1774:     pc->ops->apply = PCApply_FieldSplit;
1775:     pc->ops->view  = PCView_FieldSplit;

1777:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
1778:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);
1779:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);
1780:     PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);
1781:   }
1782:   return(0);
1783: }

1787: static PetscErrorCode  PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs)
1788: {
1789:   PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;

1792:   if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs);
1793:   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);
1794:   jac->bs = bs;
1795:   return(0);
1796: }

1800: /*@
1801:    PCFieldSplitSetType - Sets the type of fieldsplit preconditioner.

1803:    Collective on PC

1805:    Input Parameter:
1806: .  pc - the preconditioner context
1807: .  type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR

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

1812:    Level: Intermediate

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

1816: .seealso: PCCompositeSetType()

1818: @*/
1819: PetscErrorCode  PCFieldSplitSetType(PC pc,PCCompositeType type)
1820: {

1825:   PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));
1826:   return(0);
1827: }

1831: /*@
1832:   PCFieldSplitGetType - Gets the type of fieldsplit preconditioner.

1834:   Not collective

1836:   Input Parameter:
1837: . pc - the preconditioner context

1839:   Output Parameter:
1840: . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR

1842:   Level: Intermediate

1844: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
1845: .seealso: PCCompositeSetType()
1846: @*/
1847: PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
1848: {
1849:   PC_FieldSplit *jac = (PC_FieldSplit*) pc->data;

1854:   *type = jac->type;
1855:   return(0);
1856: }

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

1863:    Logically Collective

1865:    Input Parameters:
1866: +  pc   - the preconditioner context
1867: -  flg  - boolean indicating whether to use field splits defined by the DM

1869:    Options Database Key:
1870: .  -pc_fieldsplit_dm_splits

1872:    Level: Intermediate

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

1876: .seealso: PCFieldSplitGetDMSplits()

1878: @*/
1879: PetscErrorCode  PCFieldSplitSetDMSplits(PC pc,PetscBool flg)
1880: {
1881:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1882:   PetscBool      isfs;

1888:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1889:   if (isfs) {
1890:     jac->dm_splits = flg;
1891:   }
1892:   return(0);
1893: }


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

1901:    Logically Collective

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

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

1909:    Level: Intermediate

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

1913: .seealso: PCFieldSplitSetDMSplits()

1915: @*/
1916: PetscErrorCode  PCFieldSplitGetDMSplits(PC pc,PetscBool* flg)
1917: {
1918:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;
1919:   PetscBool      isfs;

1925:   PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);
1926:   if (isfs) {
1927:     if(flg) *flg = jac->dm_splits;
1928:   }
1929:   return(0);
1930: }

1932: /* -------------------------------------------------------------------------------------*/
1933: /*MC
1934:    PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
1935:                   fields or groups of fields. See the users manual section "Solving Block Matrices" for more details.

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

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

1943:    Level: intermediate

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

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

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

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

1966: $     For the Schur complement preconditioner if J = ( A00 A01 )
1967: $                                                    ( A10 A11 )
1968: $     the preconditioner using full factorization is
1969: $              ( I   -ksp(A00) A01 ) ( inv(A00)     0  ) (     I          0  )
1970: $              ( 0         I       ) (   0      ksp(S) ) ( -A10 ksp(A00)  I  )
1971:      where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_.  S is the Schur complement
1972: $              S = A11 - A10 ksp(A00) A01
1973:      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
1974:      in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0,
1975:      it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default
1976:      A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() to turn on or off this
1977:      option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc.
1978:      When option -fieldsplit_schur_precondition selfp is given, an approximation to S is assembled --
1979:      Sp = A11 - A10 inv(diag(A00)) A01, which has type AIJ and can be used with a variety of preconditioners
1980:      (e.g., -fieldsplit_1_pc_type asm).
1981:      The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above,
1982:      diag gives
1983: $              ( inv(A00)     0   )
1984: $              (   0      -ksp(S) )
1985:      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
1986: $              (  A00   0 )
1987: $              (  A10   S )
1988:      where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of
1989: $              ( A00 A01 )
1990: $              (  0   S  )
1991:      where again the inverses of A00 and S are applied using KSPs.

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

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

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

2003:    Concepts: physics based preconditioners, block preconditioners

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

2011: PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
2012: {
2014:   PC_FieldSplit  *jac;

2017:   PetscNewLog(pc,&jac);

2019:   jac->bs                 = -1;
2020:   jac->nsplits            = 0;
2021:   jac->type               = PC_COMPOSITE_MULTIPLICATIVE;
2022:   jac->schurpre           = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
2023:   jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
2024:   jac->dm_splits          = PETSC_TRUE;

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

2028:   pc->ops->apply           = PCApply_FieldSplit;
2029:   pc->ops->applytranspose  = PCApplyTranspose_FieldSplit;
2030:   pc->ops->setup           = PCSetUp_FieldSplit;
2031:   pc->ops->reset           = PCReset_FieldSplit;
2032:   pc->ops->destroy         = PCDestroy_FieldSplit;
2033:   pc->ops->setfromoptions  = PCSetFromOptions_FieldSplit;
2034:   pc->ops->view            = PCView_FieldSplit;
2035:   pc->ops->applyrichardson = 0;

2037:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);
2038:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);
2039:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);
2040:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);
2041:   PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);
2042:   return(0);
2043: }