newuoa/newuoa.c

   1 /* Copyright (c) 2004 M. J. D. Powell (mjdp@cam.ac.uk)
   2  * Copyright (c) 2007-2014 Massachusetts Institute of Technology
   3  *
   4  * Permission is hereby granted, free of charge, to any person obtaining
   5  * a copy of this software and associated documentation files (the
   6  * "Software"), to deal in the Software without restriction, including
   7  * without limitation the rights to use, copy, modify, merge, publish,
   8  * distribute, sublicense, and/or sell copies of the Software, and to
   9  * permit persons to whom the Software is furnished to do so, subject to
  10  * the following conditions:
  11  *
  12  * The above copyright notice and this permission notice shall be
  13  * included in all copies or substantial portions of the Software.
  14  *
  15  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  16  * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  17  * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  18  * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
  19  * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
  20  * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
  21  * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  22  */
  23
  24 /* NEWUOA derivative-free optimization algorithm by M. J. D. Powell.
  25    Original Fortran code by Powell (2004).  Converted via f2c, cleaned up,
  26    and incorporated into NLopt by S. G. Johnson (2008).  See README. */
  27
  28 #include <math.h>
  29 #include <stdlib.h>
  30 #include <stdio.h>
  31 #include <string.h>
  32
  33 #include "newuoa.h"
  34
  35 #define MIN2(a,b) ((a) <= (b) ? (a) : (b))
  36 #define MAX2(a,b) ((a) >= (b) ? (a) : (b))
  37
  38 /*************************************************************************/
  39 /* trsapp.f */
  40
  41 typedef struct {
  42      int npt;
  43      double *xpt, *pq, *hq, *gq, *xopt;
  44      double *hd;
  45      int iter;
  46 } quad_model_data;
  47
  48 static double quad_model(int n, const double *x, double *grad, void *data)
  49 {
  50      quad_model_data *d = (quad_model_data *) data;
  51      const double *xpt = d->xpt, *pq = d->pq, *hq = d->hq, *gq = d->gq, *xopt = d->xopt;
  52      double *hd = d->hd;
  53      int npt = d->npt;
  54      int i, j, k;
  55      double val = 0;
  56
  57      /* first, set hd to be Hessian matrix times x */
  58      memset(hd, 0, sizeof(double) * n);
  59      /* implicit Hessian terms (a sum of outer products of xpt vectors) */
  60      for (k = 0; k < npt; ++k) {
  61           double temp = 0;
  62           for (j = 0; j < n; ++j)
  63                temp += xpt[k + j * npt] * (xopt[j] + x[j]);
  64           temp *= pq[k];
  65           for (i = 0; i < n; ++i)
  66                hd[i] += temp * xpt[k + i * npt];
  67      }
  68      /* explicit Hessian terms (stored as compressed lower triangle hq) */
  69      k = 0;
  70      for (j = 0; j < n; ++j) {
  71           for (i = 0; i < j; ++i) {
  72                hd[j] += hq[k] * (xopt[i] + x[i]);
  73                hd[i] += hq[k] * (xopt[j] + x[j]);
  74                ++k;
  75           }
  76           hd[j] += hq[k++] * (xopt[j] + x[j]);
  77      }
  78
  79      for (i = 0; i < n; ++i) {
  80           val += (gq[i] + 0.5 * hd[i]) * (xopt[i] + x[i]);
  81           if (grad) grad[i] = gq[i] + hd[i];
  82      }
  83      d->iter++;
  84      return val;
  85 }
  86
  87 /* constraint function to enforce |x|^2 <= rho*rho */
  88 static double rho_constraint(int n, const double *x, double *grad, void *data)
  89 {
  90      double rho = *((double *) data);
  91      double val = - rho*rho;
  92      int i;
  93      for (i = 0; i < n; ++i)
  94           val += x[i] * x[i];
  95      if (grad) for (i = 0; i < n; ++i) grad[i] = 2*x[i];
  96      return val;
  97 }
  98
  99 static nlopt_result trsapp_(int *n, int *npt, double *xopt,
 100         double *xpt, double *gq, double *hq, double *pq,
 101         double *delta, double *step, double *d__, double *g,
 102         double *hd, double *hs, double *crvmin,
 103                     const double *xbase, const double *lb, const double *ub)
 104 {
 105     /* System generated locals */
 106     int xpt_dim1, xpt_offset, i__1, i__2;
 107     double d__1, d__2;
 108
 109     /* Local variables */
 110     int i__, j, k;
 111     double dd = 0.0, cf, dg, gg = 0.0;
 112     int ih;
 113     double ds, sg = 0.0;
 114     int iu;
 115     double ss, dhd, dhs, cth, sgk, shs = 0.0, sth, qadd, half, qbeg, qred = 0.0, qmin,
 116              temp, qsav, qnew, zero, ggbeg = 0.0, alpha, angle, reduc;
 117     int iterc;
 118     double ggsav, delsq, tempa = 0.0, tempb = 0.0;
 119     int isave;
 120     double bstep = 0.0, ratio, twopi;
 121     int itersw;
 122     double angtest;
 123     int itermax;
 124
 125
 126 /* N is the number of variables of a quadratic objective function, Q say. */
 127 /* The arguments NPT, XOPT, XPT, GQ, HQ and PQ have their usual meanings, */
 128 /*   in order to define the current quadratic model Q. */
 129 /* DELTA is the trust region radius, and has to be positive. */
 130 /* STEP will be set to the calculated trial step. */
 131 /* The arrays D, G, HD and HS will be used for working space. */
 132 /* CRVMIN will be set to the least curvature of H along the conjugate */
 133 /*   directions that occur, except that it is set to zero if STEP goes */
 134 /*   all the way to the trust region boundary. */
 135
 136 /* The calculation of STEP begins with the truncated conjugate gradient */
 137 /* method. If the boundary of the trust region is reached, then further */
 138 /* changes to STEP may be made, each one being in the 2D space spanned */
 139 /* by the current STEP and the corresponding gradient of Q. Thus STEP */
 140 /* should provide a substantial reduction to Q within the trust region. */
 141
 142 /* Initialization, which includes setting HD to H times XOPT. */
 143
 144     /* Parameter adjustments */
 145     xpt_dim1 = *npt;
 146     xpt_offset = 1 + xpt_dim1;
 147     xpt -= xpt_offset;
 148     --xopt;
 149     --gq;
 150     --hq;
 151     --pq;
 152     --step;
 153     --d__;
 154     --g;
 155     --hd;
 156     --hs;
 157
 158     if (lb && ub) {
 159          double *slb, *sub, *xtol, minf, crv;
 160          nlopt_result ret;
 161          quad_model_data qmd;
 162          qmd.npt = *npt;
 163          qmd.xpt = &xpt[1 + 1 * xpt_dim1];
 164          qmd.pq = &pq[1];
 165          qmd.hq = &hq[1];
 166          qmd.gq = &gq[1];
 167          qmd.xopt = &xopt[1];
 168          qmd.hd = &hd[1];
 169          qmd.iter = 0;
 170          slb = &g[1];
 171          sub = &hs[1];
 172          xtol = &d__[1];
 173          for (j = 0; j < *n; ++j) {
 174               slb[j] = -(sub[j] = *delta);
 175               if (slb[j] < lb[j] - xbase[j] - xopt[j+1])
 176                    slb[j] = lb[j] - xbase[j] - xopt[j+1];
 177               if (sub[j] > ub[j] - xbase[j] - xopt[j+1])
 178                    sub[j] = ub[j] - xbase[j] - xopt[j+1];
 179               if (slb[j] > 0) slb[j] = 0;
 180               if (sub[j] < 0) sub[j] = 0;
 181               xtol[j] = 1e-7 * *delta; /* absolute x tolerance */
 182          }
 183          memset(&step[1], 0, sizeof(double) * *n);
 184          ret = nlopt_minimize_constrained(NLOPT_LD_MMA, *n, quad_model, &qmd,
 185                                     1, rho_constraint, delta, sizeof(double),
 186                                     slb, sub, &step[1], &minf, -HUGE_VAL,
 187                                     0., 0., 0., xtol, 1000, 0.);
 188          if (rho_constraint(*n, &step[1], 0, delta) > -1e-6*(*delta)*(*delta))
 189               crv = 0;
 190          else {
 191               for (j = 1; j <= *n; ++j) d__[j] = step[j] - xopt[j];
 192               quad_model(*n, &d__[1], &g[1], &qmd);
 193               crv = gg = 0;
 194               for (j = 1; j <= *n; ++j) {
 195                    crv += step[j] * (g[j] - gq[j]);
 196                    gg += step[j] * step[j];
 197               }
 198               if (gg <= 1e-16 * crv)
 199                    crv = 1e16;
 200               else
 201                    crv = crv / gg;
 202          }
 203          *crvmin = crv;
 204          return ret;
 205     }
 206
 207     /* Function Body */
 208     half = .5;
 209     zero = 0.;
 210     twopi = atan(1.) * 8.;
 211     delsq = *delta * *delta;
 212     iterc = 0;
 213     itermax = *n;
 214     itersw = itermax;
 215     i__1 = *n;
 216     for (i__ = 1; i__ <= i__1; ++i__) {
 217 /* L10: */
 218         d__[i__] = xopt[i__];
 219     }
 220     goto L170;
 221
 222 /* Prepare for the first line search. */
 223
 224 L20:
 225     qred = zero;
 226     dd = zero;
 227     i__1 = *n;
 228     for (i__ = 1; i__ <= i__1; ++i__) {
 229         step[i__] = zero;
 230         hs[i__] = zero;
 231         g[i__] = gq[i__] + hd[i__];
 232         d__[i__] = -g[i__];
 233 /* L30: */
 234 /* Computing 2nd power */
 235         d__1 = d__[i__];
 236         dd += d__1 * d__1;
 237     }
 238     *crvmin = zero;
 239     if (dd == zero) {
 240         goto L160;
 241     }
 242     ds = zero;
 243     ss = zero;
 244     gg = dd;
 245     ggbeg = gg;
 246
 247 /* Calculate the step to the trust region boundary and the product HD. */
 248
 249 L40:
 250     ++iterc;
 251     temp = delsq - ss;
 252     bstep = temp / (ds + sqrt(ds * ds + dd * temp));
 253     goto L170;
 254 L50:
 255     dhd = zero;
 256     i__1 = *n;
 257     for (j = 1; j <= i__1; ++j) {
 258 /* L60: */
 259         dhd += d__[j] * hd[j];
 260     }
 261
 262 /* Update CRVMIN and set the step-length ALPHA. */
 263
 264     alpha = bstep;
 265     if (dhd > zero) {
 266         temp = dhd / dd;
 267         if (iterc == 1) {
 268             *crvmin = temp;
 269         }
 270         *crvmin = MIN2(*crvmin,temp);
 271 /* Computing MIN */
 272         d__1 = alpha, d__2 = gg / dhd;
 273         alpha = MIN2(d__1,d__2);
 274     }
 275     qadd = alpha * (gg - half * alpha * dhd);
 276     qred += qadd;
 277
 278 /* Update STEP and HS. */
 279
 280     ggsav = gg;
 281     gg = zero;
 282     i__1 = *n;
 283     for (i__ = 1; i__ <= i__1; ++i__) {
 284         step[i__] += alpha * d__[i__];
 285         hs[i__] += alpha * hd[i__];
 286 /* L70: */
 287 /* Computing 2nd power */
 288         d__1 = g[i__] + hs[i__];
 289         gg += d__1 * d__1;
 290     }
 291
 292 /* Begin another conjugate direction iteration if required. */
 293
 294     if (alpha < bstep) {
 295         if (qadd <= qred * .01) {
 296             goto L160;
 297         }
 298         if (gg <= ggbeg * 1e-4) {
 299             goto L160;
 300         }
 301         if (iterc == itermax) {
 302             goto L160;
 303         }
 304         temp = gg / ggsav;
 305         dd = zero;
 306         ds = zero;
 307         ss = zero;
 308         i__1 = *n;
 309         for (i__ = 1; i__ <= i__1; ++i__) {
 310             d__[i__] = temp * d__[i__] - g[i__] - hs[i__];
 311 /* Computing 2nd power */
 312             d__1 = d__[i__];
 313             dd += d__1 * d__1;
 314             ds += d__[i__] * step[i__];
 315 /* L80: */
 316 /* Computing 2nd power */
 317             d__1 = step[i__];
 318             ss += d__1 * d__1;
 319         }
 320         if (ds <= zero) {
 321             goto L160;
 322         }
 323         if (ss < delsq) {
 324             goto L40;
 325         }
 326     }
 327     *crvmin = zero;
 328     itersw = iterc;
 329
 330 /* Test whether an alternative iteration is required. */
 331
 332 L90:
 333     if (gg <= ggbeg * 1e-4) {
 334         goto L160;
 335     }
 336     sg = zero;
 337     shs = zero;
 338     i__1 = *n;
 339     for (i__ = 1; i__ <= i__1; ++i__) {
 340         sg += step[i__] * g[i__];
 341 /* L100: */
 342         shs += step[i__] * hs[i__];
 343     }
 344     sgk = sg + shs;
 345     angtest = sgk / sqrt(gg * delsq);
 346     if (angtest <= -.99) {
 347         goto L160;
 348     }
 349
 350 /* Begin the alternative iteration by calculating D and HD and some */
 351 /* scalar products. */
 352
 353     ++iterc;
 354     temp = sqrt(delsq * gg - sgk * sgk);
 355     tempa = delsq / temp;
 356     tempb = sgk / temp;
 357     i__1 = *n;
 358     for (i__ = 1; i__ <= i__1; ++i__) {
 359 /* L110: */
 360         d__[i__] = tempa * (g[i__] + hs[i__]) - tempb * step[i__];
 361     }
 362     goto L170;
 363 L120:
 364     dg = zero;
 365     dhd = zero;
 366     dhs = zero;
 367     i__1 = *n;
 368     for (i__ = 1; i__ <= i__1; ++i__) {
 369         dg += d__[i__] * g[i__];
 370         dhd += hd[i__] * d__[i__];
 371 /* L130: */
 372         dhs += hd[i__] * step[i__];
 373     }
 374
 375 /* Seek the value of the angle that minimizes Q. */
 376
 377     cf = half * (shs - dhd);
 378     qbeg = sg + cf;
 379     qsav = qbeg;
 380     qmin = qbeg;
 381     isave = 0;
 382     iu = 49;
 383     temp = twopi / (double) (iu + 1);
 384     i__1 = iu;
 385     for (i__ = 1; i__ <= i__1; ++i__) {
 386         angle = (double) i__ * temp;
 387         cth = cos(angle);
 388         sth = sin(angle);
 389         qnew = (sg + cf * cth) * cth + (dg + dhs * cth) * sth;
 390         if (qnew < qmin) {
 391             qmin = qnew;
 392             isave = i__;
 393             tempa = qsav;
 394         } else if (i__ == isave + 1) {
 395             tempb = qnew;
 396         }
 397 /* L140: */
 398         qsav = qnew;
 399     }
 400     if ((double) isave == zero) {
 401         tempa = qnew;
 402     }
 403     if (isave == iu) {
 404         tempb = qbeg;
 405     }
 406     angle = zero;
 407     if (tempa != tempb) {
 408         tempa -= qmin;
 409         tempb -= qmin;
 410         angle = half * (tempa - tempb) / (tempa + tempb);
 411     }
 412     angle = temp * ((double) isave + angle);
 413
 414 /* Calculate the new STEP and HS. Then test for convergence. */
 415
 416     cth = cos(angle);
 417     sth = sin(angle);
 418     reduc = qbeg - (sg + cf * cth) * cth - (dg + dhs * cth) * sth;
 419     gg = zero;
 420     i__1 = *n;
 421     for (i__ = 1; i__ <= i__1; ++i__) {
 422         step[i__] = cth * step[i__] + sth * d__[i__];
 423         hs[i__] = cth * hs[i__] + sth * hd[i__];
 424 /* L150: */
 425 /* Computing 2nd power */
 426         d__1 = g[i__] + hs[i__];
 427         gg += d__1 * d__1;
 428     }
 429     qred += reduc;
 430     ratio = reduc / qred;
 431     if (iterc < itermax && ratio > .01) {
 432         goto L90;
 433     }
 434 L160:
 435     return NLOPT_SUCCESS;
 436
 437 /* The following instructions act as a subroutine for setting the vector */
 438 /* HD to the vector D multiplied by the second derivative matrix of Q. */
 439 /* They are called from three different places, which are distinguished */
 440 /* by the value of ITERC. */
 441
 442 L170:
 443     i__1 = *n;
 444     for (i__ = 1; i__ <= i__1; ++i__) {
 445 /* L180: */
 446         hd[i__] = zero;
 447     }
 448     i__1 = *npt;
 449     for (k = 1; k <= i__1; ++k) {
 450         temp = zero;
 451         i__2 = *n;
 452         for (j = 1; j <= i__2; ++j) {
 453 /* L190: */
 454             temp += xpt[k + j * xpt_dim1] * d__[j];
 455         }
 456         temp *= pq[k];
 457         i__2 = *n;
 458         for (i__ = 1; i__ <= i__2; ++i__) {
 459 /* L200: */
 460             hd[i__] += temp * xpt[k + i__ * xpt_dim1];
 461         }
 462     }
 463     ih = 0;
 464     i__2 = *n;
 465     for (j = 1; j <= i__2; ++j) {
 466         i__1 = j;
 467         for (i__ = 1; i__ <= i__1; ++i__) {
 468             ++ih;
 469             if (i__ < j) {
 470                 hd[j] += hq[ih] * d__[i__];
 471             }
 472 /* L210: */
 473             hd[i__] += hq[ih] * d__[j];
 474         }
 475     }
 476     if (iterc == 0) {
 477         goto L20;
 478     }
 479     if (iterc <= itersw) {
 480         goto L50;
 481     }
 482     goto L120;
 483 } /* trsapp_ */
 484
 485
 486 /*************************************************************************/
 487 /* bigden.f */
 488
 489 static nlopt_result bigden_(int *n, int *npt, double *xopt,
 490                     double *xpt, double *bmat, double *zmat, int *idz,
 491                     int *ndim, int *kopt, int *knew, double *d__,
 492                     double *w, double *vlag, double *beta, double *s,
 493                     double *wvec, double *prod,
 494                     const double *xbase, const double *lb, const double *ub)
 495 {
 496     /* System generated locals */
 497     int xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1,
 498             zmat_offset, wvec_dim1, wvec_offset, prod_dim1, prod_offset, i__1,
 499              i__2;
 500     double d__1;
 501
 502     /* Local variables */
 503     int i__, j, k;
 504     double dd;
 505     int jc;
 506     double ds;
 507     int ip, iu, nw;
 508     double ss, den[9], one, par[9], tau, sum, two, diff, half, temp;
 509     int ksav;
 510     double step;
 511     int nptm;
 512     double zero, alpha, angle, denex[9];
 513     int iterc;
 514     double tempa, tempb, tempc;
 515     int isave;
 516     double ssden, dtest, quart, xoptd, twopi, xopts, denold, denmax,
 517             densav, dstemp, sumold, sstemp, xoptsq;
 518
 519
 520 /* N is the number of variables. */
 521 /* NPT is the number of interpolation equations. */
 522 /* XOPT is the best interpolation point so far. */
 523 /* XPT contains the coordinates of the current interpolation points. */
 524 /* BMAT provides the last N columns of H. */
 525 /* ZMAT and IDZ give a factorization of the first NPT by NPT submatrix of H. */
 526 /* NDIM is the first dimension of BMAT and has the value NPT+N. */
 527 /* KOPT is the index of the optimal interpolation point. */
 528 /* KNEW is the index of the interpolation point that is going to be moved. */
 529 /* D will be set to the step from XOPT to the new point, and on entry it */
 530 /*   should be the D that was calculated by the last call of BIGLAG. The */
 531 /*   length of the initial D provides a trust region bound on the final D. */
 532 /* W will be set to Wcheck for the final choice of D. */
 533 /* VLAG will be set to Theta*Wcheck+e_b for the final choice of D. */
 534 /* BETA will be set to the value that will occur in the updating formula */
 535 /*   when the KNEW-th interpolation point is moved to its new position. */
 536 /* S, WVEC, PROD and the private arrays DEN, DENEX and PAR will be used */
 537 /*   for working space. */
 538
 539 /* D is calculated in a way that should provide a denominator with a large */
 540 /* modulus in the updating formula when the KNEW-th interpolation point is */
 541 /* shifted to the new position XOPT+D. */
 542
 543 /* Set some constants. */
 544
 545     /* Parameter adjustments */
 546     zmat_dim1 = *npt;
 547     zmat_offset = 1 + zmat_dim1;
 548     zmat -= zmat_offset;
 549     xpt_dim1 = *npt;
 550     xpt_offset = 1 + xpt_dim1;
 551     xpt -= xpt_offset;
 552     --xopt;
 553     prod_dim1 = *ndim;
 554     prod_offset = 1 + prod_dim1;
 555     prod -= prod_offset;
 556     wvec_dim1 = *ndim;
 557     wvec_offset = 1 + wvec_dim1;
 558     wvec -= wvec_offset;
 559     bmat_dim1 = *ndim;
 560     bmat_offset = 1 + bmat_dim1;
 561     bmat -= bmat_offset;
 562     --d__;
 563     --w;
 564     --vlag;
 565     --s;
 566
 567     /* Function Body */
 568     half = .5;
 569     one = 1.;
 570     quart = .25;
 571     two = 2.;
 572     zero = 0.;
 573     twopi = atan(one) * 8.;
 574     nptm = *npt - *n - 1;
 575
 576 /* Store the first NPT elements of the KNEW-th column of H in W(N+1) */
 577 /* to W(N+NPT). */
 578
 579     i__1 = *npt;
 580     for (k = 1; k <= i__1; ++k) {
 581 /* L10: */
 582         w[*n + k] = zero;
 583     }
 584     i__1 = nptm;
 585     for (j = 1; j <= i__1; ++j) {
 586         temp = zmat[*knew + j * zmat_dim1];
 587         if (j < *idz) {
 588             temp = -temp;
 589         }
 590         i__2 = *npt;
 591         for (k = 1; k <= i__2; ++k) {
 592 /* L20: */
 593             w[*n + k] += temp * zmat[k + j * zmat_dim1];
 594         }
 595     }
 596     alpha = w[*n + *knew];
 597
 598 /* The initial search direction D is taken from the last call of BIGLAG, */
 599 /* and the initial S is set below, usually to the direction from X_OPT */
 600 /* to X_KNEW, but a different direction to an interpolation point may */
 601 /* be chosen, in order to prevent S from being nearly parallel to D. */
 602
 603     dd = zero;
 604     ds = zero;
 605     ss = zero;
 606     xoptsq = zero;
 607     i__2 = *n;
 608     for (i__ = 1; i__ <= i__2; ++i__) {
 609 /* Computing 2nd power */
 610         d__1 = d__[i__];
 611         dd += d__1 * d__1;
 612         s[i__] = xpt[*knew + i__ * xpt_dim1] - xopt[i__];
 613         ds += d__[i__] * s[i__];
 614 /* Computing 2nd power */
 615         d__1 = s[i__];
 616         ss += d__1 * d__1;
 617 /* L30: */
 618 /* Computing 2nd power */
 619         d__1 = xopt[i__];
 620         xoptsq += d__1 * d__1;
 621     }
 622     if (ds * ds > dd * .99 * ss) {
 623         ksav = *knew;
 624         dtest = ds * ds / ss;
 625         i__2 = *npt;
 626         for (k = 1; k <= i__2; ++k) {
 627             if (k != *kopt) {
 628                 dstemp = zero;
 629                 sstemp = zero;
 630                 i__1 = *n;
 631                 for (i__ = 1; i__ <= i__1; ++i__) {
 632                     diff = xpt[k + i__ * xpt_dim1] - xopt[i__];
 633                     dstemp += d__[i__] * diff;
 634 /* L40: */
 635                     sstemp += diff * diff;
 636                 }
 637                 if (sstemp == 0) return NLOPT_ROUNDOFF_LIMITED;
 638                 if (dstemp * dstemp / sstemp < dtest) {
 639                     ksav = k;
 640                     dtest = dstemp * dstemp / sstemp;
 641                     ds = dstemp;
 642                     ss = sstemp;
 643                 }
 644             }
 645 /* L50: */
 646         }
 647         i__2 = *n;
 648         for (i__ = 1; i__ <= i__2; ++i__) {
 649 /* L60: */
 650             s[i__] = xpt[ksav + i__ * xpt_dim1] - xopt[i__];
 651         }
 652     }
 653     ssden = dd * ss - ds * ds;
 654     iterc = 0;
 655     densav = zero;
 656
 657 /* Begin the iteration by overwriting S with a vector that has the */
 658 /* required length and direction. */
 659
 660 L70:
 661     ++iterc;
 662     if (ssden < 0) return NLOPT_ROUNDOFF_LIMITED;
 663     temp = one / sqrt(ssden);
 664     xoptd = zero;
 665     xopts = zero;
 666     i__2 = *n;
 667     for (i__ = 1; i__ <= i__2; ++i__) {
 668         s[i__] = temp * (dd * s[i__] - ds * d__[i__]);
 669         xoptd += xopt[i__] * d__[i__];
 670 /* L80: */
 671         if (nlopt_isinf(s[i__])) return NLOPT_ROUNDOFF_LIMITED;
 672         xopts += xopt[i__] * s[i__];
 673     }
 674
 675 /* Set the coefficients of the first two terms of BETA. */
 676
 677     tempa = half * xoptd * xoptd;
 678     tempb = half * xopts * xopts;
 679     den[0] = dd * (xoptsq + half * dd) + tempa + tempb;
 680     den[1] = two * xoptd * dd;
 681     den[2] = two * xopts * dd;
 682     den[3] = tempa - tempb;
 683     den[4] = xoptd * xopts;
 684     for (i__ = 6; i__ <= 9; ++i__) {
 685 /* L90: */
 686         den[i__ - 1] = zero;
 687     }
 688
 689 /* Put the coefficients of Wcheck in WVEC. */
 690
 691     i__2 = *npt;
 692     for (k = 1; k <= i__2; ++k) {
 693         tempa = zero;
 694         tempb = zero;
 695         tempc = zero;
 696         i__1 = *n;
 697         for (i__ = 1; i__ <= i__1; ++i__) {
 698             tempa += xpt[k + i__ * xpt_dim1] * d__[i__];
 699             tempb += xpt[k + i__ * xpt_dim1] * s[i__];
 700 /* L100: */
 701             tempc += xpt[k + i__ * xpt_dim1] * xopt[i__];
 702         }
 703         wvec[k + wvec_dim1] = quart * (tempa * tempa + tempb * tempb);
 704         wvec[k + (wvec_dim1 << 1)] = tempa * tempc;
 705         wvec[k + wvec_dim1 * 3] = tempb * tempc;
 706         wvec[k + (wvec_dim1 << 2)] = quart * (tempa * tempa - tempb * tempb);
 707 /* L110: */
 708         wvec[k + wvec_dim1 * 5] = half * tempa * tempb;
 709     }
 710     i__2 = *n;
 711     for (i__ = 1; i__ <= i__2; ++i__) {
 712         ip = i__ + *npt;
 713         wvec[ip + wvec_dim1] = zero;
 714         wvec[ip + (wvec_dim1 << 1)] = d__[i__];
 715         wvec[ip + wvec_dim1 * 3] = s[i__];
 716         wvec[ip + (wvec_dim1 << 2)] = zero;
 717 /* L120: */
 718         wvec[ip + wvec_dim1 * 5] = zero;
 719     }
 720
 721 /* Put the coefficents of THETA*Wcheck in PROD. */
 722
 723     for (jc = 1; jc <= 5; ++jc) {
 724         nw = *npt;
 725         if (jc == 2 || jc == 3) {
 726             nw = *ndim;
 727         }
 728         i__2 = *npt;
 729         for (k = 1; k <= i__2; ++k) {
 730 /* L130: */
 731             prod[k + jc * prod_dim1] = zero;
 732         }
 733         i__2 = nptm;
 734         for (j = 1; j <= i__2; ++j) {
 735             sum = zero;
 736             i__1 = *npt;
 737             for (k = 1; k <= i__1; ++k) {
 738 /* L140: */
 739                 sum += zmat[k + j * zmat_dim1] * wvec[k + jc * wvec_dim1];
 740             }
 741             if (j < *idz) {
 742                 sum = -sum;
 743             }
 744             i__1 = *npt;
 745             for (k = 1; k <= i__1; ++k) {
 746 /* L150: */
 747                 prod[k + jc * prod_dim1] += sum * zmat[k + j * zmat_dim1];
 748             }
 749         }
 750         if (nw == *ndim) {
 751             i__1 = *npt;
 752             for (k = 1; k <= i__1; ++k) {
 753                 sum = zero;
 754                 i__2 = *n;
 755                 for (j = 1; j <= i__2; ++j) {
 756 /* L160: */
 757                     sum += bmat[k + j * bmat_dim1] * wvec[*npt + j + jc *
 758                             wvec_dim1];
 759                 }
 760 /* L170: */
 761                 prod[k + jc * prod_dim1] += sum;
 762             }
 763         }
 764         i__1 = *n;
 765         for (j = 1; j <= i__1; ++j) {
 766             sum = zero;
 767             i__2 = nw;
 768             for (i__ = 1; i__ <= i__2; ++i__) {
 769 /* L180: */
 770                 sum += bmat[i__ + j * bmat_dim1] * wvec[i__ + jc * wvec_dim1];
 771             }
 772 /* L190: */
 773             prod[*npt + j + jc * prod_dim1] = sum;
 774         }
 775     }
 776
 777 /* Include in DEN the part of BETA that depends on THETA. */
 778
 779     i__1 = *ndim;
 780     for (k = 1; k <= i__1; ++k) {
 781         sum = zero;
 782         for (i__ = 1; i__ <= 5; ++i__) {
 783             par[i__ - 1] = half * prod[k + i__ * prod_dim1] * wvec[k + i__ *
 784                     wvec_dim1];
 785 /* L200: */
 786             sum += par[i__ - 1];
 787         }
 788         den[0] = den[0] - par[0] - sum;
 789         tempa = prod[k + prod_dim1] * wvec[k + (wvec_dim1 << 1)] + prod[k + (
 790                 prod_dim1 << 1)] * wvec[k + wvec_dim1];
 791         tempb = prod[k + (prod_dim1 << 1)] * wvec[k + (wvec_dim1 << 2)] +
 792                 prod[k + (prod_dim1 << 2)] * wvec[k + (wvec_dim1 << 1)];
 793         tempc = prod[k + prod_dim1 * 3] * wvec[k + wvec_dim1 * 5] + prod[k +
 794                 prod_dim1 * 5] * wvec[k + wvec_dim1 * 3];
 795         den[1] = den[1] - tempa - half * (tempb + tempc);
 796         den[5] -= half * (tempb - tempc);
 797         tempa = prod[k + prod_dim1] * wvec[k + wvec_dim1 * 3] + prod[k +
 798                 prod_dim1 * 3] * wvec[k + wvec_dim1];
 799         tempb = prod[k + (prod_dim1 << 1)] * wvec[k + wvec_dim1 * 5] + prod[k
 800                 + prod_dim1 * 5] * wvec[k + (wvec_dim1 << 1)];
 801         tempc = prod[k + prod_dim1 * 3] * wvec[k + (wvec_dim1 << 2)] + prod[k
 802                 + (prod_dim1 << 2)] * wvec[k + wvec_dim1 * 3];
 803         den[2] = den[2] - tempa - half * (tempb - tempc);
 804         den[6] -= half * (tempb + tempc);
 805         tempa = prod[k + prod_dim1] * wvec[k + (wvec_dim1 << 2)] + prod[k + (
 806                 prod_dim1 << 2)] * wvec[k + wvec_dim1];
 807         den[3] = den[3] - tempa - par[1] + par[2];
 808         tempa = prod[k + prod_dim1] * wvec[k + wvec_dim1 * 5] + prod[k +
 809                 prod_dim1 * 5] * wvec[k + wvec_dim1];
 810         tempb = prod[k + (prod_dim1 << 1)] * wvec[k + wvec_dim1 * 3] + prod[k
 811                 + prod_dim1 * 3] * wvec[k + (wvec_dim1 << 1)];
 812         den[4] = den[4] - tempa - half * tempb;
 813         den[7] = den[7] - par[3] + par[4];
 814         tempa = prod[k + (prod_dim1 << 2)] * wvec[k + wvec_dim1 * 5] + prod[k
 815                 + prod_dim1 * 5] * wvec[k + (wvec_dim1 << 2)];
 816 /* L210: */
 817         den[8] -= half * tempa;
 818     }
 819
 820 /* Extend DEN so that it holds all the coefficients of DENOM. */
 821
 822     sum = zero;
 823     for (i__ = 1; i__ <= 5; ++i__) {
 824 /* Computing 2nd power */
 825         d__1 = prod[*knew + i__ * prod_dim1];
 826         par[i__ - 1] = half * (d__1 * d__1);
 827 /* L220: */
 828         sum += par[i__ - 1];
 829     }
 830     denex[0] = alpha * den[0] + par[0] + sum;
 831     tempa = two * prod[*knew + prod_dim1] * prod[*knew + (prod_dim1 << 1)];
 832     tempb = prod[*knew + (prod_dim1 << 1)] * prod[*knew + (prod_dim1 << 2)];
 833     tempc = prod[*knew + prod_dim1 * 3] * prod[*knew + prod_dim1 * 5];
 834     denex[1] = alpha * den[1] + tempa + tempb + tempc;
 835     denex[5] = alpha * den[5] + tempb - tempc;
 836     tempa = two * prod[*knew + prod_dim1] * prod[*knew + prod_dim1 * 3];
 837     tempb = prod[*knew + (prod_dim1 << 1)] * prod[*knew + prod_dim1 * 5];
 838     tempc = prod[*knew + prod_dim1 * 3] * prod[*knew + (prod_dim1 << 2)];
 839     denex[2] = alpha * den[2] + tempa + tempb - tempc;
 840     denex[6] = alpha * den[6] + tempb + tempc;
 841     tempa = two * prod[*knew + prod_dim1] * prod[*knew + (prod_dim1 << 2)];
 842     denex[3] = alpha * den[3] + tempa + par[1] - par[2];
 843     tempa = two * prod[*knew + prod_dim1] * prod[*knew + prod_dim1 * 5];
 844     denex[4] = alpha * den[4] + tempa + prod[*knew + (prod_dim1 << 1)] * prod[
 845             *knew + prod_dim1 * 3];
 846     denex[7] = alpha * den[7] + par[3] - par[4];
 847     denex[8] = alpha * den[8] + prod[*knew + (prod_dim1 << 2)] * prod[*knew +
 848             prod_dim1 * 5];
 849
 850 /* Seek the value of the angle that maximizes the modulus of DENOM. */
 851
 852     sum = denex[0] + denex[1] + denex[3] + denex[5] + denex[7];
 853     denold = sum;
 854     denmax = sum;
 855     isave = 0;
 856     iu = 49;
 857     temp = twopi / (double) (iu + 1);
 858     par[0] = one;
 859     i__1 = iu;
 860     for (i__ = 1; i__ <= i__1; ++i__) {
 861         angle = (double) i__ * temp;
 862         par[1] = cos(angle);
 863         par[2] = sin(angle);
 864         for (j = 4; j <= 8; j += 2) {
 865             par[j - 1] = par[1] * par[j - 3] - par[2] * par[j - 2];
 866 /* L230: */
 867             par[j] = par[1] * par[j - 2] + par[2] * par[j - 3];
 868         }
 869         sumold = sum;
 870         sum = zero;
 871         for (j = 1; j <= 9; ++j) {
 872 /* L240: */
 873             sum += denex[j - 1] * par[j - 1];
 874         }
 875         if (fabs(sum) > fabs(denmax)) {
 876             denmax = sum;
 877             isave = i__;
 878             tempa = sumold;
 879         } else if (i__ == isave + 1) {
 880             tempb = sum;
 881         }
 882 /* L250: */
 883     }
 884     if (isave == 0) {
 885         tempa = sum;
 886     }
 887     if (isave == iu) {
 888         tempb = denold;
 889     }
 890     step = zero;
 891     if (tempa != tempb) {
 892         tempa -= denmax;
 893         tempb -= denmax;
 894         step = half * (tempa - tempb) / (tempa + tempb);
 895     }
 896     angle = temp * ((double) isave + step);
 897
 898 /* Calculate the new parameters of the denominator, the new VLAG vector */
 899 /* and the new D. Then test for convergence. */
 900
 901     par[1] = cos(angle);
 902     par[2] = sin(angle);
 903     for (j = 4; j <= 8; j += 2) {
 904         par[j - 1] = par[1] * par[j - 3] - par[2] * par[j - 2];
 905 /* L260: */
 906         par[j] = par[1] * par[j - 2] + par[2] * par[j - 3];
 907     }
 908     *beta = zero;
 909     denmax = zero;
 910     for (j = 1; j <= 9; ++j) {
 911         *beta += den[j - 1] * par[j - 1];
 912 /* L270: */
 913         denmax += denex[j - 1] * par[j - 1];
 914     }
 915     i__1 = *ndim;
 916     for (k = 1; k <= i__1; ++k) {
 917         vlag[k] = zero;
 918         for (j = 1; j <= 5; ++j) {
 919 /* L280: */
 920             vlag[k] += prod[k + j * prod_dim1] * par[j - 1];
 921         }
 922     }
 923     tau = vlag[*knew];
 924     dd = zero;
 925     tempa = zero;
 926     tempb = zero;
 927     i__1 = *n;
 928     for (i__ = 1; i__ <= i__1; ++i__) {
 929         d__[i__] = par[1] * d__[i__] + par[2] * s[i__];
 930         w[i__] = xopt[i__] + d__[i__];
 931 /* Computing 2nd power */
 932         d__1 = d__[i__];
 933         dd += d__1 * d__1;
 934         tempa += d__[i__] * w[i__];
 935 /* L290: */
 936         tempb += w[i__] * w[i__];
 937     }
 938     if (iterc >= *n) {
 939         goto L340;
 940     }
 941     if (iterc > 1) {
 942         densav = MAX2(densav,denold);
 943     }
 944     if (fabs(denmax) <= fabs(densav) * 1.1) {
 945         goto L340;
 946     }
 947     densav = denmax;
 948
 949 /* Set S to half the gradient of the denominator with respect to D. */
 950 /* Then branch for the next iteration. */
 951
 952     i__1 = *n;
 953     for (i__ = 1; i__ <= i__1; ++i__) {
 954         temp = tempa * xopt[i__] + tempb * d__[i__] - vlag[*npt + i__];
 955 /* L300: */
 956         s[i__] = tau * bmat[*knew + i__ * bmat_dim1] + alpha * temp;
 957     }
 958     i__1 = *npt;
 959     for (k = 1; k <= i__1; ++k) {
 960         sum = zero;
 961         i__2 = *n;
 962         for (j = 1; j <= i__2; ++j) {
 963 /* L310: */
 964             sum += xpt[k + j * xpt_dim1] * w[j];
 965         }
 966         if (nlopt_isinf(tau * w[*n + k]) ||
 967             nlopt_isinf(alpha * vlag[k])) return NLOPT_ROUNDOFF_LIMITED;
 968         temp = (tau * w[*n + k] - alpha * vlag[k]) * sum;
 969         i__2 = *n;
 970         for (i__ = 1; i__ <= i__2; ++i__) {
 971 /* L320: */
 972             s[i__] += temp * xpt[k + i__ * xpt_dim1];
 973         }
 974     }
 975     ss = zero;
 976     ds = zero;
 977     i__2 = *n;
 978     for (i__ = 1; i__ <= i__2; ++i__) {
 979 /* Computing 2nd power */
 980         d__1 = s[i__];
 981         ss += d__1 * d__1;
 982 /* L330: */
 983         ds += d__[i__] * s[i__];
 984     }
 985     ssden = dd * ss - ds * ds;
 986     if (ssden >= dd * 1e-8 * ss) {
 987         goto L70;
 988     }
 989
 990 /* Set the vector W before the RETURN from the subroutine. */
 991
 992 L340:
 993     /* SGJ, 2008: crude hack: truncate d to [lb,ub] bounds if given */
 994     if (lb && ub) {
 995          for (k = 1; k <= *n; ++k) {
 996               if (d__[k] > ub[k-1] - xbase[k-1] - xopt[k])
 997                    d__[k] = ub[k-1] - xbase[k-1] - xopt[k];
 998               else if (d__[k] < lb[k-1] - xbase[k-1] - xopt[k])
 999                    d__[k] = lb[k-1] - xbase[k-1] - xopt[k];
1000          }
1001     }
1002
1003     i__2 = *ndim;
1004     for (k = 1; k <= i__2; ++k) {
1005         w[k] = zero;
1006         for (j = 1; j <= 5; ++j) {
1007 /* L350: */
1008             w[k] += wvec[k + j * wvec_dim1] * par[j - 1];
1009         }
1010     }
1011     vlag[*kopt] += one;
1012     return NLOPT_SUCCESS;
1013 } /* bigden_ */
1014
1015 /*************************************************************************/
1016 /* biglag.f */
1017
1018 typedef struct {
1019      int npt, ndim, iter;
1020      double *hcol, *xpt, *bmat, *xopt;
1021      int flipsign;
1022 } lag_data;
1023
1024 /* the Lagrange function, whose absolute value biglag maximizes */
1025 static double lag(int n, const double *dx, double *grad, void *data)
1026 {
1027      lag_data *d = (lag_data *) data;
1028      int i, j, npt = d->npt, ndim = d->ndim;
1029      const double *hcol = d->hcol, *xpt = d->xpt, *bmat = d->bmat, *xopt = d->xopt;
1030      double val = 0;
1031      for (j = 0; j < n; ++j) {
1032           val += bmat[j * ndim] * (xopt[j] + dx[j]);
1033           if (grad) grad[j] = bmat[j * ndim];
1034      }
1035      for (i = 0; i < npt; ++i) {
1036           double dot = 0;
1037           for (j = 0; j < n; ++j)
1038                dot += xpt[i + j * npt] * (xopt[j] + dx[j]);
1039           val += 0.5 * hcol[i] * (dot * dot);
1040           dot *= hcol[i];
1041           if (grad) for (j = 0; j < n; ++j)
1042                          grad[j] += dot * xpt[i + j * npt];
1043      }
1044      if (d->flipsign) {
1045           val = -val;
1046           if (grad) for (j = 0; j < n; ++j) grad[j] = -grad[j];
1047      }
1048      d->iter++;
1049      return val;
1050 }
1051
1052 static nlopt_result biglag_(int *n, int *npt, double *xopt,
1053                     double *xpt, double *bmat, double *zmat, int *idz,
1054                     int *ndim, int *knew, double *delta, double *d__,
1055                     double *alpha, double *hcol, double *gc, double *gd,
1056                     double *s, double *w,
1057                     const double *xbase, const double *lb, const double *ub)
1058 {
1059     /* System generated locals */
1060     int xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1,
1061             zmat_offset, i__1, i__2;
1062     double d__1;
1063
1064     /* Local variables */
1065     int i__, j, k;
1066     double dd, gg;
1067     int iu;
1068     double sp, ss, cf1, cf2, cf3, cf4, cf5, dhd, cth, one, tau, sth, sum,
1069             half, temp, step;
1070     int nptm;
1071     double zero, angle, scale, denom;
1072     int iterc, isave;
1073     double delsq, tempa, tempb = 0.0, twopi, taubeg, tauold, taumax;
1074
1075
1076 /* N is the number of variables. */
1077 /* NPT is the number of interpolation equations. */
1078 /* XOPT is the best interpolation point so far. */
1079 /* XPT contains the coordinates of the current interpolation points. */
1080 /* BMAT provides the last N columns of H. */
1081 /* ZMAT and IDZ give a factorization of the first NPT by NPT submatrix of H. */
1082 /* NDIM is the first dimension of BMAT and has the value NPT+N. */
1083 /* KNEW is the index of the interpolation point that is going to be moved. */
1084 /* DELTA is the current trust region bound. */
1085 /* D will be set to the step from XOPT to the new point. */
1086 /* ALPHA will be set to the KNEW-th diagonal element of the H matrix. */
1087 /* HCOL, GC, GD, S and W will be used for working space. */
1088
1089 /* The step D is calculated in a way that attempts to maximize the modulus */
1090 /* of LFUNC(XOPT+D), subject to the bound ||D|| .LE. DELTA, where LFUNC is */
1091 /* the KNEW-th Lagrange function. */
1092
1093 /* Set some constants. */
1094
1095     /* Parameter adjustments */
1096     zmat_dim1 = *npt;
1097     zmat_offset = 1 + zmat_dim1;
1098     zmat -= zmat_offset;
1099     xpt_dim1 = *npt;
1100     xpt_offset = 1 + xpt_dim1;
1101     xpt -= xpt_offset;
1102     --xopt;
1103     bmat_dim1 = *ndim;
1104     bmat_offset = 1 + bmat_dim1;
1105     bmat -= bmat_offset;
1106     --d__;
1107     --hcol;
1108     --gc;
1109     --gd;
1110     --s;
1111     --w;
1112
1113     /* Function Body */
1114     half = .5;
1115     one = 1.;
1116     zero = 0.;
1117     twopi = atan(one) * 8.;
1118     delsq = *delta * *delta;
1119     nptm = *npt - *n - 1;
1120
1121 /* Set the first NPT components of HCOL to the leading elements of the */
1122 /* KNEW-th column of H. */
1123
1124     iterc = 0;
1125     i__1 = *npt;
1126     for (k = 1; k <= i__1; ++k) {
1127 /* L10: */
1128         hcol[k] = zero;
1129     }
1130     i__1 = nptm;
1131     for (j = 1; j <= i__1; ++j) {
1132         temp = zmat[*knew + j * zmat_dim1];
1133         if (j < *idz) {
1134             temp = -temp;
1135         }
1136         i__2 = *npt;
1137         for (k = 1; k <= i__2; ++k) {
1138 /* L20: */
1139             hcol[k] += temp * zmat[k + j * zmat_dim1];
1140             if (nlopt_isinf(hcol[k])) return NLOPT_ROUNDOFF_LIMITED;
1141         }
1142     }
1143     *alpha = hcol[*knew];
1144
1145 /* Set the unscaled initial direction D. Form the gradient of LFUNC at */
1146 /* XOPT, and multiply D by the second derivative matrix of LFUNC. */
1147
1148     dd = zero;
1149     i__2 = *n;
1150     for (i__ = 1; i__ <= i__2; ++i__) {
1151         d__[i__] = xpt[*knew + i__ * xpt_dim1] - xopt[i__];
1152         gc[i__] = bmat[*knew + i__ * bmat_dim1];
1153         gd[i__] = zero;
1154 /* L30: */
1155 /* Computing 2nd power */
1156         d__1 = d__[i__];
1157         dd += d__1 * d__1;
1158     }
1159     i__2 = *npt;
1160     for (k = 1; k <= i__2; ++k) {
1161         temp = zero;
1162         sum = zero;
1163         i__1 = *n;
1164         for (j = 1; j <= i__1; ++j) {
1165             temp += xpt[k + j * xpt_dim1] * xopt[j];
1166 /* L40: */
1167             sum += xpt[k + j * xpt_dim1] * d__[j];
1168         }
1169         temp = hcol[k] * temp;
1170         sum = hcol[k] * sum;
1171         i__1 = *n;
1172         for (i__ = 1; i__ <= i__1; ++i__) {
1173             gc[i__] += temp * xpt[k + i__ * xpt_dim1];
1174 /* L50: */
1175             gd[i__] += sum * xpt[k + i__ * xpt_dim1];
1176         }
1177     }
1178
1179 /* Scale D and GD, with a sign change if required. Set S to another */
1180 /* vector in the initial two dimensional subspace. */
1181
1182     gg = zero;
1183     sp = zero;
1184     dhd = zero;
1185     i__1 = *n;
1186     for (i__ = 1; i__ <= i__1; ++i__) {
1187 /* Computing 2nd power */
1188         d__1 = gc[i__];
1189         gg += d__1 * d__1;
1190         sp += d__[i__] * gc[i__];
1191 /* L60: */
1192         dhd += d__[i__] * gd[i__];
1193     }
1194     scale = *delta / sqrt(dd);
1195     if (sp * dhd < zero) {
1196         scale = -scale;
1197     }
1198     temp = zero;
1199     if (sp * sp > dd * .99 * gg) {
1200         temp = one;
1201     }
1202     tau = scale * (fabs(sp) + half * scale * fabs(dhd));
1203     if (gg * delsq < tau * .01 * tau) {
1204         temp = one;
1205     }
1206     i__1 = *n;
1207     for (i__ = 1; i__ <= i__1; ++i__) {
1208         d__[i__] = scale * d__[i__];
1209         gd[i__] = scale * gd[i__];
1210 /* L70: */
1211         s[i__] = gc[i__] + temp * gd[i__];
1212     }
1213
1214     if (lb && ub) {
1215          double minf, *dlb, *dub, *xtol;
1216          lag_data ld;
1217          ld.npt = *npt;
1218          ld.ndim = *ndim;
1219          ld.iter = 0;
1220          ld.hcol = &hcol[1];
1221          ld.xpt = &xpt[1 + 1 * xpt_dim1];
1222          ld.bmat = &bmat[*knew + 1 * bmat_dim1];
1223          ld.xopt = &xopt[1];
1224          ld.flipsign = 0;
1225          dlb = &gc[1]; dub = &gd[1]; xtol = &s[1];
1226          /* make sure rounding errors don't push initial |d| > delta */
1227          for (j = 1; j <= *n; ++j) d__[j] *= 0.99999;
1228          for (j = 0; j < *n; ++j) {
1229               dlb[j] = -(dub[j] = *delta);
1230               if (dlb[j] < lb[j] - xbase[j] - xopt[j+1])
1231                    dlb[j] = lb[j] - xbase[j] - xopt[j+1];
1232               if (dub[j] > ub[j] - xbase[j] - xopt[j+1])
1233                    dub[j] = ub[j] - xbase[j] - xopt[j+1];
1234               if (dlb[j] > 0) dlb[j] = 0;
1235               if (dub[j] < 0) dub[j] = 0;
1236               if (d__[j+1] < dlb[j]) d__[j+1] = dlb[j];
1237               else if (d__[j+1] > dub[j]) d__[j+1] = dub[j];
1238               xtol[j] = 1e-5 * *delta;
1239          }
1240          ld.flipsign = lag(*n, &d__[1], 0, &ld) > 0; /* maximize if > 0 */
1241          return nlopt_minimize_constrained(NLOPT_LD_MMA, *n, lag, &ld,
1242                                     1, rho_constraint, delta, sizeof(double),
1243                                     dlb, dub, &d__[1], &minf, -HUGE_VAL,
1244                                     0., 0., 0., xtol, 1000, 0.);
1245     }
1246
1247 /* Begin the iteration by overwriting S with a vector that has the */
1248 /* required length and direction, except that termination occurs if */
1249 /* the given D and S are nearly parallel. */
1250
1251 L80:
1252     ++iterc;
1253     dd = zero;
1254     sp = zero;
1255     ss = zero;
1256     i__1 = *n;
1257     for (i__ = 1; i__ <= i__1; ++i__) {
1258 /* Computing 2nd power */
1259         d__1 = d__[i__];
1260         dd += d__1 * d__1;
1261         sp += d__[i__] * s[i__];
1262 /* L90: */
1263 /* Computing 2nd power */
1264         d__1 = s[i__];
1265         ss += d__1 * d__1;
1266     }
1267     temp = dd * ss - sp * sp;
1268     if (temp <= dd * 1e-8 * ss) {
1269         goto L160;
1270     }
1271     denom = sqrt(temp);
1272     i__1 = *n;
1273     for (i__ = 1; i__ <= i__1; ++i__) {
1274         s[i__] = (dd * s[i__] - sp * d__[i__]) / denom;
1275 /* L100: */
1276         w[i__] = zero;
1277     }
1278
1279 /* Calculate the coefficients of the objective function on the circle, */
1280 /* beginning with the multiplication of S by the second derivative matrix. */
1281
1282     i__1 = *npt;
1283     for (k = 1; k <= i__1; ++k) {
1284         sum = zero;
1285         i__2 = *n;
1286         for (j = 1; j <= i__2; ++j) {
1287 /* L110: */
1288             sum += xpt[k + j * xpt_dim1] * s[j];
1289         }
1290         sum = hcol[k] * sum;
1291         i__2 = *n;
1292         for (i__ = 1; i__ <= i__2; ++i__) {
1293 /* L120: */
1294             w[i__] += sum * xpt[k + i__ * xpt_dim1];
1295         }
1296     }
1297     cf1 = zero;
1298     cf2 = zero;
1299     cf3 = zero;
1300     cf4 = zero;
1301     cf5 = zero;
1302     i__2 = *n;
1303     for (i__ = 1; i__ <= i__2; ++i__) {
1304         cf1 += s[i__] * w[i__];
1305         cf2 += d__[i__] * gc[i__];
1306         cf3 += s[i__] * gc[i__];
1307         cf4 += d__[i__] * gd[i__];
1308 /* L130: */
1309         cf5 += s[i__] * gd[i__];
1310     }
1311     cf1 = half * cf1;
1312     cf4 = half * cf4 - cf1;
1313
1314 /* Seek the value of the angle that maximizes the modulus of TAU. */
1315
1316     taubeg = cf1 + cf2 + cf4;
1317     taumax = taubeg;
1318     tauold = taubeg;
1319     isave = 0;
1320     iu = 49;
1321     temp = twopi / (double) (iu + 1);
1322     i__2 = iu;
1323     for (i__ = 1; i__ <= i__2; ++i__) {
1324         angle = (double) i__ * temp;
1325         cth = cos(angle);
1326         sth = sin(angle);
1327         tau = cf1 + (cf2 + cf4 * cth) * cth + (cf3 + cf5 * cth) * sth;
1328         if (fabs(tau) > fabs(taumax)) {
1329             taumax = tau;
1330             isave = i__;
1331             tempa = tauold;
1332         } else if (i__ == isave + 1) {
1333             tempb = tau;
1334         }
1335 /* L140: */
1336         tauold = tau;
1337     }
1338     if (isave == 0) {
1339         tempa = tau;
1340     }
1341     if (isave == iu) {
1342         tempb = taubeg;
1343     }
1344     step = zero;
1345     if (tempa != tempb) {
1346         tempa -= taumax;
1347         tempb -= taumax;
1348         step = half * (tempa - tempb) / (tempa + tempb);
1349     }
1350     angle = temp * ((double) isave + step);
1351
1352 /* Calculate the new D and GD. Then test for convergence. */
1353
1354     cth = cos(angle);
1355     sth = sin(angle);
1356     tau = cf1 + (cf2 + cf4 * cth) * cth + (cf3 + cf5 * cth) * sth;
1357     i__2 = *n;
1358     for (i__ = 1; i__ <= i__2; ++i__) {
1359         d__[i__] = cth * d__[i__] + sth * s[i__];
1360         gd[i__] = cth * gd[i__] + sth * w[i__];
1361 /* L150: */
1362         s[i__] = gc[i__] + gd[i__];
1363     }
1364     if (fabs(tau) <= fabs(taubeg) * 1.1) {
1365         goto L160;
1366     }
1367     if (iterc < *n) {
1368         goto L80;
1369     }
1370 L160:
1371     return NLOPT_SUCCESS;
1372 } /* biglag_ */
1373
1374
1375
1376 /*************************************************************************/
1377 /* update.f */
1378
1379 static void update_(int *n, int *npt, double *bmat,
1380         double *zmat, int *idz, int *ndim, double *vlag,
1381         double *beta, int *knew, double *w)
1382 {
1383     /* System generated locals */
1384     int bmat_dim1, bmat_offset, zmat_dim1, zmat_offset, i__1, i__2;
1385     double d__1, d__2;
1386
1387     /* Local variables */
1388     int i__, j, ja, jb, jl, jp;
1389     double one, tau, temp;
1390     int nptm;
1391     double zero;
1392     int iflag;
1393     double scala, scalb_, alpha, denom, tempa, tempb = 0.0, tausq;
1394
1395
1396 /* The arrays BMAT and ZMAT with IDZ are updated, in order to shift the */
1397 /* interpolation point that has index KNEW. On entry, VLAG contains the */
1398 /* components of the vector Theta*Wcheck+e_b of the updating formula */
1399 /* (6.11), and BETA holds the value of the parameter that has this name. */
1400 /* The vector W is used for working space. */
1401
1402 /* Set some constants. */
1403
1404     /* Parameter adjustments */
1405     zmat_dim1 = *npt;
1406     zmat_offset = 1 + zmat_dim1;
1407     zmat -= zmat_offset;
1408     bmat_dim1 = *ndim;
1409     bmat_offset = 1 + bmat_dim1;
1410     bmat -= bmat_offset;
1411     --vlag;
1412     --w;
1413
1414     /* Function Body */
1415     one = 1.;
1416     zero = 0.;
1417     nptm = *npt - *n - 1;
1418
1419 /* Apply the rotations that put zeros in the KNEW-th row of ZMAT. */
1420
1421     jl = 1;
1422     i__1 = nptm;
1423     for (j = 2; j <= i__1; ++j) {
1424         if (j == *idz) {
1425             jl = *idz;
1426         } else if (zmat[*knew + j * zmat_dim1] != zero) {
1427 /* Computing 2nd power */
1428             d__1 = zmat[*knew + jl * zmat_dim1];
1429 /* Computing 2nd power */
1430             d__2 = zmat[*knew + j * zmat_dim1];
1431             temp = sqrt(d__1 * d__1 + d__2 * d__2);
1432             tempa = zmat[*knew + jl * zmat_dim1] / temp;
1433             tempb = zmat[*knew + j * zmat_dim1] / temp;
1434             i__2 = *npt;
1435             for (i__ = 1; i__ <= i__2; ++i__) {
1436                 temp = tempa * zmat[i__ + jl * zmat_dim1] + tempb * zmat[i__
1437                         + j * zmat_dim1];
1438                 zmat[i__ + j * zmat_dim1] = tempa * zmat[i__ + j * zmat_dim1]
1439                         - tempb * zmat[i__ + jl * zmat_dim1];
1440 /* L10: */
1441                 zmat[i__ + jl * zmat_dim1] = temp;
1442             }
1443             zmat[*knew + j * zmat_dim1] = zero;
1444         }
1445 /* L20: */
1446     }
1447
1448 /* Put the first NPT components of the KNEW-th column of HLAG into W, */
1449 /* and calculate the parameters of the updating formula. */
1450
1451     tempa = zmat[*knew + zmat_dim1];
1452     if (*idz >= 2) {
1453         tempa = -tempa;
1454     }
1455     if (jl > 1) {
1456         tempb = zmat[*knew + jl * zmat_dim1];
1457     }
1458     i__1 = *npt;
1459     for (i__ = 1; i__ <= i__1; ++i__) {
1460         w[i__] = tempa * zmat[i__ + zmat_dim1];
1461         if (jl > 1) {
1462             w[i__] += tempb * zmat[i__ + jl * zmat_dim1];
1463         }
1464 /* L30: */
1465     }
1466     alpha = w[*knew];
1467     tau = vlag[*knew];
1468     tausq = tau * tau;
1469     denom = alpha * *beta + tausq;
1470     vlag[*knew] -= one;
1471
1472 /* Complete the updating of ZMAT when there is only one nonzero element */
1473 /* in the KNEW-th row of the new matrix ZMAT, but, if IFLAG is set to one, */
1474 /* then the first column of ZMAT will be exchanged with another one later. */
1475
1476     iflag = 0;
1477     if (jl == 1) {
1478         temp = sqrt((fabs(denom)));
1479         tempb = tempa / temp;
1480         tempa = tau / temp;
1481         i__1 = *npt;
1482         for (i__ = 1; i__ <= i__1; ++i__) {
1483 /* L40: */
1484             zmat[i__ + zmat_dim1] = tempa * zmat[i__ + zmat_dim1] - tempb *
1485                     vlag[i__];
1486         }
1487         if (*idz == 1 && temp < zero) {
1488             *idz = 2;
1489         }
1490         if (*idz >= 2 && temp >= zero) {
1491             iflag = 1;
1492         }
1493     } else {
1494
1495 /* Complete the updating of ZMAT in the alternative case. */
1496
1497         ja = 1;
1498         if (*beta >= zero) {
1499             ja = jl;
1500         }
1501         jb = jl + 1 - ja;
1502         temp = zmat[*knew + jb * zmat_dim1] / denom;
1503         tempa = temp * *beta;
1504         tempb = temp * tau;
1505         temp = zmat[*knew + ja * zmat_dim1];
1506         scala = one / sqrt(fabs(*beta) * temp * temp + tausq);
1507         scalb_ = scala * sqrt((fabs(denom)));
1508         i__1 = *npt;
1509         for (i__ = 1; i__ <= i__1; ++i__) {
1510             zmat[i__ + ja * zmat_dim1] = scala * (tau * zmat[i__ + ja *
1511                     zmat_dim1] - temp * vlag[i__]);
1512 /* L50: */
1513             zmat[i__ + jb * zmat_dim1] = scalb_ * (zmat[i__ + jb * zmat_dim1]
1514                     - tempa * w[i__] - tempb * vlag[i__]);
1515         }
1516         if (denom <= zero) {
1517             if (*beta < zero) {
1518                 ++(*idz);
1519             }
1520             if (*beta >= zero) {
1521                 iflag = 1;
1522             }
1523         }
1524     }
1525
1526 /* IDZ is reduced in the following case, and usually the first column */
1527 /* of ZMAT is exchanged with a later one. */
1528
1529     if (iflag == 1) {
1530         --(*idz);
1531         i__1 = *npt;
1532         for (i__ = 1; i__ <= i__1; ++i__) {
1533             temp = zmat[i__ + zmat_dim1];
1534             zmat[i__ + zmat_dim1] = zmat[i__ + *idz * zmat_dim1];
1535 /* L60: */
1536             zmat[i__ + *idz * zmat_dim1] = temp;
1537         }
1538     }
1539
1540 /* Finally, update the matrix BMAT. */
1541
1542     i__1 = *n;
1543     for (j = 1; j <= i__1; ++j) {
1544         jp = *npt + j;
1545         w[jp] = bmat[*knew + j * bmat_dim1];
1546         tempa = (alpha * vlag[jp] - tau * w[jp]) / denom;
1547         tempb = (-(*beta) * w[jp] - tau * vlag[jp]) / denom;
1548         i__2 = jp;
1549         for (i__ = 1; i__ <= i__2; ++i__) {
1550             bmat[i__ + j * bmat_dim1] = bmat[i__ + j * bmat_dim1] + tempa *
1551                     vlag[i__] + tempb * w[i__];
1552             if (i__ > *npt) {
1553                 bmat[jp + (i__ - *npt) * bmat_dim1] = bmat[i__ + j *
1554                         bmat_dim1];
1555             }
1556 /* L70: */
1557         }
1558     }
1559     return;
1560 } /* update_ */
1561
1562
1563 /*************************************************************************/
1564 /* newuob.f */
1565
1566 static nlopt_result newuob_(int *n, int *npt, double *x,
1567                             double *rhobeg,
1568                             const double *lb, const double *ub,
1569                             nlopt_stopping *stop, double *minf,
1570                             newuoa_func calfun, void *calfun_data,
1571                     double *xbase, double *xopt, double *xnew,
1572                     double *xpt, double *fval, double *gq, double *hq,
1573                     double *pq, double *bmat, double *zmat, int *ndim,
1574                     double *d__, double *vlag, double *w)
1575 {
1576     /* System generated locals */
1577     int xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1,
1578             zmat_offset, i__1, i__2, i__3;
1579     double d__1, d__2, d__3;
1580
1581     /* Local variables */
1582     double f = 0.0;
1583     int i__, j, k, ih, nf, nh, ip, jp;
1584     double dx;
1585     int np, nfm;
1586     double one;
1587     int idz;
1588     double dsq, rho = 0.0;
1589     int ipt = 0, jpt = 0;
1590     double sum, fbeg = 0.0, diff, half, beta;
1591     int nfmm;
1592     double gisq;
1593     int knew;
1594     double temp, suma, sumb, fopt = HUGE_VAL, bsum, gqsq;
1595     int kopt, nptm;
1596     double zero, xipt = 0.0, xjpt = 0.0, sumz, diffa = 0.0, diffb = 0.0, diffc = 0.0, hdiag, alpha = 0.0,
1597             delta, recip, reciq, fsave;
1598     int ksave, nfsav = 0, itemp;
1599     double dnorm = 0.0, ratio = 0.0, dstep, tenth, vquad;
1600     int ktemp;
1601     double tempq;
1602     int itest = 0;
1603     double rhosq;
1604     double detrat, crvmin = 0.0;
1605     double distsq;
1606     double xoptsq = 0.0;
1607     double rhoend;
1608     nlopt_result rc = NLOPT_SUCCESS, rc2;
1609
1610 /* SGJ, 2008: compute rhoend from NLopt stop info */
1611     rhoend = stop->xtol_rel * (*rhobeg);
1612     for (j = 0; j < *n; ++j)
1613          if (rhoend < stop->xtol_abs[j])
1614               rhoend = stop->xtol_abs[j];
1615
1616 /* The arguments N, NPT, X, RHOBEG, RHOEND, IPRINT and MAXFUN are identical */
1617 /*   to the corresponding arguments in SUBROUTINE NEWUOA. */
1618 /* XBASE will hold a shift of origin that should reduce the contributions */
1619 /*   from rounding errors to values of the model and Lagrange functions. */
1620 /* XOPT will be set to the displacement from XBASE of the vector of */
1621 /*   variables that provides the least calculated F so far. */
1622 /* XNEW will be set to the displacement from XBASE of the vector of */
1623 /*   variables for the current calculation of F. */
1624 /* XPT will contain the interpolation point coordinates relative to XBASE. */
1625 /* FVAL will hold the values of F at the interpolation points. */
1626 /* GQ will hold the gradient of the quadratic model at XBASE. */
1627 /* HQ will hold the explicit second derivatives of the quadratic model. */
1628 /* PQ will contain the parameters of the implicit second derivatives of */
1629 /*   the quadratic model. */
1630 /* BMAT will hold the last N columns of H. */
1631 /* ZMAT will hold the factorization of the leading NPT by NPT submatrix of */
1632 /*   H, this factorization being ZMAT times Diag(DZ) times ZMAT^T, where */
1633 /*   the elements of DZ are plus or minus one, as specified by IDZ. */
1634 /* NDIM is the first dimension of BMAT and has the value NPT+N. */
1635 /* D is reserved for trial steps from XOPT. */
1636 /* VLAG will contain the values of the Lagrange functions at a new point X. */
1637 /*   They are part of a product that requires VLAG to be of length NDIM. */
1638 /* The array W will be used for working space. Its length must be at least */
1639 /*   10*NDIM = 10*(NPT+N). */
1640
1641 /* Set some constants. */
1642
1643     /* Parameter adjustments */
1644     zmat_dim1 = *npt;
1645     zmat_offset = 1 + zmat_dim1;
1646     zmat -= zmat_offset;
1647     xpt_dim1 = *npt;
1648     xpt_offset = 1 + xpt_dim1;
1649     xpt -= xpt_offset;
1650     --x;
1651     --xbase;
1652     --xopt;
1653     --xnew;
1654     --fval;
1655     --gq;
1656     --hq;
1657     --pq;
1658     bmat_dim1 = *ndim;
1659     bmat_offset = 1 + bmat_dim1;
1660     bmat -= bmat_offset;
1661     --d__;
1662     --vlag;
1663     --w;
1664
1665     /* Function Body */
1666     half = .5;
1667     one = 1.;
1668     tenth = .1;
1669     zero = 0.;
1670     np = *n + 1;
1671     nh = *n * np / 2;
1672     nptm = *npt - np;
1673
1674 /* Set the initial elements of XPT, BMAT, HQ, PQ and ZMAT to zero. */
1675
1676     i__1 = *n;
1677     for (j = 1; j <= i__1; ++j) {
1678         xbase[j] = x[j];
1679         i__2 = *npt;
1680         for (k = 1; k <= i__2; ++k) {
1681 /* L10: */
1682             xpt[k + j * xpt_dim1] = zero;
1683         }
1684         i__2 = *ndim;
1685         for (i__ = 1; i__ <= i__2; ++i__) {
1686 /* L20: */
1687             bmat[i__ + j * bmat_dim1] = zero;
1688         }
1689     }
1690     i__2 = nh;
1691     for (ih = 1; ih <= i__2; ++ih) {
1692 /* L30: */
1693         hq[ih] = zero;
1694     }
1695     i__2 = *npt;
1696     for (k = 1; k <= i__2; ++k) {
1697         pq[k] = zero;
1698         i__1 = nptm;
1699         for (j = 1; j <= i__1; ++j) {
1700 /* L40: */
1701             zmat[k + j * zmat_dim1] = zero;
1702         }
1703     }
1704
1705 /* Begin the initialization procedure. NF becomes one more than the number */
1706 /* of function values so far. The coordinates of the displacement of the */
1707 /* next initial interpolation point from XBASE are set in XPT(NF,.). */
1708
1709     rhosq = *rhobeg * *rhobeg;
1710     recip = one / rhosq;
1711     reciq = sqrt(half) / rhosq;
1712     nf = 0;
1713 L50:
1714     nfm = nf;
1715     nfmm = nf - *n;
1716     ++nf;
1717     if (nfm <= *n << 1) {
1718         if (nfm >= 1 && nfm <= *n) {
1719             xpt[nf + nfm * xpt_dim1] = *rhobeg;
1720         } else if (nfm > *n) {
1721             xpt[nf + nfmm * xpt_dim1] = -(*rhobeg);
1722         }
1723     } else {
1724         itemp = (nfmm - 1) / *n;
1725         jpt = nfm - itemp * *n - *n;
1726         ipt = jpt + itemp;
1727         if (ipt > *n) {
1728             itemp = jpt;
1729             jpt = ipt - *n;
1730             ipt = itemp;
1731         }
1732         xipt = *rhobeg;
1733         if (fval[ipt + np] < fval[ipt + 1]) {
1734             xipt = -xipt;
1735         }
1736         xjpt = *rhobeg;
1737         if (fval[jpt + np] < fval[jpt + 1]) {
1738             xjpt = -xjpt;
1739         }
1740         xpt[nf + ipt * xpt_dim1] = xipt;
1741         xpt[nf + jpt * xpt_dim1] = xjpt;
1742     }
1743
1744 /* Calculate the next value of F, label 70 being reached immediately */
1745 /* after this calculation. The least function value so far and its index */
1746 /* are required. */
1747
1748     i__1 = *n;
1749     for (j = 1; j <= i__1; ++j) {
1750 /* L60: */
1751         x[j] = xpt[nf + j * xpt_dim1] + xbase[j];
1752     }
1753     if (lb && ub) { /* SGJ, 2008: make sure we are within bounds */
1754          for (j = 1; j <= i__1; ++j) {
1755               if (x[j] < lb[j-1]) x[j] = lb[j-1];
1756               else if (x[j] > ub[j-1]) x[j] = ub[j-1];
1757          }
1758     }
1759     goto L310;
1760 L70:
1761     fval[nf] = f;
1762     if (nf == 1) {
1763         fbeg = f;
1764         fopt = f;
1765         kopt = 1;
1766     } else if (f < fopt) {
1767         fopt = f;
1768         kopt = nf;
1769     }
1770
1771 /* Set the nonzero initial elements of BMAT and the quadratic model in */
1772 /* the cases when NF is at most 2*N+1. */
1773
1774     if (nfm <= *n << 1) {
1775         if (nfm >= 1 && nfm <= *n) {
1776             gq[nfm] = (f - fbeg) / *rhobeg;
1777             if (*npt < nf + *n) {
1778                 bmat[nfm * bmat_dim1 + 1] = -one / *rhobeg;
1779                 bmat[nf + nfm * bmat_dim1] = one / *rhobeg;
1780                 bmat[*npt + nfm + nfm * bmat_dim1] = -half * rhosq;
1781             }
1782         } else if (nfm > *n) {
1783             bmat[nf - *n + nfmm * bmat_dim1] = half / *rhobeg;
1784             bmat[nf + nfmm * bmat_dim1] = -half / *rhobeg;
1785             zmat[nfmm * zmat_dim1 + 1] = -reciq - reciq;
1786             zmat[nf - *n + nfmm * zmat_dim1] = reciq;
1787             zmat[nf + nfmm * zmat_dim1] = reciq;
1788             ih = nfmm * (nfmm + 1) / 2;
1789             temp = (fbeg - f) / *rhobeg;
1790             hq[ih] = (gq[nfmm] - temp) / *rhobeg;
1791             gq[nfmm] = half * (gq[nfmm] + temp);
1792         }
1793
1794 /* Set the off-diagonal second derivatives of the Lagrange functions and */
1795 /* the initial quadratic model. */
1796
1797     } else {
1798         ih = ipt * (ipt - 1) / 2 + jpt;
1799         if (xipt < zero) {
1800             ipt += *n;
1801         }
1802         if (xjpt < zero) {
1803             jpt += *n;
1804         }
1805         zmat[nfmm * zmat_dim1 + 1] = recip;
1806         zmat[nf + nfmm * zmat_dim1] = recip;
1807         zmat[ipt + 1 + nfmm * zmat_dim1] = -recip;
1808         zmat[jpt + 1 + nfmm * zmat_dim1] = -recip;
1809         hq[ih] = (fbeg - fval[ipt + 1] - fval[jpt + 1] + f) / (xipt * xjpt);
1810     }
1811     if (nf < *npt) {
1812         goto L50;
1813     }
1814
1815 /* Begin the iterative procedure, because the initial model is complete. */
1816
1817     rho = *rhobeg;
1818     delta = rho;
1819     idz = 1;
1820     diffa = zero;
1821     diffb = zero;
1822     itest = 0;
1823     xoptsq = zero;
1824     i__1 = *n;
1825     for (i__ = 1; i__ <= i__1; ++i__) {
1826         xopt[i__] = xpt[kopt + i__ * xpt_dim1];
1827 /* L80: */
1828 /* Computing 2nd power */
1829         d__1 = xopt[i__];
1830         xoptsq += d__1 * d__1;
1831     }
1832 L90:
1833     nfsav = nf;
1834
1835 /* Generate the next trust region step and test its length. Set KNEW */
1836 /* to -1 if the purpose of the next F will be to improve the model. */
1837
1838 L100:
1839     knew = 0;
1840     rc2 = trsapp_(n, npt, &xopt[1], &xpt[xpt_offset], &gq[1], &hq[1], &pq[1], &
1841             delta, &d__[1], &w[1], &w[np], &w[np + *n], &w[np + (*n << 1)], &
1842             crvmin, &xbase[1], lb, ub);
1843     if (rc2 < 0) { rc = rc2; goto L530; }
1844     dsq = zero;
1845     i__1 = *n;
1846     for (i__ = 1; i__ <= i__1; ++i__) {
1847 /* L110: */
1848 /* Computing 2nd power */
1849         d__1 = d__[i__];
1850         dsq += d__1 * d__1;
1851     }
1852 /* Computing MIN */
1853     d__1 = delta, d__2 = sqrt(dsq);
1854     dnorm = MIN2(d__1,d__2);
1855     if (dnorm < half * rho) {
1856         knew = -1;
1857         delta = tenth * delta;
1858         ratio = -1.;
1859         if (delta <= rho * 1.5) {
1860             delta = rho;
1861         }
1862         if (nf <= nfsav + 2) {
1863             goto L460;
1864         }
1865         temp = crvmin * .125 * rho * rho;
1866 /* Computing MAX */
1867         d__1 = MAX2(diffa,diffb);
1868         if (temp <= MAX2(d__1,diffc)) {
1869             goto L460;
1870         }
1871         goto L490;
1872     }
1873
1874 /* Shift XBASE if XOPT may be too far from XBASE. First make the changes */
1875 /* to BMAT that do not depend on ZMAT. */
1876
1877 L120:
1878     if (dsq <= xoptsq * .001) {
1879         tempq = xoptsq * .25;
1880         i__1 = *npt;
1881         for (k = 1; k <= i__1; ++k) {
1882             sum = zero;
1883             i__2 = *n;
1884             for (i__ = 1; i__ <= i__2; ++i__) {
1885 /* L130: */
1886                 sum += xpt[k + i__ * xpt_dim1] * xopt[i__];
1887             }
1888             temp = pq[k] * sum;
1889             sum -= half * xoptsq;
1890             w[*npt + k] = sum;
1891             i__2 = *n;
1892             for (i__ = 1; i__ <= i__2; ++i__) {
1893                 gq[i__] += temp * xpt[k + i__ * xpt_dim1];
1894                 xpt[k + i__ * xpt_dim1] -= half * xopt[i__];
1895                 vlag[i__] = bmat[k + i__ * bmat_dim1];
1896                 w[i__] = sum * xpt[k + i__ * xpt_dim1] + tempq * xopt[i__];
1897                 ip = *npt + i__;
1898                 i__3 = i__;
1899                 for (j = 1; j <= i__3; ++j) {
1900 /* L140: */
1901                     bmat[ip + j * bmat_dim1] = bmat[ip + j * bmat_dim1] +
1902                             vlag[i__] * w[j] + w[i__] * vlag[j];
1903                 }
1904             }
1905         }
1906
1907 /* Then the revisions of BMAT that depend on ZMAT are calculated. */
1908
1909         i__3 = nptm;
1910         for (k = 1; k <= i__3; ++k) {
1911             sumz = zero;
1912             i__2 = *npt;
1913             for (i__ = 1; i__ <= i__2; ++i__) {
1914                 sumz += zmat[i__ + k * zmat_dim1];
1915 /* L150: */
1916                 w[i__] = w[*npt + i__] * zmat[i__ + k * zmat_dim1];
1917             }
1918             i__2 = *n;
1919             for (j = 1; j <= i__2; ++j) {
1920                 sum = tempq * sumz * xopt[j];
1921                 i__1 = *npt;
1922                 for (i__ = 1; i__ <= i__1; ++i__) {
1923 /* L160: */
1924                     sum += w[i__] * xpt[i__ + j * xpt_dim1];
1925                 }
1926                 vlag[j] = sum;
1927                 if (k < idz) {
1928                     sum = -sum;
1929                 }
1930                 i__1 = *npt;
1931                 for (i__ = 1; i__ <= i__1; ++i__) {
1932 /* L170: */
1933                     bmat[i__ + j * bmat_dim1] += sum * zmat[i__ + k *
1934                             zmat_dim1];
1935                 }
1936             }
1937             i__1 = *n;
1938             for (i__ = 1; i__ <= i__1; ++i__) {
1939                 ip = i__ + *npt;
1940                 temp = vlag[i__];
1941                 if (k < idz) {
1942                     temp = -temp;
1943                 }
1944                 i__2 = i__;
1945                 for (j = 1; j <= i__2; ++j) {
1946 /* L180: */
1947                     bmat[ip + j * bmat_dim1] += temp * vlag[j];
1948                 }
1949             }
1950         }
1951
1952 /* The following instructions complete the shift of XBASE, including */
1953 /* the changes to the parameters of the quadratic model. */
1954
1955         ih = 0;
1956         i__2 = *n;
1957         for (j = 1; j <= i__2; ++j) {
1958             w[j] = zero;
1959             i__1 = *npt;
1960             for (k = 1; k <= i__1; ++k) {
1961                 w[j] += pq[k] * xpt[k + j * xpt_dim1];
1962 /* L190: */
1963                 xpt[k + j * xpt_dim1] -= half * xopt[j];
1964             }
1965             i__1 = j;
1966             for (i__ = 1; i__ <= i__1; ++i__) {
1967                 ++ih;
1968                 if (i__ < j) {
1969                     gq[j] += hq[ih] * xopt[i__];
1970                 }
1971                 gq[i__] += hq[ih] * xopt[j];
1972                 hq[ih] = hq[ih] + w[i__] * xopt[j] + xopt[i__] * w[j];
1973 /* L200: */
1974                 bmat[*npt + i__ + j * bmat_dim1] = bmat[*npt + j + i__ *
1975                         bmat_dim1];
1976             }
1977         }
1978         i__1 = *n;
1979         for (j = 1; j <= i__1; ++j) {
1980             xbase[j] += xopt[j];
1981 /* L210: */
1982             xopt[j] = zero;
1983         }
1984         xoptsq = zero;
1985     }
1986
1987 /* Pick the model step if KNEW is positive. A different choice of D */
1988 /* may be made later, if the choice of D by BIGLAG causes substantial */
1989 /* cancellation in DENOM. */
1990
1991     if (knew > 0) {
1992        rc2 =
1993          biglag_(n, npt, &xopt[1], &xpt[xpt_offset], &bmat[bmat_offset], &zmat[
1994                 zmat_offset], &idz, ndim, &knew, &dstep, &d__[1], &alpha, &
1995                 vlag[1], &vlag[*npt + 1], &w[1], &w[np], &w[np + *n],
1996                 &xbase[1], lb, ub);
1997        if (rc2 < 0) { rc = rc2; goto L530; }
1998     }
1999
2000 /* Calculate VLAG and BETA for the current choice of D. The first NPT */
2001 /* components of W_check will be held in W. */
2002
2003     i__1 = *npt;
2004     for (k = 1; k <= i__1; ++k) {
2005         suma = zero;
2006         sumb = zero;
2007         sum = zero;
2008         i__2 = *n;
2009         for (j = 1; j <= i__2; ++j) {
2010             suma += xpt[k + j * xpt_dim1] * d__[j];
2011             sumb += xpt[k + j * xpt_dim1] * xopt[j];
2012 /* L220: */
2013             sum += bmat[k + j * bmat_dim1] * d__[j];
2014         }
2015         w[k] = suma * (half * suma + sumb);
2016 /* L230: */
2017         vlag[k] = sum;
2018     }
2019     beta = zero;
2020     i__1 = nptm;
2021     for (k = 1; k <= i__1; ++k) {
2022         sum = zero;
2023         i__2 = *npt;
2024         for (i__ = 1; i__ <= i__2; ++i__) {
2025 /* L240: */
2026             sum += zmat[i__ + k * zmat_dim1] * w[i__];
2027         }
2028         if (k < idz) {
2029             beta += sum * sum;
2030             sum = -sum;
2031         } else {
2032             beta -= sum * sum;
2033         }
2034         i__2 = *npt;
2035         for (i__ = 1; i__ <= i__2; ++i__) {
2036 /* L250: */
2037             vlag[i__] += sum * zmat[i__ + k * zmat_dim1];
2038         }
2039     }
2040     bsum = zero;
2041     dx = zero;
2042     i__2 = *n;
2043     for (j = 1; j <= i__2; ++j) {
2044         sum = zero;
2045         i__1 = *npt;
2046         for (i__ = 1; i__ <= i__1; ++i__) {
2047 /* L260: */
2048             sum += w[i__] * bmat[i__ + j * bmat_dim1];
2049         }
2050         bsum += sum * d__[j];
2051         jp = *npt + j;
2052         i__1 = *n;
2053         for (k = 1; k <= i__1; ++k) {
2054 /* L270: */
2055             sum += bmat[jp + k * bmat_dim1] * d__[k];
2056         }
2057         vlag[jp] = sum;
2058         bsum += sum * d__[j];
2059 /* L280: */
2060         dx += d__[j] * xopt[j];
2061     }
2062     beta = dx * dx + dsq * (xoptsq + dx + dx + half * dsq) + beta - bsum;
2063     vlag[kopt] += one;
2064
2065 /* If KNEW is positive and if the cancellation in DENOM is unacceptable, */
2066 /* then BIGDEN calculates an alternative model step, XNEW being used for */
2067 /* working space. */
2068
2069     if (knew > 0) {
2070 /* Computing 2nd power */
2071         d__1 = vlag[knew];
2072         if (d__1 == 0) { rc = NLOPT_ROUNDOFF_LIMITED; goto L530; }
2073         temp = one + alpha * beta / (d__1 * d__1);
2074         if (fabs(temp) <= .8) {
2075           rc2 =
2076             bigden_(n, npt, &xopt[1], &xpt[xpt_offset], &bmat[bmat_offset], &
2077                     zmat[zmat_offset], &idz, ndim, &kopt, &knew, &d__[1], &w[
2078                     1], &vlag[1], &beta, &xnew[1], &w[*ndim + 1], &w[*ndim *
2079                     6 + 1], &xbase[1], lb, ub);
2080           if (rc2 < 0) { rc = rc2; goto L530; }
2081         }
2082     }
2083
2084 /* Calculate the next value of the objective function. */
2085
2086 L290:
2087     i__2 = *n;
2088     for (i__ = 1; i__ <= i__2; ++i__) {
2089         xnew[i__] = xopt[i__] + d__[i__];
2090 /* L300: */
2091         x[i__] = xbase[i__] + xnew[i__];
2092     }
2093     if (lb && ub) { /* SGJ, 2008: make sure we are within bounds,
2094                        since roundoff errors can push us slightly outside */
2095          for (j = 1; j <= i__1; ++j) {
2096               if (x[j] < lb[j-1]) x[j] = lb[j-1];
2097               else if (x[j] > ub[j-1]) x[j] = ub[j-1];
2098          }
2099     }
2100     ++nf;
2101 L310:
2102     if (nlopt_stop_forced(stop)) rc = NLOPT_FORCED_STOP;
2103     else if (*(stop->nevals_p) > 0) {
2104          if (nlopt_stop_evals(stop)) rc = NLOPT_MAXEVAL_REACHED;
2105          else if (nlopt_stop_time(stop)) rc = NLOPT_MAXTIME_REACHED;
2106     }
2107     if (rc != NLOPT_SUCCESS) goto L530;
2108
2109     ++ *(stop->nevals_p);
2110     f = calfun(*n, &x[1], calfun_data);
2111     if (f < stop->minf_max) {
2112        rc = NLOPT_MINF_MAX_REACHED;
2113        goto L530;
2114     }
2115
2116     if (nf <= *npt) {
2117         goto L70;
2118     }
2119     if (knew == -1) {
2120         goto L530;
2121     }
2122
2123 /* Use the quadratic model to predict the change in F due to the step D, */
2124 /* and set DIFF to the error of this prediction. */
2125
2126     vquad = zero;
2127     ih = 0;
2128     i__2 = *n;
2129     for (j = 1; j <= i__2; ++j) {
2130         vquad += d__[j] * gq[j];
2131         i__1 = j;
2132         for (i__ = 1; i__ <= i__1; ++i__) {
2133             ++ih;
2134             temp = d__[i__] * xnew[j] + d__[j] * xopt[i__];
2135             if (i__ == j) {
2136                 temp = half * temp;
2137             }
2138 /* L340: */
2139             vquad += temp * hq[ih];
2140         }
2141     }
2142     i__1 = *npt;
2143     for (k = 1; k <= i__1; ++k) {
2144 /* L350: */
2145         vquad += pq[k] * w[k];
2146     }
2147     diff = f - fopt - vquad;
2148     diffc = diffb;
2149     diffb = diffa;
2150     diffa = fabs(diff);
2151     if (dnorm > rho) {
2152         nfsav = nf;
2153     }
2154
2155 /* Update FOPT and XOPT if the new F is the least value of the objective */
2156 /* function so far. The branch when KNEW is positive occurs if D is not */
2157 /* a trust region step. */
2158
2159     fsave = fopt;
2160     if (f < fopt) {
2161         fopt = f;
2162         xoptsq = zero;
2163         i__1 = *n;
2164         for (i__ = 1; i__ <= i__1; ++i__) {
2165             xopt[i__] = xnew[i__];
2166 /* L360: */
2167 /* Computing 2nd power */
2168             d__1 = xopt[i__];
2169             xoptsq += d__1 * d__1;
2170         }
2171         if (nlopt_stop_ftol(stop, fopt, fsave)) {
2172              rc = NLOPT_FTOL_REACHED;
2173              goto L530;
2174         }
2175
2176     }
2177     ksave = knew;
2178     if (knew > 0) {
2179         goto L410;
2180     }
2181
2182 /* Pick the next value of DELTA after a trust region step. */
2183
2184     if (vquad >= zero) {
2185         goto L530;
2186     }
2187     ratio = (f - fsave) / vquad;
2188     if (ratio <= tenth) {
2189         delta = half * dnorm;
2190     } else if (ratio <= .7) {
2191 /* Computing MAX */
2192         d__1 = half * delta;
2193         delta = MAX2(d__1,dnorm);
2194     } else {
2195 /* Computing MAX */
2196         d__1 = half * delta, d__2 = dnorm + dnorm;
2197         delta = MAX2(d__1,d__2);
2198     }
2199     if (delta <= rho * 1.5) {
2200         delta = rho;
2201     }
2202
2203 /* Set KNEW to the index of the next interpolation point to be deleted. */
2204
2205 /* Computing MAX */
2206     d__2 = tenth * delta;
2207 /* Computing 2nd power */
2208     d__1 = MAX2(d__2,rho);
2209     rhosq = d__1 * d__1;
2210     ktemp = 0;
2211     detrat = zero;
2212     if (f >= fsave) {
2213         ktemp = kopt;
2214         detrat = one;
2215     }
2216     i__1 = *npt;
2217     for (k = 1; k <= i__1; ++k) {
2218         hdiag = zero;
2219         i__2 = nptm;
2220         for (j = 1; j <= i__2; ++j) {
2221             temp = one;
2222             if (j < idz) {
2223                 temp = -one;
2224             }
2225 /* L380: */
2226 /* Computing 2nd power */
2227             d__1 = zmat[k + j * zmat_dim1];
2228             hdiag += temp * (d__1 * d__1);
2229         }
2230 /* Computing 2nd power */
2231         d__2 = vlag[k];
2232         temp = (d__1 = beta * hdiag + d__2 * d__2, fabs(d__1));
2233         distsq = zero;
2234         i__2 = *n;
2235         for (j = 1; j <= i__2; ++j) {
2236 /* L390: */
2237 /* Computing 2nd power */
2238             d__1 = xpt[k + j * xpt_dim1] - xopt[j];
2239             distsq += d__1 * d__1;
2240         }
2241         if (distsq > rhosq) {
2242 /* Computing 3rd power */
2243             d__1 = distsq / rhosq;
2244             temp *= d__1 * (d__1 * d__1);
2245         }
2246         if (temp > detrat && k != ktemp) {
2247             detrat = temp;
2248             knew = k;
2249         }
2250 /* L400: */
2251     }
2252     if (knew == 0) {
2253         goto L460;
2254     }
2255
2256 /* Update BMAT, ZMAT and IDZ, so that the KNEW-th interpolation point */
2257 /* can be moved. Begin the updating of the quadratic model, starting */
2258 /* with the explicit second derivative term. */
2259
2260 L410:
2261     update_(n, npt, &bmat[bmat_offset], &zmat[zmat_offset], &idz, ndim, &vlag[
2262             1], &beta, &knew, &w[1]);
2263     fval[knew] = f;
2264     ih = 0;
2265     i__1 = *n;
2266     for (i__ = 1; i__ <= i__1; ++i__) {
2267         temp = pq[knew] * xpt[knew + i__ * xpt_dim1];
2268         i__2 = i__;
2269         for (j = 1; j <= i__2; ++j) {
2270             ++ih;
2271 /* L420: */
2272             hq[ih] += temp * xpt[knew + j * xpt_dim1];
2273         }
2274     }
2275     pq[knew] = zero;
2276
2277 /* Update the other second derivative parameters, and then the gradient */
2278 /* vector of the model. Also include the new interpolation point. */
2279
2280     i__2 = nptm;
2281     for (j = 1; j <= i__2; ++j) {
2282         temp = diff * zmat[knew + j * zmat_dim1];
2283         if (j < idz) {
2284             temp = -temp;
2285         }
2286         i__1 = *npt;
2287         for (k = 1; k <= i__1; ++k) {
2288 /* L440: */
2289             pq[k] += temp * zmat[k + j * zmat_dim1];
2290         }
2291     }
2292     gqsq = zero;
2293     i__1 = *n;
2294     for (i__ = 1; i__ <= i__1; ++i__) {
2295         gq[i__] += diff * bmat[knew + i__ * bmat_dim1];
2296 /* Computing 2nd power */
2297         d__1 = gq[i__];
2298         gqsq += d__1 * d__1;
2299 /* L450: */
2300         xpt[knew + i__ * xpt_dim1] = xnew[i__];
2301     }
2302
2303 /* If a trust region step makes a small change to the objective function, */
2304 /* then calculate the gradient of the least Frobenius norm interpolant at */
2305 /* XBASE, and store it in W, using VLAG for a vector of right hand sides. */
2306
2307     if (ksave == 0 && delta == rho) {
2308         if (fabs(ratio) > .01) {
2309             itest = 0;
2310         } else {
2311             i__1 = *npt;
2312             for (k = 1; k <= i__1; ++k) {
2313 /* L700: */
2314                 vlag[k] = fval[k] - fval[kopt];
2315             }
2316             gisq = zero;
2317             i__1 = *n;
2318             for (i__ = 1; i__ <= i__1; ++i__) {
2319                 sum = zero;
2320                 i__2 = *npt;
2321                 for (k = 1; k <= i__2; ++k) {
2322 /* L710: */
2323                     sum += bmat[k + i__ * bmat_dim1] * vlag[k];
2324                 }
2325                 gisq += sum * sum;
2326 /* L720: */
2327                 w[i__] = sum;
2328             }
2329
2330 /* Test whether to replace the new quadratic model by the least Frobenius */
2331 /* norm interpolant, making the replacement if the test is satisfied. */
2332
2333             ++itest;
2334             if (gqsq < gisq * 100.) {
2335                 itest = 0;
2336             }
2337             if (itest >= 3) {
2338                 i__1 = *n;
2339                 for (i__ = 1; i__ <= i__1; ++i__) {
2340 /* L730: */
2341                     gq[i__] = w[i__];
2342                 }
2343                 i__1 = nh;
2344                 for (ih = 1; ih <= i__1; ++ih) {
2345 /* L740: */
2346                     hq[ih] = zero;
2347                 }
2348                 i__1 = nptm;
2349                 for (j = 1; j <= i__1; ++j) {
2350                     w[j] = zero;
2351                     i__2 = *npt;
2352                     for (k = 1; k <= i__2; ++k) {
2353 /* L750: */
2354                         w[j] += vlag[k] * zmat[k + j * zmat_dim1];
2355                     }
2356 /* L760: */
2357                     if (j < idz) {
2358                         w[j] = -w[j];
2359                     }
2360                 }
2361                 i__1 = *npt;
2362                 for (k = 1; k <= i__1; ++k) {
2363                     pq[k] = zero;
2364                     i__2 = nptm;
2365                     for (j = 1; j <= i__2; ++j) {
2366 /* L770: */
2367                         pq[k] += zmat[k + j * zmat_dim1] * w[j];
2368                     }
2369                 }
2370                 itest = 0;
2371             }
2372         }
2373     }
2374     if (f < fsave) {
2375         kopt = knew;
2376     }
2377
2378 /* If a trust region step has provided a sufficient decrease in F, then */
2379 /* branch for another trust region calculation. The case KSAVE>0 occurs */
2380 /* when the new function value was calculated by a model step. */
2381
2382     if (f <= fsave + tenth * vquad) {
2383         goto L100;
2384     }
2385     if (ksave > 0) {
2386         goto L100;
2387     }
2388
2389 /* Alternatively, find out if the interpolation points are close enough */
2390 /* to the best point so far. */
2391
2392     knew = 0;
2393 L460:
2394     distsq = delta * 4. * delta;
2395     i__2 = *npt;
2396     for (k = 1; k <= i__2; ++k) {
2397         sum = zero;
2398         i__1 = *n;
2399         for (j = 1; j <= i__1; ++j) {
2400 /* L470: */
2401 /* Computing 2nd power */
2402             d__1 = xpt[k + j * xpt_dim1] - xopt[j];
2403             sum += d__1 * d__1;
2404         }
2405         if (sum > distsq) {
2406             knew = k;
2407             distsq = sum;
2408         }
2409 /* L480: */
2410     }
2411
2412 /* If KNEW is positive, then set DSTEP, and branch back for the next */
2413 /* iteration, which will generate a "model step". */
2414
2415     if (knew > 0) {
2416 /* Computing MAX */
2417 /* Computing MIN */
2418         d__2 = tenth * sqrt(distsq), d__3 = half * delta;
2419         d__1 = MIN2(d__2,d__3);
2420         dstep = MAX2(d__1,rho);
2421         dsq = dstep * dstep;
2422         goto L120;
2423     }
2424     if (ratio > zero) {
2425         goto L100;
2426     }
2427     if (MAX2(delta,dnorm) > rho) {
2428         goto L100;
2429     }
2430
2431 /* The calculations with the current value of RHO are complete. Pick the */
2432 /* next values of RHO and DELTA. */
2433
2434 L490:
2435     if (rho > rhoend) {
2436         delta = half * rho;
2437         ratio = rho / rhoend;
2438         if (ratio <= 16.) {
2439             rho = rhoend;
2440         } else if (ratio <= 250.) {
2441             rho = sqrt(ratio) * rhoend;
2442         } else {
2443             rho = tenth * rho;
2444         }
2445         delta = MAX2(delta,rho);
2446         goto L90;
2447     }
2448
2449 /* Return from the calculation, after another Newton-Raphson step, if */
2450 /* it is too short to have been tried before. */
2451
2452     if (knew == -1) {
2453         goto L290;
2454     }
2455     rc = NLOPT_XTOL_REACHED;
2456 L530:
2457     if (fopt <= f) {
2458         i__2 = *n;
2459         for (i__ = 1; i__ <= i__2; ++i__) {
2460 /* L540: */
2461             x[i__] = xbase[i__] + xopt[i__];
2462         }
2463         f = fopt;
2464     }
2465     *minf = f;
2466     return rc;
2467 } /* newuob_ */
2468
2469 /*************************************************************************/
2470 /* newuoa.f */
2471
2472 nlopt_result newuoa(int n, int npt, double *x,
2473                     const double *lb, const double *ub,
2474                     double rhobeg, nlopt_stopping *stop, double *minf,
2475                     newuoa_func calfun, void *calfun_data)
2476 {
2477     /* Local variables */
2478     int id, np, iw, igq, ihq, ixb, ifv, ipq, ivl, ixn, ixo, ixp, ndim,
2479             nptm, ibmat, izmat;
2480     nlopt_result ret;
2481     double *w;
2482
2483 /* This subroutine seeks the least value of a function of many variables, */
2484 /* by a trust region method that forms quadratic models by interpolation. */
2485 /* There can be some freedom in the interpolation conditions, which is */
2486 /* taken up by minimizing the Frobenius norm of the change to the second */
2487 /* derivative of the quadratic model, beginning with a zero matrix. The */
2488 /* arguments of the subroutine are as follows. */
2489
2490 /* N must be set to the number of variables and must be at least two. */
2491 /* NPT is the number of interpolation conditions. Its value must be in the */
2492 /*   interval [N+2,(N+1)(N+2)/2]. */
2493 /* Initial values of the variables must be set in X(1),X(2),...,X(N). They */
2494 /*   will be changed to the values that give the least calculated F. */
2495 /* RHOBEG and RHOEND must be set to the initial and final values of a trust */
2496 /*   region radius, so both must be positive with RHOEND<=RHOBEG. Typically */
2497 /*   RHOBEG should be about one tenth of the greatest expected change to a */
2498 /*   variable, and RHOEND should indicate the accuracy that is required in */
2499 /*   the final values of the variables. */
2500 /* The value of IPRINT should be set to 0, 1, 2 or 3, which controls the */
2501 /*   amount of printing. Specifically, there is no output if IPRINT=0 and */
2502 /*   there is output only at the return if IPRINT=1. Otherwise, each new */
2503 /*   value of RHO is printed, with the best vector of variables so far and */
2504 /*   the corresponding value of the objective function. Further, each new */
2505 /*   value of F with its variables are output if IPRINT=3. */
2506 /* MAXFUN must be set to an upper bound on the number of calls of CALFUN. */
2507 /* The array W will be used for working space. Its length must be at least */
2508 /* (NPT+13)*(NPT+N)+3*N*(N+3)/2. */
2509
2510 /* SUBROUTINE CALFUN (N,X,F) must be provided by the user. It must set F to */
2511 /* the value of the objective function for the variables X(1),X(2),...,X(N). */
2512
2513 /* Partition the working space array, so that different parts of it can be */
2514 /* treated separately by the subroutine that performs the main calculation. */
2515
2516     /* Parameter adjustments */
2517     --x;
2518
2519     /* Function Body */
2520     np = n + 1;
2521     nptm = npt - np;
2522     if (n < 2) {
2523         nlopt_stop_msg(stop, "dimension %d must be >= 2", n);
2524         return NLOPT_INVALID_ARGS;
2525     }
2526     if (npt < n + 2 || npt > (n + 2) * np / 2) {
2527         nlopt_stop_msg(stop, "invalid # of interpolation conditions %d", npt);
2528         return NLOPT_INVALID_ARGS;
2529     }
2530     ndim = npt + n;
2531     ixb = 1;
2532     ixo = ixb + n;
2533     ixn = ixo + n;
2534     ixp = ixn + n;
2535     ifv = ixp + n * npt;
2536     igq = ifv + npt;
2537     ihq = igq + n;
2538     ipq = ihq + n * np / 2;
2539     ibmat = ipq + npt;
2540     izmat = ibmat + ndim * n;
2541     id = izmat + npt * nptm;
2542     ivl = id + n;
2543     iw = ivl + ndim;
2544
2545     w = (double *) malloc(sizeof(double) * ((npt+13)*(npt+n) + 3*(n*(n+3))/2));
2546     if (!w) return NLOPT_OUT_OF_MEMORY;
2547     --w;
2548
2549 /* The above settings provide a partition of W for subroutine NEWUOB. */
2550 /* The partition requires the first NPT*(NPT+N)+5*N*(N+3)/2 elements of */
2551 /* W plus the space that is needed by the last array of NEWUOB. */
2552
2553     ret = newuob_(&n, &npt, &x[1], &rhobeg,
2554                   lb, ub, stop, minf, calfun, calfun_data,
2555                   &w[ixb], &w[ixo], &w[ixn], &w[ixp], &w[ifv],
2556                   &w[igq], &w[ihq], &w[ipq], &w[ibmat], &w[izmat],
2557                   &ndim, &w[id], &w[ivl], &w[iw]);
2558
2559     ++w;
2560     free(w);
2561     return ret;
2562 } /* newuoa_ */
2563
2564 /*************************************************************************/
2565 /*************************************************************************/