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;
 40:   PetscInt       max_k,k;
 41:   KSP_GMRES      *gmres = (KSP_GMRES*)ksp->data;

 44:   max_k = gmres->max_k;          /* restart size */
 45:   hh    = (max_k + 2) * (max_k + 1);
 46:   hes   = (max_k + 1) * (max_k + 1);
 47:   rs    = (max_k + 2);
 48:   cc    = (max_k + 1);

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

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

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

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

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

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

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

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

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

 96: /*
 97:     Run gmres, possibly with restart.  Return residual history if requested.
 98:     input parameters:

100: .        gmres  - structure containing parameters and work areas

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

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

122:   if (itcount) *itcount = 0;
123:   VecNormalize(VEC_VV(0),&res);
124:   KSPCheckNorm(ksp,res);

126:   /* the constant .1 is arbitrary, just some measure at how incorrect the residuals are */
127:   if ((ksp->rnorm > 0.0) && (PetscAbsReal(res-ksp->rnorm) > gmres->breakdowntol*gmres->rnorm0)) {
128:     if (ksp->errorifnotconverged) SETERRQ3(PetscObjectComm((PetscObject)ksp),PETSC_ERR_CONV_FAILED,"Residual norm computed by GMRES recursion formula %g is far from the computed residual norm %g at restart, residual norm at start of cycle %g",(double)ksp->rnorm,(double)res,(double)gmres->rnorm0);
129:     else {
130:       PetscInfo3(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);
131:       ksp->reason =  KSP_DIVERGED_BREAKDOWN;
132:       return(0);
133:     }
134:   }
135:   *GRS(0) = gmres->rnorm0 = res;

137:   /* check for the convergence */
138:   PetscObjectSAWsTakeAccess((PetscObject)ksp);
139:   ksp->rnorm = res;
140:   PetscObjectSAWsGrantAccess((PetscObject)ksp);
141:   gmres->it  = (it - 1);
142:   KSPLogResidualHistory(ksp,res);
143:   KSPLogErrorHistory(ksp);
144:   KSPMonitor(ksp,ksp->its,res);
145:   if (!res) {
146:     ksp->reason = KSP_CONVERGED_ATOL;
147:     PetscInfo(ksp,"Converged due to zero residual norm on entry\n");
148:     return(0);
149:   }

151:   (*ksp->converged)(ksp,ksp->its,res,&ksp->reason,ksp->cnvP);
152:   while (!ksp->reason && it < max_k && ksp->its < ksp->max_it) {
153:     if (it) {
154:       KSPLogResidualHistory(ksp,res);
155:       KSPLogErrorHistory(ksp);
156:       KSPMonitor(ksp,ksp->its,res);
157:     }
158:     gmres->it = (it - 1);
159:     if (gmres->vv_allocated <= it + VEC_OFFSET + 1) {
160:       KSPGMRESGetNewVectors(ksp,it+1);
161:     }
162:     KSP_PCApplyBAorAB(ksp,VEC_VV(it),VEC_VV(1+it),VEC_TEMP_MATOP);

164:     /* update hessenberg matrix and do Gram-Schmidt */
165:     (*gmres->orthog)(ksp,it);
166:     if (ksp->reason) break;

168:     /* vv(i+1) . vv(i+1) */
169:     VecNormalize(VEC_VV(it+1),&tt);
170:     KSPCheckNorm(ksp,tt);

172:     /* save the magnitude */
173:     *HH(it+1,it)  = tt;
174:     *HES(it+1,it) = tt;

176:     /* check for the happy breakdown */
177:     hapbnd = PetscAbsScalar(tt / *GRS(it));
178:     if (hapbnd > gmres->haptol) hapbnd = gmres->haptol;
179:     if (tt < hapbnd) {
180:       PetscInfo2(ksp,"Detected happy breakdown, current hapbnd = %14.12e tt = %14.12e\n",(double)hapbnd,(double)tt);
181:       hapend = PETSC_TRUE;
182:     }
183:     KSPGMRESUpdateHessenberg(ksp,it,hapend,&res);

185:     it++;
186:     gmres->it = (it-1);   /* For converged */
187:     ksp->its++;
188:     ksp->rnorm = res;
189:     if (ksp->reason) break;

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

193:     /* Catch error in happy breakdown and signal convergence and break from loop */
194:     if (hapend) {
195:       if (ksp->normtype == KSP_NORM_NONE) { /* convergence test was skipped in this case */
196:         ksp->reason = KSP_CONVERGED_HAPPY_BREAKDOWN;
197:       } else if (!ksp->reason) {
198:         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);
199:         else {
200:           ksp->reason = KSP_DIVERGED_BREAKDOWN;
201:           break;
202:         }
203:       }
204:     }
205:   }

207:   /* Monitor if we know that we will not return for a restart */
208:   if (it && (ksp->reason || ksp->its >= ksp->max_it)) {
209:     KSPLogResidualHistory(ksp,res);
210:     KSPLogErrorHistory(ksp);
211:     KSPMonitor(ksp,ksp->its,res);
212:   }

214:   if (itcount) *itcount = it;

216:   /*
217:     Down here we have to solve for the "best" coefficients of the Krylov
218:     columns, add the solution values together, and possibly unwind the
219:     preconditioning from the solution
220:    */
221:   /* Form the solution (or the solution so far) */
222:   KSPGMRESBuildSoln(GRS(0),ksp->vec_sol,ksp->vec_sol,ksp,it-1);
223:   return(0);
224: }

226: PetscErrorCode KSPSolve_GMRES(KSP ksp)
227: {
229:   PetscInt       its,itcount,i;
230:   KSP_GMRES      *gmres     = (KSP_GMRES*)ksp->data;
231:   PetscBool      guess_zero = ksp->guess_zero;
232:   PetscInt       N = gmres->max_k + 1;

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

237:   PetscObjectSAWsTakeAccess((PetscObject)ksp);
238:   ksp->its = 0;
239:   PetscObjectSAWsGrantAccess((PetscObject)ksp);

241:   itcount     = 0;
242:   gmres->fullcycle = 0;
243:   ksp->reason = KSP_CONVERGED_ITERATING;
244:   ksp->rnorm  = -1.0; /* special marker for KSPGMRESCycle() */
245:   while (!ksp->reason) {
246:     KSPInitialResidual(ksp,ksp->vec_sol,VEC_TEMP,VEC_TEMP_MATOP,VEC_VV(0),ksp->vec_rhs);
247:     KSPGMRESCycle(&its,ksp);
248:     /* Store the Hessenberg matrix and the basis vectors of the Krylov subspace
249:     if the cycle is complete for the computation of the Ritz pairs */
250:     if (its == gmres->max_k) {
251:       gmres->fullcycle++;
252:       if (ksp->calc_ritz) {
253:         if (!gmres->hes_ritz) {
254:           PetscMalloc1(N*N,&gmres->hes_ritz);
255:           PetscLogObjectMemory((PetscObject)ksp,N*N*sizeof(PetscScalar));
256:           VecDuplicateVecs(VEC_VV(0),N,&gmres->vecb);
257:         }
258:         PetscArraycpy(gmres->hes_ritz,gmres->hes_origin,N*N);
259:         for (i=0; i<gmres->max_k+1; i++) {
260:           VecCopy(VEC_VV(i),gmres->vecb[i]);
261:         }
262:       }
263:     }
264:     itcount += its;
265:     if (itcount >= ksp->max_it) {
266:       if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
267:       break;
268:     }
269:     ksp->guess_zero = PETSC_FALSE; /* every future call to KSPInitialResidual() will have nonzero guess */
270:   }
271:   ksp->guess_zero = guess_zero; /* restore if user provided nonzero initial guess */
272:   return(0);
273: }

275: PetscErrorCode KSPReset_GMRES(KSP ksp)
276: {
277:   KSP_GMRES      *gmres = (KSP_GMRES*)ksp->data;
279:   PetscInt       i;

282:   /* Free the Hessenberg matrices */
283:   PetscFree5(gmres->hh_origin,gmres->hes_origin,gmres->rs_origin,gmres->cc_origin,gmres->ss_origin);
284:   PetscFree(gmres->hes_ritz);

286:   /* free work vectors */
287:   PetscFree(gmres->vecs);
288:   for (i=0; i<gmres->nwork_alloc; i++) {
289:     VecDestroyVecs(gmres->mwork_alloc[i],&gmres->user_work[i]);
290:   }
291:   gmres->nwork_alloc = 0;
292:   if (gmres->vecb)  {
293:     VecDestroyVecs(gmres->max_k+1,&gmres->vecb);
294:   }

296:   PetscFree(gmres->user_work);
297:   PetscFree(gmres->mwork_alloc);
298:   PetscFree(gmres->nrs);
299:   VecDestroy(&gmres->sol_temp);
300:   PetscFree(gmres->Rsvd);
301:   PetscFree(gmres->Dsvd);
302:   PetscFree(gmres->orthogwork);

304:   gmres->vv_allocated   = 0;
305:   gmres->vecs_allocated = 0;
306:   gmres->sol_temp       = NULL;
307:   return(0);
308: }

310: PetscErrorCode KSPDestroy_GMRES(KSP ksp)
311: {

315:   KSPReset_GMRES(ksp);
316:   PetscFree(ksp->data);
317:   /* clear composed functions */
318:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetPreAllocateVectors_C",NULL);
319:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetOrthogonalization_C",NULL);
320:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetOrthogonalization_C",NULL);
321:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetRestart_C",NULL);
322:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetRestart_C",NULL);
323:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetHapTol_C",NULL);
324:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetBreakdownTolerance_C",NULL);
325:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetCGSRefinementType_C",NULL);
326:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetCGSRefinementType_C",NULL);
327:   return(0);
328: }
329: /*
330:     KSPGMRESBuildSoln - create the solution from the starting vector and the
331:     current iterates.

333:     Input parameters:
334:         nrs - work area of size it + 1.
335:         vs  - index of initial guess
336:         vdest - index of result.  Note that vs may == vdest (replace
337:                 guess with the solution).

339:      This is an internal routine that knows about the GMRES internals.
340:  */
341: static PetscErrorCode KSPGMRESBuildSoln(PetscScalar *nrs,Vec vs,Vec vdest,KSP ksp,PetscInt it)
342: {
343:   PetscScalar    tt;
345:   PetscInt       ii,k,j;
346:   KSP_GMRES      *gmres = (KSP_GMRES*)(ksp->data);

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

351:   /* If it is < 0, no gmres steps have been performed */
352:   if (it < 0) {
353:     VecCopy(vs,vdest); /* VecCopy() is smart, exists immediately if vguess == vdest */
354:     return(0);
355:   }
356:   if (*HH(it,it) != 0.0) {
357:     nrs[it] = *GRS(it) / *HH(it,it);
358:   } else {
359:     if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"You reached the break down in GMRES; HH(it,it) = 0");
360:     else ksp->reason = KSP_DIVERGED_BREAKDOWN;

362:     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)));
363:     return(0);
364:   }
365:   for (ii=1; ii<=it; ii++) {
366:     k  = it - ii;
367:     tt = *GRS(k);
368:     for (j=k+1; j<=it; j++) tt = tt - *HH(k,j) * nrs[j];
369:     if (*HH(k,k) == 0.0) {
370:       if (ksp->errorifnotconverged) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"Likely your matrix or preconditioner is singular. HH(k,k) is identically zero; k = %D\n",k);
371:       else {
372:         ksp->reason = KSP_DIVERGED_BREAKDOWN;
373:         PetscInfo1(ksp,"Likely your matrix or preconditioner is singular. HH(k,k) is identically zero; k = %D\n",k);
374:         return(0);
375:       }
376:     }
377:     nrs[k] = tt / *HH(k,k);
378:   }

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

384:   KSPUnwindPreconditioner(ksp,VEC_TEMP,VEC_TEMP_MATOP);
385:   /* add solution to previous solution */
386:   if (vdest != vs) {
387:     VecCopy(vs,vdest);
388:   }
389:   VecAXPY(vdest,1.0,VEC_TEMP);
390:   return(0);
391: }
392: /*
393:    Do the scalar work for the orthogonalization.  Return new residual norm.
394:  */
395: static PetscErrorCode KSPGMRESUpdateHessenberg(KSP ksp,PetscInt it,PetscBool hapend,PetscReal *res)
396: {
397:   PetscScalar *hh,*cc,*ss,tt;
398:   PetscInt    j;
399:   KSP_GMRES   *gmres = (KSP_GMRES*)(ksp->data);

402:   hh = HH(0,it);
403:   cc = CC(0);
404:   ss = SS(0);

406:   /* Apply all the previously computed plane rotations to the new column
407:      of the Hessenberg matrix */
408:   for (j=1; j<=it; j++) {
409:     tt  = *hh;
410:     *hh = PetscConj(*cc) * tt + *ss * *(hh+1);
411:     hh++;
412:     *hh = *cc++ * *hh - (*ss++ * tt);
413:   }

415:   /*
416:     compute the new plane rotation, and apply it to:
417:      1) the right-hand-side of the Hessenberg system
418:      2) the new column of the Hessenberg matrix
419:     thus obtaining the updated value of the residual
420:   */
421:   if (!hapend) {
422:     tt = PetscSqrtScalar(PetscConj(*hh) * *hh + PetscConj(*(hh+1)) * *(hh+1));
423:     if (tt == 0.0) {
424:       if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"tt == 0.0");
425:       else {
426:         ksp->reason = KSP_DIVERGED_NULL;
427:         return(0);
428:       }
429:     }
430:     *cc        = *hh / tt;
431:     *ss        = *(hh+1) / tt;
432:     *GRS(it+1) = -(*ss * *GRS(it));
433:     *GRS(it)   = PetscConj(*cc) * *GRS(it);
434:     *hh        = PetscConj(*cc) * *hh + *ss * *(hh+1);
435:     *res       = PetscAbsScalar(*GRS(it+1));
436:   } else {
437:     /* happy breakdown: HH(it+1, it) = 0, therefore we don't need to apply
438:             another rotation matrix (so RH doesn't change).  The new residual is
439:             always the new sine term times the residual from last time (GRS(it)),
440:             but now the new sine rotation would be zero...so the residual should
441:             be zero...so we will multiply "zero" by the last residual.  This might
442:             not be exactly what we want to do here -could just return "zero". */

444:     *res = 0.0;
445:   }
446:   return(0);
447: }
448: /*
449:    This routine allocates more work vectors, starting from VEC_VV(it).
450:  */
451: PetscErrorCode KSPGMRESGetNewVectors(KSP ksp,PetscInt it)
452: {
453:   KSP_GMRES      *gmres = (KSP_GMRES*)ksp->data;
455:   PetscInt       nwork = gmres->nwork_alloc,k,nalloc;

458:   nalloc = PetscMin(ksp->max_it,gmres->delta_allocate);
459:   /* Adjust the number to allocate to make sure that we don't exceed the
460:     number of available slots */
461:   if (it + VEC_OFFSET + nalloc >= gmres->vecs_allocated) {
462:     nalloc = gmres->vecs_allocated - it - VEC_OFFSET;
463:   }
464:   if (!nalloc) return(0);

466:   gmres->vv_allocated += nalloc;

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

471:   gmres->mwork_alloc[nwork] = nalloc;
472:   for (k=0; k<nalloc; k++) {
473:     gmres->vecs[it+VEC_OFFSET+k] = gmres->user_work[nwork][k];
474:   }
475:   gmres->nwork_alloc++;
476:   return(0);
477: }

479: PetscErrorCode KSPBuildSolution_GMRES(KSP ksp,Vec ptr,Vec *result)
480: {
481:   KSP_GMRES      *gmres = (KSP_GMRES*)ksp->data;

485:   if (!ptr) {
486:     if (!gmres->sol_temp) {
487:       VecDuplicate(ksp->vec_sol,&gmres->sol_temp);
488:       PetscLogObjectParent((PetscObject)ksp,(PetscObject)gmres->sol_temp);
489:     }
490:     ptr = gmres->sol_temp;
491:   }
492:   if (!gmres->nrs) {
493:     /* allocate the work area */
494:     PetscMalloc1(gmres->max_k,&gmres->nrs);
495:     PetscLogObjectMemory((PetscObject)ksp,gmres->max_k);
496:   }

498:   KSPGMRESBuildSoln(gmres->nrs,ksp->vec_sol,ptr,ksp,gmres->it);
499:   if (result) *result = ptr;
500:   return(0);
501: }

503: PetscErrorCode KSPView_GMRES(KSP ksp,PetscViewer viewer)
504: {
505:   KSP_GMRES      *gmres = (KSP_GMRES*)ksp->data;
506:   const char     *cstr;
508:   PetscBool      iascii,isstring;

511:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
512:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
513:   if (gmres->orthog == KSPGMRESClassicalGramSchmidtOrthogonalization) {
514:     switch (gmres->cgstype) {
515:     case (KSP_GMRES_CGS_REFINE_NEVER):
516:       cstr = "Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement";
517:       break;
518:     case (KSP_GMRES_CGS_REFINE_ALWAYS):
519:       cstr = "Classical (unmodified) Gram-Schmidt Orthogonalization with one step of iterative refinement";
520:       break;
521:     case (KSP_GMRES_CGS_REFINE_IFNEEDED):
522:       cstr = "Classical (unmodified) Gram-Schmidt Orthogonalization with one step of iterative refinement when needed";
523:       break;
524:     default:
525:       SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Unknown orthogonalization");
526:     }
527:   } else if (gmres->orthog == KSPGMRESModifiedGramSchmidtOrthogonalization) {
528:     cstr = "Modified Gram-Schmidt Orthogonalization";
529:   } else {
530:     cstr = "unknown orthogonalization";
531:   }
532:   if (iascii) {
533:     PetscViewerASCIIPrintf(viewer,"  restart=%D, using %s\n",gmres->max_k,cstr);
534:     PetscViewerASCIIPrintf(viewer,"  happy breakdown tolerance %g\n",(double)gmres->haptol);
535:   } else if (isstring) {
536:     PetscViewerStringSPrintf(viewer,"%s restart %D",cstr,gmres->max_k);
537:   }
538:   return(0);
539: }

541: /*@C
542:    KSPGMRESMonitorKrylov - Calls VecView() for each new direction in the GMRES accumulated Krylov space.

544:    Collective on ksp

546:    Input Parameters:
547: +  ksp - the KSP context
548: .  its - iteration number
549: .  fgnorm - 2-norm of residual (or gradient)
550: -  dummy - an collection of viewers created with KSPViewerCreate()

552:    Options Database Keys:
553: .   -ksp_gmres_kyrlov_monitor

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

559: .seealso: KSPMonitorSet(), KSPMonitorResidual(), VecView(), KSPViewersCreate(), KSPViewersDestroy()
560: @*/
561: PetscErrorCode  KSPGMRESMonitorKrylov(KSP ksp,PetscInt its,PetscReal fgnorm,void *dummy)
562: {
563:   PetscViewers   viewers = (PetscViewers)dummy;
564:   KSP_GMRES      *gmres  = (KSP_GMRES*)ksp->data;
566:   Vec            x;
567:   PetscViewer    viewer;
568:   PetscBool      flg;

571:   PetscViewersGetViewer(viewers,gmres->it+1,&viewer);
572:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&flg);
573:   if (!flg) {
574:     PetscViewerSetType(viewer,PETSCVIEWERDRAW);
575:     PetscViewerDrawSetInfo(viewer,NULL,"Krylov GMRES Monitor",PETSC_DECIDE,PETSC_DECIDE,300,300);
576:   }
577:   x    = VEC_VV(gmres->it+1);
578:   VecView(x,viewer);
579:   return(0);
580: }

582: PetscErrorCode KSPSetFromOptions_GMRES(PetscOptionItems *PetscOptionsObject,KSP ksp)
583: {
585:   PetscInt       restart;
586:   PetscReal      haptol,breakdowntol;
587:   KSP_GMRES      *gmres = (KSP_GMRES*)ksp->data;
588:   PetscBool      flg;

591:   PetscOptionsHead(PetscOptionsObject,"KSP GMRES Options");
592:   PetscOptionsInt("-ksp_gmres_restart","Number of Krylov search directions","KSPGMRESSetRestart",gmres->max_k,&restart,&flg);
593:   if (flg) { KSPGMRESSetRestart(ksp,restart); }
594:   PetscOptionsReal("-ksp_gmres_haptol","Tolerance for exact convergence (happy ending)","KSPGMRESSetHapTol",gmres->haptol,&haptol,&flg);
595:   if (flg) { KSPGMRESSetHapTol(ksp,haptol); }
596:   PetscOptionsReal("-ksp_gmres_breakdown_tolerance","Divergence breakdown tolerance during GMRES restart","KSPGMRESSetBreakdownTolerance",gmres->breakdowntol,&breakdowntol,&flg);
597:   if (flg) { KSPGMRESSetBreakdownTolerance(ksp,breakdowntol); }
598:   flg  = PETSC_FALSE;
599:   PetscOptionsBool("-ksp_gmres_preallocate","Preallocate Krylov vectors","KSPGMRESSetPreAllocateVectors",flg,&flg,NULL);
600:   if (flg) {KSPGMRESSetPreAllocateVectors(ksp);}
601:   PetscOptionsBoolGroupBegin("-ksp_gmres_classicalgramschmidt","Classical (unmodified) Gram-Schmidt (fast)","KSPGMRESSetOrthogonalization",&flg);
602:   if (flg) {KSPGMRESSetOrthogonalization(ksp,KSPGMRESClassicalGramSchmidtOrthogonalization);}
603:   PetscOptionsBoolGroupEnd("-ksp_gmres_modifiedgramschmidt","Modified Gram-Schmidt (slow,more stable)","KSPGMRESSetOrthogonalization",&flg);
604:   if (flg) {KSPGMRESSetOrthogonalization(ksp,KSPGMRESModifiedGramSchmidtOrthogonalization);}
605:   PetscOptionsEnum("-ksp_gmres_cgs_refinement_type","Type of iterative refinement for classical (unmodified) Gram-Schmidt","KSPGMRESSetCGSRefinementType",
606:                           KSPGMRESCGSRefinementTypes,(PetscEnum)gmres->cgstype,(PetscEnum*)&gmres->cgstype,&flg);
607:   flg  = PETSC_FALSE;
608:   PetscOptionsBool("-ksp_gmres_krylov_monitor","Plot the Krylov directions","KSPMonitorSet",flg,&flg,NULL);
609:   if (flg) {
610:     PetscViewers viewers;
611:     PetscViewersCreate(PetscObjectComm((PetscObject)ksp),&viewers);
612:     KSPMonitorSet(ksp,KSPGMRESMonitorKrylov,viewers,(PetscErrorCode (*)(void**))PetscViewersDestroy);
613:   }
614:   PetscOptionsTail();
615:   return(0);
616: }

618: PetscErrorCode  KSPGMRESSetHapTol_GMRES(KSP ksp,PetscReal tol)
619: {
620:   KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;

623:   if (tol < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Tolerance must be non-negative");
624:   gmres->haptol = tol;
625:   return(0);
626: }

628: PetscErrorCode  KSPGMRESSetBreakdownTolerance_GMRES(KSP ksp,PetscReal tol)
629: {
630:   KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;

633:   if (tol == PETSC_DEFAULT) {
634:     gmres->breakdowntol = 0.1;
635:     return(0);
636:   }
637:   if (tol < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Breakdown tolerance must be non-negative");
638:   gmres->breakdowntol = tol;
639:   return(0);
640: }

642: PetscErrorCode  KSPGMRESGetRestart_GMRES(KSP ksp,PetscInt *max_k)
643: {
644:   KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;

647:   *max_k = gmres->max_k;
648:   return(0);
649: }

651: PetscErrorCode  KSPGMRESSetRestart_GMRES(KSP ksp,PetscInt max_k)
652: {
653:   KSP_GMRES      *gmres = (KSP_GMRES*)ksp->data;

657:   if (max_k < 1) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Restart must be positive");
658:   if (!ksp->setupstage) {
659:     gmres->max_k = max_k;
660:   } else if (gmres->max_k != max_k) {
661:     gmres->max_k    = max_k;
662:     ksp->setupstage = KSP_SETUP_NEW;
663:     /* free the data structures, then create them again */
664:     KSPReset_GMRES(ksp);
665:   }
666:   return(0);
667: }

669: PetscErrorCode  KSPGMRESSetOrthogonalization_GMRES(KSP ksp,FCN fcn)
670: {
672:   ((KSP_GMRES*)ksp->data)->orthog = fcn;
673:   return(0);
674: }

676: PetscErrorCode  KSPGMRESGetOrthogonalization_GMRES(KSP ksp,FCN *fcn)
677: {
679:   *fcn = ((KSP_GMRES*)ksp->data)->orthog;
680:   return(0);
681: }

683: PetscErrorCode  KSPGMRESSetPreAllocateVectors_GMRES(KSP ksp)
684: {
685:   KSP_GMRES *gmres;

688:   gmres = (KSP_GMRES*)ksp->data;
689:   gmres->q_preallocate = 1;
690:   return(0);
691: }

693: PetscErrorCode  KSPGMRESSetCGSRefinementType_GMRES(KSP ksp,KSPGMRESCGSRefinementType type)
694: {
695:   KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;

698:   gmres->cgstype = type;
699:   return(0);
700: }

702: PetscErrorCode  KSPGMRESGetCGSRefinementType_GMRES(KSP ksp,KSPGMRESCGSRefinementType *type)
703: {
704:   KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;

707:   *type = gmres->cgstype;
708:   return(0);
709: }

711: /*@
712:    KSPGMRESSetCGSRefinementType - Sets the type of iterative refinement to use
713:          in the classical Gram Schmidt orthogonalization.

715:    Logically Collective on ksp

717:    Input Parameters:
718: +  ksp - the Krylov space context
719: -  type - the type of refinement

721:   Options Database:
722: .  -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always>

724:    Level: intermediate

726: .seealso: KSPGMRESSetOrthogonalization(), KSPGMRESCGSRefinementType, KSPGMRESClassicalGramSchmidtOrthogonalization(), KSPGMRESGetCGSRefinementType(),
727:           KSPGMRESGetOrthogonalization()
728: @*/
729: PetscErrorCode  KSPGMRESSetCGSRefinementType(KSP ksp,KSPGMRESCGSRefinementType type)
730: {

736:   PetscTryMethod(ksp,"KSPGMRESSetCGSRefinementType_C",(KSP,KSPGMRESCGSRefinementType),(ksp,type));
737:   return(0);
738: }

740: /*@
741:    KSPGMRESGetCGSRefinementType - Gets the type of iterative refinement to use
742:          in the classical Gram Schmidt orthogonalization.

744:    Not Collective

746:    Input Parameter:
747: .  ksp - the Krylov space context

749:    Output Parameter:
750: .  type - the type of refinement

752:   Options Database:
753: .  -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always>

755:    Level: intermediate

757: .seealso: KSPGMRESSetOrthogonalization(), KSPGMRESCGSRefinementType, KSPGMRESClassicalGramSchmidtOrthogonalization(), KSPGMRESSetCGSRefinementType(),
758:           KSPGMRESGetOrthogonalization()
759: @*/
760: PetscErrorCode  KSPGMRESGetCGSRefinementType(KSP ksp,KSPGMRESCGSRefinementType *type)
761: {

766:   PetscUseMethod(ksp,"KSPGMRESGetCGSRefinementType_C",(KSP,KSPGMRESCGSRefinementType*),(ksp,type));
767:   return(0);
768: }

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

773:    Logically Collective on ksp

775:    Input Parameters:
776: +  ksp - the Krylov space context
777: -  restart - integer restart value

779:   Options Database:
780: .  -ksp_gmres_restart <positive integer>

782:     Note: The default value is 30.

784:    Level: intermediate

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


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

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: .seealso: KSPSetTolerances(), KSPGMRESSetOrthogonalization(), KSPGMRESSetPreAllocateVectors(), KSPGMRESSetRestart()
815: @*/
816: PetscErrorCode  KSPGMRESGetRestart(KSP ksp, PetscInt *restart)
817: {

821:   PetscUseMethod(ksp,"KSPGMRESGetRestart_C",(KSP,PetscInt*),(ksp,restart));
822:   return(0);
823: }

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

828:    Logically Collective on ksp

830:    Input Parameters:
831: +  ksp - the Krylov space context
832: -  tol - the tolerance

834:   Options Database:
835: .  -ksp_gmres_haptol <positive real value>

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

841:    Level: intermediate

843: .seealso: KSPSetTolerances()
844: @*/
845: PetscErrorCode  KSPGMRESSetHapTol(KSP ksp,PetscReal tol)
846: {

851:   PetscTryMethod((ksp),"KSPGMRESSetHapTol_C",(KSP,PetscReal),((ksp),(tol)));
852:   return(0);
853: }

855: /*@
856:    KSPGMRESSetBreakdownTolerance - Sets tolerance for determining divergence breakdown in GMRES.

858:    Logically Collective on ksp

860:    Input Parameters:
861: +  ksp - the Krylov space context
862: -  tol - the tolerance

864:   Options Database:
865: .  -ksp_gmres_breakdown_tolerance <positive real value>

867:    Note: divergence breakdown occurs when GMRES residual increases significantly
868:          during restart

870:    Level: intermediate

872: .seealso: KSPSetTolerances(), KSPGMRESSetHapTol()
873: @*/
874: PetscErrorCode  KSPGMRESSetBreakdownTolerance(KSP ksp,PetscReal tol)
875: {

880:   PetscTryMethod((ksp),"KSPGMRESSetBreakdownTolerance_C",(KSP,PetscReal),(ksp,tol));
881:   return(0);
882: }

884: /*MC
885:      KSPGMRES - Implements the Generalized Minimal Residual method.
886:                 (Saad and Schultz, 1986) with restart

888:    Options Database Keys:
889: +   -ksp_gmres_restart <restart> - the number of Krylov directions to orthogonalize against
890: .   -ksp_gmres_haptol <tol> - sets the tolerance for "happy ending" (exact convergence)
891: .   -ksp_gmres_preallocate - preallocate all the Krylov search directions initially (otherwise groups of
892:                              vectors are allocated as needed)
893: .   -ksp_gmres_classicalgramschmidt - use classical (unmodified) Gram-Schmidt to orthogonalize against the Krylov space (fast) (the default)
894: .   -ksp_gmres_modifiedgramschmidt - use modified Gram-Schmidt in the orthogonalization (more stable, but slower)
895: .   -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always> - determine if iterative refinement is used to increase the
896:                                    stability of the classical Gram-Schmidt  orthogonalization.
897: -   -ksp_gmres_krylov_monitor - plot the Krylov space generated

899:    Level: beginner

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

904:    References:
905: .     1. - YOUCEF SAAD AND MARTIN H. SCHULTZ, GMRES: A GENERALIZED MINIMAL RESIDUAL ALGORITHM FOR SOLVING NONSYMMETRIC LINEAR SYSTEMS.
906:           SIAM J. ScI. STAT. COMPUT. Vo|. 7, No. 3, July 1986.

908: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPFGMRES, KSPLGMRES,
909:            KSPGMRESSetRestart(), KSPGMRESSetHapTol(), KSPGMRESSetPreAllocateVectors(), KSPGMRESSetOrthogonalization(), KSPGMRESGetOrthogonalization(),
910:            KSPGMRESClassicalGramSchmidtOrthogonalization(), KSPGMRESModifiedGramSchmidtOrthogonalization(),
911:            KSPGMRESCGSRefinementType, KSPGMRESSetCGSRefinementType(), KSPGMRESGetCGSRefinementType(), KSPGMRESMonitorKrylov(), KSPSetPCSide()

913: M*/

915: PETSC_EXTERN PetscErrorCode KSPCreate_GMRES(KSP ksp)
916: {
917:   KSP_GMRES      *gmres;

921:   PetscNewLog(ksp,&gmres);
922:   ksp->data = (void*)gmres;

924:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,4);
925:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,3);
926:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_SYMMETRIC,2);
927:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);
928:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);

930:   ksp->ops->buildsolution                = KSPBuildSolution_GMRES;
931:   ksp->ops->setup                        = KSPSetUp_GMRES;
932:   ksp->ops->solve                        = KSPSolve_GMRES;
933:   ksp->ops->reset                        = KSPReset_GMRES;
934:   ksp->ops->destroy                      = KSPDestroy_GMRES;
935:   ksp->ops->view                         = KSPView_GMRES;
936:   ksp->ops->setfromoptions               = KSPSetFromOptions_GMRES;
937:   ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_GMRES;
938:   ksp->ops->computeeigenvalues           = KSPComputeEigenvalues_GMRES;
939: #if !defined(PETSC_USE_COMPLEX)
940:   ksp->ops->computeritz                  = KSPComputeRitz_GMRES;
941: #endif
942:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetPreAllocateVectors_C",KSPGMRESSetPreAllocateVectors_GMRES);
943:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetOrthogonalization_C",KSPGMRESSetOrthogonalization_GMRES);
944:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetOrthogonalization_C",KSPGMRESGetOrthogonalization_GMRES);
945:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetRestart_C",KSPGMRESSetRestart_GMRES);
946:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetRestart_C",KSPGMRESGetRestart_GMRES);
947:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetHapTol_C",KSPGMRESSetHapTol_GMRES);
948:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetBreakdownTolerance_C",KSPGMRESSetBreakdownTolerance_GMRES);
949:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetCGSRefinementType_C",KSPGMRESSetCGSRefinementType_GMRES);
950:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetCGSRefinementType_C",KSPGMRESGetCGSRefinementType_GMRES);

952:   gmres->haptol         = 1.0e-30;
953:   gmres->breakdowntol   = 0.1;
954:   gmres->q_preallocate  = 0;
955:   gmres->delta_allocate = GMRES_DELTA_DIRECTIONS;
956:   gmres->orthog         = KSPGMRESClassicalGramSchmidtOrthogonalization;
957:   gmres->nrs            = NULL;
958:   gmres->sol_temp       = NULL;
959:   gmres->max_k          = GMRES_DEFAULT_MAXK;
960:   gmres->Rsvd           = NULL;
961:   gmres->cgstype        = KSP_GMRES_CGS_REFINE_NEVER;
962:   gmres->orthogwork     = NULL;
963:   return(0);
964: }