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: PetscCall((*ctx->convtest)(ksp, n, rnorm, reason, ctx->convctx));
39: if (*reason) PetscCall(PetscInfo(snes, "Default or user provided convergence test KSP iterations=%" PetscInt_FMT ", rnorm=%g\n", n, (double)rnorm));
40: /* Determine norm of solution */
41: PetscCall(KSPBuildSolution(ksp, NULL, &x));
42: PetscCall(VecNorm(x, neP->norm, &nrm));
43: if (nrm >= neP->delta) {
44: PetscCall(PetscInfo(snes, "Ending linear iteration early, delta=%g, length=%g\n", (double)neP->delta, (double)nrm));
45: *reason = KSP_CONVERGED_STEP_LENGTH;
46: }
47: PetscFunctionReturn(PETSC_SUCCESS);
48: }
50: static PetscErrorCode SNESTR_KSPConverged_Destroy(void *cctx)
51: {
52: SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx *)cctx;
54: PetscFunctionBegin;
55: PetscCall((*ctx->convdestroy)(ctx->convctx));
56: PetscCall(PetscFree(ctx));
57: PetscFunctionReturn(PETSC_SUCCESS);
58: }
60: static PetscErrorCode SNESTR_Converged_Private(SNES snes, PetscInt it, PetscReal xnorm, PetscReal pnorm, PetscReal fnorm, SNESConvergedReason *reason, void *dummy)
61: {
62: SNES_NEWTONTR *neP = (SNES_NEWTONTR *)snes->data;
64: PetscFunctionBegin;
65: *reason = SNES_CONVERGED_ITERATING;
66: if (neP->delta < snes->deltatol) {
67: PetscCall(PetscInfo(snes, "Diverged due to too small a trust region %g<%g\n", (double)neP->delta, (double)snes->deltatol));
68: *reason = SNES_DIVERGED_TR_DELTA;
69: } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
70: PetscCall(PetscInfo(snes, "Exceeded maximum number of function evaluations: %" PetscInt_FMT "\n", snes->max_funcs));
71: *reason = SNES_DIVERGED_FUNCTION_COUNT;
72: }
73: PetscFunctionReturn(PETSC_SUCCESS);
74: }
76: /*@
77: SNESNewtonTRSetNormType - Specify the type of norm to use for the computation of the trust region.
79: Input Parameters:
80: + snes - the nonlinear solver object
81: - norm - the norm type
83: Level: intermediate
85: .seealso: `SNESNEWTONTR`, `NormType`
86: @*/
87: PetscErrorCode SNESNewtonTRSetNormType(SNES snes, NormType norm)
88: {
89: PetscBool flg;
91: PetscFunctionBegin;
94: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
95: if (flg) {
96: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
98: tr->norm = norm;
99: }
100: PetscFunctionReturn(PETSC_SUCCESS);
101: }
103: /*@
104: SNESNewtonTRSetQNType - Specify to use a quasi-Newton model.
106: Input Parameters:
107: + snes - the nonlinear solver object
108: - use - the type of approximations to be used
110: Level: intermediate
112: Notes:
113: Options for the approximations can be set with the snes_tr_qn_ and snes_tr_qn_pre_ prefixes.
115: .seealso: `SNESNEWTONTR`, `SNESNewtonTRQNType`, `MATLMVM`
116: @*/
117: PetscErrorCode SNESNewtonTRSetQNType(SNES snes, SNESNewtonTRQNType use)
118: {
119: PetscBool flg;
121: PetscFunctionBegin;
124: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
125: if (flg) {
126: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
128: tr->qn = use;
129: }
130: PetscFunctionReturn(PETSC_SUCCESS);
131: }
133: /*@
134: SNESNewtonTRSetFallbackType - Set the type of fallback to use if the solution of the trust region subproblem is outside the radius
136: Input Parameters:
137: + snes - the nonlinear solver object
138: - ftype - the fallback type, see `SNESNewtonTRFallbackType`
140: Level: intermediate
142: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRPreCheck()`, `SNESNewtonTRGetPreCheck()`, `SNESNewtonTRSetPreCheck()`,
143: `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()`
144: @*/
145: PetscErrorCode SNESNewtonTRSetFallbackType(SNES snes, SNESNewtonTRFallbackType ftype)
146: {
147: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
148: PetscBool flg;
150: PetscFunctionBegin;
153: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
154: if (flg) tr->fallback = ftype;
155: PetscFunctionReturn(PETSC_SUCCESS);
156: }
158: /*@C
159: SNESNewtonTRSetPreCheck - Sets a user function that is called before the search step has been determined.
160: Allows the user a chance to change or override the trust region decision.
162: Logically Collective
164: Input Parameters:
165: + snes - the nonlinear solver object
166: . func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRPreCheck()`
167: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
169: Level: intermediate
171: Note:
172: This function is called BEFORE the function evaluation within the solver.
174: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRPreCheck()`, `SNESNewtonTRGetPreCheck()`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()`,
175: @*/
176: PetscErrorCode SNESNewtonTRSetPreCheck(SNES snes, PetscErrorCode (*func)(SNES, Vec, Vec, PetscBool *, void *), void *ctx)
177: {
178: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
179: PetscBool flg;
181: PetscFunctionBegin;
183: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
184: if (flg) {
185: if (func) tr->precheck = func;
186: if (ctx) tr->precheckctx = ctx;
187: }
188: PetscFunctionReturn(PETSC_SUCCESS);
189: }
191: /*@C
192: SNESNewtonTRGetPreCheck - Gets the pre-check function
194: Not Collective
196: Input Parameter:
197: . snes - the nonlinear solver context
199: Output Parameters:
200: + func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRPreCheck()`
201: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
203: Level: intermediate
205: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRPreCheck()`
206: @*/
207: PetscErrorCode SNESNewtonTRGetPreCheck(SNES snes, PetscErrorCode (**func)(SNES, Vec, Vec, PetscBool *, void *), void **ctx)
208: {
209: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
210: PetscBool flg;
212: PetscFunctionBegin;
214: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
215: PetscAssert(flg, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)snes)->type_name);
216: if (func) *func = tr->precheck;
217: if (ctx) *ctx = tr->precheckctx;
218: PetscFunctionReturn(PETSC_SUCCESS);
219: }
221: /*@C
222: SNESNewtonTRSetPostCheck - Sets a user function that is called after the search step has been determined but before the next
223: function evaluation. Allows the user a chance to change or override the internal decision of the solver
225: Logically Collective
227: Input Parameters:
228: + snes - the nonlinear solver object
229: . func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRPostCheck()`
230: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
232: Level: intermediate
234: Note:
235: This function is called BEFORE the function evaluation within the solver while the function set in
236: `SNESLineSearchSetPostCheck()` is called AFTER the function evaluation.
238: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRPostCheck()`, `SNESNewtonTRGetPostCheck()`, `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRGetPreCheck()`
239: @*/
240: PetscErrorCode SNESNewtonTRSetPostCheck(SNES snes, PetscErrorCode (*func)(SNES, Vec, Vec, Vec, PetscBool *, PetscBool *, void *), void *ctx)
241: {
242: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
243: PetscBool flg;
245: PetscFunctionBegin;
247: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
248: if (flg) {
249: if (func) tr->postcheck = func;
250: if (ctx) tr->postcheckctx = ctx;
251: }
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: /*@C
256: SNESNewtonTRGetPostCheck - Gets the post-check function
258: Not Collective
260: Input Parameter:
261: . snes - the nonlinear solver context
263: Output Parameters:
264: + func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRPostCheck()`
265: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
267: Level: intermediate
269: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRPostCheck()`
270: @*/
271: PetscErrorCode SNESNewtonTRGetPostCheck(SNES snes, PetscErrorCode (**func)(SNES, Vec, Vec, Vec, PetscBool *, PetscBool *, void *), void **ctx)
272: {
273: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
274: PetscBool flg;
276: PetscFunctionBegin;
278: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
279: PetscAssert(flg, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)snes)->type_name);
280: if (func) *func = tr->postcheck;
281: if (ctx) *ctx = tr->postcheckctx;
282: PetscFunctionReturn(PETSC_SUCCESS);
283: }
285: /*@C
286: SNESNewtonTRPreCheck - Runs the precheck routine
288: Logically Collective
290: Input Parameters:
291: + snes - the solver
292: . X - The last solution
293: - Y - The step direction
295: Output Parameter:
296: . changed_Y - Indicator that the step direction `Y` has been changed.
298: Level: intermediate
300: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRGetPreCheck()`, `SNESNewtonTRPostCheck()`
301: @*/
302: PetscErrorCode SNESNewtonTRPreCheck(SNES snes, Vec X, Vec Y, PetscBool *changed_Y)
303: {
304: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
305: PetscBool flg;
307: PetscFunctionBegin;
309: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
310: PetscAssert(flg, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)snes)->type_name);
311: *changed_Y = PETSC_FALSE;
312: if (tr->precheck) {
313: PetscCall((*tr->precheck)(snes, X, Y, changed_Y, tr->precheckctx));
315: }
316: PetscFunctionReturn(PETSC_SUCCESS);
317: }
319: /*@C
320: SNESNewtonTRPostCheck - Runs the postcheck routine
322: Logically Collective
324: Input Parameters:
325: + snes - the solver
326: . X - The last solution
327: . Y - The full step direction
328: - W - The updated solution, W = X - Y
330: Output Parameters:
331: + changed_Y - indicator if step has been changed
332: - changed_W - Indicator if the new candidate solution W has been changed.
334: Note:
335: If Y is changed then W is recomputed as X - Y
337: Level: intermediate
339: .seealso: [](ch_snes), `SNESNEWTONTR`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()`, `SNESNewtonTRPreCheck()`
340: @*/
341: PetscErrorCode SNESNewtonTRPostCheck(SNES snes, Vec X, Vec Y, Vec W, PetscBool *changed_Y, PetscBool *changed_W)
342: {
343: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
344: PetscBool flg;
346: PetscFunctionBegin;
348: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
349: PetscAssert(flg, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONG, "Not for type %s", ((PetscObject)snes)->type_name);
350: *changed_Y = PETSC_FALSE;
351: *changed_W = PETSC_FALSE;
352: if (tr->postcheck) {
353: PetscCall((*tr->postcheck)(snes, X, Y, W, changed_Y, changed_W, tr->postcheckctx));
356: }
357: PetscFunctionReturn(PETSC_SUCCESS);
358: }
360: /* stable implementation of roots of a*x^2 + b*x + c = 0 */
361: static inline void PetscQuadraticRoots(PetscReal a, PetscReal b, PetscReal c, PetscReal *xm, PetscReal *xp)
362: {
363: PetscReal temp = -0.5 * (b + PetscCopysignReal(1.0, b) * PetscSqrtReal(b * b - 4 * a * c));
364: PetscReal x1 = temp / a;
365: PetscReal x2 = c / temp;
366: *xm = PetscMin(x1, x2);
367: *xp = PetscMax(x1, x2);
368: }
370: /* Computes the quadratic model difference */
371: static PetscErrorCode SNESNewtonTRQuadraticDelta(SNES snes, Mat J, PetscBool has_objective, Vec Y, Vec GradF, Vec W, PetscReal *yTHy_, PetscReal *gTy_, PetscReal *deltaqm)
372: {
373: PetscReal yTHy, gTy;
375: PetscFunctionBegin;
376: PetscCall(MatMult(J, Y, W));
377: if (has_objective) PetscCall(VecDotRealPart(Y, W, &yTHy));
378: else PetscCall(VecDotRealPart(W, W, &yTHy)); /* Gauss-Newton approximation J^t * J */
379: PetscCall(VecDotRealPart(GradF, Y, &gTy));
380: *deltaqm = -(-(gTy) + 0.5 * (yTHy)); /* difference in quadratic model, -gTy because SNES solves it this way */
381: if (yTHy_) *yTHy_ = yTHy;
382: if (gTy_) *gTy_ = gTy;
383: PetscFunctionReturn(PETSC_SUCCESS);
384: }
386: /* Computes the new objective given X = Xk, Y = direction
387: W work vector, on output W = X - Y
388: G work vector, on output G = SNESFunction(W) */
389: static PetscErrorCode SNESNewtonTRObjective(SNES snes, PetscBool has_objective, Vec X, Vec Y, Vec W, Vec G, PetscReal *gnorm, PetscReal *fkp1)
390: {
391: PetscBool changed_y, changed_w;
393: PetscFunctionBegin;
394: /* TODO: we can add a linesearch here */
395: PetscCall(SNESNewtonTRPreCheck(snes, X, Y, &changed_y));
396: PetscCall(VecWAXPY(W, -1.0, Y, X)); /* Xkp1 */
397: PetscCall(SNESNewtonTRPostCheck(snes, X, Y, W, &changed_y, &changed_w));
398: if (changed_y && !changed_w) PetscCall(VecWAXPY(W, -1.0, Y, X));
400: PetscCall(SNESComputeFunction(snes, W, G)); /* F(Xkp1) = G */
401: PetscCall(VecNorm(G, NORM_2, gnorm));
402: if (has_objective) PetscCall(SNESComputeObjective(snes, W, fkp1));
403: else *fkp1 = 0.5 * PetscSqr(*gnorm);
404: PetscFunctionReturn(PETSC_SUCCESS);
405: }
407: static PetscErrorCode SNESSetUpQN_NEWTONTR(SNES snes)
408: {
409: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
411: PetscFunctionBegin;
412: PetscCall(MatDestroy(&tr->qnB));
413: PetscCall(MatDestroy(&tr->qnB_pre));
414: if (tr->qn) {
415: PetscInt n, N;
416: const char *optionsprefix;
417: Mat B;
419: PetscCall(MatCreate(PetscObjectComm((PetscObject)snes), &B));
420: PetscCall(SNESGetOptionsPrefix(snes, &optionsprefix));
421: PetscCall(MatSetOptionsPrefix(B, "snes_tr_qn_"));
422: PetscCall(MatAppendOptionsPrefix(B, optionsprefix));
423: PetscCall(MatSetType(B, MATLMVMBFGS));
424: PetscCall(VecGetLocalSize(snes->vec_sol, &n));
425: PetscCall(VecGetSize(snes->vec_sol, &N));
426: PetscCall(MatSetSizes(B, n, n, N, N));
427: PetscCall(MatSetUp(B));
428: PetscCall(MatSetFromOptions(B));
429: PetscCall(MatLMVMAllocate(B, snes->vec_sol, snes->vec_func));
430: tr->qnB = B;
431: if (tr->qn == SNES_TR_QN_DIFFERENT) {
432: PetscCall(MatCreate(PetscObjectComm((PetscObject)snes), &B));
433: PetscCall(SNESGetOptionsPrefix(snes, &optionsprefix));
434: PetscCall(MatSetOptionsPrefix(B, "snes_tr_qn_pre_"));
435: PetscCall(MatAppendOptionsPrefix(B, optionsprefix));
436: PetscCall(MatSetType(B, MATLMVMBFGS));
437: PetscCall(MatSetSizes(B, n, n, N, N));
438: PetscCall(MatSetUp(B));
439: PetscCall(MatSetFromOptions(B));
440: PetscCall(MatLMVMAllocate(B, snes->vec_sol, snes->vec_func));
441: tr->qnB_pre = B;
442: } else {
443: PetscCall(PetscObjectReference((PetscObject)tr->qnB));
444: tr->qnB_pre = tr->qnB;
445: }
446: }
447: PetscFunctionReturn(PETSC_SUCCESS);
448: }
450: /*
451: SNESSolve_NEWTONTR - Implements Newton's Method with trust-region subproblem and adds dogleg Cauchy
452: (Steepest Descent direction) step and direction if the trust region is not satisfied for solving system of
453: nonlinear equations
455: */
456: static PetscErrorCode SNESSolve_NEWTONTR(SNES snes)
457: {
458: SNES_NEWTONTR *neP = (SNES_NEWTONTR *)snes->data;
459: Vec X, F, Y, G, W, GradF, YU, Yc;
460: PetscInt maxits, lits;
461: PetscReal rho, fnorm, gnorm = 0.0, xnorm = 0.0, delta, ynorm;
462: PetscReal deltaM, fk, fkp1, deltaqm = 0.0, gTy = 0.0, yTHy = 0.0;
463: PetscReal auk, tauk, gfnorm, gfnorm_k, ycnorm, gTBg, objmin = 0.0, beta_k = 1.0;
464: PC pc;
465: Mat J, Jp;
466: PetscBool already_done = PETSC_FALSE, on_boundary;
467: PetscBool clear_converged_test, rho_satisfied, has_objective;
468: SNES_TR_KSPConverged_Ctx *ctx;
469: void *convctx;
471: PetscErrorCode (*convtest)(KSP, PetscInt, PetscReal, KSPConvergedReason *, void *), (*convdestroy)(void *);
472: PetscErrorCode (*objective)(SNES, Vec, PetscReal *, void *);
474: PetscFunctionBegin;
475: PetscCall(SNESGetObjective(snes, &objective, NULL));
476: has_objective = objective ? PETSC_TRUE : PETSC_FALSE;
478: maxits = snes->max_its; /* maximum number of iterations */
479: X = snes->vec_sol; /* solution vector */
480: F = snes->vec_func; /* residual vector */
481: Y = snes->vec_sol_update; /* update vector */
482: G = snes->work[0]; /* updated residual */
483: W = snes->work[1]; /* temporary vector */
484: GradF = !has_objective ? snes->work[2] : snes->vec_func; /* grad f = J^T F */
485: YU = snes->work[3]; /* work vector for dogleg method */
486: Yc = snes->work[4]; /* Cauchy point */
488: 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);
490: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
491: snes->iter = 0;
492: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
494: /* setup QN matrices if needed */
495: PetscCall(SNESSetUpQN_NEWTONTR(snes));
497: /* Set the linear stopping criteria to use the More' trick if needed */
498: clear_converged_test = PETSC_FALSE;
499: PetscCall(SNESGetKSP(snes, &snes->ksp));
500: PetscCall(KSPGetConvergenceTest(snes->ksp, &convtest, &convctx, &convdestroy));
501: if (convtest != SNESTR_KSPConverged_Private) {
502: clear_converged_test = PETSC_TRUE;
503: PetscCall(PetscNew(&ctx));
504: ctx->snes = snes;
505: PetscCall(KSPGetAndClearConvergenceTest(snes->ksp, &ctx->convtest, &ctx->convctx, &ctx->convdestroy));
506: PetscCall(KSPSetConvergenceTest(snes->ksp, SNESTR_KSPConverged_Private, ctx, SNESTR_KSPConverged_Destroy));
507: PetscCall(PetscInfo(snes, "Using Krylov convergence test SNESTR_KSPConverged_Private\n"));
508: }
510: if (!snes->vec_func_init_set) {
511: PetscCall(SNESComputeFunction(snes, X, F)); /* F(X) */
512: } else snes->vec_func_init_set = PETSC_FALSE;
514: PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- || F || */
515: SNESCheckFunctionNorm(snes, fnorm);
516: PetscCall(VecNorm(X, NORM_2, &xnorm)); /* xnorm <- || X || */
518: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
519: snes->norm = fnorm;
520: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
521: delta = neP->delta0;
522: deltaM = neP->deltaM;
523: neP->delta = delta;
524: PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));
526: /* test convergence */
527: rho_satisfied = PETSC_FALSE;
528: PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm));
529: PetscCall(SNESMonitor(snes, 0, fnorm));
530: if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
532: if (has_objective) PetscCall(SNESComputeObjective(snes, X, &fk));
533: else fk = 0.5 * PetscSqr(fnorm); /* obj(x) = 0.5 * ||F(x)||^2 */
535: /* hook state vector to BFGS preconditioner */
536: PetscCall(KSPGetPC(snes->ksp, &pc));
537: PetscCall(PCLMVMSetUpdateVec(pc, X));
539: if (neP->kmdc) PetscCall(KSPSetComputeEigenvalues(snes->ksp, PETSC_TRUE));
541: while (snes->iter < maxits) {
542: /* calculating Jacobian and GradF of minimization function only once */
543: if (!already_done) {
544: /* Call general purpose update function */
545: PetscTryTypeMethod(snes, update, snes->iter);
547: /* apply the nonlinear preconditioner */
548: if (snes->npc && snes->npcside == PC_RIGHT) {
549: SNESConvergedReason reason;
551: PetscCall(SNESSetInitialFunction(snes->npc, F));
552: PetscCall(PetscLogEventBegin(SNES_NPCSolve, snes->npc, X, snes->vec_rhs, 0));
553: PetscCall(SNESSolve(snes->npc, snes->vec_rhs, X));
554: PetscCall(PetscLogEventEnd(SNES_NPCSolve, snes->npc, X, snes->vec_rhs, 0));
555: PetscCall(SNESGetConvergedReason(snes->npc, &reason));
556: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
557: snes->reason = SNES_DIVERGED_INNER;
558: PetscFunctionReturn(PETSC_SUCCESS);
559: }
560: // XXX
561: PetscCall(SNESGetNPCFunction(snes, F, &fnorm));
562: } else if (snes->ops->update) { /* if update is present, recompute objective function and function norm */
563: PetscCall(SNESComputeFunction(snes, X, F));
564: }
566: /* Jacobian */
567: J = NULL;
568: Jp = NULL;
569: if (!neP->qnB) {
570: PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
571: J = snes->jacobian;
572: Jp = snes->jacobian_pre;
573: } else { /* QN model */
574: PetscCall(SNESComputeJacobian_MATLMVM(snes, X, neP->qnB, neP->qnB_pre, NULL));
575: J = neP->qnB;
576: Jp = neP->qnB_pre;
577: }
578: SNESCheckJacobianDomainerror(snes);
580: /* objective function */
581: PetscCall(VecNorm(F, NORM_2, &fnorm));
582: if (has_objective) PetscCall(SNESComputeObjective(snes, X, &fk));
583: else fk = 0.5 * PetscSqr(fnorm); /* obj(x) = 0.5 * ||F(x)||^2 */
585: /* GradF */
586: if (has_objective) gfnorm = fnorm;
587: else {
588: PetscCall(MatMultTranspose(J, F, GradF)); /* grad f = J^T F */
589: PetscCall(VecNorm(GradF, NORM_2, &gfnorm));
590: }
591: PetscCall(VecNorm(GradF, neP->norm, &gfnorm_k));
592: }
593: already_done = PETSC_TRUE;
595: /* solve trust-region subproblem */
597: /* first compute Cauchy Point */
598: PetscCall(MatMult(J, GradF, W));
599: if (has_objective) PetscCall(VecDotRealPart(GradF, W, &gTBg));
600: else PetscCall(VecDotRealPart(W, W, &gTBg)); /* B = J^t * J */
601: /* Eqs 4.11 and 4.12 in Nocedal and Wright 2nd Edition (4.7 and 4.8 in 1st Edition) */
602: auk = delta / gfnorm_k;
603: if (gTBg < 0.0) tauk = 1.0;
604: else tauk = PetscMin(gfnorm * gfnorm * gfnorm_k / (delta * gTBg), 1);
605: auk *= tauk;
606: ycnorm = auk * gfnorm;
607: PetscCall(VecAXPBY(Yc, auk, 0.0, GradF));
609: on_boundary = PETSC_FALSE;
610: if (tauk != 1.0) {
611: KSPConvergedReason reason;
613: /* sufficient decrease (see 6.3.27 in Conn, Gould, Toint "Trust Region Methods")
614: beta_k the largest eigenvalue of the Hessian. Here we use the previous estimated value */
615: objmin = -neP->kmdc * gnorm * PetscMin(gnorm / beta_k, delta);
616: PetscCall(KSPCGSetObjectiveTarget(snes->ksp, objmin));
618: /* specify radius if looking for Newton step and trust region norm is the l2 norm */
619: PetscCall(KSPCGSetRadius(snes->ksp, neP->fallback == SNES_TR_FALLBACK_NEWTON && neP->norm == NORM_2 ? delta : 0.0));
620: PetscCall(KSPSetOperators(snes->ksp, J, Jp));
621: PetscCall(KSPSolve(snes->ksp, F, Y));
622: SNESCheckKSPSolve(snes);
623: PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
624: PetscCall(KSPGetConvergedReason(snes->ksp, &reason));
625: on_boundary = (PetscBool)(reason == KSP_CONVERGED_STEP_LENGTH);
626: PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits));
627: if (neP->kmdc) { /* update estimated Hessian largest eigenvalue */
628: PetscReal emax, emin;
629: PetscCall(KSPComputeExtremeSingularValues(snes->ksp, &emax, &emin));
630: if (emax > 0.0) beta_k = emax + 1;
631: }
632: } else { /* Cauchy point is on the boundary, accept it */
633: on_boundary = PETSC_TRUE;
634: PetscCall(VecCopy(Yc, Y));
635: PetscCall(PetscInfo(snes, "CP evaluated on boundary. delta: %g, ycnorm: %g, gTBg: %g\n", (double)delta, (double)ycnorm, (double)gTBg));
636: }
637: PetscCall(VecNorm(Y, neP->norm, &ynorm));
639: /* decide what to do when the update is outside of trust region */
640: if (ynorm > delta || ynorm == 0.0) {
641: SNESNewtonTRFallbackType fallback = ynorm > 0.0 ? neP->fallback : SNES_TR_FALLBACK_CAUCHY;
643: PetscCheck(neP->norm == NORM_2 || fallback != SNES_TR_FALLBACK_DOGLEG, PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "DOGLEG without l2 norm not implemented");
644: switch (fallback) {
645: case SNES_TR_FALLBACK_NEWTON:
646: auk = delta / ynorm;
647: PetscCall(VecScale(Y, auk));
648: PetscCall(PetscInfo(snes, "SN evaluated. delta: %g, ynorm: %g\n", (double)delta, (double)ynorm));
649: break;
650: case SNES_TR_FALLBACK_CAUCHY:
651: case SNES_TR_FALLBACK_DOGLEG:
652: if (fallback == SNES_TR_FALLBACK_CAUCHY || gTBg <= 0.0) {
653: PetscCall(VecCopy(Yc, Y));
654: PetscCall(PetscInfo(snes, "CP evaluated. delta: %g, ynorm: %g, ycnorm: %g, gTBg: %g\n", (double)delta, (double)ynorm, (double)ycnorm, (double)gTBg));
655: } else { /* take linear combination of Cauchy and Newton direction and step */
656: auk = gfnorm * gfnorm / gTBg;
657: if (gfnorm_k * auk >= delta) { /* first leg: Cauchy point outside of trust region */
658: PetscCall(VecAXPBY(Y, delta / gfnorm_k, 0.0, GradF));
659: PetscCall(PetscInfo(snes, "CP evaluated (outside region). delta: %g, ynorm: %g, ycnorm: %g\n", (double)delta, (double)ynorm, (double)ycnorm));
660: } else { /* second leg */
661: PetscReal c0, c1, c2, tau = 0.0, tpos, tneg;
662: PetscBool noroots;
664: /* Find solutions of (Eq. 4.16 in Nocedal and Wright)
665: ||p_U + lambda * (p_B - p_U)||^2 - delta^2 = 0,
666: where p_U the Cauchy direction, p_B the Newton direction */
667: PetscCall(VecAXPBY(YU, auk, 0.0, GradF));
668: PetscCall(VecAXPY(Y, -1.0, YU));
669: PetscCall(VecNorm(Y, NORM_2, &c0));
670: PetscCall(VecDotRealPart(YU, Y, &c1));
671: c0 = PetscSqr(c0);
672: c2 = PetscSqr(ycnorm) - PetscSqr(delta);
673: PetscQuadraticRoots(c0, 2 * c1, c2, &tneg, &tpos);
675: /* In principle the DL strategy as a unique solution in [0,1]
676: here we check that for some reason we numerically failed
677: to compute it. In that case, we use the Cauchy point */
678: noroots = PetscIsInfOrNanReal(tneg);
679: if (!noroots) {
680: if (tpos > 1) {
681: if (tneg >= 0 && tneg <= 1) {
682: tau = tneg;
683: } else noroots = PETSC_TRUE;
684: } else if (tpos >= 0) {
685: tau = tpos;
686: } else noroots = PETSC_TRUE;
687: }
688: if (noroots) { /* No roots, select Cauchy point */
689: PetscCall(VecCopy(Yc, Y));
690: } else {
691: PetscCall(VecAXPBY(Y, 1.0, tau, YU));
692: }
693: 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));
694: }
695: }
696: break;
697: default:
698: SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "Unknown fallback mode");
699: break;
700: }
701: }
703: /* compute the quadratic model difference */
704: PetscCall(SNESNewtonTRQuadraticDelta(snes, J, has_objective, Y, GradF, W, &yTHy, &gTy, &deltaqm));
706: /* Compute new objective function */
707: PetscCall(SNESNewtonTRObjective(snes, has_objective, X, Y, W, G, &gnorm, &fkp1));
708: if (PetscIsInfOrNanReal(fkp1)) rho = neP->eta1;
709: else {
710: if (deltaqm > 0.0) rho = (fk - fkp1) / deltaqm; /* actual improvement over predicted improvement */
711: else rho = neP->eta1; /* no reduction in quadratic model, step must be rejected */
712: }
714: PetscCall(VecNorm(Y, neP->norm, &ynorm));
715: 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));
717: /* update the size of the trust region */
718: if (rho < neP->eta2) delta *= neP->t1; /* shrink the region */
719: else if (rho > neP->eta3 && on_boundary) delta *= neP->t2; /* expand the region */
720: delta = PetscMin(delta, deltaM); /* but not greater than deltaM */
722: /* log 2-norm of update for moniroting routines */
723: PetscCall(VecNorm(Y, NORM_2, &ynorm));
725: /* decide on new step */
726: neP->delta = delta;
727: if (rho > neP->eta1) {
728: rho_satisfied = PETSC_TRUE;
729: } else {
730: rho_satisfied = PETSC_FALSE;
731: PetscCall(PetscInfo(snes, "Trying again in smaller region\n"));
732: /* check to see if progress is hopeless */
733: PetscCall(SNESTR_Converged_Private(snes, snes->iter, xnorm, ynorm, fnorm, &snes->reason, snes->cnvP));
734: if (!snes->reason) PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm));
735: if (snes->reason == SNES_CONVERGED_SNORM_RELATIVE) snes->reason = SNES_DIVERGED_TR_DELTA;
736: snes->numFailures++;
737: /* We're not progressing, so return with the current iterate */
738: if (snes->reason) break;
739: }
740: if (rho_satisfied) {
741: /* Update function values */
742: already_done = PETSC_FALSE;
743: fnorm = gnorm;
744: fk = fkp1;
746: /* New residual and linearization point */
747: PetscCall(VecCopy(G, F));
748: PetscCall(VecCopy(W, X));
750: /* Monitor convergence */
751: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
752: snes->iter++;
753: snes->norm = fnorm;
754: snes->xnorm = xnorm;
755: snes->ynorm = ynorm;
756: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
757: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits));
759: /* Test for convergence, xnorm = || X || */
760: PetscCall(VecNorm(X, NORM_2, &xnorm));
761: PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm));
762: PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
763: if (snes->reason) break;
764: }
765: }
767: if (clear_converged_test) {
768: PetscCall(KSPGetAndClearConvergenceTest(snes->ksp, &ctx->convtest, &ctx->convctx, &ctx->convdestroy));
769: PetscCall(PetscFree(ctx));
770: PetscCall(KSPSetConvergenceTest(snes->ksp, convtest, convctx, convdestroy));
771: }
772: PetscFunctionReturn(PETSC_SUCCESS);
773: }
775: static PetscErrorCode SNESSetUp_NEWTONTR(SNES snes)
776: {
777: PetscFunctionBegin;
778: PetscCall(SNESSetWorkVecs(snes, 5));
779: PetscCall(SNESSetUpMatrices(snes));
780: PetscFunctionReturn(PETSC_SUCCESS);
781: }
783: static PetscErrorCode SNESReset_NEWTONTR(SNES snes)
784: {
785: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
787: PetscFunctionBegin;
788: PetscCall(MatDestroy(&tr->qnB));
789: PetscCall(MatDestroy(&tr->qnB_pre));
790: PetscFunctionReturn(PETSC_SUCCESS);
791: }
793: static PetscErrorCode SNESDestroy_NEWTONTR(SNES snes)
794: {
795: PetscFunctionBegin;
796: PetscCall(SNESReset_NEWTONTR(snes));
797: PetscCall(PetscFree(snes->data));
798: PetscFunctionReturn(PETSC_SUCCESS);
799: }
801: static PetscErrorCode SNESSetFromOptions_NEWTONTR(SNES snes, PetscOptionItems *PetscOptionsObject)
802: {
803: SNES_NEWTONTR *ctx = (SNES_NEWTONTR *)snes->data;
804: SNESNewtonTRQNType qn;
805: SNESNewtonTRFallbackType fallback;
806: NormType norm;
807: PetscReal deltatol;
808: PetscBool flg;
810: PetscFunctionBegin;
811: PetscOptionsHeadBegin(PetscOptionsObject, "SNES trust region options for nonlinear equations");
812: PetscCall(PetscOptionsReal("-snes_tr_eta1", "eta1", "None", ctx->eta1, &ctx->eta1, NULL));
813: PetscCall(PetscOptionsReal("-snes_tr_eta2", "eta2", "None", ctx->eta2, &ctx->eta2, NULL));
814: PetscCall(PetscOptionsReal("-snes_tr_eta3", "eta3", "None", ctx->eta3, &ctx->eta3, NULL));
815: PetscCall(PetscOptionsReal("-snes_tr_t1", "t1", "None", ctx->t1, &ctx->t1, NULL));
816: PetscCall(PetscOptionsReal("-snes_tr_t2", "t2", "None", ctx->t2, &ctx->t2, NULL));
817: PetscCall(PetscOptionsReal("-snes_tr_deltaM", "deltaM", "None", ctx->deltaM, &ctx->deltaM, NULL));
818: PetscCall(PetscOptionsReal("-snes_tr_delta0", "delta0", "None", ctx->delta0, &ctx->delta0, NULL));
819: PetscCall(PetscOptionsReal("-snes_tr_kmdc", "sufficient decrease parameter", "None", ctx->kmdc, &ctx->kmdc, NULL));
821: deltatol = snes->deltatol;
822: PetscCall(PetscOptionsReal("-snes_tr_tol", "Trust region tolerance", "SNESSetTrustRegionTolerance", deltatol, &deltatol, &flg));
823: if (flg) PetscCall(SNESSetTrustRegionTolerance(snes, deltatol));
825: fallback = ctx->fallback;
826: 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));
827: if (flg) PetscCall(SNESNewtonTRSetFallbackType(snes, fallback));
829: qn = ctx->qn;
830: PetscCall(PetscOptionsEnum("-snes_tr_qn", "Use Quasi-Newton approximations for the model", "SNESNewtonTRSetQNType", SNESNewtonTRQNTypes, (PetscEnum)qn, (PetscEnum *)&qn, &flg));
831: if (flg) PetscCall(SNESNewtonTRSetQNType(snes, qn));
833: norm = ctx->norm;
834: PetscCall(PetscOptionsEnum("-snes_tr_norm_type", "Type of norm for trust region bounds", "SNESNewtonTRSetNormType", NormTypes, (PetscEnum)norm, (PetscEnum *)&norm, &flg));
835: if (flg) PetscCall(SNESNewtonTRSetNormType(snes, norm));
837: PetscOptionsHeadEnd();
838: PetscFunctionReturn(PETSC_SUCCESS);
839: }
841: static PetscErrorCode SNESView_NEWTONTR(SNES snes, PetscViewer viewer)
842: {
843: SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
844: PetscBool iascii;
846: PetscFunctionBegin;
847: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
848: if (iascii) {
849: PetscCall(PetscViewerASCIIPrintf(viewer, " Trust region tolerance %g\n", (double)snes->deltatol));
850: PetscCall(PetscViewerASCIIPrintf(viewer, " eta1=%g, eta2=%g, eta3=%g\n", (double)tr->eta1, (double)tr->eta2, (double)tr->eta3));
851: PetscCall(PetscViewerASCIIPrintf(viewer, " delta0=%g, t1=%g, t2=%g, deltaM=%g\n", (double)tr->delta0, (double)tr->t1, (double)tr->t2, (double)tr->deltaM));
852: PetscCall(PetscViewerASCIIPrintf(viewer, " kmdc=%g\n", (double)tr->kmdc));
853: PetscCall(PetscViewerASCIIPrintf(viewer, " fallback=%s\n", SNESNewtonTRFallbackTypes[tr->fallback]));
854: if (tr->qn) PetscCall(PetscViewerASCIIPrintf(viewer, " qn=%s\n", SNESNewtonTRQNTypes[tr->qn]));
855: if (tr->norm != NORM_2) PetscCall(PetscViewerASCIIPrintf(viewer, " norm=%s\n", NormTypes[tr->norm]));
856: }
857: PetscFunctionReturn(PETSC_SUCCESS);
858: }
860: /*MC
861: SNESNEWTONTR - Newton based nonlinear solver that uses trust-region dogleg method with Cauchy direction {cite}`nocedal2006numerical`
863: Options Database Keys:
864: + -snes_tr_tol <tol> - trust region tolerance
865: . -snes_tr_eta1 <eta1> - trust region parameter eta1 <= eta2, rho > eta1 breaks out of the inner iteration (default: eta1=0.001)
866: . -snes_tr_eta2 <eta2> - trust region parameter, rho <= eta2 shrinks the trust region (default: eta2=0.25)
867: . -snes_tr_eta3 <eta3> - trust region parameter eta3 > eta2, rho >= eta3 expands the trust region (default: eta3=0.75)
868: . -snes_tr_t1 <t1> - trust region parameter, shrinking factor of trust region (default: 0.25)
869: . -snes_tr_t2 <t2> - trust region parameter, expanding factor of trust region (default: 2.0)
870: . -snes_tr_deltaM <deltaM> - trust region parameter, max size of trust region (default: MAX_REAL)
871: . -snes_tr_delta0 <delta0> - trust region parameter, initial size of trust region (default: 0.2)
872: - -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.
874: Level: beginner
876: .seealso: [](ch_snes), `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESSetTrustRegionTolerance()`,
877: `SNESNewtonTRPreCheck()`, `SNESNewtonTRGetPreCheck()`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()`,
878: `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRSetFallbackType()`, `SNESNewtonTRSetQNType()`
879: M*/
880: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONTR(SNES snes)
881: {
882: SNES_NEWTONTR *neP;
884: PetscFunctionBegin;
885: snes->ops->setup = SNESSetUp_NEWTONTR;
886: snes->ops->solve = SNESSolve_NEWTONTR;
887: snes->ops->reset = SNESReset_NEWTONTR;
888: snes->ops->destroy = SNESDestroy_NEWTONTR;
889: snes->ops->setfromoptions = SNESSetFromOptions_NEWTONTR;
890: snes->ops->view = SNESView_NEWTONTR;
892: snes->stol = 0.0;
893: snes->usesksp = PETSC_TRUE;
894: snes->npcside = PC_RIGHT;
895: snes->usesnpc = PETSC_TRUE;
897: snes->alwayscomputesfinalresidual = PETSC_TRUE;
899: PetscCall(PetscNew(&neP));
900: snes->data = (void *)neP;
901: neP->delta = 0.0;
902: neP->delta0 = 0.2;
903: neP->eta1 = 0.001;
904: neP->eta2 = 0.25;
905: neP->eta3 = 0.75;
906: neP->t1 = 0.25;
907: neP->t2 = 2.0;
908: neP->deltaM = 1.e10;
909: neP->norm = NORM_2;
910: neP->fallback = SNES_TR_FALLBACK_NEWTON;
911: neP->kmdc = 0.0; /* by default do not use sufficient decrease */
912: PetscFunctionReturn(PETSC_SUCCESS);
913: }