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 PetscErrorCode (*computeasub)(PC,Mat,Mat,Mat*,void*); 13 void *computeasub_ctx; 14 } PC_Galerkin; 15 16 static PetscErrorCode PCApply_Galerkin(PC pc,Vec x,Vec y) 17 { 18 PC_Galerkin *jac = (PC_Galerkin*)pc->data; 19 20 PetscFunctionBegin; 21 if (jac->R) { 22 PetscCall(MatRestrict(jac->R,x,jac->b)); 23 } else { 24 PetscCall(MatRestrict(jac->P,x,jac->b)); 25 } 26 PetscCall(KSPSolve(jac->ksp,jac->b,jac->x)); 27 PetscCall(KSPCheckSolve(jac->ksp,pc,jac->x)); 28 if (jac->P) { 29 PetscCall(MatInterpolate(jac->P,jac->x,y)); 30 } else { 31 PetscCall(MatInterpolate(jac->R,jac->x,y)); 32 } 33 PetscFunctionReturn(0); 34 } 35 36 static PetscErrorCode PCSetUp_Galerkin(PC pc) 37 { 38 PC_Galerkin *jac = (PC_Galerkin*)pc->data; 39 PetscBool a; 40 Vec *xx,*yy; 41 42 PetscFunctionBegin; 43 if (jac->computeasub) { 44 Mat Ap; 45 if (!pc->setupcalled) { 46 PetscCall((*jac->computeasub)(pc,pc->pmat,NULL,&Ap,jac->computeasub_ctx)); 47 PetscCall(KSPSetOperators(jac->ksp,Ap,Ap)); 48 PetscCall(MatDestroy(&Ap)); 49 } else { 50 PetscCall(KSPGetOperators(jac->ksp,NULL,&Ap)); 51 PetscCall((*jac->computeasub)(pc,pc->pmat,Ap,NULL,jac->computeasub_ctx)); 52 } 53 } 54 55 if (!jac->x) { 56 PetscCall(KSPGetOperatorsSet(jac->ksp,&a,NULL)); 57 PetscCheck(a,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set operator of PCGALERKIN KSP with PCGalerkinGetKSP()/KSPSetOperators()"); 58 PetscCall(KSPCreateVecs(jac->ksp,1,&xx,1,&yy)); 59 jac->x = *xx; 60 jac->b = *yy; 61 PetscCall(PetscFree(xx)); 62 PetscCall(PetscFree(yy)); 63 } 64 PetscCheck(jac->R || jac->P,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set restriction or interpolation of PCGALERKIN with PCGalerkinSetRestriction()/Interpolation()"); 65 /* should check here that sizes of R/P match size of a */ 66 67 PetscFunctionReturn(0); 68 } 69 70 static PetscErrorCode PCReset_Galerkin(PC pc) 71 { 72 PC_Galerkin *jac = (PC_Galerkin*)pc->data; 73 74 PetscFunctionBegin; 75 PetscCall(MatDestroy(&jac->R)); 76 PetscCall(MatDestroy(&jac->P)); 77 PetscCall(VecDestroy(&jac->x)); 78 PetscCall(VecDestroy(&jac->b)); 79 PetscCall(KSPReset(jac->ksp)); 80 PetscFunctionReturn(0); 81 } 82 83 static PetscErrorCode PCDestroy_Galerkin(PC pc) 84 { 85 PC_Galerkin *jac = (PC_Galerkin*)pc->data; 86 87 PetscFunctionBegin; 88 PetscCall(PCReset_Galerkin(pc)); 89 PetscCall(KSPDestroy(&jac->ksp)); 90 PetscCall(PetscFree(pc->data)); 91 PetscFunctionReturn(0); 92 } 93 94 static PetscErrorCode PCView_Galerkin(PC pc,PetscViewer viewer) 95 { 96 PC_Galerkin *jac = (PC_Galerkin*)pc->data; 97 PetscBool iascii; 98 99 PetscFunctionBegin; 100 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 101 if (iascii) { 102 PetscCall(PetscViewerASCIIPrintf(viewer," KSP on Galerkin follow\n")); 103 PetscCall(PetscViewerASCIIPrintf(viewer," ---------------------------------\n")); 104 } 105 PetscCall(KSPView(jac->ksp,viewer)); 106 PetscFunctionReturn(0); 107 } 108 109 static PetscErrorCode PCGalerkinGetKSP_Galerkin(PC pc,KSP *ksp) 110 { 111 PC_Galerkin *jac = (PC_Galerkin*)pc->data; 112 113 PetscFunctionBegin; 114 *ksp = jac->ksp; 115 PetscFunctionReturn(0); 116 } 117 118 static PetscErrorCode PCGalerkinSetRestriction_Galerkin(PC pc,Mat R) 119 { 120 PC_Galerkin *jac = (PC_Galerkin*)pc->data; 121 122 PetscFunctionBegin; 123 PetscCall(PetscObjectReference((PetscObject)R)); 124 PetscCall(MatDestroy(&jac->R)); 125 jac->R = R; 126 PetscFunctionReturn(0); 127 } 128 129 static PetscErrorCode PCGalerkinSetInterpolation_Galerkin(PC pc,Mat P) 130 { 131 PC_Galerkin *jac = (PC_Galerkin*)pc->data; 132 133 PetscFunctionBegin; 134 PetscCall(PetscObjectReference((PetscObject)P)); 135 PetscCall(MatDestroy(&jac->P)); 136 jac->P = P; 137 PetscFunctionReturn(0); 138 } 139 140 static PetscErrorCode PCGalerkinSetComputeSubmatrix_Galerkin(PC pc,PetscErrorCode (*computeAsub)(PC,Mat,Mat,Mat*,void*),void *ctx) 141 { 142 PC_Galerkin *jac = (PC_Galerkin*)pc->data; 143 144 PetscFunctionBegin; 145 jac->computeasub = computeAsub; 146 jac->computeasub_ctx = ctx; 147 PetscFunctionReturn(0); 148 } 149 150 /* -------------------------------------------------------------------------------- */ 151 /*@ 152 PCGalerkinSetRestriction - Sets the restriction operator for the "Galerkin-type" preconditioner 153 154 Logically Collective on PC 155 156 Input Parameters: 157 + pc - the preconditioner context 158 - R - the restriction operator 159 160 Notes: 161 Either this or PCGalerkinSetInterpolation() or both must be called 162 163 Level: Intermediate 164 165 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`, 166 `PCGalerkinSetInterpolation()`, `PCGalerkinGetKSP()` 167 168 @*/ 169 PetscErrorCode PCGalerkinSetRestriction(PC pc,Mat R) 170 { 171 PetscFunctionBegin; 172 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 173 PetscTryMethod(pc,"PCGalerkinSetRestriction_C",(PC,Mat),(pc,R)); 174 PetscFunctionReturn(0); 175 } 176 177 /*@ 178 PCGalerkinSetInterpolation - Sets the interpolation operator for the "Galerkin-type" preconditioner 179 180 Logically Collective on PC 181 182 Input Parameters: 183 + pc - the preconditioner context 184 - R - the interpolation operator 185 186 Notes: 187 Either this or PCGalerkinSetRestriction() or both must be called 188 189 Level: Intermediate 190 191 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`, 192 `PCGalerkinSetRestriction()`, `PCGalerkinGetKSP()` 193 194 @*/ 195 PetscErrorCode PCGalerkinSetInterpolation(PC pc,Mat P) 196 { 197 PetscFunctionBegin; 198 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 199 PetscTryMethod(pc,"PCGalerkinSetInterpolation_C",(PC,Mat),(pc,P)); 200 PetscFunctionReturn(0); 201 } 202 203 /*@ 204 PCGalerkinSetComputeSubmatrix - Provide a routine that will be called to compute the Galerkin submatrix 205 206 Logically Collective 207 208 Input Parameters: 209 + pc - the preconditioner context 210 . computeAsub - routine that computes the submatrix from the global matrix 211 - ctx - context used by the routine, or NULL 212 213 Calling sequence of computeAsub: 214 $ computeAsub(PC pc,Mat A, Mat Ap, Mat *cAP,void *ctx); 215 216 + PC - the Galerkin PC 217 . A - the matrix in the Galerkin PC 218 . Ap - the computed submatrix from any previous computation, if NULL it has not previously been computed 219 . cAp - the submatrix computed by this routine 220 - ctx - optional user-defined function context 221 222 Level: Intermediate 223 224 Notes: 225 Instead of providing this routine you can call PCGalerkinGetKSP() and then KSPSetOperators() to provide the submatrix, 226 but that will not work for multiple KSPSolves with different matrices unless you call it for each solve. 227 228 This routine is called each time the outer matrix is changed. In the first call the Ap argument is NULL and the routine should create the 229 matrix and computes its values in cAp. On each subsequent call the routine should up the Ap matrix. 230 231 Developer Notes: 232 If the user does not call this routine nor call PCGalerkinGetKSP() and KSPSetOperators() then PCGalerkin could 233 could automatically compute the submatrix via calls to MatGalerkin() or MatRARt() 234 235 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`, 236 `PCGalerkinSetRestriction()`, `PCGalerkinSetInterpolation()`, `PCGalerkinGetKSP()` 237 238 @*/ 239 PetscErrorCode PCGalerkinSetComputeSubmatrix(PC pc,PetscErrorCode (*computeAsub)(PC,Mat,Mat,Mat*,void*),void *ctx) 240 { 241 PetscFunctionBegin; 242 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 243 PetscTryMethod(pc,"PCGalerkinSetComputeSubmatrix_C",(PC,PetscErrorCode (*)(PC,Mat,Mat,Mat*,void*),void*),(pc,computeAsub,ctx)); 244 PetscFunctionReturn(0); 245 } 246 247 /*@ 248 PCGalerkinGetKSP - Gets the KSP object in the Galerkin PC. 249 250 Not Collective 251 252 Input Parameter: 253 . pc - the preconditioner context 254 255 Output Parameters: 256 . ksp - the KSP object 257 258 Level: Intermediate 259 260 Notes: 261 Once you have called this routine you can call KSPSetOperators() on the resulting ksp to provide the operator for the Galerkin problem, 262 an alternative is to use PCGalerkinSetComputeSubmatrix() to provide a routine that computes the submatrix as needed. 263 264 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`, 265 `PCGalerkinSetRestriction()`, `PCGalerkinSetInterpolation()`, `PCGalerkinSetComputeSubmatrix()` 266 267 @*/ 268 PetscErrorCode PCGalerkinGetKSP(PC pc,KSP *ksp) 269 { 270 PetscFunctionBegin; 271 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 272 PetscValidPointer(ksp,2); 273 PetscUseMethod(pc,"PCGalerkinGetKSP_C",(PC,KSP*),(pc,ksp)); 274 PetscFunctionReturn(0); 275 } 276 277 static PetscErrorCode PCSetFromOptions_Galerkin(PetscOptionItems *PetscOptionsObject,PC pc) 278 { 279 PC_Galerkin *jac = (PC_Galerkin*)pc->data; 280 const char *prefix; 281 PetscBool flg; 282 283 PetscFunctionBegin; 284 PetscCall(KSPGetOptionsPrefix(jac->ksp,&prefix)); 285 PetscCall(PetscStrendswith(prefix,"galerkin_",&flg)); 286 if (!flg) { 287 PetscCall(PCGetOptionsPrefix(pc,&prefix)); 288 PetscCall(KSPSetOptionsPrefix(jac->ksp,prefix)); 289 PetscCall(KSPAppendOptionsPrefix(jac->ksp,"galerkin_")); 290 } 291 292 PetscOptionsHeadBegin(PetscOptionsObject,"Galerkin options"); 293 if (jac->ksp) PetscCall(KSPSetFromOptions(jac->ksp)); 294 PetscOptionsHeadEnd(); 295 PetscFunctionReturn(0); 296 } 297 298 /* -------------------------------------------------------------------------------------------*/ 299 300 /*MC 301 PCGALERKIN - Build (part of) a preconditioner by P S R (where P is often R^T) 302 303 $ Use PCGalerkinSetRestriction(pc,R) and/or PCGalerkinSetInterpolation(pc,P) followed by 304 $ PCGalerkinGetKSP(pc,&ksp); KSPSetOperators(ksp,A,....) 305 306 Level: intermediate 307 308 Developer Note: If KSPSetOperators() has not been called on the inner KSP then PCGALERKIN could use MatRARt() or MatPtAP() to compute 309 the operators automatically. 310 Should there be a prefix for the inner KSP. 311 There is no KSPSetFromOptions_Galerkin() that calls KSPSetFromOptions() on the inner KSP 312 313 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, 314 `PCSHELL`, `PCKSP`, `PCGalerkinSetRestriction()`, `PCGalerkinSetInterpolation()`, `PCGalerkinGetKSP()` 315 316 M*/ 317 318 PETSC_EXTERN PetscErrorCode PCCreate_Galerkin(PC pc) 319 { 320 PC_Galerkin *jac; 321 322 PetscFunctionBegin; 323 PetscCall(PetscNewLog(pc,&jac)); 324 325 pc->ops->apply = PCApply_Galerkin; 326 pc->ops->setup = PCSetUp_Galerkin; 327 pc->ops->reset = PCReset_Galerkin; 328 pc->ops->destroy = PCDestroy_Galerkin; 329 pc->ops->view = PCView_Galerkin; 330 pc->ops->setfromoptions = PCSetFromOptions_Galerkin; 331 pc->ops->applyrichardson = NULL; 332 333 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc),&jac->ksp)); 334 PetscCall(KSPSetErrorIfNotConverged(jac->ksp,pc->erroriffailure)); 335 PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->ksp,(PetscObject)pc,1)); 336 337 pc->data = (void*)jac; 338 339 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinSetRestriction_C",PCGalerkinSetRestriction_Galerkin)); 340 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinSetInterpolation_C",PCGalerkinSetInterpolation_Galerkin)); 341 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinGetKSP_C",PCGalerkinGetKSP_Galerkin)); 342 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinSetComputeSubmatrix_C",PCGalerkinSetComputeSubmatrix_Galerkin)); 343 PetscFunctionReturn(0); 344 } 345