xref: /petsc/src/snes/impls/ngmres/snesngmres.c (revision fbada4e6d62d82be69a0720ca41ff3cebed5b071)
1 /* Defines the basic SNES object */
2 #include <../src/snes/impls/ngmres/snesngmres.h>
3 #include <petscblaslapack.h>
4 
5 
6 
7 
8 #undef __FUNCT__
9 #define __FUNCT__ "SNESReset_NGMRES"
10 PetscErrorCode SNESReset_NGMRES(SNES snes)
11 {
12   SNES_NGMRES   *ngmres = (SNES_NGMRES*) snes->data;
13   PetscErrorCode ierr;
14 
15   PetscFunctionBegin;
16   ierr = VecDestroyVecs(ngmres->msize, &ngmres->Fdot);CHKERRQ(ierr);
17   ierr = VecDestroyVecs(ngmres->msize, &ngmres->Xdot);CHKERRQ(ierr);
18   if (snes->work) {ierr = VecDestroyVecs(snes->nwork, &snes->work);CHKERRQ(ierr);}
19   PetscFunctionReturn(0);
20 }
21 
22 #undef __FUNCT__
23 #define __FUNCT__ "SNESDestroy_NGMRES"
24 PetscErrorCode SNESDestroy_NGMRES(SNES snes)
25 {
26   PetscErrorCode ierr;
27 
28   PetscFunctionBegin;
29   ierr = SNESReset_NGMRES(snes);CHKERRQ(ierr);
30   if (snes->data) {
31     SNES_NGMRES * ngmres = (SNES_NGMRES *)snes->data;
32     ierr = PetscFree5(ngmres->h, ngmres->beta, ngmres->xi, ngmres->fnorms, ngmres->q);CHKERRQ(ierr);
33     ierr = PetscFree(ngmres->s);CHKERRQ(ierr);
34 #if PETSC_USE_COMPLEX
35     ierr = PetscFree(ngmres->rwork);
36 #endif
37     ierr = PetscFree(ngmres->work);
38   }
39   ierr = PetscFree(snes->data);
40   PetscFunctionReturn(0);
41 }
42 
43 #undef __FUNCT__
44 #define __FUNCT__ "SNESSetUp_NGMRES"
45 PetscErrorCode SNESSetUp_NGMRES(SNES snes)
46 {
47   SNES_NGMRES    *ngmres = (SNES_NGMRES *) snes->data;
48   PetscInt       msize,hsize;
49   PetscErrorCode ierr;
50 
51   PetscFunctionBegin;
52   msize         = ngmres->msize;  /* restart size */
53   hsize         = msize * msize;
54 
55 
56   /* explicit least squares minimization solve */
57   ierr = PetscMalloc5(hsize,PetscScalar,&ngmres->h,
58                       msize,PetscScalar,&ngmres->beta,
59                       msize,PetscScalar,&ngmres->xi,
60                       msize,PetscReal,  &ngmres->fnorms,
61                       hsize,PetscScalar,&ngmres->q);CHKERRQ(ierr);
62   ngmres->nrhs = 1;
63   ngmres->lda = msize;
64   ngmres->ldb = msize;
65   ierr = PetscMalloc(msize*sizeof(PetscScalar),&ngmres->s);CHKERRQ(ierr);
66   ierr = PetscMemzero(ngmres->h,   hsize*sizeof(PetscScalar));CHKERRQ(ierr);
67   ierr = PetscMemzero(ngmres->q,   hsize*sizeof(PetscScalar));CHKERRQ(ierr);
68   ierr = PetscMemzero(ngmres->xi,  msize*sizeof(PetscScalar));CHKERRQ(ierr);
69   ierr = PetscMemzero(ngmres->beta,msize*sizeof(PetscScalar));CHKERRQ(ierr);
70   ngmres->lwork = 12*msize;
71 #if PETSC_USE_COMPLEX
72   ierr = PetscMalloc(sizeof(PetscReal)*ngmres->lwork,&ngmres->rwork);
73 #endif
74   ierr = PetscMalloc(sizeof(PetscScalar)*ngmres->lwork,&ngmres->work);
75 
76   ierr = VecDuplicateVecs(snes->vec_sol, ngmres->msize, &ngmres->Xdot);CHKERRQ(ierr);
77   ierr = VecDuplicateVecs(snes->vec_sol, ngmres->msize, &ngmres->Fdot);CHKERRQ(ierr);
78   ierr = SNESDefaultGetWork(snes, 5);CHKERRQ(ierr);
79   PetscFunctionReturn(0);
80 }
81 
82 #undef __FUNCT__
83 #define __FUNCT__ "SNESSetFromOptions_NGMRES"
84 PetscErrorCode SNESSetFromOptions_NGMRES(SNES snes)
85 {
86   SNES_NGMRES   *ngmres = (SNES_NGMRES *) snes->data;
87   PetscErrorCode ierr;
88   PetscBool      debug;
89 
90   PetscFunctionBegin;
91   ierr = PetscOptionsHead("SNES NGMRES options");CHKERRQ(ierr);
92   ierr = PetscOptionsInt("-snes_ngmres_m",        "Number of directions",               "SNES", ngmres->msize,  &ngmres->msize, PETSC_NULL);CHKERRQ(ierr);
93   ierr = PetscOptionsInt("-snes_ngmres_restart",  "Maximum iterations before restart.", "SNES", ngmres->k_rmax, &ngmres->k_rmax, PETSC_NULL);CHKERRQ(ierr);
94   ierr = PetscOptionsBool("-snes_ngmres_monitor", "Monitor actions of NGMRES",          "SNES", ngmres->monitor ? PETSC_TRUE: PETSC_FALSE, &debug, PETSC_NULL);CHKERRQ(ierr);
95   if (debug) {
96     ngmres->monitor = PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm);CHKERRQ(ierr);
97   }
98   ierr = PetscOptionsReal("-snes_ngmres_gammaA",   "Residual selection constant",   "SNES", ngmres->gammaA, &ngmres->gammaA, PETSC_NULL);CHKERRQ(ierr);
99   ierr = PetscOptionsReal("-snes_ngmres_gammaC", "  Residual restart constant",     "SNES", ngmres->gammaC, &ngmres->gammaC, PETSC_NULL);CHKERRQ(ierr);
100   ierr = PetscOptionsReal("-snes_ngmres_epsilonB", "Difference selection constant", "SNES", ngmres->epsilonB, &ngmres->epsilonB, PETSC_NULL);CHKERRQ(ierr);
101   ierr = PetscOptionsReal("-snes_ngmres_deltaB",   "Difference residual selection constant", "SNES", ngmres->deltaB, &ngmres->deltaB, PETSC_NULL);CHKERRQ(ierr);
102   ierr = PetscOptionsTail();CHKERRQ(ierr);
103   if ((ngmres->gammaA > ngmres->gammaC) && (ngmres->gammaC > 2.)) ngmres->gammaC = ngmres->gammaA;
104   PetscFunctionReturn(0);
105 }
106 
107 #undef __FUNCT__
108 #define __FUNCT__ "SNESView_NGMRES"
109 PetscErrorCode SNESView_NGMRES(SNES snes, PetscViewer viewer)
110 {
111   SNES_NGMRES   *ngmres = (SNES_NGMRES *) snes->data;
112   PetscBool      iascii;
113   const char     *cstr;
114   PetscErrorCode ierr;
115 
116   PetscFunctionBegin;
117   ierr = PetscTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
118   if (iascii) {
119     cstr = SNESLineSearchTypeName(snes->ls_type);
120     ierr = PetscViewerASCIIPrintf(viewer,"  line search variant: %s\n",cstr);CHKERRQ(ierr);
121     ierr = PetscViewerASCIIPrintf(viewer, "  Number of stored past updates: %d\n", ngmres->msize);CHKERRQ(ierr);
122     ierr = PetscViewerASCIIPrintf(viewer, "  Maximum iterations before total restart: %d\n", ngmres->k_rmax);CHKERRQ(ierr);
123     ierr = PetscViewerASCIIPrintf(viewer, "  Residual selection: gammaA=%1.0e, gammaC=%1.0e\n", ngmres->gammaA, ngmres->gammaC);CHKERRQ(ierr);
124     ierr = PetscViewerASCIIPrintf(viewer, "  Difference restart: epsilonB=%1.0e, deltaB=%1.0e\n", ngmres->epsilonB, ngmres->deltaB);CHKERRQ(ierr);
125   }
126   PetscFunctionReturn(0);
127 }
128 
129 EXTERN_C_BEGIN
130 #undef __FUNCT__
131 #define __FUNCT__ "SNESLineSearchSetType_NGMRES"
132 PetscErrorCode  SNESLineSearchSetType_NGMRES(SNES snes, SNESLineSearchType type)
133 {
134   PetscErrorCode ierr;
135   PetscFunctionBegin;
136 
137   switch (type) {
138   case SNES_LS_BASIC:
139     ierr = SNESLineSearchSet(snes,SNESLineSearchNo,PETSC_NULL);CHKERRQ(ierr);
140     break;
141   case SNES_LS_BASIC_NONORMS:
142     ierr = SNESLineSearchSet(snes,SNESLineSearchNoNorms,PETSC_NULL);CHKERRQ(ierr);
143     break;
144   case SNES_LS_QUADRATIC:
145     ierr = SNESLineSearchSet(snes,SNESLineSearchQuadraticSecant,PETSC_NULL);CHKERRQ(ierr);
146     break;
147   default:
148     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP,"Unknown line search type.");
149     break;
150   }
151   snes->ls_type = type;
152   PetscFunctionReturn(0);
153 }
154 EXTERN_C_END
155 
156 
157 #undef __FUNCT__
158 #define __FUNCT__ "SNESSolve_NGMRES"
159 
160 PetscErrorCode SNESSolve_NGMRES(SNES snes)
161 {
162   SNES_NGMRES        *ngmres = (SNES_NGMRES *) snes->data;
163   /* present solution, residual, and preconditioned residual */
164   Vec                 X, F, B, D, G, W, Y;
165 
166   /* candidate linear combination answers */
167   Vec                 XA, FA;
168 
169   /* previous iterations to construct the subspace */
170   Vec                 *Fdot = ngmres->Fdot;
171   Vec                 *Xdot = ngmres->Xdot;
172 
173   /* coefficients and RHS to the minimization problem */
174   PetscScalar         *beta = ngmres->beta;
175   PetscScalar         *xi   = ngmres->xi;
176   PetscReal           fnorm, fAnorm, gnorm, ynorm, xnorm = 0.0;
177   PetscReal           nu;
178   PetscScalar         alph_total = 0.;
179   PetscScalar         qentry;
180   PetscInt            i, j, k, k_restart, l, ivec;
181 
182   /* solution selection data */
183   PetscBool           selectA, selectRestart;
184   PetscReal           dnorm, dminnorm, dcurnorm;
185   PetscReal           fminnorm;
186 
187   SNESConvergedReason reason;
188   PetscBool           lssucceed;
189   PetscErrorCode      ierr;
190 
191   PetscFunctionBegin;
192   /* variable initialization */
193   snes->reason  = SNES_CONVERGED_ITERATING;
194   X             = snes->vec_sol;
195   F             = snes->vec_func;
196   B             = snes->vec_rhs;
197   XA            = snes->vec_sol_update;
198   FA            = snes->work[0];
199   D             = snes->work[1];
200 
201   /* work for the line search */
202   Y             = snes->work[2];
203   G             = snes->work[3];
204   W             = snes->work[4];
205 
206   ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr);
207   snes->iter = 0;
208   snes->norm = 0.;
209   ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr);
210 
211   /* initialization */
212 
213   /* r = F(x) */
214   ierr = SNESComputeFunction(snes, X, F);CHKERRQ(ierr);
215   if (snes->domainerror) {
216     snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
217     PetscFunctionReturn(0);
218   }
219 
220   /* nu = (r, r) */
221   ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr);
222   fminnorm = fnorm;
223   nu = fnorm*fnorm;
224   if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP, "Infinite or not-a-number generated in function evaluation");
225 
226   /* q_{00} = nu  */
227   Q(0,0) = nu;
228   ngmres->fnorms[0] = fnorm;
229   /* Fdot[0] = F */
230   ierr = VecCopy(X, Xdot[0]);CHKERRQ(ierr);
231   ierr = VecCopy(F, Fdot[0]);CHKERRQ(ierr);
232 
233   ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr);
234   snes->norm = fnorm;
235   ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr);
236   SNESLogConvHistory(snes, fnorm, 0);
237   ierr = SNESMonitor(snes, 0, fnorm);CHKERRQ(ierr);
238   ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
239   if (snes->reason) PetscFunctionReturn(0);
240 
241   k_restart = 1;
242   l = 1;
243   for (k=1; k<snes->max_its; k++) {
244     /* select which vector of the stored subspace will be updated */
245     ivec = k_restart % ngmres->msize; /* replace the last used part of the subspace */
246 
247     /* Computation of x^M */
248     if (!snes->pc) {
249       /* no preconditioner -- just take gradient descent with line search */
250       ierr = VecCopy(F, Y);CHKERRQ(ierr);
251       ierr = VecScale(Y, -1.0);CHKERRQ(ierr);
252       ierr = (*snes->ops->linesearch)(snes,snes->lsP,X,F,Y,fnorm,xnorm,G,W,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr);
253       if (!lssucceed) {
254         if (++snes->numFailures >= snes->maxFailures) {
255           snes->reason = SNES_DIVERGED_LINE_SEARCH;
256           PetscFunctionReturn(0);
257         }
258       }
259       fnorm = gnorm;
260       ierr = VecCopy(G, F);CHKERRQ(ierr);
261       ierr = VecCopy(W, X);CHKERRQ(ierr);
262     } else {
263       ierr = SNESSolve(snes->pc, B, X);CHKERRQ(ierr);
264       ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr);
265       if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
266         snes->reason = SNES_DIVERGED_INNER;
267         PetscFunctionReturn(0);
268       }
269       ierr = SNESComputeFunction(snes, X, F);CHKERRQ(ierr);
270       ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr);
271     }
272 
273     /* r = F(x) */
274     nu = fnorm*fnorm;
275     if (fminnorm > fnorm) fminnorm = fnorm;  /* the minimum norm is now of F^M */
276 
277     /* construct the right hand side and xi factors */
278     for (i = 0; i < l; i++) {
279       ierr = VecDot(Fdot[i], F, &xi[i]);CHKERRQ(ierr);
280       beta[i] = nu - xi[i];
281     }
282 
283     /* construct h */
284     for (j = 0; j < l; j++) {
285       for (i = 0; i < l; i++) {
286         H(i, j) = Q(i, j) - xi[i] - xi[j] + nu;
287       }
288     }
289 
290     if(l == 1) {
291       /* simply set alpha[0] = beta[0] / H[0, 0] */
292       beta[0] = beta[0] / H(0, 0);
293     } else {
294 #ifdef PETSC_MISSING_LAPACK_GELSS
295     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "NGMRES with LS requires the LAPACK GELSS routine.");
296 #else
297     ngmres->m = PetscBLASIntCast(l);
298     ngmres->n = PetscBLASIntCast(l);
299     ngmres->info = PetscBLASIntCast(0);
300     ngmres->rcond = -1.;
301 #ifdef PETSC_USE_COMPLEX
302     LAPACKgelss_(&ngmres->m,
303                  &ngmres->n,
304                  &ngmres->nrhs,
305                  ngmres->h,
306                  &ngmres->lda,
307                  ngmres->beta,
308                  &ngmres->ldb,
309                  ngmres->s,
310                  &ngmres->rcond,
311                  &ngmres->rank,
312                  ngmres->work,
313                  &ngmres->lwork,
314                  ngmres->rwork,
315                  &ngmres->info);
316 #else
317     LAPACKgelss_(&ngmres->m,
318                  &ngmres->n,
319                  &ngmres->nrhs,
320                  ngmres->h,
321                  &ngmres->lda,
322                  ngmres->beta,
323                  &ngmres->ldb,
324                  ngmres->s,
325                  &ngmres->rcond,
326                  &ngmres->rank,
327                  ngmres->work,
328                  &ngmres->lwork,
329                  &ngmres->info);
330 #endif
331     if (ngmres->info < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELSS");
332     if (ngmres->info > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SVD failed to converge");
333 #endif
334     }
335 
336     alph_total = 0.;
337     for (i = 0; i < l; i++) {
338       alph_total += beta[i];
339     }
340 
341     ierr = VecCopy(X, XA);CHKERRQ(ierr);
342     ierr = VecScale(XA, 1. - alph_total);CHKERRQ(ierr);
343 
344     for(i=0;i<l;i++){
345       ierr= VecAXPY(XA, beta[i], Xdot[i]);CHKERRQ(ierr);
346     }
347     ierr = SNESComputeFunction(snes, XA, FA);CHKERRQ(ierr);
348     ierr = VecNorm(FA, NORM_2, &fAnorm);CHKERRQ(ierr);
349 
350     selectA = PETSC_TRUE;
351     /* Conditions for choosing the accelerated answer */
352 
353     /* Criterion A -- the norm of the function isn't increased above the minimum by too much */
354     if (fAnorm >= ngmres->gammaA*fminnorm) {
355       selectA = PETSC_FALSE;
356     }
357 
358     /* Criterion B -- the choice of x^A isn't too close to some other choice */
359     ierr=VecCopy(XA, D);CHKERRQ(ierr);
360     ierr=VecAXPY(D, -1.0, X);CHKERRQ(ierr);
361     ierr=VecNorm(D, NORM_2, &dnorm);CHKERRQ(ierr);
362     dminnorm = -1.0;
363     for(i=0;i<l;i++) {
364       ierr=VecCopy(XA, D);CHKERRQ(ierr);
365       ierr=VecAXPY(D, -1.0, Xdot[i]);CHKERRQ(ierr);
366       ierr=VecNorm(D, NORM_2, &dcurnorm);CHKERRQ(ierr);
367       if((dcurnorm < dminnorm) || (dminnorm < 0.0)) dminnorm = dcurnorm;
368     }
369     if (ngmres->epsilonB*dnorm<dminnorm || sqrt(fnorm)<ngmres->deltaB*sqrt(fminnorm)) {
370     } else {
371       selectA=PETSC_FALSE;
372     }
373 
374 
375     if (selectA) {
376       if (ngmres->monitor) {
377         ierr = PetscViewerASCIIPrintf(ngmres->monitor, "picked X_A, ||F_A||_2 = %e, ||F_M||_2 = %e\n", fAnorm, fnorm);CHKERRQ(ierr);
378       }
379       /* copy it over */
380       fnorm = fAnorm;
381       nu = fnorm*fnorm;
382       ierr = VecCopy(FA, F);CHKERRQ(ierr);
383       ierr = VecCopy(XA, X);CHKERRQ(ierr);
384     } else {
385       if (ngmres->monitor) {
386         ierr = PetscViewerASCIIPrintf(ngmres->monitor, "picked X_M, ||F_A||_2 = %e, ||F_M||_2 = %e\n", fAnorm, fnorm);CHKERRQ(ierr);
387       }
388     }
389 
390     selectRestart = PETSC_FALSE;
391 
392     /* maximum iteration criterion */
393     if (k_restart > ngmres->k_rmax) {
394       selectRestart = PETSC_TRUE;
395     }
396 
397     /* difference stagnation restart */
398     if((ngmres->epsilonB*dnorm > dminnorm) && (sqrt(fAnorm) > ngmres->deltaB*sqrt(fminnorm))) {
399       if (ngmres->monitor) {
400         ierr = PetscViewerASCIIPrintf(ngmres->monitor, "difference restart: %e > %e\n", ngmres->epsilonB*dnorm, dminnorm);CHKERRQ(ierr);
401       }
402       selectRestart = PETSC_TRUE;
403     }
404     /* residual stagnation restart */
405     if (sqrt(fAnorm) > ngmres->gammaC*sqrt(fminnorm)) {
406       if (ngmres->monitor) {
407         ierr = PetscViewerASCIIPrintf(ngmres->monitor, "residual restart: %e > %e\n", sqrt(fAnorm), ngmres->gammaC*sqrt(fminnorm));CHKERRQ(ierr);
408       }
409       selectRestart = PETSC_TRUE;
410     }
411 
412     if (selectRestart) {
413       if (ngmres->monitor){
414         ierr = PetscViewerASCIIPrintf(ngmres->monitor, "Restarted at iteration %d\n", k_restart);CHKERRQ(ierr);
415       }
416       k_restart = 1;
417       l = 1;
418       /* q_{00} = nu */
419       ngmres->fnorms[0] = fnorm;
420       nu = fnorm*fnorm;
421       Q(0,0) = nu;
422       /* Fdot[0] = F */
423       ierr = VecCopy(X, Xdot[0]);CHKERRQ(ierr);
424       ierr = VecCopy(F, Fdot[0]);CHKERRQ(ierr);
425     } else {
426       /* select the current size of the subspace */
427       if (l < ngmres->msize) l++;
428       k_restart++;
429       /* place the current entry in the list of previous entries */
430       ierr = VecCopy(F, Fdot[ivec]);CHKERRQ(ierr);
431       ierr = VecCopy(X, Xdot[ivec]);CHKERRQ(ierr);
432       ngmres->fnorms[ivec] = fnorm;
433       if (fminnorm > fnorm) fminnorm = fnorm;  /* the minimum norm is now of r^A */
434       for (i = 0; i < l; i++) {
435         ierr = VecDot(F, Fdot[i], &qentry);CHKERRQ(ierr);
436         Q(i, ivec) = qentry;
437         Q(ivec, i) = qentry;
438       }
439     }
440 
441     ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr);
442     snes->iter = k;
443     snes->norm = fnorm;
444     ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr);
445     SNESLogConvHistory(snes, snes->norm, snes->iter);
446     ierr = SNESMonitor(snes, snes->iter, snes->norm);CHKERRQ(ierr);
447 
448     ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
449     if (snes->reason) PetscFunctionReturn(0);
450   }
451   snes->reason = SNES_DIVERGED_MAX_IT;
452   PetscFunctionReturn(0);
453 }
454 
455 /*MC
456   SNESNGMRES - The Nonlinear Generalized Minimum Residual (NGMRES) method of Oosterlee and Washio.
457 
458    Level: beginner
459 
460    "Krylov Subspace Acceleration of Nonlinear Multigrid with Application to Recirculating Flows", C. W. Oosterlee and T. Washio,
461    SIAM Journal on Scientific Computing, 21(5), 2000.
462 
463    This is also the same as the algorithm called Anderson acceleration introduced in "D. G. Anderson. Iterative procedures for nonlinear integral equations.
464    J. Assoc. Comput. Mach., 12:547–560, 1965."
465 
466 .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types)
467 M*/
468 
469 EXTERN_C_BEGIN
470 #undef __FUNCT__
471 #define __FUNCT__ "SNESCreate_NGMRES"
472 PetscErrorCode SNESCreate_NGMRES(SNES snes)
473 {
474   SNES_NGMRES   *ngmres;
475   PetscErrorCode ierr;
476 
477   PetscFunctionBegin;
478   snes->ops->destroy        = SNESDestroy_NGMRES;
479   snes->ops->setup          = SNESSetUp_NGMRES;
480   snes->ops->setfromoptions = SNESSetFromOptions_NGMRES;
481   snes->ops->view           = SNESView_NGMRES;
482   snes->ops->solve          = SNESSolve_NGMRES;
483   snes->ops->reset          = SNESReset_NGMRES;
484 
485   snes->usespc          = PETSC_TRUE;
486   snes->usesksp         = PETSC_FALSE;
487 
488   ierr = PetscNewLog(snes, SNES_NGMRES, &ngmres);CHKERRQ(ierr);
489   snes->data = (void*) ngmres;
490   ngmres->msize = 10;
491 
492   ngmres->gammaA   = 2.0;
493   ngmres->gammaC   = 2.0;
494   ngmres->deltaB   = 0.9;
495   ngmres->epsilonB = 0.1;
496   ngmres->k_rmax   = 200;
497   ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetType_C","SNESLineSearchSetType_NGMRES",SNESLineSearchSetType_NGMRES);CHKERRQ(ierr);
498   ierr = SNESLineSearchSetType(snes, SNES_LS_QUADRATIC);CHKERRQ(ierr);
499   PetscFunctionReturn(0);
500 }
501 EXTERN_C_END
502