Actual source code: morethuente.c

      1: #include <petsc/private/taolinesearchimpl.h>
      2: #include <../src/tao/linesearch/impls/morethuente/morethuente.h>

      4: /*
      5:    This algorithm is taken from More' and Thuente, "Line search algorithms
      6:    with guaranteed sufficient decrease", Argonne National Laboratory,
      7:    Technical Report MCS-P330-1092.
      8: */

     10: static PetscErrorCode Tao_mcstep(TaoLineSearch ls, PetscReal *stx, PetscReal *fx, PetscReal *dx, PetscReal *sty, PetscReal *fy, PetscReal *dy, PetscReal *stp, PetscReal *fp, PetscReal *dp);

     12: static PetscErrorCode TaoLineSearchDestroy_MT(TaoLineSearch ls)
     13: {
     14:   TaoLineSearch_MT *mt = (TaoLineSearch_MT *)ls->data;

     16:   PetscFunctionBegin;
     17:   PetscCall(PetscObjectDereference((PetscObject)mt->x));
     18:   PetscCall(VecDestroy(&mt->work));
     19:   PetscCall(PetscFree(ls->data));
     20:   PetscFunctionReturn(PETSC_SUCCESS);
     21: }

     23: static PetscErrorCode TaoLineSearchSetFromOptions_MT(TaoLineSearch ls, PetscOptionItems *PetscOptionsObject)
     24: {
     25:   PetscFunctionBegin;
     26:   PetscFunctionReturn(PETSC_SUCCESS);
     27: }

     29: static PetscErrorCode TaoLineSearchMonitor_MT(TaoLineSearch ls)
     30: {
     31:   TaoLineSearch_MT *mt = (TaoLineSearch_MT *)ls->data;

     33:   PetscFunctionBegin;
     34:   PetscCall(PetscViewerASCIIPrintf(ls->viewer, "stx: %g, fx: %g, dgx: %g\n", (double)mt->stx, (double)mt->fx, (double)mt->dgx));
     35:   PetscCall(PetscViewerASCIIPrintf(ls->viewer, "sty: %g, fy: %g, dgy: %g\n", (double)mt->sty, (double)mt->fy, (double)mt->dgy));
     36:   PetscFunctionReturn(PETSC_SUCCESS);
     37: }

     39: static PetscErrorCode TaoLineSearchApply_MT(TaoLineSearch ls, Vec x, PetscReal *f, Vec g, Vec s)
     40: {
     41:   TaoLineSearch_MT *mt     = (TaoLineSearch_MT *)ls->data;
     42:   PetscReal         xtrapf = 4.0;
     43:   PetscReal         finit, width, width1, dginit, fm, fxm, fym, dgm, dgxm, dgym;
     44:   PetscReal         dgx, dgy, dg, dg2, fx, fy, stx, sty, dgtest;
     45:   PetscReal         ftest1 = 0.0, ftest2 = 0.0;
     46:   PetscInt          i, stage1, n1, n2, nn1, nn2;
     47:   PetscReal         bstepmin1, bstepmin2, bstepmax, ostepmin, ostepmax;
     48:   PetscBool         g_computed = PETSC_FALSE; /* to prevent extra gradient computation */

     50:   PetscFunctionBegin;
     51:   ls->reason = TAOLINESEARCH_CONTINUE_ITERATING;
     52:   PetscCall(TaoLineSearchMonitor(ls, 0, *f, 0.0));
     53:   /* Check work vector */
     54:   if (!mt->work) {
     55:     PetscCall(VecDuplicate(x, &mt->work));
     56:     mt->x = x;
     57:     PetscCall(PetscObjectReference((PetscObject)mt->x));
     58:   } else if (x != mt->x) {
     59:     PetscCall(VecDestroy(&mt->work));
     60:     PetscCall(VecDuplicate(x, &mt->work));
     61:     PetscCall(PetscObjectDereference((PetscObject)mt->x));
     62:     mt->x = x;
     63:     PetscCall(PetscObjectReference((PetscObject)mt->x));
     64:   }

     66:   ostepmax = ls->stepmax;
     67:   ostepmin = ls->stepmin;

     69:   if (ls->bounded) {
     70:     /* Compute step length needed to make all variables equal a bound */
     71:     /* Compute the smallest steplength that will make one nonbinding variable
     72:      equal the bound */
     73:     PetscCall(VecGetLocalSize(ls->upper, &n1));
     74:     PetscCall(VecGetLocalSize(mt->x, &n2));
     75:     PetscCall(VecGetSize(ls->upper, &nn1));
     76:     PetscCall(VecGetSize(mt->x, &nn2));
     77:     PetscCheck(n1 == n2 && nn1 == nn2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Variable vector not compatible with bounds vector");
     78:     PetscCall(VecScale(s, -1.0));
     79:     PetscCall(VecBoundGradientProjection(s, x, ls->lower, ls->upper, s));
     80:     PetscCall(VecScale(s, -1.0));
     81:     PetscCall(VecStepBoundInfo(x, s, ls->lower, ls->upper, &bstepmin1, &bstepmin2, &bstepmax));
     82:     ls->stepmax = PetscMin(bstepmax, ls->stepmax);
     83:   }

     85:   PetscCall(VecDot(g, s, &dginit));
     86:   if (PetscIsInfOrNanReal(dginit)) {
     87:     PetscCall(PetscInfo(ls, "Initial Line Search step * g is Inf or Nan (%g)\n", (double)dginit));
     88:     ls->reason = TAOLINESEARCH_FAILED_INFORNAN;
     89:     PetscFunctionReturn(PETSC_SUCCESS);
     90:   }
     91:   if (dginit >= 0.0) {
     92:     PetscCall(PetscInfo(ls, "Initial Line Search step * g is not descent direction (%g)\n", (double)dginit));
     93:     ls->reason = TAOLINESEARCH_FAILED_ASCENT;
     94:     PetscFunctionReturn(PETSC_SUCCESS);
     95:   }

     97:   /* Initialization */
     98:   mt->bracket = 0;
     99:   stage1      = 1;
    100:   finit       = *f;
    101:   dgtest      = ls->ftol * dginit;
    102:   width       = ls->stepmax - ls->stepmin;
    103:   width1      = width * 2.0;
    104:   PetscCall(VecCopy(x, mt->work));
    105:   /* Variable dictionary:
    106:    stx, fx, dgx - the step, function, and derivative at the best step
    107:    sty, fy, dgy - the step, function, and derivative at the other endpoint
    108:    of the interval of uncertainty
    109:    step, f, dg - the step, function, and derivative at the current step */

    111:   stx = 0.0;
    112:   fx  = finit;
    113:   dgx = dginit;
    114:   sty = 0.0;
    115:   fy  = finit;
    116:   dgy = dginit;

    118:   ls->step = ls->initstep;
    119:   for (i = 0; i < ls->max_funcs; i++) {
    120:     /* Set min and max steps to correspond to the interval of uncertainty */
    121:     if (mt->bracket) {
    122:       ls->stepmin = PetscMin(stx, sty);
    123:       ls->stepmax = PetscMax(stx, sty);
    124:     } else {
    125:       ls->stepmin = stx;
    126:       ls->stepmax = ls->step + xtrapf * (ls->step - stx);
    127:     }

    129:     /* Force the step to be within the bounds */
    130:     ls->step = PetscMax(ls->step, ls->stepmin);
    131:     ls->step = PetscMin(ls->step, ls->stepmax);

    133:     /* If an unusual termination is to occur, then let step be the lowest
    134:      point obtained thus far */
    135:     if (stx != 0 && ((mt->bracket && (ls->step <= ls->stepmin || ls->step >= ls->stepmax)) || (mt->bracket && (ls->stepmax - ls->stepmin <= ls->rtol * ls->stepmax)) || (ls->nfeval + ls->nfgeval >= ls->max_funcs - 1) || mt->infoc == 0))
    136:       ls->step = stx;

    138:     PetscCall(VecWAXPY(mt->work, ls->step, s, x)); /* W = X + step*S */

    140:     if (ls->step == 0.0) {
    141:       PetscCall(PetscInfo(ls, "Step size is zero.\n"));
    142:       ls->reason = TAOLINESEARCH_HALTED_LOWERBOUND;
    143:       break;
    144:     }

    146:     if (ls->bounded) PetscCall(VecMedian(ls->lower, mt->work, ls->upper, mt->work));
    147:     if (ls->usegts) {
    148:       PetscCall(TaoLineSearchComputeObjectiveAndGTS(ls, mt->work, f, &dg));
    149:       g_computed = PETSC_FALSE;
    150:     } else {
    151:       PetscCall(TaoLineSearchComputeObjectiveAndGradient(ls, mt->work, f, g));
    152:       g_computed = PETSC_TRUE;
    153:       if (ls->bounded) {
    154:         PetscCall(VecDot(g, x, &dg));
    155:         PetscCall(VecDot(g, mt->work, &dg2));
    156:         dg = (dg2 - dg) / ls->step;
    157:       } else {
    158:         PetscCall(VecDot(g, s, &dg));
    159:       }
    160:     }

    162:     /* update bracketing parameters in the MT context for printouts in monitor */
    163:     mt->stx = stx;
    164:     mt->fx  = fx;
    165:     mt->dgx = dgx;
    166:     mt->sty = sty;
    167:     mt->fy  = fy;
    168:     mt->dgy = dgy;
    169:     PetscCall(TaoLineSearchMonitor(ls, i + 1, *f, ls->step));

    171:     if (i == 0) ls->f_fullstep = *f;

    173:     if (PetscIsInfOrNanReal(*f) || PetscIsInfOrNanReal(dg)) {
    174:       /* User provided compute function generated Not-a-Number, assume
    175:        domain violation and set function value and directional
    176:        derivative to infinity. */
    177:       *f = PETSC_INFINITY;
    178:       dg = PETSC_INFINITY;
    179:     }

    181:     ftest1 = finit + ls->step * dgtest;
    182:     if (ls->bounded) ftest2 = finit + ls->step * dgtest * ls->ftol;

    184:     /* Convergence testing */
    185:     if ((*f - ftest1 <= PETSC_SMALL * PetscAbsReal(finit)) && (PetscAbsReal(dg) + ls->gtol * dginit <= 0.0)) {
    186:       PetscCall(PetscInfo(ls, "Line search success: Sufficient decrease and directional deriv conditions hold\n"));
    187:       ls->reason = TAOLINESEARCH_SUCCESS;
    188:       break;
    189:     }

    191:     /* Check Armijo if beyond the first breakpoint */
    192:     if (ls->bounded && *f <= ftest2 && ls->step >= bstepmin2) {
    193:       PetscCall(PetscInfo(ls, "Line search success: Sufficient decrease.\n"));
    194:       ls->reason = TAOLINESEARCH_SUCCESS;
    195:       break;
    196:     }

    198:     /* Checks for bad cases */
    199:     if ((mt->bracket && (ls->step <= ls->stepmin || ls->step >= ls->stepmax)) || !mt->infoc) {
    200:       PetscCall(PetscInfo(ls, "Rounding errors may prevent further progress. May not be a step satisfying\nsufficient decrease and curvature conditions. Tolerances may be too small.\n"));
    201:       ls->reason = TAOLINESEARCH_HALTED_OTHER;
    202:       break;
    203:     }
    204:     if (ls->step == ls->stepmax && *f <= ftest1 && dg <= dgtest) {
    205:       PetscCall(PetscInfo(ls, "Step is at the upper bound, stepmax (%g)\n", (double)ls->stepmax));
    206:       ls->reason = TAOLINESEARCH_HALTED_UPPERBOUND;
    207:       break;
    208:     }
    209:     if (ls->step == ls->stepmin && *f >= ftest1 && dg >= dgtest) {
    210:       PetscCall(PetscInfo(ls, "Step is at the lower bound, stepmin (%g)\n", (double)ls->stepmin));
    211:       ls->reason = TAOLINESEARCH_HALTED_LOWERBOUND;
    212:       break;
    213:     }
    214:     if (mt->bracket && (ls->stepmax - ls->stepmin <= ls->rtol * ls->stepmax)) {
    215:       PetscCall(PetscInfo(ls, "Relative width of interval of uncertainty is at most rtol (%g)\n", (double)ls->rtol));
    216:       ls->reason = TAOLINESEARCH_HALTED_RTOL;
    217:       break;
    218:     }

    220:     /* In the first stage, we seek a step for which the modified function
    221:      has a nonpositive value and nonnegative derivative */
    222:     if (stage1 && *f <= ftest1 && dg >= dginit * PetscMin(ls->ftol, ls->gtol)) stage1 = 0;

    224:     /* A modified function is used to predict the step only if we
    225:      have not obtained a step for which the modified function has a
    226:      nonpositive function value and nonnegative derivative, and if a
    227:      lower function value has been obtained but the decrease is not
    228:      sufficient */

    230:     if (stage1 && *f <= fx && *f > ftest1) {
    231:       fm   = *f - ls->step * dgtest; /* Define modified function */
    232:       fxm  = fx - stx * dgtest;      /* and derivatives */
    233:       fym  = fy - sty * dgtest;
    234:       dgm  = dg - dgtest;
    235:       dgxm = dgx - dgtest;
    236:       dgym = dgy - dgtest;

    238:       /* if (dgxm * (ls->step - stx) >= 0.0) */
    239:       /* Update the interval of uncertainty and compute the new step */
    240:       PetscCall(Tao_mcstep(ls, &stx, &fxm, &dgxm, &sty, &fym, &dgym, &ls->step, &fm, &dgm));

    242:       fx  = fxm + stx * dgtest; /* Reset the function and */
    243:       fy  = fym + sty * dgtest; /* gradient values */
    244:       dgx = dgxm + dgtest;
    245:       dgy = dgym + dgtest;
    246:     } else {
    247:       /* Update the interval of uncertainty and compute the new step */
    248:       PetscCall(Tao_mcstep(ls, &stx, &fx, &dgx, &sty, &fy, &dgy, &ls->step, f, &dg));
    249:     }

    251:     /* Force a sufficient decrease in the interval of uncertainty */
    252:     if (mt->bracket) {
    253:       if (PetscAbsReal(sty - stx) >= 0.66 * width1) ls->step = stx + 0.5 * (sty - stx);
    254:       width1 = width;
    255:       width  = PetscAbsReal(sty - stx);
    256:     }
    257:   }
    258:   if (ls->nfeval + ls->nfgeval > ls->max_funcs) {
    259:     PetscCall(PetscInfo(ls, "Number of line search function evals (%" PetscInt_FMT ") > maximum (%" PetscInt_FMT ")\n", ls->nfeval + ls->nfgeval, ls->max_funcs));
    260:     ls->reason = TAOLINESEARCH_HALTED_MAXFCN;
    261:   }
    262:   ls->stepmax = ostepmax;
    263:   ls->stepmin = ostepmin;

    265:   /* Finish computations */
    266:   PetscCall(PetscInfo(ls, "%" PetscInt_FMT " function evals in line search, step = %g\n", ls->nfeval + ls->nfgeval, (double)ls->step));

    268:   /* Set new solution vector and compute gradient if needed */
    269:   PetscCall(VecCopy(mt->work, x));
    270:   if (!g_computed) PetscCall(TaoLineSearchComputeGradient(ls, mt->work, g));
    271:   PetscFunctionReturn(PETSC_SUCCESS);
    272: }

    274: /*MC
    275:    TAOLINESEARCHMT - More-Thuente line-search type with cubic interpolation that satisfies both the sufficient decrease and
    276:    curvature conditions. This method can take step lengths greater than 1, {cite}`more:92`

    278:    Options Database Key:
    279: .  -tao_ls_type more-thuente - use this line search type

    281:    Level: developer

    283: .seealso: `TaoLineSearchCreate()`, `TaoLineSearchSetType()`, `TaoLineSearchApply()`
    284: M*/
    285: PETSC_EXTERN PetscErrorCode TaoLineSearchCreate_MT(TaoLineSearch ls)
    286: {
    287:   TaoLineSearch_MT *ctx;

    289:   PetscFunctionBegin;
    291:   PetscCall(PetscNew(&ctx));
    292:   ctx->bracket            = 0;
    293:   ctx->infoc              = 1;
    294:   ls->data                = (void *)ctx;
    295:   ls->initstep            = 1.0;
    296:   ls->ops->setup          = NULL;
    297:   ls->ops->reset          = NULL;
    298:   ls->ops->apply          = TaoLineSearchApply_MT;
    299:   ls->ops->destroy        = TaoLineSearchDestroy_MT;
    300:   ls->ops->setfromoptions = TaoLineSearchSetFromOptions_MT;
    301:   ls->ops->monitor        = TaoLineSearchMonitor_MT;
    302:   PetscFunctionReturn(PETSC_SUCCESS);
    303: }

    305: /*
    306:      The subroutine mcstep is taken from the work of Jorge Nocedal.
    307:      this is a variant of More' and Thuente's routine.

    309:      subroutine mcstep

    311:      the purpose of mcstep is to compute a safeguarded step for
    312:      a linesearch and to update an interval of uncertainty for
    313:      a minimizer of the function.

    315:      the parameter stx contains the step with the least function
    316:      value. the parameter stp contains the current step. it is
    317:      assumed that the derivative at stx is negative in the
    318:      direction of the step. if bracket is set true then a
    319:      minimizer has been bracketed in an interval of uncertainty
    320:      with endpoints stx and sty.

    322:      the subroutine statement is

    324:      subroutine mcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,bracket,
    325:                        stpmin,stpmax,info)

    327:      where

    329:        stx, fx, and dx are variables which specify the step,
    330:          the function, and the derivative at the best step obtained
    331:          so far. The derivative must be negative in the direction
    332:          of the step, that is, dx and stp-stx must have opposite
    333:          signs. On output these parameters are updated appropriately.

    335:        sty, fy, and dy are variables which specify the step,
    336:          the function, and the derivative at the other endpoint of
    337:          the interval of uncertainty. On output these parameters are
    338:          updated appropriately.

    340:        stp, fp, and dp are variables which specify the step,
    341:          the function, and the derivative at the current step.
    342:          If bracket is set true then on input stp must be
    343:          between stx and sty. On output stp is set to the new step.

    345:        bracket is a logical variable which specifies if a minimizer
    346:          has been bracketed.  If the minimizer has not been bracketed
    347:          then on input bracket must be set false.  If the minimizer
    348:          is bracketed then on output bracket is set true.

    350:        stpmin and stpmax are input variables which specify lower
    351:          and upper bounds for the step.

    353:        info is an integer output variable set as follows:
    354:          if info = 1,2,3,4,5, then the step has been computed
    355:          according to one of the five cases below. otherwise
    356:          info = 0, and this indicates improper input parameters.

    358:      subprograms called

    360:        fortran-supplied ... abs,max,min,sqrt

    362:      argonne national laboratory. minpack project. june 1983
    363:      jorge j. more', david j. thuente

    365: */

    367: static PetscErrorCode Tao_mcstep(TaoLineSearch ls, PetscReal *stx, PetscReal *fx, PetscReal *dx, PetscReal *sty, PetscReal *fy, PetscReal *dy, PetscReal *stp, PetscReal *fp, PetscReal *dp)
    368: {
    369:   TaoLineSearch_MT *mtP = (TaoLineSearch_MT *)ls->data;
    370:   PetscReal         gamma1, p, q, r, s, sgnd, stpc, stpf, stpq, theta;
    371:   PetscInt          bound;

    373:   PetscFunctionBegin;
    374:   /* Check the input parameters for errors */
    375:   mtP->infoc = 0;
    376:   PetscCheck(!mtP->bracket || (*stp > PetscMin(*stx, *sty) && *stp < PetscMax(*stx, *sty)), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "bad stp in bracket");
    377:   PetscCheck(*dx * (*stp - *stx) < 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "dx * (stp-stx) >= 0.0");
    378:   PetscCheck(ls->stepmax >= ls->stepmin, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "stepmax > stepmin");

    380:   /* Determine if the derivatives have opposite sign */
    381:   sgnd = *dp * (*dx / PetscAbsReal(*dx));

    383:   if (*fp > *fx) {
    384:     /* Case 1: a higher function value.
    385:      The minimum is bracketed. If the cubic step is closer
    386:      to stx than the quadratic step, the cubic step is taken,
    387:      else the average of the cubic and quadratic steps is taken. */

    389:     mtP->infoc = 1;
    390:     bound      = 1;
    391:     theta      = 3 * (*fx - *fp) / (*stp - *stx) + *dx + *dp;
    392:     s          = PetscMax(PetscAbsReal(theta), PetscAbsReal(*dx));
    393:     s          = PetscMax(s, PetscAbsReal(*dp));
    394:     gamma1     = s * PetscSqrtScalar(PetscPowScalar(theta / s, 2.0) - (*dx / s) * (*dp / s));
    395:     if (*stp < *stx) gamma1 = -gamma1;
    396:     /* Can p be 0?  Check */
    397:     p    = (gamma1 - *dx) + theta;
    398:     q    = ((gamma1 - *dx) + gamma1) + *dp;
    399:     r    = p / q;
    400:     stpc = *stx + r * (*stp - *stx);
    401:     stpq = *stx + ((*dx / ((*fx - *fp) / (*stp - *stx) + *dx)) * 0.5) * (*stp - *stx);

    403:     if (PetscAbsReal(stpc - *stx) < PetscAbsReal(stpq - *stx)) stpf = stpc;
    404:     else stpf = stpc + 0.5 * (stpq - stpc);
    405:     mtP->bracket = 1;
    406:   } else if (sgnd < 0.0) {
    407:     /* Case 2: A lower function value and derivatives of
    408:      opposite sign. The minimum is bracketed. If the cubic
    409:      step is closer to stx than the quadratic (secant) step,
    410:      the cubic step is taken, else the quadratic step is taken. */

    412:     mtP->infoc = 2;
    413:     bound      = 0;
    414:     theta      = 3 * (*fx - *fp) / (*stp - *stx) + *dx + *dp;
    415:     s          = PetscMax(PetscAbsReal(theta), PetscAbsReal(*dx));
    416:     s          = PetscMax(s, PetscAbsReal(*dp));
    417:     gamma1     = s * PetscSqrtScalar(PetscPowScalar(theta / s, 2.0) - (*dx / s) * (*dp / s));
    418:     if (*stp > *stx) gamma1 = -gamma1;
    419:     p    = (gamma1 - *dp) + theta;
    420:     q    = ((gamma1 - *dp) + gamma1) + *dx;
    421:     r    = p / q;
    422:     stpc = *stp + r * (*stx - *stp);
    423:     stpq = *stp + (*dp / (*dp - *dx)) * (*stx - *stp);

    425:     if (PetscAbsReal(stpc - *stp) > PetscAbsReal(stpq - *stp)) stpf = stpc;
    426:     else stpf = stpq;
    427:     mtP->bracket = 1;
    428:   } else if (PetscAbsReal(*dp) < PetscAbsReal(*dx)) {
    429:     /* Case 3: A lower function value, derivatives of the
    430:      same sign, and the magnitude of the derivative decreases.
    431:      The cubic step is only used if the cubic tends to infinity
    432:      in the direction of the step or if the minimum of the cubic
    433:      is beyond stp. Otherwise the cubic step is defined to be
    434:      either stepmin or stepmax. The quadratic (secant) step is also
    435:      computed and if the minimum is bracketed then the step
    436:      closest to stx is taken, else the step farthest away is taken. */

    438:     mtP->infoc = 3;
    439:     bound      = 1;
    440:     theta      = 3 * (*fx - *fp) / (*stp - *stx) + *dx + *dp;
    441:     s          = PetscMax(PetscAbsReal(theta), PetscAbsReal(*dx));
    442:     s          = PetscMax(s, PetscAbsReal(*dp));

    444:     /* The case gamma1 = 0 only arises if the cubic does not tend
    445:        to infinity in the direction of the step. */
    446:     gamma1 = s * PetscSqrtScalar(PetscMax(0.0, PetscPowScalar(theta / s, 2.0) - (*dx / s) * (*dp / s)));
    447:     if (*stp > *stx) gamma1 = -gamma1;
    448:     p = (gamma1 - *dp) + theta;
    449:     q = (gamma1 + (*dx - *dp)) + gamma1;
    450:     r = p / q;
    451:     if (r < 0.0 && gamma1 != 0.0) stpc = *stp + r * (*stx - *stp);
    452:     else if (*stp > *stx) stpc = ls->stepmax;
    453:     else stpc = ls->stepmin;
    454:     stpq = *stp + (*dp / (*dp - *dx)) * (*stx - *stp);

    456:     if (mtP->bracket) {
    457:       if (PetscAbsReal(*stp - stpc) < PetscAbsReal(*stp - stpq)) stpf = stpc;
    458:       else stpf = stpq;
    459:     } else {
    460:       if (PetscAbsReal(*stp - stpc) > PetscAbsReal(*stp - stpq)) stpf = stpc;
    461:       else stpf = stpq;
    462:     }
    463:   } else {
    464:     /* Case 4: A lower function value, derivatives of the
    465:        same sign, and the magnitude of the derivative does
    466:        not decrease. If the minimum is not bracketed, the step
    467:        is either stpmin or stpmax, else the cubic step is taken. */

    469:     mtP->infoc = 4;
    470:     bound      = 0;
    471:     if (mtP->bracket) {
    472:       theta  = 3 * (*fp - *fy) / (*sty - *stp) + *dy + *dp;
    473:       s      = PetscMax(PetscAbsReal(theta), PetscAbsReal(*dy));
    474:       s      = PetscMax(s, PetscAbsReal(*dp));
    475:       gamma1 = s * PetscSqrtScalar(PetscPowScalar(theta / s, 2.0) - (*dy / s) * (*dp / s));
    476:       if (*stp > *sty) gamma1 = -gamma1;
    477:       p    = (gamma1 - *dp) + theta;
    478:       q    = ((gamma1 - *dp) + gamma1) + *dy;
    479:       r    = p / q;
    480:       stpc = *stp + r * (*sty - *stp);
    481:       stpf = stpc;
    482:     } else if (*stp > *stx) {
    483:       stpf = ls->stepmax;
    484:     } else {
    485:       stpf = ls->stepmin;
    486:     }
    487:   }

    489:   /* Update the interval of uncertainty.  This update does not
    490:      depend on the new step or the case analysis above. */

    492:   if (*fp > *fx) {
    493:     *sty = *stp;
    494:     *fy  = *fp;
    495:     *dy  = *dp;
    496:   } else {
    497:     if (sgnd < 0.0) {
    498:       *sty = *stx;
    499:       *fy  = *fx;
    500:       *dy  = *dx;
    501:     }
    502:     *stx = *stp;
    503:     *fx  = *fp;
    504:     *dx  = *dp;
    505:   }

    507:   /* Compute the new step and safeguard it. */
    508:   stpf = PetscMin(ls->stepmax, stpf);
    509:   stpf = PetscMax(ls->stepmin, stpf);
    510:   *stp = stpf;
    511:   if (mtP->bracket && bound) {
    512:     if (*sty > *stx) *stp = PetscMin(*stx + 0.66 * (*sty - *stx), *stp);
    513:     else *stp = PetscMax(*stx + 0.66 * (*sty - *stx), *stp);
    514:   }
    515:   PetscFunctionReturn(PETSC_SUCCESS);
    516: }