Actual source code: fgmres.c


  2: /*
  3:     This file implements FGMRES (a Generalized Minimal Residual) method.
  4:     Reference:  Saad, 1993.

  6:     Preconditioning:  If the preconditioner is constant then this fgmres
  7:     code is equivalent to RIGHT-PRECONDITIONED GMRES.
  8:     FGMRES is a modification of gmres that allows the preconditioner to change
  9:     at each iteration.

 11:     Restarts:  Restarts are basically solves with x0 not equal to zero.
 12: */

 14: #include <../src/ksp/ksp/impls/gmres/fgmres/fgmresimpl.h>
 15: #define FGMRES_DELTA_DIRECTIONS 10
 16: #define FGMRES_DEFAULT_MAXK     30
 17: static PetscErrorCode KSPFGMRESGetNewVectors(KSP, PetscInt);
 18: static PetscErrorCode KSPFGMRESUpdateHessenberg(KSP, PetscInt, PetscBool, PetscReal *);
 19: static PetscErrorCode KSPFGMRESBuildSoln(PetscScalar *, Vec, Vec, KSP, PetscInt);

 21: /*

 23:     KSPSetUp_FGMRES - Sets up the workspace needed by fgmres.

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

 28: */
 29: PetscErrorCode KSPSetUp_FGMRES(KSP ksp)
 30: {
 31:   PetscInt    max_k, k;
 32:   KSP_FGMRES *fgmres = (KSP_FGMRES *)ksp->data;

 34:   PetscFunctionBegin;
 35:   max_k = fgmres->max_k;

 37:   PetscCall(KSPSetUp_GMRES(ksp));

 39:   PetscCall(PetscMalloc1(max_k + 2, &fgmres->prevecs));
 40:   PetscCall(PetscMalloc1(max_k + 2, &fgmres->prevecs_user_work));

 42:   /* fgmres->vv_allocated includes extra work vectors, which are not used in the additional
 43:      block of vectors used to store the preconditioned directions, hence  the -VEC_OFFSET
 44:      term for this first allocation of vectors holding preconditioned directions */
 45:   PetscCall(KSPCreateVecs(ksp, fgmres->vv_allocated - VEC_OFFSET, &fgmres->prevecs_user_work[0], 0, NULL));
 46:   for (k = 0; k < fgmres->vv_allocated - VEC_OFFSET; k++) fgmres->prevecs[k] = fgmres->prevecs_user_work[0][k];
 47:   PetscFunctionReturn(PETSC_SUCCESS);
 48: }

 50: /*
 51:     KSPFGMRESResidual - This routine computes the initial residual (NOT PRECONDITIONED)
 52: */
 53: static PetscErrorCode KSPFGMRESResidual(KSP ksp)
 54: {
 55:   KSP_FGMRES *fgmres = (KSP_FGMRES *)(ksp->data);
 56:   Mat         Amat, Pmat;

 58:   PetscFunctionBegin;
 59:   PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));

 61:   /* put A*x into VEC_TEMP */
 62:   PetscCall(KSP_MatMult(ksp, Amat, ksp->vec_sol, VEC_TEMP));
 63:   /* now put residual (-A*x + f) into vec_vv(0) */
 64:   PetscCall(VecWAXPY(VEC_VV(0), -1.0, VEC_TEMP, ksp->vec_rhs));
 65:   PetscFunctionReturn(PETSC_SUCCESS);
 66: }

 68: /*

 70:     KSPFGMRESCycle - Run fgmres, possibly with restart.  Return residual
 71:                   history if requested.

 73:     input parameters:
 74: .        fgmres  - structure containing parameters and work areas

 76:     output parameters:
 77: .        itcount - number of iterations used.  If null, ignored.
 78: .        converged - 0 if not converged

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

 84:  */
 85: PetscErrorCode KSPFGMRESCycle(PetscInt *itcount, KSP ksp)
 86: {
 87:   KSP_FGMRES *fgmres = (KSP_FGMRES *)(ksp->data);
 88:   PetscReal   res_norm;
 89:   PetscReal   hapbnd, tt;
 90:   PetscBool   hapend = PETSC_FALSE;  /* indicates happy breakdown ending */
 91:   PetscInt    loc_it;                /* local count of # of dir. in Krylov space */
 92:   PetscInt    max_k = fgmres->max_k; /* max # of directions Krylov space */
 93:   Mat         Amat, Pmat;

 95:   PetscFunctionBegin;
 96:   /* Number of pseudo iterations since last restart is the number
 97:      of prestart directions */
 98:   loc_it = 0;

100:   /* note: (fgmres->it) is always set one less than (loc_it) It is used in
101:      KSPBUILDSolution_FGMRES, where it is passed to KSPFGMRESBuildSoln.
102:      Note that when KSPFGMRESBuildSoln is called from this function,
103:      (loc_it -1) is passed, so the two are equivalent */
104:   fgmres->it = (loc_it - 1);

106:   /* initial residual is in VEC_VV(0)  - compute its norm*/
107:   PetscCall(VecNorm(VEC_VV(0), NORM_2, &res_norm));
108:   KSPCheckNorm(ksp, res_norm);

110:   /* first entry in right-hand-side of hessenberg system is just
111:      the initial residual norm */
112:   *RS(0) = res_norm;

114:   ksp->rnorm = res_norm;
115:   PetscCall(KSPLogResidualHistory(ksp, res_norm));
116:   PetscCall(KSPMonitor(ksp, ksp->its, res_norm));

118:   /* check for the convergence - maybe the current guess is good enough */
119:   PetscCall((*ksp->converged)(ksp, ksp->its, res_norm, &ksp->reason, ksp->cnvP));
120:   if (ksp->reason) {
121:     if (itcount) *itcount = 0;
122:     PetscFunctionReturn(PETSC_SUCCESS);
123:   }

125:   /* scale VEC_VV (the initial residual) */
126:   PetscCall(VecScale(VEC_VV(0), 1.0 / res_norm));

128:   /* MAIN ITERATION LOOP BEGINNING*/
129:   /* keep iterating until we have converged OR generated the max number
130:      of directions OR reached the max number of iterations for the method */
131:   while (!ksp->reason && loc_it < max_k && ksp->its < ksp->max_it) {
132:     if (loc_it) {
133:       PetscCall(KSPLogResidualHistory(ksp, res_norm));
134:       PetscCall(KSPMonitor(ksp, ksp->its, res_norm));
135:     }
136:     fgmres->it = (loc_it - 1);

138:     /* see if more space is needed for work vectors */
139:     if (fgmres->vv_allocated <= loc_it + VEC_OFFSET + 1) {
140:       PetscCall(KSPFGMRESGetNewVectors(ksp, loc_it + 1));
141:       /* (loc_it+1) is passed in as number of the first vector that should
142:          be allocated */
143:     }

145:     /* CHANGE THE PRECONDITIONER? */
146:     /* ModifyPC is the callback function that can be used to
147:        change the PC or its attributes before its applied */
148:     PetscCall((*fgmres->modifypc)(ksp, ksp->its, loc_it, res_norm, fgmres->modifyctx));

150:     /* apply PRECONDITIONER to direction vector and store with
151:        preconditioned vectors in prevec */
152:     PetscCall(KSP_PCApply(ksp, VEC_VV(loc_it), PREVEC(loc_it)));

154:     PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
155:     /* Multiply preconditioned vector by operator - put in VEC_VV(loc_it+1) */
156:     PetscCall(KSP_MatMult(ksp, Amat, PREVEC(loc_it), VEC_VV(1 + loc_it)));

158:     /* update hessenberg matrix and do Gram-Schmidt - new direction is in
159:        VEC_VV(1+loc_it)*/
160:     PetscCall((*fgmres->orthog)(ksp, loc_it));

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

166:     *HH(loc_it + 1, loc_it)  = tt;
167:     *HES(loc_it + 1, loc_it) = tt;

169:     /* Happy Breakdown Check */
170:     hapbnd = PetscAbsScalar((tt) / *RS(loc_it));
171:     /* RS(loc_it) contains the res_norm from the last iteration  */
172:     hapbnd = PetscMin(fgmres->haptol, hapbnd);
173:     if (tt > hapbnd) {
174:       /* scale new direction by its norm */
175:       PetscCall(VecScale(VEC_VV(loc_it + 1), 1.0 / tt));
176:     } else {
177:       /* This happens when the solution is exactly reached. */
178:       /* So there is no new direction... */
179:       PetscCall(VecSet(VEC_TEMP, 0.0)); /* set VEC_TEMP to 0 */
180:       hapend = PETSC_TRUE;
181:     }
182:     /* note that for FGMRES we could get HES(loc_it+1, loc_it)  = 0 and the
183:        current solution would not be exact if HES was singular.  Note that
184:        HH non-singular implies that HES is no singular, and HES is guaranteed
185:        to be nonsingular when PREVECS are linearly independent and A is
186:        nonsingular (in GMRES, the nonsingularity of A implies the nonsingularity
187:        of HES). So we should really add a check to verify that HES is nonsingular.*/

189:     /* Now apply rotations to new col of hessenberg (and right side of system),
190:        calculate new rotation, and get new residual norm at the same time*/
191:     PetscCall(KSPFGMRESUpdateHessenberg(ksp, loc_it, hapend, &res_norm));
192:     if (ksp->reason) break;

194:     loc_it++;
195:     fgmres->it = (loc_it - 1); /* Add this here in case it has converged */

197:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
198:     ksp->its++;
199:     ksp->rnorm = res_norm;
200:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));

202:     PetscCall((*ksp->converged)(ksp, ksp->its, res_norm, &ksp->reason, ksp->cnvP));

204:     /* Catch error in happy breakdown and signal convergence and break from loop */
205:     if (hapend) {
206:       if (!ksp->reason) {
207:         PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "You reached the happy break down, but convergence was not indicated. Residual norm = %g", (double)res_norm);
208:         ksp->reason = KSP_DIVERGED_BREAKDOWN;
209:         break;
210:       }
211:     }
212:   }
213:   /* END OF ITERATION LOOP */
214:   PetscCall(KSPLogResidualHistory(ksp, res_norm));

216:   /*
217:      Monitor if we know that we will not return for a restart */
218:   if (loc_it && (ksp->reason || ksp->its >= ksp->max_it)) PetscCall(KSPMonitor(ksp, ksp->its, res_norm));

220:   if (itcount) *itcount = loc_it;

222:   /*
223:     Down here we have to solve for the "best" coefficients of the Krylov
224:     columns, add the solution values together, and possibly unwind the
225:     preconditioning from the solution
226:    */

228:   /* Form the solution (or the solution so far) */
229:   /* Note: must pass in (loc_it-1) for iteration count so that KSPFGMRESBuildSoln
230:      properly navigates */

232:   PetscCall(KSPFGMRESBuildSoln(RS(0), ksp->vec_sol, ksp->vec_sol, ksp, loc_it - 1));
233:   PetscFunctionReturn(PETSC_SUCCESS);
234: }

236: /*
237:     KSPSolve_FGMRES - This routine applies the FGMRES method.

239:    Input Parameter:
240: .     ksp - the Krylov space object that was set to use fgmres

242:    Output Parameter:
243: .     outits - number of iterations used

245: */

247: PetscErrorCode KSPSolve_FGMRES(KSP ksp)
248: {
249:   PetscInt    cycle_its = 0; /* iterations done in a call to KSPFGMRESCycle */
250:   KSP_FGMRES *fgmres    = (KSP_FGMRES *)ksp->data;
251:   PetscBool   diagonalscale;

253:   PetscFunctionBegin;
254:   PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
255:   PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);

257:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
258:   ksp->its = 0;
259:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));

261:   /* Compute the initial (NOT preconditioned) residual */
262:   if (!ksp->guess_zero) {
263:     PetscCall(KSPFGMRESResidual(ksp));
264:   } else { /* guess is 0 so residual is F (which is in ksp->vec_rhs) */
265:     PetscCall(VecCopy(ksp->vec_rhs, VEC_VV(0)));
266:   }
267:   /* This may be true only on a subset of MPI ranks; setting it here so it will be detected by the first norm computation in the Krylov method */
268:   if (ksp->reason == KSP_DIVERGED_PC_FAILED) PetscCall(VecSetInf(VEC_VV(0)));

270:   /* now the residual is in VEC_VV(0) - which is what
271:      KSPFGMRESCycle expects... */

273:   PetscCall(KSPFGMRESCycle(&cycle_its, ksp));
274:   while (!ksp->reason) {
275:     PetscCall(KSPFGMRESResidual(ksp));
276:     if (ksp->its >= ksp->max_it) break;
277:     PetscCall(KSPFGMRESCycle(&cycle_its, ksp));
278:   }
279:   /* mark lack of convergence */
280:   if (ksp->its >= ksp->max_it && !ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
281:   PetscFunctionReturn(PETSC_SUCCESS);
282: }

284: extern PetscErrorCode KSPReset_FGMRES(KSP);
285: /*

287:    KSPDestroy_FGMRES - Frees all memory space used by the Krylov method.

289: */
290: PetscErrorCode KSPDestroy_FGMRES(KSP ksp)
291: {
292:   PetscFunctionBegin;
293:   PetscCall(KSPReset_FGMRES(ksp));
294:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPFGMRESSetModifyPC_C", NULL));
295:   PetscCall(KSPDestroy_GMRES(ksp));
296:   PetscFunctionReturn(PETSC_SUCCESS);
297: }

299: /*
300:     KSPFGMRESBuildSoln - create the solution from the starting vector and the
301:                       current iterates.

303:     Input parameters:
304:         nrs - work area of size it + 1.
305:         vguess  - index of initial guess
306:         vdest - index of result.  Note that vguess may == vdest (replace
307:                 guess with the solution).
308:         it - HH upper triangular part is a block of size (it+1) x (it+1)

310:      This is an internal routine that knows about the FGMRES internals.
311:  */
312: static PetscErrorCode KSPFGMRESBuildSoln(PetscScalar *nrs, Vec vguess, Vec vdest, KSP ksp, PetscInt it)
313: {
314:   PetscScalar tt;
315:   PetscInt    ii, k, j;
316:   KSP_FGMRES *fgmres = (KSP_FGMRES *)(ksp->data);

318:   PetscFunctionBegin;
319:   /* Solve for solution vector that minimizes the residual */

321:   /* If it is < 0, no fgmres steps have been performed */
322:   if (it < 0) {
323:     PetscCall(VecCopy(vguess, vdest)); /* VecCopy() is smart, exists immediately if vguess == vdest */
324:     PetscFunctionReturn(PETSC_SUCCESS);
325:   }

327:   /* so fgmres steps HAVE been performed */

329:   /* solve the upper triangular system - RS is the right side and HH is
330:      the upper triangular matrix  - put soln in nrs */
331:   if (*HH(it, it) != 0.0) {
332:     nrs[it] = *RS(it) / *HH(it, it);
333:   } else {
334:     nrs[it] = 0.0;
335:   }
336:   for (ii = 1; ii <= it; ii++) {
337:     k  = it - ii;
338:     tt = *RS(k);
339:     for (j = k + 1; j <= it; j++) tt = tt - *HH(k, j) * nrs[j];
340:     nrs[k] = tt / *HH(k, k);
341:   }

343:   /* Accumulate the correction to the soln of the preconditioned prob. in
344:      VEC_TEMP - note that we use the preconditioned vectors  */
345:   PetscCall(VecSet(VEC_TEMP, 0.0)); /* set VEC_TEMP components to 0 */
346:   PetscCall(VecMAXPY(VEC_TEMP, it + 1, nrs, &PREVEC(0)));

348:   /* put updated solution into vdest.*/
349:   if (vdest != vguess) {
350:     PetscCall(VecCopy(VEC_TEMP, vdest));
351:     PetscCall(VecAXPY(vdest, 1.0, vguess));
352:   } else { /* replace guess with solution */
353:     PetscCall(VecAXPY(vdest, 1.0, VEC_TEMP));
354:   }
355:   PetscFunctionReturn(PETSC_SUCCESS);
356: }

358: /*

360:     KSPFGMRESUpdateHessenberg - Do the scalar work for the orthogonalization.
361:                             Return new residual.

363:     input parameters:

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

370:     output parameters:
371: .        res - the new residual

373:  */
374: static PetscErrorCode KSPFGMRESUpdateHessenberg(KSP ksp, PetscInt it, PetscBool hapend, PetscReal *res)
375: {
376:   PetscScalar *hh, *cc, *ss, tt;
377:   PetscInt     j;
378:   KSP_FGMRES  *fgmres = (KSP_FGMRES *)(ksp->data);

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

386:   /* Apply all the previously computed plane rotations to the new column
387:      of the Hessenberg matrix */
388:   /* Note: this uses the rotation [conj(c)  s ; -s   c], c= cos(theta), s= sin(theta),
389:      and some refs have [c   s ; -conj(s)  c] (don't be confused!) */

391:   for (j = 1; j <= it; j++) {
392:     tt  = *hh;
393:     *hh = PetscConj(*cc) * tt + *ss * *(hh + 1);
394:     hh++;
395:     *hh = *cc++ * *hh - (*ss++ * tt);
396:     /* hh, cc, and ss have all been incremented one by end of loop */
397:   }

399:   /*
400:     compute the new plane rotation, and apply it to:
401:      1) the right-hand-side of the Hessenberg system (RS)
402:         note: it affects RS(it) and RS(it+1)
403:      2) the new column of the Hessenberg matrix
404:         note: it affects HH(it,it) which is currently pointed to
405:         by hh and HH(it+1, it) (*(hh+1))
406:     thus obtaining the updated value of the residual...
407:   */

409:   /* compute new plane rotation */

411:   if (!hapend) {
412:     tt = PetscSqrtScalar(PetscConj(*hh) * *hh + PetscConj(*(hh + 1)) * *(hh + 1));
413:     if (tt == 0.0) {
414:       ksp->reason = KSP_DIVERGED_NULL;
415:       PetscFunctionReturn(PETSC_SUCCESS);
416:     }

418:     *cc = *hh / tt;       /* new cosine value */
419:     *ss = *(hh + 1) / tt; /* new sine value */

421:     /* apply to 1) and 2) */
422:     *RS(it + 1) = -(*ss * *RS(it));
423:     *RS(it)     = PetscConj(*cc) * *RS(it);
424:     *hh         = PetscConj(*cc) * *hh + *ss * *(hh + 1);

426:     /* residual is the last element (it+1) of right-hand side! */
427:     *res = PetscAbsScalar(*RS(it + 1));

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

436:     *res = 0.0;
437:   }
438:   PetscFunctionReturn(PETSC_SUCCESS);
439: }

441: /*

443:    KSPFGMRESGetNewVectors - This routine allocates more work vectors, starting from
444:                          VEC_VV(it), and more preconditioned work vectors, starting
445:                          from PREVEC(i).

447: */
448: static PetscErrorCode KSPFGMRESGetNewVectors(KSP ksp, PetscInt it)
449: {
450:   KSP_FGMRES *fgmres = (KSP_FGMRES *)ksp->data;
451:   PetscInt    nwork  = fgmres->nwork_alloc; /* number of work vector chunks allocated */
452:   PetscInt    nalloc;                       /* number to allocate */
453:   PetscInt    k;

455:   PetscFunctionBegin;
456:   nalloc = fgmres->delta_allocate; /* number of vectors to allocate
457:                                       in a single chunk */

459:   /* Adjust the number to allocate to make sure that we don't exceed the
460:      number of available slots (fgmres->vecs_allocated)*/
461:   if (it + VEC_OFFSET + nalloc >= fgmres->vecs_allocated) nalloc = fgmres->vecs_allocated - it - VEC_OFFSET;
462:   if (!nalloc) PetscFunctionReturn(PETSC_SUCCESS);

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

466:   /* work vectors */
467:   PetscCall(KSPCreateVecs(ksp, nalloc, &fgmres->user_work[nwork], 0, NULL));
468:   for (k = 0; k < nalloc; k++) fgmres->vecs[it + VEC_OFFSET + k] = fgmres->user_work[nwork][k];
469:   /* specify size of chunk allocated */
470:   fgmres->mwork_alloc[nwork] = nalloc;

472:   /* preconditioned vectors */
473:   PetscCall(KSPCreateVecs(ksp, nalloc, &fgmres->prevecs_user_work[nwork], 0, NULL));
474:   for (k = 0; k < nalloc; k++) fgmres->prevecs[it + k] = fgmres->prevecs_user_work[nwork][k];

476:   /* increment the number of work vector chunks */
477:   fgmres->nwork_alloc++;
478:   PetscFunctionReturn(PETSC_SUCCESS);
479: }

481: /*

483:    KSPBuildSolution_FGMRES

485:      Input Parameter:
486: .     ksp - the Krylov space object
487: .     ptr-

489:    Output Parameter:
490: .     result - the solution

492:    Note: this calls KSPFGMRESBuildSoln - the same function that KSPFGMRESCycle
493:    calls directly.

495: */
496: PetscErrorCode KSPBuildSolution_FGMRES(KSP ksp, Vec ptr, Vec *result)
497: {
498:   KSP_FGMRES *fgmres = (KSP_FGMRES *)ksp->data;

500:   PetscFunctionBegin;
501:   if (!ptr) {
502:     if (!fgmres->sol_temp) PetscCall(VecDuplicate(ksp->vec_sol, &fgmres->sol_temp));
503:     ptr = fgmres->sol_temp;
504:   }
505:   if (!fgmres->nrs) {
506:     /* allocate the work area */
507:     PetscCall(PetscMalloc1(fgmres->max_k, &fgmres->nrs));
508:   }

510:   PetscCall(KSPFGMRESBuildSoln(fgmres->nrs, ksp->vec_sol, ptr, ksp, fgmres->it));
511:   if (result) *result = ptr;
512:   PetscFunctionReturn(PETSC_SUCCESS);
513: }

515: PetscErrorCode KSPSetFromOptions_FGMRES(KSP ksp, PetscOptionItems *PetscOptionsObject)
516: {
517:   PetscBool flg;

519:   PetscFunctionBegin;
520:   PetscCall(KSPSetFromOptions_GMRES(ksp, PetscOptionsObject));
521:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP flexible GMRES Options");
522:   PetscCall(PetscOptionsBoolGroupBegin("-ksp_fgmres_modifypcnochange", "do not vary the preconditioner", "KSPFGMRESSetModifyPC", &flg));
523:   if (flg) PetscCall(KSPFGMRESSetModifyPC(ksp, KSPFGMRESModifyPCNoChange, NULL, NULL));
524:   PetscCall(PetscOptionsBoolGroupEnd("-ksp_fgmres_modifypcksp", "vary the KSP based preconditioner", "KSPFGMRESSetModifyPC", &flg));
525:   if (flg) PetscCall(KSPFGMRESSetModifyPC(ksp, KSPFGMRESModifyPCKSP, NULL, NULL));
526:   PetscOptionsHeadEnd();
527:   PetscFunctionReturn(PETSC_SUCCESS);
528: }

530: typedef PetscErrorCode (*FCN1)(KSP, PetscInt, PetscInt, PetscReal, void *); /* force argument to next function to not be extern C*/
531: typedef PetscErrorCode (*FCN2)(void *);

533: static PetscErrorCode KSPFGMRESSetModifyPC_FGMRES(KSP ksp, FCN1 fcn, void *ctx, FCN2 d)
534: {
535:   PetscFunctionBegin;
537:   ((KSP_FGMRES *)ksp->data)->modifypc      = fcn;
538:   ((KSP_FGMRES *)ksp->data)->modifydestroy = d;
539:   ((KSP_FGMRES *)ksp->data)->modifyctx     = ctx;
540:   PetscFunctionReturn(PETSC_SUCCESS);
541: }

543: PetscErrorCode KSPReset_FGMRES(KSP ksp)
544: {
545:   KSP_FGMRES *fgmres = (KSP_FGMRES *)ksp->data;
546:   PetscInt    i;

548:   PetscFunctionBegin;
549:   PetscCall(PetscFree(fgmres->prevecs));
550:   if (fgmres->nwork_alloc > 0) {
551:     i = 0;
552:     /* In the first allocation we allocated VEC_OFFSET fewer vectors in prevecs */
553:     PetscCall(VecDestroyVecs(fgmres->mwork_alloc[i] - VEC_OFFSET, &fgmres->prevecs_user_work[i]));
554:     for (i = 1; i < fgmres->nwork_alloc; i++) PetscCall(VecDestroyVecs(fgmres->mwork_alloc[i], &fgmres->prevecs_user_work[i]));
555:   }
556:   PetscCall(PetscFree(fgmres->prevecs_user_work));
557:   if (fgmres->modifydestroy) PetscCall((*fgmres->modifydestroy)(fgmres->modifyctx));
558:   PetscCall(KSPReset_GMRES(ksp));
559:   PetscFunctionReturn(PETSC_SUCCESS);
560: }

562: PetscErrorCode KSPGMRESSetRestart_FGMRES(KSP ksp, PetscInt max_k)
563: {
564:   KSP_FGMRES *gmres = (KSP_FGMRES *)ksp->data;

566:   PetscFunctionBegin;
567:   PetscCheck(max_k >= 1, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Restart must be positive");
568:   if (!ksp->setupstage) {
569:     gmres->max_k = max_k;
570:   } else if (gmres->max_k != max_k) {
571:     gmres->max_k    = max_k;
572:     ksp->setupstage = KSP_SETUP_NEW;
573:     /* free the data structures, then create them again */
574:     PetscCall(KSPReset_FGMRES(ksp));
575:   }
576:   PetscFunctionReturn(PETSC_SUCCESS);
577: }

579: PetscErrorCode KSPGMRESGetRestart_FGMRES(KSP ksp, PetscInt *max_k)
580: {
581:   KSP_FGMRES *gmres = (KSP_FGMRES *)ksp->data;

583:   PetscFunctionBegin;
584:   *max_k = gmres->max_k;
585:   PetscFunctionReturn(PETSC_SUCCESS);
586: }

588: /*MC
589:      KSPFGMRES - Implements the Flexible Generalized Minimal Residual method. [](sec_flexibleksp)

591:    Options Database Keys:
592: +   -ksp_gmres_restart <restart> - the number of Krylov directions to orthogonalize against
593: .   -ksp_gmres_haptol <tol> - sets the tolerance for "happy ending" (exact convergence)
594: .   -ksp_gmres_preallocate - preallocate all the Krylov search directions initially (otherwise groups of
595:                              vectors are allocated as needed)
596: .   -ksp_gmres_classicalgramschmidt - use classical (unmodified) Gram-Schmidt to orthogonalize against the Krylov space (fast) (the default)
597: .   -ksp_gmres_modifiedgramschmidt - use modified Gram-Schmidt in the orthogonalization (more stable, but slower)
598: .   -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always> - determine if iterative refinement is used to increase the
599:                                    stability of the classical Gram-Schmidt  orthogonalization.
600: .   -ksp_gmres_krylov_monitor - plot the Krylov space generated
601: .   -ksp_fgmres_modifypcnochange - do not change the preconditioner between iterations
602: -   -ksp_fgmres_modifypcksp - modify the preconditioner using `KSPFGMRESModifyPCKSP()`

604:    Level: beginner

606:     Notes:
607:     See `KSPFGMRESSetModifyPC()` for how to vary the preconditioner between iterations

609:     Only right preconditioning is supported.

611:     The following options -ksp_type fgmres -pc_type ksp -ksp_ksp_type bcgs -ksp_view -ksp_pc_type jacobi make the preconditioner (or inner solver)
612:     be bi-CG-stab with a preconditioner of Jacobi.

614:     Developer Note:
615:     This object is subclassed off of `KSPGMRES`

617:     Contributed by:
618:     Allison Baker

620: .seealso: [](ch_ksp), [](sec_flexibleksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPGMRES`, `KSPLGMRES`,
621:           `KSPGMRESSetRestart()`, `KSPGMRESSetHapTol()`, `KSPGMRESSetPreAllocateVectors()`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESGetOrthogonalization()`,
622:           `KSPGMRESClassicalGramSchmidtOrthogonalization()`, `KSPGMRESModifiedGramSchmidtOrthogonalization()`,
623:           `KSPGMRESCGSRefinementType`, `KSPGMRESSetCGSRefinementType()`, `KSPGMRESGetCGSRefinementType()`, `KSPGMRESMonitorKrylov()`, `KSPFGMRESSetModifyPC()`,
624:           `KSPFGMRESModifyPCKSP()`
625: M*/

627: PETSC_EXTERN PetscErrorCode KSPCreate_FGMRES(KSP ksp)
628: {
629:   KSP_FGMRES *fgmres;

631:   PetscFunctionBegin;
632:   PetscCall(PetscNew(&fgmres));

634:   ksp->data                              = (void *)fgmres;
635:   ksp->ops->buildsolution                = KSPBuildSolution_FGMRES;
636:   ksp->ops->setup                        = KSPSetUp_FGMRES;
637:   ksp->ops->solve                        = KSPSolve_FGMRES;
638:   ksp->ops->reset                        = KSPReset_FGMRES;
639:   ksp->ops->destroy                      = KSPDestroy_FGMRES;
640:   ksp->ops->view                         = KSPView_GMRES;
641:   ksp->ops->setfromoptions               = KSPSetFromOptions_FGMRES;
642:   ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_GMRES;
643:   ksp->ops->computeeigenvalues           = KSPComputeEigenvalues_GMRES;

645:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 3));
646:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));

648:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetPreAllocateVectors_C", KSPGMRESSetPreAllocateVectors_GMRES));
649:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetOrthogonalization_C", KSPGMRESSetOrthogonalization_GMRES));
650:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetOrthogonalization_C", KSPGMRESGetOrthogonalization_GMRES));
651:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetRestart_C", KSPGMRESSetRestart_FGMRES));
652:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetRestart_C", KSPGMRESGetRestart_FGMRES));
653:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPFGMRESSetModifyPC_C", KSPFGMRESSetModifyPC_FGMRES));
654:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetCGSRefinementType_C", KSPGMRESSetCGSRefinementType_GMRES));
655:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetCGSRefinementType_C", KSPGMRESGetCGSRefinementType_GMRES));

657:   fgmres->haptol         = 1.0e-30;
658:   fgmres->q_preallocate  = 0;
659:   fgmres->delta_allocate = FGMRES_DELTA_DIRECTIONS;
660:   fgmres->orthog         = KSPGMRESClassicalGramSchmidtOrthogonalization;
661:   fgmres->nrs            = NULL;
662:   fgmres->sol_temp       = NULL;
663:   fgmres->max_k          = FGMRES_DEFAULT_MAXK;
664:   fgmres->Rsvd           = NULL;
665:   fgmres->orthogwork     = NULL;
666:   fgmres->modifypc       = KSPFGMRESModifyPCNoChange;
667:   fgmres->modifyctx      = NULL;
668:   fgmres->modifydestroy  = NULL;
669:   fgmres->cgstype        = KSP_GMRES_CGS_REFINE_NEVER;
670:   PetscFunctionReturn(PETSC_SUCCESS);
671: }