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};
11: const char *const SNESNewtonTRQNTypes[] = {"NONE", "SAME", "DIFFERENT", "SNESNewtonTRQNType", "SNES_TR_QN_", NULL};
13: static PetscErrorCode SNESNewtonTRSetTolerances_TR(SNES snes, PetscReal delta_min, PetscReal delta_max, PetscReal delta_0)
14: {
15: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
17: PetscFunctionBegin;
18: if (delta_min == PETSC_DETERMINE) delta_min = tr->default_deltam;
19: if (delta_max == PETSC_DETERMINE) delta_max = tr->default_deltaM;
20: if (delta_0 == PETSC_DETERMINE) delta_0 = tr->default_delta0;
21: if (delta_min != PETSC_CURRENT) tr->deltam = delta_min;
22: if (delta_max != PETSC_CURRENT) tr->deltaM = delta_max;
23: if (delta_0 != PETSC_CURRENT) tr->delta0 = delta_0;
24: PetscFunctionReturn(PETSC_SUCCESS);
25: }
27: static PetscErrorCode SNESNewtonTRGetTolerances_TR(SNES snes, PetscReal *delta_min, PetscReal *delta_max, PetscReal *delta_0)
28: {
29: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
31: PetscFunctionBegin;
32: if (delta_min) *delta_min = tr->deltam;
33: if (delta_max) *delta_max = tr->deltaM;
34: if (delta_0) *delta_0 = tr->delta0;
35: PetscFunctionReturn(PETSC_SUCCESS);
36: }
38: static PetscErrorCode SNESComputeJacobian_MATLMVM(SNES snes, Vec X, Mat J, Mat B, void *dummy)
39: {
40: PetscFunctionBegin;
41: // PetscCall(MatLMVMSymBroydenSetDelta(B, _some_delta));
42: PetscCall(MatLMVMUpdate(B, X, snes->vec_func));
43: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
44: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
45: if (J != B) {
46: // PetscCall(MatLMVMSymBroydenSetDelta(J, _some_delta));
47: PetscCall(MatLMVMUpdate(J, X, snes->vec_func));
48: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
49: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
50: }
51: PetscFunctionReturn(PETSC_SUCCESS);
52: }
54: static PetscErrorCode SNESTR_KSPConverged_Private(KSP ksp, PetscInt n, PetscReal rnorm, KSPConvergedReason *reason, void *cctx)
55: {
56: SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx *)cctx;
57: SNES snes = ctx->snes;
58: SNES_NEWTONTR *neP = (SNES_NEWTONTR *)snes->data;
59: Vec x;
60: PetscReal nrm;
62: PetscFunctionBegin;
63: /* Determine norm of solution */
64: PetscCall(KSPBuildSolution(ksp, NULL, &x));
65: PetscCall(VecNorm(x, neP->norm, &nrm));
66: if (nrm >= neP->delta) {
67: PetscCall(PetscInfo(snes, "Ending linear iteration early due to exiting trust region, delta=%g, length=%g\n", (double)neP->delta, (double)nrm));
68: *reason = KSP_CONVERGED_STEP_LENGTH;
69: PetscFunctionReturn(PETSC_SUCCESS);
70: }
71: PetscCall((*ctx->convtest)(ksp, n, rnorm, reason, ctx->convctx));
72: if (*reason) PetscCall(PetscInfo(snes, "Default or user provided convergence test KSP iterations=%" PetscInt_FMT ", rnorm=%g\n", n, (double)rnorm));
73: PetscFunctionReturn(PETSC_SUCCESS);
74: }
76: static PetscErrorCode SNESTR_KSPConverged_Destroy(void *cctx)
77: {
78: SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx *)cctx;
80: PetscFunctionBegin;
81: PetscCall((*ctx->convdestroy)(ctx->convctx));
82: PetscCall(PetscFree(ctx));
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: static PetscErrorCode SNESTR_Converged_Private(SNES snes, PetscInt it, PetscReal xnorm, PetscReal pnorm, PetscReal fnorm, SNESConvergedReason *reason, void *dummy)
87: {
88: SNES_NEWTONTR *neP = (SNES_NEWTONTR *)snes->data;
90: PetscFunctionBegin;
91: *reason = SNES_CONVERGED_ITERATING;
92: if (neP->delta < neP->deltam) {
93: PetscCall(PetscInfo(snes, "Diverged due to too small a trust region %g<%g\n", (double)neP->delta, (double)neP->deltam));
94: *reason = SNES_DIVERGED_TR_DELTA;
95: } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
96: PetscCall(PetscInfo(snes, "Exceeded maximum number of function evaluations: %" PetscInt_FMT "\n", snes->max_funcs));
97: *reason = SNES_DIVERGED_FUNCTION_COUNT;
98: }
99: PetscFunctionReturn(PETSC_SUCCESS);
100: }
102: /*@
103: SNESNewtonTRSetNormType - Specify the type of norm to use for the computation of the trust region.
105: Input Parameters:
106: + snes - the nonlinear solver object
107: - norm - the norm type
109: Level: intermediate
111: .seealso: `SNESNEWTONTR`, `NormType`
112: @*/
113: PetscErrorCode SNESNewtonTRSetNormType(SNES snes, NormType norm)
114: {
115: PetscBool flg;
117: PetscFunctionBegin;
120: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
121: if (flg) {
122: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
124: tr->norm = norm;
125: }
126: PetscFunctionReturn(PETSC_SUCCESS);
127: }
129: /*@
130: SNESNewtonTRSetQNType - Specify to use a quasi-Newton model.
132: Input Parameters:
133: + snes - the nonlinear solver object
134: - use - the type of approximations to be used
136: Level: intermediate
138: Notes:
139: Options for the approximations can be set with the snes_tr_qn_ and snes_tr_qn_pre_ prefixes.
141: .seealso: `SNESNEWTONTR`, `SNESNewtonTRQNType`, `MATLMVM`
142: @*/
143: PetscErrorCode SNESNewtonTRSetQNType(SNES snes, SNESNewtonTRQNType use)
144: {
145: PetscBool flg;
147: PetscFunctionBegin;
150: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
151: if (flg) {
152: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
154: tr->qn = use;
155: }
156: PetscFunctionReturn(PETSC_SUCCESS);
157: }
159: /*@
160: SNESNewtonTRSetFallbackType - Set the type of fallback to use if the solution of the trust region subproblem is outside the radius
162: Input Parameters:
163: + snes - the nonlinear solver object
164: - ftype - the fallback type, see `SNESNewtonTRFallbackType`
166: Level: intermediate
168: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRPreCheck()`, `SNESNewtonTRGetPreCheck()`, `SNESNewtonTRSetPreCheck()`,
169: `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()`
170: @*/
171: PetscErrorCode SNESNewtonTRSetFallbackType(SNES snes, SNESNewtonTRFallbackType ftype)
172: {
173: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
174: PetscBool flg;
176: PetscFunctionBegin;
179: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
180: if (flg) tr->fallback = ftype;
181: PetscFunctionReturn(PETSC_SUCCESS);
182: }
184: /*@C
185: SNESNewtonTRSetPreCheck - Sets a user function that is called before the search step has been determined.
186: Allows the user a chance to change or override the trust region decision.
188: Logically Collective
190: Input Parameters:
191: + snes - the nonlinear solver object
192: . func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRPreCheck()`
193: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
195: Level: intermediate
197: Note:
198: This function is called BEFORE the function evaluation within the solver.
200: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRPreCheck()`, `SNESNewtonTRGetPreCheck()`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()`,
201: @*/
202: PetscErrorCode SNESNewtonTRSetPreCheck(SNES snes, PetscErrorCode (*func)(SNES, Vec, Vec, PetscBool *, void *), void *ctx)
203: {
204: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
205: PetscBool flg;
207: PetscFunctionBegin;
209: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
210: if (flg) {
211: if (func) tr->precheck = func;
212: if (ctx) tr->precheckctx = ctx;
213: }
214: PetscFunctionReturn(PETSC_SUCCESS);
215: }
217: /*@C
218: SNESNewtonTRGetPreCheck - Gets the pre-check function
220: Not Collective
222: Input Parameter:
223: . snes - the nonlinear solver context
225: Output Parameters:
226: + func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRPreCheck()`
227: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
229: Level: intermediate
231: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRPreCheck()`
232: @*/
233: PetscErrorCode SNESNewtonTRGetPreCheck(SNES snes, PetscErrorCode (**func)(SNES, Vec, Vec, PetscBool *, void *), void **ctx)
234: {
235: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
236: PetscBool flg;
238: PetscFunctionBegin;
240: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
241: PetscAssert(flg, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)snes)->type_name);
242: if (func) *func = tr->precheck;
243: if (ctx) *ctx = tr->precheckctx;
244: PetscFunctionReturn(PETSC_SUCCESS);
245: }
247: /*@C
248: SNESNewtonTRSetPostCheck - Sets a user function that is called after the search step has been determined but before the next
249: function evaluation. Allows the user a chance to change or override the internal decision of the solver
251: Logically Collective
253: Input Parameters:
254: + snes - the nonlinear solver object
255: . func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRPostCheck()`
256: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
258: Level: intermediate
260: Note:
261: This function is called BEFORE the function evaluation within the solver while the function set in
262: `SNESLineSearchSetPostCheck()` is called AFTER the function evaluation.
264: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRPostCheck()`, `SNESNewtonTRGetPostCheck()`, `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRGetPreCheck()`
265: @*/
266: PetscErrorCode SNESNewtonTRSetPostCheck(SNES snes, PetscErrorCode (*func)(SNES, Vec, Vec, Vec, PetscBool *, PetscBool *, void *), void *ctx)
267: {
268: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
269: PetscBool flg;
271: PetscFunctionBegin;
273: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
274: if (flg) {
275: if (func) tr->postcheck = func;
276: if (ctx) tr->postcheckctx = ctx;
277: }
278: PetscFunctionReturn(PETSC_SUCCESS);
279: }
281: /*@C
282: SNESNewtonTRGetPostCheck - Gets the post-check function
284: Not Collective
286: Input Parameter:
287: . snes - the nonlinear solver context
289: Output Parameters:
290: + func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRPostCheck()`
291: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
293: Level: intermediate
295: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRPostCheck()`
296: @*/
297: PetscErrorCode SNESNewtonTRGetPostCheck(SNES snes, PetscErrorCode (**func)(SNES, Vec, Vec, Vec, PetscBool *, PetscBool *, void *), void **ctx)
298: {
299: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
300: PetscBool flg;
302: PetscFunctionBegin;
304: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
305: PetscAssert(flg, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)snes)->type_name);
306: if (func) *func = tr->postcheck;
307: if (ctx) *ctx = tr->postcheckctx;
308: PetscFunctionReturn(PETSC_SUCCESS);
309: }
311: /*@C
312: SNESNewtonTRPreCheck - Runs the precheck routine
314: Logically Collective
316: Input Parameters:
317: + snes - the solver
318: . X - The last solution
319: - Y - The step direction
321: Output Parameter:
322: . changed_Y - Indicator that the step direction `Y` has been changed.
324: Level: intermediate
326: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRGetPreCheck()`, `SNESNewtonTRPostCheck()`
327: @*/
328: PetscErrorCode SNESNewtonTRPreCheck(SNES snes, Vec X, Vec Y, PetscBool *changed_Y)
329: {
330: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
331: PetscBool flg;
333: PetscFunctionBegin;
335: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
336: PetscAssert(flg, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)snes)->type_name);
337: *changed_Y = PETSC_FALSE;
338: if (tr->precheck) {
339: PetscCall((*tr->precheck)(snes, X, Y, changed_Y, tr->precheckctx));
341: }
342: PetscFunctionReturn(PETSC_SUCCESS);
343: }
345: /*@C
346: SNESNewtonTRPostCheck - Runs the postcheck routine
348: Logically Collective
350: Input Parameters:
351: + snes - the solver
352: . X - The last solution
353: . Y - The full step direction
354: - W - The updated solution, W = X - Y
356: Output Parameters:
357: + changed_Y - indicator if step has been changed
358: - changed_W - Indicator if the new candidate solution W has been changed.
360: Note:
361: If Y is changed then W is recomputed as X - Y
363: Level: intermediate
365: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()`, `SNESNewtonTRPreCheck()`
366: @*/
367: PetscErrorCode SNESNewtonTRPostCheck(SNES snes, Vec X, Vec Y, Vec W, PetscBool *changed_Y, PetscBool *changed_W)
368: {
369: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
370: PetscBool flg;
372: PetscFunctionBegin;
374: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
375: PetscAssert(flg, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)snes)->type_name);
376: *changed_Y = PETSC_FALSE;
377: *changed_W = PETSC_FALSE;
378: if (tr->postcheck) {
379: PetscCall((*tr->postcheck)(snes, X, Y, W, changed_Y, changed_W, tr->postcheckctx));
382: }
383: PetscFunctionReturn(PETSC_SUCCESS);
384: }
386: /* stable implementation of roots of a*x^2 + b*x + c = 0 */
387: static inline void PetscQuadraticRoots(PetscReal a, PetscReal b, PetscReal c, PetscReal *xm, PetscReal *xp)
388: {
389: PetscReal temp = -0.5 * (b + PetscCopysignReal(1.0, b) * PetscSqrtReal(b * b - 4 * a * c));
390: PetscReal x1 = temp / a;
391: PetscReal x2 = c / temp;
392: *xm = PetscMin(x1, x2);
393: *xp = PetscMax(x1, x2);
394: }
396: /* Computes the quadratic model difference */
397: static PetscErrorCode SNESNewtonTRQuadraticDelta(SNES snes, Mat J, PetscBool has_objective, Vec Y, Vec GradF, Vec W, PetscReal *yTHy_, PetscReal *gTy_, PetscReal *deltaqm)
398: {
399: PetscReal yTHy, gTy;
401: PetscFunctionBegin;
402: PetscCall(MatMult(J, Y, W));
403: if (has_objective) PetscCall(VecDotRealPart(Y, W, &yTHy));
404: else PetscCall(VecDotRealPart(W, W, &yTHy)); /* Gauss-Newton approximation J^t * J */
405: PetscCall(VecDotRealPart(GradF, Y, &gTy));
406: *deltaqm = -(-(gTy) + 0.5 * (yTHy)); /* difference in quadratic model, -gTy because SNES solves it this way */
407: if (yTHy_) *yTHy_ = yTHy;
408: if (gTy_) *gTy_ = gTy;
409: PetscFunctionReturn(PETSC_SUCCESS);
410: }
412: /* Computes the new objective given X = Xk, Y = direction
413: W work vector, on output W = X - Y
414: G work vector, on output G = SNESFunction(W) */
415: static PetscErrorCode SNESNewtonTRObjective(SNES snes, PetscBool has_objective, Vec X, Vec Y, Vec W, Vec G, PetscReal *gnorm, PetscReal *fkp1)
416: {
417: PetscBool changed_y, changed_w;
419: PetscFunctionBegin;
420: /* TODO: we can add a linesearch here */
421: PetscCall(SNESNewtonTRPreCheck(snes, X, Y, &changed_y));
422: PetscCall(VecWAXPY(W, -1.0, Y, X)); /* Xkp1 */
423: PetscCall(SNESNewtonTRPostCheck(snes, X, Y, W, &changed_y, &changed_w));
424: if (changed_y && !changed_w) PetscCall(VecWAXPY(W, -1.0, Y, X));
426: PetscCall(SNESComputeFunction(snes, W, G)); /* F(Xkp1) = G */
427: PetscCall(VecNorm(G, NORM_2, gnorm));
428: if (has_objective) PetscCall(SNESComputeObjective(snes, W, fkp1));
429: else *fkp1 = 0.5 * PetscSqr(*gnorm);
430: PetscFunctionReturn(PETSC_SUCCESS);
431: }
433: static PetscErrorCode SNESSetUpQN_NEWTONTR(SNES snes)
434: {
435: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
437: PetscFunctionBegin;
438: PetscCall(MatDestroy(&tr->qnB));
439: PetscCall(MatDestroy(&tr->qnB_pre));
440: if (tr->qn) {
441: PetscInt n, N;
442: const char *optionsprefix;
443: Mat B;
445: PetscCall(MatCreate(PetscObjectComm((PetscObject)snes), &B));
446: PetscCall(SNESGetOptionsPrefix(snes, &optionsprefix));
447: PetscCall(MatSetOptionsPrefix(B, "snes_tr_qn_"));
448: PetscCall(MatAppendOptionsPrefix(B, optionsprefix));
449: PetscCall(MatSetType(B, MATLMVMBFGS));
450: PetscCall(VecGetLocalSize(snes->vec_sol, &n));
451: PetscCall(VecGetSize(snes->vec_sol, &N));
452: PetscCall(MatSetSizes(B, n, n, N, N));
453: PetscCall(MatSetUp(B));
454: PetscCall(MatSetFromOptions(B));
455: PetscCall(MatLMVMAllocate(B, snes->vec_sol, snes->vec_func));
456: tr->qnB = B;
457: if (tr->qn == SNES_TR_QN_DIFFERENT) {
458: PetscCall(MatCreate(PetscObjectComm((PetscObject)snes), &B));
459: PetscCall(SNESGetOptionsPrefix(snes, &optionsprefix));
460: PetscCall(MatSetOptionsPrefix(B, "snes_tr_qn_pre_"));
461: PetscCall(MatAppendOptionsPrefix(B, optionsprefix));
462: PetscCall(MatSetType(B, MATLMVMBFGS));
463: PetscCall(MatSetSizes(B, n, n, N, N));
464: PetscCall(MatSetUp(B));
465: PetscCall(MatSetFromOptions(B));
466: PetscCall(MatLMVMAllocate(B, snes->vec_sol, snes->vec_func));
467: tr->qnB_pre = B;
468: } else {
469: PetscCall(PetscObjectReference((PetscObject)tr->qnB));
470: tr->qnB_pre = tr->qnB;
471: }
472: }
473: PetscFunctionReturn(PETSC_SUCCESS);
474: }
476: /*
477: SNESSolve_NEWTONTR - Implements Newton's Method with trust-region subproblem and adds dogleg Cauchy
478: (Steepest Descent direction) step and direction if the trust region is not satisfied for solving system of
479: nonlinear equations
481: */
482: static PetscErrorCode SNESSolve_NEWTONTR(SNES snes)
483: {
484: SNES_NEWTONTR *neP = (SNES_NEWTONTR *)snes->data;
485: Vec X, F, Y, G, W, GradF, YU, Yc;
486: PetscInt maxits, lits;
487: PetscReal rho, fnorm, gnorm = 0.0, xnorm = 0.0, delta, ynorm;
488: PetscReal deltaM, fk, fkp1, deltaqm = 0.0, gTy = 0.0, yTHy = 0.0;
489: PetscReal auk, tauk, gfnorm, gfnorm_k, ycnorm, gTBg, objmin = 0.0, beta_k = 1.0;
490: PC pc;
491: Mat J, Jp;
492: PetscBool already_done = PETSC_FALSE, on_boundary, use_cauchy;
493: PetscBool clear_converged_test, rho_satisfied, has_objective;
494: SNES_TR_KSPConverged_Ctx *ctx;
495: void *convctx;
496: SNESObjectiveFn *objective;
497: PetscErrorCode (*convtest)(KSP, PetscInt, PetscReal, KSPConvergedReason *, void *), (*convdestroy)(void *);
499: PetscFunctionBegin;
500: PetscCall(SNESGetObjective(snes, &objective, NULL));
501: has_objective = objective ? PETSC_TRUE : PETSC_FALSE;
503: maxits = snes->max_its; /* maximum number of iterations */
504: X = snes->vec_sol; /* solution vector */
505: F = snes->vec_func; /* residual vector */
506: Y = snes->vec_sol_update; /* update vector */
507: G = snes->work[0]; /* updated residual */
508: W = snes->work[1]; /* temporary vector */
509: GradF = !has_objective ? snes->work[2] : snes->vec_func; /* grad f = J^T F */
510: YU = snes->work[3]; /* work vector for dogleg method */
511: Yc = snes->work[4]; /* Cauchy point */
513: 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);
515: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
516: snes->iter = 0;
517: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
519: /* setup QN matrices if needed */
520: PetscCall(SNESSetUpQN_NEWTONTR(snes));
522: /* Set the linear stopping criteria to use the More' trick if needed */
523: clear_converged_test = PETSC_FALSE;
524: PetscCall(SNESGetKSP(snes, &snes->ksp));
525: PetscCall(KSPGetConvergenceTest(snes->ksp, &convtest, &convctx, &convdestroy));
526: if (convtest != SNESTR_KSPConverged_Private) {
527: clear_converged_test = PETSC_TRUE;
528: PetscCall(PetscNew(&ctx));
529: ctx->snes = snes;
530: PetscCall(KSPGetAndClearConvergenceTest(snes->ksp, &ctx->convtest, &ctx->convctx, &ctx->convdestroy));
531: PetscCall(KSPSetConvergenceTest(snes->ksp, SNESTR_KSPConverged_Private, ctx, SNESTR_KSPConverged_Destroy));
532: PetscCall(PetscInfo(snes, "Using Krylov convergence test SNESTR_KSPConverged_Private\n"));
533: }
535: if (!snes->vec_func_init_set) {
536: PetscCall(SNESComputeFunction(snes, X, F)); /* F(X) */
537: } else snes->vec_func_init_set = PETSC_FALSE;
539: PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- || F || */
540: SNESCheckFunctionNorm(snes, fnorm);
541: PetscCall(VecNorm(X, NORM_2, &xnorm)); /* xnorm <- || X || */
543: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
544: snes->norm = fnorm;
545: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
546: delta = neP->delta0;
547: deltaM = neP->deltaM;
548: neP->delta = delta;
549: PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));
551: /* test convergence */
552: rho_satisfied = PETSC_FALSE;
553: PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm));
554: PetscCall(SNESMonitor(snes, 0, fnorm));
555: if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
557: if (has_objective) PetscCall(SNESComputeObjective(snes, X, &fk));
558: else fk = 0.5 * PetscSqr(fnorm); /* obj(x) = 0.5 * ||F(x)||^2 */
560: /* hook state vector to BFGS preconditioner */
561: PetscCall(KSPGetPC(snes->ksp, &pc));
562: PetscCall(PCLMVMSetUpdateVec(pc, X));
564: if (neP->kmdc) PetscCall(KSPSetComputeEigenvalues(snes->ksp, PETSC_TRUE));
566: while (snes->iter < maxits) {
567: /* calculating Jacobian and GradF of minimization function only once */
568: if (!already_done) {
569: /* Call general purpose update function */
570: PetscTryTypeMethod(snes, update, snes->iter);
572: /* apply the nonlinear preconditioner */
573: if (snes->npc && snes->npcside == PC_RIGHT) {
574: SNESConvergedReason reason;
576: PetscCall(SNESSetInitialFunction(snes->npc, F));
577: PetscCall(PetscLogEventBegin(SNES_NPCSolve, snes->npc, X, snes->vec_rhs, 0));
578: PetscCall(SNESSolve(snes->npc, snes->vec_rhs, X));
579: PetscCall(PetscLogEventEnd(SNES_NPCSolve, snes->npc, X, snes->vec_rhs, 0));
580: PetscCall(SNESGetConvergedReason(snes->npc, &reason));
581: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
582: snes->reason = SNES_DIVERGED_INNER;
583: PetscFunctionReturn(PETSC_SUCCESS);
584: }
585: // XXX
586: PetscCall(SNESGetNPCFunction(snes, F, &fnorm));
587: }
589: /* Jacobian */
590: J = NULL;
591: Jp = NULL;
592: if (!neP->qnB) {
593: PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
594: J = snes->jacobian;
595: Jp = snes->jacobian_pre;
596: } else { /* QN model */
597: PetscCall(SNESComputeJacobian_MATLMVM(snes, X, neP->qnB, neP->qnB_pre, NULL));
598: J = neP->qnB;
599: Jp = neP->qnB_pre;
600: }
601: SNESCheckJacobianDomainerror(snes);
603: /* objective function */
604: PetscCall(VecNorm(F, NORM_2, &fnorm));
605: if (has_objective) PetscCall(SNESComputeObjective(snes, X, &fk));
606: else fk = 0.5 * PetscSqr(fnorm); /* obj(x) = 0.5 * ||F(x)||^2 */
608: /* GradF */
609: if (has_objective) gfnorm = fnorm;
610: else {
611: PetscCall(MatMultTranspose(J, F, GradF)); /* grad f = J^T F */
612: PetscCall(VecNorm(GradF, NORM_2, &gfnorm));
613: }
614: PetscCall(VecNorm(GradF, neP->norm, &gfnorm_k));
615: }
616: already_done = PETSC_TRUE;
618: /* solve trust-region subproblem */
620: /* first compute Cauchy Point */
621: PetscCall(MatMult(J, GradF, W));
622: if (has_objective) PetscCall(VecDotRealPart(GradF, W, &gTBg));
623: else PetscCall(VecDotRealPart(W, W, &gTBg)); /* B = J^t * J */
624: /* Eqs 4.11 and 4.12 in Nocedal and Wright 2nd Edition (4.7 and 4.8 in 1st Edition) */
625: auk = delta / gfnorm_k;
626: if (gTBg < 0.0) tauk = 1.0;
627: else tauk = PetscMin(gfnorm * gfnorm * gfnorm_k / (delta * gTBg), 1);
628: auk *= tauk;
629: ycnorm = auk * gfnorm;
630: PetscCall(VecAXPBY(Yc, auk, 0.0, GradF));
632: on_boundary = PETSC_FALSE;
633: use_cauchy = (PetscBool)(tauk == 1.0 && has_objective);
634: if (!use_cauchy) {
635: KSPConvergedReason reason;
637: /* sufficient decrease (see 6.3.27 in Conn, Gould, Toint "Trust Region Methods")
638: beta_k the largest eigenvalue of the Hessian. Here we use the previous estimated value */
639: objmin = -neP->kmdc * gnorm * PetscMin(gnorm / beta_k, delta);
640: PetscCall(KSPCGSetObjectiveTarget(snes->ksp, objmin));
642: /* specify radius if looking for Newton step and trust region norm is the l2 norm */
643: PetscCall(KSPCGSetRadius(snes->ksp, neP->fallback == SNES_TR_FALLBACK_NEWTON && neP->norm == NORM_2 ? delta : 0.0));
644: PetscCall(KSPSetOperators(snes->ksp, J, Jp));
645: PetscCall(KSPSolve(snes->ksp, F, Y));
646: SNESCheckKSPSolve(snes);
647: PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
648: PetscCall(KSPGetConvergedReason(snes->ksp, &reason));
649: on_boundary = (PetscBool)(reason == KSP_CONVERGED_STEP_LENGTH);
650: PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits));
651: if (neP->kmdc) { /* update estimated Hessian largest eigenvalue */
652: PetscReal emax, emin;
653: PetscCall(KSPComputeExtremeSingularValues(snes->ksp, &emax, &emin));
654: if (emax > 0.0) beta_k = emax + 1;
655: }
656: } else { /* Cauchy point is on the boundary, accept it */
657: on_boundary = PETSC_TRUE;
658: PetscCall(VecCopy(Yc, Y));
659: PetscCall(PetscInfo(snes, "CP evaluated on boundary. delta: %g, ycnorm: %g, gTBg: %g\n", (double)delta, (double)ycnorm, (double)gTBg));
660: }
661: PetscCall(VecNorm(Y, neP->norm, &ynorm));
663: /* decide what to do when the update is outside of trust region */
664: if (!use_cauchy && (ynorm > delta || ynorm == 0.0)) {
665: SNESNewtonTRFallbackType fallback = ynorm > 0.0 ? neP->fallback : SNES_TR_FALLBACK_CAUCHY;
667: PetscCheck(neP->norm == NORM_2 || fallback != SNES_TR_FALLBACK_DOGLEG, PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "DOGLEG without l2 norm not implemented");
668: switch (fallback) {
669: case SNES_TR_FALLBACK_NEWTON:
670: auk = delta / ynorm;
671: PetscCall(VecScale(Y, auk));
672: PetscCall(PetscInfo(snes, "SN evaluated. delta: %g, ynorm: %g\n", (double)delta, (double)ynorm));
673: break;
674: case SNES_TR_FALLBACK_CAUCHY:
675: case SNES_TR_FALLBACK_DOGLEG:
676: if (fallback == SNES_TR_FALLBACK_CAUCHY || gTBg <= 0.0) {
677: PetscCall(VecCopy(Yc, Y));
678: PetscCall(PetscInfo(snes, "CP evaluated. delta: %g, ynorm: %g, ycnorm: %g, gTBg: %g\n", (double)delta, (double)ynorm, (double)ycnorm, (double)gTBg));
679: } else { /* take linear combination of Cauchy and Newton direction and step */
680: auk = gfnorm * gfnorm / gTBg;
681: if (gfnorm_k * auk >= delta) { /* first leg: Cauchy point outside of trust region */
682: PetscCall(VecAXPBY(Y, delta / gfnorm_k, 0.0, GradF));
683: PetscCall(PetscInfo(snes, "CP evaluated (outside region). delta: %g, ynorm: %g, ycnorm: %g\n", (double)delta, (double)ynorm, (double)ycnorm));
684: } else { /* second leg */
685: PetscReal c0, c1, c2, tau = 0.0, tpos, tneg;
686: PetscBool noroots;
688: /* Find solutions of (Eq. 4.16 in Nocedal and Wright)
689: ||p_U + lambda * (p_B - p_U)||^2 - delta^2 = 0,
690: where p_U the Cauchy direction, p_B the Newton direction */
691: PetscCall(VecAXPBY(YU, auk, 0.0, GradF));
692: PetscCall(VecAXPY(Y, -1.0, YU));
693: PetscCall(VecNorm(Y, NORM_2, &c0));
694: PetscCall(VecDotRealPart(YU, Y, &c1));
695: c0 = PetscSqr(c0);
696: c2 = PetscSqr(ycnorm) - PetscSqr(delta);
697: PetscQuadraticRoots(c0, 2 * c1, c2, &tneg, &tpos);
699: /* In principle the DL strategy as a unique solution in [0,1]
700: here we check that for some reason we numerically failed
701: to compute it. In that case, we use the Cauchy point */
702: noroots = PetscIsInfOrNanReal(tneg);
703: if (!noroots) {
704: if (tpos > 1) {
705: if (tneg >= 0 && tneg <= 1) {
706: tau = tneg;
707: } else noroots = PETSC_TRUE;
708: } else if (tpos >= 0) {
709: tau = tpos;
710: } else noroots = PETSC_TRUE;
711: }
712: if (noroots) { /* No roots, select Cauchy point */
713: PetscCall(VecCopy(Yc, Y));
714: } else {
715: PetscCall(VecAXPBY(Y, 1.0, tau, YU));
716: }
717: PetscCall(PetscInfo(snes, "%s evaluated. roots: (%g, %g), tau %g, ynorm: %g, ycnorm: %g, gTBg: %g\n", noroots ? "CP" : "DL", (double)tneg, (double)tpos, (double)tau, (double)ynorm, (double)ycnorm, (double)gTBg));
718: }
719: }
720: break;
721: default:
722: SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "Unknown fallback mode");
723: break;
724: }
725: }
727: /* compute the quadratic model difference */
728: PetscCall(SNESNewtonTRQuadraticDelta(snes, J, has_objective, Y, GradF, W, &yTHy, &gTy, &deltaqm));
730: /* Compute new objective function */
731: PetscCall(SNESNewtonTRObjective(snes, has_objective, X, Y, W, G, &gnorm, &fkp1));
732: if (PetscIsInfOrNanReal(fkp1)) rho = neP->eta1;
733: else {
734: if (deltaqm > 0.0) rho = (fk - fkp1) / deltaqm; /* actual improvement over predicted improvement */
735: else rho = neP->eta1; /* no reduction in quadratic model, step must be rejected */
736: }
738: PetscCall(VecNorm(Y, neP->norm, &ynorm));
739: PetscCall(PetscInfo(snes, "rho=%g, delta=%g, fk=%g, fkp1=%g, deltaqm=%g, gTy=%g, yTHy=%g, ynormk=%g\n", (double)rho, (double)delta, (double)fk, (double)fkp1, (double)deltaqm, (double)gTy, (double)yTHy, (double)ynorm));
741: /* update the size of the trust region */
742: if (rho < neP->eta2) delta *= neP->t1; /* shrink the region */
743: else if (rho > neP->eta3 && on_boundary) delta *= neP->t2; /* expand the region */
744: delta = PetscMin(delta, deltaM); /* but not greater than deltaM */
746: /* log 2-norm of update for moniroting routines */
747: PetscCall(VecNorm(Y, NORM_2, &ynorm));
749: /* decide on new step */
750: neP->delta = delta;
751: if (rho > neP->eta1) {
752: rho_satisfied = PETSC_TRUE;
753: } else {
754: rho_satisfied = PETSC_FALSE;
755: PetscCall(PetscInfo(snes, "Trying again in smaller region\n"));
756: /* check to see if progress is hopeless */
757: PetscCall(SNESTR_Converged_Private(snes, snes->iter, xnorm, ynorm, fnorm, &snes->reason, snes->cnvP));
758: if (!snes->reason) PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm));
759: if (snes->reason == SNES_CONVERGED_SNORM_RELATIVE) snes->reason = SNES_DIVERGED_TR_DELTA;
760: snes->numFailures++;
761: /* We're not progressing, so return with the current iterate */
762: if (snes->reason) break;
763: }
764: if (rho_satisfied) {
765: /* Update function values */
766: already_done = PETSC_FALSE;
767: fnorm = gnorm;
768: fk = fkp1;
770: /* New residual and linearization point */
771: PetscCall(VecCopy(G, F));
772: PetscCall(VecCopy(W, X));
774: /* Monitor convergence */
775: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
776: snes->iter++;
777: snes->norm = fnorm;
778: snes->xnorm = xnorm;
779: snes->ynorm = ynorm;
780: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
781: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits));
783: /* Test for convergence, xnorm = || X || */
784: PetscCall(VecNorm(X, NORM_2, &xnorm));
785: PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm));
786: PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
787: if (snes->reason) break;
788: }
789: }
791: if (clear_converged_test) {
792: PetscCall(KSPGetAndClearConvergenceTest(snes->ksp, &ctx->convtest, &ctx->convctx, &ctx->convdestroy));
793: PetscCall(PetscFree(ctx));
794: PetscCall(KSPSetConvergenceTest(snes->ksp, convtest, convctx, convdestroy));
795: }
796: PetscFunctionReturn(PETSC_SUCCESS);
797: }
799: static PetscErrorCode SNESSetUp_NEWTONTR(SNES snes)
800: {
801: PetscFunctionBegin;
802: PetscCall(SNESSetWorkVecs(snes, 5));
803: PetscCall(SNESSetUpMatrices(snes));
804: PetscFunctionReturn(PETSC_SUCCESS);
805: }
807: static PetscErrorCode SNESReset_NEWTONTR(SNES snes)
808: {
809: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
811: PetscFunctionBegin;
812: PetscCall(MatDestroy(&tr->qnB));
813: PetscCall(MatDestroy(&tr->qnB_pre));
814: PetscFunctionReturn(PETSC_SUCCESS);
815: }
817: static PetscErrorCode SNESDestroy_NEWTONTR(SNES snes)
818: {
819: PetscFunctionBegin;
820: PetscCall(SNESReset_NEWTONTR(snes));
821: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonTRSetTolerances_C", NULL));
822: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonTRGetTolerances_C", NULL));
823: PetscCall(PetscFree(snes->data));
824: PetscFunctionReturn(PETSC_SUCCESS);
825: }
827: static PetscErrorCode SNESSetFromOptions_NEWTONTR(SNES snes, PetscOptionItems *PetscOptionsObject)
828: {
829: SNES_NEWTONTR *ctx = (SNES_NEWTONTR *)snes->data;
830: SNESNewtonTRQNType qn;
831: SNESNewtonTRFallbackType fallback;
832: NormType norm;
833: PetscBool flg;
835: PetscFunctionBegin;
836: PetscOptionsHeadBegin(PetscOptionsObject, "SNES trust region options for nonlinear equations");
837: PetscCall(PetscOptionsDeprecated("-snes_tr_deltaM", "-snes_tr_deltamax", "3.22", NULL));
838: PetscCall(PetscOptionsReal("-snes_tr_eta1", "eta1", "SNESNewtonTRSetUpdateParameters", ctx->eta1, &ctx->eta1, NULL));
839: PetscCall(PetscOptionsReal("-snes_tr_eta2", "eta2", "SNESNewtonTRSetUpdateParameters", ctx->eta2, &ctx->eta2, NULL));
840: PetscCall(PetscOptionsReal("-snes_tr_eta3", "eta3", "SNESNewtonTRSetUpdateParameters", ctx->eta3, &ctx->eta3, NULL));
841: PetscCall(PetscOptionsReal("-snes_tr_t1", "t1", "SNESNewtonTRSetUpdateParameters", ctx->t1, &ctx->t1, NULL));
842: PetscCall(PetscOptionsReal("-snes_tr_t2", "t2", "SNESNewtonTRSetUpdateParameters", ctx->t2, &ctx->t2, NULL));
843: PetscCall(PetscOptionsReal("-snes_tr_delta0", "Initial trust region size", "SNESNewtonTRSetTolerances", ctx->delta0, &ctx->delta0, NULL));
844: PetscCall(PetscOptionsReal("-snes_tr_deltamin", "Minimum allowed trust region size", "SNESNewtonTRSetTolerances", ctx->deltam, &ctx->deltam, NULL));
845: PetscCall(PetscOptionsReal("-snes_tr_deltamax", "Maximum allowed trust region size", "SNESNewtonTRSetTolerances", ctx->deltaM, &ctx->deltaM, NULL));
846: PetscCall(PetscOptionsReal("-snes_tr_kmdc", "sufficient decrease parameter", "None", ctx->kmdc, &ctx->kmdc, NULL));
848: fallback = ctx->fallback;
849: PetscCall(PetscOptionsEnum("-snes_tr_fallback_type", "Type of fallback if subproblem solution is outside of the trust region", "SNESNewtonTRSetFallbackType", SNESNewtonTRFallbackTypes, (PetscEnum)fallback, (PetscEnum *)&fallback, &flg));
850: if (flg) PetscCall(SNESNewtonTRSetFallbackType(snes, fallback));
852: qn = ctx->qn;
853: PetscCall(PetscOptionsEnum("-snes_tr_qn", "Use Quasi-Newton approximations for the model", "SNESNewtonTRSetQNType", SNESNewtonTRQNTypes, (PetscEnum)qn, (PetscEnum *)&qn, &flg));
854: if (flg) PetscCall(SNESNewtonTRSetQNType(snes, qn));
856: norm = ctx->norm;
857: PetscCall(PetscOptionsEnum("-snes_tr_norm_type", "Type of norm for trust region bounds", "SNESNewtonTRSetNormType", NormTypes, (PetscEnum)norm, (PetscEnum *)&norm, &flg));
858: if (flg) PetscCall(SNESNewtonTRSetNormType(snes, norm));
860: PetscOptionsHeadEnd();
861: PetscFunctionReturn(PETSC_SUCCESS);
862: }
864: static PetscErrorCode SNESView_NEWTONTR(SNES snes, PetscViewer viewer)
865: {
866: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
867: PetscBool iascii;
869: PetscFunctionBegin;
870: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
871: if (iascii) {
872: PetscCall(PetscViewerASCIIPrintf(viewer, " Trust region parameters:\n"));
873: PetscCall(PetscViewerASCIIPrintf(viewer, " eta1=%g, eta2=%g, eta3=%g\n", (double)tr->eta1, (double)tr->eta2, (double)tr->eta3));
874: PetscCall(PetscViewerASCIIPrintf(viewer, " t1=%g, t2=%g\n", (double)tr->t1, (double)tr->t2));
875: PetscCall(PetscViewerASCIIPrintf(viewer, " delta_min=%g, delta_0=%g, delta_max=%g\n", (double)tr->deltam, (double)tr->delta0, (double)tr->deltaM));
876: PetscCall(PetscViewerASCIIPrintf(viewer, " kmdc=%g\n", (double)tr->kmdc));
877: PetscCall(PetscViewerASCIIPrintf(viewer, " fallback=%s\n", SNESNewtonTRFallbackTypes[tr->fallback]));
878: if (tr->qn) PetscCall(PetscViewerASCIIPrintf(viewer, " qn=%s\n", SNESNewtonTRQNTypes[tr->qn]));
879: if (tr->norm != NORM_2) PetscCall(PetscViewerASCIIPrintf(viewer, " norm=%s\n", NormTypes[tr->norm]));
880: }
881: PetscFunctionReturn(PETSC_SUCCESS);
882: }
884: /*@
885: SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance.
887: Logically Collective
889: Input Parameters:
890: + snes - the `SNES` context
891: - tol - tolerance
893: Level: deprecated
895: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTR`, `SNESSetTolerances()`
896: @*/
897: PetscErrorCode SNESSetTrustRegionTolerance(SNES snes, PetscReal tol)
898: {
899: return SNESNewtonTRSetTolerances(snes, tol, PETSC_CURRENT, PETSC_CURRENT);
900: }
902: /*@
903: SNESNewtonTRSetTolerances - Sets the trust region parameter tolerances.
905: Logically Collective
907: Input Parameters:
908: + snes - the `SNES` context
909: . delta_min - minimum allowed trust region size
910: . delta_max - maximum allowed trust region size
911: - delta_0 - initial trust region size
913: Options Database Key:
914: + -snes_tr_deltamin <tol> - Set minimum size
915: . -snes_tr_deltamax <tol> - Set maximum size
916: - -snes_tr_delta0 <tol> - Set initial size
918: Note:
919: Use `PETSC_DETERMINE` to use the default value for the given `SNES`.
920: Use `PETSC_CURRENT` to retain a value.
922: Fortran Note:
923: Use `PETSC_DETERMINE_REAL`, `PETSC_CURRENT_REAL`
925: Level: intermediate
927: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTR`, `SNESNewtonTRGetTolerances()`
928: @*/
929: PetscErrorCode SNESNewtonTRSetTolerances(SNES snes, PetscReal delta_min, PetscReal delta_max, PetscReal delta_0)
930: {
931: PetscFunctionBegin;
936: PetscTryMethod(snes, "SNESNewtonTRSetTolerances_C", (SNES, PetscReal, PetscReal, PetscReal), (snes, delta_min, delta_max, delta_0));
937: PetscFunctionReturn(PETSC_SUCCESS);
938: }
940: /*@
941: SNESNewtonTRGetTolerances - Gets the trust region parameter tolerances.
943: Not Collective
945: Input Parameter:
946: . snes - the `SNES` context
948: Output Parameters:
949: + delta_min - minimum allowed trust region size or `NULL`
950: . delta_max - maximum allowed trust region size or `NULL`
951: - delta_0 - initial trust region size or `NULL`
953: Level: intermediate
955: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTR`, `SNESNewtonTRSetTolerances()`
956: @*/
957: PetscErrorCode SNESNewtonTRGetTolerances(SNES snes, PetscReal *delta_min, PetscReal *delta_max, PetscReal *delta_0)
958: {
959: PetscFunctionBegin;
961: if (delta_min) PetscAssertPointer(delta_min, 2);
962: if (delta_max) PetscAssertPointer(delta_max, 3);
963: if (delta_0) PetscAssertPointer(delta_0, 4);
964: PetscUseMethod(snes, "SNESNewtonTRGetTolerances_C", (SNES, PetscReal *, PetscReal *, PetscReal *), (snes, delta_min, delta_max, delta_0));
965: PetscFunctionReturn(PETSC_SUCCESS);
966: }
968: /*@
969: SNESNewtonTRSetUpdateParameters - Sets the trust region update parameters.
971: Logically Collective
973: Input Parameters:
974: + snes - the `SNES` context
975: . eta1 - acceptance tolerance
976: . eta2 - shrinking tolerance
977: . eta3 - enlarging tolerance
978: . t1 - shrink factor
979: - t2 - enlarge factor
981: Options Database Key:
982: + -snes_tr_eta1 <tol> - Set eta1
983: . -snes_tr_eta2 <tol> - Set eta2
984: . -snes_tr_eta3 <tol> - Set eta3
985: . -snes_tr_t1 <tol> - Set t1
986: - -snes_tr_t2 <tol> - Set t2
988: Notes:
989: Given the ratio $\rho = \frac{f(x_k) - f(x_k+s_k)}{m(0) - m(s_k)}$, with $x_k$ the current iterate,
990: $s_k$ the computed step, $f$ the objective function, and $m$ the quadratic model, the trust region
991: radius is modified as follows
992: $$
993: \delta =
994: \begin{cases}
995: \delta * t_1 ,& \rho < \eta_2 \\
996: \delta * t_2 ,& \rho > \eta_3 \\
997: \end{cases}
998: $$
999: The step is accepted if $\rho > \eta_1$.
1000: Use `PETSC_DETERMINE` to use the default value for the given `SNES`.
1001: Use `PETSC_CURRENT` to retain a value.
1003: Fortran Note:
1004: Use `PETSC_DETERMINE_REAL`, `PETSC_CURRENT_REAL`
1006: Level: intermediate
1008: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTR`, `SNESSetObjective()`, `SNESNewtonTRGetUpdateParameters()`
1009: @*/
1010: PetscErrorCode SNESNewtonTRSetUpdateParameters(SNES snes, PetscReal eta1, PetscReal eta2, PetscReal eta3, PetscReal t1, PetscReal t2)
1011: {
1012: PetscBool flg;
1014: PetscFunctionBegin;
1021: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
1022: if (flg) {
1023: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
1025: if (eta1 == PETSC_DETERMINE) eta1 = tr->default_eta1;
1026: if (eta2 == PETSC_DETERMINE) eta2 = tr->default_eta2;
1027: if (eta3 == PETSC_DETERMINE) eta3 = tr->default_eta3;
1028: if (t1 == PETSC_DETERMINE) t1 = tr->default_t1;
1029: if (t2 == PETSC_DETERMINE) t2 = tr->default_t2;
1030: if (eta1 != PETSC_CURRENT) tr->eta1 = eta1;
1031: if (eta2 != PETSC_CURRENT) tr->eta2 = eta2;
1032: if (eta3 != PETSC_CURRENT) tr->eta3 = eta3;
1033: if (t1 != PETSC_CURRENT) tr->t1 = t1;
1034: if (t2 != PETSC_CURRENT) tr->t2 = t2;
1035: }
1036: PetscFunctionReturn(PETSC_SUCCESS);
1037: }
1039: /*@
1040: SNESNewtonTRGetUpdateParameters - Gets the trust region update parameters.
1042: Not Collective
1044: Input Parameter:
1045: . snes - the `SNES` context
1047: Output Parameters:
1048: + eta1 - acceptance tolerance
1049: . eta2 - shrinking tolerance
1050: . eta3 - enlarging tolerance
1051: . t1 - shrink factor
1052: - t2 - enlarge factor
1054: Level: intermediate
1056: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTR`, `SNESNewtonTRSetUpdateParameters()`
1057: @*/
1058: PetscErrorCode SNESNewtonTRGetUpdateParameters(SNES snes, PetscReal *eta1, PetscReal *eta2, PetscReal *eta3, PetscReal *t1, PetscReal *t2)
1059: {
1060: SNES_NEWTONTR *tr;
1061: PetscBool flg;
1063: PetscFunctionBegin;
1065: if (eta1) PetscAssertPointer(eta1, 2);
1066: if (eta2) PetscAssertPointer(eta2, 3);
1067: if (eta3) PetscAssertPointer(eta3, 4);
1068: if (t1) PetscAssertPointer(t1, 5);
1069: if (t2) PetscAssertPointer(t2, 6);
1070: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
1071: PetscAssert(flg, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)snes)->type_name);
1072: tr = (SNES_NEWTONTR *)snes->data;
1073: if (eta1) *eta1 = tr->eta1;
1074: if (eta2) *eta2 = tr->eta2;
1075: if (eta3) *eta3 = tr->eta3;
1076: if (t1) *t1 = tr->t1;
1077: if (t2) *t2 = tr->t2;
1078: PetscFunctionReturn(PETSC_SUCCESS);
1079: }
1081: /*MC
1082: SNESNEWTONTR - Newton based nonlinear solver that uses a trust-region strategy
1084: Options Database Keys:
1085: + -snes_tr_deltamin <deltamin> - trust region parameter, minimum size of trust region
1086: . -snes_tr_deltamax <deltamax> - trust region parameter, max size of trust region (default: 1e10)
1087: . -snes_tr_delta0 <delta0> - trust region parameter, initial size of trust region (default: 0.2)
1088: . -snes_tr_eta1 <eta1> - trust region parameter eta1 <= eta2, rho > eta1 breaks out of the inner iteration (default: 0.001)
1089: . -snes_tr_eta2 <eta2> - trust region parameter, rho <= eta2 shrinks the trust region (default: 0.25)
1090: . -snes_tr_eta3 <eta3> - trust region parameter eta3 > eta2, rho >= eta3 expands the trust region (default: 0.75)
1091: . -snes_tr_t1 <t1> - trust region parameter, shrinking factor of trust region (default: 0.25)
1092: . -snes_tr_t2 <t2> - trust region parameter, expanding factor of trust region (default: 2.0)
1093: . -snes_tr_norm_type <1,2,infinity> - Type of norm for trust region bounds (default: "2")
1094: - -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.
1096: Level: beginner
1098: Notes:
1099: The code is largely based on the book {cite}`nocedal2006numerical` and supports minimizing objective functions using a quadratic model.
1100: Quasi-Newton models are also supported.
1102: Default step computation uses the Newton direction, but a dogleg type update is also supported.
1103: The 1- and infinity-norms are also supported when computing the trust region bounds.
1105: .seealso: [](ch_snes), `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESSetObjective()`,
1106: `SNESNewtonTRSetTolerances()`, `SNESNewtonTRSetUpdateParameters()`
1107: `SNESNewtonTRSetNormType()`, `SNESNewtonTRSetFallbackType()`, `SNESNewtonTRSetQNType()`
1108: `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRSetPreCheck()`,
1109: M*/
1110: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONTR(SNES snes)
1111: {
1112: SNES_NEWTONTR *neP;
1114: PetscFunctionBegin;
1115: snes->ops->setup = SNESSetUp_NEWTONTR;
1116: snes->ops->solve = SNESSolve_NEWTONTR;
1117: snes->ops->reset = SNESReset_NEWTONTR;
1118: snes->ops->destroy = SNESDestroy_NEWTONTR;
1119: snes->ops->setfromoptions = SNESSetFromOptions_NEWTONTR;
1120: snes->ops->view = SNESView_NEWTONTR;
1122: PetscCall(SNESParametersInitialize(snes));
1123: PetscObjectParameterSetDefault(snes, stol, 0.0);
1125: snes->usesksp = PETSC_TRUE;
1126: snes->npcside = PC_RIGHT;
1127: snes->usesnpc = PETSC_TRUE;
1129: snes->alwayscomputesfinalresidual = PETSC_TRUE;
1131: PetscCall(PetscNew(&neP));
1132: snes->data = (void *)neP;
1134: PetscObjectParameterSetDefault(neP, eta1, 0.001);
1135: PetscObjectParameterSetDefault(neP, eta2, 0.25);
1136: PetscObjectParameterSetDefault(neP, eta3, 0.75);
1137: PetscObjectParameterSetDefault(neP, t1, 0.25);
1138: PetscObjectParameterSetDefault(neP, t2, 2.0);
1139: PetscObjectParameterSetDefault(neP, deltam, PetscDefined(USE_REAL_SINGLE) ? 1.e-6 : 1.e-12);
1140: PetscObjectParameterSetDefault(neP, delta0, 0.2);
1141: PetscObjectParameterSetDefault(neP, deltaM, 1.e10);
1143: neP->norm = NORM_2;
1144: neP->fallback = SNES_TR_FALLBACK_NEWTON;
1145: neP->kmdc = 0.0; /* by default do not use sufficient decrease */
1147: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonTRSetTolerances_C", SNESNewtonTRSetTolerances_TR));
1148: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonTRGetTolerances_C", SNESNewtonTRGetTolerances_TR));
1149: PetscFunctionReturn(PETSC_SUCCESS);
1150: }