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 (it < pgmres->max_k + 1 || ksp->reason || ksp->its == ksp->max_it) { /* Monitor if we are done or still iterating, but not before a restart. */
 78:         PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
 79:         PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm));
 80:       }
 81:       if (ksp->reason) break;
 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));
155:   PetscFunctionReturn(PETSC_SUCCESS);
156: }

158: static PetscErrorCode KSPSolve_PGMRES(KSP ksp)
159: {
160:   PetscInt    its, itcount;
161:   KSP_PGMRES *pgmres     = (KSP_PGMRES *)ksp->data;
162:   PetscBool   guess_zero = ksp->guess_zero;

164:   PetscFunctionBegin;
165:   PetscCheck(!ksp->calc_sings || pgmres->Rsvd, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ORDER, "Must call KSPSetComputeSingularValues() before KSPSetUp() is called");
166:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
167:   ksp->its = 0;
168:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));

170:   itcount     = 0;
171:   ksp->reason = KSP_CONVERGED_ITERATING;
172:   while (!ksp->reason) {
173:     PetscCall(KSPInitialResidual(ksp, ksp->vec_sol, VEC_TEMP, VEC_TEMP_MATOP, VEC_VV(0), ksp->vec_rhs));
174:     PetscCall(KSPPGMRESCycle(&its, ksp));
175:     itcount += its;
176:     if (itcount >= ksp->max_it) {
177:       if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
178:       break;
179:     }
180:     ksp->guess_zero = PETSC_FALSE; /* every future call to KSPInitialResidual() will have nonzero guess */
181:   }
182:   ksp->guess_zero = guess_zero; /* restore if user provided nonzero initial guess */
183:   PetscFunctionReturn(PETSC_SUCCESS);
184: }

186: static PetscErrorCode KSPDestroy_PGMRES(KSP ksp)
187: {
188:   PetscFunctionBegin;
189:   PetscCall(KSPDestroy_GMRES(ksp));
190:   PetscFunctionReturn(PETSC_SUCCESS);
191: }

193: static PetscErrorCode KSPPGMRESBuildSoln(PetscScalar *nrs, Vec vguess, Vec vdest, KSP ksp, PetscInt it)
194: {
195:   PetscScalar tt;
196:   PetscInt    k, j;
197:   KSP_PGMRES *pgmres = (KSP_PGMRES *)(ksp->data);

199:   PetscFunctionBegin;
200:   /* Solve for solution vector that minimizes the residual */

202:   if (it < 0) {                        /* no pgmres steps have been performed */
203:     PetscCall(VecCopy(vguess, vdest)); /* VecCopy() is smart, exits immediately if vguess == vdest */
204:     PetscFunctionReturn(PETSC_SUCCESS);
205:   }

207:   /* solve the upper triangular system - RS is the right side and HH is
208:      the upper triangular matrix  - put soln in nrs */
209:   if (*HH(it, it) != 0.0) nrs[it] = *RS(it) / *HH(it, it);
210:   else nrs[it] = 0.0;

212:   for (k = it - 1; k >= 0; k--) {
213:     tt = *RS(k);
214:     for (j = k + 1; j <= it; j++) tt -= *HH(k, j) * nrs[j];
215:     nrs[k] = tt / *HH(k, k);
216:   }

218:   /* Accumulate the correction to the solution of the preconditioned problem in TEMP */
219:   PetscCall(VecMAXPBY(VEC_TEMP, it + 1, nrs, 0, &VEC_VV(0)));
220:   PetscCall(KSPUnwindPreconditioner(ksp, VEC_TEMP, VEC_TEMP_MATOP));
221:   /* add solution to previous solution */
222:   if (vdest == vguess) {
223:     PetscCall(VecAXPY(vdest, 1.0, VEC_TEMP));
224:   } else {
225:     PetscCall(VecWAXPY(vdest, 1.0, VEC_TEMP, vguess));
226:   }
227:   PetscFunctionReturn(PETSC_SUCCESS);
228: }

230: static PetscErrorCode KSPPGMRESUpdateHessenberg(KSP ksp, PetscInt it, PetscBool *hapend, PetscReal *res)
231: {
232:   PetscScalar *hh, *cc, *ss, *rs;
233:   PetscInt     j;
234:   PetscReal    hapbnd;
235:   KSP_PGMRES  *pgmres = (KSP_PGMRES *)(ksp->data);

237:   PetscFunctionBegin;
238:   hh = HH(0, it); /* pointer to beginning of column to update */
239:   cc = CC(0);     /* beginning of cosine rotations */
240:   ss = SS(0);     /* beginning of sine rotations */
241:   rs = RS(0);     /* right hand side of least squares system */

243:   /* The Hessenberg matrix is now correct through column it, save that form for possible spectral analysis */
244:   for (j = 0; j <= it + 1; j++) *HES(j, it) = hh[j];

246:   /* check for the happy breakdown */
247:   hapbnd = PetscMin(PetscAbsScalar(hh[it + 1] / rs[it]), pgmres->haptol);
248:   if (PetscAbsScalar(hh[it + 1]) < hapbnd) {
249:     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))));
250:     *hapend = PETSC_TRUE;
251:   }

253:   /* Apply all the previously computed plane rotations to the new column of the Hessenberg matrix */
254:   /* Note: this uses the rotation [conj(c)  s ; -s   c], c= cos(theta), s= sin(theta),
255:      and some refs have [c   s ; -conj(s)  c] (don't be confused!) */

257:   for (j = 0; j < it; j++) {
258:     PetscScalar hhj = hh[j];
259:     hh[j]           = PetscConj(cc[j]) * hhj + ss[j] * hh[j + 1];
260:     hh[j + 1]       = -ss[j] * hhj + cc[j] * hh[j + 1];
261:   }

263:   /*
264:     compute the new plane rotation, and apply it to:
265:      1) the right-hand-side of the Hessenberg system (RS)
266:         note: it affects RS(it) and RS(it+1)
267:      2) the new column of the Hessenberg matrix
268:         note: it affects HH(it,it) which is currently pointed to
269:         by hh and HH(it+1, it) (*(hh+1))
270:     thus obtaining the updated value of the residual...
271:   */

273:   /* compute new plane rotation */

275:   if (!*hapend) {
276:     PetscReal delta = PetscSqrtReal(PetscSqr(PetscAbsScalar(hh[it])) + PetscSqr(PetscAbsScalar(hh[it + 1])));
277:     if (delta == 0.0) {
278:       ksp->reason = KSP_DIVERGED_NULL;
279:       PetscFunctionReturn(PETSC_SUCCESS);
280:     }

282:     cc[it] = hh[it] / delta;     /* new cosine value */
283:     ss[it] = hh[it + 1] / delta; /* new sine value */

285:     hh[it]     = PetscConj(cc[it]) * hh[it] + ss[it] * hh[it + 1];
286:     rs[it + 1] = -ss[it] * rs[it];
287:     rs[it]     = PetscConj(cc[it]) * rs[it];
288:     *res       = PetscAbsScalar(rs[it + 1]);
289:   } else { /* happy breakdown: HH(it+1, it) = 0, therefore we don't need to apply
290:             another rotation matrix (so RH doesn't change).  The new residual is
291:             always the new sine term times the residual from last time (RS(it)),
292:             but now the new sine rotation would be zero...so the residual should
293:             be zero...so we will multiply "zero" by the last residual.  This might
294:             not be exactly what we want to do here -could just return "zero". */
295:     *res = 0.0;
296:   }
297:   PetscFunctionReturn(PETSC_SUCCESS);
298: }

300: static PetscErrorCode KSPBuildSolution_PGMRES(KSP ksp, Vec ptr, Vec *result)
301: {
302:   KSP_PGMRES *pgmres = (KSP_PGMRES *)ksp->data;

304:   PetscFunctionBegin;
305:   if (!ptr) {
306:     if (!pgmres->sol_temp) PetscCall(VecDuplicate(ksp->vec_sol, &pgmres->sol_temp));
307:     ptr = pgmres->sol_temp;
308:   }
309:   if (!pgmres->nrs) {
310:     /* allocate the work area */
311:     PetscCall(PetscMalloc1(pgmres->max_k, &pgmres->nrs));
312:   }

314:   PetscCall(KSPPGMRESBuildSoln(pgmres->nrs, ksp->vec_sol, ptr, ksp, pgmres->it));
315:   if (result) *result = ptr;
316:   PetscFunctionReturn(PETSC_SUCCESS);
317: }

319: static PetscErrorCode KSPSetFromOptions_PGMRES(KSP ksp, PetscOptionItems *PetscOptionsObject)
320: {
321:   PetscFunctionBegin;
322:   PetscCall(KSPSetFromOptions_GMRES(ksp, PetscOptionsObject));
323:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP pipelined GMRES Options");
324:   PetscOptionsHeadEnd();
325:   PetscFunctionReturn(PETSC_SUCCESS);
326: }

328: static PetscErrorCode KSPReset_PGMRES(KSP ksp)
329: {
330:   PetscFunctionBegin;
331:   PetscCall(KSPReset_GMRES(ksp));
332:   PetscFunctionReturn(PETSC_SUCCESS);
333: }

335: /*MC
336:    KSPPGMRES - Implements the Pipelined Generalized Minimal Residual method {cite}`ghyselsashbymeerbergenvanroose2013`. [](sec_pipelineksp)

338:    Options Database Keys:
339: +   -ksp_gmres_restart <restart>                                                - the number of Krylov directions to orthogonalize against
340: .   -ksp_gmres_haptol <tol>                                                     - sets the tolerance for "happy ending" (exact convergence)
341: .   -ksp_gmres_preallocate                                                      - preallocate all the Krylov search directions initially
342:                                                                                 (otherwise groups of vectors are allocated as needed)
343: .   -ksp_gmres_classicalgramschmidt                                             - use classical (unmodified) Gram-Schmidt to orthogonalize
344:                                                                                 against the Krylov space (fast) (the default)
345: .   -ksp_gmres_modifiedgramschmidt                                              - use modified Gram-Schmidt in the orthogonalization (more stable, but slower)
346: .   -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always> - determine if iterative refinement is used to increase the
347:                                                                                 stability of the classical Gram-Schmidt  orthogonalization.
348: -   -ksp_gmres_krylov_monitor                                                   - plot the Krylov space generated

350:    Level: beginner

352:    Note:
353:    MPI configuration may be necessary for reductions to make asynchronous progress, which is important for performance of pipelined methods.
354:    See [](doc_faq_pipelined)

356:    Developer Note:
357:    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

359: .seealso: [](ch_ksp), [](sec_pipelineksp), [](doc_faq_pipelined), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPGMRES`, `KSPLGMRES`, `KSPPIPECG`, `KSPPIPECR`,
360:           `KSPGMRESSetRestart()`, `KSPGMRESSetHapTol()`, `KSPGMRESSetPreAllocateVectors()`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESGetOrthogonalization()`,
361:           `KSPGMRESClassicalGramSchmidtOrthogonalization()`, `KSPGMRESModifiedGramSchmidtOrthogonalization()`,
362:           `KSPGMRESCGSRefinementType`, `KSPGMRESSetCGSRefinementType()`, `KSPGMRESGetCGSRefinementType()`, `KSPGMRESMonitorKrylov()`
363: M*/

365: PETSC_EXTERN PetscErrorCode KSPCreate_PGMRES(KSP ksp)
366: {
367:   KSP_PGMRES *pgmres;

369:   PetscFunctionBegin;
370:   PetscCall(PetscNew(&pgmres));

372:   ksp->data                              = (void *)pgmres;
373:   ksp->ops->buildsolution                = KSPBuildSolution_PGMRES;
374:   ksp->ops->setup                        = KSPSetUp_PGMRES;
375:   ksp->ops->solve                        = KSPSolve_PGMRES;
376:   ksp->ops->reset                        = KSPReset_PGMRES;
377:   ksp->ops->destroy                      = KSPDestroy_PGMRES;
378:   ksp->ops->view                         = KSPView_GMRES;
379:   ksp->ops->setfromoptions               = KSPSetFromOptions_PGMRES;
380:   ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_GMRES;
381:   ksp->ops->computeeigenvalues           = KSPComputeEigenvalues_GMRES;

383:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
384:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2));
385:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));

387:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetPreAllocateVectors_C", KSPGMRESSetPreAllocateVectors_GMRES));
388:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetOrthogonalization_C", KSPGMRESSetOrthogonalization_GMRES));
389:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetOrthogonalization_C", KSPGMRESGetOrthogonalization_GMRES));
390:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetRestart_C", KSPGMRESSetRestart_GMRES));
391:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetRestart_C", KSPGMRESGetRestart_GMRES));
392:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetCGSRefinementType_C", KSPGMRESSetCGSRefinementType_GMRES));
393:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetCGSRefinementType_C", KSPGMRESGetCGSRefinementType_GMRES));

395:   pgmres->nextra_vecs    = 1;
396:   pgmres->haptol         = 1.0e-30;
397:   pgmres->q_preallocate  = 0;
398:   pgmres->delta_allocate = PGMRES_DELTA_DIRECTIONS;
399:   pgmres->orthog         = KSPGMRESClassicalGramSchmidtOrthogonalization;
400:   pgmres->nrs            = NULL;
401:   pgmres->sol_temp       = NULL;
402:   pgmres->max_k          = PGMRES_DEFAULT_MAXK;
403:   pgmres->Rsvd           = NULL;
404:   pgmres->orthogwork     = NULL;
405:   pgmres->cgstype        = KSP_GMRES_CGS_REFINE_NEVER;
406:   PetscFunctionReturn(PETSC_SUCCESS);
407: }