xref: /petsc/src/snes/impls/ngmres/snesngmres.c (revision bfcb38ea38335faa6e7f8d97f6bc6ce9aa2a1dd1)
1 #include <../src/snes/impls/ngmres/snesngmres.h> /*I "petscsnes.h" I*/
2 #include <petscblaslapack.h>
3 #include <petscdm.h>
4 
5 const char *const SNESNGMRESRestartTypes[] = {"NONE","PERIODIC","DIFFERENCE","SNESNGMRESRestartType","SNES_NGMRES_RESTART_",0};
6 const char *const SNESNGMRESSelectTypes[] = {"NONE","DIFFERENCE","LINESEARCH","SNESNGMRESSelectType","SNES_NGMRES_SELECT_",0};
7 
8 PetscErrorCode SNESReset_NGMRES(SNES snes)
9 {
10   SNES_NGMRES    *ngmres = (SNES_NGMRES*) snes->data;
11   PetscErrorCode ierr;
12 
13   PetscFunctionBegin;
14   ierr = VecDestroyVecs(ngmres->msize,&ngmres->Fdot);CHKERRQ(ierr);
15   ierr = VecDestroyVecs(ngmres->msize,&ngmres->Xdot);CHKERRQ(ierr);
16   ierr = SNESLineSearchDestroy(&ngmres->additive_linesearch);CHKERRQ(ierr);
17   PetscFunctionReturn(0);
18 }
19 
20 PetscErrorCode SNESDestroy_NGMRES(SNES snes)
21 {
22   PetscErrorCode ierr;
23   SNES_NGMRES    *ngmres = (SNES_NGMRES*)snes->data;
24 
25   PetscFunctionBegin;
26   ierr = SNESReset_NGMRES(snes);CHKERRQ(ierr);
27   ierr = PetscFree5(ngmres->h,ngmres->beta,ngmres->xi,ngmres->fnorms,ngmres->q);CHKERRQ(ierr);
28   ierr = PetscFree(ngmres->s);CHKERRQ(ierr);
29   ierr = PetscFree(ngmres->xnorms);CHKERRQ(ierr);
30 #if defined(PETSC_USE_COMPLEX)
31   ierr = PetscFree(ngmres->rwork);CHKERRQ(ierr);
32 #endif
33   ierr = PetscFree(ngmres->work);CHKERRQ(ierr);
34   ierr = PetscFree(snes->data);CHKERRQ(ierr);
35   PetscFunctionReturn(0);
36 }
37 
38 PetscErrorCode SNESSetUp_NGMRES(SNES snes)
39 {
40   SNES_NGMRES    *ngmres = (SNES_NGMRES*) snes->data;
41   const char     *optionsprefix;
42   PetscInt       msize,hsize;
43   PetscErrorCode ierr;
44   DM             dm;
45 
46   PetscFunctionBegin;
47   if (snes->npc && snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
48     SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"SNESNGMRES does not support left preconditioning with unpreconditioned function");
49   }
50   if (snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_PRECONDITIONED;
51   ierr = SNESSetWorkVecs(snes,5);CHKERRQ(ierr);
52 
53   if (!snes->vec_sol) {
54     ierr             = SNESGetDM(snes,&dm);CHKERRQ(ierr);
55     ierr             = DMCreateGlobalVector(dm,&snes->vec_sol);CHKERRQ(ierr);
56   }
57 
58   if (!ngmres->Xdot) {ierr = VecDuplicateVecs(snes->vec_sol,ngmres->msize,&ngmres->Xdot);CHKERRQ(ierr);}
59   if (!ngmres->Fdot) {ierr = VecDuplicateVecs(snes->vec_sol,ngmres->msize,&ngmres->Fdot);CHKERRQ(ierr);}
60   if (!ngmres->setup_called) {
61     msize = ngmres->msize;          /* restart size */
62     hsize = msize * msize;
63 
64     /* explicit least squares minimization solve */
65     ierr = PetscMalloc5(hsize,&ngmres->h, msize,&ngmres->beta, msize,&ngmres->xi, msize,&ngmres->fnorms, hsize,&ngmres->q);CHKERRQ(ierr);
66     ierr = PetscMalloc1(msize,&ngmres->xnorms);CHKERRQ(ierr);
67     ngmres->nrhs  = 1;
68     ngmres->lda   = msize;
69     ngmres->ldb   = msize;
70     ierr          = PetscMalloc1(msize,&ngmres->s);CHKERRQ(ierr);
71     ierr          = PetscMemzero(ngmres->h,   hsize*sizeof(PetscScalar));CHKERRQ(ierr);
72     ierr          = PetscMemzero(ngmres->q,   hsize*sizeof(PetscScalar));CHKERRQ(ierr);
73     ierr          = PetscMemzero(ngmres->xi,  msize*sizeof(PetscScalar));CHKERRQ(ierr);
74     ierr          = PetscMemzero(ngmres->beta,msize*sizeof(PetscScalar));CHKERRQ(ierr);
75     ngmres->lwork = 12*msize;
76 #if defined(PETSC_USE_COMPLEX)
77     ierr = PetscMalloc1(ngmres->lwork,&ngmres->rwork);CHKERRQ(ierr);
78 #endif
79     ierr = PetscMalloc1(ngmres->lwork,&ngmres->work);CHKERRQ(ierr);
80   }
81 
82   /* linesearch setup */
83   ierr = SNESGetOptionsPrefix(snes,&optionsprefix);CHKERRQ(ierr);
84 
85   if (ngmres->select_type == SNES_NGMRES_SELECT_LINESEARCH) {
86     ierr = SNESLineSearchCreate(PetscObjectComm((PetscObject)snes),&ngmres->additive_linesearch);CHKERRQ(ierr);
87     ierr = SNESLineSearchSetSNES(ngmres->additive_linesearch,snes);CHKERRQ(ierr);
88     ierr = SNESLineSearchSetType(ngmres->additive_linesearch,SNESLINESEARCHL2);CHKERRQ(ierr);
89     ierr = SNESLineSearchAppendOptionsPrefix(ngmres->additive_linesearch,"additive_");CHKERRQ(ierr);
90     ierr = SNESLineSearchAppendOptionsPrefix(ngmres->additive_linesearch,optionsprefix);CHKERRQ(ierr);
91     ierr = SNESLineSearchSetFromOptions(ngmres->additive_linesearch);CHKERRQ(ierr);
92   }
93 
94   ngmres->setup_called = PETSC_TRUE;
95   PetscFunctionReturn(0);
96 }
97 
98 PetscErrorCode SNESSetFromOptions_NGMRES(PetscOptionItems *PetscOptionsObject,SNES snes)
99 {
100   SNES_NGMRES    *ngmres = (SNES_NGMRES*) snes->data;
101   PetscErrorCode ierr;
102   PetscBool      debug = PETSC_FALSE;
103 
104   PetscFunctionBegin;
105   ierr = PetscOptionsHead(PetscOptionsObject,"SNES NGMRES options");CHKERRQ(ierr);
106   ierr = PetscOptionsEnum("-snes_ngmres_select_type","Select type","SNESNGMRESSetSelectType",SNESNGMRESSelectTypes,
107                           (PetscEnum)ngmres->select_type,(PetscEnum*)&ngmres->select_type,NULL);CHKERRQ(ierr);
108   ierr = PetscOptionsEnum("-snes_ngmres_restart_type","Restart type","SNESNGMRESSetRestartType",SNESNGMRESRestartTypes,
109                           (PetscEnum)ngmres->restart_type,(PetscEnum*)&ngmres->restart_type,NULL);CHKERRQ(ierr);
110   ierr = PetscOptionsBool("-snes_ngmres_candidate", "Use candidate storage",              "SNES",ngmres->candidate,&ngmres->candidate,NULL);CHKERRQ(ierr);
111   ierr = PetscOptionsBool("-snes_ngmres_approxfunc","Linearly approximate the function", "SNES",ngmres->approxfunc,&ngmres->approxfunc,NULL);CHKERRQ(ierr);
112   ierr = PetscOptionsInt("-snes_ngmres_m",          "Number of directions",               "SNES",ngmres->msize,&ngmres->msize,NULL);CHKERRQ(ierr);
113   ierr = PetscOptionsInt("-snes_ngmres_restart",    "Iterations before forced restart",   "SNES",ngmres->restart_periodic,&ngmres->restart_periodic,NULL);CHKERRQ(ierr);
114   ierr = PetscOptionsInt("-snes_ngmres_restart_it", "Tolerance iterations before restart","SNES",ngmres->restart_it,&ngmres->restart_it,NULL);CHKERRQ(ierr);
115   ierr = PetscOptionsBool("-snes_ngmres_monitor",   "Monitor actions of NGMRES",          "SNES",ngmres->monitor ? PETSC_TRUE : PETSC_FALSE,&debug,NULL);CHKERRQ(ierr);
116   if (debug) {
117     ngmres->monitor = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes));CHKERRQ(ierr);
118   }
119   ierr = PetscOptionsReal("-snes_ngmres_gammaA",    "Residual selection constant",   "SNES",ngmres->gammaA,&ngmres->gammaA,NULL);CHKERRQ(ierr);
120   ierr = PetscOptionsReal("-snes_ngmres_gammaC",    "Residual restart constant",     "SNES",ngmres->gammaC,&ngmres->gammaC,NULL);CHKERRQ(ierr);
121   ierr = PetscOptionsReal("-snes_ngmres_epsilonB",  "Difference selection constant", "SNES",ngmres->epsilonB,&ngmres->epsilonB,NULL);CHKERRQ(ierr);
122   ierr = PetscOptionsReal("-snes_ngmres_deltaB",    "Difference residual selection constant", "SNES",ngmres->deltaB,&ngmres->deltaB,NULL);CHKERRQ(ierr);
123   ierr = PetscOptionsBool("-snes_ngmres_single_reduction", "Aggregate reductions",  "SNES",ngmres->singlereduction,&ngmres->singlereduction,NULL);CHKERRQ(ierr);
124   ierr = PetscOptionsBool("-snes_ngmres_restart_fm_rise", "Restart on F_M residual rise",  "SNESNGMRESSetRestartFmRise",ngmres->restart_fm_rise,&ngmres->restart_fm_rise,NULL);CHKERRQ(ierr);
125   ierr = PetscOptionsTail();CHKERRQ(ierr);
126   if ((ngmres->gammaA > ngmres->gammaC) && (ngmres->gammaC > 2.)) ngmres->gammaC = ngmres->gammaA;
127   PetscFunctionReturn(0);
128 }
129 
130 PetscErrorCode SNESView_NGMRES(SNES snes,PetscViewer viewer)
131 {
132   SNES_NGMRES    *ngmres = (SNES_NGMRES*) snes->data;
133   PetscBool      iascii;
134   PetscErrorCode ierr;
135 
136   PetscFunctionBegin;
137   ierr = PetscObjectTypeCompare((PetscObject) viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
138   if (iascii) {
139     ierr = PetscViewerASCIIPrintf(viewer,"  Number of stored past updates: %d\n", ngmres->msize);CHKERRQ(ierr);
140     ierr = PetscViewerASCIIPrintf(viewer,"  Residual selection: gammaA=%1.0e, gammaC=%1.0e\n",ngmres->gammaA,ngmres->gammaC);CHKERRQ(ierr);
141     ierr = PetscViewerASCIIPrintf(viewer,"  Difference restart: epsilonB=%1.0e, deltaB=%1.0e\n",ngmres->epsilonB,ngmres->deltaB);CHKERRQ(ierr);
142     ierr = PetscViewerASCIIPrintf(viewer,"  Restart on F_M residual increase: %s\n",ngmres->restart_fm_rise?"TRUE":"FALSE");CHKERRQ(ierr);
143   }
144   PetscFunctionReturn(0);
145 }
146 
147 PetscErrorCode SNESSolve_NGMRES(SNES snes)
148 {
149   SNES_NGMRES          *ngmres = (SNES_NGMRES*) snes->data;
150   /* present solution, residual, and preconditioned residual */
151   Vec                  X,F,B,D,Y;
152 
153   /* candidate linear combination answers */
154   Vec                  XA,FA,XM,FM;
155 
156   /* coefficients and RHS to the minimization problem */
157   PetscReal            fnorm,fMnorm,fAnorm;
158   PetscReal            xnorm,xMnorm,xAnorm;
159   PetscReal            ynorm,yMnorm,yAnorm;
160   PetscInt             k,k_restart,l,ivec,restart_count = 0;
161 
162   /* solution selection data */
163   PetscBool            selectRestart;
164   /*
165       These two variables are initialized to prevent compilers/analyzers from producing false warnings about these variables being passed
166       to SNESNGMRESSelect_Private() without being set when SNES_NGMRES_RESTART_DIFFERENCE, the values are not used in the subroutines in that case
167       so the code is correct as written.
168   */
169   PetscReal            dnorm = 0.0,dminnorm = 0.0;
170   PetscReal            fminnorm;
171 
172   SNESConvergedReason  reason;
173   SNESLineSearchReason lssucceed;
174   PetscErrorCode       ierr;
175 
176   PetscFunctionBegin;
177   if (snes->xl || snes->xu || snes->ops->computevariablebounds) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);
178 
179   ierr = PetscCitationsRegister(SNESCitation,&SNEScite);CHKERRQ(ierr);
180   /* variable initialization */
181   snes->reason = SNES_CONVERGED_ITERATING;
182   X            = snes->vec_sol;
183   F            = snes->vec_func;
184   B            = snes->vec_rhs;
185   XA           = snes->vec_sol_update;
186   FA           = snes->work[0];
187   D            = snes->work[1];
188 
189   /* work for the line search */
190   Y  = snes->work[2];
191   XM = snes->work[3];
192   FM = snes->work[4];
193 
194   ierr       = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
195   snes->iter = 0;
196   snes->norm = 0.;
197   ierr       = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
198 
199   /* initialization */
200 
201   if (snes->npc && snes->npcside== PC_LEFT) {
202     ierr = SNESApplyNPC(snes,X,NULL,F);CHKERRQ(ierr);
203     ierr = SNESGetConvergedReason(snes->npc,&reason);CHKERRQ(ierr);
204     if (reason < 0  && reason != SNES_DIVERGED_MAX_IT) {
205       snes->reason = SNES_DIVERGED_INNER;
206       PetscFunctionReturn(0);
207     }
208     ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr);
209   } else {
210     if (!snes->vec_func_init_set) {
211       ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr);
212     } else snes->vec_func_init_set = PETSC_FALSE;
213 
214     ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr);
215     SNESCheckFunctionNorm(snes,fnorm);
216   }
217   fminnorm = fnorm;
218 
219   ierr       = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
220   snes->norm = fnorm;
221   ierr       = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
222   ierr       = SNESLogConvergenceHistory(snes,fnorm,0);CHKERRQ(ierr);
223   ierr       = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr);
224   ierr       = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
225   if (snes->reason) PetscFunctionReturn(0);
226   SNESNGMRESUpdateSubspace_Private(snes,0,0,F,fnorm,X);
227 
228   k_restart = 1;
229   l         = 1;
230   ivec      = 0;
231   for (k=1; k < snes->max_its+1; k++) {
232     /* Computation of x^M */
233     if (snes->npc && snes->npcside== PC_RIGHT) {
234       ierr = VecCopy(X,XM);CHKERRQ(ierr);
235       ierr = SNESSetInitialFunction(snes->npc,F);CHKERRQ(ierr);
236 
237       ierr = PetscLogEventBegin(SNES_NPCSolve,snes->npc,XM,B,0);CHKERRQ(ierr);
238       ierr = SNESSolve(snes->npc,B,XM);CHKERRQ(ierr);
239       ierr = PetscLogEventEnd(SNES_NPCSolve,snes->npc,XM,B,0);CHKERRQ(ierr);
240 
241       ierr = SNESGetConvergedReason(snes->npc,&reason);CHKERRQ(ierr);
242       if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
243         snes->reason = SNES_DIVERGED_INNER;
244         PetscFunctionReturn(0);
245       }
246       ierr = SNESGetNPCFunction(snes,FM,&fMnorm);CHKERRQ(ierr);
247     } else {
248       /* no preconditioner -- just take gradient descent with line search */
249       ierr = VecCopy(F,Y);CHKERRQ(ierr);
250       ierr = VecCopy(F,FM);CHKERRQ(ierr);
251       ierr = VecCopy(X,XM);CHKERRQ(ierr);
252 
253       fMnorm = fnorm;
254 
255       ierr = SNESLineSearchApply(snes->linesearch,XM,FM,&fMnorm,Y);CHKERRQ(ierr);
256       ierr = SNESLineSearchGetReason(snes->linesearch,&lssucceed);CHKERRQ(ierr);
257       if (lssucceed) {
258         if (++snes->numFailures >= snes->maxFailures) {
259           snes->reason = SNES_DIVERGED_LINE_SEARCH;
260           PetscFunctionReturn(0);
261         }
262       }
263     }
264 
265     ierr = SNESNGMRESFormCombinedSolution_Private(snes,ivec,l,XM,FM,fMnorm,X,XA,FA);CHKERRQ(ierr);
266     /* r = F(x) */
267     if (fminnorm > fMnorm) fminnorm = fMnorm;  /* the minimum norm is now of F^M */
268 
269     /* differences for selection and restart */
270     if (ngmres->restart_type == SNES_NGMRES_RESTART_DIFFERENCE || ngmres->select_type == SNES_NGMRES_SELECT_DIFFERENCE) {
271       ierr = SNESNGMRESNorms_Private(snes,l,X,F,XM,FM,XA,FA,D,&dnorm,&dminnorm,&xMnorm,NULL,&yMnorm,&xAnorm,&fAnorm,&yAnorm);CHKERRQ(ierr);
272     } else {
273       ierr = SNESNGMRESNorms_Private(snes,l,X,F,XM,FM,XA,FA,D,NULL,NULL,&xMnorm,NULL,&yMnorm,&xAnorm,&fAnorm,&yAnorm);CHKERRQ(ierr);
274     }
275     SNESCheckFunctionNorm(snes,fnorm);
276 
277     /* combination (additive) or selection (multiplicative) of the N-GMRES solution */
278     ierr          = SNESNGMRESSelect_Private(snes,k_restart,XM,FM,xMnorm,fMnorm,yMnorm,XA,FA,xAnorm,fAnorm,yAnorm,dnorm,fminnorm,dminnorm,X,F,Y,&xnorm,&fnorm,&ynorm);CHKERRQ(ierr);
279     selectRestart = PETSC_FALSE;
280 
281     if (ngmres->restart_type == SNES_NGMRES_RESTART_DIFFERENCE) {
282       ierr = SNESNGMRESSelectRestart_Private(snes,l,fMnorm,fAnorm,dnorm,fminnorm,dminnorm,&selectRestart);CHKERRQ(ierr);
283 
284       /* if the restart conditions persist for more than restart_it iterations, restart. */
285       if (selectRestart) restart_count++;
286       else restart_count = 0;
287     } else if (ngmres->restart_type == SNES_NGMRES_RESTART_PERIODIC) {
288       if (k_restart > ngmres->restart_periodic) {
289         if (ngmres->monitor) ierr = PetscViewerASCIIPrintf(ngmres->monitor,"periodic restart after %D iterations\n",k_restart);CHKERRQ(ierr);
290         restart_count = ngmres->restart_it;
291       }
292     }
293 
294     ivec = k_restart % ngmres->msize; /* replace the last used part of the subspace */
295 
296     /* restart after restart conditions have persisted for a fixed number of iterations */
297     if (restart_count >= ngmres->restart_it) {
298       if (ngmres->monitor) {
299         ierr = PetscViewerASCIIPrintf(ngmres->monitor,"Restarted at iteration %d\n",k_restart);CHKERRQ(ierr);
300       }
301       restart_count = 0;
302       k_restart     = 1;
303       l             = 1;
304       ivec          = 0;
305       /* q_{00} = nu */
306       ierr = SNESNGMRESUpdateSubspace_Private(snes,0,0,FM,fMnorm,XM);CHKERRQ(ierr);
307     } else {
308       /* select the current size of the subspace */
309       if (l < ngmres->msize) l++;
310       k_restart++;
311       /* place the current entry in the list of previous entries */
312       if (ngmres->candidate) {
313         if (fminnorm > fMnorm) fminnorm = fMnorm;
314         ierr = SNESNGMRESUpdateSubspace_Private(snes,ivec,l,FM,fMnorm,XM);CHKERRQ(ierr);
315       } else {
316         if (fminnorm > fnorm) fminnorm = fnorm;
317         ierr = SNESNGMRESUpdateSubspace_Private(snes,ivec,l,F,fnorm,X);CHKERRQ(ierr);
318       }
319     }
320 
321     ierr       = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
322     snes->iter = k;
323     snes->norm = fnorm;
324     ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
325     ierr = SNESLogConvergenceHistory(snes,snes->norm,snes->iter);CHKERRQ(ierr);
326     ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr);
327     ierr = (*snes->ops->converged)(snes,snes->iter,0,0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
328     if (snes->reason) PetscFunctionReturn(0);
329   }
330   snes->reason = SNES_DIVERGED_MAX_IT;
331   PetscFunctionReturn(0);
332 }
333 
334 /*@
335  SNESNGMRESSetRestartFmRise - Increase the restart count if the step x_M increases the residual F_M
336 
337   Input Parameters:
338   +  snes - the SNES context.
339   -  flg  - boolean value deciding whether to use the option or not
340 
341   Options Database:
342   + -snes_ngmres_restart_fm_rise - Increase the restart count if the step x_M increases the residual F_M
343 
344   Level: intermediate
345 
346   Notes:
347   If the proposed step x_M increases the residual F_M, it might be trying to get out of a stagnation area.
348   To help the solver do that, reset the Krylov subspace whenever F_M increases.
349 
350   This option must be used with SNES_NGMRES_RESTART_DIFFERENCE
351 
352   The default is FALSE.
353   .seealso: SNES_NGMRES_RESTART_DIFFERENCE
354   @*/
355 PetscErrorCode SNESNGMRESSetRestartFmRise(SNES snes,PetscBool flg)
356 {
357     PetscErrorCode (*f)(SNES,PetscBool);
358     PetscErrorCode ierr;
359 
360     PetscFunctionBegin;
361     ierr = PetscObjectQueryFunction((PetscObject)snes,"SNESNGMRESSetRestartFmRise_C",&f);CHKERRQ(ierr);
362     if (f) {ierr = (f)(snes,flg);CHKERRQ(ierr);}
363     PetscFunctionReturn(0);
364 }
365 
366 PetscErrorCode SNESNGMRESSetRestartFmRise_NGMRES(SNES snes,PetscBool flg)
367 {
368   SNES_NGMRES *ngmres = (SNES_NGMRES*)snes->data;
369 
370   PetscFunctionBegin;
371   ngmres->restart_fm_rise = flg;
372   PetscFunctionReturn(0);
373 }
374 
375 PetscErrorCode SNESNGMRESGetRestartFmRise(SNES snes,PetscBool *flg)
376 {
377     PetscErrorCode (*f)(SNES,PetscBool*);
378     PetscErrorCode ierr;
379 
380     PetscFunctionBegin;
381     ierr = PetscObjectQueryFunction((PetscObject)snes,"SNESNGMRESGetRestartFmRise_C",&f);CHKERRQ(ierr);
382     if (f) {ierr = (f)(snes,flg);CHKERRQ(ierr);}
383     PetscFunctionReturn(0);
384 }
385 
386 PetscErrorCode SNESNGMRESGetRestartFmRise_NGMRES(SNES snes,PetscBool *flg)
387 {
388   SNES_NGMRES *ngmres = (SNES_NGMRES*)snes->data;
389 
390   PetscFunctionBegin;
391   *flg = ngmres->restart_fm_rise;
392   PetscFunctionReturn(0);
393 }
394 
395 
396 /*@
397     SNESNGMRESSetRestartType - Sets the restart type for SNESNGMRES.
398 
399     Logically Collective on SNES
400 
401     Input Parameters:
402 +   snes - the iterative context
403 -   rtype - restart type
404 
405     Options Database:
406 +   -snes_ngmres_restart_type<difference,periodic,none> - set the restart type
407 -   -snes_ngmres_restart[30] - sets the number of iterations before restart for periodic
408 
409     Level: intermediate
410 
411     SNESNGMRESRestartTypes:
412 +   SNES_NGMRES_RESTART_NONE - never restart
413 .   SNES_NGMRES_RESTART_DIFFERENCE - restart based upon difference criteria
414 -   SNES_NGMRES_RESTART_PERIODIC - restart after a fixed number of iterations
415 
416     Notes:
417     The default line search used is the L2 line search and it requires two additional function evaluations.
418 
419 @*/
420 PetscErrorCode SNESNGMRESSetRestartType(SNES snes,SNESNGMRESRestartType rtype)
421 {
422   PetscErrorCode ierr;
423 
424   PetscFunctionBegin;
425   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
426   ierr = PetscTryMethod(snes,"SNESNGMRESSetRestartType_C",(SNES,SNESNGMRESRestartType),(snes,rtype));CHKERRQ(ierr);
427   PetscFunctionReturn(0);
428 }
429 
430 /*@
431     SNESNGMRESSetSelectType - Sets the selection type for SNESNGMRES.  This determines how the candidate solution and
432     combined solution are used to create the next iterate.
433 
434     Logically Collective on SNES
435 
436     Input Parameters:
437 +   snes - the iterative context
438 -   stype - selection type
439 
440     Options Database:
441 .   -snes_ngmres_select_type<difference,none,linesearch>
442 
443     Level: intermediate
444 
445     SNESNGMRESSelectTypes:
446 +   SNES_NGMRES_SELECT_NONE - choose the combined solution all the time
447 .   SNES_NGMRES_SELECT_DIFFERENCE - choose based upon the selection criteria
448 -   SNES_NGMRES_SELECT_LINESEARCH - choose based upon line search combination
449 
450     Notes:
451     The default line search used is the L2 line search and it requires two additional function evaluations.
452 
453 @*/
454 PetscErrorCode SNESNGMRESSetSelectType(SNES snes,SNESNGMRESSelectType stype)
455 {
456   PetscErrorCode ierr;
457 
458   PetscFunctionBegin;
459   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
460   ierr = PetscTryMethod(snes,"SNESNGMRESSetSelectType_C",(SNES,SNESNGMRESSelectType),(snes,stype));CHKERRQ(ierr);
461   PetscFunctionReturn(0);
462 }
463 
464 PetscErrorCode SNESNGMRESSetSelectType_NGMRES(SNES snes,SNESNGMRESSelectType stype)
465 {
466   SNES_NGMRES *ngmres = (SNES_NGMRES*)snes->data;
467 
468   PetscFunctionBegin;
469   ngmres->select_type = stype;
470   PetscFunctionReturn(0);
471 }
472 
473 PetscErrorCode SNESNGMRESSetRestartType_NGMRES(SNES snes,SNESNGMRESRestartType rtype)
474 {
475   SNES_NGMRES *ngmres = (SNES_NGMRES*)snes->data;
476 
477   PetscFunctionBegin;
478   ngmres->restart_type = rtype;
479   PetscFunctionReturn(0);
480 }
481 
482 /*MC
483   SNESNGMRES - The Nonlinear Generalized Minimum Residual method.
484 
485    Level: beginner
486 
487    Options Database:
488 +  -snes_ngmres_select_type<difference,none,linesearch> - choose the select between candidate and combined solution
489 .  -snes_ngmres_restart_type<difference,none,periodic> - choose the restart conditions
490 .  -snes_ngmres_candidate        - Use NGMRES variant which combines candidate solutions instead of actual solutions
491 .  -snes_ngmres_m                - Number of stored previous solutions and residuals
492 .  -snes_ngmres_restart_it       - Number of iterations the restart conditions hold before restart
493 .  -snes_ngmres_gammaA           - Residual tolerance for solution select between the candidate and combination
494 .  -snes_ngmres_gammaC           - Residual tolerance for restart
495 .  -snes_ngmres_epsilonB         - Difference tolerance between subsequent solutions triggering restart
496 .  -snes_ngmres_deltaB           - Difference tolerance between residuals triggering restart
497 .  -snes_ngmres_restart_fm_rise  - Restart on residual rise from x_M step
498 .  -snes_ngmres_monitor          - Prints relevant information about the ngmres iteration
499 .  -snes_linesearch_type <basic,l2,cp> - Line search type used for the default smoother
500 -  -additive_snes_linesearch_type - linesearch type used to select between the candidate and combined solution with additive select type
501 
502    Notes:
503 
504    The N-GMRES method combines m previous solutions into a minimum-residual solution by solving a small linearized
505    optimization problem at each iteration.
506 
507    Very similar to the SNESANDERSON algorithm.
508 
509    References:
510 +  1. - C. W. Oosterlee and T. Washio, "Krylov Subspace Acceleration of Nonlinear Multigrid with Application to Recirculating Flows",
511    SIAM Journal on Scientific Computing, 21(5), 2000.
512 -  2. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
513    SIAM Review, 57(4), 2015
514 
515 
516 .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types)
517 M*/
518 
519 PETSC_EXTERN PetscErrorCode SNESCreate_NGMRES(SNES snes)
520 {
521   SNES_NGMRES    *ngmres;
522   PetscErrorCode ierr;
523   SNESLineSearch linesearch;
524 
525   PetscFunctionBegin;
526   snes->ops->destroy        = SNESDestroy_NGMRES;
527   snes->ops->setup          = SNESSetUp_NGMRES;
528   snes->ops->setfromoptions = SNESSetFromOptions_NGMRES;
529   snes->ops->view           = SNESView_NGMRES;
530   snes->ops->solve          = SNESSolve_NGMRES;
531   snes->ops->reset          = SNESReset_NGMRES;
532 
533   snes->usesnpc  = PETSC_TRUE;
534   snes->usesksp  = PETSC_FALSE;
535   snes->npcside  = PC_RIGHT;
536 
537   snes->alwayscomputesfinalresidual = PETSC_TRUE;
538 
539   ierr          = PetscNewLog(snes,&ngmres);CHKERRQ(ierr);
540   snes->data    = (void*) ngmres;
541   ngmres->msize = 30;
542 
543   if (!snes->tolerancesset) {
544     snes->max_funcs = 30000;
545     snes->max_its   = 10000;
546   }
547 
548   ngmres->candidate = PETSC_FALSE;
549 
550  ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr);
551  ierr = SNESLineSearchSetType(linesearch,SNESLINESEARCHBASIC);CHKERRQ(ierr);
552 
553   ngmres->additive_linesearch = NULL;
554   ngmres->approxfunc          = PETSC_FALSE;
555   ngmres->restart_it          = 2;
556   ngmres->restart_periodic    = 30;
557   ngmres->gammaA              = 2.0;
558   ngmres->gammaC              = 2.0;
559   ngmres->deltaB              = 0.9;
560   ngmres->epsilonB            = 0.1;
561   ngmres->restart_fm_rise     = PETSC_FALSE;
562 
563   ngmres->restart_type = SNES_NGMRES_RESTART_DIFFERENCE;
564   ngmres->select_type  = SNES_NGMRES_SELECT_DIFFERENCE;
565 
566   ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESNGMRESSetSelectType_C",SNESNGMRESSetSelectType_NGMRES);CHKERRQ(ierr);
567   ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESNGMRESSetRestartType_C",SNESNGMRESSetRestartType_NGMRES);CHKERRQ(ierr);
568   ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESNGMRESSetRestartFmRise_C",SNESNGMRESSetRestartFmRise_NGMRES);CHKERRQ(ierr);
569   ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESNGMRESGetRestartFmRise_C",SNESNGMRESGetRestartFmRise_NGMRES);CHKERRQ(ierr);
570   PetscFunctionReturn(0);
571 }
572