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) PetscInfo(snes, "Default or user provided convergence test KSP iterations=%" PetscInt_FMT ", rnorm=%g\n", n, (double)rnorm);
21: /* Determine norm of solution */
22: KSPBuildSolution(ksp, NULL, &x);
23: VecNorm(x, NORM_2, &nrm);
24: if (nrm >= neP->delta) {
25: PetscInfo(snes, "Ending linear iteration early, delta=%g, length=%g\n", (double)neP->delta, (double)nrm);
26: *reason = KSP_CONVERGED_STEP_LENGTH;
27: }
28: return 0;
29: }
31: static PetscErrorCode SNESTR_KSPConverged_Destroy(void *cctx)
32: {
33: SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx *)cctx;
35: (*ctx->convdestroy)(ctx->convctx);
36: PetscFree(ctx);
37: return 0;
38: }
40: /*
41: SNESTR_Converged_Private -test convergence JUST for
42: the trust region tolerance.
43: */
44: static PetscErrorCode SNESTR_Converged_Private(SNES snes, PetscInt it, PetscReal xnorm, PetscReal pnorm, PetscReal fnorm, SNESConvergedReason *reason, void *dummy)
45: {
46: SNES_NEWTONTR *neP = (SNES_NEWTONTR *)snes->data;
48: *reason = SNES_CONVERGED_ITERATING;
49: if (neP->delta < xnorm * snes->deltatol) {
50: PetscInfo(snes, "Converged due to trust region param %g<%g*%g\n", (double)neP->delta, (double)xnorm, (double)snes->deltatol);
51: *reason = SNES_DIVERGED_TR_DELTA;
52: } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
53: PetscInfo(snes, "Exceeded maximum number of function evaluations: %" PetscInt_FMT "\n", snes->max_funcs);
54: *reason = SNES_DIVERGED_FUNCTION_COUNT;
55: }
56: return 0;
57: }
59: /*@C
60: SNESNewtonTRSetPreCheck - Sets a user function that is called before the search step has been determined.
61: Allows the user a chance to change or override the decision of the line search routine.
63: Deprecated use `SNESNEWTONDCTRDC`
65: Logically Collective
67: Input Parameters:
68: + snes - the nonlinear solver object
69: . func - [optional] function evaluation routine, see `SNESNewtonTRPreCheck()` for the calling sequence
70: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
72: Level: intermediate
74: Note:
75: This function is called BEFORE the function evaluation within the `SNESNEWTONTR` solver.
77: .seealso: `SNESNEWTONDCTRDC`, `SNESNEWTONDCTR`, `SNESNewtonTRPreCheck()`, `SNESNewtonTRGetPreCheck()`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()`
78: @*/
79: PetscErrorCode SNESNewtonTRSetPreCheck(SNES snes, PetscErrorCode (*func)(SNES, Vec, Vec, PetscBool *, void *), void *ctx)
80: {
81: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
84: if (func) tr->precheck = func;
85: if (ctx) tr->precheckctx = ctx;
86: return 0;
87: }
89: /*@C
90: SNESNewtonTRGetPreCheck - Gets the pre-check function
92: Deprecated use `SNESNEWTONDCTRDC`
94: Not collective
96: Input Parameter:
97: . snes - the nonlinear solver context
99: Output Parameters:
100: + func - [optional] function evaluation routine, see for the calling sequence `SNESNewtonTRPreCheck()`
101: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
103: Level: intermediate
105: .seealso: `SNESNEWTONDCTRDC`, `SNESNEWTONDCTR`, `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRPreCheck()`
106: @*/
107: PetscErrorCode SNESNewtonTRGetPreCheck(SNES snes, PetscErrorCode (**func)(SNES, Vec, Vec, PetscBool *, void *), void **ctx)
108: {
109: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
112: if (func) *func = tr->precheck;
113: if (ctx) *ctx = tr->precheckctx;
114: return 0;
115: }
117: /*@C
118: SNESNewtonTRSetPostCheck - Sets a user function that is called after the search step has been determined but before the next
119: function evaluation. Allows the user a chance to change or override the decision of the line search routine
121: Deprecated use `SNESNEWTONDCTRDC`
123: Logically Collective
125: Input Parameters:
126: + snes - the nonlinear solver object
127: . func - [optional] function evaluation routine, see `SNESNewtonTRPostCheck()` for the calling sequence
128: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
130: Level: intermediate
132: Note:
133: This function is called BEFORE the function evaluation within the `SNESNEWTONTR` solver while the function set in
134: `SNESLineSearchSetPostCheck()` is called AFTER the function evaluation.
136: .seealso: `SNESNEWTONDCTRDC`, `SNESNEWTONDCTR`, `SNESNewtonTRPostCheck()`, `SNESNewtonTRGetPostCheck()`
137: @*/
138: PetscErrorCode SNESNewtonTRSetPostCheck(SNES snes, PetscErrorCode (*func)(SNES, Vec, Vec, Vec, PetscBool *, PetscBool *, void *), void *ctx)
139: {
140: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
143: if (func) tr->postcheck = func;
144: if (ctx) tr->postcheckctx = ctx;
145: return 0;
146: }
148: /*@C
149: SNESNewtonTRGetPostCheck - Gets the post-check function
151: Deprecated use `SNESNEWTONDCTRDC`
153: Not collective
155: Input Parameter:
156: . snes - the nonlinear solver context
158: Output Parameters:
159: + func - [optional] function evaluation routine, see for the calling sequence `SNESNewtonTRPostCheck()`
160: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
162: Level: intermediate
164: .seealso: `SNESNEWTONDCTRDC`, `SNESNEWTONDCTR`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRPostCheck()`
165: @*/
166: PetscErrorCode SNESNewtonTRGetPostCheck(SNES snes, PetscErrorCode (**func)(SNES, Vec, Vec, Vec, PetscBool *, PetscBool *, void *), void **ctx)
167: {
168: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
171: if (func) *func = tr->postcheck;
172: if (ctx) *ctx = tr->postcheckctx;
173: return 0;
174: }
176: /*@C
177: SNESNewtonTRPreCheck - Called before the step has been determined in `SNESNEWTONTR`
179: Deprecated use `SNESNEWTONDCTRDC`
181: Logically Collective
183: Input Parameters:
184: + snes - the solver
185: . X - The last solution
186: - Y - The step direction
188: Output Parameters:
189: . changed_Y - Indicator that the step direction Y has been changed.
191: Level: developer
193: .seealso: `SNESNEWTONDCTRDC`, `SNESNEWTONDCTR`, `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRGetPreCheck()`
194: @*/
195: static PetscErrorCode SNESNewtonTRPreCheck(SNES snes, Vec X, Vec Y, PetscBool *changed_Y)
196: {
197: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
199: *changed_Y = PETSC_FALSE;
200: if (tr->precheck) {
201: (*tr->precheck)(snes, X, Y, changed_Y, tr->precheckctx);
203: }
204: return 0;
205: }
207: /*@C
208: SNESNewtonTRPostCheck - Called after the step has been determined in `SNESNEWTONTR` but before the function evaluation
210: Deprecated use `SNESNEWTONDCTRDC`
212: Logically Collective
214: Input Parameters:
215: + snes - the solver
216: . X - The last solution
217: . Y - The full step direction
218: - W - The updated solution, W = X - Y
220: Output Parameters:
221: + changed_Y - indicator if step has been changed
222: - changed_W - Indicator if the new candidate solution W has been changed.
224: Note:
225: If Y is changed then W is recomputed as X - Y
227: Level: developer
229: .seealso: `SNESNEWTONDCTRDC`, `SNESNEWTONDCTR`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()`
230: @*/
231: static PetscErrorCode SNESNewtonTRPostCheck(SNES snes, Vec X, Vec Y, Vec W, PetscBool *changed_Y, PetscBool *changed_W)
232: {
233: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
235: *changed_Y = PETSC_FALSE;
236: *changed_W = PETSC_FALSE;
237: if (tr->postcheck) {
238: (*tr->postcheck)(snes, X, Y, W, changed_Y, changed_W, tr->postcheckctx);
241: }
242: return 0;
243: }
245: /*
246: SNESSolve_NEWTONTR - Implements Newton's Method with a very simple trust
247: region approach for solving systems of nonlinear equations.
249: */
250: static PetscErrorCode SNESSolve_NEWTONTR(SNES snes)
251: {
252: SNES_NEWTONTR *neP = (SNES_NEWTONTR *)snes->data;
253: Vec X, F, Y, G, Ytmp, W;
254: PetscInt maxits, i, lits;
255: PetscReal rho, fnorm, gnorm, gpnorm, xnorm = 0, delta, nrm, ynorm, norm1;
256: PetscScalar cnorm;
257: KSP ksp;
258: SNESConvergedReason reason = SNES_CONVERGED_ITERATING;
259: PetscBool breakout = PETSC_FALSE;
260: SNES_TR_KSPConverged_Ctx *ctx;
261: PetscErrorCode (*convtest)(KSP, PetscInt, PetscReal, KSPConvergedReason *, void *), (*convdestroy)(void *);
262: void *convctx;
266: maxits = snes->max_its; /* maximum number of iterations */
267: X = snes->vec_sol; /* solution vector */
268: F = snes->vec_func; /* residual vector */
269: Y = snes->work[0]; /* work vectors */
270: G = snes->work[1];
271: Ytmp = snes->work[2];
272: W = snes->work[3];
274: PetscObjectSAWsTakeAccess((PetscObject)snes);
275: snes->iter = 0;
276: PetscObjectSAWsGrantAccess((PetscObject)snes);
278: /* Set the linear stopping criteria to use the More' trick. */
279: SNESGetKSP(snes, &ksp);
280: KSPGetConvergenceTest(ksp, &convtest, &convctx, &convdestroy);
281: if (convtest != SNESTR_KSPConverged_Private) {
282: PetscNew(&ctx);
283: ctx->snes = snes;
284: KSPGetAndClearConvergenceTest(ksp, &ctx->convtest, &ctx->convctx, &ctx->convdestroy);
285: KSPSetConvergenceTest(ksp, SNESTR_KSPConverged_Private, ctx, SNESTR_KSPConverged_Destroy);
286: PetscInfo(snes, "Using Krylov convergence test SNESTR_KSPConverged_Private\n");
287: }
289: if (!snes->vec_func_init_set) {
290: SNESComputeFunction(snes, X, F); /* F(X) */
291: } else snes->vec_func_init_set = PETSC_FALSE;
293: VecNorm(F, NORM_2, &fnorm); /* fnorm <- || F || */
294: SNESCheckFunctionNorm(snes, fnorm);
295: VecNorm(X, NORM_2, &xnorm); /* xnorm <- || X || */
296: PetscObjectSAWsTakeAccess((PetscObject)snes);
297: snes->norm = fnorm;
298: PetscObjectSAWsGrantAccess((PetscObject)snes);
299: delta = xnorm ? neP->delta0 * xnorm : neP->delta0;
300: neP->delta = delta;
301: SNESLogConvergenceHistory(snes, fnorm, 0);
302: SNESMonitor(snes, 0, fnorm);
304: /* test convergence */
305: PetscUseTypeMethod(snes, converged, snes->iter, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
306: if (snes->reason) return 0;
308: for (i = 0; i < maxits; i++) {
309: /* Call general purpose update function */
310: PetscTryTypeMethod(snes, update, snes->iter);
312: /* Solve J Y = F, where J is Jacobian matrix */
313: SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre);
314: SNESCheckJacobianDomainerror(snes);
315: KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre);
316: KSPSolve(snes->ksp, F, Ytmp);
317: KSPGetIterationNumber(snes->ksp, &lits);
318: snes->linear_its += lits;
320: PetscInfo(snes, "iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits);
321: VecNorm(Ytmp, NORM_2, &nrm);
322: norm1 = nrm;
324: while (1) {
325: PetscBool changed_y;
326: PetscBool changed_w;
327: VecCopy(Ytmp, Y);
328: nrm = norm1;
330: /* Scale Y if need be and predict new value of F norm */
331: if (nrm >= delta) {
332: nrm = delta / nrm;
333: gpnorm = (1.0 - nrm) * fnorm;
334: cnorm = nrm;
335: PetscInfo(snes, "Scaling direction by %g\n", (double)nrm);
336: VecScale(Y, cnorm);
337: nrm = gpnorm;
338: ynorm = delta;
339: } else {
340: gpnorm = 0.0;
341: PetscInfo(snes, "Direction is in Trust Region\n");
342: ynorm = nrm;
343: }
344: /* PreCheck() allows for updates to Y prior to W <- X - Y */
345: SNESNewtonTRPreCheck(snes, X, Y, &changed_y);
346: VecWAXPY(W, -1.0, Y, X); /* W <- X - Y */
347: SNESNewtonTRPostCheck(snes, X, Y, W, &changed_y, &changed_w);
348: if (changed_y) VecWAXPY(W, -1.0, Y, X);
349: VecCopy(Y, snes->vec_sol_update);
350: SNESComputeFunction(snes, W, G); /* F(X-Y) = G */
351: VecNorm(G, NORM_2, &gnorm); /* gnorm <- || g || */
352: SNESCheckFunctionNorm(snes, gnorm);
353: if (fnorm == gpnorm) rho = 0.0;
354: else rho = (fnorm * fnorm - gnorm * gnorm) / (fnorm * fnorm - gpnorm * gpnorm);
356: /* Update size of trust region */
357: if (rho < neP->mu) delta *= neP->delta1;
358: else if (rho < neP->eta) delta *= neP->delta2;
359: else delta *= neP->delta3;
360: PetscInfo(snes, "fnorm=%g, gnorm=%g, ynorm=%g\n", (double)fnorm, (double)gnorm, (double)ynorm);
361: PetscInfo(snes, "gpred=%g, rho=%g, delta=%g\n", (double)gpnorm, (double)rho, (double)delta);
363: neP->delta = delta;
364: if (rho > neP->sigma) break;
365: PetscInfo(snes, "Trying again in smaller region\n");
367: /* check to see if progress is hopeless */
368: neP->itflag = PETSC_FALSE;
369: SNESTR_Converged_Private(snes, snes->iter, xnorm, ynorm, fnorm, &reason, snes->cnvP);
370: if (!reason) PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, fnorm, &reason, snes->cnvP);
371: if (reason == SNES_CONVERGED_SNORM_RELATIVE) reason = SNES_DIVERGED_INNER;
372: if (reason) {
373: /* We're not progressing, so return with the current iterate */
374: SNESMonitor(snes, i + 1, fnorm);
375: breakout = PETSC_TRUE;
376: break;
377: }
378: snes->numFailures++;
379: }
380: if (!breakout) {
381: /* Update function and solution vectors */
382: fnorm = gnorm;
383: VecCopy(G, F);
384: VecCopy(W, X);
385: /* Monitor convergence */
386: PetscObjectSAWsTakeAccess((PetscObject)snes);
387: snes->iter = i + 1;
388: snes->norm = fnorm;
389: snes->xnorm = xnorm;
390: snes->ynorm = ynorm;
391: PetscObjectSAWsGrantAccess((PetscObject)snes);
392: SNESLogConvergenceHistory(snes, snes->norm, lits);
393: SNESMonitor(snes, snes->iter, snes->norm);
394: /* Test for convergence, xnorm = || X || */
395: neP->itflag = PETSC_TRUE;
396: if (snes->ops->converged != SNESConvergedSkip) VecNorm(X, NORM_2, &xnorm);
397: PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, fnorm, &reason, snes->cnvP);
398: if (reason) break;
399: } else break;
400: }
402: if (i == maxits) {
403: PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", maxits);
404: if (!reason) reason = SNES_DIVERGED_MAX_IT;
405: }
406: PetscObjectSAWsTakeAccess((PetscObject)snes);
407: snes->reason = reason;
408: PetscObjectSAWsGrantAccess((PetscObject)snes);
409: if (convtest != SNESTR_KSPConverged_Private) {
410: KSPGetAndClearConvergenceTest(ksp, &ctx->convtest, &ctx->convctx, &ctx->convdestroy);
411: PetscFree(ctx);
412: KSPSetConvergenceTest(ksp, convtest, convctx, convdestroy);
413: }
414: return 0;
415: }
417: static PetscErrorCode SNESSetUp_NEWTONTR(SNES snes)
418: {
419: SNESSetWorkVecs(snes, 4);
420: SNESSetUpMatrices(snes);
421: return 0;
422: }
424: PetscErrorCode SNESReset_NEWTONTR(SNES snes)
425: {
426: return 0;
427: }
429: static PetscErrorCode SNESDestroy_NEWTONTR(SNES snes)
430: {
431: SNESReset_NEWTONTR(snes);
432: PetscFree(snes->data);
433: return 0;
434: }
436: static PetscErrorCode SNESSetFromOptions_NEWTONTR(SNES snes, PetscOptionItems *PetscOptionsObject)
437: {
438: SNES_NEWTONTR *ctx = (SNES_NEWTONTR *)snes->data;
440: PetscOptionsHeadBegin(PetscOptionsObject, "SNES trust region options for nonlinear equations");
441: PetscOptionsReal("-snes_trtol", "Trust region tolerance", "SNESSetTrustRegionTolerance", snes->deltatol, &snes->deltatol, NULL);
442: PetscOptionsReal("-snes_tr_mu", "mu", "None", ctx->mu, &ctx->mu, NULL);
443: PetscOptionsReal("-snes_tr_eta", "eta", "None", ctx->eta, &ctx->eta, NULL);
444: PetscOptionsReal("-snes_tr_sigma", "sigma", "None", ctx->sigma, &ctx->sigma, NULL);
445: PetscOptionsReal("-snes_tr_delta0", "delta0", "None", ctx->delta0, &ctx->delta0, NULL);
446: PetscOptionsReal("-snes_tr_delta1", "delta1", "None", ctx->delta1, &ctx->delta1, NULL);
447: PetscOptionsReal("-snes_tr_delta2", "delta2", "None", ctx->delta2, &ctx->delta2, NULL);
448: PetscOptionsReal("-snes_tr_delta3", "delta3", "None", ctx->delta3, &ctx->delta3, NULL);
449: PetscOptionsHeadEnd();
450: return 0;
451: }
453: static PetscErrorCode SNESView_NEWTONTR(SNES snes, PetscViewer viewer)
454: {
455: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
456: PetscBool iascii;
458: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
459: if (iascii) {
460: PetscViewerASCIIPrintf(viewer, " Trust region tolerance %g (-snes_trtol)\n", (double)snes->deltatol);
461: PetscViewerASCIIPrintf(viewer, " mu=%g, eta=%g, sigma=%g\n", (double)tr->mu, (double)tr->eta, (double)tr->sigma);
462: PetscViewerASCIIPrintf(viewer, " delta0=%g, delta1=%g, delta2=%g, delta3=%g\n", (double)tr->delta0, (double)tr->delta1, (double)tr->delta2, (double)tr->delta3);
463: }
464: return 0;
465: }
467: /*MC
468: SNESNEWTONTR - Newton based nonlinear solver that uses a trust region
470: Deprecated use `SNESNEWTONTRDC`
472: Options Database Keys:
473: + -snes_trtol <tol> - trust region tolerance
474: . -snes_tr_mu <mu> - trust region parameter
475: . -snes_tr_eta <eta> - trust region parameter
476: . -snes_tr_sigma <sigma> - trust region parameter
477: . -snes_tr_delta0 <delta0> - initial size of the trust region is delta0*norm2(x)
478: . -snes_tr_delta1 <delta1> - trust region parameter
479: . -snes_tr_delta2 <delta2> - trust region parameter
480: - -snes_tr_delta3 <delta3> - trust region parameter
482: Reference:
483: . - * "The Minpack Project", by More', Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development
484: of Mathematical Software", Wayne Cowell, editor.
486: Level: intermediate
488: .seealso: `SNESNEWTONTRDC`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESSetTrustRegionTolerance()`
489: M*/
490: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONTR(SNES snes)
491: {
492: SNES_NEWTONTR *neP;
494: snes->ops->setup = SNESSetUp_NEWTONTR;
495: snes->ops->solve = SNESSolve_NEWTONTR;
496: snes->ops->destroy = SNESDestroy_NEWTONTR;
497: snes->ops->setfromoptions = SNESSetFromOptions_NEWTONTR;
498: snes->ops->view = SNESView_NEWTONTR;
499: snes->ops->reset = SNESReset_NEWTONTR;
501: snes->usesksp = PETSC_TRUE;
502: snes->usesnpc = PETSC_FALSE;
504: snes->alwayscomputesfinalresidual = PETSC_TRUE;
506: PetscNew(&neP);
507: snes->data = (void *)neP;
508: neP->mu = 0.25;
509: neP->eta = 0.75;
510: neP->delta = 0.0;
511: neP->delta0 = 0.2;
512: neP->delta1 = 0.3;
513: neP->delta2 = 0.75;
514: neP->delta3 = 2.0;
515: neP->sigma = 0.0001;
516: neP->itflag = PETSC_FALSE;
517: neP->rnorm0 = 0.0;
518: neP->ttol = 0.0;
519: return 0;
520: }