Actual source code: composite.c

petsc-master 2014-12-19
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;

 33:   if (!next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC() or -pc_composite_pcs");
 34:   if (next->next && !jac->work2) { /* allocate second work vector */
 35:     VecDuplicate(jac->work1,&jac->work2);
 36:   }
 37:   if (pc->useAmat) mat = pc->mat;
 38:   PCApply(next->pc,x,y);                      /* y <- B x */
 39:   while (next->next) {
 40:     next = next->next;
 41:     MatMult(mat,y,jac->work1);                /* work1 <- A y */
 42:     VecWAXPY(jac->work2,-1.0,jac->work1,x);   /* work2 <- x - work1 */
 43:     VecSet(jac->work1,0.0);  /* zero since some PC's may not set all entries in the result */
 44:     PCApply(next->pc,jac->work2,jac->work1);  /* work1 <- C work2 */
 45:     VecAXPY(y,1.0,jac->work1);                /* y <- y + work1 = B x + C (x - A B x) = (B + C (1 - A B)) x */
 46:   }
 47:   if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
 48:     while (next->previous) {
 49:       next = next->previous;
 50:       MatMult(mat,y,jac->work1);
 51:       VecWAXPY(jac->work2,-1.0,jac->work1,x);
 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);
 54:       VecAXPY(y,1.0,jac->work1);
 55:     }
 56:   }
 57:   return(0);
 58: }

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

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

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

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

119:   PCApply(next->pc,x,jac->work1);
120:   PCApply(next->next->pc,jac->work1,y);
121:   return(0);
122: }

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

133:   if (!next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC() or -pc_composite_pcs");
134:   PCApply(next->pc,x,y);
135:   while (next->next) {
136:     next = next->next;
137:     VecSet(jac->work1,0.0);  /* zero since some PC's may not set all entries in the result */
138:     PCApply(next->pc,x,jac->work1);
139:     VecAXPY(y,1.0,jac->work1);
140:   }
141:   return(0);
142: }

146: static PetscErrorCode PCApplyTranspose_Composite_Additive(PC pc,Vec x,Vec y)
147: {
148:   PetscErrorCode   ierr;
149:   PC_Composite     *jac = (PC_Composite*)pc->data;
150:   PC_CompositeLink next = jac->head;

153:   if (!next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC() or -pc_composite_pcs");
154:   PCApplyTranspose(next->pc,x,y);
155:   while (next->next) {
156:     next = next->next;
157:     VecSet(jac->work1,0.0);  /* zero since some PC's may not set all entries in the result */
158:     PCApplyTranspose(next->pc,x,jac->work1);
159:     VecAXPY(y,1.0,jac->work1);
160:   }
161:   return(0);
162: }

166: static PetscErrorCode PCSetUp_Composite(PC pc)
167: {
168:   PetscErrorCode   ierr;
169:   PC_Composite     *jac = (PC_Composite*)pc->data;
170:   PC_CompositeLink next = jac->head;

173:   if (!jac->work1) {
174:     MatCreateVecs(pc->pmat,&jac->work1,0);
175:   }
176:   while (next) {
177:     PCSetOperators(next->pc,pc->mat,pc->pmat);
178:     next = next->next;
179:   }
180:   return(0);
181: }

185: static PetscErrorCode PCReset_Composite(PC pc)
186: {
187:   PC_Composite     *jac = (PC_Composite*)pc->data;
188:   PetscErrorCode   ierr;
189:   PC_CompositeLink next = jac->head;

192:   while (next) {
193:     PCReset(next->pc);
194:     next = next->next;
195:   }
196:   VecDestroy(&jac->work1);
197:   VecDestroy(&jac->work2);
198:   return(0);
199: }

203: static PetscErrorCode PCDestroy_Composite(PC pc)
204: {
205:   PC_Composite     *jac = (PC_Composite*)pc->data;
206:   PetscErrorCode   ierr;
207:   PC_CompositeLink next = jac->head,next_tmp;

210:   PCReset_Composite(pc);
211:   while (next) {
212:     PCDestroy(&next->pc);
213:     next_tmp = next;
214:     next     = next->next;
215:     PetscFree(next_tmp);
216:   }
217:   PetscFree(pc->data);
218:   return(0);
219: }

223: static PetscErrorCode PCSetFromOptions_Composite(PC pc)
224: {
225:   PC_Composite     *jac = (PC_Composite*)pc->data;
226:   PetscErrorCode   ierr;
227:   PetscInt         nmax = 8,i;
228:   PC_CompositeLink next;
229:   char             *pcs[8];
230:   PetscBool        flg;

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

247:   next = jac->head;
248:   while (next) {
249:     PCSetFromOptions(next->pc);
250:     next = next->next;
251:   }
252:   return(0);
253: }

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

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

285: /* ------------------------------------------------------------------------------*/

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: }

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: }

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

326:   *type = jac->type;
327:   return(0);
328: }

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

342:   PetscNewLog(pc,&ilink);
343:   ilink->next = 0;
344:   PCCreate(PetscObjectComm((PetscObject)pc),&ilink->pc);
345:   PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->pc);

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

372: static PetscErrorCode  PCCompositeGetPC_Composite(PC pc,PetscInt n,PC *subpc)
373: {
374:   PC_Composite     *jac;
375:   PC_CompositeLink next;
376:   PetscInt         i;

379:   jac  = (PC_Composite*)pc->data;
380:   next = jac->head;
381:   for (i=0; i<n; i++) {
382:     if (!next->next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Not enough PCs in composite preconditioner");
383:     next = next->next;
384:   }
385:   *subpc = next->pc;
386:   return(0);
387: }

389: /* -------------------------------------------------------------------------------- */
392: /*@
393:    PCCompositeSetType - Sets the type of composite preconditioner.

395:    Logically Collective on PC

397:    Input Parameters:
398: +  pc - the preconditioner context
399: -  type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL

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

404:    Level: Developer

406: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
407: @*/
408: PetscErrorCode  PCCompositeSetType(PC pc,PCCompositeType type)
409: {

415:   PetscTryMethod(pc,"PCCompositeSetType_C",(PC,PCCompositeType),(pc,type));
416:   return(0);
417: }

421: /*@
422:    PCCompositeGetType - Gets the type of composite preconditioner.

424:    Logically Collective on PC

426:    Input Parameter:
427: .  pc - the preconditioner context

429:    Output Parameter:
430: .  type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL

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

435:    Level: Developer

437: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
438: @*/
439: PetscErrorCode  PCCompositeGetType(PC pc,PCCompositeType *type)
440: {

445:   PetscUseMethod(pc,"PCCompositeGetType_C",(PC,PCCompositeType*),(pc,type));
446:   return(0);
447: }

451: /*@
452:    PCCompositeSpecialSetAlpha - Sets alpha for the special composite preconditioner
453:      for alphaI + R + S

455:    Logically Collective on PC

457:    Input Parameter:
458: +  pc - the preconditioner context
459: -  alpha - scale on identity

461:    Level: Developer

463: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
464: @*/
465: PetscErrorCode  PCCompositeSpecialSetAlpha(PC pc,PetscScalar alpha)
466: {

472:   PetscTryMethod(pc,"PCCompositeSpecialSetAlpha_C",(PC,PetscScalar),(pc,alpha));
473:   return(0);
474: }

478: /*@C
479:    PCCompositeAddPC - Adds another PC to the composite PC.

481:    Collective on PC

483:    Input Parameters:
484: +  pc - the preconditioner context
485: -  type - the type of the new preconditioner

487:    Level: Developer

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

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

503: /*@
504:    PCCompositeGetPC - Gets one of the PC objects in the composite PC.

506:    Not Collective

508:    Input Parameter:
509: +  pc - the preconditioner context
510: -  n - the number of the pc requested

512:    Output Parameters:
513: .  subpc - the PC requested

515:    Level: Developer

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

519: .seealso: PCCompositeAddPC()
520: @*/
521: PetscErrorCode  PCCompositeGetPC(PC pc,PetscInt n,PC *subpc)
522: {

528:   PetscUseMethod(pc,"PCCompositeGetPC_C",(PC,PetscInt,PC*),(pc,n,subpc));
529:   return(0);
530: }

532: /* -------------------------------------------------------------------------------------------*/

534: /*MC
535:      PCCOMPOSITE - Build a preconditioner by composing together several preconditioners

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

542:    Level: intermediate

544:    Concepts: composing solvers

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


552: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
553:            PCSHELL, PCKSP, PCCompositeSetType(), PCCompositeSpecialSetAlpha(), PCCompositeAddPC(),
554:            PCCompositeGetPC(), PCSetUseAmat()

556: M*/

560: PETSC_EXTERN PetscErrorCode PCCreate_Composite(PC pc)
561: {
563:   PC_Composite   *jac;

566:   PetscNewLog(pc,&jac);

568:   pc->ops->apply           = PCApply_Composite_Additive;
569:   pc->ops->applytranspose  = PCApplyTranspose_Composite_Additive;
570:   pc->ops->setup           = PCSetUp_Composite;
571:   pc->ops->reset           = PCReset_Composite;
572:   pc->ops->destroy         = PCDestroy_Composite;
573:   pc->ops->setfromoptions  = PCSetFromOptions_Composite;
574:   pc->ops->view            = PCView_Composite;
575:   pc->ops->applyrichardson = 0;

577:   pc->data   = (void*)jac;
578:   jac->type  = PC_COMPOSITE_ADDITIVE;
579:   jac->work1 = 0;
580:   jac->work2 = 0;
581:   jac->head  = 0;

583:   PetscObjectComposeFunction((PetscObject)pc,"PCCompositeSetType_C",PCCompositeSetType_Composite);
584:   PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetType_C",PCCompositeGetType_Composite);
585:   PetscObjectComposeFunction((PetscObject)pc,"PCCompositeAddPC_C",PCCompositeAddPC_Composite);
586:   PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetPC_C",PCCompositeGetPC_Composite);
587:   PetscObjectComposeFunction((PetscObject)pc,"PCCompositeSpecialSetAlpha_C",PCCompositeSpecialSetAlpha_Composite);
588:   return(0);
589: }