Actual source code: composite.c
petsc-3.5.4 2015-05-23
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: MatGetVecs(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 PCCompositeAddPC_Composite(PC pc,PCType type)
322: {
323: PC_Composite *jac;
324: PC_CompositeLink next,ilink;
325: PetscErrorCode ierr;
326: PetscInt cnt = 0;
327: const char *prefix;
328: char newprefix[8];
331: PetscNewLog(pc,&ilink);
332: ilink->next = 0;
333: PCCreate(PetscObjectComm((PetscObject)pc),&ilink->pc);
334: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->pc);
336: jac = (PC_Composite*)pc->data;
337: next = jac->head;
338: if (!next) {
339: jac->head = ilink;
340: ilink->previous = NULL;
341: } else {
342: cnt++;
343: while (next->next) {
344: next = next->next;
345: cnt++;
346: }
347: next->next = ilink;
348: ilink->previous = next;
349: }
350: PCGetOptionsPrefix(pc,&prefix);
351: PCSetOptionsPrefix(ilink->pc,prefix);
352: sprintf(newprefix,"sub_%d_",(int)cnt);
353: PCAppendOptionsPrefix(ilink->pc,newprefix);
354: /* type is set after prefix, because some methods may modify prefix, e.g. pcksp */
355: PCSetType(ilink->pc,type);
356: return(0);
357: }
361: static PetscErrorCode PCCompositeGetPC_Composite(PC pc,PetscInt n,PC *subpc)
362: {
363: PC_Composite *jac;
364: PC_CompositeLink next;
365: PetscInt i;
368: jac = (PC_Composite*)pc->data;
369: next = jac->head;
370: for (i=0; i<n; i++) {
371: if (!next->next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Not enough PCs in composite preconditioner");
372: next = next->next;
373: }
374: *subpc = next->pc;
375: return(0);
376: }
378: /* -------------------------------------------------------------------------------- */
381: /*@
382: PCCompositeSetType - Sets the type of composite preconditioner.
384: Logically Collective on PC
386: Input Parameter:
387: + pc - the preconditioner context
388: - type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL
390: Options Database Key:
391: . -pc_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type
393: Level: Developer
395: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
396: @*/
397: PetscErrorCode PCCompositeSetType(PC pc,PCCompositeType type)
398: {
404: PetscTryMethod(pc,"PCCompositeSetType_C",(PC,PCCompositeType),(pc,type));
405: return(0);
406: }
410: /*@
411: PCCompositeSpecialSetAlpha - Sets alpha for the special composite preconditioner
412: for alphaI + R + S
414: Logically Collective on PC
416: Input Parameter:
417: + pc - the preconditioner context
418: - alpha - scale on identity
420: Level: Developer
422: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
423: @*/
424: PetscErrorCode PCCompositeSpecialSetAlpha(PC pc,PetscScalar alpha)
425: {
431: PetscTryMethod(pc,"PCCompositeSpecialSetAlpha_C",(PC,PetscScalar),(pc,alpha));
432: return(0);
433: }
437: /*@C
438: PCCompositeAddPC - Adds another PC to the composite PC.
440: Collective on PC
442: Input Parameters:
443: + pc - the preconditioner context
444: - type - the type of the new preconditioner
446: Level: Developer
448: .keywords: PC, composite preconditioner, add
449: @*/
450: PetscErrorCode PCCompositeAddPC(PC pc,PCType type)
451: {
456: PetscTryMethod(pc,"PCCompositeAddPC_C",(PC,PCType),(pc,type));
457: return(0);
458: }
462: /*@
463: PCCompositeGetPC - Gets one of the PC objects in the composite PC.
465: Not Collective
467: Input Parameter:
468: + pc - the preconditioner context
469: - n - the number of the pc requested
471: Output Parameters:
472: . subpc - the PC requested
474: Level: Developer
476: .keywords: PC, get, composite preconditioner, sub preconditioner
478: .seealso: PCCompositeAddPC()
479: @*/
480: PetscErrorCode PCCompositeGetPC(PC pc,PetscInt n,PC *subpc)
481: {
487: PetscUseMethod(pc,"PCCompositeGetPC_C",(PC,PetscInt,PC*),(pc,n,subpc));
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, symmetric_multiplicative, special> - Sets composite preconditioner type
498: . -pc_use_amat - Activates PCSetUseAmat()
499: - -pc_composite_pcs - <pc0,pc1,...> list of PCs to compose
501: Level: intermediate
503: Concepts: composing solvers
505: Notes: To use a Krylov method inside the composite preconditioner, set the PCType of one or more
506: inner PCs to be PCKSP.
507: Using a Krylov method inside another Krylov method can be dangerous (you get divergence or
508: the incorrect answer) unless you use KSPFGMRES as the outer Krylov method
511: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
512: PCSHELL, PCKSP, PCCompositeSetType(), PCCompositeSpecialSetAlpha(), PCCompositeAddPC(),
513: PCCompositeGetPC(), PCSetUseAmat()
515: M*/
519: PETSC_EXTERN PetscErrorCode PCCreate_Composite(PC pc)
520: {
522: PC_Composite *jac;
525: PetscNewLog(pc,&jac);
527: pc->ops->apply = PCApply_Composite_Additive;
528: pc->ops->applytranspose = PCApplyTranspose_Composite_Additive;
529: pc->ops->setup = PCSetUp_Composite;
530: pc->ops->reset = PCReset_Composite;
531: pc->ops->destroy = PCDestroy_Composite;
532: pc->ops->setfromoptions = PCSetFromOptions_Composite;
533: pc->ops->view = PCView_Composite;
534: pc->ops->applyrichardson = 0;
536: pc->data = (void*)jac;
537: jac->type = PC_COMPOSITE_ADDITIVE;
538: jac->work1 = 0;
539: jac->work2 = 0;
540: jac->head = 0;
542: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeSetType_C",PCCompositeSetType_Composite);
543: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeAddPC_C",PCCompositeAddPC_Composite);
544: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetPC_C",PCCompositeGetPC_Composite);
545: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeSpecialSetAlpha_C",PCCompositeSpecialSetAlpha_Composite);
546: return(0);
547: }