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 PCGalerkinSetRestriction - Sets the restriction operator for the `PCGALERKIN` preconditioner 152 153 Logically Collective 154 155 Input Parameters: 156 + pc - the preconditioner context 157 - R - the restriction operator 158 159 Note: 160 Either this or `PCGalerkinSetInterpolation()` or both must be called 161 162 Level: Intermediate 163 164 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`, 165 `PCGalerkinSetInterpolation()`, `PCGalerkinGetKSP()` 166 @*/ 167 PetscErrorCode PCGalerkinSetRestriction(PC pc, Mat R) 168 { 169 PetscFunctionBegin; 170 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 171 PetscTryMethod(pc, "PCGalerkinSetRestriction_C", (PC, Mat), (pc, R)); 172 PetscFunctionReturn(0); 173 } 174 175 /*@ 176 PCGalerkinSetInterpolation - Sets the interpolation operator for the `PCGALERKIN` preconditioner 177 178 Logically Collective 179 180 Input Parameters: 181 + pc - the preconditioner context 182 - R - the interpolation operator 183 184 Note: 185 Either this or `PCGalerkinSetRestriction()` or both must be called 186 187 Level: Intermediate 188 189 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`, 190 `PCGalerkinSetRestriction()`, `PCGalerkinGetKSP()` 191 @*/ 192 PetscErrorCode PCGalerkinSetInterpolation(PC pc, Mat P) 193 { 194 PetscFunctionBegin; 195 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 196 PetscTryMethod(pc, "PCGalerkinSetInterpolation_C", (PC, Mat), (pc, P)); 197 PetscFunctionReturn(0); 198 } 199 200 /*@ 201 PCGalerkinSetComputeSubmatrix - Provide a routine that will be called to compute the Galerkin submatrix 202 203 Logically Collective 204 205 Input Parameters: 206 + pc - the preconditioner context 207 . computeAsub - routine that computes the submatrix from the global matrix 208 - ctx - context used by the routine, or NULL 209 210 Calling sequence of computeAsub: 211 $ computeAsub(PC pc,Mat A, Mat Ap, Mat *cAP,void *ctx); 212 213 + PC - the `PCGALERKIN` 214 . A - the matrix in the `PCGALERKIN` 215 . Ap - the computed submatrix from any previous computation, if NULL it has not previously been computed 216 . cAp - the submatrix computed by this routine 217 - ctx - optional user-defined function context 218 219 Level: Intermediate 220 221 Notes: 222 Instead of providing this routine you can call `PCGalerkinGetKSP()` and then `KSPSetOperators()` to provide the submatrix, 223 but that will not work for multiple `KSPSolve()`s with different matrices unless you call it for each solve. 224 225 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 226 matrix and computes its values in cAp. On each subsequent call the routine should up the Ap matrix. 227 228 Developer Note: 229 If the user does not call this routine nor call `PCGalerkinGetKSP()` and `KSPSetOperators()` then `PCGALERKIN` 230 could automatically compute the submatrix via calls to `MatGalerkin()` or `MatRARt()` 231 232 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`, 233 `PCGalerkinSetRestriction()`, `PCGalerkinSetInterpolation()`, `PCGalerkinGetKSP()` 234 @*/ 235 PetscErrorCode PCGalerkinSetComputeSubmatrix(PC pc, PetscErrorCode (*computeAsub)(PC, Mat, Mat, Mat *, void *), void *ctx) 236 { 237 PetscFunctionBegin; 238 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 239 PetscTryMethod(pc, "PCGalerkinSetComputeSubmatrix_C", (PC, PetscErrorCode(*)(PC, Mat, Mat, Mat *, void *), void *), (pc, computeAsub, ctx)); 240 PetscFunctionReturn(0); 241 } 242 243 /*@ 244 PCGalerkinGetKSP - Gets the `KSP` object in the `PCGALERKIN` 245 246 Not Collective 247 248 Input Parameter: 249 . pc - the preconditioner context 250 251 Output Parameter: 252 . ksp - the `KSP` object 253 254 Level: Intermediate 255 256 Note: 257 Once you have called this routine you can call `KSPSetOperators()` on the resulting ksp to provide the operator for the Galerkin problem, 258 an alternative is to use `PCGalerkinSetComputeSubmatrix()` to provide a routine that computes the submatrix as needed. 259 260 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`, 261 `PCGalerkinSetRestriction()`, `PCGalerkinSetInterpolation()`, `PCGalerkinSetComputeSubmatrix()` 262 @*/ 263 PetscErrorCode PCGalerkinGetKSP(PC pc, KSP *ksp) 264 { 265 PetscFunctionBegin; 266 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 267 PetscValidPointer(ksp, 2); 268 PetscUseMethod(pc, "PCGalerkinGetKSP_C", (PC, KSP *), (pc, ksp)); 269 PetscFunctionReturn(0); 270 } 271 272 static PetscErrorCode PCSetFromOptions_Galerkin(PC pc, PetscOptionItems *PetscOptionsObject) 273 { 274 PC_Galerkin *jac = (PC_Galerkin *)pc->data; 275 const char *prefix; 276 PetscBool flg; 277 278 PetscFunctionBegin; 279 PetscCall(KSPGetOptionsPrefix(jac->ksp, &prefix)); 280 PetscCall(PetscStrendswith(prefix, "galerkin_", &flg)); 281 if (!flg) { 282 PetscCall(PCGetOptionsPrefix(pc, &prefix)); 283 PetscCall(KSPSetOptionsPrefix(jac->ksp, prefix)); 284 PetscCall(KSPAppendOptionsPrefix(jac->ksp, "galerkin_")); 285 } 286 287 PetscOptionsHeadBegin(PetscOptionsObject, "Galerkin options"); 288 if (jac->ksp) PetscCall(KSPSetFromOptions(jac->ksp)); 289 PetscOptionsHeadEnd(); 290 PetscFunctionReturn(0); 291 } 292 293 /*MC 294 PCGALERKIN - Build (part of) a preconditioner by P S R (where P is often R^T) 295 296 Use `PCGalerkinSetRestriction`(pc,R) and/or `PCGalerkinSetInterpolation`(pc,P) followed by `PCGalerkinGetKSP`(pc,&ksp); `KSPSetOperators`(ksp,A,....) 297 298 Level: intermediate 299 300 Developer Notes: 301 If `KSPSetOperators()` has not been called on the inner `KSP` then `PCGALERKIN` could use `MatRARt()` or `MatPtAP()` to compute 302 the operators automatically. 303 304 Should there be a prefix for the inner `KSP`? 305 306 There is no `KSPSetFromOptions_Galerkin()` that calls `KSPSetFromOptions()` on the inner `KSP` 307 308 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, 309 `PCSHELL`, `PCKSP`, `PCGalerkinSetRestriction()`, `PCGalerkinSetInterpolation()`, `PCGalerkinGetKSP()` 310 M*/ 311 312 PETSC_EXTERN PetscErrorCode PCCreate_Galerkin(PC pc) 313 { 314 PC_Galerkin *jac; 315 316 PetscFunctionBegin; 317 PetscCall(PetscNew(&jac)); 318 319 pc->ops->apply = PCApply_Galerkin; 320 pc->ops->setup = PCSetUp_Galerkin; 321 pc->ops->reset = PCReset_Galerkin; 322 pc->ops->destroy = PCDestroy_Galerkin; 323 pc->ops->view = PCView_Galerkin; 324 pc->ops->setfromoptions = PCSetFromOptions_Galerkin; 325 pc->ops->applyrichardson = NULL; 326 327 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &jac->ksp)); 328 PetscCall(KSPSetErrorIfNotConverged(jac->ksp, pc->erroriffailure)); 329 PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->ksp, (PetscObject)pc, 1)); 330 331 pc->data = (void *)jac; 332 333 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGalerkinSetRestriction_C", PCGalerkinSetRestriction_Galerkin)); 334 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGalerkinSetInterpolation_C", PCGalerkinSetInterpolation_Galerkin)); 335 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGalerkinGetKSP_C", PCGalerkinGetKSP_Galerkin)); 336 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGalerkinSetComputeSubmatrix_C", PCGalerkinSetComputeSubmatrix_Galerkin)); 337 PetscFunctionReturn(0); 338 } 339