Actual source code: composite.c
petsc-3.8.4 2018-03-24
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: VecSet(jac->work1,0.0); /* zero since some PC's may not set all entries in the result */
51: PCApply(next->pc,jac->work2,jac->work1); /* work1 <- C work2 */
52: VecAXPY(y,1.0,jac->work1); /* y <- y + work1 = B x + C (x - A B x) = (B + C (1 - A B)) x */
53: }
54: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
55: while (next->previous) {
56: next = next->previous;
57: MatMult(mat,y,jac->work1);
58: VecWAXPY(jac->work2,-1.0,jac->work1,x);
59: VecSet(jac->work1,0.0); /* zero since some PC's may not set all entries in the result */
60: PCApply(next->pc,jac->work2,jac->work1);
61: VecAXPY(y,1.0,jac->work1);
62: }
63: }
64: return(0);
65: }
67: static PetscErrorCode PCApplyTranspose_Composite_Multiplicative(PC pc,Vec x,Vec y)
68: {
69: PetscErrorCode ierr;
70: PC_Composite *jac = (PC_Composite*)pc->data;
71: PC_CompositeLink next = jac->head;
72: Mat mat = pc->pmat;
75: if (!next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC() or -pc_composite_pcs");
76: if (next->next && !jac->work2) { /* allocate second work vector */
77: VecDuplicate(jac->work1,&jac->work2);
78: }
79: if (pc->useAmat) mat = pc->mat;
80: /* locate last PC */
81: while (next->next) {
82: next = next->next;
83: }
84: PCApplyTranspose(next->pc,x,y);
85: while (next->previous) {
86: next = next->previous;
87: MatMultTranspose(mat,y,jac->work1);
88: VecWAXPY(jac->work2,-1.0,jac->work1,x);
89: VecSet(jac->work1,0.0); /* zero since some PC's may not set all entries in the result */
90: PCApplyTranspose(next->pc,jac->work2,jac->work1);
91: VecAXPY(y,1.0,jac->work1);
92: }
93: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
94: next = jac->head;
95: while (next->next) {
96: next = next->next;
97: MatMultTranspose(mat,y,jac->work1);
98: VecWAXPY(jac->work2,-1.0,jac->work1,x);
99: VecSet(jac->work1,0.0); /* zero since some PC's may not set all entries in the result */
100: PCApplyTranspose(next->pc,jac->work2,jac->work1);
101: VecAXPY(y,1.0,jac->work1);
102: }
103: }
104: return(0);
105: }
107: /*
108: This is very special for a matrix of the form alpha I + R + S
109: where first preconditioner is built from alpha I + S and second from
110: alpha I + R
111: */
112: static PetscErrorCode PCApply_Composite_Special(PC pc,Vec x,Vec y)
113: {
114: PetscErrorCode ierr;
115: PC_Composite *jac = (PC_Composite*)pc->data;
116: PC_CompositeLink next = jac->head;
119: if (!next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC() or -pc_composite_pcs");
120: if (!next->next || next->next->next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Special composite preconditioners requires exactly two PCs");
122: /* Set the reuse flag on children PCs */
123: PCSetReusePreconditioner(next->pc,pc->reusepreconditioner);
124: PCSetReusePreconditioner(next->next->pc,pc->reusepreconditioner);
126: PCApply(next->pc,x,jac->work1);
127: PCApply(next->next->pc,jac->work1,y);
128: return(0);
129: }
131: static PetscErrorCode PCApply_Composite_Additive(PC pc,Vec x,Vec y)
132: {
133: PetscErrorCode ierr;
134: PC_Composite *jac = (PC_Composite*)pc->data;
135: PC_CompositeLink next = jac->head;
138: if (!next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC() or -pc_composite_pcs");
140: /* Set the reuse flag on children PCs */
141: while (next) {
142: PCSetReusePreconditioner(next->pc,pc->reusepreconditioner);
143: next = next->next;
144: }
145: next = jac->head;
147: PCApply(next->pc,x,y);
148: while (next->next) {
149: next = next->next;
150: VecSet(jac->work1,0.0); /* zero since some PC's may not set all entries in the result */
151: PCApply(next->pc,x,jac->work1);
152: VecAXPY(y,1.0,jac->work1);
153: }
154: return(0);
155: }
157: static PetscErrorCode PCApplyTranspose_Composite_Additive(PC pc,Vec x,Vec y)
158: {
159: PetscErrorCode ierr;
160: PC_Composite *jac = (PC_Composite*)pc->data;
161: PC_CompositeLink next = jac->head;
164: if (!next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC() or -pc_composite_pcs");
165: PCApplyTranspose(next->pc,x,y);
166: while (next->next) {
167: next = next->next;
168: VecSet(jac->work1,0.0); /* zero since some PC's may not set all entries in the result */
169: PCApplyTranspose(next->pc,x,jac->work1);
170: VecAXPY(y,1.0,jac->work1);
171: }
172: return(0);
173: }
175: static PetscErrorCode PCSetUp_Composite(PC pc)
176: {
177: PetscErrorCode ierr;
178: PC_Composite *jac = (PC_Composite*)pc->data;
179: PC_CompositeLink next = jac->head;
180: DM dm;
183: if (!jac->work1) {
184: MatCreateVecs(pc->pmat,&jac->work1,0);
185: }
186: PCGetDM(pc,&dm);
187: while (next) {
188: PCSetDM(next->pc,dm);
189: PCSetOperators(next->pc,pc->mat,pc->pmat);
190: next = next->next;
191: }
192: return(0);
193: }
195: static PetscErrorCode PCReset_Composite(PC pc)
196: {
197: PC_Composite *jac = (PC_Composite*)pc->data;
198: PetscErrorCode ierr;
199: PC_CompositeLink next = jac->head;
202: while (next) {
203: PCReset(next->pc);
204: next = next->next;
205: }
206: VecDestroy(&jac->work1);
207: VecDestroy(&jac->work2);
208: return(0);
209: }
211: static PetscErrorCode PCDestroy_Composite(PC pc)
212: {
213: PC_Composite *jac = (PC_Composite*)pc->data;
214: PetscErrorCode ierr;
215: PC_CompositeLink next = jac->head,next_tmp;
218: PCReset_Composite(pc);
219: while (next) {
220: PCDestroy(&next->pc);
221: next_tmp = next;
222: next = next->next;
223: PetscFree(next_tmp);
224: }
225: PetscFree(pc->data);
226: return(0);
227: }
229: static PetscErrorCode PCSetFromOptions_Composite(PetscOptionItems *PetscOptionsObject,PC pc)
230: {
231: PC_Composite *jac = (PC_Composite*)pc->data;
232: PetscErrorCode ierr;
233: PetscInt nmax = 8,i;
234: PC_CompositeLink next;
235: char *pcs[8];
236: PetscBool flg;
239: PetscOptionsHead(PetscOptionsObject,"Composite preconditioner options");
240: PetscOptionsEnum("-pc_composite_type","Type of composition","PCCompositeSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&jac->type,&flg);
241: if (flg) {
242: PCCompositeSetType(pc,jac->type);
243: }
244: PetscOptionsStringArray("-pc_composite_pcs","List of composite solvers","PCCompositeAddPC",pcs,&nmax,&flg);
245: if (flg) {
246: for (i=0; i<nmax; i++) {
247: PCCompositeAddPC(pc,pcs[i]);
248: PetscFree(pcs[i]); /* deallocate string pcs[i], which is allocated in PetscOptionsStringArray() */
249: }
250: }
251: PetscOptionsTail();
253: next = jac->head;
254: while (next) {
255: PCSetFromOptions(next->pc);
256: next = next->next;
257: }
258: return(0);
259: }
261: static PetscErrorCode PCView_Composite(PC pc,PetscViewer viewer)
262: {
263: PC_Composite *jac = (PC_Composite*)pc->data;
264: PetscErrorCode ierr;
265: PC_CompositeLink next = jac->head;
266: PetscBool iascii;
269: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
270: if (iascii) {
271: PetscViewerASCIIPrintf(viewer,"Composite PC type - %s\n",PCCompositeTypes[jac->type]);
272: PetscViewerASCIIPrintf(viewer,"PCs on composite preconditioner follow\n");
273: PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
274: }
275: if (iascii) {
276: PetscViewerASCIIPushTab(viewer);
277: }
278: while (next) {
279: PCView(next->pc,viewer);
280: next = next->next;
281: }
282: if (iascii) {
283: PetscViewerASCIIPopTab(viewer);
284: PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
285: }
286: return(0);
287: }
289: /* ------------------------------------------------------------------------------*/
291: static PetscErrorCode PCCompositeSpecialSetAlpha_Composite(PC pc,PetscScalar alpha)
292: {
293: PC_Composite *jac = (PC_Composite*)pc->data;
296: jac->alpha = alpha;
297: return(0);
298: }
300: static PetscErrorCode PCCompositeSetType_Composite(PC pc,PCCompositeType type)
301: {
302: PC_Composite *jac = (PC_Composite*)pc->data;
305: if (type == PC_COMPOSITE_ADDITIVE) {
306: pc->ops->apply = PCApply_Composite_Additive;
307: pc->ops->applytranspose = PCApplyTranspose_Composite_Additive;
308: } else if (type == PC_COMPOSITE_MULTIPLICATIVE || type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
309: pc->ops->apply = PCApply_Composite_Multiplicative;
310: pc->ops->applytranspose = PCApplyTranspose_Composite_Multiplicative;
311: } else if (type == PC_COMPOSITE_SPECIAL) {
312: pc->ops->apply = PCApply_Composite_Special;
313: pc->ops->applytranspose = NULL;
314: } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Unkown composite preconditioner type");
315: jac->type = type;
316: return(0);
317: }
319: static PetscErrorCode PCCompositeGetType_Composite(PC pc,PCCompositeType *type)
320: {
321: PC_Composite *jac = (PC_Composite*)pc->data;
324: *type = jac->type;
325: return(0);
326: }
328: static PetscErrorCode PCCompositeAddPC_Composite(PC pc,PCType type)
329: {
330: PC_Composite *jac;
331: PC_CompositeLink next,ilink;
332: PetscErrorCode ierr;
333: PetscInt cnt = 0;
334: const char *prefix;
335: char newprefix[20];
338: PetscNewLog(pc,&ilink);
339: ilink->next = 0;
340: PCCreate(PetscObjectComm((PetscObject)pc),&ilink->pc);
341: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->pc);
343: jac = (PC_Composite*)pc->data;
344: next = jac->head;
345: if (!next) {
346: jac->head = ilink;
347: ilink->previous = NULL;
348: } else {
349: cnt++;
350: while (next->next) {
351: next = next->next;
352: cnt++;
353: }
354: next->next = ilink;
355: ilink->previous = next;
356: }
357: PCGetOptionsPrefix(pc,&prefix);
358: PCSetOptionsPrefix(ilink->pc,prefix);
359: sprintf(newprefix,"sub_%d_",(int)cnt);
360: PCAppendOptionsPrefix(ilink->pc,newprefix);
361: /* type is set after prefix, because some methods may modify prefix, e.g. pcksp */
362: PCSetType(ilink->pc,type);
363: return(0);
364: }
366: static PetscErrorCode PCCompositeGetNumberPC_Composite(PC pc,PetscInt *n)
367: {
368: PC_Composite *jac;
369: PC_CompositeLink next;
372: jac = (PC_Composite*)pc->data;
373: next = jac->head;
374: *n = 0;
375: while (next) {
376: next = next->next;
377: (*n) ++;
378: }
379: return(0);
380: }
382: static PetscErrorCode PCCompositeGetPC_Composite(PC pc,PetscInt n,PC *subpc)
383: {
384: PC_Composite *jac;
385: PC_CompositeLink next;
386: PetscInt i;
389: jac = (PC_Composite*)pc->data;
390: next = jac->head;
391: for (i=0; i<n; i++) {
392: if (!next->next) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Not enough PCs in composite preconditioner");
393: next = next->next;
394: }
395: *subpc = next->pc;
396: return(0);
397: }
399: /* -------------------------------------------------------------------------------- */
400: /*@
401: PCCompositeSetType - Sets the type of composite preconditioner.
403: Logically Collective on PC
405: Input Parameters:
406: + pc - the preconditioner context
407: - type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL
409: Options Database Key:
410: . -pc_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type
412: Level: Developer
414: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
415: @*/
416: PetscErrorCode PCCompositeSetType(PC pc,PCCompositeType type)
417: {
423: PetscTryMethod(pc,"PCCompositeSetType_C",(PC,PCCompositeType),(pc,type));
424: return(0);
425: }
427: /*@
428: PCCompositeGetType - Gets the type of composite preconditioner.
430: Logically Collective on PC
432: Input Parameter:
433: . pc - the preconditioner context
435: Output Parameter:
436: . type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL
438: Options Database Key:
439: . -pc_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type
441: Level: Developer
443: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
444: @*/
445: PetscErrorCode PCCompositeGetType(PC pc,PCCompositeType *type)
446: {
451: PetscUseMethod(pc,"PCCompositeGetType_C",(PC,PCCompositeType*),(pc,type));
452: return(0);
453: }
455: /*@
456: PCCompositeSpecialSetAlpha - Sets alpha for the special composite preconditioner
457: for alphaI + R + S
459: Logically Collective on PC
461: Input Parameter:
462: + pc - the preconditioner context
463: - alpha - scale on identity
465: Level: Developer
467: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
468: @*/
469: PetscErrorCode PCCompositeSpecialSetAlpha(PC pc,PetscScalar alpha)
470: {
476: PetscTryMethod(pc,"PCCompositeSpecialSetAlpha_C",(PC,PetscScalar),(pc,alpha));
477: return(0);
478: }
480: /*@C
481: PCCompositeAddPC - Adds another PC to the composite PC.
483: Collective on PC
485: Input Parameters:
486: + pc - the preconditioner context
487: - type - the type of the new preconditioner
489: Level: Developer
491: .keywords: PC, composite preconditioner, add
492: @*/
493: PetscErrorCode PCCompositeAddPC(PC pc,PCType type)
494: {
499: PetscTryMethod(pc,"PCCompositeAddPC_C",(PC,PCType),(pc,type));
500: return(0);
501: }
503: /*@
504: PCCompositeGetNumberPC - Gets the number of PC objects in the composite PC.
506: Not Collective
508: Input Parameter:
509: . pc - the preconditioner context
511: Output Parameter:
512: . num - the number of sub pcs
514: Level: Developer
516: .keywords: PC, get, composite preconditioner, sub preconditioner
518: .seealso: PCCompositeGetPC()
519: @*/
520: PetscErrorCode PCCompositeGetNumberPC(PC pc,PetscInt *num)
521: {
527: PetscUseMethod(pc,"PCCompositeGetNumberPC_C",(PC,PetscInt*),(pc,num));
528: return(0);
529: }
531: /*@
532: PCCompositeGetPC - Gets one of the PC objects in the composite PC.
534: Not Collective
536: Input Parameter:
537: + pc - the preconditioner context
538: - n - the number of the pc requested
540: Output Parameters:
541: . subpc - the PC requested
543: Level: Developer
545: .keywords: PC, get, composite preconditioner, sub preconditioner
547: .seealso: PCCompositeAddPC(), PCCompositeGetNumberPC()
548: @*/
549: PetscErrorCode PCCompositeGetPC(PC pc,PetscInt n,PC *subpc)
550: {
556: PetscUseMethod(pc,"PCCompositeGetPC_C",(PC,PetscInt,PC*),(pc,n,subpc));
557: return(0);
558: }
560: /* -------------------------------------------------------------------------------------------*/
562: /*MC
563: PCCOMPOSITE - Build a preconditioner by composing together several preconditioners
565: Options Database Keys:
566: + -pc_composite_type <type: one of multiplicative, additive, symmetric_multiplicative, special> - Sets composite preconditioner type
567: . -pc_use_amat - Activates PCSetUseAmat()
568: - -pc_composite_pcs - <pc0,pc1,...> list of PCs to compose
570: Level: intermediate
572: Concepts: composing solvers
574: Notes: To use a Krylov method inside the composite preconditioner, set the PCType of one or more
575: inner PCs to be PCKSP.
576: Using a Krylov method inside another Krylov method can be dangerous (you get divergence or
577: the incorrect answer) unless you use KSPFGMRES as the outer Krylov method
580: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
581: PCSHELL, PCKSP, PCCompositeSetType(), PCCompositeSpecialSetAlpha(), PCCompositeAddPC(),
582: PCCompositeGetPC(), PCSetUseAmat()
584: M*/
586: PETSC_EXTERN PetscErrorCode PCCreate_Composite(PC pc)
587: {
589: PC_Composite *jac;
592: PetscNewLog(pc,&jac);
594: pc->ops->apply = PCApply_Composite_Additive;
595: pc->ops->applytranspose = PCApplyTranspose_Composite_Additive;
596: pc->ops->setup = PCSetUp_Composite;
597: pc->ops->reset = PCReset_Composite;
598: pc->ops->destroy = PCDestroy_Composite;
599: pc->ops->setfromoptions = PCSetFromOptions_Composite;
600: pc->ops->view = PCView_Composite;
601: pc->ops->applyrichardson = 0;
603: pc->data = (void*)jac;
604: jac->type = PC_COMPOSITE_ADDITIVE;
605: jac->work1 = 0;
606: jac->work2 = 0;
607: jac->head = 0;
609: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeSetType_C",PCCompositeSetType_Composite);
610: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetType_C",PCCompositeGetType_Composite);
611: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeAddPC_C",PCCompositeAddPC_Composite);
612: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetNumberPC_C",PCCompositeGetNumberPC_Composite);
613: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeGetPC_C",PCCompositeGetPC_Composite);
614: PetscObjectComposeFunction((PetscObject)pc,"PCCompositeSpecialSetAlpha_C",PCCompositeSpecialSetAlpha_Composite);
615: return(0);
616: }