Actual source code: lgmres.c

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

  3: static PetscErrorCode KSPLGMRESGetNewVectors(KSP, PetscInt);
  4: static PetscErrorCode KSPLGMRESUpdateHessenberg(KSP, PetscInt, PetscBool, PetscReal *);
  5: static PetscErrorCode KSPLGMRESBuildSoln(PetscScalar *, Vec, Vec, KSP, PetscInt);

  7: /*@
  8:   KSPLGMRESSetAugDim - Set the number of error approximations to include in the approximation space (default is 2) for `KSPLGMRES`

 10:   Collective

 12:   Input Parameters:
 13: + ksp - the `KSP` context
 14: - dim - the number of vectors to use

 16:   Options Database Key:
 17: . -ksp_lgmres_augment dim - the number of error approximations to include

 19:   Level: intermediate

 21:   Note:
 22:   If this is set to zero, then this method is equivalent to `KSPGMRES`

 24: .seealso: [](ch_ksp), `KSPLGMRES`, `KSPLGMRESSetConstant()`
 25: @*/
 26: PetscErrorCode KSPLGMRESSetAugDim(KSP ksp, PetscInt dim)
 27: {
 28:   PetscFunctionBegin;
 29:   PetscTryMethod((ksp), "KSPLGMRESSetAugDim_C", (KSP, PetscInt), (ksp, dim));
 30:   PetscFunctionReturn(PETSC_SUCCESS);
 31: }

 33: /*@
 34:   KSPLGMRESSetConstant - keep the error approximation space a constant size for every restart cycle

 36:   Collective

 38:   Input Parameters:
 39: . ksp - the `KSP` context

 41:   Options Database Key:
 42: . -ksp_lgmres_constant - set the size to be constant

 44:   Level: intermediate

 46:   Note:
 47:   This only affects the first couple of restart cycles when the total number of desired error approximations may not
 48:   be available.

 50: .seealso: [](ch_ksp), `KSPLGMRES`, `KSPLGMRESSetAugDim()`
 51: @*/
 52: PetscErrorCode KSPLGMRESSetConstant(KSP ksp)
 53: {
 54:   PetscFunctionBegin;
 55:   PetscTryMethod((ksp), "KSPLGMRESSetConstant_C", (KSP), (ksp));
 56:   PetscFunctionReturn(PETSC_SUCCESS);
 57: }

 59: static PetscErrorCode KSPSetUp_LGMRES(KSP ksp)
 60: {
 61:   PetscInt    max_k, k, aug_dim;
 62:   KSP_LGMRES *lgmres = (KSP_LGMRES *)ksp->data;

 64:   PetscFunctionBegin;
 65:   max_k   = lgmres->max_k;
 66:   aug_dim = lgmres->aug_dim;
 67:   PetscCall(KSPSetUp_GMRES(ksp));

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

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

 74:   PetscCall(PetscMalloc1(2 * aug_dim + AUG_OFFSET, &lgmres->augvecs_user_work));
 75:   PetscCall(PetscMalloc1(aug_dim, &lgmres->aug_order));

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

 82:   PetscCall(KSPCreateVecs(ksp, lgmres->aug_vv_allocated, &lgmres->augvecs_user_work[0], 0, NULL));
 83:   PetscCall(PetscMalloc1(max_k + 1, &lgmres->hwork));
 84:   for (k = 0; k < lgmres->aug_vv_allocated; k++) lgmres->augvecs[k] = lgmres->augvecs_user_work[0][k];
 85:   PetscFunctionReturn(PETSC_SUCCESS);
 86: }

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

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

110:   PetscFunctionBegin;
111:   /* Number of pseudo iterations since last restart is the number of prestart directions */
112:   loc_it = 0;

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

120:   it_total = it_arnoldi + lgmres->aug_ct;

122:   /* initial residual is in VEC_VV(0)  - compute its norm*/
123:   PetscCall(VecNorm(VEC_VV(0), NORM_2, &res_norm));
124:   KSPCheckNorm(ksp, res_norm);
125:   res = res_norm;

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

130:   /* check for the convergence */
131:   if (!res) {
132:     if (itcount) *itcount = 0;
133:     ksp->reason = KSP_CONVERGED_ATOL;
134:     PetscCall(PetscInfo(ksp, "Converged due to zero residual norm on entry\n"));
135:     PetscFunctionReturn(PETSC_SUCCESS);
136:   }

138:   /* scale VEC_VV (the initial residual) */
139:   tmp = 1.0 / res_norm;
140:   PetscCall(VecScale(VEC_VV(0), tmp));

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

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

151:   /* MAIN ITERATION LOOP BEGINNING*/

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

157:   while (!ksp->reason && loc_it < it_total && ksp->its < max_it) { /* LGMRES_MOD: changed to it_total */
158:     PetscCall(KSPLogResidualHistory(ksp, res));
159:     lgmres->it = (loc_it - 1);
160:     PetscCall(KSPMonitor(ksp, ksp->its, res));

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

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

180:       PetscCall(VecCopy(A_AUGVEC(spot), VEC_VV(1 + loc_it)));
181:       /* note: an alternate implementation choice would be to only save the AUGVECS and
182:          not A_AUGVEC and then apply the PC here to the augvec */
183:     }

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

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

191:     *HH(loc_it + 1, loc_it)  = tt;
192:     *HES(loc_it + 1, loc_it) = tt;

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

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

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

213:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
214:     ksp->its++;
215:     if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = res;
216:     else ksp->rnorm = 0.0;
217:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));

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

221:     /* Catch error in happy breakdown and signal convergence and break from loop */
222:     if (hapend) {
223:       if (!ksp->reason) {
224:         PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "Reached happy break down, but convergence was not indicated. Residual norm = %g", (double)res);
225:         ksp->reason = KSP_DIVERGED_BREAKDOWN;
226:         break;
227:       }
228:     }
229:   }
230:   /* END OF ITERATION LOOP */
231:   PetscCall(KSPLogResidualHistory(ksp, res));

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

236:   if (itcount) *itcount = loc_it;

238:   /*
239:     Solve for the "best" coefficients of the Krylov
240:     columns, add the solution values together, and possibly unwind the
241:     preconditioning from the solution
242:    */

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

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

249:   /* LGMRES_MOD collect aug vector and A*augvector for future restarts -
250:      only if we will be restarting (i.e. this cycle performed it_total iterations)  */
251:   if (!ksp->reason && ksp->its < max_it && aug_dim > 0) {
252:     /* AUG_TEMP contains the new augmentation vector (assigned in  KSPLGMRESBuildSoln) */
253:     if (!lgmres->aug_ct) {
254:       spot = 0;
255:       lgmres->aug_ct++;
256:     } else if (lgmres->aug_ct < aug_dim) {
257:       spot = lgmres->aug_ct;
258:       lgmres->aug_ct++;
259:     } else { /* truncate */
260:       for (ii = 0; ii < aug_dim; ii++) {
261:         if (lgmres->aug_order[ii] == aug_dim) spot = ii;
262:       }
263:     }

265:     PetscCall(VecCopy(AUG_TEMP, AUGVEC(spot)));
266:     /* need to normalize */
267:     PetscCall(VecNorm(AUGVEC(spot), NORM_2, &tmp_norm));

269:     inv_tmp_norm = 1.0 / tmp_norm;

271:     PetscCall(VecScale(AUGVEC(spot), inv_tmp_norm));

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

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

280:     /* do H+*y */
281:     avec = lgmres->hwork;
282:     PetscCall(PetscArrayzero(avec, it_total + 1));
283:     for (ii = 0; ii < it_total + 1; ii++) {
284:       for (jj = 0; jj <= ii + 1 && jj < it_total + 1; jj++) avec[jj] += *HES(jj, ii) * *GRS(ii);
285:     }

287:     /* multiply result by V+ */
288:     PetscCall(VecMAXPBY(VEC_TEMP, it_total + 1, avec, 0, &VEC_VV(0))); /* answer is in VEC_TEMP */

290:     /* copy answer to aug location  and scale */
291:     PetscCall(VecCopy(VEC_TEMP, A_AUGVEC(spot)));
292:     PetscCall(VecScale(A_AUGVEC(spot), inv_tmp_norm));
293:   }
294:   PetscFunctionReturn(PETSC_SUCCESS);
295: }

297: static PetscErrorCode KSPSolve_LGMRES(KSP ksp)
298: {
299:   PetscInt    itcount; /* running total of iterations, incl. those in restarts */
300:   KSP_LGMRES *lgmres     = (KSP_LGMRES *)ksp->data;
301:   PetscBool   guess_zero = ksp->guess_zero;
302:   PetscInt    ii; /* LGMRES_MOD variable */

304:   PetscFunctionBegin;
305:   PetscCheck(!ksp->calc_sings || lgmres->Rsvd, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ORDER, "Must call KSPSetComputeSingularValues() before KSPSetUp() is called");

307:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));

309:   ksp->its        = 0;
310:   lgmres->aug_ct  = 0;
311:   lgmres->matvecs = 0;

313:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));

315:   /* initialize */
316:   itcount = 0;
317:   /* LGMRES_MOD */
318:   for (ii = 0; ii < lgmres->aug_dim; ii++) lgmres->aug_order[ii] = 0;

320:   while (!ksp->reason) {
321:     PetscInt cycle_its = 0; /* iterations done in a call to KSPLGMRESCycle */
322:     /* calc residual - puts in VEC_VV(0) */
323:     PetscCall(KSPInitialResidual(ksp, ksp->vec_sol, VEC_TEMP, VEC_TEMP_MATOP, VEC_VV(0), ksp->vec_rhs));
324:     PetscCall(KSPLGMRESCycle(&cycle_its, ksp));
325:     itcount += cycle_its;
326:     if (itcount >= ksp->max_it) {
327:       if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
328:       break;
329:     }
330:     ksp->guess_zero = PETSC_FALSE; /* every future call to KSPInitialResidual() will have nonzero guess */
331:   }
332:   ksp->guess_zero = guess_zero; /* restore if user provided nonzero initial guess */
333:   PetscFunctionReturn(PETSC_SUCCESS);
334: }

336: static PetscErrorCode KSPDestroy_LGMRES(KSP ksp)
337: {
338:   KSP_LGMRES *lgmres = (KSP_LGMRES *)ksp->data;

340:   PetscFunctionBegin;
341:   PetscCall(PetscFree(lgmres->augvecs));
342:   if (lgmres->augwork_alloc) PetscCall(VecDestroyVecs(lgmres->augwork_alloc, &lgmres->augvecs_user_work[0]));
343:   PetscCall(PetscFree(lgmres->augvecs_user_work));
344:   PetscCall(PetscFree(lgmres->aug_order));
345:   PetscCall(PetscFree(lgmres->hwork));
346:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPLGMRESSetConstant_C", NULL));
347:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPLGMRESSetAugDim_C", NULL));
348:   PetscCall(KSPDestroy_GMRES(ksp));
349:   PetscFunctionReturn(PETSC_SUCCESS);
350: }

352: static PetscErrorCode KSPLGMRESBuildSoln(PetscScalar *nrs, Vec vguess, Vec vdest, KSP ksp, PetscInt it)
353: {
354:   PetscScalar tt;
355:   PetscInt    ii, k, j;
356:   KSP_LGMRES *lgmres = (KSP_LGMRES *)(ksp->data);
357:   /* LGMRES_MOD */
358:   PetscInt it_arnoldi, it_aug;
359:   PetscInt jj, spot = 0;

361:   PetscFunctionBegin;
362:   /* Solve for solution vector that minimizes the residual */

364:   /* If it is < 0, no lgmres steps have been performed */
365:   if (it < 0) {
366:     PetscCall(VecCopy(vguess, vdest)); /* VecCopy() is smart, exists immediately if vguess == vdest */
367:     PetscFunctionReturn(PETSC_SUCCESS);
368:   }

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

372:   /* LGMRES_MOD - determine if we need to use augvecs for the soln  - do not assume that
373:      this is called after the total its allowed for an approx space */
374:   if (lgmres->approx_constant) {
375:     it_arnoldi = lgmres->max_k - lgmres->aug_ct;
376:   } else {
377:     it_arnoldi = lgmres->max_k - lgmres->aug_dim;
378:   }
379:   if (it_arnoldi >= it + 1) {
380:     it_aug     = 0;
381:     it_arnoldi = it + 1;
382:   } else {
383:     it_aug = (it + 1) - it_arnoldi;
384:   }

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

389:   /* solve the upper triangular system - GRS is the right side and HH is
390:      the upper triangular matrix  - put soln in nrs */
391:   PetscCheck(*HH(it, it) != 0.0, PETSC_COMM_SELF, PETSC_ERR_CONV_FAILED, "HH(it,it) is identically zero; it = %" PetscInt_FMT " GRS(it) = %g", it, (double)PetscAbsScalar(*GRS(it)));
392:   if (*HH(it, it) != 0.0) {
393:     nrs[it] = *GRS(it) / *HH(it, it);
394:   } else {
395:     nrs[it] = 0.0;
396:   }

398:   for (ii = 1; ii <= it; ii++) {
399:     k  = it - ii;
400:     tt = *GRS(k);
401:     for (j = k + 1; j <= it; j++) tt = tt - *HH(k, j) * nrs[j];
402:     nrs[k] = tt / *HH(k, k);
403:   }

405:   /* Accumulate the correction to the soln of the preconditioned prob. in VEC_TEMP */

407:   /* LGMRES_MOD - if augmenting has happened we need to form the solution using the augvecs */
408:   if (!it_aug) { /* all its are from arnoldi */
409:     PetscCall(VecMAXPBY(VEC_TEMP, it + 1, nrs, 0, &VEC_VV(0)));
410:   } else { /* use aug vecs */
411:     /* first do regular Krylov directions */
412:     PetscCall(VecMAXPBY(VEC_TEMP, it_arnoldi, nrs, 0, &VEC_VV(0)));
413:     /* now add augmented portions - add contribution of aug vectors one at a time*/

415:     for (ii = 0; ii < it_aug; ii++) {
416:       for (jj = 0; jj < lgmres->aug_dim; jj++) {
417:         if (lgmres->aug_order[jj] == (ii + 1)) {
418:           spot = jj;
419:           break; /* must have this because there will be duplicates before aug_ct = aug_dim */
420:         }
421:       }
422:       PetscCall(VecAXPY(VEC_TEMP, nrs[it_arnoldi + ii], AUGVEC(spot)));
423:     }
424:   }
425:   /* now VEC_TEMP is what we want to keep for augmenting purposes - grab before the
426:      preconditioner is "unwound" from right-precondtioning*/
427:   PetscCall(VecCopy(VEC_TEMP, AUG_TEMP));

429:   PetscCall(KSPUnwindPreconditioner(ksp, VEC_TEMP, VEC_TEMP_MATOP));

431:   /* add solution to previous solution */
432:   /* put updated solution into vdest.*/
433:   PetscCall(VecCopy(vguess, vdest));
434:   PetscCall(VecAXPY(vdest, 1.0, VEC_TEMP));
435:   PetscFunctionReturn(PETSC_SUCCESS);
436: }

438: static PetscErrorCode KSPLGMRESUpdateHessenberg(KSP ksp, PetscInt it, PetscBool hapend, PetscReal *res)
439: {
440:   PetscScalar *hh, *cc, *ss, tt;
441:   PetscInt     j;
442:   KSP_LGMRES  *lgmres = (KSP_LGMRES *)(ksp->data);

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

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

453:   for (j = 1; j <= it; j++) {
454:     tt  = *hh;
455:     *hh = PetscConj(*cc) * tt + *ss * *(hh + 1);
456:     hh++;
457:     *hh = *cc++ * *hh - (*ss++ * tt);
458:     /* hh, cc, and ss have all been incremented one by end of loop */
459:   }

461:   /*
462:     compute the new plane rotation, and apply it to:
463:      1) the right-hand-side of the Hessenberg system (GRS)
464:         note: it affects GRS(it) and GRS(it+1)
465:      2) the new column of the Hessenberg matrix
466:         note: it affects HH(it,it) which is currently pointed to
467:         by hh and HH(it+1, it) (*(hh+1))
468:     thus obtaining the updated value of the residual...
469:   */

471:   /* compute new plane rotation */

473:   if (!hapend) {
474:     tt = PetscSqrtScalar(PetscConj(*hh) * *hh + PetscConj(*(hh + 1)) * *(hh + 1));
475:     if (tt == 0.0) {
476:       ksp->reason = KSP_DIVERGED_NULL;
477:       PetscFunctionReturn(PETSC_SUCCESS);
478:     }
479:     *cc = *hh / tt;       /* new cosine value */
480:     *ss = *(hh + 1) / tt; /* new sine value */

482:     /* apply to 1) and 2) */
483:     *GRS(it + 1) = -(*ss * *GRS(it));
484:     *GRS(it)     = PetscConj(*cc) * *GRS(it);
485:     *hh          = PetscConj(*cc) * *hh + *ss * *(hh + 1);

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

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

497:     *res = 0.0;
498:   }
499:   PetscFunctionReturn(PETSC_SUCCESS);
500: }

502: /*
503:    KSPLGMRESGetNewVectors - Allocates more work vectors, starting from VEC_VV(it)

505: */
506: static PetscErrorCode KSPLGMRESGetNewVectors(KSP ksp, PetscInt it)
507: {
508:   KSP_LGMRES *lgmres = (KSP_LGMRES *)ksp->data;
509:   PetscInt    nwork  = lgmres->nwork_alloc; /* number of work vector chunks allocated */
510:   PetscInt    nalloc;                       /* number to allocate */
511:   PetscInt    k;

513:   PetscFunctionBegin;
514:   nalloc = lgmres->delta_allocate; /* number of vectors to allocate in a single chunk */

516:   /* Adjust the number to allocate to make sure that we don't exceed the
517:      number of available slots (lgmres->vecs_allocated)*/
518:   if (it + VEC_OFFSET + nalloc >= lgmres->vecs_allocated) nalloc = lgmres->vecs_allocated - it - VEC_OFFSET;
519:   if (!nalloc) PetscFunctionReturn(PETSC_SUCCESS);

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

523:   /* work vectors */
524:   PetscCall(KSPCreateVecs(ksp, nalloc, &lgmres->user_work[nwork], 0, NULL));
525:   /* specify size of chunk allocated */
526:   lgmres->mwork_alloc[nwork] = nalloc;

528:   for (k = 0; k < nalloc; k++) lgmres->vecs[it + VEC_OFFSET + k] = lgmres->user_work[nwork][k];

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

532:   /* increment the number of work vector chunks */
533:   lgmres->nwork_alloc++;
534:   PetscFunctionReturn(PETSC_SUCCESS);
535: }

537: static PetscErrorCode KSPBuildSolution_LGMRES(KSP ksp, Vec ptr, Vec *result)
538: {
539:   KSP_LGMRES *lgmres = (KSP_LGMRES *)ksp->data;

541:   PetscFunctionBegin;
542:   if (!ptr) {
543:     if (!lgmres->sol_temp) PetscCall(VecDuplicate(ksp->vec_sol, &lgmres->sol_temp));
544:     ptr = lgmres->sol_temp;
545:   }
546:   if (!lgmres->nrs) {
547:     /* allocate the work area */
548:     PetscCall(PetscMalloc1(lgmres->max_k, &lgmres->nrs));
549:   }

551:   PetscCall(KSPLGMRESBuildSoln(lgmres->nrs, ksp->vec_sol, ptr, ksp, lgmres->it));
552:   if (result) *result = ptr;
553:   PetscFunctionReturn(PETSC_SUCCESS);
554: }

556: static PetscErrorCode KSPView_LGMRES(KSP ksp, PetscViewer viewer)
557: {
558:   KSP_LGMRES *lgmres = (KSP_LGMRES *)ksp->data;
559:   PetscBool   iascii;

561:   PetscFunctionBegin;
562:   PetscCall(KSPView_GMRES(ksp, viewer));
563:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
564:   if (iascii) {
565:     /* LGMRES_MOD */
566:     PetscCall(PetscViewerASCIIPrintf(viewer, "  aug. dimension=%" PetscInt_FMT "\n", lgmres->aug_dim));
567:     if (lgmres->approx_constant) PetscCall(PetscViewerASCIIPrintf(viewer, "  approx. space size was kept constant.\n"));
568:     PetscCall(PetscViewerASCIIPrintf(viewer, "  number of matvecs=%" PetscInt_FMT "\n", lgmres->matvecs));
569:   }
570:   PetscFunctionReturn(PETSC_SUCCESS);
571: }

573: static PetscErrorCode KSPSetFromOptions_LGMRES(KSP ksp, PetscOptionItems *PetscOptionsObject)
574: {
575:   PetscInt    aug;
576:   KSP_LGMRES *lgmres = (KSP_LGMRES *)ksp->data;
577:   PetscBool   flg    = PETSC_FALSE;

579:   PetscFunctionBegin;
580:   PetscCall(KSPSetFromOptions_GMRES(ksp, PetscOptionsObject));
581:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP LGMRES Options");
582:   PetscCall(PetscOptionsBool("-ksp_lgmres_constant", "Use constant approx. space size", "KSPGMRESSetConstant", lgmres->approx_constant, &lgmres->approx_constant, NULL));
583:   PetscCall(PetscOptionsInt("-ksp_lgmres_augment", "Number of error approximations to augment the Krylov space with", "KSPLGMRESSetAugDim", lgmres->aug_dim, &aug, &flg));
584:   if (flg) PetscCall(KSPLGMRESSetAugDim(ksp, aug));
585:   PetscOptionsHeadEnd();
586:   PetscFunctionReturn(PETSC_SUCCESS);
587: }

589: static PetscErrorCode KSPLGMRESSetConstant_LGMRES(KSP ksp)
590: {
591:   KSP_LGMRES *lgmres = (KSP_LGMRES *)ksp->data;

593:   PetscFunctionBegin;
594:   lgmres->approx_constant = PETSC_TRUE;
595:   PetscFunctionReturn(PETSC_SUCCESS);
596: }

598: static PetscErrorCode KSPLGMRESSetAugDim_LGMRES(KSP ksp, PetscInt aug_dim)
599: {
600:   KSP_LGMRES *lgmres = (KSP_LGMRES *)ksp->data;

602:   PetscFunctionBegin;
603:   PetscCheck(aug_dim >= 0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Augmentation dimension must be nonegative");
604:   PetscCheck(aug_dim <= (lgmres->max_k - 1), PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Augmentation dimension must be <= (restart size-1)");
605:   lgmres->aug_dim = aug_dim;
606:   PetscFunctionReturn(PETSC_SUCCESS);
607: }

609: /*MC
610:   KSPLGMRES - Augments the standard GMRES approximation space with approximations to the error from previous restart cycles as in {cite}`bjm2005`.

612:   Options Database Keys:
613: +   -ksp_gmres_restart <restart>                                                - total approximation space size (Krylov directions + error approximations)
614: .   -ksp_gmres_haptol <tol>                                                     - sets the tolerance for "happy ending" (exact convergence)
615: .   -ksp_gmres_preallocate                                                      - preallocate all the Krylov search directions initially
616:                                                                                 (otherwise groups of vectors are allocated as needed)
617: .   -ksp_gmres_classicalgramschmidt                                             - use classical (unmodified) Gram-Schmidt to orthogonalize against
618:                                                                                 the Krylov space (fast) (the default)
619: .   -ksp_gmres_modifiedgramschmidt                                              - use modified Gram-Schmidt in the orthogonalization (more stable, but slower)
620: .   -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always> - determine if iterative refinement is used to increase the
621:                                                                                 stability of the classical Gram-Schmidt  orthogonalization.
622: .   -ksp_gmres_krylov_monitor                                                   - plot the Krylov space generated
623: .   -ksp_lgmres_augment <k>                                                     - number of error approximations to augment the Krylov space with
624: -   -ksp_lgmres_constant                                                        - use a constant approximate space size
625:                                                                                 (only affects restart cycles < num. error approx.(k), i.e. the first k restarts)

627:   Level: beginner

629:   Notes:
630:   Supports both left and right preconditioning, but not symmetric.

632:   Augmenting with 1,2, or 3 approximations is generally optimal.

634:   This method is an accelerator for `KSPGMRES` - it is not useful for problems that stall. When gmres(m) stalls then lgmres with a size m
635:   approximation space will also generally stall.

637:   If gmres(m) converges in a small number of restart cycles, then lgmres also tends not to be very helpful.

639:   Developer Notes:
640:   To run LGMRES(m, k) as described in {cite}`bjm2005`, use\:
641: .vb
642:    -ksp_gmres_restart <m+k>
643:   -ksp_lgmres_augment <k>
644: .ve

646:   This object is subclassed off of `KSPGMRES`, see the source code in src/ksp/ksp/impls/gmres for comments on the structure of the code

648:   Contributed by:
649:   Allison Baker

651: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPFGMRES`, `KSPGMRES`,
652:           `KSPGMRESSetRestart()`, `KSPGMRESSetHapTol()`, `KSPGMRESSetPreAllocateVectors()`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESGetOrthogonalization()`,
653:           `KSPGMRESClassicalGramSchmidtOrthogonalization()`, `KSPGMRESModifiedGramSchmidtOrthogonalization()`,
654:           `KSPGMRESCGSRefinementType`, `KSPGMRESSetCGSRefinementType()`, `KSPGMRESGetCGSRefinementType()`, `KSPGMRESMonitorKrylov()`, `KSPLGMRESSetAugDim()`,
655:           `KSPGMRESSetConstant()`, `KSPLGMRESSetConstant()`, `KSPLGMRESSetAugDim()`
656: M*/

658: PETSC_EXTERN PetscErrorCode KSPCreate_LGMRES(KSP ksp)
659: {
660:   KSP_LGMRES *lgmres;

662:   PetscFunctionBegin;
663:   PetscCall(PetscNew(&lgmres));

665:   ksp->data               = (void *)lgmres;
666:   ksp->ops->buildsolution = KSPBuildSolution_LGMRES;

668:   ksp->ops->setup                        = KSPSetUp_LGMRES;
669:   ksp->ops->solve                        = KSPSolve_LGMRES;
670:   ksp->ops->destroy                      = KSPDestroy_LGMRES;
671:   ksp->ops->view                         = KSPView_LGMRES;
672:   ksp->ops->setfromoptions               = KSPSetFromOptions_LGMRES;
673:   ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_GMRES;
674:   ksp->ops->computeeigenvalues           = KSPComputeEigenvalues_GMRES;

676:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
677:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2));
678:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));

680:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetPreAllocateVectors_C", KSPGMRESSetPreAllocateVectors_GMRES));
681:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetOrthogonalization_C", KSPGMRESSetOrthogonalization_GMRES));
682:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetOrthogonalization_C", KSPGMRESGetOrthogonalization_GMRES));
683:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetRestart_C", KSPGMRESSetRestart_GMRES));
684:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetRestart_C", KSPGMRESGetRestart_GMRES));
685:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetHapTol_C", KSPGMRESSetHapTol_GMRES));
686:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetCGSRefinementType_C", KSPGMRESSetCGSRefinementType_GMRES));
687:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetCGSRefinementType_C", KSPGMRESGetCGSRefinementType_GMRES));

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

693:   /* defaults */
694:   lgmres->haptol         = 1.0e-30;
695:   lgmres->q_preallocate  = 0;
696:   lgmres->delta_allocate = LGMRES_DELTA_DIRECTIONS;
697:   lgmres->orthog         = KSPGMRESClassicalGramSchmidtOrthogonalization;
698:   lgmres->nrs            = NULL;
699:   lgmres->sol_temp       = NULL;
700:   lgmres->max_k          = LGMRES_DEFAULT_MAXK;
701:   lgmres->Rsvd           = NULL;
702:   lgmres->cgstype        = KSP_GMRES_CGS_REFINE_NEVER;
703:   lgmres->orthogwork     = NULL;

705:   /* LGMRES_MOD - new defaults */
706:   lgmres->aug_dim         = LGMRES_DEFAULT_AUGDIM;
707:   lgmres->aug_ct          = 0; /* start with no aug vectors */
708:   lgmres->approx_constant = PETSC_FALSE;
709:   lgmres->matvecs         = 0;
710:   PetscFunctionReturn(PETSC_SUCCESS);
711: }