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