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:   PetscErrorCode (*convtest)(KSP, PetscInt, PetscReal, KSPConvergedReason *, void *);
 11:   PetscErrorCode (*convdestroy)(void *);
 12:   void *convctx;
 13: } SNES_TRDC_KSPConverged_Ctx;

 15: static PetscErrorCode SNESTRDC_KSPConverged_Private(KSP ksp, PetscInt n, PetscReal rnorm, KSPConvergedReason *reason, void *cctx)
 16: {
 17:   SNES_TRDC_KSPConverged_Ctx *ctx  = (SNES_TRDC_KSPConverged_Ctx *)cctx;
 18:   SNES                        snes = ctx->snes;
 19:   SNES_NEWTONTRDC            *neP  = (SNES_NEWTONTRDC *)snes->data;
 20:   Vec                         x;
 21:   PetscReal                   nrm;

 23:   PetscFunctionBegin;
 24:   PetscCall((*ctx->convtest)(ksp, n, rnorm, reason, ctx->convctx));
 25:   if (*reason) PetscCall(PetscInfo(snes, "Default or user provided convergence test KSP iterations=%" PetscInt_FMT ", rnorm=%g\n", n, (double)rnorm));
 26:   /* Determine norm of solution */
 27:   PetscCall(KSPBuildSolution(ksp, NULL, &x));
 28:   PetscCall(VecNorm(x, NORM_2, &nrm));
 29:   if (nrm >= neP->delta) {
 30:     PetscCall(PetscInfo(snes, "Ending linear iteration early, delta=%g, length=%g\n", (double)neP->delta, (double)nrm));
 31:     *reason = KSP_CONVERGED_STEP_LENGTH;
 32:   }
 33:   PetscFunctionReturn(PETSC_SUCCESS);
 34: }

 36: static PetscErrorCode SNESTRDC_KSPConverged_Destroy(void *cctx)
 37: {
 38:   SNES_TRDC_KSPConverged_Ctx *ctx = (SNES_TRDC_KSPConverged_Ctx *)cctx;

 40:   PetscFunctionBegin;
 41:   PetscCall((*ctx->convdestroy)(ctx->convctx));
 42:   PetscCall(PetscFree(ctx));

 44:   PetscFunctionReturn(PETSC_SUCCESS);
 45: }

 47: /*
 48:    SNESTRDC_Converged_Private -test convergence JUST for
 49:    the trust region tolerance.

 51: */
 52: static PetscErrorCode SNESTRDC_Converged_Private(SNES snes, PetscInt it, PetscReal xnorm, PetscReal pnorm, PetscReal fnorm, SNESConvergedReason *reason, void *dummy)
 53: {
 54:   SNES_NEWTONTRDC *neP = (SNES_NEWTONTRDC *)snes->data;

 56:   PetscFunctionBegin;
 57:   *reason = SNES_CONVERGED_ITERATING;
 58:   if (neP->delta < xnorm * snes->deltatol) {
 59:     PetscCall(PetscInfo(snes, "Diverged due to too small a trust region %g<%g*%g\n", (double)neP->delta, (double)xnorm, (double)snes->deltatol));
 60:     *reason = SNES_DIVERGED_TR_DELTA;
 61:   } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
 62:     PetscCall(PetscInfo(snes, "Exceeded maximum number of function evaluations: %" PetscInt_FMT "\n", snes->max_funcs));
 63:     *reason = SNES_DIVERGED_FUNCTION_COUNT;
 64:   }
 65:   PetscFunctionReturn(PETSC_SUCCESS);
 66: }

 68: /*@
 69:   SNESNewtonTRDCGetRhoFlag - Get whether the current solution update is within the trust-region.

 71:   Logically Collective

 73:   Input Parameter:
 74: . snes - the nonlinear solver object

 76:   Output Parameter:
 77: . rho_flag - `PETSC_FALSE` or `PETSC_TRUE`

 79:   Level: developer

 81: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTRDC`, `SNESNewtonTRDCPreCheck()`, `SNESNewtonTRDCGetPreCheck()`, `SNESNewtonTRDCSetPreCheck()`,
 82:           `SNESNewtonTRDCSetPostCheck()`, `SNESNewtonTRDCGetPostCheck()`
 83: @*/
 84: PetscErrorCode SNESNewtonTRDCGetRhoFlag(SNES snes, PetscBool *rho_flag)
 85: {
 86:   SNES_NEWTONTRDC *tr = (SNES_NEWTONTRDC *)snes->data;

 88:   PetscFunctionBegin;
 90:   PetscAssertPointer(rho_flag, 2);
 91:   *rho_flag = tr->rho_satisfied;
 92:   PetscFunctionReturn(PETSC_SUCCESS);
 93: }

 95: /*@C
 96:   SNESNewtonTRDCSetPreCheck - Sets a user function that is called before the search step has been determined.
 97:   Allows the user a chance to change or override the trust region decision.

 99:   Logically Collective

101:   Input Parameters:
102: + snes - the nonlinear solver object
103: . func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRDCPreCheck()`
104: - ctx  - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)

106:   Level: intermediate

108:   Note:
109:   This function is called BEFORE the function evaluation within the `SNESNEWTONTRDC` solver.

111: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTRDC`, `SNESNewtonTRDCPreCheck()`, `SNESNewtonTRDCGetPreCheck()`, `SNESNewtonTRDCSetPostCheck()`, `SNESNewtonTRDCGetPostCheck()`,
112:           `SNESNewtonTRDCGetRhoFlag()`
113: @*/
114: PetscErrorCode SNESNewtonTRDCSetPreCheck(SNES snes, PetscErrorCode (*func)(SNES, Vec, Vec, PetscBool *, void *), void *ctx)
115: {
116:   SNES_NEWTONTRDC *tr = (SNES_NEWTONTRDC *)snes->data;

118:   PetscFunctionBegin;
120:   if (func) tr->precheck = func;
121:   if (ctx) tr->precheckctx = ctx;
122:   PetscFunctionReturn(PETSC_SUCCESS);
123: }

125: /*@C
126:   SNESNewtonTRDCGetPreCheck - Gets the pre-check function optionally set with `SNESNewtonTRDCSetPreCheck()`

128:   Not Collective

130:   Input Parameter:
131: . snes - the nonlinear solver context

133:   Output Parameters:
134: + func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRDCPreCheck()`
135: - ctx  - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)

137:   Level: intermediate

139: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTRDC`, `SNESNewtonTRDCSetPreCheck()`, `SNESNewtonTRDCPreCheck()`
140: @*/
141: PetscErrorCode SNESNewtonTRDCGetPreCheck(SNES snes, PetscErrorCode (**func)(SNES, Vec, Vec, PetscBool *, void *), void **ctx)
142: {
143:   SNES_NEWTONTRDC *tr = (SNES_NEWTONTRDC *)snes->data;

145:   PetscFunctionBegin;
147:   if (func) *func = tr->precheck;
148:   if (ctx) *ctx = tr->precheckctx;
149:   PetscFunctionReturn(PETSC_SUCCESS);
150: }

152: /*@C
153:   SNESNewtonTRDCSetPostCheck - Sets a user function that is called after the search step has been determined but before the next
154:   function evaluation. Allows the user a chance to change or override the decision of the line search routine

156:   Logically Collective

158:   Input Parameters:
159: + snes - the nonlinear solver object
160: . func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRDCPostCheck()`
161: - ctx  - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)

163:   Level: intermediate

165:   Note:
166:   This function is called BEFORE the function evaluation within the `SNESNEWTONTRDC` solver while the function set in
167:   `SNESLineSearchSetPostCheck()` is called AFTER the function evaluation.

169: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTRDC`, `SNESNewtonTRDCPostCheck()`, `SNESNewtonTRDCGetPostCheck()`, `SNESNewtonTRDCSetPreCheck()`, `SNESNewtonTRDCGetPreCheck()`
170: @*/
171: PetscErrorCode SNESNewtonTRDCSetPostCheck(SNES snes, PetscErrorCode (*func)(SNES, Vec, Vec, Vec, PetscBool *, PetscBool *, void *), void *ctx)
172: {
173:   SNES_NEWTONTRDC *tr = (SNES_NEWTONTRDC *)snes->data;

175:   PetscFunctionBegin;
177:   if (func) tr->postcheck = func;
178:   if (ctx) tr->postcheckctx = ctx;
179:   PetscFunctionReturn(PETSC_SUCCESS);
180: }

182: /*@C
183:   SNESNewtonTRDCGetPostCheck - Gets the post-check function optionally set with `SNESNewtonTRDCSetPostCheck()`

185:   Not Collective

187:   Input Parameter:
188: . snes - the nonlinear solver context

190:   Output Parameters:
191: + func - [optional] function evaluation routine, for the calling sequence see `SNESNewtonTRDCPostCheck()`
192: - ctx  - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)

194:   Level: intermediate

196: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTRDC`, `SNESNewtonTRDCSetPostCheck()`, `SNESNewtonTRDCPostCheck()`, `SNESNewtonTRDCSetPreCheck()`, `SNESNewtonTRDCGetPreCheck()`
197: @*/
198: PetscErrorCode SNESNewtonTRDCGetPostCheck(SNES snes, PetscErrorCode (**func)(SNES, Vec, Vec, Vec, PetscBool *, PetscBool *, void *), void **ctx)
199: {
200:   SNES_NEWTONTRDC *tr = (SNES_NEWTONTRDC *)snes->data;

202:   PetscFunctionBegin;
204:   if (func) *func = tr->postcheck;
205:   if (ctx) *ctx = tr->postcheckctx;
206:   PetscFunctionReturn(PETSC_SUCCESS);
207: }

209: // PetscClangLinter pragma disable: -fdoc-internal-linkage
210: /*@C
211:    SNESNewtonTRDCPreCheck - Called before the step has been determined in `SNESNEWTONTRDC`

213:    Logically Collective

215:    Input Parameters:
216: +  snes - the solver
217: .  X - The last solution
218: -  Y - The step direction

220:    Output Parameter:
221: .  changed_Y - Indicator that the step direction `Y` has been changed.

223:    Level: developer

225: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTRDC`, `SNESNewtonTRDCSetPreCheck()`, `SNESNewtonTRDCGetPreCheck()`, `SNESNewtonTRDCPostCheck()`
226: @*/
227: static PetscErrorCode SNESNewtonTRDCPreCheck(SNES snes, Vec X, Vec Y, PetscBool *changed_Y)
228: {
229:   SNES_NEWTONTRDC *tr = (SNES_NEWTONTRDC *)snes->data;

231:   PetscFunctionBegin;
232:   *changed_Y = PETSC_FALSE;
233:   if (tr->precheck) {
234:     PetscCall((*tr->precheck)(snes, X, Y, changed_Y, tr->precheckctx));
236:   }
237:   PetscFunctionReturn(PETSC_SUCCESS);
238: }

240: // PetscClangLinter pragma disable: -fdoc-internal-linkage
241: /*@C
242:    SNESNewtonTRDCPostCheck - Called after the step has been determined in `SNESNEWTONTRDC` but before the function evaluation at that step

244:    Logically Collective

246:    Input Parameters:
247: +  snes - the solver
248: .  X - The last solution
249: .  Y - The full step direction
250: -  W - The updated solution, W = X - Y

252:    Output Parameters:
253: +  changed_Y - indicator if step has been changed
254: -  changed_W - Indicator if the new candidate solution `W` has been changed.

256:    Level: developer

258:    Note:
259:      If `Y` is changed then `W` is recomputed as `X` - `Y`

261: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTRDC`, `SNESNEWTONTRDC`, `SNESNewtonTRDCSetPostCheck()`, `SNESNewtonTRDCGetPostCheck()`, `SNESNewtonTRDCPreCheck()
262: @*/
263: static PetscErrorCode SNESNewtonTRDCPostCheck(SNES snes, Vec X, Vec Y, Vec W, PetscBool *changed_Y, PetscBool *changed_W)
264: {
265:   SNES_NEWTONTRDC *tr = (SNES_NEWTONTRDC *)snes->data;

267:   PetscFunctionBegin;
268:   *changed_Y = PETSC_FALSE;
269:   *changed_W = PETSC_FALSE;
270:   if (tr->postcheck) {
271:     PetscCall((*tr->postcheck)(snes, X, Y, W, changed_Y, changed_W, tr->postcheckctx));
274:   }
275:   PetscFunctionReturn(PETSC_SUCCESS);
276: }

278: /*
279:    SNESSolve_NEWTONTRDC - Implements Newton's Method with trust-region subproblem and adds dogleg Cauchy
280:    (Steepest Descent direction) step and direction if the trust region is not satisfied for solving system of
281:    nonlinear equations

283: */
284: static PetscErrorCode SNESSolve_NEWTONTRDC(SNES snes)
285: {
286:   SNES_NEWTONTRDC            *neP = (SNES_NEWTONTRDC *)snes->data;
287:   Vec                         X, F, Y, G, W, GradF, YNtmp;
288:   Vec                         YCtmp;
289:   Mat                         jac;
290:   PetscInt                    maxits, i, j, lits, inner_count, bs;
291:   PetscReal                   rho, fnorm, gnorm, xnorm = 0, delta, ynorm, temp_xnorm, temp_ynorm; /* TRDC inner iteration */
292:   PetscReal                   inorms[99];                                                         /* need to make it dynamic eventually, fixed max block size of 99 for now */
293:   PetscReal                   deltaM, ynnorm, f0, mp, gTy, g, yTHy;                               /* rho calculation */
294:   PetscReal                   auk, gfnorm, ycnorm, c0, c1, c2, tau, tau_pos, tau_neg, gTBg;       /* Cauchy Point */
295:   KSP                         ksp;
296:   SNESConvergedReason         reason   = SNES_CONVERGED_ITERATING;
297:   PetscBool                   breakout = PETSC_FALSE;
298:   SNES_TRDC_KSPConverged_Ctx *ctx;
299:   PetscErrorCode (*convtest)(KSP, PetscInt, PetscReal, KSPConvergedReason *, void *), (*convdestroy)(void *);
300:   void *convctx;

302:   PetscFunctionBegin;
303:   maxits = snes->max_its;  /* maximum number of iterations */
304:   X      = snes->vec_sol;  /* solution vector */
305:   F      = snes->vec_func; /* residual vector */
306:   Y      = snes->work[0];  /* update vector */
307:   G      = snes->work[1];  /* updated residual */
308:   W      = snes->work[2];  /* temporary vector */
309:   GradF  = snes->work[3];  /* grad f = J^T F */
310:   YNtmp  = snes->work[4];  /* Newton solution */
311:   YCtmp  = snes->work[5];  /* Cauchy solution */

313:   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);

315:   PetscCall(VecGetBlockSize(YNtmp, &bs));

317:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
318:   snes->iter = 0;
319:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));

321:   /* Set the linear stopping criteria to use the More' trick. From tr.c */
322:   PetscCall(SNESGetKSP(snes, &ksp));
323:   PetscCall(KSPGetConvergenceTest(ksp, &convtest, &convctx, &convdestroy));
324:   if (convtest != SNESTRDC_KSPConverged_Private) {
325:     PetscCall(PetscNew(&ctx));
326:     ctx->snes = snes;
327:     PetscCall(KSPGetAndClearConvergenceTest(ksp, &ctx->convtest, &ctx->convctx, &ctx->convdestroy));
328:     PetscCall(KSPSetConvergenceTest(ksp, SNESTRDC_KSPConverged_Private, ctx, SNESTRDC_KSPConverged_Destroy));
329:     PetscCall(PetscInfo(snes, "Using Krylov convergence test SNESTRDC_KSPConverged_Private\n"));
330:   }

332:   if (!snes->vec_func_init_set) {
333:     PetscCall(SNESComputeFunction(snes, X, F)); /* F(X) */
334:   } else snes->vec_func_init_set = PETSC_FALSE;

336:   PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- || F || */
337:   SNESCheckFunctionNorm(snes, fnorm);
338:   PetscCall(VecNorm(X, NORM_2, &xnorm)); /* xnorm <- || X || */

340:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
341:   snes->norm = fnorm;
342:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
343:   delta      = xnorm ? neP->delta0 * xnorm : neP->delta0; /* initial trust region size scaled by xnorm */
344:   deltaM     = xnorm ? neP->deltaM * xnorm : neP->deltaM; /* maximum trust region size scaled by xnorm */
345:   neP->delta = delta;
346:   PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));
347:   PetscCall(SNESMonitor(snes, 0, fnorm));

349:   neP->rho_satisfied = PETSC_FALSE;

351:   /* test convergence */
352:   PetscUseTypeMethod(snes, converged, snes->iter, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
353:   if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);

355:   for (i = 0; i < maxits; i++) {
356:     PetscBool changed_y;
357:     PetscBool changed_w;

359:     /* dogleg method */
360:     PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
361:     SNESCheckJacobianDomainerror(snes);
362:     PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian));
363:     PetscCall(KSPSolve(snes->ksp, F, YNtmp)); /* Quasi Newton Solution */
364:     SNESCheckKSPSolve(snes);                  /* this is necessary but old tr.c did not have it*/
365:     PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
366:     PetscCall(SNESGetJacobian(snes, &jac, NULL, NULL, NULL));

368:     /* rescale Jacobian, Newton solution update, and re-calculate delta for multiphase (multivariable)
369:        for inner iteration and Cauchy direction calculation
370:     */
371:     if (bs > 1 && neP->auto_scale_multiphase) {
372:       PetscCall(VecStrideNormAll(YNtmp, NORM_INFINITY, inorms));
373:       for (j = 0; j < bs; j++) {
374:         if (neP->auto_scale_max > 1.0) {
375:           if (inorms[j] < 1.0 / neP->auto_scale_max) inorms[j] = 1.0 / neP->auto_scale_max;
376:         }
377:         PetscCall(VecStrideSet(W, j, inorms[j]));
378:         PetscCall(VecStrideScale(YNtmp, j, 1.0 / inorms[j]));
379:         PetscCall(VecStrideScale(X, j, 1.0 / inorms[j]));
380:       }
381:       PetscCall(VecNorm(X, NORM_2, &xnorm));
382:       if (i == 0) {
383:         delta = neP->delta0 * xnorm;
384:       } else {
385:         delta = neP->delta * xnorm;
386:       }
387:       deltaM = neP->deltaM * xnorm;
388:       PetscCall(MatDiagonalScale(jac, NULL, W));
389:     }

391:     /* calculating GradF of minimization function */
392:     PetscCall(MatMultTranspose(jac, F, GradF)); /* grad f = J^T F */
393:     PetscCall(VecNorm(YNtmp, NORM_2, &ynnorm)); /* ynnorm <- || Y_newton || */

395:     inner_count        = 0;
396:     neP->rho_satisfied = PETSC_FALSE;
397:     while (1) {
398:       if (ynnorm <= delta) { /* see if the Newton solution is within the trust region */
399:         PetscCall(VecCopy(YNtmp, Y));
400:       } else if (neP->use_cauchy) { /* use Cauchy direction if enabled */
401:         PetscCall(MatMult(jac, GradF, W));
402:         PetscCall(VecDotRealPart(W, W, &gTBg));     /* completes GradF^T J^T J GradF */
403:         PetscCall(VecNorm(GradF, NORM_2, &gfnorm)); /* grad f norm <- || grad f || */
404:         if (gTBg <= 0.0) {
405:           auk = PETSC_MAX_REAL;
406:         } else {
407:           auk = PetscSqr(gfnorm) / gTBg;
408:         }
409:         auk = PetscMin(delta / gfnorm, auk);
410:         PetscCall(VecCopy(GradF, YCtmp));           /* this could be improved */
411:         PetscCall(VecScale(YCtmp, auk));            /* YCtmp, Cauchy solution*/
412:         PetscCall(VecNorm(YCtmp, NORM_2, &ycnorm)); /* ycnorm <- || Y_cauchy || */
413:         if (ycnorm >= delta) {                      /* see if the Cauchy solution meets the criteria */
414:           PetscCall(VecCopy(YCtmp, Y));
415:           PetscCall(PetscInfo(snes, "DL evaluated. delta: %8.4e, ynnorm: %8.4e, ycnorm: %8.4e\n", (double)delta, (double)ynnorm, (double)ycnorm));
416:         } else {                                  /* take ratio, tau, of Cauchy and Newton direction and step */
417:           PetscCall(VecAXPY(YNtmp, -1.0, YCtmp)); /* YCtmp = A, YNtmp = B */
418:           PetscCall(VecNorm(YNtmp, NORM_2, &c0)); /* this could be improved */
419:           c0 = PetscSqr(c0);
420:           PetscCall(VecDotRealPart(YCtmp, YNtmp, &c1));
421:           c1 = 2.0 * c1;
422:           PetscCall(VecNorm(YCtmp, NORM_2, &c2)); /* this could be improved */
423:           c2      = PetscSqr(c2) - PetscSqr(delta);
424:           tau_pos = (c1 + PetscSqrtReal(PetscSqr(c1) - 4. * c0 * c2)) / (2. * c0); /* quadratic formula */
425:           tau_neg = (c1 - PetscSqrtReal(PetscSqr(c1) - 4. * c0 * c2)) / (2. * c0);
426:           tau     = PetscMax(tau_pos, tau_neg); /* can tau_neg > tau_pos? I don't think so, but just in case. */
427:           PetscCall(PetscInfo(snes, "DL evaluated. tau: %8.4e, ynnorm: %8.4e, ycnorm: %8.4e\n", (double)tau, (double)ynnorm, (double)ycnorm));
428:           PetscCall(VecWAXPY(W, tau, YNtmp, YCtmp));
429:           PetscCall(VecAXPY(W, -tau, YCtmp));
430:           PetscCall(VecCopy(W, Y)); /* this could be improved */
431:         }
432:       } else {
433:         /* if Cauchy is disabled, only use Newton direction */
434:         auk = delta / ynnorm;
435:         PetscCall(VecScale(YNtmp, auk));
436:         PetscCall(VecCopy(YNtmp, Y)); /* this could be improved (many VecCopy, VecNorm)*/
437:       }

439:       PetscCall(VecNorm(Y, NORM_2, &ynorm)); /* compute the final ynorm  */
440:       f0 = 0.5 * PetscSqr(fnorm);            /* minimizing function f(X) */
441:       PetscCall(MatMult(jac, Y, W));
442:       PetscCall(VecDotRealPart(W, W, &yTHy)); /* completes GradY^T J^T J GradY */
443:       PetscCall(VecDotRealPart(GradF, Y, &gTy));
444:       mp = f0 - gTy + 0.5 * yTHy; /* quadratic model to satisfy, -gTy because our update is X-Y*/

446:       /* scale back solution update */
447:       if (bs > 1 && neP->auto_scale_multiphase) {
448:         for (j = 0; j < bs; j++) {
449:           PetscCall(VecStrideScale(Y, j, inorms[j]));
450:           if (inner_count == 0) {
451:             /* TRDC inner algorithm does not need scaled X after calculating delta in the outer iteration */
452:             /* need to scale back X to match Y and provide proper update to the external code */
453:             PetscCall(VecStrideScale(X, j, inorms[j]));
454:           }
455:         }
456:         if (inner_count == 0) PetscCall(VecNorm(X, NORM_2, &temp_xnorm)); /* only in the first iteration */
457:         PetscCall(VecNorm(Y, NORM_2, &temp_ynorm));
458:       } else {
459:         temp_xnorm = xnorm;
460:         temp_ynorm = ynorm;
461:       }
462:       inner_count++;

464:       /* Evaluate the solution to meet the improvement ratio criteria */
465:       PetscCall(SNESNewtonTRDCPreCheck(snes, X, Y, &changed_y));
466:       PetscCall(VecWAXPY(W, -1.0, Y, X));
467:       PetscCall(SNESNewtonTRDCPostCheck(snes, X, Y, W, &changed_y, &changed_w));
468:       if (changed_y) PetscCall(VecWAXPY(W, -1.0, Y, X));
469:       PetscCall(VecCopy(Y, snes->vec_sol_update));
470:       PetscCall(SNESComputeFunction(snes, W, G)); /*  F(X-Y) = G */
471:       PetscCall(VecNorm(G, NORM_2, &gnorm));      /* gnorm <- || g || */
472:       SNESCheckFunctionNorm(snes, gnorm);
473:       g = 0.5 * PetscSqr(gnorm); /* minimizing function g(W) */
474:       if (f0 == mp) rho = 0.0;
475:       else rho = (f0 - g) / (f0 - mp); /* actual improvement over predicted improvement */

477:       if (rho < neP->eta2) {
478:         delta *= neP->t1; /* shrink the region */
479:       } else if (rho > neP->eta3) {
480:         delta = PetscMin(neP->t2 * delta, deltaM); /* expand the region, but not greater than deltaM */
481:       }

483:       neP->delta = delta;
484:       if (rho >= neP->eta1) {
485:         /* unscale delta and xnorm before going to the next outer iteration */
486:         if (bs > 1 && neP->auto_scale_multiphase) {
487:           neP->delta = delta / xnorm;
488:           xnorm      = temp_xnorm;
489:           ynorm      = temp_ynorm;
490:         }
491:         neP->rho_satisfied = PETSC_TRUE;
492:         break; /* the improvement ratio is satisfactory */
493:       }
494:       PetscCall(PetscInfo(snes, "Trying again in smaller region\n"));

496:       /* check to see if progress is hopeless */
497:       neP->itflag = PETSC_FALSE;
498:       /* both delta, ynorm, and xnorm are either scaled or unscaled */
499:       PetscCall(SNESTRDC_Converged_Private(snes, snes->iter, xnorm, ynorm, fnorm, &reason, snes->cnvP));
500:       /* if multiphase state changes, break out inner iteration */
501:       if (reason == SNES_BREAKOUT_INNER_ITER) {
502:         if (bs > 1 && neP->auto_scale_multiphase) {
503:           /* unscale delta and xnorm before going to the next outer iteration */
504:           neP->delta = delta / xnorm;
505:           xnorm      = temp_xnorm;
506:           ynorm      = temp_ynorm;
507:         }
508:         reason = SNES_CONVERGED_ITERATING;
509:         break;
510:       }
511:       if (reason == SNES_CONVERGED_SNORM_RELATIVE) reason = SNES_DIVERGED_INNER;
512:       if (reason) {
513:         if (reason < 0) {
514:           /* We're not progressing, so return with the current iterate */
515:           PetscCall(SNESMonitor(snes, i + 1, fnorm));
516:           breakout = PETSC_TRUE;
517:           break;
518:         } else if (reason > 0) {
519:           /* We're converged, so return with the current iterate and update solution */
520:           PetscCall(SNESMonitor(snes, i + 1, fnorm));
521:           breakout = PETSC_FALSE;
522:           break;
523:         }
524:       }
525:       snes->numFailures++;
526:     }
527:     if (!breakout) {
528:       /* Update function and solution vectors */
529:       fnorm = gnorm;
530:       PetscCall(VecCopy(G, F));
531:       PetscCall(VecCopy(W, X));
532:       /* Monitor convergence */
533:       PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
534:       snes->iter  = i + 1;
535:       snes->norm  = fnorm;
536:       snes->xnorm = xnorm;
537:       snes->ynorm = ynorm;
538:       PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
539:       PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits));
540:       PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
541:       /* Test for convergence, xnorm = || X || */
542:       neP->itflag = PETSC_TRUE;
543:       if (snes->ops->converged != SNESConvergedSkip) PetscCall(VecNorm(X, NORM_2, &xnorm));
544:       PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, fnorm, &reason, snes->cnvP);
545:       if (reason) break;
546:     } else break;
547:   }

549:   /* PetscCall(PetscFree(inorms)); */
550:   if (i == maxits) {
551:     PetscCall(PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", maxits));
552:     if (!reason) reason = SNES_DIVERGED_MAX_IT;
553:   }
554:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
555:   snes->reason = reason;
556:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
557:   if (convtest != SNESTRDC_KSPConverged_Private) {
558:     PetscCall(KSPGetAndClearConvergenceTest(ksp, &ctx->convtest, &ctx->convctx, &ctx->convdestroy));
559:     PetscCall(PetscFree(ctx));
560:     PetscCall(KSPSetConvergenceTest(ksp, convtest, convctx, convdestroy));
561:   }
562:   PetscFunctionReturn(PETSC_SUCCESS);
563: }

565: static PetscErrorCode SNESSetUp_NEWTONTRDC(SNES snes)
566: {
567:   PetscFunctionBegin;
568:   PetscCall(SNESSetWorkVecs(snes, 6));
569:   PetscCall(SNESSetUpMatrices(snes));
570:   PetscFunctionReturn(PETSC_SUCCESS);
571: }

573: static PetscErrorCode SNESReset_NEWTONTRDC(SNES snes)
574: {
575:   PetscFunctionBegin;
576:   PetscFunctionReturn(PETSC_SUCCESS);
577: }

579: static PetscErrorCode SNESDestroy_NEWTONTRDC(SNES snes)
580: {
581:   PetscFunctionBegin;
582:   PetscCall(SNESReset_NEWTONTRDC(snes));
583:   PetscCall(PetscFree(snes->data));
584:   PetscFunctionReturn(PETSC_SUCCESS);
585: }

587: static PetscErrorCode SNESSetFromOptions_NEWTONTRDC(SNES snes, PetscOptionItems *PetscOptionsObject)
588: {
589:   SNES_NEWTONTRDC *ctx = (SNES_NEWTONTRDC *)snes->data;

591:   PetscFunctionBegin;
592:   PetscOptionsHeadBegin(PetscOptionsObject, "SNES trust region options for nonlinear equations");
593:   PetscCall(PetscOptionsReal("-snes_trdc_tol", "Trust region tolerance", "SNESSetTrustRegionTolerance", snes->deltatol, &snes->deltatol, NULL));
594:   PetscCall(PetscOptionsReal("-snes_trdc_eta1", "eta1", "None", ctx->eta1, &ctx->eta1, NULL));
595:   PetscCall(PetscOptionsReal("-snes_trdc_eta2", "eta2", "None", ctx->eta2, &ctx->eta2, NULL));
596:   PetscCall(PetscOptionsReal("-snes_trdc_eta3", "eta3", "None", ctx->eta3, &ctx->eta3, NULL));
597:   PetscCall(PetscOptionsReal("-snes_trdc_t1", "t1", "None", ctx->t1, &ctx->t1, NULL));
598:   PetscCall(PetscOptionsReal("-snes_trdc_t2", "t2", "None", ctx->t2, &ctx->t2, NULL));
599:   PetscCall(PetscOptionsReal("-snes_trdc_deltaM", "deltaM", "None", ctx->deltaM, &ctx->deltaM, NULL));
600:   PetscCall(PetscOptionsReal("-snes_trdc_delta0", "delta0", "None", ctx->delta0, &ctx->delta0, NULL));
601:   PetscCall(PetscOptionsReal("-snes_trdc_auto_scale_max", "auto_scale_max", "None", ctx->auto_scale_max, &ctx->auto_scale_max, NULL));
602:   PetscCall(PetscOptionsBool("-snes_trdc_use_cauchy", "use_cauchy", "use Cauchy step and direction", ctx->use_cauchy, &ctx->use_cauchy, NULL));
603:   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));
604:   PetscOptionsHeadEnd();
605:   PetscFunctionReturn(PETSC_SUCCESS);
606: }

608: static PetscErrorCode SNESView_NEWTONTRDC(SNES snes, PetscViewer viewer)
609: {
610:   SNES_NEWTONTRDC *tr = (SNES_NEWTONTRDC *)snes->data;
611:   PetscBool        iascii;

613:   PetscFunctionBegin;
614:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
615:   if (iascii) {
616:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Trust region tolerance %g (-snes_trtol)\n", (double)snes->deltatol));
617:     PetscCall(PetscViewerASCIIPrintf(viewer, "  eta1=%g, eta2=%g, eta3=%g\n", (double)tr->eta1, (double)tr->eta2, (double)tr->eta3));
618:     PetscCall(PetscViewerASCIIPrintf(viewer, "  delta0=%g, t1=%g, t2=%g, deltaM=%g\n", (double)tr->delta0, (double)tr->t1, (double)tr->t2, (double)tr->deltaM));
619:   }
620:   PetscFunctionReturn(PETSC_SUCCESS);
621: }

623: /*MC
624:       SNESNEWTONTRDC - Newton based nonlinear solver that uses trust-region dogleg method with Cauchy direction

626:    Options Database Keys:
627: +   -snes_trdc_tol <tol>                                     - trust region tolerance
628: .   -snes_trdc_eta1 <eta1>                                   - trust region parameter 0.0 <= eta1 <= eta2, rho >= eta1 breaks out of the inner iteration (default: eta1=0.001)
629: .   -snes_trdc_eta2 <eta2>                                   - trust region parameter 0.0 <= eta1 <= eta2, rho <= eta2 shrinks the trust region (default: eta2=0.25)
630: .   -snes_trdc_eta3 <eta3>                                   - trust region parameter eta3 > eta2, rho >= eta3 expands the trust region (default: eta3=0.75)
631: .   -snes_trdc_t1 <t1>                                       - trust region parameter, shrinking factor of trust region (default: 0.25)
632: .   -snes_trdc_t2 <t2>                                       - trust region parameter, expanding factor of trust region (default: 2.0)
633: .   -snes_trdc_deltaM <deltaM>                               - trust region parameter, max size of trust region, $deltaM*norm2(x)$ (default: 0.5)
634: .   -snes_trdc_delta0 <delta0>                               - trust region parameter, initial size of trust region, $delta0*norm2(x)$ (default: 0.1)
635: .   -snes_trdc_auto_scale_max <auto_scale_max>               - used with auto_scale_multiphase, caps the maximum auto-scaling factor
636: .   -snes_trdc_use_cauchy <use_cauchy>                       - True uses dogleg Cauchy (Steepest Descent direction) step & direction in the trust region algorithm
637: -   -snes_trdc_auto_scale_multiphase <auto_scale_multiphase> - True turns on auto-scaling for multivariable block matrix for Cauchy and trust region

639:    Level: intermediate

641:    Note:
642:    See {cite}`park2021linear`

644: .seealso: [](ch_snes), `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESSetTrustRegionTolerance()`,
645:           `SNESNewtonTRDCPreCheck()`, `SNESNewtonTRDCGetPreCheck()`, `SNESNewtonTRDCSetPostCheck()`, `SNESNewtonTRDCGetPostCheck()`,
646:           `SNESNewtonTRDCGetRhoFlag()`, `SNESNewtonTRDCSetPreCheck()`
647: M*/
648: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONTRDC(SNES snes)
649: {
650:   SNES_NEWTONTRDC *neP;

652:   PetscFunctionBegin;
653:   snes->ops->setup          = SNESSetUp_NEWTONTRDC;
654:   snes->ops->solve          = SNESSolve_NEWTONTRDC;
655:   snes->ops->destroy        = SNESDestroy_NEWTONTRDC;
656:   snes->ops->setfromoptions = SNESSetFromOptions_NEWTONTRDC;
657:   snes->ops->view           = SNESView_NEWTONTRDC;
658:   snes->ops->reset          = SNESReset_NEWTONTRDC;

660:   snes->usesksp = PETSC_TRUE;
661:   snes->usesnpc = PETSC_FALSE;

663:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

665:   PetscCall(PetscNew(&neP));
666:   snes->data                 = (void *)neP;
667:   neP->delta                 = 0.0;
668:   neP->delta0                = 0.1;
669:   neP->eta1                  = 0.001;
670:   neP->eta2                  = 0.25;
671:   neP->eta3                  = 0.75;
672:   neP->t1                    = 0.25;
673:   neP->t2                    = 2.0;
674:   neP->deltaM                = 0.5;
675:   neP->sigma                 = 0.0001;
676:   neP->itflag                = PETSC_FALSE;
677:   neP->rnorm0                = 0.0;
678:   neP->ttol                  = 0.0;
679:   neP->use_cauchy            = PETSC_TRUE;
680:   neP->auto_scale_multiphase = PETSC_FALSE;
681:   neP->auto_scale_max        = -1.0;
682:   neP->rho_satisfied         = PETSC_FALSE;
683:   snes->deltatol             = 1.e-12;

685:   /* for multiphase (multivariable) scaling */
686:   /* may be used for dynamic allocation of inorms, but it fails snes_tutorials-ex3_13
687:      on test forced DIVERGED_JACOBIAN_DOMAIN test. I will use static array for now.
688:   PetscCall(VecGetBlockSize(snes->work[0],&neP->bs));
689:   PetscCall(PetscCalloc1(neP->bs,&neP->inorms));
690:   */

692:   PetscFunctionReturn(PETSC_SUCCESS);
693: }