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