Actual source code: tr.c
1: #include <../src/snes/impls/tr/trimpl.h>
3: typedef struct {
4: SNES snes;
5: KSPConvergenceTestFn *convtest;
6: PetscCtxDestroyFn *convdestroy;
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(PetscCtxRt 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, PetscReal xnorm, PetscReal ynorm, PetscReal fnorm)
87: {
88: SNES_NEWTONTR *neP = (SNES_NEWTONTR *)snes->data;
90: PetscFunctionBegin;
91: /* only check convergence for iter > 0 */
92: if (snes->iter) PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm));
93: if (snes->reason > 0) PetscFunctionReturn(PETSC_SUCCESS);
94: if (neP->delta < neP->deltam) {
95: PetscCall(PetscInfo(snes, "Diverged due to too small a trust region %g<%g\n", (double)neP->delta, (double)neP->deltam));
96: snes->reason = SNES_DIVERGED_TR_DELTA;
97: } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
98: PetscCall(PetscInfo(snes, "Exceeded maximum number of function evaluations: %" PetscInt_FMT "\n", snes->max_funcs));
99: snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
100: }
101: PetscFunctionReturn(PETSC_SUCCESS);
102: }
104: /*@
105: SNESNewtonTRSetNormType - Specify the type of norm to use for the computation of the trust region.
107: Input Parameters:
108: + snes - the nonlinear solver object
109: - norm - the norm type
111: Level: intermediate
113: .seealso: `SNESNEWTONTR`, `NormType`
114: @*/
115: PetscErrorCode SNESNewtonTRSetNormType(SNES snes, NormType norm)
116: {
117: PetscBool flg;
119: PetscFunctionBegin;
122: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
123: if (flg) {
124: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
126: tr->norm = norm;
127: }
128: PetscFunctionReturn(PETSC_SUCCESS);
129: }
131: /*@
132: SNESNewtonTRSetQNType - Specify to use a quasi-Newton model.
134: Input Parameters:
135: + snes - the nonlinear solver object
136: - use - the type of approximations to be used
138: Level: intermediate
140: Notes:
141: Options for the approximations can be set with the `snes_tr_qn_` and `snes_tr_qn_pre_` prefixes.
143: .seealso: `SNESNEWTONTR`, `SNESNewtonTRQNType`, `MATLMVM`
144: @*/
145: PetscErrorCode SNESNewtonTRSetQNType(SNES snes, SNESNewtonTRQNType use)
146: {
147: PetscBool flg;
149: PetscFunctionBegin;
152: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
153: if (flg) {
154: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
156: tr->qn = use;
157: }
158: PetscFunctionReturn(PETSC_SUCCESS);
159: }
161: /*@
162: SNESNewtonTRSetFallbackType - Set the type of fallback to use if the solution of the trust region subproblem is outside the radius
164: Input Parameters:
165: + snes - the nonlinear solver object
166: - ftype - the fallback type, see `SNESNewtonTRFallbackType`
168: Level: intermediate
170: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRPreCheck()`, `SNESNewtonTRGetPreCheck()`, `SNESNewtonTRSetPreCheck()`,
171: `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()`
172: @*/
173: PetscErrorCode SNESNewtonTRSetFallbackType(SNES snes, SNESNewtonTRFallbackType ftype)
174: {
175: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
176: PetscBool flg;
178: PetscFunctionBegin;
181: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
182: if (flg) tr->fallback = ftype;
183: PetscFunctionReturn(PETSC_SUCCESS);
184: }
186: /*@
187: SNESNewtonTRSetPreCheck - Sets a user function that is called before the search step has been determined.
188: Allows the user a chance to change or override the trust region decision.
190: Logically Collective
192: Input Parameters:
193: + snes - the nonlinear solver object
194: . func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRPreCheck()`
195: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
197: Calling sequence of `func`:
198: + snes - the nonlinear solver object
199: . X - the current solution value
200: . Y - the tentative update step
201: . changed - output, flag indicating `Y` has been changed by the pre-check
202: - ctx - the optional application context
204: Level: intermediate
206: Note:
207: This function is called BEFORE the function evaluation within the solver.
209: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRPreCheck()`, `SNESNewtonTRGetPreCheck()`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()`
210: @*/
211: PetscErrorCode SNESNewtonTRSetPreCheck(SNES snes, PetscErrorCode (*func)(SNES snes, Vec X, Vec Y, PetscBool *changed, PetscCtx ctx), PetscCtx ctx)
212: {
213: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
214: PetscBool flg;
216: PetscFunctionBegin;
218: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
219: if (flg) {
220: if (func) tr->precheck = func;
221: if (ctx) tr->precheckctx = ctx;
222: }
223: PetscFunctionReturn(PETSC_SUCCESS);
224: }
226: /*@
227: SNESNewtonTRGetPreCheck - Gets the pre-check function
229: Not Collective
231: Input Parameter:
232: . snes - the nonlinear solver context
234: Output Parameters:
235: + func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRPreCheck()`
236: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
238: Calling sequence of `func`:
239: + snes - the nonlinear solver object
240: . X - the current solution value
241: . Y - the tentative update step
242: . changed - output, flag indicating `Y` has been changed by the pre-check
243: - ctx - the optional application context
245: Level: intermediate
247: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRPreCheck()`
248: @*/
249: PetscErrorCode SNESNewtonTRGetPreCheck(SNES snes, PetscErrorCode (**func)(SNES snes, Vec X, Vec Y, PetscBool *changed, PetscCtx ctx), PetscCtxRt ctx)
250: {
251: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
252: PetscBool flg;
254: PetscFunctionBegin;
256: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
257: PetscAssert(flg, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)snes)->type_name);
258: if (func) *func = tr->precheck;
259: if (ctx) *(void **)ctx = tr->precheckctx;
260: PetscFunctionReturn(PETSC_SUCCESS);
261: }
263: /*@
264: SNESNewtonTRSetPostCheck - Sets a user function that is called after the search step has been determined but before the next
265: function evaluation. Allows the user a chance to change or override the internal decision of the solver
267: Logically Collective
269: Input Parameters:
270: + snes - the nonlinear solver object
271: . func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRPostCheck()`
272: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
274: Calling sequence of `func`:
275: + snes - the nonlinear solver object
276: . X - the current solution value
277: . Y - the tentative update step
278: . W - the tentative new solution value
279: . changed_Y - output, flag indicated `Y` has been changed by the post-check
280: . changed_W - output, flag indicated `W` has been changed by the post-check
281: - ctx - the optional application context
283: Level: intermediate
285: Note:
286: This function is called BEFORE the function evaluation within the solver while the function set in
287: `SNESLineSearchSetPostCheck()` is called AFTER the function evaluation.
289: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRPostCheck()`, `SNESNewtonTRGetPostCheck()`, `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRGetPreCheck()`
290: @*/
291: PetscErrorCode SNESNewtonTRSetPostCheck(SNES snes, PetscErrorCode (*func)(SNES snes, Vec X, Vec Y, Vec W, PetscBool *changed_Y, PetscBool *changed_W, PetscCtx ctx), PetscCtx ctx)
292: {
293: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
294: PetscBool flg;
296: PetscFunctionBegin;
298: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
299: if (flg) {
300: if (func) tr->postcheck = func;
301: if (ctx) tr->postcheckctx = ctx;
302: }
303: PetscFunctionReturn(PETSC_SUCCESS);
304: }
306: /*@
307: SNESNewtonTRGetPostCheck - Gets the post-check function
309: Not Collective
311: Input Parameter:
312: . snes - the nonlinear solver context
314: Output Parameters:
315: + func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRPostCheck()`
316: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
318: Calling sequence of `func`:
319: + snes - the nonlinear solver object
320: . X - the current solution value
321: . Y - the tentative update step
322: . W - the tentative new solution value
323: . changed_Y - output, flag indicated `Y` has been changed by the post-check
324: . changed_W - output, flag indicated `W` has been changed by the post-check
325: - ctx - the optional application context
327: Level: intermediate
329: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRPostCheck()`
330: @*/
331: PetscErrorCode SNESNewtonTRGetPostCheck(SNES snes, PetscErrorCode (**func)(SNES snes, Vec X, Vec Y, Vec W, PetscBool *changed_Y, PetscBool *changed_W, PetscCtx ctx), PetscCtxRt ctx)
332: {
333: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
334: PetscBool flg;
336: PetscFunctionBegin;
338: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
339: PetscAssert(flg, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)snes)->type_name);
340: if (func) *func = tr->postcheck;
341: if (ctx) *(void **)ctx = tr->postcheckctx;
342: PetscFunctionReturn(PETSC_SUCCESS);
343: }
345: /*@
346: SNESNewtonTRPreCheck - Runs the precheck routine
348: Logically Collective
350: Input Parameters:
351: + snes - the solver
352: . X - The last solution
353: - Y - The step direction
355: Output Parameter:
356: . changed_Y - Indicator that the step direction `Y` has been changed.
358: Level: intermediate
360: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRGetPreCheck()`, `SNESNewtonTRPostCheck()`
361: @*/
362: PetscErrorCode SNESNewtonTRPreCheck(SNES snes, Vec X, Vec Y, PetscBool *changed_Y)
363: {
364: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
365: PetscBool flg;
367: PetscFunctionBegin;
369: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
370: PetscAssert(flg, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)snes)->type_name);
371: *changed_Y = PETSC_FALSE;
372: if (tr->precheck) {
373: PetscCall((*tr->precheck)(snes, X, Y, changed_Y, tr->precheckctx));
375: }
376: PetscFunctionReturn(PETSC_SUCCESS);
377: }
379: /*@
380: SNESNewtonTRPostCheck - Runs the postcheck routine
382: Logically Collective
384: Input Parameters:
385: + snes - the solver
386: . X - The last solution
387: . Y - The full step direction
388: - W - The updated solution, W = X - Y
390: Output Parameters:
391: + changed_Y - indicator if step has been changed
392: - changed_W - Indicator if the new candidate solution W has been changed.
394: Note:
395: If Y is changed then W is recomputed as X - Y
397: Level: intermediate
399: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()`, `SNESNewtonTRPreCheck()`
400: @*/
401: PetscErrorCode SNESNewtonTRPostCheck(SNES snes, Vec X, Vec Y, Vec W, PetscBool *changed_Y, PetscBool *changed_W)
402: {
403: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
404: PetscBool flg;
406: PetscFunctionBegin;
408: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
409: PetscAssert(flg, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)snes)->type_name);
410: *changed_Y = PETSC_FALSE;
411: *changed_W = PETSC_FALSE;
412: if (tr->postcheck) {
413: PetscCall((*tr->postcheck)(snes, X, Y, W, changed_Y, changed_W, tr->postcheckctx));
416: }
417: PetscFunctionReturn(PETSC_SUCCESS);
418: }
420: /* stable implementation of roots of a*x^2 + b*x + c = 0 */
421: static inline void PetscQuadraticRoots(PetscReal a, PetscReal b, PetscReal c, PetscReal *xm, PetscReal *xp)
422: {
423: PetscReal temp = -0.5 * (b + PetscCopysignReal(1.0, b) * PetscSqrtReal(b * b - 4 * a * c));
424: PetscReal x1 = temp / a;
425: PetscReal x2 = c / temp;
426: *xm = PetscMin(x1, x2);
427: *xp = PetscMax(x1, x2);
428: }
430: /* Computes the quadratic model difference */
431: static PetscErrorCode SNESNewtonTRQuadraticDelta(SNES snes, Mat J, PetscBool has_objective, Vec Y, Vec GradF, Vec W, PetscReal *yTHy_, PetscReal *gTy_, PetscReal *deltaqm)
432: {
433: PetscReal yTHy, gTy;
435: PetscFunctionBegin;
436: PetscCall(MatMult(J, Y, W));
437: if (has_objective) PetscCall(VecDotRealPart(Y, W, &yTHy));
438: else PetscCall(VecDotRealPart(W, W, &yTHy)); /* Gauss-Newton approximation J^t * J */
439: PetscCall(VecDotRealPart(GradF, Y, &gTy));
440: *deltaqm = -(-gTy + 0.5 * yTHy); /* difference in quadratic model, -gTy because SNES solves it this way */
441: if (yTHy_) *yTHy_ = yTHy;
442: if (gTy_) *gTy_ = gTy;
443: PetscFunctionReturn(PETSC_SUCCESS);
444: }
446: /* Computes the new objective given X = Xk, Y = direction
447: W work vector, on output W = X - Y
448: G work vector, on output G = SNESFunction(W) */
449: static PetscErrorCode SNESNewtonTRObjective(SNES snes, PetscBool has_objective, Vec X, Vec Y, Vec W, Vec G, PetscReal *gnorm, PetscReal *fkp1)
450: {
451: PetscBool changed_Y, changed_W;
453: PetscFunctionBegin;
454: /* TODO: we can add a linesearch here */
455: PetscCall(SNESNewtonTRPreCheck(snes, X, Y, &changed_Y));
456: PetscCall(VecWAXPY(W, -1.0, Y, X)); /* Xkp1 */
457: PetscCall(SNESNewtonTRPostCheck(snes, X, Y, W, &changed_Y, &changed_W));
458: if (changed_Y && !changed_W) PetscCall(VecWAXPY(W, -1.0, Y, X));
460: PetscCall(SNESComputeFunction(snes, W, G)); /* F(Xkp1) = G */
461: PetscCall(VecNorm(G, NORM_2, gnorm));
462: SNESCheckFunctionDomainError(snes, *gnorm);
463: if (has_objective) {
464: PetscCall(SNESComputeObjective(snes, W, fkp1));
465: SNESCheckObjectiveDomainError(snes, *fkp1);
466: } else *fkp1 = 0.5 * PetscSqr(*gnorm);
467: PetscFunctionReturn(PETSC_SUCCESS);
468: }
470: static PetscErrorCode SNESSetUpQN_NEWTONTR(SNES snes)
471: {
472: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
474: PetscFunctionBegin;
475: PetscCall(MatDestroy(&tr->qnB));
476: PetscCall(MatDestroy(&tr->qnB_pre));
477: if (tr->qn) {
478: PetscInt n, N;
479: const char *optionsprefix;
480: Mat B;
482: PetscCall(MatCreate(PetscObjectComm((PetscObject)snes), &B));
483: PetscCall(SNESGetOptionsPrefix(snes, &optionsprefix));
484: PetscCall(MatSetOptionsPrefix(B, "snes_tr_qn_"));
485: PetscCall(MatAppendOptionsPrefix(B, optionsprefix));
486: PetscCall(MatSetType(B, MATLMVMBFGS));
487: PetscCall(VecGetLocalSize(snes->vec_sol, &n));
488: PetscCall(VecGetSize(snes->vec_sol, &N));
489: PetscCall(MatSetSizes(B, n, n, N, N));
490: PetscCall(MatSetUp(B));
491: PetscCall(MatSetFromOptions(B));
492: PetscCall(MatLMVMAllocate(B, snes->vec_sol, snes->vec_func));
493: tr->qnB = B;
494: if (tr->qn == SNES_TR_QN_DIFFERENT) {
495: PetscCall(MatCreate(PetscObjectComm((PetscObject)snes), &B));
496: PetscCall(SNESGetOptionsPrefix(snes, &optionsprefix));
497: PetscCall(MatSetOptionsPrefix(B, "snes_tr_qn_pre_"));
498: PetscCall(MatAppendOptionsPrefix(B, optionsprefix));
499: PetscCall(MatSetType(B, MATLMVMBFGS));
500: PetscCall(MatSetSizes(B, n, n, N, N));
501: PetscCall(MatSetUp(B));
502: PetscCall(MatSetFromOptions(B));
503: PetscCall(MatLMVMAllocate(B, snes->vec_sol, snes->vec_func));
504: tr->qnB_pre = B;
505: } else {
506: PetscCall(PetscObjectReference((PetscObject)tr->qnB));
507: tr->qnB_pre = tr->qnB;
508: }
509: }
510: PetscFunctionReturn(PETSC_SUCCESS);
511: }
513: /*
514: SNESSolve_NEWTONTR - Implements Newton's Method with trust-region subproblem and adds dogleg Cauchy
515: (Steepest Descent direction) step and direction if the trust region is not satisfied for solving system of
516: nonlinear equations
518: */
519: static PetscErrorCode SNESSolve_NEWTONTR(SNES snes)
520: {
521: SNES_NEWTONTR *neP = (SNES_NEWTONTR *)snes->data;
522: Vec X, F, Y, G, W, GradF, YU, Yc;
523: PetscInt maxits, lits;
524: PetscReal rho, fnorm, gnorm = 0.0, xnorm = 0.0, delta, ynorm;
525: PetscReal fk, fkp1, deltaqm = 0.0, gTy = 0.0, yTHy = 0.0;
526: PetscReal auk, tauk, gfnorm, gfnorm_k, ycnorm, gTBg, objmin = 0.0, beta_k = 1.0;
527: PC pc;
528: Mat J, Jp;
529: PetscBool already_done = PETSC_FALSE, on_boundary = PETSC_FALSE, use_cauchy;
530: PetscBool clear_converged_test, step_ok = PETSC_FALSE, has_objective;
531: SNES_TR_KSPConverged_Ctx *ctx;
532: void *convctx;
533: SNESObjectiveFn *objective;
534: KSPConvergenceTestFn *convtest;
535: PetscCtxDestroyFn *convdestroy;
537: PetscFunctionBegin;
538: PetscCall(SNESGetObjective(snes, &objective, NULL));
539: has_objective = objective ? PETSC_TRUE : PETSC_FALSE;
541: maxits = snes->max_its; /* maximum number of iterations */
542: X = snes->vec_sol; /* solution vector */
543: F = snes->vec_func; /* residual vector */
544: Y = snes->vec_sol_update; /* update vector */
545: G = snes->work[0]; /* updated residual */
546: W = snes->work[1]; /* temporary vector */
547: GradF = !has_objective ? snes->work[2] : snes->vec_func; /* grad f = J^T F */
548: YU = snes->work[3]; /* work vector for dogleg method */
549: Yc = snes->work[4]; /* Cauchy point */
551: 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);
553: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
554: snes->iter = 0;
555: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
557: /* setup QN matrices if needed */
558: PetscCall(SNESSetUpQN_NEWTONTR(snes));
560: /* Set the linear stopping criteria to use the More' trick if needed */
561: clear_converged_test = PETSC_FALSE;
562: PetscCall(SNESGetKSP(snes, &snes->ksp));
563: PetscCall(KSPGetConvergenceTest(snes->ksp, &convtest, &convctx, &convdestroy));
564: if (convtest != SNESTR_KSPConverged_Private) {
565: clear_converged_test = PETSC_TRUE;
566: PetscCall(PetscNew(&ctx));
567: ctx->snes = snes;
568: PetscCall(KSPGetAndClearConvergenceTest(snes->ksp, &ctx->convtest, &ctx->convctx, &ctx->convdestroy));
569: PetscCall(KSPSetConvergenceTest(snes->ksp, SNESTR_KSPConverged_Private, ctx, SNESTR_KSPConverged_Destroy));
570: PetscCall(PetscInfo(snes, "Using Krylov convergence test SNESTR_KSPConverged_Private\n"));
571: }
573: if (!snes->vec_func_init_set) {
574: PetscCall(SNESComputeFunction(snes, X, F)); /* F(X) */
575: } else snes->vec_func_init_set = PETSC_FALSE;
577: PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- || F || */
578: SNESCheckFunctionDomainError(snes, fnorm);
579: PetscCall(VecNorm(X, NORM_2, &xnorm)); /* xnorm <- || X || */
581: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
582: snes->norm = fnorm;
583: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
584: delta = neP->delta0;
585: neP->delta = delta;
586: PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));
588: /* test convergence */
589: PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm));
590: PetscCall(SNESMonitor(snes, 0, fnorm));
591: if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
593: if (has_objective) {
594: PetscCall(SNESComputeObjective(snes, X, &fk));
595: SNESCheckObjectiveDomainError(snes, fk);
596: } else fk = 0.5 * PetscSqr(fnorm); /* obj(x) = 0.5 * ||F(x)||^2 */
598: /* hook state vector to BFGS preconditioner */
599: PetscCall(KSPGetPC(snes->ksp, &pc));
600: PetscCall(PCLMVMSetUpdateVec(pc, X));
602: if (neP->kmdc) PetscCall(KSPSetComputeEigenvalues(snes->ksp, PETSC_TRUE));
604: while (snes->iter < maxits) {
605: /* calculating Jacobian and GradF of minimization function only once */
606: if (!already_done) {
607: /* Call general purpose update function */
608: PetscTryTypeMethod(snes, update, snes->iter);
610: /* apply the nonlinear preconditioner */
611: if (snes->npc && snes->npcside == PC_RIGHT) {
612: SNESConvergedReason reason;
614: PetscCall(SNESSetInitialFunction(snes->npc, F));
615: PetscCall(PetscLogEventBegin(SNES_NPCSolve, snes->npc, X, snes->vec_rhs, 0));
616: PetscCall(SNESSolve(snes->npc, snes->vec_rhs, X));
617: PetscCall(PetscLogEventEnd(SNES_NPCSolve, snes->npc, X, snes->vec_rhs, 0));
618: PetscCall(SNESGetConvergedReason(snes->npc, &reason));
619: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
620: snes->reason = SNES_DIVERGED_INNER;
621: PetscFunctionReturn(PETSC_SUCCESS);
622: }
623: // XXX
624: PetscCall(SNESGetNPCFunction(snes, F, &fnorm));
625: }
627: /* Jacobian */
628: J = NULL;
629: Jp = NULL;
630: if (!neP->qnB) {
631: PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
632: SNESCheckJacobianDomainError(snes);
633: J = snes->jacobian;
634: Jp = snes->jacobian_pre;
635: } else { /* QN model */
636: PetscCall(SNESComputeJacobian_MATLMVM(snes, X, neP->qnB, neP->qnB_pre, NULL));
637: J = neP->qnB;
638: Jp = neP->qnB_pre;
639: }
640: SNESCheckJacobianDomainError(snes);
642: /* objective function */
643: PetscCall(VecNorm(F, NORM_2, &fnorm));
644: SNESCheckFunctionDomainError(snes, fnorm);
645: if (has_objective) {
646: PetscCall(SNESComputeObjective(snes, X, &fk));
647: SNESCheckObjectiveDomainError(snes, fk);
648: } else fk = 0.5 * PetscSqr(fnorm); /* obj(x) = 0.5 * ||F(x)||^2 */
650: /* GradF */
651: if (has_objective) gfnorm = fnorm;
652: else {
653: PetscCall(MatMultTranspose(J, F, GradF)); /* grad f = J^T F */
654: PetscCall(VecNorm(GradF, NORM_2, &gfnorm));
655: }
656: PetscCall(VecNorm(GradF, neP->norm, &gfnorm_k));
657: }
658: already_done = PETSC_TRUE;
660: /* solve trust-region subproblem */
662: /* first compute Cauchy Point */
663: PetscCall(MatMult(J, GradF, W));
664: if (has_objective) PetscCall(VecDotRealPart(GradF, W, &gTBg));
665: else PetscCall(VecDotRealPart(W, W, &gTBg)); /* B = J^t * J */
666: /* Eqs 4.11 and 4.12 in Nocedal and Wright 2nd Edition (4.7 and 4.8 in 1st Edition) */
667: auk = delta / gfnorm_k;
668: if (gTBg < 0.0) tauk = 1.0;
669: else tauk = PetscMin(gfnorm * gfnorm * gfnorm_k / (delta * gTBg), 1);
670: auk *= tauk;
671: ycnorm = auk * gfnorm;
672: PetscCall(VecAXPBY(Yc, auk, 0.0, GradF));
674: use_cauchy = (PetscBool)(tauk == 1.0 && has_objective);
675: if (!use_cauchy) {
676: KSPConvergedReason reason;
678: /* sufficient decrease (see 6.3.27 in Conn, Gould, Toint "Trust Region Methods")
679: beta_k the largest eigenvalue of the Hessian. Here we use the previous estimated value */
680: objmin = -neP->kmdc * gnorm * PetscMin(gnorm / beta_k, delta);
681: PetscCall(KSPCGSetObjectiveTarget(snes->ksp, objmin));
683: /* specify radius if looking for Newton step and trust region norm is the l2 norm */
684: PetscCall(KSPCGSetRadius(snes->ksp, neP->fallback == SNES_TR_FALLBACK_NEWTON && neP->norm == NORM_2 ? delta : 0.0));
685: PetscCall(KSPSetOperators(snes->ksp, J, Jp));
686: PetscCall(KSPSolve(snes->ksp, F, Y));
687: SNESCheckKSPSolve(snes);
688: PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
689: PetscCall(KSPGetConvergedReason(snes->ksp, &reason));
690: on_boundary = (PetscBool)(reason == KSP_CONVERGED_STEP_LENGTH);
691: PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits));
692: if (neP->kmdc) { /* update estimated Hessian largest eigenvalue */
693: PetscReal emax, emin;
694: PetscCall(KSPComputeExtremeSingularValues(snes->ksp, &emax, &emin));
695: if (emax > 0.0) beta_k = emax + 1;
696: }
697: } else { /* Cauchy point is on the boundary, accept it */
698: on_boundary = PETSC_TRUE;
699: PetscCall(VecCopy(Yc, Y));
700: PetscCall(PetscInfo(snes, "CP evaluated on boundary. delta: %g, ycnorm: %g, gTBg: %g\n", (double)delta, (double)ycnorm, (double)gTBg));
701: }
702: PetscCall(VecNorm(Y, neP->norm, &ynorm));
704: /* decide what to do when the update is outside of trust region */
705: if (!use_cauchy && (ynorm > delta || ynorm == 0.0)) {
706: SNESNewtonTRFallbackType fallback = ynorm > 0.0 ? neP->fallback : SNES_TR_FALLBACK_CAUCHY;
708: PetscCheck(neP->norm == NORM_2 || fallback != SNES_TR_FALLBACK_DOGLEG, PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "DOGLEG without l2 norm not implemented");
709: switch (fallback) {
710: case SNES_TR_FALLBACK_NEWTON:
711: auk = delta / ynorm;
712: PetscCall(VecScale(Y, auk));
713: PetscCall(PetscInfo(snes, "SN evaluated. delta: %g, ynorm: %g\n", (double)delta, (double)ynorm));
714: break;
715: case SNES_TR_FALLBACK_CAUCHY:
716: case SNES_TR_FALLBACK_DOGLEG:
717: if (fallback == SNES_TR_FALLBACK_CAUCHY || gTBg <= 0.0) {
718: PetscCall(VecCopy(Yc, Y));
719: PetscCall(PetscInfo(snes, "CP evaluated. delta: %g, ynorm: %g, ycnorm: %g, gTBg: %g\n", (double)delta, (double)ynorm, (double)ycnorm, (double)gTBg));
720: } else { /* take linear combination of Cauchy and Newton direction and step */
721: auk = gfnorm * gfnorm / gTBg;
722: if (gfnorm_k * auk >= delta) { /* first leg: Cauchy point outside of trust region */
723: PetscCall(VecAXPBY(Y, delta / gfnorm_k, 0.0, GradF));
724: PetscCall(PetscInfo(snes, "CP evaluated (outside region). delta: %g, ynorm: %g, ycnorm: %g\n", (double)delta, (double)ynorm, (double)ycnorm));
725: } else { /* second leg */
726: PetscReal c0, c1, c2, tau = 0.0, tpos, tneg;
727: PetscBool noroots;
729: /* Find solutions of (Eq. 4.16 in Nocedal and Wright)
730: ||p_U + lambda * (p_B - p_U)||^2 - delta^2 = 0,
731: where p_U the Cauchy direction, p_B the Newton direction */
732: PetscCall(VecAXPBY(YU, auk, 0.0, GradF));
733: PetscCall(VecAXPY(Y, -1.0, YU));
734: PetscCall(VecNorm(Y, NORM_2, &c0));
735: PetscCall(VecDotRealPart(YU, Y, &c1));
736: c0 = PetscSqr(c0);
737: c2 = PetscSqr(ycnorm) - PetscSqr(delta);
738: PetscQuadraticRoots(c0, 2 * c1, c2, &tneg, &tpos);
740: /* In principle the DL strategy as a unique solution in [0,1]
741: here we check that for some reason we numerically failed
742: to compute it. In that case, we use the Cauchy point */
743: noroots = PetscIsInfOrNanReal(tneg);
744: if (!noroots) {
745: if (tpos > 1) {
746: if (tneg >= 0 && tneg <= 1) {
747: tau = tneg;
748: } else noroots = PETSC_TRUE;
749: } else if (tpos >= 0) {
750: tau = tpos;
751: } else noroots = PETSC_TRUE;
752: }
753: if (noroots) { /* No roots, select Cauchy point */
754: PetscCall(VecCopy(Yc, Y));
755: } else {
756: PetscCall(VecAXPBY(Y, 1.0, tau, YU));
757: }
758: 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));
759: }
760: }
761: break;
762: default:
763: SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "Unknown fallback mode");
764: break;
765: }
766: }
768: /* compute the quadratic model difference */
769: PetscCall(SNESNewtonTRQuadraticDelta(snes, J, has_objective, Y, GradF, W, &yTHy, &gTy, &deltaqm));
771: /* Compute new objective function */
772: PetscCall(SNESNewtonTRObjective(snes, has_objective, X, Y, W, G, &gnorm, &fkp1));
773: if (PetscIsInfOrNanReal(fkp1)) rho = neP->eta1;
774: else {
775: if (deltaqm > 0.0) rho = (fk - fkp1) / deltaqm; /* actual improvement over predicted improvement */
776: else rho = neP->eta1; /* no reduction in quadratic model, step must be rejected */
777: }
779: PetscCall(VecNorm(Y, neP->norm, &ynorm));
780: PetscCall(PetscInfo(snes, "rho=%g, delta=%g, fk=%g, fkp1=%g, deltaqm=%g, gTy=%g, yTHy=%g, ynormk=%g, gnorm=%g\n", (double)rho, (double)delta, (double)fk, (double)fkp1, (double)deltaqm, (double)gTy, (double)yTHy, (double)ynorm, (double)gnorm));
782: /* update the size of the trust region */
783: if (rho < neP->eta2) delta *= neP->t1; /* shrink the region */
784: else if (rho > neP->eta3 && on_boundary) delta *= neP->t2; /* expand the region */
785: delta = PetscMin(delta, neP->deltaM); /* but not greater than deltaM */
787: /* log 2-norm of update for monitoring routines */
788: PetscCall(VecNorm(Y, NORM_2, &ynorm));
790: /* decide on new step */
791: neP->delta = delta;
792: if (rho > neP->eta1) {
793: step_ok = PETSC_TRUE;
794: } else {
795: step_ok = PETSC_FALSE;
796: /* check convergence */
797: PetscCall(SNESTR_Converged_Private(snes, xnorm, ynorm, gnorm));
798: if (snes->reason < 0) {
799: snes->numFailures++;
800: break;
801: } else if (!snes->reason) {
802: snes->numFailures++;
803: PetscCall(PetscInfo(snes, "Trying again in smaller region\n"));
804: } else step_ok = PETSC_TRUE;
805: }
806: if (step_ok) {
807: /* Update function values */
808: already_done = PETSC_FALSE;
809: fnorm = gnorm;
810: fk = fkp1;
812: /* New residual and linearization point */
813: PetscCall(VecCopy(G, F));
814: PetscCall(VecCopy(W, X));
816: /* Monitor convergence */
817: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
818: snes->iter++;
819: snes->norm = fnorm;
820: snes->xnorm = xnorm;
821: snes->ynorm = ynorm;
822: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
823: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits));
825: /* Test for convergence, xnorm = || X || */
826: PetscCall(VecNorm(X, NORM_2, &xnorm));
827: PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm));
828: PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
829: if (snes->reason) break;
830: }
831: }
833: if (clear_converged_test) {
834: PetscCall(KSPGetAndClearConvergenceTest(snes->ksp, &ctx->convtest, &ctx->convctx, &ctx->convdestroy));
835: PetscCall(PetscFree(ctx));
836: PetscCall(KSPSetConvergenceTest(snes->ksp, convtest, convctx, convdestroy));
837: }
838: PetscFunctionReturn(PETSC_SUCCESS);
839: }
841: static PetscErrorCode SNESSetUp_NEWTONTR(SNES snes)
842: {
843: PetscFunctionBegin;
844: PetscCall(SNESSetWorkVecs(snes, 5));
845: PetscCall(SNESSetUpMatrices(snes));
846: PetscFunctionReturn(PETSC_SUCCESS);
847: }
849: static PetscErrorCode SNESReset_NEWTONTR(SNES snes)
850: {
851: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
853: PetscFunctionBegin;
854: PetscCall(MatDestroy(&tr->qnB));
855: PetscCall(MatDestroy(&tr->qnB_pre));
856: PetscFunctionReturn(PETSC_SUCCESS);
857: }
859: static PetscErrorCode SNESDestroy_NEWTONTR(SNES snes)
860: {
861: PetscFunctionBegin;
862: PetscCall(SNESReset_NEWTONTR(snes));
863: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonTRSetTolerances_C", NULL));
864: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonTRGetTolerances_C", NULL));
865: PetscCall(PetscFree(snes->data));
866: PetscFunctionReturn(PETSC_SUCCESS);
867: }
869: static PetscErrorCode SNESSetFromOptions_NEWTONTR(SNES snes, PetscOptionItems PetscOptionsObject)
870: {
871: SNES_NEWTONTR *ctx = (SNES_NEWTONTR *)snes->data;
872: SNESNewtonTRQNType qn;
873: SNESNewtonTRFallbackType fallback;
874: NormType norm;
875: PetscBool flg;
877: PetscFunctionBegin;
878: PetscOptionsHeadBegin(PetscOptionsObject, "SNES trust region options for nonlinear equations");
879: PetscCall(PetscOptionsDeprecated("-snes_tr_deltaM", "-snes_tr_deltamax", "3.22", NULL));
880: PetscCall(PetscOptionsReal("-snes_tr_eta1", "eta1", "SNESNewtonTRSetUpdateParameters", ctx->eta1, &ctx->eta1, NULL));
881: PetscCall(PetscOptionsReal("-snes_tr_eta2", "eta2", "SNESNewtonTRSetUpdateParameters", ctx->eta2, &ctx->eta2, NULL));
882: PetscCall(PetscOptionsReal("-snes_tr_eta3", "eta3", "SNESNewtonTRSetUpdateParameters", ctx->eta3, &ctx->eta3, NULL));
883: PetscCall(PetscOptionsReal("-snes_tr_t1", "t1", "SNESNewtonTRSetUpdateParameters", ctx->t1, &ctx->t1, NULL));
884: PetscCall(PetscOptionsReal("-snes_tr_t2", "t2", "SNESNewtonTRSetUpdateParameters", ctx->t2, &ctx->t2, NULL));
885: PetscCall(PetscOptionsReal("-snes_tr_delta0", "Initial trust region size", "SNESNewtonTRSetTolerances", ctx->delta0, &ctx->delta0, NULL));
886: PetscCall(PetscOptionsReal("-snes_tr_deltamin", "Minimum allowed trust region size", "SNESNewtonTRSetTolerances", ctx->deltam, &ctx->deltam, NULL));
887: PetscCall(PetscOptionsReal("-snes_tr_deltamax", "Maximum allowed trust region size", "SNESNewtonTRSetTolerances", ctx->deltaM, &ctx->deltaM, NULL));
888: PetscCall(PetscOptionsReal("-snes_tr_kmdc", "sufficient decrease parameter", "None", ctx->kmdc, &ctx->kmdc, NULL));
890: fallback = ctx->fallback;
891: 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));
892: if (flg) PetscCall(SNESNewtonTRSetFallbackType(snes, fallback));
894: qn = ctx->qn;
895: PetscCall(PetscOptionsEnum("-snes_tr_qn", "Use Quasi-Newton approximations for the model", "SNESNewtonTRSetQNType", SNESNewtonTRQNTypes, (PetscEnum)qn, (PetscEnum *)&qn, &flg));
896: if (flg) PetscCall(SNESNewtonTRSetQNType(snes, qn));
898: norm = ctx->norm;
899: PetscCall(PetscOptionsEnum("-snes_tr_norm_type", "Type of norm for trust region bounds", "SNESNewtonTRSetNormType", NormTypes, (PetscEnum)norm, (PetscEnum *)&norm, &flg));
900: if (flg) PetscCall(SNESNewtonTRSetNormType(snes, norm));
902: PetscOptionsHeadEnd();
903: PetscFunctionReturn(PETSC_SUCCESS);
904: }
906: static PetscErrorCode SNESView_NEWTONTR(SNES snes, PetscViewer viewer)
907: {
908: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
909: PetscBool isascii;
911: PetscFunctionBegin;
912: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
913: if (isascii) {
914: PetscCall(PetscViewerASCIIPrintf(viewer, " Trust region parameters:\n"));
915: PetscCall(PetscViewerASCIIPrintf(viewer, " eta1=%g, eta2=%g, eta3=%g\n", (double)tr->eta1, (double)tr->eta2, (double)tr->eta3));
916: PetscCall(PetscViewerASCIIPrintf(viewer, " t1=%g, t2=%g\n", (double)tr->t1, (double)tr->t2));
917: PetscCall(PetscViewerASCIIPrintf(viewer, " delta_min=%g, delta_0=%g, delta_max=%g\n", (double)tr->deltam, (double)tr->delta0, (double)tr->deltaM));
918: PetscCall(PetscViewerASCIIPrintf(viewer, " kmdc=%g\n", (double)tr->kmdc));
919: PetscCall(PetscViewerASCIIPrintf(viewer, " fallback=%s\n", SNESNewtonTRFallbackTypes[tr->fallback]));
920: if (tr->qn) PetscCall(PetscViewerASCIIPrintf(viewer, " qn=%s\n", SNESNewtonTRQNTypes[tr->qn]));
921: if (tr->norm != NORM_2) PetscCall(PetscViewerASCIIPrintf(viewer, " norm=%s\n", NormTypes[tr->norm]));
922: }
923: PetscFunctionReturn(PETSC_SUCCESS);
924: }
926: /*@
927: SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance.
929: Logically Collective
931: Input Parameters:
932: + snes - the `SNES` context
933: - tol - tolerance
935: Level: deprecated
937: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTR`, `SNESSetTolerances()`
938: @*/
939: PetscErrorCode SNESSetTrustRegionTolerance(SNES snes, PetscReal tol)
940: {
941: return SNESNewtonTRSetTolerances(snes, tol, PETSC_CURRENT, PETSC_CURRENT);
942: }
944: /*@
945: SNESNewtonTRSetTolerances - Sets the trust region parameter tolerances.
947: Logically Collective
949: Input Parameters:
950: + snes - the `SNES` context
951: . delta_min - minimum allowed trust region size
952: . delta_max - maximum allowed trust region size
953: - delta_0 - initial trust region size
955: Options Database Key:
956: + -snes_tr_deltamin tol - Set minimum size
957: . -snes_tr_deltamax tol - Set maximum size
958: - -snes_tr_delta0 tol - Set initial size
960: Note:
961: Use `PETSC_DETERMINE` to use the default value for the given `SNES`.
962: Use `PETSC_CURRENT` to retain a value.
964: Fortran Note:
965: Use `PETSC_DETERMINE_REAL`, `PETSC_CURRENT_REAL`
967: Level: intermediate
969: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTR`, `SNESNewtonTRGetTolerances()`
970: @*/
971: PetscErrorCode SNESNewtonTRSetTolerances(SNES snes, PetscReal delta_min, PetscReal delta_max, PetscReal delta_0)
972: {
973: PetscFunctionBegin;
978: PetscTryMethod(snes, "SNESNewtonTRSetTolerances_C", (SNES, PetscReal, PetscReal, PetscReal), (snes, delta_min, delta_max, delta_0));
979: PetscFunctionReturn(PETSC_SUCCESS);
980: }
982: /*@
983: SNESNewtonTRGetTolerances - Gets the trust region parameter tolerances.
985: Not Collective
987: Input Parameter:
988: . snes - the `SNES` context
990: Output Parameters:
991: + delta_min - minimum allowed trust region size or `NULL`
992: . delta_max - maximum allowed trust region size or `NULL`
993: - delta_0 - initial trust region size or `NULL`
995: Level: intermediate
997: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTR`, `SNESNewtonTRSetTolerances()`
998: @*/
999: PetscErrorCode SNESNewtonTRGetTolerances(SNES snes, PetscReal *delta_min, PetscReal *delta_max, PetscReal *delta_0)
1000: {
1001: PetscFunctionBegin;
1003: if (delta_min) PetscAssertPointer(delta_min, 2);
1004: if (delta_max) PetscAssertPointer(delta_max, 3);
1005: if (delta_0) PetscAssertPointer(delta_0, 4);
1006: PetscUseMethod(snes, "SNESNewtonTRGetTolerances_C", (SNES, PetscReal *, PetscReal *, PetscReal *), (snes, delta_min, delta_max, delta_0));
1007: PetscFunctionReturn(PETSC_SUCCESS);
1008: }
1010: /*@
1011: SNESNewtonTRSetUpdateParameters - Sets the trust region update parameters.
1013: Logically Collective
1015: Input Parameters:
1016: + snes - the `SNES` context
1017: . eta1 - acceptance tolerance
1018: . eta2 - shrinking tolerance
1019: . eta3 - enlarging tolerance
1020: . t1 - shrink factor
1021: - t2 - enlarge factor
1023: Options Database Key:
1024: + -snes_tr_eta1 tol - Set `eta1`
1025: . -snes_tr_eta2 tol - Set `eta2`
1026: . -snes_tr_eta3 tol - Set `eta3`
1027: . -snes_tr_t1 tol - Set `t1`
1028: - -snes_tr_t2 tol - Set `t2`
1030: Notes:
1031: Given the ratio $\rho = \frac{f(x_k) - f(x_k+s_k)}{m(0) - m(s_k)}$, with $x_k$ the current iterate,
1032: $s_k$ the computed step, $f$ the objective function, and $m$ the quadratic model, the trust region
1033: radius is modified as follows
1035: $
1036: \delta =
1037: \begin{cases}
1038: \delta * t_1 ,& \rho < \eta_2 \\
1039: \delta * t_2 ,& \rho > \eta_3 \\
1040: \end{cases}
1041: $
1043: The step is accepted if $\rho > \eta_1$.
1044: Use `PETSC_DETERMINE` to use the default value for the given `SNES`.
1045: Use `PETSC_CURRENT` to retain a value.
1047: Fortran Note:
1048: Use `PETSC_DETERMINE_REAL`, `PETSC_CURRENT_REAL`
1050: Level: intermediate
1052: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTR`, `SNESSetObjective()`, `SNESNewtonTRGetUpdateParameters()`
1053: @*/
1054: PetscErrorCode SNESNewtonTRSetUpdateParameters(SNES snes, PetscReal eta1, PetscReal eta2, PetscReal eta3, PetscReal t1, PetscReal t2)
1055: {
1056: PetscBool flg;
1058: PetscFunctionBegin;
1065: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
1066: if (flg) {
1067: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
1069: if (eta1 == PETSC_DETERMINE) eta1 = tr->default_eta1;
1070: if (eta2 == PETSC_DETERMINE) eta2 = tr->default_eta2;
1071: if (eta3 == PETSC_DETERMINE) eta3 = tr->default_eta3;
1072: if (t1 == PETSC_DETERMINE) t1 = tr->default_t1;
1073: if (t2 == PETSC_DETERMINE) t2 = tr->default_t2;
1074: if (eta1 != PETSC_CURRENT) tr->eta1 = eta1;
1075: if (eta2 != PETSC_CURRENT) tr->eta2 = eta2;
1076: if (eta3 != PETSC_CURRENT) tr->eta3 = eta3;
1077: if (t1 != PETSC_CURRENT) tr->t1 = t1;
1078: if (t2 != PETSC_CURRENT) tr->t2 = t2;
1079: }
1080: PetscFunctionReturn(PETSC_SUCCESS);
1081: }
1083: /*@
1084: SNESNewtonTRGetUpdateParameters - Gets the trust region update parameters.
1086: Not Collective
1088: Input Parameter:
1089: . snes - the `SNES` context
1091: Output Parameters:
1092: + eta1 - acceptance tolerance
1093: . eta2 - shrinking tolerance
1094: . eta3 - enlarging tolerance
1095: . t1 - shrink factor
1096: - t2 - enlarge factor
1098: Level: intermediate
1100: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTR`, `SNESNewtonTRSetUpdateParameters()`
1101: @*/
1102: PetscErrorCode SNESNewtonTRGetUpdateParameters(SNES snes, PetscReal *eta1, PetscReal *eta2, PetscReal *eta3, PetscReal *t1, PetscReal *t2)
1103: {
1104: SNES_NEWTONTR *tr;
1105: PetscBool flg;
1107: PetscFunctionBegin;
1109: if (eta1) PetscAssertPointer(eta1, 2);
1110: if (eta2) PetscAssertPointer(eta2, 3);
1111: if (eta3) PetscAssertPointer(eta3, 4);
1112: if (t1) PetscAssertPointer(t1, 5);
1113: if (t2) PetscAssertPointer(t2, 6);
1114: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
1115: PetscAssert(flg, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)snes)->type_name);
1116: tr = (SNES_NEWTONTR *)snes->data;
1117: if (eta1) *eta1 = tr->eta1;
1118: if (eta2) *eta2 = tr->eta2;
1119: if (eta3) *eta3 = tr->eta3;
1120: if (t1) *t1 = tr->t1;
1121: if (t2) *t2 = tr->t2;
1122: PetscFunctionReturn(PETSC_SUCCESS);
1123: }
1125: /*MC
1126: SNESNEWTONTR - Newton based nonlinear solver that uses a trust-region strategy
1128: Options Database Keys:
1129: + -snes_tr_deltamin deltamin - trust region parameter, minimum size of trust region
1130: . -snes_tr_deltamax deltamax - trust region parameter, max size of trust region (default: 1e10)
1131: . -snes_tr_delta0 delta0 - trust region parameter, initial size of trust region (default: 0.2)
1132: . -snes_tr_eta1 eta1 - trust region parameter $eta1 \le eta2$, $rho > eta1$ breaks out of the inner iteration (default: 0.001)
1133: . -snes_tr_eta2 eta2 - trust region parameter, $rho \le eta2$ shrinks the trust region (default: 0.25)
1134: . -snes_tr_eta3 eta3 - trust region parameter $eta3 > eta2$, $rho \ge eta3$ expands the trust region (default: 0.75)
1135: . -snes_tr_t1 t1 - trust region parameter, shrinking factor of trust region (default: 0.25)
1136: . -snes_tr_t2 t2 - trust region parameter, expanding factor of trust region (default: 2.0)
1137: . -snes_tr_norm_type (1|2|infinity) - Type of norm for trust region bounds (default: 2)
1138: - -snes_tr_fallback_type (newton,cauchy,dogleg) - Solution strategy to test reduction when step is outside of trust region.
1139: Can use scaled Newton direction, Cauchy point (Steepest Descent direction) or dogleg method.
1141: Level: beginner
1143: Notes:
1144: The code is largely based on the book {cite}`nocedal2006numerical` and supports minimizing objective functions using a quadratic model.
1145: Quasi-Newton models are also supported.
1147: Default step computation uses the Newton direction, but a dogleg type update is also supported.
1148: The 1- and infinity-norms are also supported, via `SNESNewtonTRSetNormType()`, when computing the trust region bounds.
1150: .seealso: [](ch_snes), `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESSetObjective()`,
1151: `SNESNewtonTRSetTolerances()`, `SNESNewtonTRSetUpdateParameters()`,
1152: `SNESNewtonTRSetNormType()`, `SNESNewtonTRSetFallbackType()`, `SNESNewtonTRSetQNType()`,
1153: `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRSetPreCheck()`
1154: M*/
1155: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONTR(SNES snes)
1156: {
1157: SNES_NEWTONTR *neP;
1159: PetscFunctionBegin;
1160: snes->ops->setup = SNESSetUp_NEWTONTR;
1161: snes->ops->solve = SNESSolve_NEWTONTR;
1162: snes->ops->reset = SNESReset_NEWTONTR;
1163: snes->ops->destroy = SNESDestroy_NEWTONTR;
1164: snes->ops->setfromoptions = SNESSetFromOptions_NEWTONTR;
1165: snes->ops->view = SNESView_NEWTONTR;
1167: PetscObjectParameterSetDefault(snes, stol, 0.0);
1169: snes->usesksp = PETSC_TRUE;
1170: snes->npcside = PC_RIGHT;
1171: snes->usesnpc = PETSC_TRUE;
1173: snes->alwayscomputesfinalresidual = PETSC_TRUE;
1175: PetscCall(PetscNew(&neP));
1176: snes->data = (void *)neP;
1178: PetscObjectParameterSetDefault(neP, eta1, 0.001);
1179: PetscObjectParameterSetDefault(neP, eta2, 0.25);
1180: PetscObjectParameterSetDefault(neP, eta3, 0.75);
1181: PetscObjectParameterSetDefault(neP, t1, 0.25);
1182: PetscObjectParameterSetDefault(neP, t2, 2.0);
1183: PetscObjectParameterSetDefault(neP, deltam, PetscDefined(USE_REAL_SINGLE) ? 1.e-6 : 1.e-12);
1184: PetscObjectParameterSetDefault(neP, delta0, 0.2);
1185: PetscObjectParameterSetDefault(neP, deltaM, 1.e10);
1187: neP->norm = NORM_2;
1188: neP->fallback = SNES_TR_FALLBACK_NEWTON;
1189: neP->kmdc = 0.0; /* by default do not use sufficient decrease */
1191: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonTRSetTolerances_C", SNESNewtonTRSetTolerances_TR));
1192: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonTRGetTolerances_C", SNESNewtonTRGetTolerances_TR));
1193: PetscFunctionReturn(PETSC_SUCCESS);
1194: }