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: PetscErrorCode ierr;
26: PC_Composite *jac = (PC_Composite*)pc->data;
27: PC_CompositeLink next = jac->head;
28: Mat mat = pc->pmat;
32: if (!next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPCType() or -pc_composite_pcs");
34: /* Set the reuse flag on children PCs */
35: while (next) {
36: PCSetReusePreconditioner(next->pc,pc->reusepreconditioner);
37: next = next->next;
38: }
39: next = jac->head;
41: if (next->next && !jac->work2) { /* allocate second work vector */
42: VecDuplicate(jac->work1,&jac->work2);
43: }
44: if (pc->useAmat) mat = pc->mat;
45: PCApply(next->pc,x,y); /* y <- B x */
46: while (next->next) {
47: next = next->next;
48: MatMult(mat,y,jac->work1); /* work1 <- A y */
49: VecWAXPY(jac->work2,-1.0,jac->work1,x); /* work2 <- x - work1 */
50: PCApply(next->pc,jac->work2,jac->work1); /* work1 <- C work2 */
51: VecAXPY(y,1.0,jac->work1); /* y <- y + work1 = B x + C (x - A B x) = (B + C (1 - A B)) x */
52: }
53: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
54: while (next->previous) {
55: next = next->previous;
56: MatMult(mat,y,jac->work1);
57: VecWAXPY(jac->work2,-1.0,jac->work1,x);
58: PCApply(next->pc,jac->work2,jac->work1);
59: VecAXPY(y,1.0,jac->work1);
60: }
61: }
62: return(0);
63: }
65: static PetscErrorCode PCApplyTranspose_Composite_Multiplicative(PC pc,Vec x,Vec y)
66: {
67: PetscErrorCode ierr;
68: PC_Composite *jac = (PC_Composite*)pc->data;
69: PC_CompositeLink next = jac->head;
70: Mat mat = pc->pmat;
73: if (!next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPCType() or -pc_composite_pcs");
74: if (next->next && !jac->work2) { /* allocate second work vector */
75: VecDuplicate(jac->work1,&jac->work2);
76: }
77: if (pc->useAmat) mat = pc->mat;
78: /* locate last PC */
79: while (next->next) {
80: next = next->next;
81: }
82: PCApplyTranspose(next->pc,x,y);
83: while (next->previous) {
84: next = next->previous;
85: MatMultTranspose(mat,y,jac->work1);
86: VecWAXPY(jac->work2,-1.0,jac->work1,x);
87: PCApplyTranspose(next->pc,jac->work2,jac->work1);
88: VecAXPY(y,1.0,jac->work1);
89: }
90: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
91: next = jac->head;
92: while (next->next) {
93: next = next->next;
94: MatMultTranspose(mat,y,jac->work1);
95: VecWAXPY(jac->work2,-1.0,jac->work1,x);
96: PCApplyTranspose(next->pc,jac->work2,jac->work1);
97: VecAXPY(y,1.0,jac->work1);
98: }
99: }
100: return(0);
101: }
103: /*
104: This is very special for a matrix of the form alpha I + R + S
105: where first preconditioner is built from alpha I + S and second from
106: alpha I + R
107: */
108: static PetscErrorCode PCApply_Composite_Special(PC pc,Vec x,Vec y)
109: {
110: PetscErrorCode ierr;
111: PC_Composite *jac = (PC_Composite*)pc->data;
112: PC_CompositeLink next = jac->head;
115: if (!next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPCType() or -pc_composite_pcs");
116: if (!next->next || next->next->next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Special composite preconditioners requires exactly two PCs");
118: /* Set the reuse flag on children PCs */
119: PCSetReusePreconditioner(next->pc,pc->reusepreconditioner);
120: PCSetReusePreconditioner(next->next->pc,pc->reusepreconditioner);
122: PCApply(next->pc,x,jac->work1);
123: PCApply(next->next->pc,jac->work1,y);
124: return(0);
125: }
127: static PetscErrorCode PCApply_Composite_Additive(PC pc,Vec x,Vec y)
128: {
129: PetscErrorCode ierr;
130: PC_Composite *jac = (PC_Composite*)pc->data;
131: PC_CompositeLink next = jac->head;
134: if (!next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPCType() or -pc_composite_pcs");
136: /* Set the reuse flag on children PCs */
137: while (next) {
138: PCSetReusePreconditioner(next->pc,pc->reusepreconditioner);
139: next = next->next;
140: }
141: next = jac->head;
143: PCApply(next->pc,x,y);
144: while (next->next) {
145: next = next->next;
146: PCApply(next->pc,x,jac->work1);
147: VecAXPY(y,1.0,jac->work1);
148: }
149: return(0);
150: }
152: static PetscErrorCode PCApplyTranspose_Composite_Additive(PC pc,Vec x,Vec y)
153: {
154: PetscErrorCode ierr;
155: PC_Composite *jac = (PC_Composite*)pc->data;
156: PC_CompositeLink next = jac->head;
159: if (!next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPCType() or -pc_composite_pcs");
160: PCApplyTranspose(next->pc,x,y);
161: while (next->next) {
162: next = next->next;
163: PCApplyTranspose(next->pc,x,jac->work1);
164: VecAXPY(y,1.0,jac->work1);
165: }
166: return(0);
167: }
169: static PetscErrorCode PCSetUp_Composite(PC pc)
170: {
171: PetscErrorCode ierr;
172: PC_Composite *jac = (PC_Composite*)pc->data;
173: PC_CompositeLink next = jac->head;
174: DM dm;
177: if (!jac->work1) {
178: MatCreateVecs(pc->pmat,&jac->work1,NULL);
179: }
180: PCGetDM(pc,&dm);
181: while (next) {
182: if (!next->pc->dm) {
183: PCSetDM(next->pc,dm);
184: }
185: if (!next->pc->mat) {
186: PCSetOperators(next->pc,pc->mat,pc->pmat);
187: }
188: next = next->next;
189: }
190: return(0);
191: }
193: static PetscErrorCode PCReset_Composite(PC pc)
194: {
195: PC_Composite *jac = (PC_Composite*)pc->data;
196: PetscErrorCode ierr;
197: PC_CompositeLink next = jac->head;
200: while (next) {
201: PCReset(next->pc);
202: next = next->next;
203: }
204: VecDestroy(&jac->work1);
205: VecDestroy(&jac->work2);
206: return(0);
207: }
209: static PetscErrorCode PCDestroy_Composite(PC pc)
210: {
211: PC_Composite *jac = (PC_Composite*)pc->data;
212: PetscErrorCode ierr;
213: PC_CompositeLink next = jac->head,next_tmp;
216: PCReset_Composite(pc);
217: while (next) {
218: PCDestroy(&next->pc);
219: next_tmp = next;
220: next = next->next;
221: PetscFree(next_tmp);
222: }
223: PetscFree(pc->data);
224: return(0);
225: }
227: static PetscErrorCode PCSetFromOptions_Composite(PetscOptionItems *PetscOptionsObject,PC pc)
228: {
229: PC_Composite *jac = (PC_Composite*)pc->data;
230: PetscErrorCode ierr;
231: PetscInt nmax = 8,i;
232: PC_CompositeLink next;
233: char *pcs[8];
234: PetscBool flg;
237: PetscOptionsHead(PetscOptionsObject,"Composite preconditioner options");
238: PetscOptionsEnum("-pc_composite_type","Type of composition","PCCompositeSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&jac->type,&flg);
239: if (flg) {
240: PCCompositeSetType(pc,jac->type);
241: }
242: PetscOptionsStringArray("-pc_composite_pcs","List of composite solvers","PCCompositeAddPCType",pcs,&nmax,&flg);
243: if (flg) {
244: for (i=0; i<nmax; i++) {
245: PCCompositeAddPCType(pc,pcs[i]);
246: PetscFree(pcs[i]); /* deallocate string pcs[i], which is allocated in PetscOptionsStringArray() */
247: }
248: }
249: PetscOptionsTail();
251: next = jac->head;
252: while (next) {
253: PCSetFromOptions(next->pc);
254: next = next->next;
255: }
256: return(0);
257: }
259: static PetscErrorCode PCView_Composite(PC pc,PetscViewer viewer)
260: {
261: PC_Composite *jac = (PC_Composite*)pc->data;
262: PetscErrorCode ierr;
263: PC_CompositeLink next = jac->head;
264: PetscBool iascii;
267: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
268: if (iascii) {
269: PetscViewerASCIIPrintf(viewer,"Composite PC type - %s\n",PCCompositeTypes[jac->type]);
270: PetscViewerASCIIPrintf(viewer,"PCs on composite preconditioner follow\n");
271: PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
272: }
273: if (iascii) {
274: PetscViewerASCIIPushTab(viewer);
275: }
276: while (next) {
277: PCView(next->pc,viewer);
278: next = next->next;
279: }
280: if (iascii) {
281: PetscViewerASCIIPopTab(viewer);
282: PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
283: }
284: return(0);
285: }
287: /* ------------------------------------------------------------------------------*/
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: }
298: static PetscErrorCode PCCompositeSetType_Composite(PC pc,PCCompositeType type)
299: {
300: PC_Composite *jac = (PC_Composite*)pc->data;
303: if (type == PC_COMPOSITE_ADDITIVE) {
304: pc->ops->apply = PCApply_Composite_Additive;
305: pc->ops->applytranspose = PCApplyTranspose_Composite_Additive;
306: } else if (type == PC_COMPOSITE_MULTIPLICATIVE || type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
307: pc->ops->apply = PCApply_Composite_Multiplicative;
308: pc->ops->applytranspose = PCApplyTranspose_Composite_Multiplicative;
309: } else if (type == PC_COMPOSITE_SPECIAL) {
310: pc->ops->apply = PCApply_Composite_Special;
311: pc->ops->applytranspose = NULL;
312: } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Unknown composite preconditioner type");
313: jac->type = type;
314: return(0);
315: }
317: static PetscErrorCode PCCompositeGetType_Composite(PC pc,PCCompositeType *type)
318: {
319: PC_Composite *jac = (PC_Composite*)pc->data;
322: *type = jac->type;
323: return(0);
324: }
326: static PetscErrorCode PCCompositeAddPC_Composite(PC pc, PC subpc)
327: {
328: PC_Composite *jac;
329: PC_CompositeLink next, ilink;
330: PetscInt cnt = 0;
331: const char *prefix;
332: char newprefix[20];
333: PetscErrorCode ierr;
336: PetscNewLog(pc, &ilink);
337: ilink->next = NULL;
338: ilink->pc = subpc;
340: jac = (PC_Composite *) pc->data;
341: next = jac->head;
342: if (!next) {
343: jac->head = ilink;
344: ilink->previous = NULL;
345: } else {
346: cnt++;
347: while (next->next) {
348: next = next->next;
349: cnt++;
350: }
351: next->next = ilink;
352: ilink->previous = next;
353: }
354: PCGetOptionsPrefix(pc, &prefix);
355: PCSetOptionsPrefix(subpc, prefix);
356: PetscSNPrintf(newprefix, 20, "sub_%d_", (int) cnt);
357: PCAppendOptionsPrefix(subpc, newprefix);
358: PetscObjectReference((PetscObject) subpc);
359: return(0);
360: }
362: static PetscErrorCode PCCompositeAddPCType_Composite(PC pc, PCType type)
363: {
364: PC subpc;
368: PCCreate(PetscObjectComm((PetscObject)pc), &subpc);
369: PetscObjectIncrementTabLevel((PetscObject) subpc, (PetscObject) pc, 1);
370: PetscLogObjectParent((PetscObject) pc, (PetscObject) subpc);
371: PCCompositeAddPC_Composite(pc, subpc);
372: /* type is set after prefix, because some methods may modify prefix, e.g. pcksp */
373: PCSetType(subpc, type);
374: PCDestroy(&subpc);
375: return(0);
376: }
378: static PetscErrorCode PCCompositeGetNumberPC_Composite(PC pc,PetscInt *n)
379: {
380: PC_Composite *jac;
381: PC_CompositeLink next;
384: jac = (PC_Composite*)pc->data;
385: next = jac->head;
386: *n = 0;
387: while (next) {
388: next = next->next;
389: (*n) ++;
390: }
391: return(0);
392: }
394: static PetscErrorCode PCCompositeGetPC_Composite(PC pc,PetscInt n,PC *subpc)
395: {
396: PC_Composite *jac;
397: PC_CompositeLink next;
398: PetscInt i;
401: jac = (PC_Composite*)pc->data;
402: next = jac->head;
403: for (i=0; i<n; i++) {
404: if (!next->next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Not enough PCs in composite preconditioner");
405: next = next->next;
406: }
407: *subpc = next->pc;
408: return(0);
409: }
411: /* -------------------------------------------------------------------------------- */
412: /*@
413: PCCompositeSetType - Sets the type of composite preconditioner.
415: Logically Collective on PC
417: Input Parameters:
418: + pc - the preconditioner context
419: - type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL
421: Options Database Key:
422: . -pc_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type
424: Level: Developer
426: @*/
427: PetscErrorCode PCCompositeSetType(PC pc,PCCompositeType type)
428: {
434: PetscTryMethod(pc,"PCCompositeSetType_C",(PC,PCCompositeType),(pc,type));
435: return(0);
436: }
438: /*@
439: PCCompositeGetType - Gets the type of composite preconditioner.
441: Logically Collective on PC
443: Input Parameter:
444: . pc - the preconditioner context
446: Output Parameter:
447: . type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL
449: Options Database Key:
450: . -pc_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type
452: Level: Developer
454: @*/
455: PetscErrorCode PCCompositeGetType(PC pc,PCCompositeType *type)
456: {
461: PetscUseMethod(pc,"PCCompositeGetType_C",(PC,PCCompositeType*),(pc,type));
462: return(0);
463: }
465: /*@
466: PCCompositeSpecialSetAlpha - Sets alpha for the special composite preconditioner
467: for alphaI + R + S
469: Logically Collective on PC
471: Input Parameters:
472: + pc - the preconditioner context
473: - alpha - scale on identity
475: Level: Developer
477: @*/
478: PetscErrorCode PCCompositeSpecialSetAlpha(PC pc,PetscScalar alpha)
479: {
485: PetscTryMethod(pc,"PCCompositeSpecialSetAlpha_C",(PC,PetscScalar),(pc,alpha));
486: return(0);
487: }
489: /*@C
490: PCCompositeAddPCType - Adds another PC of the given type to the composite PC.
492: Collective on PC
494: Input Parameters:
495: + pc - the preconditioner context
496: - type - the type of the new preconditioner
498: Level: Developer
500: .seealso: PCCompositeAddPC()
501: @*/
502: PetscErrorCode PCCompositeAddPCType(PC pc,PCType type)
503: {
508: PetscTryMethod(pc,"PCCompositeAddPCType_C",(PC,PCType),(pc,type));
509: return(0);
510: }
512: /*@
513: PCCompositeAddPC - Adds another PC to the composite PC.
515: Collective on PC
517: Input Parameters:
518: + pc - the preconditioner context
519: - subpc - the new preconditioner
521: Level: Developer
523: .seealso: PCCompositeAddPCType()
524: @*/
525: PetscErrorCode PCCompositeAddPC(PC pc, PC subpc)
526: {
532: PetscTryMethod(pc,"PCCompositeAddPC_C",(PC,PC),(pc,subpc));
533: return(0);
534: }
536: /*@
537: PCCompositeGetNumberPC - Gets the number of PC objects in the composite PC.
539: Not Collective
541: Input Parameter:
542: . pc - the preconditioner context
544: Output Parameter:
545: . num - the number of sub pcs
547: Level: Developer
549: .seealso: PCCompositeGetPC()
550: @*/
551: PetscErrorCode PCCompositeGetNumberPC(PC pc,PetscInt *num)
552: {
558: PetscUseMethod(pc,"PCCompositeGetNumberPC_C",(PC,PetscInt*),(pc,num));
559: return(0);
560: }
562: /*@
563: PCCompositeGetPC - Gets one of the PC objects in the composite PC.
565: Not Collective
567: Input Parameters:
568: + pc - the preconditioner context
569: - n - the number of the pc requested
571: Output Parameters:
572: . subpc - the PC requested
574: Level: Developer
576: Notes:
577: To use a different operator to construct one of the inner preconditioners first call PCCompositeGetPC(), then
578: call PCSetOperators() on that PC.
580: .seealso: PCCompositeAddPCType(), PCCompositeGetNumberPC(), PCSetOperators()
581: @*/
582: PetscErrorCode PCCompositeGetPC(PC pc,PetscInt n,PC *subpc)
583: {
589: PetscUseMethod(pc,"PCCompositeGetPC_C",(PC,PetscInt,PC*),(pc,n,subpc));
590: return(0);
591: }
593: /* -------------------------------------------------------------------------------------------*/
595: /*MC
596: PCCOMPOSITE - Build a preconditioner by composing together several preconditioners
598: Options Database Keys:
599: + -pc_composite_type <type: one of multiplicative, additive, symmetric_multiplicative, special> - Sets composite preconditioner type
600: . -pc_use_amat - activates PCSetUseAmat()
601: - -pc_composite_pcs - <pc0,pc1,...> list of PCs to compose
603: Level: intermediate
605: Notes:
606: To use a Krylov method inside the composite preconditioner, set the PCType of one or more
607: inner PCs to be PCKSP.
608: Using a Krylov method inside another Krylov method can be dangerous (you get divergence or
609: the incorrect answer) unless you use KSPFGMRES as the outer Krylov method
610: To use a different operator to construct one of the inner preconditioners first call PCCompositeGetPC(), then
611: call PCSetOperators() on that PC.
613: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
614: PCSHELL, PCKSP, PCCompositeSetType(), PCCompositeSpecialSetAlpha(), PCCompositeAddPCType(),
615: PCCompositeGetPC(), PCSetUseAmat()
617: M*/
619: PETSC_EXTERN PetscErrorCode PCCreate_Composite(PC pc)
620: {
622: PC_Composite *jac;
625: PetscNewLog(pc,&jac);
627: pc->ops->apply = PCApply_Composite_Additive;
628: pc->ops->applytranspose = PCApplyTranspose_Composite_Additive;
629: pc->ops->setup = PCSetUp_Composite;
630: pc->ops->reset = PCReset_Composite;
631: pc->ops->destroy = PCDestroy_Composite;
632: pc->ops->setfromoptions = PCSetFromOptions_Composite;
633: pc->ops->view = PCView_Composite;
634: pc->ops->applyrichardson = NULL;
636: pc->data = (void*)jac;
637: jac->type = PC_COMPOSITE_ADDITIVE;
638: jac->work1 = NULL;
639: jac->work2 = NULL;
640: jac->head = NULL;
642: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeSetType_C",PCCompositeSetType_Composite);
643: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetType_C",PCCompositeGetType_Composite);
644: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeAddPCType_C",PCCompositeAddPCType_Composite);
645: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeAddPC_C",PCCompositeAddPC_Composite);
646: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetNumberPC_C",PCCompositeGetNumberPC_Composite);
647: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetPC_C",PCCompositeGetPC_Composite);
648: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeSpecialSetAlpha_C",PCCompositeSpecialSetAlpha_Composite);
649: return(0);
650: }