Actual source code: composite.c
1: #define PETSCKSP_DLL
3: /*
4: Defines a preconditioner that can consist of a collection of PCs
5: */
6: #include private/pcimpl.h
7: #include petscksp.h
9: typedef struct _PC_CompositeLink *PC_CompositeLink;
10: struct _PC_CompositeLink {
11: PC pc;
12: PC_CompositeLink next;
13: PC_CompositeLink previous;
14: };
15:
16: typedef struct {
17: PC_CompositeLink head;
18: PCCompositeType type;
19: Vec work1;
20: Vec work2;
21: PetscScalar alpha;
22: PetscTruth use_true_matrix;
23: } PC_Composite;
27: static PetscErrorCode PCApply_Composite_Multiplicative(PC pc,Vec x,Vec y)
28: {
29: PetscErrorCode ierr;
30: PC_Composite *jac = (PC_Composite*)pc->data;
31: PC_CompositeLink next = jac->head;
32: Mat mat = pc->pmat;
35: if (!next) {
36: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC() or -pc_composite_pcs");
37: }
38: if (next->next && !jac->work2) { /* allocate second work vector */
39: VecDuplicate(jac->work1,&jac->work2);
40: }
41: PCApply(next->pc,x,y);
42: if (jac->use_true_matrix) mat = pc->mat;
43: while (next->next) {
44: next = next->next;
45: MatMult(mat,y,jac->work1);
46: VecWAXPY(jac->work2,-1.0,jac->work1,x);
47: VecSet(jac->work1,0.0); /* zero since some PC's may not set all entries in the result */
48: PCApply(next->pc,jac->work2,jac->work1);
49: VecAXPY(y,1.0,jac->work1);
50: }
51: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
52: while (next->previous) {
53: next = next->previous;
54: MatMult(mat,y,jac->work1);
55: VecWAXPY(jac->work2,-1.0,jac->work1,x);
56: VecSet(jac->work1,0.0); /* zero since some PC's may not set all entries in the result */
57: PCApply(next->pc,jac->work2,jac->work1);
58: VecAXPY(y,1.0,jac->work1);
59: }
60: }
61: return(0);
62: }
64: /*
65: This is very special for a matrix of the form alpha I + R + S
66: where first preconditioner is built from alpha I + S and second from
67: alpha I + R
68: */
71: static PetscErrorCode PCApply_Composite_Special(PC pc,Vec x,Vec y)
72: {
73: PetscErrorCode ierr;
74: PC_Composite *jac = (PC_Composite*)pc->data;
75: PC_CompositeLink next = jac->head;
78: if (!next) {
79: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC() or -pc_composite_pcs");
80: }
81: if (!next->next || next->next->next) {
82: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Special composite preconditioners requires exactly two PCs");
83: }
85: PCApply(next->pc,x,jac->work1);
86: PCApply(next->next->pc,jac->work1,y);
87: return(0);
88: }
92: static PetscErrorCode PCApply_Composite_Additive(PC pc,Vec x,Vec y)
93: {
94: PetscErrorCode ierr;
95: PC_Composite *jac = (PC_Composite*)pc->data;
96: PC_CompositeLink next = jac->head;
99: if (!next) {
100: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC() or -pc_composite_pcs");
101: }
102: PCApply(next->pc,x,y);
103: while (next->next) {
104: next = next->next;
105: VecSet(jac->work1,0.0); /* zero since some PC's may not set all entries in the result */
106: PCApply(next->pc,x,jac->work1);
107: VecAXPY(y,1.0,jac->work1);
108: }
109: return(0);
110: }
114: static PetscErrorCode PCSetUp_Composite(PC pc)
115: {
116: PetscErrorCode ierr;
117: PC_Composite *jac = (PC_Composite*)pc->data;
118: PC_CompositeLink next = jac->head;
121: if (!jac->work1) {
122: MatGetVecs(pc->pmat,&jac->work1,0);
123: }
124: while (next) {
125: PCSetOperators(next->pc,pc->mat,pc->pmat,pc->flag);
126: next = next->next;
127: }
128: return(0);
129: }
133: static PetscErrorCode PCDestroy_Composite(PC pc)
134: {
135: PC_Composite *jac = (PC_Composite*)pc->data;
136: PetscErrorCode ierr;
137: PC_CompositeLink next = jac->head;
140: while (next) {
141: PCDestroy(next->pc);
142: next = next->next;
143: }
145: if (jac->work1) {VecDestroy(jac->work1);}
146: if (jac->work2) {VecDestroy(jac->work2);}
147: PetscFree(jac);
148: return(0);
149: }
153: static PetscErrorCode PCSetFromOptions_Composite(PC pc)
154: {
155: PC_Composite *jac = (PC_Composite*)pc->data;
156: PetscErrorCode ierr;
157: PetscInt nmax = 8,i;
158: PC_CompositeLink next;
159: char *pcs[8];
160: PetscTruth flg;
163: PetscOptionsHead("Composite preconditioner options");
164: PetscOptionsEnum("-pc_composite_type","Type of composition","PCCompositeSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&jac->type,&flg);
165: if (flg) {
166: PCCompositeSetType(pc,jac->type);
167: }
168: PetscOptionsName("-pc_composite_true","Use true matrix for inner solves","PCCompositeSetUseTrue",&flg);
169: if (flg) {
170: PCCompositeSetUseTrue(pc);
171: }
172: PetscOptionsStringArray("-pc_composite_pcs","List of composite solvers","PCCompositeAddPC",pcs,&nmax,&flg);
173: if (flg) {
174: for (i=0; i<nmax; i++) {
175: PCCompositeAddPC(pc,pcs[i]);
176: }
177: }
178: PetscOptionsTail();
180: next = jac->head;
181: while (next) {
182: PCSetFromOptions(next->pc);
183: next = next->next;
184: }
185: return(0);
186: }
190: static PetscErrorCode PCView_Composite(PC pc,PetscViewer viewer)
191: {
192: PC_Composite *jac = (PC_Composite*)pc->data;
193: PetscErrorCode ierr;
194: PC_CompositeLink next = jac->head;
195: PetscTruth iascii;
198: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
199: if (iascii) {
200: PetscViewerASCIIPrintf(viewer,"Composite PC type - %s\n",PCCompositeTypes[jac->type]);
201: PetscViewerASCIIPrintf(viewer,"PCs on composite preconditioner follow\n");
202: PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
203: } else {
204: SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCComposite",((PetscObject)viewer)->type_name);
205: }
206: if (iascii) {
207: PetscViewerASCIIPushTab(viewer);
208: }
209: while (next) {
210: PCView(next->pc,viewer);
211: next = next->next;
212: }
213: if (iascii) {
214: PetscViewerASCIIPopTab(viewer);
215: PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
216: }
217: return(0);
218: }
220: /* ------------------------------------------------------------------------------*/
225: PetscErrorCode PCCompositeSpecialSetAlpha_Composite(PC pc,PetscScalar alpha)
226: {
227: PC_Composite *jac = (PC_Composite*)pc->data;
229: jac->alpha = alpha;
230: return(0);
231: }
237: PetscErrorCode PCCompositeSetType_Composite(PC pc,PCCompositeType type)
238: {
240: if (type == PC_COMPOSITE_ADDITIVE) {
241: pc->ops->apply = PCApply_Composite_Additive;
242: } else if (type == PC_COMPOSITE_MULTIPLICATIVE || type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
243: pc->ops->apply = PCApply_Composite_Multiplicative;
244: } else if (type == PC_COMPOSITE_SPECIAL) {
245: pc->ops->apply = PCApply_Composite_Special;
246: } else {
247: SETERRQ(PETSC_ERR_ARG_WRONG,"Unkown composite preconditioner type");
248: }
249: return(0);
250: }
256: PetscErrorCode PCCompositeAddPC_Composite(PC pc,PCType type)
257: {
258: PC_Composite *jac;
259: PC_CompositeLink next,ilink;
260: PetscErrorCode ierr;
261: PetscInt cnt = 0;
262: const char *prefix;
263: char newprefix[8];
266: PetscNew(struct _PC_CompositeLink,&ilink);
267: ilink->next = 0;
268: PCCreate(pc->comm,&ilink->pc);
270: jac = (PC_Composite*)pc->data;
271: next = jac->head;
272: if (!next) {
273: jac->head = ilink;
274: ilink->previous = PETSC_NULL;
275: } else {
276: cnt++;
277: while (next->next) {
278: next = next->next;
279: cnt++;
280: }
281: next->next = ilink;
282: ilink->previous = next;
283: }
284: PCGetOptionsPrefix(pc,&prefix);
285: PCSetOptionsPrefix(ilink->pc,prefix);
286: sprintf(newprefix,"sub_%d_",(int)cnt);
287: PCAppendOptionsPrefix(ilink->pc,newprefix);
288: /* type is set after prefix, because some methods may modify prefix, e.g. pcksp */
289: PCSetType(ilink->pc,type);
290: return(0);
291: }
297: PetscErrorCode PCCompositeGetPC_Composite(PC pc,PetscInt n,PC *subpc)
298: {
299: PC_Composite *jac;
300: PC_CompositeLink next;
301: PetscInt i;
304: jac = (PC_Composite*)pc->data;
305: next = jac->head;
306: for (i=0; i<n; i++) {
307: if (!next->next) {
308: SETERRQ(PETSC_ERR_ARG_INCOMP,"Not enough PCs in composite preconditioner");
309: }
310: next = next->next;
311: }
312: *subpc = next->pc;
313: return(0);
314: }
320: PetscErrorCode PCCompositeSetUseTrue_Composite(PC pc)
321: {
322: PC_Composite *jac;
325: jac = (PC_Composite*)pc->data;
326: jac->use_true_matrix = PETSC_TRUE;
327: return(0);
328: }
331: /* -------------------------------------------------------------------------------- */
334: /*@
335: PCCompositeSetType - Sets the type of composite preconditioner.
336:
337: Collective on PC
339: Input Parameter:
340: + pc - the preconditioner context
341: - type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL
343: Options Database Key:
344: . -pc_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type
346: Level: Developer
348: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
349: @*/
350: PetscErrorCode PCCompositeSetType(PC pc,PCCompositeType type)
351: {
352: PetscErrorCode ierr,(*f)(PC,PCCompositeType);
356: PetscObjectQueryFunction((PetscObject)pc,"PCCompositeSetType_C",(void (**)(void))&f);
357: if (f) {
358: (*f)(pc,type);
359: }
360: return(0);
361: }
365: /*@
366: PCCompositeSpecialSetAlpha - Sets alpha for the special composite preconditioner
367: for alphaI + R + S
368:
369: Collective on PC
371: Input Parameter:
372: + pc - the preconditioner context
373: - alpha - scale on identity
375: Level: Developer
377: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
378: @*/
379: PetscErrorCode PCCompositeSpecialSetAlpha(PC pc,PetscScalar alpha)
380: {
381: PetscErrorCode ierr,(*f)(PC,PetscScalar);
385: PetscObjectQueryFunction((PetscObject)pc,"PCCompositeSpecialSetAlpha_C",(void (**)(void))&f);
386: if (f) {
387: (*f)(pc,alpha);
388: }
389: return(0);
390: }
394: /*@C
395: PCCompositeAddPC - Adds another PC to the composite PC.
396:
397: Collective on PC
399: Input Parameters:
400: + pc - the preconditioner context
401: - type - the type of the new preconditioner
403: Level: Developer
405: .keywords: PC, composite preconditioner, add
406: @*/
407: PetscErrorCode PCCompositeAddPC(PC pc,PCType type)
408: {
409: PetscErrorCode ierr,(*f)(PC,PCType);
413: PetscObjectQueryFunction((PetscObject)pc,"PCCompositeAddPC_C",(void (**)(void))&f);
414: if (f) {
415: (*f)(pc,type);
416: }
417: return(0);
418: }
422: /*@
423: PCCompositeGetPC - Gets one of the PC objects in the composite PC.
424:
425: Not Collective
427: Input Parameter:
428: + pc - the preconditioner context
429: - n - the number of the pc requested
431: Output Parameters:
432: . subpc - the PC requested
434: Level: Developer
436: .keywords: PC, get, composite preconditioner, sub preconditioner
438: .seealso: PCCompositeAddPC()
439: @*/
440: PetscErrorCode PCCompositeGetPC(PC pc,PetscInt n,PC *subpc)
441: {
442: PetscErrorCode ierr,(*f)(PC,PetscInt,PC *);
447: PetscObjectQueryFunction((PetscObject)pc,"PCCompositeGetPC_C",(void (**)(void))&f);
448: if (f) {
449: (*f)(pc,n,subpc);
450: } else {
451: SETERRQ(PETSC_ERR_ARG_WRONG,"Cannot get pc, not composite type");
452: }
453: return(0);
454: }
458: /*@
459: PCCompositeSetUseTrue - Sets a flag to indicate that the true matrix (rather than
460: the matrix used to define the preconditioner) is used to compute
461: the residual when the multiplicative scheme is used.
463: Collective on PC
465: Input Parameters:
466: . pc - the preconditioner context
468: Options Database Key:
469: . -pc_composite_true - Activates PCCompositeSetUseTrue()
471: Note:
472: For the common case in which the preconditioning and linear
473: system matrices are identical, this routine is unnecessary.
475: Level: Developer
477: .keywords: PC, composite preconditioner, set, true, flag
479: .seealso: PCSetOperators(), PCBJacobiSetUseTrueLocal(), PCKSPSetUseTrue()
480: @*/
481: PetscErrorCode PCCompositeSetUseTrue(PC pc)
482: {
483: PetscErrorCode ierr,(*f)(PC);
487: PetscObjectQueryFunction((PetscObject)pc,"PCCompositeSetUseTrue_C",(void (**)(void))&f);
488: if (f) {
489: (*f)(pc);
490: }
491: return(0);
492: }
494: /* -------------------------------------------------------------------------------------------*/
496: /*MC
497: PCCOMPOSITE - Build a preconditioner by composing together several preconditioners
499: Options Database Keys:
500: + -pc_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type
501: . -pc_composite_true - Activates PCCompositeSetUseTrue()
502: - -pc_composite_pcs - <pc0,pc1,...> list of PCs to compose
504: Level: intermediate
506: Concepts: composing solvers
508: Notes: To use a Krylov method inside the composite preconditioner, set the PCType of one or more
509: inner PCs to be PCKSP.
510: Using a Krylov method inside another Krylov method can be dangerous (you get divergence or
511: the incorrect answer) unless you use KSPFGMRES as the outter Krylov method
514: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
515: PCSHELL, PCKSP, PCCompositeSetType(), PCCompositeSpecialSetAlpha(), PCCompositeAddPC(),
516: PCCompositeGetPC(), PCCompositeSetUseTrue()
518: M*/
523: PetscErrorCode PCCreate_Composite(PC pc)
524: {
526: PC_Composite *jac;
529: PetscNew(PC_Composite,&jac);
530: PetscLogObjectMemory(pc,sizeof(PC_Composite));
531: pc->ops->apply = PCApply_Composite_Additive;
532: pc->ops->setup = PCSetUp_Composite;
533: pc->ops->destroy = PCDestroy_Composite;
534: pc->ops->setfromoptions = PCSetFromOptions_Composite;
535: pc->ops->view = PCView_Composite;
536: pc->ops->applyrichardson = 0;
538: pc->data = (void*)jac;
539: jac->type = PC_COMPOSITE_ADDITIVE;
540: jac->work1 = 0;
541: jac->work2 = 0;
542: jac->head = 0;
544: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCompositeSetType_C","PCCompositeSetType_Composite",
545: PCCompositeSetType_Composite);
546: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCompositeAddPC_C","PCCompositeAddPC_Composite",
547: PCCompositeAddPC_Composite);
548: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCompositeGetPC_C","PCCompositeGetPC_Composite",
549: PCCompositeGetPC_Composite);
550: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCompositeSetUseTrue_C","PCCompositeSetUseTrue_Composite",
551: PCCompositeSetUseTrue_Composite);
552: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCompositeSpecialSetAlpha_C","PCCompositeSpecialSetAlpha_Composite",
553: PCCompositeSpecialSetAlpha_Composite);
555: return(0);
556: }