Actual source code: tr.c
1: #include <../src/snes/impls/tr/trimpl.h>
3: typedef struct {
4: SNES snes;
5: /* Information on the regular SNES convergence test; which may have been user provided */
6: PetscErrorCode (*convtest)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);
7: PetscErrorCode (*convdestroy)(void*);
8: void *convctx;
9: } SNES_TR_KSPConverged_Ctx;
11: static PetscErrorCode SNESTR_KSPConverged_Private(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *cctx)
12: {
13: SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx*)cctx;
14: SNES snes = ctx->snes;
15: SNES_NEWTONTR *neP = (SNES_NEWTONTR*)snes->data;
16: Vec x;
17: PetscReal nrm;
19: (*ctx->convtest)(ksp,n,rnorm,reason,ctx->convctx);
20: if (*reason) {
21: PetscInfo(snes,"Default or user provided convergence test KSP iterations=%D, rnorm=%g\n",n,(double)rnorm);
22: }
23: /* Determine norm of solution */
24: KSPBuildSolution(ksp,NULL,&x);
25: VecNorm(x,NORM_2,&nrm);
26: if (nrm >= neP->delta) {
27: PetscInfo(snes,"Ending linear iteration early, delta=%g, length=%g\n",(double)neP->delta,(double)nrm);
28: *reason = KSP_CONVERGED_STEP_LENGTH;
29: }
30: return 0;
31: }
33: static PetscErrorCode SNESTR_KSPConverged_Destroy(void *cctx)
34: {
35: SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx*)cctx;
37: (*ctx->convdestroy)(ctx->convctx);
38: PetscFree(ctx);
39: return 0;
40: }
42: /* ---------------------------------------------------------------- */
43: /*
44: SNESTR_Converged_Private -test convergence JUST for
45: the trust region tolerance.
47: */
48: static PetscErrorCode SNESTR_Converged_Private(SNES snes,PetscInt it,PetscReal xnorm,PetscReal pnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy)
49: {
50: SNES_NEWTONTR *neP = (SNES_NEWTONTR*)snes->data;
52: *reason = SNES_CONVERGED_ITERATING;
53: if (neP->delta < xnorm * snes->deltatol) {
54: PetscInfo(snes,"Converged due to trust region param %g<%g*%g\n",(double)neP->delta,(double)xnorm,(double)snes->deltatol);
55: *reason = SNES_DIVERGED_TR_DELTA;
56: } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
57: PetscInfo(snes,"Exceeded maximum number of function evaluations: %D\n",snes->max_funcs);
58: *reason = SNES_DIVERGED_FUNCTION_COUNT;
59: }
60: return 0;
61: }
63: /*@C
64: SNESNewtonTRSetPreCheck - Sets a user function that is called before the search step has been determined.
65: Allows the user a chance to change or override the decision of the line search routine.
67: Logically Collective on snes
69: Input Parameters:
70: + snes - the nonlinear solver object
71: . func - [optional] function evaluation routine, see SNESNewtonTRPreCheck() for the calling sequence
72: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
74: Level: intermediate
76: Note: This function is called BEFORE the function evaluation within the SNESNEWTONTR solver.
78: .seealso: SNESNewtonTRPreCheck(), SNESNewtonTRGetPreCheck(), SNESNewtonTRSetPostCheck(), SNESNewtonTRGetPostCheck()
79: @*/
80: PetscErrorCode SNESNewtonTRSetPreCheck(SNES snes, PetscErrorCode (*func)(SNES,Vec,Vec,PetscBool*,void*),void *ctx)
81: {
82: SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data;
85: if (func) tr->precheck = func;
86: if (ctx) tr->precheckctx = ctx;
87: return 0;
88: }
90: /*@C
91: SNESNewtonTRGetPreCheck - Gets the pre-check function
93: Not collective
95: Input Parameter:
96: . snes - the nonlinear solver context
98: Output Parameters:
99: + func - [optional] function evaluation routine, see for the calling sequence SNESNewtonTRPreCheck()
100: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
102: Level: intermediate
104: .seealso: SNESNewtonTRSetPreCheck(), SNESNewtonTRPreCheck()
105: @*/
106: PetscErrorCode SNESNewtonTRGetPreCheck(SNES snes, PetscErrorCode (**func)(SNES,Vec,Vec,PetscBool*,void*),void **ctx)
107: {
108: SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data;
111: if (func) *func = tr->precheck;
112: if (ctx) *ctx = tr->precheckctx;
113: return 0;
114: }
116: /*@C
117: SNESNewtonTRSetPostCheck - Sets a user function that is called after the search step has been determined but before the next
118: function evaluation. Allows the user a chance to change or override the decision of the line search routine
120: Logically Collective on snes
122: Input Parameters:
123: + snes - the nonlinear solver object
124: . func - [optional] function evaluation routine, see SNESNewtonTRPostCheck() for the calling sequence
125: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
127: Level: intermediate
129: Note: This function is called BEFORE the function evaluation within the SNESNEWTONTR solver while the function set in
130: SNESLineSearchSetPostCheck() is called AFTER the function evaluation.
132: .seealso: SNESNewtonTRPostCheck(), SNESNewtonTRGetPostCheck()
133: @*/
134: PetscErrorCode SNESNewtonTRSetPostCheck(SNES snes,PetscErrorCode (*func)(SNES,Vec,Vec,Vec,PetscBool*,PetscBool*,void*),void *ctx)
135: {
136: SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data;
139: if (func) tr->postcheck = func;
140: if (ctx) tr->postcheckctx = ctx;
141: return 0;
142: }
144: /*@C
145: SNESNewtonTRGetPostCheck - Gets the post-check function
147: Not collective
149: Input Parameter:
150: . snes - the nonlinear solver context
152: Output Parameters:
153: + func - [optional] function evaluation routine, see for the calling sequence SNESNewtonTRPostCheck()
154: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
156: Level: intermediate
158: .seealso: SNESNewtonTRSetPostCheck(), SNESNewtonTRPostCheck()
159: @*/
160: PetscErrorCode SNESNewtonTRGetPostCheck(SNES snes,PetscErrorCode (**func)(SNES,Vec,Vec,Vec,PetscBool*,PetscBool*,void*),void **ctx)
161: {
162: SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data;
165: if (func) *func = tr->postcheck;
166: if (ctx) *ctx = tr->postcheckctx;
167: return 0;
168: }
170: /*@C
171: SNESNewtonTRPreCheck - Called before the step has been determined in SNESNEWTONTR
173: Logically Collective on snes
175: Input Parameters:
176: + snes - the solver
177: . X - The last solution
178: - Y - The step direction
180: Output Parameters:
181: . changed_Y - Indicator that the step direction Y has been changed.
183: Level: developer
185: .seealso: SNESNewtonTRSetPreCheck(), SNESNewtonTRGetPreCheck()
186: @*/
187: static PetscErrorCode SNESNewtonTRPreCheck(SNES snes,Vec X,Vec Y,PetscBool *changed_Y)
188: {
189: SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data;
191: *changed_Y = PETSC_FALSE;
192: if (tr->precheck) {
193: (*tr->precheck)(snes,X,Y,changed_Y,tr->precheckctx);
195: }
196: return 0;
197: }
199: /*@C
200: SNESNewtonTRPostCheck - Called after the step has been determined in SNESNEWTONTR but before the function evaluation
202: Logically Collective on snes
204: Input Parameters:
205: + snes - the solver
206: . X - The last solution
207: . Y - The full step direction
208: - W - The updated solution, W = X - Y
210: Output Parameters:
211: + changed_Y - indicator if step has been changed
212: - changed_W - Indicator if the new candidate solution W has been changed.
214: Notes:
215: If Y is changed then W is recomputed as X - Y
217: Level: developer
219: .seealso: SNESNewtonTRSetPostCheck(), SNESNewtonTRGetPostCheck()
220: @*/
221: static PetscErrorCode SNESNewtonTRPostCheck(SNES snes,Vec X,Vec Y,Vec W,PetscBool *changed_Y,PetscBool *changed_W)
222: {
223: SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data;
225: *changed_Y = PETSC_FALSE;
226: *changed_W = PETSC_FALSE;
227: if (tr->postcheck) {
228: (*tr->postcheck)(snes,X,Y,W,changed_Y,changed_W,tr->postcheckctx);
231: }
232: return 0;
233: }
235: /*
236: SNESSolve_NEWTONTR - Implements Newton's Method with a very simple trust
237: region approach for solving systems of nonlinear equations.
239: */
240: static PetscErrorCode SNESSolve_NEWTONTR(SNES snes)
241: {
242: SNES_NEWTONTR *neP = (SNES_NEWTONTR*)snes->data;
243: Vec X,F,Y,G,Ytmp,W;
244: PetscInt maxits,i,lits;
245: PetscReal rho,fnorm,gnorm,gpnorm,xnorm=0,delta,nrm,ynorm,norm1;
246: PetscScalar cnorm;
247: KSP ksp;
248: SNESConvergedReason reason = SNES_CONVERGED_ITERATING;
249: PetscBool breakout = PETSC_FALSE;
250: SNES_TR_KSPConverged_Ctx *ctx;
251: PetscErrorCode (*convtest)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*),(*convdestroy)(void*);
252: void *convctx;
256: maxits = snes->max_its; /* maximum number of iterations */
257: X = snes->vec_sol; /* solution vector */
258: F = snes->vec_func; /* residual vector */
259: Y = snes->work[0]; /* work vectors */
260: G = snes->work[1];
261: Ytmp = snes->work[2];
262: W = snes->work[3];
264: PetscObjectSAWsTakeAccess((PetscObject)snes);
265: snes->iter = 0;
266: PetscObjectSAWsGrantAccess((PetscObject)snes);
268: /* Set the linear stopping criteria to use the More' trick. */
269: SNESGetKSP(snes,&ksp);
270: KSPGetConvergenceTest(ksp,&convtest,&convctx,&convdestroy);
271: if (convtest != SNESTR_KSPConverged_Private) {
272: PetscNew(&ctx);
273: ctx->snes = snes;
274: KSPGetAndClearConvergenceTest(ksp,&ctx->convtest,&ctx->convctx,&ctx->convdestroy);
275: KSPSetConvergenceTest(ksp,SNESTR_KSPConverged_Private,ctx,SNESTR_KSPConverged_Destroy);
276: PetscInfo(snes,"Using Krylov convergence test SNESTR_KSPConverged_Private\n");
277: }
279: if (!snes->vec_func_init_set) {
280: SNESComputeFunction(snes,X,F); /* F(X) */
281: } else snes->vec_func_init_set = PETSC_FALSE;
283: VecNorm(F,NORM_2,&fnorm); /* fnorm <- || F || */
284: SNESCheckFunctionNorm(snes,fnorm);
285: VecNorm(X,NORM_2,&xnorm); /* xnorm <- || X || */
286: PetscObjectSAWsTakeAccess((PetscObject)snes);
287: snes->norm = fnorm;
288: PetscObjectSAWsGrantAccess((PetscObject)snes);
289: delta = xnorm ? neP->delta0*xnorm : neP->delta0;
290: neP->delta = delta;
291: SNESLogConvergenceHistory(snes,fnorm,0);
292: SNESMonitor(snes,0,fnorm);
294: /* test convergence */
295: (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
296: if (snes->reason) return 0;
298: for (i=0; i<maxits; i++) {
300: /* Call general purpose update function */
301: if (snes->ops->update) {
302: (*snes->ops->update)(snes, snes->iter);
303: }
305: /* Solve J Y = F, where J is Jacobian matrix */
306: SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
307: SNESCheckJacobianDomainerror(snes);
308: KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);
309: KSPSolve(snes->ksp,F,Ytmp);
310: KSPGetIterationNumber(snes->ksp,&lits);
311: snes->linear_its += lits;
313: PetscInfo(snes,"iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n",snes->iter,lits);
314: VecNorm(Ytmp,NORM_2,&nrm);
315: norm1 = nrm;
317: while (1) {
318: PetscBool changed_y;
319: PetscBool changed_w;
320: VecCopy(Ytmp,Y);
321: nrm = norm1;
323: /* Scale Y if need be and predict new value of F norm */
324: if (nrm >= delta) {
325: nrm = delta/nrm;
326: gpnorm = (1.0 - nrm)*fnorm;
327: cnorm = nrm;
328: PetscInfo(snes,"Scaling direction by %g\n",(double)nrm);
329: VecScale(Y,cnorm);
330: nrm = gpnorm;
331: ynorm = delta;
332: } else {
333: gpnorm = 0.0;
334: PetscInfo(snes,"Direction is in Trust Region\n");
335: ynorm = nrm;
336: }
337: /* PreCheck() allows for updates to Y prior to W <- X - Y */
338: SNESNewtonTRPreCheck(snes,X,Y,&changed_y);
339: VecWAXPY(W,-1.0,Y,X); /* W <- X - Y */
340: SNESNewtonTRPostCheck(snes,X,Y,W,&changed_y,&changed_w);
341: if (changed_y) VecWAXPY(W,-1.0,Y,X);
342: VecCopy(Y,snes->vec_sol_update);
343: SNESComputeFunction(snes,W,G); /* F(X-Y) = G */
344: VecNorm(G,NORM_2,&gnorm); /* gnorm <- || g || */
345: SNESCheckFunctionNorm(snes,gnorm);
346: if (fnorm == gpnorm) rho = 0.0;
347: else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm);
349: /* Update size of trust region */
350: if (rho < neP->mu) delta *= neP->delta1;
351: else if (rho < neP->eta) delta *= neP->delta2;
352: else delta *= neP->delta3;
353: PetscInfo(snes,"fnorm=%g, gnorm=%g, ynorm=%g\n",(double)fnorm,(double)gnorm,(double)ynorm);
354: PetscInfo(snes,"gpred=%g, rho=%g, delta=%g\n",(double)gpnorm,(double)rho,(double)delta);
356: neP->delta = delta;
357: if (rho > neP->sigma) break;
358: PetscInfo(snes,"Trying again in smaller region\n");
360: /* check to see if progress is hopeless */
361: neP->itflag = PETSC_FALSE;
362: SNESTR_Converged_Private(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);
363: if (!reason) (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);
364: if (reason == SNES_CONVERGED_SNORM_RELATIVE) reason = SNES_DIVERGED_INNER;
365: if (reason) {
366: /* We're not progressing, so return with the current iterate */
367: SNESMonitor(snes,i+1,fnorm);
368: breakout = PETSC_TRUE;
369: break;
370: }
371: snes->numFailures++;
372: }
373: if (!breakout) {
374: /* Update function and solution vectors */
375: fnorm = gnorm;
376: VecCopy(G,F);
377: VecCopy(W,X);
378: /* Monitor convergence */
379: PetscObjectSAWsTakeAccess((PetscObject)snes);
380: snes->iter = i+1;
381: snes->norm = fnorm;
382: snes->xnorm = xnorm;
383: snes->ynorm = ynorm;
384: PetscObjectSAWsGrantAccess((PetscObject)snes);
385: SNESLogConvergenceHistory(snes,snes->norm,lits);
386: SNESMonitor(snes,snes->iter,snes->norm);
387: /* Test for convergence, xnorm = || X || */
388: neP->itflag = PETSC_TRUE;
389: if (snes->ops->converged != SNESConvergedSkip) VecNorm(X,NORM_2,&xnorm);
390: (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);
391: if (reason) break;
392: } else break;
393: }
395: if (i == maxits) {
396: PetscInfo(snes,"Maximum number of iterations has been reached: %" PetscInt_FMT "\n",maxits);
397: if (!reason) reason = SNES_DIVERGED_MAX_IT;
398: }
399: PetscObjectSAWsTakeAccess((PetscObject)snes);
400: snes->reason = reason;
401: PetscObjectSAWsGrantAccess((PetscObject)snes);
402: if (convtest != SNESTR_KSPConverged_Private) {
403: KSPGetAndClearConvergenceTest(ksp,&ctx->convtest,&ctx->convctx,&ctx->convdestroy);
404: PetscFree(ctx);
405: KSPSetConvergenceTest(ksp,convtest,convctx,convdestroy);
406: }
407: return 0;
408: }
410: /*------------------------------------------------------------*/
411: static PetscErrorCode SNESSetUp_NEWTONTR(SNES snes)
412: {
413: SNESSetWorkVecs(snes,4);
414: SNESSetUpMatrices(snes);
415: return 0;
416: }
418: PetscErrorCode SNESReset_NEWTONTR(SNES snes)
419: {
420: return 0;
421: }
423: static PetscErrorCode SNESDestroy_NEWTONTR(SNES snes)
424: {
425: SNESReset_NEWTONTR(snes);
426: PetscFree(snes->data);
427: return 0;
428: }
429: /*------------------------------------------------------------*/
431: static PetscErrorCode SNESSetFromOptions_NEWTONTR(PetscOptionItems *PetscOptionsObject,SNES snes)
432: {
433: SNES_NEWTONTR *ctx = (SNES_NEWTONTR*)snes->data;
435: PetscOptionsHead(PetscOptionsObject,"SNES trust region options for nonlinear equations");
436: PetscOptionsReal("-snes_trtol","Trust region tolerance","SNESSetTrustRegionTolerance",snes->deltatol,&snes->deltatol,NULL);
437: PetscOptionsReal("-snes_tr_mu","mu","None",ctx->mu,&ctx->mu,NULL);
438: PetscOptionsReal("-snes_tr_eta","eta","None",ctx->eta,&ctx->eta,NULL);
439: PetscOptionsReal("-snes_tr_sigma","sigma","None",ctx->sigma,&ctx->sigma,NULL);
440: PetscOptionsReal("-snes_tr_delta0","delta0","None",ctx->delta0,&ctx->delta0,NULL);
441: PetscOptionsReal("-snes_tr_delta1","delta1","None",ctx->delta1,&ctx->delta1,NULL);
442: PetscOptionsReal("-snes_tr_delta2","delta2","None",ctx->delta2,&ctx->delta2,NULL);
443: PetscOptionsReal("-snes_tr_delta3","delta3","None",ctx->delta3,&ctx->delta3,NULL);
444: PetscOptionsTail();
445: return 0;
446: }
448: static PetscErrorCode SNESView_NEWTONTR(SNES snes,PetscViewer viewer)
449: {
450: SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data;
451: PetscBool iascii;
453: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
454: if (iascii) {
455: PetscViewerASCIIPrintf(viewer," Trust region tolerance (-snes_trtol)\n",(double)snes->deltatol);
456: PetscViewerASCIIPrintf(viewer," mu=%g, eta=%g, sigma=%g\n",(double)tr->mu,(double)tr->eta,(double)tr->sigma);
457: PetscViewerASCIIPrintf(viewer," delta0=%g, delta1=%g, delta2=%g, delta3=%g\n",(double)tr->delta0,(double)tr->delta1,(double)tr->delta2,(double)tr->delta3);
458: }
459: return 0;
460: }
461: /* ------------------------------------------------------------ */
462: /*MC
463: SNESNEWTONTR - Newton based nonlinear solver that uses a trust region
465: Options Database:
466: + -snes_trtol <tol> - trust region tolerance
467: . -snes_tr_mu <mu> - trust region parameter
468: . -snes_tr_eta <eta> - trust region parameter
469: . -snes_tr_sigma <sigma> - trust region parameter
470: . -snes_tr_delta0 <delta0> - initial size of the trust region is delta0*norm2(x)
471: . -snes_tr_delta1 <delta1> - trust region parameter
472: . -snes_tr_delta2 <delta2> - trust region parameter
473: - -snes_tr_delta3 <delta3> - trust region parameter
475: The basic algorithm is taken from "The Minpack Project", by More',
476: Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development
477: of Mathematical Software", Wayne Cowell, editor.
479: Level: intermediate
481: .seealso: SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESSetTrustRegionTolerance()
483: M*/
484: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONTR(SNES snes)
485: {
486: SNES_NEWTONTR *neP;
488: snes->ops->setup = SNESSetUp_NEWTONTR;
489: snes->ops->solve = SNESSolve_NEWTONTR;
490: snes->ops->destroy = SNESDestroy_NEWTONTR;
491: snes->ops->setfromoptions = SNESSetFromOptions_NEWTONTR;
492: snes->ops->view = SNESView_NEWTONTR;
493: snes->ops->reset = SNESReset_NEWTONTR;
495: snes->usesksp = PETSC_TRUE;
496: snes->usesnpc = PETSC_FALSE;
498: snes->alwayscomputesfinalresidual = PETSC_TRUE;
500: PetscNewLog(snes,&neP);
501: snes->data = (void*)neP;
502: neP->mu = 0.25;
503: neP->eta = 0.75;
504: neP->delta = 0.0;
505: neP->delta0 = 0.2;
506: neP->delta1 = 0.3;
507: neP->delta2 = 0.75;
508: neP->delta3 = 2.0;
509: neP->sigma = 0.0001;
510: neP->itflag = PETSC_FALSE;
511: neP->rnorm0 = 0.0;
512: neP->ttol = 0.0;
513: return 0;
514: }