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) next = next->next;
76: PCApplyTranspose(next->pc, x, y);
77: while (next->previous) {
78: next = next->previous;
79: MatMultTranspose(mat, y, jac->work1);
80: VecWAXPY(jac->work2, -1.0, jac->work1, x);
81: PCApplyTranspose(next->pc, jac->work2, jac->work1);
82: VecAXPY(y, 1.0, jac->work1);
83: }
84: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
85: next = jac->head;
86: while (next->next) {
87: next = next->next;
88: MatMultTranspose(mat, y, jac->work1);
89: VecWAXPY(jac->work2, -1.0, jac->work1, x);
90: PCApplyTranspose(next->pc, jac->work2, jac->work1);
91: VecAXPY(y, 1.0, jac->work1);
92: }
93: }
94: return 0;
95: }
97: /*
98: This is very special for a matrix of the form alpha I + R + S
99: where first preconditioner is built from alpha I + S and second from
100: alpha I + R
101: */
102: static PetscErrorCode PCApply_Composite_Special(PC pc, Vec x, Vec y)
103: {
104: PC_Composite *jac = (PC_Composite *)pc->data;
105: PC_CompositeLink next = jac->head;
110: /* Set the reuse flag on children PCs */
111: PCSetReusePreconditioner(next->pc, pc->reusepreconditioner);
112: PCSetReusePreconditioner(next->next->pc, pc->reusepreconditioner);
114: PCApply(next->pc, x, jac->work1);
115: PCApply(next->next->pc, jac->work1, y);
116: return 0;
117: }
119: static PetscErrorCode PCApply_Composite_Additive(PC pc, Vec x, Vec y)
120: {
121: PC_Composite *jac = (PC_Composite *)pc->data;
122: PC_CompositeLink next = jac->head;
126: /* Set the reuse flag on children PCs */
127: while (next) {
128: PCSetReusePreconditioner(next->pc, pc->reusepreconditioner);
129: next = next->next;
130: }
131: next = jac->head;
133: PCApply(next->pc, x, y);
134: while (next->next) {
135: next = next->next;
136: PCApply(next->pc, x, jac->work1);
137: VecAXPY(y, 1.0, jac->work1);
138: }
139: return 0;
140: }
142: static PetscErrorCode PCApplyTranspose_Composite_Additive(PC pc, Vec x, Vec y)
143: {
144: PC_Composite *jac = (PC_Composite *)pc->data;
145: PC_CompositeLink next = jac->head;
148: PCApplyTranspose(next->pc, x, y);
149: while (next->next) {
150: next = next->next;
151: PCApplyTranspose(next->pc, x, jac->work1);
152: VecAXPY(y, 1.0, jac->work1);
153: }
154: return 0;
155: }
157: static PetscErrorCode PCSetUp_Composite(PC pc)
158: {
159: PC_Composite *jac = (PC_Composite *)pc->data;
160: PC_CompositeLink next = jac->head;
161: DM dm;
163: if (!jac->work1) MatCreateVecs(pc->pmat, &jac->work1, NULL);
164: PCGetDM(pc, &dm);
165: while (next) {
166: if (!next->pc->dm) PCSetDM(next->pc, dm);
167: if (!next->pc->mat) PCSetOperators(next->pc, pc->mat, pc->pmat);
168: next = next->next;
169: }
170: return 0;
171: }
173: static PetscErrorCode PCReset_Composite(PC pc)
174: {
175: PC_Composite *jac = (PC_Composite *)pc->data;
176: PC_CompositeLink next = jac->head;
178: while (next) {
179: PCReset(next->pc);
180: next = next->next;
181: }
182: VecDestroy(&jac->work1);
183: VecDestroy(&jac->work2);
184: return 0;
185: }
187: static PetscErrorCode PCDestroy_Composite(PC pc)
188: {
189: PC_Composite *jac = (PC_Composite *)pc->data;
190: PC_CompositeLink next = jac->head, next_tmp;
192: PCReset_Composite(pc);
193: while (next) {
194: PCDestroy(&next->pc);
195: next_tmp = next;
196: next = next->next;
197: PetscFree(next_tmp);
198: }
199: PetscObjectComposeFunction((PetscObject)pc, "PCCompositeSetType_C", NULL);
200: PetscObjectComposeFunction((PetscObject)pc, "PCCompositeGetType_C", NULL);
201: PetscObjectComposeFunction((PetscObject)pc, "PCCompositeAddPCType_C", NULL);
202: PetscObjectComposeFunction((PetscObject)pc, "PCCompositeAddPC_C", NULL);
203: PetscObjectComposeFunction((PetscObject)pc, "PCCompositeGetNumberPC_C", NULL);
204: PetscObjectComposeFunction((PetscObject)pc, "PCCompositeGetPC_C", NULL);
205: PetscObjectComposeFunction((PetscObject)pc, "PCCompositeSpecialSetAlpha_C", NULL);
206: PetscFree(pc->data);
207: return 0;
208: }
210: static PetscErrorCode PCSetFromOptions_Composite(PC pc, PetscOptionItems *PetscOptionsObject)
211: {
212: PC_Composite *jac = (PC_Composite *)pc->data;
213: PetscInt nmax = 8, i;
214: PC_CompositeLink next;
215: char *pcs[8];
216: PetscBool flg;
218: PetscOptionsHeadBegin(PetscOptionsObject, "Composite preconditioner options");
219: PetscOptionsEnum("-pc_composite_type", "Type of composition", "PCCompositeSetType", PCCompositeTypes, (PetscEnum)jac->type, (PetscEnum *)&jac->type, &flg);
220: if (flg) PCCompositeSetType(pc, jac->type);
221: PetscOptionsStringArray("-pc_composite_pcs", "List of composite solvers", "PCCompositeAddPCType", pcs, &nmax, &flg);
222: if (flg) {
223: for (i = 0; i < nmax; i++) {
224: PCCompositeAddPCType(pc, pcs[i]);
225: PetscFree(pcs[i]); /* deallocate string pcs[i], which is allocated in PetscOptionsStringArray() */
226: }
227: }
228: PetscOptionsHeadEnd();
230: next = jac->head;
231: while (next) {
232: PCSetFromOptions(next->pc);
233: next = next->next;
234: }
235: return 0;
236: }
238: static PetscErrorCode PCView_Composite(PC pc, PetscViewer viewer)
239: {
240: PC_Composite *jac = (PC_Composite *)pc->data;
241: PC_CompositeLink next = jac->head;
242: PetscBool iascii;
244: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
245: if (iascii) {
246: PetscViewerASCIIPrintf(viewer, "Composite PC type - %s\n", PCCompositeTypes[jac->type]);
247: PetscViewerASCIIPrintf(viewer, "PCs on composite preconditioner follow\n");
248: PetscViewerASCIIPrintf(viewer, "---------------------------------\n");
249: }
250: if (iascii) PetscViewerASCIIPushTab(viewer);
251: while (next) {
252: PCView(next->pc, viewer);
253: next = next->next;
254: }
255: if (iascii) {
256: PetscViewerASCIIPopTab(viewer);
257: PetscViewerASCIIPrintf(viewer, "---------------------------------\n");
258: }
259: return 0;
260: }
262: static PetscErrorCode PCCompositeSpecialSetAlpha_Composite(PC pc, PetscScalar alpha)
263: {
264: PC_Composite *jac = (PC_Composite *)pc->data;
266: jac->alpha = alpha;
267: return 0;
268: }
270: static PetscErrorCode PCCompositeSetType_Composite(PC pc, PCCompositeType type)
271: {
272: PC_Composite *jac = (PC_Composite *)pc->data;
274: if (type == PC_COMPOSITE_ADDITIVE) {
275: pc->ops->apply = PCApply_Composite_Additive;
276: pc->ops->applytranspose = PCApplyTranspose_Composite_Additive;
277: } else if (type == PC_COMPOSITE_MULTIPLICATIVE || type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
278: pc->ops->apply = PCApply_Composite_Multiplicative;
279: pc->ops->applytranspose = PCApplyTranspose_Composite_Multiplicative;
280: } else if (type == PC_COMPOSITE_SPECIAL) {
281: pc->ops->apply = PCApply_Composite_Special;
282: pc->ops->applytranspose = NULL;
283: } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Unknown composite preconditioner type");
284: jac->type = type;
285: return 0;
286: }
288: static PetscErrorCode PCCompositeGetType_Composite(PC pc, PCCompositeType *type)
289: {
290: PC_Composite *jac = (PC_Composite *)pc->data;
292: *type = jac->type;
293: return 0;
294: }
296: static PetscErrorCode PCCompositeAddPC_Composite(PC pc, PC subpc)
297: {
298: PC_Composite *jac;
299: PC_CompositeLink next, ilink;
300: PetscInt cnt = 0;
301: const char *prefix;
302: char newprefix[20];
304: PetscNew(&ilink);
305: ilink->next = NULL;
306: ilink->pc = subpc;
308: jac = (PC_Composite *)pc->data;
309: next = jac->head;
310: if (!next) {
311: jac->head = ilink;
312: ilink->previous = NULL;
313: } else {
314: cnt++;
315: while (next->next) {
316: next = next->next;
317: cnt++;
318: }
319: next->next = ilink;
320: ilink->previous = next;
321: }
322: PCGetOptionsPrefix(pc, &prefix);
323: PCSetOptionsPrefix(subpc, prefix);
324: PetscSNPrintf(newprefix, 20, "sub_%d_", (int)cnt);
325: PCAppendOptionsPrefix(subpc, newprefix);
326: PetscObjectReference((PetscObject)subpc);
327: return 0;
328: }
330: static PetscErrorCode PCCompositeAddPCType_Composite(PC pc, PCType type)
331: {
332: PC subpc;
334: PCCreate(PetscObjectComm((PetscObject)pc), &subpc);
335: PetscObjectIncrementTabLevel((PetscObject)subpc, (PetscObject)pc, 1);
336: PCCompositeAddPC_Composite(pc, subpc);
337: /* type is set after prefix, because some methods may modify prefix, e.g. pcksp */
338: PCSetType(subpc, type);
339: PCDestroy(&subpc);
340: return 0;
341: }
343: static PetscErrorCode PCCompositeGetNumberPC_Composite(PC pc, PetscInt *n)
344: {
345: PC_Composite *jac;
346: PC_CompositeLink next;
348: jac = (PC_Composite *)pc->data;
349: next = jac->head;
350: *n = 0;
351: while (next) {
352: next = next->next;
353: (*n)++;
354: }
355: return 0;
356: }
358: static PetscErrorCode PCCompositeGetPC_Composite(PC pc, PetscInt n, PC *subpc)
359: {
360: PC_Composite *jac;
361: PC_CompositeLink next;
362: PetscInt i;
364: jac = (PC_Composite *)pc->data;
365: next = jac->head;
366: for (i = 0; i < n; i++) {
368: next = next->next;
369: }
370: *subpc = next->pc;
371: return 0;
372: }
374: /*@
375: PCCompositeSetType - Sets the type of composite preconditioner.
377: Logically Collective
379: Input Parameters:
380: + pc - the preconditioner context
381: - type - `PC_COMPOSITE_ADDITIVE` (default), `PC_COMPOSITE_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`
383: Options Database Key:
384: . -pc_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type
386: Level: advanced
388: .seealso: `PCCOMPOSITE`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PCCompositeType`,
389: `PCCompositeGetType()`
390: @*/
391: PetscErrorCode PCCompositeSetType(PC pc, PCCompositeType type)
392: {
395: PetscTryMethod(pc, "PCCompositeSetType_C", (PC, PCCompositeType), (pc, type));
396: return 0;
397: }
399: /*@
400: PCCompositeGetType - Gets the type of composite preconditioner.
402: Logically Collective
404: Input Parameter:
405: . pc - the preconditioner context
407: Output Parameter:
408: . type - `PC_COMPOSITE_ADDITIVE` (default), `PC_COMPOSITE_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`
410: Level: advanced
412: .seealso: `PCCOMPOSITE`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PCCompositeType`,
413: `PCCompositeSetType()`
414: @*/
415: PetscErrorCode PCCompositeGetType(PC pc, PCCompositeType *type)
416: {
418: PetscUseMethod(pc, "PCCompositeGetType_C", (PC, PCCompositeType *), (pc, type));
419: return 0;
420: }
422: /*@
423: PCCompositeSpecialSetAlpha - Sets alpha for the special composite preconditioner, `PC_COMPOSITE_SPECIAL`,
424: for alphaI + R + S
426: Logically Collective
428: Input Parameters:
429: + pc - the preconditioner context
430: - alpha - scale on identity
432: Level: Developer
434: .seealso: `PCCOMPOSITE`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PCCompositeType`,
435: `PCCompositeSetType()`, `PCCompositeGetType()`
436: @*/
437: PetscErrorCode PCCompositeSpecialSetAlpha(PC pc, PetscScalar alpha)
438: {
441: PetscTryMethod(pc, "PCCompositeSpecialSetAlpha_C", (PC, PetscScalar), (pc, alpha));
442: return 0;
443: }
445: /*@C
446: PCCompositeAddPCType - Adds another `PC` of the given type to the composite `PC`.
448: Collective
450: Input Parameters:
451: + pc - the preconditioner context
452: - type - the type of the new preconditioner
454: Level: intermediate
456: .seealso: `PCCOMPOSITE`, `PCCompositeAddPC()`, `PCCompositeGetNumberPC()`
457: @*/
458: PetscErrorCode PCCompositeAddPCType(PC pc, PCType type)
459: {
461: PetscTryMethod(pc, "PCCompositeAddPCType_C", (PC, PCType), (pc, type));
462: return 0;
463: }
465: /*@
466: PCCompositeAddPC - Adds another `PC` to the composite `PC`.
468: Collective
470: Input Parameters:
471: + pc - the preconditioner context
472: - subpc - the new preconditioner
474: Level: intermediate
476: .seealso: `PCCOMPOSITE`, `PCCompositeAddPCType()`, `PCCompositeGetNumberPC()`
477: @*/
478: PetscErrorCode PCCompositeAddPC(PC pc, PC subpc)
479: {
482: PetscTryMethod(pc, "PCCompositeAddPC_C", (PC, PC), (pc, subpc));
483: return 0;
484: }
486: /*@
487: PCCompositeGetNumberPC - Gets the number of `PC` objects in the composite `PC`.
489: Not Collective
491: Input Parameter:
492: . pc - the preconditioner context
494: Output Parameter:
495: . num - the number of sub pcs
497: Level: Developer
499: .seealso: `PCCOMPOSITE`, `PCCompositeGetPC()`, `PCCompositeAddPC()`, `PCCompositeAddPCType()`
500: @*/
501: PetscErrorCode PCCompositeGetNumberPC(PC pc, PetscInt *num)
502: {
505: PetscUseMethod(pc, "PCCompositeGetNumberPC_C", (PC, PetscInt *), (pc, num));
506: return 0;
507: }
509: /*@
510: PCCompositeGetPC - Gets one of the `PC` objects in the composite `PC`.
512: Not Collective
514: Input Parameters:
515: + pc - the preconditioner context
516: - n - the number of the pc requested
518: Output Parameter:
519: . subpc - the PC requested
521: Level: intermediate
523: Note:
524: To use a different operator to construct one of the inner preconditioners first call `PCCompositeGetPC()`, then
525: call `PCSetOperators()` on that `PC`.
527: .seealso: `PCCOMPOSITE`, `PCCompositeAddPCType()`, `PCCompositeGetNumberPC()`, `PCSetOperators()`
528: @*/
529: PetscErrorCode PCCompositeGetPC(PC pc, PetscInt n, PC *subpc)
530: {
533: PetscUseMethod(pc, "PCCompositeGetPC_C", (PC, PetscInt, PC *), (pc, n, subpc));
534: return 0;
535: }
537: /*MC
538: PCCOMPOSITE - Build a preconditioner by composing together several preconditioners
540: Options Database Keys:
541: + -pc_composite_type <type: one of multiplicative, additive, symmetric_multiplicative, special> - Sets composite preconditioner type
542: . -pc_use_amat - activates `PCSetUseAmat()`
543: - -pc_composite_pcs - <pc0,pc1,...> list of PCs to compose
545: Level: intermediate
547: Notes:
548: To use a Krylov method inside the composite preconditioner, set the `PCType` of one or more
549: inner `PC`s to be `PCKSP`. Using a Krylov method inside another Krylov method can be dangerous (you get divergence or
550: the incorrect answer) unless you use `KSPFGMRES` as the outer Krylov method
552: To use a different operator to construct one of the inner preconditioners first call `PCCompositeGetPC()`, then
553: call `PCSetOperators()` on that `PC`.
555: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
556: `PCSHELL`, `PCKSP`, `PCCompositeSetType()`, `PCCompositeSpecialSetAlpha()`, `PCCompositeAddPCType()`,
557: `PCCompositeGetPC()`, `PCSetUseAmat()`, `PCCompositeAddPC()`, `PCCompositeGetNumberPC()`
558: M*/
560: PETSC_EXTERN PetscErrorCode PCCreate_Composite(PC pc)
561: {
562: PC_Composite *jac;
564: PetscNew(&jac);
566: pc->ops->apply = PCApply_Composite_Additive;
567: pc->ops->applytranspose = PCApplyTranspose_Composite_Additive;
568: pc->ops->setup = PCSetUp_Composite;
569: pc->ops->reset = PCReset_Composite;
570: pc->ops->destroy = PCDestroy_Composite;
571: pc->ops->setfromoptions = PCSetFromOptions_Composite;
572: pc->ops->view = PCView_Composite;
573: pc->ops->applyrichardson = NULL;
575: pc->data = (void *)jac;
576: jac->type = PC_COMPOSITE_ADDITIVE;
577: jac->work1 = NULL;
578: jac->work2 = NULL;
579: jac->head = NULL;
581: PetscObjectComposeFunction((PetscObject)pc, "PCCompositeSetType_C", PCCompositeSetType_Composite);
582: PetscObjectComposeFunction((PetscObject)pc, "PCCompositeGetType_C", PCCompositeGetType_Composite);
583: PetscObjectComposeFunction((PetscObject)pc, "PCCompositeAddPCType_C", PCCompositeAddPCType_Composite);
584: PetscObjectComposeFunction((PetscObject)pc, "PCCompositeAddPC_C", PCCompositeAddPC_Composite);
585: PetscObjectComposeFunction((PetscObject)pc, "PCCompositeGetNumberPC_C", PCCompositeGetNumberPC_Composite);
586: PetscObjectComposeFunction((PetscObject)pc, "PCCompositeGetPC_C", PCCompositeGetPC_Composite);
587: PetscObjectComposeFunction((PetscObject)pc, "PCCompositeSpecialSetAlpha_C", PCCompositeSpecialSetAlpha_Composite);
588: return 0;
589: }