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