Actual source code: tr.c
1: #include <../src/snes/impls/tr/trimpl.h>
3: typedef struct {
4: SNES snes;
5: PetscErrorCode (*convtest)(KSP, PetscInt, PetscReal, KSPConvergedReason *, void *);
6: PetscErrorCode (*convdestroy)(void *);
7: void *convctx;
8: } SNES_TR_KSPConverged_Ctx;
10: const char *const SNESNewtonTRFallbackTypes[] = {"NEWTON", "CAUCHY", "DOGLEG", "SNESNewtonTRFallbackType", "SNES_TR_FALLBACK_", NULL};
12: static PetscErrorCode SNESTR_KSPConverged_Private(KSP ksp, PetscInt n, PetscReal rnorm, KSPConvergedReason *reason, void *cctx)
13: {
14: SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx *)cctx;
15: SNES snes = ctx->snes;
16: SNES_NEWTONTR *neP = (SNES_NEWTONTR *)snes->data;
17: Vec x;
18: PetscReal nrm;
20: PetscFunctionBegin;
21: PetscCall((*ctx->convtest)(ksp, n, rnorm, reason, ctx->convctx));
22: if (*reason) PetscCall(PetscInfo(snes, "Default or user provided convergence test KSP iterations=%" PetscInt_FMT ", rnorm=%g\n", n, (double)rnorm));
23: /* Determine norm of solution */
24: PetscCall(KSPBuildSolution(ksp, NULL, &x));
25: PetscCall(VecNorm(x, NORM_2, &nrm));
26: if (nrm >= neP->delta) {
27: PetscCall(PetscInfo(snes, "Ending linear iteration early, delta=%g, length=%g\n", (double)neP->delta, (double)nrm));
28: *reason = KSP_CONVERGED_STEP_LENGTH;
29: }
30: PetscFunctionReturn(PETSC_SUCCESS);
31: }
33: static PetscErrorCode SNESTR_KSPConverged_Destroy(void *cctx)
34: {
35: SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx *)cctx;
37: PetscFunctionBegin;
38: PetscCall((*ctx->convdestroy)(ctx->convctx));
39: PetscCall(PetscFree(ctx));
40: PetscFunctionReturn(PETSC_SUCCESS);
41: }
43: static PetscErrorCode SNESTR_Converged_Private(SNES snes, PetscInt it, PetscReal xnorm, PetscReal pnorm, PetscReal fnorm, SNESConvergedReason *reason, void *dummy)
44: {
45: SNES_NEWTONTR *neP = (SNES_NEWTONTR *)snes->data;
47: PetscFunctionBegin;
48: *reason = SNES_CONVERGED_ITERATING;
49: if (neP->delta < snes->deltatol) {
50: PetscCall(PetscInfo(snes, "Diverged due to too small a trust region %g<%g\n", (double)neP->delta, (double)snes->deltatol));
51: *reason = SNES_DIVERGED_TR_DELTA;
52: } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
53: PetscCall(PetscInfo(snes, "Exceeded maximum number of function evaluations: %" PetscInt_FMT "\n", snes->max_funcs));
54: *reason = SNES_DIVERGED_FUNCTION_COUNT;
55: }
56: PetscFunctionReturn(PETSC_SUCCESS);
57: }
59: /*@
60: SNESNewtonTRSetFallbackType - Set the type of fallback to use if the solution of the trust region subproblem is outside the radius
62: Input Parameters:
63: + snes - the nonlinear solver object
64: - ftype - the fallback type, see `SNESNewtonTRFallbackType`
66: Level: intermediate
68: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRPreCheck()`, `SNESNewtonTRGetPreCheck()`, `SNESNewtonTRSetPreCheck()`,
69: `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()`
70: @*/
71: PetscErrorCode SNESNewtonTRSetFallbackType(SNES snes, SNESNewtonTRFallbackType ftype)
72: {
73: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
74: PetscBool flg;
76: PetscFunctionBegin;
79: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
80: if (flg) tr->fallback = ftype;
81: PetscFunctionReturn(PETSC_SUCCESS);
82: }
84: /*@C
85: SNESNewtonTRSetPreCheck - Sets a user function that is called before the search step has been determined.
86: Allows the user a chance to change or override the trust region decision.
88: Logically Collective
90: Input Parameters:
91: + snes - the nonlinear solver object
92: . func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRPreCheck()`
93: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
95: Level: deprecated (since 3.19)
97: Note:
98: This function is called BEFORE the function evaluation within the solver.
100: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRPreCheck()`, `SNESNewtonTRGetPreCheck()`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()`,
101: @*/
102: PetscErrorCode SNESNewtonTRSetPreCheck(SNES snes, PetscErrorCode (*func)(SNES, Vec, Vec, PetscBool *, void *), void *ctx)
103: {
104: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
105: PetscBool flg;
107: PetscFunctionBegin;
109: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
110: if (flg) {
111: if (func) tr->precheck = func;
112: if (ctx) tr->precheckctx = ctx;
113: }
114: PetscFunctionReturn(PETSC_SUCCESS);
115: }
117: /*@C
118: SNESNewtonTRGetPreCheck - Gets the pre-check function
120: Deprecated use `SNESNEWTONDCTRDC`
122: Not Collective
124: Input Parameter:
125: . snes - the nonlinear solver context
127: Output Parameters:
128: + func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRPreCheck()`
129: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
131: Level: deprecated (since 3.19)
133: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRPreCheck()`
134: @*/
135: PetscErrorCode SNESNewtonTRGetPreCheck(SNES snes, PetscErrorCode (**func)(SNES, Vec, Vec, PetscBool *, void *), void **ctx)
136: {
137: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
138: PetscBool flg;
140: PetscFunctionBegin;
142: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
143: PetscAssert(flg, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)snes)->type_name);
144: if (func) *func = tr->precheck;
145: if (ctx) *ctx = tr->precheckctx;
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: /*@C
150: SNESNewtonTRSetPostCheck - Sets a user function that is called after the search step has been determined but before the next
151: function evaluation. Allows the user a chance to change or override the internal decision of the solver
153: Logically Collective
155: Input Parameters:
156: + snes - the nonlinear solver object
157: . func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRPostCheck()`
158: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
160: Level: deprecated (since 3.19)
162: Note:
163: This function is called BEFORE the function evaluation within the solver while the function set in
164: `SNESLineSearchSetPostCheck()` is called AFTER the function evaluation.
166: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRPostCheck()`, `SNESNewtonTRGetPostCheck()`, `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRGetPreCheck()`
167: @*/
168: PetscErrorCode SNESNewtonTRSetPostCheck(SNES snes, PetscErrorCode (*func)(SNES, Vec, Vec, Vec, PetscBool *, PetscBool *, void *), void *ctx)
169: {
170: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
171: PetscBool flg;
173: PetscFunctionBegin;
175: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
176: if (flg) {
177: if (func) tr->postcheck = func;
178: if (ctx) tr->postcheckctx = ctx;
179: }
180: PetscFunctionReturn(PETSC_SUCCESS);
181: }
183: /*@C
184: SNESNewtonTRGetPostCheck - Gets the post-check function
186: Not Collective
188: Input Parameter:
189: . snes - the nonlinear solver context
191: Output Parameters:
192: + func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRPostCheck()`
193: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
195: Level: intermediate
197: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRPostCheck()`
198: @*/
199: PetscErrorCode SNESNewtonTRGetPostCheck(SNES snes, PetscErrorCode (**func)(SNES, Vec, Vec, Vec, PetscBool *, PetscBool *, void *), void **ctx)
200: {
201: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
202: PetscBool flg;
204: PetscFunctionBegin;
206: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
207: PetscAssert(flg, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)snes)->type_name);
208: if (func) *func = tr->postcheck;
209: if (ctx) *ctx = tr->postcheckctx;
210: PetscFunctionReturn(PETSC_SUCCESS);
211: }
213: /*@C
214: SNESNewtonTRPreCheck - Runs the precheck routine
216: Logically Collective
218: Input Parameters:
219: + snes - the solver
220: . X - The last solution
221: - Y - The step direction
223: Output Parameter:
224: . changed_Y - Indicator that the step direction `Y` has been changed.
226: Level: intermediate
228: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRGetPreCheck()`, `SNESNewtonTRPostCheck()`
229: @*/
230: PetscErrorCode SNESNewtonTRPreCheck(SNES snes, Vec X, Vec Y, PetscBool *changed_Y)
231: {
232: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
233: PetscBool flg;
235: PetscFunctionBegin;
237: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
238: PetscAssert(flg, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)snes)->type_name);
239: *changed_Y = PETSC_FALSE;
240: if (tr->precheck) {
241: PetscCall((*tr->precheck)(snes, X, Y, changed_Y, tr->precheckctx));
243: }
244: PetscFunctionReturn(PETSC_SUCCESS);
245: }
247: /*@C
248: SNESNewtonTRPostCheck - Runs the postcheck routine
250: Logically Collective
252: Input Parameters:
253: + snes - the solver
254: . X - The last solution
255: . Y - The full step direction
256: - W - The updated solution, W = X - Y
258: Output Parameters:
259: + changed_Y - indicator if step has been changed
260: - changed_W - Indicator if the new candidate solution W has been changed.
262: Note:
263: If Y is changed then W is recomputed as X - Y
265: Level: intermediate
267: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()`, `SNESNewtonTRPreCheck()`
268: @*/
269: PetscErrorCode SNESNewtonTRPostCheck(SNES snes, Vec X, Vec Y, Vec W, PetscBool *changed_Y, PetscBool *changed_W)
270: {
271: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
272: PetscBool flg;
274: PetscFunctionBegin;
276: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
277: PetscAssert(flg, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)snes)->type_name);
278: *changed_Y = PETSC_FALSE;
279: *changed_W = PETSC_FALSE;
280: if (tr->postcheck) {
281: PetscCall((*tr->postcheck)(snes, X, Y, W, changed_Y, changed_W, tr->postcheckctx));
284: }
285: PetscFunctionReturn(PETSC_SUCCESS);
286: }
288: static inline void PetscQuadraticRoots(PetscReal a, PetscReal b, PetscReal c, PetscReal *xm, PetscReal *xp)
289: {
290: PetscReal temp = -0.5 * (b + PetscCopysignReal(1.0, b) * PetscSqrtReal(b * b - 4 * a * c));
291: PetscReal x1 = temp / a;
292: PetscReal x2 = c / temp;
293: *xm = PetscMin(x1, x2);
294: *xp = PetscMax(x1, x2);
295: }
297: /* Computes the quadratic model difference */
298: static PetscErrorCode SNESNewtonTRQuadraticDelta(SNES snes, PetscBool has_objective, Vec Y, Vec GradF, Vec W, PetscReal *yTHy, PetscReal *gTy, PetscReal *deltaqm)
299: {
300: PetscFunctionBegin;
301: PetscCall(MatMult(snes->jacobian, Y, W));
302: if (has_objective) PetscCall(VecDotRealPart(Y, W, yTHy));
303: else PetscCall(VecDotRealPart(W, W, yTHy)); /* Gauss-Newton approximation J^t * J */
304: PetscCall(VecDotRealPart(GradF, Y, gTy));
305: *deltaqm = -(-(*gTy) + 0.5 * (*yTHy)); /* difference in quadratic model, -gTy because SNES solves it this way */
306: PetscFunctionReturn(PETSC_SUCCESS);
307: }
309: /*
310: SNESSolve_NEWTONTR - Implements Newton's Method with trust-region subproblem and adds dogleg Cauchy
311: (Steepest Descent direction) step and direction if the trust region is not satisfied for solving system of
312: nonlinear equations
314: */
315: static PetscErrorCode SNESSolve_NEWTONTR(SNES snes)
316: {
317: SNES_NEWTONTR *neP = (SNES_NEWTONTR *)snes->data;
318: Vec X, F, Y, G, W, GradF, YU;
319: PetscInt maxits, lits;
320: PetscReal rho, fnorm, gnorm = 0.0, xnorm = 0.0, delta, ynorm;
321: PetscReal deltaM, fk, fkp1, deltaqm, gTy, yTHy;
322: PetscReal auk, gfnorm, ycnorm, gTBg, objmin = 0.0;
323: KSP ksp;
324: PetscBool already_done = PETSC_FALSE;
325: PetscBool clear_converged_test, rho_satisfied, has_objective;
326: SNES_TR_KSPConverged_Ctx *ctx;
327: void *convctx;
328: PetscErrorCode (*convtest)(KSP, PetscInt, PetscReal, KSPConvergedReason *, void *), (*convdestroy)(void *);
329: PetscErrorCode (*objective)(SNES, Vec, PetscReal *, void *);
331: PetscFunctionBegin;
332: PetscCall(SNESGetObjective(snes, &objective, NULL));
333: has_objective = objective ? PETSC_TRUE : PETSC_FALSE;
335: maxits = snes->max_its; /* maximum number of iterations */
336: X = snes->vec_sol; /* solution vector */
337: F = snes->vec_func; /* residual vector */
338: Y = snes->vec_sol_update; /* update vector */
339: G = snes->work[0]; /* updated residual */
340: W = snes->work[1]; /* temporary vector */
341: GradF = !has_objective ? snes->work[2] : snes->vec_func; /* grad f = J^T F */
342: YU = snes->work[3]; /* work vector for dogleg method */
344: PetscCheck(!snes->xl && !snes->xu && !snes->ops->computevariablebounds, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);
346: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
347: snes->iter = 0;
348: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
350: /* Set the linear stopping criteria to use the More' trick if needed */
351: clear_converged_test = PETSC_FALSE;
352: PetscCall(SNESGetKSP(snes, &ksp));
353: PetscCall(KSPGetConvergenceTest(ksp, &convtest, &convctx, &convdestroy));
354: if (convtest != SNESTR_KSPConverged_Private) {
355: clear_converged_test = PETSC_TRUE;
356: PetscCall(PetscNew(&ctx));
357: ctx->snes = snes;
358: PetscCall(KSPGetAndClearConvergenceTest(ksp, &ctx->convtest, &ctx->convctx, &ctx->convdestroy));
359: PetscCall(KSPSetConvergenceTest(ksp, SNESTR_KSPConverged_Private, ctx, SNESTR_KSPConverged_Destroy));
360: PetscCall(PetscInfo(snes, "Using Krylov convergence test SNESTR_KSPConverged_Private\n"));
361: }
363: if (!snes->vec_func_init_set) {
364: PetscCall(SNESComputeFunction(snes, X, F)); /* F(X) */
365: } else snes->vec_func_init_set = PETSC_FALSE;
367: PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- || F || */
368: SNESCheckFunctionNorm(snes, fnorm);
369: PetscCall(VecNorm(X, NORM_2, &xnorm)); /* xnorm <- || X || */
371: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
372: snes->norm = fnorm;
373: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
374: delta = neP->delta0;
375: deltaM = neP->deltaM;
376: neP->delta = delta;
377: PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));
379: /* test convergence */
380: rho_satisfied = PETSC_FALSE;
381: PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm));
382: PetscCall(SNESMonitor(snes, 0, fnorm));
383: if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
385: if (has_objective) PetscCall(SNESComputeObjective(snes, X, &fk));
386: else fk = 0.5 * PetscSqr(fnorm); /* obj(x) = 0.5 * ||F(x)||^2 */
388: while (snes->iter < maxits) {
389: PetscBool changed_y;
390: PetscBool changed_w;
392: /* calculating Jacobian and GradF of minimization function only once */
393: if (!already_done) {
394: /* Call general purpose update function */
395: PetscTryTypeMethod(snes, update, snes->iter);
397: /* apply the nonlinear preconditioner */
398: if (snes->npc && snes->npcside == PC_RIGHT) {
399: SNESConvergedReason reason;
401: PetscCall(SNESSetInitialFunction(snes->npc, F));
402: PetscCall(PetscLogEventBegin(SNES_NPCSolve, snes->npc, X, snes->vec_rhs, 0));
403: PetscCall(SNESSolve(snes->npc, snes->vec_rhs, X));
404: PetscCall(PetscLogEventEnd(SNES_NPCSolve, snes->npc, X, snes->vec_rhs, 0));
405: PetscCall(SNESGetConvergedReason(snes->npc, &reason));
406: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
407: snes->reason = SNES_DIVERGED_INNER;
408: PetscFunctionReturn(PETSC_SUCCESS);
409: }
410: // XXX
411: PetscCall(SNESGetNPCFunction(snes, F, &fnorm));
412: if (has_objective) PetscCall(SNESComputeObjective(snes, X, &fk));
413: else fk = 0.5 * PetscSqr(fnorm); /* obj(x) = 0.5 * ||F(x)||^2 */
414: // XXX
415: } else if (snes->ops->update) { /* if update is present, recompute objective function and function norm */
416: PetscCall(SNESComputeFunction(snes, X, F));
417: PetscCall(VecNorm(F, NORM_2, &fnorm));
418: if (has_objective) PetscCall(SNESComputeObjective(snes, X, &fk));
419: else fk = 0.5 * PetscSqr(fnorm); /* obj(x) = 0.5 * ||F(x)||^2 */
420: }
422: /* Jacobian */
423: PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
424: SNESCheckJacobianDomainerror(snes);
426: /* GradF */
427: if (has_objective) gfnorm = fnorm;
428: else {
429: PetscCall(MatMultTranspose(snes->jacobian, F, GradF)); /* grad f = J^T F */
430: PetscCall(VecNorm(GradF, NORM_2, &gfnorm));
431: }
432: }
433: already_done = PETSC_TRUE;
435: /* solve trust-region subproblem */
437: /* sufficient decrease (see 6.3.27 in Conn, Gould, Toint "Trust Region Methods")
438: This is actually more relaxed, since they use include gnorm/beta_k, with
439: beta_k the largest eigenvalue of the Hessian */
440: objmin = -neP->kmdc * gnorm * PetscMin(gnorm, delta);
441: PetscCall(KSPCGSetObjectiveTarget(snes->ksp, objmin));
443: /* don't specify radius if not looking for Newton step only */
444: PetscCall(KSPCGSetRadius(snes->ksp, neP->fallback == SNES_TR_FALLBACK_NEWTON ? delta : 0.0));
446: PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre));
447: PetscCall(KSPSolve(snes->ksp, F, Y));
448: SNESCheckKSPSolve(snes);
449: PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
450: PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits));
452: /* decide what to do when the update is outside of trust region */
453: PetscCall(VecNorm(Y, NORM_2, &ynorm));
454: if (ynorm > delta || ynorm == 0.0) {
455: SNESNewtonTRFallbackType fallback = ynorm > 0.0 ? neP->fallback : SNES_TR_FALLBACK_CAUCHY;
457: switch (fallback) {
458: case SNES_TR_FALLBACK_NEWTON:
459: auk = delta / ynorm;
460: PetscCall(VecScale(Y, auk));
461: PetscCall(PetscInfo(snes, "SN evaluated. delta: %g, ynorm: %g\n", (double)delta, (double)ynorm));
462: break;
463: case SNES_TR_FALLBACK_CAUCHY:
464: case SNES_TR_FALLBACK_DOGLEG:
465: PetscCall(MatMult(snes->jacobian, GradF, W));
466: if (has_objective) PetscCall(VecDotRealPart(GradF, W, &gTBg));
467: else PetscCall(VecDotRealPart(W, W, &gTBg)); /* B = J^t * J */
468: /* Eqs 4.7 and 4.8 in Nocedal and Wright */
469: auk = delta / gfnorm;
470: if (gTBg > 0.0) auk *= PetscMin(gfnorm * gfnorm * gfnorm / (delta * gTBg), 1);
471: ycnorm = auk * gfnorm;
472: if (fallback == SNES_TR_FALLBACK_CAUCHY || gTBg <= 0.0) {
473: /* Cauchy solution */
474: PetscCall(VecAXPBY(Y, auk, 0.0, GradF));
475: PetscCall(PetscInfo(snes, "CP evaluated. delta: %g, ynorm: %g, ycnorm: %g, gTBg: %g\n", (double)delta, (double)ynorm, (double)ycnorm, (double)gTBg));
476: } else { /* take linear combination of Cauchy and Newton direction and step */
477: PetscReal c0, c1, c2, tau = 0.0, tpos, tneg;
478: PetscBool noroots;
480: auk = gfnorm * gfnorm / gTBg;
481: PetscCall(VecAXPBY(YU, auk, 0.0, GradF));
482: PetscCall(VecAXPY(Y, -1.0, YU));
483: PetscCall(VecNorm(Y, NORM_2, &c0));
484: PetscCall(VecDotRealPart(YU, Y, &c1));
485: c0 = PetscSqr(c0);
486: c2 = PetscSqr(ycnorm) - PetscSqr(delta);
487: PetscQuadraticRoots(c0, c1, c2, &tneg, &tpos);
489: noroots = PetscIsInfOrNanReal(tneg);
490: if (noroots) { /* No roots, select Cauchy point */
491: auk = delta / gfnorm;
492: auk *= PetscMin(gfnorm * gfnorm * gfnorm / (delta * gTBg), 1);
493: PetscCall(VecAXPBY(Y, auk, 0.0, GradF));
494: } else { /* Here roots corresponds to tau-1 in Nocedal and Wright */
495: tpos += 1.0;
496: tneg += 1.0;
497: tau = PetscClipInterval(tpos, 0.0, 2.0); /* clip to tau [0,2] */
498: if (tau < 1.0) {
499: PetscCall(VecAXPBY(Y, tau, 0.0, YU));
500: } else {
501: PetscCall(VecAXPBY(Y, 1.0, tau - 1, YU));
502: }
503: }
504: PetscCall(VecNorm(Y, NORM_2, &c0)); /* this norm will be cached and reused later */
505: PetscCall(PetscInfo(snes, "%s evaluated. roots: (%g, %g), tau %g, ynorm: %g, ycnorm: %g, ydlnorm %g, gTBg: %g\n", noroots ? "CP" : "DL", (double)tneg, (double)tpos, (double)tau, (double)ynorm, (double)ycnorm, (double)c0, (double)gTBg));
506: }
507: break;
508: default:
509: SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "Unknown fallback mode");
510: break;
511: }
512: }
514: /* Evaluate the solution to meet the improvement ratio criteria */
516: /* compute the final ynorm */
517: PetscCall(SNESNewtonTRPreCheck(snes, X, Y, &changed_y));
518: PetscCall(VecNorm(Y, NORM_2, &ynorm));
520: /* compute the quadratic model difference */
521: PetscCall(SNESNewtonTRQuadraticDelta(snes, has_objective, Y, GradF, W, &yTHy, &gTy, &deltaqm));
523: /* update */
524: PetscCall(VecWAXPY(W, -1.0, Y, X)); /* Xkp1 */
525: PetscCall(SNESNewtonTRPostCheck(snes, X, Y, W, &changed_y, &changed_w));
526: if (changed_y) {
527: /* Need to recompute the quadratic model difference */
528: PetscCall(SNESNewtonTRQuadraticDelta(snes, has_objective, Y, GradF, YU, &yTHy, &gTy, &deltaqm));
529: /* User changed Y but not W */
530: if (!changed_w) PetscCall(VecWAXPY(W, -1.0, Y, X));
531: }
533: /* Compute new objective function */
534: PetscCall(SNESComputeFunction(snes, W, G)); /* F(Xkp1) = G */
535: PetscCall(VecNorm(G, NORM_2, &gnorm));
536: if (has_objective) PetscCall(SNESComputeObjective(snes, W, &fkp1));
537: else fkp1 = 0.5 * PetscSqr(gnorm);
538: SNESCheckFunctionNorm(snes, fkp1);
540: if (deltaqm > 0.0) rho = (fk - fkp1) / deltaqm; /* actual improvement over predicted improvement */
541: else rho = -1.0; /* no reduction in quadratic model, step must be rejected */
542: PetscCall(PetscInfo(snes, "rho=%g, delta=%g, fk=%g, fkp1=%g, deltaqm=%g, gTy=%g, yTHy=%g\n", (double)rho, (double)delta, (double)fk, (double)fkp1, (double)deltaqm, (double)gTy, (double)yTHy));
544: if (rho < neP->eta2) delta *= neP->t1; /* shrink the region */
545: else if (rho > neP->eta3) delta *= neP->t2; /* expand the region */
546: delta = PetscMin(delta, deltaM); /* but not greater than deltaM */
548: neP->delta = delta;
549: if (rho >= neP->eta1) {
550: rho_satisfied = PETSC_TRUE;
551: } else {
552: rho_satisfied = PETSC_FALSE;
553: PetscCall(PetscInfo(snes, "Trying again in smaller region\n"));
554: /* check to see if progress is hopeless */
555: PetscCall(SNESTR_Converged_Private(snes, snes->iter, xnorm, ynorm, fnorm, &snes->reason, snes->cnvP));
556: if (!snes->reason) PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm));
557: if (snes->reason == SNES_CONVERGED_SNORM_RELATIVE) snes->reason = SNES_DIVERGED_TR_DELTA;
558: snes->numFailures++;
559: /* We're not progressing, so return with the current iterate */
560: if (snes->reason) break;
561: }
562: if (rho_satisfied) {
563: /* Update function values */
564: already_done = PETSC_FALSE;
565: fnorm = gnorm;
566: fk = fkp1;
568: /* New residual and linearization point */
569: PetscCall(VecCopy(G, F));
570: PetscCall(VecCopy(W, X));
572: /* Monitor convergence */
573: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
574: snes->iter++;
575: snes->norm = fnorm;
576: snes->xnorm = xnorm;
577: snes->ynorm = ynorm;
578: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
579: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits));
581: /* Test for convergence, xnorm = || X || */
582: PetscCall(VecNorm(X, NORM_2, &xnorm));
583: PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm));
584: PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
585: if (snes->reason) break;
586: }
587: }
589: if (clear_converged_test) {
590: PetscCall(KSPGetAndClearConvergenceTest(ksp, &ctx->convtest, &ctx->convctx, &ctx->convdestroy));
591: PetscCall(PetscFree(ctx));
592: PetscCall(KSPSetConvergenceTest(ksp, convtest, convctx, convdestroy));
593: }
594: PetscFunctionReturn(PETSC_SUCCESS);
595: }
597: static PetscErrorCode SNESSetUp_NEWTONTR(SNES snes)
598: {
599: PetscFunctionBegin;
600: PetscCall(SNESSetWorkVecs(snes, 4));
601: PetscCall(SNESSetUpMatrices(snes));
602: PetscFunctionReturn(PETSC_SUCCESS);
603: }
605: static PetscErrorCode SNESDestroy_NEWTONTR(SNES snes)
606: {
607: PetscFunctionBegin;
608: PetscCall(PetscFree(snes->data));
609: PetscFunctionReturn(PETSC_SUCCESS);
610: }
612: static PetscErrorCode SNESSetFromOptions_NEWTONTR(SNES snes, PetscOptionItems *PetscOptionsObject)
613: {
614: SNES_NEWTONTR *ctx = (SNES_NEWTONTR *)snes->data;
616: PetscFunctionBegin;
617: PetscOptionsHeadBegin(PetscOptionsObject, "SNES trust region options for nonlinear equations");
618: PetscCall(PetscOptionsReal("-snes_tr_tol", "Trust region tolerance", "SNESSetTrustRegionTolerance", snes->deltatol, &snes->deltatol, NULL));
619: PetscCall(PetscOptionsReal("-snes_tr_eta1", "eta1", "None", ctx->eta1, &ctx->eta1, NULL));
620: PetscCall(PetscOptionsReal("-snes_tr_eta2", "eta2", "None", ctx->eta2, &ctx->eta2, NULL));
621: PetscCall(PetscOptionsReal("-snes_tr_eta3", "eta3", "None", ctx->eta3, &ctx->eta3, NULL));
622: PetscCall(PetscOptionsReal("-snes_tr_t1", "t1", "None", ctx->t1, &ctx->t1, NULL));
623: PetscCall(PetscOptionsReal("-snes_tr_t2", "t2", "None", ctx->t2, &ctx->t2, NULL));
624: PetscCall(PetscOptionsReal("-snes_tr_deltaM", "deltaM", "None", ctx->deltaM, &ctx->deltaM, NULL));
625: PetscCall(PetscOptionsReal("-snes_tr_delta0", "delta0", "None", ctx->delta0, &ctx->delta0, NULL));
626: PetscCall(PetscOptionsReal("-snes_tr_kmdc", "sufficient decrease parameter", "None", ctx->kmdc, &ctx->kmdc, NULL));
627: PetscCall(PetscOptionsEnum("-snes_tr_fallback_type", "Type of fallback if subproblem solution is outside of the trust region", "SNESNewtonTRSetFallbackType", SNESNewtonTRFallbackTypes, (PetscEnum)ctx->fallback, (PetscEnum *)&ctx->fallback, NULL));
628: PetscOptionsHeadEnd();
629: PetscFunctionReturn(PETSC_SUCCESS);
630: }
632: static PetscErrorCode SNESView_NEWTONTR(SNES snes, PetscViewer viewer)
633: {
634: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
635: PetscBool iascii;
637: PetscFunctionBegin;
638: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
639: if (iascii) {
640: PetscCall(PetscViewerASCIIPrintf(viewer, " Trust region tolerance %g\n", (double)snes->deltatol));
641: PetscCall(PetscViewerASCIIPrintf(viewer, " eta1=%g, eta2=%g, eta3=%g\n", (double)tr->eta1, (double)tr->eta2, (double)tr->eta3));
642: PetscCall(PetscViewerASCIIPrintf(viewer, " delta0=%g, t1=%g, t2=%g, deltaM=%g\n", (double)tr->delta0, (double)tr->t1, (double)tr->t2, (double)tr->deltaM));
643: PetscCall(PetscViewerASCIIPrintf(viewer, " kmdc=%g\n", (double)tr->kmdc));
644: PetscCall(PetscViewerASCIIPrintf(viewer, " fallback=%s\n", SNESNewtonTRFallbackTypes[tr->fallback]));
645: }
646: PetscFunctionReturn(PETSC_SUCCESS);
647: }
649: /*MC
650: SNESNEWTONTR - Newton based nonlinear solver that uses trust-region dogleg method with Cauchy direction {cite}`nocedal2006numerical`
652: Options Database Keys:
653: + -snes_tr_tol <tol> - trust region tolerance
654: . -snes_tr_eta1 <eta1> - trust region parameter 0.0 <= eta1 <= eta2, rho >= eta1 breaks out of the inner iteration (default: eta1=0.001)
655: . -snes_tr_eta2 <eta2> - trust region parameter 0.0 <= eta1 <= eta2, rho <= eta2 shrinks the trust region (default: eta2=0.25)
656: . -snes_tr_eta3 <eta3> - trust region parameter eta3 > eta2, rho >= eta3 expands the trust region (default: eta3=0.75)
657: . -snes_tr_t1 <t1> - trust region parameter, shrinking factor of trust region (default: 0.25)
658: . -snes_tr_t2 <t2> - trust region parameter, expanding factor of trust region (default: 2.0)
659: . -snes_tr_deltaM <deltaM> - trust region parameter, max size of trust region (default: MAX_REAL)
660: . -snes_tr_delta0 <delta0> - trust region parameter, initial size of trust region (default: 0.2)
661: - -snes_tr_fallback_type <newton,cauchy,dogleg> - Solution strategy to test reduction when step is outside of trust region. Can use scaled Newton direction, Cauchy point (Steepest Descent direction) or dogleg method.
663: Level: deprecated (since 3.19)
665: .seealso: [](ch_snes), `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESSetTrustRegionTolerance()`,
666: `SNESNewtonTRPreCheck()`, `SNESNewtonTRGetPreCheck()`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()`,
667: `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRSetFallbackType()`
668: M*/
669: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONTR(SNES snes)
670: {
671: SNES_NEWTONTR *neP;
673: PetscFunctionBegin;
674: snes->ops->setup = SNESSetUp_NEWTONTR;
675: snes->ops->solve = SNESSolve_NEWTONTR;
676: snes->ops->destroy = SNESDestroy_NEWTONTR;
677: snes->ops->setfromoptions = SNESSetFromOptions_NEWTONTR;
678: snes->ops->view = SNESView_NEWTONTR;
680: snes->stol = 0.0;
681: snes->usesksp = PETSC_TRUE;
682: snes->npcside = PC_RIGHT;
683: snes->usesnpc = PETSC_TRUE;
685: snes->alwayscomputesfinalresidual = PETSC_TRUE;
687: PetscCall(PetscNew(&neP));
688: snes->data = (void *)neP;
689: neP->delta = 0.0;
690: neP->delta0 = 0.2;
691: neP->eta1 = 0.001;
692: neP->eta2 = 0.25;
693: neP->eta3 = 0.75;
694: neP->t1 = 0.25;
695: neP->t2 = 2.0;
696: neP->deltaM = 1.e10;
697: neP->fallback = SNES_TR_FALLBACK_NEWTON;
698: neP->kmdc = 0.0; /* by default do not use sufficient decrease */
699: PetscFunctionReturn(PETSC_SUCCESS);
700: }