Actual source code: composite.c

  1: /*$Id: composite.c,v 1.45 2001/08/07 03:03:39 balay Exp $*/
  2: /*
  3:       Defines a preconditioner that can consist of a collection of PCs
  4: */
 5:  #include src/sles/pc/pcimpl.h
 6:  #include petscsles.h

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

 23: #undef __FUNCT__  
 25: static int PCApply_Composite_Multiplicative(PC pc,Vec x,Vec y)
 26: {
 27:   int              ierr;
 28:   PC_Composite     *jac = (PC_Composite*)pc->data;
 29:   PC_CompositeLink next = jac->head;
 30:   PetscScalar      one = 1.0,mone = -1.0;
 31:   Mat              mat = pc->pmat;

 34:   if (!next) {
 35:     SETERRQ(1,"No composite preconditioners supplied via PCCompositeAddPC()");
 36:   }
 37:   if (next->next && !jac->work2) { /* allocate second work vector */
 38:     VecDuplicate(jac->work1,&jac->work2);
 39:   }
 40:   PCApply(next->pc,x,y);
 41:   if (jac->use_true_matrix) mat = pc->mat;
 42:   while (next->next) {
 43:     next = next->next;
 44:     MatMult(mat,y,jac->work1);
 45:     VecWAXPY(&mone,jac->work1,x,jac->work2);
 46:     PCApply(next->pc,jac->work2,jac->work1);
 47:     VecAXPY(&one,jac->work1,y);
 48:   }

 50:   return(0);
 51: }

 53: /*
 54:     This is very special for a matrix of the form alpha I + R + S
 55: where first preconditioner is built from alpha I + S and second from
 56: alpha I + R
 57: */
 58: #undef __FUNCT__  
 60: static int PCApply_Composite_Special(PC pc,Vec x,Vec y)
 61: {
 62:   int              ierr;
 63:   PC_Composite     *jac = (PC_Composite*)pc->data;
 64:   PC_CompositeLink next = jac->head;

 67:   if (!next) {
 68:     SETERRQ(1,"No composite preconditioners supplied via PCCompositeAddPC()");
 69:   }
 70:   if (!next->next || next->next->next) {
 71:     SETERRQ(1,"Special composite preconditioners requires exactly two PCs");
 72:   }

 74:   PCApply(next->pc,x,jac->work1);
 75:   PCApply(next->next->pc,jac->work1,y);
 76:   return(0);
 77: }

 79: #undef __FUNCT__  
 81: static int PCApply_Composite_Additive(PC pc,Vec x,Vec y)
 82: {
 83:   int              ierr;
 84:   PC_Composite     *jac = (PC_Composite*)pc->data;
 85:   PC_CompositeLink next = jac->head;
 86:   PetscScalar      one = 1.0;

 89:   if (!next) {
 90:     SETERRQ(1,"No composite preconditioners supplied via PCCompositeAddPC()");
 91:   }
 92:   PCApply(next->pc,x,y);
 93:   while (next->next) {
 94:     next = next->next;
 95:     PCApply(next->pc,x,jac->work1);
 96:     VecAXPY(&one,jac->work1,y);
 97:   }
 98:   return(0);
 99: }

101: #undef __FUNCT__  
103: static int PCSetUp_Composite(PC pc)
104: {
105:   int              ierr;
106:   PC_Composite     *jac = (PC_Composite*)pc->data;
107:   PC_CompositeLink next = jac->head;

110:   if (!jac->work1) {
111:      VecDuplicate(pc->vec,&jac->work1);
112:   }
113:   while (next) {
114:     PCSetOperators(next->pc,pc->mat,pc->pmat,pc->flag);
115:     PCSetVector(next->pc,jac->work1);
116:     next = next->next;
117:   }

119:   return(0);
120: }

122: #undef __FUNCT__  
124: static int PCDestroy_Composite(PC pc)
125: {
126:   PC_Composite     *jac = (PC_Composite*)pc->data;
127:   int              ierr;
128:   PC_CompositeLink next = jac->head;

131:   while (next) {
132:     PCDestroy(next->pc);
133:     next = next->next;
134:   }

136:   if (jac->work1) {VecDestroy(jac->work1);}
137:   if (jac->work2) {VecDestroy(jac->work2);}
138:   PetscFree(jac);
139:   return(0);
140: }

142: #undef __FUNCT__  
144: static int PCSetFromOptions_Composite(PC pc)
145: {
146:   PC_Composite     *jac = (PC_Composite*)pc->data;
147:   int              ierr,nmax = 8,i;
148:   PCCompositeType  type=PC_COMPOSITE_ADDITIVE;
149:   PC_CompositeLink next;
150:   char             *pcs[8],stype[16],*types[] = {"multiplicative","additive","special"};
151:   PetscTruth       flg;

154:   PetscOptionsHead("Composite preconditioner options");
155:     PetscOptionsEList("-pc_composite_type","Type of composition","PCCompositeSetType",types,3,"multiplicative",stype,16,&flg);
156:     if (flg) {
157:       PetscTruth ismult,isadd,isspecial;

159:       PetscStrcmp(stype,types[0],&ismult);
160:       PetscStrcmp(stype,types[1],&isadd);
161:       PetscStrcmp(stype,types[2],&isspecial);

163:       if (ismult)          type = PC_COMPOSITE_MULTIPLICATIVE;
164:       else if (isadd)      type = PC_COMPOSITE_ADDITIVE;
165:       else if (isspecial)  type = PC_COMPOSITE_SPECIAL;
166:       else SETERRQ(1,"Unknown composite type given");

168:       PCCompositeSetType(pc,type);
169:     }
170:     PetscOptionsName("-pc_composite_true","Use true matrix for inner solves","PCCompositeSetUseTrue",&flg);
171:     if (flg) {
172:       PCCompositeSetUseTrue(pc);
173:     }
174:     PetscOptionsStringArray("-pc_composite_pcs","List of composite solvers","PCCompositeAddPC",pcs,&nmax,&flg);
175:     if (flg) {
176:       for (i=0; i<nmax; i++) {
177:         PCCompositeAddPC(pc,pcs[i]);
178:       }
179:     }
180:   PetscOptionsTail();

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

190: #undef __FUNCT__  
192: static int PCView_Composite(PC pc,PetscViewer viewer)
193: {
194:   PC_Composite     *jac = (PC_Composite*)pc->data;
195:   int              ierr;
196:   PC_CompositeLink next = jac->head;
197:   PetscTruth       isascii;

200:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
201:   if (isascii) {
202:     PetscViewerASCIIPrintf(viewer,"PCs on composite preconditioner follown");
203:     PetscViewerASCIIPrintf(viewer,"---------------------------------n");
204:   } else {
205:     SETERRQ1(1,"Viewer type %s not supported for PCComposite",((PetscObject)viewer)->type_name);
206:   }
207:   if (isascii) {
208:     PetscViewerASCIIPushTab(viewer);
209:   }
210:   while (next) {
211:     PCView(next->pc,viewer);
212:     next = next->next;
213:   }
214:   if (isascii) {
215:     PetscViewerASCIIPopTab(viewer);
216:     PetscViewerASCIIPrintf(viewer,"---------------------------------n");
217:   }
218:   return(0);
219: }

221: /* ------------------------------------------------------------------------------*/

223: EXTERN_C_BEGIN
224: #undef __FUNCT__  
226: int PCCompositeSpecialSetAlpha_Composite(PC pc,PetscScalar alpha)
227: {
228:   PC_Composite *jac = (PC_Composite*)pc->data;
230:   jac->alpha = alpha;
231:   return(0);
232: }
233: EXTERN_C_END

235: EXTERN_C_BEGIN
236: #undef __FUNCT__  
238: int PCCompositeSetType_Composite(PC pc,PCCompositeType type)
239: {
241:   if (type == PC_COMPOSITE_ADDITIVE) {
242:     pc->ops->apply = PCApply_Composite_Additive;
243:   } else if (type ==  PC_COMPOSITE_MULTIPLICATIVE) {
244:     pc->ops->apply = PCApply_Composite_Multiplicative;
245:   } else if (type ==  PC_COMPOSITE_SPECIAL) {
246:     pc->ops->apply = PCApply_Composite_Special;
247:   } else {
248:     SETERRQ(1,"Unkown composite preconditioner type");
249:   }
250:   return(0);
251: }
252: EXTERN_C_END

254: EXTERN_C_BEGIN
255: #undef __FUNCT__  
257: int PCCompositeAddPC_Composite(PC pc,PCType type)
258: {
259:   PC_Composite     *jac;
260:   PC_CompositeLink next,link;
261:   int              ierr,cnt = 0;
262:   char             *prefix,newprefix[8];

265:   ierr       = PetscNew(struct _PC_CompositeLink,&link);
266:   link->next = 0;
267:   PCCreate(pc->comm,&link->pc);

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

288:   return(0);
289: }
290: EXTERN_C_END

292: EXTERN_C_BEGIN
293: #undef __FUNCT__  
295: int PCCompositeGetPC_Composite(PC pc,int n,PC *subpc)
296: {
297:   PC_Composite     *jac;
298:   PC_CompositeLink next;
299:   int              i;

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

315: EXTERN_C_BEGIN
316: #undef __FUNCT__  
318: int PCCompositeSetUseTrue_Composite(PC pc)
319: {
320:   PC_Composite   *jac;

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

329: /* -------------------------------------------------------------------------------- */
330: #undef __FUNCT__  
332: /*@C
333:    PCCompositeSetType - Sets the type of composite preconditioner.
334:    
335:    Collective on PC

337:    Input Parameter:
338: .  pc - the preconditioner context
339: .  type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL

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

344:    Level: Developer

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

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

361: #undef __FUNCT__  
363: /*@C
364:    PCCompositeSpecialSetAlpha - Sets alpha for the special composite preconditioner
365:      for alphaI + R + S
366:    
367:    Collective on PC

369:    Input Parameter:
370: +  pc - the preconditioner context
371: -  alpha - scale on identity

373:    Level: Developer

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

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

390: #undef __FUNCT__  
392: /*@C
393:    PCCompositeAddPC - Adds another PC to the composite PC.
394:    
395:    Collective on PC

397:    Input Parameters:
398: .  pc - the preconditioner context
399: .  type - the type of the new preconditioner

401:    Level: Developer

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

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

418: #undef __FUNCT__  
420: /*@C
421:    PCCompositeGetPC - Gets one of the PC objects in the composite PC.
422:    
423:    Not Collective

425:    Input Parameter:
426: .  pc - the preconditioner context
427: .  n - the number of the pc requested

429:    Output Parameters:
430: .  subpc - the PC requested

432:    Level: Developer

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

436: .seealso: PCCompositeAddPC()
437: @*/
438: int PCCompositeGetPC(PC pc,int n,PC *subpc)
439: {
440:   int ierr,(*f)(PC,int,PC *);

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

453: #undef __FUNCT__  
455: /*@
456:    PCCompositeSetUseTrue - Sets a flag to indicate that the true matrix (rather than
457:                       the matrix used to define the preconditioner) is used to compute
458:                       the residual when the multiplicative scheme is used.

460:    Collective on PC

462:    Input Parameters:
463: .  pc - the preconditioner context

465:    Options Database Key:
466: .  -pc_composite_true - Activates PCCompositeSetUseTrue()

468:    Note:
469:    For the common case in which the preconditioning and linear 
470:    system matrices are identical, this routine is unnecessary.

472:    Level: Developer

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

476: .seealso: PCSetOperators(), PCBJacobiSetUseTrueLocal(), PCSLESSetUseTrue()
477: @*/
478: int PCCompositeSetUseTrue(PC pc)
479: {
480:   int ierr,(*f)(PC);

484:   PetscObjectQueryFunction((PetscObject)pc,"PCCompositeSetUseTrue_C",(void (**)(void))&f);
485:   if (f) {
486:     (*f)(pc);
487:   }
488:   return(0);
489: }

491: /* -------------------------------------------------------------------------------------------*/

493: /*MC
494:      PCCOMPOSITE - Build a preconditioner by composing together several preconditioners 

496:    Options Database Keys:
497: .  -pc_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type
498: .  -pc_composite_true - Activates PCCompositeSetUseTrue()

500:    Level: intermediate

502:    Concepts: composing solvers

504:    Notes: To use a Krylov method inside the composite preconditioner, set the PCType of one or more
505:           inner PCs to be PCSLES. 
506:           Using a Krylov method inside another Krylov method can be dangerous (you get divergence or
507:           the incorrect answer) unless you use KSPFGMRES as the other Krylov method


510: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
511:            PCSHELL, PCSLES, PCCompositeSetType(), PCCompositeSpecialSetAlpha(), PCCompositeAddPC(),
512:            PCCompositeGetPC(), PCCompositeSetUseTrue()

514: M*/

516: EXTERN_C_BEGIN
517: #undef __FUNCT__  
519: int PCCreate_Composite(PC pc)
520: {
521:   int            ierr;
522:   PC_Composite   *jac;

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

534:   pc->data               = (void*)jac;
535:   jac->type              = PC_COMPOSITE_ADDITIVE;
536:   jac->work1             = 0;
537:   jac->work2             = 0;
538:   jac->head              = 0;

540:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCompositeSetType_C","PCCompositeSetType_Composite",
541:                     PCCompositeSetType_Composite);
542:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCompositeAddPC_C","PCCompositeAddPC_Composite",
543:                     PCCompositeAddPC_Composite);
544:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCompositeGetPC_C","PCCompositeGetPC_Composite",
545:                     PCCompositeGetPC_Composite);
546:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCompositeSetUseTrue_C","PCCompositeSetUseTrue_Composite",
547:                     PCCompositeSetUseTrue_Composite);
548:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCompositeSpecialSetAlpha_C","PCCompositeSpecialSetAlpha_Composite",
549:                     PCCompositeSpecialSetAlpha_Composite);

551:   return(0);
552: }
553: EXTERN_C_END