1 : //------------------------------------------------------------------------
2 : // Copyright (C) 1993:
3 : // J.C. Meza
4 : // Sandia National Laboratories
5 : // meza@california.sandia.gov
6 : //------------------------------------------------------------------------
7 :
8 : #include "Opt.h"
9 : using NEWMAT::ColumnVector;
10 :
11 : namespace OPTPP {
12 :
13 : int mcsrch(NLP1* nlp, ColumnVector& s, ostream *fout, double *stp,
14 : int itnmax, double ftol, double xtol, double gtol, double stpmax,
15 0 : double stpmin)
16 : {
17 : /****************************************************************************
18 : * subroutine mcsrch
19 : * Purpose
20 : * find a step which satisfies
21 : * a sufficient decrease condition and a curvature condition.
22 : *
23 : * at each stage the subroutine updates an interval of
24 : * uncertainty with endpoints stx and sty. the interval of
25 : * uncertainty is initially chosen so that it contains a
26 : * minimizer of the modified function
27 : * f(x+stp*s) - f(x) - ftol*stp*(gradf(x)'s).
28 : *
29 : * if a step is obtained for which the modified function
30 : * has a nonpositive function value and nonnegative derivative,
31 : * then the interval of uncertainty is chosen so that it
32 : * contains a minimizer of f(x+stp*s).
33 : * the algorithm is designed to find a step which satisfies
34 : * the sufficient decrease condition
35 : * f(x+stp*s) .le. f(x) + ftol*stp*(gradf(x)'s),
36 : * and the curvature condition
37 : * abs(gradf(x+stp*s)'s)) .le. gtol*abs(gradf(x)'s).
38 : * if ftol is less than gtol and if, for example, the function
39 : * is bounded below, then there is always a step which satisfies
40 : * both conditions. if no step can be found which satisfies both
41 : * conditions, then the algorithm usually stops when rounding
42 : * errors prevent further progress. in this case stp only
43 : * satisfies the sufficient decrease condition.
44 : *
45 : * Parameters
46 : * n is a positive int input variable set to the number
47 : * of variables.
48 : * x is an array of length n. on input it must contain the
49 : * base point for the line search. on output it contains
50 : * x + stp*s.
51 : * f is a variable. on input it must contain the value of f
52 : * at x. on output it contains the value of f at x + stp*s.
53 : * g is an array of length n. on input it must contain the
54 : * gradient of f at x. on output it contains the gradient
55 : * of f at x + stp*s.
56 : * s is an input array of length n which specifies the
57 : * search direction.
58 : * stp is a nonnegative variable. on input stp contains an
59 : * initial estimate of a satisfactory step. on output
60 : * stp contains the final estimate.
61 : * ftol and gtol are nonnegative input variables.
62 : * termination occurs when the sufficient decrease
63 : * condition and the directional derivative condition are
64 : * satisfied.
65 : * ftol should be smaller than 5.e-1
66 : * suggested value = 1.e-4 for newton methods
67 : * = 1.e-1 for more exact line searches
68 : * Default Value = 1.e-4
69 : *
70 : * gtol should be greater than 1.e-4
71 : * Default Value = 0.9
72 : * xtol is a nonnegative input variable. termination occurs
73 : * when the relative width of the interval of uncertainty
74 : * is at most xtol.
75 : * Default Value = 2.2e-16
76 : * stpmin and stpmax are nonnegative input variables which
77 : * specify lower and upper bounds for the step. (in this reverse
78 : * communication implementatin they are defined in a common
79 : * statement).
80 : * stpmin Default Value = 1.e-9
81 : * stpmax Default Value = 1.e3
82 : * info is an integer output variable set as follows:
83 : * info =-1 improper input parameters.
84 : * info = 1 the sufficient decrease condition and the
85 : * directional derivative condition hold.
86 : * info =-2 relative width of the interval of uncertainty
87 : * is at most xtol.
88 : * info =-4 the step is at the lower bound stpmin.
89 : * info =-5 the step is at the upper bound stpmax.
90 : * info =-6 rounding errors prevent further progress.
91 : * there may not be a step which satisfies the
92 : * sufficient decrease and curvature conditions.
93 : * tolerances may be too small.
94 : * subprograms called
95 : * mcstep
96 : * fortran-supplied...abs,max,min
97 : * argonne national laboratory. minpack project. june 1983
98 : * jorge j. more', david j. thuente
99 : *
100 : *
101 : * Recoded in C++ by Juan Meza December 1992
102 : *
103 : *
104 : *****************************************************************************/
105 :
106 : /* initialized data */
107 :
108 : static double half = .5;
109 : static double p66 = .66;
110 : static double xtrapf = 4.;
111 : static double zero = 0.;
112 :
113 : /* local variables */
114 : static double dgxm, dgym;
115 : static int j, info, infoc;
116 : static double finit, width, stmin, stmax;
117 : static bool stage1;
118 : static double width1, ftest1, dg, fm, fx, fy;
119 : static bool brackt;
120 : static double dginit, dgtest;
121 : static double dgm, dgx, dgy, fxm, fym, stx, sty;
122 :
123 : int siter;
124 : //int maxiter = itnmax;
125 0 : int maxiter = 10;
126 0 : int n = nlp->getDim();
127 :
128 : double fvalue;
129 0 : ColumnVector work(n), xc(n), grad(n);
130 :
131 0 : infoc = 1;
132 :
133 : /* check the input parameters for errors. */
134 :
135 0 : if (n <= 0 || *stp <= zero || ftol < zero || gtol < zero ||
136 : xtol < zero || stpmin < zero || stpmax < stpmin ) {
137 0 : infoc = -1;
138 0 : return infoc;
139 : }
140 :
141 : /* compute the initial gradient in the search direction */
142 : /* and check that s is a descent direction. */
143 :
144 0 : dginit = zero;
145 0 : grad = nlp->getGrad();
146 0 : for (j = 1; j <= n; ++j) {
147 0 : dginit += grad(j) * s(j);
148 : }
149 :
150 0 : if (dginit >= zero) {
151 0 : cout << "\nmcsrch: Initial search direction not a descent direction\n";
152 0 : return -1;
153 : }
154 :
155 : /* initialize local variables. */
156 :
157 0 : brackt = false;
158 0 : stage1 = true;
159 0 : finit = nlp->getF();
160 0 : dgtest = ftol * dginit;
161 0 : width = stpmax - stpmin;
162 0 : width1 = width / half;
163 :
164 0 : work = nlp->getXc();
165 :
166 : /* the variables stx, fx, dgx contain the values of the step, */
167 : /* function, and directional derivative at the best step. */
168 : /* the variables sty, fy, dgy contain the value of the step, */
169 : /* function, and derivative at the other endpoint of */
170 : /* the interval of uncertainty. */
171 : /* the variables stp, f, dg contain the values of the step, */
172 : /* function, and derivative at the current step. */
173 :
174 0 : stx = zero;
175 0 : fx = finit;
176 0 : dgx = dginit;
177 0 : sty = zero;
178 0 : fy = finit;
179 0 : dgy = dginit;
180 :
181 0 : siter = 0;
182 :
183 : /* start of iteration. */
184 :
185 0 : while (siter < maxiter) {
186 0 : siter++;
187 :
188 : /* set the minimum and maximum steps to correspond */
189 : /* to the present interval of uncertainty. */
190 :
191 0 : if (brackt) {
192 0 : stmin = min(stx,sty);
193 0 : stmax = max(stx,sty);
194 : } else {
195 0 : stmin = stx;
196 0 : stmax = *stp + xtrapf * (*stp - stx);
197 : }
198 :
199 : /* force the step to be within the bounds stpmax and stpmin. */
200 :
201 0 : *stp = max(*stp,stpmin);
202 0 : *stp = min(*stp,stpmax);
203 :
204 : /* if an unusual termination is to occur then let */
205 : /* stp be the lowest point obtained so far. */
206 :
207 0 : if (brackt && (*stp <= stmin || *stp >= stmax) ||
208 : infoc == 0 || brackt && stmax - stmin <= xtol * stmax) {
209 0 : *stp = stx;
210 : }
211 :
212 : /* evaluate the function and gradient at stp */
213 : /* and compute the directional derivative. */
214 :
215 0 : xc = work + s*(*stp);
216 0 : nlp->setX(xc);
217 0 : nlp->eval();
218 0 : fvalue = nlp->getF();
219 0 : grad = nlp->getGrad();
220 :
221 0 : info = 0;
222 0 : dg = zero;
223 0 : for (j = 1; j <= n; ++j) {
224 0 : dg += grad(j) * s(j);
225 : }
226 0 : ftest1 = finit + *stp * dgtest;
227 :
228 : /* test for convergence. */
229 :
230 0 : if (brackt && (*stp <= stmin || *stp >= stmax) || infoc == 0) {
231 0 : info = -6; // CPJW 12/10/2003 original stmt info = 6
232 : }
233 0 : if (*stp == stpmax && fvalue <= ftest1 && dg <= dgtest) {
234 0 : info = -5; // CPJW 12/10/2003 original stmt info = 5
235 : }
236 0 : if (*stp == stpmin && (fvalue > ftest1 || dg >= dgtest)) {
237 0 : info = -4; // CPJW 12/10/2003 original stmt info = 4
238 : }
239 0 : if (brackt && stmax - stmin <= xtol * stmax) {
240 0 : info = -2; // CPJW 12/10/2003 original stmt info = 2
241 : }
242 0 : if (fvalue <= ftest1 && fabs(dg) <= gtol * (-dginit)) {
243 0 : info = 1;
244 : }
245 :
246 : /* check for termination. */
247 :
248 0 : if (info < 0) {
249 0 : return(-1);
250 : }
251 0 : if (info != 0) {
252 0 : if (siter == 1) return(Newton_Step);
253 0 : else return(Backtrack_Step);
254 : }
255 :
256 : /* in the first stage we seek a step for which the modified */
257 : /* function has a nonpositive value and nonnegative derivative. */
258 :
259 0 : if (stage1 && fvalue <= ftest1 && dg >= min(ftol,gtol) * dginit) {
260 0 : stage1 = false;
261 : }
262 :
263 : /* a modified function is used to predict the step only if */
264 : /* we have not obtained a step for which the modified */
265 : /* function has a nonpositive function value and nonnegative */
266 : /* derivative, and if a lower function value has been */
267 : /* obtained but the decrease is not sufficient. */
268 :
269 0 : if (stage1 && fvalue <= fx && fvalue > ftest1) {
270 :
271 : /* define the modified function and derivative values. */
272 :
273 0 : fm = fvalue - *stp * dgtest;
274 0 : fxm = fx - stx * dgtest;
275 0 : fym = fy - sty * dgtest;
276 0 : dgm = dg - dgtest;
277 0 : dgxm = dgx - dgtest;
278 0 : dgym = dgy - dgtest;
279 :
280 : /* call cstep to update the interval of uncertainty */
281 : /* and to compute the new step. */
282 :
283 : mcstep(&stx, &fxm, &dgxm, &sty, &fym, &dgym, stp, fm,
284 0 : dgm, &brackt, stmin, stmax, &infoc);
285 :
286 : /* reset the function and gradient values for f. */
287 :
288 0 : fx = fxm + stx * dgtest;
289 0 : fy = fym + sty * dgtest;
290 0 : dgx = dgxm + dgtest;
291 0 : dgy = dgym + dgtest;
292 : } else {
293 :
294 : /* call mcstep to update the interval of uncertainty */
295 : /* and to compute the new step. */
296 :
297 : mcstep(&stx, &fx, &dgx, &sty, &fy, &dgy, stp, fvalue,
298 0 : dg, &brackt, stmin, stmax, &infoc);
299 : }
300 :
301 : /* force a sufficient decrease in the size of the */
302 : /* interval of uncertainty. */
303 :
304 0 : if (brackt) {
305 0 : if (fabs(sty - stx) >= p66 * width1) {
306 0 : *stp = stx + (sty - stx)/2.;
307 : }
308 0 : width1 = width;
309 0 : width = fabs(sty - stx);
310 : }
311 :
312 : } /* end of iteration. */
313 :
314 0 : return (-1); // too many iterations;
315 :
316 : } /* last line of subroutine mcsrch. */
317 : /****************************************************************************/
318 : int mcstep(double *stx, double *fx, double *dx, double *sty, double *fy,
319 : double *dy, double *stp, double fp, double dp, bool *brackt,
320 0 : double stpmin, double stpmax, int *info)
321 : {
322 : /* system generated locals */
323 : double d_1, d_2, d_3;
324 :
325 : /* local variables */
326 : static double sgnd, stpc, stpf, stpq, p, q, gamma, r, s, theta;
327 : static bool bound;
328 :
329 :
330 : /******************************************************************************
331 : * subroutine mcstep
332 :
333 : * the purpose of mcstep is to compute a safeguarded step for
334 : * a linesearch and to update an interval of uncertainty for
335 : * a minimizer of the function.
336 :
337 : * the parameter stx contains the step with the least function
338 : * value. the parameter stp contains the current step. it is
339 : * assumed that the derivative at stx is negative in the
340 : * direction of the step. if brackt is set true then a
341 : * minimizer has been bracketed in an interval of uncertainty
342 : * with endpoints stx and sty.
343 :
344 : * the subroutine statement is
345 :
346 : * subroutine mcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,brackt,
347 : * stpmin,stpmax,info)
348 :
349 : * where
350 :
351 : * stx, fx, and dx are variables which specify the step,
352 : * the function, and the derivative at the best step obtained
353 : * so far. the derivative must be negative in the direction
354 : * of the step, that is, dx and stp-stx must have opposite
355 : * signs. on output these parameters are updated appropriately.
356 :
357 :
358 : * sty, fy, and dy are variables which specify the step,
359 : * the function, and the derivative at the other endpoint of
360 : * the interval of uncertainty. on output these parameters are
361 : * updated appropriately.
362 :
363 : * stp, fp, and dp are variables which specify the step,
364 : * the function, and the derivative at the current step.
365 : * if brackt is set true then on input stp must be
366 : * between stx and sty. on output stp is set to the new step.
367 :
368 : * brackt is a long int variable which specifies if a minimizer
369 : * has been bracketed. if the minimizer has not been bracketed
370 : * then on input brackt must be set false. if the minimizer
371 : * is bracketed then on output brackt is set true.
372 :
373 : * stpmin and stpmax are input variables which specify lower
374 : * and upper bounds for the step.
375 :
376 : * info is an int output variable set as follows:
377 : * if info = 1,2,3,4,5, then the step has been computed
378 : * according to one of the five cases below. otherwise
379 : * info = 0, and this indicates improper input parameters.
380 :
381 : * subprograms called
382 :
383 : * fortran-supplied ... abs,max,min,sqrt
384 :
385 : * argonne national laboratory. minpack project. june 1983
386 : * jorge j. more', david j. thuente
387 :
388 : ******************************************************************************/
389 :
390 0 : *info = 0;
391 :
392 : /* check the input parameters for errors. */
393 :
394 0 : if (*brackt &&
395 : (*stp <= min(*stx,*sty) || *stp >= max(*stx,*sty)) ||
396 : *dx * (*stp - *stx) >= 0. || stpmax < stpmin) {
397 0 : return 0;
398 : }
399 :
400 : /* determine if the derivatives have opposite sign. */
401 :
402 0 : sgnd = dp * (*dx / fabs(*dx));
403 :
404 : /* first case. a higher function value. */
405 : /* the minimum is bracketed. if the cubic step is closer */
406 : /* to stx than the quadratic step, the cubic step is taken, */
407 : /* else the average of the cubic and quadratic steps is taken. */
408 :
409 0 : if (fp > *fx) {
410 0 : *info = 1;
411 0 : bound = true;
412 0 : theta = (*fx - fp) * 3 / (*stp - *stx) + *dx + dp;
413 : /* computing max */
414 0 : d_1 = fabs(theta), d_2 = fabs(*dx), d_1 = max(d_2,d_1), d_2 = fabs(dp);
415 0 : s = max(d_2,d_1);
416 : /* computing 2nd power */
417 0 : d_1 = theta / s;
418 0 : gamma = s * sqrt(d_1 * d_1 - *dx / s * (dp / s));
419 0 : if (*stp < *stx) {
420 0 : gamma = -gamma;
421 : }
422 0 : p = gamma - *dx + theta;
423 0 : q = gamma - *dx + gamma + dp;
424 0 : r = p / q;
425 0 : stpc = *stx + r * (*stp - *stx);
426 : stpq = *stx + *dx / ((*fx - fp) / (*stp - *stx) + *dx) / 2 * (*stp -
427 0 : *stx);
428 0 : if ((d_1 = stpc - *stx, fabs(d_1)) < (d_2 = stpq - *stx, fabs(d_2))) {
429 0 : stpf = stpc;
430 : } else {
431 0 : stpf = stpc + (stpq - stpc) / 2;
432 : }
433 0 : *brackt = true;
434 :
435 : /* second case. a lower function value and derivatives of */
436 : /* opposite sign. the minimum is bracketed. if the cubic */
437 : /* step is closer to stx than the quadratic (secant) step, */
438 : /* the cubic step is taken, else the quadratic step is taken. */
439 :
440 0 : } else if (sgnd < 0.) {
441 0 : *info = 2;
442 0 : bound = false;
443 0 : theta = (*fx - fp) * 3 / (*stp - *stx) + *dx + dp;
444 : /* computing max */
445 0 : d_1 = fabs(theta), d_2 = fabs(*dx), d_1 = max(d_2,d_1), d_2 = fabs(dp);
446 0 : s = max(d_2,d_1);
447 : /* computing 2nd power */
448 0 : d_1 = theta / s;
449 0 : gamma = s * sqrt(d_1 * d_1 - *dx / s * (dp / s));
450 0 : if (*stp > *stx) {
451 0 : gamma = -gamma;
452 : }
453 0 : p = gamma - dp + theta;
454 0 : q = gamma - dp + gamma + *dx;
455 0 : r = p / q;
456 0 : stpc = *stp + r * (*stx - *stp);
457 0 : stpq = *stp + dp / (dp - *dx) * (*stx - *stp);
458 0 : if ((d_1 = stpc - *stp, fabs(d_1)) > (d_2 = stpq - *stp, fabs(d_2))) {
459 0 : stpf = stpc;
460 : } else {
461 0 : stpf = stpq;
462 : }
463 0 : *brackt = true;
464 :
465 : /* third case. a lower function value, derivatives of the */
466 : /* same sign, and the magnitude of the derivative decreases. */
467 : /* the cubic step is only used if the cubic tends to infinity */
468 : /* in the direction of the step or if the minimum of the cubic */
469 : /* is beyond stp. otherwise the cubic step is defined to be */
470 : /* either stpmin or stpmax. the quadratic (secant) step is also */
471 :
472 : /* computed and if the minimum is bracketed then the the step */
473 : /* closest to stx is taken, else the step farthest away is taken.
474 : */
475 :
476 0 : } else if (fabs(dp) < fabs(*dx)) {
477 0 : *info = 3;
478 0 : bound = true;
479 0 : theta = (*fx - fp) * 3 / (*stp - *stx) + *dx + dp;
480 : /* computing max */
481 0 : d_1 = fabs(theta), d_2 = fabs(*dx), d_1 = max(d_2,d_1), d_2 = fabs(dp);
482 0 : s = max(d_2,d_1);
483 :
484 : /* the case gamma = 0 only arises if the cubic does not tend */
485 : /* to infinity in the direction of the step. */
486 :
487 : /* computing max */
488 : /* computing 2nd power */
489 0 : d_3 = theta / s;
490 0 : d_1 = 0., d_2 = d_3 * d_3 - *dx / s * (dp / s);
491 0 : gamma = s * sqrt((max(d_2,d_1)));
492 0 : if (*stp > *stx) {
493 0 : gamma = -gamma;
494 : }
495 0 : p = gamma - dp + theta;
496 0 : q = gamma + (*dx - dp) + gamma;
497 0 : r = p / q;
498 0 : if (r < 0. && gamma != 0.) {
499 0 : stpc = *stp + r * (*stx - *stp);
500 0 : } else if (*stp > *stx) {
501 0 : stpc = stpmax;
502 : } else {
503 0 : stpc = stpmin;
504 : }
505 0 : stpq = *stp + dp / (dp - *dx) * (*stx - *stp);
506 0 : if (*brackt) {
507 0 : if ((d_1 = *stp - stpc, fabs(d_1)) < (d_2 = *stp - stpq, fabs(d_2)))
508 : {
509 0 : stpf = stpc;
510 : } else {
511 0 : stpf = stpq;
512 : }
513 : } else {
514 0 : if ((d_1 = *stp - stpc, fabs(d_1)) > (d_2 = *stp - stpq, fabs(d_2)))
515 : {
516 0 : stpf = stpc;
517 : } else {
518 0 : stpf = stpq;
519 : }
520 : }
521 :
522 : /* fourth case. a lower function value, derivatives of the */
523 : /* same sign, and the magnitude of the derivative does */
524 : /* not decrease. if the minimum is not bracketed, the step */
525 : /* is either stpmin or stpmax, else the cubic step is taken. */
526 :
527 : } else {
528 0 : *info = 4;
529 0 : bound = false;
530 0 : if (*brackt) {
531 0 : theta = (fp - *fy) * 3 / (*sty - *stp) + *dy + dp;
532 : /* computing max */
533 0 : d_1 = fabs(theta), d_2 = fabs(*dy), d_1 = max(d_2,d_1), d_2 = fabs(dp);
534 0 : s = max(d_2,d_1);
535 : /* computing 2nd power */
536 0 : d_1 = theta / s;
537 0 : gamma = s * sqrt(d_1 * d_1 - *dy / s * (dp / s));
538 0 : if (*stp > *sty) {
539 0 : gamma = -gamma;
540 : }
541 0 : p = gamma - dp + theta;
542 0 : q = gamma - dp + gamma + *dy;
543 0 : r = p / q;
544 0 : stpc = *stp + r * (*sty - *stp);
545 0 : stpf = stpc;
546 0 : } else if (*stp > *stx) {
547 0 : stpf = stpmax;
548 : } else {
549 0 : stpf = stpmin;
550 : }
551 : }
552 :
553 : /* update the interval of uncertainty. this update does not */
554 : /* depend on the new step or the case analysis above. */
555 :
556 0 : if (fp > *fx) {
557 0 : *sty = *stp;
558 0 : *fy = fp;
559 0 : *dy = dp;
560 : } else {
561 0 : if (sgnd < 0.) {
562 0 : *sty = *stx;
563 0 : *fy = *fx;
564 0 : *dy = *dx;
565 : }
566 0 : *stx = *stp;
567 0 : *fx = fp;
568 0 : *dx = dp;
569 : }
570 :
571 : /* compute the new step and safeguard it. */
572 :
573 0 : stpf = min(stpmax,stpf);
574 0 : stpf = max(stpmin,stpf);
575 0 : *stp = stpf;
576 0 : if (*brackt && bound) {
577 0 : if (*sty > *stx) {
578 : /* computing max */
579 0 : d_1 = *stx + (*sty - *stx) * (float).66;
580 0 : *stp = min(*stp,d_1);
581 : } else {
582 : /* computing max */
583 0 : d_1 = *stx + (*sty - *stx) * (float).66;
584 0 : *stp = max(*stp,d_1);
585 : }
586 : }
587 0 : return 0;
588 :
589 : }
590 :
591 156 : } // namespace OPTPP
|
