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