Actual source code: composite.c

petsc-master 2017-06-21
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>

  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 PCCompositeAddPC() 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:     VecSet(jac->work1,0.0);  /* zero since some PC's may not set all entries in the result */
 51:     PCApply(next->pc,jac->work2,jac->work1);  /* work1 <- C work2 */
 52:     VecAXPY(y,1.0,jac->work1);                /* y <- y + work1 = B x + C (x - A B x) = (B + C (1 - A B)) x */
 53:   }
 54:   if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
 55:     while (next->previous) {
 56:       next = next->previous;
 57:       MatMult(mat,y,jac->work1);
 58:       VecWAXPY(jac->work2,-1.0,jac->work1,x);
 59:       VecSet(jac->work1,0.0);  /* zero since some PC's may not set all entries in the result */
 60:       PCApply(next->pc,jac->work2,jac->work1);
 61:       VecAXPY(y,1.0,jac->work1);
 62:     }
 63:   }
 64:   return(0);
 65: }

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

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

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

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

122:   /* Set the reuse flag on children PCs */
123:   PCSetReusePreconditioner(next->pc,pc->reusepreconditioner);
124:   PCSetReusePreconditioner(next->next->pc,pc->reusepreconditioner);

126:   PCApply(next->pc,x,jac->work1);
127:   PCApply(next->next->pc,jac->work1,y);
128:   return(0);
129: }

131: static PetscErrorCode PCApply_Composite_Additive(PC pc,Vec x,Vec y)
132: {
133:   PetscErrorCode   ierr;
134:   PC_Composite     *jac = (PC_Composite*)pc->data;
135:   PC_CompositeLink next = jac->head;

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

140:   /* Set the reuse flag on children PCs */
141:   while (next) {
142:     PCSetReusePreconditioner(next->pc,pc->reusepreconditioner);
143:     next = next->next;
144:   }
145:   next = jac->head;

147:   PCApply(next->pc,x,y);
148:   while (next->next) {
149:     next = next->next;
150:     VecSet(jac->work1,0.0);  /* zero since some PC's may not set all entries in the result */
151:     PCApply(next->pc,x,jac->work1);
152:     VecAXPY(y,1.0,jac->work1);
153:   }
154:   return(0);
155: }

157: static PetscErrorCode PCApplyTranspose_Composite_Additive(PC pc,Vec x,Vec y)
158: {
159:   PetscErrorCode   ierr;
160:   PC_Composite     *jac = (PC_Composite*)pc->data;
161:   PC_CompositeLink next = jac->head;

164:   if (!next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC() or -pc_composite_pcs");
165:   PCApplyTranspose(next->pc,x,y);
166:   while (next->next) {
167:     next = next->next;
168:     VecSet(jac->work1,0.0);  /* zero since some PC's may not set all entries in the result */
169:     PCApplyTranspose(next->pc,x,jac->work1);
170:     VecAXPY(y,1.0,jac->work1);
171:   }
172:   return(0);
173: }

175: static PetscErrorCode PCSetUp_Composite(PC pc)
176: {
177:   PetscErrorCode   ierr;
178:   PC_Composite     *jac = (PC_Composite*)pc->data;
179:   PC_CompositeLink next = jac->head;
180:   DM               dm;

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

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

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

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

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

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

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

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

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

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

289: /* ------------------------------------------------------------------------------*/

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

296:   jac->alpha = alpha;
297:   return(0);
298: }

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

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

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

324:   *type = jac->type;
325:   return(0);
326: }

328: static PetscErrorCode  PCCompositeAddPC_Composite(PC pc,PCType type)
329: {
330:   PC_Composite     *jac;
331:   PC_CompositeLink next,ilink;
332:   PetscErrorCode   ierr;
333:   PetscInt         cnt = 0;
334:   const char       *prefix;
335:   char             newprefix[8];

338:   PetscNewLog(pc,&ilink);
339:   ilink->next = 0;
340:   PCCreate(PetscObjectComm((PetscObject)pc),&ilink->pc);
341:   PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->pc);

343:   jac  = (PC_Composite*)pc->data;
344:   next = jac->head;
345:   if (!next) {
346:     jac->head       = ilink;
347:     ilink->previous = NULL;
348:   } else {
349:     cnt++;
350:     while (next->next) {
351:       next = next->next;
352:       cnt++;
353:     }
354:     next->next      = ilink;
355:     ilink->previous = next;
356:   }
357:   PCGetOptionsPrefix(pc,&prefix);
358:   PCSetOptionsPrefix(ilink->pc,prefix);
359:   sprintf(newprefix,"sub_%d_",(int)cnt);
360:   PCAppendOptionsPrefix(ilink->pc,newprefix);
361:   /* type is set after prefix, because some methods may modify prefix, e.g. pcksp */
362:   PCSetType(ilink->pc,type);
363:   return(0);
364: }

366: static PetscErrorCode  PCCompositeGetNumberPC_Composite(PC pc,PetscInt *n)
367: {
368:   PC_Composite     *jac;
369:   PC_CompositeLink next;

372:   jac  = (PC_Composite*)pc->data;
373:   next = jac->head;
374:   *n = 0;
375:   while (next) {
376:     next = next->next;
377:     (*n) ++;
378:   }
379:   return(0);
380: }

382: static PetscErrorCode  PCCompositeGetPC_Composite(PC pc,PetscInt n,PC *subpc)
383: {
384:   PC_Composite     *jac;
385:   PC_CompositeLink next;
386:   PetscInt         i;

389:   jac  = (PC_Composite*)pc->data;
390:   next = jac->head;
391:   for (i=0; i<n; i++) {
392:     if (!next->next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Not enough PCs in composite preconditioner");
393:     next = next->next;
394:   }
395:   *subpc = next->pc;
396:   return(0);
397: }

399: /* -------------------------------------------------------------------------------- */
400: /*@
401:    PCCompositeSetType - Sets the type of composite preconditioner.

403:    Logically Collective on PC

405:    Input Parameters:
406: +  pc - the preconditioner context
407: -  type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL

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

412:    Level: Developer

414: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
415: @*/
416: PetscErrorCode  PCCompositeSetType(PC pc,PCCompositeType type)
417: {

423:   PetscTryMethod(pc,"PCCompositeSetType_C",(PC,PCCompositeType),(pc,type));
424:   return(0);
425: }

427: /*@
428:    PCCompositeGetType - Gets the type of composite preconditioner.

430:    Logically Collective on PC

432:    Input Parameter:
433: .  pc - the preconditioner context

435:    Output Parameter:
436: .  type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL

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

441:    Level: Developer

443: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
444: @*/
445: PetscErrorCode  PCCompositeGetType(PC pc,PCCompositeType *type)
446: {

451:   PetscUseMethod(pc,"PCCompositeGetType_C",(PC,PCCompositeType*),(pc,type));
452:   return(0);
453: }

455: /*@
456:    PCCompositeSpecialSetAlpha - Sets alpha for the special composite preconditioner
457:      for alphaI + R + S

459:    Logically Collective on PC

461:    Input Parameter:
462: +  pc - the preconditioner context
463: -  alpha - scale on identity

465:    Level: Developer

467: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
468: @*/
469: PetscErrorCode  PCCompositeSpecialSetAlpha(PC pc,PetscScalar alpha)
470: {

476:   PetscTryMethod(pc,"PCCompositeSpecialSetAlpha_C",(PC,PetscScalar),(pc,alpha));
477:   return(0);
478: }

480: /*@C
481:    PCCompositeAddPC - Adds another PC to the composite PC.

483:    Collective on PC

485:    Input Parameters:
486: +  pc - the preconditioner context
487: -  type - the type of the new preconditioner

489:    Level: Developer

491: .keywords: PC, composite preconditioner, add
492: @*/
493: PetscErrorCode  PCCompositeAddPC(PC pc,PCType type)
494: {

499:   PetscTryMethod(pc,"PCCompositeAddPC_C",(PC,PCType),(pc,type));
500:   return(0);
501: }

503: /*@
504:    PCCompositeGetNumberPC - Gets the number of PC objects in the composite PC.

506:    Not Collective

508:    Input Parameter:
509: .  pc - the preconditioner context

511:    Output Parameter:
512: .  num - the number of sub pcs

514:    Level: Developer

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

518: .seealso: PCCompositeGetPC()
519: @*/
520: PetscErrorCode  PCCompositeGetNumberPC(PC pc,PetscInt *num)
521: {

527:   PetscUseMethod(pc,"PCCompositeGetNumberPC_C",(PC,PetscInt*),(pc,num));
528:   return(0);
529: }

531: /*@
532:    PCCompositeGetPC - Gets one of the PC objects in the composite PC.

534:    Not Collective

536:    Input Parameter:
537: +  pc - the preconditioner context
538: -  n - the number of the pc requested

540:    Output Parameters:
541: .  subpc - the PC requested

543:    Level: Developer

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

547: .seealso: PCCompositeAddPC(), PCCompositeGetNumberPC()
548: @*/
549: PetscErrorCode  PCCompositeGetPC(PC pc,PetscInt n,PC *subpc)
550: {

556:   PetscUseMethod(pc,"PCCompositeGetPC_C",(PC,PetscInt,PC*),(pc,n,subpc));
557:   return(0);
558: }

560: /* -------------------------------------------------------------------------------------------*/

562: /*MC
563:      PCCOMPOSITE - Build a preconditioner by composing together several preconditioners

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

570:    Level: intermediate

572:    Concepts: composing solvers

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


580: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
581:            PCSHELL, PCKSP, PCCompositeSetType(), PCCompositeSpecialSetAlpha(), PCCompositeAddPC(),
582:            PCCompositeGetPC(), PCSetUseAmat()

584: M*/

586: PETSC_EXTERN PetscErrorCode PCCreate_Composite(PC pc)
587: {
589:   PC_Composite   *jac;

592:   PetscNewLog(pc,&jac);

594:   pc->ops->apply           = PCApply_Composite_Additive;
595:   pc->ops->applytranspose  = PCApplyTranspose_Composite_Additive;
596:   pc->ops->setup           = PCSetUp_Composite;
597:   pc->ops->reset           = PCReset_Composite;
598:   pc->ops->destroy         = PCDestroy_Composite;
599:   pc->ops->setfromoptions  = PCSetFromOptions_Composite;
600:   pc->ops->view            = PCView_Composite;
601:   pc->ops->applyrichardson = 0;

603:   pc->data   = (void*)jac;
604:   jac->type  = PC_COMPOSITE_ADDITIVE;
605:   jac->work1 = 0;
606:   jac->work2 = 0;
607:   jac->head  = 0;

609:   PetscObjectComposeFunction((PetscObject)pc,"PCCompositeSetType_C",PCCompositeSetType_Composite);
610:   PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetType_C",PCCompositeGetType_Composite);
611:   PetscObjectComposeFunction((PetscObject)pc,"PCCompositeAddPC_C",PCCompositeAddPC_Composite);
612:   PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetNumberPC_C",PCCompositeGetNumberPC_Composite);
613:   PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetPC_C",PCCompositeGetPC_Composite);
614:   PetscObjectComposeFunction((PetscObject)pc,"PCCompositeSpecialSetAlpha_C",PCCompositeSpecialSetAlpha_Composite);
615:   return(0);
616: }