1 2 /* 3 Defines a preconditioner defined by R^T S R 4 */ 5 #include <petsc-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) 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 } else { 99 SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for PCGalerkin",((PetscObject)viewer)->type_name); 100 } 101 ierr = KSPView(jac->ksp,viewer);CHKERRQ(ierr); 102 PetscFunctionReturn(0); 103 } 104 105 EXTERN_C_BEGIN 106 #undef __FUNCT__ 107 #define __FUNCT__ "PCGalerkinGetKSP_Galerkin" 108 PetscErrorCode PCGalerkinGetKSP_Galerkin(PC pc,KSP *ksp) 109 { 110 PC_Galerkin *jac = (PC_Galerkin*)pc->data; 111 112 PetscFunctionBegin; 113 *ksp = jac->ksp; 114 PetscFunctionReturn(0); 115 } 116 EXTERN_C_END 117 118 EXTERN_C_BEGIN 119 #undef __FUNCT__ 120 #define __FUNCT__ "PCGalerkinSetRestriction_Galerkin" 121 PetscErrorCode PCGalerkinSetRestriction_Galerkin(PC pc,Mat R) 122 { 123 PC_Galerkin *jac = (PC_Galerkin*)pc->data; 124 PetscErrorCode ierr; 125 126 PetscFunctionBegin; 127 ierr = PetscObjectReference((PetscObject)R);CHKERRQ(ierr); 128 ierr = MatDestroy(&jac->R);CHKERRQ(ierr); 129 jac->R = R; 130 PetscFunctionReturn(0); 131 } 132 EXTERN_C_END 133 134 EXTERN_C_BEGIN 135 #undef __FUNCT__ 136 #define __FUNCT__ "PCGalerkinSetInterpolation_Galerkin" 137 PetscErrorCode PCGalerkinSetInterpolation_Galerkin(PC pc,Mat P) 138 { 139 PC_Galerkin *jac = (PC_Galerkin*)pc->data; 140 PetscErrorCode ierr; 141 142 PetscFunctionBegin; 143 ierr = PetscObjectReference((PetscObject)P);CHKERRQ(ierr); 144 ierr = MatDestroy(&jac->P);CHKERRQ(ierr); 145 jac->P = P; 146 PetscFunctionReturn(0); 147 } 148 EXTERN_C_END 149 150 /* -------------------------------------------------------------------------------- */ 151 #undef __FUNCT__ 152 #define __FUNCT__ "PCGalerkinSetRestriction" 153 /*@ 154 PCGalerkinSetRestriction - Sets the restriction operator for the "Galerkin-type" preconditioner 155 156 Logically Collective on PC 157 158 Input Parameter: 159 + pc - the preconditioner context 160 - R - the restriction operator 161 162 Notes: Either this or PCGalerkinSetInterpolation() or both must be called 163 164 Level: Intermediate 165 166 .keywords: PC, set, Galerkin preconditioner 167 168 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN, 169 PCGalerkinSetInterpolation(), PCGalerkinGetKSP() 170 171 @*/ 172 PetscErrorCode PCGalerkinSetRestriction(PC pc,Mat R) 173 { 174 PetscErrorCode ierr; 175 176 PetscFunctionBegin; 177 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 178 ierr = PetscTryMethod(pc,"PCGalerkinSetRestriction_C",(PC,Mat),(pc,R));CHKERRQ(ierr); 179 PetscFunctionReturn(0); 180 } 181 182 #undef __FUNCT__ 183 #define __FUNCT__ "PCGalerkinSetInterpolation" 184 /*@ 185 PCGalerkinSetInterpolation - Sets the interpolation operator for the "Galerkin-type" preconditioner 186 187 Logically Collective on PC 188 189 Input Parameter: 190 + pc - the preconditioner context 191 - R - the interpolation operator 192 193 Notes: Either this or PCGalerkinSetRestriction() or both must be called 194 195 Level: Intermediate 196 197 .keywords: PC, set, Galerkin preconditioner 198 199 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN, 200 PCGalerkinSetRestriction(), PCGalerkinGetKSP() 201 202 @*/ 203 PetscErrorCode PCGalerkinSetInterpolation(PC pc,Mat P) 204 { 205 PetscErrorCode ierr; 206 207 PetscFunctionBegin; 208 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 209 ierr = PetscTryMethod(pc,"PCGalerkinSetInterpolation_C",(PC,Mat),(pc,P));CHKERRQ(ierr); 210 PetscFunctionReturn(0); 211 } 212 213 #undef __FUNCT__ 214 #define __FUNCT__ "PCGalerkinGetKSP" 215 /*@ 216 PCGalerkinGetKSP - Gets the KSP object in the Galerkin PC. 217 218 Not Collective 219 220 Input Parameter: 221 . pc - the preconditioner context 222 223 Output Parameters: 224 . ksp - the KSP object 225 226 Level: Intermediate 227 228 .keywords: PC, get, Galerkin preconditioner, sub preconditioner 229 230 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN, 231 PCGalerkinSetRestriction(), PCGalerkinSetInterpolation() 232 233 @*/ 234 PetscErrorCode PCGalerkinGetKSP(PC pc,KSP *ksp) 235 { 236 PetscErrorCode ierr; 237 238 PetscFunctionBegin; 239 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 240 PetscValidPointer(ksp,2); 241 ierr = PetscUseMethod(pc,"PCGalerkinGetKSP_C",(PC,KSP *),(pc,ksp));CHKERRQ(ierr); 242 PetscFunctionReturn(0); 243 } 244 245 246 /* -------------------------------------------------------------------------------------------*/ 247 248 /*MC 249 PCGALERKIN - Build (part of) a preconditioner by P S R (where P is often R^T) 250 251 $ Use PCGalerkinSetRestriction(pc,R) and/or PCGalerkinSetInterpolation(pc,P) followed by 252 $ PCGalerkinGetKSP(pc,&ksp); KSPSetOperators(ksp,A,....) 253 254 Level: intermediate 255 256 Developer Note: If KSPSetOperators() has not been called then PCGALERKIN could use MatRARt() or MatPtAP() to compute 257 the operators automatically. 258 Should there be a prefix for the inner KSP. 259 There is no KSPSetFromOptions_Galerkin() that calls KSPSetFromOptions() on the inner KSP 260 261 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 262 PCSHELL, PCKSP, PCGalerkinSetRestriction(), PCGalerkinSetInterpolation(), PCGalerkinGetKSP() 263 264 M*/ 265 266 EXTERN_C_BEGIN 267 #undef __FUNCT__ 268 #define __FUNCT__ "PCCreate_Galerkin" 269 PetscErrorCode PCCreate_Galerkin(PC pc) 270 { 271 PetscErrorCode ierr; 272 PC_Galerkin *jac; 273 274 PetscFunctionBegin; 275 ierr = PetscNewLog(pc,PC_Galerkin,&jac);CHKERRQ(ierr); 276 pc->ops->apply = PCApply_Galerkin; 277 pc->ops->setup = PCSetUp_Galerkin; 278 pc->ops->reset = PCReset_Galerkin; 279 pc->ops->destroy = PCDestroy_Galerkin; 280 pc->ops->view = PCView_Galerkin; 281 pc->ops->applyrichardson = 0; 282 283 ierr = KSPCreate(((PetscObject)pc)->comm,&jac->ksp);CHKERRQ(ierr); 284 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 285 286 pc->data = (void*)jac; 287 288 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGalerkinSetRestriction_C","PCGalerkinSetRestriction_Galerkin", 289 PCGalerkinSetRestriction_Galerkin);CHKERRQ(ierr); 290 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGalerkinSetInterpolation_C","PCGalerkinSetInterpolation_Galerkin", 291 PCGalerkinSetInterpolation_Galerkin);CHKERRQ(ierr); 292 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGalerkinGetKSP_C","PCGalerkinGetKSP_Galerkin", 293 PCGalerkinGetKSP_Galerkin);CHKERRQ(ierr); 294 PetscFunctionReturn(0); 295 } 296 EXTERN_C_END 297 298