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