1 /*
2 This file implements PGMRES (a Pipelined Generalized Minimal Residual method)
3 */
4
5 #include <../src/ksp/ksp/impls/gmres/pgmres/pgmresimpl.h> /*I "petscksp.h" I*/
6
7 static PetscErrorCode KSPPGMRESUpdateHessenberg(KSP, PetscInt, PetscBool *, PetscReal *);
8 static PetscErrorCode KSPPGMRESBuildSoln(PetscScalar *, Vec, Vec, KSP, PetscInt);
9
KSPSetUp_PGMRES(KSP ksp)10 static PetscErrorCode KSPSetUp_PGMRES(KSP ksp)
11 {
12 PetscFunctionBegin;
13 PetscCall(KSPSetUp_GMRES(ksp));
14 PetscFunctionReturn(PETSC_SUCCESS);
15 }
16
KSPPGMRESCycle(PetscInt * itcount,KSP ksp)17 static PetscErrorCode KSPPGMRESCycle(PetscInt *itcount, KSP ksp)
18 {
19 KSP_PGMRES *pgmres = (KSP_PGMRES *)ksp->data;
20 PetscReal res_norm, res, newnorm;
21 PetscInt it = 0, j, k;
22 PetscBool hapend = PETSC_FALSE;
23
24 PetscFunctionBegin;
25 if (itcount) *itcount = 0;
26 PetscCall(VecNormalize(VEC_VV(0), &res_norm));
27 KSPCheckNorm(ksp, res_norm);
28 res = res_norm;
29 *RS(0) = res_norm;
30
31 /* check for the convergence */
32 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
33 if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = res;
34 else ksp->rnorm = 0;
35 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
36 pgmres->it = it - 2;
37 PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
38 PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm));
39 if (!res) {
40 ksp->reason = KSP_CONVERGED_ATOL;
41 PetscCall(PetscInfo(ksp, "Converged due to zero residual norm on entry\n"));
42 PetscFunctionReturn(PETSC_SUCCESS);
43 }
44
45 PetscCall((*ksp->converged)(ksp, ksp->its, ksp->rnorm, &ksp->reason, ksp->cnvP));
46 for (; !ksp->reason; it++) {
47 Vec Zcur, Znext;
48 if (pgmres->vv_allocated <= it + VEC_OFFSET + 1) PetscCall(KSPGMRESGetNewVectors(ksp, it + 1));
49 /* VEC_VV(it-1) is orthogonal, it will be normalized once the VecNorm arrives. */
50 Zcur = VEC_VV(it); /* Zcur is not yet orthogonal, but the VecMDot to orthogonalize it has been started. */
51 Znext = VEC_VV(it + 1); /* This iteration will compute Znext, update with a deferred correction once we know how
52 Zcur relates to the previous vectors, and start the reduction to orthogonalize it. */
53
54 if (it < pgmres->max_k + 1 && ksp->its + 1 < PetscMax(2, ksp->max_it)) { /* We don't know whether what we have computed is enough, so apply the matrix. */
55 PetscCall(KSP_PCApplyBAorAB(ksp, Zcur, Znext, VEC_TEMP_MATOP));
56 }
57
58 if (it > 1) { /* Complete the pending reduction */
59 PetscCall(VecNormEnd(VEC_VV(it - 1), NORM_2, &newnorm));
60 *HH(it - 1, it - 2) = newnorm;
61 }
62 if (it > 0) { /* Finish the reduction computing the latest column of H */
63 PetscCall(VecMDotEnd(Zcur, it, &(VEC_VV(0)), HH(0, it - 1)));
64 }
65
66 if (it > 1) {
67 /* normalize the base vector from two iterations ago, basis is complete up to here */
68 PetscCall(VecScale(VEC_VV(it - 1), 1. / *HH(it - 1, it - 2)));
69
70 PetscCall(KSPPGMRESUpdateHessenberg(ksp, it - 2, &hapend, &res));
71 pgmres->it = it - 2;
72 ksp->its++;
73 if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = res;
74 else ksp->rnorm = 0;
75
76 PetscCall((*ksp->converged)(ksp, ksp->its, ksp->rnorm, &ksp->reason, ksp->cnvP));
77 if (ksp->reason) break;
78 if (it < pgmres->max_k + 1) { /* Monitor if we are not done or still iterating, but not before a restart. */
79 PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
80 PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm));
81 }
82 /* Catch error in happy breakdown and signal convergence and break from loop */
83 if (hapend) {
84 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "Reached happy break down, but convergence was not indicated. Residual norm = %g", (double)res);
85 ksp->reason = KSP_DIVERGED_BREAKDOWN;
86 break;
87 }
88
89 if (!(it < pgmres->max_k + 1 && ksp->its < ksp->max_it)) break;
90
91 /* The it-2 column of H was not scaled when we computed Zcur, apply correction */
92 PetscCall(VecScale(Zcur, 1. / *HH(it - 1, it - 2)));
93 /* And Znext computed in this iteration was computed using the under-scaled Zcur */
94 PetscCall(VecScale(Znext, 1. / *HH(it - 1, it - 2)));
95
96 /* In the previous iteration, we projected an unnormalized Zcur against the Krylov basis, so we need to fix the column of H resulting from that projection. */
97 for (k = 0; k < it; k++) *HH(k, it - 1) /= *HH(it - 1, it - 2);
98 /* When Zcur was projected against the Krylov basis, VV(it-1) was still not normalized, so fix that too. This
99 * column is complete except for HH(it,it-1) which we won't know until the next iteration. */
100 *HH(it - 1, it - 1) /= *HH(it - 1, it - 2);
101 }
102
103 if (it > 0) {
104 PetscScalar *work;
105 if (!pgmres->orthogwork) PetscCall(PetscMalloc1(pgmres->max_k + 2, &pgmres->orthogwork));
106 work = pgmres->orthogwork;
107 /* Apply correction computed by the VecMDot in the last iteration to Znext. The original form is
108 *
109 * Znext -= sum_{j=0}^{i-1} Z[j+1] * H[j,i-1]
110 *
111 * where
112 *
113 * Z[j] = sum_{k=0}^j V[k] * H[k,j-1]
114 *
115 * substituting
116 *
117 * Znext -= sum_{j=0}^{i-1} sum_{k=0}^{j+1} V[k] * H[k,j] * H[j,i-1]
118 *
119 * rearranging the iteration space from row-column to column-row
120 *
121 * Znext -= sum_{k=0}^i sum_{j=k-1}^{i-1} V[k] * H[k,j] * H[j,i-1]
122 *
123 * Note that column it-1 of HH is correct. For all previous columns, we must look at HES because HH has already
124 * been transformed to upper triangular form.
125 */
126 for (k = 0; k < it + 1; k++) {
127 work[k] = 0;
128 for (j = PetscMax(0, k - 1); j < it - 1; j++) work[k] -= *HES(k, j) * *HH(j, it - 1);
129 }
130 PetscCall(VecMAXPY(Znext, it + 1, work, &VEC_VV(0)));
131 PetscCall(VecAXPY(Znext, -*HH(it - 1, it - 1), Zcur));
132
133 /* Orthogonalize Zcur against existing basis vectors. */
134 for (k = 0; k < it; k++) work[k] = -*HH(k, it - 1);
135 PetscCall(VecMAXPY(Zcur, it, work, &VEC_VV(0)));
136 /* Zcur is now orthogonal, and will be referred to as VEC_VV(it) again, though it is still not normalized. */
137 /* Begin computing the norm of the new vector, will be normalized after the MatMult in the next iteration. */
138 PetscCall(VecNormBegin(VEC_VV(it), NORM_2, &newnorm));
139 }
140
141 /* Compute column of H (to the diagonal, but not the subdiagonal) to be able to orthogonalize the newest vector. */
142 PetscCall(VecMDotBegin(Znext, it + 1, &VEC_VV(0), HH(0, it)));
143
144 /* Start an asynchronous split-mode reduction, the result of the MDot and Norm will be collected on the next iteration. */
145 PetscCall(PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)Znext)));
146 }
147 if (itcount) *itcount = it - 1; /* Number of iterations actually completed. */
148
149 /*
150 Solve for the "best" coefficients of the Krylov
151 columns, add the solution values together, and possibly unwind the preconditioning from the solution
152 */
153 /* Form the solution (or the solution so far) */
154 PetscCall(KSPPGMRESBuildSoln(RS(0), ksp->vec_sol, ksp->vec_sol, ksp, it - 2));
155
156 if (ksp->reason == KSP_CONVERGED_ITERATING && ksp->its == ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
157 if (ksp->reason) {
158 PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
159 PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm));
160 }
161 PetscFunctionReturn(PETSC_SUCCESS);
162 }
163
KSPSolve_PGMRES(KSP ksp)164 static PetscErrorCode KSPSolve_PGMRES(KSP ksp)
165 {
166 PetscInt its, itcount;
167 KSP_PGMRES *pgmres = (KSP_PGMRES *)ksp->data;
168 PetscBool guess_zero = ksp->guess_zero;
169
170 PetscFunctionBegin;
171 PetscCheck(!ksp->calc_sings || pgmres->Rsvd, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ORDER, "Must call KSPSetComputeSingularValues() before KSPSetUp() is called");
172 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
173 ksp->its = 0;
174 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
175
176 itcount = 0;
177 ksp->reason = KSP_CONVERGED_ITERATING;
178 while (!ksp->reason) {
179 PetscCall(KSPInitialResidual(ksp, ksp->vec_sol, VEC_TEMP, VEC_TEMP_MATOP, VEC_VV(0), ksp->vec_rhs));
180 PetscCall(KSPPGMRESCycle(&its, ksp));
181 itcount += its;
182 if (itcount >= ksp->max_it) {
183 if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
184 break;
185 }
186 ksp->guess_zero = PETSC_FALSE; /* every future call to KSPInitialResidual() will have nonzero guess */
187 }
188 ksp->guess_zero = guess_zero; /* restore if user provided nonzero initial guess */
189 PetscFunctionReturn(PETSC_SUCCESS);
190 }
191
KSPDestroy_PGMRES(KSP ksp)192 static PetscErrorCode KSPDestroy_PGMRES(KSP ksp)
193 {
194 PetscFunctionBegin;
195 PetscCall(KSPDestroy_GMRES(ksp));
196 PetscFunctionReturn(PETSC_SUCCESS);
197 }
198
KSPPGMRESBuildSoln(PetscScalar * nrs,Vec vguess,Vec vdest,KSP ksp,PetscInt it)199 static PetscErrorCode KSPPGMRESBuildSoln(PetscScalar *nrs, Vec vguess, Vec vdest, KSP ksp, PetscInt it)
200 {
201 PetscScalar tt;
202 PetscInt k, j;
203 KSP_PGMRES *pgmres = (KSP_PGMRES *)ksp->data;
204
205 PetscFunctionBegin;
206 /* Solve for solution vector that minimizes the residual */
207
208 if (it < 0) { /* no pgmres steps have been performed */
209 PetscCall(VecCopy(vguess, vdest)); /* VecCopy() is smart, exits immediately if vguess == vdest */
210 PetscFunctionReturn(PETSC_SUCCESS);
211 }
212
213 /* solve the upper triangular system - RS is the right side and HH is
214 the upper triangular matrix - put soln in nrs */
215 if (*HH(it, it) != 0.0) nrs[it] = *RS(it) / *HH(it, it);
216 else nrs[it] = 0.0;
217
218 for (k = it - 1; k >= 0; k--) {
219 tt = *RS(k);
220 for (j = k + 1; j <= it; j++) tt -= *HH(k, j) * nrs[j];
221 nrs[k] = tt / *HH(k, k);
222 }
223
224 /* Accumulate the correction to the solution of the preconditioned problem in TEMP */
225 PetscCall(VecMAXPBY(VEC_TEMP, it + 1, nrs, 0, &VEC_VV(0)));
226 PetscCall(KSPUnwindPreconditioner(ksp, VEC_TEMP, VEC_TEMP_MATOP));
227 /* add solution to previous solution */
228 if (vdest == vguess) {
229 PetscCall(VecAXPY(vdest, 1.0, VEC_TEMP));
230 } else {
231 PetscCall(VecWAXPY(vdest, 1.0, VEC_TEMP, vguess));
232 }
233 PetscFunctionReturn(PETSC_SUCCESS);
234 }
235
KSPPGMRESUpdateHessenberg(KSP ksp,PetscInt it,PetscBool * hapend,PetscReal * res)236 static PetscErrorCode KSPPGMRESUpdateHessenberg(KSP ksp, PetscInt it, PetscBool *hapend, PetscReal *res)
237 {
238 PetscScalar *hh, *cc, *ss, *rs;
239 PetscInt j;
240 PetscReal hapbnd;
241 KSP_PGMRES *pgmres = (KSP_PGMRES *)ksp->data;
242
243 PetscFunctionBegin;
244 hh = HH(0, it); /* pointer to beginning of column to update */
245 cc = CC(0); /* beginning of cosine rotations */
246 ss = SS(0); /* beginning of sine rotations */
247 rs = RS(0); /* right-hand side of least squares system */
248
249 /* The Hessenberg matrix is now correct through column it, save that form for possible spectral analysis */
250 for (j = 0; j <= it + 1; j++) *HES(j, it) = hh[j];
251
252 /* check for the happy breakdown */
253 hapbnd = PetscMin(PetscAbsScalar(hh[it + 1] / rs[it]), pgmres->haptol);
254 if (PetscAbsScalar(hh[it + 1]) < hapbnd) {
255 PetscCall(PetscInfo(ksp, "Detected happy breakdown, current hapbnd = %14.12e H(%" PetscInt_FMT ",%" PetscInt_FMT ") = %14.12e\n", (double)hapbnd, it + 1, it, (double)PetscAbsScalar(*HH(it + 1, it))));
256 *hapend = PETSC_TRUE;
257 }
258
259 /* Apply all the previously computed plane rotations to the new column of the Hessenberg matrix */
260 /* Note: this uses the rotation [conj(c) s ; -s c], c= cos(theta), s= sin(theta),
261 and some refs have [c s ; -conj(s) c] (don't be confused!) */
262
263 for (j = 0; j < it; j++) {
264 PetscScalar hhj = hh[j];
265 hh[j] = PetscConj(cc[j]) * hhj + ss[j] * hh[j + 1];
266 hh[j + 1] = -ss[j] * hhj + cc[j] * hh[j + 1];
267 }
268
269 /*
270 compute the new plane rotation, and apply it to:
271 1) the right-hand side of the Hessenberg system (RS)
272 note: it affects RS(it) and RS(it+1)
273 2) the new column of the Hessenberg matrix
274 note: it affects HH(it,it) which is currently pointed to
275 by hh and HH(it+1, it) (*(hh+1))
276 thus obtaining the updated value of the residual...
277 */
278
279 /* compute new plane rotation */
280
281 if (!*hapend) {
282 PetscReal delta = PetscSqrtReal(PetscSqr(PetscAbsScalar(hh[it])) + PetscSqr(PetscAbsScalar(hh[it + 1])));
283 if (delta == 0.0) {
284 ksp->reason = KSP_DIVERGED_NULL;
285 PetscFunctionReturn(PETSC_SUCCESS);
286 }
287
288 cc[it] = hh[it] / delta; /* new cosine value */
289 ss[it] = hh[it + 1] / delta; /* new sine value */
290
291 hh[it] = PetscConj(cc[it]) * hh[it] + ss[it] * hh[it + 1];
292 rs[it + 1] = -ss[it] * rs[it];
293 rs[it] = PetscConj(cc[it]) * rs[it];
294 *res = PetscAbsScalar(rs[it + 1]);
295 } else { /* happy breakdown: HH(it+1, it) = 0, therefore we don't need to apply
296 another rotation matrix (so RH doesn't change). The new residual is
297 always the new sine term times the residual from last time (RS(it)),
298 but now the new sine rotation would be zero...so the residual should
299 be zero...so we will multiply "zero" by the last residual. This might
300 not be exactly what we want to do here -could just return "zero". */
301 *res = 0.0;
302 }
303 PetscFunctionReturn(PETSC_SUCCESS);
304 }
305
KSPBuildSolution_PGMRES(KSP ksp,Vec ptr,Vec * result)306 static PetscErrorCode KSPBuildSolution_PGMRES(KSP ksp, Vec ptr, Vec *result)
307 {
308 KSP_PGMRES *pgmres = (KSP_PGMRES *)ksp->data;
309
310 PetscFunctionBegin;
311 if (!ptr) {
312 if (!pgmres->sol_temp) PetscCall(VecDuplicate(ksp->vec_sol, &pgmres->sol_temp));
313 ptr = pgmres->sol_temp;
314 }
315 if (!pgmres->nrs) {
316 /* allocate the work area */
317 PetscCall(PetscMalloc1(pgmres->max_k, &pgmres->nrs));
318 }
319
320 PetscCall(KSPPGMRESBuildSoln(pgmres->nrs, ksp->vec_sol, ptr, ksp, pgmres->it));
321 if (result) *result = ptr;
322 PetscFunctionReturn(PETSC_SUCCESS);
323 }
324
KSPSetFromOptions_PGMRES(KSP ksp,PetscOptionItems PetscOptionsObject)325 static PetscErrorCode KSPSetFromOptions_PGMRES(KSP ksp, PetscOptionItems PetscOptionsObject)
326 {
327 PetscFunctionBegin;
328 PetscCall(KSPSetFromOptions_GMRES(ksp, PetscOptionsObject));
329 PetscOptionsHeadBegin(PetscOptionsObject, "KSP pipelined GMRES Options");
330 PetscOptionsHeadEnd();
331 PetscFunctionReturn(PETSC_SUCCESS);
332 }
333
KSPReset_PGMRES(KSP ksp)334 static PetscErrorCode KSPReset_PGMRES(KSP ksp)
335 {
336 PetscFunctionBegin;
337 PetscCall(KSPReset_GMRES(ksp));
338 PetscFunctionReturn(PETSC_SUCCESS);
339 }
340
341 /*MC
342 KSPPGMRES - Implements the Pipelined Generalized Minimal Residual method {cite}`ghyselsashbymeerbergenvanroose2013`. [](sec_pipelineksp)
343
344 Options Database Keys:
345 + -ksp_gmres_restart <restart> - the number of Krylov directions to orthogonalize against
346 . -ksp_gmres_haptol <tol> - sets the tolerance for "happy ending" (exact convergence)
347 . -ksp_gmres_preallocate - preallocate all the Krylov search directions initially
348 (otherwise groups of vectors are allocated as needed)
349 . -ksp_gmres_classicalgramschmidt - use classical (unmodified) Gram-Schmidt to orthogonalize
350 against the Krylov space (fast) (the default)
351 . -ksp_gmres_modifiedgramschmidt - use modified Gram-Schmidt in the orthogonalization (more stable, but slower)
352 . -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always> - determine if iterative refinement is used to increase the
353 stability of the classical Gram-Schmidt orthogonalization.
354 - -ksp_gmres_krylov_monitor - plot the Krylov space generated
355
356 Level: beginner
357
358 Note:
359 MPI configuration may be necessary for reductions to make asynchronous progress, which is important for performance of pipelined methods.
360 See [](doc_faq_pipelined)
361
362 Developer Note:
363 This object is subclassed off of `KSPGMRES`, see the source code in src/ksp/ksp/impls/gmres for comments on the structure of the code
364
365 .seealso: [](ch_ksp), [](sec_pipelineksp), [](doc_faq_pipelined), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPGMRES`, `KSPLGMRES`, `KSPPIPECG`, `KSPPIPECR`,
366 `KSPGMRESSetRestart()`, `KSPGMRESSetHapTol()`, `KSPGMRESSetPreAllocateVectors()`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESGetOrthogonalization()`,
367 `KSPGMRESClassicalGramSchmidtOrthogonalization()`, `KSPGMRESModifiedGramSchmidtOrthogonalization()`,
368 `KSPGMRESCGSRefinementType`, `KSPGMRESSetCGSRefinementType()`, `KSPGMRESGetCGSRefinementType()`, `KSPGMRESMonitorKrylov()`
369 M*/
370
KSPCreate_PGMRES(KSP ksp)371 PETSC_EXTERN PetscErrorCode KSPCreate_PGMRES(KSP ksp)
372 {
373 KSP_PGMRES *pgmres;
374
375 PetscFunctionBegin;
376 PetscCall(PetscNew(&pgmres));
377
378 ksp->data = (void *)pgmres;
379 ksp->ops->buildsolution = KSPBuildSolution_PGMRES;
380 ksp->ops->setup = KSPSetUp_PGMRES;
381 ksp->ops->solve = KSPSolve_PGMRES;
382 ksp->ops->reset = KSPReset_PGMRES;
383 ksp->ops->destroy = KSPDestroy_PGMRES;
384 ksp->ops->view = KSPView_GMRES;
385 ksp->ops->setfromoptions = KSPSetFromOptions_PGMRES;
386 ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_GMRES;
387 ksp->ops->computeeigenvalues = KSPComputeEigenvalues_GMRES;
388
389 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
390 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2));
391 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));
392
393 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetPreAllocateVectors_C", KSPGMRESSetPreAllocateVectors_GMRES));
394 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetOrthogonalization_C", KSPGMRESSetOrthogonalization_GMRES));
395 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetOrthogonalization_C", KSPGMRESGetOrthogonalization_GMRES));
396 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetRestart_C", KSPGMRESSetRestart_GMRES));
397 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetRestart_C", KSPGMRESGetRestart_GMRES));
398 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetCGSRefinementType_C", KSPGMRESSetCGSRefinementType_GMRES));
399 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetCGSRefinementType_C", KSPGMRESGetCGSRefinementType_GMRES));
400
401 pgmres->nextra_vecs = 1;
402 pgmres->haptol = 1.0e-30;
403 pgmres->q_preallocate = PETSC_FALSE;
404 pgmres->delta_allocate = PGMRES_DELTA_DIRECTIONS;
405 pgmres->orthog = KSPGMRESClassicalGramSchmidtOrthogonalization;
406 pgmres->nrs = NULL;
407 pgmres->sol_temp = NULL;
408 pgmres->max_k = PGMRES_DEFAULT_MAXK;
409 pgmres->Rsvd = NULL;
410 pgmres->orthogwork = NULL;
411 pgmres->cgstype = KSP_GMRES_CGS_REFINE_NEVER;
412 PetscFunctionReturn(PETSC_SUCCESS);
413 }
414