Actual source code: composite.c
petsc-3.11.4 2019-09-28
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 PCCompositeAddPC() 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 PCCompositeAddPC() 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 PCCompositeAddPC() 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 PCCompositeAddPC() 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 PCCompositeAddPC() 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,0);
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","PCCompositeAddPC",pcs,&nmax,&flg);
243: if (flg) {
244: for (i=0; i<nmax; i++) {
245: PCCompositeAddPC(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,"Unkown 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,PCType type)
327: {
328: PC_Composite *jac;
329: PC_CompositeLink next,ilink;
330: PetscErrorCode ierr;
331: PetscInt cnt = 0;
332: const char *prefix;
333: char newprefix[20];
336: PetscNewLog(pc,&ilink);
337: ilink->next = 0;
338: PCCreate(PetscObjectComm((PetscObject)pc),&ilink->pc);
339: PetscObjectIncrementTabLevel((PetscObject)ilink->pc,(PetscObject)pc,1);
340: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->pc);
342: jac = (PC_Composite*)pc->data;
343: next = jac->head;
344: if (!next) {
345: jac->head = ilink;
346: ilink->previous = NULL;
347: } else {
348: cnt++;
349: while (next->next) {
350: next = next->next;
351: cnt++;
352: }
353: next->next = ilink;
354: ilink->previous = next;
355: }
356: PCGetOptionsPrefix(pc,&prefix);
357: PCSetOptionsPrefix(ilink->pc,prefix);
358: sprintf(newprefix,"sub_%d_",(int)cnt);
359: PCAppendOptionsPrefix(ilink->pc,newprefix);
360: /* type is set after prefix, because some methods may modify prefix, e.g. pcksp */
361: PCSetType(ilink->pc,type);
362: return(0);
363: }
365: static PetscErrorCode PCCompositeGetNumberPC_Composite(PC pc,PetscInt *n)
366: {
367: PC_Composite *jac;
368: PC_CompositeLink next;
371: jac = (PC_Composite*)pc->data;
372: next = jac->head;
373: *n = 0;
374: while (next) {
375: next = next->next;
376: (*n) ++;
377: }
378: return(0);
379: }
381: static PetscErrorCode PCCompositeGetPC_Composite(PC pc,PetscInt n,PC *subpc)
382: {
383: PC_Composite *jac;
384: PC_CompositeLink next;
385: PetscInt i;
388: jac = (PC_Composite*)pc->data;
389: next = jac->head;
390: for (i=0; i<n; i++) {
391: if (!next->next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Not enough PCs in composite preconditioner");
392: next = next->next;
393: }
394: *subpc = next->pc;
395: return(0);
396: }
398: /* -------------------------------------------------------------------------------- */
399: /*@
400: PCCompositeSetType - Sets the type of composite preconditioner.
402: Logically Collective on PC
404: Input Parameters:
405: + pc - the preconditioner context
406: - type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL
408: Options Database Key:
409: . -pc_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type
411: Level: Developer
413: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
414: @*/
415: PetscErrorCode PCCompositeSetType(PC pc,PCCompositeType type)
416: {
422: PetscTryMethod(pc,"PCCompositeSetType_C",(PC,PCCompositeType),(pc,type));
423: return(0);
424: }
426: /*@
427: PCCompositeGetType - Gets the type of composite preconditioner.
429: Logically Collective on PC
431: Input Parameter:
432: . pc - the preconditioner context
434: Output Parameter:
435: . type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL
437: Options Database Key:
438: . -pc_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type
440: Level: Developer
442: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
443: @*/
444: PetscErrorCode PCCompositeGetType(PC pc,PCCompositeType *type)
445: {
450: PetscUseMethod(pc,"PCCompositeGetType_C",(PC,PCCompositeType*),(pc,type));
451: return(0);
452: }
454: /*@
455: PCCompositeSpecialSetAlpha - Sets alpha for the special composite preconditioner
456: for alphaI + R + S
458: Logically Collective on PC
460: Input Parameter:
461: + pc - the preconditioner context
462: - alpha - scale on identity
464: Level: Developer
466: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
467: @*/
468: PetscErrorCode PCCompositeSpecialSetAlpha(PC pc,PetscScalar alpha)
469: {
475: PetscTryMethod(pc,"PCCompositeSpecialSetAlpha_C",(PC,PetscScalar),(pc,alpha));
476: return(0);
477: }
479: /*@C
480: PCCompositeAddPC - Adds another PC to the composite PC.
482: Collective on PC
484: Input Parameters:
485: + pc - the preconditioner context
486: - type - the type of the new preconditioner
488: Level: Developer
490: .keywords: PC, composite preconditioner, add
491: @*/
492: PetscErrorCode PCCompositeAddPC(PC pc,PCType type)
493: {
498: PetscTryMethod(pc,"PCCompositeAddPC_C",(PC,PCType),(pc,type));
499: return(0);
500: }
502: /*@
503: PCCompositeGetNumberPC - Gets the number of PC objects in the composite PC.
505: Not Collective
507: Input Parameter:
508: . pc - the preconditioner context
510: Output Parameter:
511: . num - the number of sub pcs
513: Level: Developer
515: .keywords: PC, get, composite preconditioner, sub preconditioner
517: .seealso: PCCompositeGetPC()
518: @*/
519: PetscErrorCode PCCompositeGetNumberPC(PC pc,PetscInt *num)
520: {
526: PetscUseMethod(pc,"PCCompositeGetNumberPC_C",(PC,PetscInt*),(pc,num));
527: return(0);
528: }
530: /*@
531: PCCompositeGetPC - Gets one of the PC objects in the composite PC.
533: Not Collective
535: Input Parameter:
536: + pc - the preconditioner context
537: - n - the number of the pc requested
539: Output Parameters:
540: . subpc - the PC requested
542: Level: Developer
544: Notes:
545: To use a different operator to construct one of the inner preconditioners first call PCCompositeGetPC(), then
546: call PCSetOperators() on that PC.
548: .keywords: PC, get, composite preconditioner, sub preconditioner
550: .seealso: PCCompositeAddPC(), PCCompositeGetNumberPC(), PCSetOperators()
551: @*/
552: PetscErrorCode PCCompositeGetPC(PC pc,PetscInt n,PC *subpc)
553: {
559: PetscUseMethod(pc,"PCCompositeGetPC_C",(PC,PetscInt,PC*),(pc,n,subpc));
560: return(0);
561: }
563: /* -------------------------------------------------------------------------------------------*/
565: /*MC
566: PCCOMPOSITE - Build a preconditioner by composing together several preconditioners
568: Options Database Keys:
569: + -pc_composite_type <type: one of multiplicative, additive, symmetric_multiplicative, special> - Sets composite preconditioner type
570: . -pc_use_amat - activates PCSetUseAmat()
571: - -pc_composite_pcs - <pc0,pc1,...> list of PCs to compose
573: Level: intermediate
575: Concepts: composing solvers
577: Notes:
578: To use a Krylov method inside the composite preconditioner, set the PCType of one or more
579: inner PCs to be PCKSP.
580: Using a Krylov method inside another Krylov method can be dangerous (you get divergence or
581: the incorrect answer) unless you use KSPFGMRES as the outer Krylov method
582: To use a different operator to construct one of the inner preconditioners first call PCCompositeGetPC(), then
583: call PCSetOperators() on that PC.
586: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
587: PCSHELL, PCKSP, PCCompositeSetType(), PCCompositeSpecialSetAlpha(), PCCompositeAddPC(),
588: PCCompositeGetPC(), PCSetUseAmat()
590: M*/
592: PETSC_EXTERN PetscErrorCode PCCreate_Composite(PC pc)
593: {
595: PC_Composite *jac;
598: PetscNewLog(pc,&jac);
600: pc->ops->apply = PCApply_Composite_Additive;
601: pc->ops->applytranspose = PCApplyTranspose_Composite_Additive;
602: pc->ops->setup = PCSetUp_Composite;
603: pc->ops->reset = PCReset_Composite;
604: pc->ops->destroy = PCDestroy_Composite;
605: pc->ops->setfromoptions = PCSetFromOptions_Composite;
606: pc->ops->view = PCView_Composite;
607: pc->ops->applyrichardson = 0;
609: pc->data = (void*)jac;
610: jac->type = PC_COMPOSITE_ADDITIVE;
611: jac->work1 = 0;
612: jac->work2 = 0;
613: jac->head = 0;
615: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeSetType_C",PCCompositeSetType_Composite);
616: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetType_C",PCCompositeGetType_Composite);
617: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeAddPC_C",PCCompositeAddPC_Composite);
618: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetNumberPC_C",PCCompositeGetNumberPC_Composite);
619: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetPC_C",PCCompositeGetPC_Composite);
620: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeSpecialSetAlpha_C",PCCompositeSpecialSetAlpha_Composite);
621: return(0);
622: }