Actual source code: ntrdc.c
1: #include <../src/snes/impls/ntrdc/ntrdcimpl.h>
3: typedef struct {
4: SNES snes;
5: /* Information on the regular SNES convergence test; which may have been user provided
6: Copied from tr.c (maybe able to disposed, but this is a private function) - Heeho
7: Same with SNESTR_KSPConverged_Private, SNESTR_KSPConverged_Destroy, and SNESTR_Converged_Private
8: */
10: KSPConvergenceTestFn *convtest;
11: PetscCtxDestroyFn *convdestroy;
12: void *convctx;
13: } SNES_TRDC_KSPConverged_Ctx;
15: static PetscErrorCode SNESNewtonTRSetTolerances_TRDC(SNES snes, PetscReal delta_min, PetscReal delta_max, PetscReal delta_0)
16: {
17: SNES_NEWTONTRDC *tr = (SNES_NEWTONTRDC *)snes->data;
19: PetscFunctionBegin;
20: if (delta_min == PETSC_DETERMINE) delta_min = 1.e-12;
21: if (delta_max == PETSC_DETERMINE) delta_max = 0.5;
22: if (delta_0 == PETSC_DETERMINE) delta_0 = 0.1;
23: if (delta_min != PETSC_CURRENT) tr->deltatol = delta_min;
24: if (delta_max != PETSC_CURRENT) tr->deltaM = delta_max;
25: if (delta_0 != PETSC_CURRENT) tr->delta0 = delta_0;
26: PetscFunctionReturn(PETSC_SUCCESS);
27: }
29: static PetscErrorCode SNESTRDC_KSPConverged_Private(KSP ksp, PetscInt n, PetscReal rnorm, KSPConvergedReason *reason, PetscCtx cctx)
30: {
31: SNES_TRDC_KSPConverged_Ctx *ctx = (SNES_TRDC_KSPConverged_Ctx *)cctx;
32: SNES snes = ctx->snes;
33: SNES_NEWTONTRDC *neP = (SNES_NEWTONTRDC *)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, NORM_2, &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 SNESTRDC_KSPConverged_Destroy(PetscCtxRt cctx)
51: {
52: SNES_TRDC_KSPConverged_Ctx *ctx = *(SNES_TRDC_KSPConverged_Ctx **)cctx;
54: PetscFunctionBegin;
55: PetscCall((*ctx->convdestroy)(&ctx->convctx));
56: PetscCall(PetscFree(ctx));
57: PetscFunctionReturn(PETSC_SUCCESS);
58: }
60: /*
61: SNESTRDC_Converged_Private -test convergence JUST for the trust region tolerance.
62: */
63: static PetscErrorCode SNESTRDC_Converged_Private(SNES snes, PetscInt it, PetscReal xnorm, PetscReal pnorm, PetscReal fnorm, SNESConvergedReason *reason, void *dummy)
64: {
65: SNES_NEWTONTRDC *neP = (SNES_NEWTONTRDC *)snes->data;
67: PetscFunctionBegin;
68: *reason = SNES_CONVERGED_ITERATING;
69: if (neP->delta < xnorm * neP->deltatol) {
70: PetscCall(PetscInfo(snes, "Diverged due to too small a trust region %g<%g*%g\n", (double)neP->delta, (double)xnorm, (double)neP->deltatol));
71: *reason = SNES_DIVERGED_TR_DELTA;
72: } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
73: PetscCall(PetscInfo(snes, "Exceeded maximum number of function evaluations: %" PetscInt_FMT "\n", snes->max_funcs));
74: *reason = SNES_DIVERGED_FUNCTION_COUNT;
75: }
76: PetscFunctionReturn(PETSC_SUCCESS);
77: }
79: /*@
80: SNESNewtonTRDCGetRhoFlag - Get whether the current solution update is within the trust-region.
82: Logically Collective
84: Input Parameter:
85: . snes - the nonlinear solver object
87: Output Parameter:
88: . rho_flag - `PETSC_FALSE` or `PETSC_TRUE`
90: Level: developer
92: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTRDC`, `SNESNewtonTRDCPreCheck()`, `SNESNewtonTRDCGetPreCheck()`, `SNESNewtonTRDCSetPreCheck()`,
93: `SNESNewtonTRDCSetPostCheck()`, `SNESNewtonTRDCGetPostCheck()`
94: @*/
95: PetscErrorCode SNESNewtonTRDCGetRhoFlag(SNES snes, PetscBool *rho_flag)
96: {
97: SNES_NEWTONTRDC *tr = (SNES_NEWTONTRDC *)snes->data;
99: PetscFunctionBegin;
101: PetscAssertPointer(rho_flag, 2);
102: *rho_flag = tr->rho_satisfied;
103: PetscFunctionReturn(PETSC_SUCCESS);
104: }
106: /*@C
107: SNESNewtonTRDCSetPreCheck - Sets a user function that is called before the search step has been determined.
108: Allows the user a chance to change or override the trust region decision.
110: Logically Collective
112: Input Parameters:
113: + snes - the nonlinear solver object
114: . func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRDCPreCheck()`
115: - ctx - [optional] application context for private data for the function evaluation routine (may be `NULL`)
117: Calling sequence of `func`:
118: + snes - the nonlinear solver object
119: . X - the current solution value
120: . Y - the tentative update step
121: . changed - output, flag indicating `Y` has been changed by the pre-check
122: - ctx - the optional application context
124: Level: intermediate
126: Note:
127: This function is called BEFORE the function evaluation within the `SNESNEWTONTRDC` solver.
129: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTRDC`, `SNESNewtonTRDCPreCheck()`, `SNESNewtonTRDCGetPreCheck()`, `SNESNewtonTRDCSetPostCheck()`, `SNESNewtonTRDCGetPostCheck()`,
130: `SNESNewtonTRDCGetRhoFlag()`
131: @*/
132: PetscErrorCode SNESNewtonTRDCSetPreCheck(SNES snes, PetscErrorCode (*func)(SNES snes, Vec X, Vec Y, PetscBool *changed, PetscCtx ctx), PetscCtx ctx)
133: {
134: SNES_NEWTONTRDC *tr = (SNES_NEWTONTRDC *)snes->data;
136: PetscFunctionBegin;
138: if (func) tr->precheck = func;
139: if (ctx) tr->precheckctx = ctx;
140: PetscFunctionReturn(PETSC_SUCCESS);
141: }
143: /*@C
144: SNESNewtonTRDCGetPreCheck - Gets the pre-check function optionally set with `SNESNewtonTRDCSetPreCheck()`
146: Not Collective
148: Input Parameter:
149: . snes - the nonlinear solver context
151: Output Parameters:
152: + func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRDCPreCheck()`
153: - ctx - [optional] context for private data for the function evaluation routine (may be `NULL`)
155: Calling sequence of `func`:
156: + snes - the nonlinear solver object
157: . X - the current solution value
158: . Y - the tentative update step
159: . changed - output, flag indicating `Y` has been changed by the pre-check
160: - ctx - the optional application context
162: Level: intermediate
164: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTRDC`, `SNESNewtonTRDCSetPreCheck()`, `SNESNewtonTRDCPreCheck()`
165: @*/
166: PetscErrorCode SNESNewtonTRDCGetPreCheck(SNES snes, PetscErrorCode (**func)(SNES snes, Vec X, Vec Y, PetscBool *changed, PetscCtx ctx), PetscCtxRt ctx)
167: {
168: SNES_NEWTONTRDC *tr = (SNES_NEWTONTRDC *)snes->data;
170: PetscFunctionBegin;
172: if (func) *func = tr->precheck;
173: if (ctx) *(void **)ctx = tr->precheckctx;
174: PetscFunctionReturn(PETSC_SUCCESS);
175: }
177: /*@C
178: SNESNewtonTRDCSetPostCheck - Sets a user function that is called after the search step has been determined but before the next
179: function evaluation. Allows the user a chance to change or override the decision of the line search routine
181: Logically Collective
183: Input Parameters:
184: + snes - the nonlinear solver object
185: . func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRDCPostCheck()`
186: - ctx - [optional] context for private data for the function evaluation routine (may be `NULL`)
188: Calling sequence of `func`:
189: + snes - the nonlinear solver object
190: . X - the current solution value
191: . Y - the tentative update step
192: . W - the tentative new solution value
193: . changed_y - output, flag indicated `Y` has been changed by the post-check
194: . changed_w - output, flag indicated `W` has been changed by the post-check
195: - ctx - the optional application context
197: Level: intermediate
199: Note:
200: This function is called BEFORE the function evaluation within the `SNESNEWTONTRDC` solver while the function set in
201: `SNESLineSearchSetPostCheck()` is called AFTER the function evaluation.
203: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTRDC`, `SNESNewtonTRDCPostCheck()`, `SNESNewtonTRDCGetPostCheck()`, `SNESNewtonTRDCSetPreCheck()`, `SNESNewtonTRDCGetPreCheck()`
204: @*/
205: PetscErrorCode SNESNewtonTRDCSetPostCheck(SNES snes, PetscErrorCode (*func)(SNES snes, Vec X, Vec Y, Vec W, PetscBool *changed_y, PetscBool *changed_w, PetscCtx ctx), PetscCtx ctx)
206: {
207: SNES_NEWTONTRDC *tr = (SNES_NEWTONTRDC *)snes->data;
209: PetscFunctionBegin;
211: if (func) tr->postcheck = func;
212: if (ctx) tr->postcheckctx = ctx;
213: PetscFunctionReturn(PETSC_SUCCESS);
214: }
216: /*@C
217: SNESNewtonTRDCGetPostCheck - Gets the post-check function optionally set with `SNESNewtonTRDCSetPostCheck()`
219: Not Collective
221: Input Parameter:
222: . snes - the nonlinear solver context
224: Output Parameters:
225: + func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRDCPostCheck()`
226: - ctx - [optional] context for private data for the function evaluation routine (may be `NULL`)
228: Calling sequence of `func`:
229: + snes - the nonlinear solver object
230: . X - the current solution value
231: . Y - the tentative update step
232: . W - the tentative new solution value
233: . changed_y - output, flag indicated `Y` has been changed by the post-check
234: . changed_w - output, flag indicated `W` has been changed by the post-check
235: - ctx - the optional application context
237: Level: intermediate
239: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTRDC`, `SNESNewtonTRDCSetPostCheck()`, `SNESNewtonTRDCPostCheck()`, `SNESNewtonTRDCSetPreCheck()`, `SNESNewtonTRDCGetPreCheck()`
240: @*/
241: PetscErrorCode SNESNewtonTRDCGetPostCheck(SNES snes, PetscErrorCode (**func)(SNES snes, Vec X, Vec Y, Vec W, PetscBool *changed_y, PetscBool *changed_w, PetscCtx ctx), PetscCtxRt ctx)
242: {
243: SNES_NEWTONTRDC *tr = (SNES_NEWTONTRDC *)snes->data;
245: PetscFunctionBegin;
247: if (func) *func = tr->postcheck;
248: if (ctx) *(void **)ctx = tr->postcheckctx;
249: PetscFunctionReturn(PETSC_SUCCESS);
250: }
252: // PetscClangLinter pragma disable: -fdoc-internal-linkage
253: /*@C
254: SNESNewtonTRDCPreCheck - Called before the step has been determined in `SNESNEWTONTRDC`
256: Logically Collective
258: Input Parameters:
259: + snes - the solver
260: . X - The last solution
261: - Y - The step direction
263: Output Parameter:
264: . changed_y - Indicator that the step direction `Y` has been changed.
266: Level: developer
268: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTRDC`, `SNESNewtonTRDCSetPreCheck()`, `SNESNewtonTRDCGetPreCheck()`, `SNESNewtonTRDCPostCheck()`
269: @*/
270: static PetscErrorCode SNESNewtonTRDCPreCheck(SNES snes, Vec X, Vec Y, PetscBool *changed_Y)
271: {
272: SNES_NEWTONTRDC *tr = (SNES_NEWTONTRDC *)snes->data;
274: PetscFunctionBegin;
275: *changed_Y = PETSC_FALSE;
276: if (tr->precheck) {
277: PetscCall((*tr->precheck)(snes, X, Y, changed_Y, tr->precheckctx));
279: }
280: PetscFunctionReturn(PETSC_SUCCESS);
281: }
283: // PetscClangLinter pragma disable: -fdoc-internal-linkage
284: /*@C
285: SNESNewtonTRDCPostCheck - Called after the step has been determined in `SNESNEWTONTRDC` but before the function evaluation at that step
287: Logically Collective
289: Input Parameters:
290: + snes - the solver
291: . X - The last solution
292: . Y - The full step direction
293: - W - The updated solution, $W = X - Y$
295: Output Parameters:
296: + changed_Y - indicator if the step `Y` has been changed
297: - changed_W - Indicator if the new candidate solution `W` has been changed.
299: Level: developer
301: Note:
302: If `Y` is changed then `W` is recomputed as `X` - `Y`
304: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTRDC`, `SNESNEWTONTRDC`, `SNESNewtonTRDCSetPostCheck()`, `SNESNewtonTRDCGetPostCheck()`, `SNESNewtonTRDCPreCheck()
305: @*/
306: static PetscErrorCode SNESNewtonTRDCPostCheck(SNES snes, Vec X, Vec Y, Vec W, PetscBool *changed_Y, PetscBool *changed_W)
307: {
308: SNES_NEWTONTRDC *tr = (SNES_NEWTONTRDC *)snes->data;
310: PetscFunctionBegin;
311: *changed_Y = PETSC_FALSE;
312: *changed_W = PETSC_FALSE;
313: if (tr->postcheck) {
314: PetscCall((*tr->postcheck)(snes, X, Y, W, changed_Y, changed_W, tr->postcheckctx));
317: }
318: PetscFunctionReturn(PETSC_SUCCESS);
319: }
321: /*
322: SNESSolve_NEWTONTRDC - Implements Newton's Method with trust-region subproblem and adds dogleg Cauchy
323: (Steepest Descent direction) step and direction if the trust region is not satisfied for solving system of
324: nonlinear equations
326: */
327: static PetscErrorCode SNESSolve_NEWTONTRDC(SNES snes)
328: {
329: SNES_NEWTONTRDC *neP = (SNES_NEWTONTRDC *)snes->data;
330: Vec X, F, Y, G, W, GradF, YNtmp;
331: Vec YCtmp;
332: Mat jac;
333: PetscInt maxits, i, j, lits, inner_count, bs;
334: PetscReal rho, fnorm, gnorm, xnorm = 0, delta, ynorm, temp_xnorm, temp_ynorm; /* TRDC inner iteration */
335: PetscReal inorms[99]; /* need to make it dynamic eventually, fixed max block size of 99 for now */
336: PetscReal deltaM, ynnorm, f0, mp, gTy, g, yTHy; /* rho calculation */
337: PetscReal auk, gfnorm, ycnorm, c0, c1, c2, tau, tau_pos, tau_neg, gTBg; /* Cauchy Point */
338: KSP ksp;
339: SNESConvergedReason reason = SNES_CONVERGED_ITERATING;
340: PetscBool breakout = PETSC_FALSE;
341: SNES_TRDC_KSPConverged_Ctx *ctx;
342: KSPConvergenceTestFn *convtest;
343: PetscCtxDestroyFn *convdestroy;
344: void *convctx;
346: PetscFunctionBegin;
347: maxits = snes->max_its; /* maximum number of iterations */
348: X = snes->vec_sol; /* solution vector */
349: F = snes->vec_func; /* residual vector */
350: Y = snes->work[0]; /* update vector */
351: G = snes->work[1]; /* updated residual */
352: W = snes->work[2]; /* temporary vector */
353: GradF = snes->work[3]; /* grad f = J^T F */
354: YNtmp = snes->work[4]; /* Newton solution */
355: YCtmp = snes->work[5]; /* Cauchy solution */
357: 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);
359: PetscCall(VecGetBlockSize(YNtmp, &bs));
361: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
362: snes->iter = 0;
363: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
365: /* Set the linear stopping criteria to use the More' trick. From tr.c */
366: PetscCall(SNESGetKSP(snes, &ksp));
367: PetscCall(KSPGetConvergenceTest(ksp, &convtest, &convctx, &convdestroy));
368: if (convtest != SNESTRDC_KSPConverged_Private) {
369: PetscCall(PetscNew(&ctx));
370: ctx->snes = snes;
371: PetscCall(KSPGetAndClearConvergenceTest(ksp, &ctx->convtest, &ctx->convctx, &ctx->convdestroy));
372: PetscCall(KSPSetConvergenceTest(ksp, SNESTRDC_KSPConverged_Private, ctx, SNESTRDC_KSPConverged_Destroy));
373: PetscCall(PetscInfo(snes, "Using Krylov convergence test SNESTRDC_KSPConverged_Private\n"));
374: }
376: if (!snes->vec_func_init_set) {
377: PetscCall(SNESComputeFunction(snes, X, F)); /* F(X) */
378: } else snes->vec_func_init_set = PETSC_FALSE;
380: PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- || F || */
381: SNESCheckFunctionDomainError(snes, fnorm);
382: PetscCall(VecNorm(X, NORM_2, &xnorm)); /* xnorm <- || X || */
384: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
385: snes->norm = fnorm;
386: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
387: delta = xnorm ? neP->delta0 * xnorm : neP->delta0; /* initial trust region size scaled by xnorm */
388: deltaM = xnorm ? neP->deltaM * xnorm : neP->deltaM; /* maximum trust region size scaled by xnorm */
389: neP->delta = delta;
390: PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));
391: PetscCall(SNESMonitor(snes, 0, fnorm));
393: neP->rho_satisfied = PETSC_FALSE;
395: /* test convergence */
396: PetscUseTypeMethod(snes, converged, snes->iter, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
397: if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
399: for (i = 0; i < maxits; i++) {
400: PetscBool changed_y;
401: PetscBool changed_w;
403: /* dogleg method */
404: PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
405: SNESCheckJacobianDomainError(snes);
406: PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian));
407: PetscCall(KSPSolve(snes->ksp, F, YNtmp)); /* Quasi Newton Solution */
408: SNESCheckKSPSolve(snes); /* this is necessary but old tr.c did not have it*/
409: PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
410: PetscCall(SNESGetJacobian(snes, &jac, NULL, NULL, NULL));
412: /* rescale Jacobian, Newton solution update, and re-calculate delta for multiphase (multivariable)
413: for inner iteration and Cauchy direction calculation
414: */
415: if (bs > 1 && neP->auto_scale_multiphase) {
416: PetscCall(VecStrideNormAll(YNtmp, NORM_INFINITY, inorms));
417: for (j = 0; j < bs; j++) {
418: if (neP->auto_scale_max > 1.0) {
419: if (inorms[j] < 1.0 / neP->auto_scale_max) inorms[j] = 1.0 / neP->auto_scale_max;
420: }
421: PetscCall(VecStrideSet(W, j, inorms[j]));
422: PetscCall(VecStrideScale(YNtmp, j, 1.0 / inorms[j]));
423: PetscCall(VecStrideScale(X, j, 1.0 / inorms[j]));
424: }
425: PetscCall(VecNorm(X, NORM_2, &xnorm));
426: if (i == 0) {
427: delta = neP->delta0 * xnorm;
428: } else {
429: delta = neP->delta * xnorm;
430: }
431: deltaM = neP->deltaM * xnorm;
432: PetscCall(MatDiagonalScale(jac, NULL, W));
433: }
435: /* calculating GradF of minimization function */
436: PetscCall(MatMultTranspose(jac, F, GradF)); /* grad f = J^T F */
437: PetscCall(VecNorm(YNtmp, NORM_2, &ynnorm)); /* ynnorm <- || Y_newton || */
439: inner_count = 0;
440: neP->rho_satisfied = PETSC_FALSE;
441: while (1) {
442: if (ynnorm <= delta) { /* see if the Newton solution is within the trust region */
443: PetscCall(VecCopy(YNtmp, Y));
444: } else if (neP->use_cauchy) { /* use Cauchy direction if enabled */
445: PetscCall(MatMult(jac, GradF, W));
446: PetscCall(VecDotRealPart(W, W, &gTBg)); /* completes GradF^T J^T J GradF */
447: PetscCall(VecNorm(GradF, NORM_2, &gfnorm)); /* grad f norm <- || grad f || */
448: if (gTBg <= 0.0) {
449: auk = PETSC_MAX_REAL;
450: } else {
451: auk = PetscSqr(gfnorm) / gTBg;
452: }
453: auk = PetscMin(delta / gfnorm, auk);
454: PetscCall(VecCopy(GradF, YCtmp)); /* this could be improved */
455: PetscCall(VecScale(YCtmp, auk)); /* YCtmp, Cauchy solution*/
456: PetscCall(VecNorm(YCtmp, NORM_2, &ycnorm)); /* ycnorm <- || Y_cauchy || */
457: if (ycnorm >= delta) { /* see if the Cauchy solution meets the criteria */
458: PetscCall(VecCopy(YCtmp, Y));
459: PetscCall(PetscInfo(snes, "DL evaluated. delta: %8.4e, ynnorm: %8.4e, ycnorm: %8.4e\n", (double)delta, (double)ynnorm, (double)ycnorm));
460: } else { /* take ratio, tau, of Cauchy and Newton direction and step */
461: PetscCall(VecAXPY(YNtmp, -1.0, YCtmp)); /* YCtmp = A, YNtmp = B */
462: PetscCall(VecNorm(YNtmp, NORM_2, &c0)); /* this could be improved */
463: c0 = PetscSqr(c0);
464: PetscCall(VecDotRealPart(YCtmp, YNtmp, &c1));
465: c1 = 2.0 * c1;
466: PetscCall(VecNorm(YCtmp, NORM_2, &c2)); /* this could be improved */
467: c2 = PetscSqr(c2) - PetscSqr(delta);
468: tau_pos = (c1 + PetscSqrtReal(PetscSqr(c1) - 4. * c0 * c2)) / (2. * c0); /* quadratic formula */
469: tau_neg = (c1 - PetscSqrtReal(PetscSqr(c1) - 4. * c0 * c2)) / (2. * c0);
470: tau = PetscMax(tau_pos, tau_neg); /* can tau_neg > tau_pos? I don't think so, but just in case. */
471: PetscCall(PetscInfo(snes, "DL evaluated. tau: %8.4e, ynnorm: %8.4e, ycnorm: %8.4e\n", (double)tau, (double)ynnorm, (double)ycnorm));
472: PetscCall(VecWAXPY(W, tau, YNtmp, YCtmp));
473: PetscCall(VecAXPY(W, -tau, YCtmp));
474: PetscCall(VecCopy(W, Y)); /* this could be improved */
475: }
476: } else {
477: /* if Cauchy is disabled, only use Newton direction */
478: auk = delta / ynnorm;
479: PetscCall(VecScale(YNtmp, auk));
480: PetscCall(VecCopy(YNtmp, Y)); /* this could be improved (many VecCopy, VecNorm)*/
481: }
483: PetscCall(VecNorm(Y, NORM_2, &ynorm)); /* compute the final ynorm */
484: f0 = 0.5 * PetscSqr(fnorm); /* minimizing function f(X) */
485: PetscCall(MatMult(jac, Y, W));
486: PetscCall(VecDotRealPart(W, W, &yTHy)); /* completes GradY^T J^T J GradY */
487: PetscCall(VecDotRealPart(GradF, Y, &gTy));
488: mp = f0 - gTy + 0.5 * yTHy; /* quadratic model to satisfy, -gTy because our update is X-Y*/
490: /* scale back solution update */
491: if (bs > 1 && neP->auto_scale_multiphase) {
492: for (j = 0; j < bs; j++) {
493: PetscCall(VecStrideScale(Y, j, inorms[j]));
494: if (inner_count == 0) {
495: /* TRDC inner algorithm does not need scaled X after calculating delta in the outer iteration */
496: /* need to scale back X to match Y and provide proper update to the external code */
497: PetscCall(VecStrideScale(X, j, inorms[j]));
498: }
499: }
500: if (inner_count == 0) PetscCall(VecNorm(X, NORM_2, &temp_xnorm)); /* only in the first iteration */
501: PetscCall(VecNorm(Y, NORM_2, &temp_ynorm));
502: } else {
503: temp_xnorm = xnorm;
504: temp_ynorm = ynorm;
505: }
506: inner_count++;
508: /* Evaluate the solution to meet the improvement ratio criteria */
509: PetscCall(SNESNewtonTRDCPreCheck(snes, X, Y, &changed_y));
510: PetscCall(VecWAXPY(W, -1.0, Y, X));
511: PetscCall(SNESNewtonTRDCPostCheck(snes, X, Y, W, &changed_y, &changed_w));
512: if (changed_y) PetscCall(VecWAXPY(W, -1.0, Y, X));
513: PetscCall(VecCopy(Y, snes->vec_sol_update));
514: PetscCall(SNESComputeFunction(snes, W, G)); /* F(X-Y) = G */
515: PetscCall(VecNorm(G, NORM_2, &gnorm)); /* gnorm <- || g || */
516: SNESCheckFunctionDomainError(snes, gnorm);
517: g = 0.5 * PetscSqr(gnorm); /* minimizing function g(W) */
518: if (f0 == mp) rho = 0.0;
519: else rho = (f0 - g) / (f0 - mp); /* actual improvement over predicted improvement */
521: if (rho < neP->eta2) {
522: delta *= neP->t1; /* shrink the region */
523: } else if (rho > neP->eta3) {
524: delta = PetscMin(neP->t2 * delta, deltaM); /* expand the region, but not greater than deltaM */
525: }
527: neP->delta = delta;
528: if (rho >= neP->eta1) {
529: /* unscale delta and xnorm before going to the next outer iteration */
530: if (bs > 1 && neP->auto_scale_multiphase) {
531: neP->delta = delta / xnorm;
532: xnorm = temp_xnorm;
533: ynorm = temp_ynorm;
534: }
535: neP->rho_satisfied = PETSC_TRUE;
536: break; /* the improvement ratio is satisfactory */
537: }
538: PetscCall(PetscInfo(snes, "Trying again in smaller region\n"));
540: /* check to see if progress is hopeless */
541: neP->itflag = PETSC_FALSE;
542: /* both delta, ynorm, and xnorm are either scaled or unscaled */
543: PetscCall(SNESTRDC_Converged_Private(snes, snes->iter, xnorm, ynorm, fnorm, &reason, snes->cnvP));
544: /* if multiphase state changes, break out inner iteration */
545: if (reason == SNES_BREAKOUT_INNER_ITER) {
546: if (bs > 1 && neP->auto_scale_multiphase) {
547: /* unscale delta and xnorm before going to the next outer iteration */
548: neP->delta = delta / xnorm;
549: xnorm = temp_xnorm;
550: ynorm = temp_ynorm;
551: }
552: reason = SNES_CONVERGED_ITERATING;
553: break;
554: }
555: if (reason == SNES_CONVERGED_SNORM_RELATIVE) reason = SNES_DIVERGED_INNER;
556: if (reason) {
557: if (reason < 0) {
558: /* We're not progressing, so return with the current iterate */
559: PetscCall(SNESMonitor(snes, i + 1, fnorm));
560: breakout = PETSC_TRUE;
561: break;
562: } else if (reason > 0) {
563: /* We're converged, so return with the current iterate and update solution */
564: PetscCall(SNESMonitor(snes, i + 1, fnorm));
565: breakout = PETSC_FALSE;
566: break;
567: }
568: }
569: snes->numFailures++;
570: }
571: if (!breakout) {
572: /* Update function and solution vectors */
573: fnorm = gnorm;
574: PetscCall(VecCopy(G, F));
575: PetscCall(VecCopy(W, X));
576: /* Monitor convergence */
577: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
578: snes->iter = i + 1;
579: snes->norm = fnorm;
580: snes->xnorm = xnorm;
581: snes->ynorm = ynorm;
582: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
583: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits));
584: PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
585: /* Test for convergence, xnorm = || X || */
586: neP->itflag = PETSC_TRUE;
587: if (snes->ops->converged != SNESConvergedSkip) PetscCall(VecNorm(X, NORM_2, &xnorm));
588: PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, fnorm, &reason, snes->cnvP);
589: if (reason) break;
590: } else break;
591: }
593: /* PetscCall(PetscFree(inorms)); */
594: if (i == maxits) {
595: PetscCall(PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", maxits));
596: if (!reason) reason = SNES_DIVERGED_MAX_IT;
597: }
598: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
599: snes->reason = reason;
600: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
601: if (convtest != SNESTRDC_KSPConverged_Private) {
602: PetscCall(KSPGetAndClearConvergenceTest(ksp, &ctx->convtest, &ctx->convctx, &ctx->convdestroy));
603: PetscCall(PetscFree(ctx));
604: PetscCall(KSPSetConvergenceTest(ksp, convtest, convctx, convdestroy));
605: }
606: PetscFunctionReturn(PETSC_SUCCESS);
607: }
609: static PetscErrorCode SNESSetUp_NEWTONTRDC(SNES snes)
610: {
611: PetscFunctionBegin;
612: PetscCall(SNESSetWorkVecs(snes, 6));
613: PetscCall(SNESSetUpMatrices(snes));
614: PetscFunctionReturn(PETSC_SUCCESS);
615: }
617: static PetscErrorCode SNESDestroy_NEWTONTRDC(SNES snes)
618: {
619: PetscFunctionBegin;
620: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonTRSetTolerances_C", NULL));
621: PetscCall(PetscFree(snes->data));
622: PetscFunctionReturn(PETSC_SUCCESS);
623: }
625: static PetscErrorCode SNESSetFromOptions_NEWTONTRDC(SNES snes, PetscOptionItems PetscOptionsObject)
626: {
627: SNES_NEWTONTRDC *ctx = (SNES_NEWTONTRDC *)snes->data;
629: PetscFunctionBegin;
630: PetscOptionsHeadBegin(PetscOptionsObject, "SNES trust region options for nonlinear equations");
631: PetscCall(PetscOptionsReal("-snes_trdc_tol", "Trust region tolerance", "SNESNewtonTRSetTolerances", ctx->deltatol, &ctx->deltatol, NULL));
632: PetscCall(PetscOptionsReal("-snes_trdc_eta1", "eta1", "None", ctx->eta1, &ctx->eta1, NULL));
633: PetscCall(PetscOptionsReal("-snes_trdc_eta2", "eta2", "None", ctx->eta2, &ctx->eta2, NULL));
634: PetscCall(PetscOptionsReal("-snes_trdc_eta3", "eta3", "None", ctx->eta3, &ctx->eta3, NULL));
635: PetscCall(PetscOptionsReal("-snes_trdc_t1", "t1", "None", ctx->t1, &ctx->t1, NULL));
636: PetscCall(PetscOptionsReal("-snes_trdc_t2", "t2", "None", ctx->t2, &ctx->t2, NULL));
637: PetscCall(PetscOptionsReal("-snes_trdc_deltaM", "deltaM", "None", ctx->deltaM, &ctx->deltaM, NULL));
638: PetscCall(PetscOptionsReal("-snes_trdc_delta0", "delta0", "None", ctx->delta0, &ctx->delta0, NULL));
639: PetscCall(PetscOptionsReal("-snes_trdc_auto_scale_max", "auto_scale_max", "None", ctx->auto_scale_max, &ctx->auto_scale_max, NULL));
640: PetscCall(PetscOptionsBool("-snes_trdc_use_cauchy", "use_cauchy", "use Cauchy step and direction", ctx->use_cauchy, &ctx->use_cauchy, NULL));
641: PetscCall(PetscOptionsBool("-snes_trdc_auto_scale_multiphase", "auto_scale_multiphase", "Auto scaling for proper cauchy direction", ctx->auto_scale_multiphase, &ctx->auto_scale_multiphase, NULL));
642: PetscOptionsHeadEnd();
643: PetscFunctionReturn(PETSC_SUCCESS);
644: }
646: static PetscErrorCode SNESView_NEWTONTRDC(SNES snes, PetscViewer viewer)
647: {
648: SNES_NEWTONTRDC *tr = (SNES_NEWTONTRDC *)snes->data;
649: PetscBool isascii;
651: PetscFunctionBegin;
652: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
653: if (isascii) {
654: PetscCall(PetscViewerASCIIPrintf(viewer, " Trust region tolerance %g\n", (double)tr->deltatol));
655: PetscCall(PetscViewerASCIIPrintf(viewer, " eta1=%g, eta2=%g, eta3=%g\n", (double)tr->eta1, (double)tr->eta2, (double)tr->eta3));
656: PetscCall(PetscViewerASCIIPrintf(viewer, " delta0=%g, t1=%g, t2=%g, deltaM=%g\n", (double)tr->delta0, (double)tr->t1, (double)tr->t2, (double)tr->deltaM));
657: }
658: PetscFunctionReturn(PETSC_SUCCESS);
659: }
661: /*MC
662: SNESNEWTONTRDC - Newton based nonlinear solver that uses trust-region dogleg method with Cauchy direction
664: Options Database Keys:
665: + -snes_trdc_tol <tol> - trust region tolerance
666: . -snes_trdc_eta1 <eta1> - trust region parameter 0.0 <= eta1 <= eta2, rho >= eta1 breaks out of the inner iteration (default: eta1=0.001)
667: . -snes_trdc_eta2 <eta2> - trust region parameter 0.0 <= eta1 <= eta2, rho <= eta2 shrinks the trust region (default: eta2=0.25)
668: . -snes_trdc_eta3 <eta3> - trust region parameter eta3 > eta2, rho >= eta3 expands the trust region (default: eta3=0.75)
669: . -snes_trdc_t1 <t1> - trust region parameter, shrinking factor of trust region (default: 0.25)
670: . -snes_trdc_t2 <t2> - trust region parameter, expanding factor of trust region (default: 2.0)
671: . -snes_trdc_deltaM <deltaM> - trust region parameter, max size of trust region, $deltaM*norm2(x)$ (default: 0.5)
672: . -snes_trdc_delta0 <delta0> - trust region parameter, initial size of trust region, $delta0*norm2(x)$ (default: 0.1)
673: . -snes_trdc_auto_scale_max <auto_scale_max> - used with auto_scale_multiphase, caps the maximum auto-scaling factor
674: . -snes_trdc_use_cauchy <use_cauchy> - True uses dogleg Cauchy (Steepest Descent direction) step & direction in the trust region algorithm
675: - -snes_trdc_auto_scale_multiphase <auto_scale_multiphase> - True turns on auto-scaling for multivariable block matrix for Cauchy and trust region
677: Level: advanced
679: Notes:
680: `SNESNEWTONTRDC` only works for root-finding problems and does not support objective functions.
681: The main difference with respect to `SNESNEWTONTR` is that `SNESNEWTONTRDC` scales the trust region by the norm of the current linearization point.
682: Future version may extend the `SNESNEWTONTR` code and deprecate `SNESNEWTONTRDC`.
684: For details, see {cite}`park2021linear`
686: .seealso: [](ch_snes), `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESNewtonTRSetTolerances()`,
687: `SNESNewtonTRDCPreCheck()`, `SNESNewtonTRDCGetPreCheck()`, `SNESNewtonTRDCSetPostCheck()`, `SNESNewtonTRDCGetPostCheck()`,
688: `SNESNewtonTRDCGetRhoFlag()`, `SNESNewtonTRDCSetPreCheck()`
689: M*/
690: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONTRDC(SNES snes)
691: {
692: SNES_NEWTONTRDC *neP;
694: PetscFunctionBegin;
695: snes->ops->setup = SNESSetUp_NEWTONTRDC;
696: snes->ops->solve = SNESSolve_NEWTONTRDC;
697: snes->ops->destroy = SNESDestroy_NEWTONTRDC;
698: snes->ops->setfromoptions = SNESSetFromOptions_NEWTONTRDC;
699: snes->ops->view = SNESView_NEWTONTRDC;
701: snes->usesksp = PETSC_TRUE;
702: snes->usesnpc = PETSC_FALSE;
704: snes->alwayscomputesfinalresidual = PETSC_TRUE;
706: PetscCall(SNESParametersInitialize(snes));
708: PetscCall(PetscNew(&neP));
709: snes->data = (void *)neP;
710: neP->eta1 = 0.001;
711: neP->eta2 = 0.25;
712: neP->eta3 = 0.75;
713: neP->t1 = 0.25;
714: neP->t2 = 2.0;
715: neP->sigma = 0.0001;
716: neP->itflag = PETSC_FALSE;
717: neP->rnorm0 = 0.0;
718: neP->ttol = 0.0;
719: neP->use_cauchy = PETSC_TRUE;
720: neP->auto_scale_multiphase = PETSC_FALSE;
721: neP->auto_scale_max = -1.0;
722: neP->rho_satisfied = PETSC_FALSE;
723: neP->delta = 0.0;
724: neP->deltaM = 0.5;
725: neP->delta0 = 0.1;
726: neP->deltatol = 1.e-12;
728: /* for multiphase (multivariable) scaling */
729: /* may be used for dynamic allocation of inorms, but it fails snes_tutorials-ex3_13
730: on test forced DIVERGED_JACOBIAN_DOMAIN test. I will use static array for now.
731: PetscCall(VecGetBlockSize(snes->work[0],&neP->bs));
732: PetscCall(PetscCalloc1(neP->bs,&neP->inorms));
733: */
734: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNewtonTRSetTolerances_C", SNESNewtonTRSetTolerances_TRDC));
735: PetscFunctionReturn(PETSC_SUCCESS);
736: }