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