Actual source code: composite.c

  1: #define PETSCKSP_DLL

  3: /*
  4:       Defines a preconditioner that can consist of a collection of PCs
  5: */
 6:  #include private/pcimpl.h
 7:  #include petscksp.h

  9: typedef struct _PC_CompositeLink *PC_CompositeLink;
 10: struct _PC_CompositeLink {
 11:   PC               pc;
 12:   PC_CompositeLink next;
 13:   PC_CompositeLink previous;
 14: };
 15: 
 16: typedef struct {
 17:   PC_CompositeLink head;
 18:   PCCompositeType  type;
 19:   Vec              work1;
 20:   Vec              work2;
 21:   PetscScalar      alpha;
 22:   PetscTruth       use_true_matrix;
 23: } PC_Composite;

 27: static PetscErrorCode PCApply_Composite_Multiplicative(PC pc,Vec x,Vec y)
 28: {
 29:   PetscErrorCode   ierr;
 30:   PC_Composite     *jac = (PC_Composite*)pc->data;
 31:   PC_CompositeLink next = jac->head;
 32:   Mat              mat = pc->pmat;

 35:   if (!next) {
 36:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC() or -pc_composite_pcs");
 37:   }
 38:   if (next->next && !jac->work2) { /* allocate second work vector */
 39:     VecDuplicate(jac->work1,&jac->work2);
 40:   }
 41:   PCApply(next->pc,x,y);
 42:   if (jac->use_true_matrix) mat = pc->mat;
 43:   while (next->next) {
 44:     next = next->next;
 45:     MatMult(mat,y,jac->work1);
 46:     VecWAXPY(jac->work2,-1.0,jac->work1,x);
 47:     VecSet(jac->work1,0.0);  /* zero since some PC's may not set all entries in the result */
 48:     PCApply(next->pc,jac->work2,jac->work1);
 49:     VecAXPY(y,1.0,jac->work1);
 50:   }
 51:   if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
 52:     while (next->previous) {
 53:       next = next->previous;
 54:       MatMult(mat,y,jac->work1);
 55:       VecWAXPY(jac->work2,-1.0,jac->work1,x);
 56:       VecSet(jac->work1,0.0);  /* zero since some PC's may not set all entries in the result */
 57:       PCApply(next->pc,jac->work2,jac->work1);
 58:       VecAXPY(y,1.0,jac->work1);
 59:     }
 60:   }
 61:   return(0);
 62: }

 64: /*
 65:     This is very special for a matrix of the form alpha I + R + S
 66: where first preconditioner is built from alpha I + S and second from
 67: alpha I + R
 68: */
 71: static PetscErrorCode PCApply_Composite_Special(PC pc,Vec x,Vec y)
 72: {
 73:   PetscErrorCode   ierr;
 74:   PC_Composite     *jac = (PC_Composite*)pc->data;
 75:   PC_CompositeLink next = jac->head;

 78:   if (!next) {
 79:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC() or -pc_composite_pcs");
 80:   }
 81:   if (!next->next || next->next->next) {
 82:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Special composite preconditioners requires exactly two PCs");
 83:   }

 85:   PCApply(next->pc,x,jac->work1);
 86:   PCApply(next->next->pc,jac->work1,y);
 87:   return(0);
 88: }

 92: static PetscErrorCode PCApply_Composite_Additive(PC pc,Vec x,Vec y)
 93: {
 94:   PetscErrorCode   ierr;
 95:   PC_Composite     *jac = (PC_Composite*)pc->data;
 96:   PC_CompositeLink next = jac->head;

 99:   if (!next) {
100:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC() or -pc_composite_pcs");
101:   }
102:   PCApply(next->pc,x,y);
103:   while (next->next) {
104:     next = next->next;
105:     VecSet(jac->work1,0.0);  /* zero since some PC's may not set all entries in the result */
106:     PCApply(next->pc,x,jac->work1);
107:     VecAXPY(y,1.0,jac->work1);
108:   }
109:   return(0);
110: }

114: static PetscErrorCode PCSetUp_Composite(PC pc)
115: {
116:   PetscErrorCode   ierr;
117:   PC_Composite     *jac = (PC_Composite*)pc->data;
118:   PC_CompositeLink next = jac->head;

121:   if (!jac->work1) {
122:    MatGetVecs(pc->pmat,&jac->work1,0);
123:   }
124:   while (next) {
125:     PCSetOperators(next->pc,pc->mat,pc->pmat,pc->flag);
126:     next = next->next;
127:   }
128:   return(0);
129: }

133: static PetscErrorCode PCDestroy_Composite(PC pc)
134: {
135:   PC_Composite     *jac = (PC_Composite*)pc->data;
136:   PetscErrorCode   ierr;
137:   PC_CompositeLink next = jac->head;

140:   while (next) {
141:     PCDestroy(next->pc);
142:     next = next->next;
143:   }

145:   if (jac->work1) {VecDestroy(jac->work1);}
146:   if (jac->work2) {VecDestroy(jac->work2);}
147:   PetscFree(jac);
148:   return(0);
149: }

153: static PetscErrorCode PCSetFromOptions_Composite(PC pc)
154: {
155:   PC_Composite     *jac = (PC_Composite*)pc->data;
156:   PetscErrorCode   ierr;
157:   PetscInt         nmax = 8,i;
158:   PC_CompositeLink next;
159:   char             *pcs[8];
160:   PetscTruth       flg;

163:   PetscOptionsHead("Composite preconditioner options");
164:     PetscOptionsEnum("-pc_composite_type","Type of composition","PCCompositeSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&jac->type,&flg);
165:     if (flg) {
166:       PCCompositeSetType(pc,jac->type);
167:     }
168:     PetscOptionsName("-pc_composite_true","Use true matrix for inner solves","PCCompositeSetUseTrue",&flg);
169:     if (flg) {
170:       PCCompositeSetUseTrue(pc);
171:     }
172:     PetscOptionsStringArray("-pc_composite_pcs","List of composite solvers","PCCompositeAddPC",pcs,&nmax,&flg);
173:     if (flg) {
174:       for (i=0; i<nmax; i++) {
175:         PCCompositeAddPC(pc,pcs[i]);
176:       }
177:     }
178:   PetscOptionsTail();

180:   next = jac->head;
181:   while (next) {
182:     PCSetFromOptions(next->pc);
183:     next = next->next;
184:   }
185:   return(0);
186: }

190: static PetscErrorCode PCView_Composite(PC pc,PetscViewer viewer)
191: {
192:   PC_Composite     *jac = (PC_Composite*)pc->data;
193:   PetscErrorCode   ierr;
194:   PC_CompositeLink next = jac->head;
195:   PetscTruth       iascii;

198:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
199:   if (iascii) {
200:     PetscViewerASCIIPrintf(viewer,"Composite PC type - %s\n",PCCompositeTypes[jac->type]);
201:     PetscViewerASCIIPrintf(viewer,"PCs on composite preconditioner follow\n");
202:     PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
203:   } else {
204:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCComposite",((PetscObject)viewer)->type_name);
205:   }
206:   if (iascii) {
207:     PetscViewerASCIIPushTab(viewer);
208:   }
209:   while (next) {
210:     PCView(next->pc,viewer);
211:     next = next->next;
212:   }
213:   if (iascii) {
214:     PetscViewerASCIIPopTab(viewer);
215:     PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
216:   }
217:   return(0);
218: }

220: /* ------------------------------------------------------------------------------*/

225: PetscErrorCode  PCCompositeSpecialSetAlpha_Composite(PC pc,PetscScalar alpha)
226: {
227:   PC_Composite *jac = (PC_Composite*)pc->data;
229:   jac->alpha = alpha;
230:   return(0);
231: }

237: PetscErrorCode  PCCompositeSetType_Composite(PC pc,PCCompositeType type)
238: {
240:   if (type == PC_COMPOSITE_ADDITIVE) {
241:     pc->ops->apply = PCApply_Composite_Additive;
242:   } else if (type ==  PC_COMPOSITE_MULTIPLICATIVE || type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
243:     pc->ops->apply = PCApply_Composite_Multiplicative;
244:   } else if (type ==  PC_COMPOSITE_SPECIAL) {
245:     pc->ops->apply = PCApply_Composite_Special;
246:   } else {
247:     SETERRQ(PETSC_ERR_ARG_WRONG,"Unkown composite preconditioner type");
248:   }
249:   return(0);
250: }

256: PetscErrorCode  PCCompositeAddPC_Composite(PC pc,PCType type)
257: {
258:   PC_Composite     *jac;
259:   PC_CompositeLink next,ilink;
260:   PetscErrorCode   ierr;
261:   PetscInt         cnt = 0;
262:   const char       *prefix;
263:   char             newprefix[8];

266:   PetscNew(struct _PC_CompositeLink,&ilink);
267:   ilink->next = 0;
268:   PCCreate(pc->comm,&ilink->pc);

270:   jac  = (PC_Composite*)pc->data;
271:   next = jac->head;
272:   if (!next) {
273:     jac->head       = ilink;
274:     ilink->previous = PETSC_NULL;
275:   } else {
276:     cnt++;
277:     while (next->next) {
278:       next = next->next;
279:       cnt++;
280:     }
281:     next->next      = ilink;
282:     ilink->previous = next;
283:   }
284:   PCGetOptionsPrefix(pc,&prefix);
285:   PCSetOptionsPrefix(ilink->pc,prefix);
286:   sprintf(newprefix,"sub_%d_",(int)cnt);
287:   PCAppendOptionsPrefix(ilink->pc,newprefix);
288:   /* type is set after prefix, because some methods may modify prefix, e.g. pcksp */
289:   PCSetType(ilink->pc,type);
290:   return(0);
291: }

297: PetscErrorCode  PCCompositeGetPC_Composite(PC pc,PetscInt n,PC *subpc)
298: {
299:   PC_Composite     *jac;
300:   PC_CompositeLink next;
301:   PetscInt         i;

304:   jac  = (PC_Composite*)pc->data;
305:   next = jac->head;
306:   for (i=0; i<n; i++) {
307:     if (!next->next) {
308:       SETERRQ(PETSC_ERR_ARG_INCOMP,"Not enough PCs in composite preconditioner");
309:     }
310:     next = next->next;
311:   }
312:   *subpc = next->pc;
313:   return(0);
314: }

320: PetscErrorCode  PCCompositeSetUseTrue_Composite(PC pc)
321: {
322:   PC_Composite   *jac;

325:   jac                  = (PC_Composite*)pc->data;
326:   jac->use_true_matrix = PETSC_TRUE;
327:   return(0);
328: }

331: /* -------------------------------------------------------------------------------- */
334: /*@
335:    PCCompositeSetType - Sets the type of composite preconditioner.
336:    
337:    Collective on PC

339:    Input Parameter:
340: +  pc - the preconditioner context
341: -  type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL

343:    Options Database Key:
344: .  -pc_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type

346:    Level: Developer

348: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
349: @*/
350: PetscErrorCode  PCCompositeSetType(PC pc,PCCompositeType type)
351: {
352:   PetscErrorCode ierr,(*f)(PC,PCCompositeType);

356:   PetscObjectQueryFunction((PetscObject)pc,"PCCompositeSetType_C",(void (**)(void))&f);
357:   if (f) {
358:     (*f)(pc,type);
359:   }
360:   return(0);
361: }

365: /*@
366:    PCCompositeSpecialSetAlpha - Sets alpha for the special composite preconditioner
367:      for alphaI + R + S
368:    
369:    Collective on PC

371:    Input Parameter:
372: +  pc - the preconditioner context
373: -  alpha - scale on identity

375:    Level: Developer

377: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
378: @*/
379: PetscErrorCode  PCCompositeSpecialSetAlpha(PC pc,PetscScalar alpha)
380: {
381:   PetscErrorCode ierr,(*f)(PC,PetscScalar);

385:   PetscObjectQueryFunction((PetscObject)pc,"PCCompositeSpecialSetAlpha_C",(void (**)(void))&f);
386:   if (f) {
387:     (*f)(pc,alpha);
388:   }
389:   return(0);
390: }

394: /*@C
395:    PCCompositeAddPC - Adds another PC to the composite PC.
396:    
397:    Collective on PC

399:    Input Parameters:
400: +  pc - the preconditioner context
401: -  type - the type of the new preconditioner

403:    Level: Developer

405: .keywords: PC, composite preconditioner, add
406: @*/
407: PetscErrorCode  PCCompositeAddPC(PC pc,PCType type)
408: {
409:   PetscErrorCode ierr,(*f)(PC,PCType);

413:   PetscObjectQueryFunction((PetscObject)pc,"PCCompositeAddPC_C",(void (**)(void))&f);
414:   if (f) {
415:     (*f)(pc,type);
416:   }
417:   return(0);
418: }

422: /*@
423:    PCCompositeGetPC - Gets one of the PC objects in the composite PC.
424:    
425:    Not Collective

427:    Input Parameter:
428: +  pc - the preconditioner context
429: -  n - the number of the pc requested

431:    Output Parameters:
432: .  subpc - the PC requested

434:    Level: Developer

436: .keywords: PC, get, composite preconditioner, sub preconditioner

438: .seealso: PCCompositeAddPC()
439: @*/
440: PetscErrorCode  PCCompositeGetPC(PC pc,PetscInt n,PC *subpc)
441: {
442:   PetscErrorCode ierr,(*f)(PC,PetscInt,PC *);

447:   PetscObjectQueryFunction((PetscObject)pc,"PCCompositeGetPC_C",(void (**)(void))&f);
448:   if (f) {
449:     (*f)(pc,n,subpc);
450:   } else {
451:     SETERRQ(PETSC_ERR_ARG_WRONG,"Cannot get pc, not composite type");
452:   }
453:   return(0);
454: }

458: /*@
459:    PCCompositeSetUseTrue - Sets a flag to indicate that the true matrix (rather than
460:                       the matrix used to define the preconditioner) is used to compute
461:                       the residual when the multiplicative scheme is used.

463:    Collective on PC

465:    Input Parameters:
466: .  pc - the preconditioner context

468:    Options Database Key:
469: .  -pc_composite_true - Activates PCCompositeSetUseTrue()

471:    Note:
472:    For the common case in which the preconditioning and linear 
473:    system matrices are identical, this routine is unnecessary.

475:    Level: Developer

477: .keywords: PC, composite preconditioner, set, true, flag

479: .seealso: PCSetOperators(), PCBJacobiSetUseTrueLocal(), PCKSPSetUseTrue()
480: @*/
481: PetscErrorCode  PCCompositeSetUseTrue(PC pc)
482: {
483:   PetscErrorCode ierr,(*f)(PC);

487:   PetscObjectQueryFunction((PetscObject)pc,"PCCompositeSetUseTrue_C",(void (**)(void))&f);
488:   if (f) {
489:     (*f)(pc);
490:   }
491:   return(0);
492: }

494: /* -------------------------------------------------------------------------------------------*/

496: /*MC
497:      PCCOMPOSITE - Build a preconditioner by composing together several preconditioners 

499:    Options Database Keys:
500: +  -pc_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type
501: .  -pc_composite_true - Activates PCCompositeSetUseTrue()
502: -  -pc_composite_pcs - <pc0,pc1,...> list of PCs to compose

504:    Level: intermediate

506:    Concepts: composing solvers

508:    Notes: To use a Krylov method inside the composite preconditioner, set the PCType of one or more
509:           inner PCs to be PCKSP. 
510:           Using a Krylov method inside another Krylov method can be dangerous (you get divergence or
511:           the incorrect answer) unless you use KSPFGMRES as the outter Krylov method


514: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
515:            PCSHELL, PCKSP, PCCompositeSetType(), PCCompositeSpecialSetAlpha(), PCCompositeAddPC(),
516:            PCCompositeGetPC(), PCCompositeSetUseTrue()

518: M*/

523: PetscErrorCode  PCCreate_Composite(PC pc)
524: {
526:   PC_Composite   *jac;

529:   PetscNew(PC_Composite,&jac);
530:   PetscLogObjectMemory(pc,sizeof(PC_Composite));
531:   pc->ops->apply              = PCApply_Composite_Additive;
532:   pc->ops->setup              = PCSetUp_Composite;
533:   pc->ops->destroy            = PCDestroy_Composite;
534:   pc->ops->setfromoptions     = PCSetFromOptions_Composite;
535:   pc->ops->view               = PCView_Composite;
536:   pc->ops->applyrichardson    = 0;

538:   pc->data               = (void*)jac;
539:   jac->type              = PC_COMPOSITE_ADDITIVE;
540:   jac->work1             = 0;
541:   jac->work2             = 0;
542:   jac->head              = 0;

544:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCompositeSetType_C","PCCompositeSetType_Composite",
545:                     PCCompositeSetType_Composite);
546:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCompositeAddPC_C","PCCompositeAddPC_Composite",
547:                     PCCompositeAddPC_Composite);
548:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCompositeGetPC_C","PCCompositeGetPC_Composite",
549:                     PCCompositeGetPC_Composite);
550:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCompositeSetUseTrue_C","PCCompositeSetUseTrue_Composite",
551:                     PCCompositeSetUseTrue_Composite);
552:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCompositeSpecialSetAlpha_C","PCCompositeSpecialSetAlpha_Composite",
553:                     PCCompositeSpecialSetAlpha_Composite);

555:   return(0);
556: }