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