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 SNESComputeJacobian_MATLMVM(SNES snes, Vec X, Mat J, Mat B, void *dummy)
14: {
15: PetscFunctionBegin;
16: // PetscCall(MatLMVMSymBroydenSetDelta(B, _some_delta));
17: PetscCall(MatLMVMUpdate(B, X, snes->vec_func));
18: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
19: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
20: if (J != B) {
21: // PetscCall(MatLMVMSymBroydenSetDelta(J, _some_delta));
22: PetscCall(MatLMVMUpdate(J, X, snes->vec_func));
23: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
24: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
25: }
26: PetscFunctionReturn(PETSC_SUCCESS);
27: }
29: static PetscErrorCode SNESTR_KSPConverged_Private(KSP ksp, PetscInt n, PetscReal rnorm, KSPConvergedReason *reason, void *cctx)
30: {
31: SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx *)cctx;
32: SNES snes = ctx->snes;
33: SNES_NEWTONTR *neP = (SNES_NEWTONTR *)snes->data;
34: Vec x;
35: PetscReal nrm;
37: PetscFunctionBegin;
38: /* Determine norm of solution */
39: PetscCall(KSPBuildSolution(ksp, NULL, &x));
40: PetscCall(VecNorm(x, neP->norm, &nrm));
41: if (nrm >= neP->delta) {
42: PetscCall(PetscInfo(snes, "Ending linear iteration early due to exiting trust region, delta=%g, length=%g\n", (double)neP->delta, (double)nrm));
43: *reason = KSP_CONVERGED_STEP_LENGTH;
44: PetscFunctionReturn(PETSC_SUCCESS);
45: }
46: PetscCall((*ctx->convtest)(ksp, n, rnorm, reason, ctx->convctx));
47: if (*reason) PetscCall(PetscInfo(snes, "Default or user provided convergence test KSP iterations=%" PetscInt_FMT ", rnorm=%g\n", n, (double)rnorm));
48: PetscFunctionReturn(PETSC_SUCCESS);
49: }
51: static PetscErrorCode SNESTR_KSPConverged_Destroy(void *cctx)
52: {
53: SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx *)cctx;
55: PetscFunctionBegin;
56: PetscCall((*ctx->convdestroy)(ctx->convctx));
57: PetscCall(PetscFree(ctx));
58: PetscFunctionReturn(PETSC_SUCCESS);
59: }
61: static PetscErrorCode SNESTR_Converged_Private(SNES snes, PetscInt it, PetscReal xnorm, PetscReal pnorm, PetscReal fnorm, SNESConvergedReason *reason, void *dummy)
62: {
63: SNES_NEWTONTR *neP = (SNES_NEWTONTR *)snes->data;
65: PetscFunctionBegin;
66: *reason = SNES_CONVERGED_ITERATING;
67: if (neP->delta < snes->deltatol) {
68: PetscCall(PetscInfo(snes, "Diverged due to too small a trust region %g<%g\n", (double)neP->delta, (double)snes->deltatol));
69: *reason = SNES_DIVERGED_TR_DELTA;
70: } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
71: PetscCall(PetscInfo(snes, "Exceeded maximum number of function evaluations: %" PetscInt_FMT "\n", snes->max_funcs));
72: *reason = SNES_DIVERGED_FUNCTION_COUNT;
73: }
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
77: /*@
78: SNESNewtonTRSetNormType - Specify the type of norm to use for the computation of the trust region.
80: Input Parameters:
81: + snes - the nonlinear solver object
82: - norm - the norm type
84: Level: intermediate
86: .seealso: `SNESNEWTONTR`, `NormType`
87: @*/
88: PetscErrorCode SNESNewtonTRSetNormType(SNES snes, NormType norm)
89: {
90: PetscBool flg;
92: PetscFunctionBegin;
95: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
96: if (flg) {
97: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
99: tr->norm = norm;
100: }
101: PetscFunctionReturn(PETSC_SUCCESS);
102: }
104: /*@
105: SNESNewtonTRSetQNType - Specify to use a quasi-Newton model.
107: Input Parameters:
108: + snes - the nonlinear solver object
109: - use - the type of approximations to be used
111: Level: intermediate
113: Notes:
114: Options for the approximations can be set with the snes_tr_qn_ and snes_tr_qn_pre_ prefixes.
116: .seealso: `SNESNEWTONTR`, `SNESNewtonTRQNType`, `MATLMVM`
117: @*/
118: PetscErrorCode SNESNewtonTRSetQNType(SNES snes, SNESNewtonTRQNType use)
119: {
120: PetscBool flg;
122: PetscFunctionBegin;
125: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
126: if (flg) {
127: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
129: tr->qn = use;
130: }
131: PetscFunctionReturn(PETSC_SUCCESS);
132: }
134: /*@
135: SNESNewtonTRSetFallbackType - Set the type of fallback to use if the solution of the trust region subproblem is outside the radius
137: Input Parameters:
138: + snes - the nonlinear solver object
139: - ftype - the fallback type, see `SNESNewtonTRFallbackType`
141: Level: intermediate
143: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRPreCheck()`, `SNESNewtonTRGetPreCheck()`, `SNESNewtonTRSetPreCheck()`,
144: `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()`
145: @*/
146: PetscErrorCode SNESNewtonTRSetFallbackType(SNES snes, SNESNewtonTRFallbackType ftype)
147: {
148: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
149: PetscBool flg;
151: PetscFunctionBegin;
154: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
155: if (flg) tr->fallback = ftype;
156: PetscFunctionReturn(PETSC_SUCCESS);
157: }
159: /*@C
160: SNESNewtonTRSetPreCheck - Sets a user function that is called before the search step has been determined.
161: Allows the user a chance to change or override the trust region decision.
163: Logically Collective
165: Input Parameters:
166: + snes - the nonlinear solver object
167: . func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRPreCheck()`
168: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
170: Level: intermediate
172: Note:
173: This function is called BEFORE the function evaluation within the solver.
175: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRPreCheck()`, `SNESNewtonTRGetPreCheck()`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()`,
176: @*/
177: PetscErrorCode SNESNewtonTRSetPreCheck(SNES snes, PetscErrorCode (*func)(SNES, Vec, Vec, PetscBool *, void *), void *ctx)
178: {
179: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
180: PetscBool flg;
182: PetscFunctionBegin;
184: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
185: if (flg) {
186: if (func) tr->precheck = func;
187: if (ctx) tr->precheckctx = ctx;
188: }
189: PetscFunctionReturn(PETSC_SUCCESS);
190: }
192: /*@C
193: SNESNewtonTRGetPreCheck - Gets the pre-check function
195: Not Collective
197: Input Parameter:
198: . snes - the nonlinear solver context
200: Output Parameters:
201: + func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRPreCheck()`
202: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
204: Level: intermediate
206: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRPreCheck()`
207: @*/
208: PetscErrorCode SNESNewtonTRGetPreCheck(SNES snes, PetscErrorCode (**func)(SNES, Vec, Vec, PetscBool *, void *), void **ctx)
209: {
210: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
211: PetscBool flg;
213: PetscFunctionBegin;
215: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
216: PetscAssert(flg, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)snes)->type_name);
217: if (func) *func = tr->precheck;
218: if (ctx) *ctx = tr->precheckctx;
219: PetscFunctionReturn(PETSC_SUCCESS);
220: }
222: /*@C
223: SNESNewtonTRSetPostCheck - Sets a user function that is called after the search step has been determined but before the next
224: function evaluation. Allows the user a chance to change or override the internal decision of the solver
226: Logically Collective
228: Input Parameters:
229: + snes - the nonlinear solver object
230: . func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRPostCheck()`
231: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
233: Level: intermediate
235: Note:
236: This function is called BEFORE the function evaluation within the solver while the function set in
237: `SNESLineSearchSetPostCheck()` is called AFTER the function evaluation.
239: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRPostCheck()`, `SNESNewtonTRGetPostCheck()`, `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRGetPreCheck()`
240: @*/
241: PetscErrorCode SNESNewtonTRSetPostCheck(SNES snes, PetscErrorCode (*func)(SNES, Vec, Vec, Vec, PetscBool *, PetscBool *, void *), void *ctx)
242: {
243: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
244: PetscBool flg;
246: PetscFunctionBegin;
248: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
249: if (flg) {
250: if (func) tr->postcheck = func;
251: if (ctx) tr->postcheckctx = ctx;
252: }
253: PetscFunctionReturn(PETSC_SUCCESS);
254: }
256: /*@C
257: SNESNewtonTRGetPostCheck - Gets the post-check function
259: Not Collective
261: Input Parameter:
262: . snes - the nonlinear solver context
264: Output Parameters:
265: + func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRPostCheck()`
266: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
268: Level: intermediate
270: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRPostCheck()`
271: @*/
272: PetscErrorCode SNESNewtonTRGetPostCheck(SNES snes, PetscErrorCode (**func)(SNES, Vec, Vec, Vec, PetscBool *, PetscBool *, void *), void **ctx)
273: {
274: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
275: PetscBool flg;
277: PetscFunctionBegin;
279: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
280: PetscAssert(flg, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)snes)->type_name);
281: if (func) *func = tr->postcheck;
282: if (ctx) *ctx = tr->postcheckctx;
283: PetscFunctionReturn(PETSC_SUCCESS);
284: }
286: /*@C
287: SNESNewtonTRPreCheck - Runs the precheck routine
289: Logically Collective
291: Input Parameters:
292: + snes - the solver
293: . X - The last solution
294: - Y - The step direction
296: Output Parameter:
297: . changed_Y - Indicator that the step direction `Y` has been changed.
299: Level: intermediate
301: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRGetPreCheck()`, `SNESNewtonTRPostCheck()`
302: @*/
303: PetscErrorCode SNESNewtonTRPreCheck(SNES snes, Vec X, Vec Y, PetscBool *changed_Y)
304: {
305: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
306: PetscBool flg;
308: PetscFunctionBegin;
310: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
311: PetscAssert(flg, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)snes)->type_name);
312: *changed_Y = PETSC_FALSE;
313: if (tr->precheck) {
314: PetscCall((*tr->precheck)(snes, X, Y, changed_Y, tr->precheckctx));
316: }
317: PetscFunctionReturn(PETSC_SUCCESS);
318: }
320: /*@C
321: SNESNewtonTRPostCheck - Runs the postcheck routine
323: Logically Collective
325: Input Parameters:
326: + snes - the solver
327: . X - The last solution
328: . Y - The full step direction
329: - W - The updated solution, W = X - Y
331: Output Parameters:
332: + changed_Y - indicator if step has been changed
333: - changed_W - Indicator if the new candidate solution W has been changed.
335: Note:
336: If Y is changed then W is recomputed as X - Y
338: Level: intermediate
340: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()`, `SNESNewtonTRPreCheck()`
341: @*/
342: PetscErrorCode SNESNewtonTRPostCheck(SNES snes, Vec X, Vec Y, Vec W, PetscBool *changed_Y, PetscBool *changed_W)
343: {
344: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
345: PetscBool flg;
347: PetscFunctionBegin;
349: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
350: PetscAssert(flg, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)snes)->type_name);
351: *changed_Y = PETSC_FALSE;
352: *changed_W = PETSC_FALSE;
353: if (tr->postcheck) {
354: PetscCall((*tr->postcheck)(snes, X, Y, W, changed_Y, changed_W, tr->postcheckctx));
357: }
358: PetscFunctionReturn(PETSC_SUCCESS);
359: }
361: /* stable implementation of roots of a*x^2 + b*x + c = 0 */
362: static inline void PetscQuadraticRoots(PetscReal a, PetscReal b, PetscReal c, PetscReal *xm, PetscReal *xp)
363: {
364: PetscReal temp = -0.5 * (b + PetscCopysignReal(1.0, b) * PetscSqrtReal(b * b - 4 * a * c));
365: PetscReal x1 = temp / a;
366: PetscReal x2 = c / temp;
367: *xm = PetscMin(x1, x2);
368: *xp = PetscMax(x1, x2);
369: }
371: /* Computes the quadratic model difference */
372: static PetscErrorCode SNESNewtonTRQuadraticDelta(SNES snes, Mat J, PetscBool has_objective, Vec Y, Vec GradF, Vec W, PetscReal *yTHy_, PetscReal *gTy_, PetscReal *deltaqm)
373: {
374: PetscReal yTHy, gTy;
376: PetscFunctionBegin;
377: PetscCall(MatMult(J, Y, W));
378: if (has_objective) PetscCall(VecDotRealPart(Y, W, &yTHy));
379: else PetscCall(VecDotRealPart(W, W, &yTHy)); /* Gauss-Newton approximation J^t * J */
380: PetscCall(VecDotRealPart(GradF, Y, &gTy));
381: *deltaqm = -(-(gTy) + 0.5 * (yTHy)); /* difference in quadratic model, -gTy because SNES solves it this way */
382: if (yTHy_) *yTHy_ = yTHy;
383: if (gTy_) *gTy_ = gTy;
384: PetscFunctionReturn(PETSC_SUCCESS);
385: }
387: /* Computes the new objective given X = Xk, Y = direction
388: W work vector, on output W = X - Y
389: G work vector, on output G = SNESFunction(W) */
390: static PetscErrorCode SNESNewtonTRObjective(SNES snes, PetscBool has_objective, Vec X, Vec Y, Vec W, Vec G, PetscReal *gnorm, PetscReal *fkp1)
391: {
392: PetscBool changed_y, changed_w;
394: PetscFunctionBegin;
395: /* TODO: we can add a linesearch here */
396: PetscCall(SNESNewtonTRPreCheck(snes, X, Y, &changed_y));
397: PetscCall(VecWAXPY(W, -1.0, Y, X)); /* Xkp1 */
398: PetscCall(SNESNewtonTRPostCheck(snes, X, Y, W, &changed_y, &changed_w));
399: if (changed_y && !changed_w) PetscCall(VecWAXPY(W, -1.0, Y, X));
401: PetscCall(SNESComputeFunction(snes, W, G)); /* F(Xkp1) = G */
402: PetscCall(VecNorm(G, NORM_2, gnorm));
403: if (has_objective) PetscCall(SNESComputeObjective(snes, W, fkp1));
404: else *fkp1 = 0.5 * PetscSqr(*gnorm);
405: PetscFunctionReturn(PETSC_SUCCESS);
406: }
408: static PetscErrorCode SNESSetUpQN_NEWTONTR(SNES snes)
409: {
410: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
412: PetscFunctionBegin;
413: PetscCall(MatDestroy(&tr->qnB));
414: PetscCall(MatDestroy(&tr->qnB_pre));
415: if (tr->qn) {
416: PetscInt n, N;
417: const char *optionsprefix;
418: Mat B;
420: PetscCall(MatCreate(PetscObjectComm((PetscObject)snes), &B));
421: PetscCall(SNESGetOptionsPrefix(snes, &optionsprefix));
422: PetscCall(MatSetOptionsPrefix(B, "snes_tr_qn_"));
423: PetscCall(MatAppendOptionsPrefix(B, optionsprefix));
424: PetscCall(MatSetType(B, MATLMVMBFGS));
425: PetscCall(VecGetLocalSize(snes->vec_sol, &n));
426: PetscCall(VecGetSize(snes->vec_sol, &N));
427: PetscCall(MatSetSizes(B, n, n, N, N));
428: PetscCall(MatSetUp(B));
429: PetscCall(MatSetFromOptions(B));
430: PetscCall(MatLMVMAllocate(B, snes->vec_sol, snes->vec_func));
431: tr->qnB = B;
432: if (tr->qn == SNES_TR_QN_DIFFERENT) {
433: PetscCall(MatCreate(PetscObjectComm((PetscObject)snes), &B));
434: PetscCall(SNESGetOptionsPrefix(snes, &optionsprefix));
435: PetscCall(MatSetOptionsPrefix(B, "snes_tr_qn_pre_"));
436: PetscCall(MatAppendOptionsPrefix(B, optionsprefix));
437: PetscCall(MatSetType(B, MATLMVMBFGS));
438: PetscCall(MatSetSizes(B, n, n, N, N));
439: PetscCall(MatSetUp(B));
440: PetscCall(MatSetFromOptions(B));
441: PetscCall(MatLMVMAllocate(B, snes->vec_sol, snes->vec_func));
442: tr->qnB_pre = B;
443: } else {
444: PetscCall(PetscObjectReference((PetscObject)tr->qnB));
445: tr->qnB_pre = tr->qnB;
446: }
447: }
448: PetscFunctionReturn(PETSC_SUCCESS);
449: }
451: /*
452: SNESSolve_NEWTONTR - Implements Newton's Method with trust-region subproblem and adds dogleg Cauchy
453: (Steepest Descent direction) step and direction if the trust region is not satisfied for solving system of
454: nonlinear equations
456: */
457: static PetscErrorCode SNESSolve_NEWTONTR(SNES snes)
458: {
459: SNES_NEWTONTR *neP = (SNES_NEWTONTR *)snes->data;
460: Vec X, F, Y, G, W, GradF, YU, Yc;
461: PetscInt maxits, lits;
462: PetscReal rho, fnorm, gnorm = 0.0, xnorm = 0.0, delta, ynorm;
463: PetscReal deltaM, fk, fkp1, deltaqm = 0.0, gTy = 0.0, yTHy = 0.0;
464: PetscReal auk, tauk, gfnorm, gfnorm_k, ycnorm, gTBg, objmin = 0.0, beta_k = 1.0;
465: PC pc;
466: Mat J, Jp;
467: PetscBool already_done = PETSC_FALSE, on_boundary;
468: PetscBool clear_converged_test, rho_satisfied, has_objective;
469: SNES_TR_KSPConverged_Ctx *ctx;
470: void *convctx;
472: PetscErrorCode (*convtest)(KSP, PetscInt, PetscReal, KSPConvergedReason *, void *), (*convdestroy)(void *);
473: PetscErrorCode (*objective)(SNES, Vec, PetscReal *, void *);
475: PetscFunctionBegin;
476: PetscCall(SNESGetObjective(snes, &objective, NULL));
477: has_objective = objective ? PETSC_TRUE : PETSC_FALSE;
479: maxits = snes->max_its; /* maximum number of iterations */
480: X = snes->vec_sol; /* solution vector */
481: F = snes->vec_func; /* residual vector */
482: Y = snes->vec_sol_update; /* update vector */
483: G = snes->work[0]; /* updated residual */
484: W = snes->work[1]; /* temporary vector */
485: GradF = !has_objective ? snes->work[2] : snes->vec_func; /* grad f = J^T F */
486: YU = snes->work[3]; /* work vector for dogleg method */
487: Yc = snes->work[4]; /* Cauchy point */
489: 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);
491: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
492: snes->iter = 0;
493: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
495: /* setup QN matrices if needed */
496: PetscCall(SNESSetUpQN_NEWTONTR(snes));
498: /* Set the linear stopping criteria to use the More' trick if needed */
499: clear_converged_test = PETSC_FALSE;
500: PetscCall(SNESGetKSP(snes, &snes->ksp));
501: PetscCall(KSPGetConvergenceTest(snes->ksp, &convtest, &convctx, &convdestroy));
502: if (convtest != SNESTR_KSPConverged_Private) {
503: clear_converged_test = PETSC_TRUE;
504: PetscCall(PetscNew(&ctx));
505: ctx->snes = snes;
506: PetscCall(KSPGetAndClearConvergenceTest(snes->ksp, &ctx->convtest, &ctx->convctx, &ctx->convdestroy));
507: PetscCall(KSPSetConvergenceTest(snes->ksp, SNESTR_KSPConverged_Private, ctx, SNESTR_KSPConverged_Destroy));
508: PetscCall(PetscInfo(snes, "Using Krylov convergence test SNESTR_KSPConverged_Private\n"));
509: }
511: if (!snes->vec_func_init_set) {
512: PetscCall(SNESComputeFunction(snes, X, F)); /* F(X) */
513: } else snes->vec_func_init_set = PETSC_FALSE;
515: PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- || F || */
516: SNESCheckFunctionNorm(snes, fnorm);
517: PetscCall(VecNorm(X, NORM_2, &xnorm)); /* xnorm <- || X || */
519: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
520: snes->norm = fnorm;
521: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
522: delta = neP->delta0;
523: deltaM = neP->deltaM;
524: neP->delta = delta;
525: PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));
527: /* test convergence */
528: rho_satisfied = PETSC_FALSE;
529: PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm));
530: PetscCall(SNESMonitor(snes, 0, fnorm));
531: if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
533: if (has_objective) PetscCall(SNESComputeObjective(snes, X, &fk));
534: else fk = 0.5 * PetscSqr(fnorm); /* obj(x) = 0.5 * ||F(x)||^2 */
536: /* hook state vector to BFGS preconditioner */
537: PetscCall(KSPGetPC(snes->ksp, &pc));
538: PetscCall(PCLMVMSetUpdateVec(pc, X));
540: if (neP->kmdc) PetscCall(KSPSetComputeEigenvalues(snes->ksp, PETSC_TRUE));
542: while (snes->iter < maxits) {
543: /* calculating Jacobian and GradF of minimization function only once */
544: if (!already_done) {
545: /* Call general purpose update function */
546: PetscTryTypeMethod(snes, update, snes->iter);
548: /* apply the nonlinear preconditioner */
549: if (snes->npc && snes->npcside == PC_RIGHT) {
550: SNESConvergedReason reason;
552: PetscCall(SNESSetInitialFunction(snes->npc, F));
553: PetscCall(PetscLogEventBegin(SNES_NPCSolve, snes->npc, X, snes->vec_rhs, 0));
554: PetscCall(SNESSolve(snes->npc, snes->vec_rhs, X));
555: PetscCall(PetscLogEventEnd(SNES_NPCSolve, snes->npc, X, snes->vec_rhs, 0));
556: PetscCall(SNESGetConvergedReason(snes->npc, &reason));
557: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
558: snes->reason = SNES_DIVERGED_INNER;
559: PetscFunctionReturn(PETSC_SUCCESS);
560: }
561: // XXX
562: PetscCall(SNESGetNPCFunction(snes, F, &fnorm));
563: } else if (snes->ops->update) { /* if update is present, recompute objective function and function norm */
564: PetscCall(SNESComputeFunction(snes, X, F));
565: }
567: /* Jacobian */
568: J = NULL;
569: Jp = NULL;
570: if (!neP->qnB) {
571: PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
572: J = snes->jacobian;
573: Jp = snes->jacobian_pre;
574: } else { /* QN model */
575: PetscCall(SNESComputeJacobian_MATLMVM(snes, X, neP->qnB, neP->qnB_pre, NULL));
576: J = neP->qnB;
577: Jp = neP->qnB_pre;
578: }
579: SNESCheckJacobianDomainerror(snes);
581: /* objective function */
582: PetscCall(VecNorm(F, NORM_2, &fnorm));
583: if (has_objective) PetscCall(SNESComputeObjective(snes, X, &fk));
584: else fk = 0.5 * PetscSqr(fnorm); /* obj(x) = 0.5 * ||F(x)||^2 */
586: /* GradF */
587: if (has_objective) gfnorm = fnorm;
588: else {
589: PetscCall(MatMultTranspose(J, F, GradF)); /* grad f = J^T F */
590: PetscCall(VecNorm(GradF, NORM_2, &gfnorm));
591: }
592: PetscCall(VecNorm(GradF, neP->norm, &gfnorm_k));
593: }
594: already_done = PETSC_TRUE;
596: /* solve trust-region subproblem */
598: /* first compute Cauchy Point */
599: PetscCall(MatMult(J, GradF, W));
600: if (has_objective) PetscCall(VecDotRealPart(GradF, W, &gTBg));
601: else PetscCall(VecDotRealPart(W, W, &gTBg)); /* B = J^t * J */
602: /* Eqs 4.11 and 4.12 in Nocedal and Wright 2nd Edition (4.7 and 4.8 in 1st Edition) */
603: auk = delta / gfnorm_k;
604: if (gTBg < 0.0) tauk = 1.0;
605: else tauk = PetscMin(gfnorm * gfnorm * gfnorm_k / (delta * gTBg), 1);
606: auk *= tauk;
607: ycnorm = auk * gfnorm;
608: PetscCall(VecAXPBY(Yc, auk, 0.0, GradF));
610: on_boundary = PETSC_FALSE;
611: if (tauk != 1.0) {
612: KSPConvergedReason reason;
614: /* sufficient decrease (see 6.3.27 in Conn, Gould, Toint "Trust Region Methods")
615: beta_k the largest eigenvalue of the Hessian. Here we use the previous estimated value */
616: objmin = -neP->kmdc * gnorm * PetscMin(gnorm / beta_k, delta);
617: PetscCall(KSPCGSetObjectiveTarget(snes->ksp, objmin));
619: /* specify radius if looking for Newton step and trust region norm is the l2 norm */
620: PetscCall(KSPCGSetRadius(snes->ksp, neP->fallback == SNES_TR_FALLBACK_NEWTON && neP->norm == NORM_2 ? delta : 0.0));
621: PetscCall(KSPSetOperators(snes->ksp, J, Jp));
622: PetscCall(KSPSolve(snes->ksp, F, Y));
623: SNESCheckKSPSolve(snes);
624: PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
625: PetscCall(KSPGetConvergedReason(snes->ksp, &reason));
626: on_boundary = (PetscBool)(reason == KSP_CONVERGED_STEP_LENGTH);
627: PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits));
628: if (neP->kmdc) { /* update estimated Hessian largest eigenvalue */
629: PetscReal emax, emin;
630: PetscCall(KSPComputeExtremeSingularValues(snes->ksp, &emax, &emin));
631: if (emax > 0.0) beta_k = emax + 1;
632: }
633: } else { /* Cauchy point is on the boundary, accept it */
634: on_boundary = PETSC_TRUE;
635: PetscCall(VecCopy(Yc, Y));
636: PetscCall(PetscInfo(snes, "CP evaluated on boundary. delta: %g, ycnorm: %g, gTBg: %g\n", (double)delta, (double)ycnorm, (double)gTBg));
637: }
638: PetscCall(VecNorm(Y, neP->norm, &ynorm));
640: /* decide what to do when the update is outside of trust region */
641: if (ynorm > delta || ynorm == 0.0) {
642: SNESNewtonTRFallbackType fallback = ynorm > 0.0 ? neP->fallback : SNES_TR_FALLBACK_CAUCHY;
644: PetscCheck(neP->norm == NORM_2 || fallback != SNES_TR_FALLBACK_DOGLEG, PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "DOGLEG without l2 norm not implemented");
645: switch (fallback) {
646: case SNES_TR_FALLBACK_NEWTON:
647: auk = delta / ynorm;
648: PetscCall(VecScale(Y, auk));
649: PetscCall(PetscInfo(snes, "SN evaluated. delta: %g, ynorm: %g\n", (double)delta, (double)ynorm));
650: break;
651: case SNES_TR_FALLBACK_CAUCHY:
652: case SNES_TR_FALLBACK_DOGLEG:
653: if (fallback == SNES_TR_FALLBACK_CAUCHY || gTBg <= 0.0) {
654: PetscCall(VecCopy(Yc, Y));
655: PetscCall(PetscInfo(snes, "CP evaluated. delta: %g, ynorm: %g, ycnorm: %g, gTBg: %g\n", (double)delta, (double)ynorm, (double)ycnorm, (double)gTBg));
656: } else { /* take linear combination of Cauchy and Newton direction and step */
657: auk = gfnorm * gfnorm / gTBg;
658: if (gfnorm_k * auk >= delta) { /* first leg: Cauchy point outside of trust region */
659: PetscCall(VecAXPBY(Y, delta / gfnorm_k, 0.0, GradF));
660: PetscCall(PetscInfo(snes, "CP evaluated (outside region). delta: %g, ynorm: %g, ycnorm: %g\n", (double)delta, (double)ynorm, (double)ycnorm));
661: } else { /* second leg */
662: PetscReal c0, c1, c2, tau = 0.0, tpos, tneg;
663: PetscBool noroots;
665: /* Find solutions of (Eq. 4.16 in Nocedal and Wright)
666: ||p_U + lambda * (p_B - p_U)||^2 - delta^2 = 0,
667: where p_U the Cauchy direction, p_B the Newton direction */
668: PetscCall(VecAXPBY(YU, auk, 0.0, GradF));
669: PetscCall(VecAXPY(Y, -1.0, YU));
670: PetscCall(VecNorm(Y, NORM_2, &c0));
671: PetscCall(VecDotRealPart(YU, Y, &c1));
672: c0 = PetscSqr(c0);
673: c2 = PetscSqr(ycnorm) - PetscSqr(delta);
674: PetscQuadraticRoots(c0, 2 * c1, c2, &tneg, &tpos);
676: /* In principle the DL strategy as a unique solution in [0,1]
677: here we check that for some reason we numerically failed
678: to compute it. In that case, we use the Cauchy point */
679: noroots = PetscIsInfOrNanReal(tneg);
680: if (!noroots) {
681: if (tpos > 1) {
682: if (tneg >= 0 && tneg <= 1) {
683: tau = tneg;
684: } else noroots = PETSC_TRUE;
685: } else if (tpos >= 0) {
686: tau = tpos;
687: } else noroots = PETSC_TRUE;
688: }
689: if (noroots) { /* No roots, select Cauchy point */
690: PetscCall(VecCopy(Yc, Y));
691: } else {
692: PetscCall(VecAXPBY(Y, 1.0, tau, YU));
693: }
694: 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));
695: }
696: }
697: break;
698: default:
699: SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "Unknown fallback mode");
700: break;
701: }
702: }
704: /* compute the quadratic model difference */
705: PetscCall(SNESNewtonTRQuadraticDelta(snes, J, has_objective, Y, GradF, W, &yTHy, &gTy, &deltaqm));
707: /* Compute new objective function */
708: PetscCall(SNESNewtonTRObjective(snes, has_objective, X, Y, W, G, &gnorm, &fkp1));
709: if (PetscIsInfOrNanReal(fkp1)) rho = neP->eta1;
710: else {
711: if (deltaqm > 0.0) rho = (fk - fkp1) / deltaqm; /* actual improvement over predicted improvement */
712: else rho = neP->eta1; /* no reduction in quadratic model, step must be rejected */
713: }
715: PetscCall(VecNorm(Y, neP->norm, &ynorm));
716: 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));
718: /* update the size of the trust region */
719: if (rho < neP->eta2) delta *= neP->t1; /* shrink the region */
720: else if (rho > neP->eta3 && on_boundary) delta *= neP->t2; /* expand the region */
721: delta = PetscMin(delta, deltaM); /* but not greater than deltaM */
723: /* log 2-norm of update for moniroting routines */
724: PetscCall(VecNorm(Y, NORM_2, &ynorm));
726: /* decide on new step */
727: neP->delta = delta;
728: if (rho > neP->eta1) {
729: rho_satisfied = PETSC_TRUE;
730: } else {
731: rho_satisfied = PETSC_FALSE;
732: PetscCall(PetscInfo(snes, "Trying again in smaller region\n"));
733: /* check to see if progress is hopeless */
734: PetscCall(SNESTR_Converged_Private(snes, snes->iter, xnorm, ynorm, fnorm, &snes->reason, snes->cnvP));
735: if (!snes->reason) PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm));
736: if (snes->reason == SNES_CONVERGED_SNORM_RELATIVE) snes->reason = SNES_DIVERGED_TR_DELTA;
737: snes->numFailures++;
738: /* We're not progressing, so return with the current iterate */
739: if (snes->reason) break;
740: }
741: if (rho_satisfied) {
742: /* Update function values */
743: already_done = PETSC_FALSE;
744: fnorm = gnorm;
745: fk = fkp1;
747: /* New residual and linearization point */
748: PetscCall(VecCopy(G, F));
749: PetscCall(VecCopy(W, X));
751: /* Monitor convergence */
752: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
753: snes->iter++;
754: snes->norm = fnorm;
755: snes->xnorm = xnorm;
756: snes->ynorm = ynorm;
757: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
758: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits));
760: /* Test for convergence, xnorm = || X || */
761: PetscCall(VecNorm(X, NORM_2, &xnorm));
762: PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm));
763: PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
764: if (snes->reason) break;
765: }
766: }
768: if (clear_converged_test) {
769: PetscCall(KSPGetAndClearConvergenceTest(snes->ksp, &ctx->convtest, &ctx->convctx, &ctx->convdestroy));
770: PetscCall(PetscFree(ctx));
771: PetscCall(KSPSetConvergenceTest(snes->ksp, convtest, convctx, convdestroy));
772: }
773: PetscFunctionReturn(PETSC_SUCCESS);
774: }
776: static PetscErrorCode SNESSetUp_NEWTONTR(SNES snes)
777: {
778: PetscFunctionBegin;
779: PetscCall(SNESSetWorkVecs(snes, 5));
780: PetscCall(SNESSetUpMatrices(snes));
781: PetscFunctionReturn(PETSC_SUCCESS);
782: }
784: static PetscErrorCode SNESReset_NEWTONTR(SNES snes)
785: {
786: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
788: PetscFunctionBegin;
789: PetscCall(MatDestroy(&tr->qnB));
790: PetscCall(MatDestroy(&tr->qnB_pre));
791: PetscFunctionReturn(PETSC_SUCCESS);
792: }
794: static PetscErrorCode SNESDestroy_NEWTONTR(SNES snes)
795: {
796: PetscFunctionBegin;
797: PetscCall(SNESReset_NEWTONTR(snes));
798: PetscCall(PetscFree(snes->data));
799: PetscFunctionReturn(PETSC_SUCCESS);
800: }
802: static PetscErrorCode SNESSetFromOptions_NEWTONTR(SNES snes, PetscOptionItems *PetscOptionsObject)
803: {
804: SNES_NEWTONTR *ctx = (SNES_NEWTONTR *)snes->data;
805: SNESNewtonTRQNType qn;
806: SNESNewtonTRFallbackType fallback;
807: NormType norm;
808: PetscReal deltatol;
809: PetscBool flg;
811: PetscFunctionBegin;
812: PetscOptionsHeadBegin(PetscOptionsObject, "SNES trust region options for nonlinear equations");
813: PetscCall(PetscOptionsReal("-snes_tr_eta1", "eta1", "None", ctx->eta1, &ctx->eta1, NULL));
814: PetscCall(PetscOptionsReal("-snes_tr_eta2", "eta2", "None", ctx->eta2, &ctx->eta2, NULL));
815: PetscCall(PetscOptionsReal("-snes_tr_eta3", "eta3", "None", ctx->eta3, &ctx->eta3, NULL));
816: PetscCall(PetscOptionsReal("-snes_tr_t1", "t1", "None", ctx->t1, &ctx->t1, NULL));
817: PetscCall(PetscOptionsReal("-snes_tr_t2", "t2", "None", ctx->t2, &ctx->t2, NULL));
818: PetscCall(PetscOptionsReal("-snes_tr_deltaM", "deltaM", "None", ctx->deltaM, &ctx->deltaM, NULL));
819: PetscCall(PetscOptionsReal("-snes_tr_delta0", "delta0", "None", ctx->delta0, &ctx->delta0, NULL));
820: PetscCall(PetscOptionsReal("-snes_tr_kmdc", "sufficient decrease parameter", "None", ctx->kmdc, &ctx->kmdc, NULL));
822: deltatol = snes->deltatol;
823: PetscCall(PetscOptionsReal("-snes_tr_tol", "Trust region tolerance", "SNESSetTrustRegionTolerance", deltatol, &deltatol, &flg));
824: if (flg) PetscCall(SNESSetTrustRegionTolerance(snes, deltatol));
826: fallback = ctx->fallback;
827: 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));
828: if (flg) PetscCall(SNESNewtonTRSetFallbackType(snes, fallback));
830: qn = ctx->qn;
831: PetscCall(PetscOptionsEnum("-snes_tr_qn", "Use Quasi-Newton approximations for the model", "SNESNewtonTRSetQNType", SNESNewtonTRQNTypes, (PetscEnum)qn, (PetscEnum *)&qn, &flg));
832: if (flg) PetscCall(SNESNewtonTRSetQNType(snes, qn));
834: norm = ctx->norm;
835: PetscCall(PetscOptionsEnum("-snes_tr_norm_type", "Type of norm for trust region bounds", "SNESNewtonTRSetNormType", NormTypes, (PetscEnum)norm, (PetscEnum *)&norm, &flg));
836: if (flg) PetscCall(SNESNewtonTRSetNormType(snes, norm));
838: PetscOptionsHeadEnd();
839: PetscFunctionReturn(PETSC_SUCCESS);
840: }
842: static PetscErrorCode SNESView_NEWTONTR(SNES snes, PetscViewer viewer)
843: {
844: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
845: PetscBool iascii;
847: PetscFunctionBegin;
848: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
849: if (iascii) {
850: PetscCall(PetscViewerASCIIPrintf(viewer, " Trust region tolerance %g\n", (double)snes->deltatol));
851: PetscCall(PetscViewerASCIIPrintf(viewer, " eta1=%g, eta2=%g, eta3=%g\n", (double)tr->eta1, (double)tr->eta2, (double)tr->eta3));
852: PetscCall(PetscViewerASCIIPrintf(viewer, " delta0=%g, t1=%g, t2=%g, deltaM=%g\n", (double)tr->delta0, (double)tr->t1, (double)tr->t2, (double)tr->deltaM));
853: PetscCall(PetscViewerASCIIPrintf(viewer, " kmdc=%g\n", (double)tr->kmdc));
854: PetscCall(PetscViewerASCIIPrintf(viewer, " fallback=%s\n", SNESNewtonTRFallbackTypes[tr->fallback]));
855: if (tr->qn) PetscCall(PetscViewerASCIIPrintf(viewer, " qn=%s\n", SNESNewtonTRQNTypes[tr->qn]));
856: if (tr->norm != NORM_2) PetscCall(PetscViewerASCIIPrintf(viewer, " norm=%s\n", NormTypes[tr->norm]));
857: }
858: PetscFunctionReturn(PETSC_SUCCESS);
859: }
861: /*MC
862: SNESNEWTONTR - Newton based nonlinear solver that uses trust-region dogleg method with Cauchy direction {cite}`nocedal2006numerical`
864: Options Database Keys:
865: + -snes_tr_tol <tol> - trust region tolerance
866: . -snes_tr_eta1 <eta1> - trust region parameter eta1 <= eta2, rho > eta1 breaks out of the inner iteration (default: eta1=0.001)
867: . -snes_tr_eta2 <eta2> - trust region parameter, rho <= eta2 shrinks the trust region (default: eta2=0.25)
868: . -snes_tr_eta3 <eta3> - trust region parameter eta3 > eta2, rho >= eta3 expands the trust region (default: eta3=0.75)
869: . -snes_tr_t1 <t1> - trust region parameter, shrinking factor of trust region (default: 0.25)
870: . -snes_tr_t2 <t2> - trust region parameter, expanding factor of trust region (default: 2.0)
871: . -snes_tr_deltaM <deltaM> - trust region parameter, max size of trust region (default: MAX_REAL)
872: . -snes_tr_delta0 <delta0> - trust region parameter, initial size of trust region (default: 0.2)
873: - -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.
875: Level: beginner
877: .seealso: [](ch_snes), `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESSetTrustRegionTolerance()`,
878: `SNESNewtonTRPreCheck()`, `SNESNewtonTRGetPreCheck()`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()`,
879: `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRSetFallbackType()`, `SNESNewtonTRSetQNType()`
880: M*/
881: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONTR(SNES snes)
882: {
883: SNES_NEWTONTR *neP;
885: PetscFunctionBegin;
886: snes->ops->setup = SNESSetUp_NEWTONTR;
887: snes->ops->solve = SNESSolve_NEWTONTR;
888: snes->ops->reset = SNESReset_NEWTONTR;
889: snes->ops->destroy = SNESDestroy_NEWTONTR;
890: snes->ops->setfromoptions = SNESSetFromOptions_NEWTONTR;
891: snes->ops->view = SNESView_NEWTONTR;
893: snes->stol = 0.0;
894: snes->usesksp = PETSC_TRUE;
895: snes->npcside = PC_RIGHT;
896: snes->usesnpc = PETSC_TRUE;
898: snes->alwayscomputesfinalresidual = PETSC_TRUE;
900: PetscCall(PetscNew(&neP));
901: snes->data = (void *)neP;
902: neP->delta = 0.0;
903: neP->delta0 = 0.2;
904: neP->eta1 = 0.001;
905: neP->eta2 = 0.25;
906: neP->eta3 = 0.75;
907: neP->t1 = 0.25;
908: neP->t2 = 2.0;
909: neP->deltaM = 1.e10;
910: neP->norm = NORM_2;
911: neP->fallback = SNES_TR_FALLBACK_NEWTON;
912: neP->kmdc = 0.0; /* by default do not use sufficient decrease */
913: PetscFunctionReturn(PETSC_SUCCESS);
914: }