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 #define min(a,b) ((a) <= (b) ? (a) : (b))
  13 #define max(a,b) ((a) >= (b) ? (a) : (b))
  14 #define iabs(x) ((x) < 0 ? -(x) : (x))
  15
  16 static void update_(int *n, int *npt, double *bmat,
  17         double *zmat, int *ndim, double *vlag, double *beta,
  18         double *denom, int *knew, double *w)
  19 {
  20     /* System generated locals */
  21     int bmat_dim1, bmat_offset, zmat_dim1, zmat_offset, i__1, i__2;
  22     double d__1, d__2, d__3;
  23
  24     /* Local variables */
  25     int i__, j, k, jl, jp;
  26     double one, tau, temp;
  27     int nptm;
  28     double zero, alpha, tempa, tempb, ztest;
  29
  30
  31 /*     The arrays BMAT and ZMAT are updated, as required by the new position */
  32 /*     of the interpolation point that has the index KNEW. The vector VLAG has */
  33 /*     N+NPT components, set on entry to the first NPT and last N components */
  34 /*     of the product Hw in equation (4.11) of the Powell (2006) paper on */
  35 /*     NEWUOA. Further, BETA is set on entry to the value of the parameter */
  36 /*     with that name, and DENOM is set to the denominator of the updating */
  37 /*     formula. Elements of ZMAT may be treated as zero if their moduli are */
  38 /*     at most ZTEST. The first NDIM elements of W are used for working space. */
  39
  40 /*     Set some constants. */
  41
  42     /* Parameter adjustments */
  43     zmat_dim1 = *npt;
  44     zmat_offset = 1 + zmat_dim1;
  45     zmat -= zmat_offset;
  46     bmat_dim1 = *ndim;
  47     bmat_offset = 1 + bmat_dim1;
  48     bmat -= bmat_offset;
  49     --vlag;
  50     --w;
  51
  52     /* Function Body */
  53     one = 1.;
  54     zero = 0.;
  55     nptm = *npt - *n - 1;
  56     ztest = zero;
  57     i__1 = *npt;
  58     for (k = 1; k <= i__1; ++k) {
  59         i__2 = nptm;
  60         for (j = 1; j <= i__2; ++j) {
  61 /* L10: */
  62 /* Computing MAX */
  63             d__2 = ztest, d__3 = (d__1 = zmat[k + j * zmat_dim1], fabs(d__1));
  64             ztest = max(d__2,d__3);
  65         }
  66     }
  67     ztest *= 1e-20;
  68
  69 /*     Apply the rotations that put zeros in the KNEW-th row of ZMAT. */
  70
  71     jl = 1;
  72     i__2 = nptm;
  73     for (j = 2; j <= i__2; ++j) {
  74         if ((d__1 = zmat[*knew + j * zmat_dim1], fabs(d__1)) > ztest) {
  75 /* Computing 2nd power */
  76             d__1 = zmat[*knew + zmat_dim1];
  77 /* Computing 2nd power */
  78             d__2 = zmat[*knew + j * zmat_dim1];
  79             temp = sqrt(d__1 * d__1 + d__2 * d__2);
  80             tempa = zmat[*knew + zmat_dim1] / temp;
  81             tempb = zmat[*knew + j * zmat_dim1] / temp;
  82             i__1 = *npt;
  83             for (i__ = 1; i__ <= i__1; ++i__) {
  84                 temp = tempa * zmat[i__ + zmat_dim1] + tempb * zmat[i__ + j *
  85                         zmat_dim1];
  86                 zmat[i__ + j * zmat_dim1] = tempa * zmat[i__ + j * zmat_dim1]
  87                         - tempb * zmat[i__ + zmat_dim1];
  88 /* L20: */
  89                 zmat[i__ + zmat_dim1] = temp;
  90             }
  91         }
  92         zmat[*knew + j * zmat_dim1] = zero;
  93 /* L30: */
  94     }
  95
  96 /*     Put the first NPT components of the KNEW-th column of HLAG into W, */
  97 /*     and calculate the parameters of the updating formula. */
  98
  99     i__2 = *npt;
 100     for (i__ = 1; i__ <= i__2; ++i__) {
 101         w[i__] = zmat[*knew + zmat_dim1] * zmat[i__ + zmat_dim1];
 102 /* L40: */
 103     }
 104     alpha = w[*knew];
 105     tau = vlag[*knew];
 106     vlag[*knew] -= one;
 107
 108 /*     Complete the updating of ZMAT. */
 109
 110     temp = sqrt(*denom);
 111     tempb = zmat[*knew + zmat_dim1] / temp;
 112     tempa = tau / temp;
 113     i__2 = *npt;
 114     for (i__ = 1; i__ <= i__2; ++i__) {
 115 /* L50: */
 116         zmat[i__ + zmat_dim1] = tempa * zmat[i__ + zmat_dim1] - tempb * vlag[
 117                 i__];
 118     }
 119
 120 /*     Finally, update the matrix BMAT. */
 121
 122     i__2 = *n;
 123     for (j = 1; j <= i__2; ++j) {
 124         jp = *npt + j;
 125         w[jp] = bmat[*knew + j * bmat_dim1];
 126         tempa = (alpha * vlag[jp] - tau * w[jp]) / *denom;
 127         tempb = (-(*beta) * w[jp] - tau * vlag[jp]) / *denom;
 128         i__1 = jp;
 129         for (i__ = 1; i__ <= i__1; ++i__) {
 130             bmat[i__ + j * bmat_dim1] = bmat[i__ + j * bmat_dim1] + tempa *
 131                     vlag[i__] + tempb * w[i__];
 132             if (i__ > *npt) {
 133                 bmat[jp + (i__ - *npt) * bmat_dim1] = bmat[i__ + j *
 134                         bmat_dim1];
 135             }
 136 /* L60: */
 137         }
 138     }
 139 } /* update_ */
 140
 141 static nlopt_result rescue_(int *n, int *npt, const double *xl, const double *xu,
 142                     /* int *maxfun */
 143                     nlopt_stopping *stop,
 144                     bobyqa_func calfun, void *calfun_data,
 145         double *xbase, double *xpt, double *fval, double *xopt, double *gopt,
 146          double *hq, double *pq, double *bmat, double *zmat,
 147         int *ndim, double *sl, double *su, /* int *nf,  */
 148         double *delta, int *kopt, double *vlag, double *
 149         ptsaux, double *ptsid, double *w)
 150 {
 151     /* System generated locals */
 152     int xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1,
 153             zmat_offset, i__1, i__2, i__3;
 154     double d__1, d__2, d__3, d__4;
 155
 156     /* Local variables */
 157     double f;
 158     int i__, j, k, ih, jp, ip, iq, np, iw;
 159     double xp, xq, den;
 160     int ihp;
 161     double one;
 162     int ihq, jpn, kpt;
 163     double sum, diff, half, beta;
 164     int kold;
 165     double winc;
 166     int nrem, knew;
 167     double temp, bsum;
 168     int nptm;
 169     double zero, hdiag, fbase, sfrac, denom, vquad, sumpq;
 170     double dsqmin, distsq, vlmxsq;
 171
 172
 173 /*     The arguments N, NPT, XL, XU, MAXFUN, XBASE, XPT, FVAL, XOPT, */
 174 /*       GOPT, HQ, PQ, BMAT, ZMAT, NDIM, SL and SU have the same meanings as */
 175 /*       the corresponding arguments of BOBYQB on the entry to RESCUE. */
 176 /*     NF is maintained as the number of calls of CALFUN so far, except that */
 177 /*       NF is set to -1 if the value of MAXFUN prevents further progress. */
 178 /*     KOPT is maintained so that FVAL(KOPT) is the least calculated function */
 179 /*       value. Its correct value must be given on entry. It is updated if a */
 180 /*       new least function value is found, but the corresponding changes to */
 181 /*       XOPT and GOPT have to be made later by the calling program. */
 182 /*     DELTA is the current trust region radius. */
 183 /*     VLAG is a working space vector that will be used for the values of the */
 184 /*       provisional Lagrange functions at each of the interpolation points. */
 185 /*       They are part of a product that requires VLAG to be of length NDIM. */
 186 /*     PTSAUX is also a working space array. For J=1,2,...,N, PTSAUX(1,J) and */
 187 /*       PTSAUX(2,J) specify the two positions of provisional interpolation */
 188 /*       points when a nonzero step is taken along e_J (the J-th coordinate */
 189 /*       direction) through XBASE+XOPT, as specified below. Usually these */
 190 /*       steps have length DELTA, but other lengths are chosen if necessary */
 191 /*       in order to satisfy the given bounds on the variables. */
 192 /*     PTSID is also a working space array. It has NPT components that denote */
 193 /*       provisional new positions of the original interpolation points, in */
 194 /*       case changes are needed to restore the linear independence of the */
 195 /*       interpolation conditions. The K-th point is a candidate for change */
 196 /*       if and only if PTSID(K) is nonzero. In this case let p and q be the */
 197 /*       int parts of PTSID(K) and (PTSID(K)-p) multiplied by N+1. If p */
 198 /*       and q are both positive, the step from XBASE+XOPT to the new K-th */
 199 /*       interpolation point is PTSAUX(1,p)*e_p + PTSAUX(1,q)*e_q. Otherwise */
 200 /*       the step is PTSAUX(1,p)*e_p or PTSAUX(2,q)*e_q in the cases q=0 or */
 201 /*       p=0, respectively. */
 202 /*     The first NDIM+NPT elements of the array W are used for working space. */
 203 /*     The final elements of BMAT and ZMAT are set in a well-conditioned way */
 204 /*       to the values that are appropriate for the new interpolation points. */
 205 /*     The elements of GOPT, HQ and PQ are also revised to the values that are */
 206 /*       appropriate to the final quadratic model. */
 207
 208 /*     Set some constants. */
 209
 210     /* Parameter adjustments */
 211     zmat_dim1 = *npt;
 212     zmat_offset = 1 + zmat_dim1;
 213     zmat -= zmat_offset;
 214     xpt_dim1 = *npt;
 215     xpt_offset = 1 + xpt_dim1;
 216     xpt -= xpt_offset;
 217     --xl;
 218     --xu;
 219     --xbase;
 220     --fval;
 221     --xopt;
 222     --gopt;
 223     --hq;
 224     --pq;
 225     bmat_dim1 = *ndim;
 226     bmat_offset = 1 + bmat_dim1;
 227     bmat -= bmat_offset;
 228     --sl;
 229     --su;
 230     --vlag;
 231     ptsaux -= 3;
 232     --ptsid;
 233     --w;
 234
 235     /* Function Body */
 236     half = .5;
 237     one = 1.;
 238     zero = 0.;
 239     np = *n + 1;
 240     sfrac = half / (double) np;
 241     nptm = *npt - np;
 242
 243 /*     Shift the interpolation points so that XOPT becomes the origin, and set */
 244 /*     the elements of ZMAT to zero. The value of SUMPQ is required in the */
 245 /*     updating of HQ below. The squares of the distances from XOPT to the */
 246 /*     other interpolation points are set at the end of W. Increments of WINC */
 247 /*     may be added later to these squares to balance the consideration of */
 248 /*     the choice of point that is going to become current. */
 249
 250     sumpq = zero;
 251     winc = zero;
 252     i__1 = *npt;
 253     for (k = 1; k <= i__1; ++k) {
 254         distsq = zero;
 255         i__2 = *n;
 256         for (j = 1; j <= i__2; ++j) {
 257             xpt[k + j * xpt_dim1] -= xopt[j];
 258 /* L10: */
 259 /* Computing 2nd power */
 260             d__1 = xpt[k + j * xpt_dim1];
 261             distsq += d__1 * d__1;
 262         }
 263         sumpq += pq[k];
 264         w[*ndim + k] = distsq;
 265         winc = max(winc,distsq);
 266         i__2 = nptm;
 267         for (j = 1; j <= i__2; ++j) {
 268 /* L20: */
 269             zmat[k + j * zmat_dim1] = zero;
 270         }
 271     }
 272
 273 /*     Update HQ so that HQ and PQ define the second derivatives of the model */
 274 /*     after XBASE has been shifted to the trust region centre. */
 275
 276     ih = 0;
 277     i__2 = *n;
 278     for (j = 1; j <= i__2; ++j) {
 279         w[j] = half * sumpq * xopt[j];
 280         i__1 = *npt;
 281         for (k = 1; k <= i__1; ++k) {
 282 /* L30: */
 283             w[j] += pq[k] * xpt[k + j * xpt_dim1];
 284         }
 285         i__1 = j;
 286         for (i__ = 1; i__ <= i__1; ++i__) {
 287             ++ih;
 288 /* L40: */
 289             hq[ih] = hq[ih] + w[i__] * xopt[j] + w[j] * xopt[i__];
 290         }
 291     }
 292
 293 /*     Shift XBASE, SL, SU and XOPT. Set the elements of BMAT to zero, and */
 294 /*     also set the elements of PTSAUX. */
 295
 296     i__1 = *n;
 297     for (j = 1; j <= i__1; ++j) {
 298         xbase[j] += xopt[j];
 299         sl[j] -= xopt[j];
 300         su[j] -= xopt[j];
 301         xopt[j] = zero;
 302 /* Computing MIN */
 303         d__1 = *delta, d__2 = su[j];
 304         ptsaux[(j << 1) + 1] = min(d__1,d__2);
 305 /* Computing MAX */
 306         d__1 = -(*delta), d__2 = sl[j];
 307         ptsaux[(j << 1) + 2] = max(d__1,d__2);
 308         if (ptsaux[(j << 1) + 1] + ptsaux[(j << 1) + 2] < zero) {
 309             temp = ptsaux[(j << 1) + 1];
 310             ptsaux[(j << 1) + 1] = ptsaux[(j << 1) + 2];
 311             ptsaux[(j << 1) + 2] = temp;
 312         }
 313         if ((d__2 = ptsaux[(j << 1) + 2], fabs(d__2)) < half * (d__1 = ptsaux[(
 314                 j << 1) + 1], fabs(d__1))) {
 315             ptsaux[(j << 1) + 2] = half * ptsaux[(j << 1) + 1];
 316         }
 317         i__2 = *ndim;
 318         for (i__ = 1; i__ <= i__2; ++i__) {
 319 /* L50: */
 320             bmat[i__ + j * bmat_dim1] = zero;
 321         }
 322     }
 323     fbase = fval[*kopt];
 324
 325 /*     Set the identifiers of the artificial interpolation points that are */
 326 /*     along a coordinate direction from XOPT, and set the corresponding */
 327 /*     nonzero elements of BMAT and ZMAT. */
 328
 329     ptsid[1] = sfrac;
 330     i__2 = *n;
 331     for (j = 1; j <= i__2; ++j) {
 332         jp = j + 1;
 333         jpn = jp + *n;
 334         ptsid[jp] = (double) j + sfrac;
 335         if (jpn <= *npt) {
 336             ptsid[jpn] = (double) j / (double) np + sfrac;
 337             temp = one / (ptsaux[(j << 1) + 1] - ptsaux[(j << 1) + 2]);
 338             bmat[jp + j * bmat_dim1] = -temp + one / ptsaux[(j << 1) + 1];
 339             bmat[jpn + j * bmat_dim1] = temp + one / ptsaux[(j << 1) + 2];
 340             bmat[j * bmat_dim1 + 1] = -bmat[jp + j * bmat_dim1] - bmat[jpn +
 341                     j * bmat_dim1];
 342             zmat[j * zmat_dim1 + 1] = sqrt(2.) / (d__1 = ptsaux[(j << 1) + 1]
 343                     * ptsaux[(j << 1) + 2], fabs(d__1));
 344             zmat[jp + j * zmat_dim1] = zmat[j * zmat_dim1 + 1] * ptsaux[(j <<
 345                     1) + 2] * temp;
 346             zmat[jpn + j * zmat_dim1] = -zmat[j * zmat_dim1 + 1] * ptsaux[(j
 347                     << 1) + 1] * temp;
 348         } else {
 349             bmat[j * bmat_dim1 + 1] = -one / ptsaux[(j << 1) + 1];
 350             bmat[jp + j * bmat_dim1] = one / ptsaux[(j << 1) + 1];
 351 /* Computing 2nd power */
 352             d__1 = ptsaux[(j << 1) + 1];
 353             bmat[j + *npt + j * bmat_dim1] = -half * (d__1 * d__1);
 354         }
 355 /* L60: */
 356     }
 357
 358 /*     Set any remaining identifiers with their nonzero elements of ZMAT. */
 359
 360     if (*npt >= *n + np) {
 361         i__2 = *npt;
 362         for (k = np << 1; k <= i__2; ++k) {
 363             iw = (int) (((double) (k - np) - half) / (double) (*n)
 364                     );
 365             ip = k - np - iw * *n;
 366             iq = ip + iw;
 367             if (iq > *n) {
 368                 iq -= *n;
 369             }
 370             ptsid[k] = (double) ip + (double) iq / (double) np +
 371                     sfrac;
 372             temp = one / (ptsaux[(ip << 1) + 1] * ptsaux[(iq << 1) + 1]);
 373             zmat[(k - np) * zmat_dim1 + 1] = temp;
 374             zmat[ip + 1 + (k - np) * zmat_dim1] = -temp;
 375             zmat[iq + 1 + (k - np) * zmat_dim1] = -temp;
 376 /* L70: */
 377             zmat[k + (k - np) * zmat_dim1] = temp;
 378         }
 379     }
 380     nrem = *npt;
 381     kold = 1;
 382     knew = *kopt;
 383
 384 /*     Reorder the provisional points in the way that exchanges PTSID(KOLD) */
 385 /*     with PTSID(KNEW). */
 386
 387 L80:
 388     i__2 = *n;
 389     for (j = 1; j <= i__2; ++j) {
 390         temp = bmat[kold + j * bmat_dim1];
 391         bmat[kold + j * bmat_dim1] = bmat[knew + j * bmat_dim1];
 392 /* L90: */
 393         bmat[knew + j * bmat_dim1] = temp;
 394     }
 395     i__2 = nptm;
 396     for (j = 1; j <= i__2; ++j) {
 397         temp = zmat[kold + j * zmat_dim1];
 398         zmat[kold + j * zmat_dim1] = zmat[knew + j * zmat_dim1];
 399 /* L100: */
 400         zmat[knew + j * zmat_dim1] = temp;
 401     }
 402     ptsid[kold] = ptsid[knew];
 403     ptsid[knew] = zero;
 404     w[*ndim + knew] = zero;
 405     --nrem;
 406     if (knew != *kopt) {
 407         temp = vlag[kold];
 408         vlag[kold] = vlag[knew];
 409         vlag[knew] = temp;
 410
 411 /*     Update the BMAT and ZMAT matrices so that the status of the KNEW-th */
 412 /*     interpolation point can be changed from provisional to original. The */
 413 /*     branch to label 350 occurs if all the original points are reinstated. */
 414 /*     The nonnegative values of W(NDIM+K) are required in the search below. */
 415
 416         update_(n, npt, &bmat[bmat_offset], &zmat[zmat_offset], ndim, &vlag[1]
 417                 , &beta, &denom, &knew, &w[1]);
 418         if (nrem == 0) {
 419             goto L350;
 420         }
 421         i__2 = *npt;
 422         for (k = 1; k <= i__2; ++k) {
 423 /* L110: */
 424             w[*ndim + k] = (d__1 = w[*ndim + k], fabs(d__1));
 425         }
 426     }
 427
 428 /*     Pick the index KNEW of an original interpolation point that has not */
 429 /*     yet replaced one of the provisional interpolation points, giving */
 430 /*     attention to the closeness to XOPT and to previous tries with KNEW. */
 431
 432 L120:
 433     dsqmin = zero;
 434     i__2 = *npt;
 435     for (k = 1; k <= i__2; ++k) {
 436         if (w[*ndim + k] > zero) {
 437             if (dsqmin == zero || w[*ndim + k] < dsqmin) {
 438                 knew = k;
 439                 dsqmin = w[*ndim + k];
 440             }
 441         }
 442 /* L130: */
 443     }
 444     if (dsqmin == zero) {
 445         goto L260;
 446     }
 447
 448 /*     Form the W-vector of the chosen original interpolation point. */
 449
 450     i__2 = *n;
 451     for (j = 1; j <= i__2; ++j) {
 452 /* L140: */
 453         w[*npt + j] = xpt[knew + j * xpt_dim1];
 454     }
 455     i__2 = *npt;
 456     for (k = 1; k <= i__2; ++k) {
 457         sum = zero;
 458         if (k == *kopt) {
 459         } else if (ptsid[k] == zero) {
 460             i__1 = *n;
 461             for (j = 1; j <= i__1; ++j) {
 462 /* L150: */
 463                 sum += w[*npt + j] * xpt[k + j * xpt_dim1];
 464             }
 465         } else {
 466             ip = (int) ptsid[k];
 467             if (ip > 0) {
 468                 sum = w[*npt + ip] * ptsaux[(ip << 1) + 1];
 469             }
 470             iq = (int) ((double) np * ptsid[k] - (double) (ip *
 471                     np));
 472             if (iq > 0) {
 473                 iw = 1;
 474                 if (ip == 0) {
 475                     iw = 2;
 476                 }
 477                 sum += w[*npt + iq] * ptsaux[iw + (iq << 1)];
 478             }
 479         }
 480 /* L160: */
 481         w[k] = half * sum * sum;
 482     }
 483
 484 /*     Calculate VLAG and BETA for the required updating of the H matrix if */
 485 /*     XPT(KNEW,.) is reinstated in the set of interpolation points. */
 486
 487     i__2 = *npt;
 488     for (k = 1; k <= i__2; ++k) {
 489         sum = zero;
 490         i__1 = *n;
 491         for (j = 1; j <= i__1; ++j) {
 492 /* L170: */
 493             sum += bmat[k + j * bmat_dim1] * w[*npt + j];
 494         }
 495 /* L180: */
 496         vlag[k] = sum;
 497     }
 498     beta = zero;
 499     i__2 = nptm;
 500     for (j = 1; j <= i__2; ++j) {
 501         sum = zero;
 502         i__1 = *npt;
 503         for (k = 1; k <= i__1; ++k) {
 504 /* L190: */
 505             sum += zmat[k + j * zmat_dim1] * w[k];
 506         }
 507         beta -= sum * sum;
 508         i__1 = *npt;
 509         for (k = 1; k <= i__1; ++k) {
 510 /* L200: */
 511             vlag[k] += sum * zmat[k + j * zmat_dim1];
 512         }
 513     }
 514     bsum = zero;
 515     distsq = zero;
 516     i__1 = *n;
 517     for (j = 1; j <= i__1; ++j) {
 518         sum = zero;
 519         i__2 = *npt;
 520         for (k = 1; k <= i__2; ++k) {
 521 /* L210: */
 522             sum += bmat[k + j * bmat_dim1] * w[k];
 523         }
 524         jp = j + *npt;
 525         bsum += sum * w[jp];
 526         i__2 = *ndim;
 527         for (ip = *npt + 1; ip <= i__2; ++ip) {
 528 /* L220: */
 529             sum += bmat[ip + j * bmat_dim1] * w[ip];
 530         }
 531         bsum += sum * w[jp];
 532         vlag[jp] = sum;
 533 /* L230: */
 534 /* Computing 2nd power */
 535         d__1 = xpt[knew + j * xpt_dim1];
 536         distsq += d__1 * d__1;
 537     }
 538     beta = half * distsq * distsq + beta - bsum;
 539     vlag[*kopt] += one;
 540
 541 /*     KOLD is set to the index of the provisional interpolation point that is */
 542 /*     going to be deleted to make way for the KNEW-th original interpolation */
 543 /*     point. The choice of KOLD is governed by the avoidance of a small value */
 544 /*     of the denominator in the updating calculation of UPDATE. */
 545
 546     denom = zero;
 547     vlmxsq = zero;
 548     i__1 = *npt;
 549     for (k = 1; k <= i__1; ++k) {
 550         if (ptsid[k] != zero) {
 551             hdiag = zero;
 552             i__2 = nptm;
 553             for (j = 1; j <= i__2; ++j) {
 554 /* L240: */
 555 /* Computing 2nd power */
 556                 d__1 = zmat[k + j * zmat_dim1];
 557                 hdiag += d__1 * d__1;
 558             }
 559 /* Computing 2nd power */
 560             d__1 = vlag[k];
 561             den = beta * hdiag + d__1 * d__1;
 562             if (den > denom) {
 563                 kold = k;
 564                 denom = den;
 565             }
 566         }
 567 /* L250: */
 568 /* Computing MAX */
 569 /* Computing 2nd power */
 570         d__3 = vlag[k];
 571         d__1 = vlmxsq, d__2 = d__3 * d__3;
 572         vlmxsq = max(d__1,d__2);
 573     }
 574     if (denom <= vlmxsq * .01) {
 575         w[*ndim + knew] = -w[*ndim + knew] - winc;
 576         goto L120;
 577     }
 578     goto L80;
 579
 580 /*     When label 260 is reached, all the final positions of the interpolation */
 581 /*     points have been chosen although any changes have not been included yet */
 582 /*     in XPT. Also the final BMAT and ZMAT matrices are complete, but, apart */
 583 /*     from the shift of XBASE, the updating of the quadratic model remains to */
 584 /*     be done. The following cycle through the new interpolation points begins */
 585 /*     by putting the new point in XPT(KPT,.) and by setting PQ(KPT) to zero, */
 586 /*     except that a RETURN occurs if MAXFUN prohibits another value of F. */
 587
 588 L260:
 589     i__1 = *npt;
 590     for (kpt = 1; kpt <= i__1; ++kpt) {
 591         if (ptsid[kpt] == zero) {
 592             goto L340;
 593         }
 594
 595         if (nlopt_stop_forced(stop)) return NLOPT_FORCE_STOP;
 596         else if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
 597         else if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
 598
 599         ih = 0;
 600         i__2 = *n;
 601         for (j = 1; j <= i__2; ++j) {
 602             w[j] = xpt[kpt + j * xpt_dim1];
 603             xpt[kpt + j * xpt_dim1] = zero;
 604             temp = pq[kpt] * w[j];
 605             i__3 = j;
 606             for (i__ = 1; i__ <= i__3; ++i__) {
 607                 ++ih;
 608 /* L270: */
 609                 hq[ih] += temp * w[i__];
 610             }
 611         }
 612         pq[kpt] = zero;
 613         ip = (int) ptsid[kpt];
 614         iq = (int) ((double) np * ptsid[kpt] - (double) (ip * np))
 615                 ;
 616         if (ip > 0) {
 617             xp = ptsaux[(ip << 1) + 1];
 618             xpt[kpt + ip * xpt_dim1] = xp;
 619         }
 620         if (iq > 0) {
 621             xq = ptsaux[(iq << 1) + 1];
 622             if (ip == 0) {
 623                 xq = ptsaux[(iq << 1) + 2];
 624             }
 625             xpt[kpt + iq * xpt_dim1] = xq;
 626         }
 627
 628 /*     Set VQUAD to the value of the current model at the new point. */
 629
 630         vquad = fbase;
 631         if (ip > 0) {
 632             ihp = (ip + ip * ip) / 2;
 633             vquad += xp * (gopt[ip] + half * xp * hq[ihp]);
 634         }
 635         if (iq > 0) {
 636             ihq = (iq + iq * iq) / 2;
 637             vquad += xq * (gopt[iq] + half * xq * hq[ihq]);
 638             if (ip > 0) {
 639                 iw = max(ihp,ihq) - (i__3 = ip - iq, iabs(i__3));
 640                 vquad += xp * xq * hq[iw];
 641             }
 642         }
 643         i__3 = *npt;
 644         for (k = 1; k <= i__3; ++k) {
 645             temp = zero;
 646             if (ip > 0) {
 647                 temp += xp * xpt[k + ip * xpt_dim1];
 648             }
 649             if (iq > 0) {
 650                 temp += xq * xpt[k + iq * xpt_dim1];
 651             }
 652 /* L280: */
 653             vquad += half * pq[k] * temp * temp;
 654         }
 655
 656 /*     Calculate F at the new interpolation point, and set DIFF to the factor */
 657 /*     that is going to multiply the KPT-th Lagrange function when the model */
 658 /*     is updated to provide interpolation to the new function value. */
 659
 660         i__3 = *n;
 661         for (i__ = 1; i__ <= i__3; ++i__) {
 662 /* Computing MIN */
 663 /* Computing MAX */
 664             d__3 = xl[i__], d__4 = xbase[i__] + xpt[kpt + i__ * xpt_dim1];
 665             d__1 = max(d__3,d__4), d__2 = xu[i__];
 666             w[i__] = min(d__1,d__2);
 667             if (xpt[kpt + i__ * xpt_dim1] == sl[i__]) {
 668                 w[i__] = xl[i__];
 669             }
 670             if (xpt[kpt + i__ * xpt_dim1] == su[i__]) {
 671                 w[i__] = xu[i__];
 672             }
 673 /* L290: */
 674         }
 675
 676         stop->nevals++;
 677         f = calfun(*n, &w[1], calfun_data);
 678         fval[kpt] = f;
 679         if (f < fval[*kopt]) {
 680             *kopt = kpt;
 681         }
 682         if (nlopt_stop_forced(stop)) return NLOPT_FORCE_STOP;
 683         else if (f < stop->minf_max) return NLOPT_MINF_MAX_REACHED;
 684         else if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
 685         else if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
 686
 687         diff = f - vquad;
 688
 689 /*     Update the quadratic model. The RETURN from the subroutine occurs when */
 690 /*     all the new interpolation points are included in the model. */
 691
 692         i__3 = *n;
 693         for (i__ = 1; i__ <= i__3; ++i__) {
 694 /* L310: */
 695             gopt[i__] += diff * bmat[kpt + i__ * bmat_dim1];
 696         }
 697         i__3 = *npt;
 698         for (k = 1; k <= i__3; ++k) {
 699             sum = zero;
 700             i__2 = nptm;
 701             for (j = 1; j <= i__2; ++j) {
 702 /* L320: */
 703                 sum += zmat[k + j * zmat_dim1] * zmat[kpt + j * zmat_dim1];
 704             }
 705             temp = diff * sum;
 706             if (ptsid[k] == zero) {
 707                 pq[k] += temp;
 708             } else {
 709                 ip = (int) ptsid[k];
 710                 iq = (int) ((double) np * ptsid[k] - (double) (ip
 711                         * np));
 712                 ihq = (iq * iq + iq) / 2;
 713                 if (ip == 0) {
 714 /* Computing 2nd power */
 715                     d__1 = ptsaux[(iq << 1) + 2];
 716                     hq[ihq] += temp * (d__1 * d__1);
 717                 } else {
 718                     ihp = (ip * ip + ip) / 2;
 719 /* Computing 2nd power */
 720                     d__1 = ptsaux[(ip << 1) + 1];
 721                     hq[ihp] += temp * (d__1 * d__1);
 722                     if (iq > 0) {
 723 /* Computing 2nd power */
 724                         d__1 = ptsaux[(iq << 1) + 1];
 725                         hq[ihq] += temp * (d__1 * d__1);
 726                         iw = max(ihp,ihq) - (i__2 = iq - ip, iabs(i__2));
 727                         hq[iw] += temp * ptsaux[(ip << 1) + 1] * ptsaux[(iq <<
 728                                  1) + 1];
 729                     }
 730                 }
 731             }
 732 /* L330: */
 733         }
 734         ptsid[kpt] = zero;
 735 L340:
 736         ;
 737     }
 738 L350:
 739     return NLOPT_SUCCESS;
 740 } /* rescue_ */
 741
 742 static void altmov_(int *n, int *npt, double *xpt,
 743         double *xopt, double *bmat, double *zmat, int *ndim,
 744         double *sl, double *su, int *kopt, int *knew,
 745         double *adelt, double *xnew, double *xalt, double *
 746         alpha, double *cauchy, double *glag, double *hcol,
 747         double *w)
 748 {
 749     /* System generated locals */
 750     int xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1,
 751             zmat_offset, i__1, i__2;
 752     double d__1, d__2, d__3, d__4;
 753
 754     /* Local variables */
 755     int i__, j, k;
 756     double ha, gw, one, diff, half;
 757     int ilbd, isbd;
 758     double slbd;
 759     int iubd;
 760     double vlag, subd, temp;
 761     int ksav;
 762     double step, zero, curv;
 763     int iflag;
 764     double scale, csave, tempa, tempb, tempd, const__, sumin, ggfree;
 765     int ibdsav;
 766     double dderiv, bigstp, predsq, presav, distsq, stpsav, wfixsq, wsqsav;
 767
 768
 769 /*     The arguments N, NPT, XPT, XOPT, BMAT, ZMAT, NDIM, SL and SU all have */
 770 /*       the same meanings as the corresponding arguments of BOBYQB. */
 771 /*     KOPT is the index of the optimal interpolation point. */
 772 /*     KNEW is the index of the interpolation point that is going to be moved. */
 773 /*     ADELT is the current trust region bound. */
 774 /*     XNEW will be set to a suitable new position for the interpolation point */
 775 /*       XPT(KNEW,.). Specifically, it satisfies the SL, SU and trust region */
 776 /*       bounds and it should provide a large denominator in the next call of */
 777 /*       UPDATE. The step XNEW-XOPT from XOPT is restricted to moves along the */
 778 /*       straight lines through XOPT and another interpolation point. */
 779 /*     XALT also provides a large value of the modulus of the KNEW-th Lagrange */
 780 /*       function subject to the constraints that have been mentioned, its main */
 781 /*       difference from XNEW being that XALT-XOPT is a constrained version of */
 782 /*       the Cauchy step within the trust region. An exception is that XALT is */
 783 /*       not calculated if all components of GLAG (see below) are zero. */
 784 /*     ALPHA will be set to the KNEW-th diagonal element of the H matrix. */
 785 /*     CAUCHY will be set to the square of the KNEW-th Lagrange function at */
 786 /*       the step XALT-XOPT from XOPT for the vector XALT that is returned, */
 787 /*       except that CAUCHY is set to zero if XALT is not calculated. */
 788 /*     GLAG is a working space vector of length N for the gradient of the */
 789 /*       KNEW-th Lagrange function at XOPT. */
 790 /*     HCOL is a working space vector of length NPT for the second derivative */
 791 /*       coefficients of the KNEW-th Lagrange function. */
 792 /*     W is a working space vector of length 2N that is going to hold the */
 793 /*       constrained Cauchy step from XOPT of the Lagrange function, followed */
 794 /*       by the downhill version of XALT when the uphill step is calculated. */
 795
 796 /*     Set the first NPT components of W to the leading elements of the */
 797 /*     KNEW-th column of the H matrix. */
 798
 799     /* Parameter adjustments */
 800     zmat_dim1 = *npt;
 801     zmat_offset = 1 + zmat_dim1;
 802     zmat -= zmat_offset;
 803     xpt_dim1 = *npt;
 804     xpt_offset = 1 + xpt_dim1;
 805     xpt -= xpt_offset;
 806     --xopt;
 807     bmat_dim1 = *ndim;
 808     bmat_offset = 1 + bmat_dim1;
 809     bmat -= bmat_offset;
 810     --sl;
 811     --su;
 812     --xnew;
 813     --xalt;
 814     --glag;
 815     --hcol;
 816     --w;
 817
 818     /* Function Body */
 819     half = .5;
 820     one = 1.;
 821     zero = 0.;
 822     const__ = one + sqrt(2.);
 823     i__1 = *npt;
 824     for (k = 1; k <= i__1; ++k) {
 825 /* L10: */
 826         hcol[k] = zero;
 827     }
 828     i__1 = *npt - *n - 1;
 829     for (j = 1; j <= i__1; ++j) {
 830         temp = zmat[*knew + j * zmat_dim1];
 831         i__2 = *npt;
 832         for (k = 1; k <= i__2; ++k) {
 833 /* L20: */
 834             hcol[k] += temp * zmat[k + j * zmat_dim1];
 835         }
 836     }
 837     *alpha = hcol[*knew];
 838     ha = half * *alpha;
 839
 840 /*     Calculate the gradient of the KNEW-th Lagrange function at XOPT. */
 841
 842     i__2 = *n;
 843     for (i__ = 1; i__ <= i__2; ++i__) {
 844 /* L30: */
 845         glag[i__] = bmat[*knew + i__ * bmat_dim1];
 846     }
 847     i__2 = *npt;
 848     for (k = 1; k <= i__2; ++k) {
 849         temp = zero;
 850         i__1 = *n;
 851         for (j = 1; j <= i__1; ++j) {
 852 /* L40: */
 853             temp += xpt[k + j * xpt_dim1] * xopt[j];
 854         }
 855         temp = hcol[k] * temp;
 856         i__1 = *n;
 857         for (i__ = 1; i__ <= i__1; ++i__) {
 858 /* L50: */
 859             glag[i__] += temp * xpt[k + i__ * xpt_dim1];
 860         }
 861     }
 862
 863 /*     Search for a large denominator along the straight lines through XOPT */
 864 /*     and another interpolation point. SLBD and SUBD will be lower and upper */
 865 /*     bounds on the step along each of these lines in turn. PREDSQ will be */
 866 /*     set to the square of the predicted denominator for each line. PRESAV */
 867 /*     will be set to the largest admissible value of PREDSQ that occurs. */
 868
 869     presav = zero;
 870     i__1 = *npt;
 871     for (k = 1; k <= i__1; ++k) {
 872         if (k == *kopt) {
 873             goto L80;
 874         }
 875         dderiv = zero;
 876         distsq = zero;
 877         i__2 = *n;
 878         for (i__ = 1; i__ <= i__2; ++i__) {
 879             temp = xpt[k + i__ * xpt_dim1] - xopt[i__];
 880             dderiv += glag[i__] * temp;
 881 /* L60: */
 882             distsq += temp * temp;
 883         }
 884         subd = *adelt / sqrt(distsq);
 885         slbd = -subd;
 886         ilbd = 0;
 887         iubd = 0;
 888         sumin = min(one,subd);
 889
 890 /*     Revise SLBD and SUBD if necessary because of the bounds in SL and SU. */
 891
 892         i__2 = *n;
 893         for (i__ = 1; i__ <= i__2; ++i__) {
 894             temp = xpt[k + i__ * xpt_dim1] - xopt[i__];
 895             if (temp > zero) {
 896                 if (slbd * temp < sl[i__] - xopt[i__]) {
 897                     slbd = (sl[i__] - xopt[i__]) / temp;
 898                     ilbd = -i__;
 899                 }
 900                 if (subd * temp > su[i__] - xopt[i__]) {
 901 /* Computing MAX */
 902                     d__1 = sumin, d__2 = (su[i__] - xopt[i__]) / temp;
 903                     subd = max(d__1,d__2);
 904                     iubd = i__;
 905                 }
 906             } else if (temp < zero) {
 907                 if (slbd * temp > su[i__] - xopt[i__]) {
 908                     slbd = (su[i__] - xopt[i__]) / temp;
 909                     ilbd = i__;
 910                 }
 911                 if (subd * temp < sl[i__] - xopt[i__]) {
 912 /* Computing MAX */
 913                     d__1 = sumin, d__2 = (sl[i__] - xopt[i__]) / temp;
 914                     subd = max(d__1,d__2);
 915                     iubd = -i__;
 916                 }
 917             }
 918 /* L70: */
 919         }
 920
 921 /*     Seek a large modulus of the KNEW-th Lagrange function when the index */
 922 /*     of the other interpolation point on the line through XOPT is KNEW. */
 923
 924         if (k == *knew) {
 925             diff = dderiv - one;
 926             step = slbd;
 927             vlag = slbd * (dderiv - slbd * diff);
 928             isbd = ilbd;
 929             temp = subd * (dderiv - subd * diff);
 930             if (fabs(temp) > fabs(vlag)) {
 931                 step = subd;
 932                 vlag = temp;
 933                 isbd = iubd;
 934             }
 935             tempd = half * dderiv;
 936             tempa = tempd - diff * slbd;
 937             tempb = tempd - diff * subd;
 938             if (tempa * tempb < zero) {
 939                 temp = tempd * tempd / diff;
 940                 if (fabs(temp) > fabs(vlag)) {
 941                     step = tempd / diff;
 942                     vlag = temp;
 943                     isbd = 0;
 944                 }
 945             }
 946
 947 /*     Search along each of the other lines through XOPT and another point. */
 948
 949         } else {
 950             step = slbd;
 951             vlag = slbd * (one - slbd);
 952             isbd = ilbd;
 953             temp = subd * (one - subd);
 954             if (fabs(temp) > fabs(vlag)) {
 955                 step = subd;
 956                 vlag = temp;
 957                 isbd = iubd;
 958             }
 959             if (subd > half) {
 960                 if (fabs(vlag) < .25) {
 961                     step = half;
 962                     vlag = .25;
 963                     isbd = 0;
 964                 }
 965             }
 966             vlag *= dderiv;
 967         }
 968
 969 /*     Calculate PREDSQ for the current line search and maintain PRESAV. */
 970
 971         temp = step * (one - step) * distsq;
 972         predsq = vlag * vlag * (vlag * vlag + ha * temp * temp);
 973         if (predsq > presav) {
 974             presav = predsq;
 975             ksav = k;
 976             stpsav = step;
 977             ibdsav = isbd;
 978         }
 979 L80:
 980         ;
 981     }
 982
 983 /*     Construct XNEW in a way that satisfies the bound constraints exactly. */
 984
 985     i__1 = *n;
 986     for (i__ = 1; i__ <= i__1; ++i__) {
 987         temp = xopt[i__] + stpsav * (xpt[ksav + i__ * xpt_dim1] - xopt[i__]);
 988 /* L90: */
 989 /* Computing MAX */
 990 /* Computing MIN */
 991         d__3 = su[i__];
 992         d__1 = sl[i__], d__2 = min(d__3,temp);
 993         xnew[i__] = max(d__1,d__2);
 994     }
 995     if (ibdsav < 0) {
 996         xnew[-ibdsav] = sl[-ibdsav];
 997     }
 998     if (ibdsav > 0) {
 999         xnew[ibdsav] = su[ibdsav];
1000     }
1001
1002 /*     Prepare for the iterative method that assembles the constrained Cauchy */
1003 /*     step in W. The sum of squares of the fixed components of W is formed in */
1004 /*     WFIXSQ, and the free components of W are set to BIGSTP. */
1005
1006     bigstp = *adelt + *adelt;
1007     iflag = 0;
1008 L100:
1009     wfixsq = zero;
1010     ggfree = zero;
1011     i__1 = *n;
1012     for (i__ = 1; i__ <= i__1; ++i__) {
1013         w[i__] = zero;
1014 /* Computing MIN */
1015         d__1 = xopt[i__] - sl[i__], d__2 = glag[i__];
1016         tempa = min(d__1,d__2);
1017 /* Computing MAX */
1018         d__1 = xopt[i__] - su[i__], d__2 = glag[i__];
1019         tempb = max(d__1,d__2);
1020         if (tempa > zero || tempb < zero) {
1021             w[i__] = bigstp;
1022 /* Computing 2nd power */
1023             d__1 = glag[i__];
1024             ggfree += d__1 * d__1;
1025         }
1026 /* L110: */
1027     }
1028     if (ggfree == zero) {
1029         *cauchy = zero;
1030         goto L200;
1031     }
1032
1033 /*     Investigate whether more components of W can be fixed. */
1034
1035 L120:
1036     temp = *adelt * *adelt - wfixsq;
1037     if (temp > zero) {
1038         wsqsav = wfixsq;
1039         step = sqrt(temp / ggfree);
1040         ggfree = zero;
1041         i__1 = *n;
1042         for (i__ = 1; i__ <= i__1; ++i__) {
1043             if (w[i__] == bigstp) {
1044                 temp = xopt[i__] - step * glag[i__];
1045                 if (temp <= sl[i__]) {
1046                     w[i__] = sl[i__] - xopt[i__];
1047 /* Computing 2nd power */
1048                     d__1 = w[i__];
1049                     wfixsq += d__1 * d__1;
1050                 } else if (temp >= su[i__]) {
1051                     w[i__] = su[i__] - xopt[i__];
1052 /* Computing 2nd power */
1053                     d__1 = w[i__];
1054                     wfixsq += d__1 * d__1;
1055                 } else {
1056 /* Computing 2nd power */
1057                     d__1 = glag[i__];
1058                     ggfree += d__1 * d__1;
1059                 }
1060             }
1061 /* L130: */
1062         }
1063         if (wfixsq > wsqsav && ggfree > zero) {
1064             goto L120;
1065         }
1066     }
1067
1068 /*     Set the remaining free components of W and all components of XALT, */
1069 /*     except that W may be scaled later. */
1070
1071     gw = zero;
1072     i__1 = *n;
1073     for (i__ = 1; i__ <= i__1; ++i__) {
1074         if (w[i__] == bigstp) {
1075             w[i__] = -step * glag[i__];
1076 /* Computing MAX */
1077 /* Computing MIN */
1078             d__3 = su[i__], d__4 = xopt[i__] + w[i__];
1079             d__1 = sl[i__], d__2 = min(d__3,d__4);
1080             xalt[i__] = max(d__1,d__2);
1081         } else if (w[i__] == zero) {
1082             xalt[i__] = xopt[i__];
1083         } else if (glag[i__] > zero) {
1084             xalt[i__] = sl[i__];
1085         } else {
1086             xalt[i__] = su[i__];
1087         }
1088 /* L140: */
1089         gw += glag[i__] * w[i__];
1090     }
1091
1092 /*     Set CURV to the curvature of the KNEW-th Lagrange function along W. */
1093 /*     Scale W by a factor less than one if that can reduce the modulus of */
1094 /*     the Lagrange function at XOPT+W. Set CAUCHY to the final value of */
1095 /*     the square of this function. */
1096
1097     curv = zero;
1098     i__1 = *npt;
1099     for (k = 1; k <= i__1; ++k) {
1100         temp = zero;
1101         i__2 = *n;
1102         for (j = 1; j <= i__2; ++j) {
1103 /* L150: */
1104             temp += xpt[k + j * xpt_dim1] * w[j];
1105         }
1106 /* L160: */
1107         curv += hcol[k] * temp * temp;
1108     }
1109     if (iflag == 1) {
1110         curv = -curv;
1111     }
1112     if (curv > -gw && curv < -const__ * gw) {
1113         scale = -gw / curv;
1114         i__1 = *n;
1115         for (i__ = 1; i__ <= i__1; ++i__) {
1116             temp = xopt[i__] + scale * w[i__];
1117 /* L170: */
1118 /* Computing MAX */
1119 /* Computing MIN */
1120             d__3 = su[i__];
1121             d__1 = sl[i__], d__2 = min(d__3,temp);
1122             xalt[i__] = max(d__1,d__2);
1123         }
1124 /* Computing 2nd power */
1125         d__1 = half * gw * scale;
1126         *cauchy = d__1 * d__1;
1127     } else {
1128 /* Computing 2nd power */
1129         d__1 = gw + half * curv;
1130         *cauchy = d__1 * d__1;
1131     }
1132
1133 /*     If IFLAG is zero, then XALT is calculated as before after reversing */
1134 /*     the sign of GLAG. Thus two XALT vectors become available. The one that */
1135 /*     is chosen is the one that gives the larger value of CAUCHY. */
1136
1137     if (iflag == 0) {
1138         i__1 = *n;
1139         for (i__ = 1; i__ <= i__1; ++i__) {
1140             glag[i__] = -glag[i__];
1141 /* L180: */
1142             w[*n + i__] = xalt[i__];
1143         }
1144         csave = *cauchy;
1145         iflag = 1;
1146         goto L100;
1147     }
1148     if (csave > *cauchy) {
1149         i__1 = *n;
1150         for (i__ = 1; i__ <= i__1; ++i__) {
1151 /* L190: */
1152             xalt[i__] = w[*n + i__];
1153         }
1154         *cauchy = csave;
1155     }
1156 L200:
1157     return;
1158 } /* altmov_ */
1159
1160 static void trsbox_(int *n, int *npt, double *xpt,
1161         double *xopt, double *gopt, double *hq, double *pq,
1162         double *sl, double *su, double *delta, double *xnew,
1163         double *d__, double *gnew, double *xbdi, double *s,
1164         double *hs, double *hred, double *dsq, double *crvmin)
1165 {
1166     /* System generated locals */
1167     int xpt_dim1, xpt_offset, i__1, i__2;
1168     double d__1, d__2, d__3, d__4;
1169
1170     /* Local variables */
1171     int i__, j, k, ih;
1172     double ds;
1173     int iu;
1174     double dhd, dhs, cth, one, shs, sth, ssq, half, beta, sdec, blen;
1175     int iact, nact;
1176     double angt, qred;
1177     int isav;
1178     double temp, zero, xsav, xsum, angbd, dredg, sredg;
1179     int iterc;
1180     double resid, delsq, ggsav, tempa, tempb, ratio, sqstp, redmax,
1181             dredsq, redsav, onemin, gredsq, rednew;
1182     int itcsav;
1183     double rdprev, rdnext, stplen, stepsq;
1184     int itermax;
1185
1186
1187 /*     The arguments N, NPT, XPT, XOPT, GOPT, HQ, PQ, SL and SU have the same */
1188 /*       meanings as the corresponding arguments of BOBYQB. */
1189 /*     DELTA is the trust region radius for the present calculation, which */
1190 /*       seeks a small value of the quadratic model within distance DELTA of */
1191 /*       XOPT subject to the bounds on the variables. */
1192 /*     XNEW will be set to a new vector of variables that is approximately */
1193 /*       the one that minimizes the quadratic model within the trust region */
1194 /*       subject to the SL and SU constraints on the variables. It satisfies */
1195 /*       as equations the bounds that become active during the calculation. */
1196 /*     D is the calculated trial step from XOPT, generated iteratively from an */
1197 /*       initial value of zero. Thus XNEW is XOPT+D after the final iteration. */
1198 /*     GNEW holds the gradient of the quadratic model at XOPT+D. It is updated */
1199 /*       when D is updated. */
1200 /*     XBDI is a working space vector. For I=1,2,...,N, the element XBDI(I) is */
1201 /*       set to -1.0, 0.0, or 1.0, the value being nonzero if and only if the */
1202 /*       I-th variable has become fixed at a bound, the bound being SL(I) or */
1203 /*       SU(I) in the case XBDI(I)=-1.0 or XBDI(I)=1.0, respectively. This */
1204 /*       information is accumulated during the construction of XNEW. */
1205 /*     The arrays S, HS and HRED are also used for working space. They hold the */
1206 /*       current search direction, and the changes in the gradient of Q along S */
1207 /*       and the reduced D, respectively, where the reduced D is the same as D, */
1208 /*       except that the components of the fixed variables are zero. */
1209 /*     DSQ will be set to the square of the length of XNEW-XOPT. */
1210 /*     CRVMIN is set to zero if D reaches the trust region boundary. Otherwise */
1211 /*       it is set to the least curvature of H that occurs in the conjugate */
1212 /*       gradient searches that are not restricted by any constraints. The */
1213 /*       value CRVMIN=-1.0D0 is set, however, if all of these searches are */
1214 /*       constrained. */
1215
1216 /*     A version of the truncated conjugate gradient is applied. If a line */
1217 /*     search is restricted by a constraint, then the procedure is restarted, */
1218 /*     the values of the variables that are at their bounds being fixed. If */
1219 /*     the trust region boundary is reached, then further changes may be made */
1220 /*     to D, each one being in the two dimensional space that is spanned */
1221 /*     by the current D and the gradient of Q at XOPT+D, staying on the trust */
1222 /*     region boundary. Termination occurs when the reduction in Q seems to */
1223 /*     be close to the greatest reduction that can be achieved. */
1224
1225 /*     Set some constants. */
1226
1227     /* Parameter adjustments */
1228     xpt_dim1 = *npt;
1229     xpt_offset = 1 + xpt_dim1;
1230     xpt -= xpt_offset;
1231     --xopt;
1232     --gopt;
1233     --hq;
1234     --pq;
1235     --sl;
1236     --su;
1237     --xnew;
1238     --d__;
1239     --gnew;
1240     --xbdi;
1241     --s;
1242     --hs;
1243     --hred;
1244
1245     /* Function Body */
1246     half = .5;
1247     one = 1.;
1248     onemin = -1.;
1249     zero = 0.;
1250
1251 /*     The sign of GOPT(I) gives the sign of the change to the I-th variable */
1252 /*     that will reduce Q from its value at XOPT. Thus XBDI(I) shows whether */
1253 /*     or not to fix the I-th variable at one of its bounds initially, with */
1254 /*     NACT being set to the number of fixed variables. D and GNEW are also */
1255 /*     set for the first iteration. DELSQ is the upper bound on the sum of */
1256 /*     squares of the free variables. QRED is the reduction in Q so far. */
1257
1258     iterc = 0;
1259     nact = 0;
1260     sqstp = zero;
1261     i__1 = *n;
1262     for (i__ = 1; i__ <= i__1; ++i__) {
1263         xbdi[i__] = zero;
1264         if (xopt[i__] <= sl[i__]) {
1265             if (gopt[i__] >= zero) {
1266                 xbdi[i__] = onemin;
1267             }
1268         } else if (xopt[i__] >= su[i__]) {
1269             if (gopt[i__] <= zero) {
1270                 xbdi[i__] = one;
1271             }
1272         }
1273         if (xbdi[i__] != zero) {
1274             ++nact;
1275         }
1276         d__[i__] = zero;
1277 /* L10: */
1278         gnew[i__] = gopt[i__];
1279     }
1280     delsq = *delta * *delta;
1281     qred = zero;
1282     *crvmin = onemin;
1283
1284 /*     Set the next search direction of the conjugate gradient method. It is */
1285 /*     the steepest descent direction initially and when the iterations are */
1286 /*     restarted because a variable has just been fixed by a bound, and of */
1287 /*     course the components of the fixed variables are zero. ITERMAX is an */
1288 /*     upper bound on the indices of the conjugate gradient iterations. */
1289
1290 L20:
1291     beta = zero;
1292 L30:
1293     stepsq = zero;
1294     i__1 = *n;
1295     for (i__ = 1; i__ <= i__1; ++i__) {
1296         if (xbdi[i__] != zero) {
1297             s[i__] = zero;
1298         } else if (beta == zero) {
1299             s[i__] = -gnew[i__];
1300         } else {
1301             s[i__] = beta * s[i__] - gnew[i__];
1302         }
1303 /* L40: */
1304 /* Computing 2nd power */
1305         d__1 = s[i__];
1306         stepsq += d__1 * d__1;
1307     }
1308     if (stepsq == zero) {
1309         goto L190;
1310     }
1311     if (beta == zero) {
1312         gredsq = stepsq;
1313         itermax = iterc + *n - nact;
1314     }
1315     if (gredsq * delsq <= qred * 1e-4 * qred) {
1316         goto L190;
1317     }
1318
1319 /*     Multiply the search direction by the second derivative matrix of Q and */
1320 /*     calculate some scalars for the choice of steplength. Then set BLEN to */
1321 /*     the length of the the step to the trust region boundary and STPLEN to */
1322 /*     the steplength, ignoring the simple bounds. */
1323
1324     goto L210;
1325 L50:
1326     resid = delsq;
1327     ds = zero;
1328     shs = zero;
1329     i__1 = *n;
1330     for (i__ = 1; i__ <= i__1; ++i__) {
1331         if (xbdi[i__] == zero) {
1332 /* Computing 2nd power */
1333             d__1 = d__[i__];
1334             resid -= d__1 * d__1;
1335             ds += s[i__] * d__[i__];
1336             shs += s[i__] * hs[i__];
1337         }
1338 /* L60: */
1339     }
1340     if (resid <= zero) {
1341         goto L90;
1342     }
1343     temp = sqrt(stepsq * resid + ds * ds);
1344     if (ds < zero) {
1345         blen = (temp - ds) / stepsq;
1346     } else {
1347         blen = resid / (temp + ds);
1348     }
1349     stplen = blen;
1350     if (shs > zero) {
1351 /* Computing MIN */
1352         d__1 = blen, d__2 = gredsq / shs;
1353         stplen = min(d__1,d__2);
1354     }
1355
1356 /*     Reduce STPLEN if necessary in order to preserve the simple bounds, */
1357 /*     letting IACT be the index of the new constrained variable. */
1358
1359     iact = 0;
1360     i__1 = *n;
1361     for (i__ = 1; i__ <= i__1; ++i__) {
1362         if (s[i__] != zero) {
1363             xsum = xopt[i__] + d__[i__];
1364             if (s[i__] > zero) {
1365                 temp = (su[i__] - xsum) / s[i__];
1366             } else {
1367                 temp = (sl[i__] - xsum) / s[i__];
1368             }
1369             if (temp < stplen) {
1370                 stplen = temp;
1371                 iact = i__;
1372             }
1373         }
1374 /* L70: */
1375     }
1376
1377 /*     Update CRVMIN, GNEW and D. Set SDEC to the decrease that occurs in Q. */
1378
1379     sdec = zero;
1380     if (stplen > zero) {
1381         ++iterc;
1382         temp = shs / stepsq;
1383         if (iact == 0 && temp > zero) {
1384             *crvmin = min(*crvmin,temp);
1385             if (*crvmin == onemin) {
1386                 *crvmin = temp;
1387             }
1388         }
1389         ggsav = gredsq;
1390         gredsq = zero;
1391         i__1 = *n;
1392         for (i__ = 1; i__ <= i__1; ++i__) {
1393             gnew[i__] += stplen * hs[i__];
1394             if (xbdi[i__] == zero) {
1395 /* Computing 2nd power */
1396                 d__1 = gnew[i__];
1397                 gredsq += d__1 * d__1;
1398             }
1399 /* L80: */
1400             d__[i__] += stplen * s[i__];
1401         }
1402 /* Computing MAX */
1403         d__1 = stplen * (ggsav - half * stplen * shs);
1404         sdec = max(d__1,zero);
1405         qred += sdec;
1406     }
1407
1408 /*     Restart the conjugate gradient method if it has hit a new bound. */
1409
1410     if (iact > 0) {
1411         ++nact;
1412         xbdi[iact] = one;
1413         if (s[iact] < zero) {
1414             xbdi[iact] = onemin;
1415         }
1416 /* Computing 2nd power */
1417         d__1 = d__[iact];
1418         delsq -= d__1 * d__1;
1419         if (delsq <= zero) {
1420             goto L90;
1421         }
1422         goto L20;
1423     }
1424
1425 /*     If STPLEN is less than BLEN, then either apply another conjugate */
1426 /*     gradient iteration or RETURN. */
1427
1428     if (stplen < blen) {
1429         if (iterc == itermax) {
1430             goto L190;
1431         }
1432         if (sdec <= qred * .01) {
1433             goto L190;
1434         }
1435         beta = gredsq / ggsav;
1436         goto L30;
1437     }
1438 L90:
1439     *crvmin = zero;
1440
1441 /*     Prepare for the alternative iteration by calculating some scalars */
1442 /*     and by multiplying the reduced D by the second derivative matrix of */
1443 /*     Q, where S holds the reduced D in the call of GGMULT. */
1444
1445 L100:
1446     if (nact >= *n - 1) {
1447         goto L190;
1448     }
1449     dredsq = zero;
1450     dredg = zero;
1451     gredsq = zero;
1452     i__1 = *n;
1453     for (i__ = 1; i__ <= i__1; ++i__) {
1454         if (xbdi[i__] == zero) {
1455 /* Computing 2nd power */
1456             d__1 = d__[i__];
1457             dredsq += d__1 * d__1;
1458             dredg += d__[i__] * gnew[i__];
1459 /* Computing 2nd power */
1460             d__1 = gnew[i__];
1461             gredsq += d__1 * d__1;
1462             s[i__] = d__[i__];
1463         } else {
1464             s[i__] = zero;
1465         }
1466 /* L110: */
1467     }
1468     itcsav = iterc;
1469     goto L210;
1470
1471 /*     Let the search direction S be a linear combination of the reduced D */
1472 /*     and the reduced G that is orthogonal to the reduced D. */
1473
1474 L120:
1475     ++iterc;
1476     temp = gredsq * dredsq - dredg * dredg;
1477     if (temp <= qred * 1e-4 * qred) {
1478         goto L190;
1479     }
1480     temp = sqrt(temp);
1481     i__1 = *n;
1482     for (i__ = 1; i__ <= i__1; ++i__) {
1483         if (xbdi[i__] == zero) {
1484             s[i__] = (dredg * d__[i__] - dredsq * gnew[i__]) / temp;
1485         } else {
1486             s[i__] = zero;
1487         }
1488 /* L130: */
1489     }
1490     sredg = -temp;
1491
1492 /*     By considering the simple bounds on the variables, calculate an upper */
1493 /*     bound on the tangent of half the angle of the alternative iteration, */
1494 /*     namely ANGBD, except that, if already a free variable has reached a */
1495 /*     bound, there is a branch back to label 100 after fixing that variable. */
1496
1497     angbd = one;
1498     iact = 0;
1499     i__1 = *n;
1500     for (i__ = 1; i__ <= i__1; ++i__) {
1501         if (xbdi[i__] == zero) {
1502             tempa = xopt[i__] + d__[i__] - sl[i__];
1503             tempb = su[i__] - xopt[i__] - d__[i__];
1504             if (tempa <= zero) {
1505                 ++nact;
1506                 xbdi[i__] = onemin;
1507                 goto L100;
1508             } else if (tempb <= zero) {
1509                 ++nact;
1510                 xbdi[i__] = one;
1511                 goto L100;
1512             }
1513             ratio = one;
1514 /* Computing 2nd power */
1515             d__1 = d__[i__];
1516 /* Computing 2nd power */
1517             d__2 = s[i__];
1518             ssq = d__1 * d__1 + d__2 * d__2;
1519 /* Computing 2nd power */
1520             d__1 = xopt[i__] - sl[i__];
1521             temp = ssq - d__1 * d__1;
1522             if (temp > zero) {
1523                 temp = sqrt(temp) - s[i__];
1524                 if (angbd * temp > tempa) {
1525                     angbd = tempa / temp;
1526                     iact = i__;
1527                     xsav = onemin;
1528                 }
1529             }
1530 /* Computing 2nd power */
1531             d__1 = su[i__] - xopt[i__];
1532             temp = ssq - d__1 * d__1;
1533             if (temp > zero) {
1534                 temp = sqrt(temp) + s[i__];
1535                 if (angbd * temp > tempb) {
1536                     angbd = tempb / temp;
1537                     iact = i__;
1538                     xsav = one;
1539                 }
1540             }
1541         }
1542 /* L140: */
1543     }
1544
1545 /*     Calculate HHD and some curvatures for the alternative iteration. */
1546
1547     goto L210;
1548 L150:
1549     shs = zero;
1550     dhs = zero;
1551     dhd = zero;
1552     i__1 = *n;
1553     for (i__ = 1; i__ <= i__1; ++i__) {
1554         if (xbdi[i__] == zero) {
1555             shs += s[i__] * hs[i__];
1556             dhs += d__[i__] * hs[i__];
1557             dhd += d__[i__] * hred[i__];
1558         }
1559 /* L160: */
1560     }
1561
1562 /*     Seek the greatest reduction in Q for a range of equally spaced values */
1563 /*     of ANGT in [0,ANGBD], where ANGT is the tangent of half the angle of */
1564 /*     the alternative iteration. */
1565
1566     redmax = zero;
1567     isav = 0;
1568     redsav = zero;
1569     iu = (int) (angbd * 17. + 3.1);
1570     i__1 = iu;
1571     for (i__ = 1; i__ <= i__1; ++i__) {
1572         angt = angbd * (double) i__ / (double) iu;
1573         sth = (angt + angt) / (one + angt * angt);
1574         temp = shs + angt * (angt * dhd - dhs - dhs);
1575         rednew = sth * (angt * dredg - sredg - half * sth * temp);
1576         if (rednew > redmax) {
1577             redmax = rednew;
1578             isav = i__;
1579             rdprev = redsav;
1580         } else if (i__ == isav + 1) {
1581             rdnext = rednew;
1582         }
1583 /* L170: */
1584         redsav = rednew;
1585     }
1586
1587 /*     Return if the reduction is zero. Otherwise, set the sine and cosine */
1588 /*     of the angle of the alternative iteration, and calculate SDEC. */
1589
1590     if (isav == 0) {
1591         goto L190;
1592     }
1593     if (isav < iu) {
1594         temp = (rdnext - rdprev) / (redmax + redmax - rdprev - rdnext);
1595         angt = angbd * ((double) isav + half * temp) / (double) iu;
1596     }
1597     cth = (one - angt * angt) / (one + angt * angt);
1598     sth = (angt + angt) / (one + angt * angt);
1599     temp = shs + angt * (angt * dhd - dhs - dhs);
1600     sdec = sth * (angt * dredg - sredg - half * sth * temp);
1601     if (sdec <= zero) {
1602         goto L190;
1603     }
1604
1605 /*     Update GNEW, D and HRED. If the angle of the alternative iteration */
1606 /*     is restricted by a bound on a free variable, that variable is fixed */
1607 /*     at the bound. */
1608
1609     dredg = zero;
1610     gredsq = zero;
1611     i__1 = *n;
1612     for (i__ = 1; i__ <= i__1; ++i__) {
1613         gnew[i__] = gnew[i__] + (cth - one) * hred[i__] + sth * hs[i__];
1614         if (xbdi[i__] == zero) {
1615             d__[i__] = cth * d__[i__] + sth * s[i__];
1616             dredg += d__[i__] * gnew[i__];
1617 /* Computing 2nd power */
1618             d__1 = gnew[i__];
1619             gredsq += d__1 * d__1;
1620         }
1621 /* L180: */
1622         hred[i__] = cth * hred[i__] + sth * hs[i__];
1623     }
1624     qred += sdec;
1625     if (iact > 0 && isav == iu) {
1626         ++nact;
1627         xbdi[iact] = xsav;
1628         goto L100;
1629     }
1630
1631 /*     If SDEC is sufficiently small, then RETURN after setting XNEW to */
1632 /*     XOPT+D, giving careful attention to the bounds. */
1633
1634     if (sdec > qred * .01) {
1635         goto L120;
1636     }
1637 L190:
1638     *dsq = zero;
1639     i__1 = *n;
1640     for (i__ = 1; i__ <= i__1; ++i__) {
1641 /* Computing MAX */
1642 /* Computing MIN */
1643         d__3 = xopt[i__] + d__[i__], d__4 = su[i__];
1644         d__1 = min(d__3,d__4), d__2 = sl[i__];
1645         xnew[i__] = max(d__1,d__2);
1646         if (xbdi[i__] == onemin) {
1647             xnew[i__] = sl[i__];
1648         }
1649         if (xbdi[i__] == one) {
1650             xnew[i__] = su[i__];
1651         }
1652         d__[i__] = xnew[i__] - xopt[i__];
1653 /* L200: */
1654 /* Computing 2nd power */
1655         d__1 = d__[i__];
1656         *dsq += d__1 * d__1;
1657     }
1658     return;
1659 /*     The following instructions multiply the current S-vector by the second */
1660 /*     derivative matrix of the quadratic model, putting the product in HS. */
1661 /*     They are reached from three different parts of the software above and */
1662 /*     they can be regarded as an external subroutine. */
1663
1664 L210:
1665     ih = 0;
1666     i__1 = *n;
1667     for (j = 1; j <= i__1; ++j) {
1668         hs[j] = zero;
1669         i__2 = j;
1670         for (i__ = 1; i__ <= i__2; ++i__) {
1671             ++ih;
1672             if (i__ < j) {
1673                 hs[j] += hq[ih] * s[i__];
1674             }
1675 /* L220: */
1676             hs[i__] += hq[ih] * s[j];
1677         }
1678     }
1679     i__2 = *npt;
1680     for (k = 1; k <= i__2; ++k) {
1681         if (pq[k] != zero) {
1682             temp = zero;
1683             i__1 = *n;
1684             for (j = 1; j <= i__1; ++j) {
1685 /* L230: */
1686                 temp += xpt[k + j * xpt_dim1] * s[j];
1687             }
1688             temp *= pq[k];
1689             i__1 = *n;
1690             for (i__ = 1; i__ <= i__1; ++i__) {
1691 /* L240: */
1692                 hs[i__] += temp * xpt[k + i__ * xpt_dim1];
1693             }
1694         }
1695 /* L250: */
1696     }
1697     if (*crvmin != zero) {
1698         goto L50;
1699     }
1700     if (iterc > itcsav) {
1701         goto L150;
1702     }
1703     i__2 = *n;
1704     for (i__ = 1; i__ <= i__2; ++i__) {
1705 /* L260: */
1706         hred[i__] = hs[i__];
1707     }
1708     goto L120;
1709 } /* trsbox_ */
1710
1711 static nlopt_result prelim_(int *n, int *npt, double *x,
1712         const double *xl, const double *xu, double *rhobeg,
1713                     nlopt_stopping *stop,
1714                     bobyqa_func calfun, void *calfun_data,
1715          double *xbase, double *xpt, double *fval,
1716          double *gopt, double *hq, double *pq, double *bmat,
1717         double *zmat, int *ndim, double *sl, double *su,
1718                     int *kopt)
1719 {
1720     /* System generated locals */
1721     int xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1,
1722             zmat_offset, i__1, i__2;
1723     double d__1, d__2, d__3, d__4;
1724
1725     /* Local variables */
1726     double f;
1727     int i__, j, k, ih, np, nfm;
1728     double one;
1729     int nfx, ipt, jpt;
1730     double two, fbeg, diff, half, temp, zero, recip, stepa, stepb;
1731     int itemp;
1732     double rhosq;
1733
1734     int nf;
1735
1736 /*     The arguments N, NPT, X, XL, XU, RHOBEG, and MAXFUN are the */
1737 /*       same as the corresponding arguments in SUBROUTINE BOBYQA. */
1738 /*     The arguments XBASE, XPT, FVAL, HQ, PQ, BMAT, ZMAT, NDIM, SL and SU */
1739 /*       are the same as the corresponding arguments in BOBYQB, the elements */
1740 /*       of SL and SU being set in BOBYQA. */
1741 /*     GOPT is usually the gradient of the quadratic model at XOPT+XBASE, but */
1742 /*       it is set by PRELIM to the gradient of the quadratic model at XBASE. */
1743 /*       If XOPT is nonzero, BOBYQB will change it to its usual value later. */
1744 /*     NF is maintaned as the number of calls of CALFUN so far. */
1745 /*     KOPT will be such that the least calculated value of F so far is at */
1746 /*       the point XPT(KOPT,.)+XBASE in the space of the variables. */
1747
1748 /*     SUBROUTINE PRELIM sets the elements of XBASE, XPT, FVAL, GOPT, HQ, PQ, */
1749 /*     BMAT and ZMAT for the first iteration, and it maintains the values of */
1750 /*     NF and KOPT. The vector X is also changed by PRELIM. */
1751
1752 /*     Set some constants. */
1753
1754     /* Parameter adjustments */
1755     zmat_dim1 = *npt;
1756     zmat_offset = 1 + zmat_dim1;
1757     zmat -= zmat_offset;
1758     xpt_dim1 = *npt;
1759     xpt_offset = 1 + xpt_dim1;
1760     xpt -= xpt_offset;
1761     --x;
1762     --xl;
1763     --xu;
1764     --xbase;
1765     --fval;
1766     --gopt;
1767     --hq;
1768     --pq;
1769     bmat_dim1 = *ndim;
1770     bmat_offset = 1 + bmat_dim1;
1771     bmat -= bmat_offset;
1772     --sl;
1773     --su;
1774
1775     /* Function Body */
1776     half = .5;
1777     one = 1.;
1778     two = 2.;
1779     zero = 0.;
1780     rhosq = *rhobeg * *rhobeg;
1781     recip = one / rhosq;
1782     np = *n + 1;
1783
1784 /*     Set XBASE to the initial vector of variables, and set the initial */
1785 /*     elements of XPT, BMAT, HQ, PQ and ZMAT to zero. */
1786
1787     i__1 = *n;
1788     for (j = 1; j <= i__1; ++j) {
1789         xbase[j] = x[j];
1790         i__2 = *npt;
1791         for (k = 1; k <= i__2; ++k) {
1792 /* L10: */
1793             xpt[k + j * xpt_dim1] = zero;
1794         }
1795         i__2 = *ndim;
1796         for (i__ = 1; i__ <= i__2; ++i__) {
1797 /* L20: */
1798             bmat[i__ + j * bmat_dim1] = zero;
1799         }
1800     }
1801     i__2 = *n * np / 2;
1802     for (ih = 1; ih <= i__2; ++ih) {
1803 /* L30: */
1804         hq[ih] = zero;
1805     }
1806     i__2 = *npt;
1807     for (k = 1; k <= i__2; ++k) {
1808         pq[k] = zero;
1809         i__1 = *npt - np;
1810         for (j = 1; j <= i__1; ++j) {
1811 /* L40: */
1812             zmat[k + j * zmat_dim1] = zero;
1813         }
1814     }
1815
1816 /*     Begin the initialization procedure. NF becomes one more than the number */
1817 /*     of function values so far. The coordinates of the displacement of the */
1818 /*     next initial interpolation point from XBASE are set in XPT(NF+1,.). */
1819
1820     nf = 0;
1821 L50:
1822     nfm = nf;
1823     nfx = nf - *n;
1824     ++(nf);
1825     if (nfm <= *n << 1) {
1826         if (nfm >= 1 && nfm <= *n) {
1827             stepa = *rhobeg;
1828             if (su[nfm] == zero) {
1829                 stepa = -stepa;
1830             }
1831             xpt[nf + nfm * xpt_dim1] = stepa;
1832         } else if (nfm > *n) {
1833             stepa = xpt[nf - *n + nfx * xpt_dim1];
1834             stepb = -(*rhobeg);
1835             if (sl[nfx] == zero) {
1836 /* Computing MIN */
1837                 d__1 = two * *rhobeg, d__2 = su[nfx];
1838                 stepb = min(d__1,d__2);
1839             }
1840             if (su[nfx] == zero) {
1841 /* Computing MAX */
1842                 d__1 = -two * *rhobeg, d__2 = sl[nfx];
1843                 stepb = max(d__1,d__2);
1844             }
1845             xpt[nf + nfx * xpt_dim1] = stepb;
1846         }
1847     } else {
1848         itemp = (nfm - np) / *n;
1849         jpt = nfm - itemp * *n - *n;
1850         ipt = jpt + itemp;
1851         if (ipt > *n) {
1852             itemp = jpt;
1853             jpt = ipt - *n;
1854             ipt = itemp;
1855         }
1856         xpt[nf + ipt * xpt_dim1] = xpt[ipt + 1 + ipt * xpt_dim1];
1857         xpt[nf + jpt * xpt_dim1] = xpt[jpt + 1 + jpt * xpt_dim1];
1858     }
1859
1860 /*     Calculate the next value of F. The least function value so far and */
1861 /*     its index are required. */
1862
1863     i__1 = *n;
1864     for (j = 1; j <= i__1; ++j) {
1865 /* Computing MIN */
1866 /* Computing MAX */
1867         d__3 = xl[j], d__4 = xbase[j] + xpt[nf + j * xpt_dim1];
1868         d__1 = max(d__3,d__4), d__2 = xu[j];
1869         x[j] = min(d__1,d__2);
1870         if (xpt[nf + j * xpt_dim1] == sl[j]) {
1871             x[j] = xl[j];
1872         }
1873         if (xpt[nf + j * xpt_dim1] == su[j]) {
1874             x[j] = xu[j];
1875         }
1876 /* L60: */
1877     }
1878     stop->nevals++;
1879     f = calfun(*n, &x[1], calfun_data);
1880     fval[nf] = f;
1881     if (nf == 1) {
1882         fbeg = f;
1883         *kopt = 1;
1884     } else if (f < fval[*kopt]) {
1885         *kopt = nf;
1886     }
1887
1888 /*     Set the nonzero initial elements of BMAT and the quadratic model in the */
1889 /*     cases when NF is at most 2*N+1. If NF exceeds N+1, then the positions */
1890 /*     of the NF-th and (NF-N)-th interpolation points may be switched, in */
1891 /*     order that the function value at the first of them contributes to the */
1892 /*     off-diagonal second derivative terms of the initial quadratic model. */
1893
1894     if (nf <= (*n << 1) + 1) {
1895         if (nf >= 2 && nf <= *n + 1) {
1896             gopt[nfm] = (f - fbeg) / stepa;
1897             if (*npt < nf + *n) {
1898                 bmat[nfm * bmat_dim1 + 1] = -one / stepa;
1899                 bmat[nf + nfm * bmat_dim1] = one / stepa;
1900                 bmat[*npt + nfm + nfm * bmat_dim1] = -half * rhosq;
1901             }
1902         } else if (nf >= *n + 2) {
1903             ih = nfx * (nfx + 1) / 2;
1904             temp = (f - fbeg) / stepb;
1905             diff = stepb - stepa;
1906             hq[ih] = two * (temp - gopt[nfx]) / diff;
1907             gopt[nfx] = (gopt[nfx] * stepb - temp * stepa) / diff;
1908             if (stepa * stepb < zero) {
1909                 if (f < fval[nf - *n]) {
1910                     fval[nf] = fval[nf - *n];
1911                     fval[nf - *n] = f;
1912                     if (*kopt == nf) {
1913                         *kopt = nf - *n;
1914                     }
1915                     xpt[nf - *n + nfx * xpt_dim1] = stepb;
1916                     xpt[nf + nfx * xpt_dim1] = stepa;
1917                 }
1918             }
1919             bmat[nfx * bmat_dim1 + 1] = -(stepa + stepb) / (stepa * stepb);
1920             bmat[nf + nfx * bmat_dim1] = -half / xpt[nf - *n + nfx *
1921                     xpt_dim1];
1922             bmat[nf - *n + nfx * bmat_dim1] = -bmat[nfx * bmat_dim1 + 1] -
1923                     bmat[nf + nfx * bmat_dim1];
1924             zmat[nfx * zmat_dim1 + 1] = sqrt(two) / (stepa * stepb);
1925             zmat[nf + nfx * zmat_dim1] = sqrt(half) / rhosq;
1926             zmat[nf - *n + nfx * zmat_dim1] = -zmat[nfx * zmat_dim1 + 1] -
1927                     zmat[nf + nfx * zmat_dim1];
1928         }
1929
1930 /*     Set the off-diagonal second derivatives of the Lagrange functions and */
1931 /*     the initial quadratic model. */
1932
1933     } else {
1934         ih = ipt * (ipt - 1) / 2 + jpt;
1935         zmat[nfx * zmat_dim1 + 1] = recip;
1936         zmat[nf + nfx * zmat_dim1] = recip;
1937         zmat[ipt + 1 + nfx * zmat_dim1] = -recip;
1938         zmat[jpt + 1 + nfx * zmat_dim1] = -recip;
1939         temp = xpt[nf + ipt * xpt_dim1] * xpt[nf + jpt * xpt_dim1];
1940         hq[ih] = (fbeg - fval[ipt + 1] - fval[jpt + 1] + f) / temp;
1941     }
1942     if (nlopt_stop_forced(stop)) return NLOPT_FORCE_STOP;
1943     else if (f < stop->minf_max) return NLOPT_MINF_MAX_REACHED;
1944     else if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
1945     else if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
1946     if (nf < *npt) {
1947         goto L50;
1948     }
1949     return NLOPT_SUCCESS;
1950 } /* prelim_ */
1951
1952 static nlopt_result bobyqb_(int *n, int *npt, double *x,
1953         const double *xl, const double *xu, double *rhobeg, double *
1954         rhoend,
1955                             nlopt_stopping *stop,
1956                             bobyqa_func calfun, void *calfun_data,
1957                             double *minf,
1958         double *xbase,
1959         double *xpt, double *fval, double *xopt, double *gopt,
1960          double *hq, double *pq, double *bmat, double *zmat,
1961         int *ndim, double *sl, double *su, double *xnew,
1962         double *xalt, double *d__, double *vlag, double *w)
1963 {
1964     /* System generated locals */
1965     int xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1,
1966             zmat_offset, i__1, i__2, i__3;
1967     double d__1, d__2, d__3, d__4;
1968
1969     /* Local variables */
1970     double f;
1971     int i__, j, k, ih, jj, nh, ip, jp;
1972     double dx;
1973     int np;
1974     double den, one, ten, dsq, rho, sum, two, diff, half, beta, gisq;
1975     int knew;
1976     double temp, suma, sumb, bsum, fopt;
1977     int kopt, nptm;
1978     double zero, curv;
1979     int ksav;
1980     double gqsq, dist, sumw, sumz, diffa, diffb, diffc, hdiag;
1981     int kbase;
1982     double alpha, delta, adelt, denom, fsave, bdtol, delsq;
1983     int nresc, nfsav;
1984     double ratio, dnorm, vquad, pqold, tenth;
1985     int itest;
1986     double sumpq, scaden;
1987     double errbig, cauchy, fracsq, biglsq, densav;
1988     double bdtest;
1989     double crvmin, frhosq;
1990     double distsq;
1991     int ntrits;
1992     double xoptsq;
1993
1994     nlopt_result rc = NLOPT_SUCCESS, rc2;
1995
1996 /*     The arguments N, NPT, X, XL, XU, RHOBEG, RHOEND, and MAXFUN */
1997 /*       are identical to the corresponding arguments in SUBROUTINE BOBYQA. */
1998 /*     XBASE holds a shift of origin that should reduce the contributions */
1999 /*       from rounding errors to values of the model and Lagrange functions. */
2000 /*     XPT is a two-dimensional array that holds the coordinates of the */
2001 /*       interpolation points relative to XBASE. */
2002 /*     FVAL holds the values of F at the interpolation points. */
2003 /*     XOPT is set to the displacement from XBASE of the trust region centre. */
2004 /*     GOPT holds the gradient of the quadratic model at XBASE+XOPT. */
2005 /*     HQ holds the explicit second derivatives of the quadratic model. */
2006 /*     PQ contains the parameters of the implicit second derivatives of the */
2007 /*       quadratic model. */
2008 /*     BMAT holds the last N columns of H. */
2009 /*     ZMAT holds the factorization of the leading NPT by NPT submatrix of H, */
2010 /*       this factorization being ZMAT times ZMAT^T, which provides both the */
2011 /*       correct rank and positive semi-definiteness. */
2012 /*     NDIM is the first dimension of BMAT and has the value NPT+N. */
2013 /*     SL and SU hold the differences XL-XBASE and XU-XBASE, respectively. */
2014 /*       All the components of every XOPT are going to satisfy the bounds */
2015 /*       SL(I) .LEQ. XOPT(I) .LEQ. SU(I), with appropriate equalities when */
2016 /*       XOPT is on a constraint boundary. */
2017 /*     XNEW is chosen by SUBROUTINE TRSBOX or ALTMOV. Usually XBASE+XNEW is the */
2018 /*       vector of variables for the next call of CALFUN. XNEW also satisfies */
2019 /*       the SL and SU constraints in the way that has just been mentioned. */
2020 /*     XALT is an alternative to XNEW, chosen by ALTMOV, that may replace XNEW */
2021 /*       in order to increase the denominator in the updating of UPDATE. */
2022 /*     D is reserved for a trial step from XOPT, which is usually XNEW-XOPT. */
2023 /*     VLAG contains the values of the Lagrange functions at a new point X. */
2024 /*       They are part of a product that requires VLAG to be of length NDIM. */
2025 /*     W is a one-dimensional array that is used for working space. Its length */
2026 /*       must be at least 3*NDIM = 3*(NPT+N). */
2027
2028 /*     Set some constants. */
2029
2030     /* Parameter adjustments */
2031     zmat_dim1 = *npt;
2032     zmat_offset = 1 + zmat_dim1;
2033     zmat -= zmat_offset;
2034     xpt_dim1 = *npt;
2035     xpt_offset = 1 + xpt_dim1;
2036     xpt -= xpt_offset;
2037     --x;
2038     --xl;
2039     --xu;
2040     --xbase;
2041     --fval;
2042     --xopt;
2043     --gopt;
2044     --hq;
2045     --pq;
2046     bmat_dim1 = *ndim;
2047     bmat_offset = 1 + bmat_dim1;
2048     bmat -= bmat_offset;
2049     --sl;
2050     --su;
2051     --xnew;
2052     --xalt;
2053     --d__;
2054     --vlag;
2055     --w;
2056
2057     /* Function Body */
2058     half = .5;
2059     one = 1.;
2060     ten = 10.;
2061     tenth = .1;
2062     two = 2.;
2063     zero = 0.;
2064     np = *n + 1;
2065     nptm = *npt - np;
2066     nh = *n * np / 2;
2067
2068 /*     The call of PRELIM sets the elements of XBASE, XPT, FVAL, GOPT, HQ, PQ, */
2069 /*     BMAT and ZMAT for the first iteration, with the corresponding values of */
2070 /*     of NF and KOPT, which are the number of calls of CALFUN so far and the */
2071 /*     index of the interpolation point at the trust region centre. Then the */
2072 /*     initial XOPT is set too. The branch to label 720 occurs if MAXFUN is */
2073 /*     less than NPT. GOPT will be updated if KOPT is different from KBASE. */
2074
2075     rc2 = prelim_(n, npt, &x[1], &xl[1], &xu[1], rhobeg,
2076                   stop, calfun, calfun_data,
2077             &xbase[1], &xpt[xpt_offset], &fval[1], &gopt[1], &hq[1], &pq[1], &bmat[
2078             bmat_offset], &zmat[zmat_offset], ndim, &sl[1], &su[1], &kopt);
2079     xoptsq = zero;
2080     i__1 = *n;
2081     for (i__ = 1; i__ <= i__1; ++i__) {
2082         xopt[i__] = xpt[kopt + i__ * xpt_dim1];
2083 /* L10: */
2084 /* Computing 2nd power */
2085         d__1 = xopt[i__];
2086         xoptsq += d__1 * d__1;
2087     }
2088     fsave = fval[1];
2089     if (rc2 != NLOPT_SUCCESS) {
2090       rc = rc2;
2091       goto L720;
2092     }
2093     kbase = 1;
2094
2095 /*     Complete the settings that are required for the iterative procedure. */
2096
2097     rho = *rhobeg;
2098     delta = rho;
2099     nresc = stop->nevals;
2100     ntrits = 0;
2101     diffa = zero;
2102     diffb = zero;
2103     itest = 0;
2104     nfsav = stop->nevals;
2105
2106 /*     Update GOPT if necessary before the first iteration and after each */
2107 /*     call of RESCUE that makes a call of CALFUN. */
2108
2109 L20:
2110     if (kopt != kbase) {
2111         ih = 0;
2112         i__1 = *n;
2113         for (j = 1; j <= i__1; ++j) {
2114             i__2 = j;
2115             for (i__ = 1; i__ <= i__2; ++i__) {
2116                 ++ih;
2117                 if (i__ < j) {
2118                     gopt[j] += hq[ih] * xopt[i__];
2119                 }
2120 /* L30: */
2121                 gopt[i__] += hq[ih] * xopt[j];
2122             }
2123         }
2124         if (stop->nevals > *npt) {
2125             i__2 = *npt;
2126             for (k = 1; k <= i__2; ++k) {
2127                 temp = zero;
2128                 i__1 = *n;
2129                 for (j = 1; j <= i__1; ++j) {
2130 /* L40: */
2131                     temp += xpt[k + j * xpt_dim1] * xopt[j];
2132                 }
2133                 temp = pq[k] * temp;
2134                 i__1 = *n;
2135                 for (i__ = 1; i__ <= i__1; ++i__) {
2136 /* L50: */
2137                     gopt[i__] += temp * xpt[k + i__ * xpt_dim1];
2138                 }
2139             }
2140         }
2141     }
2142
2143 /*     Generate the next point in the trust region that provides a small value */
2144 /*     of the quadratic model subject to the constraints on the variables. */
2145 /*     The int NTRITS is set to the number "trust region" iterations that */
2146 /*     have occurred since the last "alternative" iteration. If the length */
2147 /*     of XNEW-XOPT is less than HALF*RHO, however, then there is a branch to */
2148 /*     label 650 or 680 with NTRITS=-1, instead of calculating F at XNEW. */
2149
2150 L60:
2151     trsbox_(n, npt, &xpt[xpt_offset], &xopt[1], &gopt[1], &hq[1], &pq[1], &sl[
2152             1], &su[1], &delta, &xnew[1], &d__[1], &w[1], &w[np], &w[np + *n],
2153              &w[np + (*n << 1)], &w[np + *n * 3], &dsq, &crvmin);
2154 /* Computing MIN */
2155     d__1 = delta, d__2 = sqrt(dsq);
2156     dnorm = min(d__1,d__2);
2157     if (dnorm < half * rho) {
2158         ntrits = -1;
2159 /* Computing 2nd power */
2160         d__1 = ten * rho;
2161         distsq = d__1 * d__1;
2162         if (stop->nevals <= nfsav + 2) {
2163             goto L650;
2164         }
2165
2166 /*     The following choice between labels 650 and 680 depends on whether or */
2167 /*     not our work with the current RHO seems to be complete. Either RHO is */
2168 /*     decreased or termination occurs if the errors in the quadratic model at */
2169 /*     the last three interpolation points compare favourably with predictions */
2170 /*     of likely improvements to the model within distance HALF*RHO of XOPT. */
2171
2172 /* Computing MAX */
2173         d__1 = max(diffa,diffb);
2174         errbig = max(d__1,diffc);
2175         frhosq = rho * .125 * rho;
2176         if (crvmin > zero && errbig > frhosq * crvmin) {
2177             goto L650;
2178         }
2179         bdtol = errbig / rho;
2180         i__1 = *n;
2181         for (j = 1; j <= i__1; ++j) {
2182             bdtest = bdtol;
2183             if (xnew[j] == sl[j]) {
2184                 bdtest = w[j];
2185             }
2186             if (xnew[j] == su[j]) {
2187                 bdtest = -w[j];
2188             }
2189             if (bdtest < bdtol) {
2190                 curv = hq[(j + j * j) / 2];
2191                 i__2 = *npt;
2192                 for (k = 1; k <= i__2; ++k) {
2193 /* L70: */
2194 /* Computing 2nd power */
2195                     d__1 = xpt[k + j * xpt_dim1];
2196                     curv += pq[k] * (d__1 * d__1);
2197                 }
2198                 bdtest += half * curv * rho;
2199                 if (bdtest < bdtol) {
2200                     goto L650;
2201                 }
2202             }
2203 /* L80: */
2204         }
2205         goto L680;
2206     }
2207     ++ntrits;
2208
2209 /*     Severe cancellation is likely to occur if XOPT is too far from XBASE. */
2210 /*     If the following test holds, then XBASE is shifted so that XOPT becomes */
2211 /*     zero. The appropriate changes are made to BMAT and to the second */
2212 /*     derivatives of the current model, beginning with the changes to BMAT */
2213 /*     that do not depend on ZMAT. VLAG is used temporarily for working space. */
2214
2215 L90:
2216     if (dsq <= xoptsq * .001) {
2217         fracsq = xoptsq * .25;
2218         sumpq = zero;
2219         i__1 = *npt;
2220         for (k = 1; k <= i__1; ++k) {
2221             sumpq += pq[k];
2222             sum = -half * xoptsq;
2223             i__2 = *n;
2224             for (i__ = 1; i__ <= i__2; ++i__) {
2225 /* L100: */
2226                 sum += xpt[k + i__ * xpt_dim1] * xopt[i__];
2227             }
2228             w[*npt + k] = sum;
2229             temp = fracsq - half * sum;
2230             i__2 = *n;
2231             for (i__ = 1; i__ <= i__2; ++i__) {
2232                 w[i__] = bmat[k + i__ * bmat_dim1];
2233                 vlag[i__] = sum * xpt[k + i__ * xpt_dim1] + temp * xopt[i__];
2234                 ip = *npt + i__;
2235                 i__3 = i__;
2236                 for (j = 1; j <= i__3; ++j) {
2237 /* L110: */
2238                     bmat[ip + j * bmat_dim1] = bmat[ip + j * bmat_dim1] + w[
2239                             i__] * vlag[j] + vlag[i__] * w[j];
2240                 }
2241             }
2242         }
2243
2244 /*     Then the revisions of BMAT that depend on ZMAT are calculated. */
2245
2246         i__3 = nptm;
2247         for (jj = 1; jj <= i__3; ++jj) {
2248             sumz = zero;
2249             sumw = zero;
2250             i__2 = *npt;
2251             for (k = 1; k <= i__2; ++k) {
2252                 sumz += zmat[k + jj * zmat_dim1];
2253                 vlag[k] = w[*npt + k] * zmat[k + jj * zmat_dim1];
2254 /* L120: */
2255                 sumw += vlag[k];
2256             }
2257             i__2 = *n;
2258             for (j = 1; j <= i__2; ++j) {
2259                 sum = (fracsq * sumz - half * sumw) * xopt[j];
2260                 i__1 = *npt;
2261                 for (k = 1; k <= i__1; ++k) {
2262 /* L130: */
2263                     sum += vlag[k] * xpt[k + j * xpt_dim1];
2264                 }
2265                 w[j] = sum;
2266                 i__1 = *npt;
2267                 for (k = 1; k <= i__1; ++k) {
2268 /* L140: */
2269                     bmat[k + j * bmat_dim1] += sum * zmat[k + jj * zmat_dim1];
2270                 }
2271             }
2272             i__1 = *n;
2273             for (i__ = 1; i__ <= i__1; ++i__) {
2274                 ip = i__ + *npt;
2275                 temp = w[i__];
2276                 i__2 = i__;
2277                 for (j = 1; j <= i__2; ++j) {
2278 /* L150: */
2279                     bmat[ip + j * bmat_dim1] += temp * w[j];
2280                 }
2281             }
2282         }
2283
2284 /*     The following instructions complete the shift, including the changes */
2285 /*     to the second derivative parameters of the quadratic model. */
2286
2287         ih = 0;
2288         i__2 = *n;
2289         for (j = 1; j <= i__2; ++j) {
2290             w[j] = -half * sumpq * xopt[j];
2291             i__1 = *npt;
2292             for (k = 1; k <= i__1; ++k) {
2293                 w[j] += pq[k] * xpt[k + j * xpt_dim1];
2294 /* L160: */
2295                 xpt[k + j * xpt_dim1] -= xopt[j];
2296             }
2297             i__1 = j;
2298             for (i__ = 1; i__ <= i__1; ++i__) {
2299                 ++ih;
2300                 hq[ih] = hq[ih] + w[i__] * xopt[j] + xopt[i__] * w[j];
2301 /* L170: */
2302                 bmat[*npt + i__ + j * bmat_dim1] = bmat[*npt + j + i__ *
2303                         bmat_dim1];
2304             }
2305         }
2306         i__1 = *n;
2307         for (i__ = 1; i__ <= i__1; ++i__) {
2308             xbase[i__] += xopt[i__];
2309             xnew[i__] -= xopt[i__];
2310             sl[i__] -= xopt[i__];
2311             su[i__] -= xopt[i__];
2312 /* L180: */
2313             xopt[i__] = zero;
2314         }
2315         xoptsq = zero;
2316     }
2317     if (ntrits == 0) {
2318         goto L210;
2319     }
2320     goto L230;
2321
2322 /*     XBASE is also moved to XOPT by a call of RESCUE. This calculation is */
2323 /*     more expensive than the previous shift, because new matrices BMAT and */
2324 /*     ZMAT are generated from scratch, which may include the replacement of */
2325 /*     interpolation points whose positions seem to be causing near linear */
2326 /*     dependence in the interpolation conditions. Therefore RESCUE is called */
2327 /*     only if rounding errors have reduced by at least a factor of two the */
2328 /*     denominator of the formula for updating the H matrix. It provides a */
2329 /*     useful safeguard, but is not invoked in most applications of BOBYQA. */
2330
2331 L190:
2332     nfsav = stop->nevals;
2333     kbase = kopt;
2334     rc2 = rescue_(n, npt, &xl[1], &xu[1],
2335                   stop, calfun, calfun_data,
2336                   &xbase[1], &xpt[xpt_offset], &fval[1], &xopt[1], &gopt[1],
2337                   &hq[1], &pq[1], &bmat[bmat_offset], &zmat[zmat_offset], ndim,
2338                   &sl[1], &su[1], &delta, &kopt, &vlag[1],
2339                   &w[1], &w[*n + np], &w[*ndim + np]);
2340
2341 /*     XOPT is updated now in case the branch below to label 720 is taken. */
2342 /*     Any updating of GOPT occurs after the branch below to label 20, which */
2343 /*     leads to a trust region iteration as does the branch to label 60. */
2344
2345     xoptsq = zero;
2346     if (kopt != kbase) {
2347         i__1 = *n;
2348         for (i__ = 1; i__ <= i__1; ++i__) {
2349             xopt[i__] = xpt[kopt + i__ * xpt_dim1];
2350 /* L200: */
2351 /* Computing 2nd power */
2352             d__1 = xopt[i__];
2353             xoptsq += d__1 * d__1;
2354         }
2355     }
2356     if (rc2 != NLOPT_SUCCESS) {
2357       rc = rc2;
2358       goto L720;
2359     }
2360     nresc = stop->nevals;
2361     if (nfsav < stop->nevals) {
2362         nfsav = stop->nevals;
2363         goto L20;
2364     }
2365     if (ntrits > 0) {
2366         goto L60;
2367     }
2368
2369 /*     Pick two alternative vectors of variables, relative to XBASE, that */
2370 /*     are suitable as new positions of the KNEW-th interpolation point. */
2371 /*     Firstly, XNEW is set to the point on a line through XOPT and another */
2372 /*     interpolation point that minimizes the predicted value of the next */
2373 /*     denominator, subject to ||XNEW - XOPT|| .LEQ. ADELT and to the SL */
2374 /*     and SU bounds. Secondly, XALT is set to the best feasible point on */
2375 /*     a constrained version of the Cauchy step of the KNEW-th Lagrange */
2376 /*     function, the corresponding value of the square of this function */
2377 /*     being returned in CAUCHY. The choice between these alternatives is */
2378 /*     going to be made when the denominator is calculated. */
2379
2380 L210:
2381     altmov_(n, npt, &xpt[xpt_offset], &xopt[1], &bmat[bmat_offset], &zmat[
2382             zmat_offset], ndim, &sl[1], &su[1], &kopt, &knew, &adelt, &xnew[1]
2383             , &xalt[1], &alpha, &cauchy, &w[1], &w[np], &w[*ndim + 1]);
2384     i__1 = *n;
2385     for (i__ = 1; i__ <= i__1; ++i__) {
2386 /* L220: */
2387         d__[i__] = xnew[i__] - xopt[i__];
2388     }
2389
2390 /*     Calculate VLAG and BETA for the current choice of D. The scalar */
2391 /*     product of D with XPT(K,.) is going to be held in W(NPT+K) for */
2392 /*     use when VQUAD is calculated. */
2393
2394 L230:
2395     i__1 = *npt;
2396     for (k = 1; k <= i__1; ++k) {
2397         suma = zero;
2398         sumb = zero;
2399         sum = zero;
2400         i__2 = *n;
2401         for (j = 1; j <= i__2; ++j) {
2402             suma += xpt[k + j * xpt_dim1] * d__[j];
2403             sumb += xpt[k + j * xpt_dim1] * xopt[j];
2404 /* L240: */
2405             sum += bmat[k + j * bmat_dim1] * d__[j];
2406         }
2407         w[k] = suma * (half * suma + sumb);
2408         vlag[k] = sum;
2409 /* L250: */
2410         w[*npt + k] = suma;
2411     }
2412     beta = zero;
2413     i__1 = nptm;
2414     for (jj = 1; jj <= i__1; ++jj) {
2415         sum = zero;
2416         i__2 = *npt;
2417         for (k = 1; k <= i__2; ++k) {
2418 /* L260: */
2419             sum += zmat[k + jj * zmat_dim1] * w[k];
2420         }
2421         beta -= sum * sum;
2422         i__2 = *npt;
2423         for (k = 1; k <= i__2; ++k) {
2424 /* L270: */
2425             vlag[k] += sum * zmat[k + jj * zmat_dim1];
2426         }
2427     }
2428     dsq = zero;
2429     bsum = zero;
2430     dx = zero;
2431     i__2 = *n;
2432     for (j = 1; j <= i__2; ++j) {
2433 /* Computing 2nd power */
2434         d__1 = d__[j];
2435         dsq += d__1 * d__1;
2436         sum = zero;
2437         i__1 = *npt;
2438         for (k = 1; k <= i__1; ++k) {
2439 /* L280: */
2440             sum += w[k] * bmat[k + j * bmat_dim1];
2441         }
2442         bsum += sum * d__[j];
2443         jp = *npt + j;
2444         i__1 = *n;
2445         for (i__ = 1; i__ <= i__1; ++i__) {
2446 /* L290: */
2447             sum += bmat[jp + i__ * bmat_dim1] * d__[i__];
2448         }
2449         vlag[jp] = sum;
2450         bsum += sum * d__[j];
2451 /* L300: */
2452         dx += d__[j] * xopt[j];
2453     }
2454     beta = dx * dx + dsq * (xoptsq + dx + dx + half * dsq) + beta - bsum;
2455     vlag[kopt] += one;
2456
2457 /*     If NTRITS is zero, the denominator may be increased by replacing */
2458 /*     the step D of ALTMOV by a Cauchy step. Then RESCUE may be called if */
2459 /*     rounding errors have damaged the chosen denominator. */
2460
2461     if (ntrits == 0) {
2462 /* Computing 2nd power */
2463         d__1 = vlag[knew];
2464         denom = d__1 * d__1 + alpha * beta;
2465         if (denom < cauchy && cauchy > zero) {
2466             i__2 = *n;
2467             for (i__ = 1; i__ <= i__2; ++i__) {
2468                 xnew[i__] = xalt[i__];
2469 /* L310: */
2470                 d__[i__] = xnew[i__] - xopt[i__];
2471             }
2472             cauchy = zero;
2473             goto L230;
2474         }
2475 /* Computing 2nd power */
2476         d__1 = vlag[knew];
2477         if (denom <= half * (d__1 * d__1)) {
2478             if (stop->nevals > nresc) {
2479                 goto L190;
2480             }
2481             /* Return from BOBYQA because of much cancellation in a
2482                denominator. */
2483             rc = NLOPT_ROUNDOFF_LIMITED;
2484             goto L720;
2485         }
2486
2487 /*     Alternatively, if NTRITS is positive, then set KNEW to the index of */
2488 /*     the next interpolation point to be deleted to make room for a trust */
2489 /*     region step. Again RESCUE may be called if rounding errors have damaged */
2490 /*     the chosen denominator, which is the reason for attempting to select */
2491 /*     KNEW before calculating the next value of the objective function. */
2492
2493     } else {
2494         delsq = delta * delta;
2495         scaden = zero;
2496         biglsq = zero;
2497         knew = 0;
2498         i__2 = *npt;
2499         for (k = 1; k <= i__2; ++k) {
2500             if (k == kopt) {
2501                 goto L350;
2502             }
2503             hdiag = zero;
2504             i__1 = nptm;
2505             for (jj = 1; jj <= i__1; ++jj) {
2506 /* L330: */
2507 /* Computing 2nd power */
2508                 d__1 = zmat[k + jj * zmat_dim1];
2509                 hdiag += d__1 * d__1;
2510             }
2511 /* Computing 2nd power */
2512             d__1 = vlag[k];
2513             den = beta * hdiag + d__1 * d__1;
2514             distsq = zero;
2515             i__1 = *n;
2516             for (j = 1; j <= i__1; ++j) {
2517 /* L340: */
2518 /* Computing 2nd power */
2519                 d__1 = xpt[k + j * xpt_dim1] - xopt[j];
2520                 distsq += d__1 * d__1;
2521             }
2522 /* Computing MAX */
2523 /* Computing 2nd power */
2524             d__3 = distsq / delsq;
2525             d__1 = one, d__2 = d__3 * d__3;
2526             temp = max(d__1,d__2);
2527             if (temp * den > scaden) {
2528                 scaden = temp * den;
2529                 knew = k;
2530                 denom = den;
2531             }
2532 /* Computing MAX */
2533 /* Computing 2nd power */
2534             d__3 = vlag[k];
2535             d__1 = biglsq, d__2 = temp * (d__3 * d__3);
2536             biglsq = max(d__1,d__2);
2537 L350:
2538             ;
2539         }
2540         if (scaden <= half * biglsq) {
2541             if (stop->nevals > nresc) {
2542                 goto L190;
2543             }
2544             /* Return from BOBYQA because of much cancellation in a
2545                denominator. */
2546             rc = NLOPT_ROUNDOFF_LIMITED;
2547             goto L720;
2548         }
2549     }
2550
2551 /*     Put the variables for the next calculation of the objective function */
2552 /*       in XNEW, with any adjustments for the bounds. */
2553
2554
2555 /*     Calculate the value of the objective function at XBASE+XNEW, unless */
2556 /*       the limit on the number of calculations of F has been reached. */
2557
2558 L360:
2559     i__2 = *n;
2560     for (i__ = 1; i__ <= i__2; ++i__) {
2561 /* Computing MIN */
2562 /* Computing MAX */
2563         d__3 = xl[i__], d__4 = xbase[i__] + xnew[i__];
2564         d__1 = max(d__3,d__4), d__2 = xu[i__];
2565         x[i__] = min(d__1,d__2);
2566         if (xnew[i__] == sl[i__]) {
2567             x[i__] = xl[i__];
2568         }
2569         if (xnew[i__] == su[i__]) {
2570             x[i__] = xu[i__];
2571         }
2572 /* L380: */
2573     }
2574
2575     if (nlopt_stop_forced(stop)) rc = NLOPT_FORCE_STOP;
2576     else if (nlopt_stop_evals(stop)) rc = NLOPT_MAXEVAL_REACHED;
2577     else if (nlopt_stop_time(stop)) rc = NLOPT_MAXTIME_REACHED;
2578     if (rc != NLOPT_SUCCESS) goto L720;
2579
2580     stop->nevals++;
2581     f = calfun(*n, &x[1], calfun_data);
2582     if (ntrits == -1) {
2583         fsave = f;
2584         rc = NLOPT_XTOL_REACHED;
2585         if (fsave < fval[kopt]) { *minf = f; return rc; }
2586         goto L720;
2587     }
2588
2589     if (f < stop->minf_max) {
2590       *minf = f;
2591       return NLOPT_MINF_MAX_REACHED;
2592     }
2593
2594 /*     Use the quadratic model to predict the change in F due to the step D, */
2595 /*       and set DIFF to the error of this prediction. */
2596
2597     fopt = fval[kopt];
2598     vquad = zero;
2599     ih = 0;
2600     i__2 = *n;
2601     for (j = 1; j <= i__2; ++j) {
2602         vquad += d__[j] * gopt[j];
2603         i__1 = j;
2604         for (i__ = 1; i__ <= i__1; ++i__) {
2605             ++ih;
2606             temp = d__[i__] * d__[j];
2607             if (i__ == j) {
2608                 temp = half * temp;
2609             }
2610 /* L410: */
2611             vquad += hq[ih] * temp;
2612         }
2613     }
2614     i__1 = *npt;
2615     for (k = 1; k <= i__1; ++k) {
2616 /* L420: */
2617 /* Computing 2nd power */
2618         d__1 = w[*npt + k];
2619         vquad += half * pq[k] * (d__1 * d__1);
2620     }
2621     diff = f - fopt - vquad;
2622     diffc = diffb;
2623     diffb = diffa;
2624     diffa = fabs(diff);
2625     if (dnorm > rho) {
2626         nfsav = stop->nevals;
2627     }
2628
2629 /*     Pick the next value of DELTA after a trust region step. */
2630
2631     if (ntrits > 0) {
2632         if (vquad >= zero) {
2633           /* Return from BOBYQA because a trust region step has failed
2634              to reduce Q. */
2635           rc = NLOPT_ROUNDOFF_LIMITED; /* or FTOL_REACHED? */
2636           goto L720;
2637         }
2638         ratio = (f - fopt) / vquad;
2639         if (ratio <= tenth) {
2640 /* Computing MIN */
2641             d__1 = half * delta;
2642             delta = min(d__1,dnorm);
2643         } else if (ratio <= .7) {
2644 /* Computing MAX */
2645             d__1 = half * delta;
2646             delta = max(d__1,dnorm);
2647         } else {
2648 /* Computing MAX */
2649             d__1 = half * delta, d__2 = dnorm + dnorm;
2650             delta = max(d__1,d__2);
2651         }
2652         if (delta <= rho * 1.5) {
2653             delta = rho;
2654         }
2655
2656 /*     Recalculate KNEW and DENOM if the new F is less than FOPT. */
2657
2658         if (f < fopt) {
2659             ksav = knew;
2660             densav = denom;
2661             delsq = delta * delta;
2662             scaden = zero;
2663             biglsq = zero;
2664             knew = 0;
2665             i__1 = *npt;
2666             for (k = 1; k <= i__1; ++k) {
2667                 hdiag = zero;
2668                 i__2 = nptm;
2669                 for (jj = 1; jj <= i__2; ++jj) {
2670 /* L440: */
2671 /* Computing 2nd power */
2672                     d__1 = zmat[k + jj * zmat_dim1];
2673                     hdiag += d__1 * d__1;
2674                 }
2675 /* Computing 2nd power */
2676                 d__1 = vlag[k];
2677                 den = beta * hdiag + d__1 * d__1;
2678                 distsq = zero;
2679                 i__2 = *n;
2680                 for (j = 1; j <= i__2; ++j) {
2681 /* L450: */
2682 /* Computing 2nd power */
2683                     d__1 = xpt[k + j * xpt_dim1] - xnew[j];
2684                     distsq += d__1 * d__1;
2685                 }
2686 /* Computing MAX */
2687 /* Computing 2nd power */
2688                 d__3 = distsq / delsq;
2689                 d__1 = one, d__2 = d__3 * d__3;
2690                 temp = max(d__1,d__2);
2691                 if (temp * den > scaden) {
2692                     scaden = temp * den;
2693                     knew = k;
2694                     denom = den;
2695                 }
2696 /* L460: */
2697 /* Computing MAX */
2698 /* Computing 2nd power */
2699                 d__3 = vlag[k];
2700                 d__1 = biglsq, d__2 = temp * (d__3 * d__3);
2701                 biglsq = max(d__1,d__2);
2702             }
2703             if (scaden <= half * biglsq) {
2704                 knew = ksav;
2705                 denom = densav;
2706             }
2707         }
2708     }
2709
2710 /*     Update BMAT and ZMAT, so that the KNEW-th interpolation point can be */
2711 /*     moved. Also update the second derivative terms of the model. */
2712
2713     update_(n, npt, &bmat[bmat_offset], &zmat[zmat_offset], ndim, &vlag[1], &
2714             beta, &denom, &knew, &w[1]);
2715     ih = 0;
2716     pqold = pq[knew];
2717     pq[knew] = zero;
2718     i__1 = *n;
2719     for (i__ = 1; i__ <= i__1; ++i__) {
2720         temp = pqold * xpt[knew + i__ * xpt_dim1];
2721         i__2 = i__;
2722         for (j = 1; j <= i__2; ++j) {
2723             ++ih;
2724 /* L470: */
2725             hq[ih] += temp * xpt[knew + j * xpt_dim1];
2726         }
2727     }
2728     i__2 = nptm;
2729     for (jj = 1; jj <= i__2; ++jj) {
2730         temp = diff * zmat[knew + jj * zmat_dim1];
2731         i__1 = *npt;
2732         for (k = 1; k <= i__1; ++k) {
2733 /* L480: */
2734             pq[k] += temp * zmat[k + jj * zmat_dim1];
2735         }
2736     }
2737
2738 /*     Include the new interpolation point, and make the changes to GOPT at */
2739 /*     the old XOPT that are caused by the updating of the quadratic model. */
2740
2741     fval[knew] = f;
2742     i__1 = *n;
2743     for (i__ = 1; i__ <= i__1; ++i__) {
2744         xpt[knew + i__ * xpt_dim1] = xnew[i__];
2745 /* L490: */
2746         w[i__] = bmat[knew + i__ * bmat_dim1];
2747     }
2748     i__1 = *npt;
2749     for (k = 1; k <= i__1; ++k) {
2750         suma = zero;
2751         i__2 = nptm;
2752         for (jj = 1; jj <= i__2; ++jj) {
2753 /* L500: */
2754             suma += zmat[knew + jj * zmat_dim1] * zmat[k + jj * zmat_dim1];
2755         }
2756         if (nlopt_isinf(suma)) {
2757           /* SGJ: detect singularity here (happend if we run
2758              for too many iterations) ... is there another way to recover? */
2759           rc = NLOPT_ROUNDOFF_LIMITED;
2760           goto L720;
2761         }
2762         sumb = zero;
2763         i__2 = *n;
2764         for (j = 1; j <= i__2; ++j) {
2765 /* L510: */
2766             sumb += xpt[k + j * xpt_dim1] * xopt[j];
2767         }
2768         temp = suma * sumb;
2769         i__2 = *n;
2770         for (i__ = 1; i__ <= i__2; ++i__) {
2771 /* L520: */
2772             w[i__] += temp * xpt[k + i__ * xpt_dim1];
2773         }
2774     }
2775     i__2 = *n;
2776     for (i__ = 1; i__ <= i__2; ++i__) {
2777 /* L530: */
2778         gopt[i__] += diff * w[i__];
2779     }
2780
2781 /*     Update XOPT, GOPT and KOPT if the new calculated F is less than FOPT. */
2782
2783     if (f < fopt) {
2784         kopt = knew;
2785         xoptsq = zero;
2786         ih = 0;
2787         i__2 = *n;
2788         for (j = 1; j <= i__2; ++j) {
2789             xopt[j] = xnew[j];
2790 /* Computing 2nd power */
2791             d__1 = xopt[j];
2792             xoptsq += d__1 * d__1;
2793             i__1 = j;
2794             for (i__ = 1; i__ <= i__1; ++i__) {
2795                 ++ih;
2796                 if (i__ < j) {
2797                     gopt[j] += hq[ih] * d__[i__];
2798                 }
2799 /* L540: */
2800                 gopt[i__] += hq[ih] * d__[j];
2801             }
2802         }
2803         i__1 = *npt;
2804         for (k = 1; k <= i__1; ++k) {
2805             temp = zero;
2806             i__2 = *n;
2807             for (j = 1; j <= i__2; ++j) {
2808 /* L550: */
2809                 temp += xpt[k + j * xpt_dim1] * d__[j];
2810             }
2811             temp = pq[k] * temp;
2812             i__2 = *n;
2813             for (i__ = 1; i__ <= i__2; ++i__) {
2814 /* L560: */
2815                 gopt[i__] += temp * xpt[k + i__ * xpt_dim1];
2816             }
2817         }
2818         if (nlopt_stop_ftol(stop, f, fopt)) {
2819              rc = NLOPT_FTOL_REACHED;
2820              goto L720;
2821         }
2822     }
2823
2824 /*     Calculate the parameters of the least Frobenius norm interpolant to */
2825 /*     the current data, the gradient of this interpolant at XOPT being put */
2826 /*     into VLAG(NPT+I), I=1,2,...,N. */
2827
2828     if (ntrits > 0) {
2829         i__2 = *npt;
2830         for (k = 1; k <= i__2; ++k) {
2831             vlag[k] = fval[k] - fval[kopt];
2832 /* L570: */
2833             w[k] = zero;
2834         }
2835         i__2 = nptm;
2836         for (j = 1; j <= i__2; ++j) {
2837             sum = zero;
2838             i__1 = *npt;
2839             for (k = 1; k <= i__1; ++k) {
2840 /* L580: */
2841                 sum += zmat[k + j * zmat_dim1] * vlag[k];
2842             }
2843             i__1 = *npt;
2844             for (k = 1; k <= i__1; ++k) {
2845 /* L590: */
2846                 w[k] += sum * zmat[k + j * zmat_dim1];
2847             }
2848         }
2849         i__1 = *npt;
2850         for (k = 1; k <= i__1; ++k) {
2851             sum = zero;
2852             i__2 = *n;
2853             for (j = 1; j <= i__2; ++j) {
2854 /* L600: */
2855                 sum += xpt[k + j * xpt_dim1] * xopt[j];
2856             }
2857             w[k + *npt] = w[k];
2858 /* L610: */
2859             w[k] = sum * w[k];
2860         }
2861         gqsq = zero;
2862         gisq = zero;
2863         i__1 = *n;
2864         for (i__ = 1; i__ <= i__1; ++i__) {
2865             sum = zero;
2866             i__2 = *npt;
2867             for (k = 1; k <= i__2; ++k) {
2868 /* L620: */
2869                 sum = sum + bmat[k + i__ * bmat_dim1] * vlag[k] + xpt[k + i__
2870                         * xpt_dim1] * w[k];
2871             }
2872             if (xopt[i__] == sl[i__]) {
2873 /* Computing MIN */
2874                 d__2 = zero, d__3 = gopt[i__];
2875 /* Computing 2nd power */
2876                 d__1 = min(d__2,d__3);
2877                 gqsq += d__1 * d__1;
2878 /* Computing 2nd power */
2879                 d__1 = min(zero,sum);
2880                 gisq += d__1 * d__1;
2881             } else if (xopt[i__] == su[i__]) {
2882 /* Computing MAX */
2883                 d__2 = zero, d__3 = gopt[i__];
2884 /* Computing 2nd power */
2885                 d__1 = max(d__2,d__3);
2886                 gqsq += d__1 * d__1;
2887 /* Computing 2nd power */
2888                 d__1 = max(zero,sum);
2889                 gisq += d__1 * d__1;
2890             } else {
2891 /* Computing 2nd power */
2892                 d__1 = gopt[i__];
2893                 gqsq += d__1 * d__1;
2894                 gisq += sum * sum;
2895             }
2896 /* L630: */
2897             vlag[*npt + i__] = sum;
2898         }
2899
2900 /*     Test whether to replace the new quadratic model by the least Frobenius */
2901 /*     norm interpolant, making the replacement if the test is satisfied. */
2902
2903         ++itest;
2904         if (gqsq < ten * gisq) {
2905             itest = 0;
2906         }
2907         if (itest >= 3) {
2908             i__1 = max(*npt,nh);
2909             for (i__ = 1; i__ <= i__1; ++i__) {
2910                 if (i__ <= *n) {
2911                     gopt[i__] = vlag[*npt + i__];
2912                 }
2913                 if (i__ <= *npt) {
2914                     pq[i__] = w[*npt + i__];
2915                 }
2916                 if (i__ <= nh) {
2917                     hq[i__] = zero;
2918                 }
2919                 itest = 0;
2920 /* L640: */
2921             }
2922         }
2923     }
2924
2925 /*     If a trust region step has provided a sufficient decrease in F, then */
2926 /*     branch for another trust region calculation. The case NTRITS=0 occurs */
2927 /*     when the new interpolation point was reached by an alternative step. */
2928
2929     if (ntrits == 0) {
2930         goto L60;
2931     }
2932     if (f <= fopt + tenth * vquad) {
2933         goto L60;
2934     }
2935
2936 /*     Alternatively, find out if the interpolation points are close enough */
2937 /*       to the best point so far. */
2938
2939 /* Computing MAX */
2940 /* Computing 2nd power */
2941     d__3 = two * delta;
2942 /* Computing 2nd power */
2943     d__4 = ten * rho;
2944     d__1 = d__3 * d__3, d__2 = d__4 * d__4;
2945     distsq = max(d__1,d__2);
2946 L650:
2947     knew = 0;
2948     i__1 = *npt;
2949     for (k = 1; k <= i__1; ++k) {
2950         sum = zero;
2951         i__2 = *n;
2952         for (j = 1; j <= i__2; ++j) {
2953 /* L660: */
2954 /* Computing 2nd power */
2955             d__1 = xpt[k + j * xpt_dim1] - xopt[j];
2956             sum += d__1 * d__1;
2957         }
2958         if (sum > distsq) {
2959             knew = k;
2960             distsq = sum;
2961         }
2962 /* L670: */
2963     }
2964
2965 /*     If KNEW is positive, then ALTMOV finds alternative new positions for */
2966 /*     the KNEW-th interpolation point within distance ADELT of XOPT. It is */
2967 /*     reached via label 90. Otherwise, there is a branch to label 60 for */
2968 /*     another trust region iteration, unless the calculations with the */
2969 /*     current RHO are complete. */
2970
2971     if (knew > 0) {
2972         dist = sqrt(distsq);
2973         if (ntrits == -1) {
2974 /* Computing MIN */
2975             d__1 = tenth * delta, d__2 = half * dist;
2976             delta = min(d__1,d__2);
2977             if (delta <= rho * 1.5) {
2978                 delta = rho;
2979             }
2980         }
2981         ntrits = 0;
2982 /* Computing MAX */
2983 /* Computing MIN */
2984         d__2 = tenth * dist;
2985         d__1 = min(d__2,delta);
2986         adelt = max(d__1,rho);
2987         dsq = adelt * adelt;
2988         goto L90;
2989     }
2990     if (ntrits == -1) {
2991         goto L680;
2992     }
2993     if (ratio > zero) {
2994         goto L60;
2995     }
2996     if (max(delta,dnorm) > rho) {
2997         goto L60;
2998     }
2999
3000 /*     The calculations with the current value of RHO are complete. Pick the */
3001 /*       next values of RHO and DELTA. */
3002
3003 L680:
3004     if (rho > *rhoend) {
3005         delta = half * rho;
3006         ratio = rho / *rhoend;
3007         if (ratio <= 16.) {
3008             rho = *rhoend;
3009         } else if (ratio <= 250.) {
3010             rho = sqrt(ratio) * *rhoend;
3011         } else {
3012             rho = tenth * rho;
3013         }
3014         delta = max(delta,rho);
3015         ntrits = 0;
3016         nfsav = stop->nevals;
3017         goto L60;
3018     }
3019
3020 /*     Return from the calculation, after another Newton-Raphson step, if */
3021 /*       it is too short to have been tried before. */
3022
3023     if (ntrits == -1) {
3024         goto L360;
3025     }
3026 L720:
3027     /* originally: if (fval[kopt] <= fsave) -- changed by SGJ, since
3028        this seems like a slight optimization to not update x[]
3029        unnecessarily, at the expense of possibly not returning the
3030        best x[] found so far if the algorithm is stopped suddenly
3031        (e.g. runs out of time) ... it seems safer to execute this
3032        unconditionally, and the efficiency loss seems negligible. */
3033     {
3034         i__1 = *n;
3035         for (i__ = 1; i__ <= i__1; ++i__) {
3036 /* Computing MIN */
3037 /* Computing MAX */
3038             d__3 = xl[i__], d__4 = xbase[i__] + xopt[i__];
3039             d__1 = max(d__3,d__4), d__2 = xu[i__];
3040             x[i__] = min(d__1,d__2);
3041             if (xopt[i__] == sl[i__]) {
3042                 x[i__] = xl[i__];
3043             }
3044             if (xopt[i__] == su[i__]) {
3045                 x[i__] = xu[i__];
3046             }
3047 /* L730: */
3048         }
3049         f = fval[kopt];
3050     }
3051     *minf = f;
3052     return rc;
3053 } /* bobyqb_ */
3054
3055 /**************************************************************************/
3056
3057 nlopt_result bobyqa(int n, int npt, double *x,
3058                     const double *xl, const double *xu, double rhobeg,
3059                     nlopt_stopping *stop, double *minf,
3060                     bobyqa_func calfun, void *calfun_data)
3061 {
3062     /* System generated locals */
3063     int i__1;
3064     double d__1, d__2;
3065
3066     /* Local variables */
3067     int j, id, np, iw, igo, ihq, ixb, ixa, ifv, isl, jsl, ipq, ivl, ixn,
3068             ixo, ixp, isu, jsu, ndim;
3069     double temp, zero;
3070     int ibmat, izmat;
3071
3072     double rhoend;
3073     double *w;
3074     nlopt_result ret;
3075
3076 /* SGJ, 2009: compute rhoend from NLopt stop info */
3077     rhoend = stop->xtol_rel * (rhobeg);
3078     for (j = 0; j < n; ++j)
3079          if (rhoend < stop->xtol_abs[j])
3080               rhoend = stop->xtol_abs[j];
3081
3082
3083 /*     This subroutine seeks the least value of a function of many variables, */
3084 /*     by applying a trust region method that forms quadratic models by */
3085 /*     interpolation. There is usually some freedom in the interpolation */
3086 /*     conditions, which is taken up by minimizing the Frobenius norm of */
3087 /*     the change to the second derivative of the model, beginning with the */
3088 /*     zero matrix. The values of the variables are constrained by upper and */
3089 /*     lower bounds. The arguments of the subroutine are as follows. */
3090
3091 /*     N must be set to the number of variables and must be at least two. */
3092 /*     NPT is the number of interpolation conditions. Its value must be in */
3093 /*       the interval [N+2,(N+1)(N+2)/2]. Choices that exceed 2*N+1 are not */
3094 /*       recommended. */
3095 /*     Initial values of the variables must be set in X(1),X(2),...,X(N). They */
3096 /*       will be changed to the values that give the least calculated F. */
3097 /*     For I=1,2,...,N, XL(I) and XU(I) must provide the lower and upper */
3098 /*       bounds, respectively, on X(I). The construction of quadratic models */
3099 /*       requires XL(I) to be strictly less than XU(I) for each I. Further, */
3100 /*       the contribution to a model from changes to the I-th variable is */
3101 /*       damaged severely by rounding errors if XU(I)-XL(I) is too small. */
3102 /*     RHOBEG and RHOEND must be set to the initial and final values of a trust */
3103 /*       region radius, so both must be positive with RHOEND no greater than */
3104 /*       RHOBEG. Typically, RHOBEG should be about one tenth of the greatest */
3105 /*       expected change to a variable, while RHOEND should indicate the */
3106 /*       accuracy that is required in the final values of the variables. An */
3107 /*       error return occurs if any of the differences XU(I)-XL(I), I=1,...,N, */
3108 /*       is less than 2*RHOBEG. */
3109 /*     MAXFUN must be set to an upper bound on the number of calls of CALFUN. */
3110 /*     The array W will be used for working space. Its length must be at least */
3111 /*       (NPT+5)*(NPT+N)+3*N*(N+5)/2. */
3112
3113 /*     SUBROUTINE CALFUN (N,X,F) has to be provided by the user. It must set */
3114 /*     F to the value of the objective function for the current values of the */
3115 /*     variables X(1),X(2),...,X(N), which are generated automatically in a */
3116 /*     way that satisfies the bounds given in XL and XU. */
3117
3118 /*     Return if the value of NPT is unacceptable. */
3119
3120     /* Parameter adjustments */
3121     --xu;
3122     --xl;
3123     --x;
3124
3125     /* Function Body */
3126     np = n + 1;
3127     if (npt < n + 2 || npt > (n + 2) * np / 2) {
3128       /* Return from BOBYQA because NPT is not in the required interval */
3129       return NLOPT_INVALID_ARGS;
3130     }
3131
3132 /*     Partition the working space array, so that different parts of it can */
3133 /*     be treated separately during the calculation of BOBYQB. The partition */
3134 /*     requires the first (NPT+2)*(NPT+N)+3*N*(N+5)/2 elements of W plus the */
3135 /*     space that is taken by the last array in the argument list of BOBYQB. */
3136
3137     ndim = npt + n;
3138     ixb = 1;
3139     ixp = ixb + n;
3140     ifv = ixp + n * npt;
3141     ixo = ifv + npt;
3142     igo = ixo + n;
3143     ihq = igo + n;
3144     ipq = ihq + n * np / 2;
3145     ibmat = ipq + npt;
3146     izmat = ibmat + ndim * n;
3147     isl = izmat + npt * (npt - np);
3148     isu = isl + n;
3149     ixn = isu + n;
3150     ixa = ixn + n;
3151     id = ixa + n;
3152     ivl = id + n;
3153     iw = ivl + ndim;
3154
3155     w = (double *) malloc(sizeof(double) * ((npt+5)*(npt+n)+3*n*(n+5)/2));
3156     if (!w) return NLOPT_OUT_OF_MEMORY;
3157     --w;
3158
3159 /*   Return if there is insufficient space between the bounds. Modify the */
3160 /*   initial X if necessary in order to avoid conflicts between the bounds */
3161 /*   and the construction of the first quadratic model. The lower and upper */
3162 /*   bounds on moves from the updated X are set now, in the ISL and ISU */
3163 /*   partitions of W, in order to provide useful and exact information about */
3164 /*   components of X that become within distance RHOBEG from their bounds. */
3165
3166     zero = 0.;
3167     i__1 = n;
3168     for (j = 1; j <= i__1; ++j) {
3169         temp = xu[j] - xl[j];
3170         if (temp < rhobeg + rhobeg) {
3171           /* Return from BOBYQA because one of the differences
3172              XU(I)-XL(I)s is less than 2*RHOBEG. */
3173           free(w+1);
3174           return NLOPT_INVALID_ARGS;
3175         }
3176         jsl = isl + j - 1;
3177         jsu = jsl + n;
3178         w[jsl] = xl[j] - x[j];
3179         w[jsu] = xu[j] - x[j];
3180         if (w[jsl] >= -(rhobeg)) {
3181             if (w[jsl] >= zero) {
3182                 x[j] = xl[j];
3183                 w[jsl] = zero;
3184                 w[jsu] = temp;
3185             } else {
3186                 x[j] = xl[j] + rhobeg;
3187                 w[jsl] = -(rhobeg);
3188 /* Computing MAX */
3189                 d__1 = xu[j] - x[j];
3190                 w[jsu] = max(d__1,rhobeg);
3191             }
3192         } else if (w[jsu] <= rhobeg) {
3193             if (w[jsu] <= zero) {
3194                 x[j] = xu[j];
3195                 w[jsl] = -temp;
3196                 w[jsu] = zero;
3197             } else {
3198                 x[j] = xu[j] - rhobeg;
3199 /* Computing MIN */
3200                 d__1 = xl[j] - x[j], d__2 = -(rhobeg);
3201                 w[jsl] = min(d__1,d__2);
3202                 w[jsu] = rhobeg;
3203             }
3204         }
3205 /* L30: */
3206     }
3207
3208 /*     Make the call of BOBYQB. */
3209
3210     ret = bobyqb_(&n, &npt, &x[1], &xl[1], &xu[1], &rhobeg, &rhoend,
3211                   stop, calfun, calfun_data, minf,
3212                   &w[ixb], &w[ixp], &w[ifv], &w[ixo], &w[igo], &w[ihq], &w[ipq],
3213                   &w[ibmat], &w[izmat], &ndim, &w[isl], &w[isu], &w[ixn], &w[ixa],
3214                   &w[id], &w[ivl], &w[iw]);
3215     free(w+1);
3216     return ret;
3217 } /* bobyqa_ */
3218