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