Actual source code: galerkin.c
petsc-3.5.4 2015-05-23
2: /*
3: Defines a preconditioner defined by R^T S R
4: */
5: #include <petsc-private/pcimpl.h>
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: {
19: PC_Galerkin *jac = (PC_Galerkin*)pc->data;
22: if (jac->R) {
23: MatRestrict(jac->R,x,jac->b);
24: } else {
25: MatRestrict(jac->P,x,jac->b);
26: }
27: KSPSolve(jac->ksp,jac->b,jac->x);
28: if (jac->P) {
29: MatInterpolate(jac->P,jac->x,y);
30: } else {
31: MatInterpolate(jac->R,jac->x,y);
32: }
33: return(0);
34: }
38: static PetscErrorCode PCSetUp_Galerkin(PC pc)
39: {
41: PC_Galerkin *jac = (PC_Galerkin*)pc->data;
42: PetscBool a;
43: Vec *xx,*yy;
46: if (!jac->x) {
47: KSPGetOperatorsSet(jac->ksp,&a,NULL);
48: if (!a) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set operator of PCGALERKIN KSP with PCGalerkinGetKSP()/KSPSetOperators()");
49: KSPGetVecs(jac->ksp,1,&xx,1,&yy);
50: jac->x = *xx;
51: jac->b = *yy;
52: PetscFree(xx);
53: PetscFree(yy);
54: }
55: if (!jac->R && !jac->P) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set restriction or interpolation of PCGALERKIN with PCGalerkinSetRestriction/Interpolation()");
56: /* should check here that sizes of R/P match size of a */
57: return(0);
58: }
62: static PetscErrorCode PCReset_Galerkin(PC pc)
63: {
64: PC_Galerkin *jac = (PC_Galerkin*)pc->data;
68: MatDestroy(&jac->R);
69: MatDestroy(&jac->P);
70: VecDestroy(&jac->x);
71: VecDestroy(&jac->b);
72: KSPReset(jac->ksp);
73: return(0);
74: }
78: static PetscErrorCode PCDestroy_Galerkin(PC pc)
79: {
80: PC_Galerkin *jac = (PC_Galerkin*)pc->data;
84: PCReset_Galerkin(pc);
85: KSPDestroy(&jac->ksp);
86: PetscFree(pc->data);
87: return(0);
88: }
92: static PetscErrorCode PCView_Galerkin(PC pc,PetscViewer viewer)
93: {
94: PC_Galerkin *jac = (PC_Galerkin*)pc->data;
96: PetscBool iascii;
99: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
100: if (iascii) {
101: PetscViewerASCIIPrintf(viewer,"Galerkin PC\n");
102: PetscViewerASCIIPrintf(viewer,"KSP on Galerkin follow\n");
103: PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
104: }
105: KSPView(jac->ksp,viewer);
106: return(0);
107: }
111: static PetscErrorCode PCGalerkinGetKSP_Galerkin(PC pc,KSP *ksp)
112: {
113: PC_Galerkin *jac = (PC_Galerkin*)pc->data;
116: *ksp = jac->ksp;
117: return(0);
118: }
122: static PetscErrorCode PCGalerkinSetRestriction_Galerkin(PC pc,Mat R)
123: {
124: PC_Galerkin *jac = (PC_Galerkin*)pc->data;
128: PetscObjectReference((PetscObject)R);
129: MatDestroy(&jac->R);
130: jac->R = R;
131: return(0);
132: }
136: static PetscErrorCode PCGalerkinSetInterpolation_Galerkin(PC pc,Mat P)
137: {
138: PC_Galerkin *jac = (PC_Galerkin*)pc->data;
142: PetscObjectReference((PetscObject)P);
143: MatDestroy(&jac->P);
144: jac->P = P;
145: return(0);
146: }
148: /* -------------------------------------------------------------------------------- */
151: /*@
152: PCGalerkinSetRestriction - Sets the restriction operator for the "Galerkin-type" preconditioner
154: Logically Collective on PC
156: Input Parameter:
157: + pc - the preconditioner context
158: - R - the restriction operator
160: Notes: Either this or PCGalerkinSetInterpolation() or both must be called
162: Level: Intermediate
164: .keywords: PC, set, Galerkin preconditioner
166: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN,
167: PCGalerkinSetInterpolation(), PCGalerkinGetKSP()
169: @*/
170: PetscErrorCode PCGalerkinSetRestriction(PC pc,Mat R)
171: {
176: PetscTryMethod(pc,"PCGalerkinSetRestriction_C",(PC,Mat),(pc,R));
177: return(0);
178: }
182: /*@
183: PCGalerkinSetInterpolation - Sets the interpolation operator for the "Galerkin-type" preconditioner
185: Logically Collective on PC
187: Input Parameter:
188: + pc - the preconditioner context
189: - R - the interpolation operator
191: Notes: Either this or PCGalerkinSetRestriction() or both must be called
193: Level: Intermediate
195: .keywords: PC, set, Galerkin preconditioner
197: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN,
198: PCGalerkinSetRestriction(), PCGalerkinGetKSP()
200: @*/
201: PetscErrorCode PCGalerkinSetInterpolation(PC pc,Mat P)
202: {
207: PetscTryMethod(pc,"PCGalerkinSetInterpolation_C",(PC,Mat),(pc,P));
208: return(0);
209: }
213: /*@
214: PCGalerkinGetKSP - Gets the KSP object in the Galerkin PC.
216: Not Collective
218: Input Parameter:
219: . pc - the preconditioner context
221: Output Parameters:
222: . ksp - the KSP object
224: Level: Intermediate
226: .keywords: PC, get, Galerkin preconditioner, sub preconditioner
228: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN,
229: PCGalerkinSetRestriction(), PCGalerkinSetInterpolation()
231: @*/
232: PetscErrorCode PCGalerkinGetKSP(PC pc,KSP *ksp)
233: {
239: PetscUseMethod(pc,"PCGalerkinGetKSP_C",(PC,KSP*),(pc,ksp));
240: return(0);
241: }
244: /* -------------------------------------------------------------------------------------------*/
246: /*MC
247: PCGALERKIN - Build (part of) a preconditioner by P S R (where P is often R^T)
249: $ Use PCGalerkinSetRestriction(pc,R) and/or PCGalerkinSetInterpolation(pc,P) followed by
250: $ PCGalerkinGetKSP(pc,&ksp); KSPSetOperators(ksp,A,....)
252: Level: intermediate
254: Developer Note: If KSPSetOperators() has not been called then PCGALERKIN could use MatRARt() or MatPtAP() to compute
255: the operators automatically.
256: Should there be a prefix for the inner KSP.
257: There is no KSPSetFromOptions_Galerkin() that calls KSPSetFromOptions() on the inner KSP
259: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
260: PCSHELL, PCKSP, PCGalerkinSetRestriction(), PCGalerkinSetInterpolation(), PCGalerkinGetKSP()
262: M*/
266: PETSC_EXTERN PetscErrorCode PCCreate_Galerkin(PC pc)
267: {
269: PC_Galerkin *jac;
272: PetscNewLog(pc,&jac);
274: pc->ops->apply = PCApply_Galerkin;
275: pc->ops->setup = PCSetUp_Galerkin;
276: pc->ops->reset = PCReset_Galerkin;
277: pc->ops->destroy = PCDestroy_Galerkin;
278: pc->ops->view = PCView_Galerkin;
279: pc->ops->applyrichardson = 0;
281: KSPCreate(PetscObjectComm((PetscObject)pc),&jac->ksp);
282: PetscObjectIncrementTabLevel((PetscObject)jac->ksp,(PetscObject)pc,1);
284: pc->data = (void*)jac;
286: PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinSetRestriction_C",PCGalerkinSetRestriction_Galerkin);
287: PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinSetInterpolation_C",PCGalerkinSetInterpolation_Galerkin);
288: PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinGetKSP_C",PCGalerkinGetKSP_Galerkin);
289: return(0);
290: }