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: }