Actual source code: composite.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  2: /*
  3:       Defines a preconditioner that can consist of a collection of PCs
  4: */
  5: #include <petsc/private/pcimpl.h>
  6: #include <petscksp.h>            /*I "petscksp.h" I*/

  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;

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


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

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

 43:   if (next->next && !jac->work2) { /* allocate second work vector */
 44:     VecDuplicate(jac->work1,&jac->work2);
 45:   }
 46:   if (pc->useAmat) mat = pc->mat;
 47:   PCApply(next->pc,x,y);                      /* y <- B x */
 48:   while (next->next) {
 49:     next = next->next;
 50:     MatMult(mat,y,jac->work1);                /* work1 <- A y */
 51:     VecWAXPY(jac->work2,-1.0,jac->work1,x);   /* work2 <- x - work1 */
 52:     VecSet(jac->work1,0.0);  /* zero since some PC's may not set all entries in the result */
 53:     PCApply(next->pc,jac->work2,jac->work1);  /* work1 <- C work2 */
 54:     VecAXPY(y,1.0,jac->work1);                /* y <- y + work1 = B x + C (x - A B x) = (B + C (1 - A B)) x */
 55:   }
 56:   if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
 57:     while (next->previous) {
 58:       next = next->previous;
 59:       MatMult(mat,y,jac->work1);
 60:       VecWAXPY(jac->work2,-1.0,jac->work1,x);
 61:       VecSet(jac->work1,0.0);  /* zero since some PC's may not set all entries in the result */
 62:       PCApply(next->pc,jac->work2,jac->work1);
 63:       VecAXPY(y,1.0,jac->work1);
 64:     }
 65:   }
 66:   return(0);
 67: }

 71: static PetscErrorCode PCApplyTranspose_Composite_Multiplicative(PC pc,Vec x,Vec y)
 72: {
 73:   PetscErrorCode   ierr;
 74:   PC_Composite     *jac = (PC_Composite*)pc->data;
 75:   PC_CompositeLink next = jac->head;
 76:   Mat              mat  = pc->pmat;

 79:   if (!next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC() or -pc_composite_pcs");
 80:   if (next->next && !jac->work2) { /* allocate second work vector */
 81:     VecDuplicate(jac->work1,&jac->work2);
 82:   }
 83:   if (pc->useAmat) mat = pc->mat;
 84:   /* locate last PC */
 85:   while (next->next) {
 86:     next = next->next;
 87:   }
 88:   PCApplyTranspose(next->pc,x,y);
 89:   while (next->previous) {
 90:     next = next->previous;
 91:     MatMultTranspose(mat,y,jac->work1);
 92:     VecWAXPY(jac->work2,-1.0,jac->work1,x);
 93:     VecSet(jac->work1,0.0);  /* zero since some PC's may not set all entries in the result */
 94:     PCApplyTranspose(next->pc,jac->work2,jac->work1);
 95:     VecAXPY(y,1.0,jac->work1);
 96:   }
 97:   if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
 98:     next = jac->head;
 99:     while (next->next) {
100:       next = next->next;
101:       MatMultTranspose(mat,y,jac->work1);
102:       VecWAXPY(jac->work2,-1.0,jac->work1,x);
103:       VecSet(jac->work1,0.0);  /* zero since some PC's may not set all entries in the result */
104:       PCApplyTranspose(next->pc,jac->work2,jac->work1);
105:       VecAXPY(y,1.0,jac->work1);
106:     }
107:   }
108:   return(0);
109: }

111: /*
112:     This is very special for a matrix of the form alpha I + R + S
113: where first preconditioner is built from alpha I + S and second from
114: alpha I + R
115: */
118: static PetscErrorCode PCApply_Composite_Special(PC pc,Vec x,Vec y)
119: {
120:   PetscErrorCode   ierr;
121:   PC_Composite     *jac = (PC_Composite*)pc->data;
122:   PC_CompositeLink next = jac->head;

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

128:   /* Set the reuse flag on children PCs */
129:   PCSetReusePreconditioner(next->pc,pc->reusepreconditioner);
130:   PCSetReusePreconditioner(next->next->pc,pc->reusepreconditioner);

132:   PCApply(next->pc,x,jac->work1);
133:   PCApply(next->next->pc,jac->work1,y);
134:   return(0);
135: }

139: static PetscErrorCode PCApply_Composite_Additive(PC pc,Vec x,Vec y)
140: {
141:   PetscErrorCode   ierr;
142:   PC_Composite     *jac = (PC_Composite*)pc->data;
143:   PC_CompositeLink next = jac->head;

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

148:   /* Set the reuse flag on children PCs */
149:   while (next) {
150:     PCSetReusePreconditioner(next->pc,pc->reusepreconditioner);
151:     next = next->next;
152:   }
153:   next = jac->head;

155:   PCApply(next->pc,x,y);
156:   while (next->next) {
157:     next = next->next;
158:     VecSet(jac->work1,0.0);  /* zero since some PC's may not set all entries in the result */
159:     PCApply(next->pc,x,jac->work1);
160:     VecAXPY(y,1.0,jac->work1);
161:   }
162:   return(0);
163: }

167: static PetscErrorCode PCApplyTranspose_Composite_Additive(PC pc,Vec x,Vec y)
168: {
169:   PetscErrorCode   ierr;
170:   PC_Composite     *jac = (PC_Composite*)pc->data;
171:   PC_CompositeLink next = jac->head;

174:   if (!next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC() or -pc_composite_pcs");
175:   PCApplyTranspose(next->pc,x,y);
176:   while (next->next) {
177:     next = next->next;
178:     VecSet(jac->work1,0.0);  /* zero since some PC's may not set all entries in the result */
179:     PCApplyTranspose(next->pc,x,jac->work1);
180:     VecAXPY(y,1.0,jac->work1);
181:   }
182:   return(0);
183: }

187: static PetscErrorCode PCSetUp_Composite(PC pc)
188: {
189:   PetscErrorCode   ierr;
190:   PC_Composite     *jac = (PC_Composite*)pc->data;
191:   PC_CompositeLink next = jac->head;
192:   DM               dm;

195:   if (!jac->work1) {
196:     MatCreateVecs(pc->pmat,&jac->work1,0);
197:   }
198:   PCGetDM(pc,&dm);
199:   while (next) {
200:     PCSetDM(next->pc,dm);
201:     PCSetOperators(next->pc,pc->mat,pc->pmat);
202:     next = next->next;
203:   }
204:   return(0);
205: }

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

216:   while (next) {
217:     PCReset(next->pc);
218:     next = next->next;
219:   }
220:   VecDestroy(&jac->work1);
221:   VecDestroy(&jac->work2);
222:   return(0);
223: }

227: static PetscErrorCode PCDestroy_Composite(PC pc)
228: {
229:   PC_Composite     *jac = (PC_Composite*)pc->data;
230:   PetscErrorCode   ierr;
231:   PC_CompositeLink next = jac->head,next_tmp;

234:   PCReset_Composite(pc);
235:   while (next) {
236:     PCDestroy(&next->pc);
237:     next_tmp = next;
238:     next     = next->next;
239:     PetscFree(next_tmp);
240:   }
241:   PetscFree(pc->data);
242:   return(0);
243: }

247: static PetscErrorCode PCSetFromOptions_Composite(PetscOptionItems *PetscOptionsObject,PC pc)
248: {
249:   PC_Composite     *jac = (PC_Composite*)pc->data;
250:   PetscErrorCode   ierr;
251:   PetscInt         nmax = 8,i;
252:   PC_CompositeLink next;
253:   char             *pcs[8];
254:   PetscBool        flg;

257:   PetscOptionsHead(PetscOptionsObject,"Composite preconditioner options");
258:   PetscOptionsEnum("-pc_composite_type","Type of composition","PCCompositeSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&jac->type,&flg);
259:   if (flg) {
260:     PCCompositeSetType(pc,jac->type);
261:   }
262:   PetscOptionsStringArray("-pc_composite_pcs","List of composite solvers","PCCompositeAddPC",pcs,&nmax,&flg);
263:   if (flg) {
264:     for (i=0; i<nmax; i++) {
265:       PCCompositeAddPC(pc,pcs[i]);
266:       PetscFree(pcs[i]);   /* deallocate string pcs[i], which is allocated in PetscOptionsStringArray() */
267:     }
268:   }
269:   PetscOptionsTail();

271:   next = jac->head;
272:   while (next) {
273:     PCSetFromOptions(next->pc);
274:     next = next->next;
275:   }
276:   return(0);
277: }

281: static PetscErrorCode PCView_Composite(PC pc,PetscViewer viewer)
282: {
283:   PC_Composite     *jac = (PC_Composite*)pc->data;
284:   PetscErrorCode   ierr;
285:   PC_CompositeLink next = jac->head;
286:   PetscBool        iascii;

289:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
290:   if (iascii) {
291:     PetscViewerASCIIPrintf(viewer,"Composite PC type - %s\n",PCCompositeTypes[jac->type]);
292:     PetscViewerASCIIPrintf(viewer,"PCs on composite preconditioner follow\n");
293:     PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
294:   }
295:   if (iascii) {
296:     PetscViewerASCIIPushTab(viewer);
297:   }
298:   while (next) {
299:     PCView(next->pc,viewer);
300:     next = next->next;
301:   }
302:   if (iascii) {
303:     PetscViewerASCIIPopTab(viewer);
304:     PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
305:   }
306:   return(0);
307: }

309: /* ------------------------------------------------------------------------------*/

313: static PetscErrorCode  PCCompositeSpecialSetAlpha_Composite(PC pc,PetscScalar alpha)
314: {
315:   PC_Composite *jac = (PC_Composite*)pc->data;

318:   jac->alpha = alpha;
319:   return(0);
320: }

324: static PetscErrorCode  PCCompositeSetType_Composite(PC pc,PCCompositeType type)
325: {
326:   PC_Composite *jac = (PC_Composite*)pc->data;

329:   if (type == PC_COMPOSITE_ADDITIVE) {
330:     pc->ops->apply          = PCApply_Composite_Additive;
331:     pc->ops->applytranspose = PCApplyTranspose_Composite_Additive;
332:   } else if (type ==  PC_COMPOSITE_MULTIPLICATIVE || type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
333:     pc->ops->apply          = PCApply_Composite_Multiplicative;
334:     pc->ops->applytranspose = PCApplyTranspose_Composite_Multiplicative;
335:   } else if (type ==  PC_COMPOSITE_SPECIAL) {
336:     pc->ops->apply          = PCApply_Composite_Special;
337:     pc->ops->applytranspose = NULL;
338:   } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Unkown composite preconditioner type");
339:   jac->type = type;
340:   return(0);
341: }

345: static PetscErrorCode  PCCompositeGetType_Composite(PC pc,PCCompositeType *type)
346: {
347:   PC_Composite *jac = (PC_Composite*)pc->data;

350:   *type = jac->type;
351:   return(0);
352: }

356: static PetscErrorCode  PCCompositeAddPC_Composite(PC pc,PCType type)
357: {
358:   PC_Composite     *jac;
359:   PC_CompositeLink next,ilink;
360:   PetscErrorCode   ierr;
361:   PetscInt         cnt = 0;
362:   const char       *prefix;
363:   char             newprefix[8];

366:   PetscNewLog(pc,&ilink);
367:   ilink->next = 0;
368:   PCCreate(PetscObjectComm((PetscObject)pc),&ilink->pc);
369:   PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->pc);

371:   jac  = (PC_Composite*)pc->data;
372:   next = jac->head;
373:   if (!next) {
374:     jac->head       = ilink;
375:     ilink->previous = NULL;
376:   } else {
377:     cnt++;
378:     while (next->next) {
379:       next = next->next;
380:       cnt++;
381:     }
382:     next->next      = ilink;
383:     ilink->previous = next;
384:   }
385:   PCGetOptionsPrefix(pc,&prefix);
386:   PCSetOptionsPrefix(ilink->pc,prefix);
387:   sprintf(newprefix,"sub_%d_",(int)cnt);
388:   PCAppendOptionsPrefix(ilink->pc,newprefix);
389:   /* type is set after prefix, because some methods may modify prefix, e.g. pcksp */
390:   PCSetType(ilink->pc,type);
391:   return(0);
392: }

396: static PetscErrorCode  PCCompositeGetNumberPC_Composite(PC pc,PetscInt *n)
397: {
398:   PC_Composite     *jac;
399:   PC_CompositeLink next;

402:   jac  = (PC_Composite*)pc->data;
403:   next = jac->head;
404:   *n = 0;
405:   while (next) {
406:     next = next->next;
407:     (*n) ++;
408:   }
409:   return(0);
410: }

414: static PetscErrorCode  PCCompositeGetPC_Composite(PC pc,PetscInt n,PC *subpc)
415: {
416:   PC_Composite     *jac;
417:   PC_CompositeLink next;
418:   PetscInt         i;

421:   jac  = (PC_Composite*)pc->data;
422:   next = jac->head;
423:   for (i=0; i<n; i++) {
424:     if (!next->next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Not enough PCs in composite preconditioner");
425:     next = next->next;
426:   }
427:   *subpc = next->pc;
428:   return(0);
429: }

431: /* -------------------------------------------------------------------------------- */
434: /*@
435:    PCCompositeSetType - Sets the type of composite preconditioner.

437:    Logically Collective on PC

439:    Input Parameters:
440: +  pc - the preconditioner context
441: -  type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL

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

446:    Level: Developer

448: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
449: @*/
450: PetscErrorCode  PCCompositeSetType(PC pc,PCCompositeType type)
451: {

457:   PetscTryMethod(pc,"PCCompositeSetType_C",(PC,PCCompositeType),(pc,type));
458:   return(0);
459: }

463: /*@
464:    PCCompositeGetType - Gets the type of composite preconditioner.

466:    Logically Collective on PC

468:    Input Parameter:
469: .  pc - the preconditioner context

471:    Output Parameter:
472: .  type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL

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

477:    Level: Developer

479: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
480: @*/
481: PetscErrorCode  PCCompositeGetType(PC pc,PCCompositeType *type)
482: {

487:   PetscUseMethod(pc,"PCCompositeGetType_C",(PC,PCCompositeType*),(pc,type));
488:   return(0);
489: }

493: /*@
494:    PCCompositeSpecialSetAlpha - Sets alpha for the special composite preconditioner
495:      for alphaI + R + S

497:    Logically Collective on PC

499:    Input Parameter:
500: +  pc - the preconditioner context
501: -  alpha - scale on identity

503:    Level: Developer

505: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
506: @*/
507: PetscErrorCode  PCCompositeSpecialSetAlpha(PC pc,PetscScalar alpha)
508: {

514:   PetscTryMethod(pc,"PCCompositeSpecialSetAlpha_C",(PC,PetscScalar),(pc,alpha));
515:   return(0);
516: }

520: /*@C
521:    PCCompositeAddPC - Adds another PC to the composite PC.

523:    Collective on PC

525:    Input Parameters:
526: +  pc - the preconditioner context
527: -  type - the type of the new preconditioner

529:    Level: Developer

531: .keywords: PC, composite preconditioner, add
532: @*/
533: PetscErrorCode  PCCompositeAddPC(PC pc,PCType type)
534: {

539:   PetscTryMethod(pc,"PCCompositeAddPC_C",(PC,PCType),(pc,type));
540:   return(0);
541: }

545: /*@
546:    PCCompositeGetNumberPC - Gets the number of PC objects in the composite PC.

548:    Not Collective

550:    Input Parameter:
551: .  pc - the preconditioner context

553:    Output Parameter:
554: .  num - the number of sub pcs

556:    Level: Developer

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

560: .seealso: PCCompositeGetPC()
561: @*/
562: PetscErrorCode  PCCompositeGetNumberPC(PC pc,PetscInt *num)
563: {

569:   PetscUseMethod(pc,"PCCompositeGetNumberPC_C",(PC,PetscInt*),(pc,num));
570:   return(0);
571: }

575: /*@
576:    PCCompositeGetPC - Gets one of the PC objects in the composite PC.

578:    Not Collective

580:    Input Parameter:
581: +  pc - the preconditioner context
582: -  n - the number of the pc requested

584:    Output Parameters:
585: .  subpc - the PC requested

587:    Level: Developer

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

591: .seealso: PCCompositeAddPC(), PCCompositeGetNumberPC()
592: @*/
593: PetscErrorCode  PCCompositeGetPC(PC pc,PetscInt n,PC *subpc)
594: {

600:   PetscUseMethod(pc,"PCCompositeGetPC_C",(PC,PetscInt,PC*),(pc,n,subpc));
601:   return(0);
602: }

604: /* -------------------------------------------------------------------------------------------*/

606: /*MC
607:      PCCOMPOSITE - Build a preconditioner by composing together several preconditioners

609:    Options Database Keys:
610: +  -pc_composite_type <type: one of multiplicative, additive, symmetric_multiplicative, special> - Sets composite preconditioner type
611: .  -pc_use_amat - Activates PCSetUseAmat()
612: -  -pc_composite_pcs - <pc0,pc1,...> list of PCs to compose

614:    Level: intermediate

616:    Concepts: composing solvers

618:    Notes: To use a Krylov method inside the composite preconditioner, set the PCType of one or more
619:           inner PCs to be PCKSP.
620:           Using a Krylov method inside another Krylov method can be dangerous (you get divergence or
621:           the incorrect answer) unless you use KSPFGMRES as the outer Krylov method


624: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
625:            PCSHELL, PCKSP, PCCompositeSetType(), PCCompositeSpecialSetAlpha(), PCCompositeAddPC(),
626:            PCCompositeGetPC(), PCSetUseAmat()

628: M*/

632: PETSC_EXTERN PetscErrorCode PCCreate_Composite(PC pc)
633: {
635:   PC_Composite   *jac;

638:   PetscNewLog(pc,&jac);

640:   pc->ops->apply           = PCApply_Composite_Additive;
641:   pc->ops->applytranspose  = PCApplyTranspose_Composite_Additive;
642:   pc->ops->setup           = PCSetUp_Composite;
643:   pc->ops->reset           = PCReset_Composite;
644:   pc->ops->destroy         = PCDestroy_Composite;
645:   pc->ops->setfromoptions  = PCSetFromOptions_Composite;
646:   pc->ops->view            = PCView_Composite;
647:   pc->ops->applyrichardson = 0;

649:   pc->data   = (void*)jac;
650:   jac->type  = PC_COMPOSITE_ADDITIVE;
651:   jac->work1 = 0;
652:   jac->work2 = 0;
653:   jac->head  = 0;

655:   PetscObjectComposeFunction((PetscObject)pc,"PCCompositeSetType_C",PCCompositeSetType_Composite);
656:   PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetType_C",PCCompositeGetType_Composite);
657:   PetscObjectComposeFunction((PetscObject)pc,"PCCompositeAddPC_C",PCCompositeAddPC_Composite);
658:   PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetNumberPC_C",PCCompositeGetNumberPC_Composite);
659:   PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetPC_C",PCCompositeGetPC_Composite);
660:   PetscObjectComposeFunction((PetscObject)pc,"PCCompositeSpecialSetAlpha_C",PCCompositeSpecialSetAlpha_Composite);
661:   return(0);
662: }