Actual source code: gmres.c

petsc-3.6.1 2015-07-22
Report Typos and Errors
  2: /*
  3:     This file implements GMRES (a Generalized Minimal Residual) method.
  4:     Reference:  Saad and Schultz, 1986.


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

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

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

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

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

 47:   max_k = gmres->max_k;          /* restart size */
 48:   hh    = (max_k + 2) * (max_k + 1);
 49:   hes   = (max_k + 1) * (max_k + 1);
 50:   rs    = (max_k + 2);
 51:   cc    = (max_k + 1);

 53:   PetscCalloc5(hh,&gmres->hh_origin,hes,&gmres->hes_origin,rs,&gmres->rs_origin,cc,&gmres->cc_origin,cc,&gmres->ss_origin);
 54:   PetscLogObjectMemory((PetscObject)ksp,(hh + hes + rs + 2*cc)*sizeof(PetscScalar));

 56:   if (ksp->calc_sings) {
 57:     /* Allocate workspace to hold Hessenberg matrix needed by lapack */
 58:     PetscMalloc1((max_k + 3)*(max_k + 9),&gmres->Rsvd);
 59:     PetscLogObjectMemory((PetscObject)ksp,(max_k + 3)*(max_k + 9)*sizeof(PetscScalar));
 60:     PetscMalloc1(6*(max_k+2),&gmres->Dsvd);
 61:     PetscLogObjectMemory((PetscObject)ksp,6*(max_k+2)*sizeof(PetscReal));
 62:   }

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

 68:   PetscMalloc1(gmres->vecs_allocated,&gmres->vecs);
 69:   PetscMalloc1(VEC_OFFSET+2+max_k,&gmres->user_work);
 70:   PetscMalloc1(VEC_OFFSET+2+max_k,&gmres->mwork_alloc);
 71:   PetscLogObjectMemory((PetscObject)ksp,(VEC_OFFSET+2+max_k)*(sizeof(Vec*)+sizeof(PetscInt)) + gmres->vecs_allocated*sizeof(Vec));

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

 76:     KSPCreateVecs(ksp,gmres->vv_allocated,&gmres->user_work[0],0,NULL);
 77:     PetscLogObjectParents(ksp,gmres->vv_allocated,gmres->user_work[0]);

 79:     gmres->mwork_alloc[0] = gmres->vv_allocated;
 80:     gmres->nwork_alloc    = 1;
 81:     for (k=0; k<gmres->vv_allocated; k++) {
 82:       gmres->vecs[k] = gmres->user_work[0][k];
 83:     }
 84:   } else {
 85:     gmres->vv_allocated = 5;

 87:     KSPCreateVecs(ksp,5,&gmres->user_work[0],0,NULL);
 88:     PetscLogObjectParents(ksp,5,gmres->user_work[0]);

 90:     gmres->mwork_alloc[0] = 5;
 91:     gmres->nwork_alloc    = 1;
 92:     for (k=0; k<gmres->vv_allocated; k++) {
 93:       gmres->vecs[k] = gmres->user_work[0][k];
 94:     }
 95:   }
 96:   return(0);
 97: }

 99: /*
100:     Run gmres, possibly with restart.  Return residual history if requested.
101:     input parameters:

103: .        gmres  - structure containing parameters and work areas

105:     output parameters:
106: .        nres    - residuals (from preconditioned system) at each step.
107:                   If restarting, consider passing nres+it.  If null,
108:                   ignored
109: .        itcount - number of iterations used.  nres[0] to nres[itcount]
110:                   are defined.  If null, ignored.

112:     Notes:
113:     On entry, the value in vector VEC_VV(0) should be the initial residual
114:     (this allows shortcuts where the initial preconditioned residual is 0).
115:  */
118: PetscErrorCode KSPGMRESCycle(PetscInt *itcount,KSP ksp)
119: {
120:   KSP_GMRES      *gmres = (KSP_GMRES*)(ksp->data);
121:   PetscReal      res_norm,res,hapbnd,tt;
123:   PetscInt       it     = 0, max_k = gmres->max_k;
124:   PetscBool      hapend = PETSC_FALSE;

127:   VecNormalize(VEC_VV(0),&res_norm);
128:   res     = res_norm;
129:   *GRS(0) = res_norm;

131:   /* check for the convergence */
132:   PetscObjectSAWsTakeAccess((PetscObject)ksp);
133:   ksp->rnorm = res;
134:   PetscObjectSAWsGrantAccess((PetscObject)ksp);
135:   gmres->it  = (it - 1);
136:   KSPLogResidualHistory(ksp,res);
137:   KSPMonitor(ksp,ksp->its,res);
138:   if (!res) {
139:     if (itcount) *itcount = 0;
140:     ksp->reason = KSP_CONVERGED_ATOL;
141:     PetscInfo(ksp,"Converged due to zero residual norm on entry\n");
142:     return(0);
143:   }

145:   (*ksp->converged)(ksp,ksp->its,res,&ksp->reason,ksp->cnvP);
146:   while (!ksp->reason && it < max_k && ksp->its < ksp->max_it) {
147:     if (it) {
148:       KSPLogResidualHistory(ksp,res);
149:       KSPMonitor(ksp,ksp->its,res);
150:     }
151:     gmres->it = (it - 1);
152:     if (gmres->vv_allocated <= it + VEC_OFFSET + 1) {
153:       KSPGMRESGetNewVectors(ksp,it+1);
154:     }
155:     KSP_PCApplyBAorAB(ksp,VEC_VV(it),VEC_VV(1+it),VEC_TEMP_MATOP);

157:     /* update hessenberg matrix and do Gram-Schmidt */
158:     (*gmres->orthog)(ksp,it);
159:     if (ksp->reason) break;

161:     /* vv(i+1) . vv(i+1) */
162:     VecNormalize(VEC_VV(it+1),&tt);

164:     /* save the magnitude */
165:     *HH(it+1,it)  = tt;
166:     *HES(it+1,it) = tt;

168:     /* check for the happy breakdown */
169:     hapbnd = PetscAbsScalar(tt / *GRS(it));
170:     if (hapbnd > gmres->haptol) hapbnd = gmres->haptol;
171:     if (tt < hapbnd) {
172:       PetscInfo2(ksp,"Detected happy breakdown, current hapbnd = %14.12e tt = %14.12e\n",(double)hapbnd,(double)tt);
173:       hapend = PETSC_TRUE;
174:     }
175:     KSPGMRESUpdateHessenberg(ksp,it,hapend,&res);

177:     it++;
178:     gmres->it = (it-1);   /* For converged */
179:     ksp->its++;
180:     ksp->rnorm = res;
181:     if (ksp->reason) break;

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

185:     /* Catch error in happy breakdown and signal convergence and break from loop */
186:     if (hapend) {
187:       if (!ksp->reason) {
188:         if (ksp->errorifnotconverged) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"You reached the happy break down, but convergence was not indicated. Residual norm = %g",(double)res);
189:         else {
190:           ksp->reason = KSP_DIVERGED_BREAKDOWN;
191:           break;
192:         }
193:       }
194:     }
195:   }

197:   /* Monitor if we know that we will not return for a restart */
198:   if (it && (ksp->reason || ksp->its >= ksp->max_it)) {
199:     KSPLogResidualHistory(ksp,res);
200:     KSPMonitor(ksp,ksp->its,res);
201:   }

203:   if (itcount) *itcount = it;


206:   /*
207:     Down here we have to solve for the "best" coefficients of the Krylov
208:     columns, add the solution values together, and possibly unwind the
209:     preconditioning from the solution
210:    */
211:   /* Form the solution (or the solution so far) */
212:   KSPGMRESBuildSoln(GRS(0),ksp->vec_sol,ksp->vec_sol,ksp,it-1);
213:   return(0);
214: }

218: PetscErrorCode KSPSolve_GMRES(KSP ksp)
219: {
221:   PetscInt       its,itcount;
222:   KSP_GMRES      *gmres     = (KSP_GMRES*)ksp->data;
223:   PetscBool      guess_zero = ksp->guess_zero;

226:   if (ksp->calc_sings && !gmres->Rsvd) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ORDER,"Must call KSPSetComputeSingularValues() before KSPSetUp() is called");

228:   PetscObjectSAWsTakeAccess((PetscObject)ksp);
229:   ksp->its = 0;
230:   PetscObjectSAWsGrantAccess((PetscObject)ksp);

232:   itcount     = 0;
233:   ksp->reason = KSP_CONVERGED_ITERATING;
234:   while (!ksp->reason) {
235:     KSPInitialResidual(ksp,ksp->vec_sol,VEC_TEMP,VEC_TEMP_MATOP,VEC_VV(0),ksp->vec_rhs);
236:     KSPGMRESCycle(&its,ksp);
237:     itcount += its;
238:     if (itcount >= ksp->max_it) {
239:       if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
240:       break;
241:     }
242:     ksp->guess_zero = PETSC_FALSE; /* every future call to KSPInitialResidual() will have nonzero guess */
243:   }
244:   ksp->guess_zero = guess_zero; /* restore if user provided nonzero initial guess */
245:   return(0);
246: }

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

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

260:   /* free work vectors */
261:   PetscFree(gmres->vecs);
262:   for (i=0; i<gmres->nwork_alloc; i++) {
263:     VecDestroyVecs(gmres->mwork_alloc[i],&gmres->user_work[i]);
264:   }
265:   gmres->nwork_alloc = 0;

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

275:   gmres->sol_temp       = 0;
276:   gmres->vv_allocated   = 0;
277:   gmres->vecs_allocated = 0;
278:   gmres->sol_temp       = 0;
279:   return(0);
280: }

284: PetscErrorCode KSPDestroy_GMRES(KSP ksp)
285: {

289:   KSPReset_GMRES(ksp);
290:   PetscFree(ksp->data);
291:   /* clear composed functions */
292:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetPreAllocateVectors_C",NULL);
293:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetOrthogonalization_C",NULL);
294:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetOrthogonalization_C",NULL);
295:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetRestart_C",NULL);
296:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetRestart_C",NULL);
297:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetHapTol_C",NULL);
298:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetCGSRefinementType_C",NULL);
299:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetCGSRefinementType_C",NULL);
300:   return(0);
301: }
302: /*
303:     KSPGMRESBuildSoln - create the solution from the starting vector and the
304:     current iterates.

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

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

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

326:   /* If it is < 0, no gmres steps have been performed */
327:   if (it < 0) {
328:     VecCopy(vs,vdest); /* VecCopy() is smart, exists immediately if vguess == vdest */
329:     return(0);
330:   }
331:   if (*HH(it,it) != 0.0) {
332:     nrs[it] = *GRS(it) / *HH(it,it);
333:   } else {
334:     ksp->reason = KSP_DIVERGED_BREAKDOWN;

336:     PetscInfo2(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)));
337:     return(0);
338:   }
339:   for (ii=1; ii<=it; ii++) {
340:     k  = it - ii;
341:     tt = *GRS(k);
342:     for (j=k+1; j<=it; j++) tt = tt - *HH(k,j) * nrs[j];
343:     if (*HH(k,k) == 0.0) {
344:       ksp->reason = KSP_DIVERGED_BREAKDOWN;

346:       PetscInfo1(ksp,"Likely your matrix or preconditioner is singular. HH(k,k) is identically zero; k = %D\n",k);
347:       return(0);
348:     }
349:     nrs[k] = tt / *HH(k,k);
350:   }

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

356:   KSPUnwindPreconditioner(ksp,VEC_TEMP,VEC_TEMP_MATOP);
357:   /* add solution to previous solution */
358:   if (vdest != vs) {
359:     VecCopy(vs,vdest);
360:   }
361:   VecAXPY(vdest,1.0,VEC_TEMP);
362:   return(0);
363: }
364: /*
365:    Do the scalar work for the orthogonalization.  Return new residual norm.
366:  */
369: static PetscErrorCode KSPGMRESUpdateHessenberg(KSP ksp,PetscInt it,PetscBool hapend,PetscReal *res)
370: {
371:   PetscScalar *hh,*cc,*ss,tt;
372:   PetscInt    j;
373:   KSP_GMRES   *gmres = (KSP_GMRES*)(ksp->data);

376:   hh = HH(0,it);
377:   cc = CC(0);
378:   ss = SS(0);

380:   /* Apply all the previously computed plane rotations to the new column
381:      of the Hessenberg matrix */
382:   for (j=1; j<=it; j++) {
383:     tt  = *hh;
384:     *hh = PetscConj(*cc) * tt + *ss * *(hh+1);
385:     hh++;
386:     *hh = *cc++ * *hh - (*ss++ * tt);
387:   }

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

415:     *res = 0.0;
416:   }
417:   return(0);
418: }
419: /*
420:    This routine allocates more work vectors, starting from VEC_VV(it).
421:  */
424: PetscErrorCode KSPGMRESGetNewVectors(KSP ksp,PetscInt it)
425: {
426:   KSP_GMRES      *gmres = (KSP_GMRES*)ksp->data;
428:   PetscInt       nwork = gmres->nwork_alloc,k,nalloc;

431:   nalloc = PetscMin(ksp->max_it,gmres->delta_allocate);
432:   /* Adjust the number to allocate to make sure that we don't exceed the
433:     number of available slots */
434:   if (it + VEC_OFFSET + nalloc >= gmres->vecs_allocated) {
435:     nalloc = gmres->vecs_allocated - it - VEC_OFFSET;
436:   }
437:   if (!nalloc) return(0);

439:   gmres->vv_allocated += nalloc;

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

444:   gmres->mwork_alloc[nwork] = nalloc;
445:   for (k=0; k<nalloc; k++) {
446:     gmres->vecs[it+VEC_OFFSET+k] = gmres->user_work[nwork][k];
447:   }
448:   gmres->nwork_alloc++;
449:   return(0);
450: }

454: PetscErrorCode KSPBuildSolution_GMRES(KSP ksp,Vec ptr,Vec *result)
455: {
456:   KSP_GMRES      *gmres = (KSP_GMRES*)ksp->data;

460:   if (!ptr) {
461:     if (!gmres->sol_temp) {
462:       VecDuplicate(ksp->vec_sol,&gmres->sol_temp);
463:       PetscLogObjectParent((PetscObject)ksp,(PetscObject)gmres->sol_temp);
464:     }
465:     ptr = gmres->sol_temp;
466:   }
467:   if (!gmres->nrs) {
468:     /* allocate the work area */
469:     PetscMalloc1(gmres->max_k,&gmres->nrs);
470:     PetscLogObjectMemory((PetscObject)ksp,gmres->max_k*sizeof(PetscScalar));
471:   }

473:   KSPGMRESBuildSoln(gmres->nrs,ksp->vec_sol,ptr,ksp,gmres->it);
474:   if (result) *result = ptr;
475:   return(0);
476: }

480: PetscErrorCode KSPView_GMRES(KSP ksp,PetscViewer viewer)
481: {
482:   KSP_GMRES      *gmres = (KSP_GMRES*)ksp->data;
483:   const char     *cstr;
485:   PetscBool      iascii,isstring;

488:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
489:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
490:   if (gmres->orthog == KSPGMRESClassicalGramSchmidtOrthogonalization) {
491:     switch (gmres->cgstype) {
492:     case (KSP_GMRES_CGS_REFINE_NEVER):
493:       cstr = "Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement";
494:       break;
495:     case (KSP_GMRES_CGS_REFINE_ALWAYS):
496:       cstr = "Classical (unmodified) Gram-Schmidt Orthogonalization with one step of iterative refinement";
497:       break;
498:     case (KSP_GMRES_CGS_REFINE_IFNEEDED):
499:       cstr = "Classical (unmodified) Gram-Schmidt Orthogonalization with one step of iterative refinement when needed";
500:       break;
501:     default:
502:       SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Unknown orthogonalization");
503:     }
504:   } else if (gmres->orthog == KSPGMRESModifiedGramSchmidtOrthogonalization) {
505:     cstr = "Modified Gram-Schmidt Orthogonalization";
506:   } else {
507:     cstr = "unknown orthogonalization";
508:   }
509:   if (iascii) {
510:     PetscViewerASCIIPrintf(viewer,"  GMRES: restart=%D, using %s\n",gmres->max_k,cstr);
511:     PetscViewerASCIIPrintf(viewer,"  GMRES: happy breakdown tolerance %g\n",(double)gmres->haptol);
512:   } else if (isstring) {
513:     PetscViewerStringSPrintf(viewer,"%s restart %D",cstr,gmres->max_k);
514:   }
515:   return(0);
516: }

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_kyrlov_monitor

534:    Notes: A new PETSCVIEWERDRAW is created for each Krylov vector so they can all be simultaneously viewed
535:    Level: intermediate

537: .keywords: KSP, nonlinear, vector, monitor, view, Krylov space

539: .seealso: KSPMonitorSet(), KSPMonitorDefault(), VecView(), KSPViewersCreate(), KSPViewersDestroy()
540: @*/
541: PetscErrorCode  KSPGMRESMonitorKrylov(KSP ksp,PetscInt its,PetscReal fgnorm,void *dummy)
542: {
543:   PetscViewers   viewers = (PetscViewers)dummy;
544:   KSP_GMRES      *gmres  = (KSP_GMRES*)ksp->data;
546:   Vec            x;
547:   PetscViewer    viewer;
548:   PetscBool      flg;

551:   PetscViewersGetViewer(viewers,gmres->it+1,&viewer);
552:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&flg);
553:   if (!flg) {
554:     PetscViewerSetType(viewer,PETSCVIEWERDRAW);
555:     PetscViewerDrawSetInfo(viewer,NULL,"Krylov GMRES Monitor",PETSC_DECIDE,PETSC_DECIDE,300,300);
556:   }
557:   x    = VEC_VV(gmres->it+1);
558:   VecView(x,viewer);
559:   return(0);
560: }

564: PetscErrorCode KSPSetFromOptions_GMRES(PetscOptions *PetscOptionsObject,KSP ksp)
565: {
567:   PetscInt       restart;
568:   PetscReal      haptol;
569:   KSP_GMRES      *gmres = (KSP_GMRES*)ksp->data;
570:   PetscBool      flg;

573:   PetscOptionsHead(PetscOptionsObject,"KSP GMRES Options");
574:   PetscOptionsInt("-ksp_gmres_restart","Number of Krylov search directions","KSPGMRESSetRestart",gmres->max_k,&restart,&flg);
575:   if (flg) { KSPGMRESSetRestart(ksp,restart); }
576:   PetscOptionsReal("-ksp_gmres_haptol","Tolerance for exact convergence (happy ending)","KSPGMRESSetHapTol",gmres->haptol,&haptol,&flg);
577:   if (flg) { KSPGMRESSetHapTol(ksp,haptol); }
578:   flg  = PETSC_FALSE;
579:   PetscOptionsBool("-ksp_gmres_preallocate","Preallocate Krylov vectors","KSPGMRESSetPreAllocateVectors",flg,&flg,NULL);
580:   if (flg) {KSPGMRESSetPreAllocateVectors(ksp);}
581:   PetscOptionsBoolGroupBegin("-ksp_gmres_classicalgramschmidt","Classical (unmodified) Gram-Schmidt (fast)","KSPGMRESSetOrthogonalization",&flg);
582:   if (flg) {KSPGMRESSetOrthogonalization(ksp,KSPGMRESClassicalGramSchmidtOrthogonalization);}
583:   PetscOptionsBoolGroupEnd("-ksp_gmres_modifiedgramschmidt","Modified Gram-Schmidt (slow,more stable)","KSPGMRESSetOrthogonalization",&flg);
584:   if (flg) {KSPGMRESSetOrthogonalization(ksp,KSPGMRESModifiedGramSchmidtOrthogonalization);}
585:   PetscOptionsEnum("-ksp_gmres_cgs_refinement_type","Type of iterative refinement for classical (unmodified) Gram-Schmidt","KSPGMRESSetCGSRefinementType",
586:                           KSPGMRESCGSRefinementTypes,(PetscEnum)gmres->cgstype,(PetscEnum*)&gmres->cgstype,&flg);
587:   flg  = PETSC_FALSE;
588:   PetscOptionsBool("-ksp_gmres_krylov_monitor","Plot the Krylov directions","KSPMonitorSet",flg,&flg,NULL);
589:   if (flg) {
590:     PetscViewers viewers;
591:     PetscViewersCreate(PetscObjectComm((PetscObject)ksp),&viewers);
592:     KSPMonitorSet(ksp,KSPGMRESMonitorKrylov,viewers,(PetscErrorCode (*)(void**))PetscViewersDestroy);
593:   }
594:   PetscOptionsTail();
595:   return(0);
596: }

598: extern PetscErrorCode KSPComputeExtremeSingularValues_GMRES(KSP,PetscReal*,PetscReal*);
599: extern PetscErrorCode KSPComputeEigenvalues_GMRES(KSP,PetscInt,PetscReal*,PetscReal*,PetscInt*);

603: PetscErrorCode  KSPGMRESSetHapTol_GMRES(KSP ksp,PetscReal tol)
604: {
605:   KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;

608:   if (tol < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Tolerance must be non-negative");
609:   gmres->haptol = tol;
610:   return(0);
611: }

615: PetscErrorCode  KSPGMRESGetRestart_GMRES(KSP ksp,PetscInt *max_k)
616: {
617:   KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;

620:   *max_k = gmres->max_k;
621:   return(0);
622: }

626: PetscErrorCode  KSPGMRESSetRestart_GMRES(KSP ksp,PetscInt max_k)
627: {
628:   KSP_GMRES      *gmres = (KSP_GMRES*)ksp->data;

632:   if (max_k < 1) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Restart must be positive");
633:   if (!ksp->setupstage) {
634:     gmres->max_k = max_k;
635:   } else if (gmres->max_k != max_k) {
636:     gmres->max_k    = max_k;
637:     ksp->setupstage = KSP_SETUP_NEW;
638:     /* free the data structures, then create them again */
639:     KSPReset_GMRES(ksp);
640:   }
641:   return(0);
642: }

646: PetscErrorCode  KSPGMRESSetOrthogonalization_GMRES(KSP ksp,FCN fcn)
647: {
649:   ((KSP_GMRES*)ksp->data)->orthog = fcn;
650:   return(0);
651: }

655: PetscErrorCode  KSPGMRESGetOrthogonalization_GMRES(KSP ksp,FCN *fcn)
656: {
658:   *fcn = ((KSP_GMRES*)ksp->data)->orthog;
659:   return(0);
660: }

664: PetscErrorCode  KSPGMRESSetPreAllocateVectors_GMRES(KSP ksp)
665: {
666:   KSP_GMRES *gmres;

669:   gmres = (KSP_GMRES*)ksp->data;
670:   gmres->q_preallocate = 1;
671:   return(0);
672: }

676: PetscErrorCode  KSPGMRESSetCGSRefinementType_GMRES(KSP ksp,KSPGMRESCGSRefinementType type)
677: {
678:   KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;

681:   gmres->cgstype = type;
682:   return(0);
683: }

687: PetscErrorCode  KSPGMRESGetCGSRefinementType_GMRES(KSP ksp,KSPGMRESCGSRefinementType *type)
688: {
689:   KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;

692:   *type = gmres->cgstype;
693:   return(0);
694: }

698: /*@
699:    KSPGMRESSetCGSRefinementType - Sets the type of iterative refinement to use
700:          in the classical Gram Schmidt orthogonalization.

702:    Logically Collective on KSP

704:    Input Parameters:
705: +  ksp - the Krylov space context
706: -  type - the type of refinement

708:   Options Database:
709: .  -ksp_gmres_cgs_refinement_type <never,ifneeded,always>

711:    Level: intermediate

713: .keywords: KSP, GMRES, iterative refinement

715: .seealso: KSPGMRESSetOrthogonalization(), KSPGMRESCGSRefinementType, KSPGMRESClassicalGramSchmidtOrthogonalization(), KSPGMRESGetCGSRefinementType(),
716:           KSPGMRESGetOrthogonalization()
717: @*/
718: PetscErrorCode  KSPGMRESSetCGSRefinementType(KSP ksp,KSPGMRESCGSRefinementType type)
719: {

725:   PetscTryMethod(ksp,"KSPGMRESSetCGSRefinementType_C",(KSP,KSPGMRESCGSRefinementType),(ksp,type));
726:   return(0);
727: }

731: /*@
732:    KSPGMRESGetCGSRefinementType - Gets the type of iterative refinement to use
733:          in the classical Gram Schmidt orthogonalization.

735:    Not Collective

737:    Input Parameter:
738: .  ksp - the Krylov space context

740:    Output Parameter:
741: .  type - the type of refinement

743:   Options Database:
744: .  -ksp_gmres_cgs_refinement_type <never,ifneeded,always>

746:    Level: intermediate

748: .keywords: KSP, GMRES, iterative refinement

750: .seealso: KSPGMRESSetOrthogonalization(), KSPGMRESCGSRefinementType, KSPGMRESClassicalGramSchmidtOrthogonalization(), KSPGMRESSetCGSRefinementType(),
751:           KSPGMRESGetOrthogonalization()
752: @*/
753: PetscErrorCode  KSPGMRESGetCGSRefinementType(KSP ksp,KSPGMRESCGSRefinementType *type)
754: {

759:   PetscUseMethod(ksp,"KSPGMRESGetCGSRefinementType_C",(KSP,KSPGMRESCGSRefinementType*),(ksp,type));
760:   return(0);
761: }


766: /*@
767:    KSPGMRESSetRestart - Sets number of iterations at which GMRES, FGMRES and LGMRES restarts.

769:    Logically Collective on KSP

771:    Input Parameters:
772: +  ksp - the Krylov space context
773: -  restart - integer restart value

775:   Options Database:
776: .  -ksp_gmres_restart <positive integer>

778:     Note: The default value is 30.

780:    Level: intermediate

782: .keywords: KSP, GMRES, restart, iterations

784: .seealso: KSPSetTolerances(), KSPGMRESSetOrthogonalization(), KSPGMRESSetPreAllocateVectors(), KSPGMRESGetRestart()
785: @*/
786: PetscErrorCode  KSPGMRESSetRestart(KSP ksp, PetscInt restart)
787: {


793:   PetscTryMethod(ksp,"KSPGMRESSetRestart_C",(KSP,PetscInt),(ksp,restart));
794:   return(0);
795: }

799: /*@
800:    KSPGMRESGetRestart - Gets number of iterations at which GMRES, FGMRES and LGMRES restarts.

802:    Not Collective

804:    Input Parameter:
805: .  ksp - the Krylov space context

807:    Output Parameter:
808: .   restart - integer restart value

810:     Note: The default value is 30.

812:    Level: intermediate

814: .keywords: KSP, GMRES, restart, iterations

816: .seealso: KSPSetTolerances(), KSPGMRESSetOrthogonalization(), KSPGMRESSetPreAllocateVectors(), KSPGMRESSetRestart()
817: @*/
818: PetscErrorCode  KSPGMRESGetRestart(KSP ksp, PetscInt *restart)
819: {

823:   PetscTryMethod(ksp,"KSPGMRESGetRestart_C",(KSP,PetscInt*),(ksp,restart));
824:   return(0);
825: }

829: /*@
830:    KSPGMRESSetHapTol - Sets tolerance for determining happy breakdown in GMRES, FGMRES and LGMRES.

832:    Logically Collective on KSP

834:    Input Parameters:
835: +  ksp - the Krylov space context
836: -  tol - the tolerance

838:   Options Database:
839: .  -ksp_gmres_haptol <positive real value>

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

845:    Level: intermediate

847: .keywords: KSP, GMRES, tolerance

849: .seealso: KSPSetTolerances()
850: @*/
851: PetscErrorCode  KSPGMRESSetHapTol(KSP ksp,PetscReal tol)
852: {

857:   PetscTryMethod((ksp),"KSPGMRESSetHapTol_C",(KSP,PetscReal),((ksp),(tol)));
858:   return(0);
859: }

861: /*MC
862:      KSPGMRES - Implements the Generalized Minimal Residual method.
863:                 (Saad and Schultz, 1986) with restart


866:    Options Database Keys:
867: +   -ksp_gmres_restart <restart> - the number of Krylov directions to orthogonalize against
868: .   -ksp_gmres_haptol <tol> - sets the tolerance for "happy ending" (exact convergence)
869: .   -ksp_gmres_preallocate - preallocate all the Krylov search directions initially (otherwise groups of
870:                              vectors are allocated as needed)
871: .   -ksp_gmres_classicalgramschmidt - use classical (unmodified) Gram-Schmidt to orthogonalize against the Krylov space (fast) (the default)
872: .   -ksp_gmres_modifiedgramschmidt - use modified Gram-Schmidt in the orthogonalization (more stable, but slower)
873: .   -ksp_gmres_cgs_refinement_type <never,ifneeded,always> - determine if iterative refinement is used to increase the
874:                                    stability of the classical Gram-Schmidt  orthogonalization.
875: -   -ksp_gmres_krylov_monitor - plot the Krylov space generated

877:    Level: beginner

879:    Notes: Left and right preconditioning are supported, but not symmetric preconditioning.

881:    References:
882:      GMRES: A GENERALIZED MINIMAL RESIDUAL ALGORITHM FOR SOLVING NONSYMMETRIC LINEAR SYSTEMS. YOUCEF SAAD AND MARTIN H. SCHULTZ,
883:           SIAM J. ScI. STAT. COMPUT. Vo|. 7, No. 3, July 1986, pp. 856--869.

885: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPFGMRES, KSPLGMRES,
886:            KSPGMRESSetRestart(), KSPGMRESSetHapTol(), KSPGMRESSetPreAllocateVectors(), KSPGMRESSetOrthogonalization(), KSPGMRESGetOrthogonalization(),
887:            KSPGMRESClassicalGramSchmidtOrthogonalization(), KSPGMRESModifiedGramSchmidtOrthogonalization(),
888:            KSPGMRESCGSRefinementType, KSPGMRESSetCGSRefinementType(), KSPGMRESGetCGSRefinementType(), KSPGMRESMonitorKrylov(), KSPSetPCSide()

890: M*/

894: PETSC_EXTERN PetscErrorCode KSPCreate_GMRES(KSP ksp)
895: {
896:   KSP_GMRES      *gmres;

900:   PetscNewLog(ksp,&gmres);
901:   ksp->data = (void*)gmres;

903:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
904:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);

906:   ksp->ops->buildsolution                = KSPBuildSolution_GMRES;
907:   ksp->ops->setup                        = KSPSetUp_GMRES;
908:   ksp->ops->solve                        = KSPSolve_GMRES;
909:   ksp->ops->reset                        = KSPReset_GMRES;
910:   ksp->ops->destroy                      = KSPDestroy_GMRES;
911:   ksp->ops->view                         = KSPView_GMRES;
912:   ksp->ops->setfromoptions               = KSPSetFromOptions_GMRES;
913:   ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_GMRES;
914:   ksp->ops->computeeigenvalues           = KSPComputeEigenvalues_GMRES;

916:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetPreAllocateVectors_C",KSPGMRESSetPreAllocateVectors_GMRES);
917:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetOrthogonalization_C",KSPGMRESSetOrthogonalization_GMRES);
918:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetOrthogonalization_C",KSPGMRESGetOrthogonalization_GMRES);
919:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetRestart_C",KSPGMRESSetRestart_GMRES);
920:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetRestart_C",KSPGMRESGetRestart_GMRES);
921:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetHapTol_C",KSPGMRESSetHapTol_GMRES);
922:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetCGSRefinementType_C",KSPGMRESSetCGSRefinementType_GMRES);
923:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetCGSRefinementType_C",KSPGMRESGetCGSRefinementType_GMRES);

925:   gmres->haptol         = 1.0e-30;
926:   gmres->q_preallocate  = 0;
927:   gmres->delta_allocate = GMRES_DELTA_DIRECTIONS;
928:   gmres->orthog         = KSPGMRESClassicalGramSchmidtOrthogonalization;
929:   gmres->nrs            = 0;
930:   gmres->sol_temp       = 0;
931:   gmres->max_k          = GMRES_DEFAULT_MAXK;
932:   gmres->Rsvd           = 0;
933:   gmres->cgstype        = KSP_GMRES_CGS_REFINE_NEVER;
934:   gmres->orthogwork     = 0;
935:   return(0);
936: }