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 PetscBool 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(((PetscObject)pc)->comm,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(((PetscObject)pc)->comm,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(pc->data);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 PetscBool iascii; 81 82 PetscFunctionBegin; 83 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&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(((PetscObject)pc)->comm,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 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 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 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 Logically 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 PCGalerkinSetRestriction(PC pc,Mat R) 163 { 164 PetscErrorCode ierr; 165 166 PetscFunctionBegin; 167 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 168 ierr = PetscTryMethod(pc,"PCGalerkinSetRestriction_C",(PC,Mat),(pc,R));CHKERRQ(ierr); 169 PetscFunctionReturn(0); 170 } 171 172 #undef __FUNCT__ 173 #define __FUNCT__ "PCGalerkinSetRestriction" 174 /*@ 175 PCGalerkinSetInterpolation - Sets the interpolation operator for the "Galerkin-type" preconditioner 176 177 Logically Collective on PC 178 179 Input Parameter: 180 + pc - the preconditioner context 181 - R - the interpolation operator 182 183 Notes: Either this or PCGalerkinSetRestriction() or both must be called 184 185 Level: Intermediate 186 187 .keywords: PC, set, Galerkin preconditioner 188 189 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN, 190 PCGalerkinSetRestriction(), PCGalerkinGetKSP() 191 192 @*/ 193 PetscErrorCode PCGalerkinSetInterpolation(PC pc,Mat P) 194 { 195 PetscErrorCode ierr; 196 197 PetscFunctionBegin; 198 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 199 ierr = PetscTryMethod(pc,"PCGalerkinSetInterpolation_C",(PC,Mat),(pc,P));CHKERRQ(ierr); 200 PetscFunctionReturn(0); 201 } 202 203 #undef __FUNCT__ 204 #define __FUNCT__ "PCGalerkinGetKSP" 205 /*@ 206 PCGalerkinGetKSP - Gets the KSP object in the Galerkin PC. 207 208 Not Collective 209 210 Input Parameter: 211 . pc - the preconditioner context 212 213 Output Parameters: 214 . ksp - the KSP object 215 216 Level: Intermediate 217 218 .keywords: PC, get, Galerkin preconditioner, sub preconditioner 219 220 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN, 221 PCGalerkinSetRestriction(), PCGalerkinSetInterpolation() 222 223 @*/ 224 PetscErrorCode PCGalerkinGetKSP(PC pc,KSP *ksp) 225 { 226 PetscErrorCode ierr; 227 228 PetscFunctionBegin; 229 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 230 PetscValidPointer(ksp,2); 231 ierr = PetscUseMethod(pc,"PCGalerkinGetKSP_C",(PC,KSP *),(pc,ksp));CHKERRQ(ierr); 232 PetscFunctionReturn(0); 233 } 234 235 236 /* -------------------------------------------------------------------------------------------*/ 237 238 /*MC 239 PCGALERKIN - Build (part of) a preconditioner by P S R (where P is often R^T) 240 241 $ Use PCGalerkinSetRestriction(pc,R) and/or PCGalerkinSetInterpolation(pc,P) followed by 242 $ PCGalerkinGetKSP(pc,&ksp); KSPSetOperations(ksp,A,....) 243 244 Level: intermediate 245 246 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 247 PCSHELL, PCKSP, PCGalerkinSetRestriction(), PCGalerkinSetInterpolation(), PCGalerkinGetKSP() 248 249 M*/ 250 251 EXTERN_C_BEGIN 252 #undef __FUNCT__ 253 #define __FUNCT__ "PCCreate_Galerkin" 254 PetscErrorCode PCCreate_Galerkin(PC pc) 255 { 256 PetscErrorCode ierr; 257 PC_Galerkin *jac; 258 259 PetscFunctionBegin; 260 ierr = PetscNewLog(pc,PC_Galerkin,&jac);CHKERRQ(ierr); 261 pc->ops->apply = PCApply_Galerkin; 262 pc->ops->setup = PCSetUp_Galerkin; 263 pc->ops->destroy = PCDestroy_Galerkin; 264 pc->ops->view = PCView_Galerkin; 265 pc->ops->applyrichardson = 0; 266 267 ierr = KSPCreate(((PetscObject)pc)->comm,&jac->ksp);CHKERRQ(ierr); 268 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 269 270 pc->data = (void*)jac; 271 272 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGalerkinSetRestriction_C","PCGalerkinSetRestriction_Galerkin", 273 PCGalerkinSetRestriction_Galerkin);CHKERRQ(ierr); 274 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGalerkinSetInterpolation_C","PCGalerkinSetInterpolation_Galerkin", 275 PCGalerkinSetInterpolation_Galerkin);CHKERRQ(ierr); 276 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGalerkinGetKSP_C","PCGalerkinGetKSP_Galerkin", 277 PCGalerkinGetKSP_Galerkin);CHKERRQ(ierr); 278 PetscFunctionReturn(0); 279 } 280 EXTERN_C_END 281 282