Actual source code: composite.c
petsc-3.3-p7 2013-05-11
2: /*
3: Defines a preconditioner that can consist of a collection of PCs
4: */
5: #include <petsc-private/pcimpl.h> /*I "petscpc.h" I*/
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: };
14:
15: typedef struct {
16: PC_CompositeLink head;
17: PCCompositeType type;
18: Vec work1;
19: Vec work2;
20: PetscScalar alpha;
21: PetscBool use_true_matrix;
22: } PC_Composite;
26: static PetscErrorCode PCApply_Composite_Multiplicative(PC pc,Vec x,Vec y)
27: {
28: PetscErrorCode ierr;
29: PC_Composite *jac = (PC_Composite*)pc->data;
30: PC_CompositeLink next = jac->head;
31: Mat mat = pc->pmat;
34: if (!next) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC() or -pc_composite_pcs");
35: if (next->next && !jac->work2) { /* allocate second work vector */
36: VecDuplicate(jac->work1,&jac->work2);
37: }
38: if (jac->use_true_matrix) mat = pc->mat;
39: PCApply(next->pc,x,y);
40: while (next->next) {
41: next = next->next;
42: MatMult(mat,y,jac->work1);
43: VecWAXPY(jac->work2,-1.0,jac->work1,x);
44: VecSet(jac->work1,0.0); /* zero since some PC's may not set all entries in the result */
45: PCApply(next->pc,jac->work2,jac->work1);
46: VecAXPY(y,1.0,jac->work1);
47: }
48: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
49: while (next->previous) {
50: next = next->previous;
51: MatMult(mat,y,jac->work1);
52: VecWAXPY(jac->work2,-1.0,jac->work1,x);
53: VecSet(jac->work1,0.0); /* zero since some PC's may not set all entries in the result */
54: PCApply(next->pc,jac->work2,jac->work1);
55: VecAXPY(y,1.0,jac->work1);
56: }
57: }
58: return(0);
59: }
63: static PetscErrorCode PCApplyTranspose_Composite_Multiplicative(PC pc,Vec x,Vec y)
64: {
65: PetscErrorCode ierr;
66: PC_Composite *jac = (PC_Composite*)pc->data;
67: PC_CompositeLink next = jac->head;
68: Mat mat = pc->pmat;
71: if (!next) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC() or -pc_composite_pcs");
72: if (next->next && !jac->work2) { /* allocate second work vector */
73: VecDuplicate(jac->work1,&jac->work2);
74: }
75: if (jac->use_true_matrix) mat = pc->mat;
76: /* locate last PC */
77: while (next->next) {
78: next = next->next;
79: }
80: PCApplyTranspose(next->pc,x,y);
81: while (next->previous) {
82: next = next->previous;
83: MatMultTranspose(mat,y,jac->work1);
84: VecWAXPY(jac->work2,-1.0,jac->work1,x);
85: VecSet(jac->work1,0.0); /* zero since some PC's may not set all entries in the result */
86: PCApplyTranspose(next->pc,jac->work2,jac->work1);
87: VecAXPY(y,1.0,jac->work1);
88: }
89: if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
90: next = jac->head;
91: while (next->next) {
92: next = next->next;
93: MatMultTranspose(mat,y,jac->work1);
94: VecWAXPY(jac->work2,-1.0,jac->work1,x);
95: VecSet(jac->work1,0.0); /* zero since some PC's may not set all entries in the result */
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: */
110: static PetscErrorCode PCApply_Composite_Special(PC pc,Vec x,Vec y)
111: {
112: PetscErrorCode ierr;
113: PC_Composite *jac = (PC_Composite*)pc->data;
114: PC_CompositeLink next = jac->head;
117: if (!next) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC() or -pc_composite_pcs");
118: if (!next->next || next->next->next) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Special composite preconditioners requires exactly two PCs");
120: PCApply(next->pc,x,jac->work1);
121: PCApply(next->next->pc,jac->work1,y);
122: return(0);
123: }
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(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC() or -pc_composite_pcs");
135: PCApply(next->pc,x,y);
136: while (next->next) {
137: next = next->next;
138: VecSet(jac->work1,0.0); /* zero since some PC's may not set all entries in the result */
139: PCApply(next->pc,x,jac->work1);
140: VecAXPY(y,1.0,jac->work1);
141: }
142: return(0);
143: }
147: static PetscErrorCode PCApplyTranspose_Composite_Additive(PC pc,Vec x,Vec y)
148: {
149: PetscErrorCode ierr;
150: PC_Composite *jac = (PC_Composite*)pc->data;
151: PC_CompositeLink next = jac->head;
154: if (!next) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC() or -pc_composite_pcs");
155: PCApplyTranspose(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: PCApplyTranspose(next->pc,x,jac->work1);
160: VecAXPY(y,1.0,jac->work1);
161: }
162: return(0);
163: }
167: static PetscErrorCode PCSetUp_Composite(PC pc)
168: {
169: PetscErrorCode ierr;
170: PC_Composite *jac = (PC_Composite*)pc->data;
171: PC_CompositeLink next = jac->head;
174: if (!jac->work1) {
175: MatGetVecs(pc->pmat,&jac->work1,0);
176: }
177: while (next) {
178: PCSetOperators(next->pc,pc->mat,pc->pmat,pc->flag);
179: next = next->next;
180: }
181: return(0);
182: }
186: static PetscErrorCode PCReset_Composite(PC pc)
187: {
188: PC_Composite *jac = (PC_Composite*)pc->data;
189: PetscErrorCode ierr;
190: PC_CompositeLink next = jac->head;
193: while (next) {
194: PCReset(next->pc);
195: next = next->next;
196: }
197: VecDestroy(&jac->work1);
198: VecDestroy(&jac->work2);
199: return(0);
200: }
204: static PetscErrorCode PCDestroy_Composite(PC pc)
205: {
206: PC_Composite *jac = (PC_Composite*)pc->data;
207: PetscErrorCode ierr;
208: PC_CompositeLink next = jac->head,next_tmp;
211: PCReset_Composite(pc);
212: while (next) {
213: PCDestroy(&next->pc);
214: next_tmp = next;
215: next = next->next;
216: PetscFree(next_tmp);
217: }
218: PetscFree(pc->data);
219: return(0);
220: }
224: static PetscErrorCode PCSetFromOptions_Composite(PC pc)
225: {
226: PC_Composite *jac = (PC_Composite*)pc->data;
227: PetscErrorCode ierr;
228: PetscInt nmax = 8,i;
229: PC_CompositeLink next;
230: char *pcs[8];
231: PetscBool flg;
234: PetscOptionsHead("Composite preconditioner options");
235: PetscOptionsEnum("-pc_composite_type","Type of composition","PCCompositeSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&jac->type,&flg);
236: if (flg) {
237: PCCompositeSetType(pc,jac->type);
238: }
239: flg = PETSC_FALSE;
240: PetscOptionsBool("-pc_composite_true","Use true matrix for inner solves","PCCompositeSetUseTrue",flg,&flg,PETSC_NULL);
241: if (flg) {
242: PCCompositeSetUseTrue(pc);
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: }
263: static PetscErrorCode PCView_Composite(PC pc,PetscViewer viewer)
264: {
265: PC_Composite *jac = (PC_Composite*)pc->data;
266: PetscErrorCode ierr;
267: PC_CompositeLink next = jac->head;
268: PetscBool iascii;
271: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
272: if (iascii) {
273: PetscViewerASCIIPrintf(viewer,"Composite PC type - %s\n",PCCompositeTypes[jac->type]);
274: PetscViewerASCIIPrintf(viewer,"PCs on composite preconditioner follow\n");
275: PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
276: } else {
277: SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for PCComposite",((PetscObject)viewer)->type_name);
278: }
279: if (iascii) {
280: PetscViewerASCIIPushTab(viewer);
281: }
282: while (next) {
283: PCView(next->pc,viewer);
284: next = next->next;
285: }
286: if (iascii) {
287: PetscViewerASCIIPopTab(viewer);
288: PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
289: }
290: return(0);
291: }
293: /* ------------------------------------------------------------------------------*/
295: EXTERN_C_BEGIN
298: PetscErrorCode PCCompositeSpecialSetAlpha_Composite(PC pc,PetscScalar alpha)
299: {
300: PC_Composite *jac = (PC_Composite*)pc->data;
302: jac->alpha = alpha;
303: return(0);
304: }
305: EXTERN_C_END
307: EXTERN_C_BEGIN
310: PetscErrorCode PCCompositeSetType_Composite(PC pc,PCCompositeType type)
311: {
312: PC_Composite *jac = (PC_Composite*)pc->data;
315: if (type == PC_COMPOSITE_ADDITIVE) {
316: pc->ops->apply = PCApply_Composite_Additive;
317: pc->ops->applytranspose = PCApplyTranspose_Composite_Additive;
318: } else if (type == PC_COMPOSITE_MULTIPLICATIVE || type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
319: pc->ops->apply = PCApply_Composite_Multiplicative;
320: pc->ops->applytranspose = PCApplyTranspose_Composite_Multiplicative;
321: } else if (type == PC_COMPOSITE_SPECIAL) {
322: pc->ops->apply = PCApply_Composite_Special;
323: pc->ops->applytranspose = PETSC_NULL;
324: } else SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG,"Unkown composite preconditioner type");
325: jac->type = type;
326: return(0);
327: }
328: EXTERN_C_END
330: EXTERN_C_BEGIN
333: PetscErrorCode PCCompositeAddPC_Composite(PC pc,const PCType type)
334: {
335: PC_Composite *jac;
336: PC_CompositeLink next,ilink;
337: PetscErrorCode ierr;
338: PetscInt cnt = 0;
339: const char *prefix;
340: char newprefix[8];
343: PetscNewLog(pc,struct _PC_CompositeLink,&ilink);
344: ilink->next = 0;
345: PCCreate(((PetscObject)pc)->comm,&ilink->pc);
346: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->pc);
348: jac = (PC_Composite*)pc->data;
349: next = jac->head;
350: if (!next) {
351: jac->head = ilink;
352: ilink->previous = PETSC_NULL;
353: } else {
354: cnt++;
355: while (next->next) {
356: next = next->next;
357: cnt++;
358: }
359: next->next = ilink;
360: ilink->previous = next;
361: }
362: PCGetOptionsPrefix(pc,&prefix);
363: PCSetOptionsPrefix(ilink->pc,prefix);
364: sprintf(newprefix,"sub_%d_",(int)cnt);
365: PCAppendOptionsPrefix(ilink->pc,newprefix);
366: /* type is set after prefix, because some methods may modify prefix, e.g. pcksp */
367: PCSetType(ilink->pc,type);
368: return(0);
369: }
370: EXTERN_C_END
372: EXTERN_C_BEGIN
375: PetscErrorCode PCCompositeGetPC_Composite(PC pc,PetscInt n,PC *subpc)
376: {
377: PC_Composite *jac;
378: PC_CompositeLink next;
379: PetscInt i;
382: jac = (PC_Composite*)pc->data;
383: next = jac->head;
384: for (i=0; i<n; i++) {
385: if (!next->next) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_INCOMP,"Not enough PCs in composite preconditioner");
386: next = next->next;
387: }
388: *subpc = next->pc;
389: return(0);
390: }
391: EXTERN_C_END
393: EXTERN_C_BEGIN
396: PetscErrorCode PCCompositeSetUseTrue_Composite(PC pc)
397: {
398: PC_Composite *jac;
401: jac = (PC_Composite*)pc->data;
402: jac->use_true_matrix = PETSC_TRUE;
403: return(0);
404: }
405: EXTERN_C_END
407: /* -------------------------------------------------------------------------------- */
410: /*@
411: PCCompositeSetType - Sets the type of composite preconditioner.
412:
413: Logically Collective on PC
415: Input Parameter:
416: + pc - the preconditioner context
417: - type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL
419: Options Database Key:
420: . -pc_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type
422: Level: Developer
424: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
425: @*/
426: PetscErrorCode PCCompositeSetType(PC pc,PCCompositeType type)
427: {
433: PetscTryMethod(pc,"PCCompositeSetType_C",(PC,PCCompositeType),(pc,type));
434: return(0);
435: }
439: /*@
440: PCCompositeSpecialSetAlpha - Sets alpha for the special composite preconditioner
441: for alphaI + R + S
442:
443: Logically Collective on PC
445: Input Parameter:
446: + pc - the preconditioner context
447: - alpha - scale on identity
449: Level: Developer
451: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
452: @*/
453: PetscErrorCode PCCompositeSpecialSetAlpha(PC pc,PetscScalar alpha)
454: {
460: PetscTryMethod(pc,"PCCompositeSpecialSetAlpha_C",(PC,PetscScalar),(pc,alpha));
461: return(0);
462: }
466: /*@C
467: PCCompositeAddPC - Adds another PC to the composite PC.
468:
469: Collective on PC
471: Input Parameters:
472: + pc - the preconditioner context
473: - type - the type of the new preconditioner
475: Level: Developer
477: .keywords: PC, composite preconditioner, add
478: @*/
479: PetscErrorCode PCCompositeAddPC(PC pc,const PCType type)
480: {
485: PetscTryMethod(pc,"PCCompositeAddPC_C",(PC,const PCType),(pc,type));
486: return(0);
487: }
491: /*@
492: PCCompositeGetPC - Gets one of the PC objects in the composite PC.
493:
494: Not Collective
496: Input Parameter:
497: + pc - the preconditioner context
498: - n - the number of the pc requested
500: Output Parameters:
501: . subpc - the PC requested
503: Level: Developer
505: .keywords: PC, get, composite preconditioner, sub preconditioner
507: .seealso: PCCompositeAddPC()
508: @*/
509: PetscErrorCode PCCompositeGetPC(PC pc,PetscInt n,PC *subpc)
510: {
516: PetscUseMethod(pc,"PCCompositeGetPC_C",(PC,PetscInt,PC *),(pc,n,subpc));
517: return(0);
518: }
522: /*@
523: PCCompositeSetUseTrue - Sets a flag to indicate that the true matrix (rather than
524: the matrix used to define the preconditioner) is used to compute
525: the residual when the multiplicative scheme is used.
527: Logically Collective on PC
529: Input Parameters:
530: . pc - the preconditioner context
532: Options Database Key:
533: . -pc_composite_true - Activates PCCompositeSetUseTrue()
535: Note:
536: For the common case in which the preconditioning and linear
537: system matrices are identical, this routine is unnecessary.
539: Level: Developer
541: .keywords: PC, composite preconditioner, set, true, flag
543: .seealso: PCSetOperators(), PCBJacobiSetUseTrueLocal(), PCKSPSetUseTrue()
544: @*/
545: PetscErrorCode PCCompositeSetUseTrue(PC pc)
546: {
551: PetscTryMethod(pc,"PCCompositeSetUseTrue_C",(PC),(pc));
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_composite_true - Activates PCCompositeSetUseTrue()
563: - -pc_composite_pcs - <pc0,pc1,...> list of PCs to compose
565: Level: intermediate
567: Concepts: composing solvers
569: Notes: To use a Krylov method inside the composite preconditioner, set the PCType of one or more
570: inner PCs to be PCKSP.
571: Using a Krylov method inside another Krylov method can be dangerous (you get divergence or
572: the incorrect answer) unless you use KSPFGMRES as the outer Krylov method
575: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
576: PCSHELL, PCKSP, PCCompositeSetType(), PCCompositeSpecialSetAlpha(), PCCompositeAddPC(),
577: PCCompositeGetPC(), PCCompositeSetUseTrue()
579: M*/
581: EXTERN_C_BEGIN
584: PetscErrorCode PCCreate_Composite(PC pc)
585: {
587: PC_Composite *jac;
590: PetscNewLog(pc,PC_Composite,&jac);
591: pc->ops->apply = PCApply_Composite_Additive;
592: pc->ops->applytranspose = PCApplyTranspose_Composite_Additive;
593: pc->ops->setup = PCSetUp_Composite;
594: pc->ops->reset = PCReset_Composite;
595: pc->ops->destroy = PCDestroy_Composite;
596: pc->ops->setfromoptions = PCSetFromOptions_Composite;
597: pc->ops->view = PCView_Composite;
598: pc->ops->applyrichardson = 0;
600: pc->data = (void*)jac;
601: jac->type = PC_COMPOSITE_ADDITIVE;
602: jac->work1 = 0;
603: jac->work2 = 0;
604: jac->head = 0;
606: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCompositeSetType_C","PCCompositeSetType_Composite",
607: PCCompositeSetType_Composite);
608: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCompositeAddPC_C","PCCompositeAddPC_Composite",
609: PCCompositeAddPC_Composite);
610: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCompositeGetPC_C","PCCompositeGetPC_Composite",
611: PCCompositeGetPC_Composite);
612: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCompositeSetUseTrue_C","PCCompositeSetUseTrue_Composite",
613: PCCompositeSetUseTrue_Composite);
614: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCompositeSpecialSetAlpha_C","PCCompositeSpecialSetAlpha_Composite",
615: PCCompositeSpecialSetAlpha_Composite);
617: return(0);
618: }
619: EXTERN_C_END