Actual source code: galerkin.c
petsc-3.3-p7 2013-05-11
2: /*
3: Defines a preconditioner defined by R^T S R
4: */
5: #include <petsc-private/pcimpl.h> /*I "petscpc.h" I*/
6: #include <petscksp.h> /*I "petscksp.h" I*/
8: typedef struct {
9: KSP ksp;
10: Mat R,P;
11: Vec b,x;
12: } PC_Galerkin;
16: static PetscErrorCode PCApply_Galerkin(PC pc,Vec x,Vec y)
17: {
18: PetscErrorCode ierr;
19: PC_Galerkin *jac = (PC_Galerkin*)pc->data;
22: if (jac->R) {MatRestrict(jac->R,x,jac->b);}
23: else {MatRestrict(jac->P,x,jac->b);}
24: KSPSolve(jac->ksp,jac->b,jac->x);
25: if (jac->P) {MatInterpolate(jac->P,jac->x,y);}
26: else {MatInterpolate(jac->R,jac->x,y);}
27: return(0);
28: }
32: static PetscErrorCode PCSetUp_Galerkin(PC pc)
33: {
34: PetscErrorCode ierr;
35: PC_Galerkin *jac = (PC_Galerkin*)pc->data;
36: PetscBool a;
37: Vec *xx,*yy;
40: if (!jac->x) {
41: KSPGetOperatorsSet(jac->ksp,&a,PETSC_NULL);
42: if (!a) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set operator of PCGALERKIN KSP with PCGalerkinGetKSP()/KSPSetOperators()");
43: KSPGetVecs(jac->ksp,1,&xx,1,&yy);
44: jac->x = *xx;
45: jac->b = *yy;
46: PetscFree(xx);
47: PetscFree(yy);
48: }
49: if (!jac->R && !jac->P) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set restriction or interpolation of PCGALERKIN with PCGalerkinSetRestriction/Interpolation()");
50: /* should check here that sizes of R/P match size of a */
51: return(0);
52: }
56: static PetscErrorCode PCReset_Galerkin(PC pc)
57: {
58: PC_Galerkin *jac = (PC_Galerkin*)pc->data;
59: PetscErrorCode ierr;
62: MatDestroy(&jac->R);
63: MatDestroy(&jac->P);
64: VecDestroy(&jac->x);
65: VecDestroy(&jac->b);
66: KSPReset(jac->ksp);
67: return(0);
68: }
72: static PetscErrorCode PCDestroy_Galerkin(PC pc)
73: {
74: PC_Galerkin *jac = (PC_Galerkin*)pc->data;
75: PetscErrorCode ierr;
78: PCReset_Galerkin(pc);
79: KSPDestroy(&jac->ksp);
80: PetscFree(pc->data);
81: return(0);
82: }
86: static PetscErrorCode PCView_Galerkin(PC pc,PetscViewer viewer)
87: {
88: PC_Galerkin *jac = (PC_Galerkin*)pc->data;
89: PetscErrorCode ierr;
90: PetscBool iascii;
93: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
94: if (iascii) {
95: PetscViewerASCIIPrintf(viewer,"Galerkin PC\n");
96: PetscViewerASCIIPrintf(viewer,"KSP on Galerkin follow\n");
97: PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
98: } else {
99: SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for PCGalerkin",((PetscObject)viewer)->type_name);
100: }
101: KSPView(jac->ksp,viewer);
102: return(0);
103: }
105: EXTERN_C_BEGIN
108: PetscErrorCode PCGalerkinGetKSP_Galerkin(PC pc,KSP *ksp)
109: {
110: PC_Galerkin *jac = (PC_Galerkin*)pc->data;
113: *ksp = jac->ksp;
114: return(0);
115: }
116: EXTERN_C_END
118: EXTERN_C_BEGIN
121: PetscErrorCode PCGalerkinSetRestriction_Galerkin(PC pc,Mat R)
122: {
123: PC_Galerkin *jac = (PC_Galerkin*)pc->data;
124: PetscErrorCode ierr;
127: PetscObjectReference((PetscObject)R);
128: MatDestroy(&jac->R);
129: jac->R = R;
130: return(0);
131: }
132: EXTERN_C_END
134: EXTERN_C_BEGIN
137: PetscErrorCode PCGalerkinSetInterpolation_Galerkin(PC pc,Mat P)
138: {
139: PC_Galerkin *jac = (PC_Galerkin*)pc->data;
140: PetscErrorCode ierr;
143: PetscObjectReference((PetscObject)P);
144: MatDestroy(&jac->P);
145: jac->P = P;
146: return(0);
147: }
148: EXTERN_C_END
150: /* -------------------------------------------------------------------------------- */
153: /*@
154: PCGalerkinSetRestriction - Sets the restriction operator for the "Galerkin-type" preconditioner
155:
156: Logically Collective on PC
158: Input Parameter:
159: + pc - the preconditioner context
160: - R - the restriction operator
162: Notes: Either this or PCGalerkinSetInterpolation() or both must be called
164: Level: Intermediate
166: .keywords: PC, set, Galerkin preconditioner
168: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN,
169: PCGalerkinSetInterpolation(), PCGalerkinGetKSP()
171: @*/
172: PetscErrorCode PCGalerkinSetRestriction(PC pc,Mat R)
173: {
178: PetscTryMethod(pc,"PCGalerkinSetRestriction_C",(PC,Mat),(pc,R));
179: return(0);
180: }
184: /*@
185: PCGalerkinSetInterpolation - Sets the interpolation operator for the "Galerkin-type" preconditioner
186:
187: Logically Collective on PC
189: Input Parameter:
190: + pc - the preconditioner context
191: - R - the interpolation operator
193: Notes: Either this or PCGalerkinSetRestriction() or both must be called
195: Level: Intermediate
197: .keywords: PC, set, Galerkin preconditioner
199: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN,
200: PCGalerkinSetRestriction(), PCGalerkinGetKSP()
202: @*/
203: PetscErrorCode PCGalerkinSetInterpolation(PC pc,Mat P)
204: {
209: PetscTryMethod(pc,"PCGalerkinSetInterpolation_C",(PC,Mat),(pc,P));
210: return(0);
211: }
215: /*@
216: PCGalerkinGetKSP - Gets the KSP object in the Galerkin PC.
217:
218: Not Collective
220: Input Parameter:
221: . pc - the preconditioner context
223: Output Parameters:
224: . ksp - the KSP object
226: Level: Intermediate
228: .keywords: PC, get, Galerkin preconditioner, sub preconditioner
230: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN,
231: PCGalerkinSetRestriction(), PCGalerkinSetInterpolation()
233: @*/
234: PetscErrorCode PCGalerkinGetKSP(PC pc,KSP *ksp)
235: {
241: PetscUseMethod(pc,"PCGalerkinGetKSP_C",(PC,KSP *),(pc,ksp));
242: return(0);
243: }
246: /* -------------------------------------------------------------------------------------------*/
248: /*MC
249: PCGALERKIN - Build (part of) a preconditioner by P S R (where P is often R^T)
251: $ Use PCGalerkinSetRestriction(pc,R) and/or PCGalerkinSetInterpolation(pc,P) followed by
252: $ PCGalerkinGetKSP(pc,&ksp); KSPSetOperators(ksp,A,....)
254: Level: intermediate
256: Developer Note: If KSPSetOperators() has not been called then PCGALERKIN could use MatRARt() or MatPtAP() to compute
257: the operators automatically.
258: Should there be a prefix for the inner KSP.
259: There is no KSPSetFromOptions_Galerkin() that calls KSPSetFromOptions() on the inner KSP
261: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
262: PCSHELL, PCKSP, PCGalerkinSetRestriction(), PCGalerkinSetInterpolation(), PCGalerkinGetKSP()
264: M*/
266: EXTERN_C_BEGIN
269: PetscErrorCode PCCreate_Galerkin(PC pc)
270: {
272: PC_Galerkin *jac;
275: PetscNewLog(pc,PC_Galerkin,&jac);
276: pc->ops->apply = PCApply_Galerkin;
277: pc->ops->setup = PCSetUp_Galerkin;
278: pc->ops->reset = PCReset_Galerkin;
279: pc->ops->destroy = PCDestroy_Galerkin;
280: pc->ops->view = PCView_Galerkin;
281: pc->ops->applyrichardson = 0;
283: KSPCreate(((PetscObject)pc)->comm,&jac->ksp);
284: PetscObjectIncrementTabLevel((PetscObject)jac->ksp,(PetscObject)pc,1);
286: pc->data = (void*)jac;
288: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGalerkinSetRestriction_C","PCGalerkinSetRestriction_Galerkin",
289: PCGalerkinSetRestriction_Galerkin);
290: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGalerkinSetInterpolation_C","PCGalerkinSetInterpolation_Galerkin",
291: PCGalerkinSetInterpolation_Galerkin);
292: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGalerkinGetKSP_C","PCGalerkinGetKSP_Galerkin",
293: PCGalerkinGetKSP_Galerkin);
294: return(0);
295: }
296: EXTERN_C_END