1 /* BOBYQA derivative-free optimization code by M. J. D. Powell.
2 Original Fortran code by Powell (2009). Converted via v2c,
3 cleaned up, and incorporated into NLopt by S. G. Johnson (2009).
12 typedef double (*bobyqa_func)(int n, const double *x, void *func_data);
14 #define MIN2(a,b) ((a) <= (b) ? (a) : (b))
15 #define MAX2(a,b) ((a) >= (b) ? (a) : (b))
16 #define IABS(x) ((x) < 0 ? -(x) : (x))
18 static void update_(int *n, int *npt, double *bmat,
19 double *zmat, int *ndim, double *vlag, double *beta,
20 double *denom, int *knew, double *w)
22 /* System generated locals */
23 int bmat_dim1, bmat_offset, zmat_dim1, zmat_offset, i__1, i__2;
24 double d__1, d__2, d__3;
27 int i__, j, k, jl, jp;
28 double one, tau, temp;
30 double zero, alpha, tempa, tempb, ztest;
33 /* The arrays BMAT and ZMAT are updated, as required by the new position */
34 /* of the interpolation point that has the index KNEW. The vector VLAG has */
35 /* N+NPT components, set on entry to the first NPT and last N components */
36 /* of the product Hw in equation (4.11) of the Powell (2006) paper on */
37 /* NEWUOA. Further, BETA is set on entry to the value of the parameter */
38 /* with that name, and DENOM is set to the denominator of the updating */
39 /* formula. Elements of ZMAT may be treated as zero if their moduli are */
40 /* at most ZTEST. The first NDIM elements of W are used for working space. */
42 /* Set some constants. */
44 /* Parameter adjustments */
46 zmat_offset = 1 + zmat_dim1;
49 bmat_offset = 1 + bmat_dim1;
60 for (k = 1; k <= i__1; ++k) {
62 for (j = 1; j <= i__2; ++j) {
65 d__2 = ztest, d__3 = (d__1 = zmat[k + j * zmat_dim1], fabs(d__1));
66 ztest = MAX2(d__2,d__3);
71 /* Apply the rotations that put zeros in the KNEW-th row of ZMAT. */
75 for (j = 2; j <= i__2; ++j) {
76 if ((d__1 = zmat[*knew + j * zmat_dim1], fabs(d__1)) > ztest) {
77 /* Computing 2nd power */
78 d__1 = zmat[*knew + zmat_dim1];
79 /* Computing 2nd power */
80 d__2 = zmat[*knew + j * zmat_dim1];
81 temp = sqrt(d__1 * d__1 + d__2 * d__2);
82 tempa = zmat[*knew + zmat_dim1] / temp;
83 tempb = zmat[*knew + j * zmat_dim1] / temp;
85 for (i__ = 1; i__ <= i__1; ++i__) {
86 temp = tempa * zmat[i__ + zmat_dim1] + tempb * zmat[i__ + j *
88 zmat[i__ + j * zmat_dim1] = tempa * zmat[i__ + j * zmat_dim1]
89 - tempb * zmat[i__ + zmat_dim1];
91 zmat[i__ + zmat_dim1] = temp;
94 zmat[*knew + j * zmat_dim1] = zero;
98 /* Put the first NPT components of the KNEW-th column of HLAG into W, */
99 /* and calculate the parameters of the updating formula. */
102 for (i__ = 1; i__ <= i__2; ++i__) {
103 w[i__] = zmat[*knew + zmat_dim1] * zmat[i__ + zmat_dim1];
110 /* Complete the updating of ZMAT. */
113 tempb = zmat[*knew + zmat_dim1] / temp;
116 for (i__ = 1; i__ <= i__2; ++i__) {
118 zmat[i__ + zmat_dim1] = tempa * zmat[i__ + zmat_dim1] - tempb * vlag[
122 /* Finally, update the matrix BMAT. */
125 for (j = 1; j <= i__2; ++j) {
127 w[jp] = bmat[*knew + j * bmat_dim1];
128 tempa = (alpha * vlag[jp] - tau * w[jp]) / *denom;
129 tempb = (-(*beta) * w[jp] - tau * vlag[jp]) / *denom;
131 for (i__ = 1; i__ <= i__1; ++i__) {
132 bmat[i__ + j * bmat_dim1] = bmat[i__ + j * bmat_dim1] + tempa *
133 vlag[i__] + tempb * w[i__];
135 bmat[jp + (i__ - *npt) * bmat_dim1] = bmat[i__ + j *
143 static nlopt_result rescue_(int *n, int *npt, const double *xl, const double *xu,
145 nlopt_stopping *stop,
146 bobyqa_func calfun, void *calfun_data,
147 double *xbase, double *xpt, double *fval, double *xopt, double *gopt,
148 double *hq, double *pq, double *bmat, double *zmat,
149 int *ndim, double *sl, double *su, /* int *nf, */
150 double *delta, int *kopt, double *vlag, double *
151 ptsaux, double *ptsid, double *w)
153 /* System generated locals */
154 int xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1,
155 zmat_offset, i__1, i__2, i__3;
156 double d__1, d__2, d__3, d__4;
158 /* Local variables */
160 int i__, j, k, ih, jp, ip, iq, np, iw;
165 double sum, diff, half, beta;
171 double zero, hdiag, fbase, sfrac, denom, vquad, sumpq;
172 double dsqmin, distsq, vlmxsq;
175 /* The arguments N, NPT, XL, XU, MAXFUN, XBASE, XPT, FVAL, XOPT, */
176 /* GOPT, HQ, PQ, BMAT, ZMAT, NDIM, SL and SU have the same meanings as */
177 /* the corresponding arguments of BOBYQB on the entry to RESCUE. */
178 /* NF is maintained as the number of calls of CALFUN so far, except that */
179 /* NF is set to -1 if the value of MAXFUN prevents further progress. */
180 /* KOPT is maintained so that FVAL(KOPT) is the least calculated function */
181 /* value. Its correct value must be given on entry. It is updated if a */
182 /* new least function value is found, but the corresponding changes to */
183 /* XOPT and GOPT have to be made later by the calling program. */
184 /* DELTA is the current trust region radius. */
185 /* VLAG is a working space vector that will be used for the values of the */
186 /* provisional Lagrange functions at each of the interpolation points. */
187 /* They are part of a product that requires VLAG to be of length NDIM. */
188 /* PTSAUX is also a working space array. For J=1,2,...,N, PTSAUX(1,J) and */
189 /* PTSAUX(2,J) specify the two positions of provisional interpolation */
190 /* points when a nonzero step is taken along e_J (the J-th coordinate */
191 /* direction) through XBASE+XOPT, as specified below. Usually these */
192 /* steps have length DELTA, but other lengths are chosen if necessary */
193 /* in order to satisfy the given bounds on the variables. */
194 /* PTSID is also a working space array. It has NPT components that denote */
195 /* provisional new positions of the original interpolation points, in */
196 /* case changes are needed to restore the linear independence of the */
197 /* interpolation conditions. The K-th point is a candidate for change */
198 /* if and only if PTSID(K) is nonzero. In this case let p and q be the */
199 /* int parts of PTSID(K) and (PTSID(K)-p) multiplied by N+1. If p */
200 /* and q are both positive, the step from XBASE+XOPT to the new K-th */
201 /* interpolation point is PTSAUX(1,p)*e_p + PTSAUX(1,q)*e_q. Otherwise */
202 /* the step is PTSAUX(1,p)*e_p or PTSAUX(2,q)*e_q in the cases q=0 or */
203 /* p=0, respectively. */
204 /* The first NDIM+NPT elements of the array W are used for working space. */
205 /* The final elements of BMAT and ZMAT are set in a well-conditioned way */
206 /* to the values that are appropriate for the new interpolation points. */
207 /* The elements of GOPT, HQ and PQ are also revised to the values that are */
208 /* appropriate to the final quadratic model. */
210 /* Set some constants. */
212 /* Parameter adjustments */
214 zmat_offset = 1 + zmat_dim1;
217 xpt_offset = 1 + xpt_dim1;
228 bmat_offset = 1 + bmat_dim1;
242 sfrac = half / (double) np;
245 /* Shift the interpolation points so that XOPT becomes the origin, and set */
246 /* the elements of ZMAT to zero. The value of SUMPQ is required in the */
247 /* updating of HQ below. The squares of the distances from XOPT to the */
248 /* other interpolation points are set at the end of W. Increments of WINC */
249 /* may be added later to these squares to balance the consideration of */
250 /* the choice of point that is going to become current. */
255 for (k = 1; k <= i__1; ++k) {
258 for (j = 1; j <= i__2; ++j) {
259 xpt[k + j * xpt_dim1] -= xopt[j];
261 /* Computing 2nd power */
262 d__1 = xpt[k + j * xpt_dim1];
263 distsq += d__1 * d__1;
266 w[*ndim + k] = distsq;
267 winc = MAX2(winc,distsq);
269 for (j = 1; j <= i__2; ++j) {
271 zmat[k + j * zmat_dim1] = zero;
275 /* Update HQ so that HQ and PQ define the second derivatives of the model */
276 /* after XBASE has been shifted to the trust region centre. */
280 for (j = 1; j <= i__2; ++j) {
281 w[j] = half * sumpq * xopt[j];
283 for (k = 1; k <= i__1; ++k) {
285 w[j] += pq[k] * xpt[k + j * xpt_dim1];
288 for (i__ = 1; i__ <= i__1; ++i__) {
291 hq[ih] = hq[ih] + w[i__] * xopt[j] + w[j] * xopt[i__];
295 /* Shift XBASE, SL, SU and XOPT. Set the elements of BMAT to zero, and */
296 /* also set the elements of PTSAUX. */
299 for (j = 1; j <= i__1; ++j) {
305 d__1 = *delta, d__2 = su[j];
306 ptsaux[(j << 1) + 1] = MIN2(d__1,d__2);
308 d__1 = -(*delta), d__2 = sl[j];
309 ptsaux[(j << 1) + 2] = MAX2(d__1,d__2);
310 if (ptsaux[(j << 1) + 1] + ptsaux[(j << 1) + 2] < zero) {
311 temp = ptsaux[(j << 1) + 1];
312 ptsaux[(j << 1) + 1] = ptsaux[(j << 1) + 2];
313 ptsaux[(j << 1) + 2] = temp;
315 if ((d__2 = ptsaux[(j << 1) + 2], fabs(d__2)) < half * (d__1 = ptsaux[(
316 j << 1) + 1], fabs(d__1))) {
317 ptsaux[(j << 1) + 2] = half * ptsaux[(j << 1) + 1];
320 for (i__ = 1; i__ <= i__2; ++i__) {
322 bmat[i__ + j * bmat_dim1] = zero;
327 /* Set the identifiers of the artificial interpolation points that are */
328 /* along a coordinate direction from XOPT, and set the corresponding */
329 /* nonzero elements of BMAT and ZMAT. */
333 for (j = 1; j <= i__2; ++j) {
336 ptsid[jp] = (double) j + sfrac;
338 ptsid[jpn] = (double) j / (double) np + sfrac;
339 temp = one / (ptsaux[(j << 1) + 1] - ptsaux[(j << 1) + 2]);
340 bmat[jp + j * bmat_dim1] = -temp + one / ptsaux[(j << 1) + 1];
341 bmat[jpn + j * bmat_dim1] = temp + one / ptsaux[(j << 1) + 2];
342 bmat[j * bmat_dim1 + 1] = -bmat[jp + j * bmat_dim1] - bmat[jpn +
344 zmat[j * zmat_dim1 + 1] = sqrt(2.) / (d__1 = ptsaux[(j << 1) + 1]
345 * ptsaux[(j << 1) + 2], fabs(d__1));
346 zmat[jp + j * zmat_dim1] = zmat[j * zmat_dim1 + 1] * ptsaux[(j <<
348 zmat[jpn + j * zmat_dim1] = -zmat[j * zmat_dim1 + 1] * ptsaux[(j
351 bmat[j * bmat_dim1 + 1] = -one / ptsaux[(j << 1) + 1];
352 bmat[jp + j * bmat_dim1] = one / ptsaux[(j << 1) + 1];
353 /* Computing 2nd power */
354 d__1 = ptsaux[(j << 1) + 1];
355 bmat[j + *npt + j * bmat_dim1] = -half * (d__1 * d__1);
360 /* Set any remaining identifiers with their nonzero elements of ZMAT. */
362 if (*npt >= *n + np) {
364 for (k = np << 1; k <= i__2; ++k) {
365 iw = (int) (((double) (k - np) - half) / (double) (*n)
367 ip = k - np - iw * *n;
372 ptsid[k] = (double) ip + (double) iq / (double) np +
374 temp = one / (ptsaux[(ip << 1) + 1] * ptsaux[(iq << 1) + 1]);
375 zmat[(k - np) * zmat_dim1 + 1] = temp;
376 zmat[ip + 1 + (k - np) * zmat_dim1] = -temp;
377 zmat[iq + 1 + (k - np) * zmat_dim1] = -temp;
379 zmat[k + (k - np) * zmat_dim1] = temp;
386 /* Reorder the provisional points in the way that exchanges PTSID(KOLD) */
387 /* with PTSID(KNEW). */
391 for (j = 1; j <= i__2; ++j) {
392 temp = bmat[kold + j * bmat_dim1];
393 bmat[kold + j * bmat_dim1] = bmat[knew + j * bmat_dim1];
395 bmat[knew + j * bmat_dim1] = temp;
398 for (j = 1; j <= i__2; ++j) {
399 temp = zmat[kold + j * zmat_dim1];
400 zmat[kold + j * zmat_dim1] = zmat[knew + j * zmat_dim1];
402 zmat[knew + j * zmat_dim1] = temp;
404 ptsid[kold] = ptsid[knew];
406 w[*ndim + knew] = zero;
410 vlag[kold] = vlag[knew];
413 /* Update the BMAT and ZMAT matrices so that the status of the KNEW-th */
414 /* interpolation point can be changed from provisional to original. The */
415 /* branch to label 350 occurs if all the original points are reinstated. */
416 /* The nonnegative values of W(NDIM+K) are required in the search below. */
418 update_(n, npt, &bmat[bmat_offset], &zmat[zmat_offset], ndim, &vlag[1]
419 , &beta, &denom, &knew, &w[1]);
424 for (k = 1; k <= i__2; ++k) {
426 w[*ndim + k] = (d__1 = w[*ndim + k], fabs(d__1));
430 /* Pick the index KNEW of an original interpolation point that has not */
431 /* yet replaced one of the provisional interpolation points, giving */
432 /* attention to the closeness to XOPT and to previous tries with KNEW. */
437 for (k = 1; k <= i__2; ++k) {
438 if (w[*ndim + k] > zero) {
439 if (dsqmin == zero || w[*ndim + k] < dsqmin) {
441 dsqmin = w[*ndim + k];
446 if (dsqmin == zero) {
450 /* Form the W-vector of the chosen original interpolation point. */
453 for (j = 1; j <= i__2; ++j) {
455 w[*npt + j] = xpt[knew + j * xpt_dim1];
458 for (k = 1; k <= i__2; ++k) {
461 } else if (ptsid[k] == zero) {
463 for (j = 1; j <= i__1; ++j) {
465 sum += w[*npt + j] * xpt[k + j * xpt_dim1];
470 sum = w[*npt + ip] * ptsaux[(ip << 1) + 1];
472 iq = (int) ((double) np * ptsid[k] - (double) (ip *
479 sum += w[*npt + iq] * ptsaux[iw + (iq << 1)];
483 w[k] = half * sum * sum;
486 /* Calculate VLAG and BETA for the required updating of the H matrix if */
487 /* XPT(KNEW,.) is reinstated in the set of interpolation points. */
490 for (k = 1; k <= i__2; ++k) {
493 for (j = 1; j <= i__1; ++j) {
495 sum += bmat[k + j * bmat_dim1] * w[*npt + j];
502 for (j = 1; j <= i__2; ++j) {
505 for (k = 1; k <= i__1; ++k) {
507 sum += zmat[k + j * zmat_dim1] * w[k];
511 for (k = 1; k <= i__1; ++k) {
513 vlag[k] += sum * zmat[k + j * zmat_dim1];
519 for (j = 1; j <= i__1; ++j) {
522 for (k = 1; k <= i__2; ++k) {
524 sum += bmat[k + j * bmat_dim1] * w[k];
529 for (ip = *npt + 1; ip <= i__2; ++ip) {
531 sum += bmat[ip + j * bmat_dim1] * w[ip];
536 /* Computing 2nd power */
537 d__1 = xpt[knew + j * xpt_dim1];
538 distsq += d__1 * d__1;
540 beta = half * distsq * distsq + beta - bsum;
543 /* KOLD is set to the index of the provisional interpolation point that is */
544 /* going to be deleted to make way for the KNEW-th original interpolation */
545 /* point. The choice of KOLD is governed by the avoidance of a small value */
546 /* of the denominator in the updating calculation of UPDATE. */
551 for (k = 1; k <= i__1; ++k) {
552 if (ptsid[k] != zero) {
555 for (j = 1; j <= i__2; ++j) {
557 /* Computing 2nd power */
558 d__1 = zmat[k + j * zmat_dim1];
559 hdiag += d__1 * d__1;
561 /* Computing 2nd power */
563 den = beta * hdiag + d__1 * d__1;
571 /* Computing 2nd power */
573 d__1 = vlmxsq, d__2 = d__3 * d__3;
574 vlmxsq = MAX2(d__1,d__2);
576 if (denom <= vlmxsq * .01) {
577 w[*ndim + knew] = -w[*ndim + knew] - winc;
582 /* When label 260 is reached, all the final positions of the interpolation */
583 /* points have been chosen although any changes have not been included yet */
584 /* in XPT. Also the final BMAT and ZMAT matrices are complete, but, apart */
585 /* from the shift of XBASE, the updating of the quadratic model remains to */
586 /* be done. The following cycle through the new interpolation points begins */
587 /* by putting the new point in XPT(KPT,.) and by setting PQ(KPT) to zero, */
588 /* except that a RETURN occurs if MAXFUN prohibits another value of F. */
592 for (kpt = 1; kpt <= i__1; ++kpt) {
593 if (ptsid[kpt] == zero) {
597 if (nlopt_stop_forced(stop)) return NLOPT_FORCED_STOP;
598 else if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
599 else if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
603 for (j = 1; j <= i__2; ++j) {
604 w[j] = xpt[kpt + j * xpt_dim1];
605 xpt[kpt + j * xpt_dim1] = zero;
606 temp = pq[kpt] * w[j];
608 for (i__ = 1; i__ <= i__3; ++i__) {
611 hq[ih] += temp * w[i__];
615 ip = (int) ptsid[kpt];
616 iq = (int) ((double) np * ptsid[kpt] - (double) (ip * np))
619 xp = ptsaux[(ip << 1) + 1];
620 xpt[kpt + ip * xpt_dim1] = xp;
623 xq = ptsaux[(iq << 1) + 1];
625 xq = ptsaux[(iq << 1) + 2];
627 xpt[kpt + iq * xpt_dim1] = xq;
630 /* Set VQUAD to the value of the current model at the new point. */
634 ihp = (ip + ip * ip) / 2;
635 vquad += xp * (gopt[ip] + half * xp * hq[ihp]);
638 ihq = (iq + iq * iq) / 2;
639 vquad += xq * (gopt[iq] + half * xq * hq[ihq]);
641 iw = MAX2(ihp,ihq) - (i__3 = ip - iq, IABS(i__3));
642 vquad += xp * xq * hq[iw];
646 for (k = 1; k <= i__3; ++k) {
649 temp += xp * xpt[k + ip * xpt_dim1];
652 temp += xq * xpt[k + iq * xpt_dim1];
655 vquad += half * pq[k] * temp * temp;
658 /* Calculate F at the new interpolation point, and set DIFF to the factor */
659 /* that is going to multiply the KPT-th Lagrange function when the model */
660 /* is updated to provide interpolation to the new function value. */
663 for (i__ = 1; i__ <= i__3; ++i__) {
666 d__3 = xl[i__], d__4 = xbase[i__] + xpt[kpt + i__ * xpt_dim1];
667 d__1 = MAX2(d__3,d__4), d__2 = xu[i__];
668 w[i__] = MIN2(d__1,d__2);
669 if (xpt[kpt + i__ * xpt_dim1] == sl[i__]) {
672 if (xpt[kpt + i__ * xpt_dim1] == su[i__]) {
679 f = calfun(*n, &w[1], calfun_data);
681 if (f < fval[*kopt]) {
684 if (nlopt_stop_forced(stop)) return NLOPT_FORCED_STOP;
685 else if (f < stop->minf_max) return NLOPT_MINF_MAX_REACHED;
686 else if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
687 else if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
691 /* Update the quadratic model. The RETURN from the subroutine occurs when */
692 /* all the new interpolation points are included in the model. */
695 for (i__ = 1; i__ <= i__3; ++i__) {
697 gopt[i__] += diff * bmat[kpt + i__ * bmat_dim1];
700 for (k = 1; k <= i__3; ++k) {
703 for (j = 1; j <= i__2; ++j) {
705 sum += zmat[k + j * zmat_dim1] * zmat[kpt + j * zmat_dim1];
708 if (ptsid[k] == zero) {
712 iq = (int) ((double) np * ptsid[k] - (double) (ip
714 ihq = (iq * iq + iq) / 2;
716 /* Computing 2nd power */
717 d__1 = ptsaux[(iq << 1) + 2];
718 hq[ihq] += temp * (d__1 * d__1);
720 ihp = (ip * ip + ip) / 2;
721 /* Computing 2nd power */
722 d__1 = ptsaux[(ip << 1) + 1];
723 hq[ihp] += temp * (d__1 * d__1);
725 /* Computing 2nd power */
726 d__1 = ptsaux[(iq << 1) + 1];
727 hq[ihq] += temp * (d__1 * d__1);
728 iw = MAX2(ihp,ihq) - (i__2 = iq - ip, IABS(i__2));
729 hq[iw] += temp * ptsaux[(ip << 1) + 1] * ptsaux[(iq <<
741 return NLOPT_SUCCESS;
744 static void altmov_(int *n, int *npt, double *xpt,
745 double *xopt, double *bmat, double *zmat, int *ndim,
746 double *sl, double *su, int *kopt, int *knew,
747 double *adelt, double *xnew, double *xalt, double *
748 alpha, double *cauchy, double *glag, double *hcol,
751 /* System generated locals */
752 int xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1,
753 zmat_offset, i__1, i__2;
754 double d__1, d__2, d__3, d__4;
756 /* Local variables */
758 double ha, gw, one, diff, half;
762 double vlag, subd, temp;
764 double step, zero, curv;
766 double scale, csave, tempa, tempb, tempd, const__, sumin, ggfree;
768 double dderiv, bigstp, predsq, presav, distsq, stpsav, wfixsq, wsqsav;
771 /* The arguments N, NPT, XPT, XOPT, BMAT, ZMAT, NDIM, SL and SU all have */
772 /* the same meanings as the corresponding arguments of BOBYQB. */
773 /* KOPT is the index of the optimal interpolation point. */
774 /* KNEW is the index of the interpolation point that is going to be moved. */
775 /* ADELT is the current trust region bound. */
776 /* XNEW will be set to a suitable new position for the interpolation point */
777 /* XPT(KNEW,.). Specifically, it satisfies the SL, SU and trust region */
778 /* bounds and it should provide a large denominator in the next call of */
779 /* UPDATE. The step XNEW-XOPT from XOPT is restricted to moves along the */
780 /* straight lines through XOPT and another interpolation point. */
781 /* XALT also provides a large value of the modulus of the KNEW-th Lagrange */
782 /* function subject to the constraints that have been mentioned, its main */
783 /* difference from XNEW being that XALT-XOPT is a constrained version of */
784 /* the Cauchy step within the trust region. An exception is that XALT is */
785 /* not calculated if all components of GLAG (see below) are zero. */
786 /* ALPHA will be set to the KNEW-th diagonal element of the H matrix. */
787 /* CAUCHY will be set to the square of the KNEW-th Lagrange function at */
788 /* the step XALT-XOPT from XOPT for the vector XALT that is returned, */
789 /* except that CAUCHY is set to zero if XALT is not calculated. */
790 /* GLAG is a working space vector of length N for the gradient of the */
791 /* KNEW-th Lagrange function at XOPT. */
792 /* HCOL is a working space vector of length NPT for the second derivative */
793 /* coefficients of the KNEW-th Lagrange function. */
794 /* W is a working space vector of length 2N that is going to hold the */
795 /* constrained Cauchy step from XOPT of the Lagrange function, followed */
796 /* by the downhill version of XALT when the uphill step is calculated. */
798 /* Set the first NPT components of W to the leading elements of the */
799 /* KNEW-th column of the H matrix. */
801 /* Parameter adjustments */
803 zmat_offset = 1 + zmat_dim1;
806 xpt_offset = 1 + xpt_dim1;
810 bmat_offset = 1 + bmat_dim1;
824 const__ = one + sqrt(2.);
826 for (k = 1; k <= i__1; ++k) {
830 i__1 = *npt - *n - 1;
831 for (j = 1; j <= i__1; ++j) {
832 temp = zmat[*knew + j * zmat_dim1];
834 for (k = 1; k <= i__2; ++k) {
836 hcol[k] += temp * zmat[k + j * zmat_dim1];
839 *alpha = hcol[*knew];
842 /* Calculate the gradient of the KNEW-th Lagrange function at XOPT. */
845 for (i__ = 1; i__ <= i__2; ++i__) {
847 glag[i__] = bmat[*knew + i__ * bmat_dim1];
850 for (k = 1; k <= i__2; ++k) {
853 for (j = 1; j <= i__1; ++j) {
855 temp += xpt[k + j * xpt_dim1] * xopt[j];
857 temp = hcol[k] * temp;
859 for (i__ = 1; i__ <= i__1; ++i__) {
861 glag[i__] += temp * xpt[k + i__ * xpt_dim1];
865 /* Search for a large denominator along the straight lines through XOPT */
866 /* and another interpolation point. SLBD and SUBD will be lower and upper */
867 /* bounds on the step along each of these lines in turn. PREDSQ will be */
868 /* set to the square of the predicted denominator for each line. PRESAV */
869 /* will be set to the largest admissible value of PREDSQ that occurs. */
873 for (k = 1; k <= i__1; ++k) {
880 for (i__ = 1; i__ <= i__2; ++i__) {
881 temp = xpt[k + i__ * xpt_dim1] - xopt[i__];
882 dderiv += glag[i__] * temp;
884 distsq += temp * temp;
886 subd = *adelt / sqrt(distsq);
890 sumin = MIN2(one,subd);
892 /* Revise SLBD and SUBD if necessary because of the bounds in SL and SU. */
895 for (i__ = 1; i__ <= i__2; ++i__) {
896 temp = xpt[k + i__ * xpt_dim1] - xopt[i__];
898 if (slbd * temp < sl[i__] - xopt[i__]) {
899 slbd = (sl[i__] - xopt[i__]) / temp;
902 if (subd * temp > su[i__] - xopt[i__]) {
904 d__1 = sumin, d__2 = (su[i__] - xopt[i__]) / temp;
905 subd = MAX2(d__1,d__2);
908 } else if (temp < zero) {
909 if (slbd * temp > su[i__] - xopt[i__]) {
910 slbd = (su[i__] - xopt[i__]) / temp;
913 if (subd * temp < sl[i__] - xopt[i__]) {
915 d__1 = sumin, d__2 = (sl[i__] - xopt[i__]) / temp;
916 subd = MAX2(d__1,d__2);
923 /* Seek a large modulus of the KNEW-th Lagrange function when the index */
924 /* of the other interpolation point on the line through XOPT is KNEW. */
929 vlag = slbd * (dderiv - slbd * diff);
931 temp = subd * (dderiv - subd * diff);
932 if (fabs(temp) > fabs(vlag)) {
937 tempd = half * dderiv;
938 tempa = tempd - diff * slbd;
939 tempb = tempd - diff * subd;
940 if (tempa * tempb < zero) {
941 temp = tempd * tempd / diff;
942 if (fabs(temp) > fabs(vlag)) {
949 /* Search along each of the other lines through XOPT and another point. */
953 vlag = slbd * (one - slbd);
955 temp = subd * (one - subd);
956 if (fabs(temp) > fabs(vlag)) {
962 if (fabs(vlag) < .25) {
971 /* Calculate PREDSQ for the current line search and maintain PRESAV. */
973 temp = step * (one - step) * distsq;
974 predsq = vlag * vlag * (vlag * vlag + ha * temp * temp);
975 if (predsq > presav) {
985 /* Construct XNEW in a way that satisfies the bound constraints exactly. */
988 for (i__ = 1; i__ <= i__1; ++i__) {
989 temp = xopt[i__] + stpsav * (xpt[ksav + i__ * xpt_dim1] - xopt[i__]);
994 d__1 = sl[i__], d__2 = MIN2(d__3,temp);
995 xnew[i__] = MAX2(d__1,d__2);
998 xnew[-ibdsav] = sl[-ibdsav];
1001 xnew[ibdsav] = su[ibdsav];
1004 /* Prepare for the iterative method that assembles the constrained Cauchy */
1005 /* step in W. The sum of squares of the fixed components of W is formed in */
1006 /* WFIXSQ, and the free components of W are set to BIGSTP. */
1008 bigstp = *adelt + *adelt;
1014 for (i__ = 1; i__ <= i__1; ++i__) {
1017 d__1 = xopt[i__] - sl[i__], d__2 = glag[i__];
1018 tempa = MIN2(d__1,d__2);
1020 d__1 = xopt[i__] - su[i__], d__2 = glag[i__];
1021 tempb = MAX2(d__1,d__2);
1022 if (tempa > zero || tempb < zero) {
1024 /* Computing 2nd power */
1026 ggfree += d__1 * d__1;
1030 if (ggfree == zero) {
1035 /* Investigate whether more components of W can be fixed. */
1038 temp = *adelt * *adelt - wfixsq;
1041 step = sqrt(temp / ggfree);
1044 for (i__ = 1; i__ <= i__1; ++i__) {
1045 if (w[i__] == bigstp) {
1046 temp = xopt[i__] - step * glag[i__];
1047 if (temp <= sl[i__]) {
1048 w[i__] = sl[i__] - xopt[i__];
1049 /* Computing 2nd power */
1051 wfixsq += d__1 * d__1;
1052 } else if (temp >= su[i__]) {
1053 w[i__] = su[i__] - xopt[i__];
1054 /* Computing 2nd power */
1056 wfixsq += d__1 * d__1;
1058 /* Computing 2nd power */
1060 ggfree += d__1 * d__1;
1065 if (wfixsq > wsqsav && ggfree > zero) {
1070 /* Set the remaining free components of W and all components of XALT, */
1071 /* except that W may be scaled later. */
1075 for (i__ = 1; i__ <= i__1; ++i__) {
1076 if (w[i__] == bigstp) {
1077 w[i__] = -step * glag[i__];
1080 d__3 = su[i__], d__4 = xopt[i__] + w[i__];
1081 d__1 = sl[i__], d__2 = MIN2(d__3,d__4);
1082 xalt[i__] = MAX2(d__1,d__2);
1083 } else if (w[i__] == zero) {
1084 xalt[i__] = xopt[i__];
1085 } else if (glag[i__] > zero) {
1086 xalt[i__] = sl[i__];
1088 xalt[i__] = su[i__];
1091 gw += glag[i__] * w[i__];
1094 /* Set CURV to the curvature of the KNEW-th Lagrange function along W. */
1095 /* Scale W by a factor less than one if that can reduce the modulus of */
1096 /* the Lagrange function at XOPT+W. Set CAUCHY to the final value of */
1097 /* the square of this function. */
1101 for (k = 1; k <= i__1; ++k) {
1104 for (j = 1; j <= i__2; ++j) {
1106 temp += xpt[k + j * xpt_dim1] * w[j];
1109 curv += hcol[k] * temp * temp;
1114 if (curv > -gw && curv < -const__ * gw) {
1117 for (i__ = 1; i__ <= i__1; ++i__) {
1118 temp = xopt[i__] + scale * w[i__];
1123 d__1 = sl[i__], d__2 = MIN2(d__3,temp);
1124 xalt[i__] = MAX2(d__1,d__2);
1126 /* Computing 2nd power */
1127 d__1 = half * gw * scale;
1128 *cauchy = d__1 * d__1;
1130 /* Computing 2nd power */
1131 d__1 = gw + half * curv;
1132 *cauchy = d__1 * d__1;
1135 /* If IFLAG is zero, then XALT is calculated as before after reversing */
1136 /* the sign of GLAG. Thus two XALT vectors become available. The one that */
1137 /* is chosen is the one that gives the larger value of CAUCHY. */
1141 for (i__ = 1; i__ <= i__1; ++i__) {
1142 glag[i__] = -glag[i__];
1144 w[*n + i__] = xalt[i__];
1150 if (csave > *cauchy) {
1152 for (i__ = 1; i__ <= i__1; ++i__) {
1154 xalt[i__] = w[*n + i__];
1162 static void trsbox_(int *n, int *npt, double *xpt,
1163 double *xopt, double *gopt, double *hq, double *pq,
1164 double *sl, double *su, double *delta, double *xnew,
1165 double *d__, double *gnew, double *xbdi, double *s,
1166 double *hs, double *hred, double *dsq, double *crvmin)
1168 /* System generated locals */
1169 int xpt_dim1, xpt_offset, i__1, i__2;
1170 double d__1, d__2, d__3, d__4;
1172 /* Local variables */
1176 double dhd, dhs, cth, one, shs, sth, ssq, half, beta, sdec, blen;
1180 double temp, zero, xsav, xsum, angbd, dredg, sredg;
1182 double resid, delsq, ggsav, tempa, tempb, ratio, sqstp, redmax,
1183 dredsq, redsav, onemin, gredsq, rednew;
1185 double rdprev, rdnext, stplen, stepsq;
1189 /* The arguments N, NPT, XPT, XOPT, GOPT, HQ, PQ, SL and SU have the same */
1190 /* meanings as the corresponding arguments of BOBYQB. */
1191 /* DELTA is the trust region radius for the present calculation, which */
1192 /* seeks a small value of the quadratic model within distance DELTA of */
1193 /* XOPT subject to the bounds on the variables. */
1194 /* XNEW will be set to a new vector of variables that is approximately */
1195 /* the one that minimizes the quadratic model within the trust region */
1196 /* subject to the SL and SU constraints on the variables. It satisfies */
1197 /* as equations the bounds that become active during the calculation. */
1198 /* D is the calculated trial step from XOPT, generated iteratively from an */
1199 /* initial value of zero. Thus XNEW is XOPT+D after the final iteration. */
1200 /* GNEW holds the gradient of the quadratic model at XOPT+D. It is updated */
1201 /* when D is updated. */
1202 /* XBDI is a working space vector. For I=1,2,...,N, the element XBDI(I) is */
1203 /* set to -1.0, 0.0, or 1.0, the value being nonzero if and only if the */
1204 /* I-th variable has become fixed at a bound, the bound being SL(I) or */
1205 /* SU(I) in the case XBDI(I)=-1.0 or XBDI(I)=1.0, respectively. This */
1206 /* information is accumulated during the construction of XNEW. */
1207 /* The arrays S, HS and HRED are also used for working space. They hold the */
1208 /* current search direction, and the changes in the gradient of Q along S */
1209 /* and the reduced D, respectively, where the reduced D is the same as D, */
1210 /* except that the components of the fixed variables are zero. */
1211 /* DSQ will be set to the square of the length of XNEW-XOPT. */
1212 /* CRVMIN is set to zero if D reaches the trust region boundary. Otherwise */
1213 /* it is set to the least curvature of H that occurs in the conjugate */
1214 /* gradient searches that are not restricted by any constraints. The */
1215 /* value CRVMIN=-1.0D0 is set, however, if all of these searches are */
1218 /* A version of the truncated conjugate gradient is applied. If a line */
1219 /* search is restricted by a constraint, then the procedure is restarted, */
1220 /* the values of the variables that are at their bounds being fixed. If */
1221 /* the trust region boundary is reached, then further changes may be made */
1222 /* to D, each one being in the two dimensional space that is spanned */
1223 /* by the current D and the gradient of Q at XOPT+D, staying on the trust */
1224 /* region boundary. Termination occurs when the reduction in Q seems to */
1225 /* be close to the greatest reduction that can be achieved. */
1227 /* Set some constants. */
1229 /* Parameter adjustments */
1231 xpt_offset = 1 + xpt_dim1;
1253 /* The sign of GOPT(I) gives the sign of the change to the I-th variable */
1254 /* that will reduce Q from its value at XOPT. Thus XBDI(I) shows whether */
1255 /* or not to fix the I-th variable at one of its bounds initially, with */
1256 /* NACT being set to the number of fixed variables. D and GNEW are also */
1257 /* set for the first iteration. DELSQ is the upper bound on the sum of */
1258 /* squares of the free variables. QRED is the reduction in Q so far. */
1264 for (i__ = 1; i__ <= i__1; ++i__) {
1266 if (xopt[i__] <= sl[i__]) {
1267 if (gopt[i__] >= zero) {
1270 } else if (xopt[i__] >= su[i__]) {
1271 if (gopt[i__] <= zero) {
1275 if (xbdi[i__] != zero) {
1280 gnew[i__] = gopt[i__];
1282 delsq = *delta * *delta;
1286 /* Set the next search direction of the conjugate gradient method. It is */
1287 /* the steepest descent direction initially and when the iterations are */
1288 /* restarted because a variable has just been fixed by a bound, and of */
1289 /* course the components of the fixed variables are zero. ITERMAX is an */
1290 /* upper bound on the indices of the conjugate gradient iterations. */
1297 for (i__ = 1; i__ <= i__1; ++i__) {
1298 if (xbdi[i__] != zero) {
1300 } else if (beta == zero) {
1301 s[i__] = -gnew[i__];
1303 s[i__] = beta * s[i__] - gnew[i__];
1306 /* Computing 2nd power */
1308 stepsq += d__1 * d__1;
1310 if (stepsq == zero) {
1315 itermax = iterc + *n - nact;
1317 if (gredsq * delsq <= qred * 1e-4 * qred) {
1321 /* Multiply the search direction by the second derivative matrix of Q and */
1322 /* calculate some scalars for the choice of steplength. Then set BLEN to */
1323 /* the length of the the step to the trust region boundary and STPLEN to */
1324 /* the steplength, ignoring the simple bounds. */
1332 for (i__ = 1; i__ <= i__1; ++i__) {
1333 if (xbdi[i__] == zero) {
1334 /* Computing 2nd power */
1336 resid -= d__1 * d__1;
1337 ds += s[i__] * d__[i__];
1338 shs += s[i__] * hs[i__];
1342 if (resid <= zero) {
1345 temp = sqrt(stepsq * resid + ds * ds);
1347 blen = (temp - ds) / stepsq;
1349 blen = resid / (temp + ds);
1354 d__1 = blen, d__2 = gredsq / shs;
1355 stplen = MIN2(d__1,d__2);
1358 /* Reduce STPLEN if necessary in order to preserve the simple bounds, */
1359 /* letting IACT be the index of the new constrained variable. */
1363 for (i__ = 1; i__ <= i__1; ++i__) {
1364 if (s[i__] != zero) {
1365 xsum = xopt[i__] + d__[i__];
1366 if (s[i__] > zero) {
1367 temp = (su[i__] - xsum) / s[i__];
1369 temp = (sl[i__] - xsum) / s[i__];
1371 if (temp < stplen) {
1379 /* Update CRVMIN, GNEW and D. Set SDEC to the decrease that occurs in Q. */
1382 if (stplen > zero) {
1384 temp = shs / stepsq;
1385 if (iact == 0 && temp > zero) {
1386 *crvmin = MIN2(*crvmin,temp);
1387 if (*crvmin == onemin) {
1394 for (i__ = 1; i__ <= i__1; ++i__) {
1395 gnew[i__] += stplen * hs[i__];
1396 if (xbdi[i__] == zero) {
1397 /* Computing 2nd power */
1399 gredsq += d__1 * d__1;
1402 d__[i__] += stplen * s[i__];
1405 d__1 = stplen * (ggsav - half * stplen * shs);
1406 sdec = MAX2(d__1,zero);
1410 /* Restart the conjugate gradient method if it has hit a new bound. */
1415 if (s[iact] < zero) {
1416 xbdi[iact] = onemin;
1418 /* Computing 2nd power */
1420 delsq -= d__1 * d__1;
1421 if (delsq <= zero) {
1427 /* If STPLEN is less than BLEN, then either apply another conjugate */
1428 /* gradient iteration or RETURN. */
1430 if (stplen < blen) {
1431 if (iterc == itermax) {
1434 if (sdec <= qred * .01) {
1437 beta = gredsq / ggsav;
1443 /* Prepare for the alternative iteration by calculating some scalars */
1444 /* and by multiplying the reduced D by the second derivative matrix of */
1445 /* Q, where S holds the reduced D in the call of GGMULT. */
1448 if (nact >= *n - 1) {
1455 for (i__ = 1; i__ <= i__1; ++i__) {
1456 if (xbdi[i__] == zero) {
1457 /* Computing 2nd power */
1459 dredsq += d__1 * d__1;
1460 dredg += d__[i__] * gnew[i__];
1461 /* Computing 2nd power */
1463 gredsq += d__1 * d__1;
1473 /* Let the search direction S be a linear combination of the reduced D */
1474 /* and the reduced G that is orthogonal to the reduced D. */
1478 temp = gredsq * dredsq - dredg * dredg;
1479 if (temp <= qred * 1e-4 * qred) {
1484 for (i__ = 1; i__ <= i__1; ++i__) {
1485 if (xbdi[i__] == zero) {
1486 s[i__] = (dredg * d__[i__] - dredsq * gnew[i__]) / temp;
1494 /* By considering the simple bounds on the variables, calculate an upper */
1495 /* bound on the tangent of half the angle of the alternative iteration, */
1496 /* namely ANGBD, except that, if already a free variable has reached a */
1497 /* bound, there is a branch back to label 100 after fixing that variable. */
1502 for (i__ = 1; i__ <= i__1; ++i__) {
1503 if (xbdi[i__] == zero) {
1504 tempa = xopt[i__] + d__[i__] - sl[i__];
1505 tempb = su[i__] - xopt[i__] - d__[i__];
1506 if (tempa <= zero) {
1510 } else if (tempb <= zero) {
1516 /* Computing 2nd power */
1518 /* Computing 2nd power */
1520 ssq = d__1 * d__1 + d__2 * d__2;
1521 /* Computing 2nd power */
1522 d__1 = xopt[i__] - sl[i__];
1523 temp = ssq - d__1 * d__1;
1525 temp = sqrt(temp) - s[i__];
1526 if (angbd * temp > tempa) {
1527 angbd = tempa / temp;
1532 /* Computing 2nd power */
1533 d__1 = su[i__] - xopt[i__];
1534 temp = ssq - d__1 * d__1;
1536 temp = sqrt(temp) + s[i__];
1537 if (angbd * temp > tempb) {
1538 angbd = tempb / temp;
1547 /* Calculate HHD and some curvatures for the alternative iteration. */
1555 for (i__ = 1; i__ <= i__1; ++i__) {
1556 if (xbdi[i__] == zero) {
1557 shs += s[i__] * hs[i__];
1558 dhs += d__[i__] * hs[i__];
1559 dhd += d__[i__] * hred[i__];
1564 /* Seek the greatest reduction in Q for a range of equally spaced values */
1565 /* of ANGT in [0,ANGBD], where ANGT is the tangent of half the angle of */
1566 /* the alternative iteration. */
1571 iu = (int) (angbd * 17. + 3.1);
1573 for (i__ = 1; i__ <= i__1; ++i__) {
1574 angt = angbd * (double) i__ / (double) iu;
1575 sth = (angt + angt) / (one + angt * angt);
1576 temp = shs + angt * (angt * dhd - dhs - dhs);
1577 rednew = sth * (angt * dredg - sredg - half * sth * temp);
1578 if (rednew > redmax) {
1582 } else if (i__ == isav + 1) {
1589 /* Return if the reduction is zero. Otherwise, set the sine and cosine */
1590 /* of the angle of the alternative iteration, and calculate SDEC. */
1596 temp = (rdnext - rdprev) / (redmax + redmax - rdprev - rdnext);
1597 angt = angbd * ((double) isav + half * temp) / (double) iu;
1599 cth = (one - angt * angt) / (one + angt * angt);
1600 sth = (angt + angt) / (one + angt * angt);
1601 temp = shs + angt * (angt * dhd - dhs - dhs);
1602 sdec = sth * (angt * dredg - sredg - half * sth * temp);
1607 /* Update GNEW, D and HRED. If the angle of the alternative iteration */
1608 /* is restricted by a bound on a free variable, that variable is fixed */
1614 for (i__ = 1; i__ <= i__1; ++i__) {
1615 gnew[i__] = gnew[i__] + (cth - one) * hred[i__] + sth * hs[i__];
1616 if (xbdi[i__] == zero) {
1617 d__[i__] = cth * d__[i__] + sth * s[i__];
1618 dredg += d__[i__] * gnew[i__];
1619 /* Computing 2nd power */
1621 gredsq += d__1 * d__1;
1624 hred[i__] = cth * hred[i__] + sth * hs[i__];
1627 if (iact > 0 && isav == iu) {
1633 /* If SDEC is sufficiently small, then RETURN after setting XNEW to */
1634 /* XOPT+D, giving careful attention to the bounds. */
1636 if (sdec > qred * .01) {
1642 for (i__ = 1; i__ <= i__1; ++i__) {
1645 d__3 = xopt[i__] + d__[i__], d__4 = su[i__];
1646 d__1 = MIN2(d__3,d__4), d__2 = sl[i__];
1647 xnew[i__] = MAX2(d__1,d__2);
1648 if (xbdi[i__] == onemin) {
1649 xnew[i__] = sl[i__];
1651 if (xbdi[i__] == one) {
1652 xnew[i__] = su[i__];
1654 d__[i__] = xnew[i__] - xopt[i__];
1656 /* Computing 2nd power */
1658 *dsq += d__1 * d__1;
1661 /* The following instructions multiply the current S-vector by the second */
1662 /* derivative matrix of the quadratic model, putting the product in HS. */
1663 /* They are reached from three different parts of the software above and */
1664 /* they can be regarded as an external subroutine. */
1669 for (j = 1; j <= i__1; ++j) {
1672 for (i__ = 1; i__ <= i__2; ++i__) {
1675 hs[j] += hq[ih] * s[i__];
1678 hs[i__] += hq[ih] * s[j];
1682 for (k = 1; k <= i__2; ++k) {
1683 if (pq[k] != zero) {
1686 for (j = 1; j <= i__1; ++j) {
1688 temp += xpt[k + j * xpt_dim1] * s[j];
1692 for (i__ = 1; i__ <= i__1; ++i__) {
1694 hs[i__] += temp * xpt[k + i__ * xpt_dim1];
1699 if (*crvmin != zero) {
1702 if (iterc > itcsav) {
1706 for (i__ = 1; i__ <= i__2; ++i__) {
1708 hred[i__] = hs[i__];
1713 static nlopt_result prelim_(int *n, int *npt, double *x,
1714 const double *xl, const double *xu, double *rhobeg,
1715 nlopt_stopping *stop,
1716 bobyqa_func calfun, void *calfun_data,
1717 double *xbase, double *xpt, double *fval,
1718 double *gopt, double *hq, double *pq, double *bmat,
1719 double *zmat, int *ndim, double *sl, double *su,
1722 /* System generated locals */
1723 int xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1,
1724 zmat_offset, i__1, i__2;
1725 double d__1, d__2, d__3, d__4;
1727 /* Local variables */
1729 int i__, j, k, ih, np, nfm;
1732 double two, fbeg, diff, half, temp, zero, recip, stepa, stepb;
1738 /* The arguments N, NPT, X, XL, XU, RHOBEG, and MAXFUN are the */
1739 /* same as the corresponding arguments in SUBROUTINE BOBYQA. */
1740 /* The arguments XBASE, XPT, FVAL, HQ, PQ, BMAT, ZMAT, NDIM, SL and SU */
1741 /* are the same as the corresponding arguments in BOBYQB, the elements */
1742 /* of SL and SU being set in BOBYQA. */
1743 /* GOPT is usually the gradient of the quadratic model at XOPT+XBASE, but */
1744 /* it is set by PRELIM to the gradient of the quadratic model at XBASE. */
1745 /* If XOPT is nonzero, BOBYQB will change it to its usual value later. */
1746 /* NF is maintaned as the number of calls of CALFUN so far. */
1747 /* KOPT will be such that the least calculated value of F so far is at */
1748 /* the point XPT(KOPT,.)+XBASE in the space of the variables. */
1750 /* SUBROUTINE PRELIM sets the elements of XBASE, XPT, FVAL, GOPT, HQ, PQ, */
1751 /* BMAT and ZMAT for the first iteration, and it maintains the values of */
1752 /* NF and KOPT. The vector X is also changed by PRELIM. */
1754 /* Set some constants. */
1756 /* Parameter adjustments */
1758 zmat_offset = 1 + zmat_dim1;
1759 zmat -= zmat_offset;
1761 xpt_offset = 1 + xpt_dim1;
1772 bmat_offset = 1 + bmat_dim1;
1773 bmat -= bmat_offset;
1782 rhosq = *rhobeg * *rhobeg;
1783 recip = one / rhosq;
1786 /* Set XBASE to the initial vector of variables, and set the initial */
1787 /* elements of XPT, BMAT, HQ, PQ and ZMAT to zero. */
1790 for (j = 1; j <= i__1; ++j) {
1793 for (k = 1; k <= i__2; ++k) {
1795 xpt[k + j * xpt_dim1] = zero;
1798 for (i__ = 1; i__ <= i__2; ++i__) {
1800 bmat[i__ + j * bmat_dim1] = zero;
1804 for (ih = 1; ih <= i__2; ++ih) {
1809 for (k = 1; k <= i__2; ++k) {
1812 for (j = 1; j <= i__1; ++j) {
1814 zmat[k + j * zmat_dim1] = zero;
1818 /* Begin the initialization procedure. NF becomes one more than the number */
1819 /* of function values so far. The coordinates of the displacement of the */
1820 /* next initial interpolation point from XBASE are set in XPT(NF+1,.). */
1827 if (nfm <= *n << 1) {
1828 if (nfm >= 1 && nfm <= *n) {
1830 if (su[nfm] == zero) {
1833 xpt[nf + nfm * xpt_dim1] = stepa;
1834 } else if (nfm > *n) {
1835 stepa = xpt[nf - *n + nfx * xpt_dim1];
1837 if (sl[nfx] == zero) {
1839 d__1 = two * *rhobeg, d__2 = su[nfx];
1840 stepb = MIN2(d__1,d__2);
1842 if (su[nfx] == zero) {
1844 d__1 = -two * *rhobeg, d__2 = sl[nfx];
1845 stepb = MAX2(d__1,d__2);
1847 xpt[nf + nfx * xpt_dim1] = stepb;
1850 itemp = (nfm - np) / *n;
1851 jpt = nfm - itemp * *n - *n;
1858 xpt[nf + ipt * xpt_dim1] = xpt[ipt + 1 + ipt * xpt_dim1];
1859 xpt[nf + jpt * xpt_dim1] = xpt[jpt + 1 + jpt * xpt_dim1];
1862 /* Calculate the next value of F. The least function value so far and */
1863 /* its index are required. */
1866 for (j = 1; j <= i__1; ++j) {
1869 d__3 = xl[j], d__4 = xbase[j] + xpt[nf + j * xpt_dim1];
1870 d__1 = MAX2(d__3,d__4), d__2 = xu[j];
1871 x[j] = MIN2(d__1,d__2);
1872 if (xpt[nf + j * xpt_dim1] == sl[j]) {
1875 if (xpt[nf + j * xpt_dim1] == su[j]) {
1881 f = calfun(*n, &x[1], calfun_data);
1886 } else if (f < fval[*kopt]) {
1890 /* Set the nonzero initial elements of BMAT and the quadratic model in the */
1891 /* cases when NF is at most 2*N+1. If NF exceeds N+1, then the positions */
1892 /* of the NF-th and (NF-N)-th interpolation points may be switched, in */
1893 /* order that the function value at the first of them contributes to the */
1894 /* off-diagonal second derivative terms of the initial quadratic model. */
1896 if (nf <= (*n << 1) + 1) {
1897 if (nf >= 2 && nf <= *n + 1) {
1898 gopt[nfm] = (f - fbeg) / stepa;
1899 if (*npt < nf + *n) {
1900 bmat[nfm * bmat_dim1 + 1] = -one / stepa;
1901 bmat[nf + nfm * bmat_dim1] = one / stepa;
1902 bmat[*npt + nfm + nfm * bmat_dim1] = -half * rhosq;
1904 } else if (nf >= *n + 2) {
1905 ih = nfx * (nfx + 1) / 2;
1906 temp = (f - fbeg) / stepb;
1907 diff = stepb - stepa;
1908 hq[ih] = two * (temp - gopt[nfx]) / diff;
1909 gopt[nfx] = (gopt[nfx] * stepb - temp * stepa) / diff;
1910 if (stepa * stepb < zero) {
1911 if (f < fval[nf - *n]) {
1912 fval[nf] = fval[nf - *n];
1917 xpt[nf - *n + nfx * xpt_dim1] = stepb;
1918 xpt[nf + nfx * xpt_dim1] = stepa;
1921 bmat[nfx * bmat_dim1 + 1] = -(stepa + stepb) / (stepa * stepb);
1922 bmat[nf + nfx * bmat_dim1] = -half / xpt[nf - *n + nfx *
1924 bmat[nf - *n + nfx * bmat_dim1] = -bmat[nfx * bmat_dim1 + 1] -
1925 bmat[nf + nfx * bmat_dim1];
1926 zmat[nfx * zmat_dim1 + 1] = sqrt(two) / (stepa * stepb);
1927 zmat[nf + nfx * zmat_dim1] = sqrt(half) / rhosq;
1928 zmat[nf - *n + nfx * zmat_dim1] = -zmat[nfx * zmat_dim1 + 1] -
1929 zmat[nf + nfx * zmat_dim1];
1932 /* Set the off-diagonal second derivatives of the Lagrange functions and */
1933 /* the initial quadratic model. */
1936 ih = ipt * (ipt - 1) / 2 + jpt;
1937 zmat[nfx * zmat_dim1 + 1] = recip;
1938 zmat[nf + nfx * zmat_dim1] = recip;
1939 zmat[ipt + 1 + nfx * zmat_dim1] = -recip;
1940 zmat[jpt + 1 + nfx * zmat_dim1] = -recip;
1941 temp = xpt[nf + ipt * xpt_dim1] * xpt[nf + jpt * xpt_dim1];
1942 hq[ih] = (fbeg - fval[ipt + 1] - fval[jpt + 1] + f) / temp;
1944 if (nlopt_stop_forced(stop)) return NLOPT_FORCED_STOP;
1945 else if (f < stop->minf_max) return NLOPT_MINF_MAX_REACHED;
1946 else if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
1947 else if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
1951 return NLOPT_SUCCESS;
1954 static nlopt_result bobyqb_(int *n, int *npt, double *x,
1955 const double *xl, const double *xu, double *rhobeg, double *
1957 nlopt_stopping *stop,
1958 bobyqa_func calfun, void *calfun_data,
1961 double *xpt, double *fval, double *xopt, double *gopt,
1962 double *hq, double *pq, double *bmat, double *zmat,
1963 int *ndim, double *sl, double *su, double *xnew,
1964 double *xalt, double *d__, double *vlag, double *w)
1966 /* System generated locals */
1967 int xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1,
1968 zmat_offset, i__1, i__2, i__3;
1969 double d__1, d__2, d__3, d__4;
1971 /* Local variables */
1973 int i__, j, k, ih, jj, nh, ip, jp;
1976 double den, one, ten, dsq, rho, sum, two, diff, half, beta, gisq;
1978 double temp, suma, sumb, bsum, fopt;
1982 double gqsq, dist, sumw, sumz, diffa, diffb, diffc, hdiag;
1984 double alpha, delta, adelt, denom, fsave, bdtol, delsq;
1986 double ratio, dnorm, vquad, pqold, tenth;
1988 double sumpq, scaden;
1989 double errbig, cauchy, fracsq, biglsq, densav;
1991 double crvmin, frhosq;
1996 nlopt_result rc = NLOPT_SUCCESS, rc2;
1998 /* The arguments N, NPT, X, XL, XU, RHOBEG, RHOEND, and MAXFUN */
1999 /* are identical to the corresponding arguments in SUBROUTINE BOBYQA. */
2000 /* XBASE holds a shift of origin that should reduce the contributions */
2001 /* from rounding errors to values of the model and Lagrange functions. */
2002 /* XPT is a two-dimensional array that holds the coordinates of the */
2003 /* interpolation points relative to XBASE. */
2004 /* FVAL holds the values of F at the interpolation points. */
2005 /* XOPT is set to the displacement from XBASE of the trust region centre. */
2006 /* GOPT holds the gradient of the quadratic model at XBASE+XOPT. */
2007 /* HQ holds the explicit second derivatives of the quadratic model. */
2008 /* PQ contains the parameters of the implicit second derivatives of the */
2009 /* quadratic model. */
2010 /* BMAT holds the last N columns of H. */
2011 /* ZMAT holds the factorization of the leading NPT by NPT submatrix of H, */
2012 /* this factorization being ZMAT times ZMAT^T, which provides both the */
2013 /* correct rank and positive semi-definiteness. */
2014 /* NDIM is the first dimension of BMAT and has the value NPT+N. */
2015 /* SL and SU hold the differences XL-XBASE and XU-XBASE, respectively. */
2016 /* All the components of every XOPT are going to satisfy the bounds */
2017 /* SL(I) .LEQ. XOPT(I) .LEQ. SU(I), with appropriate equalities when */
2018 /* XOPT is on a constraint boundary. */
2019 /* XNEW is chosen by SUBROUTINE TRSBOX or ALTMOV. Usually XBASE+XNEW is the */
2020 /* vector of variables for the next call of CALFUN. XNEW also satisfies */
2021 /* the SL and SU constraints in the way that has just been mentioned. */
2022 /* XALT is an alternative to XNEW, chosen by ALTMOV, that may replace XNEW */
2023 /* in order to increase the denominator in the updating of UPDATE. */
2024 /* D is reserved for a trial step from XOPT, which is usually XNEW-XOPT. */
2025 /* VLAG contains the values of the Lagrange functions at a new point X. */
2026 /* They are part of a product that requires VLAG to be of length NDIM. */
2027 /* W is a one-dimensional array that is used for working space. Its length */
2028 /* must be at least 3*NDIM = 3*(NPT+N). */
2030 /* Set some constants. */
2032 /* Parameter adjustments */
2034 zmat_offset = 1 + zmat_dim1;
2035 zmat -= zmat_offset;
2037 xpt_offset = 1 + xpt_dim1;
2049 bmat_offset = 1 + bmat_dim1;
2050 bmat -= bmat_offset;
2070 /* The call of PRELIM sets the elements of XBASE, XPT, FVAL, GOPT, HQ, PQ, */
2071 /* BMAT and ZMAT for the first iteration, with the corresponding values of */
2072 /* of NF and KOPT, which are the number of calls of CALFUN so far and the */
2073 /* index of the interpolation point at the trust region centre. Then the */
2074 /* initial XOPT is set too. The branch to label 720 occurs if MAXFUN is */
2075 /* less than NPT. GOPT will be updated if KOPT is different from KBASE. */
2077 rc2 = prelim_(n, npt, &x[1], &xl[1], &xu[1], rhobeg,
2078 stop, calfun, calfun_data,
2079 &xbase[1], &xpt[xpt_offset], &fval[1], &gopt[1], &hq[1], &pq[1], &bmat[
2080 bmat_offset], &zmat[zmat_offset], ndim, &sl[1], &su[1], &kopt);
2083 for (i__ = 1; i__ <= i__1; ++i__) {
2084 xopt[i__] = xpt[kopt + i__ * xpt_dim1];
2086 /* Computing 2nd power */
2088 xoptsq += d__1 * d__1;
2091 if (rc2 != NLOPT_SUCCESS) {
2097 /* Complete the settings that are required for the iterative procedure. */
2101 nresc = stop->nevals;
2106 nfsav = stop->nevals;
2108 /* Update GOPT if necessary before the first iteration and after each */
2109 /* call of RESCUE that makes a call of CALFUN. */
2112 if (kopt != kbase) {
2115 for (j = 1; j <= i__1; ++j) {
2117 for (i__ = 1; i__ <= i__2; ++i__) {
2120 gopt[j] += hq[ih] * xopt[i__];
2123 gopt[i__] += hq[ih] * xopt[j];
2126 if (stop->nevals > *npt) {
2128 for (k = 1; k <= i__2; ++k) {
2131 for (j = 1; j <= i__1; ++j) {
2133 temp += xpt[k + j * xpt_dim1] * xopt[j];
2135 temp = pq[k] * temp;
2137 for (i__ = 1; i__ <= i__1; ++i__) {
2139 gopt[i__] += temp * xpt[k + i__ * xpt_dim1];
2145 /* Generate the next point in the trust region that provides a small value */
2146 /* of the quadratic model subject to the constraints on the variables. */
2147 /* The int NTRITS is set to the number "trust region" iterations that */
2148 /* have occurred since the last "alternative" iteration. If the length */
2149 /* of XNEW-XOPT is less than HALF*RHO, however, then there is a branch to */
2150 /* label 650 or 680 with NTRITS=-1, instead of calculating F at XNEW. */
2153 trsbox_(n, npt, &xpt[xpt_offset], &xopt[1], &gopt[1], &hq[1], &pq[1], &sl[
2154 1], &su[1], &delta, &xnew[1], &d__[1], &w[1], &w[np], &w[np + *n],
2155 &w[np + (*n << 1)], &w[np + *n * 3], &dsq, &crvmin);
2157 d__1 = delta, d__2 = sqrt(dsq);
2158 dnorm = MIN2(d__1,d__2);
2159 if (dnorm < half * rho) {
2161 /* Computing 2nd power */
2163 distsq = d__1 * d__1;
2164 if (stop->nevals <= nfsav + 2) {
2168 /* The following choice between labels 650 and 680 depends on whether or */
2169 /* not our work with the current RHO seems to be complete. Either RHO is */
2170 /* decreased or termination occurs if the errors in the quadratic model at */
2171 /* the last three interpolation points compare favourably with predictions */
2172 /* of likely improvements to the model within distance HALF*RHO of XOPT. */
2175 d__1 = MAX2(diffa,diffb);
2176 errbig = MAX2(d__1,diffc);
2177 frhosq = rho * .125 * rho;
2178 if (crvmin > zero && errbig > frhosq * crvmin) {
2181 bdtol = errbig / rho;
2183 for (j = 1; j <= i__1; ++j) {
2185 if (xnew[j] == sl[j]) {
2188 if (xnew[j] == su[j]) {
2191 if (bdtest < bdtol) {
2192 curv = hq[(j + j * j) / 2];
2194 for (k = 1; k <= i__2; ++k) {
2196 /* Computing 2nd power */
2197 d__1 = xpt[k + j * xpt_dim1];
2198 curv += pq[k] * (d__1 * d__1);
2200 bdtest += half * curv * rho;
2201 if (bdtest < bdtol) {
2211 /* Severe cancellation is likely to occur if XOPT is too far from XBASE. */
2212 /* If the following test holds, then XBASE is shifted so that XOPT becomes */
2213 /* zero. The appropriate changes are made to BMAT and to the second */
2214 /* derivatives of the current model, beginning with the changes to BMAT */
2215 /* that do not depend on ZMAT. VLAG is used temporarily for working space. */
2218 if (dsq <= xoptsq * .001) {
2219 fracsq = xoptsq * .25;
2222 for (k = 1; k <= i__1; ++k) {
2224 sum = -half * xoptsq;
2226 for (i__ = 1; i__ <= i__2; ++i__) {
2228 sum += xpt[k + i__ * xpt_dim1] * xopt[i__];
2231 temp = fracsq - half * sum;
2233 for (i__ = 1; i__ <= i__2; ++i__) {
2234 w[i__] = bmat[k + i__ * bmat_dim1];
2235 vlag[i__] = sum * xpt[k + i__ * xpt_dim1] + temp * xopt[i__];
2238 for (j = 1; j <= i__3; ++j) {
2240 bmat[ip + j * bmat_dim1] = bmat[ip + j * bmat_dim1] + w[
2241 i__] * vlag[j] + vlag[i__] * w[j];
2246 /* Then the revisions of BMAT that depend on ZMAT are calculated. */
2249 for (jj = 1; jj <= i__3; ++jj) {
2253 for (k = 1; k <= i__2; ++k) {
2254 sumz += zmat[k + jj * zmat_dim1];
2255 vlag[k] = w[*npt + k] * zmat[k + jj * zmat_dim1];
2260 for (j = 1; j <= i__2; ++j) {
2261 sum = (fracsq * sumz - half * sumw) * xopt[j];
2263 for (k = 1; k <= i__1; ++k) {
2265 sum += vlag[k] * xpt[k + j * xpt_dim1];
2269 for (k = 1; k <= i__1; ++k) {
2271 bmat[k + j * bmat_dim1] += sum * zmat[k + jj * zmat_dim1];
2275 for (i__ = 1; i__ <= i__1; ++i__) {
2279 for (j = 1; j <= i__2; ++j) {
2281 bmat[ip + j * bmat_dim1] += temp * w[j];
2286 /* The following instructions complete the shift, including the changes */
2287 /* to the second derivative parameters of the quadratic model. */
2291 for (j = 1; j <= i__2; ++j) {
2292 w[j] = -half * sumpq * xopt[j];
2294 for (k = 1; k <= i__1; ++k) {
2295 w[j] += pq[k] * xpt[k + j * xpt_dim1];
2297 xpt[k + j * xpt_dim1] -= xopt[j];
2300 for (i__ = 1; i__ <= i__1; ++i__) {
2302 hq[ih] = hq[ih] + w[i__] * xopt[j] + xopt[i__] * w[j];
2304 bmat[*npt + i__ + j * bmat_dim1] = bmat[*npt + j + i__ *
2309 for (i__ = 1; i__ <= i__1; ++i__) {
2310 xbase[i__] += xopt[i__];
2311 xnew[i__] -= xopt[i__];
2312 sl[i__] -= xopt[i__];
2313 su[i__] -= xopt[i__];
2324 /* XBASE is also moved to XOPT by a call of RESCUE. This calculation is */
2325 /* more expensive than the previous shift, because new matrices BMAT and */
2326 /* ZMAT are generated from scratch, which may include the replacement of */
2327 /* interpolation points whose positions seem to be causing near linear */
2328 /* dependence in the interpolation conditions. Therefore RESCUE is called */
2329 /* only if rounding errors have reduced by at least a factor of two the */
2330 /* denominator of the formula for updating the H matrix. It provides a */
2331 /* useful safeguard, but is not invoked in most applications of BOBYQA. */
2334 nfsav = stop->nevals;
2336 rc2 = rescue_(n, npt, &xl[1], &xu[1],
2337 stop, calfun, calfun_data,
2338 &xbase[1], &xpt[xpt_offset], &fval[1], &xopt[1], &gopt[1],
2339 &hq[1], &pq[1], &bmat[bmat_offset], &zmat[zmat_offset], ndim,
2340 &sl[1], &su[1], &delta, &kopt, &vlag[1],
2341 &w[1], &w[*n + np], &w[*ndim + np]);
2343 /* XOPT is updated now in case the branch below to label 720 is taken. */
2344 /* Any updating of GOPT occurs after the branch below to label 20, which */
2345 /* leads to a trust region iteration as does the branch to label 60. */
2348 if (kopt != kbase) {
2350 for (i__ = 1; i__ <= i__1; ++i__) {
2351 xopt[i__] = xpt[kopt + i__ * xpt_dim1];
2353 /* Computing 2nd power */
2355 xoptsq += d__1 * d__1;
2358 if (rc2 != NLOPT_SUCCESS) {
2362 nresc = stop->nevals;
2363 if (nfsav < stop->nevals) {
2364 nfsav = stop->nevals;
2371 /* Pick two alternative vectors of variables, relative to XBASE, that */
2372 /* are suitable as new positions of the KNEW-th interpolation point. */
2373 /* Firstly, XNEW is set to the point on a line through XOPT and another */
2374 /* interpolation point that minimizes the predicted value of the next */
2375 /* denominator, subject to ||XNEW - XOPT|| .LEQ. ADELT and to the SL */
2376 /* and SU bounds. Secondly, XALT is set to the best feasible point on */
2377 /* a constrained version of the Cauchy step of the KNEW-th Lagrange */
2378 /* function, the corresponding value of the square of this function */
2379 /* being returned in CAUCHY. The choice between these alternatives is */
2380 /* going to be made when the denominator is calculated. */
2383 altmov_(n, npt, &xpt[xpt_offset], &xopt[1], &bmat[bmat_offset], &zmat[
2384 zmat_offset], ndim, &sl[1], &su[1], &kopt, &knew, &adelt, &xnew[1]
2385 , &xalt[1], &alpha, &cauchy, &w[1], &w[np], &w[*ndim + 1]);
2387 for (i__ = 1; i__ <= i__1; ++i__) {
2389 d__[i__] = xnew[i__] - xopt[i__];
2392 /* Calculate VLAG and BETA for the current choice of D. The scalar */
2393 /* product of D with XPT(K,.) is going to be held in W(NPT+K) for */
2394 /* use when VQUAD is calculated. */
2398 for (k = 1; k <= i__1; ++k) {
2403 for (j = 1; j <= i__2; ++j) {
2404 suma += xpt[k + j * xpt_dim1] * d__[j];
2405 sumb += xpt[k + j * xpt_dim1] * xopt[j];
2407 sum += bmat[k + j * bmat_dim1] * d__[j];
2409 w[k] = suma * (half * suma + sumb);
2416 for (jj = 1; jj <= i__1; ++jj) {
2419 for (k = 1; k <= i__2; ++k) {
2421 sum += zmat[k + jj * zmat_dim1] * w[k];
2425 for (k = 1; k <= i__2; ++k) {
2427 vlag[k] += sum * zmat[k + jj * zmat_dim1];
2434 for (j = 1; j <= i__2; ++j) {
2435 /* Computing 2nd power */
2440 for (k = 1; k <= i__1; ++k) {
2442 sum += w[k] * bmat[k + j * bmat_dim1];
2444 bsum += sum * d__[j];
2447 for (i__ = 1; i__ <= i__1; ++i__) {
2449 sum += bmat[jp + i__ * bmat_dim1] * d__[i__];
2452 bsum += sum * d__[j];
2454 dx += d__[j] * xopt[j];
2456 beta = dx * dx + dsq * (xoptsq + dx + dx + half * dsq) + beta - bsum;
2459 /* If NTRITS is zero, the denominator may be increased by replacing */
2460 /* the step D of ALTMOV by a Cauchy step. Then RESCUE may be called if */
2461 /* rounding errors have damaged the chosen denominator. */
2464 /* Computing 2nd power */
2466 denom = d__1 * d__1 + alpha * beta;
2467 if (denom < cauchy && cauchy > zero) {
2469 for (i__ = 1; i__ <= i__2; ++i__) {
2470 xnew[i__] = xalt[i__];
2472 d__[i__] = xnew[i__] - xopt[i__];
2477 /* Computing 2nd power */
2479 if (denom <= half * (d__1 * d__1)) {
2480 if (stop->nevals > nresc) {
2483 /* Return from BOBYQA because of much cancellation in a
2485 rc = NLOPT_ROUNDOFF_LIMITED;
2489 /* Alternatively, if NTRITS is positive, then set KNEW to the index of */
2490 /* the next interpolation point to be deleted to make room for a trust */
2491 /* region step. Again RESCUE may be called if rounding errors have damaged */
2492 /* the chosen denominator, which is the reason for attempting to select */
2493 /* KNEW before calculating the next value of the objective function. */
2496 delsq = delta * delta;
2501 for (k = 1; k <= i__2; ++k) {
2507 for (jj = 1; jj <= i__1; ++jj) {
2509 /* Computing 2nd power */
2510 d__1 = zmat[k + jj * zmat_dim1];
2511 hdiag += d__1 * d__1;
2513 /* Computing 2nd power */
2515 den = beta * hdiag + d__1 * d__1;
2518 for (j = 1; j <= i__1; ++j) {
2520 /* Computing 2nd power */
2521 d__1 = xpt[k + j * xpt_dim1] - xopt[j];
2522 distsq += d__1 * d__1;
2525 /* Computing 2nd power */
2526 d__3 = distsq / delsq;
2527 d__1 = one, d__2 = d__3 * d__3;
2528 temp = MAX2(d__1,d__2);
2529 if (temp * den > scaden) {
2530 scaden = temp * den;
2535 /* Computing 2nd power */
2537 d__1 = biglsq, d__2 = temp * (d__3 * d__3);
2538 biglsq = MAX2(d__1,d__2);
2542 if (scaden <= half * biglsq) {
2543 if (stop->nevals > nresc) {
2546 /* Return from BOBYQA because of much cancellation in a
2548 rc = NLOPT_ROUNDOFF_LIMITED;
2553 /* Put the variables for the next calculation of the objective function */
2554 /* in XNEW, with any adjustments for the bounds. */
2557 /* Calculate the value of the objective function at XBASE+XNEW, unless */
2558 /* the limit on the number of calculations of F has been reached. */
2562 for (i__ = 1; i__ <= i__2; ++i__) {
2565 d__3 = xl[i__], d__4 = xbase[i__] + xnew[i__];
2566 d__1 = MAX2(d__3,d__4), d__2 = xu[i__];
2567 x[i__] = MIN2(d__1,d__2);
2568 if (xnew[i__] == sl[i__]) {
2571 if (xnew[i__] == su[i__]) {
2577 if (nlopt_stop_forced(stop)) rc = NLOPT_FORCED_STOP;
2578 else if (nlopt_stop_evals(stop)) rc = NLOPT_MAXEVAL_REACHED;
2579 else if (nlopt_stop_time(stop)) rc = NLOPT_MAXTIME_REACHED;
2580 if (rc != NLOPT_SUCCESS) goto L720;
2583 f = calfun(*n, &x[1], calfun_data);
2586 rc = NLOPT_XTOL_REACHED;
2587 if (fsave < fval[kopt]) { *minf = f; return rc; }
2591 if (f < stop->minf_max) {
2593 return NLOPT_MINF_MAX_REACHED;
2596 /* Use the quadratic model to predict the change in F due to the step D, */
2597 /* and set DIFF to the error of this prediction. */
2603 for (j = 1; j <= i__2; ++j) {
2604 vquad += d__[j] * gopt[j];
2606 for (i__ = 1; i__ <= i__1; ++i__) {
2608 temp = d__[i__] * d__[j];
2613 vquad += hq[ih] * temp;
2617 for (k = 1; k <= i__1; ++k) {
2619 /* Computing 2nd power */
2621 vquad += half * pq[k] * (d__1 * d__1);
2623 diff = f - fopt - vquad;
2628 nfsav = stop->nevals;
2631 /* Pick the next value of DELTA after a trust region step. */
2634 if (vquad >= zero) {
2635 /* Return from BOBYQA because a trust region step has failed
2637 rc = NLOPT_ROUNDOFF_LIMITED; /* or FTOL_REACHED? */
2640 ratio = (f - fopt) / vquad;
2641 if (ratio <= tenth) {
2643 d__1 = half * delta;
2644 delta = MIN2(d__1,dnorm);
2645 } else if (ratio <= .7) {
2647 d__1 = half * delta;
2648 delta = MAX2(d__1,dnorm);
2651 d__1 = half * delta, d__2 = dnorm + dnorm;
2652 delta = MAX2(d__1,d__2);
2654 if (delta <= rho * 1.5) {
2658 /* Recalculate KNEW and DENOM if the new F is less than FOPT. */
2663 delsq = delta * delta;
2668 for (k = 1; k <= i__1; ++k) {
2671 for (jj = 1; jj <= i__2; ++jj) {
2673 /* Computing 2nd power */
2674 d__1 = zmat[k + jj * zmat_dim1];
2675 hdiag += d__1 * d__1;
2677 /* Computing 2nd power */
2679 den = beta * hdiag + d__1 * d__1;
2682 for (j = 1; j <= i__2; ++j) {
2684 /* Computing 2nd power */
2685 d__1 = xpt[k + j * xpt_dim1] - xnew[j];
2686 distsq += d__1 * d__1;
2689 /* Computing 2nd power */
2690 d__3 = distsq / delsq;
2691 d__1 = one, d__2 = d__3 * d__3;
2692 temp = MAX2(d__1,d__2);
2693 if (temp * den > scaden) {
2694 scaden = temp * den;
2700 /* Computing 2nd power */
2702 d__1 = biglsq, d__2 = temp * (d__3 * d__3);
2703 biglsq = MAX2(d__1,d__2);
2705 if (scaden <= half * biglsq) {
2712 /* Update BMAT and ZMAT, so that the KNEW-th interpolation point can be */
2713 /* moved. Also update the second derivative terms of the model. */
2715 update_(n, npt, &bmat[bmat_offset], &zmat[zmat_offset], ndim, &vlag[1], &
2716 beta, &denom, &knew, &w[1]);
2721 for (i__ = 1; i__ <= i__1; ++i__) {
2722 temp = pqold * xpt[knew + i__ * xpt_dim1];
2724 for (j = 1; j <= i__2; ++j) {
2727 hq[ih] += temp * xpt[knew + j * xpt_dim1];
2731 for (jj = 1; jj <= i__2; ++jj) {
2732 temp = diff * zmat[knew + jj * zmat_dim1];
2734 for (k = 1; k <= i__1; ++k) {
2736 pq[k] += temp * zmat[k + jj * zmat_dim1];
2740 /* Include the new interpolation point, and make the changes to GOPT at */
2741 /* the old XOPT that are caused by the updating of the quadratic model. */
2745 for (i__ = 1; i__ <= i__1; ++i__) {
2746 xpt[knew + i__ * xpt_dim1] = xnew[i__];
2748 w[i__] = bmat[knew + i__ * bmat_dim1];
2751 for (k = 1; k <= i__1; ++k) {
2754 for (jj = 1; jj <= i__2; ++jj) {
2756 suma += zmat[knew + jj * zmat_dim1] * zmat[k + jj * zmat_dim1];
2758 if (nlopt_isinf(suma)) {
2759 /* SGJ: detect singularity here (happend if we run
2760 for too many iterations) ... is there another way to recover? */
2761 rc = NLOPT_ROUNDOFF_LIMITED;
2766 for (j = 1; j <= i__2; ++j) {
2768 sumb += xpt[k + j * xpt_dim1] * xopt[j];
2772 for (i__ = 1; i__ <= i__2; ++i__) {
2774 w[i__] += temp * xpt[k + i__ * xpt_dim1];
2778 for (i__ = 1; i__ <= i__2; ++i__) {
2780 gopt[i__] += diff * w[i__];
2783 /* Update XOPT, GOPT and KOPT if the new calculated F is less than FOPT. */
2790 for (j = 1; j <= i__2; ++j) {
2792 /* Computing 2nd power */
2794 xoptsq += d__1 * d__1;
2796 for (i__ = 1; i__ <= i__1; ++i__) {
2799 gopt[j] += hq[ih] * d__[i__];
2802 gopt[i__] += hq[ih] * d__[j];
2806 for (k = 1; k <= i__1; ++k) {
2809 for (j = 1; j <= i__2; ++j) {
2811 temp += xpt[k + j * xpt_dim1] * d__[j];
2813 temp = pq[k] * temp;
2815 for (i__ = 1; i__ <= i__2; ++i__) {
2817 gopt[i__] += temp * xpt[k + i__ * xpt_dim1];
2820 if (nlopt_stop_ftol(stop, f, fopt)) {
2821 rc = NLOPT_FTOL_REACHED;
2826 /* Calculate the parameters of the least Frobenius norm interpolant to */
2827 /* the current data, the gradient of this interpolant at XOPT being put */
2828 /* into VLAG(NPT+I), I=1,2,...,N. */
2832 for (k = 1; k <= i__2; ++k) {
2833 vlag[k] = fval[k] - fval[kopt];
2838 for (j = 1; j <= i__2; ++j) {
2841 for (k = 1; k <= i__1; ++k) {
2843 sum += zmat[k + j * zmat_dim1] * vlag[k];
2846 for (k = 1; k <= i__1; ++k) {
2848 w[k] += sum * zmat[k + j * zmat_dim1];
2852 for (k = 1; k <= i__1; ++k) {
2855 for (j = 1; j <= i__2; ++j) {
2857 sum += xpt[k + j * xpt_dim1] * xopt[j];
2866 for (i__ = 1; i__ <= i__1; ++i__) {
2869 for (k = 1; k <= i__2; ++k) {
2871 sum = sum + bmat[k + i__ * bmat_dim1] * vlag[k] + xpt[k + i__
2874 if (xopt[i__] == sl[i__]) {
2876 d__2 = zero, d__3 = gopt[i__];
2877 /* Computing 2nd power */
2878 d__1 = MIN2(d__2,d__3);
2879 gqsq += d__1 * d__1;
2880 /* Computing 2nd power */
2881 d__1 = MIN2(zero,sum);
2882 gisq += d__1 * d__1;
2883 } else if (xopt[i__] == su[i__]) {
2885 d__2 = zero, d__3 = gopt[i__];
2886 /* Computing 2nd power */
2887 d__1 = MAX2(d__2,d__3);
2888 gqsq += d__1 * d__1;
2889 /* Computing 2nd power */
2890 d__1 = MAX2(zero,sum);
2891 gisq += d__1 * d__1;
2893 /* Computing 2nd power */
2895 gqsq += d__1 * d__1;
2899 vlag[*npt + i__] = sum;
2902 /* Test whether to replace the new quadratic model by the least Frobenius */
2903 /* norm interpolant, making the replacement if the test is satisfied. */
2906 if (gqsq < ten * gisq) {
2910 i__1 = MAX2(*npt,nh);
2911 for (i__ = 1; i__ <= i__1; ++i__) {
2913 gopt[i__] = vlag[*npt + i__];
2916 pq[i__] = w[*npt + i__];
2927 /* If a trust region step has provided a sufficient decrease in F, then */
2928 /* branch for another trust region calculation. The case NTRITS=0 occurs */
2929 /* when the new interpolation point was reached by an alternative step. */
2934 if (f <= fopt + tenth * vquad) {
2938 /* Alternatively, find out if the interpolation points are close enough */
2939 /* to the best point so far. */
2942 /* Computing 2nd power */
2944 /* Computing 2nd power */
2946 d__1 = d__3 * d__3, d__2 = d__4 * d__4;
2947 distsq = MAX2(d__1,d__2);
2951 for (k = 1; k <= i__1; ++k) {
2954 for (j = 1; j <= i__2; ++j) {
2956 /* Computing 2nd power */
2957 d__1 = xpt[k + j * xpt_dim1] - xopt[j];
2967 /* If KNEW is positive, then ALTMOV finds alternative new positions for */
2968 /* the KNEW-th interpolation point within distance ADELT of XOPT. It is */
2969 /* reached via label 90. Otherwise, there is a branch to label 60 for */
2970 /* another trust region iteration, unless the calculations with the */
2971 /* current RHO are complete. */
2974 dist = sqrt(distsq);
2977 d__1 = tenth * delta, d__2 = half * dist;
2978 delta = MIN2(d__1,d__2);
2979 if (delta <= rho * 1.5) {
2986 d__2 = tenth * dist;
2987 d__1 = MIN2(d__2,delta);
2988 adelt = MAX2(d__1,rho);
2989 dsq = adelt * adelt;
2998 if (MAX2(delta,dnorm) > rho) {
3002 /* The calculations with the current value of RHO are complete. Pick the */
3003 /* next values of RHO and DELTA. */
3006 if (rho > *rhoend) {
3008 ratio = rho / *rhoend;
3011 } else if (ratio <= 250.) {
3012 rho = sqrt(ratio) * *rhoend;
3016 delta = MAX2(delta,rho);
3018 nfsav = stop->nevals;
3022 /* Return from the calculation, after another Newton-Raphson step, if */
3023 /* it is too short to have been tried before. */
3029 /* originally: if (fval[kopt] <= fsave) -- changed by SGJ, since
3030 this seems like a slight optimization to not update x[]
3031 unnecessarily, at the expense of possibly not returning the
3032 best x[] found so far if the algorithm is stopped suddenly
3033 (e.g. runs out of time) ... it seems safer to execute this
3034 unconditionally, and the efficiency loss seems negligible. */
3037 for (i__ = 1; i__ <= i__1; ++i__) {
3040 d__3 = xl[i__], d__4 = xbase[i__] + xopt[i__];
3041 d__1 = MAX2(d__3,d__4), d__2 = xu[i__];
3042 x[i__] = MIN2(d__1,d__2);
3043 if (xopt[i__] == sl[i__]) {
3046 if (xopt[i__] == su[i__]) {
3057 /**************************************************************************/
3059 #define U(n) ((unsigned) (n))
3063 nlopt_func f; void *f_data;
3066 static double rescale_fun(int n, const double *x, void *d_)
3068 rescale_fun_data *d = (rescale_fun_data*) d_;
3069 nlopt_unscale(U(n), d->s, x, d->xs);
3070 return d->f(U(n), d->xs, NULL, d->f_data);
3073 nlopt_result bobyqa(int n, int npt, double *x,
3074 const double *xl, const double *xu,
3076 nlopt_stopping *stop, double *minf,
3077 nlopt_func f, void *f_data)
3079 /* System generated locals */
3083 /* Local variables */
3084 int j, id, np, iw, igo, ihq, ixb, ixa, ifv, isl, jsl, ipq, ivl, ixn,
3085 ixo, ixp, isu, jsu, ndim;
3089 double rhobeg, rhoend;
3090 double *w0 = NULL, *w;
3092 double *s = NULL, *sxl = NULL, *sxu = NULL, *xs = NULL;
3093 rescale_fun_data calfun_data;
3095 /* SGJ 2010: rescale parameters to make the initial step sizes dx
3096 equal in all directions */
3097 s = nlopt_compute_rescaling(U(n), dx);
3098 if (!s) return NLOPT_OUT_OF_MEMORY;
3099 for (j = 0; j < n; ++j)
3100 if (s[j] == 0 || !nlopt_isfinite(s[j])) {
3101 nlopt_stop_msg(stop, "invalid scaling %g of dimension %d: possible over/underflow?", s[j], j);
3102 ret = NLOPT_INVALID_ARGS; goto done;
3105 /* this statement must go before goto done, so that --x occurs */
3106 nlopt_rescale(U(n), s, x, x); --x;
3108 xs = (double *) malloc(sizeof(double) * (U(n)));
3109 if (!xs) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
3111 sxl = nlopt_new_rescaled(U(n), s, xl);
3112 if (!sxl) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
3114 sxu = nlopt_new_rescaled(U(n), s, xu);
3115 if (!sxu) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
3117 nlopt_reorder_bounds(n, sxl, sxu);
3119 rhobeg = fabs(dx[0] / s[0]); /* equals all other dx[i] after rescaling */
3122 calfun_data.xs = xs;
3124 calfun_data.f_data = f_data;
3126 /* SGJ, 2009: compute rhoend from NLopt stop info */
3127 rhoend = stop->xtol_rel * (rhobeg);
3128 for (j = 0; j < n; ++j)
3129 if (rhoend < stop->xtol_abs[j] / fabs(s[j]))
3130 rhoend = stop->xtol_abs[j] / fabs(s[j]);
3133 /* This subroutine seeks the least value of a function of many variables, */
3134 /* by applying a trust region method that forms quadratic models by */
3135 /* interpolation. There is usually some freedom in the interpolation */
3136 /* conditions, which is taken up by minimizing the Frobenius norm of */
3137 /* the change to the second derivative of the model, beginning with the */
3138 /* zero matrix. The values of the variables are constrained by upper and */
3139 /* lower bounds. The arguments of the subroutine are as follows. */
3141 /* N must be set to the number of variables and must be at least two. */
3142 /* NPT is the number of interpolation conditions. Its value must be in */
3143 /* the interval [N+2,(N+1)(N+2)/2]. Choices that exceed 2*N+1 are not */
3145 /* Initial values of the variables must be set in X(1),X(2),...,X(N). They */
3146 /* will be changed to the values that give the least calculated F. */
3147 /* For I=1,2,...,N, XL(I) and XU(I) must provide the lower and upper */
3148 /* bounds, respectively, on X(I). The construction of quadratic models */
3149 /* requires XL(I) to be strictly less than XU(I) for each I. Further, */
3150 /* the contribution to a model from changes to the I-th variable is */
3151 /* damaged severely by rounding errors if XU(I)-XL(I) is too small. */
3152 /* RHOBEG and RHOEND must be set to the initial and final values of a trust */
3153 /* region radius, so both must be positive with RHOEND no greater than */
3154 /* RHOBEG. Typically, RHOBEG should be about one tenth of the greatest */
3155 /* expected change to a variable, while RHOEND should indicate the */
3156 /* accuracy that is required in the final values of the variables. An */
3157 /* error return occurs if any of the differences XU(I)-XL(I), I=1,...,N, */
3158 /* is less than 2*RHOBEG. */
3159 /* MAXFUN must be set to an upper bound on the number of calls of CALFUN. */
3160 /* The array W will be used for working space. Its length must be at least */
3161 /* (NPT+5)*(NPT+N)+3*N*(N+5)/2. */
3163 /* SUBROUTINE CALFUN (N,X,F) has to be provided by the user. It must set */
3164 /* F to the value of the objective function for the current values of the */
3165 /* variables X(1),X(2),...,X(N), which are generated automatically in a */
3166 /* way that satisfies the bounds given in XL and XU. */
3168 /* Return if the value of NPT is unacceptable. */
3170 /* Parameter adjustments */
3176 if (npt < n + 2 || npt > (n + 2) * np / 2) {
3177 /* Return from BOBYQA because NPT is not in the required interval */
3178 ret = NLOPT_INVALID_ARGS;
3179 nlopt_stop_msg(stop, "invalid number of sampling points");
3183 /* Partition the working space array, so that different parts of it can */
3184 /* be treated separately during the calculation of BOBYQB. The partition */
3185 /* requires the first (NPT+2)*(NPT+N)+3*N*(N+5)/2 elements of W plus the */
3186 /* space that is taken by the last array in the argument list of BOBYQB. */
3191 ifv = ixp + n * npt;
3195 ipq = ihq + n * np / 2;
3197 izmat = ibmat + ndim * n;
3198 isl = izmat + npt * (npt - np);
3206 w0 = (double *) malloc(sizeof(double) * U((npt+5)*(npt+n)+3*n*(n+5)/2));
3207 if (!w0) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
3210 /* Return if there is insufficient space between the bounds. Modify the */
3211 /* initial X if necessary in order to avoid conflicts between the bounds */
3212 /* and the construction of the first quadratic model. The lower and upper */
3213 /* bounds on moves from the updated X are set now, in the ISL and ISU */
3214 /* partitions of W, in order to provide useful and exact information about */
3215 /* components of X that become within distance RHOBEG from their bounds. */
3219 for (j = 1; j <= i__1; ++j) {
3220 temp = xu[j] - xl[j];
3221 if (temp < rhobeg + rhobeg) {
3222 /* Return from BOBYQA because one of the differences
3223 XU(I)-XL(I)s is less than 2*RHOBEG. */
3224 ret = NLOPT_INVALID_ARGS;
3225 nlopt_stop_msg(stop, "insufficient space between the bounds: %g - %g < %g",
3226 xu[j], xl[j], rhobeg+rhobeg);
3231 w[jsl] = xl[j] - x[j];
3232 w[jsu] = xu[j] - x[j];
3233 if (w[jsl] >= -(rhobeg)) {
3234 if (w[jsl] >= zero) {
3239 x[j] = xl[j] + rhobeg;
3242 d__1 = xu[j] - x[j];
3243 w[jsu] = MAX2(d__1,rhobeg);
3245 } else if (w[jsu] <= rhobeg) {
3246 if (w[jsu] <= zero) {
3251 x[j] = xu[j] - rhobeg;
3253 d__1 = xl[j] - x[j], d__2 = -(rhobeg);
3254 w[jsl] = MIN2(d__1,d__2);
3261 /* Make the call of BOBYQB. */
3263 ret = bobyqb_(&n, &npt, &x[1], &xl[1], &xu[1], &rhobeg, &rhoend,
3264 stop, rescale_fun, &calfun_data, minf,
3265 &w[ixb], &w[ixp], &w[ifv], &w[ixo], &w[igo], &w[ihq], &w[ipq],
3266 &w[ibmat], &w[izmat], &ndim, &w[isl], &w[isu], &w[ixn], &w[ixa],
3267 &w[id], &w[ivl], &w[iw]);
3274 ++x; nlopt_unscale(U(n), s, x, x);