Actual source code: composite.c
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: PC_Composite *jac = (PC_Composite*)pc->data;
26: PC_CompositeLink next = jac->head;
27: Mat mat = pc->pmat;
32: /* Set the reuse flag on children PCs */
33: while (next) {
34: PCSetReusePreconditioner(next->pc,pc->reusepreconditioner);
35: next = next->next;
36: }
37: next = jac->head;
39: if (next->next && !jac->work2) { /* allocate second work vector */
40: VecDuplicate(jac->work1,&jac->work2);
41: }
42: if (pc->useAmat) mat = pc->mat;
43: PCApply(next->pc,x,y); /* y <- B x */
44: while (next->next) {
45: next = next->next;
46: MatMult(mat,y,jac->work1); /* work1 <- A y */
47: VecWAXPY(jac->work2,-1.0,jac->work1,x); /* work2 <- x - work1 */
48: PCApply(next->pc,jac->work2,jac->work1); /* work1 <- C work2 */
49: VecAXPY(y,1.0,jac->work1); /* y <- y + work1 = B x + C (x - A B x) = (B + C (1 - A B)) x */
50: }
51: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
52: while (next->previous) {
53: next = next->previous;
54: MatMult(mat,y,jac->work1);
55: VecWAXPY(jac->work2,-1.0,jac->work1,x);
56: PCApply(next->pc,jac->work2,jac->work1);
57: VecAXPY(y,1.0,jac->work1);
58: }
59: }
60: return 0;
61: }
63: static PetscErrorCode PCApplyTranspose_Composite_Multiplicative(PC pc,Vec x,Vec y)
64: {
65: PC_Composite *jac = (PC_Composite*)pc->data;
66: PC_CompositeLink next = jac->head;
67: Mat mat = pc->pmat;
70: if (next->next && !jac->work2) { /* allocate second work vector */
71: VecDuplicate(jac->work1,&jac->work2);
72: }
73: if (pc->useAmat) mat = pc->mat;
74: /* locate last PC */
75: while (next->next) {
76: next = next->next;
77: }
78: PCApplyTranspose(next->pc,x,y);
79: while (next->previous) {
80: next = next->previous;
81: MatMultTranspose(mat,y,jac->work1);
82: VecWAXPY(jac->work2,-1.0,jac->work1,x);
83: PCApplyTranspose(next->pc,jac->work2,jac->work1);
84: VecAXPY(y,1.0,jac->work1);
85: }
86: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
87: next = jac->head;
88: while (next->next) {
89: next = next->next;
90: MatMultTranspose(mat,y,jac->work1);
91: VecWAXPY(jac->work2,-1.0,jac->work1,x);
92: PCApplyTranspose(next->pc,jac->work2,jac->work1);
93: VecAXPY(y,1.0,jac->work1);
94: }
95: }
96: return 0;
97: }
99: /*
100: This is very special for a matrix of the form alpha I + R + S
101: where first preconditioner is built from alpha I + S and second from
102: alpha I + R
103: */
104: static PetscErrorCode PCApply_Composite_Special(PC pc,Vec x,Vec y)
105: {
106: PC_Composite *jac = (PC_Composite*)pc->data;
107: PC_CompositeLink next = jac->head;
112: /* Set the reuse flag on children PCs */
113: PCSetReusePreconditioner(next->pc,pc->reusepreconditioner);
114: PCSetReusePreconditioner(next->next->pc,pc->reusepreconditioner);
116: PCApply(next->pc,x,jac->work1);
117: PCApply(next->next->pc,jac->work1,y);
118: return 0;
119: }
121: static PetscErrorCode PCApply_Composite_Additive(PC pc,Vec x,Vec y)
122: {
123: PC_Composite *jac = (PC_Composite*)pc->data;
124: PC_CompositeLink next = jac->head;
128: /* Set the reuse flag on children PCs */
129: while (next) {
130: PCSetReusePreconditioner(next->pc,pc->reusepreconditioner);
131: next = next->next;
132: }
133: next = jac->head;
135: PCApply(next->pc,x,y);
136: while (next->next) {
137: next = next->next;
138: PCApply(next->pc,x,jac->work1);
139: VecAXPY(y,1.0,jac->work1);
140: }
141: return 0;
142: }
144: static PetscErrorCode PCApplyTranspose_Composite_Additive(PC pc,Vec x,Vec y)
145: {
146: PC_Composite *jac = (PC_Composite*)pc->data;
147: PC_CompositeLink next = jac->head;
150: PCApplyTranspose(next->pc,x,y);
151: while (next->next) {
152: next = next->next;
153: PCApplyTranspose(next->pc,x,jac->work1);
154: VecAXPY(y,1.0,jac->work1);
155: }
156: return 0;
157: }
159: static PetscErrorCode PCSetUp_Composite(PC pc)
160: {
161: PC_Composite *jac = (PC_Composite*)pc->data;
162: PC_CompositeLink next = jac->head;
163: DM dm;
165: if (!jac->work1) {
166: MatCreateVecs(pc->pmat,&jac->work1,NULL);
167: }
168: PCGetDM(pc,&dm);
169: while (next) {
170: if (!next->pc->dm) {
171: PCSetDM(next->pc,dm);
172: }
173: if (!next->pc->mat) {
174: PCSetOperators(next->pc,pc->mat,pc->pmat);
175: }
176: next = next->next;
177: }
178: return 0;
179: }
181: static PetscErrorCode PCReset_Composite(PC pc)
182: {
183: PC_Composite *jac = (PC_Composite*)pc->data;
184: PC_CompositeLink next = jac->head;
186: while (next) {
187: PCReset(next->pc);
188: next = next->next;
189: }
190: VecDestroy(&jac->work1);
191: VecDestroy(&jac->work2);
192: return 0;
193: }
195: static PetscErrorCode PCDestroy_Composite(PC pc)
196: {
197: PC_Composite *jac = (PC_Composite*)pc->data;
198: PC_CompositeLink next = jac->head,next_tmp;
200: PCReset_Composite(pc);
201: while (next) {
202: PCDestroy(&next->pc);
203: next_tmp = next;
204: next = next->next;
205: PetscFree(next_tmp);
206: }
207: PetscFree(pc->data);
208: return 0;
209: }
211: static PetscErrorCode PCSetFromOptions_Composite(PetscOptionItems *PetscOptionsObject,PC pc)
212: {
213: PC_Composite *jac = (PC_Composite*)pc->data;
214: PetscInt nmax = 8,i;
215: PC_CompositeLink next;
216: char *pcs[8];
217: PetscBool flg;
219: PetscOptionsHead(PetscOptionsObject,"Composite preconditioner options");
220: PetscOptionsEnum("-pc_composite_type","Type of composition","PCCompositeSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&jac->type,&flg);
221: if (flg) {
222: PCCompositeSetType(pc,jac->type);
223: }
224: PetscOptionsStringArray("-pc_composite_pcs","List of composite solvers","PCCompositeAddPCType",pcs,&nmax,&flg);
225: if (flg) {
226: for (i=0; i<nmax; i++) {
227: PCCompositeAddPCType(pc,pcs[i]);
228: PetscFree(pcs[i]); /* deallocate string pcs[i], which is allocated in PetscOptionsStringArray() */
229: }
230: }
231: PetscOptionsTail();
233: next = jac->head;
234: while (next) {
235: PCSetFromOptions(next->pc);
236: next = next->next;
237: }
238: return 0;
239: }
241: static PetscErrorCode PCView_Composite(PC pc,PetscViewer viewer)
242: {
243: PC_Composite *jac = (PC_Composite*)pc->data;
244: PC_CompositeLink next = jac->head;
245: PetscBool iascii;
247: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
248: if (iascii) {
249: PetscViewerASCIIPrintf(viewer,"Composite PC type - %s\n",PCCompositeTypes[jac->type]);
250: PetscViewerASCIIPrintf(viewer,"PCs on composite preconditioner follow\n");
251: PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
252: }
253: if (iascii) {
254: PetscViewerASCIIPushTab(viewer);
255: }
256: while (next) {
257: PCView(next->pc,viewer);
258: next = next->next;
259: }
260: if (iascii) {
261: PetscViewerASCIIPopTab(viewer);
262: PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
263: }
264: return 0;
265: }
267: /* ------------------------------------------------------------------------------*/
269: static PetscErrorCode PCCompositeSpecialSetAlpha_Composite(PC pc,PetscScalar alpha)
270: {
271: PC_Composite *jac = (PC_Composite*)pc->data;
273: jac->alpha = alpha;
274: return 0;
275: }
277: static PetscErrorCode PCCompositeSetType_Composite(PC pc,PCCompositeType type)
278: {
279: PC_Composite *jac = (PC_Composite*)pc->data;
281: if (type == PC_COMPOSITE_ADDITIVE) {
282: pc->ops->apply = PCApply_Composite_Additive;
283: pc->ops->applytranspose = PCApplyTranspose_Composite_Additive;
284: } else if (type == PC_COMPOSITE_MULTIPLICATIVE || type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
285: pc->ops->apply = PCApply_Composite_Multiplicative;
286: pc->ops->applytranspose = PCApplyTranspose_Composite_Multiplicative;
287: } else if (type == PC_COMPOSITE_SPECIAL) {
288: pc->ops->apply = PCApply_Composite_Special;
289: pc->ops->applytranspose = NULL;
290: } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Unknown composite preconditioner type");
291: jac->type = type;
292: return 0;
293: }
295: static PetscErrorCode PCCompositeGetType_Composite(PC pc,PCCompositeType *type)
296: {
297: PC_Composite *jac = (PC_Composite*)pc->data;
299: *type = jac->type;
300: return 0;
301: }
303: static PetscErrorCode PCCompositeAddPC_Composite(PC pc, PC subpc)
304: {
305: PC_Composite *jac;
306: PC_CompositeLink next, ilink;
307: PetscInt cnt = 0;
308: const char *prefix;
309: char newprefix[20];
311: PetscNewLog(pc, &ilink);
312: ilink->next = NULL;
313: ilink->pc = subpc;
315: jac = (PC_Composite *) pc->data;
316: next = jac->head;
317: if (!next) {
318: jac->head = ilink;
319: ilink->previous = NULL;
320: } else {
321: cnt++;
322: while (next->next) {
323: next = next->next;
324: cnt++;
325: }
326: next->next = ilink;
327: ilink->previous = next;
328: }
329: PCGetOptionsPrefix(pc, &prefix);
330: PCSetOptionsPrefix(subpc, prefix);
331: PetscSNPrintf(newprefix, 20, "sub_%d_", (int) cnt);
332: PCAppendOptionsPrefix(subpc, newprefix);
333: PetscObjectReference((PetscObject) subpc);
334: return 0;
335: }
337: static PetscErrorCode PCCompositeAddPCType_Composite(PC pc, PCType type)
338: {
339: PC subpc;
341: PCCreate(PetscObjectComm((PetscObject)pc), &subpc);
342: PetscObjectIncrementTabLevel((PetscObject) subpc, (PetscObject) pc, 1);
343: PetscLogObjectParent((PetscObject) pc, (PetscObject) subpc);
344: PCCompositeAddPC_Composite(pc, subpc);
345: /* type is set after prefix, because some methods may modify prefix, e.g. pcksp */
346: PCSetType(subpc, type);
347: PCDestroy(&subpc);
348: return 0;
349: }
351: static PetscErrorCode PCCompositeGetNumberPC_Composite(PC pc,PetscInt *n)
352: {
353: PC_Composite *jac;
354: PC_CompositeLink next;
356: jac = (PC_Composite*)pc->data;
357: next = jac->head;
358: *n = 0;
359: while (next) {
360: next = next->next;
361: (*n) ++;
362: }
363: return 0;
364: }
366: static PetscErrorCode PCCompositeGetPC_Composite(PC pc,PetscInt n,PC *subpc)
367: {
368: PC_Composite *jac;
369: PC_CompositeLink next;
370: PetscInt i;
372: jac = (PC_Composite*)pc->data;
373: next = jac->head;
374: for (i=0; i<n; i++) {
376: next = next->next;
377: }
378: *subpc = next->pc;
379: return 0;
380: }
382: /* -------------------------------------------------------------------------------- */
383: /*@
384: PCCompositeSetType - Sets the type of composite preconditioner.
386: Logically Collective on PC
388: Input Parameters:
389: + pc - the preconditioner context
390: - type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL
392: Options Database Key:
393: . -pc_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type
395: Level: Developer
397: @*/
398: PetscErrorCode PCCompositeSetType(PC pc,PCCompositeType type)
399: {
402: PetscTryMethod(pc,"PCCompositeSetType_C",(PC,PCCompositeType),(pc,type));
403: return 0;
404: }
406: /*@
407: PCCompositeGetType - Gets the type of composite preconditioner.
409: Logically Collective on PC
411: Input Parameter:
412: . pc - the preconditioner context
414: Output Parameter:
415: . type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL
417: Options Database Key:
418: . -pc_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type
420: Level: Developer
422: @*/
423: PetscErrorCode PCCompositeGetType(PC pc,PCCompositeType *type)
424: {
426: PetscUseMethod(pc,"PCCompositeGetType_C",(PC,PCCompositeType*),(pc,type));
427: return 0;
428: }
430: /*@
431: PCCompositeSpecialSetAlpha - Sets alpha for the special composite preconditioner
432: for alphaI + R + S
434: Logically Collective on PC
436: Input Parameters:
437: + pc - the preconditioner context
438: - alpha - scale on identity
440: Level: Developer
442: @*/
443: PetscErrorCode PCCompositeSpecialSetAlpha(PC pc,PetscScalar alpha)
444: {
447: PetscTryMethod(pc,"PCCompositeSpecialSetAlpha_C",(PC,PetscScalar),(pc,alpha));
448: return 0;
449: }
451: /*@C
452: PCCompositeAddPCType - Adds another PC of the given type to the composite PC.
454: Collective on PC
456: Input Parameters:
457: + pc - the preconditioner context
458: - type - the type of the new preconditioner
460: Level: Developer
462: .seealso: PCCompositeAddPC()
463: @*/
464: PetscErrorCode PCCompositeAddPCType(PC pc,PCType type)
465: {
467: PetscTryMethod(pc,"PCCompositeAddPCType_C",(PC,PCType),(pc,type));
468: return 0;
469: }
471: /*@
472: PCCompositeAddPC - Adds another PC to the composite PC.
474: Collective on PC
476: Input Parameters:
477: + pc - the preconditioner context
478: - subpc - the new preconditioner
480: Level: Developer
482: .seealso: PCCompositeAddPCType()
483: @*/
484: PetscErrorCode PCCompositeAddPC(PC pc, PC subpc)
485: {
488: PetscTryMethod(pc,"PCCompositeAddPC_C",(PC,PC),(pc,subpc));
489: return 0;
490: }
492: /*@
493: PCCompositeGetNumberPC - Gets the number of PC objects in the composite PC.
495: Not Collective
497: Input Parameter:
498: . pc - the preconditioner context
500: Output Parameter:
501: . num - the number of sub pcs
503: Level: Developer
505: .seealso: PCCompositeGetPC()
506: @*/
507: PetscErrorCode PCCompositeGetNumberPC(PC pc,PetscInt *num)
508: {
511: PetscUseMethod(pc,"PCCompositeGetNumberPC_C",(PC,PetscInt*),(pc,num));
512: return 0;
513: }
515: /*@
516: PCCompositeGetPC - Gets one of the PC objects in the composite PC.
518: Not Collective
520: Input Parameters:
521: + pc - the preconditioner context
522: - n - the number of the pc requested
524: Output Parameters:
525: . subpc - the PC requested
527: Level: Developer
529: Notes:
530: To use a different operator to construct one of the inner preconditioners first call PCCompositeGetPC(), then
531: call PCSetOperators() on that PC.
533: .seealso: PCCompositeAddPCType(), PCCompositeGetNumberPC(), PCSetOperators()
534: @*/
535: PetscErrorCode PCCompositeGetPC(PC pc,PetscInt n,PC *subpc)
536: {
539: PetscUseMethod(pc,"PCCompositeGetPC_C",(PC,PetscInt,PC*),(pc,n,subpc));
540: return 0;
541: }
543: /* -------------------------------------------------------------------------------------------*/
545: /*MC
546: PCCOMPOSITE - Build a preconditioner by composing together several preconditioners
548: Options Database Keys:
549: + -pc_composite_type <type: one of multiplicative, additive, symmetric_multiplicative, special> - Sets composite preconditioner type
550: . -pc_use_amat - activates PCSetUseAmat()
551: - -pc_composite_pcs - <pc0,pc1,...> list of PCs to compose
553: Level: intermediate
555: Notes:
556: To use a Krylov method inside the composite preconditioner, set the PCType of one or more
557: inner PCs to be PCKSP.
558: Using a Krylov method inside another Krylov method can be dangerous (you get divergence or
559: the incorrect answer) unless you use KSPFGMRES as the outer Krylov method
560: To use a different operator to construct one of the inner preconditioners first call PCCompositeGetPC(), then
561: call PCSetOperators() on that PC.
563: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
564: PCSHELL, PCKSP, PCCompositeSetType(), PCCompositeSpecialSetAlpha(), PCCompositeAddPCType(),
565: PCCompositeGetPC(), PCSetUseAmat()
567: M*/
569: PETSC_EXTERN PetscErrorCode PCCreate_Composite(PC pc)
570: {
571: PC_Composite *jac;
573: PetscNewLog(pc,&jac);
575: pc->ops->apply = PCApply_Composite_Additive;
576: pc->ops->applytranspose = PCApplyTranspose_Composite_Additive;
577: pc->ops->setup = PCSetUp_Composite;
578: pc->ops->reset = PCReset_Composite;
579: pc->ops->destroy = PCDestroy_Composite;
580: pc->ops->setfromoptions = PCSetFromOptions_Composite;
581: pc->ops->view = PCView_Composite;
582: pc->ops->applyrichardson = NULL;
584: pc->data = (void*)jac;
585: jac->type = PC_COMPOSITE_ADDITIVE;
586: jac->work1 = NULL;
587: jac->work2 = NULL;
588: jac->head = NULL;
590: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeSetType_C",PCCompositeSetType_Composite);
591: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetType_C",PCCompositeGetType_Composite);
592: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeAddPCType_C",PCCompositeAddPCType_Composite);
593: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeAddPC_C",PCCompositeAddPC_Composite);
594: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetNumberPC_C",PCCompositeGetNumberPC_Composite);
595: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetPC_C",PCCompositeGetPC_Composite);
596: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeSpecialSetAlpha_C",PCCompositeSpecialSetAlpha_Composite);
597: return 0;
598: }