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 #define min(a,b) ((a) <= (b) ? (a) : (b))
13 #define max(a,b) ((a) >= (b) ? (a) : (b))
14 #define iabs(x) ((x) < 0 ? -(x) : (x))
16 static void update_(int *n, int *npt, double *bmat,
17 double *zmat, int *ndim, double *vlag, double *beta,
18 double *denom, int *knew, double *w)
20 /* System generated locals */
21 int bmat_dim1, bmat_offset, zmat_dim1, zmat_offset, i__1, i__2;
22 double d__1, d__2, d__3;
25 int i__, j, k, jl, jp;
26 double one, tau, temp;
28 double zero, alpha, tempa, tempb, ztest;
31 /* The arrays BMAT and ZMAT are updated, as required by the new position */
32 /* of the interpolation point that has the index KNEW. The vector VLAG has */
33 /* N+NPT components, set on entry to the first NPT and last N components */
34 /* of the product Hw in equation (4.11) of the Powell (2006) paper on */
35 /* NEWUOA. Further, BETA is set on entry to the value of the parameter */
36 /* with that name, and DENOM is set to the denominator of the updating */
37 /* formula. Elements of ZMAT may be treated as zero if their moduli are */
38 /* at most ZTEST. The first NDIM elements of W are used for working space. */
40 /* Set some constants. */
42 /* Parameter adjustments */
44 zmat_offset = 1 + zmat_dim1;
47 bmat_offset = 1 + bmat_dim1;
58 for (k = 1; k <= i__1; ++k) {
60 for (j = 1; j <= i__2; ++j) {
63 d__2 = ztest, d__3 = (d__1 = zmat[k + j * zmat_dim1], fabs(d__1));
64 ztest = max(d__2,d__3);
69 /* Apply the rotations that put zeros in the KNEW-th row of ZMAT. */
73 for (j = 2; j <= i__2; ++j) {
74 if ((d__1 = zmat[*knew + j * zmat_dim1], fabs(d__1)) > ztest) {
75 /* Computing 2nd power */
76 d__1 = zmat[*knew + zmat_dim1];
77 /* Computing 2nd power */
78 d__2 = zmat[*knew + j * zmat_dim1];
79 temp = sqrt(d__1 * d__1 + d__2 * d__2);
80 tempa = zmat[*knew + zmat_dim1] / temp;
81 tempb = zmat[*knew + j * zmat_dim1] / temp;
83 for (i__ = 1; i__ <= i__1; ++i__) {
84 temp = tempa * zmat[i__ + zmat_dim1] + tempb * zmat[i__ + j *
86 zmat[i__ + j * zmat_dim1] = tempa * zmat[i__ + j * zmat_dim1]
87 - tempb * zmat[i__ + zmat_dim1];
89 zmat[i__ + zmat_dim1] = temp;
92 zmat[*knew + j * zmat_dim1] = zero;
96 /* Put the first NPT components of the KNEW-th column of HLAG into W, */
97 /* and calculate the parameters of the updating formula. */
100 for (i__ = 1; i__ <= i__2; ++i__) {
101 w[i__] = zmat[*knew + zmat_dim1] * zmat[i__ + zmat_dim1];
108 /* Complete the updating of ZMAT. */
111 tempb = zmat[*knew + zmat_dim1] / temp;
114 for (i__ = 1; i__ <= i__2; ++i__) {
116 zmat[i__ + zmat_dim1] = tempa * zmat[i__ + zmat_dim1] - tempb * vlag[
120 /* Finally, update the matrix BMAT. */
123 for (j = 1; j <= i__2; ++j) {
125 w[jp] = bmat[*knew + j * bmat_dim1];
126 tempa = (alpha * vlag[jp] - tau * w[jp]) / *denom;
127 tempb = (-(*beta) * w[jp] - tau * vlag[jp]) / *denom;
129 for (i__ = 1; i__ <= i__1; ++i__) {
130 bmat[i__ + j * bmat_dim1] = bmat[i__ + j * bmat_dim1] + tempa *
131 vlag[i__] + tempb * w[i__];
133 bmat[jp + (i__ - *npt) * bmat_dim1] = bmat[i__ + j *
141 static nlopt_result rescue_(int *n, int *npt, const double *xl, const double *xu,
143 nlopt_stopping *stop,
144 bobyqa_func calfun, void *calfun_data,
145 double *xbase, double *xpt, double *fval, double *xopt, double *gopt,
146 double *hq, double *pq, double *bmat, double *zmat,
147 int *ndim, double *sl, double *su, /* int *nf, */
148 double *delta, int *kopt, double *vlag, double *
149 ptsaux, double *ptsid, double *w)
151 /* System generated locals */
152 int xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1,
153 zmat_offset, i__1, i__2, i__3;
154 double d__1, d__2, d__3, d__4;
156 /* Local variables */
158 int i__, j, k, ih, jp, ip, iq, np, iw;
163 double sum, diff, half, beta;
169 double zero, hdiag, fbase, sfrac, denom, vquad, sumpq;
170 double dsqmin, distsq, vlmxsq;
173 /* The arguments N, NPT, XL, XU, MAXFUN, XBASE, XPT, FVAL, XOPT, */
174 /* GOPT, HQ, PQ, BMAT, ZMAT, NDIM, SL and SU have the same meanings as */
175 /* the corresponding arguments of BOBYQB on the entry to RESCUE. */
176 /* NF is maintained as the number of calls of CALFUN so far, except that */
177 /* NF is set to -1 if the value of MAXFUN prevents further progress. */
178 /* KOPT is maintained so that FVAL(KOPT) is the least calculated function */
179 /* value. Its correct value must be given on entry. It is updated if a */
180 /* new least function value is found, but the corresponding changes to */
181 /* XOPT and GOPT have to be made later by the calling program. */
182 /* DELTA is the current trust region radius. */
183 /* VLAG is a working space vector that will be used for the values of the */
184 /* provisional Lagrange functions at each of the interpolation points. */
185 /* They are part of a product that requires VLAG to be of length NDIM. */
186 /* PTSAUX is also a working space array. For J=1,2,...,N, PTSAUX(1,J) and */
187 /* PTSAUX(2,J) specify the two positions of provisional interpolation */
188 /* points when a nonzero step is taken along e_J (the J-th coordinate */
189 /* direction) through XBASE+XOPT, as specified below. Usually these */
190 /* steps have length DELTA, but other lengths are chosen if necessary */
191 /* in order to satisfy the given bounds on the variables. */
192 /* PTSID is also a working space array. It has NPT components that denote */
193 /* provisional new positions of the original interpolation points, in */
194 /* case changes are needed to restore the linear independence of the */
195 /* interpolation conditions. The K-th point is a candidate for change */
196 /* if and only if PTSID(K) is nonzero. In this case let p and q be the */
197 /* int parts of PTSID(K) and (PTSID(K)-p) multiplied by N+1. If p */
198 /* and q are both positive, the step from XBASE+XOPT to the new K-th */
199 /* interpolation point is PTSAUX(1,p)*e_p + PTSAUX(1,q)*e_q. Otherwise */
200 /* the step is PTSAUX(1,p)*e_p or PTSAUX(2,q)*e_q in the cases q=0 or */
201 /* p=0, respectively. */
202 /* The first NDIM+NPT elements of the array W are used for working space. */
203 /* The final elements of BMAT and ZMAT are set in a well-conditioned way */
204 /* to the values that are appropriate for the new interpolation points. */
205 /* The elements of GOPT, HQ and PQ are also revised to the values that are */
206 /* appropriate to the final quadratic model. */
208 /* Set some constants. */
210 /* Parameter adjustments */
212 zmat_offset = 1 + zmat_dim1;
215 xpt_offset = 1 + xpt_dim1;
226 bmat_offset = 1 + bmat_dim1;
240 sfrac = half / (double) np;
243 /* Shift the interpolation points so that XOPT becomes the origin, and set */
244 /* the elements of ZMAT to zero. The value of SUMPQ is required in the */
245 /* updating of HQ below. The squares of the distances from XOPT to the */
246 /* other interpolation points are set at the end of W. Increments of WINC */
247 /* may be added later to these squares to balance the consideration of */
248 /* the choice of point that is going to become current. */
253 for (k = 1; k <= i__1; ++k) {
256 for (j = 1; j <= i__2; ++j) {
257 xpt[k + j * xpt_dim1] -= xopt[j];
259 /* Computing 2nd power */
260 d__1 = xpt[k + j * xpt_dim1];
261 distsq += d__1 * d__1;
264 w[*ndim + k] = distsq;
265 winc = max(winc,distsq);
267 for (j = 1; j <= i__2; ++j) {
269 zmat[k + j * zmat_dim1] = zero;
273 /* Update HQ so that HQ and PQ define the second derivatives of the model */
274 /* after XBASE has been shifted to the trust region centre. */
278 for (j = 1; j <= i__2; ++j) {
279 w[j] = half * sumpq * xopt[j];
281 for (k = 1; k <= i__1; ++k) {
283 w[j] += pq[k] * xpt[k + j * xpt_dim1];
286 for (i__ = 1; i__ <= i__1; ++i__) {
289 hq[ih] = hq[ih] + w[i__] * xopt[j] + w[j] * xopt[i__];
293 /* Shift XBASE, SL, SU and XOPT. Set the elements of BMAT to zero, and */
294 /* also set the elements of PTSAUX. */
297 for (j = 1; j <= i__1; ++j) {
303 d__1 = *delta, d__2 = su[j];
304 ptsaux[(j << 1) + 1] = min(d__1,d__2);
306 d__1 = -(*delta), d__2 = sl[j];
307 ptsaux[(j << 1) + 2] = max(d__1,d__2);
308 if (ptsaux[(j << 1) + 1] + ptsaux[(j << 1) + 2] < zero) {
309 temp = ptsaux[(j << 1) + 1];
310 ptsaux[(j << 1) + 1] = ptsaux[(j << 1) + 2];
311 ptsaux[(j << 1) + 2] = temp;
313 if ((d__2 = ptsaux[(j << 1) + 2], fabs(d__2)) < half * (d__1 = ptsaux[(
314 j << 1) + 1], fabs(d__1))) {
315 ptsaux[(j << 1) + 2] = half * ptsaux[(j << 1) + 1];
318 for (i__ = 1; i__ <= i__2; ++i__) {
320 bmat[i__ + j * bmat_dim1] = zero;
325 /* Set the identifiers of the artificial interpolation points that are */
326 /* along a coordinate direction from XOPT, and set the corresponding */
327 /* nonzero elements of BMAT and ZMAT. */
331 for (j = 1; j <= i__2; ++j) {
334 ptsid[jp] = (double) j + sfrac;
336 ptsid[jpn] = (double) j / (double) np + sfrac;
337 temp = one / (ptsaux[(j << 1) + 1] - ptsaux[(j << 1) + 2]);
338 bmat[jp + j * bmat_dim1] = -temp + one / ptsaux[(j << 1) + 1];
339 bmat[jpn + j * bmat_dim1] = temp + one / ptsaux[(j << 1) + 2];
340 bmat[j * bmat_dim1 + 1] = -bmat[jp + j * bmat_dim1] - bmat[jpn +
342 zmat[j * zmat_dim1 + 1] = sqrt(2.) / (d__1 = ptsaux[(j << 1) + 1]
343 * ptsaux[(j << 1) + 2], fabs(d__1));
344 zmat[jp + j * zmat_dim1] = zmat[j * zmat_dim1 + 1] * ptsaux[(j <<
346 zmat[jpn + j * zmat_dim1] = -zmat[j * zmat_dim1 + 1] * ptsaux[(j
349 bmat[j * bmat_dim1 + 1] = -one / ptsaux[(j << 1) + 1];
350 bmat[jp + j * bmat_dim1] = one / ptsaux[(j << 1) + 1];
351 /* Computing 2nd power */
352 d__1 = ptsaux[(j << 1) + 1];
353 bmat[j + *npt + j * bmat_dim1] = -half * (d__1 * d__1);
358 /* Set any remaining identifiers with their nonzero elements of ZMAT. */
360 if (*npt >= *n + np) {
362 for (k = np << 1; k <= i__2; ++k) {
363 iw = (int) (((double) (k - np) - half) / (double) (*n)
365 ip = k - np - iw * *n;
370 ptsid[k] = (double) ip + (double) iq / (double) np +
372 temp = one / (ptsaux[(ip << 1) + 1] * ptsaux[(iq << 1) + 1]);
373 zmat[(k - np) * zmat_dim1 + 1] = temp;
374 zmat[ip + 1 + (k - np) * zmat_dim1] = -temp;
375 zmat[iq + 1 + (k - np) * zmat_dim1] = -temp;
377 zmat[k + (k - np) * zmat_dim1] = temp;
384 /* Reorder the provisional points in the way that exchanges PTSID(KOLD) */
385 /* with PTSID(KNEW). */
389 for (j = 1; j <= i__2; ++j) {
390 temp = bmat[kold + j * bmat_dim1];
391 bmat[kold + j * bmat_dim1] = bmat[knew + j * bmat_dim1];
393 bmat[knew + j * bmat_dim1] = temp;
396 for (j = 1; j <= i__2; ++j) {
397 temp = zmat[kold + j * zmat_dim1];
398 zmat[kold + j * zmat_dim1] = zmat[knew + j * zmat_dim1];
400 zmat[knew + j * zmat_dim1] = temp;
402 ptsid[kold] = ptsid[knew];
404 w[*ndim + knew] = zero;
408 vlag[kold] = vlag[knew];
411 /* Update the BMAT and ZMAT matrices so that the status of the KNEW-th */
412 /* interpolation point can be changed from provisional to original. The */
413 /* branch to label 350 occurs if all the original points are reinstated. */
414 /* The nonnegative values of W(NDIM+K) are required in the search below. */
416 update_(n, npt, &bmat[bmat_offset], &zmat[zmat_offset], ndim, &vlag[1]
417 , &beta, &denom, &knew, &w[1]);
422 for (k = 1; k <= i__2; ++k) {
424 w[*ndim + k] = (d__1 = w[*ndim + k], fabs(d__1));
428 /* Pick the index KNEW of an original interpolation point that has not */
429 /* yet replaced one of the provisional interpolation points, giving */
430 /* attention to the closeness to XOPT and to previous tries with KNEW. */
435 for (k = 1; k <= i__2; ++k) {
436 if (w[*ndim + k] > zero) {
437 if (dsqmin == zero || w[*ndim + k] < dsqmin) {
439 dsqmin = w[*ndim + k];
444 if (dsqmin == zero) {
448 /* Form the W-vector of the chosen original interpolation point. */
451 for (j = 1; j <= i__2; ++j) {
453 w[*npt + j] = xpt[knew + j * xpt_dim1];
456 for (k = 1; k <= i__2; ++k) {
459 } else if (ptsid[k] == zero) {
461 for (j = 1; j <= i__1; ++j) {
463 sum += w[*npt + j] * xpt[k + j * xpt_dim1];
468 sum = w[*npt + ip] * ptsaux[(ip << 1) + 1];
470 iq = (int) ((double) np * ptsid[k] - (double) (ip *
477 sum += w[*npt + iq] * ptsaux[iw + (iq << 1)];
481 w[k] = half * sum * sum;
484 /* Calculate VLAG and BETA for the required updating of the H matrix if */
485 /* XPT(KNEW,.) is reinstated in the set of interpolation points. */
488 for (k = 1; k <= i__2; ++k) {
491 for (j = 1; j <= i__1; ++j) {
493 sum += bmat[k + j * bmat_dim1] * w[*npt + j];
500 for (j = 1; j <= i__2; ++j) {
503 for (k = 1; k <= i__1; ++k) {
505 sum += zmat[k + j * zmat_dim1] * w[k];
509 for (k = 1; k <= i__1; ++k) {
511 vlag[k] += sum * zmat[k + j * zmat_dim1];
517 for (j = 1; j <= i__1; ++j) {
520 for (k = 1; k <= i__2; ++k) {
522 sum += bmat[k + j * bmat_dim1] * w[k];
527 for (ip = *npt + 1; ip <= i__2; ++ip) {
529 sum += bmat[ip + j * bmat_dim1] * w[ip];
534 /* Computing 2nd power */
535 d__1 = xpt[knew + j * xpt_dim1];
536 distsq += d__1 * d__1;
538 beta = half * distsq * distsq + beta - bsum;
541 /* KOLD is set to the index of the provisional interpolation point that is */
542 /* going to be deleted to make way for the KNEW-th original interpolation */
543 /* point. The choice of KOLD is governed by the avoidance of a small value */
544 /* of the denominator in the updating calculation of UPDATE. */
549 for (k = 1; k <= i__1; ++k) {
550 if (ptsid[k] != zero) {
553 for (j = 1; j <= i__2; ++j) {
555 /* Computing 2nd power */
556 d__1 = zmat[k + j * zmat_dim1];
557 hdiag += d__1 * d__1;
559 /* Computing 2nd power */
561 den = beta * hdiag + d__1 * d__1;
569 /* Computing 2nd power */
571 d__1 = vlmxsq, d__2 = d__3 * d__3;
572 vlmxsq = max(d__1,d__2);
574 if (denom <= vlmxsq * .01) {
575 w[*ndim + knew] = -w[*ndim + knew] - winc;
580 /* When label 260 is reached, all the final positions of the interpolation */
581 /* points have been chosen although any changes have not been included yet */
582 /* in XPT. Also the final BMAT and ZMAT matrices are complete, but, apart */
583 /* from the shift of XBASE, the updating of the quadratic model remains to */
584 /* be done. The following cycle through the new interpolation points begins */
585 /* by putting the new point in XPT(KPT,.) and by setting PQ(KPT) to zero, */
586 /* except that a RETURN occurs if MAXFUN prohibits another value of F. */
590 for (kpt = 1; kpt <= i__1; ++kpt) {
591 if (ptsid[kpt] == zero) {
595 if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
596 else if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
600 for (j = 1; j <= i__2; ++j) {
601 w[j] = xpt[kpt + j * xpt_dim1];
602 xpt[kpt + j * xpt_dim1] = zero;
603 temp = pq[kpt] * w[j];
605 for (i__ = 1; i__ <= i__3; ++i__) {
608 hq[ih] += temp * w[i__];
612 ip = (int) ptsid[kpt];
613 iq = (int) ((double) np * ptsid[kpt] - (double) (ip * np))
616 xp = ptsaux[(ip << 1) + 1];
617 xpt[kpt + ip * xpt_dim1] = xp;
620 xq = ptsaux[(iq << 1) + 1];
622 xq = ptsaux[(iq << 1) + 2];
624 xpt[kpt + iq * xpt_dim1] = xq;
627 /* Set VQUAD to the value of the current model at the new point. */
631 ihp = (ip + ip * ip) / 2;
632 vquad += xp * (gopt[ip] + half * xp * hq[ihp]);
635 ihq = (iq + iq * iq) / 2;
636 vquad += xq * (gopt[iq] + half * xq * hq[ihq]);
638 iw = max(ihp,ihq) - (i__3 = ip - iq, iabs(i__3));
639 vquad += xp * xq * hq[iw];
643 for (k = 1; k <= i__3; ++k) {
646 temp += xp * xpt[k + ip * xpt_dim1];
649 temp += xq * xpt[k + iq * xpt_dim1];
652 vquad += half * pq[k] * temp * temp;
655 /* Calculate F at the new interpolation point, and set DIFF to the factor */
656 /* that is going to multiply the KPT-th Lagrange function when the model */
657 /* is updated to provide interpolation to the new function value. */
660 for (i__ = 1; i__ <= i__3; ++i__) {
663 d__3 = xl[i__], d__4 = xbase[i__] + xpt[kpt + i__ * xpt_dim1];
664 d__1 = max(d__3,d__4), d__2 = xu[i__];
665 w[i__] = min(d__1,d__2);
666 if (xpt[kpt + i__ * xpt_dim1] == sl[i__]) {
669 if (xpt[kpt + i__ * xpt_dim1] == su[i__]) {
676 f = calfun(*n, &w[1], calfun_data);
678 if (f < fval[*kopt]) {
681 if (f < stop->minf_max) return NLOPT_MINF_MAX_REACHED;
682 else if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
683 else if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
687 /* Update the quadratic model. The RETURN from the subroutine occurs when */
688 /* all the new interpolation points are included in the model. */
691 for (i__ = 1; i__ <= i__3; ++i__) {
693 gopt[i__] += diff * bmat[kpt + i__ * bmat_dim1];
696 for (k = 1; k <= i__3; ++k) {
699 for (j = 1; j <= i__2; ++j) {
701 sum += zmat[k + j * zmat_dim1] * zmat[kpt + j * zmat_dim1];
704 if (ptsid[k] == zero) {
708 iq = (int) ((double) np * ptsid[k] - (double) (ip
710 ihq = (iq * iq + iq) / 2;
712 /* Computing 2nd power */
713 d__1 = ptsaux[(iq << 1) + 2];
714 hq[ihq] += temp * (d__1 * d__1);
716 ihp = (ip * ip + ip) / 2;
717 /* Computing 2nd power */
718 d__1 = ptsaux[(ip << 1) + 1];
719 hq[ihp] += temp * (d__1 * d__1);
721 /* Computing 2nd power */
722 d__1 = ptsaux[(iq << 1) + 1];
723 hq[ihq] += temp * (d__1 * d__1);
724 iw = max(ihp,ihq) - (i__2 = iq - ip, iabs(i__2));
725 hq[iw] += temp * ptsaux[(ip << 1) + 1] * ptsaux[(iq <<
737 return NLOPT_SUCCESS;
740 static void altmov_(int *n, int *npt, double *xpt,
741 double *xopt, double *bmat, double *zmat, int *ndim,
742 double *sl, double *su, int *kopt, int *knew,
743 double *adelt, double *xnew, double *xalt, double *
744 alpha, double *cauchy, double *glag, double *hcol,
747 /* System generated locals */
748 int xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1,
749 zmat_offset, i__1, i__2;
750 double d__1, d__2, d__3, d__4;
752 /* Local variables */
754 double ha, gw, one, diff, half;
758 double vlag, subd, temp;
760 double step, zero, curv;
762 double scale, csave, tempa, tempb, tempd, const__, sumin, ggfree;
764 double dderiv, bigstp, predsq, presav, distsq, stpsav, wfixsq, wsqsav;
767 /* The arguments N, NPT, XPT, XOPT, BMAT, ZMAT, NDIM, SL and SU all have */
768 /* the same meanings as the corresponding arguments of BOBYQB. */
769 /* KOPT is the index of the optimal interpolation point. */
770 /* KNEW is the index of the interpolation point that is going to be moved. */
771 /* ADELT is the current trust region bound. */
772 /* XNEW will be set to a suitable new position for the interpolation point */
773 /* XPT(KNEW,.). Specifically, it satisfies the SL, SU and trust region */
774 /* bounds and it should provide a large denominator in the next call of */
775 /* UPDATE. The step XNEW-XOPT from XOPT is restricted to moves along the */
776 /* straight lines through XOPT and another interpolation point. */
777 /* XALT also provides a large value of the modulus of the KNEW-th Lagrange */
778 /* function subject to the constraints that have been mentioned, its main */
779 /* difference from XNEW being that XALT-XOPT is a constrained version of */
780 /* the Cauchy step within the trust region. An exception is that XALT is */
781 /* not calculated if all components of GLAG (see below) are zero. */
782 /* ALPHA will be set to the KNEW-th diagonal element of the H matrix. */
783 /* CAUCHY will be set to the square of the KNEW-th Lagrange function at */
784 /* the step XALT-XOPT from XOPT for the vector XALT that is returned, */
785 /* except that CAUCHY is set to zero if XALT is not calculated. */
786 /* GLAG is a working space vector of length N for the gradient of the */
787 /* KNEW-th Lagrange function at XOPT. */
788 /* HCOL is a working space vector of length NPT for the second derivative */
789 /* coefficients of the KNEW-th Lagrange function. */
790 /* W is a working space vector of length 2N that is going to hold the */
791 /* constrained Cauchy step from XOPT of the Lagrange function, followed */
792 /* by the downhill version of XALT when the uphill step is calculated. */
794 /* Set the first NPT components of W to the leading elements of the */
795 /* KNEW-th column of the H matrix. */
797 /* Parameter adjustments */
799 zmat_offset = 1 + zmat_dim1;
802 xpt_offset = 1 + xpt_dim1;
806 bmat_offset = 1 + bmat_dim1;
820 const__ = one + sqrt(2.);
822 for (k = 1; k <= i__1; ++k) {
826 i__1 = *npt - *n - 1;
827 for (j = 1; j <= i__1; ++j) {
828 temp = zmat[*knew + j * zmat_dim1];
830 for (k = 1; k <= i__2; ++k) {
832 hcol[k] += temp * zmat[k + j * zmat_dim1];
835 *alpha = hcol[*knew];
838 /* Calculate the gradient of the KNEW-th Lagrange function at XOPT. */
841 for (i__ = 1; i__ <= i__2; ++i__) {
843 glag[i__] = bmat[*knew + i__ * bmat_dim1];
846 for (k = 1; k <= i__2; ++k) {
849 for (j = 1; j <= i__1; ++j) {
851 temp += xpt[k + j * xpt_dim1] * xopt[j];
853 temp = hcol[k] * temp;
855 for (i__ = 1; i__ <= i__1; ++i__) {
857 glag[i__] += temp * xpt[k + i__ * xpt_dim1];
861 /* Search for a large denominator along the straight lines through XOPT */
862 /* and another interpolation point. SLBD and SUBD will be lower and upper */
863 /* bounds on the step along each of these lines in turn. PREDSQ will be */
864 /* set to the square of the predicted denominator for each line. PRESAV */
865 /* will be set to the largest admissible value of PREDSQ that occurs. */
869 for (k = 1; k <= i__1; ++k) {
876 for (i__ = 1; i__ <= i__2; ++i__) {
877 temp = xpt[k + i__ * xpt_dim1] - xopt[i__];
878 dderiv += glag[i__] * temp;
880 distsq += temp * temp;
882 subd = *adelt / sqrt(distsq);
886 sumin = min(one,subd);
888 /* Revise SLBD and SUBD if necessary because of the bounds in SL and SU. */
891 for (i__ = 1; i__ <= i__2; ++i__) {
892 temp = xpt[k + i__ * xpt_dim1] - xopt[i__];
894 if (slbd * temp < sl[i__] - xopt[i__]) {
895 slbd = (sl[i__] - xopt[i__]) / temp;
898 if (subd * temp > su[i__] - xopt[i__]) {
900 d__1 = sumin, d__2 = (su[i__] - xopt[i__]) / temp;
901 subd = max(d__1,d__2);
904 } else if (temp < zero) {
905 if (slbd * temp > su[i__] - xopt[i__]) {
906 slbd = (su[i__] - xopt[i__]) / temp;
909 if (subd * temp < sl[i__] - xopt[i__]) {
911 d__1 = sumin, d__2 = (sl[i__] - xopt[i__]) / temp;
912 subd = max(d__1,d__2);
919 /* Seek a large modulus of the KNEW-th Lagrange function when the index */
920 /* of the other interpolation point on the line through XOPT is KNEW. */
925 vlag = slbd * (dderiv - slbd * diff);
927 temp = subd * (dderiv - subd * diff);
928 if (fabs(temp) > fabs(vlag)) {
933 tempd = half * dderiv;
934 tempa = tempd - diff * slbd;
935 tempb = tempd - diff * subd;
936 if (tempa * tempb < zero) {
937 temp = tempd * tempd / diff;
938 if (fabs(temp) > fabs(vlag)) {
945 /* Search along each of the other lines through XOPT and another point. */
949 vlag = slbd * (one - slbd);
951 temp = subd * (one - subd);
952 if (fabs(temp) > fabs(vlag)) {
958 if (fabs(vlag) < .25) {
967 /* Calculate PREDSQ for the current line search and maintain PRESAV. */
969 temp = step * (one - step) * distsq;
970 predsq = vlag * vlag * (vlag * vlag + ha * temp * temp);
971 if (predsq > presav) {
981 /* Construct XNEW in a way that satisfies the bound constraints exactly. */
984 for (i__ = 1; i__ <= i__1; ++i__) {
985 temp = xopt[i__] + stpsav * (xpt[ksav + i__ * xpt_dim1] - xopt[i__]);
990 d__1 = sl[i__], d__2 = min(d__3,temp);
991 xnew[i__] = max(d__1,d__2);
994 xnew[-ibdsav] = sl[-ibdsav];
997 xnew[ibdsav] = su[ibdsav];
1000 /* Prepare for the iterative method that assembles the constrained Cauchy */
1001 /* step in W. The sum of squares of the fixed components of W is formed in */
1002 /* WFIXSQ, and the free components of W are set to BIGSTP. */
1004 bigstp = *adelt + *adelt;
1010 for (i__ = 1; i__ <= i__1; ++i__) {
1013 d__1 = xopt[i__] - sl[i__], d__2 = glag[i__];
1014 tempa = min(d__1,d__2);
1016 d__1 = xopt[i__] - su[i__], d__2 = glag[i__];
1017 tempb = max(d__1,d__2);
1018 if (tempa > zero || tempb < zero) {
1020 /* Computing 2nd power */
1022 ggfree += d__1 * d__1;
1026 if (ggfree == zero) {
1031 /* Investigate whether more components of W can be fixed. */
1034 temp = *adelt * *adelt - wfixsq;
1037 step = sqrt(temp / ggfree);
1040 for (i__ = 1; i__ <= i__1; ++i__) {
1041 if (w[i__] == bigstp) {
1042 temp = xopt[i__] - step * glag[i__];
1043 if (temp <= sl[i__]) {
1044 w[i__] = sl[i__] - xopt[i__];
1045 /* Computing 2nd power */
1047 wfixsq += d__1 * d__1;
1048 } else if (temp >= su[i__]) {
1049 w[i__] = su[i__] - xopt[i__];
1050 /* Computing 2nd power */
1052 wfixsq += d__1 * d__1;
1054 /* Computing 2nd power */
1056 ggfree += d__1 * d__1;
1061 if (wfixsq > wsqsav && ggfree > zero) {
1066 /* Set the remaining free components of W and all components of XALT, */
1067 /* except that W may be scaled later. */
1071 for (i__ = 1; i__ <= i__1; ++i__) {
1072 if (w[i__] == bigstp) {
1073 w[i__] = -step * glag[i__];
1076 d__3 = su[i__], d__4 = xopt[i__] + w[i__];
1077 d__1 = sl[i__], d__2 = min(d__3,d__4);
1078 xalt[i__] = max(d__1,d__2);
1079 } else if (w[i__] == zero) {
1080 xalt[i__] = xopt[i__];
1081 } else if (glag[i__] > zero) {
1082 xalt[i__] = sl[i__];
1084 xalt[i__] = su[i__];
1087 gw += glag[i__] * w[i__];
1090 /* Set CURV to the curvature of the KNEW-th Lagrange function along W. */
1091 /* Scale W by a factor less than one if that can reduce the modulus of */
1092 /* the Lagrange function at XOPT+W. Set CAUCHY to the final value of */
1093 /* the square of this function. */
1097 for (k = 1; k <= i__1; ++k) {
1100 for (j = 1; j <= i__2; ++j) {
1102 temp += xpt[k + j * xpt_dim1] * w[j];
1105 curv += hcol[k] * temp * temp;
1110 if (curv > -gw && curv < -const__ * gw) {
1113 for (i__ = 1; i__ <= i__1; ++i__) {
1114 temp = xopt[i__] + scale * w[i__];
1119 d__1 = sl[i__], d__2 = min(d__3,temp);
1120 xalt[i__] = max(d__1,d__2);
1122 /* Computing 2nd power */
1123 d__1 = half * gw * scale;
1124 *cauchy = d__1 * d__1;
1126 /* Computing 2nd power */
1127 d__1 = gw + half * curv;
1128 *cauchy = d__1 * d__1;
1131 /* If IFLAG is zero, then XALT is calculated as before after reversing */
1132 /* the sign of GLAG. Thus two XALT vectors become available. The one that */
1133 /* is chosen is the one that gives the larger value of CAUCHY. */
1137 for (i__ = 1; i__ <= i__1; ++i__) {
1138 glag[i__] = -glag[i__];
1140 w[*n + i__] = xalt[i__];
1146 if (csave > *cauchy) {
1148 for (i__ = 1; i__ <= i__1; ++i__) {
1150 xalt[i__] = w[*n + i__];
1158 static void trsbox_(int *n, int *npt, double *xpt,
1159 double *xopt, double *gopt, double *hq, double *pq,
1160 double *sl, double *su, double *delta, double *xnew,
1161 double *d__, double *gnew, double *xbdi, double *s,
1162 double *hs, double *hred, double *dsq, double *crvmin)
1164 /* System generated locals */
1165 int xpt_dim1, xpt_offset, i__1, i__2;
1166 double d__1, d__2, d__3, d__4;
1168 /* Local variables */
1172 double dhd, dhs, cth, one, shs, sth, ssq, half, beta, sdec, blen;
1176 double temp, zero, xsav, xsum, angbd, dredg, sredg;
1178 double resid, delsq, ggsav, tempa, tempb, ratio, sqstp, redmax,
1179 dredsq, redsav, onemin, gredsq, rednew;
1181 double rdprev, rdnext, stplen, stepsq;
1185 /* The arguments N, NPT, XPT, XOPT, GOPT, HQ, PQ, SL and SU have the same */
1186 /* meanings as the corresponding arguments of BOBYQB. */
1187 /* DELTA is the trust region radius for the present calculation, which */
1188 /* seeks a small value of the quadratic model within distance DELTA of */
1189 /* XOPT subject to the bounds on the variables. */
1190 /* XNEW will be set to a new vector of variables that is approximately */
1191 /* the one that minimizes the quadratic model within the trust region */
1192 /* subject to the SL and SU constraints on the variables. It satisfies */
1193 /* as equations the bounds that become active during the calculation. */
1194 /* D is the calculated trial step from XOPT, generated iteratively from an */
1195 /* initial value of zero. Thus XNEW is XOPT+D after the final iteration. */
1196 /* GNEW holds the gradient of the quadratic model at XOPT+D. It is updated */
1197 /* when D is updated. */
1198 /* XBDI is a working space vector. For I=1,2,...,N, the element XBDI(I) is */
1199 /* set to -1.0, 0.0, or 1.0, the value being nonzero if and only if the */
1200 /* I-th variable has become fixed at a bound, the bound being SL(I) or */
1201 /* SU(I) in the case XBDI(I)=-1.0 or XBDI(I)=1.0, respectively. This */
1202 /* information is accumulated during the construction of XNEW. */
1203 /* The arrays S, HS and HRED are also used for working space. They hold the */
1204 /* current search direction, and the changes in the gradient of Q along S */
1205 /* and the reduced D, respectively, where the reduced D is the same as D, */
1206 /* except that the components of the fixed variables are zero. */
1207 /* DSQ will be set to the square of the length of XNEW-XOPT. */
1208 /* CRVMIN is set to zero if D reaches the trust region boundary. Otherwise */
1209 /* it is set to the least curvature of H that occurs in the conjugate */
1210 /* gradient searches that are not restricted by any constraints. The */
1211 /* value CRVMIN=-1.0D0 is set, however, if all of these searches are */
1214 /* A version of the truncated conjugate gradient is applied. If a line */
1215 /* search is restricted by a constraint, then the procedure is restarted, */
1216 /* the values of the variables that are at their bounds being fixed. If */
1217 /* the trust region boundary is reached, then further changes may be made */
1218 /* to D, each one being in the two dimensional space that is spanned */
1219 /* by the current D and the gradient of Q at XOPT+D, staying on the trust */
1220 /* region boundary. Termination occurs when the reduction in Q seems to */
1221 /* be close to the greatest reduction that can be achieved. */
1223 /* Set some constants. */
1225 /* Parameter adjustments */
1227 xpt_offset = 1 + xpt_dim1;
1249 /* The sign of GOPT(I) gives the sign of the change to the I-th variable */
1250 /* that will reduce Q from its value at XOPT. Thus XBDI(I) shows whether */
1251 /* or not to fix the I-th variable at one of its bounds initially, with */
1252 /* NACT being set to the number of fixed variables. D and GNEW are also */
1253 /* set for the first iteration. DELSQ is the upper bound on the sum of */
1254 /* squares of the free variables. QRED is the reduction in Q so far. */
1260 for (i__ = 1; i__ <= i__1; ++i__) {
1262 if (xopt[i__] <= sl[i__]) {
1263 if (gopt[i__] >= zero) {
1266 } else if (xopt[i__] >= su[i__]) {
1267 if (gopt[i__] <= zero) {
1271 if (xbdi[i__] != zero) {
1276 gnew[i__] = gopt[i__];
1278 delsq = *delta * *delta;
1282 /* Set the next search direction of the conjugate gradient method. It is */
1283 /* the steepest descent direction initially and when the iterations are */
1284 /* restarted because a variable has just been fixed by a bound, and of */
1285 /* course the components of the fixed variables are zero. ITERMAX is an */
1286 /* upper bound on the indices of the conjugate gradient iterations. */
1293 for (i__ = 1; i__ <= i__1; ++i__) {
1294 if (xbdi[i__] != zero) {
1296 } else if (beta == zero) {
1297 s[i__] = -gnew[i__];
1299 s[i__] = beta * s[i__] - gnew[i__];
1302 /* Computing 2nd power */
1304 stepsq += d__1 * d__1;
1306 if (stepsq == zero) {
1311 itermax = iterc + *n - nact;
1313 if (gredsq * delsq <= qred * 1e-4 * qred) {
1317 /* Multiply the search direction by the second derivative matrix of Q and */
1318 /* calculate some scalars for the choice of steplength. Then set BLEN to */
1319 /* the length of the the step to the trust region boundary and STPLEN to */
1320 /* the steplength, ignoring the simple bounds. */
1328 for (i__ = 1; i__ <= i__1; ++i__) {
1329 if (xbdi[i__] == zero) {
1330 /* Computing 2nd power */
1332 resid -= d__1 * d__1;
1333 ds += s[i__] * d__[i__];
1334 shs += s[i__] * hs[i__];
1338 if (resid <= zero) {
1341 temp = sqrt(stepsq * resid + ds * ds);
1343 blen = (temp - ds) / stepsq;
1345 blen = resid / (temp + ds);
1350 d__1 = blen, d__2 = gredsq / shs;
1351 stplen = min(d__1,d__2);
1354 /* Reduce STPLEN if necessary in order to preserve the simple bounds, */
1355 /* letting IACT be the index of the new constrained variable. */
1359 for (i__ = 1; i__ <= i__1; ++i__) {
1360 if (s[i__] != zero) {
1361 xsum = xopt[i__] + d__[i__];
1362 if (s[i__] > zero) {
1363 temp = (su[i__] - xsum) / s[i__];
1365 temp = (sl[i__] - xsum) / s[i__];
1367 if (temp < stplen) {
1375 /* Update CRVMIN, GNEW and D. Set SDEC to the decrease that occurs in Q. */
1378 if (stplen > zero) {
1380 temp = shs / stepsq;
1381 if (iact == 0 && temp > zero) {
1382 *crvmin = min(*crvmin,temp);
1383 if (*crvmin == onemin) {
1390 for (i__ = 1; i__ <= i__1; ++i__) {
1391 gnew[i__] += stplen * hs[i__];
1392 if (xbdi[i__] == zero) {
1393 /* Computing 2nd power */
1395 gredsq += d__1 * d__1;
1398 d__[i__] += stplen * s[i__];
1401 d__1 = stplen * (ggsav - half * stplen * shs);
1402 sdec = max(d__1,zero);
1406 /* Restart the conjugate gradient method if it has hit a new bound. */
1411 if (s[iact] < zero) {
1412 xbdi[iact] = onemin;
1414 /* Computing 2nd power */
1416 delsq -= d__1 * d__1;
1417 if (delsq <= zero) {
1423 /* If STPLEN is less than BLEN, then either apply another conjugate */
1424 /* gradient iteration or RETURN. */
1426 if (stplen < blen) {
1427 if (iterc == itermax) {
1430 if (sdec <= qred * .01) {
1433 beta = gredsq / ggsav;
1439 /* Prepare for the alternative iteration by calculating some scalars */
1440 /* and by multiplying the reduced D by the second derivative matrix of */
1441 /* Q, where S holds the reduced D in the call of GGMULT. */
1444 if (nact >= *n - 1) {
1451 for (i__ = 1; i__ <= i__1; ++i__) {
1452 if (xbdi[i__] == zero) {
1453 /* Computing 2nd power */
1455 dredsq += d__1 * d__1;
1456 dredg += d__[i__] * gnew[i__];
1457 /* Computing 2nd power */
1459 gredsq += d__1 * d__1;
1469 /* Let the search direction S be a linear combination of the reduced D */
1470 /* and the reduced G that is orthogonal to the reduced D. */
1474 temp = gredsq * dredsq - dredg * dredg;
1475 if (temp <= qred * 1e-4 * qred) {
1480 for (i__ = 1; i__ <= i__1; ++i__) {
1481 if (xbdi[i__] == zero) {
1482 s[i__] = (dredg * d__[i__] - dredsq * gnew[i__]) / temp;
1490 /* By considering the simple bounds on the variables, calculate an upper */
1491 /* bound on the tangent of half the angle of the alternative iteration, */
1492 /* namely ANGBD, except that, if already a free variable has reached a */
1493 /* bound, there is a branch back to label 100 after fixing that variable. */
1498 for (i__ = 1; i__ <= i__1; ++i__) {
1499 if (xbdi[i__] == zero) {
1500 tempa = xopt[i__] + d__[i__] - sl[i__];
1501 tempb = su[i__] - xopt[i__] - d__[i__];
1502 if (tempa <= zero) {
1506 } else if (tempb <= zero) {
1512 /* Computing 2nd power */
1514 /* Computing 2nd power */
1516 ssq = d__1 * d__1 + d__2 * d__2;
1517 /* Computing 2nd power */
1518 d__1 = xopt[i__] - sl[i__];
1519 temp = ssq - d__1 * d__1;
1521 temp = sqrt(temp) - s[i__];
1522 if (angbd * temp > tempa) {
1523 angbd = tempa / temp;
1528 /* Computing 2nd power */
1529 d__1 = su[i__] - xopt[i__];
1530 temp = ssq - d__1 * d__1;
1532 temp = sqrt(temp) + s[i__];
1533 if (angbd * temp > tempb) {
1534 angbd = tempb / temp;
1543 /* Calculate HHD and some curvatures for the alternative iteration. */
1551 for (i__ = 1; i__ <= i__1; ++i__) {
1552 if (xbdi[i__] == zero) {
1553 shs += s[i__] * hs[i__];
1554 dhs += d__[i__] * hs[i__];
1555 dhd += d__[i__] * hred[i__];
1560 /* Seek the greatest reduction in Q for a range of equally spaced values */
1561 /* of ANGT in [0,ANGBD], where ANGT is the tangent of half the angle of */
1562 /* the alternative iteration. */
1567 iu = (int) (angbd * 17. + 3.1);
1569 for (i__ = 1; i__ <= i__1; ++i__) {
1570 angt = angbd * (double) i__ / (double) iu;
1571 sth = (angt + angt) / (one + angt * angt);
1572 temp = shs + angt * (angt * dhd - dhs - dhs);
1573 rednew = sth * (angt * dredg - sredg - half * sth * temp);
1574 if (rednew > redmax) {
1578 } else if (i__ == isav + 1) {
1585 /* Return if the reduction is zero. Otherwise, set the sine and cosine */
1586 /* of the angle of the alternative iteration, and calculate SDEC. */
1592 temp = (rdnext - rdprev) / (redmax + redmax - rdprev - rdnext);
1593 angt = angbd * ((double) isav + half * temp) / (double) iu;
1595 cth = (one - angt * angt) / (one + angt * angt);
1596 sth = (angt + angt) / (one + angt * angt);
1597 temp = shs + angt * (angt * dhd - dhs - dhs);
1598 sdec = sth * (angt * dredg - sredg - half * sth * temp);
1603 /* Update GNEW, D and HRED. If the angle of the alternative iteration */
1604 /* is restricted by a bound on a free variable, that variable is fixed */
1610 for (i__ = 1; i__ <= i__1; ++i__) {
1611 gnew[i__] = gnew[i__] + (cth - one) * hred[i__] + sth * hs[i__];
1612 if (xbdi[i__] == zero) {
1613 d__[i__] = cth * d__[i__] + sth * s[i__];
1614 dredg += d__[i__] * gnew[i__];
1615 /* Computing 2nd power */
1617 gredsq += d__1 * d__1;
1620 hred[i__] = cth * hred[i__] + sth * hs[i__];
1623 if (iact > 0 && isav == iu) {
1629 /* If SDEC is sufficiently small, then RETURN after setting XNEW to */
1630 /* XOPT+D, giving careful attention to the bounds. */
1632 if (sdec > qred * .01) {
1638 for (i__ = 1; i__ <= i__1; ++i__) {
1641 d__3 = xopt[i__] + d__[i__], d__4 = su[i__];
1642 d__1 = min(d__3,d__4), d__2 = sl[i__];
1643 xnew[i__] = max(d__1,d__2);
1644 if (xbdi[i__] == onemin) {
1645 xnew[i__] = sl[i__];
1647 if (xbdi[i__] == one) {
1648 xnew[i__] = su[i__];
1650 d__[i__] = xnew[i__] - xopt[i__];
1652 /* Computing 2nd power */
1654 *dsq += d__1 * d__1;
1657 /* The following instructions multiply the current S-vector by the second */
1658 /* derivative matrix of the quadratic model, putting the product in HS. */
1659 /* They are reached from three different parts of the software above and */
1660 /* they can be regarded as an external subroutine. */
1665 for (j = 1; j <= i__1; ++j) {
1668 for (i__ = 1; i__ <= i__2; ++i__) {
1671 hs[j] += hq[ih] * s[i__];
1674 hs[i__] += hq[ih] * s[j];
1678 for (k = 1; k <= i__2; ++k) {
1679 if (pq[k] != zero) {
1682 for (j = 1; j <= i__1; ++j) {
1684 temp += xpt[k + j * xpt_dim1] * s[j];
1688 for (i__ = 1; i__ <= i__1; ++i__) {
1690 hs[i__] += temp * xpt[k + i__ * xpt_dim1];
1695 if (*crvmin != zero) {
1698 if (iterc > itcsav) {
1702 for (i__ = 1; i__ <= i__2; ++i__) {
1704 hred[i__] = hs[i__];
1709 static nlopt_result prelim_(int *n, int *npt, double *x,
1710 const double *xl, const double *xu, double *rhobeg,
1711 nlopt_stopping *stop,
1712 bobyqa_func calfun, void *calfun_data,
1713 double *xbase, double *xpt, double *fval,
1714 double *gopt, double *hq, double *pq, double *bmat,
1715 double *zmat, int *ndim, double *sl, double *su,
1718 /* System generated locals */
1719 int xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1,
1720 zmat_offset, i__1, i__2;
1721 double d__1, d__2, d__3, d__4;
1723 /* Local variables */
1725 int i__, j, k, ih, np, nfm;
1728 double two, fbeg, diff, half, temp, zero, recip, stepa, stepb;
1734 /* The arguments N, NPT, X, XL, XU, RHOBEG, and MAXFUN are the */
1735 /* same as the corresponding arguments in SUBROUTINE BOBYQA. */
1736 /* The arguments XBASE, XPT, FVAL, HQ, PQ, BMAT, ZMAT, NDIM, SL and SU */
1737 /* are the same as the corresponding arguments in BOBYQB, the elements */
1738 /* of SL and SU being set in BOBYQA. */
1739 /* GOPT is usually the gradient of the quadratic model at XOPT+XBASE, but */
1740 /* it is set by PRELIM to the gradient of the quadratic model at XBASE. */
1741 /* If XOPT is nonzero, BOBYQB will change it to its usual value later. */
1742 /* NF is maintaned as the number of calls of CALFUN so far. */
1743 /* KOPT will be such that the least calculated value of F so far is at */
1744 /* the point XPT(KOPT,.)+XBASE in the space of the variables. */
1746 /* SUBROUTINE PRELIM sets the elements of XBASE, XPT, FVAL, GOPT, HQ, PQ, */
1747 /* BMAT and ZMAT for the first iteration, and it maintains the values of */
1748 /* NF and KOPT. The vector X is also changed by PRELIM. */
1750 /* Set some constants. */
1752 /* Parameter adjustments */
1754 zmat_offset = 1 + zmat_dim1;
1755 zmat -= zmat_offset;
1757 xpt_offset = 1 + xpt_dim1;
1768 bmat_offset = 1 + bmat_dim1;
1769 bmat -= bmat_offset;
1778 rhosq = *rhobeg * *rhobeg;
1779 recip = one / rhosq;
1782 /* Set XBASE to the initial vector of variables, and set the initial */
1783 /* elements of XPT, BMAT, HQ, PQ and ZMAT to zero. */
1786 for (j = 1; j <= i__1; ++j) {
1789 for (k = 1; k <= i__2; ++k) {
1791 xpt[k + j * xpt_dim1] = zero;
1794 for (i__ = 1; i__ <= i__2; ++i__) {
1796 bmat[i__ + j * bmat_dim1] = zero;
1800 for (ih = 1; ih <= i__2; ++ih) {
1805 for (k = 1; k <= i__2; ++k) {
1808 for (j = 1; j <= i__1; ++j) {
1810 zmat[k + j * zmat_dim1] = zero;
1814 /* Begin the initialization procedure. NF becomes one more than the number */
1815 /* of function values so far. The coordinates of the displacement of the */
1816 /* next initial interpolation point from XBASE are set in XPT(NF+1,.). */
1823 if (nfm <= *n << 1) {
1824 if (nfm >= 1 && nfm <= *n) {
1826 if (su[nfm] == zero) {
1829 xpt[nf + nfm * xpt_dim1] = stepa;
1830 } else if (nfm > *n) {
1831 stepa = xpt[nf - *n + nfx * xpt_dim1];
1833 if (sl[nfx] == zero) {
1835 d__1 = two * *rhobeg, d__2 = su[nfx];
1836 stepb = min(d__1,d__2);
1838 if (su[nfx] == zero) {
1840 d__1 = -two * *rhobeg, d__2 = sl[nfx];
1841 stepb = max(d__1,d__2);
1843 xpt[nf + nfx * xpt_dim1] = stepb;
1846 itemp = (nfm - np) / *n;
1847 jpt = nfm - itemp * *n - *n;
1854 xpt[nf + ipt * xpt_dim1] = xpt[ipt + 1 + ipt * xpt_dim1];
1855 xpt[nf + jpt * xpt_dim1] = xpt[jpt + 1 + jpt * xpt_dim1];
1858 /* Calculate the next value of F. The least function value so far and */
1859 /* its index are required. */
1862 for (j = 1; j <= i__1; ++j) {
1865 d__3 = xl[j], d__4 = xbase[j] + xpt[nf + j * xpt_dim1];
1866 d__1 = max(d__3,d__4), d__2 = xu[j];
1867 x[j] = min(d__1,d__2);
1868 if (xpt[nf + j * xpt_dim1] == sl[j]) {
1871 if (xpt[nf + j * xpt_dim1] == su[j]) {
1877 f = calfun(*n, &x[1], calfun_data);
1882 } else if (f < fval[*kopt]) {
1886 /* Set the nonzero initial elements of BMAT and the quadratic model in the */
1887 /* cases when NF is at most 2*N+1. If NF exceeds N+1, then the positions */
1888 /* of the NF-th and (NF-N)-th interpolation points may be switched, in */
1889 /* order that the function value at the first of them contributes to the */
1890 /* off-diagonal second derivative terms of the initial quadratic model. */
1892 if (nf <= (*n << 1) + 1) {
1893 if (nf >= 2 && nf <= *n + 1) {
1894 gopt[nfm] = (f - fbeg) / stepa;
1895 if (*npt < nf + *n) {
1896 bmat[nfm * bmat_dim1 + 1] = -one / stepa;
1897 bmat[nf + nfm * bmat_dim1] = one / stepa;
1898 bmat[*npt + nfm + nfm * bmat_dim1] = -half * rhosq;
1900 } else if (nf >= *n + 2) {
1901 ih = nfx * (nfx + 1) / 2;
1902 temp = (f - fbeg) / stepb;
1903 diff = stepb - stepa;
1904 hq[ih] = two * (temp - gopt[nfx]) / diff;
1905 gopt[nfx] = (gopt[nfx] * stepb - temp * stepa) / diff;
1906 if (stepa * stepb < zero) {
1907 if (f < fval[nf - *n]) {
1908 fval[nf] = fval[nf - *n];
1913 xpt[nf - *n + nfx * xpt_dim1] = stepb;
1914 xpt[nf + nfx * xpt_dim1] = stepa;
1917 bmat[nfx * bmat_dim1 + 1] = -(stepa + stepb) / (stepa * stepb);
1918 bmat[nf + nfx * bmat_dim1] = -half / xpt[nf - *n + nfx *
1920 bmat[nf - *n + nfx * bmat_dim1] = -bmat[nfx * bmat_dim1 + 1] -
1921 bmat[nf + nfx * bmat_dim1];
1922 zmat[nfx * zmat_dim1 + 1] = sqrt(two) / (stepa * stepb);
1923 zmat[nf + nfx * zmat_dim1] = sqrt(half) / rhosq;
1924 zmat[nf - *n + nfx * zmat_dim1] = -zmat[nfx * zmat_dim1 + 1] -
1925 zmat[nf + nfx * zmat_dim1];
1928 /* Set the off-diagonal second derivatives of the Lagrange functions and */
1929 /* the initial quadratic model. */
1932 ih = ipt * (ipt - 1) / 2 + jpt;
1933 zmat[nfx * zmat_dim1 + 1] = recip;
1934 zmat[nf + nfx * zmat_dim1] = recip;
1935 zmat[ipt + 1 + nfx * zmat_dim1] = -recip;
1936 zmat[jpt + 1 + nfx * zmat_dim1] = -recip;
1937 temp = xpt[nf + ipt * xpt_dim1] * xpt[nf + jpt * xpt_dim1];
1938 hq[ih] = (fbeg - fval[ipt + 1] - fval[jpt + 1] + f) / temp;
1940 if (f < stop->minf_max) return NLOPT_MINF_MAX_REACHED;
1941 else if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
1942 else if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
1946 return NLOPT_SUCCESS;
1949 static nlopt_result bobyqb_(int *n, int *npt, double *x,
1950 const double *xl, const double *xu, double *rhobeg, double *
1952 nlopt_stopping *stop,
1953 bobyqa_func calfun, void *calfun_data,
1956 double *xpt, double *fval, double *xopt, double *gopt,
1957 double *hq, double *pq, double *bmat, double *zmat,
1958 int *ndim, double *sl, double *su, double *xnew,
1959 double *xalt, double *d__, double *vlag, double *w)
1961 /* System generated locals */
1962 int xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1,
1963 zmat_offset, i__1, i__2, i__3;
1964 double d__1, d__2, d__3, d__4;
1966 /* Local variables */
1968 int i__, j, k, ih, jj, nh, ip, jp;
1971 double den, one, ten, dsq, rho, sum, two, diff, half, beta, gisq;
1973 double temp, suma, sumb, bsum, fopt;
1977 double gqsq, dist, sumw, sumz, diffa, diffb, diffc, hdiag;
1979 double alpha, delta, adelt, denom, fsave, bdtol, delsq;
1981 double ratio, dnorm, vquad, pqold, tenth;
1983 double sumpq, scaden;
1984 double errbig, cauchy, fracsq, biglsq, densav;
1986 double crvmin, frhosq;
1991 nlopt_result rc = NLOPT_SUCCESS, rc2;
1993 /* The arguments N, NPT, X, XL, XU, RHOBEG, RHOEND, and MAXFUN */
1994 /* are identical to the corresponding arguments in SUBROUTINE BOBYQA. */
1995 /* XBASE holds a shift of origin that should reduce the contributions */
1996 /* from rounding errors to values of the model and Lagrange functions. */
1997 /* XPT is a two-dimensional array that holds the coordinates of the */
1998 /* interpolation points relative to XBASE. */
1999 /* FVAL holds the values of F at the interpolation points. */
2000 /* XOPT is set to the displacement from XBASE of the trust region centre. */
2001 /* GOPT holds the gradient of the quadratic model at XBASE+XOPT. */
2002 /* HQ holds the explicit second derivatives of the quadratic model. */
2003 /* PQ contains the parameters of the implicit second derivatives of the */
2004 /* quadratic model. */
2005 /* BMAT holds the last N columns of H. */
2006 /* ZMAT holds the factorization of the leading NPT by NPT submatrix of H, */
2007 /* this factorization being ZMAT times ZMAT^T, which provides both the */
2008 /* correct rank and positive semi-definiteness. */
2009 /* NDIM is the first dimension of BMAT and has the value NPT+N. */
2010 /* SL and SU hold the differences XL-XBASE and XU-XBASE, respectively. */
2011 /* All the components of every XOPT are going to satisfy the bounds */
2012 /* SL(I) .LEQ. XOPT(I) .LEQ. SU(I), with appropriate equalities when */
2013 /* XOPT is on a constraint boundary. */
2014 /* XNEW is chosen by SUBROUTINE TRSBOX or ALTMOV. Usually XBASE+XNEW is the */
2015 /* vector of variables for the next call of CALFUN. XNEW also satisfies */
2016 /* the SL and SU constraints in the way that has just been mentioned. */
2017 /* XALT is an alternative to XNEW, chosen by ALTMOV, that may replace XNEW */
2018 /* in order to increase the denominator in the updating of UPDATE. */
2019 /* D is reserved for a trial step from XOPT, which is usually XNEW-XOPT. */
2020 /* VLAG contains the values of the Lagrange functions at a new point X. */
2021 /* They are part of a product that requires VLAG to be of length NDIM. */
2022 /* W is a one-dimensional array that is used for working space. Its length */
2023 /* must be at least 3*NDIM = 3*(NPT+N). */
2025 /* Set some constants. */
2027 /* Parameter adjustments */
2029 zmat_offset = 1 + zmat_dim1;
2030 zmat -= zmat_offset;
2032 xpt_offset = 1 + xpt_dim1;
2044 bmat_offset = 1 + bmat_dim1;
2045 bmat -= bmat_offset;
2065 /* The call of PRELIM sets the elements of XBASE, XPT, FVAL, GOPT, HQ, PQ, */
2066 /* BMAT and ZMAT for the first iteration, with the corresponding values of */
2067 /* of NF and KOPT, which are the number of calls of CALFUN so far and the */
2068 /* index of the interpolation point at the trust region centre. Then the */
2069 /* initial XOPT is set too. The branch to label 720 occurs if MAXFUN is */
2070 /* less than NPT. GOPT will be updated if KOPT is different from KBASE. */
2072 rc2 = prelim_(n, npt, &x[1], &xl[1], &xu[1], rhobeg,
2073 stop, calfun, calfun_data,
2074 &xbase[1], &xpt[xpt_offset], &fval[1], &gopt[1], &hq[1], &pq[1], &bmat[
2075 bmat_offset], &zmat[zmat_offset], ndim, &sl[1], &su[1], &kopt);
2078 for (i__ = 1; i__ <= i__1; ++i__) {
2079 xopt[i__] = xpt[kopt + i__ * xpt_dim1];
2081 /* Computing 2nd power */
2083 xoptsq += d__1 * d__1;
2086 if (rc2 != NLOPT_SUCCESS) {
2092 /* Complete the settings that are required for the iterative procedure. */
2096 nresc = stop->nevals;
2101 nfsav = stop->nevals;
2103 /* Update GOPT if necessary before the first iteration and after each */
2104 /* call of RESCUE that makes a call of CALFUN. */
2107 if (kopt != kbase) {
2110 for (j = 1; j <= i__1; ++j) {
2112 for (i__ = 1; i__ <= i__2; ++i__) {
2115 gopt[j] += hq[ih] * xopt[i__];
2118 gopt[i__] += hq[ih] * xopt[j];
2121 if (stop->nevals > *npt) {
2123 for (k = 1; k <= i__2; ++k) {
2126 for (j = 1; j <= i__1; ++j) {
2128 temp += xpt[k + j * xpt_dim1] * xopt[j];
2130 temp = pq[k] * temp;
2132 for (i__ = 1; i__ <= i__1; ++i__) {
2134 gopt[i__] += temp * xpt[k + i__ * xpt_dim1];
2140 /* Generate the next point in the trust region that provides a small value */
2141 /* of the quadratic model subject to the constraints on the variables. */
2142 /* The int NTRITS is set to the number "trust region" iterations that */
2143 /* have occurred since the last "alternative" iteration. If the length */
2144 /* of XNEW-XOPT is less than HALF*RHO, however, then there is a branch to */
2145 /* label 650 or 680 with NTRITS=-1, instead of calculating F at XNEW. */
2148 trsbox_(n, npt, &xpt[xpt_offset], &xopt[1], &gopt[1], &hq[1], &pq[1], &sl[
2149 1], &su[1], &delta, &xnew[1], &d__[1], &w[1], &w[np], &w[np + *n],
2150 &w[np + (*n << 1)], &w[np + *n * 3], &dsq, &crvmin);
2152 d__1 = delta, d__2 = sqrt(dsq);
2153 dnorm = min(d__1,d__2);
2154 if (dnorm < half * rho) {
2156 /* Computing 2nd power */
2158 distsq = d__1 * d__1;
2159 if (stop->nevals <= nfsav + 2) {
2163 /* The following choice between labels 650 and 680 depends on whether or */
2164 /* not our work with the current RHO seems to be complete. Either RHO is */
2165 /* decreased or termination occurs if the errors in the quadratic model at */
2166 /* the last three interpolation points compare favourably with predictions */
2167 /* of likely improvements to the model within distance HALF*RHO of XOPT. */
2170 d__1 = max(diffa,diffb);
2171 errbig = max(d__1,diffc);
2172 frhosq = rho * .125 * rho;
2173 if (crvmin > zero && errbig > frhosq * crvmin) {
2176 bdtol = errbig / rho;
2178 for (j = 1; j <= i__1; ++j) {
2180 if (xnew[j] == sl[j]) {
2183 if (xnew[j] == su[j]) {
2186 if (bdtest < bdtol) {
2187 curv = hq[(j + j * j) / 2];
2189 for (k = 1; k <= i__2; ++k) {
2191 /* Computing 2nd power */
2192 d__1 = xpt[k + j * xpt_dim1];
2193 curv += pq[k] * (d__1 * d__1);
2195 bdtest += half * curv * rho;
2196 if (bdtest < bdtol) {
2206 /* Severe cancellation is likely to occur if XOPT is too far from XBASE. */
2207 /* If the following test holds, then XBASE is shifted so that XOPT becomes */
2208 /* zero. The appropriate changes are made to BMAT and to the second */
2209 /* derivatives of the current model, beginning with the changes to BMAT */
2210 /* that do not depend on ZMAT. VLAG is used temporarily for working space. */
2213 if (dsq <= xoptsq * .001) {
2214 fracsq = xoptsq * .25;
2217 for (k = 1; k <= i__1; ++k) {
2219 sum = -half * xoptsq;
2221 for (i__ = 1; i__ <= i__2; ++i__) {
2223 sum += xpt[k + i__ * xpt_dim1] * xopt[i__];
2226 temp = fracsq - half * sum;
2228 for (i__ = 1; i__ <= i__2; ++i__) {
2229 w[i__] = bmat[k + i__ * bmat_dim1];
2230 vlag[i__] = sum * xpt[k + i__ * xpt_dim1] + temp * xopt[i__];
2233 for (j = 1; j <= i__3; ++j) {
2235 bmat[ip + j * bmat_dim1] = bmat[ip + j * bmat_dim1] + w[
2236 i__] * vlag[j] + vlag[i__] * w[j];
2241 /* Then the revisions of BMAT that depend on ZMAT are calculated. */
2244 for (jj = 1; jj <= i__3; ++jj) {
2248 for (k = 1; k <= i__2; ++k) {
2249 sumz += zmat[k + jj * zmat_dim1];
2250 vlag[k] = w[*npt + k] * zmat[k + jj * zmat_dim1];
2255 for (j = 1; j <= i__2; ++j) {
2256 sum = (fracsq * sumz - half * sumw) * xopt[j];
2258 for (k = 1; k <= i__1; ++k) {
2260 sum += vlag[k] * xpt[k + j * xpt_dim1];
2264 for (k = 1; k <= i__1; ++k) {
2266 bmat[k + j * bmat_dim1] += sum * zmat[k + jj * zmat_dim1];
2270 for (i__ = 1; i__ <= i__1; ++i__) {
2274 for (j = 1; j <= i__2; ++j) {
2276 bmat[ip + j * bmat_dim1] += temp * w[j];
2281 /* The following instructions complete the shift, including the changes */
2282 /* to the second derivative parameters of the quadratic model. */
2286 for (j = 1; j <= i__2; ++j) {
2287 w[j] = -half * sumpq * xopt[j];
2289 for (k = 1; k <= i__1; ++k) {
2290 w[j] += pq[k] * xpt[k + j * xpt_dim1];
2292 xpt[k + j * xpt_dim1] -= xopt[j];
2295 for (i__ = 1; i__ <= i__1; ++i__) {
2297 hq[ih] = hq[ih] + w[i__] * xopt[j] + xopt[i__] * w[j];
2299 bmat[*npt + i__ + j * bmat_dim1] = bmat[*npt + j + i__ *
2304 for (i__ = 1; i__ <= i__1; ++i__) {
2305 xbase[i__] += xopt[i__];
2306 xnew[i__] -= xopt[i__];
2307 sl[i__] -= xopt[i__];
2308 su[i__] -= xopt[i__];
2319 /* XBASE is also moved to XOPT by a call of RESCUE. This calculation is */
2320 /* more expensive than the previous shift, because new matrices BMAT and */
2321 /* ZMAT are generated from scratch, which may include the replacement of */
2322 /* interpolation points whose positions seem to be causing near linear */
2323 /* dependence in the interpolation conditions. Therefore RESCUE is called */
2324 /* only if rounding errors have reduced by at least a factor of two the */
2325 /* denominator of the formula for updating the H matrix. It provides a */
2326 /* useful safeguard, but is not invoked in most applications of BOBYQA. */
2329 nfsav = stop->nevals;
2331 rc2 = rescue_(n, npt, &xl[1], &xu[1],
2332 stop, calfun, calfun_data,
2333 &xbase[1], &xpt[xpt_offset], &fval[1], &xopt[1], &gopt[1],
2334 &hq[1], &pq[1], &bmat[bmat_offset], &zmat[zmat_offset], ndim,
2335 &sl[1], &su[1], &delta, &kopt, &vlag[1],
2336 &w[1], &w[*n + np], &w[*ndim + np]);
2338 /* XOPT is updated now in case the branch below to label 720 is taken. */
2339 /* Any updating of GOPT occurs after the branch below to label 20, which */
2340 /* leads to a trust region iteration as does the branch to label 60. */
2343 if (kopt != kbase) {
2345 for (i__ = 1; i__ <= i__1; ++i__) {
2346 xopt[i__] = xpt[kopt + i__ * xpt_dim1];
2348 /* Computing 2nd power */
2350 xoptsq += d__1 * d__1;
2353 if (rc2 != NLOPT_SUCCESS) {
2357 nresc = stop->nevals;
2358 if (nfsav < stop->nevals) {
2359 nfsav = stop->nevals;
2366 /* Pick two alternative vectors of variables, relative to XBASE, that */
2367 /* are suitable as new positions of the KNEW-th interpolation point. */
2368 /* Firstly, XNEW is set to the point on a line through XOPT and another */
2369 /* interpolation point that minimizes the predicted value of the next */
2370 /* denominator, subject to ||XNEW - XOPT|| .LEQ. ADELT and to the SL */
2371 /* and SU bounds. Secondly, XALT is set to the best feasible point on */
2372 /* a constrained version of the Cauchy step of the KNEW-th Lagrange */
2373 /* function, the corresponding value of the square of this function */
2374 /* being returned in CAUCHY. The choice between these alternatives is */
2375 /* going to be made when the denominator is calculated. */
2378 altmov_(n, npt, &xpt[xpt_offset], &xopt[1], &bmat[bmat_offset], &zmat[
2379 zmat_offset], ndim, &sl[1], &su[1], &kopt, &knew, &adelt, &xnew[1]
2380 , &xalt[1], &alpha, &cauchy, &w[1], &w[np], &w[*ndim + 1]);
2382 for (i__ = 1; i__ <= i__1; ++i__) {
2384 d__[i__] = xnew[i__] - xopt[i__];
2387 /* Calculate VLAG and BETA for the current choice of D. The scalar */
2388 /* product of D with XPT(K,.) is going to be held in W(NPT+K) for */
2389 /* use when VQUAD is calculated. */
2393 for (k = 1; k <= i__1; ++k) {
2398 for (j = 1; j <= i__2; ++j) {
2399 suma += xpt[k + j * xpt_dim1] * d__[j];
2400 sumb += xpt[k + j * xpt_dim1] * xopt[j];
2402 sum += bmat[k + j * bmat_dim1] * d__[j];
2404 w[k] = suma * (half * suma + sumb);
2411 for (jj = 1; jj <= i__1; ++jj) {
2414 for (k = 1; k <= i__2; ++k) {
2416 sum += zmat[k + jj * zmat_dim1] * w[k];
2420 for (k = 1; k <= i__2; ++k) {
2422 vlag[k] += sum * zmat[k + jj * zmat_dim1];
2429 for (j = 1; j <= i__2; ++j) {
2430 /* Computing 2nd power */
2435 for (k = 1; k <= i__1; ++k) {
2437 sum += w[k] * bmat[k + j * bmat_dim1];
2439 bsum += sum * d__[j];
2442 for (i__ = 1; i__ <= i__1; ++i__) {
2444 sum += bmat[jp + i__ * bmat_dim1] * d__[i__];
2447 bsum += sum * d__[j];
2449 dx += d__[j] * xopt[j];
2451 beta = dx * dx + dsq * (xoptsq + dx + dx + half * dsq) + beta - bsum;
2454 /* If NTRITS is zero, the denominator may be increased by replacing */
2455 /* the step D of ALTMOV by a Cauchy step. Then RESCUE may be called if */
2456 /* rounding errors have damaged the chosen denominator. */
2459 /* Computing 2nd power */
2461 denom = d__1 * d__1 + alpha * beta;
2462 if (denom < cauchy && cauchy > zero) {
2464 for (i__ = 1; i__ <= i__2; ++i__) {
2465 xnew[i__] = xalt[i__];
2467 d__[i__] = xnew[i__] - xopt[i__];
2472 /* Computing 2nd power */
2474 if (denom <= half * (d__1 * d__1)) {
2475 if (stop->nevals > nresc) {
2478 /* Return from BOBYQA because of much cancellation in a
2480 rc = NLOPT_ROUNDOFF_LIMITED;
2484 /* Alternatively, if NTRITS is positive, then set KNEW to the index of */
2485 /* the next interpolation point to be deleted to make room for a trust */
2486 /* region step. Again RESCUE may be called if rounding errors have damaged */
2487 /* the chosen denominator, which is the reason for attempting to select */
2488 /* KNEW before calculating the next value of the objective function. */
2491 delsq = delta * delta;
2496 for (k = 1; k <= i__2; ++k) {
2502 for (jj = 1; jj <= i__1; ++jj) {
2504 /* Computing 2nd power */
2505 d__1 = zmat[k + jj * zmat_dim1];
2506 hdiag += d__1 * d__1;
2508 /* Computing 2nd power */
2510 den = beta * hdiag + d__1 * d__1;
2513 for (j = 1; j <= i__1; ++j) {
2515 /* Computing 2nd power */
2516 d__1 = xpt[k + j * xpt_dim1] - xopt[j];
2517 distsq += d__1 * d__1;
2520 /* Computing 2nd power */
2521 d__3 = distsq / delsq;
2522 d__1 = one, d__2 = d__3 * d__3;
2523 temp = max(d__1,d__2);
2524 if (temp * den > scaden) {
2525 scaden = temp * den;
2530 /* Computing 2nd power */
2532 d__1 = biglsq, d__2 = temp * (d__3 * d__3);
2533 biglsq = max(d__1,d__2);
2537 if (scaden <= half * biglsq) {
2538 if (stop->nevals > nresc) {
2541 /* Return from BOBYQA because of much cancellation in a
2543 rc = NLOPT_ROUNDOFF_LIMITED;
2548 /* Put the variables for the next calculation of the objective function */
2549 /* in XNEW, with any adjustments for the bounds. */
2552 /* Calculate the value of the objective function at XBASE+XNEW, unless */
2553 /* the limit on the number of calculations of F has been reached. */
2557 for (i__ = 1; i__ <= i__2; ++i__) {
2560 d__3 = xl[i__], d__4 = xbase[i__] + xnew[i__];
2561 d__1 = max(d__3,d__4), d__2 = xu[i__];
2562 x[i__] = min(d__1,d__2);
2563 if (xnew[i__] == sl[i__]) {
2566 if (xnew[i__] == su[i__]) {
2572 if (nlopt_stop_evals(stop)) rc = NLOPT_MAXEVAL_REACHED;
2573 else if (nlopt_stop_time(stop)) rc = NLOPT_MAXTIME_REACHED;
2574 if (rc != NLOPT_SUCCESS) goto L720;
2577 f = calfun(*n, &x[1], calfun_data);
2580 rc = NLOPT_XTOL_REACHED;
2584 if (f < stop->minf_max) {
2585 return NLOPT_MINF_MAX_REACHED;
2588 /* Use the quadratic model to predict the change in F due to the step D, */
2589 /* and set DIFF to the error of this prediction. */
2595 for (j = 1; j <= i__2; ++j) {
2596 vquad += d__[j] * gopt[j];
2598 for (i__ = 1; i__ <= i__1; ++i__) {
2600 temp = d__[i__] * d__[j];
2605 vquad += hq[ih] * temp;
2609 for (k = 1; k <= i__1; ++k) {
2611 /* Computing 2nd power */
2613 vquad += half * pq[k] * (d__1 * d__1);
2615 diff = f - fopt - vquad;
2620 nfsav = stop->nevals;
2623 /* Pick the next value of DELTA after a trust region step. */
2626 if (vquad >= zero) {
2627 /* Return from BOBYQA because a trust region step has failed
2629 rc = NLOPT_ROUNDOFF_LIMITED; /* or FTOL_REACHED? */
2632 ratio = (f - fopt) / vquad;
2633 if (ratio <= tenth) {
2635 d__1 = half * delta;
2636 delta = min(d__1,dnorm);
2637 } else if (ratio <= .7) {
2639 d__1 = half * delta;
2640 delta = max(d__1,dnorm);
2643 d__1 = half * delta, d__2 = dnorm + dnorm;
2644 delta = max(d__1,d__2);
2646 if (delta <= rho * 1.5) {
2650 /* Recalculate KNEW and DENOM if the new F is less than FOPT. */
2655 delsq = delta * delta;
2660 for (k = 1; k <= i__1; ++k) {
2663 for (jj = 1; jj <= i__2; ++jj) {
2665 /* Computing 2nd power */
2666 d__1 = zmat[k + jj * zmat_dim1];
2667 hdiag += d__1 * d__1;
2669 /* Computing 2nd power */
2671 den = beta * hdiag + d__1 * d__1;
2674 for (j = 1; j <= i__2; ++j) {
2676 /* Computing 2nd power */
2677 d__1 = xpt[k + j * xpt_dim1] - xnew[j];
2678 distsq += d__1 * d__1;
2681 /* Computing 2nd power */
2682 d__3 = distsq / delsq;
2683 d__1 = one, d__2 = d__3 * d__3;
2684 temp = max(d__1,d__2);
2685 if (temp * den > scaden) {
2686 scaden = temp * den;
2692 /* Computing 2nd power */
2694 d__1 = biglsq, d__2 = temp * (d__3 * d__3);
2695 biglsq = max(d__1,d__2);
2697 if (scaden <= half * biglsq) {
2704 /* Update BMAT and ZMAT, so that the KNEW-th interpolation point can be */
2705 /* moved. Also update the second derivative terms of the model. */
2707 update_(n, npt, &bmat[bmat_offset], &zmat[zmat_offset], ndim, &vlag[1], &
2708 beta, &denom, &knew, &w[1]);
2713 for (i__ = 1; i__ <= i__1; ++i__) {
2714 temp = pqold * xpt[knew + i__ * xpt_dim1];
2716 for (j = 1; j <= i__2; ++j) {
2719 hq[ih] += temp * xpt[knew + j * xpt_dim1];
2723 for (jj = 1; jj <= i__2; ++jj) {
2724 temp = diff * zmat[knew + jj * zmat_dim1];
2726 for (k = 1; k <= i__1; ++k) {
2728 pq[k] += temp * zmat[k + jj * zmat_dim1];
2732 /* Include the new interpolation point, and make the changes to GOPT at */
2733 /* the old XOPT that are caused by the updating of the quadratic model. */
2737 for (i__ = 1; i__ <= i__1; ++i__) {
2738 xpt[knew + i__ * xpt_dim1] = xnew[i__];
2740 w[i__] = bmat[knew + i__ * bmat_dim1];
2743 for (k = 1; k <= i__1; ++k) {
2746 for (jj = 1; jj <= i__2; ++jj) {
2748 suma += zmat[knew + jj * zmat_dim1] * zmat[k + jj * zmat_dim1];
2750 if (nlopt_isinf(suma)) {
2751 /* SGJ: detect singularity here (happend if we run
2752 for too many iterations) ... is there another way to recover? */
2753 rc = NLOPT_ROUNDOFF_LIMITED;
2758 for (j = 1; j <= i__2; ++j) {
2760 sumb += xpt[k + j * xpt_dim1] * xopt[j];
2764 for (i__ = 1; i__ <= i__2; ++i__) {
2766 w[i__] += temp * xpt[k + i__ * xpt_dim1];
2770 for (i__ = 1; i__ <= i__2; ++i__) {
2772 gopt[i__] += diff * w[i__];
2775 /* Update XOPT, GOPT and KOPT if the new calculated F is less than FOPT. */
2782 for (j = 1; j <= i__2; ++j) {
2784 /* Computing 2nd power */
2786 xoptsq += d__1 * d__1;
2788 for (i__ = 1; i__ <= i__1; ++i__) {
2791 gopt[j] += hq[ih] * d__[i__];
2794 gopt[i__] += hq[ih] * d__[j];
2798 for (k = 1; k <= i__1; ++k) {
2801 for (j = 1; j <= i__2; ++j) {
2803 temp += xpt[k + j * xpt_dim1] * d__[j];
2805 temp = pq[k] * temp;
2807 for (i__ = 1; i__ <= i__2; ++i__) {
2809 gopt[i__] += temp * xpt[k + i__ * xpt_dim1];
2814 /* Calculate the parameters of the least Frobenius norm interpolant to */
2815 /* the current data, the gradient of this interpolant at XOPT being put */
2816 /* into VLAG(NPT+I), I=1,2,...,N. */
2820 for (k = 1; k <= i__2; ++k) {
2821 vlag[k] = fval[k] - fval[kopt];
2826 for (j = 1; j <= i__2; ++j) {
2829 for (k = 1; k <= i__1; ++k) {
2831 sum += zmat[k + j * zmat_dim1] * vlag[k];
2834 for (k = 1; k <= i__1; ++k) {
2836 w[k] += sum * zmat[k + j * zmat_dim1];
2840 for (k = 1; k <= i__1; ++k) {
2843 for (j = 1; j <= i__2; ++j) {
2845 sum += xpt[k + j * xpt_dim1] * xopt[j];
2854 for (i__ = 1; i__ <= i__1; ++i__) {
2857 for (k = 1; k <= i__2; ++k) {
2859 sum = sum + bmat[k + i__ * bmat_dim1] * vlag[k] + xpt[k + i__
2862 if (xopt[i__] == sl[i__]) {
2864 d__2 = zero, d__3 = gopt[i__];
2865 /* Computing 2nd power */
2866 d__1 = min(d__2,d__3);
2867 gqsq += d__1 * d__1;
2868 /* Computing 2nd power */
2869 d__1 = min(zero,sum);
2870 gisq += d__1 * d__1;
2871 } else if (xopt[i__] == su[i__]) {
2873 d__2 = zero, d__3 = gopt[i__];
2874 /* Computing 2nd power */
2875 d__1 = max(d__2,d__3);
2876 gqsq += d__1 * d__1;
2877 /* Computing 2nd power */
2878 d__1 = max(zero,sum);
2879 gisq += d__1 * d__1;
2881 /* Computing 2nd power */
2883 gqsq += d__1 * d__1;
2887 vlag[*npt + i__] = sum;
2890 /* Test whether to replace the new quadratic model by the least Frobenius */
2891 /* norm interpolant, making the replacement if the test is satisfied. */
2894 if (gqsq < ten * gisq) {
2898 i__1 = max(*npt,nh);
2899 for (i__ = 1; i__ <= i__1; ++i__) {
2901 gopt[i__] = vlag[*npt + i__];
2904 pq[i__] = w[*npt + i__];
2915 /* If a trust region step has provided a sufficient decrease in F, then */
2916 /* branch for another trust region calculation. The case NTRITS=0 occurs */
2917 /* when the new interpolation point was reached by an alternative step. */
2922 if (f <= fopt + tenth * vquad) {
2926 /* Alternatively, find out if the interpolation points are close enough */
2927 /* to the best point so far. */
2930 /* Computing 2nd power */
2932 /* Computing 2nd power */
2934 d__1 = d__3 * d__3, d__2 = d__4 * d__4;
2935 distsq = max(d__1,d__2);
2939 for (k = 1; k <= i__1; ++k) {
2942 for (j = 1; j <= i__2; ++j) {
2944 /* Computing 2nd power */
2945 d__1 = xpt[k + j * xpt_dim1] - xopt[j];
2955 /* If KNEW is positive, then ALTMOV finds alternative new positions for */
2956 /* the KNEW-th interpolation point within distance ADELT of XOPT. It is */
2957 /* reached via label 90. Otherwise, there is a branch to label 60 for */
2958 /* another trust region iteration, unless the calculations with the */
2959 /* current RHO are complete. */
2962 dist = sqrt(distsq);
2965 d__1 = tenth * delta, d__2 = half * dist;
2966 delta = min(d__1,d__2);
2967 if (delta <= rho * 1.5) {
2974 d__2 = tenth * dist;
2975 d__1 = min(d__2,delta);
2976 adelt = max(d__1,rho);
2977 dsq = adelt * adelt;
2986 if (max(delta,dnorm) > rho) {
2990 /* The calculations with the current value of RHO are complete. Pick the */
2991 /* next values of RHO and DELTA. */
2994 if (rho > *rhoend) {
2996 ratio = rho / *rhoend;
2999 } else if (ratio <= 250.) {
3000 rho = sqrt(ratio) * *rhoend;
3004 delta = max(delta,rho);
3006 nfsav = stop->nevals;
3010 /* Return from the calculation, after another Newton-Raphson step, if */
3011 /* it is too short to have been tried before. */
3017 /* originally: if (fval[kopt] <= fsave) -- changed by SGJ, since
3018 this seems like a slight optimization to not update x[]
3019 unnecessarily, at the expense of possibly not returning the
3020 best x[] found so far if the algorithm is stopped suddenly
3021 (e.g. runs out of time) ... it seems safer to execute this
3022 unconditionally, and the efficiency loss seems negligible. */
3025 for (i__ = 1; i__ <= i__1; ++i__) {
3028 d__3 = xl[i__], d__4 = xbase[i__] + xopt[i__];
3029 d__1 = max(d__3,d__4), d__2 = xu[i__];
3030 x[i__] = min(d__1,d__2);
3031 if (xopt[i__] == sl[i__]) {
3034 if (xopt[i__] == su[i__]) {
3045 /**************************************************************************/
3047 nlopt_result bobyqa(int n, int npt, double *x,
3048 const double *xl, const double *xu, double rhobeg,
3049 nlopt_stopping *stop, double *minf,
3050 bobyqa_func calfun, void *calfun_data)
3052 /* System generated locals */
3056 /* Local variables */
3057 int j, id, np, iw, igo, ihq, ixb, ixa, ifv, isl, jsl, ipq, ivl, ixn,
3058 ixo, ixp, isu, jsu, ndim;
3066 /* SGJ, 2009: compute rhoend from NLopt stop info */
3067 rhoend = stop->xtol_rel * (rhobeg);
3068 for (j = 0; j < n; ++j)
3069 if (rhoend < stop->xtol_abs[j])
3070 rhoend = stop->xtol_abs[j];
3073 /* This subroutine seeks the least value of a function of many variables, */
3074 /* by applying a trust region method that forms quadratic models by */
3075 /* interpolation. There is usually some freedom in the interpolation */
3076 /* conditions, which is taken up by minimizing the Frobenius norm of */
3077 /* the change to the second derivative of the model, beginning with the */
3078 /* zero matrix. The values of the variables are constrained by upper and */
3079 /* lower bounds. The arguments of the subroutine are as follows. */
3081 /* N must be set to the number of variables and must be at least two. */
3082 /* NPT is the number of interpolation conditions. Its value must be in */
3083 /* the interval [N+2,(N+1)(N+2)/2]. Choices that exceed 2*N+1 are not */
3085 /* Initial values of the variables must be set in X(1),X(2),...,X(N). They */
3086 /* will be changed to the values that give the least calculated F. */
3087 /* For I=1,2,...,N, XL(I) and XU(I) must provide the lower and upper */
3088 /* bounds, respectively, on X(I). The construction of quadratic models */
3089 /* requires XL(I) to be strictly less than XU(I) for each I. Further, */
3090 /* the contribution to a model from changes to the I-th variable is */
3091 /* damaged severely by rounding errors if XU(I)-XL(I) is too small. */
3092 /* RHOBEG and RHOEND must be set to the initial and final values of a trust */
3093 /* region radius, so both must be positive with RHOEND no greater than */
3094 /* RHOBEG. Typically, RHOBEG should be about one tenth of the greatest */
3095 /* expected change to a variable, while RHOEND should indicate the */
3096 /* accuracy that is required in the final values of the variables. An */
3097 /* error return occurs if any of the differences XU(I)-XL(I), I=1,...,N, */
3098 /* is less than 2*RHOBEG. */
3099 /* MAXFUN must be set to an upper bound on the number of calls of CALFUN. */
3100 /* The array W will be used for working space. Its length must be at least */
3101 /* (NPT+5)*(NPT+N)+3*N*(N+5)/2. */
3103 /* SUBROUTINE CALFUN (N,X,F) has to be provided by the user. It must set */
3104 /* F to the value of the objective function for the current values of the */
3105 /* variables X(1),X(2),...,X(N), which are generated automatically in a */
3106 /* way that satisfies the bounds given in XL and XU. */
3108 /* Return if the value of NPT is unacceptable. */
3110 /* Parameter adjustments */
3117 if (npt < n + 2 || npt > (n + 2) * np / 2) {
3118 /* Return from BOBYQA because NPT is not in the required interval */
3119 return NLOPT_INVALID_ARGS;
3122 /* Partition the working space array, so that different parts of it can */
3123 /* be treated separately during the calculation of BOBYQB. The partition */
3124 /* requires the first (NPT+2)*(NPT+N)+3*N*(N+5)/2 elements of W plus the */
3125 /* space that is taken by the last array in the argument list of BOBYQB. */
3130 ifv = ixp + n * npt;
3134 ipq = ihq + n * np / 2;
3136 izmat = ibmat + ndim * n;
3137 isl = izmat + npt * (npt - np);
3145 w = (double *) malloc(sizeof(double) * ((npt+5)*(npt+n)+3*n*(n+5)/2));
3146 if (!w) return NLOPT_OUT_OF_MEMORY;
3149 /* Return if there is insufficient space between the bounds. Modify the */
3150 /* initial X if necessary in order to avoid conflicts between the bounds */
3151 /* and the construction of the first quadratic model. The lower and upper */
3152 /* bounds on moves from the updated X are set now, in the ISL and ISU */
3153 /* partitions of W, in order to provide useful and exact information about */
3154 /* components of X that become within distance RHOBEG from their bounds. */
3158 for (j = 1; j <= i__1; ++j) {
3159 temp = xu[j] - xl[j];
3160 if (temp < rhobeg + rhobeg) {
3161 /* Return from BOBYQA because one of the differences
3162 XU(I)-XL(I)s is less than 2*RHOBEG. */
3164 return NLOPT_INVALID_ARGS;
3168 w[jsl] = xl[j] - x[j];
3169 w[jsu] = xu[j] - x[j];
3170 if (w[jsl] >= -(rhobeg)) {
3171 if (w[jsl] >= zero) {
3176 x[j] = xl[j] + rhobeg;
3179 d__1 = xu[j] - x[j];
3180 w[jsu] = max(d__1,rhobeg);
3182 } else if (w[jsu] <= rhobeg) {
3183 if (w[jsu] <= zero) {
3188 x[j] = xu[j] - rhobeg;
3190 d__1 = xl[j] - x[j], d__2 = -(rhobeg);
3191 w[jsl] = min(d__1,d__2);
3198 /* Make the call of BOBYQB. */
3200 ret = bobyqb_(&n, &npt, &x[1], &xl[1], &xu[1], &rhobeg, &rhoend,
3201 stop, calfun, calfun_data, minf,
3202 &w[ixb], &w[ixp], &w[ifv], &w[ixo], &w[igo], &w[ihq], &w[ipq],
3203 &w[ibmat], &w[izmat], &ndim, &w[isl], &w[isu], &w[ixn], &w[ixa],
3204 &w[id], &w[ivl], &w[iw]);