Actual source code: gmres.c

```  1: /*
2:     This file implements GMRES (a Generalized Minimal Residual) method.
3:     Reference:  Saad and Schultz, 1986.

5:     Some comments on left vs. right preconditioning, and restarts.
6:     Left and right preconditioning.
7:     If right preconditioning is chosen, then the problem being solved
8:     by gmres is actually
9:        My =  AB^-1 y = f
10:     so the initial residual is
11:           r = f - Mx
12:     Note that B^-1 y = x or y = B x, and if x is non-zero, the initial
13:     residual is
14:           r = f - A x
15:     The final solution is then
16:           x = B^-1 y

18:     If left preconditioning is chosen, then the problem being solved is
19:        My = B^-1 A x = B^-1 f,
20:     and the initial residual is
21:        r  = B^-1(f - Ax)

23:     Restarts:  Restarts are basically solves with x0 not equal to zero.
24:     Note that we can eliminate an extra application of B^-1 between
25:     restarts as long as we don't require that the solution at the end
26:     of an unsuccessful gmres iteration always be the solution x.
27:  */

29: #include <../src/ksp/ksp/impls/gmres/gmresimpl.h>
30: #define GMRES_DELTA_DIRECTIONS 10
31: #define GMRES_DEFAULT_MAXK     30
32: static PetscErrorCode KSPGMRESUpdateHessenberg(KSP, PetscInt, PetscBool, PetscReal *);
33: static PetscErrorCode KSPGMRESBuildSoln(PetscScalar *, Vec, Vec, KSP, PetscInt);

35: PetscErrorCode KSPSetUp_GMRES(KSP ksp)
36: {
37:   PetscInt   hh, hes, rs, cc;
38:   PetscInt   max_k, k;
39:   KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;

41:   PetscFunctionBegin;
42:   max_k = gmres->max_k; /* restart size */
43:   hh    = (max_k + 2) * (max_k + 1);
44:   hes   = (max_k + 1) * (max_k + 1);
45:   rs    = (max_k + 2);
46:   cc    = (max_k + 1);

48:   PetscCall(PetscCalloc5(hh, &gmres->hh_origin, hes, &gmres->hes_origin, rs, &gmres->rs_origin, cc, &gmres->cc_origin, cc, &gmres->ss_origin));

50:   if (ksp->calc_sings) {
51:     /* Allocate workspace to hold Hessenberg matrix needed by lapack */
52:     PetscCall(PetscMalloc1((max_k + 3) * (max_k + 9), &gmres->Rsvd));
53:     PetscCall(PetscMalloc1(6 * (max_k + 2), &gmres->Dsvd));
54:   }

56:   /* Allocate array to hold pointers to user vectors.  Note that we need
57:    4 + max_k + 1 (since we need it+1 vectors, and it <= max_k) */
58:   gmres->vecs_allocated = VEC_OFFSET + 2 + max_k + gmres->nextra_vecs;

60:   PetscCall(PetscMalloc1(gmres->vecs_allocated, &gmres->vecs));
61:   PetscCall(PetscMalloc1(VEC_OFFSET + 2 + max_k, &gmres->user_work));
62:   PetscCall(PetscMalloc1(VEC_OFFSET + 2 + max_k, &gmres->mwork_alloc));

64:   if (gmres->q_preallocate) {
65:     gmres->vv_allocated = VEC_OFFSET + 2 + max_k;

67:     PetscCall(KSPCreateVecs(ksp, gmres->vv_allocated, &gmres->user_work[0], 0, NULL));

69:     gmres->mwork_alloc[0] = gmres->vv_allocated;
70:     gmres->nwork_alloc    = 1;
71:     for (k = 0; k < gmres->vv_allocated; k++) gmres->vecs[k] = gmres->user_work[0][k];
72:   } else {
73:     gmres->vv_allocated = 5;

75:     PetscCall(KSPCreateVecs(ksp, 5, &gmres->user_work[0], 0, NULL));

77:     gmres->mwork_alloc[0] = 5;
78:     gmres->nwork_alloc    = 1;
79:     for (k = 0; k < gmres->vv_allocated; k++) gmres->vecs[k] = gmres->user_work[0][k];
80:   }
81:   PetscFunctionReturn(PETSC_SUCCESS);
82: }

84: /*
85:     Run gmres, possibly with restart.  Return residual history if requested.
86:     input parameters:

88: .        gmres  - structure containing parameters and work areas

90:     output parameters:
91: .        nres    - residuals (from preconditioned system) at each step.
92:                   If restarting, consider passing nres+it.  If null,
93:                   ignored
94: .        itcount - number of iterations used.  nres[0] to nres[itcount]
95:                   are defined.  If null, ignored.

97:     Notes:
98:     On entry, the value in vector VEC_VV(0) should be the initial residual
99:     (this allows shortcuts where the initial preconditioned residual is 0).
100:  */
101: static PetscErrorCode KSPGMRESCycle(PetscInt *itcount, KSP ksp)
102: {
103:   KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;
104:   PetscReal  res, hapbnd, tt;
105:   PetscInt   it = 0, max_k = gmres->max_k;
106:   PetscBool  hapend = PETSC_FALSE;

108:   PetscFunctionBegin;
109:   if (itcount) *itcount = 0;
110:   PetscCall(VecNormalize(VEC_VV(0), &res));
111:   KSPCheckNorm(ksp, res);

113:   /* the constant .1 is arbitrary, just some measure at how incorrect the residuals are */
114:   if ((ksp->rnorm > 0.0) && (PetscAbsReal(res - ksp->rnorm) > gmres->breakdowntol * gmres->rnorm0)) {
115:     PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_CONV_FAILED, "Residual norm computed by GMRES recursion formula %g is far from the computed residual norm %g at restart, residual norm at start of cycle %g",
116:                (double)ksp->rnorm, (double)res, (double)gmres->rnorm0);
117:     PetscCall(PetscInfo(ksp, "Residual norm computed by GMRES recursion formula %g is far from the computed residual norm %g at restart, residual norm at start of cycle %g\n", (double)ksp->rnorm, (double)res, (double)gmres->rnorm0));
118:     ksp->reason = KSP_DIVERGED_BREAKDOWN;
119:     PetscFunctionReturn(PETSC_SUCCESS);
120:   }
121:   *GRS(0) = gmres->rnorm0 = res;

123:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
124:   ksp->rnorm = res;
125:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
126:   gmres->it = (it - 1);
127:   PetscCall(KSPLogResidualHistory(ksp, res));
128:   PetscCall(KSPLogErrorHistory(ksp));
129:   PetscCall(KSPMonitor(ksp, ksp->its, res));
130:   if (!res) {
131:     ksp->reason = KSP_CONVERGED_ATOL;
132:     PetscCall(PetscInfo(ksp, "Converged due to zero residual norm on entry\n"));
133:     PetscFunctionReturn(PETSC_SUCCESS);
134:   }

136:   /* check for the convergence */
137:   PetscCall((*ksp->converged)(ksp, ksp->its, res, &ksp->reason, ksp->cnvP));
138:   while (!ksp->reason && it < max_k && ksp->its < ksp->max_it) {
139:     if (it) {
140:       PetscCall(KSPLogResidualHistory(ksp, res));
141:       PetscCall(KSPLogErrorHistory(ksp));
142:       PetscCall(KSPMonitor(ksp, ksp->its, res));
143:     }
144:     gmres->it = (it - 1);
145:     if (gmres->vv_allocated <= it + VEC_OFFSET + 1) PetscCall(KSPGMRESGetNewVectors(ksp, it + 1));
146:     PetscCall(KSP_PCApplyBAorAB(ksp, VEC_VV(it), VEC_VV(1 + it), VEC_TEMP_MATOP));

148:     /* update Hessenberg matrix and do Gram-Schmidt */
149:     PetscCall((*gmres->orthog)(ksp, it));
150:     if (ksp->reason) break;

152:     /* vv(i+1) . vv(i+1) */
153:     PetscCall(VecNormalize(VEC_VV(it + 1), &tt));
154:     KSPCheckNorm(ksp, tt);

156:     /* save the magnitude */
157:     *HH(it + 1, it)  = tt;
158:     *HES(it + 1, it) = tt;

160:     /* check for the happy breakdown */
161:     hapbnd = PetscAbsScalar(tt / *GRS(it));
162:     if (hapbnd > gmres->haptol) hapbnd = gmres->haptol;
163:     if (tt < hapbnd) {
164:       PetscCall(PetscInfo(ksp, "Detected happy breakdown, current hapbnd = %14.12e tt = %14.12e\n", (double)hapbnd, (double)tt));
165:       hapend = PETSC_TRUE;
166:     }
167:     PetscCall(KSPGMRESUpdateHessenberg(ksp, it, hapend, &res));

169:     it++;
170:     gmres->it = (it - 1); /* For converged */
171:     ksp->its++;
172:     ksp->rnorm = res;
173:     if (ksp->reason) break;

175:     PetscCall((*ksp->converged)(ksp, ksp->its, res, &ksp->reason, ksp->cnvP));

177:     /* Catch error in happy breakdown and signal convergence and break from loop */
178:     if (hapend) {
179:       if (ksp->normtype == KSP_NORM_NONE) { /* convergence test was skipped in this case */
180:         ksp->reason = KSP_CONVERGED_HAPPY_BREAKDOWN;
181:       } else if (!ksp->reason) {
182:         PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "Reached happy break down, but convergence was not indicated. Residual norm = %g", (double)res);
183:         ksp->reason = KSP_DIVERGED_BREAKDOWN;
184:         break;
185:       }
186:     }
187:   }

189:   if (itcount) *itcount = it;

191:   /*
192:     Down here we have to solve for the "best" coefficients of the Krylov
193:     columns, add the solution values together, and possibly unwind the
194:     preconditioning from the solution
195:    */
196:   /* Form the solution (or the solution so far) */
197:   PetscCall(KSPGMRESBuildSoln(GRS(0), ksp->vec_sol, ksp->vec_sol, ksp, it - 1));

199:   /* Monitor if we know that we will not return for a restart */
200:   if (it && (ksp->reason || ksp->its >= ksp->max_it)) {
201:     PetscCall(KSPLogResidualHistory(ksp, res));
202:     PetscCall(KSPLogErrorHistory(ksp));
203:     PetscCall(KSPMonitor(ksp, ksp->its, res));
204:   }
205:   PetscFunctionReturn(PETSC_SUCCESS);
206: }

208: static PetscErrorCode KSPSolve_GMRES(KSP ksp)
209: {
210:   PetscInt   its, itcount, i;
211:   KSP_GMRES *gmres      = (KSP_GMRES *)ksp->data;
212:   PetscBool  guess_zero = ksp->guess_zero;
213:   PetscInt   N          = gmres->max_k + 1;

215:   PetscFunctionBegin;
216:   PetscCheck(!ksp->calc_sings || gmres->Rsvd, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ORDER, "Must call KSPSetComputeSingularValues() before KSPSetUp() is called");

218:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
219:   ksp->its = 0;
220:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));

222:   itcount          = 0;
223:   gmres->fullcycle = 0;
224:   ksp->rnorm       = -1.0; /* special marker for KSPGMRESCycle() */
225:   while (!ksp->reason || (ksp->rnorm == -1 && ksp->reason == KSP_DIVERGED_PC_FAILED)) {
226:     PetscCall(KSPInitialResidual(ksp, ksp->vec_sol, VEC_TEMP, VEC_TEMP_MATOP, VEC_VV(0), ksp->vec_rhs));
227:     PetscCall(KSPGMRESCycle(&its, ksp));
228:     /* Store the Hessenberg matrix and the basis vectors of the Krylov subspace
229:     if the cycle is complete for the computation of the Ritz pairs */
230:     if (its == gmres->max_k) {
231:       gmres->fullcycle++;
232:       if (ksp->calc_ritz) {
233:         if (!gmres->hes_ritz) {
234:           PetscCall(PetscMalloc1(N * N, &gmres->hes_ritz));
235:           PetscCall(VecDuplicateVecs(VEC_VV(0), N, &gmres->vecb));
236:         }
237:         PetscCall(PetscArraycpy(gmres->hes_ritz, gmres->hes_origin, N * N));
238:         for (i = 0; i < gmres->max_k + 1; i++) PetscCall(VecCopy(VEC_VV(i), gmres->vecb[i]));
239:       }
240:     }
241:     itcount += its;
242:     if (itcount >= ksp->max_it) {
243:       if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
244:       break;
245:     }
246:     ksp->guess_zero = PETSC_FALSE; /* every future call to KSPInitialResidual() will have nonzero guess */
247:   }
248:   ksp->guess_zero = guess_zero; /* restore if user provided nonzero initial guess */
249:   PetscFunctionReturn(PETSC_SUCCESS);
250: }

252: PetscErrorCode KSPReset_GMRES(KSP ksp)
253: {
254:   KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;
255:   PetscInt   i;

257:   PetscFunctionBegin;
258:   /* Free the Hessenberg matrices */
259:   PetscCall(PetscFree5(gmres->hh_origin, gmres->hes_origin, gmres->rs_origin, gmres->cc_origin, gmres->ss_origin));
260:   PetscCall(PetscFree(gmres->hes_ritz));

262:   /* free work vectors */
263:   PetscCall(PetscFree(gmres->vecs));
264:   for (i = 0; i < gmres->nwork_alloc; i++) PetscCall(VecDestroyVecs(gmres->mwork_alloc[i], &gmres->user_work[i]));
265:   gmres->nwork_alloc = 0;
266:   if (gmres->vecb) PetscCall(VecDestroyVecs(gmres->max_k + 1, &gmres->vecb));

268:   PetscCall(PetscFree(gmres->user_work));
269:   PetscCall(PetscFree(gmres->mwork_alloc));
270:   PetscCall(PetscFree(gmres->nrs));
271:   PetscCall(VecDestroy(&gmres->sol_temp));
272:   PetscCall(PetscFree(gmres->Rsvd));
273:   PetscCall(PetscFree(gmres->Dsvd));
274:   PetscCall(PetscFree(gmres->orthogwork));

276:   gmres->vv_allocated   = 0;
277:   gmres->vecs_allocated = 0;
278:   gmres->sol_temp       = NULL;
279:   PetscFunctionReturn(PETSC_SUCCESS);
280: }

282: PetscErrorCode KSPDestroy_GMRES(KSP ksp)
283: {
284:   PetscFunctionBegin;
285:   PetscCall(KSPReset_GMRES(ksp));
286:   PetscCall(PetscFree(ksp->data));
287:   /* clear composed functions */
288:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetPreAllocateVectors_C", NULL));
289:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetOrthogonalization_C", NULL));
290:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetOrthogonalization_C", NULL));
291:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetRestart_C", NULL));
292:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetRestart_C", NULL));
293:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetHapTol_C", NULL));
294:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetBreakdownTolerance_C", NULL));
295:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetCGSRefinementType_C", NULL));
296:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetCGSRefinementType_C", NULL));
297:   PetscFunctionReturn(PETSC_SUCCESS);
298: }
299: /*
300:     KSPGMRESBuildSoln - create the solution from the starting vector and the
301:     current iterates.

303:     Input parameters:
304:         nrs - work area of size it + 1.
305:         vs  - index of initial guess
306:         vdest - index of result.  Note that vs may == vdest (replace
307:                 guess with the solution).

309:      This is an internal routine that knows about the GMRES internals.
310:  */
311: static PetscErrorCode KSPGMRESBuildSoln(PetscScalar *nrs, Vec vs, Vec vdest, KSP ksp, PetscInt it)
312: {
313:   PetscScalar tt;
314:   PetscInt    ii, k, j;
315:   KSP_GMRES  *gmres = (KSP_GMRES *)ksp->data;

317:   PetscFunctionBegin;
318:   /* Solve for solution vector that minimizes the residual */

320:   /* If it is < 0, no gmres steps have been performed */
321:   if (it < 0) {
322:     PetscCall(VecCopy(vs, vdest)); /* VecCopy() is smart, exists immediately if vguess == vdest */
323:     PetscFunctionReturn(PETSC_SUCCESS);
324:   }
325:   if (*HH(it, it) != 0.0) {
326:     nrs[it] = *GRS(it) / *HH(it, it);
327:   } else {
328:     PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "You reached the break down in GMRES; HH(it,it) = 0");
329:     ksp->reason = KSP_DIVERGED_BREAKDOWN;

331:     PetscCall(PetscInfo(ksp, "Likely your matrix or preconditioner is singular. HH(it,it) is identically zero; it = %" PetscInt_FMT " GRS(it) = %g\n", it, (double)PetscAbsScalar(*GRS(it))));
332:     PetscFunctionReturn(PETSC_SUCCESS);
333:   }
334:   for (ii = 1; ii <= it; ii++) {
335:     k  = it - ii;
336:     tt = *GRS(k);
337:     for (j = k + 1; j <= it; j++) tt = tt - *HH(k, j) * nrs[j];
338:     if (*HH(k, k) == 0.0) {
339:       PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "Likely your matrix or preconditioner is singular. HH(k,k) is identically zero; k = %" PetscInt_FMT, k);
340:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
341:       PetscCall(PetscInfo(ksp, "Likely your matrix or preconditioner is singular. HH(k,k) is identically zero; k = %" PetscInt_FMT "\n", k));
342:       PetscFunctionReturn(PETSC_SUCCESS);
343:     }
344:     nrs[k] = tt / *HH(k, k);
345:   }

347:   /* Accumulate the correction to the solution of the preconditioned problem in TEMP */
348:   PetscCall(VecMAXPBY(VEC_TEMP, it + 1, nrs, 0, &VEC_VV(0)));

350:   PetscCall(KSPUnwindPreconditioner(ksp, VEC_TEMP, VEC_TEMP_MATOP));
351:   /* add solution to previous solution */
352:   if (vdest != vs) PetscCall(VecCopy(vs, vdest));
353:   PetscCall(VecAXPY(vdest, 1.0, VEC_TEMP));
354:   PetscFunctionReturn(PETSC_SUCCESS);
355: }
356: /*
357:    Do the scalar work for the orthogonalization.  Return new residual norm.
358:  */
359: static PetscErrorCode KSPGMRESUpdateHessenberg(KSP ksp, PetscInt it, PetscBool hapend, PetscReal *res)
360: {
361:   PetscScalar *hh, *cc, *ss, tt;
362:   PetscInt     j;
363:   KSP_GMRES   *gmres = (KSP_GMRES *)ksp->data;

365:   PetscFunctionBegin;
366:   hh = HH(0, it);
367:   cc = CC(0);
368:   ss = SS(0);

370:   /* Apply all the previously computed plane rotations to the new column
371:      of the Hessenberg matrix */
372:   for (j = 1; j <= it; j++) {
373:     tt  = *hh;
374:     *hh = PetscConj(*cc) * tt + *ss * *(hh + 1);
375:     hh++;
376:     *hh = *cc++ * *hh - (*ss++ * tt);
377:   }

379:   /*
380:     compute the new plane rotation, and apply it to:
381:      1) the right-hand side of the Hessenberg system
382:      2) the new column of the Hessenberg matrix
383:     thus obtaining the updated value of the residual
384:   */
385:   if (!hapend) {
386:     tt = PetscSqrtScalar(PetscConj(*hh) * *hh + PetscConj(*(hh + 1)) * *(hh + 1));
387:     if (tt == 0.0) {
388:       PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "tt == 0.0");
389:       ksp->reason = KSP_DIVERGED_NULL;
390:       PetscFunctionReturn(PETSC_SUCCESS);
391:     }
392:     *cc          = *hh / tt;
393:     *ss          = *(hh + 1) / tt;
394:     *GRS(it + 1) = -(*ss * *GRS(it));
395:     *GRS(it)     = PetscConj(*cc) * *GRS(it);
396:     *hh          = PetscConj(*cc) * *hh + *ss * *(hh + 1);
397:     *res         = PetscAbsScalar(*GRS(it + 1));
398:   } else {
399:     /* happy breakdown: HH(it+1, it) = 0, therefore we don't need to apply
400:             another rotation matrix (so RH doesn't change).  The new residual is
401:             always the new sine term times the residual from last time (GRS(it)),
402:             but now the new sine rotation would be zero...so the residual should
403:             be zero...so we will multiply "zero" by the last residual.  This might
404:             not be exactly what we want to do here -could just return "zero". */

406:     *res = 0.0;
407:   }
408:   PetscFunctionReturn(PETSC_SUCCESS);
409: }
410: /*
411:    This routine allocates more work vectors, starting from VEC_VV(it).
412:  */
413: PetscErrorCode KSPGMRESGetNewVectors(KSP ksp, PetscInt it)
414: {
415:   KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;
416:   PetscInt   nwork = gmres->nwork_alloc, k, nalloc;

418:   PetscFunctionBegin;
419:   nalloc = PetscMin(ksp->max_it, gmres->delta_allocate);
420:   /* Adjust the number to allocate to make sure that we don't exceed the
421:     number of available slots */
422:   if (it + VEC_OFFSET + nalloc >= gmres->vecs_allocated) nalloc = gmres->vecs_allocated - it - VEC_OFFSET;
423:   if (!nalloc) PetscFunctionReturn(PETSC_SUCCESS);

425:   gmres->vv_allocated += nalloc;

427:   PetscCall(KSPCreateVecs(ksp, nalloc, &gmres->user_work[nwork], 0, NULL));

429:   gmres->mwork_alloc[nwork] = nalloc;
430:   for (k = 0; k < nalloc; k++) gmres->vecs[it + VEC_OFFSET + k] = gmres->user_work[nwork][k];
431:   gmres->nwork_alloc++;
432:   PetscFunctionReturn(PETSC_SUCCESS);
433: }

435: static PetscErrorCode KSPBuildSolution_GMRES(KSP ksp, Vec ptr, Vec *result)
436: {
437:   KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;

439:   PetscFunctionBegin;
440:   if (!ptr) {
441:     if (!gmres->sol_temp) PetscCall(VecDuplicate(ksp->vec_sol, &gmres->sol_temp));
442:     ptr = gmres->sol_temp;
443:   }
444:   if (!gmres->nrs) {
445:     /* allocate the work area */
446:     PetscCall(PetscMalloc1(gmres->max_k, &gmres->nrs));
447:   }

449:   PetscCall(KSPGMRESBuildSoln(gmres->nrs, ksp->vec_sol, ptr, ksp, gmres->it));
450:   if (result) *result = ptr;
451:   PetscFunctionReturn(PETSC_SUCCESS);
452: }

454: PetscErrorCode KSPView_GMRES(KSP ksp, PetscViewer viewer)
455: {
456:   KSP_GMRES  *gmres = (KSP_GMRES *)ksp->data;
457:   const char *cstr;
458:   PetscBool   iascii, isstring;

460:   PetscFunctionBegin;
461:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
462:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
463:   if (gmres->orthog == KSPGMRESClassicalGramSchmidtOrthogonalization) {
464:     switch (gmres->cgstype) {
465:     case (KSP_GMRES_CGS_REFINE_NEVER):
466:       cstr = "Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement";
467:       break;
468:     case (KSP_GMRES_CGS_REFINE_ALWAYS):
469:       cstr = "Classical (unmodified) Gram-Schmidt Orthogonalization with one step of iterative refinement";
470:       break;
471:     case (KSP_GMRES_CGS_REFINE_IFNEEDED):
472:       cstr = "Classical (unmodified) Gram-Schmidt Orthogonalization with one step of iterative refinement when needed";
473:       break;
474:     default:
475:       SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Unknown orthogonalization");
476:     }
477:   } else if (gmres->orthog == KSPGMRESModifiedGramSchmidtOrthogonalization) {
478:     cstr = "Modified Gram-Schmidt Orthogonalization";
479:   } else {
480:     cstr = "unknown orthogonalization";
481:   }
482:   if (iascii) {
483:     PetscCall(PetscViewerASCIIPrintf(viewer, "  restart=%" PetscInt_FMT ", using %s\n", gmres->max_k, cstr));
484:     PetscCall(PetscViewerASCIIPrintf(viewer, "  happy breakdown tolerance %g\n", (double)gmres->haptol));
485:   } else if (isstring) {
486:     PetscCall(PetscViewerStringSPrintf(viewer, "%s restart %" PetscInt_FMT, cstr, gmres->max_k));
487:   }
488:   PetscFunctionReturn(PETSC_SUCCESS);
489: }

491: /*@C
492:   KSPGMRESMonitorKrylov - Calls `VecView()` for each new direction in the `KSPGMRES` accumulated Krylov space.

494:   Collective

496:   Input Parameters:
497: + ksp    - the `KSP` context
498: . its    - iteration number
499: . fgnorm - 2-norm of residual (or gradient)
500: - dummy  - a collection of viewers created with `PetscViewersCreate()`

502:   Options Database Key:
503: . -ksp_gmres_krylov_monitor <bool> - Plot the Krylov directions

505:   Level: intermediate

507:   Note:
508:   A new `PETSCVIEWERDRAW` is created for each Krylov vector so they can all be simultaneously viewed

510: .seealso: [](ch_ksp), `KSPGMRES`, `KSPMonitorSet()`, `KSPMonitorResidual()`, `VecView()`, `PetscViewersCreate()`, `PetscViewersDestroy()`
511: @*/
512: PetscErrorCode KSPGMRESMonitorKrylov(KSP ksp, PetscInt its, PetscReal fgnorm, void *dummy)
513: {
514:   PetscViewers viewers = (PetscViewers)dummy;
515:   KSP_GMRES   *gmres   = (KSP_GMRES *)ksp->data;
516:   Vec          x;
517:   PetscViewer  viewer;
518:   PetscBool    flg;

520:   PetscFunctionBegin;
521:   PetscCall(PetscViewersGetViewer(viewers, gmres->it + 1, &viewer));
522:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &flg));
523:   if (!flg) {
524:     PetscCall(PetscViewerSetType(viewer, PETSCVIEWERDRAW));
525:     PetscCall(PetscViewerDrawSetInfo(viewer, NULL, "Krylov GMRES Monitor", PETSC_DECIDE, PETSC_DECIDE, 300, 300));
526:   }
527:   x = VEC_VV(gmres->it + 1);
528:   PetscCall(VecView(x, viewer));
529:   PetscFunctionReturn(PETSC_SUCCESS);
530: }

532: PetscErrorCode KSPSetFromOptions_GMRES(KSP ksp, PetscOptionItems *PetscOptionsObject)
533: {
534:   PetscInt   restart;
535:   PetscReal  haptol, breakdowntol;
536:   KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;
537:   PetscBool  flg;

539:   PetscFunctionBegin;
540:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP GMRES Options");
541:   PetscCall(PetscOptionsInt("-ksp_gmres_restart", "Number of Krylov search directions", "KSPGMRESSetRestart", gmres->max_k, &restart, &flg));
542:   if (flg) PetscCall(KSPGMRESSetRestart(ksp, restart));
543:   PetscCall(PetscOptionsReal("-ksp_gmres_haptol", "Tolerance for exact convergence (happy ending)", "KSPGMRESSetHapTol", gmres->haptol, &haptol, &flg));
544:   if (flg) PetscCall(KSPGMRESSetHapTol(ksp, haptol));
545:   PetscCall(PetscOptionsReal("-ksp_gmres_breakdown_tolerance", "Divergence breakdown tolerance during GMRES restart", "KSPGMRESSetBreakdownTolerance", gmres->breakdowntol, &breakdowntol, &flg));
546:   if (flg) PetscCall(KSPGMRESSetBreakdownTolerance(ksp, breakdowntol));
547:   flg = PETSC_FALSE;
548:   PetscCall(PetscOptionsBool("-ksp_gmres_preallocate", "Preallocate Krylov vectors", "KSPGMRESSetPreAllocateVectors", flg, &flg, NULL));
549:   if (flg) PetscCall(KSPGMRESSetPreAllocateVectors(ksp));
550:   PetscCall(PetscOptionsBoolGroupBegin("-ksp_gmres_classicalgramschmidt", "Classical (unmodified) Gram-Schmidt (fast)", "KSPGMRESSetOrthogonalization", &flg));
551:   if (flg) PetscCall(KSPGMRESSetOrthogonalization(ksp, KSPGMRESClassicalGramSchmidtOrthogonalization));
552:   PetscCall(PetscOptionsBoolGroupEnd("-ksp_gmres_modifiedgramschmidt", "Modified Gram-Schmidt (slow,more stable)", "KSPGMRESSetOrthogonalization", &flg));
553:   if (flg) PetscCall(KSPGMRESSetOrthogonalization(ksp, KSPGMRESModifiedGramSchmidtOrthogonalization));
554:   PetscCall(PetscOptionsEnum("-ksp_gmres_cgs_refinement_type", "Type of iterative refinement for classical (unmodified) Gram-Schmidt", "KSPGMRESSetCGSRefinementType", KSPGMRESCGSRefinementTypes, (PetscEnum)gmres->cgstype, (PetscEnum *)&gmres->cgstype, &flg));
555:   flg = PETSC_FALSE;
556:   PetscCall(PetscOptionsBool("-ksp_gmres_krylov_monitor", "Plot the Krylov directions", "KSPMonitorSet", flg, &flg, NULL));
557:   if (flg) {
558:     PetscViewers viewers;
559:     PetscCall(PetscViewersCreate(PetscObjectComm((PetscObject)ksp), &viewers));
560:     PetscCall(KSPMonitorSet(ksp, KSPGMRESMonitorKrylov, viewers, (PetscErrorCode(*)(void **))PetscViewersDestroy));
561:   }
563:   PetscFunctionReturn(PETSC_SUCCESS);
564: }

566: PetscErrorCode KSPGMRESSetHapTol_GMRES(KSP ksp, PetscReal tol)
567: {
568:   KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;

570:   PetscFunctionBegin;
571:   PetscCheck(tol >= 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Tolerance must be non-negative");
572:   gmres->haptol = tol;
573:   PetscFunctionReturn(PETSC_SUCCESS);
574: }

576: static PetscErrorCode KSPGMRESSetBreakdownTolerance_GMRES(KSP ksp, PetscReal tol)
577: {
578:   KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;

580:   PetscFunctionBegin;
581:   if (tol == (PetscReal)PETSC_DEFAULT) {
582:     gmres->breakdowntol = 0.1;
583:     PetscFunctionReturn(PETSC_SUCCESS);
584:   }
585:   PetscCheck(tol >= 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Breakdown tolerance must be non-negative");
586:   gmres->breakdowntol = tol;
587:   PetscFunctionReturn(PETSC_SUCCESS);
588: }

590: PetscErrorCode KSPGMRESGetRestart_GMRES(KSP ksp, PetscInt *max_k)
591: {
592:   KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;

594:   PetscFunctionBegin;
595:   *max_k = gmres->max_k;
596:   PetscFunctionReturn(PETSC_SUCCESS);
597: }

599: PetscErrorCode KSPGMRESSetRestart_GMRES(KSP ksp, PetscInt max_k)
600: {
601:   KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;

603:   PetscFunctionBegin;
604:   PetscCheck(max_k >= 1, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Restart must be positive");
605:   if (!ksp->setupstage) {
606:     gmres->max_k = max_k;
607:   } else if (gmres->max_k != max_k) {
608:     gmres->max_k    = max_k;
609:     ksp->setupstage = KSP_SETUP_NEW;
610:     /* free the data structures, then create them again */
611:     PetscCall(KSPReset_GMRES(ksp));
612:   }
613:   PetscFunctionReturn(PETSC_SUCCESS);
614: }

616: PetscErrorCode KSPGMRESSetOrthogonalization_GMRES(KSP ksp, FCN fcn)
617: {
618:   PetscFunctionBegin;
619:   ((KSP_GMRES *)ksp->data)->orthog = fcn;
620:   PetscFunctionReturn(PETSC_SUCCESS);
621: }

623: PetscErrorCode KSPGMRESGetOrthogonalization_GMRES(KSP ksp, FCN *fcn)
624: {
625:   PetscFunctionBegin;
626:   *fcn = ((KSP_GMRES *)ksp->data)->orthog;
627:   PetscFunctionReturn(PETSC_SUCCESS);
628: }

630: PetscErrorCode KSPGMRESSetPreAllocateVectors_GMRES(KSP ksp)
631: {
632:   KSP_GMRES *gmres;

634:   PetscFunctionBegin;
635:   gmres                = (KSP_GMRES *)ksp->data;
636:   gmres->q_preallocate = 1;
637:   PetscFunctionReturn(PETSC_SUCCESS);
638: }

640: PetscErrorCode KSPGMRESSetCGSRefinementType_GMRES(KSP ksp, KSPGMRESCGSRefinementType type)
641: {
642:   KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;

644:   PetscFunctionBegin;
645:   gmres->cgstype = type;
646:   PetscFunctionReturn(PETSC_SUCCESS);
647: }

649: PetscErrorCode KSPGMRESGetCGSRefinementType_GMRES(KSP ksp, KSPGMRESCGSRefinementType *type)
650: {
651:   KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;

653:   PetscFunctionBegin;
654:   *type = gmres->cgstype;
655:   PetscFunctionReturn(PETSC_SUCCESS);
656: }

658: /*@
659:   KSPGMRESSetCGSRefinementType - Sets the type of iterative refinement to use
660:   in the classical Gram-Schmidt orthogonalization.

662:   Logically Collective

664:   Input Parameters:
665: + ksp  - the Krylov space context
666: - type - the type of refinement
667: .vb
668:   KSP_GMRES_CGS_REFINE_NEVER
669:   KSP_GMRES_CGS_REFINE_IFNEEDED
670:   KSP_GMRES_CGS_REFINE_ALWAYS
671: .ve

673:   Options Database Key:
674: . -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always> - refinement type

676:   Level: intermediate

678: .seealso: [](ch_ksp), `KSPGMRES`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESCGSRefinementType`, `KSPGMRESClassicalGramSchmidtOrthogonalization()`, `KSPGMRESGetCGSRefinementType()`,
679:           `KSPGMRESGetOrthogonalization()`
680: @*/
681: PetscErrorCode KSPGMRESSetCGSRefinementType(KSP ksp, KSPGMRESCGSRefinementType type)
682: {
683:   PetscFunctionBegin;
686:   PetscTryMethod(ksp, "KSPGMRESSetCGSRefinementType_C", (KSP, KSPGMRESCGSRefinementType), (ksp, type));
687:   PetscFunctionReturn(PETSC_SUCCESS);
688: }

690: /*@
691:   KSPGMRESGetCGSRefinementType - Gets the type of iterative refinement to use
692:   in the classical Gram Schmidt orthogonalization.

694:   Not Collective

696:   Input Parameter:
697: . ksp - the Krylov space context

699:   Output Parameter:
700: . type - the type of refinement

702:   Level: intermediate

704: .seealso: [](ch_ksp), `KSPGMRES`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESCGSRefinementType`, `KSPGMRESClassicalGramSchmidtOrthogonalization()`, `KSPGMRESSetCGSRefinementType()`,
705:           `KSPGMRESGetOrthogonalization()`
706: @*/
707: PetscErrorCode KSPGMRESGetCGSRefinementType(KSP ksp, KSPGMRESCGSRefinementType *type)
708: {
709:   PetscFunctionBegin;
711:   PetscUseMethod(ksp, "KSPGMRESGetCGSRefinementType_C", (KSP, KSPGMRESCGSRefinementType *), (ksp, type));
712:   PetscFunctionReturn(PETSC_SUCCESS);
713: }

715: /*@
716:   KSPGMRESSetRestart - Sets number of iterations at which `KSPGMRES`, `KSPFGMRES` and `KSPLGMRES` restarts.

718:   Logically Collective

720:   Input Parameters:
721: + ksp     - the Krylov space context
722: - restart - integer restart value

724:   Options Database Key:
725: . -ksp_gmres_restart <positive integer> - integer restart value

727:   Level: intermediate

729:   Note:
730:   The default value is 30.

732: .seealso: [](ch_ksp), `KSPGMRES`, `KSPSetTolerances()`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESSetPreAllocateVectors()`, `KSPGMRESGetRestart()`
733: @*/
734: PetscErrorCode KSPGMRESSetRestart(KSP ksp, PetscInt restart)
735: {
736:   PetscFunctionBegin;

739:   PetscTryMethod(ksp, "KSPGMRESSetRestart_C", (KSP, PetscInt), (ksp, restart));
740:   PetscFunctionReturn(PETSC_SUCCESS);
741: }

743: /*@
744:   KSPGMRESGetRestart - Gets number of iterations at which `KSPGMRES`, `KSPFGMRES` and `KSPLGMRES` restarts.

746:   Not Collective

748:   Input Parameter:
749: . ksp - the Krylov space context

751:   Output Parameter:
752: . restart - integer restart value

754:   Level: intermediate

756: .seealso: [](ch_ksp), `KSPGMRES`, `KSPSetTolerances()`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESSetPreAllocateVectors()`, `KSPGMRESSetRestart()`
757: @*/
758: PetscErrorCode KSPGMRESGetRestart(KSP ksp, PetscInt *restart)
759: {
760:   PetscFunctionBegin;
761:   PetscUseMethod(ksp, "KSPGMRESGetRestart_C", (KSP, PetscInt *), (ksp, restart));
762:   PetscFunctionReturn(PETSC_SUCCESS);
763: }

765: /*@
766:   KSPGMRESSetHapTol - Sets tolerance for determining happy breakdown in `KSPGMRES`, `KSPFGMRES` and `KSPLGMRES`

768:   Logically Collective

770:   Input Parameters:
771: + ksp - the Krylov space context
772: - tol - the tolerance

774:   Options Database Key:
775: . -ksp_gmres_haptol <positive real value> - set tolerance for determining happy breakdown

777:   Level: intermediate

779:   Note:
780:   Happy breakdown is the rare case in `KSPGMRES` where an 'exact' solution is obtained after
781:   a certain number of iterations. If you attempt more iterations after this point unstable
782:   things can happen hence very occasionally you may need to set this value to detect this condition

784: .seealso: [](ch_ksp), `KSPGMRES`, `KSPSetTolerances()`
785: @*/
786: PetscErrorCode KSPGMRESSetHapTol(KSP ksp, PetscReal tol)
787: {
788:   PetscFunctionBegin;
790:   PetscTryMethod((ksp), "KSPGMRESSetHapTol_C", (KSP, PetscReal), ((ksp), (tol)));
791:   PetscFunctionReturn(PETSC_SUCCESS);
792: }

794: /*@
795:   KSPGMRESSetBreakdownTolerance - Sets tolerance for determining divergence breakdown in `KSPGMRES`.

797:   Logically Collective

799:   Input Parameters:
800: + ksp - the Krylov space context
801: - tol - the tolerance

803:   Options Database Key:
804: . -ksp_gmres_breakdown_tolerance <positive real value> - set tolerance for determining divergence breakdown

806:   Level: intermediate

808:   Note:
809:   Divergence breakdown occurs when GMRES residual increases significantly during restart

811: .seealso: [](ch_ksp), `KSPGMRES`, `KSPSetTolerances()`, `KSPGMRESSetHapTol()`
812: @*/
813: PetscErrorCode KSPGMRESSetBreakdownTolerance(KSP ksp, PetscReal tol)
814: {
815:   PetscFunctionBegin;
817:   PetscTryMethod((ksp), "KSPGMRESSetBreakdownTolerance_C", (KSP, PetscReal), (ksp, tol));
818:   PetscFunctionReturn(PETSC_SUCCESS);
819: }

821: /*MC
822:      KSPGMRES - Implements the Generalized Minimal Residual method {cite}`saad.schultz:gmres` with restart

824:    Options Database Keys:
825: +   -ksp_gmres_restart <restart> - the number of Krylov directions to orthogonalize against
826: .   -ksp_gmres_haptol <tol> - sets the tolerance for "happy ending" (exact convergence)
827: .   -ksp_gmres_preallocate - preallocate all the Krylov search directions initially (otherwise groups of
828:                              vectors are allocated as needed)
829: .   -ksp_gmres_classicalgramschmidt - use classical (unmodified) Gram-Schmidt to orthogonalize against the Krylov space (fast) (the default)
830: .   -ksp_gmres_modifiedgramschmidt - use modified Gram-Schmidt in the orthogonalization (more stable, but slower)
831: .   -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always> - determine if iterative refinement is used to increase the
832:                                    stability of the classical Gram-Schmidt  orthogonalization.
833: -   -ksp_gmres_krylov_monitor - plot the Krylov space generated

835:    Level: beginner

837:    Note:
838:    Left and right preconditioning are supported, but not symmetric preconditioning.

840: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPFGMRES`, `KSPLGMRES`,
841:           `KSPGMRESSetRestart()`, `KSPGMRESSetHapTol()`, `KSPGMRESSetPreAllocateVectors()`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESGetOrthogonalization()`,
842:           `KSPGMRESClassicalGramSchmidtOrthogonalization()`, `KSPGMRESModifiedGramSchmidtOrthogonalization()`,
843:           `KSPGMRESCGSRefinementType`, `KSPGMRESSetCGSRefinementType()`, `KSPGMRESGetCGSRefinementType()`, `KSPGMRESMonitorKrylov()`, `KSPSetPCSide()`
844: M*/

846: PETSC_EXTERN PetscErrorCode KSPCreate_GMRES(KSP ksp)
847: {
848:   KSP_GMRES *gmres;

850:   PetscFunctionBegin;
851:   PetscCall(PetscNew(&gmres));
852:   ksp->data = (void *)gmres;

854:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 4));
855:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 3));
856:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_SYMMETRIC, 2));
857:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));
858:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));

860:   ksp->ops->buildsolution                = KSPBuildSolution_GMRES;
861:   ksp->ops->setup                        = KSPSetUp_GMRES;
862:   ksp->ops->solve                        = KSPSolve_GMRES;
863:   ksp->ops->reset                        = KSPReset_GMRES;
864:   ksp->ops->destroy                      = KSPDestroy_GMRES;
865:   ksp->ops->view                         = KSPView_GMRES;
866:   ksp->ops->setfromoptions               = KSPSetFromOptions_GMRES;
867:   ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_GMRES;
868:   ksp->ops->computeeigenvalues           = KSPComputeEigenvalues_GMRES;
869:   ksp->ops->computeritz                  = KSPComputeRitz_GMRES;
870:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetPreAllocateVectors_C", KSPGMRESSetPreAllocateVectors_GMRES));
871:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetOrthogonalization_C", KSPGMRESSetOrthogonalization_GMRES));
872:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetOrthogonalization_C", KSPGMRESGetOrthogonalization_GMRES));
873:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetRestart_C", KSPGMRESSetRestart_GMRES));
874:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetRestart_C", KSPGMRESGetRestart_GMRES));
875:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetHapTol_C", KSPGMRESSetHapTol_GMRES));
876:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetBreakdownTolerance_C", KSPGMRESSetBreakdownTolerance_GMRES));
877:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetCGSRefinementType_C", KSPGMRESSetCGSRefinementType_GMRES));
878:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetCGSRefinementType_C", KSPGMRESGetCGSRefinementType_GMRES));

880:   gmres->haptol         = 1.0e-30;
881:   gmres->breakdowntol   = 0.1;
882:   gmres->q_preallocate  = 0;
883:   gmres->delta_allocate = GMRES_DELTA_DIRECTIONS;
884:   gmres->orthog         = KSPGMRESClassicalGramSchmidtOrthogonalization;
885:   gmres->nrs            = NULL;
886:   gmres->sol_temp       = NULL;
887:   gmres->max_k          = GMRES_DEFAULT_MAXK;
888:   gmres->Rsvd           = NULL;
889:   gmres->cgstype        = KSP_GMRES_CGS_REFINE_NEVER;
890:   gmres->orthogwork     = NULL;
891:   PetscFunctionReturn(PETSC_SUCCESS);
892: }
```