Actual source code: composite.c
petsc-3.7.7 2017-09-25
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;
34: if (!next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC() or -pc_composite_pcs");
36: /* Set the reuse flag on children PCs */
37: while (next) {
38: PCSetReusePreconditioner(next->pc,pc->reusepreconditioner);
39: next = next->next;
40: }
41: next = jac->head;
43: if (next->next && !jac->work2) { /* allocate second work vector */
44: VecDuplicate(jac->work1,&jac->work2);
45: }
46: if (pc->useAmat) mat = pc->mat;
47: PCApply(next->pc,x,y); /* y <- B x */
48: while (next->next) {
49: next = next->next;
50: MatMult(mat,y,jac->work1); /* work1 <- A y */
51: VecWAXPY(jac->work2,-1.0,jac->work1,x); /* work2 <- x - work1 */
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); /* work1 <- C work2 */
54: VecAXPY(y,1.0,jac->work1); /* y <- y + work1 = B x + C (x - A B x) = (B + C (1 - A B)) x */
55: }
56: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
57: while (next->previous) {
58: next = next->previous;
59: MatMult(mat,y,jac->work1);
60: VecWAXPY(jac->work2,-1.0,jac->work1,x);
61: VecSet(jac->work1,0.0); /* zero since some PC's may not set all entries in the result */
62: PCApply(next->pc,jac->work2,jac->work1);
63: VecAXPY(y,1.0,jac->work1);
64: }
65: }
66: return(0);
67: }
71: static PetscErrorCode PCApplyTranspose_Composite_Multiplicative(PC pc,Vec x,Vec y)
72: {
73: PetscErrorCode ierr;
74: PC_Composite *jac = (PC_Composite*)pc->data;
75: PC_CompositeLink next = jac->head;
76: Mat mat = pc->pmat;
79: if (!next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC() or -pc_composite_pcs");
80: if (next->next && !jac->work2) { /* allocate second work vector */
81: VecDuplicate(jac->work1,&jac->work2);
82: }
83: if (pc->useAmat) mat = pc->mat;
84: /* locate last PC */
85: while (next->next) {
86: next = next->next;
87: }
88: PCApplyTranspose(next->pc,x,y);
89: while (next->previous) {
90: next = next->previous;
91: MatMultTranspose(mat,y,jac->work1);
92: VecWAXPY(jac->work2,-1.0,jac->work1,x);
93: VecSet(jac->work1,0.0); /* zero since some PC's may not set all entries in the result */
94: PCApplyTranspose(next->pc,jac->work2,jac->work1);
95: VecAXPY(y,1.0,jac->work1);
96: }
97: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
98: next = jac->head;
99: while (next->next) {
100: next = next->next;
101: MatMultTranspose(mat,y,jac->work1);
102: VecWAXPY(jac->work2,-1.0,jac->work1,x);
103: VecSet(jac->work1,0.0); /* zero since some PC's may not set all entries in the result */
104: PCApplyTranspose(next->pc,jac->work2,jac->work1);
105: VecAXPY(y,1.0,jac->work1);
106: }
107: }
108: return(0);
109: }
111: /*
112: This is very special for a matrix of the form alpha I + R + S
113: where first preconditioner is built from alpha I + S and second from
114: alpha I + R
115: */
118: static PetscErrorCode PCApply_Composite_Special(PC pc,Vec x,Vec y)
119: {
120: PetscErrorCode ierr;
121: PC_Composite *jac = (PC_Composite*)pc->data;
122: PC_CompositeLink next = jac->head;
125: if (!next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC() or -pc_composite_pcs");
126: if (!next->next || next->next->next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Special composite preconditioners requires exactly two PCs");
128: /* Set the reuse flag on children PCs */
129: PCSetReusePreconditioner(next->pc,pc->reusepreconditioner);
130: PCSetReusePreconditioner(next->next->pc,pc->reusepreconditioner);
132: PCApply(next->pc,x,jac->work1);
133: PCApply(next->next->pc,jac->work1,y);
134: return(0);
135: }
139: static PetscErrorCode PCApply_Composite_Additive(PC pc,Vec x,Vec y)
140: {
141: PetscErrorCode ierr;
142: PC_Composite *jac = (PC_Composite*)pc->data;
143: PC_CompositeLink next = jac->head;
146: if (!next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC() or -pc_composite_pcs");
148: /* Set the reuse flag on children PCs */
149: while (next) {
150: PCSetReusePreconditioner(next->pc,pc->reusepreconditioner);
151: next = next->next;
152: }
153: next = jac->head;
155: PCApply(next->pc,x,y);
156: while (next->next) {
157: next = next->next;
158: VecSet(jac->work1,0.0); /* zero since some PC's may not set all entries in the result */
159: PCApply(next->pc,x,jac->work1);
160: VecAXPY(y,1.0,jac->work1);
161: }
162: return(0);
163: }
167: static PetscErrorCode PCApplyTranspose_Composite_Additive(PC pc,Vec x,Vec y)
168: {
169: PetscErrorCode ierr;
170: PC_Composite *jac = (PC_Composite*)pc->data;
171: PC_CompositeLink next = jac->head;
174: if (!next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC() or -pc_composite_pcs");
175: PCApplyTranspose(next->pc,x,y);
176: while (next->next) {
177: next = next->next;
178: VecSet(jac->work1,0.0); /* zero since some PC's may not set all entries in the result */
179: PCApplyTranspose(next->pc,x,jac->work1);
180: VecAXPY(y,1.0,jac->work1);
181: }
182: return(0);
183: }
187: static PetscErrorCode PCSetUp_Composite(PC pc)
188: {
189: PetscErrorCode ierr;
190: PC_Composite *jac = (PC_Composite*)pc->data;
191: PC_CompositeLink next = jac->head;
192: DM dm;
195: if (!jac->work1) {
196: MatCreateVecs(pc->pmat,&jac->work1,0);
197: }
198: PCGetDM(pc,&dm);
199: while (next) {
200: PCSetDM(next->pc,dm);
201: PCSetOperators(next->pc,pc->mat,pc->pmat);
202: next = next->next;
203: }
204: return(0);
205: }
209: static PetscErrorCode PCReset_Composite(PC pc)
210: {
211: PC_Composite *jac = (PC_Composite*)pc->data;
212: PetscErrorCode ierr;
213: PC_CompositeLink next = jac->head;
216: while (next) {
217: PCReset(next->pc);
218: next = next->next;
219: }
220: VecDestroy(&jac->work1);
221: VecDestroy(&jac->work2);
222: return(0);
223: }
227: static PetscErrorCode PCDestroy_Composite(PC pc)
228: {
229: PC_Composite *jac = (PC_Composite*)pc->data;
230: PetscErrorCode ierr;
231: PC_CompositeLink next = jac->head,next_tmp;
234: PCReset_Composite(pc);
235: while (next) {
236: PCDestroy(&next->pc);
237: next_tmp = next;
238: next = next->next;
239: PetscFree(next_tmp);
240: }
241: PetscFree(pc->data);
242: return(0);
243: }
247: static PetscErrorCode PCSetFromOptions_Composite(PetscOptionItems *PetscOptionsObject,PC pc)
248: {
249: PC_Composite *jac = (PC_Composite*)pc->data;
250: PetscErrorCode ierr;
251: PetscInt nmax = 8,i;
252: PC_CompositeLink next;
253: char *pcs[8];
254: PetscBool flg;
257: PetscOptionsHead(PetscOptionsObject,"Composite preconditioner options");
258: PetscOptionsEnum("-pc_composite_type","Type of composition","PCCompositeSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&jac->type,&flg);
259: if (flg) {
260: PCCompositeSetType(pc,jac->type);
261: }
262: PetscOptionsStringArray("-pc_composite_pcs","List of composite solvers","PCCompositeAddPC",pcs,&nmax,&flg);
263: if (flg) {
264: for (i=0; i<nmax; i++) {
265: PCCompositeAddPC(pc,pcs[i]);
266: PetscFree(pcs[i]); /* deallocate string pcs[i], which is allocated in PetscOptionsStringArray() */
267: }
268: }
269: PetscOptionsTail();
271: next = jac->head;
272: while (next) {
273: PCSetFromOptions(next->pc);
274: next = next->next;
275: }
276: return(0);
277: }
281: static PetscErrorCode PCView_Composite(PC pc,PetscViewer viewer)
282: {
283: PC_Composite *jac = (PC_Composite*)pc->data;
284: PetscErrorCode ierr;
285: PC_CompositeLink next = jac->head;
286: PetscBool iascii;
289: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
290: if (iascii) {
291: PetscViewerASCIIPrintf(viewer,"Composite PC type - %s\n",PCCompositeTypes[jac->type]);
292: PetscViewerASCIIPrintf(viewer,"PCs on composite preconditioner follow\n");
293: PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
294: }
295: if (iascii) {
296: PetscViewerASCIIPushTab(viewer);
297: }
298: while (next) {
299: PCView(next->pc,viewer);
300: next = next->next;
301: }
302: if (iascii) {
303: PetscViewerASCIIPopTab(viewer);
304: PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
305: }
306: return(0);
307: }
309: /* ------------------------------------------------------------------------------*/
313: static PetscErrorCode PCCompositeSpecialSetAlpha_Composite(PC pc,PetscScalar alpha)
314: {
315: PC_Composite *jac = (PC_Composite*)pc->data;
318: jac->alpha = alpha;
319: return(0);
320: }
324: static PetscErrorCode PCCompositeSetType_Composite(PC pc,PCCompositeType type)
325: {
326: PC_Composite *jac = (PC_Composite*)pc->data;
329: if (type == PC_COMPOSITE_ADDITIVE) {
330: pc->ops->apply = PCApply_Composite_Additive;
331: pc->ops->applytranspose = PCApplyTranspose_Composite_Additive;
332: } else if (type == PC_COMPOSITE_MULTIPLICATIVE || type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
333: pc->ops->apply = PCApply_Composite_Multiplicative;
334: pc->ops->applytranspose = PCApplyTranspose_Composite_Multiplicative;
335: } else if (type == PC_COMPOSITE_SPECIAL) {
336: pc->ops->apply = PCApply_Composite_Special;
337: pc->ops->applytranspose = NULL;
338: } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Unkown composite preconditioner type");
339: jac->type = type;
340: return(0);
341: }
345: static PetscErrorCode PCCompositeGetType_Composite(PC pc,PCCompositeType *type)
346: {
347: PC_Composite *jac = (PC_Composite*)pc->data;
350: *type = jac->type;
351: return(0);
352: }
356: static PetscErrorCode PCCompositeAddPC_Composite(PC pc,PCType type)
357: {
358: PC_Composite *jac;
359: PC_CompositeLink next,ilink;
360: PetscErrorCode ierr;
361: PetscInt cnt = 0;
362: const char *prefix;
363: char newprefix[8];
366: PetscNewLog(pc,&ilink);
367: ilink->next = 0;
368: PCCreate(PetscObjectComm((PetscObject)pc),&ilink->pc);
369: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->pc);
371: jac = (PC_Composite*)pc->data;
372: next = jac->head;
373: if (!next) {
374: jac->head = ilink;
375: ilink->previous = NULL;
376: } else {
377: cnt++;
378: while (next->next) {
379: next = next->next;
380: cnt++;
381: }
382: next->next = ilink;
383: ilink->previous = next;
384: }
385: PCGetOptionsPrefix(pc,&prefix);
386: PCSetOptionsPrefix(ilink->pc,prefix);
387: sprintf(newprefix,"sub_%d_",(int)cnt);
388: PCAppendOptionsPrefix(ilink->pc,newprefix);
389: /* type is set after prefix, because some methods may modify prefix, e.g. pcksp */
390: PCSetType(ilink->pc,type);
391: return(0);
392: }
396: static PetscErrorCode PCCompositeGetNumberPC_Composite(PC pc,PetscInt *n)
397: {
398: PC_Composite *jac;
399: PC_CompositeLink next;
402: jac = (PC_Composite*)pc->data;
403: next = jac->head;
404: *n = 0;
405: while (next) {
406: next = next->next;
407: (*n) ++;
408: }
409: return(0);
410: }
414: static PetscErrorCode PCCompositeGetPC_Composite(PC pc,PetscInt n,PC *subpc)
415: {
416: PC_Composite *jac;
417: PC_CompositeLink next;
418: PetscInt i;
421: jac = (PC_Composite*)pc->data;
422: next = jac->head;
423: for (i=0; i<n; i++) {
424: if (!next->next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Not enough PCs in composite preconditioner");
425: next = next->next;
426: }
427: *subpc = next->pc;
428: return(0);
429: }
431: /* -------------------------------------------------------------------------------- */
434: /*@
435: PCCompositeSetType - Sets the type of composite preconditioner.
437: Logically Collective on PC
439: Input Parameters:
440: + pc - the preconditioner context
441: - type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL
443: Options Database Key:
444: . -pc_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type
446: Level: Developer
448: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
449: @*/
450: PetscErrorCode PCCompositeSetType(PC pc,PCCompositeType type)
451: {
457: PetscTryMethod(pc,"PCCompositeSetType_C",(PC,PCCompositeType),(pc,type));
458: return(0);
459: }
463: /*@
464: PCCompositeGetType - Gets the type of composite preconditioner.
466: Logically Collective on PC
468: Input Parameter:
469: . pc - the preconditioner context
471: Output Parameter:
472: . type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL
474: Options Database Key:
475: . -pc_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type
477: Level: Developer
479: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
480: @*/
481: PetscErrorCode PCCompositeGetType(PC pc,PCCompositeType *type)
482: {
487: PetscUseMethod(pc,"PCCompositeGetType_C",(PC,PCCompositeType*),(pc,type));
488: return(0);
489: }
493: /*@
494: PCCompositeSpecialSetAlpha - Sets alpha for the special composite preconditioner
495: for alphaI + R + S
497: Logically Collective on PC
499: Input Parameter:
500: + pc - the preconditioner context
501: - alpha - scale on identity
503: Level: Developer
505: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
506: @*/
507: PetscErrorCode PCCompositeSpecialSetAlpha(PC pc,PetscScalar alpha)
508: {
514: PetscTryMethod(pc,"PCCompositeSpecialSetAlpha_C",(PC,PetscScalar),(pc,alpha));
515: return(0);
516: }
520: /*@C
521: PCCompositeAddPC - Adds another PC to the composite PC.
523: Collective on PC
525: Input Parameters:
526: + pc - the preconditioner context
527: - type - the type of the new preconditioner
529: Level: Developer
531: .keywords: PC, composite preconditioner, add
532: @*/
533: PetscErrorCode PCCompositeAddPC(PC pc,PCType type)
534: {
539: PetscTryMethod(pc,"PCCompositeAddPC_C",(PC,PCType),(pc,type));
540: return(0);
541: }
545: /*@
546: PCCompositeGetNumberPC - Gets the number of PC objects in the composite PC.
548: Not Collective
550: Input Parameter:
551: . pc - the preconditioner context
553: Output Parameter:
554: . num - the number of sub pcs
556: Level: Developer
558: .keywords: PC, get, composite preconditioner, sub preconditioner
560: .seealso: PCCompositeGetPC()
561: @*/
562: PetscErrorCode PCCompositeGetNumberPC(PC pc,PetscInt *num)
563: {
569: PetscUseMethod(pc,"PCCompositeGetNumberPC_C",(PC,PetscInt*),(pc,num));
570: return(0);
571: }
575: /*@
576: PCCompositeGetPC - Gets one of the PC objects in the composite PC.
578: Not Collective
580: Input Parameter:
581: + pc - the preconditioner context
582: - n - the number of the pc requested
584: Output Parameters:
585: . subpc - the PC requested
587: Level: Developer
589: .keywords: PC, get, composite preconditioner, sub preconditioner
591: .seealso: PCCompositeAddPC(), PCCompositeGetNumberPC()
592: @*/
593: PetscErrorCode PCCompositeGetPC(PC pc,PetscInt n,PC *subpc)
594: {
600: PetscUseMethod(pc,"PCCompositeGetPC_C",(PC,PetscInt,PC*),(pc,n,subpc));
601: return(0);
602: }
604: /* -------------------------------------------------------------------------------------------*/
606: /*MC
607: PCCOMPOSITE - Build a preconditioner by composing together several preconditioners
609: Options Database Keys:
610: + -pc_composite_type <type: one of multiplicative, additive, symmetric_multiplicative, special> - Sets composite preconditioner type
611: . -pc_use_amat - Activates PCSetUseAmat()
612: - -pc_composite_pcs - <pc0,pc1,...> list of PCs to compose
614: Level: intermediate
616: Concepts: composing solvers
618: Notes: To use a Krylov method inside the composite preconditioner, set the PCType of one or more
619: inner PCs to be PCKSP.
620: Using a Krylov method inside another Krylov method can be dangerous (you get divergence or
621: the incorrect answer) unless you use KSPFGMRES as the outer Krylov method
624: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
625: PCSHELL, PCKSP, PCCompositeSetType(), PCCompositeSpecialSetAlpha(), PCCompositeAddPC(),
626: PCCompositeGetPC(), PCSetUseAmat()
628: M*/
632: PETSC_EXTERN PetscErrorCode PCCreate_Composite(PC pc)
633: {
635: PC_Composite *jac;
638: PetscNewLog(pc,&jac);
640: pc->ops->apply = PCApply_Composite_Additive;
641: pc->ops->applytranspose = PCApplyTranspose_Composite_Additive;
642: pc->ops->setup = PCSetUp_Composite;
643: pc->ops->reset = PCReset_Composite;
644: pc->ops->destroy = PCDestroy_Composite;
645: pc->ops->setfromoptions = PCSetFromOptions_Composite;
646: pc->ops->view = PCView_Composite;
647: pc->ops->applyrichardson = 0;
649: pc->data = (void*)jac;
650: jac->type = PC_COMPOSITE_ADDITIVE;
651: jac->work1 = 0;
652: jac->work2 = 0;
653: jac->head = 0;
655: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeSetType_C",PCCompositeSetType_Composite);
656: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetType_C",PCCompositeGetType_Composite);
657: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeAddPC_C",PCCompositeAddPC_Composite);
658: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetNumberPC_C",PCCompositeGetNumberPC_Composite);
659: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetPC_C",PCCompositeGetPC_Composite);
660: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeSpecialSetAlpha_C",PCCompositeSpecialSetAlpha_Composite);
661: return(0);
662: }