xref: /petsc/src/snes/impls/ngmres/snesngmres.c (revision fc8a9adeb7fcdc98711d755fa2dc544ddccf0f3e)
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 .keywords: SNES, SNESNGMRES, restart, type, set SNESLineSearch
420 @*/
421 PetscErrorCode SNESNGMRESSetRestartType(SNES snes,SNESNGMRESRestartType rtype)
422 {
423   PetscErrorCode ierr;
424 
425   PetscFunctionBegin;
426   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
427   ierr = PetscTryMethod(snes,"SNESNGMRESSetRestartType_C",(SNES,SNESNGMRESRestartType),(snes,rtype));CHKERRQ(ierr);
428   PetscFunctionReturn(0);
429 }
430 
431 /*@
432     SNESNGMRESSetSelectType - Sets the selection type for SNESNGMRES.  This determines how the candidate solution and
433     combined solution are used to create the next iterate.
434 
435     Logically Collective on SNES
436 
437     Input Parameters:
438 +   snes - the iterative context
439 -   stype - selection type
440 
441     Options Database:
442 .   -snes_ngmres_select_type<difference,none,linesearch>
443 
444     Level: intermediate
445 
446     SNESNGMRESSelectTypes:
447 +   SNES_NGMRES_SELECT_NONE - choose the combined solution all the time
448 .   SNES_NGMRES_SELECT_DIFFERENCE - choose based upon the selection criteria
449 -   SNES_NGMRES_SELECT_LINESEARCH - choose based upon line search combination
450 
451     Notes:
452     The default line search used is the L2 line search and it requires two additional function evaluations.
453 
454 .keywords: SNES, SNESNGMRES, selection, type, set SNESLineSearch
455 @*/
456 PetscErrorCode SNESNGMRESSetSelectType(SNES snes,SNESNGMRESSelectType stype)
457 {
458   PetscErrorCode ierr;
459 
460   PetscFunctionBegin;
461   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
462   ierr = PetscTryMethod(snes,"SNESNGMRESSetSelectType_C",(SNES,SNESNGMRESSelectType),(snes,stype));CHKERRQ(ierr);
463   PetscFunctionReturn(0);
464 }
465 
466 PetscErrorCode SNESNGMRESSetSelectType_NGMRES(SNES snes,SNESNGMRESSelectType stype)
467 {
468   SNES_NGMRES *ngmres = (SNES_NGMRES*)snes->data;
469 
470   PetscFunctionBegin;
471   ngmres->select_type = stype;
472   PetscFunctionReturn(0);
473 }
474 
475 PetscErrorCode SNESNGMRESSetRestartType_NGMRES(SNES snes,SNESNGMRESRestartType rtype)
476 {
477   SNES_NGMRES *ngmres = (SNES_NGMRES*)snes->data;
478 
479   PetscFunctionBegin;
480   ngmres->restart_type = rtype;
481   PetscFunctionReturn(0);
482 }
483 
484 /*MC
485   SNESNGMRES - The Nonlinear Generalized Minimum Residual method.
486 
487    Level: beginner
488 
489    Options Database:
490 +  -snes_ngmres_select_type<difference,none,linesearch> - choose the select between candidate and combined solution
491 .  -snes_ngmres_restart_type<difference,none,periodic> - choose the restart conditions
492 .  -snes_ngmres_candidate        - Use NGMRES variant which combines candidate solutions instead of actual solutions
493 .  -snes_ngmres_m                - Number of stored previous solutions and residuals
494 .  -snes_ngmres_restart_it       - Number of iterations the restart conditions hold before restart
495 .  -snes_ngmres_gammaA           - Residual tolerance for solution select between the candidate and combination
496 .  -snes_ngmres_gammaC           - Residual tolerance for restart
497 .  -snes_ngmres_epsilonB         - Difference tolerance between subsequent solutions triggering restart
498 .  -snes_ngmres_deltaB           - Difference tolerance between residuals triggering restart
499 .  -snes_ngmres_restart_fm_rise  - Restart on residual rise from x_M step
500 .  -snes_ngmres_monitor          - Prints relevant information about the ngmres iteration
501 .  -snes_linesearch_type <basic,l2,cp> - Line search type used for the default smoother
502 -  -additive_snes_linesearch_type - linesearch type used to select between the candidate and combined solution with additive select type
503 
504    Notes:
505 
506    The N-GMRES method combines m previous solutions into a minimum-residual solution by solving a small linearized
507    optimization problem at each iteration.
508 
509    Very similar to the SNESANDERSON algorithm.
510 
511    References:
512 +  1. - C. W. Oosterlee and T. Washio, "Krylov Subspace Acceleration of Nonlinear Multigrid with Application to Recirculating Flows",
513    SIAM Journal on Scientific Computing, 21(5), 2000.
514 -  2. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
515    SIAM Review, 57(4), 2015
516 
517 
518 .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types)
519 M*/
520 
521 PETSC_EXTERN PetscErrorCode SNESCreate_NGMRES(SNES snes)
522 {
523   SNES_NGMRES    *ngmres;
524   PetscErrorCode ierr;
525   SNESLineSearch linesearch;
526 
527   PetscFunctionBegin;
528   snes->ops->destroy        = SNESDestroy_NGMRES;
529   snes->ops->setup          = SNESSetUp_NGMRES;
530   snes->ops->setfromoptions = SNESSetFromOptions_NGMRES;
531   snes->ops->view           = SNESView_NGMRES;
532   snes->ops->solve          = SNESSolve_NGMRES;
533   snes->ops->reset          = SNESReset_NGMRES;
534 
535   snes->usesnpc  = PETSC_TRUE;
536   snes->usesksp  = PETSC_FALSE;
537   snes->npcside  = PC_RIGHT;
538 
539   snes->alwayscomputesfinalresidual = PETSC_TRUE;
540 
541   ierr          = PetscNewLog(snes,&ngmres);CHKERRQ(ierr);
542   snes->data    = (void*) ngmres;
543   ngmres->msize = 30;
544 
545   if (!snes->tolerancesset) {
546     snes->max_funcs = 30000;
547     snes->max_its   = 10000;
548   }
549 
550   ngmres->candidate = PETSC_FALSE;
551 
552  ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr);
553  ierr = SNESLineSearchSetType(linesearch,SNESLINESEARCHBASIC);CHKERRQ(ierr);
554 
555   ngmres->additive_linesearch = NULL;
556   ngmres->approxfunc          = PETSC_FALSE;
557   ngmres->restart_it          = 2;
558   ngmres->restart_periodic    = 30;
559   ngmres->gammaA              = 2.0;
560   ngmres->gammaC              = 2.0;
561   ngmres->deltaB              = 0.9;
562   ngmres->epsilonB            = 0.1;
563   ngmres->restart_fm_rise     = PETSC_FALSE;
564 
565   ngmres->restart_type = SNES_NGMRES_RESTART_DIFFERENCE;
566   ngmres->select_type  = SNES_NGMRES_SELECT_DIFFERENCE;
567 
568   ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESNGMRESSetSelectType_C",SNESNGMRESSetSelectType_NGMRES);CHKERRQ(ierr);
569   ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESNGMRESSetRestartType_C",SNESNGMRESSetRestartType_NGMRES);CHKERRQ(ierr);
570   ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESNGMRESSetRestartFmRise_C",SNESNGMRESSetRestartFmRise_NGMRES);CHKERRQ(ierr);
571   ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESNGMRESGetRestartFmRise_C",SNESNGMRESGetRestartFmRise_NGMRES);CHKERRQ(ierr);
572   PetscFunctionReturn(0);
573 }
574