1 #include <petsc/private/taolinesearchimpl.h>
2 #include <../src/tao/linesearch/impls/morethuente/morethuente.h>
3
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 */
9
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);
11
TaoLineSearchDestroy_MT(TaoLineSearch ls)12 static PetscErrorCode TaoLineSearchDestroy_MT(TaoLineSearch ls)
13 {
14 TaoLineSearch_MT *mt = (TaoLineSearch_MT *)ls->data;
15
16 PetscFunctionBegin;
17 PetscCall(PetscObjectDereference((PetscObject)mt->x));
18 PetscCall(VecDestroy(&mt->work));
19 PetscCall(PetscFree(ls->data));
20 PetscFunctionReturn(PETSC_SUCCESS);
21 }
22
TaoLineSearchMonitor_MT(TaoLineSearch ls)23 static PetscErrorCode TaoLineSearchMonitor_MT(TaoLineSearch ls)
24 {
25 TaoLineSearch_MT *mt = (TaoLineSearch_MT *)ls->data;
26
27 PetscFunctionBegin;
28 PetscCall(PetscViewerASCIIPrintf(ls->viewer, "stx: %g, fx: %g, dgx: %g\n", (double)mt->stx, (double)mt->fx, (double)mt->dgx));
29 PetscCall(PetscViewerASCIIPrintf(ls->viewer, "sty: %g, fy: %g, dgy: %g\n", (double)mt->sty, (double)mt->fy, (double)mt->dgy));
30 PetscFunctionReturn(PETSC_SUCCESS);
31 }
32
TaoLineSearchApply_MT(TaoLineSearch ls,Vec x,PetscReal * f,Vec g,Vec s)33 static PetscErrorCode TaoLineSearchApply_MT(TaoLineSearch ls, Vec x, PetscReal *f, Vec g, Vec s)
34 {
35 TaoLineSearch_MT *mt = (TaoLineSearch_MT *)ls->data;
36 PetscReal xtrapf = 4.0;
37 PetscReal finit, width, width1, dginit, fm, fxm, fym, dgm, dgxm, dgym;
38 PetscReal dgx, dgy, dg, dg2, fx, fy, stx, sty, dgtest;
39 PetscReal ftest1 = 0.0, ftest2 = 0.0;
40 PetscInt i, stage1, n1, n2, nn1, nn2;
41 PetscReal bstepmin1, bstepmin2, bstepmax, ostepmin, ostepmax;
42 PetscBool g_computed = PETSC_FALSE; /* to prevent extra gradient computation */
43
44 PetscFunctionBegin;
45 ls->reason = TAOLINESEARCH_CONTINUE_ITERATING;
46 PetscCall(TaoLineSearchMonitor(ls, 0, *f, 0.0));
47 /* Check work vector */
48 if (!mt->work) {
49 PetscCall(VecDuplicate(x, &mt->work));
50 mt->x = x;
51 PetscCall(PetscObjectReference((PetscObject)mt->x));
52 } else if (x != mt->x) {
53 PetscCall(VecDestroy(&mt->work));
54 PetscCall(VecDuplicate(x, &mt->work));
55 PetscCall(PetscObjectDereference((PetscObject)mt->x));
56 mt->x = x;
57 PetscCall(PetscObjectReference((PetscObject)mt->x));
58 }
59
60 ostepmax = ls->stepmax;
61 ostepmin = ls->stepmin;
62
63 if (ls->bounded) {
64 /* Compute step length needed to make all variables equal a bound */
65 /* Compute the smallest steplength that will make one nonbinding variable
66 equal the bound */
67 PetscCall(VecGetLocalSize(ls->upper, &n1));
68 PetscCall(VecGetLocalSize(mt->x, &n2));
69 PetscCall(VecGetSize(ls->upper, &nn1));
70 PetscCall(VecGetSize(mt->x, &nn2));
71 PetscCheck(n1 == n2 && nn1 == nn2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Variable vector not compatible with bounds vector");
72 PetscCall(VecScale(s, -1.0));
73 PetscCall(VecBoundGradientProjection(s, x, ls->lower, ls->upper, s));
74 PetscCall(VecScale(s, -1.0));
75 PetscCall(VecStepBoundInfo(x, s, ls->lower, ls->upper, &bstepmin1, &bstepmin2, &bstepmax));
76 ls->stepmax = PetscMin(bstepmax, ls->stepmax);
77 }
78
79 PetscCall(VecDot(g, s, &dginit));
80 if (PetscIsInfOrNanReal(dginit)) {
81 PetscCall(PetscInfo(ls, "Initial Line Search step * g is infinity or NaN (%g)\n", (double)dginit));
82 ls->reason = TAOLINESEARCH_FAILED_INFORNAN;
83 PetscFunctionReturn(PETSC_SUCCESS);
84 }
85 if (dginit >= 0.0) {
86 PetscCall(PetscInfo(ls, "Initial Line Search step * g is not descent direction (%g)\n", (double)dginit));
87 ls->reason = TAOLINESEARCH_FAILED_ASCENT;
88 PetscFunctionReturn(PETSC_SUCCESS);
89 }
90
91 /* Initialization */
92 mt->bracket = 0;
93 stage1 = 1;
94 finit = *f;
95 dgtest = ls->ftol * dginit;
96 width = ls->stepmax - ls->stepmin;
97 width1 = width * 2.0;
98 PetscCall(VecCopy(x, mt->work));
99 /* Variable dictionary:
100 stx, fx, dgx - the step, function, and derivative at the best step
101 sty, fy, dgy - the step, function, and derivative at the other endpoint
102 of the interval of uncertainty
103 step, f, dg - the step, function, and derivative at the current step */
104
105 stx = 0.0;
106 fx = finit;
107 dgx = dginit;
108 sty = 0.0;
109 fy = finit;
110 dgy = dginit;
111
112 ls->step = ls->initstep;
113 for (i = 0; i < ls->max_funcs; i++) {
114 /* Set min and max steps to correspond to the interval of uncertainty */
115 if (mt->bracket) {
116 ls->stepmin = PetscMin(stx, sty);
117 ls->stepmax = PetscMax(stx, sty);
118 } else {
119 ls->stepmin = stx;
120 ls->stepmax = ls->step + xtrapf * (ls->step - stx);
121 }
122
123 /* Force the step to be within the bounds */
124 ls->step = PetscMax(ls->step, ls->stepmin);
125 ls->step = PetscMin(ls->step, ls->stepmax);
126
127 /* If an unusual termination is to occur, then let step be the lowest
128 point obtained thus far */
129 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))
130 ls->step = stx;
131
132 PetscCall(VecWAXPY(mt->work, ls->step, s, x)); /* W = X + step*S */
133
134 if (ls->step == 0.0) {
135 PetscCall(PetscInfo(ls, "Step size is zero.\n"));
136 ls->reason = TAOLINESEARCH_HALTED_LOWERBOUND;
137 break;
138 }
139
140 if (ls->bounded) PetscCall(VecMedian(ls->lower, mt->work, ls->upper, mt->work));
141 /* Make sure user code doesn't mess with the non-updated solution */
142 PetscCall(VecLockReadPush(x));
143 if (ls->usegts) {
144 PetscCall(TaoLineSearchComputeObjectiveAndGTS(ls, mt->work, f, &dg));
145 g_computed = PETSC_FALSE;
146 } else {
147 PetscCall(TaoLineSearchComputeObjectiveAndGradient(ls, mt->work, f, g));
148 g_computed = PETSC_TRUE;
149 if (ls->bounded) {
150 PetscCall(VecDot(g, x, &dg));
151 PetscCall(VecDot(g, mt->work, &dg2));
152 dg = (dg2 - dg) / ls->step;
153 } else {
154 PetscCall(VecDot(g, s, &dg));
155 }
156 }
157 PetscCall(VecLockReadPop(x));
158
159 /* update bracketing parameters in the MT context for printouts in monitor */
160 mt->stx = stx;
161 mt->fx = fx;
162 mt->dgx = dgx;
163 mt->sty = sty;
164 mt->fy = fy;
165 mt->dgy = dgy;
166 PetscCall(TaoLineSearchMonitor(ls, i + 1, *f, ls->step));
167
168 if (i == 0) ls->f_fullstep = *f;
169
170 if (PetscIsInfOrNanReal(*f) || PetscIsInfOrNanReal(dg)) {
171 /* User provided compute function generated Not-a-Number, assume
172 domain violation and set function value and directional
173 derivative to infinity. */
174 *f = PETSC_INFINITY;
175 dg = PETSC_INFINITY;
176 }
177
178 ftest1 = finit + ls->step * dgtest;
179 if (ls->bounded) ftest2 = finit + ls->step * dgtest * ls->ftol;
180
181 /* Convergence testing */
182 if ((*f - ftest1 <= PETSC_SMALL * PetscAbsReal(finit)) && (PetscAbsReal(dg) + ls->gtol * dginit <= 0.0)) {
183 PetscCall(PetscInfo(ls, "Line search success: Sufficient decrease and directional deriv conditions hold\n"));
184 ls->reason = TAOLINESEARCH_SUCCESS;
185 break;
186 }
187
188 /* Check Armijo if beyond the first breakpoint */
189 if (ls->bounded && *f <= ftest2 && ls->step >= bstepmin2) {
190 PetscCall(PetscInfo(ls, "Line search success: Sufficient decrease.\n"));
191 ls->reason = TAOLINESEARCH_SUCCESS;
192 break;
193 }
194
195 /* Checks for bad cases */
196 if ((mt->bracket && (ls->step <= ls->stepmin || ls->step >= ls->stepmax)) || !mt->infoc) {
197 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"));
198 ls->reason = TAOLINESEARCH_HALTED_OTHER;
199 break;
200 }
201 if (ls->step == ls->stepmax && *f <= ftest1 && dg <= dgtest) {
202 PetscCall(PetscInfo(ls, "Step is at the upper bound, stepmax (%g)\n", (double)ls->stepmax));
203 ls->reason = TAOLINESEARCH_HALTED_UPPERBOUND;
204 break;
205 }
206 if (ls->step == ls->stepmin && *f >= ftest1 && dg >= dgtest) {
207 PetscCall(PetscInfo(ls, "Step is at the lower bound, stepmin (%g)\n", (double)ls->stepmin));
208 ls->reason = TAOLINESEARCH_HALTED_LOWERBOUND;
209 break;
210 }
211 if (mt->bracket && (ls->stepmax - ls->stepmin <= ls->rtol * ls->stepmax)) {
212 PetscCall(PetscInfo(ls, "Relative width of interval of uncertainty is at most rtol (%g)\n", (double)ls->rtol));
213 ls->reason = TAOLINESEARCH_HALTED_RTOL;
214 break;
215 }
216
217 /* In the first stage, we seek a step for which the modified function
218 has a nonpositive value and nonnegative derivative */
219 if (stage1 && *f <= ftest1 && dg >= dginit * PetscMin(ls->ftol, ls->gtol)) stage1 = 0;
220
221 /* A modified function is used to predict the step only if we
222 have not obtained a step for which the modified function has a
223 nonpositive function value and nonnegative derivative, and if a
224 lower function value has been obtained but the decrease is not
225 sufficient */
226
227 if (stage1 && *f <= fx && *f > ftest1) {
228 fm = *f - ls->step * dgtest; /* Define modified function */
229 fxm = fx - stx * dgtest; /* and derivatives */
230 fym = fy - sty * dgtest;
231 dgm = dg - dgtest;
232 dgxm = dgx - dgtest;
233 dgym = dgy - dgtest;
234
235 /* if (dgxm * (ls->step - stx) >= 0.0) */
236 /* Update the interval of uncertainty and compute the new step */
237 PetscCall(Tao_mcstep(ls, &stx, &fxm, &dgxm, &sty, &fym, &dgym, &ls->step, &fm, &dgm));
238
239 fx = fxm + stx * dgtest; /* Reset the function and */
240 fy = fym + sty * dgtest; /* gradient values */
241 dgx = dgxm + dgtest;
242 dgy = dgym + dgtest;
243 } else {
244 /* Update the interval of uncertainty and compute the new step */
245 PetscCall(Tao_mcstep(ls, &stx, &fx, &dgx, &sty, &fy, &dgy, &ls->step, f, &dg));
246 }
247
248 /* Force a sufficient decrease in the interval of uncertainty */
249 if (mt->bracket) {
250 if (PetscAbsReal(sty - stx) >= 0.66 * width1) ls->step = stx + 0.5 * (sty - stx);
251 width1 = width;
252 width = PetscAbsReal(sty - stx);
253 }
254 }
255 if (ls->nfeval + ls->nfgeval > ls->max_funcs) {
256 PetscCall(PetscInfo(ls, "Number of line search function evals (%" PetscInt_FMT ") > maximum (%" PetscInt_FMT ")\n", ls->nfeval + ls->nfgeval, ls->max_funcs));
257 ls->reason = TAOLINESEARCH_HALTED_MAXFCN;
258 }
259 ls->stepmax = ostepmax;
260 ls->stepmin = ostepmin;
261
262 /* Finish computations */
263 PetscCall(PetscInfo(ls, "%" PetscInt_FMT " function evals in line search, step = %g\n", ls->nfeval + ls->nfgeval, (double)ls->step));
264
265 /* Set new solution vector and compute gradient if needed */
266 PetscCall(VecCopy(mt->work, x));
267 if (!g_computed) PetscCall(TaoLineSearchComputeGradient(ls, x, g));
268 PetscFunctionReturn(PETSC_SUCCESS);
269 }
270
271 /*MC
272 TAOLINESEARCHMT - More-Thuente line-search type with cubic interpolation that satisfies both the sufficient decrease and
273 curvature conditions. This method can take step lengths greater than 1, {cite}`more:92`
274
275 Options Database Key:
276 . -tao_ls_type more-thuente - use this line search type
277
278 Level: developer
279
280 .seealso: `TaoLineSearchCreate()`, `TaoLineSearchSetType()`, `TaoLineSearchApply()`
281 M*/
TaoLineSearchCreate_MT(TaoLineSearch ls)282 PETSC_EXTERN PetscErrorCode TaoLineSearchCreate_MT(TaoLineSearch ls)
283 {
284 TaoLineSearch_MT *ctx;
285
286 PetscFunctionBegin;
287 PetscValidHeaderSpecific(ls, TAOLINESEARCH_CLASSID, 1);
288 PetscCall(PetscNew(&ctx));
289 ctx->bracket = 0;
290 ctx->infoc = 1;
291 ls->data = (void *)ctx;
292 ls->initstep = 1.0;
293 ls->ops->setup = NULL;
294 ls->ops->reset = NULL;
295 ls->ops->apply = TaoLineSearchApply_MT;
296 ls->ops->destroy = TaoLineSearchDestroy_MT;
297 ls->ops->monitor = TaoLineSearchMonitor_MT;
298 PetscFunctionReturn(PETSC_SUCCESS);
299 }
300
301 /*
302 The subroutine mcstep is taken from the work of Jorge Nocedal.
303 this is a variant of More' and Thuente's routine.
304
305 subroutine mcstep
306
307 the purpose of mcstep is to compute a safeguarded step for
308 a linesearch and to update an interval of uncertainty for
309 a minimizer of the function.
310
311 the parameter stx contains the step with the least function
312 value. the parameter stp contains the current step. it is
313 assumed that the derivative at stx is negative in the
314 direction of the step. if bracket is set true then a
315 minimizer has been bracketed in an interval of uncertainty
316 with endpoints stx and sty.
317
318 the subroutine statement is
319
320 subroutine mcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,bracket,
321 stpmin,stpmax,info)
322
323 where
324
325 stx, fx, and dx are variables which specify the step,
326 the function, and the derivative at the best step obtained
327 so far. The derivative must be negative in the direction
328 of the step, that is, dx and stp-stx must have opposite
329 signs. On output these parameters are updated appropriately.
330
331 sty, fy, and dy are variables which specify the step,
332 the function, and the derivative at the other endpoint of
333 the interval of uncertainty. On output these parameters are
334 updated appropriately.
335
336 stp, fp, and dp are variables which specify the step,
337 the function, and the derivative at the current step.
338 If bracket is set true then on input stp must be
339 between stx and sty. On output stp is set to the new step.
340
341 bracket is a logical variable which specifies if a minimizer
342 has been bracketed. If the minimizer has not been bracketed
343 then on input bracket must be set false. If the minimizer
344 is bracketed then on output bracket is set true.
345
346 stpmin and stpmax are input variables which specify lower
347 and upper bounds for the step.
348
349 info is an integer output variable set as follows:
350 if info = 1,2,3,4,5, then the step has been computed
351 according to one of the five cases below. otherwise
352 info = 0, and this indicates improper input parameters.
353
354 subprograms called
355
356 fortran-supplied ... abs,max,min,sqrt
357
358 argonne national laboratory. minpack project. june 1983
359 jorge j. more', david j. thuente
360
361 */
362
Tao_mcstep(TaoLineSearch ls,PetscReal * stx,PetscReal * fx,PetscReal * dx,PetscReal * sty,PetscReal * fy,PetscReal * dy,PetscReal * stp,PetscReal * fp,PetscReal * dp)363 static PetscErrorCode Tao_mcstep(TaoLineSearch ls, PetscReal *stx, PetscReal *fx, PetscReal *dx, PetscReal *sty, PetscReal *fy, PetscReal *dy, PetscReal *stp, PetscReal *fp, PetscReal *dp)
364 {
365 TaoLineSearch_MT *mtP = (TaoLineSearch_MT *)ls->data;
366 PetscReal gamma1, p, q, r, s, sgnd, stpc, stpf, stpq, theta;
367 PetscInt bound;
368
369 PetscFunctionBegin;
370 /* Check the input parameters for errors */
371 mtP->infoc = 0;
372 PetscCheck(!mtP->bracket || (*stp > PetscMin(*stx, *sty) && *stp < PetscMax(*stx, *sty)), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "bad stp in bracket");
373 PetscCheck(*dx * (*stp - *stx) < 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "dx * (stp-stx) >= 0.0");
374 PetscCheck(ls->stepmax >= ls->stepmin, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "stepmax > stepmin");
375
376 /* Determine if the derivatives have opposite sign */
377 sgnd = *dp * (*dx / PetscAbsReal(*dx));
378
379 if (*fp > *fx) {
380 /* Case 1: a higher function value.
381 The minimum is bracketed. If the cubic step is closer
382 to stx than the quadratic step, the cubic step is taken,
383 else the average of the cubic and quadratic steps is taken. */
384
385 mtP->infoc = 1;
386 bound = 1;
387 theta = 3 * (*fx - *fp) / (*stp - *stx) + *dx + *dp;
388 s = PetscMax(PetscAbsReal(theta), PetscAbsReal(*dx));
389 s = PetscMax(s, PetscAbsReal(*dp));
390 gamma1 = s * PetscSqrtScalar(PetscPowScalar(theta / s, 2.0) - (*dx / s) * (*dp / s));
391 if (*stp < *stx) gamma1 = -gamma1;
392 /* Can p be 0? Check */
393 p = (gamma1 - *dx) + theta;
394 q = ((gamma1 - *dx) + gamma1) + *dp;
395 r = p / q;
396 stpc = *stx + r * (*stp - *stx);
397 stpq = *stx + ((*dx / ((*fx - *fp) / (*stp - *stx) + *dx)) * 0.5) * (*stp - *stx);
398
399 if (PetscAbsReal(stpc - *stx) < PetscAbsReal(stpq - *stx)) stpf = stpc;
400 else stpf = stpc + 0.5 * (stpq - stpc);
401 mtP->bracket = 1;
402 } else if (sgnd < 0.0) {
403 /* Case 2: A lower function value and derivatives of
404 opposite sign. The minimum is bracketed. If the cubic
405 step is closer to stx than the quadratic (secant) step,
406 the cubic step is taken, else the quadratic step is taken. */
407
408 mtP->infoc = 2;
409 bound = 0;
410 theta = 3 * (*fx - *fp) / (*stp - *stx) + *dx + *dp;
411 s = PetscMax(PetscAbsReal(theta), PetscAbsReal(*dx));
412 s = PetscMax(s, PetscAbsReal(*dp));
413 gamma1 = s * PetscSqrtScalar(PetscPowScalar(theta / s, 2.0) - (*dx / s) * (*dp / s));
414 if (*stp > *stx) gamma1 = -gamma1;
415 p = (gamma1 - *dp) + theta;
416 q = ((gamma1 - *dp) + gamma1) + *dx;
417 r = p / q;
418 stpc = *stp + r * (*stx - *stp);
419 stpq = *stp + (*dp / (*dp - *dx)) * (*stx - *stp);
420
421 if (PetscAbsReal(stpc - *stp) > PetscAbsReal(stpq - *stp)) stpf = stpc;
422 else stpf = stpq;
423 mtP->bracket = 1;
424 } else if (PetscAbsReal(*dp) < PetscAbsReal(*dx)) {
425 /* Case 3: A lower function value, derivatives of the
426 same sign, and the magnitude of the derivative decreases.
427 The cubic step is only used if the cubic tends to infinity
428 in the direction of the step or if the minimum of the cubic
429 is beyond stp. Otherwise the cubic step is defined to be
430 either stepmin or stepmax. The quadratic (secant) step is also
431 computed and if the minimum is bracketed then the step
432 closest to stx is taken, else the step farthest away is taken. */
433
434 mtP->infoc = 3;
435 bound = 1;
436 theta = 3 * (*fx - *fp) / (*stp - *stx) + *dx + *dp;
437 s = PetscMax(PetscAbsReal(theta), PetscAbsReal(*dx));
438 s = PetscMax(s, PetscAbsReal(*dp));
439
440 /* The case gamma1 = 0 only arises if the cubic does not tend
441 to infinity in the direction of the step. */
442 gamma1 = s * PetscSqrtScalar(PetscMax(0.0, PetscPowScalar(theta / s, 2.0) - (*dx / s) * (*dp / s)));
443 if (*stp > *stx) gamma1 = -gamma1;
444 p = (gamma1 - *dp) + theta;
445 q = (gamma1 + (*dx - *dp)) + gamma1;
446 r = p / q;
447 if (r < 0.0 && gamma1 != 0.0) stpc = *stp + r * (*stx - *stp);
448 else if (*stp > *stx) stpc = ls->stepmax;
449 else stpc = ls->stepmin;
450 stpq = *stp + (*dp / (*dp - *dx)) * (*stx - *stp);
451
452 if (mtP->bracket) {
453 if (PetscAbsReal(*stp - stpc) < PetscAbsReal(*stp - stpq)) stpf = stpc;
454 else stpf = stpq;
455 } else {
456 if (PetscAbsReal(*stp - stpc) > PetscAbsReal(*stp - stpq)) stpf = stpc;
457 else stpf = stpq;
458 }
459 } else {
460 /* Case 4: A lower function value, derivatives of the
461 same sign, and the magnitude of the derivative does
462 not decrease. If the minimum is not bracketed, the step
463 is either stpmin or stpmax, else the cubic step is taken. */
464
465 mtP->infoc = 4;
466 bound = 0;
467 if (mtP->bracket) {
468 theta = 3 * (*fp - *fy) / (*sty - *stp) + *dy + *dp;
469 s = PetscMax(PetscAbsReal(theta), PetscAbsReal(*dy));
470 s = PetscMax(s, PetscAbsReal(*dp));
471 gamma1 = s * PetscSqrtScalar(PetscPowScalar(theta / s, 2.0) - (*dy / s) * (*dp / s));
472 if (*stp > *sty) gamma1 = -gamma1;
473 p = (gamma1 - *dp) + theta;
474 q = ((gamma1 - *dp) + gamma1) + *dy;
475 r = p / q;
476 stpc = *stp + r * (*sty - *stp);
477 stpf = stpc;
478 } else if (*stp > *stx) {
479 stpf = ls->stepmax;
480 } else {
481 stpf = ls->stepmin;
482 }
483 }
484
485 /* Update the interval of uncertainty. This update does not
486 depend on the new step or the case analysis above. */
487
488 if (*fp > *fx) {
489 *sty = *stp;
490 *fy = *fp;
491 *dy = *dp;
492 } else {
493 if (sgnd < 0.0) {
494 *sty = *stx;
495 *fy = *fx;
496 *dy = *dx;
497 }
498 *stx = *stp;
499 *fx = *fp;
500 *dx = *dp;
501 }
502
503 /* Compute the new step and safeguard it. */
504 stpf = PetscMin(ls->stepmax, stpf);
505 stpf = PetscMax(ls->stepmin, stpf);
506 *stp = stpf;
507 if (mtP->bracket && bound) {
508 if (*sty > *stx) *stp = PetscMin(*stx + 0.66 * (*sty - *stx), *stp);
509 else *stp = PetscMax(*stx + 0.66 * (*sty - *stx), *stp);
510 }
511 PetscFunctionReturn(PETSC_SUCCESS);
512 }
513