1 2 /* 3 Defines a preconditioner defined by R^T S R 4 */ 5 #include "private/pcimpl.h" /*I "petscpc.h" I*/ 6 #include "petscksp.h" /*I "petscksp.h" I*/ 7 8 typedef struct { 9 KSP ksp; 10 Mat R,P; 11 Vec b,x; 12 } PC_Galerkin; 13 14 #undef __FUNCT__ 15 #define __FUNCT__ "PCApply_Galerkin" 16 static PetscErrorCode PCApply_Galerkin(PC pc,Vec x,Vec y) 17 { 18 PetscErrorCode ierr; 19 PC_Galerkin *jac = (PC_Galerkin*)pc->data; 20 21 PetscFunctionBegin; 22 if (jac->R) {ierr = MatRestrict(jac->R,x,jac->b);CHKERRQ(ierr);} 23 else {ierr = MatRestrict(jac->P,x,jac->b);CHKERRQ(ierr);} 24 ierr = KSPSolve(jac->ksp,jac->b,jac->x);CHKERRQ(ierr); 25 if (jac->P) {ierr = MatInterpolate(jac->P,jac->x,y);CHKERRQ(ierr);} 26 else {ierr = MatInterpolate(jac->R,jac->x,y);CHKERRQ(ierr);} 27 PetscFunctionReturn(0); 28 } 29 30 #undef __FUNCT__ 31 #define __FUNCT__ "PCSetUp_Galerkin" 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; 38 39 PetscFunctionBegin; 40 if (!jac->x) { 41 ierr = KSPGetOperatorsSet(jac->ksp,&a,PETSC_NULL);CHKERRQ(ierr); 42 if (!a) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set operator of PCGALERKIN KSP with PCGalerkinGetKSP()/KSPSetOperators()"); 43 ierr = KSPGetVecs(jac->ksp,1,&xx,1,&yy);CHKERRQ(ierr); 44 jac->x = *xx; 45 jac->b = *yy; 46 ierr = PetscFree(xx);CHKERRQ(ierr); 47 ierr = PetscFree(yy);CHKERRQ(ierr); 48 } 49 if (!jac->R && !jac->P) { 50 SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set restriction or interpolation of PCGALERKIN with PCGalerkinSetRestriction/Interpolation()"); 51 } 52 /* should check here that sizes of R/P match size of a */ 53 PetscFunctionReturn(0); 54 } 55 56 #undef __FUNCT__ 57 #define __FUNCT__ "PCDestroy_Galerkin" 58 static PetscErrorCode PCDestroy_Galerkin(PC pc) 59 { 60 PC_Galerkin *jac = (PC_Galerkin*)pc->data; 61 PetscErrorCode ierr; 62 63 PetscFunctionBegin; 64 if (jac->R) {ierr = MatDestroy(jac->R);CHKERRQ(ierr);} 65 if (jac->P) {ierr = MatDestroy(jac->P);CHKERRQ(ierr);} 66 if (jac->x) {ierr = VecDestroy(jac->x);CHKERRQ(ierr);} 67 if (jac->b) {ierr = VecDestroy(jac->b);CHKERRQ(ierr);} 68 ierr = KSPDestroy(jac->ksp);CHKERRQ(ierr); 69 ierr = PetscFree(pc->data);CHKERRQ(ierr); 70 PetscFunctionReturn(0); 71 } 72 73 #undef __FUNCT__ 74 #define __FUNCT__ "PCView_Galerkin" 75 static PetscErrorCode PCView_Galerkin(PC pc,PetscViewer viewer) 76 { 77 PC_Galerkin *jac = (PC_Galerkin*)pc->data; 78 PetscErrorCode ierr; 79 PetscBool iascii; 80 81 PetscFunctionBegin; 82 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 83 if (iascii) { 84 ierr = PetscViewerASCIIPrintf(viewer,"Galerkin PC\n");CHKERRQ(ierr); 85 ierr = PetscViewerASCIIPrintf(viewer,"KSP on Galerkin follow\n");CHKERRQ(ierr); 86 ierr = PetscViewerASCIIPrintf(viewer,"---------------------------------\n");CHKERRQ(ierr); 87 } else { 88 SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for PCGalerkin",((PetscObject)viewer)->type_name); 89 } 90 ierr = KSPView(jac->ksp,viewer);CHKERRQ(ierr); 91 PetscFunctionReturn(0); 92 } 93 94 EXTERN_C_BEGIN 95 #undef __FUNCT__ 96 #define __FUNCT__ "PCGalerkinGetKSP_Galerkin" 97 PetscErrorCode PCGalerkinGetKSP_Galerkin(PC pc,KSP *ksp) 98 { 99 PC_Galerkin *jac = (PC_Galerkin*)pc->data; 100 101 PetscFunctionBegin; 102 *ksp = jac->ksp; 103 PetscFunctionReturn(0); 104 } 105 EXTERN_C_END 106 107 EXTERN_C_BEGIN 108 #undef __FUNCT__ 109 #define __FUNCT__ "PCGalerkinSetRestriction_Galerkin" 110 PetscErrorCode PCGalerkinSetRestriction_Galerkin(PC pc,Mat R) 111 { 112 PC_Galerkin *jac = (PC_Galerkin*)pc->data; 113 PetscErrorCode ierr; 114 115 PetscFunctionBegin; 116 ierr = PetscObjectReference((PetscObject)R);CHKERRQ(ierr); 117 if (jac->R) {ierr = MatDestroy(jac->R);CHKERRQ(ierr);} 118 jac->R = R; 119 PetscFunctionReturn(0); 120 } 121 EXTERN_C_END 122 123 EXTERN_C_BEGIN 124 #undef __FUNCT__ 125 #define __FUNCT__ "PCGalerkinSetInterpolation_Galerkin" 126 PetscErrorCode PCGalerkinSetInterpolation_Galerkin(PC pc,Mat P) 127 { 128 PC_Galerkin *jac = (PC_Galerkin*)pc->data; 129 PetscErrorCode ierr; 130 131 PetscFunctionBegin; 132 ierr = PetscObjectReference((PetscObject)P);CHKERRQ(ierr); 133 if (jac->P) {ierr = MatDestroy(jac->P);CHKERRQ(ierr);} 134 jac->P = P; 135 PetscFunctionReturn(0); 136 } 137 EXTERN_C_END 138 139 /* -------------------------------------------------------------------------------- */ 140 #undef __FUNCT__ 141 #define __FUNCT__ "PCGalerkinSetRestriction" 142 /*@ 143 PCGalerkinSetRestriction - Sets the restriction operator for the "Galerkin-type" preconditioner 144 145 Logically Collective on PC 146 147 Input Parameter: 148 + pc - the preconditioner context 149 - R - the restriction operator 150 151 Notes: Either this or PCGalerkinSetInterpolation() or both must be called 152 153 Level: Intermediate 154 155 .keywords: PC, set, Galerkin preconditioner 156 157 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN, 158 PCGalerkinSetInterpolation(), PCGalerkinGetKSP() 159 160 @*/ 161 PetscErrorCode PCGalerkinSetRestriction(PC pc,Mat R) 162 { 163 PetscErrorCode ierr; 164 165 PetscFunctionBegin; 166 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 167 ierr = PetscTryMethod(pc,"PCGalerkinSetRestriction_C",(PC,Mat),(pc,R));CHKERRQ(ierr); 168 PetscFunctionReturn(0); 169 } 170 171 #undef __FUNCT__ 172 #define __FUNCT__ "PCGalerkinSetRestriction" 173 /*@ 174 PCGalerkinSetInterpolation - Sets the interpolation operator for the "Galerkin-type" preconditioner 175 176 Logically Collective on PC 177 178 Input Parameter: 179 + pc - the preconditioner context 180 - R - the interpolation operator 181 182 Notes: Either this or PCGalerkinSetRestriction() or both must be called 183 184 Level: Intermediate 185 186 .keywords: PC, set, Galerkin preconditioner 187 188 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN, 189 PCGalerkinSetRestriction(), PCGalerkinGetKSP() 190 191 @*/ 192 PetscErrorCode PCGalerkinSetInterpolation(PC pc,Mat P) 193 { 194 PetscErrorCode ierr; 195 196 PetscFunctionBegin; 197 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 198 ierr = PetscTryMethod(pc,"PCGalerkinSetInterpolation_C",(PC,Mat),(pc,P));CHKERRQ(ierr); 199 PetscFunctionReturn(0); 200 } 201 202 #undef __FUNCT__ 203 #define __FUNCT__ "PCGalerkinGetKSP" 204 /*@ 205 PCGalerkinGetKSP - Gets the KSP object in the Galerkin PC. 206 207 Not Collective 208 209 Input Parameter: 210 . pc - the preconditioner context 211 212 Output Parameters: 213 . ksp - the KSP object 214 215 Level: Intermediate 216 217 .keywords: PC, get, Galerkin preconditioner, sub preconditioner 218 219 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN, 220 PCGalerkinSetRestriction(), PCGalerkinSetInterpolation() 221 222 @*/ 223 PetscErrorCode PCGalerkinGetKSP(PC pc,KSP *ksp) 224 { 225 PetscErrorCode ierr; 226 227 PetscFunctionBegin; 228 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 229 PetscValidPointer(ksp,2); 230 ierr = PetscUseMethod(pc,"PCGalerkinGetKSP_C",(PC,KSP *),(pc,ksp));CHKERRQ(ierr); 231 PetscFunctionReturn(0); 232 } 233 234 235 /* -------------------------------------------------------------------------------------------*/ 236 237 /*MC 238 PCGALERKIN - Build (part of) a preconditioner by P S R (where P is often R^T) 239 240 $ Use PCGalerkinSetRestriction(pc,R) and/or PCGalerkinSetInterpolation(pc,P) followed by 241 $ PCGalerkinGetKSP(pc,&ksp); KSPSetOperations(ksp,A,....) 242 243 Level: intermediate 244 245 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 246 PCSHELL, PCKSP, PCGalerkinSetRestriction(), PCGalerkinSetInterpolation(), PCGalerkinGetKSP() 247 248 M*/ 249 250 EXTERN_C_BEGIN 251 #undef __FUNCT__ 252 #define __FUNCT__ "PCCreate_Galerkin" 253 PetscErrorCode PCCreate_Galerkin(PC pc) 254 { 255 PetscErrorCode ierr; 256 PC_Galerkin *jac; 257 258 PetscFunctionBegin; 259 ierr = PetscNewLog(pc,PC_Galerkin,&jac);CHKERRQ(ierr); 260 pc->ops->apply = PCApply_Galerkin; 261 pc->ops->setup = PCSetUp_Galerkin; 262 pc->ops->destroy = PCDestroy_Galerkin; 263 pc->ops->view = PCView_Galerkin; 264 pc->ops->applyrichardson = 0; 265 266 ierr = KSPCreate(((PetscObject)pc)->comm,&jac->ksp);CHKERRQ(ierr); 267 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 268 269 pc->data = (void*)jac; 270 271 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGalerkinSetRestriction_C","PCGalerkinSetRestriction_Galerkin", 272 PCGalerkinSetRestriction_Galerkin);CHKERRQ(ierr); 273 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGalerkinSetInterpolation_C","PCGalerkinSetInterpolation_Galerkin", 274 PCGalerkinSetInterpolation_Galerkin);CHKERRQ(ierr); 275 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGalerkinGetKSP_C","PCGalerkinGetKSP_Galerkin", 276 PCGalerkinGetKSP_Galerkin);CHKERRQ(ierr); 277 PetscFunctionReturn(0); 278 } 279 EXTERN_C_END 280 281