Actual source code: gmres.c


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

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

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

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

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

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

 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:   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:     PetscMalloc1((max_k + 3) * (max_k + 9), &gmres->Rsvd);
 53:     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:   PetscMalloc1(gmres->vecs_allocated, &gmres->vecs);
 61:   PetscMalloc1(VEC_OFFSET + 2 + max_k, &gmres->user_work);
 62:   PetscMalloc1(VEC_OFFSET + 2 + max_k, &gmres->mwork_alloc);

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

 67:     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:     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:   return 0;
 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: 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:   if (itcount) *itcount = 0;
109:   VecNormalize(VEC_VV(0), &res);
110:   KSPCheckNorm(ksp, res);

112:   /* the constant .1 is arbitrary, just some measure at how incorrect the residuals are */
113:   if ((ksp->rnorm > 0.0) && (PetscAbsReal(res - ksp->rnorm) > gmres->breakdowntol * gmres->rnorm0)) {
115:                (double)ksp->rnorm, (double)res, (double)gmres->rnorm0);
116:     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", (double)ksp->rnorm, (double)res, (double)gmres->rnorm0);
117:     ksp->reason = KSP_DIVERGED_BREAKDOWN;
118:     return 0;
119:   }
120:   *GRS(0) = gmres->rnorm0 = res;

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

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

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

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

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

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

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

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

176:     /* Catch error in happy breakdown and signal convergence and break from loop */
177:     if (hapend) {
178:       if (ksp->normtype == KSP_NORM_NONE) { /* convergence test was skipped in this case */
179:         ksp->reason = KSP_CONVERGED_HAPPY_BREAKDOWN;
180:       } else if (!ksp->reason) {
182:         ksp->reason = KSP_DIVERGED_BREAKDOWN;
183:         break;
184:       }
185:     }
186:   }

188:   /* Monitor if we know that we will not return for a restart */
189:   if (it && (ksp->reason || ksp->its >= ksp->max_it)) {
190:     KSPLogResidualHistory(ksp, res);
191:     KSPLogErrorHistory(ksp);
192:     KSPMonitor(ksp, ksp->its, res);
193:   }

195:   if (itcount) *itcount = it;

197:   /*
198:     Down here we have to solve for the "best" coefficients of the Krylov
199:     columns, add the solution values together, and possibly unwind the
200:     preconditioning from the solution
201:    */
202:   /* Form the solution (or the solution so far) */
203:   KSPGMRESBuildSoln(GRS(0), ksp->vec_sol, ksp->vec_sol, ksp, it - 1);
204:   return 0;
205: }

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


216:   PetscObjectSAWsTakeAccess((PetscObject)ksp);
217:   ksp->its = 0;
218:   PetscObjectSAWsGrantAccess((PetscObject)ksp);

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

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

255:   /* Free the Hessenberg matrices */
256:   PetscFree5(gmres->hh_origin, gmres->hes_origin, gmres->rs_origin, gmres->cc_origin, gmres->ss_origin);
257:   PetscFree(gmres->hes_ritz);

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

265:   PetscFree(gmres->user_work);
266:   PetscFree(gmres->mwork_alloc);
267:   PetscFree(gmres->nrs);
268:   VecDestroy(&gmres->sol_temp);
269:   PetscFree(gmres->Rsvd);
270:   PetscFree(gmres->Dsvd);
271:   PetscFree(gmres->orthogwork);

273:   gmres->vv_allocated   = 0;
274:   gmres->vecs_allocated = 0;
275:   gmres->sol_temp       = NULL;
276:   return 0;
277: }

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

299:     Input parameters:
300:         nrs - work area of size it + 1.
301:         vs  - index of initial guess
302:         vdest - index of result.  Note that vs may == vdest (replace
303:                 guess with the solution).

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

313:   /* Solve for solution vector that minimizes the residual */

315:   /* If it is < 0, no gmres steps have been performed */
316:   if (it < 0) {
317:     VecCopy(vs, vdest); /* VecCopy() is smart, exists immediately if vguess == vdest */
318:     return 0;
319:   }
320:   if (*HH(it, it) != 0.0) {
321:     nrs[it] = *GRS(it) / *HH(it, it);
322:   } else {
324:     ksp->reason = KSP_DIVERGED_BREAKDOWN;

326:     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)));
327:     return 0;
328:   }
329:   for (ii = 1; ii <= it; ii++) {
330:     k  = it - ii;
331:     tt = *GRS(k);
332:     for (j = k + 1; j <= it; j++) tt = tt - *HH(k, j) * nrs[j];
333:     if (*HH(k, k) == 0.0) {
335:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
336:       PetscInfo(ksp, "Likely your matrix or preconditioner is singular. HH(k,k) is identically zero; k = %" PetscInt_FMT "\n", k);
337:       return 0;
338:     }
339:     nrs[k] = tt / *HH(k, k);
340:   }

342:   /* Accumulate the correction to the solution of the preconditioned problem in TEMP */
343:   VecSet(VEC_TEMP, 0.0);
344:   VecMAXPY(VEC_TEMP, it + 1, nrs, &VEC_VV(0));

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

361:   hh = HH(0, it);
362:   cc = CC(0);
363:   ss = SS(0);

365:   /* Apply all the previously computed plane rotations to the new column
366:      of the Hessenberg matrix */
367:   for (j = 1; j <= it; j++) {
368:     tt  = *hh;
369:     *hh = PetscConj(*cc) * tt + *ss * *(hh + 1);
370:     hh++;
371:     *hh = *cc++ * *hh - (*ss++ * tt);
372:   }

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

401:     *res = 0.0;
402:   }
403:   return 0;
404: }
405: /*
406:    This routine allocates more work vectors, starting from VEC_VV(it).
407:  */
408: PetscErrorCode KSPGMRESGetNewVectors(KSP ksp, PetscInt it)
409: {
410:   KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;
411:   PetscInt   nwork = gmres->nwork_alloc, k, nalloc;

413:   nalloc = PetscMin(ksp->max_it, gmres->delta_allocate);
414:   /* Adjust the number to allocate to make sure that we don't exceed the
415:     number of available slots */
416:   if (it + VEC_OFFSET + nalloc >= gmres->vecs_allocated) nalloc = gmres->vecs_allocated - it - VEC_OFFSET;
417:   if (!nalloc) return 0;

419:   gmres->vv_allocated += nalloc;

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

423:   gmres->mwork_alloc[nwork] = nalloc;
424:   for (k = 0; k < nalloc; k++) gmres->vecs[it + VEC_OFFSET + k] = gmres->user_work[nwork][k];
425:   gmres->nwork_alloc++;
426:   return 0;
427: }

429: PetscErrorCode KSPBuildSolution_GMRES(KSP ksp, Vec ptr, Vec *result)
430: {
431:   KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;

433:   if (!ptr) {
434:     if (!gmres->sol_temp) { VecDuplicate(ksp->vec_sol, &gmres->sol_temp); }
435:     ptr = gmres->sol_temp;
436:   }
437:   if (!gmres->nrs) {
438:     /* allocate the work area */
439:     PetscMalloc1(gmres->max_k, &gmres->nrs);
440:   }

442:   KSPGMRESBuildSoln(gmres->nrs, ksp->vec_sol, ptr, ksp, gmres->it);
443:   if (result) *result = ptr;
444:   return 0;
445: }

447: PetscErrorCode KSPView_GMRES(KSP ksp, PetscViewer viewer)
448: {
449:   KSP_GMRES  *gmres = (KSP_GMRES *)ksp->data;
450:   const char *cstr;
451:   PetscBool   iascii, isstring;

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

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

486:    Collective

488:    Input Parameters:
489: +  ksp - the `KSP` context
490: .  its - iteration number
491: .  fgnorm - 2-norm of residual (or gradient)
492: -  dummy - an collection of viewers created with `KSPViewerCreate()`

494:    Options Database Key:
495: .   -ksp_gmres_krylov_monitor <bool> - Plot the Krylov directions

497:    Level: intermediate

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

502: .seealso: [](chapter_ksp), `KSPGMRES`, `KSPMonitorSet()`, `KSPMonitorResidual()`, `VecView()`, `KSPViewersCreate()`, `KSPViewersDestroy()`
503: @*/
504: PetscErrorCode KSPGMRESMonitorKrylov(KSP ksp, PetscInt its, PetscReal fgnorm, void *dummy)
505: {
506:   PetscViewers viewers = (PetscViewers)dummy;
507:   KSP_GMRES   *gmres   = (KSP_GMRES *)ksp->data;
508:   Vec          x;
509:   PetscViewer  viewer;
510:   PetscBool    flg;

512:   PetscViewersGetViewer(viewers, gmres->it + 1, &viewer);
513:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &flg);
514:   if (!flg) {
515:     PetscViewerSetType(viewer, PETSCVIEWERDRAW);
516:     PetscViewerDrawSetInfo(viewer, NULL, "Krylov GMRES Monitor", PETSC_DECIDE, PETSC_DECIDE, 300, 300);
517:   }
518:   x = VEC_VV(gmres->it + 1);
519:   VecView(x, viewer);
520:   return 0;
521: }

523: PetscErrorCode KSPSetFromOptions_GMRES(KSP ksp, PetscOptionItems *PetscOptionsObject)
524: {
525:   PetscInt   restart;
526:   PetscReal  haptol, breakdowntol;
527:   KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;
528:   PetscBool  flg;

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

556: PetscErrorCode KSPGMRESSetHapTol_GMRES(KSP ksp, PetscReal tol)
557: {
558:   KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;

561:   gmres->haptol = tol;
562:   return 0;
563: }

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

569:   if (tol == PETSC_DEFAULT) {
570:     gmres->breakdowntol = 0.1;
571:     return 0;
572:   }
574:   gmres->breakdowntol = tol;
575:   return 0;
576: }

578: PetscErrorCode KSPGMRESGetRestart_GMRES(KSP ksp, PetscInt *max_k)
579: {
580:   KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;

582:   *max_k = gmres->max_k;
583:   return 0;
584: }

586: PetscErrorCode KSPGMRESSetRestart_GMRES(KSP ksp, PetscInt max_k)
587: {
588:   KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;

591:   if (!ksp->setupstage) {
592:     gmres->max_k = max_k;
593:   } else if (gmres->max_k != max_k) {
594:     gmres->max_k    = max_k;
595:     ksp->setupstage = KSP_SETUP_NEW;
596:     /* free the data structures, then create them again */
597:     KSPReset_GMRES(ksp);
598:   }
599:   return 0;
600: }

602: PetscErrorCode KSPGMRESSetOrthogonalization_GMRES(KSP ksp, FCN fcn)
603: {
604:   ((KSP_GMRES *)ksp->data)->orthog = fcn;
605:   return 0;
606: }

608: PetscErrorCode KSPGMRESGetOrthogonalization_GMRES(KSP ksp, FCN *fcn)
609: {
610:   *fcn = ((KSP_GMRES *)ksp->data)->orthog;
611:   return 0;
612: }

614: PetscErrorCode KSPGMRESSetPreAllocateVectors_GMRES(KSP ksp)
615: {
616:   KSP_GMRES *gmres;

618:   gmres                = (KSP_GMRES *)ksp->data;
619:   gmres->q_preallocate = 1;
620:   return 0;
621: }

623: PetscErrorCode KSPGMRESSetCGSRefinementType_GMRES(KSP ksp, KSPGMRESCGSRefinementType type)
624: {
625:   KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;

627:   gmres->cgstype = type;
628:   return 0;
629: }

631: PetscErrorCode KSPGMRESGetCGSRefinementType_GMRES(KSP ksp, KSPGMRESCGSRefinementType *type)
632: {
633:   KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;

635:   *type = gmres->cgstype;
636:   return 0;
637: }

639: /*@
640:    KSPGMRESSetCGSRefinementType - Sets the type of iterative refinement to use
641:          in the classical Gram Schmidt orthogonalization.

643:    Logically Collective

645:    Input Parameters:
646: +  ksp - the Krylov space context
647: -  type - the type of refinement
648: .vb
649:   KSP_GMRES_CGS_REFINE_NEVER
650:   KSP_GMRES_CGS_REFINE_IFNEEDED
651:   KSP_GMRES_CGS_REFINE_ALWAYS
652: .ve

654:   Options Database Key:
655: .  -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always> - refinement type

657:    Level: intermediate

659: .seealso: [](chapter_ksp), `KSPGMRES`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESCGSRefinementType`, `KSPGMRESClassicalGramSchmidtOrthogonalization()`, `KSPGMRESGetCGSRefinementType()`,
660:           `KSPGMRESGetOrthogonalization()`
661: @*/
662: PetscErrorCode KSPGMRESSetCGSRefinementType(KSP ksp, KSPGMRESCGSRefinementType type)
663: {
666:   PetscTryMethod(ksp, "KSPGMRESSetCGSRefinementType_C", (KSP, KSPGMRESCGSRefinementType), (ksp, type));
667:   return 0;
668: }

670: /*@
671:    KSPGMRESGetCGSRefinementType - Gets the type of iterative refinement to use
672:          in the classical Gram Schmidt orthogonalization.

674:    Not Collective

676:    Input Parameter:
677: .  ksp - the Krylov space context

679:    Output Parameter:
680: .  type - the type of refinement

682:   Options Database Key:
683: .  -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always> - type of refinement

685:    Level: intermediate

687: .seealso: [](chapter_ksp), `KSPGMRES`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESCGSRefinementType`, `KSPGMRESClassicalGramSchmidtOrthogonalization()`, `KSPGMRESSetCGSRefinementType()`,
688:           `KSPGMRESGetOrthogonalization()`
689: @*/
690: PetscErrorCode KSPGMRESGetCGSRefinementType(KSP ksp, KSPGMRESCGSRefinementType *type)
691: {
693:   PetscUseMethod(ksp, "KSPGMRESGetCGSRefinementType_C", (KSP, KSPGMRESCGSRefinementType *), (ksp, type));
694:   return 0;
695: }

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

700:    Logically Collective

702:    Input Parameters:
703: +  ksp - the Krylov space context
704: -  restart - integer restart value

706:   Options Database Key:
707: .  -ksp_gmres_restart <positive integer> - integer restart value

709:    Level: intermediate

711:     Note:
712:     The default value is 30.

714: .seealso: [](chapter_ksp), `KSPGMRES`, `KSPSetTolerances()`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESSetPreAllocateVectors()`, `KSPGMRESGetRestart()`
715: @*/
716: PetscErrorCode KSPGMRESSetRestart(KSP ksp, PetscInt restart)
717: {

720:   PetscTryMethod(ksp, "KSPGMRESSetRestart_C", (KSP, PetscInt), (ksp, restart));
721:   return 0;
722: }

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

727:    Not Collective

729:    Input Parameter:
730: .  ksp - the Krylov space context

732:    Output Parameter:
733: .   restart - integer restart value

735:    Level: intermediate

737: .seealso: [](chapter_ksp), `KSPGMRES`, `KSPSetTolerances()`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESSetPreAllocateVectors()`, `KSPGMRESSetRestart()`
738: @*/
739: PetscErrorCode KSPGMRESGetRestart(KSP ksp, PetscInt *restart)
740: {
741:   PetscUseMethod(ksp, "KSPGMRESGetRestart_C", (KSP, PetscInt *), (ksp, restart));
742:   return 0;
743: }

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

748:    Logically Collective

750:    Input Parameters:
751: +  ksp - the Krylov space context
752: -  tol - the tolerance

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

757:    Level: intermediate

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

764: .seealso: [](chapter_ksp), `KSPGMRES`, `KSPSetTolerances()`
765: @*/
766: PetscErrorCode KSPGMRESSetHapTol(KSP ksp, PetscReal tol)
767: {
769:   PetscTryMethod((ksp), "KSPGMRESSetHapTol_C", (KSP, PetscReal), ((ksp), (tol)));
770:   return 0;
771: }

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

776:    Logically Collective

778:    Input Parameters:
779: +  ksp - the Krylov space context
780: -  tol - the tolerance

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

785:    Level: intermediate

787:    Note:
788:    Divergence breakdown occurs when GMRES residual increases significantly during restart

790: .seealso: [](chapter_ksp), `KSPGMRES`, `KSPSetTolerances()`, `KSPGMRESSetHapTol()`
791: @*/
792: PetscErrorCode KSPGMRESSetBreakdownTolerance(KSP ksp, PetscReal tol)
793: {
795:   PetscTryMethod((ksp), "KSPGMRESSetBreakdownTolerance_C", (KSP, PetscReal), (ksp, tol));
796:   return 0;
797: }

799: /*MC
800:      KSPGMRES - Implements the Generalized Minimal Residual method [1] with restart

802:    Options Database Keys:
803: +   -ksp_gmres_restart <restart> - the number of Krylov directions to orthogonalize against
804: .   -ksp_gmres_haptol <tol> - sets the tolerance for "happy ending" (exact convergence)
805: .   -ksp_gmres_preallocate - preallocate all the Krylov search directions initially (otherwise groups of
806:                              vectors are allocated as needed)
807: .   -ksp_gmres_classicalgramschmidt - use classical (unmodified) Gram-Schmidt to orthogonalize against the Krylov space (fast) (the default)
808: .   -ksp_gmres_modifiedgramschmidt - use modified Gram-Schmidt in the orthogonalization (more stable, but slower)
809: .   -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always> - determine if iterative refinement is used to increase the
810:                                    stability of the classical Gram-Schmidt  orthogonalization.
811: -   -ksp_gmres_krylov_monitor - plot the Krylov space generated

813:    Level: beginner

815:    Note:
816:     Left and right preconditioning are supported, but not symmetric preconditioning.

818:    Reference:
819: .  [1] - YOUCEF SAAD AND MARTIN H. SCHULTZ, GMRES: A GENERALIZED MINIMAL RESIDUAL ALGORITHM FOR SOLVING NONSYMMETRIC LINEAR SYSTEMS.
820:           SIAM J. ScI. STAT. COMPUT. Vo|. 7, No. 3, July 1986.

822: .seealso: [](chapter_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPFGMRES`, `KSPLGMRES`,
823:           `KSPGMRESSetRestart()`, `KSPGMRESSetHapTol()`, `KSPGMRESSetPreAllocateVectors()`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESGetOrthogonalization()`,
824:           `KSPGMRESClassicalGramSchmidtOrthogonalization()`, `KSPGMRESModifiedGramSchmidtOrthogonalization()`,
825:           `KSPGMRESCGSRefinementType`, `KSPGMRESSetCGSRefinementType()`, `KSPGMRESGetCGSRefinementType()`, `KSPGMRESMonitorKrylov()`, `KSPSetPCSide()`
826: M*/

828: PETSC_EXTERN PetscErrorCode KSPCreate_GMRES(KSP ksp)
829: {
830:   KSP_GMRES *gmres;

832:   PetscNew(&gmres);
833:   ksp->data = (void *)gmres;

835:   KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 4);
836:   KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 3);
837:   KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_SYMMETRIC, 2);
838:   KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1);
839:   KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1);

841:   ksp->ops->buildsolution                = KSPBuildSolution_GMRES;
842:   ksp->ops->setup                        = KSPSetUp_GMRES;
843:   ksp->ops->solve                        = KSPSolve_GMRES;
844:   ksp->ops->reset                        = KSPReset_GMRES;
845:   ksp->ops->destroy                      = KSPDestroy_GMRES;
846:   ksp->ops->view                         = KSPView_GMRES;
847:   ksp->ops->setfromoptions               = KSPSetFromOptions_GMRES;
848:   ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_GMRES;
849:   ksp->ops->computeeigenvalues           = KSPComputeEigenvalues_GMRES;
850:   ksp->ops->computeritz                  = KSPComputeRitz_GMRES;
851:   PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetPreAllocateVectors_C", KSPGMRESSetPreAllocateVectors_GMRES);
852:   PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetOrthogonalization_C", KSPGMRESSetOrthogonalization_GMRES);
853:   PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetOrthogonalization_C", KSPGMRESGetOrthogonalization_GMRES);
854:   PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetRestart_C", KSPGMRESSetRestart_GMRES);
855:   PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetRestart_C", KSPGMRESGetRestart_GMRES);
856:   PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetHapTol_C", KSPGMRESSetHapTol_GMRES);
857:   PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetBreakdownTolerance_C", KSPGMRESSetBreakdownTolerance_GMRES);
858:   PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetCGSRefinementType_C", KSPGMRESSetCGSRefinementType_GMRES);
859:   PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetCGSRefinementType_C", KSPGMRESGetCGSRefinementType_GMRES);

861:   gmres->haptol         = 1.0e-30;
862:   gmres->breakdowntol   = 0.1;
863:   gmres->q_preallocate  = 0;
864:   gmres->delta_allocate = GMRES_DELTA_DIRECTIONS;
865:   gmres->orthog         = KSPGMRESClassicalGramSchmidtOrthogonalization;
866:   gmres->nrs            = NULL;
867:   gmres->sol_temp       = NULL;
868:   gmres->max_k          = GMRES_DEFAULT_MAXK;
869:   gmres->Rsvd           = NULL;
870:   gmres->cgstype        = KSP_GMRES_CGS_REFINE_NEVER;
871:   gmres->orthogwork     = NULL;
872:   return 0;
873: }