Actual source code: composite.c
petsc-3.13.6 2020-09-29
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,"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,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: @*/
414: PetscErrorCode PCCompositeSetType(PC pc,PCCompositeType type)
415: {
421: PetscTryMethod(pc,"PCCompositeSetType_C",(PC,PCCompositeType),(pc,type));
422: return(0);
423: }
425: /*@
426: PCCompositeGetType - Gets the type of composite preconditioner.
428: Logically Collective on PC
430: Input Parameter:
431: . pc - the preconditioner context
433: Output Parameter:
434: . type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL
436: Options Database Key:
437: . -pc_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type
439: Level: Developer
441: @*/
442: PetscErrorCode PCCompositeGetType(PC pc,PCCompositeType *type)
443: {
448: PetscUseMethod(pc,"PCCompositeGetType_C",(PC,PCCompositeType*),(pc,type));
449: return(0);
450: }
452: /*@
453: PCCompositeSpecialSetAlpha - Sets alpha for the special composite preconditioner
454: for alphaI + R + S
456: Logically Collective on PC
458: Input Parameter:
459: + pc - the preconditioner context
460: - alpha - scale on identity
462: Level: Developer
464: @*/
465: PetscErrorCode PCCompositeSpecialSetAlpha(PC pc,PetscScalar alpha)
466: {
472: PetscTryMethod(pc,"PCCompositeSpecialSetAlpha_C",(PC,PetscScalar),(pc,alpha));
473: return(0);
474: }
476: /*@C
477: PCCompositeAddPC - Adds another PC to the composite PC.
479: Collective on PC
481: Input Parameters:
482: + pc - the preconditioner context
483: - type - the type of the new preconditioner
485: Level: Developer
487: @*/
488: PetscErrorCode PCCompositeAddPC(PC pc,PCType type)
489: {
494: PetscTryMethod(pc,"PCCompositeAddPC_C",(PC,PCType),(pc,type));
495: return(0);
496: }
498: /*@
499: PCCompositeGetNumberPC - Gets the number of PC objects in the composite PC.
501: Not Collective
503: Input Parameter:
504: . pc - the preconditioner context
506: Output Parameter:
507: . num - the number of sub pcs
509: Level: Developer
511: .seealso: PCCompositeGetPC()
512: @*/
513: PetscErrorCode PCCompositeGetNumberPC(PC pc,PetscInt *num)
514: {
520: PetscUseMethod(pc,"PCCompositeGetNumberPC_C",(PC,PetscInt*),(pc,num));
521: return(0);
522: }
524: /*@
525: PCCompositeGetPC - Gets one of the PC objects in the composite PC.
527: Not Collective
529: Input Parameter:
530: + pc - the preconditioner context
531: - n - the number of the pc requested
533: Output Parameters:
534: . subpc - the PC requested
536: Level: Developer
538: Notes:
539: To use a different operator to construct one of the inner preconditioners first call PCCompositeGetPC(), then
540: call PCSetOperators() on that PC.
542: .seealso: PCCompositeAddPC(), PCCompositeGetNumberPC(), PCSetOperators()
543: @*/
544: PetscErrorCode PCCompositeGetPC(PC pc,PetscInt n,PC *subpc)
545: {
551: PetscUseMethod(pc,"PCCompositeGetPC_C",(PC,PetscInt,PC*),(pc,n,subpc));
552: return(0);
553: }
555: /* -------------------------------------------------------------------------------------------*/
557: /*MC
558: PCCOMPOSITE - Build a preconditioner by composing together several preconditioners
560: Options Database Keys:
561: + -pc_composite_type <type: one of multiplicative, additive, symmetric_multiplicative, special> - Sets composite preconditioner type
562: . -pc_use_amat - activates PCSetUseAmat()
563: - -pc_composite_pcs - <pc0,pc1,...> list of PCs to compose
565: Level: intermediate
567: Notes:
568: To use a Krylov method inside the composite preconditioner, set the PCType of one or more
569: inner PCs to be PCKSP.
570: Using a Krylov method inside another Krylov method can be dangerous (you get divergence or
571: the incorrect answer) unless you use KSPFGMRES as the outer Krylov method
572: To use a different operator to construct one of the inner preconditioners first call PCCompositeGetPC(), then
573: call PCSetOperators() on that PC.
576: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
577: PCSHELL, PCKSP, PCCompositeSetType(), PCCompositeSpecialSetAlpha(), PCCompositeAddPC(),
578: PCCompositeGetPC(), PCSetUseAmat()
580: M*/
582: PETSC_EXTERN PetscErrorCode PCCreate_Composite(PC pc)
583: {
585: PC_Composite *jac;
588: PetscNewLog(pc,&jac);
590: pc->ops->apply = PCApply_Composite_Additive;
591: pc->ops->applytranspose = PCApplyTranspose_Composite_Additive;
592: pc->ops->setup = PCSetUp_Composite;
593: pc->ops->reset = PCReset_Composite;
594: pc->ops->destroy = PCDestroy_Composite;
595: pc->ops->setfromoptions = PCSetFromOptions_Composite;
596: pc->ops->view = PCView_Composite;
597: pc->ops->applyrichardson = 0;
599: pc->data = (void*)jac;
600: jac->type = PC_COMPOSITE_ADDITIVE;
601: jac->work1 = 0;
602: jac->work2 = 0;
603: jac->head = 0;
605: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeSetType_C",PCCompositeSetType_Composite);
606: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetType_C",PCCompositeGetType_Composite);
607: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeAddPC_C",PCCompositeAddPC_Composite);
608: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetNumberPC_C",PCCompositeGetNumberPC_Composite);
609: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetPC_C",PCCompositeGetPC_Composite);
610: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeSpecialSetAlpha_C",PCCompositeSpecialSetAlpha_Composite);
611: return(0);
612: }