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