1 #include <petsc/private/kspimpl.h> /*I "petscksp.h" I*/
2
3 typedef struct {
4 PetscInt restart;
5 PetscInt n_restarts;
6 PetscScalar *val;
7 Vec *VV, *SS;
8 Vec R;
9
10 KSPFlexibleModifyPCFn *modifypc; /* function to modify the preconditioner*/
11 PetscCtxDestroyFn *modifypc_destroy; /* function to destroy the user context for the modifypc function */
12
13 void *modifypc_ctx; /* user defined data for the modifypc function */
14 } KSP_GCR;
15
KSPSolve_GCR_cycle(KSP ksp)16 static PetscErrorCode KSPSolve_GCR_cycle(KSP ksp)
17 {
18 KSP_GCR *ctx = (KSP_GCR *)ksp->data;
19 PetscScalar r_dot_v;
20 Mat A, B;
21 PC pc;
22 Vec s, v, r;
23 /*
24 The residual norm will not be computed when ksp->its > ksp->chknorm hence need to initialize norm_r with some dummy value
25 */
26 PetscReal norm_r = 0.0, nrm;
27 PetscInt k, i, restart;
28 Vec x;
29
30 PetscFunctionBegin;
31 restart = ctx->restart;
32 PetscCall(KSPGetPC(ksp, &pc));
33 PetscCall(KSPGetOperators(ksp, &A, &B));
34
35 x = ksp->vec_sol;
36 r = ctx->R;
37
38 for (k = 0; k < restart; k++) {
39 v = ctx->VV[k];
40 s = ctx->SS[k];
41 if (ctx->modifypc) PetscCall((*ctx->modifypc)(ksp, ksp->its, k, ksp->rnorm, ctx->modifypc_ctx));
42
43 PetscCall(KSP_PCApply(ksp, r, s)); /* s = B^{-1} r */
44 PetscCall(KSP_MatMult(ksp, A, s, v)); /* v = A s */
45
46 PetscCall(VecMDot(v, k, ctx->VV, ctx->val));
47 for (i = 0; i < k; i++) ctx->val[i] = -ctx->val[i];
48 PetscCall(VecMAXPY(v, k, ctx->val, ctx->VV)); /* v = v - sum_{i=0}^{k-1} alpha_i v_i */
49 PetscCall(VecMAXPY(s, k, ctx->val, ctx->SS)); /* s = s - sum_{i=0}^{k-1} alpha_i s_i */
50
51 PetscCall(VecDotNorm2(r, v, &r_dot_v, &nrm));
52 nrm = PetscSqrtReal(nrm);
53 r_dot_v = r_dot_v / nrm;
54 PetscCall(VecScale(v, 1.0 / nrm));
55 PetscCall(VecScale(s, 1.0 / nrm));
56 PetscCall(VecAXPY(x, r_dot_v, s));
57 PetscCall(VecAXPY(r, -r_dot_v, v));
58 if (ksp->its > ksp->chknorm && ksp->normtype != KSP_NORM_NONE) {
59 PetscCall(VecNorm(r, NORM_2, &norm_r));
60 KSPCheckNorm(ksp, norm_r);
61 }
62 /* update the local counter and the global counter */
63 ksp->its++;
64 ksp->rnorm = norm_r;
65
66 PetscCall(KSPLogResidualHistory(ksp, norm_r));
67 PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
68
69 if (ksp->its - 1 > ksp->chknorm) {
70 PetscCall((*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP));
71 if (ksp->reason) break;
72 }
73
74 if (ksp->its >= ksp->max_it) {
75 ksp->reason = KSP_CONVERGED_ITS;
76 break;
77 }
78 }
79 ctx->n_restarts++;
80 PetscFunctionReturn(PETSC_SUCCESS);
81 }
82
KSPSolve_GCR(KSP ksp)83 static PetscErrorCode KSPSolve_GCR(KSP ksp)
84 {
85 KSP_GCR *ctx = (KSP_GCR *)ksp->data;
86 Mat A, B;
87 Vec r, b, x;
88 PetscReal norm_r = 0.0;
89
90 PetscFunctionBegin;
91 PetscCall(KSPGetOperators(ksp, &A, &B));
92 x = ksp->vec_sol;
93 b = ksp->vec_rhs;
94 r = ctx->R;
95
96 /* compute initial residual */
97 PetscCall(KSP_MatMult(ksp, A, x, r));
98 PetscCall(VecAYPX(r, -1.0, b)); /* r = b - A x */
99 if (ksp->normtype != KSP_NORM_NONE) {
100 PetscCall(VecNorm(r, NORM_2, &norm_r));
101 KSPCheckNorm(ksp, norm_r);
102 }
103 ksp->its = 0;
104 ksp->rnorm0 = norm_r;
105
106 PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm0));
107 PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm0));
108 PetscCall((*ksp->converged)(ksp, ksp->its, ksp->rnorm0, &ksp->reason, ksp->cnvP));
109 if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
110
111 do {
112 PetscCall(KSPSolve_GCR_cycle(ksp));
113 if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS); /* catch case when convergence occurs inside the cycle */
114 } while (ksp->its < ksp->max_it);
115
116 if (ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
117 PetscFunctionReturn(PETSC_SUCCESS);
118 }
119
KSPView_GCR(KSP ksp,PetscViewer viewer)120 static PetscErrorCode KSPView_GCR(KSP ksp, PetscViewer viewer)
121 {
122 KSP_GCR *ctx = (KSP_GCR *)ksp->data;
123 PetscBool isascii;
124
125 PetscFunctionBegin;
126 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
127 if (isascii) {
128 PetscCall(PetscViewerASCIIPrintf(viewer, " restart = %" PetscInt_FMT " \n", ctx->restart));
129 PetscCall(PetscViewerASCIIPrintf(viewer, " restarts performed = %" PetscInt_FMT " \n", ctx->n_restarts));
130 }
131 PetscFunctionReturn(PETSC_SUCCESS);
132 }
133
KSPSetUp_GCR(KSP ksp)134 static PetscErrorCode KSPSetUp_GCR(KSP ksp)
135 {
136 KSP_GCR *ctx = (KSP_GCR *)ksp->data;
137 Mat A;
138 PetscBool diagonalscale;
139
140 PetscFunctionBegin;
141 PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
142 PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
143
144 PetscCall(KSPGetOperators(ksp, &A, NULL));
145 PetscCall(MatCreateVecs(A, &ctx->R, NULL));
146 PetscCall(VecDuplicateVecs(ctx->R, ctx->restart, &ctx->VV));
147 PetscCall(VecDuplicateVecs(ctx->R, ctx->restart, &ctx->SS));
148
149 PetscCall(PetscMalloc1(ctx->restart, &ctx->val));
150 PetscFunctionReturn(PETSC_SUCCESS);
151 }
152
KSPReset_GCR(KSP ksp)153 static PetscErrorCode KSPReset_GCR(KSP ksp)
154 {
155 KSP_GCR *ctx = (KSP_GCR *)ksp->data;
156
157 PetscFunctionBegin;
158 PetscCall(VecDestroy(&ctx->R));
159 PetscCall(VecDestroyVecs(ctx->restart, &ctx->VV));
160 PetscCall(VecDestroyVecs(ctx->restart, &ctx->SS));
161 if (ctx->modifypc_destroy) PetscCall((*ctx->modifypc_destroy)(&ctx->modifypc_ctx));
162 PetscCall(PetscFree(ctx->val));
163 PetscFunctionReturn(PETSC_SUCCESS);
164 }
165
KSPDestroy_GCR(KSP ksp)166 static PetscErrorCode KSPDestroy_GCR(KSP ksp)
167 {
168 PetscFunctionBegin;
169 PetscCall(KSPReset_GCR(ksp));
170 PetscCall(KSPDestroyDefault(ksp));
171 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGCRSetRestart_C", NULL));
172 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGCRGetRestart_C", NULL));
173 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPFlexibleSetModifyPC_C", NULL));
174 PetscFunctionReturn(PETSC_SUCCESS);
175 }
176
KSPSetFromOptions_GCR(KSP ksp,PetscOptionItems PetscOptionsObject)177 static PetscErrorCode KSPSetFromOptions_GCR(KSP ksp, PetscOptionItems PetscOptionsObject)
178 {
179 KSP_GCR *ctx = (KSP_GCR *)ksp->data;
180 PetscInt restart;
181 PetscBool flg;
182
183 PetscFunctionBegin;
184 PetscOptionsHeadBegin(PetscOptionsObject, "KSP GCR options");
185 PetscCall(PetscOptionsInt("-ksp_gcr_restart", "Number of Krylov search directions", "KSPGCRSetRestart", ctx->restart, &restart, &flg));
186 if (flg) PetscCall(KSPGCRSetRestart(ksp, restart));
187 PetscOptionsHeadEnd();
188 PetscFunctionReturn(PETSC_SUCCESS);
189 }
190
KSPFlexibleSetModifyPC_GCR(KSP ksp,KSPFlexibleModifyPCFn * function,PetscCtx ctx,PetscCtxDestroyFn * destroy)191 static PetscErrorCode KSPFlexibleSetModifyPC_GCR(KSP ksp, KSPFlexibleModifyPCFn *function, PetscCtx ctx, PetscCtxDestroyFn *destroy)
192 {
193 KSP_GCR *gcr = (KSP_GCR *)ksp->data;
194
195 PetscFunctionBegin;
196 PetscValidHeaderSpecific(ksp, KSP_CLASSID, 1);
197 gcr->modifypc = function;
198 gcr->modifypc_destroy = destroy;
199 gcr->modifypc_ctx = ctx;
200 PetscFunctionReturn(PETSC_SUCCESS);
201 }
202
KSPGCRSetRestart_GCR(KSP ksp,PetscInt restart)203 static PetscErrorCode KSPGCRSetRestart_GCR(KSP ksp, PetscInt restart)
204 {
205 KSP_GCR *ctx;
206
207 PetscFunctionBegin;
208 ctx = (KSP_GCR *)ksp->data;
209 ctx->restart = restart;
210 PetscFunctionReturn(PETSC_SUCCESS);
211 }
212
KSPGCRGetRestart_GCR(KSP ksp,PetscInt * restart)213 static PetscErrorCode KSPGCRGetRestart_GCR(KSP ksp, PetscInt *restart)
214 {
215 KSP_GCR *ctx;
216
217 PetscFunctionBegin;
218 ctx = (KSP_GCR *)ksp->data;
219 *restart = ctx->restart;
220 PetscFunctionReturn(PETSC_SUCCESS);
221 }
222
223 /*@
224 KSPGCRSetRestart - Sets number of iterations at which `KSPGCR` restarts.
225
226 Not Collective
227
228 Input Parameters:
229 + ksp - the Krylov space context
230 - restart - integer restart value
231
232 Options Database Key:
233 . -ksp_gcr_restart <restart> - the number of stored vectors to orthogonalize against
234
235 Level: intermediate
236
237 Note:
238 The default value is 30.
239
240 Developer Note:
241 The API could be made uniform for all `KSP` methods that have a restart.
242
243 .seealso: [](ch_ksp), `KSPGCR`, `KSPSetTolerances()`, `KSPGCRGetRestart()`, `KSPGMRESSetRestart()`
244 @*/
KSPGCRSetRestart(KSP ksp,PetscInt restart)245 PetscErrorCode KSPGCRSetRestart(KSP ksp, PetscInt restart)
246 {
247 PetscFunctionBegin;
248 PetscTryMethod(ksp, "KSPGCRSetRestart_C", (KSP, PetscInt), (ksp, restart));
249 PetscFunctionReturn(PETSC_SUCCESS);
250 }
251
252 /*@
253 KSPGCRGetRestart - Gets number of iterations at which `KSPGCR` restarts.
254
255 Not Collective
256
257 Input Parameter:
258 . ksp - the Krylov space context
259
260 Output Parameter:
261 . restart - integer restart value
262
263 Level: intermediate
264
265 .seealso: [](ch_ksp), `KSPGCR`, `KSPSetTolerances()`, `KSPGCRSetRestart()`, `KSPGMRESGetRestart()`
266 @*/
KSPGCRGetRestart(KSP ksp,PetscInt * restart)267 PetscErrorCode KSPGCRGetRestart(KSP ksp, PetscInt *restart)
268 {
269 PetscFunctionBegin;
270 PetscTryMethod(ksp, "KSPGCRGetRestart_C", (KSP, PetscInt *), (ksp, restart));
271 PetscFunctionReturn(PETSC_SUCCESS);
272 }
273
KSPBuildSolution_GCR(KSP ksp,Vec v,Vec * V)274 static PetscErrorCode KSPBuildSolution_GCR(KSP ksp, Vec v, Vec *V)
275 {
276 Vec x;
277
278 PetscFunctionBegin;
279 x = ksp->vec_sol;
280 if (v) {
281 PetscCall(VecCopy(x, v));
282 if (V) *V = v;
283 } else if (V) {
284 *V = ksp->vec_sol;
285 }
286 PetscFunctionReturn(PETSC_SUCCESS);
287 }
288
KSPBuildResidual_GCR(KSP ksp,Vec t,Vec v,Vec * V)289 static PetscErrorCode KSPBuildResidual_GCR(KSP ksp, Vec t, Vec v, Vec *V)
290 {
291 KSP_GCR *ctx;
292
293 PetscFunctionBegin;
294 ctx = (KSP_GCR *)ksp->data;
295 if (v) {
296 PetscCall(VecCopy(ctx->R, v));
297 if (V) *V = v;
298 } else if (V) {
299 *V = ctx->R;
300 }
301 PetscFunctionReturn(PETSC_SUCCESS);
302 }
303
304 /*MC
305 KSPGCR - Implements the preconditioned flexible Generalized Conjugate Residual (GCR) method {cite}`eisenstat1983variational`. [](sec_flexibleksp)
306
307 Options Database Key:
308 . -ksp_gcr_restart <restart> - the number of stored vectors to orthogonalize against
309
310 Level: beginner
311
312 Notes:
313 This method supports non-symmetric matrices and permits the use of a preconditioner
314 which may vary from one iteration to the next.
315
316 Users can define a method to vary the preconditioner between iterates via `KSPFlexibleSetModifyPC()`.
317
318 When a restart occurs, the initial starting solution is given by the current estimate for `x`.
319
320 Unlike `KSPGMRES` and `KSPFGMRES`, when using GCR, the solution and residual vector can be directly accessed at any iterate,
321 with zero computational cost, via a call to `KSPBuildSolution()` and `KSPBuildResidual()` respectively.
322
323 The stopping condition test is only applied after the iteration count exceeds the value set by
324 `KSPSetCheckNormIteration()` (also available as `-ksp_check_norm_iteration`). The residual norm
325 reported by the monitor and stored in the residual history will be listed as 0.0 before that
326 iteration; the norm is not actually zero, it is simply not computed until then.
327
328 The method implemented requires the storage of 2 x restart + 1 vectors, twice as much as `KSPGMRES`.
329
330 Supports only right preconditioning.
331
332 Contributed by:
333 Dave May
334
335 .seealso: [](ch_ksp), [](sec_flexibleksp), `KSPFCG`, `KSPPIPEGCR`, `KSPPIPEFCG`, `KSPFGMRES`, `KSPCG`, `KSPCreate()`, `KSPSetType()`, `KSPType`,
336 `KSP`, `KSPGCRSetRestart()`, `KSPGCRGetRestart()`, `KSPFlexibleSetModifyPC()`, `KSPGMRES`
337 M*/
KSPCreate_GCR(KSP ksp)338 PETSC_EXTERN PetscErrorCode KSPCreate_GCR(KSP ksp)
339 {
340 KSP_GCR *ctx;
341
342 PetscFunctionBegin;
343 PetscCall(PetscNew(&ctx));
344
345 ctx->restart = 30;
346 ctx->n_restarts = 0;
347 ksp->data = (void *)ctx;
348
349 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));
350 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 3));
351
352 ksp->ops->setup = KSPSetUp_GCR;
353 ksp->ops->solve = KSPSolve_GCR;
354 ksp->ops->reset = KSPReset_GCR;
355 ksp->ops->destroy = KSPDestroy_GCR;
356 ksp->ops->view = KSPView_GCR;
357 ksp->ops->setfromoptions = KSPSetFromOptions_GCR;
358 ksp->ops->buildsolution = KSPBuildSolution_GCR;
359 ksp->ops->buildresidual = KSPBuildResidual_GCR;
360
361 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGCRSetRestart_C", KSPGCRSetRestart_GCR));
362 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGCRGetRestart_C", KSPGCRGetRestart_GCR));
363 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPFlexibleSetModifyPC_C", KSPFlexibleSetModifyPC_GCR));
364 PetscFunctionReturn(PETSC_SUCCESS);
365 }
366