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: 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", (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: /* check for the convergence */
124: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
125: ksp->rnorm = res;
126: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
127: gmres->it = (it - 1);
128: PetscCall(KSPLogResidualHistory(ksp, res));
129: PetscCall(KSPLogErrorHistory(ksp));
130: PetscCall(KSPMonitor(ksp, ksp->its, res));
131: if (!res) {
132: ksp->reason = KSP_CONVERGED_ATOL;
133: PetscCall(PetscInfo(ksp, "Converged due to zero residual norm on entry\n"));
134: PetscFunctionReturn(PETSC_SUCCESS);
135: }
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, "You reached the 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: /* Monitor if we know that we will not return for a restart */
190: if (it && (ksp->reason || ksp->its >= ksp->max_it)) {
191: PetscCall(KSPLogResidualHistory(ksp, res));
192: PetscCall(KSPLogErrorHistory(ksp));
193: PetscCall(KSPMonitor(ksp, ksp->its, res));
194: }
196: if (itcount) *itcount = it;
198: /*
199: Down here we have to solve for the "best" coefficients of the Krylov
200: columns, add the solution values together, and possibly unwind the
201: preconditioning from the solution
202: */
203: /* Form the solution (or the solution so far) */
204: PetscCall(KSPGMRESBuildSoln(GRS(0), ksp->vec_sol, ksp->vec_sol, ksp, it - 1));
205: PetscFunctionReturn(PETSC_SUCCESS);
206: }
208: 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(VecSet(VEC_TEMP, 0.0));
349: PetscCall(VecMAXPY(VEC_TEMP, it + 1, nrs, &VEC_VV(0)));
351: PetscCall(KSPUnwindPreconditioner(ksp, VEC_TEMP, VEC_TEMP_MATOP));
352: /* add solution to previous solution */
353: if (vdest != vs) PetscCall(VecCopy(vs, vdest));
354: PetscCall(VecAXPY(vdest, 1.0, VEC_TEMP));
355: PetscFunctionReturn(PETSC_SUCCESS);
356: }
357: /*
358: Do the scalar work for the orthogonalization. Return new residual norm.
359: */
360: static PetscErrorCode KSPGMRESUpdateHessenberg(KSP ksp, PetscInt it, PetscBool hapend, PetscReal *res)
361: {
362: PetscScalar *hh, *cc, *ss, tt;
363: PetscInt j;
364: KSP_GMRES *gmres = (KSP_GMRES *)(ksp->data);
366: PetscFunctionBegin;
367: hh = HH(0, it);
368: cc = CC(0);
369: ss = SS(0);
371: /* Apply all the previously computed plane rotations to the new column
372: of the Hessenberg matrix */
373: for (j = 1; j <= it; j++) {
374: tt = *hh;
375: *hh = PetscConj(*cc) * tt + *ss * *(hh + 1);
376: hh++;
377: *hh = *cc++ * *hh - (*ss++ * tt);
378: }
380: /*
381: compute the new plane rotation, and apply it to:
382: 1) the right-hand-side of the Hessenberg system
383: 2) the new column of the Hessenberg matrix
384: thus obtaining the updated value of the residual
385: */
386: if (!hapend) {
387: tt = PetscSqrtScalar(PetscConj(*hh) * *hh + PetscConj(*(hh + 1)) * *(hh + 1));
388: if (tt == 0.0) {
389: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "tt == 0.0");
390: ksp->reason = KSP_DIVERGED_NULL;
391: PetscFunctionReturn(PETSC_SUCCESS);
392: }
393: *cc = *hh / tt;
394: *ss = *(hh + 1) / tt;
395: *GRS(it + 1) = -(*ss * *GRS(it));
396: *GRS(it) = PetscConj(*cc) * *GRS(it);
397: *hh = PetscConj(*cc) * *hh + *ss * *(hh + 1);
398: *res = PetscAbsScalar(*GRS(it + 1));
399: } else {
400: /* happy breakdown: HH(it+1, it) = 0, therefore we don't need to apply
401: another rotation matrix (so RH doesn't change). The new residual is
402: always the new sine term times the residual from last time (GRS(it)),
403: but now the new sine rotation would be zero...so the residual should
404: be zero...so we will multiply "zero" by the last residual. This might
405: not be exactly what we want to do here -could just return "zero". */
407: *res = 0.0;
408: }
409: PetscFunctionReturn(PETSC_SUCCESS);
410: }
411: /*
412: This routine allocates more work vectors, starting from VEC_VV(it).
413: */
414: PetscErrorCode KSPGMRESGetNewVectors(KSP ksp, PetscInt it)
415: {
416: KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;
417: PetscInt nwork = gmres->nwork_alloc, k, nalloc;
419: PetscFunctionBegin;
420: nalloc = PetscMin(ksp->max_it, gmres->delta_allocate);
421: /* Adjust the number to allocate to make sure that we don't exceed the
422: number of available slots */
423: if (it + VEC_OFFSET + nalloc >= gmres->vecs_allocated) nalloc = gmres->vecs_allocated - it - VEC_OFFSET;
424: if (!nalloc) PetscFunctionReturn(PETSC_SUCCESS);
426: gmres->vv_allocated += nalloc;
428: PetscCall(KSPCreateVecs(ksp, nalloc, &gmres->user_work[nwork], 0, NULL));
430: gmres->mwork_alloc[nwork] = nalloc;
431: for (k = 0; k < nalloc; k++) gmres->vecs[it + VEC_OFFSET + k] = gmres->user_work[nwork][k];
432: gmres->nwork_alloc++;
433: PetscFunctionReturn(PETSC_SUCCESS);
434: }
436: PetscErrorCode KSPBuildSolution_GMRES(KSP ksp, Vec ptr, Vec *result)
437: {
438: KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;
440: PetscFunctionBegin;
441: if (!ptr) {
442: if (!gmres->sol_temp) PetscCall(VecDuplicate(ksp->vec_sol, &gmres->sol_temp));
443: ptr = gmres->sol_temp;
444: }
445: if (!gmres->nrs) {
446: /* allocate the work area */
447: PetscCall(PetscMalloc1(gmres->max_k, &gmres->nrs));
448: }
450: PetscCall(KSPGMRESBuildSoln(gmres->nrs, ksp->vec_sol, ptr, ksp, gmres->it));
451: if (result) *result = ptr;
452: PetscFunctionReturn(PETSC_SUCCESS);
453: }
455: PetscErrorCode KSPView_GMRES(KSP ksp, PetscViewer viewer)
456: {
457: KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;
458: const char *cstr;
459: PetscBool iascii, isstring;
461: PetscFunctionBegin;
462: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
463: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
464: if (gmres->orthog == KSPGMRESClassicalGramSchmidtOrthogonalization) {
465: switch (gmres->cgstype) {
466: case (KSP_GMRES_CGS_REFINE_NEVER):
467: cstr = "Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement";
468: break;
469: case (KSP_GMRES_CGS_REFINE_ALWAYS):
470: cstr = "Classical (unmodified) Gram-Schmidt Orthogonalization with one step of iterative refinement";
471: break;
472: case (KSP_GMRES_CGS_REFINE_IFNEEDED):
473: cstr = "Classical (unmodified) Gram-Schmidt Orthogonalization with one step of iterative refinement when needed";
474: break;
475: default:
476: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Unknown orthogonalization");
477: }
478: } else if (gmres->orthog == KSPGMRESModifiedGramSchmidtOrthogonalization) {
479: cstr = "Modified Gram-Schmidt Orthogonalization";
480: } else {
481: cstr = "unknown orthogonalization";
482: }
483: if (iascii) {
484: PetscCall(PetscViewerASCIIPrintf(viewer, " restart=%" PetscInt_FMT ", using %s\n", gmres->max_k, cstr));
485: PetscCall(PetscViewerASCIIPrintf(viewer, " happy breakdown tolerance %g\n", (double)gmres->haptol));
486: } else if (isstring) {
487: PetscCall(PetscViewerStringSPrintf(viewer, "%s restart %" PetscInt_FMT, cstr, gmres->max_k));
488: }
489: PetscFunctionReturn(PETSC_SUCCESS);
490: }
492: /*@C
493: KSPGMRESMonitorKrylov - Calls `VecView()` for each new direction in the `KSPGMRES` accumulated Krylov space.
495: Collective
497: Input Parameters:
498: + ksp - the `KSP` context
499: . its - iteration number
500: . fgnorm - 2-norm of residual (or gradient)
501: - dummy - a collection of viewers created with `PetscViewersCreate()`
503: Options Database Key:
504: . -ksp_gmres_krylov_monitor <bool> - Plot the Krylov directions
506: Level: intermediate
508: Note:
509: A new `PETSCVIEWERDRAW` is created for each Krylov vector so they can all be simultaneously viewed
511: .seealso: [](ch_ksp), `KSPGMRES`, `KSPMonitorSet()`, `KSPMonitorResidual()`, `VecView()`, `PetscViewersCreate()`, `PetscViewersDestroy()`
512: @*/
513: PetscErrorCode KSPGMRESMonitorKrylov(KSP ksp, PetscInt its, PetscReal fgnorm, void *dummy)
514: {
515: PetscViewers viewers = (PetscViewers)dummy;
516: KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;
517: Vec x;
518: PetscViewer viewer;
519: PetscBool flg;
521: PetscFunctionBegin;
522: PetscCall(PetscViewersGetViewer(viewers, gmres->it + 1, &viewer));
523: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &flg));
524: if (!flg) {
525: PetscCall(PetscViewerSetType(viewer, PETSCVIEWERDRAW));
526: PetscCall(PetscViewerDrawSetInfo(viewer, NULL, "Krylov GMRES Monitor", PETSC_DECIDE, PETSC_DECIDE, 300, 300));
527: }
528: x = VEC_VV(gmres->it + 1);
529: PetscCall(VecView(x, viewer));
530: PetscFunctionReturn(PETSC_SUCCESS);
531: }
533: PetscErrorCode KSPSetFromOptions_GMRES(KSP ksp, PetscOptionItems *PetscOptionsObject)
534: {
535: PetscInt restart;
536: PetscReal haptol, breakdowntol;
537: KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;
538: PetscBool flg;
540: PetscFunctionBegin;
541: PetscOptionsHeadBegin(PetscOptionsObject, "KSP GMRES Options");
542: PetscCall(PetscOptionsInt("-ksp_gmres_restart", "Number of Krylov search directions", "KSPGMRESSetRestart", gmres->max_k, &restart, &flg));
543: if (flg) PetscCall(KSPGMRESSetRestart(ksp, restart));
544: PetscCall(PetscOptionsReal("-ksp_gmres_haptol", "Tolerance for exact convergence (happy ending)", "KSPGMRESSetHapTol", gmres->haptol, &haptol, &flg));
545: if (flg) PetscCall(KSPGMRESSetHapTol(ksp, haptol));
546: PetscCall(PetscOptionsReal("-ksp_gmres_breakdown_tolerance", "Divergence breakdown tolerance during GMRES restart", "KSPGMRESSetBreakdownTolerance", gmres->breakdowntol, &breakdowntol, &flg));
547: if (flg) PetscCall(KSPGMRESSetBreakdownTolerance(ksp, breakdowntol));
548: flg = PETSC_FALSE;
549: PetscCall(PetscOptionsBool("-ksp_gmres_preallocate", "Preallocate Krylov vectors", "KSPGMRESSetPreAllocateVectors", flg, &flg, NULL));
550: if (flg) PetscCall(KSPGMRESSetPreAllocateVectors(ksp));
551: PetscCall(PetscOptionsBoolGroupBegin("-ksp_gmres_classicalgramschmidt", "Classical (unmodified) Gram-Schmidt (fast)", "KSPGMRESSetOrthogonalization", &flg));
552: if (flg) PetscCall(KSPGMRESSetOrthogonalization(ksp, KSPGMRESClassicalGramSchmidtOrthogonalization));
553: PetscCall(PetscOptionsBoolGroupEnd("-ksp_gmres_modifiedgramschmidt", "Modified Gram-Schmidt (slow,more stable)", "KSPGMRESSetOrthogonalization", &flg));
554: if (flg) PetscCall(KSPGMRESSetOrthogonalization(ksp, KSPGMRESModifiedGramSchmidtOrthogonalization));
555: 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));
556: flg = PETSC_FALSE;
557: PetscCall(PetscOptionsBool("-ksp_gmres_krylov_monitor", "Plot the Krylov directions", "KSPMonitorSet", flg, &flg, NULL));
558: if (flg) {
559: PetscViewers viewers;
560: PetscCall(PetscViewersCreate(PetscObjectComm((PetscObject)ksp), &viewers));
561: PetscCall(KSPMonitorSet(ksp, KSPGMRESMonitorKrylov, viewers, (PetscErrorCode(*)(void **))PetscViewersDestroy));
562: }
563: PetscOptionsHeadEnd();
564: PetscFunctionReturn(PETSC_SUCCESS);
565: }
567: PetscErrorCode KSPGMRESSetHapTol_GMRES(KSP ksp, PetscReal tol)
568: {
569: KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;
571: PetscFunctionBegin;
572: PetscCheck(tol >= 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Tolerance must be non-negative");
573: gmres->haptol = tol;
574: PetscFunctionReturn(PETSC_SUCCESS);
575: }
577: PetscErrorCode KSPGMRESSetBreakdownTolerance_GMRES(KSP ksp, PetscReal tol)
578: {
579: KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;
581: PetscFunctionBegin;
582: if (tol == (PetscReal)PETSC_DEFAULT) {
583: gmres->breakdowntol = 0.1;
584: PetscFunctionReturn(PETSC_SUCCESS);
585: }
586: PetscCheck(tol >= 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Breakdown tolerance must be non-negative");
587: gmres->breakdowntol = tol;
588: PetscFunctionReturn(PETSC_SUCCESS);
589: }
591: PetscErrorCode KSPGMRESGetRestart_GMRES(KSP ksp, PetscInt *max_k)
592: {
593: KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;
595: PetscFunctionBegin;
596: *max_k = gmres->max_k;
597: PetscFunctionReturn(PETSC_SUCCESS);
598: }
600: PetscErrorCode KSPGMRESSetRestart_GMRES(KSP ksp, PetscInt max_k)
601: {
602: KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;
604: PetscFunctionBegin;
605: PetscCheck(max_k >= 1, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Restart must be positive");
606: if (!ksp->setupstage) {
607: gmres->max_k = max_k;
608: } else if (gmres->max_k != max_k) {
609: gmres->max_k = max_k;
610: ksp->setupstage = KSP_SETUP_NEW;
611: /* free the data structures, then create them again */
612: PetscCall(KSPReset_GMRES(ksp));
613: }
614: PetscFunctionReturn(PETSC_SUCCESS);
615: }
617: PetscErrorCode KSPGMRESSetOrthogonalization_GMRES(KSP ksp, FCN fcn)
618: {
619: PetscFunctionBegin;
620: ((KSP_GMRES *)ksp->data)->orthog = fcn;
621: PetscFunctionReturn(PETSC_SUCCESS);
622: }
624: PetscErrorCode KSPGMRESGetOrthogonalization_GMRES(KSP ksp, FCN *fcn)
625: {
626: PetscFunctionBegin;
627: *fcn = ((KSP_GMRES *)ksp->data)->orthog;
628: PetscFunctionReturn(PETSC_SUCCESS);
629: }
631: PetscErrorCode KSPGMRESSetPreAllocateVectors_GMRES(KSP ksp)
632: {
633: KSP_GMRES *gmres;
635: PetscFunctionBegin;
636: gmres = (KSP_GMRES *)ksp->data;
637: gmres->q_preallocate = 1;
638: PetscFunctionReturn(PETSC_SUCCESS);
639: }
641: PetscErrorCode KSPGMRESSetCGSRefinementType_GMRES(KSP ksp, KSPGMRESCGSRefinementType type)
642: {
643: KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;
645: PetscFunctionBegin;
646: gmres->cgstype = type;
647: PetscFunctionReturn(PETSC_SUCCESS);
648: }
650: PetscErrorCode KSPGMRESGetCGSRefinementType_GMRES(KSP ksp, KSPGMRESCGSRefinementType *type)
651: {
652: KSP_GMRES *gmres = (KSP_GMRES *)ksp->data;
654: PetscFunctionBegin;
655: *type = gmres->cgstype;
656: PetscFunctionReturn(PETSC_SUCCESS);
657: }
659: /*@
660: KSPGMRESSetCGSRefinementType - Sets the type of iterative refinement to use
661: in the classical Gram Schmidt orthogonalization.
663: Logically Collective
665: Input Parameters:
666: + ksp - the Krylov space context
667: - type - the type of refinement
668: .vb
669: KSP_GMRES_CGS_REFINE_NEVER
670: KSP_GMRES_CGS_REFINE_IFNEEDED
671: KSP_GMRES_CGS_REFINE_ALWAYS
672: .ve
674: Options Database Key:
675: . -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always> - refinement type
677: Level: intermediate
679: .seealso: [](ch_ksp), `KSPGMRES`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESCGSRefinementType`, `KSPGMRESClassicalGramSchmidtOrthogonalization()`, `KSPGMRESGetCGSRefinementType()`,
680: `KSPGMRESGetOrthogonalization()`
681: @*/
682: PetscErrorCode KSPGMRESSetCGSRefinementType(KSP ksp, KSPGMRESCGSRefinementType type)
683: {
684: PetscFunctionBegin;
687: PetscTryMethod(ksp, "KSPGMRESSetCGSRefinementType_C", (KSP, KSPGMRESCGSRefinementType), (ksp, type));
688: PetscFunctionReturn(PETSC_SUCCESS);
689: }
691: /*@
692: KSPGMRESGetCGSRefinementType - Gets the type of iterative refinement to use
693: in the classical Gram Schmidt orthogonalization.
695: Not Collective
697: Input Parameter:
698: . ksp - the Krylov space context
700: Output Parameter:
701: . type - the type of refinement
703: Options Database Key:
704: . -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always> - type of refinement
706: Level: intermediate
708: .seealso: [](ch_ksp), `KSPGMRES`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESCGSRefinementType`, `KSPGMRESClassicalGramSchmidtOrthogonalization()`, `KSPGMRESSetCGSRefinementType()`,
709: `KSPGMRESGetOrthogonalization()`
710: @*/
711: PetscErrorCode KSPGMRESGetCGSRefinementType(KSP ksp, KSPGMRESCGSRefinementType *type)
712: {
713: PetscFunctionBegin;
715: PetscUseMethod(ksp, "KSPGMRESGetCGSRefinementType_C", (KSP, KSPGMRESCGSRefinementType *), (ksp, type));
716: PetscFunctionReturn(PETSC_SUCCESS);
717: }
719: /*@
720: KSPGMRESSetRestart - Sets number of iterations at which `KSPGMRES`, `KSPFGMRES` and `KSPLGMRES` restarts.
722: Logically Collective
724: Input Parameters:
725: + ksp - the Krylov space context
726: - restart - integer restart value
728: Options Database Key:
729: . -ksp_gmres_restart <positive integer> - integer restart value
731: Level: intermediate
733: Note:
734: The default value is 30.
736: .seealso: [](ch_ksp), `KSPGMRES`, `KSPSetTolerances()`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESSetPreAllocateVectors()`, `KSPGMRESGetRestart()`
737: @*/
738: PetscErrorCode KSPGMRESSetRestart(KSP ksp, PetscInt restart)
739: {
740: PetscFunctionBegin;
743: PetscTryMethod(ksp, "KSPGMRESSetRestart_C", (KSP, PetscInt), (ksp, restart));
744: PetscFunctionReturn(PETSC_SUCCESS);
745: }
747: /*@
748: KSPGMRESGetRestart - Gets number of iterations at which `KSPGMRES`, `KSPFGMRES` and `KSPLGMRES` restarts.
750: Not Collective
752: Input Parameter:
753: . ksp - the Krylov space context
755: Output Parameter:
756: . restart - integer restart value
758: Level: intermediate
760: .seealso: [](ch_ksp), `KSPGMRES`, `KSPSetTolerances()`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESSetPreAllocateVectors()`, `KSPGMRESSetRestart()`
761: @*/
762: PetscErrorCode KSPGMRESGetRestart(KSP ksp, PetscInt *restart)
763: {
764: PetscFunctionBegin;
765: PetscUseMethod(ksp, "KSPGMRESGetRestart_C", (KSP, PetscInt *), (ksp, restart));
766: PetscFunctionReturn(PETSC_SUCCESS);
767: }
769: /*@
770: KSPGMRESSetHapTol - Sets tolerance for determining happy breakdown in `KSPGMRES`, `KSPFGMRES` and `KSPLGMRES`
772: Logically Collective
774: Input Parameters:
775: + ksp - the Krylov space context
776: - tol - the tolerance
778: Options Database Key:
779: . -ksp_gmres_haptol <positive real value> - set tolerance for determining happy breakdown
781: Level: intermediate
783: Note:
784: Happy breakdown is the rare case in `KSPGMRES` where an 'exact' solution is obtained after
785: a certain number of iterations. If you attempt more iterations after this point unstable
786: things can happen hence very occasionally you may need to set this value to detect this condition
788: .seealso: [](ch_ksp), `KSPGMRES`, `KSPSetTolerances()`
789: @*/
790: PetscErrorCode KSPGMRESSetHapTol(KSP ksp, PetscReal tol)
791: {
792: PetscFunctionBegin;
794: PetscTryMethod((ksp), "KSPGMRESSetHapTol_C", (KSP, PetscReal), ((ksp), (tol)));
795: PetscFunctionReturn(PETSC_SUCCESS);
796: }
798: /*@
799: KSPGMRESSetBreakdownTolerance - Sets tolerance for determining divergence breakdown in `KSPGMRES`.
801: Logically Collective
803: Input Parameters:
804: + ksp - the Krylov space context
805: - tol - the tolerance
807: Options Database Key:
808: . -ksp_gmres_breakdown_tolerance <positive real value> - set tolerance for determining divergence breakdown
810: Level: intermediate
812: Note:
813: Divergence breakdown occurs when GMRES residual increases significantly during restart
815: .seealso: [](ch_ksp), `KSPGMRES`, `KSPSetTolerances()`, `KSPGMRESSetHapTol()`
816: @*/
817: PetscErrorCode KSPGMRESSetBreakdownTolerance(KSP ksp, PetscReal tol)
818: {
819: PetscFunctionBegin;
821: PetscTryMethod((ksp), "KSPGMRESSetBreakdownTolerance_C", (KSP, PetscReal), (ksp, tol));
822: PetscFunctionReturn(PETSC_SUCCESS);
823: }
825: /*MC
826: KSPGMRES - Implements the Generalized Minimal Residual method [1] with restart
828: Options Database Keys:
829: + -ksp_gmres_restart <restart> - the number of Krylov directions to orthogonalize against
830: . -ksp_gmres_haptol <tol> - sets the tolerance for "happy ending" (exact convergence)
831: . -ksp_gmres_preallocate - preallocate all the Krylov search directions initially (otherwise groups of
832: vectors are allocated as needed)
833: . -ksp_gmres_classicalgramschmidt - use classical (unmodified) Gram-Schmidt to orthogonalize against the Krylov space (fast) (the default)
834: . -ksp_gmres_modifiedgramschmidt - use modified Gram-Schmidt in the orthogonalization (more stable, but slower)
835: . -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always> - determine if iterative refinement is used to increase the
836: stability of the classical Gram-Schmidt orthogonalization.
837: - -ksp_gmres_krylov_monitor - plot the Krylov space generated
839: Level: beginner
841: Note:
842: Left and right preconditioning are supported, but not symmetric preconditioning.
844: Reference:
845: . [1] - YOUCEF SAAD AND MARTIN H. SCHULTZ, GMRES: A GENERALIZED MINIMAL RESIDUAL ALGORITHM FOR SOLVING NONSYMMETRIC LINEAR SYSTEMS.
846: SIAM J. ScI. STAT. COMPUT. Vo|. 7, No. 3, July 1986.
848: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPFGMRES`, `KSPLGMRES`,
849: `KSPGMRESSetRestart()`, `KSPGMRESSetHapTol()`, `KSPGMRESSetPreAllocateVectors()`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESGetOrthogonalization()`,
850: `KSPGMRESClassicalGramSchmidtOrthogonalization()`, `KSPGMRESModifiedGramSchmidtOrthogonalization()`,
851: `KSPGMRESCGSRefinementType`, `KSPGMRESSetCGSRefinementType()`, `KSPGMRESGetCGSRefinementType()`, `KSPGMRESMonitorKrylov()`, `KSPSetPCSide()`
852: M*/
854: PETSC_EXTERN PetscErrorCode KSPCreate_GMRES(KSP ksp)
855: {
856: KSP_GMRES *gmres;
858: PetscFunctionBegin;
859: PetscCall(PetscNew(&gmres));
860: ksp->data = (void *)gmres;
862: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 4));
863: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 3));
864: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_SYMMETRIC, 2));
865: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));
866: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
868: ksp->ops->buildsolution = KSPBuildSolution_GMRES;
869: ksp->ops->setup = KSPSetUp_GMRES;
870: ksp->ops->solve = KSPSolve_GMRES;
871: ksp->ops->reset = KSPReset_GMRES;
872: ksp->ops->destroy = KSPDestroy_GMRES;
873: ksp->ops->view = KSPView_GMRES;
874: ksp->ops->setfromoptions = KSPSetFromOptions_GMRES;
875: ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_GMRES;
876: ksp->ops->computeeigenvalues = KSPComputeEigenvalues_GMRES;
877: ksp->ops->computeritz = KSPComputeRitz_GMRES;
878: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetPreAllocateVectors_C", KSPGMRESSetPreAllocateVectors_GMRES));
879: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetOrthogonalization_C", KSPGMRESSetOrthogonalization_GMRES));
880: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetOrthogonalization_C", KSPGMRESGetOrthogonalization_GMRES));
881: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetRestart_C", KSPGMRESSetRestart_GMRES));
882: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetRestart_C", KSPGMRESGetRestart_GMRES));
883: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetHapTol_C", KSPGMRESSetHapTol_GMRES));
884: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetBreakdownTolerance_C", KSPGMRESSetBreakdownTolerance_GMRES));
885: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetCGSRefinementType_C", KSPGMRESSetCGSRefinementType_GMRES));
886: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetCGSRefinementType_C", KSPGMRESGetCGSRefinementType_GMRES));
888: gmres->haptol = 1.0e-30;
889: gmres->breakdowntol = 0.1;
890: gmres->q_preallocate = 0;
891: gmres->delta_allocate = GMRES_DELTA_DIRECTIONS;
892: gmres->orthog = KSPGMRESClassicalGramSchmidtOrthogonalization;
893: gmres->nrs = NULL;
894: gmres->sol_temp = NULL;
895: gmres->max_k = GMRES_DEFAULT_MAXK;
896: gmres->Rsvd = NULL;
897: gmres->cgstype = KSP_GMRES_CGS_REFINE_NEVER;
898: gmres->orthogwork = NULL;
899: PetscFunctionReturn(PETSC_SUCCESS);
900: }