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, 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;

     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(void **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] user-defined context for private data for the function evaluation routine (may be `NULL`)

    117:   Level: intermediate

    119:   Note:
    120:   This function is called BEFORE the function evaluation within the `SNESNEWTONTRDC` solver.

    122: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTRDC`, `SNESNewtonTRDCPreCheck()`, `SNESNewtonTRDCGetPreCheck()`, `SNESNewtonTRDCSetPostCheck()`, `SNESNewtonTRDCGetPostCheck()`,
    123:           `SNESNewtonTRDCGetRhoFlag()`
    124: @*/
    125: PetscErrorCode SNESNewtonTRDCSetPreCheck(SNES snes, PetscErrorCode (*func)(SNES, Vec, Vec, PetscBool *, void *), void *ctx)
    126: {
    127:   SNES_NEWTONTRDC *tr = (SNES_NEWTONTRDC *)snes->data;

    129:   PetscFunctionBegin;
    131:   if (func) tr->precheck = func;
    132:   if (ctx) tr->precheckctx = ctx;
    133:   PetscFunctionReturn(PETSC_SUCCESS);
    134: }

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

    139:   Not Collective

    141:   Input Parameter:
    142: . snes - the nonlinear solver context

    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`)

    148:   Level: intermediate

    150: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTRDC`, `SNESNewtonTRDCSetPreCheck()`, `SNESNewtonTRDCPreCheck()`
    151: @*/
    152: PetscErrorCode SNESNewtonTRDCGetPreCheck(SNES snes, PetscErrorCode (**func)(SNES, Vec, Vec, PetscBool *, void *), void **ctx)
    153: {
    154:   SNES_NEWTONTRDC *tr = (SNES_NEWTONTRDC *)snes->data;

    156:   PetscFunctionBegin;
    158:   if (func) *func = tr->precheck;
    159:   if (ctx) *ctx = tr->precheckctx;
    160:   PetscFunctionReturn(PETSC_SUCCESS);
    161: }

    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

    167:   Logically Collective

    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`)

    174:   Level: intermediate

    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.

    180: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTRDC`, `SNESNewtonTRDCPostCheck()`, `SNESNewtonTRDCGetPostCheck()`, `SNESNewtonTRDCSetPreCheck()`, `SNESNewtonTRDCGetPreCheck()`
    181: @*/
    182: PetscErrorCode SNESNewtonTRDCSetPostCheck(SNES snes, PetscErrorCode (*func)(SNES, Vec, Vec, Vec, PetscBool *, PetscBool *, void *), void *ctx)
    183: {
    184:   SNES_NEWTONTRDC *tr = (SNES_NEWTONTRDC *)snes->data;

    186:   PetscFunctionBegin;
    188:   if (func) tr->postcheck = func;
    189:   if (ctx) tr->postcheckctx = ctx;
    190:   PetscFunctionReturn(PETSC_SUCCESS);
    191: }

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

    196:   Not Collective

    198:   Input Parameter:
    199: . snes - the nonlinear solver context

    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`)

    205:   Level: intermediate

    207: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTRDC`, `SNESNewtonTRDCSetPostCheck()`, `SNESNewtonTRDCPostCheck()`, `SNESNewtonTRDCSetPreCheck()`, `SNESNewtonTRDCGetPreCheck()`
    208: @*/
    209: PetscErrorCode SNESNewtonTRDCGetPostCheck(SNES snes, PetscErrorCode (**func)(SNES, Vec, Vec, Vec, PetscBool *, PetscBool *, void *), void **ctx)
    210: {
    211:   SNES_NEWTONTRDC *tr = (SNES_NEWTONTRDC *)snes->data;

    213:   PetscFunctionBegin;
    215:   if (func) *func = tr->postcheck;
    216:   if (ctx) *ctx = tr->postcheckctx;
    217:   PetscFunctionReturn(PETSC_SUCCESS);
    218: }

    220: // PetscClangLinter pragma disable: -fdoc-internal-linkage
    221: /*@C
    222:    SNESNewtonTRDCPreCheck - Called before the step has been determined in `SNESNEWTONTRDC`

    224:    Logically Collective

    226:    Input Parameters:
    227: +  snes - the solver
    228: .  X - The last solution
    229: -  Y - The step direction

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

    234:    Level: developer

    236: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTRDC`, `SNESNewtonTRDCSetPreCheck()`, `SNESNewtonTRDCGetPreCheck()`, `SNESNewtonTRDCPostCheck()`
    237: @*/
    238: static PetscErrorCode SNESNewtonTRDCPreCheck(SNES snes, Vec X, Vec Y, PetscBool *changed_Y)
    239: {
    240:   SNES_NEWTONTRDC *tr = (SNES_NEWTONTRDC *)snes->data;

    242:   PetscFunctionBegin;
    243:   *changed_Y = PETSC_FALSE;
    244:   if (tr->precheck) {
    245:     PetscCall((*tr->precheck)(snes, X, Y, changed_Y, tr->precheckctx));
    247:   }
    248:   PetscFunctionReturn(PETSC_SUCCESS);
    249: }

    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

    255:    Logically Collective

    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

    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.

    267:    Level: developer

    269:    Note:
    270:      If `Y` is changed then `W` is recomputed as `X` - `Y`

    272: .seealso: [](ch_snes), `SNES`, `SNESNEWTONTRDC`, `SNESNEWTONTRDC`, `SNESNewtonTRDCSetPostCheck()`, `SNESNewtonTRDCGetPostCheck()`, `SNESNewtonTRDCPreCheck()
    273: @*/
    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;

    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));
    285:   }
    286:   PetscFunctionReturn(PETSC_SUCCESS);
    287: }

    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

    294: */
    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;

    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 */

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

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

    329:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
    330:   snes->iter = 0;
    331:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));

    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:   }

    344:   if (!snes->vec_func_init_set) {
    345:     PetscCall(SNESComputeFunction(snes, X, F)); /* F(X) */
    346:   } else snes->vec_func_init_set = PETSC_FALSE;

    348:   PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- || F || */
    349:   SNESCheckFunctionNorm(snes, fnorm);
    350:   PetscCall(VecNorm(X, NORM_2, &xnorm)); /* xnorm <- || X || */

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

    361:   neP->rho_satisfied = PETSC_FALSE;

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

    367:   for (i = 0; i < maxits; i++) {
    368:     PetscBool changed_y;
    369:     PetscBool changed_w;

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

    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:     }

    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 || */

    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:       }

    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*/

    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++;

    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:       SNESCheckFunctionNorm(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 */

    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:       }

    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"));

    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:   }

    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: }

    577: static PetscErrorCode SNESSetUp_NEWTONTRDC(SNES snes)
    578: {
    579:   PetscFunctionBegin;
    580:   PetscCall(SNESSetWorkVecs(snes, 6));
    581:   PetscCall(SNESSetUpMatrices(snes));
    582:   PetscFunctionReturn(PETSC_SUCCESS);
    583: }

    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: }

    593: static PetscErrorCode SNESSetFromOptions_NEWTONTRDC(SNES snes, PetscOptionItems PetscOptionsObject)
    594: {
    595:   SNES_NEWTONTRDC *ctx = (SNES_NEWTONTRDC *)snes->data;

    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: }

    614: static PetscErrorCode SNESView_NEWTONTRDC(SNES snes, PetscViewer viewer)
    615: {
    616:   SNES_NEWTONTRDC *tr = (SNES_NEWTONTRDC *)snes->data;
    617:   PetscBool        isascii;

    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: }

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

    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

    645:    Level: advanced

    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`.

    652:    For details, see {cite}`park2021linear`

    654: .seealso: [](ch_snes), `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESNewtonTRSetTolerances()`,
    655:           `SNESNewtonTRDCPreCheck()`, `SNESNewtonTRDCGetPreCheck()`, `SNESNewtonTRDCSetPostCheck()`, `SNESNewtonTRDCGetPostCheck()`,
    656:           `SNESNewtonTRDCGetRhoFlag()`, `SNESNewtonTRDCSetPreCheck()`
    657: M*/
    658: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONTRDC(SNES snes)
    659: {
    660:   SNES_NEWTONTRDC *neP;

    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;

    669:   snes->usesksp = PETSC_TRUE;
    670:   snes->usesnpc = PETSC_FALSE;

    672:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

    674:   PetscCall(SNESParametersInitialize(snes));

    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;

    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: }