xref: /petsc/src/snes/impls/ngmres/snesngmres.c (revision 794bee337f0d5af71ccfe93f8f7f88922f7fb150)
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_it","Tolerance iterations before restart","SNES", ngmres->restart_it,  &ngmres->restart_it, 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, "  Residual selection: gammaA=%1.0e, gammaC=%1.0e\n", ngmres->gammaA, ngmres->gammaC);CHKERRQ(ierr);
123     ierr = PetscViewerASCIIPrintf(viewer, "  Difference restart: epsilonB=%1.0e, deltaB=%1.0e\n", ngmres->epsilonB, ngmres->deltaB);CHKERRQ(ierr);
124   }
125   PetscFunctionReturn(0);
126 }
127 
128 EXTERN_C_BEGIN
129 #undef __FUNCT__
130 #define __FUNCT__ "SNESLineSearchSetType_NGMRES"
131 PetscErrorCode  SNESLineSearchSetType_NGMRES(SNES snes, SNESLineSearchType type)
132 {
133   PetscErrorCode ierr;
134   PetscFunctionBegin;
135 
136   switch (type) {
137   case SNES_LS_BASIC:
138     ierr = SNESLineSearchSet(snes,SNESLineSearchNo,PETSC_NULL);CHKERRQ(ierr);
139     break;
140   case SNES_LS_BASIC_NONORMS:
141     ierr = SNESLineSearchSet(snes,SNESLineSearchNoNorms,PETSC_NULL);CHKERRQ(ierr);
142     break;
143   case SNES_LS_QUADRATIC:
144     ierr = SNESLineSearchSet(snes,SNESLineSearchQuadraticSecant,PETSC_NULL);CHKERRQ(ierr);
145     break;
146   default:
147     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP,"Unknown line search type.");
148     break;
149   }
150   snes->ls_type = type;
151   PetscFunctionReturn(0);
152 }
153 EXTERN_C_END
154 
155 
156 #undef __FUNCT__
157 #define __FUNCT__ "SNESSolve_NGMRES"
158 
159 PetscErrorCode SNESSolve_NGMRES(SNES snes)
160 {
161   SNES_NGMRES        *ngmres = (SNES_NGMRES *) snes->data;
162   /* present solution, residual, and preconditioned residual */
163   Vec                 X, F, B, D, G, W, Y;
164 
165   /* candidate linear combination answers */
166   Vec                 XA, FA;
167 
168   /* previous iterations to construct the subspace */
169   Vec                 *Fdot = ngmres->Fdot;
170   Vec                 *Xdot = ngmres->Xdot;
171 
172   /* coefficients and RHS to the minimization problem */
173   PetscScalar         *beta = ngmres->beta;
174   PetscScalar         *xi   = ngmres->xi;
175   PetscReal           fnorm, fAnorm, gnorm, ynorm, xnorm = 0.0;
176   PetscReal           nu;
177   PetscScalar         alph_total = 0.;
178   PetscScalar         qentry;
179   PetscInt            i, j, k, k_restart, l, ivec, restart_count = 0;
180 
181   /* solution selection data */
182   PetscBool           selectA, selectRestart;
183   PetscReal           dnorm, dminnorm, dcurnorm;
184   PetscReal           fminnorm;
185 
186   SNESConvergedReason reason;
187   PetscBool           lssucceed;
188   PetscErrorCode      ierr;
189 
190   PetscFunctionBegin;
191   /* variable initialization */
192   snes->reason  = SNES_CONVERGED_ITERATING;
193   X             = snes->vec_sol;
194   F             = snes->vec_func;
195   B             = snes->vec_rhs;
196   XA            = snes->vec_sol_update;
197   FA            = snes->work[0];
198   D             = snes->work[1];
199 
200   /* work for the line search */
201   Y             = snes->work[2];
202   G             = snes->work[3];
203   W             = snes->work[4];
204 
205   ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr);
206   snes->iter = 0;
207   snes->norm = 0.;
208   ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr);
209 
210   /* initialization */
211 
212   /* r = F(x) */
213   ierr = SNESComputeFunction(snes, X, F);CHKERRQ(ierr);
214   if (snes->domainerror) {
215     snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
216     PetscFunctionReturn(0);
217   }
218 
219   /* nu = (r, r) */
220   ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr);
221   fminnorm = fnorm;
222   nu = fnorm*fnorm;
223   if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP, "Infinite or not-a-number generated in function evaluation");
224 
225   /* q_{00} = nu  */
226   Q(0,0) = nu;
227   ngmres->fnorms[0] = fnorm;
228   /* Fdot[0] = F */
229   ierr = VecCopy(X, Xdot[0]);CHKERRQ(ierr);
230   ierr = VecCopy(F, Fdot[0]);CHKERRQ(ierr);
231 
232   ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr);
233   snes->norm = fnorm;
234   ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr);
235   SNESLogConvHistory(snes, fnorm, 0);
236   ierr = SNESMonitor(snes, 0, fnorm);CHKERRQ(ierr);
237   ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
238   if (snes->reason) PetscFunctionReturn(0);
239 
240   k_restart = 1;
241   l = 1;
242   for (k=1; k<snes->max_its; k++) {
243     /* select which vector of the stored subspace will be updated */
244     ivec = k_restart % ngmres->msize; /* replace the last used part of the subspace */
245 
246     /* Computation of x^M */
247     if (!snes->pc) {
248       /* no preconditioner -- just take gradient descent with line search */
249       ierr = VecCopy(F, Y);CHKERRQ(ierr);
250       ierr = VecScale(Y, -1.0);CHKERRQ(ierr);
251       ierr = (*snes->ops->linesearch)(snes,snes->lsP,X,F,Y,fnorm,xnorm,G,W,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr);
252       if (!lssucceed) {
253         if (++snes->numFailures >= snes->maxFailures) {
254           snes->reason = SNES_DIVERGED_LINE_SEARCH;
255           PetscFunctionReturn(0);
256         }
257       }
258       fnorm = gnorm;
259       ierr = VecCopy(G, F);CHKERRQ(ierr);
260       ierr = VecCopy(W, X);CHKERRQ(ierr);
261     } else {
262       ierr = SNESSolve(snes->pc, B, X);CHKERRQ(ierr);
263       ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr);
264       if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
265         snes->reason = SNES_DIVERGED_INNER;
266         PetscFunctionReturn(0);
267       }
268       ierr = SNESComputeFunction(snes, X, F);CHKERRQ(ierr);
269       ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr);
270     }
271 
272     /* r = F(x) */
273     nu = fnorm*fnorm;
274     if (fminnorm > fnorm) fminnorm = fnorm;  /* the minimum norm is now of F^M */
275 
276     /* construct the right hand side and xi factors */
277     for (i = 0; i < l; i++) {
278       ierr = VecDot(Fdot[i], F, &xi[i]);CHKERRQ(ierr);
279       beta[i] = nu - xi[i];
280     }
281 
282     /* construct h */
283     for (j = 0; j < l; j++) {
284       for (i = 0; i < l; i++) {
285         H(i, j) = Q(i, j) - xi[i] - xi[j] + nu;
286       }
287     }
288 
289     if(l == 1) {
290       /* simply set alpha[0] = beta[0] / H[0, 0] */
291       beta[0] = beta[0] / H(0, 0);
292     } else {
293 #ifdef PETSC_MISSING_LAPACK_GELSS
294     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "NGMRES with LS requires the LAPACK GELSS routine.");
295 #else
296     ngmres->m = PetscBLASIntCast(l);
297     ngmres->n = PetscBLASIntCast(l);
298     ngmres->info = PetscBLASIntCast(0);
299     ngmres->rcond = -1.;
300 #ifdef PETSC_USE_COMPLEX
301     LAPACKgelss_(&ngmres->m,
302                  &ngmres->n,
303                  &ngmres->nrhs,
304                  ngmres->h,
305                  &ngmres->lda,
306                  ngmres->beta,
307                  &ngmres->ldb,
308                  ngmres->s,
309                  &ngmres->rcond,
310                  &ngmres->rank,
311                  ngmres->work,
312                  &ngmres->lwork,
313                  ngmres->rwork,
314                  &ngmres->info);
315 #else
316     LAPACKgelss_(&ngmres->m,
317                  &ngmres->n,
318                  &ngmres->nrhs,
319                  ngmres->h,
320                  &ngmres->lda,
321                  ngmres->beta,
322                  &ngmres->ldb,
323                  ngmres->s,
324                  &ngmres->rcond,
325                  &ngmres->rank,
326                  ngmres->work,
327                  &ngmres->lwork,
328                  &ngmres->info);
329 #endif
330     if (ngmres->info < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELSS");
331     if (ngmres->info > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SVD failed to converge");
332 #endif
333     }
334 
335     alph_total = 0.;
336     for (i = 0; i < l; i++) {
337       alph_total += beta[i];
338     }
339 
340     ierr = VecCopy(X, XA);CHKERRQ(ierr);
341     ierr = VecScale(XA, 1. - alph_total);CHKERRQ(ierr);
342 
343     for(i=0;i<l;i++){
344       ierr= VecAXPY(XA, beta[i], Xdot[i]);CHKERRQ(ierr);
345     }
346     ierr = SNESComputeFunction(snes, XA, FA);CHKERRQ(ierr);
347     ierr = VecNorm(FA, NORM_2, &fAnorm);CHKERRQ(ierr);
348 
349     selectA = PETSC_TRUE;
350     /* Conditions for choosing the accelerated answer */
351 
352     /* Criterion A -- the norm of the function isn't increased above the minimum by too much */
353     if (fAnorm >= ngmres->gammaA*fminnorm) {
354       selectA = PETSC_FALSE;
355     }
356 
357     /* Criterion B -- the choice of x^A isn't too close to some other choice */
358     ierr=VecCopy(XA, D);CHKERRQ(ierr);
359     ierr=VecAXPY(D, -1.0, X);CHKERRQ(ierr);
360     ierr=VecNorm(D, NORM_2, &dnorm);CHKERRQ(ierr);
361     dminnorm = -1.0;
362     for(i=0;i<l;i++) {
363       ierr=VecCopy(XA, D);CHKERRQ(ierr);
364       ierr=VecAXPY(D, -1.0, Xdot[i]);CHKERRQ(ierr);
365       ierr=VecNorm(D, NORM_2, &dcurnorm);CHKERRQ(ierr);
366       if((dcurnorm < dminnorm) || (dminnorm < 0.0)) dminnorm = dcurnorm;
367     }
368     if (ngmres->epsilonB*dnorm<dminnorm || PetscSqrtReal(fnorm)<ngmres->deltaB*PetscSqrtReal(fminnorm)) {
369     } else {
370       selectA=PETSC_FALSE;
371     }
372 
373 
374     if (selectA) {
375       if (ngmres->monitor) {
376         ierr = PetscViewerASCIIPrintf(ngmres->monitor, "picked X_A, ||F_A||_2 = %e, ||F_M||_2 = %e\n", fAnorm, fnorm);CHKERRQ(ierr);
377       }
378       /* copy it over */
379       fnorm = fAnorm;
380       nu = fnorm*fnorm;
381       ierr = VecCopy(FA, F);CHKERRQ(ierr);
382       ierr = VecCopy(XA, X);CHKERRQ(ierr);
383     } else {
384       if (ngmres->monitor) {
385         ierr = PetscViewerASCIIPrintf(ngmres->monitor, "picked X_M, ||F_A||_2 = %e, ||F_M||_2 = %e\n", fAnorm, fnorm);CHKERRQ(ierr);
386       }
387     }
388 
389     selectRestart = PETSC_FALSE;
390 
391     /* difference stagnation restart */
392     if((ngmres->epsilonB*dnorm > dminnorm) && (PetscSqrtReal(fAnorm) > ngmres->deltaB*PetscSqrtReal(fminnorm))) {
393       if (ngmres->monitor) {
394         ierr = PetscViewerASCIIPrintf(ngmres->monitor, "difference restart: %e > %e\n", ngmres->epsilonB*dnorm, dminnorm);CHKERRQ(ierr);
395       }
396       selectRestart = PETSC_TRUE;
397     }
398     /* residual stagnation restart */
399     if (PetscSqrtReal(fAnorm) > ngmres->gammaC*PetscSqrtReal(fminnorm)) {
400       if (ngmres->monitor) {
401         ierr = PetscViewerASCIIPrintf(ngmres->monitor, "residual restart: %e > %e\n", PetscSqrtReal(fAnorm), ngmres->gammaC*PetscSqrtReal(fminnorm));CHKERRQ(ierr);
402       }
403       selectRestart = PETSC_TRUE;
404     }
405 
406     /* if the restart conditions persist for more than restart_it iterations, restart. */
407     if (selectRestart) {
408       restart_count++;
409     } else {
410       restart_count = 0;
411     }
412 
413     /* restart after restart conditions have persisted for a fixed number of iterations */
414     if (restart_count >= ngmres->restart_it) {
415       if (ngmres->monitor){
416         ierr = PetscViewerASCIIPrintf(ngmres->monitor, "Restarted at iteration %d\n", k_restart);CHKERRQ(ierr);
417       }
418       restart_count = 0;
419       k_restart = 1;
420       l = 1;
421       /* q_{00} = nu */
422       ngmres->fnorms[0] = fnorm;
423       nu = fnorm*fnorm;
424       Q(0,0) = nu;
425       /* Fdot[0] = F */
426       ierr = VecCopy(X, Xdot[0]);CHKERRQ(ierr);
427       ierr = VecCopy(F, Fdot[0]);CHKERRQ(ierr);
428     } else {
429       /* select the current size of the subspace */
430       if (l < ngmres->msize) l++;
431       k_restart++;
432       /* place the current entry in the list of previous entries */
433       ierr = VecCopy(F, Fdot[ivec]);CHKERRQ(ierr);
434       ierr = VecCopy(X, Xdot[ivec]);CHKERRQ(ierr);
435       ngmres->fnorms[ivec] = fnorm;
436       if (fminnorm > fnorm) fminnorm = fnorm;  /* the minimum norm is now of r^A */
437       for (i = 0; i < l; i++) {
438         ierr = VecDot(F, Fdot[i], &qentry);CHKERRQ(ierr);
439         Q(i, ivec) = qentry;
440         Q(ivec, i) = qentry;
441       }
442     }
443 
444     ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr);
445     snes->iter = k;
446     snes->norm = fnorm;
447     ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr);
448     SNESLogConvHistory(snes, snes->norm, snes->iter);
449     ierr = SNESMonitor(snes, snes->iter, snes->norm);CHKERRQ(ierr);
450 
451     ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
452     if (snes->reason) PetscFunctionReturn(0);
453   }
454   snes->reason = SNES_DIVERGED_MAX_IT;
455   PetscFunctionReturn(0);
456 }
457 
458 /*MC
459   SNESNGMRES - The Nonlinear Generalized Minimum Residual (NGMRES) method of Oosterlee and Washio.
460 
461    Level: beginner
462 
463    "Krylov Subspace Acceleration of Nonlinear Multigrid with Application to Recirculating Flows", C. W. Oosterlee and T. Washio,
464    SIAM Journal on Scientific Computing, 21(5), 2000.
465 
466    This is also the same as the algorithm called Anderson acceleration introduced in "D. G. Anderson. Iterative procedures for nonlinear integral equations.
467    J. Assoc. Comput. Mach., 12:547–560, 1965."
468 
469 .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types)
470 M*/
471 
472 EXTERN_C_BEGIN
473 #undef __FUNCT__
474 #define __FUNCT__ "SNESCreate_NGMRES"
475 PetscErrorCode SNESCreate_NGMRES(SNES snes)
476 {
477   SNES_NGMRES   *ngmres;
478   PetscErrorCode ierr;
479 
480   PetscFunctionBegin;
481   snes->ops->destroy        = SNESDestroy_NGMRES;
482   snes->ops->setup          = SNESSetUp_NGMRES;
483   snes->ops->setfromoptions = SNESSetFromOptions_NGMRES;
484   snes->ops->view           = SNESView_NGMRES;
485   snes->ops->solve          = SNESSolve_NGMRES;
486   snes->ops->reset          = SNESReset_NGMRES;
487 
488   snes->usespc          = PETSC_TRUE;
489   snes->usesksp         = PETSC_FALSE;
490 
491   ierr = PetscNewLog(snes, SNES_NGMRES, &ngmres);CHKERRQ(ierr);
492   snes->data = (void*) ngmres;
493   ngmres->msize = 10;
494 
495   ngmres->restart_it = 2;
496   ngmres->gammaA     = 2.0;
497   ngmres->gammaC     = 2.0;
498   ngmres->deltaB     = 0.9;
499   ngmres->epsilonB   = 0.1;
500   ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetType_C","SNESLineSearchSetType_NGMRES",SNESLineSearchSetType_NGMRES);CHKERRQ(ierr);
501   ierr = SNESLineSearchSetType(snes, SNES_LS_QUADRATIC);CHKERRQ(ierr);
502   PetscFunctionReturn(0);
503 }
504 EXTERN_C_END
505