Actual source code: pgmres.c
1: /*
2: This file implements PGMRES (a Pipelined Generalized Minimal Residual method)
3: */
5: #include <../src/ksp/ksp/impls/gmres/pgmres/pgmresimpl.h>
7: static PetscErrorCode KSPPGMRESUpdateHessenberg(KSP, PetscInt, PetscBool *, PetscReal *);
8: static PetscErrorCode KSPPGMRESBuildSoln(PetscScalar *, Vec, Vec, KSP, PetscInt);
10: static PetscErrorCode KSPSetUp_PGMRES(KSP ksp)
11: {
12: PetscFunctionBegin;
13: PetscCall(KSPSetUp_GMRES(ksp));
14: PetscFunctionReturn(PETSC_SUCCESS);
15: }
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;
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;
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: }
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. */
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: }
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: }
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)));
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;
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: }
89: if (!(it < pgmres->max_k + 1 && ksp->its < ksp->max_it)) break;
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)));
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: }
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));
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: }
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)));
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. */
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));
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: }
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;
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));
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: }
192: static PetscErrorCode KSPDestroy_PGMRES(KSP ksp)
193: {
194: PetscFunctionBegin;
195: PetscCall(KSPDestroy_GMRES(ksp));
196: PetscFunctionReturn(PETSC_SUCCESS);
197: }
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;
205: PetscFunctionBegin;
206: /* Solve for solution vector that minimizes the residual */
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: }
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;
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: }
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) PetscCall(VecAXPY(vdest, 1.0, VEC_TEMP));
229: else PetscCall(VecWAXPY(vdest, 1.0, VEC_TEMP, vguess));
230: PetscFunctionReturn(PETSC_SUCCESS);
231: }
233: static PetscErrorCode KSPPGMRESUpdateHessenberg(KSP ksp, PetscInt it, PetscBool *hapend, PetscReal *res)
234: {
235: PetscScalar *hh, *cc, *ss, *rs;
236: PetscInt j;
237: PetscReal hapbnd;
238: KSP_PGMRES *pgmres = (KSP_PGMRES *)ksp->data;
240: PetscFunctionBegin;
241: hh = HH(0, it); /* pointer to beginning of column to update */
242: cc = CC(0); /* beginning of cosine rotations */
243: ss = SS(0); /* beginning of sine rotations */
244: rs = RS(0); /* right-hand side of least squares system */
246: /* The Hessenberg matrix is now correct through column it, save that form for possible spectral analysis */
247: for (j = 0; j <= it + 1; j++) *HES(j, it) = hh[j];
249: /* check for the happy breakdown */
250: hapbnd = PetscMin(PetscAbsScalar(hh[it + 1] / rs[it]), pgmres->haptol);
251: if (PetscAbsScalar(hh[it + 1]) < hapbnd) {
252: 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))));
253: *hapend = PETSC_TRUE;
254: }
256: /* Apply all the previously computed plane rotations to the new column of the Hessenberg matrix */
257: /* Note: this uses the rotation [conj(c) s ; -s c], c= cos(theta), s= sin(theta),
258: and some refs have [c s ; -conj(s) c] (don't be confused!) */
260: for (j = 0; j < it; j++) {
261: PetscScalar hhj = hh[j];
262: hh[j] = PetscConj(cc[j]) * hhj + ss[j] * hh[j + 1];
263: hh[j + 1] = -ss[j] * hhj + cc[j] * hh[j + 1];
264: }
266: /*
267: compute the new plane rotation, and apply it to:
268: 1) the right-hand side of the Hessenberg system (RS)
269: note: it affects RS(it) and RS(it+1)
270: 2) the new column of the Hessenberg matrix
271: note: it affects HH(it,it) which is currently pointed to
272: by hh and HH(it+1, it) (*(hh+1))
273: thus obtaining the updated value of the residual...
274: */
276: /* compute new plane rotation */
278: if (!*hapend) {
279: PetscReal delta = PetscSqrtReal(PetscSqr(PetscAbsScalar(hh[it])) + PetscSqr(PetscAbsScalar(hh[it + 1])));
280: if (delta == 0.0) {
281: ksp->reason = KSP_DIVERGED_NULL;
282: PetscFunctionReturn(PETSC_SUCCESS);
283: }
285: cc[it] = hh[it] / delta; /* new cosine value */
286: ss[it] = hh[it + 1] / delta; /* new sine value */
288: hh[it] = PetscConj(cc[it]) * hh[it] + ss[it] * hh[it + 1];
289: rs[it + 1] = -ss[it] * rs[it];
290: rs[it] = PetscConj(cc[it]) * rs[it];
291: *res = PetscAbsScalar(rs[it + 1]);
292: } else { /* happy breakdown: HH(it+1, it) = 0, therefore we don't need to apply
293: another rotation matrix (so RH doesn't change). The new residual is
294: always the new sine term times the residual from last time (RS(it)),
295: but now the new sine rotation would be zero...so the residual should
296: be zero...so we will multiply "zero" by the last residual. This might
297: not be exactly what we want to do here -could just return "zero". */
298: *res = 0.0;
299: }
300: PetscFunctionReturn(PETSC_SUCCESS);
301: }
303: static PetscErrorCode KSPBuildSolution_PGMRES(KSP ksp, Vec ptr, Vec *result)
304: {
305: KSP_PGMRES *pgmres = (KSP_PGMRES *)ksp->data;
307: PetscFunctionBegin;
308: if (!ptr) {
309: if (!pgmres->sol_temp) PetscCall(VecDuplicate(ksp->vec_sol, &pgmres->sol_temp));
310: ptr = pgmres->sol_temp;
311: }
312: if (!pgmres->nrs) {
313: /* allocate the work area */
314: PetscCall(PetscMalloc1(pgmres->max_k, &pgmres->nrs));
315: }
317: PetscCall(KSPPGMRESBuildSoln(pgmres->nrs, ksp->vec_sol, ptr, ksp, pgmres->it));
318: if (result) *result = ptr;
319: PetscFunctionReturn(PETSC_SUCCESS);
320: }
322: static PetscErrorCode KSPSetFromOptions_PGMRES(KSP ksp, PetscOptionItems PetscOptionsObject)
323: {
324: PetscFunctionBegin;
325: PetscCall(KSPSetFromOptions_GMRES(ksp, PetscOptionsObject));
326: PetscOptionsHeadBegin(PetscOptionsObject, "KSP pipelined GMRES Options");
327: PetscOptionsHeadEnd();
328: PetscFunctionReturn(PETSC_SUCCESS);
329: }
331: static PetscErrorCode KSPReset_PGMRES(KSP ksp)
332: {
333: PetscFunctionBegin;
334: PetscCall(KSPReset_GMRES(ksp));
335: PetscFunctionReturn(PETSC_SUCCESS);
336: }
338: /*MC
339: KSPPGMRES - Implements the Pipelined Generalized Minimal Residual method {cite}`ghyselsashbymeerbergenvanroose2013`. [](sec_pipelineksp)
341: Options Database Keys:
342: + -ksp_gmres_restart restart - the number of Krylov directions to orthogonalize against
343: . -ksp_gmres_haptol tol - sets the tolerance for "happy breakdown" (exact convergence)
344: . -ksp_gmres_preallocate (true|false) - preallocate all the Krylov search directions initially
345: (otherwise groups of vectors are allocated as needed)
346: . -ksp_gmres_classicalgramschmidt (true|false) - use classical (unmodified) Gram-Schmidt to orthogonalize
347: against the Krylov space (fast) (the default)
348: . -ksp_gmres_modifiedgramschmidt (true|false) - use modified Gram-Schmidt in the orthogonalization (more stable, but slower)
349: . -ksp_gmres_cgs_refinement_type (refine_never|refine_ifneeded|refine_always) - determine if iterative refinement is used to increase the
350: stability of the classical Gram-Schmidt orthogonalization.
351: - -ksp_gmres_krylov_monitor - plot the Krylov space generated
353: Level: beginner
355: Note:
356: MPI configuration may be necessary for reductions to make asynchronous progress, which is important for performance of pipelined methods.
357: See [](doc_faq_pipelined)
359: Developer Note:
360: 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
362: .seealso: [](ch_ksp), [](sec_pipelineksp), [](doc_faq_pipelined), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPGMRES`, `KSPLGMRES`, `KSPPIPECG`, `KSPPIPECR`,
363: `KSPGMRESSetRestart()`, `KSPGMRESSetHapTol()`, `KSPGMRESSetPreAllocateVectors()`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESGetOrthogonalization()`,
364: `KSPGMRESClassicalGramSchmidtOrthogonalization()`, `KSPGMRESModifiedGramSchmidtOrthogonalization()`,
365: `KSPGMRESCGSRefinementType`, `KSPGMRESSetCGSRefinementType()`, `KSPGMRESGetCGSRefinementType()`, `KSPGMRESMonitorKrylov()`
366: M*/
368: PETSC_EXTERN PetscErrorCode KSPCreate_PGMRES(KSP ksp)
369: {
370: KSP_PGMRES *pgmres;
372: PetscFunctionBegin;
373: PetscCall(PetscNew(&pgmres));
375: ksp->data = (void *)pgmres;
376: ksp->ops->buildsolution = KSPBuildSolution_PGMRES;
377: ksp->ops->setup = KSPSetUp_PGMRES;
378: ksp->ops->solve = KSPSolve_PGMRES;
379: ksp->ops->reset = KSPReset_PGMRES;
380: ksp->ops->destroy = KSPDestroy_PGMRES;
381: ksp->ops->view = KSPView_GMRES;
382: ksp->ops->setfromoptions = KSPSetFromOptions_PGMRES;
383: ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_GMRES;
384: ksp->ops->computeeigenvalues = KSPComputeEigenvalues_GMRES;
386: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
387: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2));
388: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));
390: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetPreAllocateVectors_C", KSPGMRESSetPreAllocateVectors_GMRES));
391: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetOrthogonalization_C", KSPGMRESSetOrthogonalization_GMRES));
392: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetOrthogonalization_C", KSPGMRESGetOrthogonalization_GMRES));
393: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetRestart_C", KSPGMRESSetRestart_GMRES));
394: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetRestart_C", KSPGMRESGetRestart_GMRES));
395: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetCGSRefinementType_C", KSPGMRESSetCGSRefinementType_GMRES));
396: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetCGSRefinementType_C", KSPGMRESGetCGSRefinementType_GMRES));
398: pgmres->nextra_vecs = 1;
399: pgmres->haptol = 1.0e-30;
400: pgmres->q_preallocate = PETSC_FALSE;
401: pgmres->delta_allocate = PGMRES_DELTA_DIRECTIONS;
402: pgmres->orthog = KSPGMRESClassicalGramSchmidtOrthogonalization;
403: pgmres->nrs = NULL;
404: pgmres->sol_temp = NULL;
405: pgmres->max_k = PGMRES_DEFAULT_MAXK;
406: pgmres->Rsvd = NULL;
407: pgmres->orthogwork = NULL;
408: pgmres->cgstype = KSP_GMRES_CGS_REFINE_NEVER;
409: PetscFunctionReturn(PETSC_SUCCESS);
410: }