xref: /petsc/src/snes/impls/ngmres/snesngmres.c (revision 8cdbd7576f676fd2f0cf4ce91c9fbf7f9b0d4574)
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 = PetscOptionsBool("-snes_ngmres_ngs",      "Use nonlinear Gauss-Seidel",    "SNES", ngmres->useGS, &ngmres->useGS, PETSC_NULL);CHKERRQ(ierr);
99   ierr = PetscOptionsReal("-snes_ngmres_gammaA",   "Residual selection constant",   "SNES", ngmres->gammaA, &ngmres->gammaA, PETSC_NULL);CHKERRQ(ierr);
100   ierr = PetscOptionsReal("-snes_ngmres_gammaC", "  Residual restart constant",     "SNES", ngmres->gammaC, &ngmres->gammaC, PETSC_NULL);CHKERRQ(ierr);
101   ierr = PetscOptionsReal("-snes_ngmres_epsilonB", "Difference selection constant", "SNES", ngmres->epsilonB, &ngmres->epsilonB, PETSC_NULL);CHKERRQ(ierr);
102   ierr = PetscOptionsReal("-snes_ngmres_deltaB",   "Difference residual selection constant", "SNES", ngmres->deltaB, &ngmres->deltaB, PETSC_NULL);CHKERRQ(ierr);
103   ierr = PetscOptionsTail();CHKERRQ(ierr);
104   if ((ngmres->gammaA > ngmres->gammaC) && (ngmres->gammaC > 2.)) ngmres->gammaC = ngmres->gammaA;
105   PetscFunctionReturn(0);
106 }
107 
108 #undef __FUNCT__
109 #define __FUNCT__ "SNESView_NGMRES"
110 PetscErrorCode SNESView_NGMRES(SNES snes, PetscViewer viewer)
111 {
112   SNES_NGMRES   *ngmres = (SNES_NGMRES *) snes->data;
113   PetscBool      iascii;
114   const char     *cstr;
115   PetscErrorCode ierr;
116 
117   PetscFunctionBegin;
118   ierr = PetscTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
119   if (iascii) {
120     cstr = SNESLineSearchTypeName(snes->ls_type);
121     ierr = PetscViewerASCIIPrintf(viewer, "  line search variant: %s\n",cstr);CHKERRQ(ierr);
122     ierr = PetscViewerASCIIPrintf(viewer, "  Number of stored past updates: %d\n", ngmres->msize);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, restart_count = 0;
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       ierr = SNESSolve(snes->pc, B, X);CHKERRQ(ierr);
250       ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr);
251       if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
252         snes->reason = SNES_DIVERGED_INNER;
253         PetscFunctionReturn(0);
254       }
255       ierr = SNESComputeFunction(snes, X, F);CHKERRQ(ierr);
256       ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr);
257     } else if (ngmres->useGS && snes->ops->computegs) {
258       /* compute the update using the supplied Gauss-Seidel routine */
259       ierr = SNESComputeGS(snes, B, X);CHKERRQ(ierr);
260       ierr = SNESComputeFunction(snes, X, F);CHKERRQ(ierr);
261       ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr);
262     } else {
263       /* no preconditioner -- just take gradient descent with line search */
264       ierr = VecCopy(F, Y);CHKERRQ(ierr);
265       ierr = VecScale(Y, -1.0);CHKERRQ(ierr);
266       ierr = (*snes->ops->linesearch)(snes,snes->lsP,X,F,Y,fnorm,xnorm,G,W,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr);
267       if (!lssucceed) {
268         if (++snes->numFailures >= snes->maxFailures) {
269           snes->reason = SNES_DIVERGED_LINE_SEARCH;
270           PetscFunctionReturn(0);
271         }
272       }
273       fnorm = gnorm;
274       ierr = VecCopy(G, F);CHKERRQ(ierr);
275       ierr = VecCopy(W, X);CHKERRQ(ierr);
276     }
277 
278     /* r = F(x) */
279     nu = fnorm*fnorm;
280     if (fminnorm > fnorm) fminnorm = fnorm;  /* the minimum norm is now of F^M */
281 
282     /* construct the right hand side and xi factors */
283     for (i = 0; i < l; i++) {
284       ierr = VecDot(Fdot[i], F, &xi[i]);CHKERRQ(ierr);
285       beta[i] = nu - xi[i];
286     }
287 
288     /* construct h */
289     for (j = 0; j < l; j++) {
290       for (i = 0; i < l; i++) {
291         H(i, j) = Q(i, j) - xi[i] - xi[j] + nu;
292       }
293     }
294 
295     if(l == 1) {
296       /* simply set alpha[0] = beta[0] / H[0, 0] */
297       beta[0] = beta[0] / H(0, 0);
298     } else {
299 #ifdef PETSC_MISSING_LAPACK_GELSS
300     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "NGMRES with LS requires the LAPACK GELSS routine.");
301 #else
302     ngmres->m = PetscBLASIntCast(l);
303     ngmres->n = PetscBLASIntCast(l);
304     ngmres->info = PetscBLASIntCast(0);
305     ngmres->rcond = -1.;
306 #ifdef PETSC_USE_COMPLEX
307     LAPACKgelss_(&ngmres->m,
308                  &ngmres->n,
309                  &ngmres->nrhs,
310                  ngmres->h,
311                  &ngmres->lda,
312                  ngmres->beta,
313                  &ngmres->ldb,
314                  ngmres->s,
315                  &ngmres->rcond,
316                  &ngmres->rank,
317                  ngmres->work,
318                  &ngmres->lwork,
319                  ngmres->rwork,
320                  &ngmres->info);
321 #else
322     LAPACKgelss_(&ngmres->m,
323                  &ngmres->n,
324                  &ngmres->nrhs,
325                  ngmres->h,
326                  &ngmres->lda,
327                  ngmres->beta,
328                  &ngmres->ldb,
329                  ngmres->s,
330                  &ngmres->rcond,
331                  &ngmres->rank,
332                  ngmres->work,
333                  &ngmres->lwork,
334                  &ngmres->info);
335 #endif
336     if (ngmres->info < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELSS");
337     if (ngmres->info > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SVD failed to converge");
338 #endif
339     }
340 
341     alph_total = 0.;
342     for (i = 0; i < l; i++) {
343       alph_total += beta[i];
344     }
345 
346     ierr = VecCopy(X, XA);CHKERRQ(ierr);
347     ierr = VecScale(XA, 1. - alph_total);CHKERRQ(ierr);
348 
349     for(i=0;i<l;i++){
350       ierr= VecAXPY(XA, beta[i], Xdot[i]);CHKERRQ(ierr);
351     }
352     ierr = SNESComputeFunction(snes, XA, FA);CHKERRQ(ierr);
353     ierr = VecNorm(FA, NORM_2, &fAnorm);CHKERRQ(ierr);
354 
355     selectA = PETSC_TRUE;
356     /* Conditions for choosing the accelerated answer */
357 
358     /* Criterion A -- the norm of the function isn't increased above the minimum by too much */
359     if (fAnorm >= ngmres->gammaA*fminnorm) {
360       selectA = PETSC_FALSE;
361     }
362 
363     /* Criterion B -- the choice of x^A isn't too close to some other choice */
364     ierr=VecCopy(XA, D);CHKERRQ(ierr);
365     ierr=VecAXPY(D, -1.0, X);CHKERRQ(ierr);
366     ierr=VecNorm(D, NORM_2, &dnorm);CHKERRQ(ierr);
367     dminnorm = -1.0;
368     for(i=0;i<l;i++) {
369       ierr=VecCopy(XA, D);CHKERRQ(ierr);
370       ierr=VecAXPY(D, -1.0, Xdot[i]);CHKERRQ(ierr);
371       ierr=VecNorm(D, NORM_2, &dcurnorm);CHKERRQ(ierr);
372       if((dcurnorm < dminnorm) || (dminnorm < 0.0)) dminnorm = dcurnorm;
373     }
374     if (ngmres->epsilonB*dnorm<dminnorm || PetscSqrtReal(fnorm)<ngmres->deltaB*PetscSqrtReal(fminnorm)) {
375     } else {
376       selectA=PETSC_FALSE;
377     }
378 
379 
380     if (selectA) {
381       if (ngmres->monitor) {
382         ierr = PetscViewerASCIIPrintf(ngmres->monitor, "picked X_A, ||F_A||_2 = %e, ||F_M||_2 = %e\n", fAnorm, fnorm);CHKERRQ(ierr);
383       }
384       /* copy it over */
385       fnorm = fAnorm;
386       nu = fnorm*fnorm;
387       ierr = VecCopy(FA, F);CHKERRQ(ierr);
388       ierr = VecCopy(XA, X);CHKERRQ(ierr);
389     } else {
390       if (ngmres->monitor) {
391         ierr = PetscViewerASCIIPrintf(ngmres->monitor, "picked X_M, ||F_A||_2 = %e, ||F_M||_2 = %e\n", fAnorm, fnorm);CHKERRQ(ierr);
392       }
393     }
394 
395     selectRestart = PETSC_FALSE;
396 
397     /* difference stagnation restart */
398     if((ngmres->epsilonB*dnorm > dminnorm) && (PetscSqrtReal(fAnorm) > ngmres->deltaB*PetscSqrtReal(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 (PetscSqrtReal(fAnorm) > ngmres->gammaC*PetscSqrtReal(fminnorm)) {
406       if (ngmres->monitor) {
407         ierr = PetscViewerASCIIPrintf(ngmres->monitor, "residual restart: %e > %e\n", PetscSqrtReal(fAnorm), ngmres->gammaC*PetscSqrtReal(fminnorm));CHKERRQ(ierr);
408       }
409       selectRestart = PETSC_TRUE;
410     }
411 
412     /* if the restart conditions persist for more than restart_it iterations, restart. */
413     if (selectRestart) {
414       restart_count++;
415     } else {
416       restart_count = 0;
417     }
418 
419     /* restart after restart conditions have persisted for a fixed number of iterations */
420     if (restart_count >= ngmres->restart_it) {
421       if (ngmres->monitor){
422         ierr = PetscViewerASCIIPrintf(ngmres->monitor, "Restarted at iteration %d\n", k_restart);CHKERRQ(ierr);
423       }
424       restart_count = 0;
425       k_restart = 1;
426       l = 1;
427       /* q_{00} = nu */
428       ngmres->fnorms[0] = fnorm;
429       nu = fnorm*fnorm;
430       Q(0,0) = nu;
431       /* Fdot[0] = F */
432       ierr = VecCopy(X, Xdot[0]);CHKERRQ(ierr);
433       ierr = VecCopy(F, Fdot[0]);CHKERRQ(ierr);
434     } else {
435       /* select the current size of the subspace */
436       if (l < ngmres->msize) l++;
437       k_restart++;
438       /* place the current entry in the list of previous entries */
439       ierr = VecCopy(F, Fdot[ivec]);CHKERRQ(ierr);
440       ierr = VecCopy(X, Xdot[ivec]);CHKERRQ(ierr);
441       ngmres->fnorms[ivec] = fnorm;
442       if (fminnorm > fnorm) fminnorm = fnorm;  /* the minimum norm is now of r^A */
443       for (i = 0; i < l; i++) {
444         ierr = VecDot(F, Fdot[i], &qentry);CHKERRQ(ierr);
445         Q(i, ivec) = qentry;
446         Q(ivec, i) = qentry;
447       }
448     }
449 
450     ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr);
451     snes->iter = k;
452     snes->norm = fnorm;
453     ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr);
454     SNESLogConvHistory(snes, snes->norm, snes->iter);
455     ierr = SNESMonitor(snes, snes->iter, snes->norm);CHKERRQ(ierr);
456 
457     ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
458     if (snes->reason) PetscFunctionReturn(0);
459   }
460   snes->reason = SNES_DIVERGED_MAX_IT;
461   PetscFunctionReturn(0);
462 }
463 
464 /*MC
465   SNESNGMRES - The Nonlinear Generalized Minimum Residual (NGMRES) method of Oosterlee and Washio.
466 
467    Level: beginner
468 
469    "Krylov Subspace Acceleration of Nonlinear Multigrid with Application to Recirculating Flows", C. W. Oosterlee and T. Washio,
470    SIAM Journal on Scientific Computing, 21(5), 2000.
471 
472    This is also the same as the algorithm called Anderson acceleration introduced in "D. G. Anderson. Iterative procedures for nonlinear integral equations.
473    J. Assoc. Comput. Mach., 12:547–560, 1965."
474 
475 .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types)
476 M*/
477 
478 EXTERN_C_BEGIN
479 #undef __FUNCT__
480 #define __FUNCT__ "SNESCreate_NGMRES"
481 PetscErrorCode SNESCreate_NGMRES(SNES snes)
482 {
483   SNES_NGMRES   *ngmres;
484   PetscErrorCode ierr;
485 
486   PetscFunctionBegin;
487   snes->ops->destroy        = SNESDestroy_NGMRES;
488   snes->ops->setup          = SNESSetUp_NGMRES;
489   snes->ops->setfromoptions = SNESSetFromOptions_NGMRES;
490   snes->ops->view           = SNESView_NGMRES;
491   snes->ops->solve          = SNESSolve_NGMRES;
492   snes->ops->reset          = SNESReset_NGMRES;
493 
494   snes->usespc          = PETSC_TRUE;
495   snes->usesksp         = PETSC_FALSE;
496 
497   ierr = PetscNewLog(snes, SNES_NGMRES, &ngmres);CHKERRQ(ierr);
498   snes->data = (void*) ngmres;
499   ngmres->msize = 10;
500   ngmres->useGS = PETSC_FALSE;
501 
502   ngmres->restart_it = 2;
503   ngmres->gammaA     = 2.0;
504   ngmres->gammaC     = 2.0;
505   ngmres->deltaB     = 0.9;
506   ngmres->epsilonB   = 0.1;
507   ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetType_C","SNESLineSearchSetType_NGMRES",SNESLineSearchSetType_NGMRES);CHKERRQ(ierr);
508   ierr = SNESLineSearchSetType(snes, SNES_LS_QUADRATIC);CHKERRQ(ierr);
509   PetscFunctionReturn(0);
510 }
511 EXTERN_C_END
512