Actual source code: composite.c


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

  8: typedef struct _PC_CompositeLink *PC_CompositeLink;
  9: struct _PC_CompositeLink {
 10:   PC               pc;
 11:   PC_CompositeLink next;
 12:   PC_CompositeLink previous;
 13: };

 15: typedef struct {
 16:   PC_CompositeLink head;
 17:   PCCompositeType  type;
 18:   Vec              work1;
 19:   Vec              work2;
 20:   PetscScalar      alpha;
 21: } PC_Composite;

 23: static PetscErrorCode PCApply_Composite_Multiplicative(PC pc,Vec x,Vec y)
 24: {
 25:   PetscErrorCode   ierr;
 26:   PC_Composite     *jac = (PC_Composite*)pc->data;
 27:   PC_CompositeLink next = jac->head;
 28:   Mat              mat  = pc->pmat;


 32:   if (!next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPCType() or -pc_composite_pcs");

 34:   /* Set the reuse flag on children PCs */
 35:   while (next) {
 36:     PCSetReusePreconditioner(next->pc,pc->reusepreconditioner);
 37:     next = next->next;
 38:   }
 39:   next = jac->head;

 41:   if (next->next && !jac->work2) { /* allocate second work vector */
 42:     VecDuplicate(jac->work1,&jac->work2);
 43:   }
 44:   if (pc->useAmat) mat = pc->mat;
 45:   PCApply(next->pc,x,y);                      /* y <- B x */
 46:   while (next->next) {
 47:     next = next->next;
 48:     MatMult(mat,y,jac->work1);                /* work1 <- A y */
 49:     VecWAXPY(jac->work2,-1.0,jac->work1,x);   /* work2 <- x - work1 */
 50:     PCApply(next->pc,jac->work2,jac->work1);  /* work1 <- C work2 */
 51:     VecAXPY(y,1.0,jac->work1);                /* y <- y + work1 = B x + C (x - A B x) = (B + C (1 - A B)) x */
 52:   }
 53:   if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
 54:     while (next->previous) {
 55:       next = next->previous;
 56:       MatMult(mat,y,jac->work1);
 57:       VecWAXPY(jac->work2,-1.0,jac->work1,x);
 58:       PCApply(next->pc,jac->work2,jac->work1);
 59:       VecAXPY(y,1.0,jac->work1);
 60:     }
 61:   }
 62:   return(0);
 63: }

 65: static PetscErrorCode PCApplyTranspose_Composite_Multiplicative(PC pc,Vec x,Vec y)
 66: {
 67:   PetscErrorCode   ierr;
 68:   PC_Composite     *jac = (PC_Composite*)pc->data;
 69:   PC_CompositeLink next = jac->head;
 70:   Mat              mat  = pc->pmat;

 73:   if (!next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPCType() or -pc_composite_pcs");
 74:   if (next->next && !jac->work2) { /* allocate second work vector */
 75:     VecDuplicate(jac->work1,&jac->work2);
 76:   }
 77:   if (pc->useAmat) mat = pc->mat;
 78:   /* locate last PC */
 79:   while (next->next) {
 80:     next = next->next;
 81:   }
 82:   PCApplyTranspose(next->pc,x,y);
 83:   while (next->previous) {
 84:     next = next->previous;
 85:     MatMultTranspose(mat,y,jac->work1);
 86:     VecWAXPY(jac->work2,-1.0,jac->work1,x);
 87:     PCApplyTranspose(next->pc,jac->work2,jac->work1);
 88:     VecAXPY(y,1.0,jac->work1);
 89:   }
 90:   if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
 91:     next = jac->head;
 92:     while (next->next) {
 93:       next = next->next;
 94:       MatMultTranspose(mat,y,jac->work1);
 95:       VecWAXPY(jac->work2,-1.0,jac->work1,x);
 96:       PCApplyTranspose(next->pc,jac->work2,jac->work1);
 97:       VecAXPY(y,1.0,jac->work1);
 98:     }
 99:   }
100:   return(0);
101: }

103: /*
104:     This is very special for a matrix of the form alpha I + R + S
105: where first preconditioner is built from alpha I + S and second from
106: alpha I + R
107: */
108: static PetscErrorCode PCApply_Composite_Special(PC pc,Vec x,Vec y)
109: {
110:   PetscErrorCode   ierr;
111:   PC_Composite     *jac = (PC_Composite*)pc->data;
112:   PC_CompositeLink next = jac->head;

115:   if (!next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPCType() or -pc_composite_pcs");
116:   if (!next->next || next->next->next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Special composite preconditioners requires exactly two PCs");

118:   /* Set the reuse flag on children PCs */
119:   PCSetReusePreconditioner(next->pc,pc->reusepreconditioner);
120:   PCSetReusePreconditioner(next->next->pc,pc->reusepreconditioner);

122:   PCApply(next->pc,x,jac->work1);
123:   PCApply(next->next->pc,jac->work1,y);
124:   return(0);
125: }

127: static PetscErrorCode PCApply_Composite_Additive(PC pc,Vec x,Vec y)
128: {
129:   PetscErrorCode   ierr;
130:   PC_Composite     *jac = (PC_Composite*)pc->data;
131:   PC_CompositeLink next = jac->head;

134:   if (!next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPCType() or -pc_composite_pcs");

136:   /* Set the reuse flag on children PCs */
137:   while (next) {
138:     PCSetReusePreconditioner(next->pc,pc->reusepreconditioner);
139:     next = next->next;
140:   }
141:   next = jac->head;

143:   PCApply(next->pc,x,y);
144:   while (next->next) {
145:     next = next->next;
146:     PCApply(next->pc,x,jac->work1);
147:     VecAXPY(y,1.0,jac->work1);
148:   }
149:   return(0);
150: }

152: static PetscErrorCode PCApplyTranspose_Composite_Additive(PC pc,Vec x,Vec y)
153: {
154:   PetscErrorCode   ierr;
155:   PC_Composite     *jac = (PC_Composite*)pc->data;
156:   PC_CompositeLink next = jac->head;

159:   if (!next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPCType() or -pc_composite_pcs");
160:   PCApplyTranspose(next->pc,x,y);
161:   while (next->next) {
162:     next = next->next;
163:     PCApplyTranspose(next->pc,x,jac->work1);
164:     VecAXPY(y,1.0,jac->work1);
165:   }
166:   return(0);
167: }

169: static PetscErrorCode PCSetUp_Composite(PC pc)
170: {
171:   PetscErrorCode   ierr;
172:   PC_Composite     *jac = (PC_Composite*)pc->data;
173:   PC_CompositeLink next = jac->head;
174:   DM               dm;

177:   if (!jac->work1) {
178:     MatCreateVecs(pc->pmat,&jac->work1,NULL);
179:   }
180:   PCGetDM(pc,&dm);
181:   while (next) {
182:     if (!next->pc->dm) {
183:       PCSetDM(next->pc,dm);
184:     }
185:     if (!next->pc->mat) {
186:       PCSetOperators(next->pc,pc->mat,pc->pmat);
187:     }
188:     next = next->next;
189:   }
190:   return(0);
191: }

193: static PetscErrorCode PCReset_Composite(PC pc)
194: {
195:   PC_Composite     *jac = (PC_Composite*)pc->data;
196:   PetscErrorCode   ierr;
197:   PC_CompositeLink next = jac->head;

200:   while (next) {
201:     PCReset(next->pc);
202:     next = next->next;
203:   }
204:   VecDestroy(&jac->work1);
205:   VecDestroy(&jac->work2);
206:   return(0);
207: }

209: static PetscErrorCode PCDestroy_Composite(PC pc)
210: {
211:   PC_Composite     *jac = (PC_Composite*)pc->data;
212:   PetscErrorCode   ierr;
213:   PC_CompositeLink next = jac->head,next_tmp;

216:   PCReset_Composite(pc);
217:   while (next) {
218:     PCDestroy(&next->pc);
219:     next_tmp = next;
220:     next     = next->next;
221:     PetscFree(next_tmp);
222:   }
223:   PetscFree(pc->data);
224:   return(0);
225: }

227: static PetscErrorCode PCSetFromOptions_Composite(PetscOptionItems *PetscOptionsObject,PC pc)
228: {
229:   PC_Composite     *jac = (PC_Composite*)pc->data;
230:   PetscErrorCode   ierr;
231:   PetscInt         nmax = 8,i;
232:   PC_CompositeLink next;
233:   char             *pcs[8];
234:   PetscBool        flg;

237:   PetscOptionsHead(PetscOptionsObject,"Composite preconditioner options");
238:   PetscOptionsEnum("-pc_composite_type","Type of composition","PCCompositeSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&jac->type,&flg);
239:   if (flg) {
240:     PCCompositeSetType(pc,jac->type);
241:   }
242:   PetscOptionsStringArray("-pc_composite_pcs","List of composite solvers","PCCompositeAddPCType",pcs,&nmax,&flg);
243:   if (flg) {
244:     for (i=0; i<nmax; i++) {
245:       PCCompositeAddPCType(pc,pcs[i]);
246:       PetscFree(pcs[i]);   /* deallocate string pcs[i], which is allocated in PetscOptionsStringArray() */
247:     }
248:   }
249:   PetscOptionsTail();

251:   next = jac->head;
252:   while (next) {
253:     PCSetFromOptions(next->pc);
254:     next = next->next;
255:   }
256:   return(0);
257: }

259: static PetscErrorCode PCView_Composite(PC pc,PetscViewer viewer)
260: {
261:   PC_Composite     *jac = (PC_Composite*)pc->data;
262:   PetscErrorCode   ierr;
263:   PC_CompositeLink next = jac->head;
264:   PetscBool        iascii;

267:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
268:   if (iascii) {
269:     PetscViewerASCIIPrintf(viewer,"Composite PC type - %s\n",PCCompositeTypes[jac->type]);
270:     PetscViewerASCIIPrintf(viewer,"PCs on composite preconditioner follow\n");
271:     PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
272:   }
273:   if (iascii) {
274:     PetscViewerASCIIPushTab(viewer);
275:   }
276:   while (next) {
277:     PCView(next->pc,viewer);
278:     next = next->next;
279:   }
280:   if (iascii) {
281:     PetscViewerASCIIPopTab(viewer);
282:     PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
283:   }
284:   return(0);
285: }

287: /* ------------------------------------------------------------------------------*/

289: static PetscErrorCode  PCCompositeSpecialSetAlpha_Composite(PC pc,PetscScalar alpha)
290: {
291:   PC_Composite *jac = (PC_Composite*)pc->data;

294:   jac->alpha = alpha;
295:   return(0);
296: }

298: static PetscErrorCode  PCCompositeSetType_Composite(PC pc,PCCompositeType type)
299: {
300:   PC_Composite *jac = (PC_Composite*)pc->data;

303:   if (type == PC_COMPOSITE_ADDITIVE) {
304:     pc->ops->apply          = PCApply_Composite_Additive;
305:     pc->ops->applytranspose = PCApplyTranspose_Composite_Additive;
306:   } else if (type ==  PC_COMPOSITE_MULTIPLICATIVE || type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
307:     pc->ops->apply          = PCApply_Composite_Multiplicative;
308:     pc->ops->applytranspose = PCApplyTranspose_Composite_Multiplicative;
309:   } else if (type ==  PC_COMPOSITE_SPECIAL) {
310:     pc->ops->apply          = PCApply_Composite_Special;
311:     pc->ops->applytranspose = NULL;
312:   } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Unknown composite preconditioner type");
313:   jac->type = type;
314:   return(0);
315: }

317: static PetscErrorCode  PCCompositeGetType_Composite(PC pc,PCCompositeType *type)
318: {
319:   PC_Composite *jac = (PC_Composite*)pc->data;

322:   *type = jac->type;
323:   return(0);
324: }

326: static PetscErrorCode PCCompositeAddPC_Composite(PC pc, PC subpc)
327: {
328:   PC_Composite    *jac;
329:   PC_CompositeLink next, ilink;
330:   PetscInt         cnt = 0;
331:   const char      *prefix;
332:   char             newprefix[20];
333:   PetscErrorCode   ierr;

336:   PetscNewLog(pc, &ilink);
337:   ilink->next = NULL;
338:   ilink->pc   = subpc;

340:   jac  = (PC_Composite *) pc->data;
341:   next = jac->head;
342:   if (!next) {
343:     jac->head       = ilink;
344:     ilink->previous = NULL;
345:   } else {
346:     cnt++;
347:     while (next->next) {
348:       next = next->next;
349:       cnt++;
350:     }
351:     next->next      = ilink;
352:     ilink->previous = next;
353:   }
354:   PCGetOptionsPrefix(pc, &prefix);
355:   PCSetOptionsPrefix(subpc, prefix);
356:   PetscSNPrintf(newprefix, 20, "sub_%d_", (int) cnt);
357:   PCAppendOptionsPrefix(subpc, newprefix);
358:   PetscObjectReference((PetscObject) subpc);
359:   return(0);
360: }

362: static PetscErrorCode PCCompositeAddPCType_Composite(PC pc, PCType type)
363: {
364:   PC             subpc;

368:   PCCreate(PetscObjectComm((PetscObject)pc), &subpc);
369:   PetscObjectIncrementTabLevel((PetscObject) subpc, (PetscObject) pc, 1);
370:   PetscLogObjectParent((PetscObject) pc, (PetscObject) subpc);
371:   PCCompositeAddPC_Composite(pc, subpc);
372:   /* type is set after prefix, because some methods may modify prefix, e.g. pcksp */
373:   PCSetType(subpc, type);
374:   PCDestroy(&subpc);
375:   return(0);
376: }

378: static PetscErrorCode  PCCompositeGetNumberPC_Composite(PC pc,PetscInt *n)
379: {
380:   PC_Composite     *jac;
381:   PC_CompositeLink next;

384:   jac  = (PC_Composite*)pc->data;
385:   next = jac->head;
386:   *n = 0;
387:   while (next) {
388:     next = next->next;
389:     (*n) ++;
390:   }
391:   return(0);
392: }

394: static PetscErrorCode  PCCompositeGetPC_Composite(PC pc,PetscInt n,PC *subpc)
395: {
396:   PC_Composite     *jac;
397:   PC_CompositeLink next;
398:   PetscInt         i;

401:   jac  = (PC_Composite*)pc->data;
402:   next = jac->head;
403:   for (i=0; i<n; i++) {
404:     if (!next->next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Not enough PCs in composite preconditioner");
405:     next = next->next;
406:   }
407:   *subpc = next->pc;
408:   return(0);
409: }

411: /* -------------------------------------------------------------------------------- */
412: /*@
413:    PCCompositeSetType - Sets the type of composite preconditioner.

415:    Logically Collective on PC

417:    Input Parameters:
418: +  pc - the preconditioner context
419: -  type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL

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

424:    Level: Developer

426: @*/
427: PetscErrorCode  PCCompositeSetType(PC pc,PCCompositeType type)
428: {

434:   PetscTryMethod(pc,"PCCompositeSetType_C",(PC,PCCompositeType),(pc,type));
435:   return(0);
436: }

438: /*@
439:    PCCompositeGetType - Gets the type of composite preconditioner.

441:    Logically Collective on PC

443:    Input Parameter:
444: .  pc - the preconditioner context

446:    Output Parameter:
447: .  type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL

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

452:    Level: Developer

454: @*/
455: PetscErrorCode  PCCompositeGetType(PC pc,PCCompositeType *type)
456: {

461:   PetscUseMethod(pc,"PCCompositeGetType_C",(PC,PCCompositeType*),(pc,type));
462:   return(0);
463: }

465: /*@
466:    PCCompositeSpecialSetAlpha - Sets alpha for the special composite preconditioner
467:      for alphaI + R + S

469:    Logically Collective on PC

471:    Input Parameter:
472: +  pc - the preconditioner context
473: -  alpha - scale on identity

475:    Level: Developer

477: @*/
478: PetscErrorCode  PCCompositeSpecialSetAlpha(PC pc,PetscScalar alpha)
479: {

485:   PetscTryMethod(pc,"PCCompositeSpecialSetAlpha_C",(PC,PetscScalar),(pc,alpha));
486:   return(0);
487: }

489: /*@C
490:   PCCompositeAddPCType - Adds another PC of the given type to the composite PC.

492:   Collective on PC

494:   Input Parameters:
495: + pc - the preconditioner context
496: - type - the type of the new preconditioner

498:   Level: Developer

500: .seealso: PCCompositeAddPC()
501: @*/
502: PetscErrorCode  PCCompositeAddPCType(PC pc,PCType type)
503: {

508:   PetscTryMethod(pc,"PCCompositeAddPCType_C",(PC,PCType),(pc,type));
509:   return(0);
510: }

512: /*@
513:   PCCompositeAddPC - Adds another PC to the composite PC.

515:   Collective on PC

517:   Input Parameters:
518: + pc    - the preconditioner context
519: - subpc - the new preconditioner

521:    Level: Developer

523: .seealso: PCCompositeAddPCType()
524: @*/
525: PetscErrorCode PCCompositeAddPC(PC pc, PC subpc)
526: {

532:   PetscTryMethod(pc,"PCCompositeAddPC_C",(PC,PC),(pc,subpc));
533:   return(0);
534: }

536: /*@
537:    PCCompositeGetNumberPC - Gets the number of PC objects in the composite PC.

539:    Not Collective

541:    Input Parameter:
542: .  pc - the preconditioner context

544:    Output Parameter:
545: .  num - the number of sub pcs

547:    Level: Developer

549: .seealso: PCCompositeGetPC()
550: @*/
551: PetscErrorCode  PCCompositeGetNumberPC(PC pc,PetscInt *num)
552: {

558:   PetscUseMethod(pc,"PCCompositeGetNumberPC_C",(PC,PetscInt*),(pc,num));
559:   return(0);
560: }

562: /*@
563:    PCCompositeGetPC - Gets one of the PC objects in the composite PC.

565:    Not Collective

567:    Input Parameter:
568: +  pc - the preconditioner context
569: -  n - the number of the pc requested

571:    Output Parameters:
572: .  subpc - the PC requested

574:    Level: Developer

576:     Notes:
577:     To use a different operator to construct one of the inner preconditioners first call PCCompositeGetPC(), then
578:             call PCSetOperators() on that PC.

580: .seealso: PCCompositeAddPCType(), PCCompositeGetNumberPC(), PCSetOperators()
581: @*/
582: PetscErrorCode PCCompositeGetPC(PC pc,PetscInt n,PC *subpc)
583: {

589:   PetscUseMethod(pc,"PCCompositeGetPC_C",(PC,PetscInt,PC*),(pc,n,subpc));
590:   return(0);
591: }

593: /* -------------------------------------------------------------------------------------------*/

595: /*MC
596:      PCCOMPOSITE - Build a preconditioner by composing together several preconditioners

598:    Options Database Keys:
599: +  -pc_composite_type <type: one of multiplicative, additive, symmetric_multiplicative, special> - Sets composite preconditioner type
600: .  -pc_use_amat - activates PCSetUseAmat()
601: -  -pc_composite_pcs - <pc0,pc1,...> list of PCs to compose

603:    Level: intermediate

605:    Notes:
606:     To use a Krylov method inside the composite preconditioner, set the PCType of one or more
607:           inner PCs to be PCKSP.
608:           Using a Krylov method inside another Krylov method can be dangerous (you get divergence or
609:           the incorrect answer) unless you use KSPFGMRES as the outer Krylov method
610:           To use a different operator to construct one of the inner preconditioners first call PCCompositeGetPC(), then
611:           call PCSetOperators() on that PC.


614: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
615:            PCSHELL, PCKSP, PCCompositeSetType(), PCCompositeSpecialSetAlpha(), PCCompositeAddPCType(),
616:            PCCompositeGetPC(), PCSetUseAmat()

618: M*/

620: PETSC_EXTERN PetscErrorCode PCCreate_Composite(PC pc)
621: {
623:   PC_Composite   *jac;

626:   PetscNewLog(pc,&jac);

628:   pc->ops->apply           = PCApply_Composite_Additive;
629:   pc->ops->applytranspose  = PCApplyTranspose_Composite_Additive;
630:   pc->ops->setup           = PCSetUp_Composite;
631:   pc->ops->reset           = PCReset_Composite;
632:   pc->ops->destroy         = PCDestroy_Composite;
633:   pc->ops->setfromoptions  = PCSetFromOptions_Composite;
634:   pc->ops->view            = PCView_Composite;
635:   pc->ops->applyrichardson = NULL;

637:   pc->data   = (void*)jac;
638:   jac->type  = PC_COMPOSITE_ADDITIVE;
639:   jac->work1 = NULL;
640:   jac->work2 = NULL;
641:   jac->head  = NULL;

643:   PetscObjectComposeFunction((PetscObject)pc,"PCCompositeSetType_C",PCCompositeSetType_Composite);
644:   PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetType_C",PCCompositeGetType_Composite);
645:   PetscObjectComposeFunction((PetscObject)pc,"PCCompositeAddPCType_C",PCCompositeAddPCType_Composite);
646:   PetscObjectComposeFunction((PetscObject)pc,"PCCompositeAddPC_C",PCCompositeAddPC_Composite);
647:   PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetNumberPC_C",PCCompositeGetNumberPC_Composite);
648:   PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetPC_C",PCCompositeGetPC_Composite);
649:   PetscObjectComposeFunction((PetscObject)pc,"PCCompositeSpecialSetAlpha_C",PCCompositeSpecialSetAlpha_Composite);
650:   return(0);
651: }