1
2 /*
3 This file implements GMRES (a Generalized Minimal Residual) method.
4 Reference: Saad and Schultz, 1986.
5
6
7 Some comments on left vs. right preconditioning, and restarts.
8 Left and right preconditioning.
9 If right preconditioning is chosen, then the problem being solved
10 by gmres is actually
11 My = AB^-1 y = f
12 so the initial residual is
13 r = f - Mx
14 Note that B^-1 y = x or y = B x, and if x is non-zero, the initial
15 residual is
16 r = f - A x
17 The final solution is then
18 x = B^-1 y
19
20 If left preconditioning is chosen, then the problem being solved is
21 My = B^-1 A x = B^-1 f,
22 and the initial residual is
23 r = B^-1(f - Ax)
24
25 Restarts: Restarts are basically solves with x0 not equal to zero.
26 Note that we can eliminate an extra application of B^-1 between
27 restarts as long as we don't require that the solution at the end
28 of an unsuccessful gmres iteration always be the solution x.
29 */
30
31 #include <../src/ksp/ksp/impls/gmres/gmresimpl.h> /*I "petscksp.h" I*/
32 #define GMRES_DELTA_DIRECTIONS 10
33 #define GMRES_DEFAULT_MAXK 30
34 static PetscErrorCode KSPGMRESUpdateHessenberg(KSP,PetscInt,PetscBool,PetscReal*);
35 static PetscErrorCode KSPGMRESBuildSoln(PetscScalar*,Vec,Vec,KSP,PetscInt);
36
KSPSetUp_GMRES(KSP ksp)37 PetscErrorCode KSPSetUp_GMRES(KSP ksp)
38 {
39 PetscInt hh,hes,rs,cc;
40 PetscErrorCode ierr;
41 PetscInt max_k,k;
42 KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
43
44 PetscFunctionBegin;
45 max_k = gmres->max_k; /* restart size */
46 hh = (max_k + 2) * (max_k + 1);
47 hes = (max_k + 1) * (max_k + 1);
48 rs = (max_k + 2);
49 cc = (max_k + 1);
50
51 ierr = PetscCalloc5(hh,&gmres->hh_origin,hes,&gmres->hes_origin,rs,&gmres->rs_origin,cc,&gmres->cc_origin,cc,&gmres->ss_origin);CHKERRQ(ierr);
52 ierr = PetscLogObjectMemory((PetscObject)ksp,(hh + hes + rs + 2*cc)*sizeof(PetscScalar));CHKERRQ(ierr);
53
54 if (ksp->calc_sings) {
55 /* Allocate workspace to hold Hessenberg matrix needed by lapack */
56 ierr = PetscMalloc1((max_k + 3)*(max_k + 9),&gmres->Rsvd);CHKERRQ(ierr);
57 ierr = PetscLogObjectMemory((PetscObject)ksp,(max_k + 3)*(max_k + 9)*sizeof(PetscScalar));CHKERRQ(ierr);
58 ierr = PetscMalloc1(6*(max_k+2),&gmres->Dsvd);CHKERRQ(ierr);
59 ierr = PetscLogObjectMemory((PetscObject)ksp,6*(max_k+2)*sizeof(PetscReal));CHKERRQ(ierr);
60 }
61
62 /* Allocate array to hold pointers to user vectors. Note that we need
63 4 + max_k + 1 (since we need it+1 vectors, and it <= max_k) */
64 gmres->vecs_allocated = VEC_OFFSET + 2 + max_k + gmres->nextra_vecs;
65
66 ierr = PetscMalloc1(gmres->vecs_allocated,&gmres->vecs);CHKERRQ(ierr);
67 ierr = PetscMalloc1(VEC_OFFSET+2+max_k,&gmres->user_work);CHKERRQ(ierr);
68 ierr = PetscMalloc1(VEC_OFFSET+2+max_k,&gmres->mwork_alloc);CHKERRQ(ierr);
69 ierr = PetscLogObjectMemory((PetscObject)ksp,(VEC_OFFSET+2+max_k)*(sizeof(Vec*)+sizeof(PetscInt)) + gmres->vecs_allocated*sizeof(Vec));CHKERRQ(ierr);
70
71 if (gmres->q_preallocate) {
72 gmres->vv_allocated = VEC_OFFSET + 2 + max_k;
73
74 ierr = KSPCreateVecs(ksp,gmres->vv_allocated,&gmres->user_work[0],0,NULL);CHKERRQ(ierr);
75 ierr = PetscLogObjectParents(ksp,gmres->vv_allocated,gmres->user_work[0]);CHKERRQ(ierr);
76
77 gmres->mwork_alloc[0] = gmres->vv_allocated;
78 gmres->nwork_alloc = 1;
79 for (k=0; k<gmres->vv_allocated; k++) {
80 gmres->vecs[k] = gmres->user_work[0][k];
81 }
82 } else {
83 gmres->vv_allocated = 5;
84
85 ierr = KSPCreateVecs(ksp,5,&gmres->user_work[0],0,NULL);CHKERRQ(ierr);
86 ierr = PetscLogObjectParents(ksp,5,gmres->user_work[0]);CHKERRQ(ierr);
87
88 gmres->mwork_alloc[0] = 5;
89 gmres->nwork_alloc = 1;
90 for (k=0; k<gmres->vv_allocated; k++) {
91 gmres->vecs[k] = gmres->user_work[0][k];
92 }
93 }
94 PetscFunctionReturn(0);
95 }
96
97 /*
98 Run gmres, possibly with restart. Return residual history if requested.
99 input parameters:
100
101 . gmres - structure containing parameters and work areas
102
103 output parameters:
104 . nres - residuals (from preconditioned system) at each step.
105 If restarting, consider passing nres+it. If null,
106 ignored
107 . itcount - number of iterations used. nres[0] to nres[itcount]
108 are defined. If null, ignored.
109
110 Notes:
111 On entry, the value in vector VEC_VV(0) should be the initial residual
112 (this allows shortcuts where the initial preconditioned residual is 0).
113 */
KSPGMRESCycle(PetscInt * itcount,KSP ksp)114 PetscErrorCode KSPGMRESCycle(PetscInt *itcount,KSP ksp)
115 {
116 KSP_GMRES *gmres = (KSP_GMRES*)(ksp->data);
117 PetscReal res,hapbnd,tt;
118 PetscErrorCode ierr;
119 PetscInt it = 0, max_k = gmres->max_k;
120 PetscBool hapend = PETSC_FALSE;
121
122 PetscFunctionBegin;
123 if (itcount) *itcount = 0;
124 ierr = VecNormalize(VEC_VV(0),&res);CHKERRQ(ierr);
125 KSPCheckNorm(ksp,res);
126
127 /* the constant .1 is arbitrary, just some measure at how incorrect the residuals are */
128 if ((ksp->rnorm > 0.0) && (PetscAbsReal(res-ksp->rnorm) > .1*gmres->rnorm0)) {
129 if (ksp->errorifnotconverged) SETERRQ3(PetscObjectComm((PetscObject)ksp),PETSC_ERR_CONV_FAILED,"Residual norm computed by GMRES recursion formula %g is far from the computed residual norm %g at restart, residual norm at start of cycle %g",(double)ksp->rnorm,(double)res,(double)gmres->rnorm0);
130 else {
131 ierr = PetscInfo3(ksp,"Residual norm computed by GMRES recursion formula %g is far from the computed residual norm %g at restart, residual norm at start of cycle %g",(double)ksp->rnorm,(double)res,(double)gmres->rnorm0);
132 ksp->reason = KSP_DIVERGED_BREAKDOWN;
133 PetscFunctionReturn(0);
134 }
135 }
136 *GRS(0) = gmres->rnorm0 = res;
137
138 /* check for the convergence */
139 ierr = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr);
140 ksp->rnorm = res;
141 ierr = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr);
142 gmres->it = (it - 1);
143 ierr = KSPLogResidualHistory(ksp,res);CHKERRQ(ierr);
144 ierr = KSPMonitor(ksp,ksp->its,res);CHKERRQ(ierr);
145 if (!res) {
146 ksp->reason = KSP_CONVERGED_ATOL;
147 ierr = PetscInfo(ksp,"Converged due to zero residual norm on entry\n");CHKERRQ(ierr);
148 PetscFunctionReturn(0);
149 }
150
151 ierr = (*ksp->converged)(ksp,ksp->its,res,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
152 while (!ksp->reason && it < max_k && ksp->its < ksp->max_it) {
153 if (it) {
154 ierr = KSPLogResidualHistory(ksp,res);CHKERRQ(ierr);
155 ierr = KSPMonitor(ksp,ksp->its,res);CHKERRQ(ierr);
156 }
157 gmres->it = (it - 1);
158 if (gmres->vv_allocated <= it + VEC_OFFSET + 1) {
159 ierr = KSPGMRESGetNewVectors(ksp,it+1);CHKERRQ(ierr);
160 }
161 ierr = KSP_PCApplyBAorAB(ksp,VEC_VV(it),VEC_VV(1+it),VEC_TEMP_MATOP);CHKERRQ(ierr);
162
163 /* update hessenberg matrix and do Gram-Schmidt */
164 ierr = (*gmres->orthog)(ksp,it);CHKERRQ(ierr);
165 if (ksp->reason) break;
166
167 /* vv(i+1) . vv(i+1) */
168 ierr = VecNormalize(VEC_VV(it+1),&tt);CHKERRQ(ierr);
169 KSPCheckNorm(ksp,tt);
170
171 /* save the magnitude */
172 *HH(it+1,it) = tt;
173 *HES(it+1,it) = tt;
174
175 /* check for the happy breakdown */
176 hapbnd = PetscAbsScalar(tt / *GRS(it));
177 if (hapbnd > gmres->haptol) hapbnd = gmres->haptol;
178 if (tt < hapbnd) {
179 ierr = PetscInfo2(ksp,"Detected happy breakdown, current hapbnd = %14.12e tt = %14.12e\n",(double)hapbnd,(double)tt);CHKERRQ(ierr);
180 hapend = PETSC_TRUE;
181 }
182 ierr = KSPGMRESUpdateHessenberg(ksp,it,hapend,&res);CHKERRQ(ierr);
183
184 it++;
185 gmres->it = (it-1); /* For converged */
186 ksp->its++;
187 ksp->rnorm = res;
188 if (ksp->reason) break;
189
190 ierr = (*ksp->converged)(ksp,ksp->its,res,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
191
192 /* Catch error in happy breakdown and signal convergence and break from loop */
193 if (hapend) {
194 if (ksp->normtype == KSP_NORM_NONE) { /* convergence test was skipped in this case */
195 ksp->reason = KSP_CONVERGED_HAPPY_BREAKDOWN;
196 } else if (!ksp->reason) {
197 if (ksp->errorifnotconverged) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"You reached the happy break down, but convergence was not indicated. Residual norm = %g",(double)res);
198 else {
199 ksp->reason = KSP_DIVERGED_BREAKDOWN;
200 break;
201 }
202 }
203 }
204 }
205
206 /* Monitor if we know that we will not return for a restart */
207 if (it && (ksp->reason || ksp->its >= ksp->max_it)) {
208 ierr = KSPLogResidualHistory(ksp,res);CHKERRQ(ierr);
209 ierr = KSPMonitor(ksp,ksp->its,res);CHKERRQ(ierr);
210 }
211
212 if (itcount) *itcount = it;
213
214
215 /*
216 Down here we have to solve for the "best" coefficients of the Krylov
217 columns, add the solution values together, and possibly unwind the
218 preconditioning from the solution
219 */
220 /* Form the solution (or the solution so far) */
221 ierr = KSPGMRESBuildSoln(GRS(0),ksp->vec_sol,ksp->vec_sol,ksp,it-1);CHKERRQ(ierr);
222 PetscFunctionReturn(0);
223 }
224
KSPSolve_GMRES(KSP ksp)225 PetscErrorCode KSPSolve_GMRES(KSP ksp)
226 {
227 PetscErrorCode ierr;
228 PetscInt its,itcount,i;
229 KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
230 PetscBool guess_zero = ksp->guess_zero;
231 PetscInt N = gmres->max_k + 1;
232
233 PetscFunctionBegin;
234 if (ksp->calc_sings && !gmres->Rsvd) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ORDER,"Must call KSPSetComputeSingularValues() before KSPSetUp() is called");
235
236 ierr = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr);
237 ksp->its = 0;
238 ierr = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr);
239
240 itcount = 0;
241 gmres->fullcycle = 0;
242 ksp->reason = KSP_CONVERGED_ITERATING;
243 ksp->rnorm = -1.0; /* special marker for KSPGMRESCycle() */
244 while (!ksp->reason) {
245 ierr = KSPInitialResidual(ksp,ksp->vec_sol,VEC_TEMP,VEC_TEMP_MATOP,VEC_VV(0),ksp->vec_rhs);CHKERRQ(ierr);
246 ierr = KSPGMRESCycle(&its,ksp);CHKERRQ(ierr);
247 /* Store the Hessenberg matrix and the basis vectors of the Krylov subspace
248 if the cycle is complete for the computation of the Ritz pairs */
249 if (its == gmres->max_k) {
250 gmres->fullcycle++;
251 if (ksp->calc_ritz) {
252 if (!gmres->hes_ritz) {
253 ierr = PetscMalloc1(N*N,&gmres->hes_ritz);CHKERRQ(ierr);
254 ierr = PetscLogObjectMemory((PetscObject)ksp,N*N*sizeof(PetscScalar));CHKERRQ(ierr);
255 ierr = VecDuplicateVecs(VEC_VV(0),N,&gmres->vecb);CHKERRQ(ierr);
256 }
257 ierr = PetscArraycpy(gmres->hes_ritz,gmres->hes_origin,N*N);CHKERRQ(ierr);
258 for (i=0; i<gmres->max_k+1; i++) {
259 ierr = VecCopy(VEC_VV(i),gmres->vecb[i]);CHKERRQ(ierr);
260 }
261 }
262 }
263 itcount += its;
264 if (itcount >= ksp->max_it) {
265 if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
266 break;
267 }
268 ksp->guess_zero = PETSC_FALSE; /* every future call to KSPInitialResidual() will have nonzero guess */
269 }
270 ksp->guess_zero = guess_zero; /* restore if user provided nonzero initial guess */
271 PetscFunctionReturn(0);
272 }
273
KSPReset_GMRES(KSP ksp)274 PetscErrorCode KSPReset_GMRES(KSP ksp)
275 {
276 KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
277 PetscErrorCode ierr;
278 PetscInt i;
279
280 PetscFunctionBegin;
281 /* Free the Hessenberg matrices */
282 ierr = PetscFree5(gmres->hh_origin,gmres->hes_origin,gmres->rs_origin,gmres->cc_origin,gmres->ss_origin);CHKERRQ(ierr);
283 ierr = PetscFree(gmres->hes_ritz);CHKERRQ(ierr);
284
285 /* free work vectors */
286 ierr = PetscFree(gmres->vecs);CHKERRQ(ierr);
287 for (i=0; i<gmres->nwork_alloc; i++) {
288 ierr = VecDestroyVecs(gmres->mwork_alloc[i],&gmres->user_work[i]);CHKERRQ(ierr);
289 }
290 gmres->nwork_alloc = 0;
291 if (gmres->vecb) {
292 ierr = VecDestroyVecs(gmres->max_k+1,&gmres->vecb);CHKERRQ(ierr);
293 }
294
295 ierr = PetscFree(gmres->user_work);CHKERRQ(ierr);
296 ierr = PetscFree(gmres->mwork_alloc);CHKERRQ(ierr);
297 ierr = PetscFree(gmres->nrs);CHKERRQ(ierr);
298 ierr = VecDestroy(&gmres->sol_temp);CHKERRQ(ierr);
299 ierr = PetscFree(gmres->Rsvd);CHKERRQ(ierr);
300 ierr = PetscFree(gmres->Dsvd);CHKERRQ(ierr);
301 ierr = PetscFree(gmres->orthogwork);CHKERRQ(ierr);
302
303 gmres->vv_allocated = 0;
304 gmres->vecs_allocated = 0;
305 gmres->sol_temp = NULL;
306 PetscFunctionReturn(0);
307 }
308
KSPDestroy_GMRES(KSP ksp)309 PetscErrorCode KSPDestroy_GMRES(KSP ksp)
310 {
311 PetscErrorCode ierr;
312
313 PetscFunctionBegin;
314 ierr = KSPReset_GMRES(ksp);CHKERRQ(ierr);
315 ierr = PetscFree(ksp->data);CHKERRQ(ierr);
316 /* clear composed functions */
317 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetPreAllocateVectors_C",NULL);CHKERRQ(ierr);
318 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetOrthogonalization_C",NULL);CHKERRQ(ierr);
319 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetOrthogonalization_C",NULL);CHKERRQ(ierr);
320 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetRestart_C",NULL);CHKERRQ(ierr);
321 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetRestart_C",NULL);CHKERRQ(ierr);
322 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetHapTol_C",NULL);CHKERRQ(ierr);
323 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetCGSRefinementType_C",NULL);CHKERRQ(ierr);
324 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetCGSRefinementType_C",NULL);CHKERRQ(ierr);
325 PetscFunctionReturn(0);
326 }
327 /*
328 KSPGMRESBuildSoln - create the solution from the starting vector and the
329 current iterates.
330
331 Input parameters:
332 nrs - work area of size it + 1.
333 vs - index of initial guess
334 vdest - index of result. Note that vs may == vdest (replace
335 guess with the solution).
336
337 This is an internal routine that knows about the GMRES internals.
338 */
KSPGMRESBuildSoln(PetscScalar * nrs,Vec vs,Vec vdest,KSP ksp,PetscInt it)339 static PetscErrorCode KSPGMRESBuildSoln(PetscScalar *nrs,Vec vs,Vec vdest,KSP ksp,PetscInt it)
340 {
341 PetscScalar tt;
342 PetscErrorCode ierr;
343 PetscInt ii,k,j;
344 KSP_GMRES *gmres = (KSP_GMRES*)(ksp->data);
345
346 PetscFunctionBegin;
347 /* Solve for solution vector that minimizes the residual */
348
349 /* If it is < 0, no gmres steps have been performed */
350 if (it < 0) {
351 ierr = VecCopy(vs,vdest);CHKERRQ(ierr); /* VecCopy() is smart, exists immediately if vguess == vdest */
352 PetscFunctionReturn(0);
353 }
354 if (*HH(it,it) != 0.0) {
355 nrs[it] = *GRS(it) / *HH(it,it);
356 } else {
357 if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"You reached the break down in GMRES; HH(it,it) = 0");
358 else ksp->reason = KSP_DIVERGED_BREAKDOWN;
359
360 ierr = PetscInfo2(ksp,"Likely your matrix or preconditioner is singular. HH(it,it) is identically zero; it = %D GRS(it) = %g\n",it,(double)PetscAbsScalar(*GRS(it)));CHKERRQ(ierr);
361 PetscFunctionReturn(0);
362 }
363 for (ii=1; ii<=it; ii++) {
364 k = it - ii;
365 tt = *GRS(k);
366 for (j=k+1; j<=it; j++) tt = tt - *HH(k,j) * nrs[j];
367 if (*HH(k,k) == 0.0) {
368 if (ksp->errorifnotconverged) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"Likely your matrix or preconditioner is singular. HH(k,k) is identically zero; k = %D\n",k);
369 else {
370 ksp->reason = KSP_DIVERGED_BREAKDOWN;
371 ierr = PetscInfo1(ksp,"Likely your matrix or preconditioner is singular. HH(k,k) is identically zero; k = %D\n",k);CHKERRQ(ierr);
372 PetscFunctionReturn(0);
373 }
374 }
375 nrs[k] = tt / *HH(k,k);
376 }
377
378 /* Accumulate the correction to the solution of the preconditioned problem in TEMP */
379 ierr = VecSet(VEC_TEMP,0.0);CHKERRQ(ierr);
380 ierr = VecMAXPY(VEC_TEMP,it+1,nrs,&VEC_VV(0));CHKERRQ(ierr);
381
382 ierr = KSPUnwindPreconditioner(ksp,VEC_TEMP,VEC_TEMP_MATOP);CHKERRQ(ierr);
383 /* add solution to previous solution */
384 if (vdest != vs) {
385 ierr = VecCopy(vs,vdest);CHKERRQ(ierr);
386 }
387 ierr = VecAXPY(vdest,1.0,VEC_TEMP);CHKERRQ(ierr);
388 PetscFunctionReturn(0);
389 }
390 /*
391 Do the scalar work for the orthogonalization. Return new residual norm.
392 */
KSPGMRESUpdateHessenberg(KSP ksp,PetscInt it,PetscBool hapend,PetscReal * res)393 static PetscErrorCode KSPGMRESUpdateHessenberg(KSP ksp,PetscInt it,PetscBool hapend,PetscReal *res)
394 {
395 PetscScalar *hh,*cc,*ss,tt;
396 PetscInt j;
397 KSP_GMRES *gmres = (KSP_GMRES*)(ksp->data);
398
399 PetscFunctionBegin;
400 hh = HH(0,it);
401 cc = CC(0);
402 ss = SS(0);
403
404 /* Apply all the previously computed plane rotations to the new column
405 of the Hessenberg matrix */
406 for (j=1; j<=it; j++) {
407 tt = *hh;
408 *hh = PetscConj(*cc) * tt + *ss * *(hh+1);
409 hh++;
410 *hh = *cc++ * *hh - (*ss++ * tt);
411 }
412
413 /*
414 compute the new plane rotation, and apply it to:
415 1) the right-hand-side of the Hessenberg system
416 2) the new column of the Hessenberg matrix
417 thus obtaining the updated value of the residual
418 */
419 if (!hapend) {
420 tt = PetscSqrtScalar(PetscConj(*hh) * *hh + PetscConj(*(hh+1)) * *(hh+1));
421 if (tt == 0.0) {
422 if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"tt == 0.0");
423 else {
424 ksp->reason = KSP_DIVERGED_NULL;
425 PetscFunctionReturn(0);
426 }
427 }
428 *cc = *hh / tt;
429 *ss = *(hh+1) / tt;
430 *GRS(it+1) = -(*ss * *GRS(it));
431 *GRS(it) = PetscConj(*cc) * *GRS(it);
432 *hh = PetscConj(*cc) * *hh + *ss * *(hh+1);
433 *res = PetscAbsScalar(*GRS(it+1));
434 } else {
435 /* happy breakdown: HH(it+1, it) = 0, therfore we don't need to apply
436 another rotation matrix (so RH doesn't change). The new residual is
437 always the new sine term times the residual from last time (GRS(it)),
438 but now the new sine rotation would be zero...so the residual should
439 be zero...so we will multiply "zero" by the last residual. This might
440 not be exactly what we want to do here -could just return "zero". */
441
442 *res = 0.0;
443 }
444 PetscFunctionReturn(0);
445 }
446 /*
447 This routine allocates more work vectors, starting from VEC_VV(it).
448 */
KSPGMRESGetNewVectors(KSP ksp,PetscInt it)449 PetscErrorCode KSPGMRESGetNewVectors(KSP ksp,PetscInt it)
450 {
451 KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
452 PetscErrorCode ierr;
453 PetscInt nwork = gmres->nwork_alloc,k,nalloc;
454
455 PetscFunctionBegin;
456 nalloc = PetscMin(ksp->max_it,gmres->delta_allocate);
457 /* Adjust the number to allocate to make sure that we don't exceed the
458 number of available slots */
459 if (it + VEC_OFFSET + nalloc >= gmres->vecs_allocated) {
460 nalloc = gmres->vecs_allocated - it - VEC_OFFSET;
461 }
462 if (!nalloc) PetscFunctionReturn(0);
463
464 gmres->vv_allocated += nalloc;
465
466 ierr = KSPCreateVecs(ksp,nalloc,&gmres->user_work[nwork],0,NULL);CHKERRQ(ierr);
467 ierr = PetscLogObjectParents(ksp,nalloc,gmres->user_work[nwork]);CHKERRQ(ierr);
468
469 gmres->mwork_alloc[nwork] = nalloc;
470 for (k=0; k<nalloc; k++) {
471 gmres->vecs[it+VEC_OFFSET+k] = gmres->user_work[nwork][k];
472 }
473 gmres->nwork_alloc++;
474 PetscFunctionReturn(0);
475 }
476
KSPBuildSolution_GMRES(KSP ksp,Vec ptr,Vec * result)477 PetscErrorCode KSPBuildSolution_GMRES(KSP ksp,Vec ptr,Vec *result)
478 {
479 KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
480 PetscErrorCode ierr;
481
482 PetscFunctionBegin;
483 if (!ptr) {
484 if (!gmres->sol_temp) {
485 ierr = VecDuplicate(ksp->vec_sol,&gmres->sol_temp);CHKERRQ(ierr);
486 ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)gmres->sol_temp);CHKERRQ(ierr);
487 }
488 ptr = gmres->sol_temp;
489 }
490 if (!gmres->nrs) {
491 /* allocate the work area */
492 ierr = PetscMalloc1(gmres->max_k,&gmres->nrs);CHKERRQ(ierr);
493 ierr = PetscLogObjectMemory((PetscObject)ksp,gmres->max_k);CHKERRQ(ierr);
494 }
495
496 ierr = KSPGMRESBuildSoln(gmres->nrs,ksp->vec_sol,ptr,ksp,gmres->it);CHKERRQ(ierr);
497 if (result) *result = ptr;
498 PetscFunctionReturn(0);
499 }
500
KSPView_GMRES(KSP ksp,PetscViewer viewer)501 PetscErrorCode KSPView_GMRES(KSP ksp,PetscViewer viewer)
502 {
503 KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
504 const char *cstr;
505 PetscErrorCode ierr;
506 PetscBool iascii,isstring;
507
508 PetscFunctionBegin;
509 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
510 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr);
511 if (gmres->orthog == KSPGMRESClassicalGramSchmidtOrthogonalization) {
512 switch (gmres->cgstype) {
513 case (KSP_GMRES_CGS_REFINE_NEVER):
514 cstr = "Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement";
515 break;
516 case (KSP_GMRES_CGS_REFINE_ALWAYS):
517 cstr = "Classical (unmodified) Gram-Schmidt Orthogonalization with one step of iterative refinement";
518 break;
519 case (KSP_GMRES_CGS_REFINE_IFNEEDED):
520 cstr = "Classical (unmodified) Gram-Schmidt Orthogonalization with one step of iterative refinement when needed";
521 break;
522 default:
523 SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Unknown orthogonalization");
524 }
525 } else if (gmres->orthog == KSPGMRESModifiedGramSchmidtOrthogonalization) {
526 cstr = "Modified Gram-Schmidt Orthogonalization";
527 } else {
528 cstr = "unknown orthogonalization";
529 }
530 if (iascii) {
531 ierr = PetscViewerASCIIPrintf(viewer," restart=%D, using %s\n",gmres->max_k,cstr);CHKERRQ(ierr);
532 ierr = PetscViewerASCIIPrintf(viewer," happy breakdown tolerance %g\n",(double)gmres->haptol);CHKERRQ(ierr);
533 } else if (isstring) {
534 ierr = PetscViewerStringSPrintf(viewer,"%s restart %D",cstr,gmres->max_k);CHKERRQ(ierr);
535 }
536 PetscFunctionReturn(0);
537 }
538
539 /*@C
540 KSPGMRESMonitorKrylov - Calls VecView() for each new direction in the GMRES accumulated Krylov space.
541
542 Collective on ksp
543
544 Input Parameters:
545 + ksp - the KSP context
546 . its - iteration number
547 . fgnorm - 2-norm of residual (or gradient)
548 - dummy - an collection of viewers created with KSPViewerCreate()
549
550 Options Database Keys:
551 . -ksp_gmres_kyrlov_monitor
552
553 Notes:
554 A new PETSCVIEWERDRAW is created for each Krylov vector so they can all be simultaneously viewed
555 Level: intermediate
556
557 .seealso: KSPMonitorSet(), KSPMonitorDefault(), VecView(), KSPViewersCreate(), KSPViewersDestroy()
558 @*/
KSPGMRESMonitorKrylov(KSP ksp,PetscInt its,PetscReal fgnorm,void * dummy)559 PetscErrorCode KSPGMRESMonitorKrylov(KSP ksp,PetscInt its,PetscReal fgnorm,void *dummy)
560 {
561 PetscViewers viewers = (PetscViewers)dummy;
562 KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
563 PetscErrorCode ierr;
564 Vec x;
565 PetscViewer viewer;
566 PetscBool flg;
567
568 PetscFunctionBegin;
569 ierr = PetscViewersGetViewer(viewers,gmres->it+1,&viewer);CHKERRQ(ierr);
570 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&flg);CHKERRQ(ierr);
571 if (!flg) {
572 ierr = PetscViewerSetType(viewer,PETSCVIEWERDRAW);CHKERRQ(ierr);
573 ierr = PetscViewerDrawSetInfo(viewer,NULL,"Krylov GMRES Monitor",PETSC_DECIDE,PETSC_DECIDE,300,300);CHKERRQ(ierr);
574 }
575 x = VEC_VV(gmres->it+1);
576 ierr = VecView(x,viewer);CHKERRQ(ierr);
577 PetscFunctionReturn(0);
578 }
579
KSPSetFromOptions_GMRES(PetscOptionItems * PetscOptionsObject,KSP ksp)580 PetscErrorCode KSPSetFromOptions_GMRES(PetscOptionItems *PetscOptionsObject,KSP ksp)
581 {
582 PetscErrorCode ierr;
583 PetscInt restart;
584 PetscReal haptol;
585 KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
586 PetscBool flg;
587
588 PetscFunctionBegin;
589 ierr = PetscOptionsHead(PetscOptionsObject,"KSP GMRES Options");CHKERRQ(ierr);
590 ierr = PetscOptionsInt("-ksp_gmres_restart","Number of Krylov search directions","KSPGMRESSetRestart",gmres->max_k,&restart,&flg);CHKERRQ(ierr);
591 if (flg) { ierr = KSPGMRESSetRestart(ksp,restart);CHKERRQ(ierr); }
592 ierr = PetscOptionsReal("-ksp_gmres_haptol","Tolerance for exact convergence (happy ending)","KSPGMRESSetHapTol",gmres->haptol,&haptol,&flg);CHKERRQ(ierr);
593 if (flg) { ierr = KSPGMRESSetHapTol(ksp,haptol);CHKERRQ(ierr); }
594 flg = PETSC_FALSE;
595 ierr = PetscOptionsBool("-ksp_gmres_preallocate","Preallocate Krylov vectors","KSPGMRESSetPreAllocateVectors",flg,&flg,NULL);CHKERRQ(ierr);
596 if (flg) {ierr = KSPGMRESSetPreAllocateVectors(ksp);CHKERRQ(ierr);}
597 ierr = PetscOptionsBoolGroupBegin("-ksp_gmres_classicalgramschmidt","Classical (unmodified) Gram-Schmidt (fast)","KSPGMRESSetOrthogonalization",&flg);CHKERRQ(ierr);
598 if (flg) {ierr = KSPGMRESSetOrthogonalization(ksp,KSPGMRESClassicalGramSchmidtOrthogonalization);CHKERRQ(ierr);}
599 ierr = PetscOptionsBoolGroupEnd("-ksp_gmres_modifiedgramschmidt","Modified Gram-Schmidt (slow,more stable)","KSPGMRESSetOrthogonalization",&flg);CHKERRQ(ierr);
600 if (flg) {ierr = KSPGMRESSetOrthogonalization(ksp,KSPGMRESModifiedGramSchmidtOrthogonalization);CHKERRQ(ierr);}
601 ierr = PetscOptionsEnum("-ksp_gmres_cgs_refinement_type","Type of iterative refinement for classical (unmodified) Gram-Schmidt","KSPGMRESSetCGSRefinementType",
602 KSPGMRESCGSRefinementTypes,(PetscEnum)gmres->cgstype,(PetscEnum*)&gmres->cgstype,&flg);CHKERRQ(ierr);
603 flg = PETSC_FALSE;
604 ierr = PetscOptionsBool("-ksp_gmres_krylov_monitor","Plot the Krylov directions","KSPMonitorSet",flg,&flg,NULL);CHKERRQ(ierr);
605 if (flg) {
606 PetscViewers viewers;
607 ierr = PetscViewersCreate(PetscObjectComm((PetscObject)ksp),&viewers);CHKERRQ(ierr);
608 ierr = KSPMonitorSet(ksp,KSPGMRESMonitorKrylov,viewers,(PetscErrorCode (*)(void**))PetscViewersDestroy);CHKERRQ(ierr);
609 }
610 ierr = PetscOptionsTail();CHKERRQ(ierr);
611 PetscFunctionReturn(0);
612 }
613
KSPGMRESSetHapTol_GMRES(KSP ksp,PetscReal tol)614 PetscErrorCode KSPGMRESSetHapTol_GMRES(KSP ksp,PetscReal tol)
615 {
616 KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
617
618 PetscFunctionBegin;
619 if (tol < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Tolerance must be non-negative");
620 gmres->haptol = tol;
621 PetscFunctionReturn(0);
622 }
623
KSPGMRESGetRestart_GMRES(KSP ksp,PetscInt * max_k)624 PetscErrorCode KSPGMRESGetRestart_GMRES(KSP ksp,PetscInt *max_k)
625 {
626 KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
627
628 PetscFunctionBegin;
629 *max_k = gmres->max_k;
630 PetscFunctionReturn(0);
631 }
632
KSPGMRESSetRestart_GMRES(KSP ksp,PetscInt max_k)633 PetscErrorCode KSPGMRESSetRestart_GMRES(KSP ksp,PetscInt max_k)
634 {
635 KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
636 PetscErrorCode ierr;
637
638 PetscFunctionBegin;
639 if (max_k < 1) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Restart must be positive");
640 if (!ksp->setupstage) {
641 gmres->max_k = max_k;
642 } else if (gmres->max_k != max_k) {
643 gmres->max_k = max_k;
644 ksp->setupstage = KSP_SETUP_NEW;
645 /* free the data structures, then create them again */
646 ierr = KSPReset_GMRES(ksp);CHKERRQ(ierr);
647 }
648 PetscFunctionReturn(0);
649 }
650
KSPGMRESSetOrthogonalization_GMRES(KSP ksp,FCN fcn)651 PetscErrorCode KSPGMRESSetOrthogonalization_GMRES(KSP ksp,FCN fcn)
652 {
653 PetscFunctionBegin;
654 ((KSP_GMRES*)ksp->data)->orthog = fcn;
655 PetscFunctionReturn(0);
656 }
657
KSPGMRESGetOrthogonalization_GMRES(KSP ksp,FCN * fcn)658 PetscErrorCode KSPGMRESGetOrthogonalization_GMRES(KSP ksp,FCN *fcn)
659 {
660 PetscFunctionBegin;
661 *fcn = ((KSP_GMRES*)ksp->data)->orthog;
662 PetscFunctionReturn(0);
663 }
664
KSPGMRESSetPreAllocateVectors_GMRES(KSP ksp)665 PetscErrorCode KSPGMRESSetPreAllocateVectors_GMRES(KSP ksp)
666 {
667 KSP_GMRES *gmres;
668
669 PetscFunctionBegin;
670 gmres = (KSP_GMRES*)ksp->data;
671 gmres->q_preallocate = 1;
672 PetscFunctionReturn(0);
673 }
674
KSPGMRESSetCGSRefinementType_GMRES(KSP ksp,KSPGMRESCGSRefinementType type)675 PetscErrorCode KSPGMRESSetCGSRefinementType_GMRES(KSP ksp,KSPGMRESCGSRefinementType type)
676 {
677 KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
678
679 PetscFunctionBegin;
680 gmres->cgstype = type;
681 PetscFunctionReturn(0);
682 }
683
KSPGMRESGetCGSRefinementType_GMRES(KSP ksp,KSPGMRESCGSRefinementType * type)684 PetscErrorCode KSPGMRESGetCGSRefinementType_GMRES(KSP ksp,KSPGMRESCGSRefinementType *type)
685 {
686 KSP_GMRES *gmres = (KSP_GMRES*)ksp->data;
687
688 PetscFunctionBegin;
689 *type = gmres->cgstype;
690 PetscFunctionReturn(0);
691 }
692
693 /*@
694 KSPGMRESSetCGSRefinementType - Sets the type of iterative refinement to use
695 in the classical Gram Schmidt orthogonalization.
696
697 Logically Collective on ksp
698
699 Input Parameters:
700 + ksp - the Krylov space context
701 - type - the type of refinement
702
703 Options Database:
704 . -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always>
705
706 Level: intermediate
707
708 .seealso: KSPGMRESSetOrthogonalization(), KSPGMRESCGSRefinementType, KSPGMRESClassicalGramSchmidtOrthogonalization(), KSPGMRESGetCGSRefinementType(),
709 KSPGMRESGetOrthogonalization()
710 @*/
KSPGMRESSetCGSRefinementType(KSP ksp,KSPGMRESCGSRefinementType type)711 PetscErrorCode KSPGMRESSetCGSRefinementType(KSP ksp,KSPGMRESCGSRefinementType type)
712 {
713 PetscErrorCode ierr;
714
715 PetscFunctionBegin;
716 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
717 PetscValidLogicalCollectiveEnum(ksp,type,2);
718 ierr = PetscTryMethod(ksp,"KSPGMRESSetCGSRefinementType_C",(KSP,KSPGMRESCGSRefinementType),(ksp,type));CHKERRQ(ierr);
719 PetscFunctionReturn(0);
720 }
721
722 /*@
723 KSPGMRESGetCGSRefinementType - Gets the type of iterative refinement to use
724 in the classical Gram Schmidt orthogonalization.
725
726 Not Collective
727
728 Input Parameter:
729 . ksp - the Krylov space context
730
731 Output Parameter:
732 . type - the type of refinement
733
734 Options Database:
735 . -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always>
736
737 Level: intermediate
738
739 .seealso: KSPGMRESSetOrthogonalization(), KSPGMRESCGSRefinementType, KSPGMRESClassicalGramSchmidtOrthogonalization(), KSPGMRESSetCGSRefinementType(),
740 KSPGMRESGetOrthogonalization()
741 @*/
KSPGMRESGetCGSRefinementType(KSP ksp,KSPGMRESCGSRefinementType * type)742 PetscErrorCode KSPGMRESGetCGSRefinementType(KSP ksp,KSPGMRESCGSRefinementType *type)
743 {
744 PetscErrorCode ierr;
745
746 PetscFunctionBegin;
747 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
748 ierr = PetscUseMethod(ksp,"KSPGMRESGetCGSRefinementType_C",(KSP,KSPGMRESCGSRefinementType*),(ksp,type));CHKERRQ(ierr);
749 PetscFunctionReturn(0);
750 }
751
752
753 /*@
754 KSPGMRESSetRestart - Sets number of iterations at which GMRES, FGMRES and LGMRES restarts.
755
756 Logically Collective on ksp
757
758 Input Parameters:
759 + ksp - the Krylov space context
760 - restart - integer restart value
761
762 Options Database:
763 . -ksp_gmres_restart <positive integer>
764
765 Note: The default value is 30.
766
767 Level: intermediate
768
769 .seealso: KSPSetTolerances(), KSPGMRESSetOrthogonalization(), KSPGMRESSetPreAllocateVectors(), KSPGMRESGetRestart()
770 @*/
KSPGMRESSetRestart(KSP ksp,PetscInt restart)771 PetscErrorCode KSPGMRESSetRestart(KSP ksp, PetscInt restart)
772 {
773 PetscErrorCode ierr;
774
775 PetscFunctionBegin;
776 PetscValidLogicalCollectiveInt(ksp,restart,2);
777
778 ierr = PetscTryMethod(ksp,"KSPGMRESSetRestart_C",(KSP,PetscInt),(ksp,restart));CHKERRQ(ierr);
779 PetscFunctionReturn(0);
780 }
781
782 /*@
783 KSPGMRESGetRestart - Gets number of iterations at which GMRES, FGMRES and LGMRES restarts.
784
785 Not Collective
786
787 Input Parameter:
788 . ksp - the Krylov space context
789
790 Output Parameter:
791 . restart - integer restart value
792
793 Note: The default value is 30.
794
795 Level: intermediate
796
797 .seealso: KSPSetTolerances(), KSPGMRESSetOrthogonalization(), KSPGMRESSetPreAllocateVectors(), KSPGMRESSetRestart()
798 @*/
KSPGMRESGetRestart(KSP ksp,PetscInt * restart)799 PetscErrorCode KSPGMRESGetRestart(KSP ksp, PetscInt *restart)
800 {
801 PetscErrorCode ierr;
802
803 PetscFunctionBegin;
804 ierr = PetscUseMethod(ksp,"KSPGMRESGetRestart_C",(KSP,PetscInt*),(ksp,restart));CHKERRQ(ierr);
805 PetscFunctionReturn(0);
806 }
807
808 /*@
809 KSPGMRESSetHapTol - Sets tolerance for determining happy breakdown in GMRES, FGMRES and LGMRES.
810
811 Logically Collective on ksp
812
813 Input Parameters:
814 + ksp - the Krylov space context
815 - tol - the tolerance
816
817 Options Database:
818 . -ksp_gmres_haptol <positive real value>
819
820 Note: Happy breakdown is the rare case in GMRES where an 'exact' solution is obtained after
821 a certain number of iterations. If you attempt more iterations after this point unstable
822 things can happen hence very occasionally you may need to set this value to detect this condition
823
824 Level: intermediate
825
826 .seealso: KSPSetTolerances()
827 @*/
KSPGMRESSetHapTol(KSP ksp,PetscReal tol)828 PetscErrorCode KSPGMRESSetHapTol(KSP ksp,PetscReal tol)
829 {
830 PetscErrorCode ierr;
831
832 PetscFunctionBegin;
833 PetscValidLogicalCollectiveReal(ksp,tol,2);
834 ierr = PetscTryMethod((ksp),"KSPGMRESSetHapTol_C",(KSP,PetscReal),((ksp),(tol)));CHKERRQ(ierr);
835 PetscFunctionReturn(0);
836 }
837
838 /*MC
839 KSPGMRES - Implements the Generalized Minimal Residual method.
840 (Saad and Schultz, 1986) with restart
841
842
843 Options Database Keys:
844 + -ksp_gmres_restart <restart> - the number of Krylov directions to orthogonalize against
845 . -ksp_gmres_haptol <tol> - sets the tolerance for "happy ending" (exact convergence)
846 . -ksp_gmres_preallocate - preallocate all the Krylov search directions initially (otherwise groups of
847 vectors are allocated as needed)
848 . -ksp_gmres_classicalgramschmidt - use classical (unmodified) Gram-Schmidt to orthogonalize against the Krylov space (fast) (the default)
849 . -ksp_gmres_modifiedgramschmidt - use modified Gram-Schmidt in the orthogonalization (more stable, but slower)
850 . -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always> - determine if iterative refinement is used to increase the
851 stability of the classical Gram-Schmidt orthogonalization.
852 - -ksp_gmres_krylov_monitor - plot the Krylov space generated
853
854 Level: beginner
855
856 Notes:
857 Left and right preconditioning are supported, but not symmetric preconditioning.
858
859 References:
860 . 1. - YOUCEF SAAD AND MARTIN H. SCHULTZ, GMRES: A GENERALIZED MINIMAL RESIDUAL ALGORITHM FOR SOLVING NONSYMMETRIC LINEAR SYSTEMS.
861 SIAM J. ScI. STAT. COMPUT. Vo|. 7, No. 3, July 1986.
862
863 .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPFGMRES, KSPLGMRES,
864 KSPGMRESSetRestart(), KSPGMRESSetHapTol(), KSPGMRESSetPreAllocateVectors(), KSPGMRESSetOrthogonalization(), KSPGMRESGetOrthogonalization(),
865 KSPGMRESClassicalGramSchmidtOrthogonalization(), KSPGMRESModifiedGramSchmidtOrthogonalization(),
866 KSPGMRESCGSRefinementType, KSPGMRESSetCGSRefinementType(), KSPGMRESGetCGSRefinementType(), KSPGMRESMonitorKrylov(), KSPSetPCSide()
867
868 M*/
869
KSPCreate_GMRES(KSP ksp)870 PETSC_EXTERN PetscErrorCode KSPCreate_GMRES(KSP ksp)
871 {
872 KSP_GMRES *gmres;
873 PetscErrorCode ierr;
874
875 PetscFunctionBegin;
876 ierr = PetscNewLog(ksp,&gmres);CHKERRQ(ierr);
877 ksp->data = (void*)gmres;
878
879 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,4);CHKERRQ(ierr);
880 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,3);CHKERRQ(ierr);
881 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_SYMMETRIC,2);CHKERRQ(ierr);
882 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);CHKERRQ(ierr);
883 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);CHKERRQ(ierr);
884
885 ksp->ops->buildsolution = KSPBuildSolution_GMRES;
886 ksp->ops->setup = KSPSetUp_GMRES;
887 ksp->ops->solve = KSPSolve_GMRES;
888 ksp->ops->reset = KSPReset_GMRES;
889 ksp->ops->destroy = KSPDestroy_GMRES;
890 ksp->ops->view = KSPView_GMRES;
891 ksp->ops->setfromoptions = KSPSetFromOptions_GMRES;
892 ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_GMRES;
893 ksp->ops->computeeigenvalues = KSPComputeEigenvalues_GMRES;
894 #if !defined(PETSC_USE_COMPLEX) && !defined(PETSC_HAVE_ESSL)
895 ksp->ops->computeritz = KSPComputeRitz_GMRES;
896 #endif
897 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetPreAllocateVectors_C",KSPGMRESSetPreAllocateVectors_GMRES);CHKERRQ(ierr);
898 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetOrthogonalization_C",KSPGMRESSetOrthogonalization_GMRES);CHKERRQ(ierr);
899 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetOrthogonalization_C",KSPGMRESGetOrthogonalization_GMRES);CHKERRQ(ierr);
900 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetRestart_C",KSPGMRESSetRestart_GMRES);CHKERRQ(ierr);
901 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetRestart_C",KSPGMRESGetRestart_GMRES);CHKERRQ(ierr);
902 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetHapTol_C",KSPGMRESSetHapTol_GMRES);CHKERRQ(ierr);
903 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetCGSRefinementType_C",KSPGMRESSetCGSRefinementType_GMRES);CHKERRQ(ierr);
904 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetCGSRefinementType_C",KSPGMRESGetCGSRefinementType_GMRES);CHKERRQ(ierr);
905
906 gmres->haptol = 1.0e-30;
907 gmres->q_preallocate = 0;
908 gmres->delta_allocate = GMRES_DELTA_DIRECTIONS;
909 gmres->orthog = KSPGMRESClassicalGramSchmidtOrthogonalization;
910 gmres->nrs = NULL;
911 gmres->sol_temp = NULL;
912 gmres->max_k = GMRES_DEFAULT_MAXK;
913 gmres->Rsvd = NULL;
914 gmres->cgstype = KSP_GMRES_CGS_REFINE_NEVER;
915 gmres->orthogwork = NULL;
916 PetscFunctionReturn(0);
917 }
918
919