Actual source code: lgmres.c


  2: #include <../src/ksp/ksp/impls/gmres/lgmres/lgmresimpl.h>

  4: #define LGMRES_DELTA_DIRECTIONS 10
  5: #define LGMRES_DEFAULT_MAXK     30
  6: #define LGMRES_DEFAULT_AUGDIM   2 /*default number of augmentation vectors */
  7: static PetscErrorCode    KSPLGMRESGetNewVectors(KSP,PetscInt);
  8: static PetscErrorCode    KSPLGMRESUpdateHessenberg(KSP,PetscInt,PetscBool,PetscReal*);
  9: static PetscErrorCode    KSPLGMRESBuildSoln(PetscScalar*,Vec,Vec,KSP,PetscInt);

 11: PetscErrorCode  KSPLGMRESSetAugDim(KSP ksp, PetscInt dim)
 12: {
 13:   PetscTryMethod((ksp),"KSPLGMRESSetAugDim_C",(KSP,PetscInt),(ksp,dim));
 14:   return 0;
 15: }

 17: PetscErrorCode  KSPLGMRESSetConstant(KSP ksp)
 18: {
 19:   PetscTryMethod((ksp),"KSPLGMRESSetConstant_C",(KSP),(ksp));
 20:   return 0;
 21: }

 23: /*
 24:     KSPSetUp_LGMRES - Sets up the workspace needed by lgmres.

 26:     This is called once, usually automatically by KSPSolve() or KSPSetUp(),
 27:     but can be called directly by KSPSetUp().

 29: */
 30: PetscErrorCode    KSPSetUp_LGMRES(KSP ksp)
 31: {
 32:   PetscInt       max_k,k, aug_dim;
 33:   KSP_LGMRES     *lgmres = (KSP_LGMRES*)ksp->data;

 35:   max_k   = lgmres->max_k;
 36:   aug_dim = lgmres->aug_dim;
 37:   KSPSetUp_GMRES(ksp);

 39:   /* need array of pointers to augvecs*/
 40:   PetscMalloc1(2*aug_dim + AUG_OFFSET,&lgmres->augvecs);

 42:   lgmres->aug_vecs_allocated = 2 *aug_dim + AUG_OFFSET;

 44:   PetscMalloc1(2*aug_dim + AUG_OFFSET,&lgmres->augvecs_user_work);
 45:   PetscMalloc1(aug_dim,&lgmres->aug_order);
 46:   PetscLogObjectMemory((PetscObject)ksp,(aug_dim)*(4*sizeof(void*) + sizeof(PetscInt)) + AUG_OFFSET*2*sizeof(void*));

 48:   /*  for now we will preallocate the augvecs - because aug_dim << restart
 49:      ... also keep in mind that we need to keep augvecs from cycle to cycle*/
 50:   lgmres->aug_vv_allocated = 2* aug_dim + AUG_OFFSET;
 51:   lgmres->augwork_alloc    =  2* aug_dim + AUG_OFFSET;

 53:   KSPCreateVecs(ksp,lgmres->aug_vv_allocated,&lgmres->augvecs_user_work[0],0,NULL);
 54:   PetscMalloc1(max_k+1,&lgmres->hwork);
 55:   PetscLogObjectParents(ksp,lgmres->aug_vv_allocated,lgmres->augvecs_user_work[0]);
 56:   for (k=0; k<lgmres->aug_vv_allocated; k++) {
 57:     lgmres->augvecs[k] = lgmres->augvecs_user_work[0][k];
 58:   }
 59:   return 0;
 60: }

 62: /*

 64:     KSPLGMRESCycle - Run lgmres, possibly with restart.  Return residual
 65:                   history if requested.

 67:     input parameters:
 68: .        lgmres  - structure containing parameters and work areas

 70:     output parameters:
 71: .        nres    - residuals (from preconditioned system) at each step.
 72:                   If restarting, consider passing nres+it.  If null,
 73:                   ignored
 74: .        itcount - number of iterations used.   nres[0] to nres[itcount]
 75:                   are defined.  If null, ignored.  If null, ignored.
 76: .        converged - 0 if not converged

 78:     Notes:
 79:     On entry, the value in vector VEC_VV(0) should be
 80:     the initial residual.

 82:  */
 83: PetscErrorCode KSPLGMRESCycle(PetscInt *itcount,KSP ksp)
 84: {
 85:   KSP_LGMRES     *lgmres = (KSP_LGMRES*)(ksp->data);
 86:   PetscReal      res_norm, res;
 87:   PetscReal      hapbnd, tt;
 88:   PetscScalar    tmp;
 89:   PetscBool      hapend = PETSC_FALSE;  /* indicates happy breakdown ending */
 90:   PetscInt       loc_it;                /* local count of # of dir. in Krylov space */
 91:   PetscInt       max_k  = lgmres->max_k; /* max approx space size */
 92:   PetscInt       max_it = ksp->max_it;  /* max # of overall iterations for the method */

 94:   /* LGMRES_MOD - new variables*/
 95:   PetscInt    aug_dim = lgmres->aug_dim;
 96:   PetscInt    spot    = 0;
 97:   PetscInt    order   = 0;
 98:   PetscInt    it_arnoldi;                /* number of arnoldi steps to take */
 99:   PetscInt    it_total;                  /* total number of its to take (=approx space size)*/
100:   PetscInt    ii, jj;
101:   PetscReal   tmp_norm;
102:   PetscScalar inv_tmp_norm;
103:   PetscScalar *avec;

105:   /* Number of pseudo iterations since last restart is the number
106:      of prestart directions */
107:   loc_it = 0;

109:   /* LGMRES_MOD: determine number of arnoldi steps to take */
110:   /* if approx_constant then we keep the space the same size even if
111:      we don't have the full number of aug vectors yet*/
112:   if (lgmres->approx_constant) it_arnoldi = max_k - lgmres->aug_ct;
113:   else it_arnoldi = max_k - aug_dim;

115:   it_total =  it_arnoldi + lgmres->aug_ct;

117:   /* initial residual is in VEC_VV(0)  - compute its norm*/
118:   VecNorm(VEC_VV(0),NORM_2,&res_norm);
119:   KSPCheckNorm(ksp,res_norm);
120:   res  = res_norm;

122:   /* first entry in right-hand-side of hessenberg system is just
123:      the initial residual norm */
124:   *GRS(0) = res_norm;

126:   /* check for the convergence */
127:   if (!res) {
128:     if (itcount) *itcount = 0;
129:     ksp->reason = KSP_CONVERGED_ATOL;
130:     PetscInfo(ksp,"Converged due to zero residual norm on entry\n");
131:     return 0;
132:   }

134:   /* scale VEC_VV (the initial residual) */
135:   tmp = 1.0/res_norm; VecScale(VEC_VV(0),tmp);

137:   if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = res;
138:   else ksp->rnorm = 0.0;

140:   /* note: (lgmres->it) is always set one less than (loc_it) It is used in
141:      KSPBUILDSolution_LGMRES, where it is passed to KSPLGMRESBuildSoln.
142:      Note that when KSPLGMRESBuildSoln is called from this function,
143:      (loc_it -1) is passed, so the two are equivalent */
144:   lgmres->it = (loc_it - 1);

146:   /* MAIN ITERATION LOOP BEGINNING*/

148:   /* keep iterating until we have converged OR generated the max number
149:      of directions OR reached the max number of iterations for the method */
150:   (*ksp->converged)(ksp,ksp->its,res,&ksp->reason,ksp->cnvP);

152:   while (!ksp->reason && loc_it < it_total && ksp->its < max_it) { /* LGMRES_MOD: changed to it_total */
153:     KSPLogResidualHistory(ksp,res);
154:     lgmres->it = (loc_it - 1);
155:     KSPMonitor(ksp,ksp->its,res);

157:     /* see if more space is needed for work vectors */
158:     if (lgmres->vv_allocated <= loc_it + VEC_OFFSET + 1) {
159:       KSPLGMRESGetNewVectors(ksp,loc_it+1);
160:       /* (loc_it+1) is passed in as number of the first vector that should
161:           be allocated */
162:     }

164:     /*LGMRES_MOD: decide whether this is an arnoldi step or an aug step */
165:     if (loc_it < it_arnoldi) { /* Arnoldi */
166:       KSP_PCApplyBAorAB(ksp,VEC_VV(loc_it),VEC_VV(1+loc_it),VEC_TEMP_MATOP);
167:     } else { /*aug step */
168:       order = loc_it - it_arnoldi + 1; /* which aug step */
169:       for (ii=0; ii<aug_dim; ii++) {
170:         if (lgmres->aug_order[ii] == order) {
171:           spot = ii;
172:           break; /* must have this because there will be duplicates before aug_ct = aug_dim */
173:         }
174:       }

176:       VecCopy(A_AUGVEC(spot), VEC_VV(1+loc_it));
177:       /*note: an alternate implementation choice would be to only save the AUGVECS and
178:         not A_AUGVEC and then apply the PC here to the augvec */
179:     }

181:     /* update hessenberg matrix and do Gram-Schmidt - new direction is in
182:        VEC_VV(1+loc_it)*/
183:     (*lgmres->orthog)(ksp,loc_it);

185:     /* new entry in hessenburg is the 2-norm of our new direction */
186:     VecNorm(VEC_VV(loc_it+1),NORM_2,&tt);

188:     *HH(loc_it+1,loc_it)  = tt;
189:     *HES(loc_it+1,loc_it) = tt;

191:     /* check for the happy breakdown */
192:     hapbnd = PetscAbsScalar(tt / *GRS(loc_it)); /* GRS(loc_it) contains the res_norm from the last iteration  */
193:     if (hapbnd > lgmres->haptol) hapbnd = lgmres->haptol;
194:     if (tt > hapbnd) {
195:       tmp  = 1.0/tt;
196:       VecScale(VEC_VV(loc_it+1),tmp); /* scale new direction by its norm */
197:     } else {
198:       PetscInfo(ksp,"Detected happy breakdown, current hapbnd = %g tt = %g\n",(double)hapbnd,(double)tt);
199:       hapend = PETSC_TRUE;
200:     }

202:     /* Now apply rotations to new col of hessenberg (and right side of system),
203:        calculate new rotation, and get new residual norm at the same time*/
204:     KSPLGMRESUpdateHessenberg(ksp,loc_it,hapend,&res);
205:     if (ksp->reason) break;

207:     loc_it++;
208:     lgmres->it = (loc_it-1);   /* Add this here in case it has converged */

210:     PetscObjectSAWsTakeAccess((PetscObject)ksp);
211:     ksp->its++;
212:     if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = res;
213:     else ksp->rnorm = 0.0;
214:     PetscObjectSAWsGrantAccess((PetscObject)ksp);

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

218:     /* Catch error in happy breakdown and signal convergence and break from loop */
219:     if (hapend) {
220:       if (!ksp->reason) {
222:         else {
223:           ksp->reason = KSP_DIVERGED_BREAKDOWN;
224:           break;
225:         }
226:       }
227:     }
228:   }
229:   /* END OF ITERATION LOOP */
230:   KSPLogResidualHistory(ksp,res);

232:   /* Monitor if we know that we will not return for a restart */
233:   if (ksp->reason || ksp->its >= max_it) {
234:     KSPMonitor(ksp, ksp->its, res);
235:   }

237:   if (itcount) *itcount = loc_it;

239:   /*
240:     Down here we have to solve for the "best" coefficients of the Krylov
241:     columns, add the solution values together, and possibly unwind the
242:     preconditioning from the solution
243:    */

245:   /* Form the solution (or the solution so far) */
246:   /* Note: must pass in (loc_it-1) for iteration count so that KSPLGMRESBuildSoln
247:      properly navigates */

249:   KSPLGMRESBuildSoln(GRS(0),ksp->vec_sol,ksp->vec_sol,ksp,loc_it-1);

251:   /* LGMRES_MOD collect aug vector and A*augvector for future restarts -
252:      only if we will be restarting (i.e. this cycle performed it_total
253:      iterations)  */
254:   if (!ksp->reason && ksp->its < max_it && aug_dim > 0) {

256:     /*AUG_TEMP contains the new augmentation vector (assigned in  KSPLGMRESBuildSoln) */
257:     if (!lgmres->aug_ct) {
258:       spot = 0;
259:       lgmres->aug_ct++;
260:     } else if (lgmres->aug_ct < aug_dim) {
261:       spot = lgmres->aug_ct;
262:       lgmres->aug_ct++;
263:     } else { /* truncate */
264:       for (ii=0; ii<aug_dim; ii++) {
265:         if (lgmres->aug_order[ii] == aug_dim) spot = ii;
266:       }
267:     }

269:     VecCopy(AUG_TEMP, AUGVEC(spot));
270:     /*need to normalize */
271:     VecNorm(AUGVEC(spot), NORM_2, &tmp_norm);

273:     inv_tmp_norm = 1.0/tmp_norm;

275:     VecScale(AUGVEC(spot),inv_tmp_norm);

277:     /*set new aug vector to order 1  - move all others back one */
278:     for (ii=0; ii < aug_dim; ii++) AUG_ORDER(ii)++;
279:     AUG_ORDER(spot) = 1;

281:     /*now add the A*aug vector to A_AUGVEC(spot)  - this is independ. of preconditioning type*/
282:     /* want V*H*y - y is in GRS, V is in VEC_VV and H is in HES */

284:     /* first do H+*y */
285:     avec = lgmres->hwork;
286:     PetscArrayzero(avec,it_total+1);
287:     for (ii=0; ii < it_total + 1; ii++) {
288:       for (jj=0; jj <= ii+1 && jj < it_total+1; jj++) {
289:         avec[jj] += *HES(jj ,ii) * *GRS(ii);
290:       }
291:     }

293:     /*now multiply result by V+ */
294:     VecSet(VEC_TEMP,0.0);
295:     VecMAXPY(VEC_TEMP, it_total+1, avec, &VEC_VV(0)); /*answer is in VEC_TEMP*/

297:     /*copy answer to aug location  and scale*/
298:     VecCopy(VEC_TEMP,  A_AUGVEC(spot));
299:     VecScale(A_AUGVEC(spot),inv_tmp_norm);
300:   }
301:   return 0;
302: }

304: /*
305:     KSPSolve_LGMRES - This routine applies the LGMRES method.

307:    Input Parameter:
308: .     ksp - the Krylov space object that was set to use lgmres

310:    Output Parameter:
311: .     outits - number of iterations used

313: */

315: PetscErrorCode KSPSolve_LGMRES(KSP ksp)
316: {
317:   PetscInt       cycle_its; /* iterations done in a call to KSPLGMRESCycle */
318:   PetscInt       itcount;   /* running total of iterations, incl. those in restarts */
319:   KSP_LGMRES     *lgmres    = (KSP_LGMRES*)ksp->data;
320:   PetscBool      guess_zero = ksp->guess_zero;
321:   PetscInt       ii;        /*LGMRES_MOD variable */


325:   PetscObjectSAWsTakeAccess((PetscObject)ksp);

327:   ksp->its        = 0;
328:   lgmres->aug_ct  = 0;
329:   lgmres->matvecs = 0;

331:   PetscObjectSAWsGrantAccess((PetscObject)ksp);

333:   /* initialize */
334:   itcount     = 0;
335:   /*LGMRES_MOD*/
336:   for (ii=0; ii<lgmres->aug_dim; ii++) lgmres->aug_order[ii] = 0;

338:   while (!ksp->reason) {
339:     /* calc residual - puts in VEC_VV(0) */
340:     KSPInitialResidual(ksp,ksp->vec_sol,VEC_TEMP,VEC_TEMP_MATOP,VEC_VV(0),ksp->vec_rhs);
341:     KSPLGMRESCycle(&cycle_its,ksp);
342:     itcount += cycle_its;
343:     if (itcount >= ksp->max_it) {
344:       if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
345:       break;
346:     }
347:     ksp->guess_zero = PETSC_FALSE; /* every future call to KSPInitialResidual() will have nonzero guess */
348:   }
349:   ksp->guess_zero = guess_zero; /* restore if user provided nonzero initial guess */
350:   return 0;
351: }

353: /*

355:    KSPDestroy_LGMRES - Frees all memory space used by the Krylov method.

357: */
358: PetscErrorCode KSPDestroy_LGMRES(KSP ksp)
359: {
360:   KSP_LGMRES     *lgmres = (KSP_LGMRES*)ksp->data;

362:   PetscFree(lgmres->augvecs);
363:   if (lgmres->augwork_alloc) {
364:     VecDestroyVecs(lgmres->augwork_alloc,&lgmres->augvecs_user_work[0]);
365:   }
366:   PetscFree(lgmres->augvecs_user_work);
367:   PetscFree(lgmres->aug_order);
368:   PetscFree(lgmres->hwork);
369:   KSPDestroy_GMRES(ksp);
370:   return 0;
371: }

373: /*
374:     KSPLGMRESBuildSoln - create the solution from the starting vector and the
375:                       current iterates.

377:     Input parameters:
378:         nrs - work area of size it + 1.
379:         vguess  - index of initial guess
380:         vdest - index of result.  Note that vguess may == vdest (replace
381:                 guess with the solution).
382:         it - HH upper triangular part is a block of size (it+1) x (it+1)

384:      This is an internal routine that knows about the LGMRES internals.
385:  */
386: static PetscErrorCode KSPLGMRESBuildSoln(PetscScalar *nrs,Vec vguess,Vec vdest,KSP ksp,PetscInt it)
387: {
388:   PetscScalar    tt;
389:   PetscInt       ii,k,j;
390:   KSP_LGMRES     *lgmres = (KSP_LGMRES*)(ksp->data);
391:   /*LGMRES_MOD */
392:   PetscInt it_arnoldi, it_aug;
393:   PetscInt jj, spot = 0;

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

397:   /* If it is < 0, no lgmres steps have been performed */
398:   if (it < 0) {
399:     VecCopy(vguess,vdest); /* VecCopy() is smart, exists immediately if vguess == vdest */
400:     return 0;
401:   }

403:   /* so (it+1) lgmres steps HAVE been performed */

405:   /* LGMRES_MOD - determine if we need to use augvecs for the soln  - do not assume that
406:      this is called after the total its allowed for an approx space */
407:   if (lgmres->approx_constant) {
408:     it_arnoldi = lgmres->max_k - lgmres->aug_ct;
409:   } else {
410:     it_arnoldi = lgmres->max_k - lgmres->aug_dim;
411:   }
412:   if (it_arnoldi >= it +1) {
413:     it_aug     = 0;
414:     it_arnoldi = it+1;
415:   } else {
416:     it_aug = (it + 1) - it_arnoldi;
417:   }

419:   /* now it_arnoldi indicates the number of matvecs that took place */
420:   lgmres->matvecs += it_arnoldi;

422:   /* solve the upper triangular system - GRS is the right side and HH is
423:      the upper triangular matrix  - put soln in nrs */
425:   if (*HH(it,it) != 0.0) {
426:     nrs[it] = *GRS(it) / *HH(it,it);
427:   } else {
428:     nrs[it] = 0.0;
429:   }

431:   for (ii=1; ii<=it; ii++) {
432:     k  = it - ii;
433:     tt = *GRS(k);
434:     for (j=k+1; j<=it; j++) tt = tt - *HH(k,j) * nrs[j];
435:     nrs[k] = tt / *HH(k,k);
436:   }

438:   /* Accumulate the correction to the soln of the preconditioned prob. in VEC_TEMP */
439:   VecSet(VEC_TEMP,0.0); /* set VEC_TEMP components to 0 */

441:   /*LGMRES_MOD - if augmenting has happened we need to form the solution
442:     using the augvecs */
443:   if (!it_aug) { /* all its are from arnoldi */
444:     VecMAXPY(VEC_TEMP,it+1,nrs,&VEC_VV(0));
445:   } else { /*use aug vecs */
446:     /*first do regular krylov directions */
447:     VecMAXPY(VEC_TEMP,it_arnoldi,nrs,&VEC_VV(0));
448:     /*now add augmented portions - add contribution of aug vectors one at a time*/

450:     for (ii=0; ii<it_aug; ii++) {
451:       for (jj=0; jj<lgmres->aug_dim; jj++) {
452:         if (lgmres->aug_order[jj] == (ii+1)) {
453:           spot = jj;
454:           break; /* must have this because there will be duplicates before aug_ct = aug_dim */
455:         }
456:       }
457:       VecAXPY(VEC_TEMP,nrs[it_arnoldi+ii],AUGVEC(spot));
458:     }
459:   }
460:   /* now VEC_TEMP is what we want to keep for augmenting purposes - grab before the
461:      preconditioner is "unwound" from right-precondtioning*/
462:   VecCopy(VEC_TEMP, AUG_TEMP);

464:   KSPUnwindPreconditioner(ksp,VEC_TEMP,VEC_TEMP_MATOP);

466:   /* add solution to previous solution */
467:   /* put updated solution into vdest.*/
468:   VecCopy(vguess,vdest);
469:   VecAXPY(vdest,1.0,VEC_TEMP);
470:   return 0;
471: }

473: /*

475:     KSPLGMRESUpdateHessenberg - Do the scalar work for the orthogonalization.
476:                             Return new residual.

478:     input parameters:

480: .        ksp -    Krylov space object
481: .        it  -    plane rotations are applied to the (it+1)th column of the
482:                   modified hessenberg (i.e. HH(:,it))
483: .        hapend - PETSC_FALSE not happy breakdown ending.

485:     output parameters:
486: .        res - the new residual

488:  */
489: static PetscErrorCode KSPLGMRESUpdateHessenberg(KSP ksp,PetscInt it,PetscBool hapend,PetscReal *res)
490: {
491:   PetscScalar *hh,*cc,*ss,tt;
492:   PetscInt    j;
493:   KSP_LGMRES  *lgmres = (KSP_LGMRES*)(ksp->data);

495:   hh = HH(0,it);   /* pointer to beginning of column to update - so
496:                       incrementing hh "steps down" the (it+1)th col of HH*/
497:   cc = CC(0);      /* beginning of cosine rotations */
498:   ss = SS(0);      /* beginning of sine rotations */

500:   /* Apply all the previously computed plane rotations to the new column
501:      of the Hessenberg matrix */
502:   /* Note: this uses the rotation [conj(c)  s ; -s   c], c= cos(theta), s= sin(theta) */

504:   for (j=1; j<=it; j++) {
505:     tt  = *hh;
506:     *hh = PetscConj(*cc) * tt + *ss * *(hh+1);
507:     hh++;
508:     *hh = *cc++ * *hh - (*ss++ * tt);
509:     /* hh, cc, and ss have all been incremented one by end of loop */
510:   }

512:   /*
513:     compute the new plane rotation, and apply it to:
514:      1) the right-hand-side of the Hessenberg system (GRS)
515:         note: it affects GRS(it) and GRS(it+1)
516:      2) the new column of the Hessenberg matrix
517:         note: it affects HH(it,it) which is currently pointed to
518:         by hh and HH(it+1, it) (*(hh+1))
519:     thus obtaining the updated value of the residual...
520:   */

522:   /* compute new plane rotation */

524:   if (!hapend) {
525:     tt = PetscSqrtScalar(PetscConj(*hh) * *hh + PetscConj(*(hh+1)) * *(hh+1));
526:     if (tt == 0.0) {
527:       ksp->reason = KSP_DIVERGED_NULL;
528:       return 0;
529:     }
530:     *cc = *hh / tt;         /* new cosine value */
531:     *ss = *(hh+1) / tt;        /* new sine value */

533:     /* apply to 1) and 2) */
534:     *GRS(it+1) = -(*ss * *GRS(it));
535:     *GRS(it)   = PetscConj(*cc) * *GRS(it);
536:     *hh        = PetscConj(*cc) * *hh + *ss * *(hh+1);

538:     /* residual is the last element (it+1) of right-hand side! */
539:     *res = PetscAbsScalar(*GRS(it+1));

541:   } else { /* happy breakdown: HH(it+1, it) = 0, therefore we don't need to apply
542:             another rotation matrix (so RH doesn't change).  The new residual is
543:             always the new sine term times the residual from last time (GRS(it)),
544:             but now the new sine rotation would be zero...so the residual should
545:             be zero...so we will multiply "zero" by the last residual.  This might
546:             not be exactly what we want to do here -could just return "zero". */

548:     *res = 0.0;
549:   }
550:   return 0;
551: }

553: /*

555:    KSPLGMRESGetNewVectors - This routine allocates more work vectors, starting from
556:                          VEC_VV(it)

558: */
559: static PetscErrorCode KSPLGMRESGetNewVectors(KSP ksp,PetscInt it)
560: {
561:   KSP_LGMRES     *lgmres = (KSP_LGMRES*)ksp->data;
562:   PetscInt       nwork   = lgmres->nwork_alloc; /* number of work vector chunks allocated */
563:   PetscInt       nalloc;                      /* number to allocate */
564:   PetscInt       k;

566:   nalloc = lgmres->delta_allocate; /* number of vectors to allocate
567:                                       in a single chunk */

569:   /* Adjust the number to allocate to make sure that we don't exceed the
570:      number of available slots (lgmres->vecs_allocated)*/
571:   if (it + VEC_OFFSET + nalloc >= lgmres->vecs_allocated) {
572:     nalloc = lgmres->vecs_allocated - it - VEC_OFFSET;
573:   }
574:   if (!nalloc) return 0;

576:   lgmres->vv_allocated += nalloc; /* vv_allocated is the number of vectors allocated */

578:   /* work vectors */
579:   KSPCreateVecs(ksp,nalloc,&lgmres->user_work[nwork],0,NULL);
580:   PetscLogObjectParents(ksp,nalloc,lgmres->user_work[nwork]);
581:   /* specify size of chunk allocated */
582:   lgmres->mwork_alloc[nwork] = nalloc;

584:   for (k=0; k < nalloc; k++) {
585:     lgmres->vecs[it+VEC_OFFSET+k] = lgmres->user_work[nwork][k];
586:   }

588:   /* LGMRES_MOD - for now we are preallocating the augmentation vectors */

590:   /* increment the number of work vector chunks */
591:   lgmres->nwork_alloc++;
592:   return 0;
593: }

595: /*

597:    KSPBuildSolution_LGMRES

599:      Input Parameter:
600: .     ksp - the Krylov space object
601: .     ptr-

603:    Output Parameter:
604: .     result - the solution

606:    Note: this calls KSPLGMRESBuildSoln - the same function that KSPLGMRESCycle
607:    calls directly.

609: */
610: PetscErrorCode KSPBuildSolution_LGMRES(KSP ksp,Vec ptr,Vec *result)
611: {
612:   KSP_LGMRES     *lgmres = (KSP_LGMRES*)ksp->data;

614:   if (!ptr) {
615:     if (!lgmres->sol_temp) {
616:       VecDuplicate(ksp->vec_sol,&lgmres->sol_temp);
617:       PetscLogObjectParent((PetscObject)ksp,(PetscObject)lgmres->sol_temp);
618:     }
619:     ptr = lgmres->sol_temp;
620:   }
621:   if (!lgmres->nrs) {
622:     /* allocate the work area */
623:     PetscMalloc1(lgmres->max_k,&lgmres->nrs);
624:     PetscLogObjectMemory((PetscObject)ksp,lgmres->max_k*sizeof(PetscScalar));
625:   }

627:   KSPLGMRESBuildSoln(lgmres->nrs,ksp->vec_sol,ptr,ksp,lgmres->it);
628:   if (result) *result = ptr;
629:   return 0;
630: }

632: PetscErrorCode KSPView_LGMRES(KSP ksp,PetscViewer viewer)
633: {
634:   KSP_LGMRES     *lgmres = (KSP_LGMRES*)ksp->data;
635:   PetscBool      iascii;

637:   KSPView_GMRES(ksp,viewer);
638:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
639:   if (iascii) {
640:     /*LGMRES_MOD */
641:     PetscViewerASCIIPrintf(viewer,"  aug. dimension=%D\n",lgmres->aug_dim);
642:     if (lgmres->approx_constant) {
643:       PetscViewerASCIIPrintf(viewer,"  approx. space size was kept constant.\n");
644:     }
645:     PetscViewerASCIIPrintf(viewer,"  number of matvecs=%D\n",lgmres->matvecs);
646:   }
647:   return 0;
648: }

650: PetscErrorCode KSPSetFromOptions_LGMRES(PetscOptionItems *PetscOptionsObject,KSP ksp)
651: {
652:   PetscInt       aug;
653:   KSP_LGMRES     *lgmres = (KSP_LGMRES*) ksp->data;
654:   PetscBool      flg     = PETSC_FALSE;

656:   KSPSetFromOptions_GMRES(PetscOptionsObject,ksp);
657:   PetscOptionsHead(PetscOptionsObject,"KSP LGMRES Options");
658:   PetscOptionsBool("-ksp_lgmres_constant","Use constant approx. space size","KSPGMRESSetConstant",lgmres->approx_constant,&lgmres->approx_constant,NULL);
659:   PetscOptionsInt("-ksp_lgmres_augment","Number of error approximations to augment the Krylov space with","KSPLGMRESSetAugDim",lgmres->aug_dim,&aug,&flg);
660:   if (flg) KSPLGMRESSetAugDim(ksp,aug);
661:   PetscOptionsTail();
662:   return 0;
663: }

665: /*functions for extra lgmres options here*/
666: static PetscErrorCode  KSPLGMRESSetConstant_LGMRES(KSP ksp)
667: {
668:   KSP_LGMRES *lgmres = (KSP_LGMRES*)ksp->data;

670:   lgmres->approx_constant = PETSC_TRUE;
671:   return 0;
672: }

674: static PetscErrorCode  KSPLGMRESSetAugDim_LGMRES(KSP ksp,PetscInt aug_dim)
675: {
676:   KSP_LGMRES *lgmres = (KSP_LGMRES*)ksp->data;

680:   lgmres->aug_dim = aug_dim;
681:   return 0;
682: }

684: /* end new lgmres functions */

686: /*MC
687:     KSPLGMRES - Augments the standard GMRES approximation space with approximations to
688:                 the error from previous restart cycles.

690:   Options Database Keys:
691: +   -ksp_gmres_restart <restart> - total approximation space size (Krylov directions + error approximations)
692: .   -ksp_gmres_haptol <tol> - sets the tolerance for "happy ending" (exact convergence)
693: .   -ksp_gmres_preallocate - preallocate all the Krylov search directions initially (otherwise groups of
694:                             vectors are allocated as needed)
695: .   -ksp_gmres_classicalgramschmidt - use classical (unmodified) Gram-Schmidt to orthogonalize against the Krylov space (fast) (the default)
696: .   -ksp_gmres_modifiedgramschmidt - use modified Gram-Schmidt in the orthogonalization (more stable, but slower)
697: .   -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always> - determine if iterative refinement is used to increase the
698:                                   stability of the classical Gram-Schmidt  orthogonalization.
699: .   -ksp_gmres_krylov_monitor - plot the Krylov space generated
700: .   -ksp_lgmres_augment <k> - number of error approximations to augment the Krylov space with
701: -   -ksp_lgmres_constant - use a constant approx. space size (only affects restart cycles < num. error approx.(k), i.e. the first k restarts)

703:     To run LGMRES(m, k) as described in the above paper, use:
704:        -ksp_gmres_restart <m+k>
705:        -ksp_lgmres_augment <k>

707:   Level: beginner

709:    Notes:
710:     Supports both left and right preconditioning, but not symmetric.

712:    References:
713: .  * - A. H. Baker, E.R. Jessup, and T.A. Manteuffel. A technique for accelerating the convergence of restarted GMRES. SIAM Journal on Matrix Analysis and Applications, 26 (2005).

715:   Developer Notes:
716:     This object is subclassed off of KSPGMRES

718:   Contributed by: Allison Baker

720: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPFGMRES, KSPGMRES,
721:           KSPGMRESSetRestart(), KSPGMRESSetHapTol(), KSPGMRESSetPreAllocateVectors(), KSPGMRESSetOrthogonalization(), KSPGMRESGetOrthogonalization(),
722:           KSPGMRESClassicalGramSchmidtOrthogonalization(), KSPGMRESModifiedGramSchmidtOrthogonalization(),
723:           KSPGMRESCGSRefinementType, KSPGMRESSetCGSRefinementType(), KSPGMRESGetCGSRefinementType(), KSPGMRESMonitorKrylov(), KSPLGMRESSetAugDim(),
724:           KSPGMRESSetConstant()

726: M*/

728: PETSC_EXTERN PetscErrorCode KSPCreate_LGMRES(KSP ksp)
729: {
730:   KSP_LGMRES     *lgmres;

732:   PetscNewLog(ksp,&lgmres);

734:   ksp->data               = (void*)lgmres;
735:   ksp->ops->buildsolution = KSPBuildSolution_LGMRES;

737:   ksp->ops->setup                        = KSPSetUp_LGMRES;
738:   ksp->ops->solve                        = KSPSolve_LGMRES;
739:   ksp->ops->destroy                      = KSPDestroy_LGMRES;
740:   ksp->ops->view                         = KSPView_LGMRES;
741:   ksp->ops->setfromoptions               = KSPSetFromOptions_LGMRES;
742:   ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_GMRES;
743:   ksp->ops->computeeigenvalues           = KSPComputeEigenvalues_GMRES;

745:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
746:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);
747:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);

749:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetPreAllocateVectors_C",KSPGMRESSetPreAllocateVectors_GMRES);
750:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetOrthogonalization_C",KSPGMRESSetOrthogonalization_GMRES);
751:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetOrthogonalization_C",KSPGMRESGetOrthogonalization_GMRES);
752:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetRestart_C",KSPGMRESSetRestart_GMRES);
753:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetRestart_C",KSPGMRESGetRestart_GMRES);
754:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetHapTol_C",KSPGMRESSetHapTol_GMRES);
755:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetCGSRefinementType_C",KSPGMRESSetCGSRefinementType_GMRES);
756:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetCGSRefinementType_C",KSPGMRESGetCGSRefinementType_GMRES);

758:   /*LGMRES_MOD add extra functions here - like the one to set num of aug vectors */
759:   PetscObjectComposeFunction((PetscObject)ksp,"KSPLGMRESSetConstant_C",KSPLGMRESSetConstant_LGMRES);
760:   PetscObjectComposeFunction((PetscObject)ksp,"KSPLGMRESSetAugDim_C",KSPLGMRESSetAugDim_LGMRES);

762:   /*defaults */
763:   lgmres->haptol         = 1.0e-30;
764:   lgmres->q_preallocate  = 0;
765:   lgmres->delta_allocate = LGMRES_DELTA_DIRECTIONS;
766:   lgmres->orthog         = KSPGMRESClassicalGramSchmidtOrthogonalization;
767:   lgmres->nrs            = NULL;
768:   lgmres->sol_temp       = NULL;
769:   lgmres->max_k          = LGMRES_DEFAULT_MAXK;
770:   lgmres->Rsvd           = NULL;
771:   lgmres->cgstype        = KSP_GMRES_CGS_REFINE_NEVER;
772:   lgmres->orthogwork     = NULL;

774:   /*LGMRES_MOD - new defaults */
775:   lgmres->aug_dim         = LGMRES_DEFAULT_AUGDIM;
776:   lgmres->aug_ct          = 0;     /* start with no aug vectors */
777:   lgmres->approx_constant = PETSC_FALSE;
778:   lgmres->matvecs         = 0;
779:   return 0;
780: }