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