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