xref: /petsc/src/snes/impls/ncg/snesncg.c (revision 65f8aed5f7eaa1e2ef2ddeffe666264e0669c876)
1 #include <../src/snes/impls/ncg/snesncgimpl.h> /*I "petscsnes.h" I*/
2 const char *const SNESNCGTypes[] = {"FR","PRP","HS","DY","CD","SNESNCGType","SNES_NCG_",0};
3 
4 PetscErrorCode SNESReset_NCG(SNES snes)
5 {
6   PetscFunctionBegin;
7   PetscFunctionReturn(0);
8 }
9 
10 #define SNESLINESEARCHNCGLINEAR "ncglinear"
11 
12 /*
13   SNESDestroy_NCG - Destroys the private SNES_NCG context that was created with SNESCreate_NCG().
14 
15   Input Parameter:
16 . snes - the SNES context
17 
18   Application Interface Routine: SNESDestroy()
19 */
20 PetscErrorCode SNESDestroy_NCG(SNES snes)
21 {
22   PetscErrorCode ierr;
23 
24   PetscFunctionBegin;
25   ierr = PetscFree(snes->data);CHKERRQ(ierr);
26   PetscFunctionReturn(0);
27 }
28 
29 /*
30    SNESSetUp_NCG - Sets up the internal data structures for the later use
31    of the SNESNCG nonlinear solver.
32 
33    Input Parameters:
34 +  snes - the SNES context
35 -  x - the solution vector
36 
37    Application Interface Routine: SNESSetUp()
38  */
39 
40 PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_NCGLinear(SNESLineSearch);
41 
42 PetscErrorCode SNESSetUp_NCG(SNES snes)
43 {
44   PetscErrorCode ierr;
45 
46   PetscFunctionBegin;
47   ierr = SNESSetWorkVecs(snes,2);CHKERRQ(ierr);
48   if (snes->npcside== PC_RIGHT) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNESNCG only supports left preconditioning");
49   if (snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_UNPRECONDITIONED;
50   PetscFunctionReturn(0);
51 }
52 /*
53   SNESSetFromOptions_NCG - Sets various parameters for the SNESNCG method.
54 
55   Input Parameter:
56 . snes - the SNES context
57 
58   Application Interface Routine: SNESSetFromOptions()
59 */
60 static PetscErrorCode SNESSetFromOptions_NCG(PetscOptionItems *PetscOptionsObject,SNES snes)
61 {
62   SNES_NCG       *ncg = (SNES_NCG*)snes->data;
63   PetscErrorCode ierr;
64   PetscBool      debug = PETSC_FALSE;
65   SNESLineSearch linesearch;
66   SNESNCGType    ncgtype=ncg->type;
67 
68   PetscFunctionBegin;
69   ierr = PetscOptionsHead(PetscOptionsObject,"SNES NCG options");CHKERRQ(ierr);
70   ierr = PetscOptionsBool("-snes_ncg_monitor","Monitor NCG iterations","SNES",ncg->monitor ? PETSC_TRUE : PETSC_FALSE, &debug, NULL);CHKERRQ(ierr);
71   if (debug) {
72     ncg->monitor = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes));CHKERRQ(ierr);
73   }
74   ierr = PetscOptionsEnum("-snes_ncg_type","NCG Beta type used","SNESNCGSetType",SNESNCGTypes,(PetscEnum)ncg->type,(PetscEnum*)&ncgtype,NULL);CHKERRQ(ierr);
75   ierr = SNESNCGSetType(snes, ncgtype);CHKERRQ(ierr);
76   ierr = PetscOptionsTail();CHKERRQ(ierr);
77   if (!snes->linesearch) {
78     ierr = SNESGetLineSearch(snes, &linesearch);CHKERRQ(ierr);
79     if (!snes->npc) {
80       ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHCP);CHKERRQ(ierr);
81     } else {
82       ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);CHKERRQ(ierr);
83     }
84   }
85   ierr = SNESLineSearchRegister(SNESLINESEARCHNCGLINEAR, SNESLineSearchCreate_NCGLinear);CHKERRQ(ierr);
86   PetscFunctionReturn(0);
87 }
88 
89 /*
90   SNESView_NCG - Prints info from the SNESNCG data structure.
91 
92   Input Parameters:
93 + SNES - the SNES context
94 - viewer - visualization context
95 
96   Application Interface Routine: SNESView()
97 */
98 static PetscErrorCode SNESView_NCG(SNES snes, PetscViewer viewer)
99 {
100   SNES_NCG      *ncg = (SNES_NCG *) snes->data;
101   PetscBool      iascii;
102   PetscErrorCode ierr;
103 
104   PetscFunctionBegin;
105   ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
106   if (iascii) {
107     ierr = PetscViewerASCIIPrintf(viewer, "  type: %s\n", SNESNCGTypes[ncg->type]);CHKERRQ(ierr);
108   }
109   PetscFunctionReturn(0);
110 }
111 
112 PetscErrorCode SNESLineSearchApply_NCGLinear(SNESLineSearch linesearch)
113 {
114   PetscScalar    alpha, ptAp;
115   Vec            X, Y, F, W;
116   SNES           snes;
117   PetscErrorCode ierr;
118   PetscReal      *fnorm, *xnorm, *ynorm;
119 
120   PetscFunctionBegin;
121   ierr  = SNESLineSearchGetSNES(linesearch, &snes);CHKERRQ(ierr);
122   X     = linesearch->vec_sol;
123   W     = linesearch->vec_sol_new;
124   F     = linesearch->vec_func;
125   Y     = linesearch->vec_update;
126   fnorm = &linesearch->fnorm;
127   xnorm = &linesearch->xnorm;
128   ynorm = &linesearch->ynorm;
129 
130   if (!snes->jacobian) {
131     ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr);
132   }
133 
134   /*
135 
136    The exact step size for unpreconditioned linear CG is just:
137    alpha = (r, r) / (p, Ap) = (f, f) / (y, Jy)
138    */
139   ierr  = SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre);CHKERRQ(ierr);
140   ierr  = VecDot(F, F, &alpha);CHKERRQ(ierr);
141   ierr  = MatMult(snes->jacobian, Y, W);CHKERRQ(ierr);
142   ierr  = VecDot(Y, W, &ptAp);CHKERRQ(ierr);
143   alpha = alpha / ptAp;
144   ierr  = VecAXPY(X,-alpha,Y);CHKERRQ(ierr);
145   ierr  = SNESComputeFunction(snes,X,F);CHKERRQ(ierr);
146 
147   ierr = VecNorm(F,NORM_2,fnorm);CHKERRQ(ierr);
148   ierr = VecNorm(X,NORM_2,xnorm);CHKERRQ(ierr);
149   ierr = VecNorm(Y,NORM_2,ynorm);CHKERRQ(ierr);
150   PetscFunctionReturn(0);
151 }
152 
153 PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_NCGLinear(SNESLineSearch linesearch)
154 {
155   PetscFunctionBegin;
156   linesearch->ops->apply          = SNESLineSearchApply_NCGLinear;
157   linesearch->ops->destroy        = NULL;
158   linesearch->ops->setfromoptions = NULL;
159   linesearch->ops->reset          = NULL;
160   linesearch->ops->view           = NULL;
161   linesearch->ops->setup          = NULL;
162   PetscFunctionReturn(0);
163 }
164 
165 /*
166 
167  Assuming F = SNESComputeFunction(X) compute Y^tJ^tF using a simple secant approximation of the jacobian.
168 
169  */
170 PetscErrorCode SNESNCGComputeYtJtF_Private(SNES snes, Vec X, Vec F, Vec Y, Vec W, Vec G, PetscScalar * ytJtf)
171 {
172   PetscErrorCode ierr;
173   PetscScalar    ftf, ftg, fty, h;
174 
175   PetscFunctionBegin;
176   ierr   = VecDot(F, F, &ftf);CHKERRQ(ierr);
177   ierr   = VecDot(F, Y, &fty);CHKERRQ(ierr);
178   h      = 1e-5*fty / fty;
179   ierr   = VecCopy(X, W);CHKERRQ(ierr);
180   ierr   = VecAXPY(W, -h, Y);CHKERRQ(ierr);          /* this is arbitrary */
181   ierr   = SNESComputeFunction(snes, W, G);CHKERRQ(ierr);
182   ierr   = VecDot(G, F, &ftg);CHKERRQ(ierr);
183   *ytJtf = (ftg - ftf) / h;
184   PetscFunctionReturn(0);
185 }
186 
187 /*@
188     SNESNCGSetType - Sets the conjugate update type for SNESNCG.
189 
190     Logically Collective on SNES
191 
192     Input Parameters:
193 +   snes - the iterative context
194 -   btype - update type
195 
196     Options Database:
197 .   -snes_ncg_type<prp,fr,hs,dy,cd>
198 
199     Level: intermediate
200 
201     SNESNCGSelectTypes:
202 +   SNES_NCG_FR - Fletcher-Reeves update
203 .   SNES_NCG_PRP - Polak-Ribiere-Polyak update
204 .   SNES_NCG_HS - Hestenes-Steifel update
205 .   SNES_NCG_DY - Dai-Yuan update
206 -   SNES_NCG_CD - Conjugate Descent update
207 
208    Notes:
209    PRP is the default, and the only one that tolerates generalized search directions.
210 
211 .keywords: SNES, SNESNCG, selection, type, set
212 @*/
213 PetscErrorCode SNESNCGSetType(SNES snes, SNESNCGType btype)
214 {
215   PetscErrorCode ierr;
216 
217   PetscFunctionBegin;
218   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
219   ierr = PetscTryMethod(snes,"SNESNCGSetType_C",(SNES,SNESNCGType),(snes,btype));CHKERRQ(ierr);
220   PetscFunctionReturn(0);
221 }
222 
223 PetscErrorCode SNESNCGSetType_NCG(SNES snes, SNESNCGType btype)
224 {
225   SNES_NCG *ncg = (SNES_NCG*)snes->data;
226 
227   PetscFunctionBegin;
228   ncg->type = btype;
229   PetscFunctionReturn(0);
230 }
231 
232 /*
233   SNESSolve_NCG - Solves a nonlinear system with the Nonlinear Conjugate Gradient method.
234 
235   Input Parameters:
236 . snes - the SNES context
237 
238   Output Parameter:
239 . outits - number of iterations until termination
240 
241   Application Interface Routine: SNESSolve()
242 */
243 PetscErrorCode SNESSolve_NCG(SNES snes)
244 {
245   SNES_NCG             *ncg = (SNES_NCG*)snes->data;
246   Vec                  X,dX,lX,F,dXold;
247   PetscReal            fnorm, ynorm, xnorm, beta = 0.0;
248   PetscScalar          dXdotdX, dXolddotdXold, dXdotdXold, lXdotdX, lXdotdXold;
249   PetscInt             maxits, i;
250   PetscErrorCode       ierr;
251   SNESLineSearchReason lsresult = SNES_LINESEARCH_SUCCEEDED;
252   SNESLineSearch       linesearch;
253   SNESConvergedReason  reason;
254 
255   PetscFunctionBegin;
256   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);
257 
258   ierr = PetscCitationsRegister(SNESCitation,&SNEScite);CHKERRQ(ierr);
259   snes->reason = SNES_CONVERGED_ITERATING;
260 
261   maxits = snes->max_its;            /* maximum number of iterations */
262   X      = snes->vec_sol;            /* X^n */
263   dXold  = snes->work[0];            /* The previous iterate of X */
264   dX     = snes->work[1];            /* the preconditioned direction */
265   lX     = snes->vec_sol_update;     /* the conjugate direction */
266   F      = snes->vec_func;           /* residual vector */
267 
268   ierr = SNESGetLineSearch(snes, &linesearch);CHKERRQ(ierr);
269 
270   ierr       = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
271   snes->iter = 0;
272   snes->norm = 0.;
273   ierr       = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
274 
275   /* compute the initial function and preconditioned update dX */
276 
277   if (snes->npc && snes->functype == SNES_FUNCTION_PRECONDITIONED) {
278     ierr = SNESApplyNPC(snes,X,NULL,dX);CHKERRQ(ierr);
279     ierr = SNESGetConvergedReason(snes->npc,&reason);CHKERRQ(ierr);
280     if (reason < 0  && reason != SNES_DIVERGED_MAX_IT) {
281       snes->reason = SNES_DIVERGED_INNER;
282       PetscFunctionReturn(0);
283     }
284     ierr = VecCopy(dX,F);CHKERRQ(ierr);
285     ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr);
286   } else {
287     if (!snes->vec_func_init_set) {
288       ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr);
289     } else snes->vec_func_init_set = PETSC_FALSE;
290 
291     /* convergence test */
292     ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr);
293     SNESCheckFunctionNorm(snes,fnorm);
294     ierr = VecCopy(F,dX);CHKERRQ(ierr);
295   }
296   if (snes->npc) {
297     if (snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
298       ierr = SNESApplyNPC(snes,X,F,dX);CHKERRQ(ierr);
299       ierr = SNESGetConvergedReason(snes->npc,&reason);CHKERRQ(ierr);
300       if (reason < 0  && reason != SNES_DIVERGED_MAX_IT) {
301         snes->reason = SNES_DIVERGED_INNER;
302         PetscFunctionReturn(0);
303       }
304     }
305   }
306   ierr = VecCopy(dX,lX);CHKERRQ(ierr);
307   ierr = VecDot(dX, dX, &dXdotdX);CHKERRQ(ierr);
308 
309 
310   ierr       = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
311   snes->norm = fnorm;
312   ierr       = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
313   ierr       = SNESLogConvergenceHistory(snes,fnorm,0);CHKERRQ(ierr);
314   ierr       = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr);
315 
316   /* test convergence */
317   ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
318   if (snes->reason) PetscFunctionReturn(0);
319 
320   /* Call general purpose update function */
321   if (snes->ops->update) {
322     ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr);
323   }
324 
325   /* first update -- just use the (preconditioned) residual direction for the initial conjugate direction */
326 
327   for (i = 1; i < maxits + 1; i++) {
328     /* some update types require the old update direction or conjugate direction */
329     if (ncg->type != SNES_NCG_FR) {
330       ierr = VecCopy(dX, dXold);CHKERRQ(ierr);
331     }
332     ierr = SNESLineSearchApply(linesearch,X,F,&fnorm,lX);CHKERRQ(ierr);
333     ierr = SNESLineSearchGetReason(linesearch, &lsresult);CHKERRQ(ierr);
334     ierr = SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm);CHKERRQ(ierr);
335     if (lsresult) {
336       if (++snes->numFailures >= snes->maxFailures) {
337         snes->reason = SNES_DIVERGED_LINE_SEARCH;
338         PetscFunctionReturn(0);
339       }
340     }
341     if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
342       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
343       PetscFunctionReturn(0);
344     }
345     /* Monitor convergence */
346     ierr       = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
347     snes->iter = i;
348     snes->norm = fnorm;
349     ierr       = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
350     ierr       = SNESLogConvergenceHistory(snes,snes->norm,0);CHKERRQ(ierr);
351     ierr       = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr);
352 
353     /* Test for convergence */
354     ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
355     if (snes->reason) PetscFunctionReturn(0);
356 
357     /* Call general purpose update function */
358     if (snes->ops->update) {
359       ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr);
360     }
361     if (snes->npc) {
362       if (snes->functype == SNES_FUNCTION_PRECONDITIONED) {
363         ierr = SNESApplyNPC(snes,X,NULL,dX);CHKERRQ(ierr);
364         ierr = SNESGetConvergedReason(snes->npc,&reason);CHKERRQ(ierr);
365         if (reason < 0  && reason != SNES_DIVERGED_MAX_IT) {
366           snes->reason = SNES_DIVERGED_INNER;
367           PetscFunctionReturn(0);
368         }
369         ierr = VecCopy(dX,F);CHKERRQ(ierr);
370       } else {
371         ierr = SNESApplyNPC(snes,X,F,dX);CHKERRQ(ierr);
372         ierr = SNESGetConvergedReason(snes->npc,&reason);CHKERRQ(ierr);
373         if (reason < 0  && reason != SNES_DIVERGED_MAX_IT) {
374           snes->reason = SNES_DIVERGED_INNER;
375           PetscFunctionReturn(0);
376         }
377       }
378     } else {
379       ierr = VecCopy(F,dX);CHKERRQ(ierr);
380     }
381 
382     /* compute the conjugate direction lX = dX + beta*lX with beta = ((dX, dX) / (dX_old, dX_old) (Fletcher-Reeves update)*/
383     switch (ncg->type) {
384     case SNES_NCG_FR: /* Fletcher-Reeves */
385       dXolddotdXold= dXdotdX;
386       ierr         = VecDot(dX, dX, &dXdotdX);CHKERRQ(ierr);
387       beta         = PetscRealPart(dXdotdX / dXolddotdXold);
388       break;
389     case SNES_NCG_PRP: /* Polak-Ribiere-Poylak */
390       dXolddotdXold= dXdotdX;
391       ierr         = VecDotBegin(dX, dX, &dXdotdXold);CHKERRQ(ierr);
392       ierr         = VecDotBegin(dXold, dX, &dXdotdXold);CHKERRQ(ierr);
393       ierr         = VecDotEnd(dX, dX, &dXdotdX);CHKERRQ(ierr);
394       ierr         = VecDotEnd(dXold, dX, &dXdotdXold);CHKERRQ(ierr);
395       beta         = PetscRealPart(((dXdotdX - dXdotdXold) / dXolddotdXold));
396       if (beta < 0.0) beta = 0.0; /* restart */
397       break;
398     case SNES_NCG_HS: /* Hestenes-Stiefel */
399       ierr = VecDotBegin(dX, dX, &dXdotdX);CHKERRQ(ierr);
400       ierr = VecDotBegin(dX, dXold, &dXdotdXold);CHKERRQ(ierr);
401       ierr = VecDotBegin(lX, dX, &lXdotdX);CHKERRQ(ierr);
402       ierr = VecDotBegin(lX, dXold, &lXdotdXold);CHKERRQ(ierr);
403       ierr = VecDotEnd(dX, dX, &dXdotdX);CHKERRQ(ierr);
404       ierr = VecDotEnd(dX, dXold, &dXdotdXold);CHKERRQ(ierr);
405       ierr = VecDotEnd(lX, dX, &lXdotdX);CHKERRQ(ierr);
406       ierr = VecDotEnd(lX, dXold, &lXdotdXold);CHKERRQ(ierr);
407       beta = PetscRealPart((dXdotdX - dXdotdXold) / (lXdotdX - lXdotdXold));
408       break;
409     case SNES_NCG_DY: /* Dai-Yuan */
410       ierr = VecDotBegin(dX, dX, &dXdotdX);CHKERRQ(ierr);
411       ierr = VecDotBegin(lX, dX, &lXdotdX);CHKERRQ(ierr);
412       ierr = VecDotBegin(lX, dXold, &lXdotdXold);CHKERRQ(ierr);
413       ierr = VecDotEnd(dX, dX, &dXdotdX);CHKERRQ(ierr);
414       ierr = VecDotEnd(lX, dX, &lXdotdX);CHKERRQ(ierr);
415       ierr = VecDotEnd(lX, dXold, &lXdotdXold);CHKERRQ(ierr);
416       beta = PetscRealPart(dXdotdX / (lXdotdXold - lXdotdX));CHKERRQ(ierr);
417       break;
418     case SNES_NCG_CD: /* Conjugate Descent */
419       ierr = VecDotBegin(dX, dX, &dXdotdX);CHKERRQ(ierr);
420       ierr = VecDotBegin(lX, dXold, &lXdotdXold);CHKERRQ(ierr);
421       ierr = VecDotEnd(dX, dX, &dXdotdX);CHKERRQ(ierr);
422       ierr = VecDotEnd(lX, dXold, &lXdotdXold);CHKERRQ(ierr);
423       beta = PetscRealPart(dXdotdX / lXdotdXold);CHKERRQ(ierr);
424       break;
425     }
426     if (ncg->monitor) {
427       ierr = PetscViewerASCIIPrintf(ncg->monitor, "beta = %e\n", (double)beta);CHKERRQ(ierr);
428     }
429     ierr = VecAYPX(lX, beta, dX);CHKERRQ(ierr);
430   }
431   ierr = PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", maxits);CHKERRQ(ierr);
432   if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
433   PetscFunctionReturn(0);
434 }
435 
436 
437 
438 /*MC
439   SNESNCG - Nonlinear Conjugate-Gradient method for the solution of general nonlinear systems.
440 
441   Level: beginner
442 
443   Options Database:
444 +   -snes_ncg_type <fr, prp, dy, hs, cd> - Choice of conjugate-gradient update parameter, default is prp.
445 .   -snes_linesearch_type <cp,l2,basic> - Line search type.
446 -   -snes_ncg_monitor - Print relevant information about the ncg iteration.
447 
448    Notes:
449     This solves the nonlinear system of equations F(x) = 0 using the nonlinear generalization of the conjugate
450           gradient method.  This may be used with a nonlinear preconditioner used to pick the new search directions, but otherwise
451           chooses the initial search direction as F(x) for the initial guess x.
452 
453           Only supports left non-linear preconditioning.
454 
455    References:
456 .  1. -  Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu,"Composing Scalable Nonlinear Algebraic Solvers",
457    SIAM Review, 57(4), 2015
458 
459 
460 .seealso:  SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESNEWTONTR, SNESNGMRES, SNESQN
461 M*/
462 PETSC_EXTERN PetscErrorCode SNESCreate_NCG(SNES snes)
463 {
464   PetscErrorCode ierr;
465   SNES_NCG       * neP;
466 
467   PetscFunctionBegin;
468   snes->ops->destroy        = SNESDestroy_NCG;
469   snes->ops->setup          = SNESSetUp_NCG;
470   snes->ops->setfromoptions = SNESSetFromOptions_NCG;
471   snes->ops->view           = SNESView_NCG;
472   snes->ops->solve          = SNESSolve_NCG;
473   snes->ops->reset          = SNESReset_NCG;
474 
475   snes->usesksp = PETSC_FALSE;
476   snes->usesnpc = PETSC_TRUE;
477   snes->npcside = PC_LEFT;
478 
479   snes->alwayscomputesfinalresidual = PETSC_TRUE;
480 
481   if (!snes->tolerancesset) {
482     snes->max_funcs = 30000;
483     snes->max_its   = 10000;
484     snes->stol      = 1e-20;
485   }
486 
487   ierr         = PetscNewLog(snes,&neP);CHKERRQ(ierr);
488   snes->data   = (void*) neP;
489   neP->monitor = NULL;
490   neP->type    = SNES_NCG_PRP;
491   ierr         = PetscObjectComposeFunction((PetscObject)snes,"SNESNCGSetType_C", SNESNCGSetType_NCG);CHKERRQ(ierr);
492   PetscFunctionReturn(0);
493 }
494