Actual source code: pgmres.c


  2: /*
  3:     This file implements PGMRES (a Pipelined Generalized Minimal Residual method)
  4: */

  6: #include <../src/ksp/ksp/impls/gmres/pgmres/pgmresimpl.h>
  7: #define PGMRES_DELTA_DIRECTIONS 10
  8: #define PGMRES_DEFAULT_MAXK     30

 10: static PetscErrorCode KSPPGMRESUpdateHessenberg(KSP, PetscInt, PetscBool *, PetscReal *);
 11: static PetscErrorCode KSPPGMRESBuildSoln(PetscScalar *, Vec, Vec, KSP, PetscInt);

 13: /*

 15:     KSPSetUp_PGMRES - Sets up the workspace needed by pgmres.

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

 20: */
 21: static PetscErrorCode KSPSetUp_PGMRES(KSP ksp)
 22: {
 23:   PetscFunctionBegin;
 24:   PetscCall(KSPSetUp_GMRES(ksp));
 25:   PetscFunctionReturn(PETSC_SUCCESS);
 26: }

 28: /*

 30:     KSPPGMRESCycle - Run pgmres, possibly with restart.  Return residual
 31:                   history if requested.

 33:     input parameters:
 34: .        pgmres  - structure containing parameters and work areas

 36:     output parameters:
 37: .        itcount - number of iterations used.  If null, ignored.
 38: .        converged - 0 if not converged

 40:     Notes:
 41:     On entry, the value in vector VEC_VV(0) should be
 42:     the initial residual.

 44:  */
 45: static PetscErrorCode KSPPGMRESCycle(PetscInt *itcount, KSP ksp)
 46: {
 47:   KSP_PGMRES *pgmres = (KSP_PGMRES *)(ksp->data);
 48:   PetscReal   res_norm, res, newnorm;
 49:   PetscInt    it     = 0, j, k;
 50:   PetscBool   hapend = PETSC_FALSE;

 52:   PetscFunctionBegin;
 53:   if (itcount) *itcount = 0;
 54:   PetscCall(VecNormalize(VEC_VV(0), &res_norm));
 55:   KSPCheckNorm(ksp, res_norm);
 56:   res    = res_norm;
 57:   *RS(0) = res_norm;

 59:   /* check for the convergence */
 60:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
 61:   if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = res;
 62:   else ksp->rnorm = 0;
 63:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
 64:   pgmres->it = it - 2;
 65:   PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
 66:   PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm));
 67:   if (!res) {
 68:     ksp->reason = KSP_CONVERGED_ATOL;
 69:     PetscCall(PetscInfo(ksp, "Converged due to zero residual norm on entry\n"));
 70:     PetscFunctionReturn(PETSC_SUCCESS);
 71:   }

 73:   PetscCall((*ksp->converged)(ksp, ksp->its, ksp->rnorm, &ksp->reason, ksp->cnvP));
 74:   for (; !ksp->reason; it++) {
 75:     Vec Zcur, Znext;
 76:     if (pgmres->vv_allocated <= it + VEC_OFFSET + 1) PetscCall(KSPGMRESGetNewVectors(ksp, it + 1));
 77:     /* VEC_VV(it-1) is orthogonal, it will be normalized once the VecNorm arrives. */
 78:     Zcur  = VEC_VV(it);     /* Zcur is not yet orthogonal, but the VecMDot to orthogonalize it has been started. */
 79:     Znext = VEC_VV(it + 1); /* This iteration will compute Znext, update with a deferred correction once we know how
 80:                                  * Zcur relates to the previous vectors, and start the reduction to orthogonalize it. */

 82:     if (it < pgmres->max_k + 1 && ksp->its + 1 < PetscMax(2, ksp->max_it)) { /* We don't know whether what we have computed is enough, so apply the matrix. */
 83:       PetscCall(KSP_PCApplyBAorAB(ksp, Zcur, Znext, VEC_TEMP_MATOP));
 84:     }

 86:     if (it > 1) { /* Complete the pending reduction */
 87:       PetscCall(VecNormEnd(VEC_VV(it - 1), NORM_2, &newnorm));
 88:       *HH(it - 1, it - 2) = newnorm;
 89:     }
 90:     if (it > 0) { /* Finish the reduction computing the latest column of H */
 91:       PetscCall(VecMDotEnd(Zcur, it, &(VEC_VV(0)), HH(0, it - 1)));
 92:     }

 94:     if (it > 1) {
 95:       /* normalize the base vector from two iterations ago, basis is complete up to here */
 96:       PetscCall(VecScale(VEC_VV(it - 1), 1. / *HH(it - 1, it - 2)));

 98:       PetscCall(KSPPGMRESUpdateHessenberg(ksp, it - 2, &hapend, &res));
 99:       pgmres->it = it - 2;
100:       ksp->its++;
101:       if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = res;
102:       else ksp->rnorm = 0;

104:       PetscCall((*ksp->converged)(ksp, ksp->its, ksp->rnorm, &ksp->reason, ksp->cnvP));
105:       if (it < pgmres->max_k + 1 || ksp->reason || ksp->its == ksp->max_it) { /* Monitor if we are done or still iterating, but not before a restart. */
106:         PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
107:         PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm));
108:       }
109:       if (ksp->reason) break;
110:       /* Catch error in happy breakdown and signal convergence and break from loop */
111:       if (hapend) {
112:         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);
113:         ksp->reason = KSP_DIVERGED_BREAKDOWN;
114:         break;
115:       }

117:       if (!(it < pgmres->max_k + 1 && ksp->its < ksp->max_it)) break;

119:       /* The it-2 column of H was not scaled when we computed Zcur, apply correction */
120:       PetscCall(VecScale(Zcur, 1. / *HH(it - 1, it - 2)));
121:       /* And Znext computed in this iteration was computed using the under-scaled Zcur */
122:       PetscCall(VecScale(Znext, 1. / *HH(it - 1, it - 2)));

124:       /* In the previous iteration, we projected an unnormalized Zcur against the Krylov basis, so we need to fix the column of H resulting from that projection. */
125:       for (k = 0; k < it; k++) *HH(k, it - 1) /= *HH(it - 1, it - 2);
126:       /* When Zcur was projected against the Krylov basis, VV(it-1) was still not normalized, so fix that too. This
127:        * column is complete except for HH(it,it-1) which we won't know until the next iteration. */
128:       *HH(it - 1, it - 1) /= *HH(it - 1, it - 2);
129:     }

131:     if (it > 0) {
132:       PetscScalar *work;
133:       if (!pgmres->orthogwork) PetscCall(PetscMalloc1(pgmres->max_k + 2, &pgmres->orthogwork));
134:       work = pgmres->orthogwork;
135:       /* Apply correction computed by the VecMDot in the last iteration to Znext. The original form is
136:        *
137:        *   Znext -= sum_{j=0}^{i-1} Z[j+1] * H[j,i-1]
138:        *
139:        * where
140:        *
141:        *   Z[j] = sum_{k=0}^j V[k] * H[k,j-1]
142:        *
143:        * substituting
144:        *
145:        *   Znext -= sum_{j=0}^{i-1} sum_{k=0}^{j+1} V[k] * H[k,j] * H[j,i-1]
146:        *
147:        * rearranging the iteration space from row-column to column-row
148:        *
149:        *   Znext -= sum_{k=0}^i sum_{j=k-1}^{i-1} V[k] * H[k,j] * H[j,i-1]
150:        *
151:        * Note that column it-1 of HH is correct. For all previous columns, we must look at HES because HH has already
152:        * been transformed to upper triangular form.
153:        */
154:       for (k = 0; k < it + 1; k++) {
155:         work[k] = 0;
156:         for (j = PetscMax(0, k - 1); j < it - 1; j++) work[k] -= *HES(k, j) * *HH(j, it - 1);
157:       }
158:       PetscCall(VecMAXPY(Znext, it + 1, work, &VEC_VV(0)));
159:       PetscCall(VecAXPY(Znext, -*HH(it - 1, it - 1), Zcur));

161:       /* Orthogonalize Zcur against existing basis vectors. */
162:       for (k = 0; k < it; k++) work[k] = -*HH(k, it - 1);
163:       PetscCall(VecMAXPY(Zcur, it, work, &VEC_VV(0)));
164:       /* Zcur is now orthogonal, and will be referred to as VEC_VV(it) again, though it is still not normalized. */
165:       /* Begin computing the norm of the new vector, will be normalized after the MatMult in the next iteration. */
166:       PetscCall(VecNormBegin(VEC_VV(it), NORM_2, &newnorm));
167:     }

169:     /* Compute column of H (to the diagonal, but not the subdiagonal) to be able to orthogonalize the newest vector. */
170:     PetscCall(VecMDotBegin(Znext, it + 1, &VEC_VV(0), HH(0, it)));

172:     /* Start an asynchronous split-mode reduction, the result of the MDot and Norm will be collected on the next iteration. */
173:     PetscCall(PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)Znext)));
174:   }

176:   if (itcount) *itcount = it - 1; /* Number of iterations actually completed. */

178:   /*
179:     Down here we have to solve for the "best" coefficients of the Krylov
180:     columns, add the solution values together, and possibly unwind the
181:     preconditioning from the solution
182:    */
183:   /* Form the solution (or the solution so far) */
184:   PetscCall(KSPPGMRESBuildSoln(RS(0), ksp->vec_sol, ksp->vec_sol, ksp, it - 2));
185:   PetscFunctionReturn(PETSC_SUCCESS);
186: }

188: /*
189:     KSPSolve_PGMRES - This routine applies the PGMRES method.

191:    Input Parameter:
192: .     ksp - the Krylov space object that was set to use pgmres

194:    Output Parameter:
195: .     outits - number of iterations used

197: */
198: static PetscErrorCode KSPSolve_PGMRES(KSP ksp)
199: {
200:   PetscInt    its, itcount;
201:   KSP_PGMRES *pgmres     = (KSP_PGMRES *)ksp->data;
202:   PetscBool   guess_zero = ksp->guess_zero;

204:   PetscFunctionBegin;
205:   PetscCheck(!ksp->calc_sings || pgmres->Rsvd, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ORDER, "Must call KSPSetComputeSingularValues() before KSPSetUp() is called");
206:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
207:   ksp->its = 0;
208:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));

210:   itcount     = 0;
211:   ksp->reason = KSP_CONVERGED_ITERATING;
212:   while (!ksp->reason) {
213:     PetscCall(KSPInitialResidual(ksp, ksp->vec_sol, VEC_TEMP, VEC_TEMP_MATOP, VEC_VV(0), ksp->vec_rhs));
214:     PetscCall(KSPPGMRESCycle(&its, ksp));
215:     itcount += its;
216:     if (itcount >= ksp->max_it) {
217:       if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
218:       break;
219:     }
220:     ksp->guess_zero = PETSC_FALSE; /* every future call to KSPInitialResidual() will have nonzero guess */
221:   }
222:   ksp->guess_zero = guess_zero; /* restore if user provided nonzero initial guess */
223:   PetscFunctionReturn(PETSC_SUCCESS);
224: }

226: static PetscErrorCode KSPDestroy_PGMRES(KSP ksp)
227: {
228:   PetscFunctionBegin;
229:   PetscCall(KSPDestroy_GMRES(ksp));
230:   PetscFunctionReturn(PETSC_SUCCESS);
231: }

233: /*
234:     KSPPGMRESBuildSoln - create the solution from the starting vector and the
235:                       current iterates.

237:     Input parameters:
238:         nrs - work area of size it + 1.
239:         vguess  - index of initial guess
240:         vdest - index of result.  Note that vguess may == vdest (replace
241:                 guess with the solution).
242:         it - HH upper triangular part is a block of size (it+1) x (it+1)

244:      This is an internal routine that knows about the PGMRES internals.
245:  */
246: static PetscErrorCode KSPPGMRESBuildSoln(PetscScalar *nrs, Vec vguess, Vec vdest, KSP ksp, PetscInt it)
247: {
248:   PetscScalar tt;
249:   PetscInt    k, j;
250:   KSP_PGMRES *pgmres = (KSP_PGMRES *)(ksp->data);

252:   PetscFunctionBegin;
253:   /* Solve for solution vector that minimizes the residual */

255:   if (it < 0) {                        /* no pgmres steps have been performed */
256:     PetscCall(VecCopy(vguess, vdest)); /* VecCopy() is smart, exits immediately if vguess == vdest */
257:     PetscFunctionReturn(PETSC_SUCCESS);
258:   }

260:   /* solve the upper triangular system - RS is the right side and HH is
261:      the upper triangular matrix  - put soln in nrs */
262:   if (*HH(it, it) != 0.0) nrs[it] = *RS(it) / *HH(it, it);
263:   else nrs[it] = 0.0;

265:   for (k = it - 1; k >= 0; k--) {
266:     tt = *RS(k);
267:     for (j = k + 1; j <= it; j++) tt -= *HH(k, j) * nrs[j];
268:     nrs[k] = tt / *HH(k, k);
269:   }

271:   /* Accumulate the correction to the solution of the preconditioned problem in TEMP */
272:   PetscCall(VecZeroEntries(VEC_TEMP));
273:   PetscCall(VecMAXPY(VEC_TEMP, it + 1, nrs, &VEC_VV(0)));
274:   PetscCall(KSPUnwindPreconditioner(ksp, VEC_TEMP, VEC_TEMP_MATOP));
275:   /* add solution to previous solution */
276:   if (vdest == vguess) {
277:     PetscCall(VecAXPY(vdest, 1.0, VEC_TEMP));
278:   } else {
279:     PetscCall(VecWAXPY(vdest, 1.0, VEC_TEMP, vguess));
280:   }
281:   PetscFunctionReturn(PETSC_SUCCESS);
282: }

284: /*

286:     KSPPGMRESUpdateHessenberg - Do the scalar work for the orthogonalization.
287:                             Return new residual.

289:     input parameters:

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

296:     output parameters:
297: .        res - the new residual

299:  */
300: /*
301: .  it - column of the Hessenberg that is complete, PGMRES is actually computing two columns ahead of this
302:  */
303: static PetscErrorCode KSPPGMRESUpdateHessenberg(KSP ksp, PetscInt it, PetscBool *hapend, PetscReal *res)
304: {
305:   PetscScalar *hh, *cc, *ss, *rs;
306:   PetscInt     j;
307:   PetscReal    hapbnd;
308:   KSP_PGMRES  *pgmres = (KSP_PGMRES *)(ksp->data);

310:   PetscFunctionBegin;
311:   hh = HH(0, it); /* pointer to beginning of column to update */
312:   cc = CC(0);     /* beginning of cosine rotations */
313:   ss = SS(0);     /* beginning of sine rotations */
314:   rs = RS(0);     /* right hand side of least squares system */

316:   /* The Hessenberg matrix is now correct through column it, save that form for possible spectral analysis */
317:   for (j = 0; j <= it + 1; j++) *HES(j, it) = hh[j];

319:   /* check for the happy breakdown */
320:   hapbnd = PetscMin(PetscAbsScalar(hh[it + 1] / rs[it]), pgmres->haptol);
321:   if (PetscAbsScalar(hh[it + 1]) < hapbnd) {
322:     PetscCall(PetscInfo(ksp, "Detected happy breakdown, current hapbnd = %14.12e H(%" PetscInt_FMT ",%" PetscInt_FMT ") = %14.12e\n", (double)hapbnd, it + 1, it, (double)PetscAbsScalar(*HH(it + 1, it))));
323:     *hapend = PETSC_TRUE;
324:   }

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

331:   for (j = 0; j < it; j++) {
332:     PetscScalar hhj = hh[j];
333:     hh[j]           = PetscConj(cc[j]) * hhj + ss[j] * hh[j + 1];
334:     hh[j + 1]       = -ss[j] * hhj + cc[j] * hh[j + 1];
335:   }

337:   /*
338:     compute the new plane rotation, and apply it to:
339:      1) the right-hand-side of the Hessenberg system (RS)
340:         note: it affects RS(it) and RS(it+1)
341:      2) the new column of the Hessenberg matrix
342:         note: it affects HH(it,it) which is currently pointed to
343:         by hh and HH(it+1, it) (*(hh+1))
344:     thus obtaining the updated value of the residual...
345:   */

347:   /* compute new plane rotation */

349:   if (!*hapend) {
350:     PetscReal delta = PetscSqrtReal(PetscSqr(PetscAbsScalar(hh[it])) + PetscSqr(PetscAbsScalar(hh[it + 1])));
351:     if (delta == 0.0) {
352:       ksp->reason = KSP_DIVERGED_NULL;
353:       PetscFunctionReturn(PETSC_SUCCESS);
354:     }

356:     cc[it] = hh[it] / delta;     /* new cosine value */
357:     ss[it] = hh[it + 1] / delta; /* new sine value */

359:     hh[it]     = PetscConj(cc[it]) * hh[it] + ss[it] * hh[it + 1];
360:     rs[it + 1] = -ss[it] * rs[it];
361:     rs[it]     = PetscConj(cc[it]) * rs[it];
362:     *res       = PetscAbsScalar(rs[it + 1]);
363:   } else { /* happy breakdown: HH(it+1, it) = 0, therefore we don't need to apply
364:             another rotation matrix (so RH doesn't change).  The new residual is
365:             always the new sine term times the residual from last time (RS(it)),
366:             but now the new sine rotation would be zero...so the residual should
367:             be zero...so we will multiply "zero" by the last residual.  This might
368:             not be exactly what we want to do here -could just return "zero". */

370:     *res = 0.0;
371:   }
372:   PetscFunctionReturn(PETSC_SUCCESS);
373: }

375: /*
376:    KSPBuildSolution_PGMRES

378:      Input Parameter:
379: .     ksp - the Krylov space object
380: .     ptr-

382:    Output Parameter:
383: .     result - the solution

385:    Note: this calls KSPPGMRESBuildSoln - the same function that KSPPGMRESCycle
386:    calls directly.

388: */
389: PetscErrorCode KSPBuildSolution_PGMRES(KSP ksp, Vec ptr, Vec *result)
390: {
391:   KSP_PGMRES *pgmres = (KSP_PGMRES *)ksp->data;

393:   PetscFunctionBegin;
394:   if (!ptr) {
395:     if (!pgmres->sol_temp) PetscCall(VecDuplicate(ksp->vec_sol, &pgmres->sol_temp));
396:     ptr = pgmres->sol_temp;
397:   }
398:   if (!pgmres->nrs) {
399:     /* allocate the work area */
400:     PetscCall(PetscMalloc1(pgmres->max_k, &pgmres->nrs));
401:   }

403:   PetscCall(KSPPGMRESBuildSoln(pgmres->nrs, ksp->vec_sol, ptr, ksp, pgmres->it));
404:   if (result) *result = ptr;
405:   PetscFunctionReturn(PETSC_SUCCESS);
406: }

408: PetscErrorCode KSPSetFromOptions_PGMRES(KSP ksp, PetscOptionItems *PetscOptionsObject)
409: {
410:   PetscFunctionBegin;
411:   PetscCall(KSPSetFromOptions_GMRES(ksp, PetscOptionsObject));
412:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP pipelined GMRES Options");
413:   PetscOptionsHeadEnd();
414:   PetscFunctionReturn(PETSC_SUCCESS);
415: }

417: PetscErrorCode KSPReset_PGMRES(KSP ksp)
418: {
419:   PetscFunctionBegin;
420:   PetscCall(KSPReset_GMRES(ksp));
421:   PetscFunctionReturn(PETSC_SUCCESS);
422: }

424: /*MC
425:      KSPPGMRES - Implements the Pipelined Generalized Minimal Residual method. [](sec_pipelineksp)

427:    Options Database Keys:
428: +   -ksp_gmres_restart <restart> - the number of Krylov directions to orthogonalize against
429: .   -ksp_gmres_haptol <tol> - sets the tolerance for "happy ending" (exact convergence)
430: .   -ksp_gmres_preallocate - preallocate all the Krylov search directions initially (otherwise groups of
431:                              vectors are allocated as needed)
432: .   -ksp_gmres_classicalgramschmidt - use classical (unmodified) Gram-Schmidt to orthogonalize against the Krylov space (fast) (the default)
433: .   -ksp_gmres_modifiedgramschmidt - use modified Gram-Schmidt in the orthogonalization (more stable, but slower)
434: .   -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always> - determine if iterative refinement is used to increase the
435:                                    stability of the classical Gram-Schmidt  orthogonalization.
436: -   -ksp_gmres_krylov_monitor - plot the Krylov space generated

438:    Level: beginner

440:    Note:
441:    MPI configuration may be necessary for reductions to make asynchronous progress, which is important for performance of pipelined methods.
442:    See [](doc_faq_pipelined)

444:    Reference:
445:    Ghysels, Ashby, Meerbergen, Vanroose, Hiding global communication latencies in the GMRES algorithm on massively parallel machines, 2012.

447:    Developer Note:
448:     This object is subclassed off of `KSPGMRES`

450: .seealso: [](ch_ksp), [](sec_pipelineksp), [](doc_faq_pipelined), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPGMRES`, `KSPLGMRES`, `KSPPIPECG`, `KSPPIPECR`,
451:           `KSPGMRESSetRestart()`, `KSPGMRESSetHapTol()`, `KSPGMRESSetPreAllocateVectors()`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESGetOrthogonalization()`,
452:           `KSPGMRESClassicalGramSchmidtOrthogonalization()`, `KSPGMRESModifiedGramSchmidtOrthogonalization()`,
453:           `KSPGMRESCGSRefinementType`, `KSPGMRESSetCGSRefinementType()`, `KSPGMRESGetCGSRefinementType()`, `KSPGMRESMonitorKrylov()`
454: M*/

456: PETSC_EXTERN PetscErrorCode KSPCreate_PGMRES(KSP ksp)
457: {
458:   KSP_PGMRES *pgmres;

460:   PetscFunctionBegin;
461:   PetscCall(PetscNew(&pgmres));

463:   ksp->data                              = (void *)pgmres;
464:   ksp->ops->buildsolution                = KSPBuildSolution_PGMRES;
465:   ksp->ops->setup                        = KSPSetUp_PGMRES;
466:   ksp->ops->solve                        = KSPSolve_PGMRES;
467:   ksp->ops->reset                        = KSPReset_PGMRES;
468:   ksp->ops->destroy                      = KSPDestroy_PGMRES;
469:   ksp->ops->view                         = KSPView_GMRES;
470:   ksp->ops->setfromoptions               = KSPSetFromOptions_PGMRES;
471:   ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_GMRES;
472:   ksp->ops->computeeigenvalues           = KSPComputeEigenvalues_GMRES;

474:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
475:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2));
476:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));

478:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetPreAllocateVectors_C", KSPGMRESSetPreAllocateVectors_GMRES));
479:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetOrthogonalization_C", KSPGMRESSetOrthogonalization_GMRES));
480:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetOrthogonalization_C", KSPGMRESGetOrthogonalization_GMRES));
481:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetRestart_C", KSPGMRESSetRestart_GMRES));
482:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetRestart_C", KSPGMRESGetRestart_GMRES));
483:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetCGSRefinementType_C", KSPGMRESSetCGSRefinementType_GMRES));
484:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetCGSRefinementType_C", KSPGMRESGetCGSRefinementType_GMRES));

486:   pgmres->nextra_vecs    = 1;
487:   pgmres->haptol         = 1.0e-30;
488:   pgmres->q_preallocate  = 0;
489:   pgmres->delta_allocate = PGMRES_DELTA_DIRECTIONS;
490:   pgmres->orthog         = KSPGMRESClassicalGramSchmidtOrthogonalization;
491:   pgmres->nrs            = NULL;
492:   pgmres->sol_temp       = NULL;
493:   pgmres->max_k          = PGMRES_DEFAULT_MAXK;
494:   pgmres->Rsvd           = NULL;
495:   pgmres->orthogwork     = NULL;
496:   pgmres->cgstype        = KSP_GMRES_CGS_REFINE_NEVER;
497:   PetscFunctionReturn(PETSC_SUCCESS);
498: }