xref: /petsc/src/ksp/ksp/impls/gmres/pgmres/pgmres.c (revision 6479e8359b1ee10dd95e9e66ea2aeaf5277fb0c9)
1 
2 /*
3     This file implements PGMRES (a Pipelined Generalized Minimal Residual method)
4 */
5 
6 #include <../src/ksp/ksp/impls/gmres/pgmres/pgmresimpl.h> /*I  "petscksp.h"  I*/
7 #define PGMRES_DELTA_DIRECTIONS 10
8 #define PGMRES_DEFAULT_MAXK     30
9 
10 static PetscErrorCode KSPPGMRESUpdateHessenberg(KSP, PetscInt, PetscBool *, PetscReal *);
11 static PetscErrorCode KSPPGMRESBuildSoln(PetscScalar *, Vec, Vec, KSP, PetscInt);
12 
13 /*
14 
15     KSPSetUp_PGMRES - Sets up the workspace needed by pgmres.
16 
17     This is called once, usually automatically by KSPSolve() or KSPSetUp(),
18     but can be called directly by KSPSetUp().
19 
20 */
21 static PetscErrorCode KSPSetUp_PGMRES(KSP ksp)
22 {
23   PetscFunctionBegin;
24   PetscCall(KSPSetUp_GMRES(ksp));
25   PetscFunctionReturn(PETSC_SUCCESS);
26 }
27 
28 /*
29 
30     KSPPGMRESCycle - Run pgmres, possibly with restart.  Return residual
31                   history if requested.
32 
33     input parameters:
34 .        pgmres  - structure containing parameters and work areas
35 
36     output parameters:
37 .        itcount - number of iterations used.  If null, ignored.
38 .        converged - 0 if not converged
39 
40     Notes:
41     On entry, the value in vector VEC_VV(0) should be
42     the initial residual.
43 
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;
51 
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;
58 
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   }
72 
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. */
81 
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     }
85 
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     }
93 
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)));
97 
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;
103 
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       }
116 
117       if (!(it < pgmres->max_k + 1 && ksp->its < ksp->max_it)) break;
118 
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)));
123 
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     }
130 
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));
160 
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     }
168 
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)));
171 
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   }
175 
176   if (itcount) *itcount = it - 1; /* Number of iterations actually completed. */
177 
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 }
187 
188 /*
189     KSPSolve_PGMRES - This routine applies the PGMRES method.
190 
191    Input Parameter:
192 .     ksp - the Krylov space object that was set to use pgmres
193 
194    Output Parameter:
195 .     outits - number of iterations used
196 
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;
203 
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));
209 
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 }
225 
226 static PetscErrorCode KSPDestroy_PGMRES(KSP ksp)
227 {
228   PetscFunctionBegin;
229   PetscCall(KSPDestroy_GMRES(ksp));
230   PetscFunctionReturn(PETSC_SUCCESS);
231 }
232 
233 /*
234     KSPPGMRESBuildSoln - create the solution from the starting vector and the
235                       current iterates.
236 
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)
243 
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);
251 
252   PetscFunctionBegin;
253   /* Solve for solution vector that minimizes the residual */
254 
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   }
259 
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;
264 
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   }
270 
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 }
283 
284 /*
285 
286     KSPPGMRESUpdateHessenberg - Do the scalar work for the orthogonalization.
287                             Return new residual.
288 
289     input parameters:
290 
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.
295 
296     output parameters:
297 .        res - the new residual
298 
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);
309 
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 */
315 
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];
318 
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   }
325 
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!) */
330 
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   }
336 
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   */
346 
347   /* compute new plane rotation */
348 
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     }
355 
356     cc[it] = hh[it] / delta;     /* new cosine value */
357     ss[it] = hh[it + 1] / delta; /* new sine value */
358 
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". */
369 
370     *res = 0.0;
371   }
372   PetscFunctionReturn(PETSC_SUCCESS);
373 }
374 
375 /*
376    KSPBuildSolution_PGMRES
377 
378      Input Parameter:
379 .     ksp - the Krylov space object
380 .     ptr-
381 
382    Output Parameter:
383 .     result - the solution
384 
385    Note: this calls KSPPGMRESBuildSoln - the same function that KSPPGMRESCycle
386    calls directly.
387 
388 */
389 PetscErrorCode KSPBuildSolution_PGMRES(KSP ksp, Vec ptr, Vec *result)
390 {
391   KSP_PGMRES *pgmres = (KSP_PGMRES *)ksp->data;
392 
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   }
402 
403   PetscCall(KSPPGMRESBuildSoln(pgmres->nrs, ksp->vec_sol, ptr, ksp, pgmres->it));
404   if (result) *result = ptr;
405   PetscFunctionReturn(PETSC_SUCCESS);
406 }
407 
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 }
416 
417 PetscErrorCode KSPReset_PGMRES(KSP ksp)
418 {
419   PetscFunctionBegin;
420   PetscCall(KSPReset_GMRES(ksp));
421   PetscFunctionReturn(PETSC_SUCCESS);
422 }
423 
424 /*MC
425      KSPPGMRES - Implements the Pipelined Generalized Minimal Residual method. [](sec_pipelineksp)
426 
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
437 
438    Level: beginner
439 
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)
443 
444    Reference:
445    Ghysels, Ashby, Meerbergen, Vanroose, Hiding global communication latencies in the GMRES algorithm on massively parallel machines, 2012.
446 
447    Developer Note:
448     This object is subclassed off of `KSPGMRES`
449 
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*/
455 
456 PETSC_EXTERN PetscErrorCode KSPCreate_PGMRES(KSP ksp)
457 {
458   KSP_PGMRES *pgmres;
459 
460   PetscFunctionBegin;
461   PetscCall(PetscNew(&pgmres));
462 
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;
473 
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));
477 
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));
485 
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 }
499