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