bobyqa/bobyqa.c

   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).
   4    See README. */
   5
   6 #include <math.h>
   7 #include <stdlib.h>
   8 #include <string.h>
   9
  10 #include "bobyqa.h"
  11
  12 typedef double (*bobyqa_func)(int n, const double *x, void *func_data);
  13
  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))
  17
  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)
  21 {
  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;
  25
  26     /* Local variables */
  27     int i__, j, k, jl, jp;
  28     double one, tau, temp;
  29     int nptm;
  30     double zero, alpha, tempa, tempb, ztest;
  31
  32
  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. */
  41
  42 /*     Set some constants. */
  43
  44     /* Parameter adjustments */
  45     zmat_dim1 = *npt;
  46     zmat_offset = 1 + zmat_dim1;
  47     zmat -= zmat_offset;
  48     bmat_dim1 = *ndim;
  49     bmat_offset = 1 + bmat_dim1;
  50     bmat -= bmat_offset;
  51     --vlag;
  52     --w;
  53
  54     /* Function Body */
  55     one = 1.;
  56     zero = 0.;
  57     nptm = *npt - *n - 1;
  58     ztest = zero;
  59     i__1 = *npt;
  60     for (k = 1; k <= i__1; ++k) {
  61         i__2 = nptm;
  62         for (j = 1; j <= i__2; ++j) {
  63 /* L10: */
  64 /* Computing MAX */
  65             d__2 = ztest, d__3 = (d__1 = zmat[k + j * zmat_dim1], fabs(d__1));
  66             ztest = MAX2(d__2,d__3);
  67         }
  68     }
  69     ztest *= 1e-20;
  70
  71 /*     Apply the rotations that put zeros in the KNEW-th row of ZMAT. */
  72
  73     jl = 1;
  74     i__2 = nptm;
  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;
  84             i__1 = *npt;
  85             for (i__ = 1; i__ <= i__1; ++i__) {
  86                 temp = tempa * zmat[i__ + zmat_dim1] + tempb * zmat[i__ + j *
  87                         zmat_dim1];
  88                 zmat[i__ + j * zmat_dim1] = tempa * zmat[i__ + j * zmat_dim1]
  89                         - tempb * zmat[i__ + zmat_dim1];
  90 /* L20: */
  91                 zmat[i__ + zmat_dim1] = temp;
  92             }
  93         }
  94         zmat[*knew + j * zmat_dim1] = zero;
  95 /* L30: */
  96     }
  97
  98 /*     Put the first NPT components of the KNEW-th column of HLAG into W, */
  99 /*     and calculate the parameters of the updating formula. */
 100
 101     i__2 = *npt;
 102     for (i__ = 1; i__ <= i__2; ++i__) {
 103         w[i__] = zmat[*knew + zmat_dim1] * zmat[i__ + zmat_dim1];
 104 /* L40: */
 105     }
 106     alpha = w[*knew];
 107     tau = vlag[*knew];
 108     vlag[*knew] -= one;
 109
 110 /*     Complete the updating of ZMAT. */
 111
 112     temp = sqrt(*denom);
 113     tempb = zmat[*knew + zmat_dim1] / temp;
 114     tempa = tau / temp;
 115     i__2 = *npt;
 116     for (i__ = 1; i__ <= i__2; ++i__) {
 117 /* L50: */
 118         zmat[i__ + zmat_dim1] = tempa * zmat[i__ + zmat_dim1] - tempb * vlag[
 119                 i__];
 120     }
 121
 122 /*     Finally, update the matrix BMAT. */
 123
 124     i__2 = *n;
 125     for (j = 1; j <= i__2; ++j) {
 126         jp = *npt + 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;
 130         i__1 = jp;
 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__];
 134             if (i__ > *npt) {
 135                 bmat[jp + (i__ - *npt) * bmat_dim1] = bmat[i__ + j *
 136                         bmat_dim1];
 137             }
 138 /* L60: */
 139         }
 140     }
 141 } /* update_ */
 142
 143 static nlopt_result rescue_(int *n, int *npt, const double *xl, const double *xu,
 144                     /* int *maxfun */
 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)
 152 {
 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;
 157
 158     /* Local variables */
 159     double f;
 160     int i__, j, k, ih, jp, ip, iq, np, iw;
 161     double xp = 0.0, xq, den;
 162     int ihp = 0;
 163     double one;
 164     int ihq, jpn, kpt;
 165     double sum, diff, half, beta;
 166     int kold;
 167     double winc;
 168     int nrem, knew;
 169     double temp, bsum;
 170     int nptm;
 171     double zero, hdiag, fbase, sfrac, denom, vquad, sumpq;
 172     double dsqmin, distsq, vlmxsq;
 173
 174
 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. */
 209
 210 /*     Set some constants. */
 211
 212     /* Parameter adjustments */
 213     zmat_dim1 = *npt;
 214     zmat_offset = 1 + zmat_dim1;
 215     zmat -= zmat_offset;
 216     xpt_dim1 = *npt;
 217     xpt_offset = 1 + xpt_dim1;
 218     xpt -= xpt_offset;
 219     --xl;
 220     --xu;
 221     --xbase;
 222     --fval;
 223     --xopt;
 224     --gopt;
 225     --hq;
 226     --pq;
 227     bmat_dim1 = *ndim;
 228     bmat_offset = 1 + bmat_dim1;
 229     bmat -= bmat_offset;
 230     --sl;
 231     --su;
 232     --vlag;
 233     ptsaux -= 3;
 234     --ptsid;
 235     --w;
 236
 237     /* Function Body */
 238     half = .5;
 239     one = 1.;
 240     zero = 0.;
 241     np = *n + 1;
 242     sfrac = half / (double) np;
 243     nptm = *npt - np;
 244
 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. */
 251
 252     sumpq = zero;
 253     winc = zero;
 254     i__1 = *npt;
 255     for (k = 1; k <= i__1; ++k) {
 256         distsq = zero;
 257         i__2 = *n;
 258         for (j = 1; j <= i__2; ++j) {
 259             xpt[k + j * xpt_dim1] -= xopt[j];
 260 /* L10: */
 261 /* Computing 2nd power */
 262             d__1 = xpt[k + j * xpt_dim1];
 263             distsq += d__1 * d__1;
 264         }
 265         sumpq += pq[k];
 266         w[*ndim + k] = distsq;
 267         winc = MAX2(winc,distsq);
 268         i__2 = nptm;
 269         for (j = 1; j <= i__2; ++j) {
 270 /* L20: */
 271             zmat[k + j * zmat_dim1] = zero;
 272         }
 273     }
 274
 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. */
 277
 278     ih = 0;
 279     i__2 = *n;
 280     for (j = 1; j <= i__2; ++j) {
 281         w[j] = half * sumpq * xopt[j];
 282         i__1 = *npt;
 283         for (k = 1; k <= i__1; ++k) {
 284 /* L30: */
 285             w[j] += pq[k] * xpt[k + j * xpt_dim1];
 286         }
 287         i__1 = j;
 288         for (i__ = 1; i__ <= i__1; ++i__) {
 289             ++ih;
 290 /* L40: */
 291             hq[ih] = hq[ih] + w[i__] * xopt[j] + w[j] * xopt[i__];
 292         }
 293     }
 294
 295 /*     Shift XBASE, SL, SU and XOPT. Set the elements of BMAT to zero, and */
 296 /*     also set the elements of PTSAUX. */
 297
 298     i__1 = *n;
 299     for (j = 1; j <= i__1; ++j) {
 300         xbase[j] += xopt[j];
 301         sl[j] -= xopt[j];
 302         su[j] -= xopt[j];
 303         xopt[j] = zero;
 304 /* Computing MIN */
 305         d__1 = *delta, d__2 = su[j];
 306         ptsaux[(j << 1) + 1] = MIN2(d__1,d__2);
 307 /* Computing MAX */
 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;
 314         }
 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];
 318         }
 319         i__2 = *ndim;
 320         for (i__ = 1; i__ <= i__2; ++i__) {
 321 /* L50: */
 322             bmat[i__ + j * bmat_dim1] = zero;
 323         }
 324     }
 325     fbase = fval[*kopt];
 326
 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. */
 330
 331     ptsid[1] = sfrac;
 332     i__2 = *n;
 333     for (j = 1; j <= i__2; ++j) {
 334         jp = j + 1;
 335         jpn = jp + *n;
 336         ptsid[jp] = (double) j + sfrac;
 337         if (jpn <= *npt) {
 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 +
 343                     j * bmat_dim1];
 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 <<
 347                     1) + 2] * temp;
 348             zmat[jpn + j * zmat_dim1] = -zmat[j * zmat_dim1 + 1] * ptsaux[(j
 349                     << 1) + 1] * temp;
 350         } else {
 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);
 356         }
 357 /* L60: */
 358     }
 359
 360 /*     Set any remaining identifiers with their nonzero elements of ZMAT. */
 361
 362     if (*npt >= *n + np) {
 363         i__2 = *npt;
 364         for (k = np << 1; k <= i__2; ++k) {
 365             iw = (int) (((double) (k - np) - half) / (double) (*n)
 366                     );
 367             ip = k - np - iw * *n;
 368             iq = ip + iw;
 369             if (iq > *n) {
 370                 iq -= *n;
 371             }
 372             ptsid[k] = (double) ip + (double) iq / (double) np +
 373                     sfrac;
 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;
 378 /* L70: */
 379             zmat[k + (k - np) * zmat_dim1] = temp;
 380         }
 381     }
 382     nrem = *npt;
 383     kold = 1;
 384     knew = *kopt;
 385
 386 /*     Reorder the provisional points in the way that exchanges PTSID(KOLD) */
 387 /*     with PTSID(KNEW). */
 388
 389 L80:
 390     i__2 = *n;
 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];
 394 /* L90: */
 395         bmat[knew + j * bmat_dim1] = temp;
 396     }
 397     i__2 = nptm;
 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];
 401 /* L100: */
 402         zmat[knew + j * zmat_dim1] = temp;
 403     }
 404     ptsid[kold] = ptsid[knew];
 405     ptsid[knew] = zero;
 406     w[*ndim + knew] = zero;
 407     --nrem;
 408     if (knew != *kopt) {
 409         temp = vlag[kold];
 410         vlag[kold] = vlag[knew];
 411         vlag[knew] = temp;
 412
 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. */
 417
 418         update_(n, npt, &bmat[bmat_offset], &zmat[zmat_offset], ndim, &vlag[1]
 419                 , &beta, &denom, &knew, &w[1]);
 420         if (nrem == 0) {
 421             goto L350;
 422         }
 423         i__2 = *npt;
 424         for (k = 1; k <= i__2; ++k) {
 425 /* L110: */
 426             w[*ndim + k] = (d__1 = w[*ndim + k], fabs(d__1));
 427         }
 428     }
 429
 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. */
 433
 434 L120:
 435     dsqmin = zero;
 436     i__2 = *npt;
 437     for (k = 1; k <= i__2; ++k) {
 438         if (w[*ndim + k] > zero) {
 439             if (dsqmin == zero || w[*ndim + k] < dsqmin) {
 440                 knew = k;
 441                 dsqmin = w[*ndim + k];
 442             }
 443         }
 444 /* L130: */
 445     }
 446     if (dsqmin == zero) {
 447         goto L260;
 448     }
 449
 450 /*     Form the W-vector of the chosen original interpolation point. */
 451
 452     i__2 = *n;
 453     for (j = 1; j <= i__2; ++j) {
 454 /* L140: */
 455         w[*npt + j] = xpt[knew + j * xpt_dim1];
 456     }
 457     i__2 = *npt;
 458     for (k = 1; k <= i__2; ++k) {
 459         sum = zero;
 460         if (k == *kopt) {
 461         } else if (ptsid[k] == zero) {
 462             i__1 = *n;
 463             for (j = 1; j <= i__1; ++j) {
 464 /* L150: */
 465                 sum += w[*npt + j] * xpt[k + j * xpt_dim1];
 466             }
 467         } else {
 468             ip = (int) ptsid[k];
 469             if (ip > 0) {
 470                 sum = w[*npt + ip] * ptsaux[(ip << 1) + 1];
 471             }
 472             iq = (int) ((double) np * ptsid[k] - (double) (ip *
 473                     np));
 474             if (iq > 0) {
 475                 iw = 1;
 476                 if (ip == 0) {
 477                     iw = 2;
 478                 }
 479                 sum += w[*npt + iq] * ptsaux[iw + (iq << 1)];
 480             }
 481         }
 482 /* L160: */
 483         w[k] = half * sum * sum;
 484     }
 485
 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. */
 488
 489     i__2 = *npt;
 490     for (k = 1; k <= i__2; ++k) {
 491         sum = zero;
 492         i__1 = *n;
 493         for (j = 1; j <= i__1; ++j) {
 494 /* L170: */
 495             sum += bmat[k + j * bmat_dim1] * w[*npt + j];
 496         }
 497 /* L180: */
 498         vlag[k] = sum;
 499     }
 500     beta = zero;
 501     i__2 = nptm;
 502     for (j = 1; j <= i__2; ++j) {
 503         sum = zero;
 504         i__1 = *npt;
 505         for (k = 1; k <= i__1; ++k) {
 506 /* L190: */
 507             sum += zmat[k + j * zmat_dim1] * w[k];
 508         }
 509         beta -= sum * sum;
 510         i__1 = *npt;
 511         for (k = 1; k <= i__1; ++k) {
 512 /* L200: */
 513             vlag[k] += sum * zmat[k + j * zmat_dim1];
 514         }
 515     }
 516     bsum = zero;
 517     distsq = zero;
 518     i__1 = *n;
 519     for (j = 1; j <= i__1; ++j) {
 520         sum = zero;
 521         i__2 = *npt;
 522         for (k = 1; k <= i__2; ++k) {
 523 /* L210: */
 524             sum += bmat[k + j * bmat_dim1] * w[k];
 525         }
 526         jp = j + *npt;
 527         bsum += sum * w[jp];
 528         i__2 = *ndim;
 529         for (ip = *npt + 1; ip <= i__2; ++ip) {
 530 /* L220: */
 531             sum += bmat[ip + j * bmat_dim1] * w[ip];
 532         }
 533         bsum += sum * w[jp];
 534         vlag[jp] = sum;
 535 /* L230: */
 536 /* Computing 2nd power */
 537         d__1 = xpt[knew + j * xpt_dim1];
 538         distsq += d__1 * d__1;
 539     }
 540     beta = half * distsq * distsq + beta - bsum;
 541     vlag[*kopt] += one;
 542
 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. */
 547
 548     denom = zero;
 549     vlmxsq = zero;
 550     i__1 = *npt;
 551     for (k = 1; k <= i__1; ++k) {
 552         if (ptsid[k] != zero) {
 553             hdiag = zero;
 554             i__2 = nptm;
 555             for (j = 1; j <= i__2; ++j) {
 556 /* L240: */
 557 /* Computing 2nd power */
 558                 d__1 = zmat[k + j * zmat_dim1];
 559                 hdiag += d__1 * d__1;
 560             }
 561 /* Computing 2nd power */
 562             d__1 = vlag[k];
 563             den = beta * hdiag + d__1 * d__1;
 564             if (den > denom) {
 565                 kold = k;
 566                 denom = den;
 567             }
 568         }
 569 /* L250: */
 570 /* Computing MAX */
 571 /* Computing 2nd power */
 572         d__3 = vlag[k];
 573         d__1 = vlmxsq, d__2 = d__3 * d__3;
 574         vlmxsq = MAX2(d__1,d__2);
 575     }
 576     if (denom <= vlmxsq * .01) {
 577         w[*ndim + knew] = -w[*ndim + knew] - winc;
 578         goto L120;
 579     }
 580     goto L80;
 581
 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. */
 589
 590 L260:
 591     i__1 = *npt;
 592     for (kpt = 1; kpt <= i__1; ++kpt) {
 593         if (ptsid[kpt] == zero) {
 594             goto L340;
 595         }
 596
 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;
 600
 601         ih = 0;
 602         i__2 = *n;
 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];
 607             i__3 = j;
 608             for (i__ = 1; i__ <= i__3; ++i__) {
 609                 ++ih;
 610 /* L270: */
 611                 hq[ih] += temp * w[i__];
 612             }
 613         }
 614         pq[kpt] = zero;
 615         ip = (int) ptsid[kpt];
 616         iq = (int) ((double) np * ptsid[kpt] - (double) (ip * np))
 617                 ;
 618         if (ip > 0) {
 619             xp = ptsaux[(ip << 1) + 1];
 620             xpt[kpt + ip * xpt_dim1] = xp;
 621         }
 622         if (iq > 0) {
 623             xq = ptsaux[(iq << 1) + 1];
 624             if (ip == 0) {
 625                 xq = ptsaux[(iq << 1) + 2];
 626             }
 627             xpt[kpt + iq * xpt_dim1] = xq;
 628         }
 629
 630 /*     Set VQUAD to the value of the current model at the new point. */
 631
 632         vquad = fbase;
 633         if (ip > 0) {
 634             ihp = (ip + ip * ip) / 2;
 635             vquad += xp * (gopt[ip] + half * xp * hq[ihp]);
 636         }
 637         if (iq > 0) {
 638             ihq = (iq + iq * iq) / 2;
 639             vquad += xq * (gopt[iq] + half * xq * hq[ihq]);
 640             if (ip > 0) {
 641                 iw = MAX2(ihp,ihq) - (i__3 = ip - iq, IABS(i__3));
 642                 vquad += xp * xq * hq[iw];
 643             }
 644         }
 645         i__3 = *npt;
 646         for (k = 1; k <= i__3; ++k) {
 647             temp = zero;
 648             if (ip > 0) {
 649                 temp += xp * xpt[k + ip * xpt_dim1];
 650             }
 651             if (iq > 0) {
 652                 temp += xq * xpt[k + iq * xpt_dim1];
 653             }
 654 /* L280: */
 655             vquad += half * pq[k] * temp * temp;
 656         }
 657
 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. */
 661
 662         i__3 = *n;
 663         for (i__ = 1; i__ <= i__3; ++i__) {
 664 /* Computing MIN */
 665 /* Computing MAX */
 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__]) {
 670                 w[i__] = xl[i__];
 671             }
 672             if (xpt[kpt + i__ * xpt_dim1] == su[i__]) {
 673                 w[i__] = xu[i__];
 674             }
 675 /* L290: */
 676         }
 677
 678         ++ *(stop->nevals_p);
 679         f = calfun(*n, &w[1], calfun_data);
 680         fval[kpt] = f;
 681         if (f < fval[*kopt]) {
 682             *kopt = kpt;
 683         }
 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;
 688
 689         diff = f - vquad;
 690
 691 /*     Update the quadratic model. The RETURN from the subroutine occurs when */
 692 /*     all the new interpolation points are included in the model. */
 693
 694         i__3 = *n;
 695         for (i__ = 1; i__ <= i__3; ++i__) {
 696 /* L310: */
 697             gopt[i__] += diff * bmat[kpt + i__ * bmat_dim1];
 698         }
 699         i__3 = *npt;
 700         for (k = 1; k <= i__3; ++k) {
 701             sum = zero;
 702             i__2 = nptm;
 703             for (j = 1; j <= i__2; ++j) {
 704 /* L320: */
 705                 sum += zmat[k + j * zmat_dim1] * zmat[kpt + j * zmat_dim1];
 706             }
 707             temp = diff * sum;
 708             if (ptsid[k] == zero) {
 709                 pq[k] += temp;
 710             } else {
 711                 ip = (int) ptsid[k];
 712                 iq = (int) ((double) np * ptsid[k] - (double) (ip
 713                         * np));
 714                 ihq = (iq * iq + iq) / 2;
 715                 if (ip == 0) {
 716 /* Computing 2nd power */
 717                     d__1 = ptsaux[(iq << 1) + 2];
 718                     hq[ihq] += temp * (d__1 * d__1);
 719                 } else {
 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);
 724                     if (iq > 0) {
 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 <<
 730                                  1) + 1];
 731                     }
 732                 }
 733             }
 734 /* L330: */
 735         }
 736         ptsid[kpt] = zero;
 737 L340:
 738         ;
 739     }
 740 L350:
 741     return NLOPT_SUCCESS;
 742 } /* rescue_ */
 743
 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,
 749         double *w)
 750 {
 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;
 755
 756     /* Local variables */
 757     int i__, j, k;
 758     double ha, gw, one, diff, half;
 759     int ilbd, isbd;
 760     double slbd;
 761     int iubd;
 762     double vlag, subd, temp;
 763     int ksav = 0;
 764     double step = 0.0, zero, curv;
 765     int iflag;
 766     double scale, csave = 0.0, tempa, tempb, tempd, const__, sumin, ggfree;
 767     int ibdsav = 0;
 768     double dderiv, bigstp, predsq, presav, distsq, stpsav = 0.0, wfixsq, wsqsav;
 769
 770
 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. */
 797
 798 /*     Set the first NPT components of W to the leading elements of the */
 799 /*     KNEW-th column of the H matrix. */
 800
 801     /* Parameter adjustments */
 802     zmat_dim1 = *npt;
 803     zmat_offset = 1 + zmat_dim1;
 804     zmat -= zmat_offset;
 805     xpt_dim1 = *npt;
 806     xpt_offset = 1 + xpt_dim1;
 807     xpt -= xpt_offset;
 808     --xopt;
 809     bmat_dim1 = *ndim;
 810     bmat_offset = 1 + bmat_dim1;
 811     bmat -= bmat_offset;
 812     --sl;
 813     --su;
 814     --xnew;
 815     --xalt;
 816     --glag;
 817     --hcol;
 818     --w;
 819
 820     /* Function Body */
 821     half = .5;
 822     one = 1.;
 823     zero = 0.;
 824     const__ = one + sqrt(2.);
 825     i__1 = *npt;
 826     for (k = 1; k <= i__1; ++k) {
 827 /* L10: */
 828         hcol[k] = zero;
 829     }
 830     i__1 = *npt - *n - 1;
 831     for (j = 1; j <= i__1; ++j) {
 832         temp = zmat[*knew + j * zmat_dim1];
 833         i__2 = *npt;
 834         for (k = 1; k <= i__2; ++k) {
 835 /* L20: */
 836             hcol[k] += temp * zmat[k + j * zmat_dim1];
 837         }
 838     }
 839     *alpha = hcol[*knew];
 840     ha = half * *alpha;
 841
 842 /*     Calculate the gradient of the KNEW-th Lagrange function at XOPT. */
 843
 844     i__2 = *n;
 845     for (i__ = 1; i__ <= i__2; ++i__) {
 846 /* L30: */
 847         glag[i__] = bmat[*knew + i__ * bmat_dim1];
 848     }
 849     i__2 = *npt;
 850     for (k = 1; k <= i__2; ++k) {
 851         temp = zero;
 852         i__1 = *n;
 853         for (j = 1; j <= i__1; ++j) {
 854 /* L40: */
 855             temp += xpt[k + j * xpt_dim1] * xopt[j];
 856         }
 857         temp = hcol[k] * temp;
 858         i__1 = *n;
 859         for (i__ = 1; i__ <= i__1; ++i__) {
 860 /* L50: */
 861             glag[i__] += temp * xpt[k + i__ * xpt_dim1];
 862         }
 863     }
 864
 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. */
 870
 871     presav = zero;
 872     i__1 = *npt;
 873     for (k = 1; k <= i__1; ++k) {
 874         if (k == *kopt) {
 875             goto L80;
 876         }
 877         dderiv = zero;
 878         distsq = zero;
 879         i__2 = *n;
 880         for (i__ = 1; i__ <= i__2; ++i__) {
 881             temp = xpt[k + i__ * xpt_dim1] - xopt[i__];
 882             dderiv += glag[i__] * temp;
 883 /* L60: */
 884             distsq += temp * temp;
 885         }
 886         subd = *adelt / sqrt(distsq);
 887         slbd = -subd;
 888         ilbd = 0;
 889         iubd = 0;
 890         sumin = MIN2(one,subd);
 891
 892 /*     Revise SLBD and SUBD if necessary because of the bounds in SL and SU. */
 893
 894         i__2 = *n;
 895         for (i__ = 1; i__ <= i__2; ++i__) {
 896             temp = xpt[k + i__ * xpt_dim1] - xopt[i__];
 897             if (temp > zero) {
 898                 if (slbd * temp < sl[i__] - xopt[i__]) {
 899                     slbd = (sl[i__] - xopt[i__]) / temp;
 900                     ilbd = -i__;
 901                 }
 902                 if (subd * temp > su[i__] - xopt[i__]) {
 903 /* Computing MAX */
 904                     d__1 = sumin, d__2 = (su[i__] - xopt[i__]) / temp;
 905                     subd = MAX2(d__1,d__2);
 906                     iubd = i__;
 907                 }
 908             } else if (temp < zero) {
 909                 if (slbd * temp > su[i__] - xopt[i__]) {
 910                     slbd = (su[i__] - xopt[i__]) / temp;
 911                     ilbd = i__;
 912                 }
 913                 if (subd * temp < sl[i__] - xopt[i__]) {
 914 /* Computing MAX */
 915                     d__1 = sumin, d__2 = (sl[i__] - xopt[i__]) / temp;
 916                     subd = MAX2(d__1,d__2);
 917                     iubd = -i__;
 918                 }
 919             }
 920 /* L70: */
 921         }
 922
 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. */
 925
 926         if (k == *knew) {
 927             diff = dderiv - one;
 928             step = slbd;
 929             vlag = slbd * (dderiv - slbd * diff);
 930             isbd = ilbd;
 931             temp = subd * (dderiv - subd * diff);
 932             if (fabs(temp) > fabs(vlag)) {
 933                 step = subd;
 934                 vlag = temp;
 935                 isbd = iubd;
 936             }
 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)) {
 943                     step = tempd / diff;
 944                     vlag = temp;
 945                     isbd = 0;
 946                 }
 947             }
 948
 949 /*     Search along each of the other lines through XOPT and another point. */
 950
 951         } else {
 952             step = slbd;
 953             vlag = slbd * (one - slbd);
 954             isbd = ilbd;
 955             temp = subd * (one - subd);
 956             if (fabs(temp) > fabs(vlag)) {
 957                 step = subd;
 958                 vlag = temp;
 959                 isbd = iubd;
 960             }
 961             if (subd > half) {
 962                 if (fabs(vlag) < .25) {
 963                     step = half;
 964                     vlag = .25;
 965                     isbd = 0;
 966                 }
 967             }
 968             vlag *= dderiv;
 969         }
 970
 971 /*     Calculate PREDSQ for the current line search and maintain PRESAV. */
 972
 973         temp = step * (one - step) * distsq;
 974         predsq = vlag * vlag * (vlag * vlag + ha * temp * temp);
 975         if (predsq > presav) {
 976             presav = predsq;
 977             ksav = k;
 978             stpsav = step;
 979             ibdsav = isbd;
 980         }
 981 L80:
 982         ;
 983     }
 984
 985 /*     Construct XNEW in a way that satisfies the bound constraints exactly. */
 986
 987     i__1 = *n;
 988     for (i__ = 1; i__ <= i__1; ++i__) {
 989         temp = xopt[i__] + stpsav * (xpt[ksav + i__ * xpt_dim1] - xopt[i__]);
 990 /* L90: */
 991 /* Computing MAX */
 992 /* Computing MIN */
 993         d__3 = su[i__];
 994         d__1 = sl[i__], d__2 = MIN2(d__3,temp);
 995         xnew[i__] = MAX2(d__1,d__2);
 996     }
 997     if (ibdsav < 0) {
 998         xnew[-ibdsav] = sl[-ibdsav];
 999     }
1000     if (ibdsav > 0) {
1001         xnew[ibdsav] = su[ibdsav];
1002     }
1003
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. */
1007
1008     bigstp = *adelt + *adelt;
1009     iflag = 0;
1010 L100:
1011     wfixsq = zero;
1012     ggfree = zero;
1013     i__1 = *n;
1014     for (i__ = 1; i__ <= i__1; ++i__) {
1015         w[i__] = zero;
1016 /* Computing MIN */
1017         d__1 = xopt[i__] - sl[i__], d__2 = glag[i__];
1018         tempa = MIN2(d__1,d__2);
1019 /* Computing MAX */
1020         d__1 = xopt[i__] - su[i__], d__2 = glag[i__];
1021         tempb = MAX2(d__1,d__2);
1022         if (tempa > zero || tempb < zero) {
1023             w[i__] = bigstp;
1024 /* Computing 2nd power */
1025             d__1 = glag[i__];
1026             ggfree += d__1 * d__1;
1027         }
1028 /* L110: */
1029     }
1030     if (ggfree == zero) {
1031         *cauchy = zero;
1032         goto L200;
1033     }
1034
1035 /*     Investigate whether more components of W can be fixed. */
1036
1037 L120:
1038     temp = *adelt * *adelt - wfixsq;
1039     if (temp > zero) {
1040         wsqsav = wfixsq;
1041         step = sqrt(temp / ggfree);
1042         ggfree = zero;
1043         i__1 = *n;
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 */
1050                     d__1 = w[i__];
1051                     wfixsq += d__1 * d__1;
1052                 } else if (temp >= su[i__]) {
1053                     w[i__] = su[i__] - xopt[i__];
1054 /* Computing 2nd power */
1055                     d__1 = w[i__];
1056                     wfixsq += d__1 * d__1;
1057                 } else {
1058 /* Computing 2nd power */
1059                     d__1 = glag[i__];
1060                     ggfree += d__1 * d__1;
1061                 }
1062             }
1063 /* L130: */
1064         }
1065         if (wfixsq > wsqsav && ggfree > zero) {
1066             goto L120;
1067         }
1068     }
1069
1070 /*     Set the remaining free components of W and all components of XALT, */
1071 /*     except that W may be scaled later. */
1072
1073     gw = zero;
1074     i__1 = *n;
1075     for (i__ = 1; i__ <= i__1; ++i__) {
1076         if (w[i__] == bigstp) {
1077             w[i__] = -step * glag[i__];
1078 /* Computing MAX */
1079 /* Computing MIN */
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__];
1087         } else {
1088             xalt[i__] = su[i__];
1089         }
1090 /* L140: */
1091         gw += glag[i__] * w[i__];
1092     }
1093
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. */
1098
1099     curv = zero;
1100     i__1 = *npt;
1101     for (k = 1; k <= i__1; ++k) {
1102         temp = zero;
1103         i__2 = *n;
1104         for (j = 1; j <= i__2; ++j) {
1105 /* L150: */
1106             temp += xpt[k + j * xpt_dim1] * w[j];
1107         }
1108 /* L160: */
1109         curv += hcol[k] * temp * temp;
1110     }
1111     if (iflag == 1) {
1112         curv = -curv;
1113     }
1114     if (curv > -gw && curv < -const__ * gw) {
1115         scale = -gw / curv;
1116         i__1 = *n;
1117         for (i__ = 1; i__ <= i__1; ++i__) {
1118             temp = xopt[i__] + scale * w[i__];
1119 /* L170: */
1120 /* Computing MAX */
1121 /* Computing MIN */
1122             d__3 = su[i__];
1123             d__1 = sl[i__], d__2 = MIN2(d__3,temp);
1124             xalt[i__] = MAX2(d__1,d__2);
1125         }
1126 /* Computing 2nd power */
1127         d__1 = half * gw * scale;
1128         *cauchy = d__1 * d__1;
1129     } else {
1130 /* Computing 2nd power */
1131         d__1 = gw + half * curv;
1132         *cauchy = d__1 * d__1;
1133     }
1134
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. */
1138
1139     if (iflag == 0) {
1140         i__1 = *n;
1141         for (i__ = 1; i__ <= i__1; ++i__) {
1142             glag[i__] = -glag[i__];
1143 /* L180: */
1144             w[*n + i__] = xalt[i__];
1145         }
1146         csave = *cauchy;
1147         iflag = 1;
1148         goto L100;
1149     }
1150     if (csave > *cauchy) {
1151         i__1 = *n;
1152         for (i__ = 1; i__ <= i__1; ++i__) {
1153 /* L190: */
1154             xalt[i__] = w[*n + i__];
1155         }
1156         *cauchy = csave;
1157     }
1158 L200:
1159     return;
1160 } /* altmov_ */
1161
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)
1167 {
1168     /* System generated locals */
1169     int xpt_dim1, xpt_offset, i__1, i__2;
1170     double d__1, d__2, d__3, d__4;
1171
1172     /* Local variables */
1173     int i__, j, k, ih;
1174     double ds;
1175     int iu;
1176     double dhd, dhs, cth, one, shs, sth, ssq, half, beta, sdec, blen;
1177     int iact = 0, nact;
1178     double angt, qred;
1179     int isav;
1180     double temp, zero, xsav = 0.0, xsum, angbd = 0.0, dredg = 0.0, sredg = 0.0;
1181     int iterc;
1182     double resid, delsq, ggsav = 0.0, tempa, tempb, ratio, sqstp, redmax,
1183             dredsq = 0.0, redsav, onemin, gredsq = 0.0, rednew;
1184     int itcsav = 0;
1185     double rdprev, rdnext = 0.0, stplen, stepsq;
1186     int itermax = 0;
1187
1188
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 */
1216 /*       constrained. */
1217
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. */
1226
1227 /*     Set some constants. */
1228
1229     /* Parameter adjustments */
1230     xpt_dim1 = *npt;
1231     xpt_offset = 1 + xpt_dim1;
1232     xpt -= xpt_offset;
1233     --xopt;
1234     --gopt;
1235     --hq;
1236     --pq;
1237     --sl;
1238     --su;
1239     --xnew;
1240     --d__;
1241     --gnew;
1242     --xbdi;
1243     --s;
1244     --hs;
1245     --hred;
1246
1247     /* Function Body */
1248     half = .5;
1249     one = 1.;
1250     onemin = -1.;
1251     zero = 0.;
1252
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. */
1259
1260     iterc = 0;
1261     nact = 0;
1262     sqstp = zero;
1263     i__1 = *n;
1264     for (i__ = 1; i__ <= i__1; ++i__) {
1265         xbdi[i__] = zero;
1266         if (xopt[i__] <= sl[i__]) {
1267             if (gopt[i__] >= zero) {
1268                 xbdi[i__] = onemin;
1269             }
1270         } else if (xopt[i__] >= su[i__]) {
1271             if (gopt[i__] <= zero) {
1272                 xbdi[i__] = one;
1273             }
1274         }
1275         if (xbdi[i__] != zero) {
1276             ++nact;
1277         }
1278         d__[i__] = zero;
1279 /* L10: */
1280         gnew[i__] = gopt[i__];
1281     }
1282     delsq = *delta * *delta;
1283     qred = zero;
1284     *crvmin = onemin;
1285
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. */
1291
1292 L20:
1293     beta = zero;
1294 L30:
1295     stepsq = zero;
1296     i__1 = *n;
1297     for (i__ = 1; i__ <= i__1; ++i__) {
1298         if (xbdi[i__] != zero) {
1299             s[i__] = zero;
1300         } else if (beta == zero) {
1301             s[i__] = -gnew[i__];
1302         } else {
1303             s[i__] = beta * s[i__] - gnew[i__];
1304         }
1305 /* L40: */
1306 /* Computing 2nd power */
1307         d__1 = s[i__];
1308         stepsq += d__1 * d__1;
1309     }
1310     if (stepsq == zero) {
1311         goto L190;
1312     }
1313     if (beta == zero) {
1314         gredsq = stepsq;
1315         itermax = iterc + *n - nact;
1316     }
1317     if (gredsq * delsq <= qred * 1e-4 * qred) {
1318         goto L190;
1319     }
1320
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. */
1325
1326     goto L210;
1327 L50:
1328     resid = delsq;
1329     ds = zero;
1330     shs = zero;
1331     i__1 = *n;
1332     for (i__ = 1; i__ <= i__1; ++i__) {
1333         if (xbdi[i__] == zero) {
1334 /* Computing 2nd power */
1335             d__1 = d__[i__];
1336             resid -= d__1 * d__1;
1337             ds += s[i__] * d__[i__];
1338             shs += s[i__] * hs[i__];
1339         }
1340 /* L60: */
1341     }
1342     if (resid <= zero) {
1343         goto L90;
1344     }
1345     temp = sqrt(stepsq * resid + ds * ds);
1346     if (ds < zero) {
1347         blen = (temp - ds) / stepsq;
1348     } else {
1349         blen = resid / (temp + ds);
1350     }
1351     stplen = blen;
1352     if (shs > zero) {
1353 /* Computing MIN */
1354         d__1 = blen, d__2 = gredsq / shs;
1355         stplen = MIN2(d__1,d__2);
1356     }
1357
1358 /*     Reduce STPLEN if necessary in order to preserve the simple bounds, */
1359 /*     letting IACT be the index of the new constrained variable. */
1360
1361     iact = 0;
1362     i__1 = *n;
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__];
1368             } else {
1369                 temp = (sl[i__] - xsum) / s[i__];
1370             }
1371             if (temp < stplen) {
1372                 stplen = temp;
1373                 iact = i__;
1374             }
1375         }
1376 /* L70: */
1377     }
1378
1379 /*     Update CRVMIN, GNEW and D. Set SDEC to the decrease that occurs in Q. */
1380
1381     sdec = zero;
1382     if (stplen > zero) {
1383         ++iterc;
1384         temp = shs / stepsq;
1385         if (iact == 0 && temp > zero) {
1386             *crvmin = MIN2(*crvmin,temp);
1387             if (*crvmin == onemin) {
1388                 *crvmin = temp;
1389             }
1390         }
1391         ggsav = gredsq;
1392         gredsq = zero;
1393         i__1 = *n;
1394         for (i__ = 1; i__ <= i__1; ++i__) {
1395             gnew[i__] += stplen * hs[i__];
1396             if (xbdi[i__] == zero) {
1397 /* Computing 2nd power */
1398                 d__1 = gnew[i__];
1399                 gredsq += d__1 * d__1;
1400             }
1401 /* L80: */
1402             d__[i__] += stplen * s[i__];
1403         }
1404 /* Computing MAX */
1405         d__1 = stplen * (ggsav - half * stplen * shs);
1406         sdec = MAX2(d__1,zero);
1407         qred += sdec;
1408     }
1409
1410 /*     Restart the conjugate gradient method if it has hit a new bound. */
1411
1412     if (iact > 0) {
1413         ++nact;
1414         xbdi[iact] = one;
1415         if (s[iact] < zero) {
1416             xbdi[iact] = onemin;
1417         }
1418 /* Computing 2nd power */
1419         d__1 = d__[iact];
1420         delsq -= d__1 * d__1;
1421         if (delsq <= zero) {
1422             goto L90;
1423         }
1424         goto L20;
1425     }
1426
1427 /*     If STPLEN is less than BLEN, then either apply another conjugate */
1428 /*     gradient iteration or RETURN. */
1429
1430     if (stplen < blen) {
1431         if (iterc == itermax) {
1432             goto L190;
1433         }
1434         if (sdec <= qred * .01) {
1435             goto L190;
1436         }
1437         beta = gredsq / ggsav;
1438         goto L30;
1439     }
1440 L90:
1441     *crvmin = zero;
1442
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. */
1446
1447 L100:
1448     if (nact >= *n - 1) {
1449         goto L190;
1450     }
1451     dredsq = zero;
1452     dredg = zero;
1453     gredsq = zero;
1454     i__1 = *n;
1455     for (i__ = 1; i__ <= i__1; ++i__) {
1456         if (xbdi[i__] == zero) {
1457 /* Computing 2nd power */
1458             d__1 = d__[i__];
1459             dredsq += d__1 * d__1;
1460             dredg += d__[i__] * gnew[i__];
1461 /* Computing 2nd power */
1462             d__1 = gnew[i__];
1463             gredsq += d__1 * d__1;
1464             s[i__] = d__[i__];
1465         } else {
1466             s[i__] = zero;
1467         }
1468 /* L110: */
1469     }
1470     itcsav = iterc;
1471     goto L210;
1472
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. */
1475
1476 L120:
1477     ++iterc;
1478     temp = gredsq * dredsq - dredg * dredg;
1479     if (temp <= qred * 1e-4 * qred) {
1480         goto L190;
1481     }
1482     temp = sqrt(temp);
1483     i__1 = *n;
1484     for (i__ = 1; i__ <= i__1; ++i__) {
1485         if (xbdi[i__] == zero) {
1486             s[i__] = (dredg * d__[i__] - dredsq * gnew[i__]) / temp;
1487         } else {
1488             s[i__] = zero;
1489         }
1490 /* L130: */
1491     }
1492     sredg = -temp;
1493
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. */
1498
1499     angbd = one;
1500     iact = 0;
1501     i__1 = *n;
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) {
1507                 ++nact;
1508                 xbdi[i__] = onemin;
1509                 goto L100;
1510             } else if (tempb <= zero) {
1511                 ++nact;
1512                 xbdi[i__] = one;
1513                 goto L100;
1514             }
1515             ratio = one;
1516 /* Computing 2nd power */
1517             d__1 = d__[i__];
1518 /* Computing 2nd power */
1519             d__2 = s[i__];
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;
1524             if (temp > zero) {
1525                 temp = sqrt(temp) - s[i__];
1526                 if (angbd * temp > tempa) {
1527                     angbd = tempa / temp;
1528                     iact = i__;
1529                     xsav = onemin;
1530                 }
1531             }
1532 /* Computing 2nd power */
1533             d__1 = su[i__] - xopt[i__];
1534             temp = ssq - d__1 * d__1;
1535             if (temp > zero) {
1536                 temp = sqrt(temp) + s[i__];
1537                 if (angbd * temp > tempb) {
1538                     angbd = tempb / temp;
1539                     iact = i__;
1540                     xsav = one;
1541                 }
1542             }
1543         }
1544 /* L140: */
1545     }
1546
1547 /*     Calculate HHD and some curvatures for the alternative iteration. */
1548
1549     goto L210;
1550 L150:
1551     shs = zero;
1552     dhs = zero;
1553     dhd = zero;
1554     i__1 = *n;
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__];
1560         }
1561 /* L160: */
1562     }
1563
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. */
1567
1568     redmax = zero;
1569     isav = 0;
1570     redsav = zero;
1571     iu = (int) (angbd * 17. + 3.1);
1572     i__1 = iu;
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) {
1579             redmax = rednew;
1580             isav = i__;
1581             rdprev = redsav;
1582         } else if (i__ == isav + 1) {
1583             rdnext = rednew;
1584         }
1585 /* L170: */
1586         redsav = rednew;
1587     }
1588
1589 /*     Return if the reduction is zero. Otherwise, set the sine and cosine */
1590 /*     of the angle of the alternative iteration, and calculate SDEC. */
1591
1592     if (isav == 0) {
1593         goto L190;
1594     }
1595     if (isav < iu) {
1596         temp = (rdnext - rdprev) / (redmax + redmax - rdprev - rdnext);
1597         angt = angbd * ((double) isav + half * temp) / (double) iu;
1598     }
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);
1603     if (sdec <= zero) {
1604         goto L190;
1605     }
1606
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 */
1609 /*     at the bound. */
1610
1611     dredg = zero;
1612     gredsq = zero;
1613     i__1 = *n;
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 */
1620             d__1 = gnew[i__];
1621             gredsq += d__1 * d__1;
1622         }
1623 /* L180: */
1624         hred[i__] = cth * hred[i__] + sth * hs[i__];
1625     }
1626     qred += sdec;
1627     if (iact > 0 && isav == iu) {
1628         ++nact;
1629         xbdi[iact] = xsav;
1630         goto L100;
1631     }
1632
1633 /*     If SDEC is sufficiently small, then RETURN after setting XNEW to */
1634 /*     XOPT+D, giving careful attention to the bounds. */
1635
1636     if (sdec > qred * .01) {
1637         goto L120;
1638     }
1639 L190:
1640     *dsq = zero;
1641     i__1 = *n;
1642     for (i__ = 1; i__ <= i__1; ++i__) {
1643 /* Computing MAX */
1644 /* Computing MIN */
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__];
1650         }
1651         if (xbdi[i__] == one) {
1652             xnew[i__] = su[i__];
1653         }
1654         d__[i__] = xnew[i__] - xopt[i__];
1655 /* L200: */
1656 /* Computing 2nd power */
1657         d__1 = d__[i__];
1658         *dsq += d__1 * d__1;
1659     }
1660     return;
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. */
1665
1666 L210:
1667     ih = 0;
1668     i__1 = *n;
1669     for (j = 1; j <= i__1; ++j) {
1670         hs[j] = zero;
1671         i__2 = j;
1672         for (i__ = 1; i__ <= i__2; ++i__) {
1673             ++ih;
1674             if (i__ < j) {
1675                 hs[j] += hq[ih] * s[i__];
1676             }
1677 /* L220: */
1678             hs[i__] += hq[ih] * s[j];
1679         }
1680     }
1681     i__2 = *npt;
1682     for (k = 1; k <= i__2; ++k) {
1683         if (pq[k] != zero) {
1684             temp = zero;
1685             i__1 = *n;
1686             for (j = 1; j <= i__1; ++j) {
1687 /* L230: */
1688                 temp += xpt[k + j * xpt_dim1] * s[j];
1689             }
1690             temp *= pq[k];
1691             i__1 = *n;
1692             for (i__ = 1; i__ <= i__1; ++i__) {
1693 /* L240: */
1694                 hs[i__] += temp * xpt[k + i__ * xpt_dim1];
1695             }
1696         }
1697 /* L250: */
1698     }
1699     if (*crvmin != zero) {
1700         goto L50;
1701     }
1702     if (iterc > itcsav) {
1703         goto L150;
1704     }
1705     i__2 = *n;
1706     for (i__ = 1; i__ <= i__2; ++i__) {
1707 /* L260: */
1708         hred[i__] = hs[i__];
1709     }
1710     goto L120;
1711 } /* trsbox_ */
1712
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,
1720                     int *kopt)
1721 {
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;
1726
1727     /* Local variables */
1728     double f;
1729     int i__, j, k, ih, np, nfm;
1730     double one;
1731     int nfx, ipt = 0, jpt = 0;
1732     double two, fbeg, diff, half, temp, zero, recip, stepa = 0.0, stepb = 0.0;
1733     int itemp;
1734     double rhosq;
1735
1736     int nf;
1737
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. */
1749
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. */
1753
1754 /*     Set some constants. */
1755
1756     /* Parameter adjustments */
1757     zmat_dim1 = *npt;
1758     zmat_offset = 1 + zmat_dim1;
1759     zmat -= zmat_offset;
1760     xpt_dim1 = *npt;
1761     xpt_offset = 1 + xpt_dim1;
1762     xpt -= xpt_offset;
1763     --x;
1764     --xl;
1765     --xu;
1766     --xbase;
1767     --fval;
1768     --gopt;
1769     --hq;
1770     --pq;
1771     bmat_dim1 = *ndim;
1772     bmat_offset = 1 + bmat_dim1;
1773     bmat -= bmat_offset;
1774     --sl;
1775     --su;
1776
1777     /* Function Body */
1778     half = .5;
1779     one = 1.;
1780     two = 2.;
1781     zero = 0.;
1782     rhosq = *rhobeg * *rhobeg;
1783     recip = one / rhosq;
1784     np = *n + 1;
1785
1786 /*     Set XBASE to the initial vector of variables, and set the initial */
1787 /*     elements of XPT, BMAT, HQ, PQ and ZMAT to zero. */
1788
1789     i__1 = *n;
1790     for (j = 1; j <= i__1; ++j) {
1791         xbase[j] = x[j];
1792         i__2 = *npt;
1793         for (k = 1; k <= i__2; ++k) {
1794 /* L10: */
1795             xpt[k + j * xpt_dim1] = zero;
1796         }
1797         i__2 = *ndim;
1798         for (i__ = 1; i__ <= i__2; ++i__) {
1799 /* L20: */
1800             bmat[i__ + j * bmat_dim1] = zero;
1801         }
1802     }
1803     i__2 = *n * np / 2;
1804     for (ih = 1; ih <= i__2; ++ih) {
1805 /* L30: */
1806         hq[ih] = zero;
1807     }
1808     i__2 = *npt;
1809     for (k = 1; k <= i__2; ++k) {
1810         pq[k] = zero;
1811         i__1 = *npt - np;
1812         for (j = 1; j <= i__1; ++j) {
1813 /* L40: */
1814             zmat[k + j * zmat_dim1] = zero;
1815         }
1816     }
1817
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,.). */
1821
1822     nf = 0;
1823 L50:
1824     nfm = nf;
1825     nfx = nf - *n;
1826     ++(nf);
1827     if (nfm <= *n << 1) {
1828         if (nfm >= 1 && nfm <= *n) {
1829             stepa = *rhobeg;
1830             if (su[nfm] == zero) {
1831                 stepa = -stepa;
1832             }
1833             xpt[nf + nfm * xpt_dim1] = stepa;
1834         } else if (nfm > *n) {
1835             stepa = xpt[nf - *n + nfx * xpt_dim1];
1836             stepb = -(*rhobeg);
1837             if (sl[nfx] == zero) {
1838 /* Computing MIN */
1839                 d__1 = two * *rhobeg, d__2 = su[nfx];
1840                 stepb = MIN2(d__1,d__2);
1841             }
1842             if (su[nfx] == zero) {
1843 /* Computing MAX */
1844                 d__1 = -two * *rhobeg, d__2 = sl[nfx];
1845                 stepb = MAX2(d__1,d__2);
1846             }
1847             xpt[nf + nfx * xpt_dim1] = stepb;
1848         }
1849     } else {
1850         itemp = (nfm - np) / *n;
1851         jpt = nfm - itemp * *n - *n;
1852         ipt = jpt + itemp;
1853         if (ipt > *n) {
1854             itemp = jpt;
1855             jpt = ipt - *n;
1856             ipt = itemp;
1857         }
1858         xpt[nf + ipt * xpt_dim1] = xpt[ipt + 1 + ipt * xpt_dim1];
1859         xpt[nf + jpt * xpt_dim1] = xpt[jpt + 1 + jpt * xpt_dim1];
1860     }
1861
1862 /*     Calculate the next value of F. The least function value so far and */
1863 /*     its index are required. */
1864
1865     i__1 = *n;
1866     for (j = 1; j <= i__1; ++j) {
1867 /* Computing MIN */
1868 /* Computing MAX */
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]) {
1873             x[j] = xl[j];
1874         }
1875         if (xpt[nf + j * xpt_dim1] == su[j]) {
1876             x[j] = xu[j];
1877         }
1878 /* L60: */
1879     }
1880     ++ *(stop->nevals_p);
1881     f = calfun(*n, &x[1], calfun_data);
1882     fval[nf] = f;
1883     if (nf == 1) {
1884         fbeg = f;
1885         *kopt = 1;
1886     } else if (f < fval[*kopt]) {
1887         *kopt = nf;
1888     }
1889
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. */
1895
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;
1903             }
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];
1913                     fval[nf - *n] = f;
1914                     if (*kopt == nf) {
1915                         *kopt = nf - *n;
1916                     }
1917                     xpt[nf - *n + nfx * xpt_dim1] = stepb;
1918                     xpt[nf + nfx * xpt_dim1] = stepa;
1919                 }
1920             }
1921             bmat[nfx * bmat_dim1 + 1] = -(stepa + stepb) / (stepa * stepb);
1922             bmat[nf + nfx * bmat_dim1] = -half / xpt[nf - *n + nfx *
1923                     xpt_dim1];
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];
1930         }
1931
1932 /*     Set the off-diagonal second derivatives of the Lagrange functions and */
1933 /*     the initial quadratic model. */
1934
1935     } else {
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;
1943     }
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;
1948     if (nf < *npt) {
1949         goto L50;
1950     }
1951     return NLOPT_SUCCESS;
1952 } /* prelim_ */
1953
1954 static nlopt_result bobyqb_(int *n, int *npt, double *x,
1955         const double *xl, const double *xu, double *rhobeg, double *
1956         rhoend,
1957                             nlopt_stopping *stop,
1958                             bobyqa_func calfun, void *calfun_data,
1959                             double *minf,
1960         double *xbase,
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)
1965 {
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;
1970
1971     /* Local variables */
1972     double f;
1973     int i__, j, k, ih, jj, nh, ip, jp;
1974     double dx;
1975     int np;
1976     double den, one, ten, dsq, rho, sum, two, diff, half, beta, gisq;
1977     int knew;
1978     double temp, suma, sumb, bsum, fopt;
1979     int kopt, nptm;
1980     double zero, curv;
1981     int ksav;
1982     double gqsq, dist, sumw, sumz, diffa, diffb, diffc = 0.0, hdiag;
1983     int kbase;
1984     double alpha, delta, adelt, denom, fsave, bdtol, delsq;
1985     int nresc, nfsav;
1986     double ratio, dnorm, vquad, pqold, tenth;
1987     int itest;
1988     double sumpq, scaden;
1989     double errbig, cauchy, fracsq, biglsq, densav;
1990     double bdtest;
1991     double crvmin, frhosq;
1992     double distsq;
1993     int ntrits;
1994     double xoptsq;
1995
1996     nlopt_result rc = NLOPT_SUCCESS, rc2;
1997
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). */
2029
2030 /*     Set some constants. */
2031
2032     /* Parameter adjustments */
2033     zmat_dim1 = *npt;
2034     zmat_offset = 1 + zmat_dim1;
2035     zmat -= zmat_offset;
2036     xpt_dim1 = *npt;
2037     xpt_offset = 1 + xpt_dim1;
2038     xpt -= xpt_offset;
2039     --x;
2040     --xl;
2041     --xu;
2042     --xbase;
2043     --fval;
2044     --xopt;
2045     --gopt;
2046     --hq;
2047     --pq;
2048     bmat_dim1 = *ndim;
2049     bmat_offset = 1 + bmat_dim1;
2050     bmat -= bmat_offset;
2051     --sl;
2052     --su;
2053     --xnew;
2054     --xalt;
2055     --d__;
2056     --vlag;
2057     --w;
2058
2059     /* Function Body */
2060     half = .5;
2061     one = 1.;
2062     ten = 10.;
2063     tenth = .1;
2064     two = 2.;
2065     zero = 0.;
2066     np = *n + 1;
2067     nptm = *npt - np;
2068     nh = *n * np / 2;
2069
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. */
2076
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);
2081     xoptsq = zero;
2082     i__1 = *n;
2083     for (i__ = 1; i__ <= i__1; ++i__) {
2084         xopt[i__] = xpt[kopt + i__ * xpt_dim1];
2085 /* L10: */
2086 /* Computing 2nd power */
2087         d__1 = xopt[i__];
2088         xoptsq += d__1 * d__1;
2089     }
2090     fsave = fval[1];
2091     if (rc2 != NLOPT_SUCCESS) {
2092       rc = rc2;
2093       goto L720;
2094     }
2095     kbase = 1;
2096
2097 /*     Complete the settings that are required for the iterative procedure. */
2098
2099     rho = *rhobeg;
2100     delta = rho;
2101     nresc = *(stop->nevals_p);
2102     ntrits = 0;
2103     diffa = zero;
2104     diffb = zero;
2105     itest = 0;
2106     nfsav = *(stop->nevals_p);
2107
2108 /*     Update GOPT if necessary before the first iteration and after each */
2109 /*     call of RESCUE that makes a call of CALFUN. */
2110
2111 L20:
2112     if (kopt != kbase) {
2113         ih = 0;
2114         i__1 = *n;
2115         for (j = 1; j <= i__1; ++j) {
2116             i__2 = j;
2117             for (i__ = 1; i__ <= i__2; ++i__) {
2118                 ++ih;
2119                 if (i__ < j) {
2120                     gopt[j] += hq[ih] * xopt[i__];
2121                 }
2122 /* L30: */
2123                 gopt[i__] += hq[ih] * xopt[j];
2124             }
2125         }
2126         if (*(stop->nevals_p) > *npt) {
2127             i__2 = *npt;
2128             for (k = 1; k <= i__2; ++k) {
2129                 temp = zero;
2130                 i__1 = *n;
2131                 for (j = 1; j <= i__1; ++j) {
2132 /* L40: */
2133                     temp += xpt[k + j * xpt_dim1] * xopt[j];
2134                 }
2135                 temp = pq[k] * temp;
2136                 i__1 = *n;
2137                 for (i__ = 1; i__ <= i__1; ++i__) {
2138 /* L50: */
2139                     gopt[i__] += temp * xpt[k + i__ * xpt_dim1];
2140                 }
2141             }
2142         }
2143     }
2144
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. */
2151
2152 L60:
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);
2156 /* Computing MIN */
2157     d__1 = delta, d__2 = sqrt(dsq);
2158     dnorm = MIN2(d__1,d__2);
2159     if (dnorm < half * rho) {
2160         ntrits = -1;
2161 /* Computing 2nd power */
2162         d__1 = ten * rho;
2163         distsq = d__1 * d__1;
2164         if (*(stop->nevals_p) <= nfsav + 2) {
2165             goto L650;
2166         }
2167
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. */
2173
2174 /* Computing MAX */
2175         d__1 = MAX2(diffa,diffb);
2176         errbig = MAX2(d__1,diffc);
2177         frhosq = rho * .125 * rho;
2178         if (crvmin > zero && errbig > frhosq * crvmin) {
2179             goto L650;
2180         }
2181         bdtol = errbig / rho;
2182         i__1 = *n;
2183         for (j = 1; j <= i__1; ++j) {
2184             bdtest = bdtol;
2185             if (xnew[j] == sl[j]) {
2186                 bdtest = w[j];
2187             }
2188             if (xnew[j] == su[j]) {
2189                 bdtest = -w[j];
2190             }
2191             if (bdtest < bdtol) {
2192                 curv = hq[(j + j * j) / 2];
2193                 i__2 = *npt;
2194                 for (k = 1; k <= i__2; ++k) {
2195 /* L70: */
2196 /* Computing 2nd power */
2197                     d__1 = xpt[k + j * xpt_dim1];
2198                     curv += pq[k] * (d__1 * d__1);
2199                 }
2200                 bdtest += half * curv * rho;
2201                 if (bdtest < bdtol) {
2202                     goto L650;
2203                 }
2204             }
2205 /* L80: */
2206         }
2207         goto L680;
2208     }
2209     ++ntrits;
2210
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. */
2216
2217 L90:
2218     if (dsq <= xoptsq * .001) {
2219         fracsq = xoptsq * .25;
2220         sumpq = zero;
2221         i__1 = *npt;
2222         for (k = 1; k <= i__1; ++k) {
2223             sumpq += pq[k];
2224             sum = -half * xoptsq;
2225             i__2 = *n;
2226             for (i__ = 1; i__ <= i__2; ++i__) {
2227 /* L100: */
2228                 sum += xpt[k + i__ * xpt_dim1] * xopt[i__];
2229             }
2230             w[*npt + k] = sum;
2231             temp = fracsq - half * sum;
2232             i__2 = *n;
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__];
2236                 ip = *npt + i__;
2237                 i__3 = i__;
2238                 for (j = 1; j <= i__3; ++j) {
2239 /* L110: */
2240                     bmat[ip + j * bmat_dim1] = bmat[ip + j * bmat_dim1] + w[
2241                             i__] * vlag[j] + vlag[i__] * w[j];
2242                 }
2243             }
2244         }
2245
2246 /*     Then the revisions of BMAT that depend on ZMAT are calculated. */
2247
2248         i__3 = nptm;
2249         for (jj = 1; jj <= i__3; ++jj) {
2250             sumz = zero;
2251             sumw = zero;
2252             i__2 = *npt;
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];
2256 /* L120: */
2257                 sumw += vlag[k];
2258             }
2259             i__2 = *n;
2260             for (j = 1; j <= i__2; ++j) {
2261                 sum = (fracsq * sumz - half * sumw) * xopt[j];
2262                 i__1 = *npt;
2263                 for (k = 1; k <= i__1; ++k) {
2264 /* L130: */
2265                     sum += vlag[k] * xpt[k + j * xpt_dim1];
2266                 }
2267                 w[j] = sum;
2268                 i__1 = *npt;
2269                 for (k = 1; k <= i__1; ++k) {
2270 /* L140: */
2271                     bmat[k + j * bmat_dim1] += sum * zmat[k + jj * zmat_dim1];
2272                 }
2273             }
2274             i__1 = *n;
2275             for (i__ = 1; i__ <= i__1; ++i__) {
2276                 ip = i__ + *npt;
2277                 temp = w[i__];
2278                 i__2 = i__;
2279                 for (j = 1; j <= i__2; ++j) {
2280 /* L150: */
2281                     bmat[ip + j * bmat_dim1] += temp * w[j];
2282                 }
2283             }
2284         }
2285
2286 /*     The following instructions complete the shift, including the changes */
2287 /*     to the second derivative parameters of the quadratic model. */
2288
2289         ih = 0;
2290         i__2 = *n;
2291         for (j = 1; j <= i__2; ++j) {
2292             w[j] = -half * sumpq * xopt[j];
2293             i__1 = *npt;
2294             for (k = 1; k <= i__1; ++k) {
2295                 w[j] += pq[k] * xpt[k + j * xpt_dim1];
2296 /* L160: */
2297                 xpt[k + j * xpt_dim1] -= xopt[j];
2298             }
2299             i__1 = j;
2300             for (i__ = 1; i__ <= i__1; ++i__) {
2301                 ++ih;
2302                 hq[ih] = hq[ih] + w[i__] * xopt[j] + xopt[i__] * w[j];
2303 /* L170: */
2304                 bmat[*npt + i__ + j * bmat_dim1] = bmat[*npt + j + i__ *
2305                         bmat_dim1];
2306             }
2307         }
2308         i__1 = *n;
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__];
2314 /* L180: */
2315             xopt[i__] = zero;
2316         }
2317         xoptsq = zero;
2318     }
2319     if (ntrits == 0) {
2320         goto L210;
2321     }
2322     goto L230;
2323
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. */
2332
2333 L190:
2334     nfsav = *(stop->nevals_p);
2335     kbase = kopt;
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]);
2342
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. */
2346
2347     xoptsq = zero;
2348     if (kopt != kbase) {
2349         i__1 = *n;
2350         for (i__ = 1; i__ <= i__1; ++i__) {
2351             xopt[i__] = xpt[kopt + i__ * xpt_dim1];
2352 /* L200: */
2353 /* Computing 2nd power */
2354             d__1 = xopt[i__];
2355             xoptsq += d__1 * d__1;
2356         }
2357     }
2358     if (rc2 != NLOPT_SUCCESS) {
2359       rc = rc2;
2360       goto L720;
2361     }
2362     nresc = *(stop->nevals_p);
2363     if (nfsav < *(stop->nevals_p)) {
2364         nfsav = *(stop->nevals_p);
2365         goto L20;
2366     }
2367     if (ntrits > 0) {
2368         goto L60;
2369     }
2370
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. */
2381
2382 L210:
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]);
2386     i__1 = *n;
2387     for (i__ = 1; i__ <= i__1; ++i__) {
2388 /* L220: */
2389         d__[i__] = xnew[i__] - xopt[i__];
2390     }
2391
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. */
2395
2396 L230:
2397     i__1 = *npt;
2398     for (k = 1; k <= i__1; ++k) {
2399         suma = zero;
2400         sumb = zero;
2401         sum = zero;
2402         i__2 = *n;
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];
2406 /* L240: */
2407             sum += bmat[k + j * bmat_dim1] * d__[j];
2408         }
2409         w[k] = suma * (half * suma + sumb);
2410         vlag[k] = sum;
2411 /* L250: */
2412         w[*npt + k] = suma;
2413     }
2414     beta = zero;
2415     i__1 = nptm;
2416     for (jj = 1; jj <= i__1; ++jj) {
2417         sum = zero;
2418         i__2 = *npt;
2419         for (k = 1; k <= i__2; ++k) {
2420 /* L260: */
2421             sum += zmat[k + jj * zmat_dim1] * w[k];
2422         }
2423         beta -= sum * sum;
2424         i__2 = *npt;
2425         for (k = 1; k <= i__2; ++k) {
2426 /* L270: */
2427             vlag[k] += sum * zmat[k + jj * zmat_dim1];
2428         }
2429     }
2430     dsq = zero;
2431     bsum = zero;
2432     dx = zero;
2433     i__2 = *n;
2434     for (j = 1; j <= i__2; ++j) {
2435 /* Computing 2nd power */
2436         d__1 = d__[j];
2437         dsq += d__1 * d__1;
2438         sum = zero;
2439         i__1 = *npt;
2440         for (k = 1; k <= i__1; ++k) {
2441 /* L280: */
2442             sum += w[k] * bmat[k + j * bmat_dim1];
2443         }
2444         bsum += sum * d__[j];
2445         jp = *npt + j;
2446         i__1 = *n;
2447         for (i__ = 1; i__ <= i__1; ++i__) {
2448 /* L290: */
2449             sum += bmat[jp + i__ * bmat_dim1] * d__[i__];
2450         }
2451         vlag[jp] = sum;
2452         bsum += sum * d__[j];
2453 /* L300: */
2454         dx += d__[j] * xopt[j];
2455     }
2456     beta = dx * dx + dsq * (xoptsq + dx + dx + half * dsq) + beta - bsum;
2457     vlag[kopt] += one;
2458
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. */
2462
2463     if (ntrits == 0) {
2464 /* Computing 2nd power */
2465         d__1 = vlag[knew];
2466         denom = d__1 * d__1 + alpha * beta;
2467         if (denom < cauchy && cauchy > zero) {
2468             i__2 = *n;
2469             for (i__ = 1; i__ <= i__2; ++i__) {
2470                 xnew[i__] = xalt[i__];
2471 /* L310: */
2472                 d__[i__] = xnew[i__] - xopt[i__];
2473             }
2474             cauchy = zero;
2475             goto L230;
2476         }
2477 /* Computing 2nd power */
2478         d__1 = vlag[knew];
2479         if (denom <= half * (d__1 * d__1)) {
2480             if (*(stop->nevals_p) > nresc) {
2481                 goto L190;
2482             }
2483             /* Return from BOBYQA because of much cancellation in a
2484                denominator. */
2485             rc = NLOPT_ROUNDOFF_LIMITED;
2486             goto L720;
2487         }
2488
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. */
2494
2495     } else {
2496         delsq = delta * delta;
2497         scaden = zero;
2498         biglsq = zero;
2499         knew = 0;
2500         i__2 = *npt;
2501         for (k = 1; k <= i__2; ++k) {
2502             if (k == kopt) {
2503                 goto L350;
2504             }
2505             hdiag = zero;
2506             i__1 = nptm;
2507             for (jj = 1; jj <= i__1; ++jj) {
2508 /* L330: */
2509 /* Computing 2nd power */
2510                 d__1 = zmat[k + jj * zmat_dim1];
2511                 hdiag += d__1 * d__1;
2512             }
2513 /* Computing 2nd power */
2514             d__1 = vlag[k];
2515             den = beta * hdiag + d__1 * d__1;
2516             distsq = zero;
2517             i__1 = *n;
2518             for (j = 1; j <= i__1; ++j) {
2519 /* L340: */
2520 /* Computing 2nd power */
2521                 d__1 = xpt[k + j * xpt_dim1] - xopt[j];
2522                 distsq += d__1 * d__1;
2523             }
2524 /* Computing MAX */
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;
2531                 knew = k;
2532                 denom = den;
2533             }
2534 /* Computing MAX */
2535 /* Computing 2nd power */
2536             d__3 = vlag[k];
2537             d__1 = biglsq, d__2 = temp * (d__3 * d__3);
2538             biglsq = MAX2(d__1,d__2);
2539 L350:
2540             ;
2541         }
2542         if (scaden <= half * biglsq) {
2543             if (*(stop->nevals_p) > nresc) {
2544                 goto L190;
2545             }
2546             /* Return from BOBYQA because of much cancellation in a
2547                denominator. */
2548             rc = NLOPT_ROUNDOFF_LIMITED;
2549             goto L720;
2550         }
2551     }
2552
2553 /*     Put the variables for the next calculation of the objective function */
2554 /*       in XNEW, with any adjustments for the bounds. */
2555
2556
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. */
2559
2560 L360:
2561     i__2 = *n;
2562     for (i__ = 1; i__ <= i__2; ++i__) {
2563 /* Computing MIN */
2564 /* Computing MAX */
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__]) {
2569             x[i__] = xl[i__];
2570         }
2571         if (xnew[i__] == su[i__]) {
2572             x[i__] = xu[i__];
2573         }
2574 /* L380: */
2575     }
2576
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;
2581
2582     ++ *(stop->nevals_p);
2583     f = calfun(*n, &x[1], calfun_data);
2584     if (ntrits == -1) {
2585         fsave = f;
2586         rc = NLOPT_XTOL_REACHED;
2587         if (fsave < fval[kopt]) { *minf = f; return rc; }
2588         goto L720;
2589     }
2590
2591     if (f < stop->minf_max) {
2592       *minf = f;
2593       return NLOPT_MINF_MAX_REACHED;
2594     }
2595
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. */
2598
2599     fopt = fval[kopt];
2600     vquad = zero;
2601     ih = 0;
2602     i__2 = *n;
2603     for (j = 1; j <= i__2; ++j) {
2604         vquad += d__[j] * gopt[j];
2605         i__1 = j;
2606         for (i__ = 1; i__ <= i__1; ++i__) {
2607             ++ih;
2608             temp = d__[i__] * d__[j];
2609             if (i__ == j) {
2610                 temp = half * temp;
2611             }
2612 /* L410: */
2613             vquad += hq[ih] * temp;
2614         }
2615     }
2616     i__1 = *npt;
2617     for (k = 1; k <= i__1; ++k) {
2618 /* L420: */
2619 /* Computing 2nd power */
2620         d__1 = w[*npt + k];
2621         vquad += half * pq[k] * (d__1 * d__1);
2622     }
2623     diff = f - fopt - vquad;
2624     diffc = diffb;
2625     diffb = diffa;
2626     diffa = fabs(diff);
2627     if (dnorm > rho) {
2628         nfsav = *(stop->nevals_p);
2629     }
2630
2631 /*     Pick the next value of DELTA after a trust region step. */
2632
2633     if (ntrits > 0) {
2634         if (vquad >= zero) {
2635           /* Return from BOBYQA because a trust region step has failed
2636              to reduce Q. */
2637           rc = NLOPT_ROUNDOFF_LIMITED; /* or FTOL_REACHED? */
2638           goto L720;
2639         }
2640         ratio = (f - fopt) / vquad;
2641         if (ratio <= tenth) {
2642 /* Computing MIN */
2643             d__1 = half * delta;
2644             delta = MIN2(d__1,dnorm);
2645         } else if (ratio <= .7) {
2646 /* Computing MAX */
2647             d__1 = half * delta;
2648             delta = MAX2(d__1,dnorm);
2649         } else {
2650 /* Computing MAX */
2651             d__1 = half * delta, d__2 = dnorm + dnorm;
2652             delta = MAX2(d__1,d__2);
2653         }
2654         if (delta <= rho * 1.5) {
2655             delta = rho;
2656         }
2657
2658 /*     Recalculate KNEW and DENOM if the new F is less than FOPT. */
2659
2660         if (f < fopt) {
2661             ksav = knew;
2662             densav = denom;
2663             delsq = delta * delta;
2664             scaden = zero;
2665             biglsq = zero;
2666             knew = 0;
2667             i__1 = *npt;
2668             for (k = 1; k <= i__1; ++k) {
2669                 hdiag = zero;
2670                 i__2 = nptm;
2671                 for (jj = 1; jj <= i__2; ++jj) {
2672 /* L440: */
2673 /* Computing 2nd power */
2674                     d__1 = zmat[k + jj * zmat_dim1];
2675                     hdiag += d__1 * d__1;
2676                 }
2677 /* Computing 2nd power */
2678                 d__1 = vlag[k];
2679                 den = beta * hdiag + d__1 * d__1;
2680                 distsq = zero;
2681                 i__2 = *n;
2682                 for (j = 1; j <= i__2; ++j) {
2683 /* L450: */
2684 /* Computing 2nd power */
2685                     d__1 = xpt[k + j * xpt_dim1] - xnew[j];
2686                     distsq += d__1 * d__1;
2687                 }
2688 /* Computing MAX */
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;
2695                     knew = k;
2696                     denom = den;
2697                 }
2698 /* L460: */
2699 /* Computing MAX */
2700 /* Computing 2nd power */
2701                 d__3 = vlag[k];
2702                 d__1 = biglsq, d__2 = temp * (d__3 * d__3);
2703                 biglsq = MAX2(d__1,d__2);
2704             }
2705             if (scaden <= half * biglsq) {
2706                 knew = ksav;
2707                 denom = densav;
2708             }
2709         }
2710     }
2711
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. */
2714
2715     update_(n, npt, &bmat[bmat_offset], &zmat[zmat_offset], ndim, &vlag[1], &
2716             beta, &denom, &knew, &w[1]);
2717     ih = 0;
2718     pqold = pq[knew];
2719     pq[knew] = zero;
2720     i__1 = *n;
2721     for (i__ = 1; i__ <= i__1; ++i__) {
2722         temp = pqold * xpt[knew + i__ * xpt_dim1];
2723         i__2 = i__;
2724         for (j = 1; j <= i__2; ++j) {
2725             ++ih;
2726 /* L470: */
2727             hq[ih] += temp * xpt[knew + j * xpt_dim1];
2728         }
2729     }
2730     i__2 = nptm;
2731     for (jj = 1; jj <= i__2; ++jj) {
2732         temp = diff * zmat[knew + jj * zmat_dim1];
2733         i__1 = *npt;
2734         for (k = 1; k <= i__1; ++k) {
2735 /* L480: */
2736             pq[k] += temp * zmat[k + jj * zmat_dim1];
2737         }
2738     }
2739
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. */
2742
2743     fval[knew] = f;
2744     i__1 = *n;
2745     for (i__ = 1; i__ <= i__1; ++i__) {
2746         xpt[knew + i__ * xpt_dim1] = xnew[i__];
2747 /* L490: */
2748         w[i__] = bmat[knew + i__ * bmat_dim1];
2749     }
2750     i__1 = *npt;
2751     for (k = 1; k <= i__1; ++k) {
2752         suma = zero;
2753         i__2 = nptm;
2754         for (jj = 1; jj <= i__2; ++jj) {
2755 /* L500: */
2756             suma += zmat[knew + jj * zmat_dim1] * zmat[k + jj * zmat_dim1];
2757         }
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;
2762           goto L720;
2763         }
2764         sumb = zero;
2765         i__2 = *n;
2766         for (j = 1; j <= i__2; ++j) {
2767 /* L510: */
2768             sumb += xpt[k + j * xpt_dim1] * xopt[j];
2769         }
2770         temp = suma * sumb;
2771         i__2 = *n;
2772         for (i__ = 1; i__ <= i__2; ++i__) {
2773 /* L520: */
2774             w[i__] += temp * xpt[k + i__ * xpt_dim1];
2775         }
2776     }
2777     i__2 = *n;
2778     for (i__ = 1; i__ <= i__2; ++i__) {
2779 /* L530: */
2780         gopt[i__] += diff * w[i__];
2781     }
2782
2783 /*     Update XOPT, GOPT and KOPT if the new calculated F is less than FOPT. */
2784
2785     if (f < fopt) {
2786         kopt = knew;
2787         xoptsq = zero;
2788         ih = 0;
2789         i__2 = *n;
2790         for (j = 1; j <= i__2; ++j) {
2791             xopt[j] = xnew[j];
2792 /* Computing 2nd power */
2793             d__1 = xopt[j];
2794             xoptsq += d__1 * d__1;
2795             i__1 = j;
2796             for (i__ = 1; i__ <= i__1; ++i__) {
2797                 ++ih;
2798                 if (i__ < j) {
2799                     gopt[j] += hq[ih] * d__[i__];
2800                 }
2801 /* L540: */
2802                 gopt[i__] += hq[ih] * d__[j];
2803             }
2804         }
2805         i__1 = *npt;
2806         for (k = 1; k <= i__1; ++k) {
2807             temp = zero;
2808             i__2 = *n;
2809             for (j = 1; j <= i__2; ++j) {
2810 /* L550: */
2811                 temp += xpt[k + j * xpt_dim1] * d__[j];
2812             }
2813             temp = pq[k] * temp;
2814             i__2 = *n;
2815             for (i__ = 1; i__ <= i__2; ++i__) {
2816 /* L560: */
2817                 gopt[i__] += temp * xpt[k + i__ * xpt_dim1];
2818             }
2819         }
2820         if (nlopt_stop_ftol(stop, f, fopt)) {
2821              rc = NLOPT_FTOL_REACHED;
2822              goto L720;
2823         }
2824     }
2825
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. */
2829
2830     if (ntrits > 0) {
2831         i__2 = *npt;
2832         for (k = 1; k <= i__2; ++k) {
2833             vlag[k] = fval[k] - fval[kopt];
2834 /* L570: */
2835             w[k] = zero;
2836         }
2837         i__2 = nptm;
2838         for (j = 1; j <= i__2; ++j) {
2839             sum = zero;
2840             i__1 = *npt;
2841             for (k = 1; k <= i__1; ++k) {
2842 /* L580: */
2843                 sum += zmat[k + j * zmat_dim1] * vlag[k];
2844             }
2845             i__1 = *npt;
2846             for (k = 1; k <= i__1; ++k) {
2847 /* L590: */
2848                 w[k] += sum * zmat[k + j * zmat_dim1];
2849             }
2850         }
2851         i__1 = *npt;
2852         for (k = 1; k <= i__1; ++k) {
2853             sum = zero;
2854             i__2 = *n;
2855             for (j = 1; j <= i__2; ++j) {
2856 /* L600: */
2857                 sum += xpt[k + j * xpt_dim1] * xopt[j];
2858             }
2859             w[k + *npt] = w[k];
2860 /* L610: */
2861             w[k] = sum * w[k];
2862         }
2863         gqsq = zero;
2864         gisq = zero;
2865         i__1 = *n;
2866         for (i__ = 1; i__ <= i__1; ++i__) {
2867             sum = zero;
2868             i__2 = *npt;
2869             for (k = 1; k <= i__2; ++k) {
2870 /* L620: */
2871                 sum = sum + bmat[k + i__ * bmat_dim1] * vlag[k] + xpt[k + i__
2872                         * xpt_dim1] * w[k];
2873             }
2874             if (xopt[i__] == sl[i__]) {
2875 /* Computing MIN */
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__]) {
2884 /* Computing MAX */
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;
2892             } else {
2893 /* Computing 2nd power */
2894                 d__1 = gopt[i__];
2895                 gqsq += d__1 * d__1;
2896                 gisq += sum * sum;
2897             }
2898 /* L630: */
2899             vlag[*npt + i__] = sum;
2900         }
2901
2902 /*     Test whether to replace the new quadratic model by the least Frobenius */
2903 /*     norm interpolant, making the replacement if the test is satisfied. */
2904
2905         ++itest;
2906         if (gqsq < ten * gisq) {
2907             itest = 0;
2908         }
2909         if (itest >= 3) {
2910             i__1 = MAX2(*npt,nh);
2911             for (i__ = 1; i__ <= i__1; ++i__) {
2912                 if (i__ <= *n) {
2913                     gopt[i__] = vlag[*npt + i__];
2914                 }
2915                 if (i__ <= *npt) {
2916                     pq[i__] = w[*npt + i__];
2917                 }
2918                 if (i__ <= nh) {
2919                     hq[i__] = zero;
2920                 }
2921                 itest = 0;
2922 /* L640: */
2923             }
2924         }
2925     }
2926
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. */
2930
2931     if (ntrits == 0) {
2932         goto L60;
2933     }
2934     if (f <= fopt + tenth * vquad) {
2935         goto L60;
2936     }
2937
2938 /*     Alternatively, find out if the interpolation points are close enough */
2939 /*       to the best point so far. */
2940
2941 /* Computing MAX */
2942 /* Computing 2nd power */
2943     d__3 = two * delta;
2944 /* Computing 2nd power */
2945     d__4 = ten * rho;
2946     d__1 = d__3 * d__3, d__2 = d__4 * d__4;
2947     distsq = MAX2(d__1,d__2);
2948 L650:
2949     knew = 0;
2950     i__1 = *npt;
2951     for (k = 1; k <= i__1; ++k) {
2952         sum = zero;
2953         i__2 = *n;
2954         for (j = 1; j <= i__2; ++j) {
2955 /* L660: */
2956 /* Computing 2nd power */
2957             d__1 = xpt[k + j * xpt_dim1] - xopt[j];
2958             sum += d__1 * d__1;
2959         }
2960         if (sum > distsq) {
2961             knew = k;
2962             distsq = sum;
2963         }
2964 /* L670: */
2965     }
2966
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. */
2972
2973     if (knew > 0) {
2974         dist = sqrt(distsq);
2975         if (ntrits == -1) {
2976 /* Computing MIN */
2977             d__1 = tenth * delta, d__2 = half * dist;
2978             delta = MIN2(d__1,d__2);
2979             if (delta <= rho * 1.5) {
2980                 delta = rho;
2981             }
2982         }
2983         ntrits = 0;
2984 /* Computing MAX */
2985 /* Computing MIN */
2986         d__2 = tenth * dist;
2987         d__1 = MIN2(d__2,delta);
2988         adelt = MAX2(d__1,rho);
2989         dsq = adelt * adelt;
2990         goto L90;
2991     }
2992     if (ntrits == -1) {
2993         goto L680;
2994     }
2995     if (ratio > zero) {
2996         goto L60;
2997     }
2998     if (MAX2(delta,dnorm) > rho) {
2999         goto L60;
3000     }
3001
3002 /*     The calculations with the current value of RHO are complete. Pick the */
3003 /*       next values of RHO and DELTA. */
3004
3005 L680:
3006     if (rho > *rhoend) {
3007         delta = half * rho;
3008         ratio = rho / *rhoend;
3009         if (ratio <= 16.) {
3010             rho = *rhoend;
3011         } else if (ratio <= 250.) {
3012             rho = sqrt(ratio) * *rhoend;
3013         } else {
3014             rho = tenth * rho;
3015         }
3016         delta = MAX2(delta,rho);
3017         ntrits = 0;
3018         nfsav = *(stop->nevals_p);
3019         goto L60;
3020     }
3021
3022 /*     Return from the calculation, after another Newton-Raphson step, if */
3023 /*       it is too short to have been tried before. */
3024
3025     if (ntrits == -1) {
3026         goto L360;
3027     }
3028 L720:
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. */
3035     {
3036         i__1 = *n;
3037         for (i__ = 1; i__ <= i__1; ++i__) {
3038 /* Computing MIN */
3039 /* Computing MAX */
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__]) {
3044                 x[i__] = xl[i__];
3045             }
3046             if (xopt[i__] == su[i__]) {
3047                 x[i__] = xu[i__];
3048             }
3049 /* L730: */
3050         }
3051         f = fval[kopt];
3052     }
3053     *minf = f;
3054     return rc;
3055 } /* bobyqb_ */
3056
3057 /**************************************************************************/
3058
3059 #define U(n) ((unsigned) (n))
3060
3061 typedef struct {
3062      double *s, *xs;
3063      nlopt_func f; void *f_data;
3064 } rescale_fun_data;
3065
3066 static double rescale_fun(int n, const double *x, void *d_)
3067 {
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);
3071 }
3072
3073 nlopt_result bobyqa(int n, int npt, double *x,
3074                     const double *xl, const double *xu,
3075                     const double *dx,
3076                     nlopt_stopping *stop, double *minf,
3077                     nlopt_func f, void *f_data)
3078 {
3079     /* System generated locals */
3080     int i__1;
3081     double d__1, d__2;
3082
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;
3086     double temp, zero;
3087     int ibmat, izmat;
3088
3089     double rhobeg, rhoend;
3090     double *w0 = NULL, *w;
3091     nlopt_result ret;
3092     double *s = NULL, *sxl = NULL, *sxu = NULL, *xs = NULL;
3093     rescale_fun_data calfun_data;
3094
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;
3103         }
3104
3105     /* this statement must go before goto done, so that --x occurs */
3106     nlopt_rescale(U(n), s, x, x); --x;
3107
3108     xs = (double *) malloc(sizeof(double) * (U(n)));
3109     if (!xs) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
3110
3111     sxl = nlopt_new_rescaled(U(n), s, xl);
3112     if (!sxl) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
3113     xl = sxl;
3114     sxu = nlopt_new_rescaled(U(n), s, xu);
3115     if (!sxu) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
3116     xu = sxu;
3117     nlopt_reorder_bounds(n, sxl, sxu);
3118
3119     rhobeg = fabs(dx[0] / s[0]); /* equals all other dx[i] after rescaling */
3120
3121     calfun_data.s = s;
3122     calfun_data.xs = xs;
3123     calfun_data.f = f;
3124     calfun_data.f_data = f_data;
3125
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]);
3131
3132
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. */
3140
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 */
3144 /*       recommended. */
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. */
3162
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. */
3167
3168 /*     Return if the value of NPT is unacceptable. */
3169
3170     /* Parameter adjustments */
3171     --xu;
3172     --xl;
3173
3174     /* Function Body */
3175     np = n + 1;
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");
3180       goto done;
3181     }
3182
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. */
3187
3188     ndim = npt + n;
3189     ixb = 1;
3190     ixp = ixb + n;
3191     ifv = ixp + n * npt;
3192     ixo = ifv + npt;
3193     igo = ixo + n;
3194     ihq = igo + n;
3195     ipq = ihq + n * np / 2;
3196     ibmat = ipq + npt;
3197     izmat = ibmat + ndim * n;
3198     isl = izmat + npt * (npt - np);
3199     isu = isl + n;
3200     ixn = isu + n;
3201     ixa = ixn + n;
3202     id = ixa + n;
3203     ivl = id + n;
3204     iw = ivl + ndim;
3205
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; }
3208     w = w0 - 1;
3209
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. */
3216
3217     zero = 0.;
3218     i__1 = n;
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);
3227              goto done;
3228         }
3229         jsl = isl + j - 1;
3230         jsu = jsl + n;
3231         w[jsl] = xl[j] - x[j];
3232         w[jsu] = xu[j] - x[j];
3233         if (w[jsl] >= -(rhobeg)) {
3234             if (w[jsl] >= zero) {
3235                 x[j] = xl[j];
3236                 w[jsl] = zero;
3237                 w[jsu] = temp;
3238             } else {
3239                 x[j] = xl[j] + rhobeg;
3240                 w[jsl] = -(rhobeg);
3241 /* Computing MAX */
3242                 d__1 = xu[j] - x[j];
3243                 w[jsu] = MAX2(d__1,rhobeg);
3244             }
3245         } else if (w[jsu] <= rhobeg) {
3246             if (w[jsu] <= zero) {
3247                 x[j] = xu[j];
3248                 w[jsl] = -temp;
3249                 w[jsu] = zero;
3250             } else {
3251                 x[j] = xu[j] - rhobeg;
3252 /* Computing MIN */
3253                 d__1 = xl[j] - x[j], d__2 = -(rhobeg);
3254                 w[jsl] = MIN2(d__1,d__2);
3255                 w[jsu] = rhobeg;
3256             }
3257         }
3258 /* L30: */
3259     }
3260
3261 /*     Make the call of BOBYQB. */
3262
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]);
3268
3269 done:
3270     free(w0);
3271     free(sxl);
3272     free(sxu);
3273     free(xs);
3274     ++x; nlopt_unscale(U(n), s, x, x);
3275     free(s);
3276     return ret;
3277 } /* bobyqa_ */
3278