Actual source code: galerkin.c
1: #define PETSCKSP_DLL
3: /*
4: Defines a preconditioner defined by R^T S R
5: */
6: #include private/pcimpl.h
7: #include petscksp.h
9: typedef struct {
10: KSP ksp;
11: Mat R,P;
12: Vec b,x;
13: } PC_Galerkin;
17: static PetscErrorCode PCApply_Galerkin(PC pc,Vec x,Vec y)
18: {
19: PetscErrorCode ierr;
20: PC_Galerkin *jac = (PC_Galerkin*)pc->data;
23: if (jac->R) {MatRestrict(jac->R,x,jac->b);}
24: else {MatRestrict(jac->P,x,jac->b);}
25: KSPSolve(jac->ksp,jac->b,jac->x);
26: if (jac->P) {MatInterpolate(jac->P,jac->x,y);}
27: else {MatInterpolate(jac->R,jac->x,y);}
28: return(0);
29: }
33: static PetscErrorCode PCSetUp_Galerkin(PC pc)
34: {
35: PetscErrorCode ierr;
36: PC_Galerkin *jac = (PC_Galerkin*)pc->data;
37: PetscTruth a;
38: Vec *xx,*yy;
41: if (!jac->x) {
42: KSPGetOperatorsSet(jac->ksp,&a,PETSC_NULL);
43: if (!a) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must set operator of PCGALERKIN KSP with PCGalerkinGetKSP()/KSPSetOperators()");
44: KSPGetVecs(jac->ksp,1,&xx,1,&yy);
45: jac->x = *xx;
46: jac->b = *yy;
47: PetscFree(xx);
48: PetscFree(yy);
49: }
50: if (!jac->R && !jac->P) {
51: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must set restriction or interpolation of PCGALERKIN with PCGalerkinSetRestriction/Interpolation()");
52: }
53: /* should check here that sizes of R/P match size of a */
54: return(0);
55: }
59: static PetscErrorCode PCDestroy_Galerkin(PC pc)
60: {
61: PC_Galerkin *jac = (PC_Galerkin*)pc->data;
62: PetscErrorCode ierr;
65: if (jac->R) {MatDestroy(jac->R);}
66: if (jac->P) {MatDestroy(jac->P);}
67: if (jac->x) {VecDestroy(jac->x);}
68: if (jac->b) {VecDestroy(jac->b);}
69: KSPDestroy(jac->ksp);
70: PetscFree(jac);
71: return(0);
72: }
76: static PetscErrorCode PCView_Galerkin(PC pc,PetscViewer viewer)
77: {
78: PC_Galerkin *jac = (PC_Galerkin*)pc->data;
79: PetscErrorCode ierr;
80: PetscTruth iascii;
83: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
84: if (iascii) {
85: PetscViewerASCIIPrintf(viewer,"Galerkin PC\n");
86: PetscViewerASCIIPrintf(viewer,"KSP on Galerkin follow\n");
87: PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
88: } else {
89: SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCGalerkin",((PetscObject)viewer)->type_name);
90: }
91: KSPView(jac->ksp,viewer);
92: return(0);
93: }
98: PetscErrorCode PCGalerkinGetKSP_Galerkin(PC pc,KSP *ksp)
99: {
100: PC_Galerkin *jac = (PC_Galerkin*)pc->data;
103: *ksp = jac->ksp;
104: return(0);
105: }
111: PetscErrorCode PCGalerkinSetRestriction_Galerkin(PC pc,Mat R)
112: {
113: PC_Galerkin *jac = (PC_Galerkin*)pc->data;
114: PetscErrorCode ierr;
117: if (jac->R) {MatDestroy(jac->R);}
118: jac->R = R;
119: PetscObjectReference((PetscObject)R);
120: return(0);
121: }
127: PetscErrorCode PCGalerkinSetInterpolation_Galerkin(PC pc,Mat P)
128: {
129: PC_Galerkin *jac = (PC_Galerkin*)pc->data;
130: PetscErrorCode ierr;
133: if (jac->P) {MatDestroy(jac->P);}
134: jac->P = P;
135: PetscObjectReference((PetscObject)P);
136: return(0);
137: }
140: /* -------------------------------------------------------------------------------- */
143: /*@
144: PCGalerkinSetRestriction - Sets the restriction operator for the "Galerkin-type" preconditioner
145:
146: Collective on PC
148: Input Parameter:
149: + pc - the preconditioner context
150: - R - the restriction operator
152: Notes: Either this or PCGalerkinSetInterpolation() or both must be called
154: Level: Intermediate
156: .keywords: PC, set, Galerkin preconditioner
158: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN,
159: PCGalerkinSetInterpolation(), PCGalerkinGetKSP()
161: @*/
162: PetscErrorCode PCGalerkinSetRestriction(PC pc,Mat R)
163: {
164: PetscErrorCode ierr,(*f)(PC,Mat);
168: PetscObjectQueryFunction((PetscObject)pc,"PCGalerkinSetRestriction_C",(void (**)(void))&f);
169: if (f) {
170: (*f)(pc,R);
171: }
172: return(0);
173: }
177: /*@
178: PCGalerkinSetInterpolation - Sets the interpolation operator for the "Galerkin-type" preconditioner
179:
180: Collective on PC
182: Input Parameter:
183: + pc - the preconditioner context
184: - R - the interpolation operator
186: Notes: Either this or PCGalerkinSetRestriction() or both must be called
188: Level: Intermediate
190: .keywords: PC, set, Galerkin preconditioner
192: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN,
193: PCGalerkinSetRestriction(), PCGalerkinGetKSP()
195: @*/
196: PetscErrorCode PCGalerkinSetInterpolation(PC pc,Mat P)
197: {
198: PetscErrorCode ierr,(*f)(PC,Mat);
202: PetscObjectQueryFunction((PetscObject)pc,"PCGalerkinSetInterpolation_C",(void (**)(void))&f);
203: if (f) {
204: (*f)(pc,P);
205: }
206: return(0);
207: }
211: /*@
212: PCGalerkinGetKSP - Gets the KSP object in the Galerkin PC.
213:
214: Not Collective
216: Input Parameter:
217: . pc - the preconditioner context
219: Output Parameters:
220: . ksp - the KSP object
222: Level: Intermediate
224: .keywords: PC, get, Galerkin preconditioner, sub preconditioner
226: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN,
227: PCGalerkinSetRestriction(), PCGalerkinSetInterpolation()
229: @*/
230: PetscErrorCode PCGalerkinGetKSP(PC pc,KSP *ksp)
231: {
232: PetscErrorCode ierr,(*f)(PC,KSP *);
237: PetscObjectQueryFunction((PetscObject)pc,"PCGalerkinGetKSP_C",(void (**)(void))&f);
238: if (f) {
239: (*f)(pc,ksp);
240: } else {
241: SETERRQ(PETSC_ERR_ARG_WRONG,"Cannot get KSP, not Galerkin type");
242: }
243: return(0);
244: }
247: /* -------------------------------------------------------------------------------------------*/
249: /*MC
250: PCGALERKIN - Build (part of) a preconditioner by P S R (where P is often R^T)
252: $ Use PCGalerkinSetRestriction(pc,R) and/or PCGalerkinSetInterpolation(pc,P) followed by
253: $ PCGalerkinGetKSP(pc,&ksp); KSPSetOperations(ksp,A,....)
255: Level: intermediate
257: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
258: PCSHELL, PCKSP, PCGalerkinSetRestriction(), PCGalerkinSetInterpolation(), PCGalerkinGetKSP()
260: M*/
265: PetscErrorCode PCCreate_Galerkin(PC pc)
266: {
268: PC_Galerkin *jac;
271: PetscNew(PC_Galerkin,&jac);
272: pc->ops->apply = PCApply_Galerkin;
273: pc->ops->setup = PCSetUp_Galerkin;
274: pc->ops->destroy = PCDestroy_Galerkin;
275: pc->ops->view = PCView_Galerkin;
276: pc->ops->applyrichardson = 0;
278: KSPCreate(pc->comm,&jac->ksp);
280: pc->data = (void*)jac;
282: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGalerkinSetRestriction_C","PCGalerkinSetRestriction_Galerkin",
283: PCGalerkinSetRestriction_Galerkin);
284: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGalerkinSetInterpolation_C","PCGalerkinSetInterpolation_Galerkin",
285: PCGalerkinSetInterpolation_Galerkin);
286: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGalerkinGetKSP_C","PCGalerkinGetKSP_Galerkin",
287: PCGalerkinGetKSP_Galerkin);
288: return(0);
289: }