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);
49: PetscLogObjectMemory((PetscObject)ksp,(hh + hes + rs + 2*cc)*sizeof(PetscScalar));
51: if (ksp->calc_sings) {
52: /* Allocate workspace to hold Hessenberg matrix needed by lapack */
53: PetscMalloc1((max_k + 3)*(max_k + 9),&gmres->Rsvd);
54: PetscLogObjectMemory((PetscObject)ksp,(max_k + 3)*(max_k + 9)*sizeof(PetscScalar));
55: PetscMalloc1(6*(max_k+2),&gmres->Dsvd);
56: PetscLogObjectMemory((PetscObject)ksp,6*(max_k+2)*sizeof(PetscReal));
57: }
59: /* Allocate array to hold pointers to user vectors. Note that we need
60: 4 + max_k + 1 (since we need it+1 vectors, and it <= max_k) */
61: gmres->vecs_allocated = VEC_OFFSET + 2 + max_k + gmres->nextra_vecs;
63: PetscMalloc1(gmres->vecs_allocated,&gmres->vecs);
64: PetscMalloc1(VEC_OFFSET+2+max_k,&gmres->user_work);
65: PetscMalloc1(VEC_OFFSET+2+max_k,&gmres->mwork_alloc);
66: PetscLogObjectMemory((PetscObject)ksp,(VEC_OFFSET+2+max_k)*(sizeof(Vec*)+sizeof(PetscInt)) + gmres->vecs_allocated*sizeof(Vec));
68: if (gmres->q_preallocate) {
69: gmres->vv_allocated = VEC_OFFSET + 2 + max_k;
71: KSPCreateVecs(ksp,gmres->vv_allocated,&gmres->user_work[0],0,NULL);
72: PetscLogObjectParents(ksp,gmres->vv_allocated,gmres->user_work[0]);
74: gmres->mwork_alloc[0] = gmres->vv_allocated;
75: gmres->nwork_alloc = 1;
76: for (k=0; k<gmres->vv_allocated; k++) {
77: gmres->vecs[k] = gmres->user_work[0][k];
78: }
79: } else {
80: gmres->vv_allocated = 5;
82: KSPCreateVecs(ksp,5,&gmres->user_work[0],0,NULL);
83: PetscLogObjectParents(ksp,5,gmres->user_work[0]);
85: gmres->mwork_alloc[0] = 5;
86: gmres->nwork_alloc = 1;
87: for (k=0; k<gmres->vv_allocated; k++) {
88: gmres->vecs[k] = gmres->user_work[0][k];
89: }
90: }
91: return 0;
92: }
94: /*
95: Run gmres, possibly with restart. Return residual history if requested.
96: input parameters:
98: . gmres - structure containing parameters and work areas
100: output parameters:
101: . nres - residuals (from preconditioned system) at each step.
102: If restarting, consider passing nres+it. If null,
103: ignored
104: . itcount - number of iterations used. nres[0] to nres[itcount]
105: are defined. If null, ignored.
107: Notes:
108: On entry, the value in vector VEC_VV(0) should be the initial residual
109: (this allows shortcuts where the initial preconditioned residual is 0).
110: */
111: PetscErrorCode KSPGMRESCycle(PetscInt *itcount,KSP ksp)
112: {
113: KSP_GMRES *gmres = (KSP_GMRES*)(ksp->data);
114: PetscReal res,hapbnd,tt;
115: PetscInt it = 0, max_k = gmres->max_k;
116: PetscBool hapend = PETSC_FALSE;
118: if (itcount) *itcount = 0;
119: VecNormalize(VEC_VV(0),&res);
120: KSPCheckNorm(ksp,res);
122: /* the constant .1 is arbitrary, just some measure at how incorrect the residuals are */
123: if ((ksp->rnorm > 0.0) && (PetscAbsReal(res-ksp->rnorm) > gmres->breakdowntol*gmres->rnorm0)) {
125: else {
126: 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);
127: ksp->reason = KSP_DIVERGED_BREAKDOWN;
128: return 0;
129: }
130: }
131: *GRS(0) = gmres->rnorm0 = res;
133: /* check for the convergence */
134: PetscObjectSAWsTakeAccess((PetscObject)ksp);
135: ksp->rnorm = res;
136: PetscObjectSAWsGrantAccess((PetscObject)ksp);
137: gmres->it = (it - 1);
138: KSPLogResidualHistory(ksp,res);
139: KSPLogErrorHistory(ksp);
140: KSPMonitor(ksp,ksp->its,res);
141: if (!res) {
142: ksp->reason = KSP_CONVERGED_ATOL;
143: PetscInfo(ksp,"Converged due to zero residual norm on entry\n");
144: return 0;
145: }
147: (*ksp->converged)(ksp,ksp->its,res,&ksp->reason,ksp->cnvP);
148: while (!ksp->reason && it < max_k && ksp->its < ksp->max_it) {
149: if (it) {
150: KSPLogResidualHistory(ksp,res);
151: KSPLogErrorHistory(ksp);
152: KSPMonitor(ksp,ksp->its,res);
153: }
154: gmres->it = (it - 1);
155: if (gmres->vv_allocated <= it + VEC_OFFSET + 1) {
156: KSPGMRESGetNewVectors(ksp,it+1);
157: }
158: KSP_PCApplyBAorAB(ksp,VEC_VV(it),VEC_VV(1+it),VEC_TEMP_MATOP);
160: /* update hessenberg matrix and do Gram-Schmidt */
161: (*gmres->orthog)(ksp,it);
162: if (ksp->reason) break;
164: /* vv(i+1) . vv(i+1) */
165: VecNormalize(VEC_VV(it+1),&tt);
166: KSPCheckNorm(ksp,tt);
168: /* save the magnitude */
169: *HH(it+1,it) = tt;
170: *HES(it+1,it) = tt;
172: /* check for the happy breakdown */
173: hapbnd = PetscAbsScalar(tt / *GRS(it));
174: if (hapbnd > gmres->haptol) hapbnd = gmres->haptol;
175: if (tt < hapbnd) {
176: PetscInfo(ksp,"Detected happy breakdown, current hapbnd = %14.12e tt = %14.12e\n",(double)hapbnd,(double)tt);
177: hapend = PETSC_TRUE;
178: }
179: KSPGMRESUpdateHessenberg(ksp,it,hapend,&res);
181: it++;
182: gmres->it = (it-1); /* For converged */
183: ksp->its++;
184: ksp->rnorm = res;
185: if (ksp->reason) break;
187: (*ksp->converged)(ksp,ksp->its,res,&ksp->reason,ksp->cnvP);
189: /* Catch error in happy breakdown and signal convergence and break from loop */
190: if (hapend) {
191: if (ksp->normtype == KSP_NORM_NONE) { /* convergence test was skipped in this case */
192: ksp->reason = KSP_CONVERGED_HAPPY_BREAKDOWN;
193: } else if (!ksp->reason) {
195: else {
196: ksp->reason = KSP_DIVERGED_BREAKDOWN;
197: break;
198: }
199: }
200: }
201: }
203: /* Monitor if we know that we will not return for a restart */
204: if (it && (ksp->reason || ksp->its >= ksp->max_it)) {
205: KSPLogResidualHistory(ksp,res);
206: KSPLogErrorHistory(ksp);
207: KSPMonitor(ksp,ksp->its,res);
208: }
210: if (itcount) *itcount = it;
212: /*
213: Down here we have to solve for the "best" coefficients of the Krylov
214: columns, add the solution values together, and possibly unwind the
215: preconditioning from the solution
216: */
217: /* Form the solution (or the solution so far) */
218: KSPGMRESBuildSoln(GRS(0),ksp->vec_sol,ksp->vec_sol,ksp,it-1);
219: return 0;
220: }
222: PetscErrorCode KSPSolve_GMRES(KSP ksp)
223: {
224: PetscInt its,itcount,i;
225: KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
226: PetscBool guess_zero = ksp->guess_zero;
227: PetscInt N = gmres->max_k + 1;
231: PetscObjectSAWsTakeAccess((PetscObject)ksp);
232: ksp->its = 0;
233: PetscObjectSAWsGrantAccess((PetscObject)ksp);
235: itcount = 0;
236: gmres->fullcycle = 0;
237: ksp->rnorm = -1.0; /* special marker for KSPGMRESCycle() */
238: while (!ksp->reason || (ksp->rnorm == -1 && ksp->reason == KSP_DIVERGED_PC_FAILED)) {
239: KSPInitialResidual(ksp,ksp->vec_sol,VEC_TEMP,VEC_TEMP_MATOP,VEC_VV(0),ksp->vec_rhs);
240: KSPGMRESCycle(&its,ksp);
241: /* Store the Hessenberg matrix and the basis vectors of the Krylov subspace
242: if the cycle is complete for the computation of the Ritz pairs */
243: if (its == gmres->max_k) {
244: gmres->fullcycle++;
245: if (ksp->calc_ritz) {
246: if (!gmres->hes_ritz) {
247: PetscMalloc1(N*N,&gmres->hes_ritz);
248: PetscLogObjectMemory((PetscObject)ksp,N*N*sizeof(PetscScalar));
249: VecDuplicateVecs(VEC_VV(0),N,&gmres->vecb);
250: }
251: PetscArraycpy(gmres->hes_ritz,gmres->hes_origin,N*N);
252: for (i=0; i<gmres->max_k+1; i++) {
253: VecCopy(VEC_VV(i),gmres->vecb[i]);
254: }
255: }
256: }
257: itcount += its;
258: if (itcount >= ksp->max_it) {
259: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
260: break;
261: }
262: ksp->guess_zero = PETSC_FALSE; /* every future call to KSPInitialResidual() will have nonzero guess */
263: }
264: ksp->guess_zero = guess_zero; /* restore if user provided nonzero initial guess */
265: return 0;
266: }
268: PetscErrorCode KSPReset_GMRES(KSP ksp)
269: {
270: KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
271: PetscInt i;
273: /* Free the Hessenberg matrices */
274: PetscFree5(gmres->hh_origin,gmres->hes_origin,gmres->rs_origin,gmres->cc_origin,gmres->ss_origin);
275: PetscFree(gmres->hes_ritz);
277: /* free work vectors */
278: PetscFree(gmres->vecs);
279: for (i=0; i<gmres->nwork_alloc; i++) {
280: VecDestroyVecs(gmres->mwork_alloc[i],&gmres->user_work[i]);
281: }
282: gmres->nwork_alloc = 0;
283: if (gmres->vecb) {
284: VecDestroyVecs(gmres->max_k+1,&gmres->vecb);
285: }
287: PetscFree(gmres->user_work);
288: PetscFree(gmres->mwork_alloc);
289: PetscFree(gmres->nrs);
290: VecDestroy(&gmres->sol_temp);
291: PetscFree(gmres->Rsvd);
292: PetscFree(gmres->Dsvd);
293: PetscFree(gmres->orthogwork);
295: gmres->vv_allocated = 0;
296: gmres->vecs_allocated = 0;
297: gmres->sol_temp = NULL;
298: return 0;
299: }
301: PetscErrorCode KSPDestroy_GMRES(KSP ksp)
302: {
303: KSPReset_GMRES(ksp);
304: PetscFree(ksp->data);
305: /* clear composed functions */
306: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetPreAllocateVectors_C",NULL);
307: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetOrthogonalization_C",NULL);
308: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetOrthogonalization_C",NULL);
309: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetRestart_C",NULL);
310: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetRestart_C",NULL);
311: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetHapTol_C",NULL);
312: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetBreakdownTolerance_C",NULL);
313: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetCGSRefinementType_C",NULL);
314: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetCGSRefinementType_C",NULL);
315: return 0;
316: }
317: /*
318: KSPGMRESBuildSoln - create the solution from the starting vector and the
319: current iterates.
321: Input parameters:
322: nrs - work area of size it + 1.
323: vs - index of initial guess
324: vdest - index of result. Note that vs may == vdest (replace
325: guess with the solution).
327: This is an internal routine that knows about the GMRES internals.
328: */
329: static PetscErrorCode KSPGMRESBuildSoln(PetscScalar *nrs,Vec vs,Vec vdest,KSP ksp,PetscInt it)
330: {
331: PetscScalar tt;
332: PetscInt ii,k,j;
333: KSP_GMRES *gmres = (KSP_GMRES*)(ksp->data);
335: /* Solve for solution vector that minimizes the residual */
337: /* If it is < 0, no gmres steps have been performed */
338: if (it < 0) {
339: VecCopy(vs,vdest); /* VecCopy() is smart, exists immediately if vguess == vdest */
340: return 0;
341: }
342: if (*HH(it,it) != 0.0) {
343: nrs[it] = *GRS(it) / *HH(it,it);
344: } else {
346: else ksp->reason = KSP_DIVERGED_BREAKDOWN;
348: PetscInfo(ksp,"Likely your matrix or preconditioner is singular. HH(it,it) is identically zero; it = %D GRS(it) = %g\n",it,(double)PetscAbsScalar(*GRS(it)));
349: return 0;
350: }
351: for (ii=1; ii<=it; ii++) {
352: k = it - ii;
353: tt = *GRS(k);
354: for (j=k+1; j<=it; j++) tt = tt - *HH(k,j) * nrs[j];
355: if (*HH(k,k) == 0.0) {
357: else {
358: ksp->reason = KSP_DIVERGED_BREAKDOWN;
359: PetscInfo(ksp,"Likely your matrix or preconditioner is singular. HH(k,k) is identically zero; k = %D\n",k);
360: return 0;
361: }
362: }
363: nrs[k] = tt / *HH(k,k);
364: }
366: /* Accumulate the correction to the solution of the preconditioned problem in TEMP */
367: VecSet(VEC_TEMP,0.0);
368: VecMAXPY(VEC_TEMP,it+1,nrs,&VEC_VV(0));
370: KSPUnwindPreconditioner(ksp,VEC_TEMP,VEC_TEMP_MATOP);
371: /* add solution to previous solution */
372: if (vdest != vs) {
373: VecCopy(vs,vdest);
374: }
375: VecAXPY(vdest,1.0,VEC_TEMP);
376: return 0;
377: }
378: /*
379: Do the scalar work for the orthogonalization. Return new residual norm.
380: */
381: static PetscErrorCode KSPGMRESUpdateHessenberg(KSP ksp,PetscInt it,PetscBool hapend,PetscReal *res)
382: {
383: PetscScalar *hh,*cc,*ss,tt;
384: PetscInt j;
385: KSP_GMRES *gmres = (KSP_GMRES*)(ksp->data);
387: hh = HH(0,it);
388: cc = CC(0);
389: ss = SS(0);
391: /* Apply all the previously computed plane rotations to the new column
392: of the Hessenberg matrix */
393: for (j=1; j<=it; j++) {
394: tt = *hh;
395: *hh = PetscConj(*cc) * tt + *ss * *(hh+1);
396: hh++;
397: *hh = *cc++ * *hh - (*ss++ * tt);
398: }
400: /*
401: compute the new plane rotation, and apply it to:
402: 1) the right-hand-side of the Hessenberg system
403: 2) the new column of the Hessenberg matrix
404: thus obtaining the updated value of the residual
405: */
406: if (!hapend) {
407: tt = PetscSqrtScalar(PetscConj(*hh) * *hh + PetscConj(*(hh+1)) * *(hh+1));
408: if (tt == 0.0) {
410: else {
411: ksp->reason = KSP_DIVERGED_NULL;
412: return 0;
413: }
414: }
415: *cc = *hh / tt;
416: *ss = *(hh+1) / tt;
417: *GRS(it+1) = -(*ss * *GRS(it));
418: *GRS(it) = PetscConj(*cc) * *GRS(it);
419: *hh = PetscConj(*cc) * *hh + *ss * *(hh+1);
420: *res = PetscAbsScalar(*GRS(it+1));
421: } else {
422: /* happy breakdown: HH(it+1, it) = 0, therefore we don't need to apply
423: another rotation matrix (so RH doesn't change). The new residual is
424: always the new sine term times the residual from last time (GRS(it)),
425: but now the new sine rotation would be zero...so the residual should
426: be zero...so we will multiply "zero" by the last residual. This might
427: not be exactly what we want to do here -could just return "zero". */
429: *res = 0.0;
430: }
431: return 0;
432: }
433: /*
434: This routine allocates more work vectors, starting from VEC_VV(it).
435: */
436: PetscErrorCode KSPGMRESGetNewVectors(KSP ksp,PetscInt it)
437: {
438: KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
439: PetscInt nwork = gmres->nwork_alloc,k,nalloc;
441: nalloc = PetscMin(ksp->max_it,gmres->delta_allocate);
442: /* Adjust the number to allocate to make sure that we don't exceed the
443: number of available slots */
444: if (it + VEC_OFFSET + nalloc >= gmres->vecs_allocated) {
445: nalloc = gmres->vecs_allocated - it - VEC_OFFSET;
446: }
447: if (!nalloc) return 0;
449: gmres->vv_allocated += nalloc;
451: KSPCreateVecs(ksp,nalloc,&gmres->user_work[nwork],0,NULL);
452: PetscLogObjectParents(ksp,nalloc,gmres->user_work[nwork]);
454: gmres->mwork_alloc[nwork] = nalloc;
455: for (k=0; k<nalloc; k++) {
456: gmres->vecs[it+VEC_OFFSET+k] = gmres->user_work[nwork][k];
457: }
458: gmres->nwork_alloc++;
459: return 0;
460: }
462: PetscErrorCode KSPBuildSolution_GMRES(KSP ksp,Vec ptr,Vec *result)
463: {
464: KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
466: if (!ptr) {
467: if (!gmres->sol_temp) {
468: VecDuplicate(ksp->vec_sol,&gmres->sol_temp);
469: PetscLogObjectParent((PetscObject)ksp,(PetscObject)gmres->sol_temp);
470: }
471: ptr = gmres->sol_temp;
472: }
473: if (!gmres->nrs) {
474: /* allocate the work area */
475: PetscMalloc1(gmres->max_k,&gmres->nrs);
476: PetscLogObjectMemory((PetscObject)ksp,gmres->max_k);
477: }
479: KSPGMRESBuildSoln(gmres->nrs,ksp->vec_sol,ptr,ksp,gmres->it);
480: if (result) *result = ptr;
481: return 0;
482: }
484: PetscErrorCode KSPView_GMRES(KSP ksp,PetscViewer viewer)
485: {
486: KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
487: const char *cstr;
488: PetscBool iascii,isstring;
490: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
491: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
492: if (gmres->orthog == KSPGMRESClassicalGramSchmidtOrthogonalization) {
493: switch (gmres->cgstype) {
494: case (KSP_GMRES_CGS_REFINE_NEVER):
495: cstr = "Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement";
496: break;
497: case (KSP_GMRES_CGS_REFINE_ALWAYS):
498: cstr = "Classical (unmodified) Gram-Schmidt Orthogonalization with one step of iterative refinement";
499: break;
500: case (KSP_GMRES_CGS_REFINE_IFNEEDED):
501: cstr = "Classical (unmodified) Gram-Schmidt Orthogonalization with one step of iterative refinement when needed";
502: break;
503: default:
504: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Unknown orthogonalization");
505: }
506: } else if (gmres->orthog == KSPGMRESModifiedGramSchmidtOrthogonalization) {
507: cstr = "Modified Gram-Schmidt Orthogonalization";
508: } else {
509: cstr = "unknown orthogonalization";
510: }
511: if (iascii) {
512: PetscViewerASCIIPrintf(viewer," restart=%D, using %s\n",gmres->max_k,cstr);
513: PetscViewerASCIIPrintf(viewer," happy breakdown tolerance %g\n",(double)gmres->haptol);
514: } else if (isstring) {
515: PetscViewerStringSPrintf(viewer,"%s restart %D",cstr,gmres->max_k);
516: }
517: return 0;
518: }
520: /*@C
521: KSPGMRESMonitorKrylov - Calls VecView() for each new direction in the GMRES accumulated Krylov space.
523: Collective on ksp
525: Input Parameters:
526: + ksp - the KSP context
527: . its - iteration number
528: . fgnorm - 2-norm of residual (or gradient)
529: - dummy - an collection of viewers created with KSPViewerCreate()
531: Options Database Keys:
532: . -ksp_gmres_krylov_monitor <bool> - Plot the Krylov directions
534: Notes:
535: A new PETSCVIEWERDRAW is created for each Krylov vector so they can all be simultaneously viewed
536: Level: intermediate
538: .seealso: KSPMonitorSet(), KSPMonitorResidual(), VecView(), KSPViewersCreate(), KSPViewersDestroy()
539: @*/
540: PetscErrorCode KSPGMRESMonitorKrylov(KSP ksp,PetscInt its,PetscReal fgnorm,void *dummy)
541: {
542: PetscViewers viewers = (PetscViewers)dummy;
543: KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
544: Vec x;
545: PetscViewer viewer;
546: PetscBool flg;
548: PetscViewersGetViewer(viewers,gmres->it+1,&viewer);
549: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&flg);
550: if (!flg) {
551: PetscViewerSetType(viewer,PETSCVIEWERDRAW);
552: PetscViewerDrawSetInfo(viewer,NULL,"Krylov GMRES Monitor",PETSC_DECIDE,PETSC_DECIDE,300,300);
553: }
554: x = VEC_VV(gmres->it+1);
555: VecView(x,viewer);
556: return 0;
557: }
559: PetscErrorCode KSPSetFromOptions_GMRES(PetscOptionItems *PetscOptionsObject,KSP ksp)
560: {
562: PetscInt restart;
563: PetscReal haptol,breakdowntol;
564: KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
565: PetscBool flg;
567: PetscOptionsHead(PetscOptionsObject,"KSP GMRES Options");
568: PetscOptionsInt("-ksp_gmres_restart","Number of Krylov search directions","KSPGMRESSetRestart",gmres->max_k,&restart,&flg);
569: if (flg) KSPGMRESSetRestart(ksp,restart);
570: PetscOptionsReal("-ksp_gmres_haptol","Tolerance for exact convergence (happy ending)","KSPGMRESSetHapTol",gmres->haptol,&haptol,&flg);
571: if (flg) KSPGMRESSetHapTol(ksp,haptol);
572: PetscOptionsReal("-ksp_gmres_breakdown_tolerance","Divergence breakdown tolerance during GMRES restart","KSPGMRESSetBreakdownTolerance",gmres->breakdowntol,&breakdowntol,&flg);
573: if (flg) KSPGMRESSetBreakdownTolerance(ksp,breakdowntol);
574: flg = PETSC_FALSE;
575: PetscOptionsBool("-ksp_gmres_preallocate","Preallocate Krylov vectors","KSPGMRESSetPreAllocateVectors",flg,&flg,NULL);
576: if (flg) KSPGMRESSetPreAllocateVectors(ksp);
577: PetscOptionsBoolGroupBegin("-ksp_gmres_classicalgramschmidt","Classical (unmodified) Gram-Schmidt (fast)","KSPGMRESSetOrthogonalization",&flg);
578: if (flg) KSPGMRESSetOrthogonalization(ksp,KSPGMRESClassicalGramSchmidtOrthogonalization);
579: PetscOptionsBoolGroupEnd("-ksp_gmres_modifiedgramschmidt","Modified Gram-Schmidt (slow,more stable)","KSPGMRESSetOrthogonalization",&flg);
580: if (flg) KSPGMRESSetOrthogonalization(ksp,KSPGMRESModifiedGramSchmidtOrthogonalization);
581: PetscOptionsEnum("-ksp_gmres_cgs_refinement_type","Type of iterative refinement for classical (unmodified) Gram-Schmidt","KSPGMRESSetCGSRefinementType",
582: KSPGMRESCGSRefinementTypes,(PetscEnum)gmres->cgstype,(PetscEnum*)&gmres->cgstype,&flg);
583: flg = PETSC_FALSE;
584: PetscOptionsBool("-ksp_gmres_krylov_monitor","Plot the Krylov directions","KSPMonitorSet",flg,&flg,NULL);
585: if (flg) {
586: PetscViewers viewers;
587: PetscViewersCreate(PetscObjectComm((PetscObject)ksp),&viewers);
588: KSPMonitorSet(ksp,KSPGMRESMonitorKrylov,viewers,(PetscErrorCode (*)(void**))PetscViewersDestroy);
589: }
590: PetscOptionsTail();
591: return 0;
592: }
594: PetscErrorCode KSPGMRESSetHapTol_GMRES(KSP ksp,PetscReal tol)
595: {
596: KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
599: gmres->haptol = tol;
600: return 0;
601: }
603: PetscErrorCode KSPGMRESSetBreakdownTolerance_GMRES(KSP ksp,PetscReal tol)
604: {
605: KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
607: if (tol == PETSC_DEFAULT) {
608: gmres->breakdowntol = 0.1;
609: return 0;
610: }
612: gmres->breakdowntol = tol;
613: return 0;
614: }
616: PetscErrorCode KSPGMRESGetRestart_GMRES(KSP ksp,PetscInt *max_k)
617: {
618: KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
620: *max_k = gmres->max_k;
621: return 0;
622: }
624: PetscErrorCode KSPGMRESSetRestart_GMRES(KSP ksp,PetscInt max_k)
625: {
626: KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
629: if (!ksp->setupstage) {
630: gmres->max_k = max_k;
631: } else if (gmres->max_k != max_k) {
632: gmres->max_k = max_k;
633: ksp->setupstage = KSP_SETUP_NEW;
634: /* free the data structures, then create them again */
635: KSPReset_GMRES(ksp);
636: }
637: return 0;
638: }
640: PetscErrorCode KSPGMRESSetOrthogonalization_GMRES(KSP ksp,FCN fcn)
641: {
642: ((KSP_GMRES*)ksp->data)->orthog = fcn;
643: return 0;
644: }
646: PetscErrorCode KSPGMRESGetOrthogonalization_GMRES(KSP ksp,FCN *fcn)
647: {
648: *fcn = ((KSP_GMRES*)ksp->data)->orthog;
649: return 0;
650: }
652: PetscErrorCode KSPGMRESSetPreAllocateVectors_GMRES(KSP ksp)
653: {
654: KSP_GMRES *gmres;
656: gmres = (KSP_GMRES*)ksp->data;
657: gmres->q_preallocate = 1;
658: return 0;
659: }
661: PetscErrorCode KSPGMRESSetCGSRefinementType_GMRES(KSP ksp,KSPGMRESCGSRefinementType type)
662: {
663: KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
665: gmres->cgstype = type;
666: return 0;
667: }
669: PetscErrorCode KSPGMRESGetCGSRefinementType_GMRES(KSP ksp,KSPGMRESCGSRefinementType *type)
670: {
671: KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
673: *type = gmres->cgstype;
674: return 0;
675: }
677: /*@
678: KSPGMRESSetCGSRefinementType - Sets the type of iterative refinement to use
679: in the classical Gram Schmidt orthogonalization.
681: Logically Collective on ksp
683: Input Parameters:
684: + ksp - the Krylov space context
685: - type - the type of refinement
687: Options Database:
688: . -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always> - refinement type
690: Level: intermediate
692: .seealso: KSPGMRESSetOrthogonalization(), KSPGMRESCGSRefinementType, KSPGMRESClassicalGramSchmidtOrthogonalization(), KSPGMRESGetCGSRefinementType(),
693: KSPGMRESGetOrthogonalization()
694: @*/
695: PetscErrorCode KSPGMRESSetCGSRefinementType(KSP ksp,KSPGMRESCGSRefinementType type)
696: {
699: PetscTryMethod(ksp,"KSPGMRESSetCGSRefinementType_C",(KSP,KSPGMRESCGSRefinementType),(ksp,type));
700: return 0;
701: }
703: /*@
704: KSPGMRESGetCGSRefinementType - Gets the type of iterative refinement to use
705: in the classical Gram Schmidt orthogonalization.
707: Not Collective
709: Input Parameter:
710: . ksp - the Krylov space context
712: Output Parameter:
713: . type - the type of refinement
715: Options Database:
716: . -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always> - type of refinement
718: Level: intermediate
720: .seealso: KSPGMRESSetOrthogonalization(), KSPGMRESCGSRefinementType, KSPGMRESClassicalGramSchmidtOrthogonalization(), KSPGMRESSetCGSRefinementType(),
721: KSPGMRESGetOrthogonalization()
722: @*/
723: PetscErrorCode KSPGMRESGetCGSRefinementType(KSP ksp,KSPGMRESCGSRefinementType *type)
724: {
726: PetscUseMethod(ksp,"KSPGMRESGetCGSRefinementType_C",(KSP,KSPGMRESCGSRefinementType*),(ksp,type));
727: return 0;
728: }
730: /*@
731: KSPGMRESSetRestart - Sets number of iterations at which GMRES, FGMRES and LGMRES restarts.
733: Logically Collective on ksp
735: Input Parameters:
736: + ksp - the Krylov space context
737: - restart - integer restart value
739: Options Database:
740: . -ksp_gmres_restart <positive integer> - integer restart value
742: Note: The default value is 30.
744: Level: intermediate
746: .seealso: KSPSetTolerances(), KSPGMRESSetOrthogonalization(), KSPGMRESSetPreAllocateVectors(), KSPGMRESGetRestart()
747: @*/
748: PetscErrorCode KSPGMRESSetRestart(KSP ksp, PetscInt restart)
749: {
752: PetscTryMethod(ksp,"KSPGMRESSetRestart_C",(KSP,PetscInt),(ksp,restart));
753: return 0;
754: }
756: /*@
757: KSPGMRESGetRestart - Gets number of iterations at which GMRES, FGMRES and LGMRES restarts.
759: Not Collective
761: Input Parameter:
762: . ksp - the Krylov space context
764: Output Parameter:
765: . restart - integer restart value
767: Note: The default value is 30.
769: Level: intermediate
771: .seealso: KSPSetTolerances(), KSPGMRESSetOrthogonalization(), KSPGMRESSetPreAllocateVectors(), KSPGMRESSetRestart()
772: @*/
773: PetscErrorCode KSPGMRESGetRestart(KSP ksp, PetscInt *restart)
774: {
775: PetscUseMethod(ksp,"KSPGMRESGetRestart_C",(KSP,PetscInt*),(ksp,restart));
776: return 0;
777: }
779: /*@
780: KSPGMRESSetHapTol - Sets tolerance for determining happy breakdown in GMRES, FGMRES and LGMRES.
782: Logically Collective on ksp
784: Input Parameters:
785: + ksp - the Krylov space context
786: - tol - the tolerance
788: Options Database:
789: . -ksp_gmres_haptol <positive real value> - set tolerance for determining happy breakdown
791: Note: Happy breakdown is the rare case in GMRES where an 'exact' solution is obtained after
792: a certain number of iterations. If you attempt more iterations after this point unstable
793: things can happen hence very occasionally you may need to set this value to detect this condition
795: Level: intermediate
797: .seealso: KSPSetTolerances()
798: @*/
799: PetscErrorCode KSPGMRESSetHapTol(KSP ksp,PetscReal tol)
800: {
802: PetscTryMethod((ksp),"KSPGMRESSetHapTol_C",(KSP,PetscReal),((ksp),(tol)));
803: return 0;
804: }
806: /*@
807: KSPGMRESSetBreakdownTolerance - Sets tolerance for determining divergence breakdown in GMRES.
809: Logically Collective on ksp
811: Input Parameters:
812: + ksp - the Krylov space context
813: - tol - the tolerance
815: Options Database:
816: . -ksp_gmres_breakdown_tolerance <positive real value> - set tolerance for determining divergence breakdown
818: Note: divergence breakdown occurs when GMRES residual increases significantly
819: during restart
821: Level: intermediate
823: .seealso: KSPSetTolerances(), KSPGMRESSetHapTol()
824: @*/
825: PetscErrorCode KSPGMRESSetBreakdownTolerance(KSP ksp,PetscReal tol)
826: {
828: PetscTryMethod((ksp),"KSPGMRESSetBreakdownTolerance_C",(KSP,PetscReal),(ksp,tol));
829: return 0;
830: }
832: /*MC
833: KSPGMRES - Implements the Generalized Minimal Residual method.
834: (Saad and Schultz, 1986) with restart
836: Options Database Keys:
837: + -ksp_gmres_restart <restart> - the number of Krylov directions to orthogonalize against
838: . -ksp_gmres_haptol <tol> - sets the tolerance for "happy ending" (exact convergence)
839: . -ksp_gmres_preallocate - preallocate all the Krylov search directions initially (otherwise groups of
840: vectors are allocated as needed)
841: . -ksp_gmres_classicalgramschmidt - use classical (unmodified) Gram-Schmidt to orthogonalize against the Krylov space (fast) (the default)
842: . -ksp_gmres_modifiedgramschmidt - use modified Gram-Schmidt in the orthogonalization (more stable, but slower)
843: . -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always> - determine if iterative refinement is used to increase the
844: stability of the classical Gram-Schmidt orthogonalization.
845: - -ksp_gmres_krylov_monitor - plot the Krylov space generated
847: Level: beginner
849: Notes:
850: Left and right preconditioning are supported, but not symmetric preconditioning.
852: References:
853: . * - YOUCEF SAAD AND MARTIN H. SCHULTZ, GMRES: A GENERALIZED MINIMAL RESIDUAL ALGORITHM FOR SOLVING NONSYMMETRIC LINEAR SYSTEMS.
854: SIAM J. ScI. STAT. COMPUT. Vo|. 7, No. 3, July 1986.
856: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPFGMRES, KSPLGMRES,
857: KSPGMRESSetRestart(), KSPGMRESSetHapTol(), KSPGMRESSetPreAllocateVectors(), KSPGMRESSetOrthogonalization(), KSPGMRESGetOrthogonalization(),
858: KSPGMRESClassicalGramSchmidtOrthogonalization(), KSPGMRESModifiedGramSchmidtOrthogonalization(),
859: KSPGMRESCGSRefinementType, KSPGMRESSetCGSRefinementType(), KSPGMRESGetCGSRefinementType(), KSPGMRESMonitorKrylov(), KSPSetPCSide()
861: M*/
863: PETSC_EXTERN PetscErrorCode KSPCreate_GMRES(KSP ksp)
864: {
865: KSP_GMRES *gmres;
867: PetscNewLog(ksp,&gmres);
868: ksp->data = (void*)gmres;
870: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,4);
871: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,3);
872: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_SYMMETRIC,2);
873: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);
874: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);
876: ksp->ops->buildsolution = KSPBuildSolution_GMRES;
877: ksp->ops->setup = KSPSetUp_GMRES;
878: ksp->ops->solve = KSPSolve_GMRES;
879: ksp->ops->reset = KSPReset_GMRES;
880: ksp->ops->destroy = KSPDestroy_GMRES;
881: ksp->ops->view = KSPView_GMRES;
882: ksp->ops->setfromoptions = KSPSetFromOptions_GMRES;
883: ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_GMRES;
884: ksp->ops->computeeigenvalues = KSPComputeEigenvalues_GMRES;
885: #if !defined(PETSC_USE_COMPLEX)
886: ksp->ops->computeritz = KSPComputeRitz_GMRES;
887: #endif
888: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetPreAllocateVectors_C",KSPGMRESSetPreAllocateVectors_GMRES);
889: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetOrthogonalization_C",KSPGMRESSetOrthogonalization_GMRES);
890: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetOrthogonalization_C",KSPGMRESGetOrthogonalization_GMRES);
891: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetRestart_C",KSPGMRESSetRestart_GMRES);
892: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetRestart_C",KSPGMRESGetRestart_GMRES);
893: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetHapTol_C",KSPGMRESSetHapTol_GMRES);
894: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetBreakdownTolerance_C",KSPGMRESSetBreakdownTolerance_GMRES);
895: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetCGSRefinementType_C",KSPGMRESSetCGSRefinementType_GMRES);
896: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetCGSRefinementType_C",KSPGMRESGetCGSRefinementType_GMRES);
898: gmres->haptol = 1.0e-30;
899: gmres->breakdowntol = 0.1;
900: gmres->q_preallocate = 0;
901: gmres->delta_allocate = GMRES_DELTA_DIRECTIONS;
902: gmres->orthog = KSPGMRESClassicalGramSchmidtOrthogonalization;
903: gmres->nrs = NULL;
904: gmres->sol_temp = NULL;
905: gmres->max_k = GMRES_DEFAULT_MAXK;
906: gmres->Rsvd = NULL;
907: gmres->cgstype = KSP_GMRES_CGS_REFINE_NEVER;
908: gmres->orthogwork = NULL;
909: return 0;
910: }