Actual source code: fieldsplit.c

  1: #define PETSCKSP_DLL

  3: /*

  5: */
 6:  #include private/pcimpl.h

  8: typedef struct _PC_FieldSplitLink *PC_FieldSplitLink;
  9: struct _PC_FieldSplitLink {
 10:   KSP               ksp;
 11:   Vec               x,y;
 12:   PetscInt          nfields;
 13:   PetscInt          *fields;
 14:   VecScatter        sctx;
 15:   PC_FieldSplitLink next,previous;
 16: };

 18: typedef struct {
 19:   PCCompositeType   type;              /* additive or multiplicative */
 20:   PetscTruth        defaultsplit;
 21:   PetscInt          bs;
 22:   PetscInt          nsplits,*csize;
 23:   Vec               *x,*y,w1,w2;
 24:   Mat               *pmat;
 25:   IS                *is,*cis;
 26:   PC_FieldSplitLink head;
 27: } PC_FieldSplit;

 31: static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer)
 32: {
 33:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
 34:   PetscErrorCode    ierr;
 35:   PetscTruth        iascii;
 36:   PetscInt          i,j;
 37:   PC_FieldSplitLink ilink = jac->head;

 40:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
 41:   if (iascii) {
 42:     PetscViewerASCIIPrintf(viewer,"  FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);
 43:     PetscViewerASCIIPrintf(viewer,"  Solver info for each split is in the following KSP objects:\n");
 44:     PetscViewerASCIIPushTab(viewer);
 45:     for (i=0; i<jac->nsplits; i++) {
 46:       PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);
 47:       PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
 48:       for (j=0; j<ilink->nfields; j++) {
 49:         if (j > 0) {
 50:           PetscViewerASCIIPrintf(viewer,",");
 51:         }
 52:         PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);
 53:       }
 54:       PetscViewerASCIIPrintf(viewer,"\n");
 55:       PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
 56:       KSPView(ilink->ksp,viewer);
 57:       ilink = ilink->next;
 58:     }
 59:     PetscViewerASCIIPopTab(viewer);
 60:   } else {
 61:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCFieldSplit",((PetscObject)viewer)->type_name);
 62:   }
 63:   return(0);
 64: }

 68: static PetscErrorCode PCFieldSplitSetDefaults(PC pc)
 69: {
 70:   PC_FieldSplit     *jac  = (PC_FieldSplit*)pc->data;
 71:   PetscErrorCode    ierr;
 72:   PC_FieldSplitLink ilink = jac->head;
 73:   PetscInt          i;
 74:   PetscTruth        flg = PETSC_FALSE,*fields;

 77:   PetscOptionsGetTruth(pc->prefix,"-pc_fieldsplit_default",&flg,PETSC_NULL);
 78:   if (!ilink || flg) {
 79:     PetscInfo(pc,"Using default splitting of fields\n");
 80:     if (jac->bs <= 0) {
 81:       MatGetBlockSize(pc->pmat,&jac->bs);
 82:     }
 83:     PetscMalloc(jac->bs*sizeof(PetscTruth),&fields);
 84:     PetscMemzero(fields,jac->bs*sizeof(PetscTruth));
 85:     while (ilink) {
 86:       for (i=0; i<ilink->nfields; i++) {
 87:         fields[ilink->fields[i]] = PETSC_TRUE;
 88:       }
 89:       ilink = ilink->next;
 90:     }
 91:     jac->defaultsplit = PETSC_TRUE;
 92:     for (i=0; i<jac->bs; i++) {
 93:       if (!fields[i]) {
 94:         PCFieldSplitSetFields(pc,1,&i);
 95:       } else {
 96:         jac->defaultsplit = PETSC_FALSE;
 97:       }
 98:     }
 99:   }
100:   return(0);
101: }


106: static PetscErrorCode PCSetUp_FieldSplit(PC pc)
107: {
108:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
109:   PetscErrorCode    ierr;
110:   PC_FieldSplitLink ilink;
111:   PetscInt          i,nsplit,ccsize;
112:   MatStructure      flag = pc->flag;
113:   PetscTruth        sorted;

116:   PCFieldSplitSetDefaults(pc);
117:   nsplit = jac->nsplits;
118:   ilink  = jac->head;

120:   /* get the matrices for each split */
121:   if (!jac->is) {
122:     PetscInt rstart,rend,nslots,bs;

124:     bs     = jac->bs;
125:     MatGetOwnershipRange(pc->pmat,&rstart,&rend);
126:     MatGetLocalSize(pc->pmat,PETSC_NULL,&ccsize);
127:     nslots = (rend - rstart)/bs;
128:     PetscMalloc(nsplit*sizeof(IS),&jac->is);
129:     PetscMalloc(nsplit*sizeof(IS),&jac->cis);
130:     PetscMalloc(nsplit*sizeof(PetscInt),&jac->csize);
131:     for (i=0; i<nsplit; i++) {
132:       if (jac->defaultsplit) {
133:         ISCreateStride(pc->comm,nslots,rstart+i,nsplit,&jac->is[i]);
134:         jac->csize[i] = ccsize/nsplit;
135:       } else {
136:         if (ilink->nfields > 1) {
137:           PetscInt   *ii,j,k,nfields = ilink->nfields,*fields = ilink->fields;
138:           PetscMalloc(ilink->nfields*nslots*sizeof(PetscInt),&ii);
139:           for (j=0; j<nslots; j++) {
140:             for (k=0; k<nfields; k++) {
141:               ii[nfields*j + k] = rstart + bs*j + fields[k];
142:             }
143:           }
144:           ISCreateGeneral(pc->comm,nslots*nfields,ii,&jac->is[i]);
145:           PetscFree(ii);
146:         } else {
147:           ISCreateStride(pc->comm,nslots,ilink->fields[0],bs,&jac->is[i]);
148:         }
149:         jac->csize[i] = (ccsize/bs)*ilink->nfields;
150:         ISSorted(jac->is[i],&sorted);
151:         if (!sorted) SETERRQ(PETSC_ERR_USER,"Fields must be sorted when creating split");
152:         ilink = ilink->next;
153:       }
154:       ISAllGather(jac->is[i],&jac->cis[i]);
155:     }
156:   }
157: 
158:   if (!jac->pmat) {
159:     PetscMalloc(nsplit*sizeof(Mat),&jac->pmat);
160:     for (i=0; i<nsplit; i++) {
161:       MatGetSubMatrix(pc->pmat,jac->is[i],jac->cis[i],jac->csize[i],MAT_INITIAL_MATRIX,&jac->pmat[i]);
162:     }
163:   } else {
164:     for (i=0; i<nsplit; i++) {
165:       MatGetSubMatrix(pc->pmat,jac->is[i],jac->cis[i],jac->csize[i],MAT_REUSE_MATRIX,&jac->pmat[i]);
166:     }
167:   }

169:   /* set up the individual PCs */
170:   i    = 0;
171:   ilink = jac->head;
172:   while (ilink) {
173:     KSPSetOperators(ilink->ksp,jac->pmat[i],jac->pmat[i],flag);
174:     KSPSetFromOptions(ilink->ksp);
175:     KSPSetUp(ilink->ksp);
176:     i++;
177:     ilink = ilink->next;
178:   }
179: 
180:   /* create work vectors for each split */
181:   if (!jac->x) {
182:     Vec xtmp;
183:     PetscMalloc2(nsplit,Vec,&jac->x,nsplit,Vec,&jac->y);
184:     ilink = jac->head;
185:     for (i=0; i<nsplit; i++) {
186:       Vec *vl,*vr;

188:       KSPGetVecs(ilink->ksp,1,&vr,1,&vl);
189:       ilink->x  = *vr;
190:       ilink->y  = *vl;
191:       PetscFree(vr);
192:       PetscFree(vl);
193:       jac->x[i] = ilink->x;
194:       jac->y[i] = ilink->y;
195:       ilink     = ilink->next;
196:     }
197:     /* compute scatter contexts needed by multiplicative versions and non-default splits */
198: 
199:     ilink = jac->head;
200:     MatGetVecs(pc->pmat,&xtmp,PETSC_NULL);
201:     for (i=0; i<nsplit; i++) {
202:       VecScatterCreate(xtmp,jac->is[i],jac->x[i],PETSC_NULL,&ilink->sctx);
203:       ilink = ilink->next;
204:     }
205:     VecDestroy(xtmp);
206:   }
207:   return(0);
208: }

210: #define FieldSplitSplitSolveAdd(ilink,xx,yy) \
211:     (VecScatterBegin(xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD,ilink->sctx) || \
212:      VecScatterEnd(xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD,ilink->sctx) || \
213:      KSPSolve(ilink->ksp,ilink->x,ilink->y) || \
214:      VecScatterBegin(ilink->y,yy,ADD_VALUES,SCATTER_REVERSE,ilink->sctx) || \
215:      VecScatterEnd(ilink->y,yy,ADD_VALUES,SCATTER_REVERSE,ilink->sctx))

219: static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y)
220: {
221:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
222:   PetscErrorCode    ierr;
223:   PC_FieldSplitLink ilink = jac->head;
224:   PetscInt          bs;

227:   CHKMEMQ;
228:   VecGetBlockSize(x,&bs);
229:   VecSetBlockSize(x,jac->bs);
230:   VecSetBlockSize(y,jac->bs);

232:   if (jac->type == PC_COMPOSITE_ADDITIVE) {
233:     if (jac->defaultsplit) {
234:       VecStrideGatherAll(x,jac->x,INSERT_VALUES);
235:       while (ilink) {
236:         KSPSolve(ilink->ksp,ilink->x,ilink->y);
237:         ilink = ilink->next;
238:       }
239:       VecStrideScatterAll(jac->y,y,INSERT_VALUES);
240:     } else {
241:       VecSet(y,0.0);
242:       while (ilink) {
243:         FieldSplitSplitSolveAdd(ilink,x,y);
244:         ilink = ilink->next;
245:       }
246:     }
247:   } else {
248:     if (!jac->w1) {
249:       VecDuplicate(x,&jac->w1);
250:       VecDuplicate(x,&jac->w2);
251:     }
252:     VecSet(y,0.0);
253:     FieldSplitSplitSolveAdd(ilink,x,y);
254:     while (ilink->next) {
255:       ilink = ilink->next;
256:       MatMult(pc->pmat,y,jac->w1);
257:       VecWAXPY(jac->w2,-1.0,jac->w1,x);
258:       FieldSplitSplitSolveAdd(ilink,jac->w2,y);
259:     }
260:     if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
261:       while (ilink->previous) {
262:         ilink = ilink->previous;
263:         MatMult(pc->pmat,y,jac->w1);
264:         VecWAXPY(jac->w2,-1.0,jac->w1,x);
265:         FieldSplitSplitSolveAdd(ilink,jac->w2,y);
266:       }
267:     }
268:   }
269:   CHKMEMQ;
270:   return(0);
271: }

273: #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \
274:     (VecScatterBegin(xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD,ilink->sctx) || \
275:      VecScatterEnd(xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD,ilink->sctx) || \
276:      KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \
277:      VecScatterBegin(ilink->x,yy,ADD_VALUES,SCATTER_REVERSE,ilink->sctx) || \
278:      VecScatterEnd(ilink->x,yy,ADD_VALUES,SCATTER_REVERSE,ilink->sctx))

282: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y)
283: {
284:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
285:   PetscErrorCode    ierr;
286:   PC_FieldSplitLink ilink = jac->head;
287:   PetscInt          bs;

290:   CHKMEMQ;
291:   VecGetBlockSize(x,&bs);
292:   VecSetBlockSize(x,jac->bs);
293:   VecSetBlockSize(y,jac->bs);

295:   if (jac->type == PC_COMPOSITE_ADDITIVE) {
296:     if (jac->defaultsplit) {
297:       VecStrideGatherAll(x,jac->x,INSERT_VALUES);
298:       while (ilink) {
299:         KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);
300:         ilink = ilink->next;
301:       }
302:       VecStrideScatterAll(jac->y,y,INSERT_VALUES);
303:     } else {
304:       VecSet(y,0.0);
305:       while (ilink) {
306:         FieldSplitSplitSolveAddTranspose(ilink,x,y);
307:         ilink = ilink->next;
308:       }
309:     }
310:   } else {
311:     if (!jac->w1) {
312:       VecDuplicate(x,&jac->w1);
313:       VecDuplicate(x,&jac->w2);
314:     }
315:     VecSet(y,0.0);
316:     if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
317:       FieldSplitSplitSolveAddTranspose(ilink,x,y);
318:       while (ilink->next) {
319:         ilink = ilink->next;
320:         MatMultTranspose(pc->pmat,y,jac->w1);
321:         VecWAXPY(jac->w2,-1.0,jac->w1,x);
322:         FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
323:       }
324:       while (ilink->previous) {
325:         ilink = ilink->previous;
326:         MatMultTranspose(pc->pmat,y,jac->w1);
327:         VecWAXPY(jac->w2,-1.0,jac->w1,x);
328:         FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
329:       }
330:     } else {
331:       while (ilink->next) {   /* get to last entry in linked list */
332:         ilink = ilink->next;
333:       }
334:       FieldSplitSplitSolveAddTranspose(ilink,x,y);
335:       while (ilink->previous) {
336:         ilink = ilink->previous;
337:         MatMultTranspose(pc->pmat,y,jac->w1);
338:         VecWAXPY(jac->w2,-1.0,jac->w1,x);
339:         FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);
340:       }
341:     }
342:   }
343:   CHKMEMQ;
344:   return(0);
345: }

349: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
350: {
351:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
352:   PetscErrorCode    ierr;
353:   PC_FieldSplitLink ilink = jac->head,next;

356:   while (ilink) {
357:     KSPDestroy(ilink->ksp);
358:     if (ilink->x) {VecDestroy(ilink->x);}
359:     if (ilink->y) {VecDestroy(ilink->y);}
360:     if (ilink->sctx) {VecScatterDestroy(ilink->sctx);}
361:     next = ilink->next;
362:     PetscFree2(ilink,ilink->fields);
363:     ilink = next;
364:   }
365:   PetscFree2(jac->x,jac->y);
366:   if (jac->pmat) {MatDestroyMatrices(jac->nsplits,&jac->pmat);}
367:   if (jac->is) {
368:     PetscInt i;
369:     for (i=0; i<jac->nsplits; i++) {ISDestroy(jac->is[i]);}
370:     PetscFree(jac->is);
371:   }
372:   if (jac->cis) {
373:     PetscInt i;
374:     for (i=0; i<jac->nsplits; i++) {ISDestroy(jac->cis[i]);}
375:     PetscFree(jac->cis);
376:   }
377:   if (jac->w1) {VecDestroy(jac->w1);}
378:   if (jac->w2) {VecDestroy(jac->w2);}
379:   PetscFree(jac->csize);
380:   PetscFree(jac);
381:   return(0);
382: }

386: static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc)
387: {
389:   PetscInt       i = 0,nfields,*fields,bs;
390:   PetscTruth     flg;
391:   char           optionname[128];
392:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

395:   PetscOptionsHead("FieldSplit options");
396:   PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);
397:   if (flg) {
398:     PCFieldSplitSetBlockSize(pc,bs);
399:   }
400:   if (jac->bs <= 0) {
401:     PCFieldSplitSetBlockSize(pc,1);
402:   }
403:   PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&jac->type,&flg);
404:   PetscMalloc(jac->bs*sizeof(PetscInt),&fields);
405:   while (PETSC_TRUE) {
406:     sprintf(optionname,"-pc_fieldsplit_%d_fields",(int)i++);
407:     nfields = jac->bs;
408:     PetscOptionsIntArray(optionname,"Fields in this split","PCFieldSplitSetFields",fields,&nfields,&flg);
409:     if (!flg) break;
410:     if (!nfields) SETERRQ(PETSC_ERR_USER,"Cannot list zero fields");
411:     PCFieldSplitSetFields(pc,nfields,fields);
412:   }
413:   PetscFree(fields);
414:   PetscOptionsTail();
415:   return(0);
416: }

418: /*------------------------------------------------------------------------------------*/

423: PetscErrorCode  PCFieldSplitSetFields_FieldSplit(PC pc,PetscInt n,PetscInt *fields)
424: {
425:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
426:   PetscErrorCode    ierr;
427:   PC_FieldSplitLink ilink,next = jac->head;
428:   char              prefix[128];
429:   PetscInt          i;

432:   if (n <= 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative number of fields requested");
433:   for (i=0; i<n; i++) {
434:     if (fields[i] >= jac->bs) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs);
435:     if (fields[i] < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]);
436:   }
437:   PetscMalloc2(1,struct _PC_FieldSplitLink,&ilink,n,PetscInt,&ilink->fields);
438:   PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));
439:   ilink->nfields = n;
440:   ilink->next    = PETSC_NULL;
441:   KSPCreate(pc->comm,&ilink->ksp);
442:   KSPSetType(ilink->ksp,KSPPREONLY);

444:   if (pc->prefix) {
445:     sprintf(prefix,"%sfieldsplit_%d_",pc->prefix,(int)jac->nsplits);
446:   } else {
447:     sprintf(prefix,"fieldsplit_%d_",(int)jac->nsplits);
448:   }
449:   KSPSetOptionsPrefix(ilink->ksp,prefix);

451:   if (!next) {
452:     jac->head       = ilink;
453:     ilink->previous = PETSC_NULL;
454:   } else {
455:     while (next->next) {
456:       next = next->next;
457:     }
458:     next->next      = ilink;
459:     ilink->previous = next;
460:   }
461:   jac->nsplits++;
462:   return(0);
463: }


470: PetscErrorCode  PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
471: {
472:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
473:   PetscErrorCode    ierr;
474:   PetscInt          cnt = 0;
475:   PC_FieldSplitLink ilink = jac->head;

478:   PetscMalloc(jac->nsplits*sizeof(KSP*),subksp);
479:   while (ilink) {
480:     (*subksp)[cnt++] = ilink->ksp;
481:     ilink = ilink->next;
482:   }
483:   if (cnt != jac->nsplits) SETERRQ2(PETSC_ERR_PLIB,"Corrupt PCFIELDSPLIT object: number splits in linked list %D in object %D",cnt,jac->nsplits);
484:   *n = jac->nsplits;
485:   return(0);
486: }

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

494:     Collective on PC

496:     Input Parameters:
497: +   pc  - the preconditioner context
498: .   n - the number of fields in this split
499: .   fields - the fields in this split

501:     Level: intermediate

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

505: @*/
506: PetscErrorCode  PCFieldSplitSetFields(PC pc,PetscInt n, PetscInt *fields)
507: {
508:   PetscErrorCode ierr,(*f)(PC,PetscInt,PetscInt *);

512:   PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetFields_C",(void (**)(void))&f);
513:   if (f) {
514:     (*f)(pc,n,fields);
515:   }
516:   return(0);
517: }

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

525:     Collective on PC

527:     Input Parameters:
528: +   pc  - the preconditioner context
529: -   bs - the block size

531:     Level: intermediate

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

535: @*/
536: PetscErrorCode  PCFieldSplitSetBlockSize(PC pc,PetscInt bs)
537: {
538:   PetscErrorCode ierr,(*f)(PC,PetscInt);

542:   PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",(void (**)(void))&f);
543:   if (f) {
544:     (*f)(pc,bs);
545:   }
546:   return(0);
547: }

551: /*@C
552:    PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits
553:    
554:    Collective on KSP

556:    Input Parameter:
557: .  pc - the preconditioner context

559:    Output Parameters:
560: +  n - the number of split
561: -  pc - the array of KSP contexts

563:    Note:  
564:    After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed

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

568:    Level: advanced

570: .seealso: PCFIELDSPLIT
571: @*/
572: PetscErrorCode  PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
573: {
574:   PetscErrorCode ierr,(*f)(PC,PetscInt*,KSP **);

579:   PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",(void (**)(void))&f);
580:   if (f) {
581:     (*f)(pc,n,subksp);
582:   } else {
583:     SETERRQ(PETSC_ERR_ARG_WRONG,"Cannot get subksp for this type of PC");
584:   }
585:   return(0);
586: }

591: PetscErrorCode  PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type)
592: {
593:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

596:   jac->type = type;
597:   return(0);
598: }

604: PetscErrorCode  PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs)
605: {
606:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

609:   if (bs < 1) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs);
610:   if (jac->bs > 0 && jac->bs != bs) SETERRQ2(PETSC_ERR_ARG_WRONGSTATE,"Cannot change fieldsplit blocksize from %D to %D after it has been set",jac->bs,bs);
611:   jac->bs = bs;
612:   return(0);
613: }

618: /*@C
619:    PCFieldSplitSetType - Sets the type of fieldsplit preconditioner.
620:    
621:    Collective on PC

623:    Input Parameter:
624: .  pc - the preconditioner context
625: .  type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE

627:    Options Database Key:
628: .  -pc_fieldsplit_type <type: one of multiplicative, additive, special> - Sets fieldsplit preconditioner type

630:    Level: Developer

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

634: .seealso: PCCompositeSetType()

636: @*/
637: PetscErrorCode  PCFieldSplitSetType(PC pc,PCCompositeType type)
638: {
639:   PetscErrorCode ierr,(*f)(PC,PCCompositeType);

643:   PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetType_C",(void (**)(void))&f);
644:   if (f) {
645:     (*f)(pc,type);
646:   }
647:   return(0);
648: }

650: /* -------------------------------------------------------------------------------------*/
651: /*MC
652:    PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
653:                   fields or groups of fields


656:      To set options on the solvers for each block append -sub_ to all the PC
657:         options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1
658:         
659:      To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP()
660:          and set the options directly on the resulting KSP object

662:    Level: intermediate

664:    Options Database Keys:
665: +   -pc_splitfield_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split
666: .   -pc_splitfield_default - automatically add any fields to additional splits that have not
667:                               been supplied explicitly by -pc_splitfield_%d_fields
668: .   -pc_splitfield_block_size <bs> - size of block that defines fields (i.e. there are bs fields)
669: -   -pc_splitfield_type <additive,multiplicative>

671:    Concepts: physics based preconditioners

673: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
674:            PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(),PCFieldSplitSetType()
675: M*/

680: PetscErrorCode  PCCreate_FieldSplit(PC pc)
681: {
683:   PC_FieldSplit  *jac;

686:   PetscNew(PC_FieldSplit,&jac);
687:   PetscLogObjectMemory(pc,sizeof(PC_FieldSplit));
688:   jac->bs        = -1;
689:   jac->nsplits   = 0;
690:   jac->type      = PC_COMPOSITE_ADDITIVE;

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

694:   pc->ops->apply             = PCApply_FieldSplit;
695:   pc->ops->applytranspose    = PCApplyTranspose_FieldSplit;
696:   pc->ops->setup             = PCSetUp_FieldSplit;
697:   pc->ops->destroy           = PCDestroy_FieldSplit;
698:   pc->ops->setfromoptions    = PCSetFromOptions_FieldSplit;
699:   pc->ops->view              = PCView_FieldSplit;
700:   pc->ops->applyrichardson   = 0;

702:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit",
703:                     PCFieldSplitGetSubKSP_FieldSplit);
704:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","PCFieldSplitSetFields_FieldSplit",
705:                     PCFieldSplitSetFields_FieldSplit);
706:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","PCFieldSplitSetType_FieldSplit",
707:                     PCFieldSplitSetType_FieldSplit);
708:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","PCFieldSplitSetBlockSize_FieldSplit",
709:                     PCFieldSplitSetBlockSize_FieldSplit);
710:   return(0);
711: }