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