newuoa/newuoa.c

   1 /* Copyright (c) 2004 M. J. D. Powell (mjdp@cam.ac.uk)
   2  * Copyright (c) 2007-2008 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 min(a,b) ((a) <= (b) ? (a) : (b))
  36 #define max(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, cf, dg, gg;
 112     int ih;
 113     double ds, sg;
 114     int iu;
 115     double ss, dhd, dhs, cth, sgk, shs, sth, qadd, half, qbeg, qred, qmin,
 116              temp, qsav, qnew, zero, ggbeg, alpha, angle, reduc;
 117     int iterc;
 118     double ggsav, delsq, tempa, tempb;
 119     int isave;
 120     double bstep, 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 = min(*crvmin,temp);
 271 /* Computing MIN */
 272         d__1 = alpha, d__2 = gg / dhd;
 273         alpha = min(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 void 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 (dstemp * dstemp / sstemp < dtest) {
 638                     ksav = k;
 639                     dtest = dstemp * dstemp / sstemp;
 640                     ds = dstemp;
 641                     ss = sstemp;
 642                 }
 643             }
 644 /* L50: */
 645         }
 646         i__2 = *n;
 647         for (i__ = 1; i__ <= i__2; ++i__) {
 648 /* L60: */
 649             s[i__] = xpt[ksav + i__ * xpt_dim1] - xopt[i__];
 650         }
 651     }
 652     ssden = dd * ss - ds * ds;
 653     iterc = 0;
 654     densav = zero;
 655
 656 /* Begin the iteration by overwriting S with a vector that has the */
 657 /* required length and direction. */
 658
 659 L70:
 660     ++iterc;
 661     temp = one / sqrt(ssden);
 662     xoptd = zero;
 663     xopts = zero;
 664     i__2 = *n;
 665     for (i__ = 1; i__ <= i__2; ++i__) {
 666         s[i__] = temp * (dd * s[i__] - ds * d__[i__]);
 667         xoptd += xopt[i__] * d__[i__];
 668 /* L80: */
 669         xopts += xopt[i__] * s[i__];
 670     }
 671
 672 /* Set the coefficients of the first two terms of BETA. */
 673
 674     tempa = half * xoptd * xoptd;
 675     tempb = half * xopts * xopts;
 676     den[0] = dd * (xoptsq + half * dd) + tempa + tempb;
 677     den[1] = two * xoptd * dd;
 678     den[2] = two * xopts * dd;
 679     den[3] = tempa - tempb;
 680     den[4] = xoptd * xopts;
 681     for (i__ = 6; i__ <= 9; ++i__) {
 682 /* L90: */
 683         den[i__ - 1] = zero;
 684     }
 685
 686 /* Put the coefficients of Wcheck in WVEC. */
 687
 688     i__2 = *npt;
 689     for (k = 1; k <= i__2; ++k) {
 690         tempa = zero;
 691         tempb = zero;
 692         tempc = zero;
 693         i__1 = *n;
 694         for (i__ = 1; i__ <= i__1; ++i__) {
 695             tempa += xpt[k + i__ * xpt_dim1] * d__[i__];
 696             tempb += xpt[k + i__ * xpt_dim1] * s[i__];
 697 /* L100: */
 698             tempc += xpt[k + i__ * xpt_dim1] * xopt[i__];
 699         }
 700         wvec[k + wvec_dim1] = quart * (tempa * tempa + tempb * tempb);
 701         wvec[k + (wvec_dim1 << 1)] = tempa * tempc;
 702         wvec[k + wvec_dim1 * 3] = tempb * tempc;
 703         wvec[k + (wvec_dim1 << 2)] = quart * (tempa * tempa - tempb * tempb);
 704 /* L110: */
 705         wvec[k + wvec_dim1 * 5] = half * tempa * tempb;
 706     }
 707     i__2 = *n;
 708     for (i__ = 1; i__ <= i__2; ++i__) {
 709         ip = i__ + *npt;
 710         wvec[ip + wvec_dim1] = zero;
 711         wvec[ip + (wvec_dim1 << 1)] = d__[i__];
 712         wvec[ip + wvec_dim1 * 3] = s[i__];
 713         wvec[ip + (wvec_dim1 << 2)] = zero;
 714 /* L120: */
 715         wvec[ip + wvec_dim1 * 5] = zero;
 716     }
 717
 718 /* Put the coefficents of THETA*Wcheck in PROD. */
 719
 720     for (jc = 1; jc <= 5; ++jc) {
 721         nw = *npt;
 722         if (jc == 2 || jc == 3) {
 723             nw = *ndim;
 724         }
 725         i__2 = *npt;
 726         for (k = 1; k <= i__2; ++k) {
 727 /* L130: */
 728             prod[k + jc * prod_dim1] = zero;
 729         }
 730         i__2 = nptm;
 731         for (j = 1; j <= i__2; ++j) {
 732             sum = zero;
 733             i__1 = *npt;
 734             for (k = 1; k <= i__1; ++k) {
 735 /* L140: */
 736                 sum += zmat[k + j * zmat_dim1] * wvec[k + jc * wvec_dim1];
 737             }
 738             if (j < *idz) {
 739                 sum = -sum;
 740             }
 741             i__1 = *npt;
 742             for (k = 1; k <= i__1; ++k) {
 743 /* L150: */
 744                 prod[k + jc * prod_dim1] += sum * zmat[k + j * zmat_dim1];
 745             }
 746         }
 747         if (nw == *ndim) {
 748             i__1 = *npt;
 749             for (k = 1; k <= i__1; ++k) {
 750                 sum = zero;
 751                 i__2 = *n;
 752                 for (j = 1; j <= i__2; ++j) {
 753 /* L160: */
 754                     sum += bmat[k + j * bmat_dim1] * wvec[*npt + j + jc *
 755                             wvec_dim1];
 756                 }
 757 /* L170: */
 758                 prod[k + jc * prod_dim1] += sum;
 759             }
 760         }
 761         i__1 = *n;
 762         for (j = 1; j <= i__1; ++j) {
 763             sum = zero;
 764             i__2 = nw;
 765             for (i__ = 1; i__ <= i__2; ++i__) {
 766 /* L180: */
 767                 sum += bmat[i__ + j * bmat_dim1] * wvec[i__ + jc * wvec_dim1];
 768             }
 769 /* L190: */
 770             prod[*npt + j + jc * prod_dim1] = sum;
 771         }
 772     }
 773
 774 /* Include in DEN the part of BETA that depends on THETA. */
 775
 776     i__1 = *ndim;
 777     for (k = 1; k <= i__1; ++k) {
 778         sum = zero;
 779         for (i__ = 1; i__ <= 5; ++i__) {
 780             par[i__ - 1] = half * prod[k + i__ * prod_dim1] * wvec[k + i__ *
 781                     wvec_dim1];
 782 /* L200: */
 783             sum += par[i__ - 1];
 784         }
 785         den[0] = den[0] - par[0] - sum;
 786         tempa = prod[k + prod_dim1] * wvec[k + (wvec_dim1 << 1)] + prod[k + (
 787                 prod_dim1 << 1)] * wvec[k + wvec_dim1];
 788         tempb = prod[k + (prod_dim1 << 1)] * wvec[k + (wvec_dim1 << 2)] +
 789                 prod[k + (prod_dim1 << 2)] * wvec[k + (wvec_dim1 << 1)];
 790         tempc = prod[k + prod_dim1 * 3] * wvec[k + wvec_dim1 * 5] + prod[k +
 791                 prod_dim1 * 5] * wvec[k + wvec_dim1 * 3];
 792         den[1] = den[1] - tempa - half * (tempb + tempc);
 793         den[5] -= half * (tempb - tempc);
 794         tempa = prod[k + prod_dim1] * wvec[k + wvec_dim1 * 3] + prod[k +
 795                 prod_dim1 * 3] * wvec[k + wvec_dim1];
 796         tempb = prod[k + (prod_dim1 << 1)] * wvec[k + wvec_dim1 * 5] + prod[k
 797                 + prod_dim1 * 5] * wvec[k + (wvec_dim1 << 1)];
 798         tempc = prod[k + prod_dim1 * 3] * wvec[k + (wvec_dim1 << 2)] + prod[k
 799                 + (prod_dim1 << 2)] * wvec[k + wvec_dim1 * 3];
 800         den[2] = den[2] - tempa - half * (tempb - tempc);
 801         den[6] -= half * (tempb + tempc);
 802         tempa = prod[k + prod_dim1] * wvec[k + (wvec_dim1 << 2)] + prod[k + (
 803                 prod_dim1 << 2)] * wvec[k + wvec_dim1];
 804         den[3] = den[3] - tempa - par[1] + par[2];
 805         tempa = prod[k + prod_dim1] * wvec[k + wvec_dim1 * 5] + prod[k +
 806                 prod_dim1 * 5] * wvec[k + wvec_dim1];
 807         tempb = prod[k + (prod_dim1 << 1)] * wvec[k + wvec_dim1 * 3] + prod[k
 808                 + prod_dim1 * 3] * wvec[k + (wvec_dim1 << 1)];
 809         den[4] = den[4] - tempa - half * tempb;
 810         den[7] = den[7] - par[3] + par[4];
 811         tempa = prod[k + (prod_dim1 << 2)] * wvec[k + wvec_dim1 * 5] + prod[k
 812                 + prod_dim1 * 5] * wvec[k + (wvec_dim1 << 2)];
 813 /* L210: */
 814         den[8] -= half * tempa;
 815     }
 816
 817 /* Extend DEN so that it holds all the coefficients of DENOM. */
 818
 819     sum = zero;
 820     for (i__ = 1; i__ <= 5; ++i__) {
 821 /* Computing 2nd power */
 822         d__1 = prod[*knew + i__ * prod_dim1];
 823         par[i__ - 1] = half * (d__1 * d__1);
 824 /* L220: */
 825         sum += par[i__ - 1];
 826     }
 827     denex[0] = alpha * den[0] + par[0] + sum;
 828     tempa = two * prod[*knew + prod_dim1] * prod[*knew + (prod_dim1 << 1)];
 829     tempb = prod[*knew + (prod_dim1 << 1)] * prod[*knew + (prod_dim1 << 2)];
 830     tempc = prod[*knew + prod_dim1 * 3] * prod[*knew + prod_dim1 * 5];
 831     denex[1] = alpha * den[1] + tempa + tempb + tempc;
 832     denex[5] = alpha * den[5] + tempb - tempc;
 833     tempa = two * prod[*knew + prod_dim1] * prod[*knew + prod_dim1 * 3];
 834     tempb = prod[*knew + (prod_dim1 << 1)] * prod[*knew + prod_dim1 * 5];
 835     tempc = prod[*knew + prod_dim1 * 3] * prod[*knew + (prod_dim1 << 2)];
 836     denex[2] = alpha * den[2] + tempa + tempb - tempc;
 837     denex[6] = alpha * den[6] + tempb + tempc;
 838     tempa = two * prod[*knew + prod_dim1] * prod[*knew + (prod_dim1 << 2)];
 839     denex[3] = alpha * den[3] + tempa + par[1] - par[2];
 840     tempa = two * prod[*knew + prod_dim1] * prod[*knew + prod_dim1 * 5];
 841     denex[4] = alpha * den[4] + tempa + prod[*knew + (prod_dim1 << 1)] * prod[
 842             *knew + prod_dim1 * 3];
 843     denex[7] = alpha * den[7] + par[3] - par[4];
 844     denex[8] = alpha * den[8] + prod[*knew + (prod_dim1 << 2)] * prod[*knew +
 845             prod_dim1 * 5];
 846
 847 /* Seek the value of the angle that maximizes the modulus of DENOM. */
 848
 849     sum = denex[0] + denex[1] + denex[3] + denex[5] + denex[7];
 850     denold = sum;
 851     denmax = sum;
 852     isave = 0;
 853     iu = 49;
 854     temp = twopi / (double) (iu + 1);
 855     par[0] = one;
 856     i__1 = iu;
 857     for (i__ = 1; i__ <= i__1; ++i__) {
 858         angle = (double) i__ * temp;
 859         par[1] = cos(angle);
 860         par[2] = sin(angle);
 861         for (j = 4; j <= 8; j += 2) {
 862             par[j - 1] = par[1] * par[j - 3] - par[2] * par[j - 2];
 863 /* L230: */
 864             par[j] = par[1] * par[j - 2] + par[2] * par[j - 3];
 865         }
 866         sumold = sum;
 867         sum = zero;
 868         for (j = 1; j <= 9; ++j) {
 869 /* L240: */
 870             sum += denex[j - 1] * par[j - 1];
 871         }
 872         if (fabs(sum) > fabs(denmax)) {
 873             denmax = sum;
 874             isave = i__;
 875             tempa = sumold;
 876         } else if (i__ == isave + 1) {
 877             tempb = sum;
 878         }
 879 /* L250: */
 880     }
 881     if (isave == 0) {
 882         tempa = sum;
 883     }
 884     if (isave == iu) {
 885         tempb = denold;
 886     }
 887     step = zero;
 888     if (tempa != tempb) {
 889         tempa -= denmax;
 890         tempb -= denmax;
 891         step = half * (tempa - tempb) / (tempa + tempb);
 892     }
 893     angle = temp * ((double) isave + step);
 894
 895 /* Calculate the new parameters of the denominator, the new VLAG vector */
 896 /* and the new D. Then test for convergence. */
 897
 898     par[1] = cos(angle);
 899     par[2] = sin(angle);
 900     for (j = 4; j <= 8; j += 2) {
 901         par[j - 1] = par[1] * par[j - 3] - par[2] * par[j - 2];
 902 /* L260: */
 903         par[j] = par[1] * par[j - 2] + par[2] * par[j - 3];
 904     }
 905     *beta = zero;
 906     denmax = zero;
 907     for (j = 1; j <= 9; ++j) {
 908         *beta += den[j - 1] * par[j - 1];
 909 /* L270: */
 910         denmax += denex[j - 1] * par[j - 1];
 911     }
 912     i__1 = *ndim;
 913     for (k = 1; k <= i__1; ++k) {
 914         vlag[k] = zero;
 915         for (j = 1; j <= 5; ++j) {
 916 /* L280: */
 917             vlag[k] += prod[k + j * prod_dim1] * par[j - 1];
 918         }
 919     }
 920     tau = vlag[*knew];
 921     dd = zero;
 922     tempa = zero;
 923     tempb = zero;
 924     i__1 = *n;
 925     for (i__ = 1; i__ <= i__1; ++i__) {
 926         d__[i__] = par[1] * d__[i__] + par[2] * s[i__];
 927         w[i__] = xopt[i__] + d__[i__];
 928 /* Computing 2nd power */
 929         d__1 = d__[i__];
 930         dd += d__1 * d__1;
 931         tempa += d__[i__] * w[i__];
 932 /* L290: */
 933         tempb += w[i__] * w[i__];
 934     }
 935     if (iterc >= *n) {
 936         goto L340;
 937     }
 938     if (iterc > 1) {
 939         densav = max(densav,denold);
 940     }
 941     if (fabs(denmax) <= fabs(densav) * 1.1) {
 942         goto L340;
 943     }
 944     densav = denmax;
 945
 946 /* Set S to half the gradient of the denominator with respect to D. */
 947 /* Then branch for the next iteration. */
 948
 949     i__1 = *n;
 950     for (i__ = 1; i__ <= i__1; ++i__) {
 951         temp = tempa * xopt[i__] + tempb * d__[i__] - vlag[*npt + i__];
 952 /* L300: */
 953         s[i__] = tau * bmat[*knew + i__ * bmat_dim1] + alpha * temp;
 954     }
 955     i__1 = *npt;
 956     for (k = 1; k <= i__1; ++k) {
 957         sum = zero;
 958         i__2 = *n;
 959         for (j = 1; j <= i__2; ++j) {
 960 /* L310: */
 961             sum += xpt[k + j * xpt_dim1] * w[j];
 962         }
 963         temp = (tau * w[*n + k] - alpha * vlag[k]) * sum;
 964         i__2 = *n;
 965         for (i__ = 1; i__ <= i__2; ++i__) {
 966 /* L320: */
 967             s[i__] += temp * xpt[k + i__ * xpt_dim1];
 968         }
 969     }
 970     ss = zero;
 971     ds = zero;
 972     i__2 = *n;
 973     for (i__ = 1; i__ <= i__2; ++i__) {
 974 /* Computing 2nd power */
 975         d__1 = s[i__];
 976         ss += d__1 * d__1;
 977 /* L330: */
 978         ds += d__[i__] * s[i__];
 979     }
 980     ssden = dd * ss - ds * ds;
 981     if (ssden >= dd * 1e-8 * ss) {
 982         goto L70;
 983     }
 984
 985 /* Set the vector W before the RETURN from the subroutine. */
 986
 987 L340:
 988     /* SGJ, 2008: crude hack: truncate d to [lb,ub] bounds if given */
 989     if (lb && ub) {
 990          for (k = 1; k <= *n; ++k) {
 991               if (d__[k] > ub[k-1] - xbase[k-1] - xopt[k])
 992                    d__[k] = ub[k-1] - xbase[k-1] - xopt[k];
 993               else if (d__[k] < lb[k-1] - xbase[k-1] - xopt[k])
 994                    d__[k] = lb[k-1] - xbase[k-1] - xopt[k];
 995          }
 996     }
 997
 998     i__2 = *ndim;
 999     for (k = 1; k <= i__2; ++k) {
1000         w[k] = zero;
1001         for (j = 1; j <= 5; ++j) {
1002 /* L350: */
1003             w[k] += wvec[k + j * wvec_dim1] * par[j - 1];
1004         }
1005     }
1006     vlag[*kopt] += one;
1007     return;
1008 } /* bigden_ */
1009
1010 /*************************************************************************/
1011 /* biglag.f */
1012
1013 typedef struct {
1014      int npt, ndim, iter;
1015      double *hcol, *xpt, *bmat, *xopt;
1016      int flipsign;
1017 } lag_data;
1018
1019 /* the Lagrange function, whose absolute value biglag maximizes */
1020 static double lag(int n, const double *dx, double *grad, void *data)
1021 {
1022      lag_data *d = (lag_data *) data;
1023      int i, j, npt = d->npt, ndim = d->ndim;
1024      const double *hcol = d->hcol, *xpt = d->xpt, *bmat = d->bmat, *xopt = d->xopt;
1025      double val = 0;
1026      for (j = 0; j < n; ++j) {
1027           val += bmat[j * ndim] * (xopt[j] + dx[j]);
1028           if (grad) grad[j] = bmat[j * ndim];
1029      }
1030      for (i = 0; i < npt; ++i) {
1031           double dot = 0;
1032           for (j = 0; j < n; ++j)
1033                dot += xpt[i + j * npt] * (xopt[j] + dx[j]);
1034           val += 0.5 * hcol[i] * (dot * dot);
1035           dot *= hcol[i];
1036           if (grad) for (j = 0; j < n; ++j)
1037                          grad[j] += dot * xpt[i + j * npt];
1038      }
1039      if (d->flipsign) {
1040           val = -val;
1041           if (grad) for (j = 0; j < n; ++j) grad[j] = -grad[j];
1042      }
1043      d->iter++;
1044      return val;
1045 }
1046
1047 static nlopt_result biglag_(int *n, int *npt, double *xopt,
1048                     double *xpt, double *bmat, double *zmat, int *idz,
1049                     int *ndim, int *knew, double *delta, double *d__,
1050                     double *alpha, double *hcol, double *gc, double *gd,
1051                     double *s, double *w,
1052                     const double *xbase, const double *lb, const double *ub)
1053 {
1054     /* System generated locals */
1055     int xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1,
1056             zmat_offset, i__1, i__2;
1057     double d__1;
1058
1059     /* Local variables */
1060     int i__, j, k;
1061     double dd, gg;
1062     int iu;
1063     double sp, ss, cf1, cf2, cf3, cf4, cf5, dhd, cth, one, tau, sth, sum,
1064             half, temp, step;
1065     int nptm;
1066     double zero, angle, scale, denom;
1067     int iterc, isave;
1068     double delsq, tempa, tempb, twopi, taubeg, tauold, taumax;
1069
1070
1071 /* N is the number of variables. */
1072 /* NPT is the number of interpolation equations. */
1073 /* XOPT is the best interpolation point so far. */
1074 /* XPT contains the coordinates of the current interpolation points. */
1075 /* BMAT provides the last N columns of H. */
1076 /* ZMAT and IDZ give a factorization of the first NPT by NPT submatrix of H. */
1077 /* NDIM is the first dimension of BMAT and has the value NPT+N. */
1078 /* KNEW is the index of the interpolation point that is going to be moved. */
1079 /* DELTA is the current trust region bound. */
1080 /* D will be set to the step from XOPT to the new point. */
1081 /* ALPHA will be set to the KNEW-th diagonal element of the H matrix. */
1082 /* HCOL, GC, GD, S and W will be used for working space. */
1083
1084 /* The step D is calculated in a way that attempts to maximize the modulus */
1085 /* of LFUNC(XOPT+D), subject to the bound ||D|| .LE. DELTA, where LFUNC is */
1086 /* the KNEW-th Lagrange function. */
1087
1088 /* Set some constants. */
1089
1090     /* Parameter adjustments */
1091     zmat_dim1 = *npt;
1092     zmat_offset = 1 + zmat_dim1;
1093     zmat -= zmat_offset;
1094     xpt_dim1 = *npt;
1095     xpt_offset = 1 + xpt_dim1;
1096     xpt -= xpt_offset;
1097     --xopt;
1098     bmat_dim1 = *ndim;
1099     bmat_offset = 1 + bmat_dim1;
1100     bmat -= bmat_offset;
1101     --d__;
1102     --hcol;
1103     --gc;
1104     --gd;
1105     --s;
1106     --w;
1107
1108     /* Function Body */
1109     half = .5;
1110     one = 1.;
1111     zero = 0.;
1112     twopi = atan(one) * 8.;
1113     delsq = *delta * *delta;
1114     nptm = *npt - *n - 1;
1115
1116 /* Set the first NPT components of HCOL to the leading elements of the */
1117 /* KNEW-th column of H. */
1118
1119     iterc = 0;
1120     i__1 = *npt;
1121     for (k = 1; k <= i__1; ++k) {
1122 /* L10: */
1123         hcol[k] = zero;
1124     }
1125     i__1 = nptm;
1126     for (j = 1; j <= i__1; ++j) {
1127         temp = zmat[*knew + j * zmat_dim1];
1128         if (j < *idz) {
1129             temp = -temp;
1130         }
1131         i__2 = *npt;
1132         for (k = 1; k <= i__2; ++k) {
1133 /* L20: */
1134             hcol[k] += temp * zmat[k + j * zmat_dim1];
1135         }
1136     }
1137     *alpha = hcol[*knew];
1138
1139 /* Set the unscaled initial direction D. Form the gradient of LFUNC at */
1140 /* XOPT, and multiply D by the second derivative matrix of LFUNC. */
1141
1142     dd = zero;
1143     i__2 = *n;
1144     for (i__ = 1; i__ <= i__2; ++i__) {
1145         d__[i__] = xpt[*knew + i__ * xpt_dim1] - xopt[i__];
1146         gc[i__] = bmat[*knew + i__ * bmat_dim1];
1147         gd[i__] = zero;
1148 /* L30: */
1149 /* Computing 2nd power */
1150         d__1 = d__[i__];
1151         dd += d__1 * d__1;
1152     }
1153     i__2 = *npt;
1154     for (k = 1; k <= i__2; ++k) {
1155         temp = zero;
1156         sum = zero;
1157         i__1 = *n;
1158         for (j = 1; j <= i__1; ++j) {
1159             temp += xpt[k + j * xpt_dim1] * xopt[j];
1160 /* L40: */
1161             sum += xpt[k + j * xpt_dim1] * d__[j];
1162         }
1163         temp = hcol[k] * temp;
1164         sum = hcol[k] * sum;
1165         i__1 = *n;
1166         for (i__ = 1; i__ <= i__1; ++i__) {
1167             gc[i__] += temp * xpt[k + i__ * xpt_dim1];
1168 /* L50: */
1169             gd[i__] += sum * xpt[k + i__ * xpt_dim1];
1170         }
1171     }
1172
1173 /* Scale D and GD, with a sign change if required. Set S to another */
1174 /* vector in the initial two dimensional subspace. */
1175
1176     gg = zero;
1177     sp = zero;
1178     dhd = zero;
1179     i__1 = *n;
1180     for (i__ = 1; i__ <= i__1; ++i__) {
1181 /* Computing 2nd power */
1182         d__1 = gc[i__];
1183         gg += d__1 * d__1;
1184         sp += d__[i__] * gc[i__];
1185 /* L60: */
1186         dhd += d__[i__] * gd[i__];
1187     }
1188     scale = *delta / sqrt(dd);
1189     if (sp * dhd < zero) {
1190         scale = -scale;
1191     }
1192     temp = zero;
1193     if (sp * sp > dd * .99 * gg) {
1194         temp = one;
1195     }
1196     tau = scale * (fabs(sp) + half * scale * fabs(dhd));
1197     if (gg * delsq < tau * .01 * tau) {
1198         temp = one;
1199     }
1200     i__1 = *n;
1201     for (i__ = 1; i__ <= i__1; ++i__) {
1202         d__[i__] = scale * d__[i__];
1203         gd[i__] = scale * gd[i__];
1204 /* L70: */
1205         s[i__] = gc[i__] + temp * gd[i__];
1206     }
1207
1208     if (lb && ub) {
1209          double minf, *dlb, *dub, *xtol;
1210          lag_data ld;
1211          ld.npt = *npt;
1212          ld.ndim = *ndim;
1213          ld.iter = 0;
1214          ld.hcol = &hcol[1];
1215          ld.xpt = &xpt[1 + 1 * xpt_dim1];
1216          ld.bmat = &bmat[*knew + 1 * bmat_dim1];
1217          ld.xopt = &xopt[1];
1218          ld.flipsign = 0;
1219          dlb = &gc[1]; dub = &gd[1]; xtol = &s[1];
1220          /* make sure rounding errors don't push initial |d| > delta */
1221          for (j = 1; j <= *n; ++j) d__[j] *= 0.99999;
1222          for (j = 0; j < *n; ++j) {
1223               dlb[j] = -(dub[j] = *delta);
1224               if (dlb[j] < lb[j] - xbase[j] - xopt[j+1])
1225                    dlb[j] = lb[j] - xbase[j] - xopt[j+1];
1226               if (dub[j] > ub[j] - xbase[j] - xopt[j+1])
1227                    dub[j] = ub[j] - xbase[j] - xopt[j+1];
1228               if (dlb[j] > 0) dlb[j] = 0;
1229               if (dub[j] < 0) dub[j] = 0;
1230               if (d__[j+1] < dlb[j]) d__[j+1] = dlb[j];
1231               else if (d__[j+1] > dub[j]) d__[j+1] = dub[j];
1232               xtol[j] = 1e-5 * *delta;
1233          }
1234          ld.flipsign = lag(*n, &d__[1], 0, &ld) > 0; /* maximize if > 0 */
1235          return nlopt_minimize_constrained(NLOPT_LD_MMA, *n, lag, &ld,
1236                                     1, rho_constraint, delta, sizeof(double),
1237                                     dlb, dub, &d__[1], &minf, -HUGE_VAL,
1238                                     0., 0., 0., xtol, 1000, 0.);
1239     }
1240
1241 /* Begin the iteration by overwriting S with a vector that has the */
1242 /* required length and direction, except that termination occurs if */
1243 /* the given D and S are nearly parallel. */
1244
1245 L80:
1246     ++iterc;
1247     dd = zero;
1248     sp = zero;
1249     ss = zero;
1250     i__1 = *n;
1251     for (i__ = 1; i__ <= i__1; ++i__) {
1252 /* Computing 2nd power */
1253         d__1 = d__[i__];
1254         dd += d__1 * d__1;
1255         sp += d__[i__] * s[i__];
1256 /* L90: */
1257 /* Computing 2nd power */
1258         d__1 = s[i__];
1259         ss += d__1 * d__1;
1260     }
1261     temp = dd * ss - sp * sp;
1262     if (temp <= dd * 1e-8 * ss) {
1263         goto L160;
1264     }
1265     denom = sqrt(temp);
1266     i__1 = *n;
1267     for (i__ = 1; i__ <= i__1; ++i__) {
1268         s[i__] = (dd * s[i__] - sp * d__[i__]) / denom;
1269 /* L100: */
1270         w[i__] = zero;
1271     }
1272
1273 /* Calculate the coefficients of the objective function on the circle, */
1274 /* beginning with the multiplication of S by the second derivative matrix. */
1275
1276     i__1 = *npt;
1277     for (k = 1; k <= i__1; ++k) {
1278         sum = zero;
1279         i__2 = *n;
1280         for (j = 1; j <= i__2; ++j) {
1281 /* L110: */
1282             sum += xpt[k + j * xpt_dim1] * s[j];
1283         }
1284         sum = hcol[k] * sum;
1285         i__2 = *n;
1286         for (i__ = 1; i__ <= i__2; ++i__) {
1287 /* L120: */
1288             w[i__] += sum * xpt[k + i__ * xpt_dim1];
1289         }
1290     }
1291     cf1 = zero;
1292     cf2 = zero;
1293     cf3 = zero;
1294     cf4 = zero;
1295     cf5 = zero;
1296     i__2 = *n;
1297     for (i__ = 1; i__ <= i__2; ++i__) {
1298         cf1 += s[i__] * w[i__];
1299         cf2 += d__[i__] * gc[i__];
1300         cf3 += s[i__] * gc[i__];
1301         cf4 += d__[i__] * gd[i__];
1302 /* L130: */
1303         cf5 += s[i__] * gd[i__];
1304     }
1305     cf1 = half * cf1;
1306     cf4 = half * cf4 - cf1;
1307
1308 /* Seek the value of the angle that maximizes the modulus of TAU. */
1309
1310     taubeg = cf1 + cf2 + cf4;
1311     taumax = taubeg;
1312     tauold = taubeg;
1313     isave = 0;
1314     iu = 49;
1315     temp = twopi / (double) (iu + 1);
1316     i__2 = iu;
1317     for (i__ = 1; i__ <= i__2; ++i__) {
1318         angle = (double) i__ * temp;
1319         cth = cos(angle);
1320         sth = sin(angle);
1321         tau = cf1 + (cf2 + cf4 * cth) * cth + (cf3 + cf5 * cth) * sth;
1322         if (fabs(tau) > fabs(taumax)) {
1323             taumax = tau;
1324             isave = i__;
1325             tempa = tauold;
1326         } else if (i__ == isave + 1) {
1327             tempb = tau;
1328         }
1329 /* L140: */
1330         tauold = tau;
1331     }
1332     if (isave == 0) {
1333         tempa = tau;
1334     }
1335     if (isave == iu) {
1336         tempb = taubeg;
1337     }
1338     step = zero;
1339     if (tempa != tempb) {
1340         tempa -= taumax;
1341         tempb -= taumax;
1342         step = half * (tempa - tempb) / (tempa + tempb);
1343     }
1344     angle = temp * ((double) isave + step);
1345
1346 /* Calculate the new D and GD. Then test for convergence. */
1347
1348     cth = cos(angle);
1349     sth = sin(angle);
1350     tau = cf1 + (cf2 + cf4 * cth) * cth + (cf3 + cf5 * cth) * sth;
1351     i__2 = *n;
1352     for (i__ = 1; i__ <= i__2; ++i__) {
1353         d__[i__] = cth * d__[i__] + sth * s[i__];
1354         gd[i__] = cth * gd[i__] + sth * w[i__];
1355 /* L150: */
1356         s[i__] = gc[i__] + gd[i__];
1357     }
1358     if (fabs(tau) <= fabs(taubeg) * 1.1) {
1359         goto L160;
1360     }
1361     if (iterc < *n) {
1362         goto L80;
1363     }
1364 L160:
1365     return NLOPT_SUCCESS;
1366 } /* biglag_ */
1367
1368
1369
1370 /*************************************************************************/
1371 /* update.f */
1372
1373 static void update_(int *n, int *npt, double *bmat,
1374         double *zmat, int *idz, int *ndim, double *vlag,
1375         double *beta, int *knew, double *w)
1376 {
1377     /* System generated locals */
1378     int bmat_dim1, bmat_offset, zmat_dim1, zmat_offset, i__1, i__2;
1379     double d__1, d__2;
1380
1381     /* Local variables */
1382     int i__, j, ja, jb, jl, jp;
1383     double one, tau, temp;
1384     int nptm;
1385     double zero;
1386     int iflag;
1387     double scala, scalb_, alpha, denom, tempa, tempb, tausq;
1388
1389
1390 /* The arrays BMAT and ZMAT with IDZ are updated, in order to shift the */
1391 /* interpolation point that has index KNEW. On entry, VLAG contains the */
1392 /* components of the vector Theta*Wcheck+e_b of the updating formula */
1393 /* (6.11), and BETA holds the value of the parameter that has this name. */
1394 /* The vector W is used for working space. */
1395
1396 /* Set some constants. */
1397
1398     /* Parameter adjustments */
1399     zmat_dim1 = *npt;
1400     zmat_offset = 1 + zmat_dim1;
1401     zmat -= zmat_offset;
1402     bmat_dim1 = *ndim;
1403     bmat_offset = 1 + bmat_dim1;
1404     bmat -= bmat_offset;
1405     --vlag;
1406     --w;
1407
1408     /* Function Body */
1409     one = 1.;
1410     zero = 0.;
1411     nptm = *npt - *n - 1;
1412
1413 /* Apply the rotations that put zeros in the KNEW-th row of ZMAT. */
1414
1415     jl = 1;
1416     i__1 = nptm;
1417     for (j = 2; j <= i__1; ++j) {
1418         if (j == *idz) {
1419             jl = *idz;
1420         } else if (zmat[*knew + j * zmat_dim1] != zero) {
1421 /* Computing 2nd power */
1422             d__1 = zmat[*knew + jl * zmat_dim1];
1423 /* Computing 2nd power */
1424             d__2 = zmat[*knew + j * zmat_dim1];
1425             temp = sqrt(d__1 * d__1 + d__2 * d__2);
1426             tempa = zmat[*knew + jl * zmat_dim1] / temp;
1427             tempb = zmat[*knew + j * zmat_dim1] / temp;
1428             i__2 = *npt;
1429             for (i__ = 1; i__ <= i__2; ++i__) {
1430                 temp = tempa * zmat[i__ + jl * zmat_dim1] + tempb * zmat[i__
1431                         + j * zmat_dim1];
1432                 zmat[i__ + j * zmat_dim1] = tempa * zmat[i__ + j * zmat_dim1]
1433                         - tempb * zmat[i__ + jl * zmat_dim1];
1434 /* L10: */
1435                 zmat[i__ + jl * zmat_dim1] = temp;
1436             }
1437             zmat[*knew + j * zmat_dim1] = zero;
1438         }
1439 /* L20: */
1440     }
1441
1442 /* Put the first NPT components of the KNEW-th column of HLAG into W, */
1443 /* and calculate the parameters of the updating formula. */
1444
1445     tempa = zmat[*knew + zmat_dim1];
1446     if (*idz >= 2) {
1447         tempa = -tempa;
1448     }
1449     if (jl > 1) {
1450         tempb = zmat[*knew + jl * zmat_dim1];
1451     }
1452     i__1 = *npt;
1453     for (i__ = 1; i__ <= i__1; ++i__) {
1454         w[i__] = tempa * zmat[i__ + zmat_dim1];
1455         if (jl > 1) {
1456             w[i__] += tempb * zmat[i__ + jl * zmat_dim1];
1457         }
1458 /* L30: */
1459     }
1460     alpha = w[*knew];
1461     tau = vlag[*knew];
1462     tausq = tau * tau;
1463     denom = alpha * *beta + tausq;
1464     vlag[*knew] -= one;
1465
1466 /* Complete the updating of ZMAT when there is only one nonzero element */
1467 /* in the KNEW-th row of the new matrix ZMAT, but, if IFLAG is set to one, */
1468 /* then the first column of ZMAT will be exchanged with another one later. */
1469
1470     iflag = 0;
1471     if (jl == 1) {
1472         temp = sqrt((fabs(denom)));
1473         tempb = tempa / temp;
1474         tempa = tau / temp;
1475         i__1 = *npt;
1476         for (i__ = 1; i__ <= i__1; ++i__) {
1477 /* L40: */
1478             zmat[i__ + zmat_dim1] = tempa * zmat[i__ + zmat_dim1] - tempb *
1479                     vlag[i__];
1480         }
1481         if (*idz == 1 && temp < zero) {
1482             *idz = 2;
1483         }
1484         if (*idz >= 2 && temp >= zero) {
1485             iflag = 1;
1486         }
1487     } else {
1488
1489 /* Complete the updating of ZMAT in the alternative case. */
1490
1491         ja = 1;
1492         if (*beta >= zero) {
1493             ja = jl;
1494         }
1495         jb = jl + 1 - ja;
1496         temp = zmat[*knew + jb * zmat_dim1] / denom;
1497         tempa = temp * *beta;
1498         tempb = temp * tau;
1499         temp = zmat[*knew + ja * zmat_dim1];
1500         scala = one / sqrt(fabs(*beta) * temp * temp + tausq);
1501         scalb_ = scala * sqrt((fabs(denom)));
1502         i__1 = *npt;
1503         for (i__ = 1; i__ <= i__1; ++i__) {
1504             zmat[i__ + ja * zmat_dim1] = scala * (tau * zmat[i__ + ja *
1505                     zmat_dim1] - temp * vlag[i__]);
1506 /* L50: */
1507             zmat[i__ + jb * zmat_dim1] = scalb_ * (zmat[i__ + jb * zmat_dim1]
1508                     - tempa * w[i__] - tempb * vlag[i__]);
1509         }
1510         if (denom <= zero) {
1511             if (*beta < zero) {
1512                 ++(*idz);
1513             }
1514             if (*beta >= zero) {
1515                 iflag = 1;
1516             }
1517         }
1518     }
1519
1520 /* IDZ is reduced in the following case, and usually the first column */
1521 /* of ZMAT is exchanged with a later one. */
1522
1523     if (iflag == 1) {
1524         --(*idz);
1525         i__1 = *npt;
1526         for (i__ = 1; i__ <= i__1; ++i__) {
1527             temp = zmat[i__ + zmat_dim1];
1528             zmat[i__ + zmat_dim1] = zmat[i__ + *idz * zmat_dim1];
1529 /* L60: */
1530             zmat[i__ + *idz * zmat_dim1] = temp;
1531         }
1532     }
1533
1534 /* Finally, update the matrix BMAT. */
1535
1536     i__1 = *n;
1537     for (j = 1; j <= i__1; ++j) {
1538         jp = *npt + j;
1539         w[jp] = bmat[*knew + j * bmat_dim1];
1540         tempa = (alpha * vlag[jp] - tau * w[jp]) / denom;
1541         tempb = (-(*beta) * w[jp] - tau * vlag[jp]) / denom;
1542         i__2 = jp;
1543         for (i__ = 1; i__ <= i__2; ++i__) {
1544             bmat[i__ + j * bmat_dim1] = bmat[i__ + j * bmat_dim1] + tempa *
1545                     vlag[i__] + tempb * w[i__];
1546             if (i__ > *npt) {
1547                 bmat[jp + (i__ - *npt) * bmat_dim1] = bmat[i__ + j *
1548                         bmat_dim1];
1549             }
1550 /* L70: */
1551         }
1552     }
1553     return;
1554 } /* update_ */
1555
1556
1557 /*************************************************************************/
1558 /* newuob.f */
1559
1560 static nlopt_result newuob_(int *n, int *npt, double *x,
1561                             double *rhobeg,
1562                             const double *lb, const double *ub,
1563                             nlopt_stopping *stop, double *minf,
1564                             newuoa_func calfun, void *calfun_data,
1565                     double *xbase, double *xopt, double *xnew,
1566                     double *xpt, double *fval, double *gq, double *hq,
1567                     double *pq, double *bmat, double *zmat, int *ndim,
1568                     double *d__, double *vlag, double *w)
1569 {
1570     /* System generated locals */
1571     int xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1,
1572             zmat_offset, i__1, i__2, i__3;
1573     double d__1, d__2, d__3;
1574
1575     /* Local variables */
1576     double f;
1577     int i__, j, k, ih, nf, nh, ip, jp;
1578     double dx;
1579     int np, nfm;
1580     double one;
1581     int idz;
1582     double dsq, rho;
1583     int ipt, jpt;
1584     double sum, fbeg, diff, half, beta;
1585     int nfmm;
1586     double gisq;
1587     int knew;
1588     double temp, suma, sumb, fopt = HUGE_VAL, bsum, gqsq;
1589     int kopt, nptm;
1590     double zero, xipt, xjpt, sumz, diffa, diffb, diffc, hdiag, alpha,
1591             delta, recip, reciq, fsave;
1592     int ksave, nfsav, itemp;
1593     double dnorm, ratio, dstep, tenth, vquad;
1594     int ktemp;
1595     double tempq;
1596     int itest;
1597     double rhosq;
1598     double detrat, crvmin;
1599     double distsq;
1600     double xoptsq;
1601     double rhoend;
1602     nlopt_result rc = NLOPT_SUCCESS, rc2;
1603
1604 /* SGJ, 2008: compute rhoend from NLopt stop info */
1605     rhoend = stop->xtol_rel * (*rhobeg);
1606     for (j = 0; j < *n; ++j)
1607          if (rhoend < stop->xtol_abs[j])
1608               rhoend = stop->xtol_abs[j];
1609
1610 /* The arguments N, NPT, X, RHOBEG, RHOEND, IPRINT and MAXFUN are identical */
1611 /*   to the corresponding arguments in SUBROUTINE NEWUOA. */
1612 /* XBASE will hold a shift of origin that should reduce the contributions */
1613 /*   from rounding errors to values of the model and Lagrange functions. */
1614 /* XOPT will be set to the displacement from XBASE of the vector of */
1615 /*   variables that provides the least calculated F so far. */
1616 /* XNEW will be set to the displacement from XBASE of the vector of */
1617 /*   variables for the current calculation of F. */
1618 /* XPT will contain the interpolation point coordinates relative to XBASE. */
1619 /* FVAL will hold the values of F at the interpolation points. */
1620 /* GQ will hold the gradient of the quadratic model at XBASE. */
1621 /* HQ will hold the explicit second derivatives of the quadratic model. */
1622 /* PQ will contain the parameters of the implicit second derivatives of */
1623 /*   the quadratic model. */
1624 /* BMAT will hold the last N columns of H. */
1625 /* ZMAT will hold the factorization of the leading NPT by NPT submatrix of */
1626 /*   H, this factorization being ZMAT times Diag(DZ) times ZMAT^T, where */
1627 /*   the elements of DZ are plus or minus one, as specified by IDZ. */
1628 /* NDIM is the first dimension of BMAT and has the value NPT+N. */
1629 /* D is reserved for trial steps from XOPT. */
1630 /* VLAG will contain the values of the Lagrange functions at a new point X. */
1631 /*   They are part of a product that requires VLAG to be of length NDIM. */
1632 /* The array W will be used for working space. Its length must be at least */
1633 /*   10*NDIM = 10*(NPT+N). */
1634
1635 /* Set some constants. */
1636
1637     /* Parameter adjustments */
1638     zmat_dim1 = *npt;
1639     zmat_offset = 1 + zmat_dim1;
1640     zmat -= zmat_offset;
1641     xpt_dim1 = *npt;
1642     xpt_offset = 1 + xpt_dim1;
1643     xpt -= xpt_offset;
1644     --x;
1645     --xbase;
1646     --xopt;
1647     --xnew;
1648     --fval;
1649     --gq;
1650     --hq;
1651     --pq;
1652     bmat_dim1 = *ndim;
1653     bmat_offset = 1 + bmat_dim1;
1654     bmat -= bmat_offset;
1655     --d__;
1656     --vlag;
1657     --w;
1658
1659     /* Function Body */
1660     half = .5;
1661     one = 1.;
1662     tenth = .1;
1663     zero = 0.;
1664     np = *n + 1;
1665     nh = *n * np / 2;
1666     nptm = *npt - np;
1667
1668 /* Set the initial elements of XPT, BMAT, HQ, PQ and ZMAT to zero. */
1669
1670     i__1 = *n;
1671     for (j = 1; j <= i__1; ++j) {
1672         xbase[j] = x[j];
1673         i__2 = *npt;
1674         for (k = 1; k <= i__2; ++k) {
1675 /* L10: */
1676             xpt[k + j * xpt_dim1] = zero;
1677         }
1678         i__2 = *ndim;
1679         for (i__ = 1; i__ <= i__2; ++i__) {
1680 /* L20: */
1681             bmat[i__ + j * bmat_dim1] = zero;
1682         }
1683     }
1684     i__2 = nh;
1685     for (ih = 1; ih <= i__2; ++ih) {
1686 /* L30: */
1687         hq[ih] = zero;
1688     }
1689     i__2 = *npt;
1690     for (k = 1; k <= i__2; ++k) {
1691         pq[k] = zero;
1692         i__1 = nptm;
1693         for (j = 1; j <= i__1; ++j) {
1694 /* L40: */
1695             zmat[k + j * zmat_dim1] = zero;
1696         }
1697     }
1698
1699 /* Begin the initialization procedure. NF becomes one more than the number */
1700 /* of function values so far. The coordinates of the displacement of the */
1701 /* next initial interpolation point from XBASE are set in XPT(NF,.). */
1702
1703     rhosq = *rhobeg * *rhobeg;
1704     recip = one / rhosq;
1705     reciq = sqrt(half) / rhosq;
1706     nf = 0;
1707 L50:
1708     nfm = nf;
1709     nfmm = nf - *n;
1710     ++nf;
1711     if (nfm <= *n << 1) {
1712         if (nfm >= 1 && nfm <= *n) {
1713             xpt[nf + nfm * xpt_dim1] = *rhobeg;
1714         } else if (nfm > *n) {
1715             xpt[nf + nfmm * xpt_dim1] = -(*rhobeg);
1716         }
1717     } else {
1718         itemp = (nfmm - 1) / *n;
1719         jpt = nfm - itemp * *n - *n;
1720         ipt = jpt + itemp;
1721         if (ipt > *n) {
1722             itemp = jpt;
1723             jpt = ipt - *n;
1724             ipt = itemp;
1725         }
1726         xipt = *rhobeg;
1727         if (fval[ipt + np] < fval[ipt + 1]) {
1728             xipt = -xipt;
1729         }
1730         xjpt = *rhobeg;
1731         if (fval[jpt + np] < fval[jpt + 1]) {
1732             xjpt = -xjpt;
1733         }
1734         xpt[nf + ipt * xpt_dim1] = xipt;
1735         xpt[nf + jpt * xpt_dim1] = xjpt;
1736     }
1737
1738 /* Calculate the next value of F, label 70 being reached immediately */
1739 /* after this calculation. The least function value so far and its index */
1740 /* are required. */
1741
1742     i__1 = *n;
1743     for (j = 1; j <= i__1; ++j) {
1744 /* L60: */
1745         x[j] = xpt[nf + j * xpt_dim1] + xbase[j];
1746     }
1747     goto L310;
1748 L70:
1749     fval[nf] = f;
1750     if (nf == 1) {
1751         fbeg = f;
1752         fopt = f;
1753         kopt = 1;
1754     } else if (f < fopt) {
1755         fopt = f;
1756         kopt = nf;
1757     }
1758
1759 /* Set the nonzero initial elements of BMAT and the quadratic model in */
1760 /* the cases when NF is at most 2*N+1. */
1761
1762     if (nfm <= *n << 1) {
1763         if (nfm >= 1 && nfm <= *n) {
1764             gq[nfm] = (f - fbeg) / *rhobeg;
1765             if (*npt < nf + *n) {
1766                 bmat[nfm * bmat_dim1 + 1] = -one / *rhobeg;
1767                 bmat[nf + nfm * bmat_dim1] = one / *rhobeg;
1768                 bmat[*npt + nfm + nfm * bmat_dim1] = -half * rhosq;
1769             }
1770         } else if (nfm > *n) {
1771             bmat[nf - *n + nfmm * bmat_dim1] = half / *rhobeg;
1772             bmat[nf + nfmm * bmat_dim1] = -half / *rhobeg;
1773             zmat[nfmm * zmat_dim1 + 1] = -reciq - reciq;
1774             zmat[nf - *n + nfmm * zmat_dim1] = reciq;
1775             zmat[nf + nfmm * zmat_dim1] = reciq;
1776             ih = nfmm * (nfmm + 1) / 2;
1777             temp = (fbeg - f) / *rhobeg;
1778             hq[ih] = (gq[nfmm] - temp) / *rhobeg;
1779             gq[nfmm] = half * (gq[nfmm] + temp);
1780         }
1781
1782 /* Set the off-diagonal second derivatives of the Lagrange functions and */
1783 /* the initial quadratic model. */
1784
1785     } else {
1786         ih = ipt * (ipt - 1) / 2 + jpt;
1787         if (xipt < zero) {
1788             ipt += *n;
1789         }
1790         if (xjpt < zero) {
1791             jpt += *n;
1792         }
1793         zmat[nfmm * zmat_dim1 + 1] = recip;
1794         zmat[nf + nfmm * zmat_dim1] = recip;
1795         zmat[ipt + 1 + nfmm * zmat_dim1] = -recip;
1796         zmat[jpt + 1 + nfmm * zmat_dim1] = -recip;
1797         hq[ih] = (fbeg - fval[ipt + 1] - fval[jpt + 1] + f) / (xipt * xjpt);
1798     }
1799     if (nf < *npt) {
1800         goto L50;
1801     }
1802
1803 /* Begin the iterative procedure, because the initial model is complete. */
1804
1805     rho = *rhobeg;
1806     delta = rho;
1807     idz = 1;
1808     diffa = zero;
1809     diffb = zero;
1810     itest = 0;
1811     xoptsq = zero;
1812     i__1 = *n;
1813     for (i__ = 1; i__ <= i__1; ++i__) {
1814         xopt[i__] = xpt[kopt + i__ * xpt_dim1];
1815 /* L80: */
1816 /* Computing 2nd power */
1817         d__1 = xopt[i__];
1818         xoptsq += d__1 * d__1;
1819     }
1820 L90:
1821     nfsav = nf;
1822
1823 /* Generate the next trust region step and test its length. Set KNEW */
1824 /* to -1 if the purpose of the next F will be to improve the model. */
1825
1826 L100:
1827     knew = 0;
1828     rc2 = trsapp_(n, npt, &xopt[1], &xpt[xpt_offset], &gq[1], &hq[1], &pq[1], &
1829             delta, &d__[1], &w[1], &w[np], &w[np + *n], &w[np + (*n << 1)], &
1830             crvmin, &xbase[1], lb, ub);
1831     if (rc2 < 0) { rc = rc2; goto L530; }
1832     dsq = zero;
1833     i__1 = *n;
1834     for (i__ = 1; i__ <= i__1; ++i__) {
1835 /* L110: */
1836 /* Computing 2nd power */
1837         d__1 = d__[i__];
1838         dsq += d__1 * d__1;
1839     }
1840 /* Computing MIN */
1841     d__1 = delta, d__2 = sqrt(dsq);
1842     dnorm = min(d__1,d__2);
1843     if (dnorm < half * rho) {
1844         knew = -1;
1845         delta = tenth * delta;
1846         ratio = -1.;
1847         if (delta <= rho * 1.5) {
1848             delta = rho;
1849         }
1850         if (nf <= nfsav + 2) {
1851             goto L460;
1852         }
1853         temp = crvmin * .125 * rho * rho;
1854 /* Computing MAX */
1855         d__1 = max(diffa,diffb);
1856         if (temp <= max(d__1,diffc)) {
1857             goto L460;
1858         }
1859         goto L490;
1860     }
1861
1862 /* Shift XBASE if XOPT may be too far from XBASE. First make the changes */
1863 /* to BMAT that do not depend on ZMAT. */
1864
1865 L120:
1866     if (dsq <= xoptsq * .001) {
1867         tempq = xoptsq * .25;
1868         i__1 = *npt;
1869         for (k = 1; k <= i__1; ++k) {
1870             sum = zero;
1871             i__2 = *n;
1872             for (i__ = 1; i__ <= i__2; ++i__) {
1873 /* L130: */
1874                 sum += xpt[k + i__ * xpt_dim1] * xopt[i__];
1875             }
1876             temp = pq[k] * sum;
1877             sum -= half * xoptsq;
1878             w[*npt + k] = sum;
1879             i__2 = *n;
1880             for (i__ = 1; i__ <= i__2; ++i__) {
1881                 gq[i__] += temp * xpt[k + i__ * xpt_dim1];
1882                 xpt[k + i__ * xpt_dim1] -= half * xopt[i__];
1883                 vlag[i__] = bmat[k + i__ * bmat_dim1];
1884                 w[i__] = sum * xpt[k + i__ * xpt_dim1] + tempq * xopt[i__];
1885                 ip = *npt + i__;
1886                 i__3 = i__;
1887                 for (j = 1; j <= i__3; ++j) {
1888 /* L140: */
1889                     bmat[ip + j * bmat_dim1] = bmat[ip + j * bmat_dim1] +
1890                             vlag[i__] * w[j] + w[i__] * vlag[j];
1891                 }
1892             }
1893         }
1894
1895 /* Then the revisions of BMAT that depend on ZMAT are calculated. */
1896
1897         i__3 = nptm;
1898         for (k = 1; k <= i__3; ++k) {
1899             sumz = zero;
1900             i__2 = *npt;
1901             for (i__ = 1; i__ <= i__2; ++i__) {
1902                 sumz += zmat[i__ + k * zmat_dim1];
1903 /* L150: */
1904                 w[i__] = w[*npt + i__] * zmat[i__ + k * zmat_dim1];
1905             }
1906             i__2 = *n;
1907             for (j = 1; j <= i__2; ++j) {
1908                 sum = tempq * sumz * xopt[j];
1909                 i__1 = *npt;
1910                 for (i__ = 1; i__ <= i__1; ++i__) {
1911 /* L160: */
1912                     sum += w[i__] * xpt[i__ + j * xpt_dim1];
1913                 }
1914                 vlag[j] = sum;
1915                 if (k < idz) {
1916                     sum = -sum;
1917                 }
1918                 i__1 = *npt;
1919                 for (i__ = 1; i__ <= i__1; ++i__) {
1920 /* L170: */
1921                     bmat[i__ + j * bmat_dim1] += sum * zmat[i__ + k *
1922                             zmat_dim1];
1923                 }
1924             }
1925             i__1 = *n;
1926             for (i__ = 1; i__ <= i__1; ++i__) {
1927                 ip = i__ + *npt;
1928                 temp = vlag[i__];
1929                 if (k < idz) {
1930                     temp = -temp;
1931                 }
1932                 i__2 = i__;
1933                 for (j = 1; j <= i__2; ++j) {
1934 /* L180: */
1935                     bmat[ip + j * bmat_dim1] += temp * vlag[j];
1936                 }
1937             }
1938         }
1939
1940 /* The following instructions complete the shift of XBASE, including */
1941 /* the changes to the parameters of the quadratic model. */
1942
1943         ih = 0;
1944         i__2 = *n;
1945         for (j = 1; j <= i__2; ++j) {
1946             w[j] = zero;
1947             i__1 = *npt;
1948             for (k = 1; k <= i__1; ++k) {
1949                 w[j] += pq[k] * xpt[k + j * xpt_dim1];
1950 /* L190: */
1951                 xpt[k + j * xpt_dim1] -= half * xopt[j];
1952             }
1953             i__1 = j;
1954             for (i__ = 1; i__ <= i__1; ++i__) {
1955                 ++ih;
1956                 if (i__ < j) {
1957                     gq[j] += hq[ih] * xopt[i__];
1958                 }
1959                 gq[i__] += hq[ih] * xopt[j];
1960                 hq[ih] = hq[ih] + w[i__] * xopt[j] + xopt[i__] * w[j];
1961 /* L200: */
1962                 bmat[*npt + i__ + j * bmat_dim1] = bmat[*npt + j + i__ *
1963                         bmat_dim1];
1964             }
1965         }
1966         i__1 = *n;
1967         for (j = 1; j <= i__1; ++j) {
1968             xbase[j] += xopt[j];
1969 /* L210: */
1970             xopt[j] = zero;
1971         }
1972         xoptsq = zero;
1973     }
1974
1975 /* Pick the model step if KNEW is positive. A different choice of D */
1976 /* may be made later, if the choice of D by BIGLAG causes substantial */
1977 /* cancellation in DENOM. */
1978
1979     if (knew > 0) {
1980        rc2 =
1981          biglag_(n, npt, &xopt[1], &xpt[xpt_offset], &bmat[bmat_offset], &zmat[
1982                 zmat_offset], &idz, ndim, &knew, &dstep, &d__[1], &alpha, &
1983                 vlag[1], &vlag[*npt + 1], &w[1], &w[np], &w[np + *n],
1984                 &xbase[1], lb, ub);
1985        if (rc2 < 0) { rc = rc2; goto L530; }
1986     }
1987
1988 /* Calculate VLAG and BETA for the current choice of D. The first NPT */
1989 /* components of W_check will be held in W. */
1990
1991     i__1 = *npt;
1992     for (k = 1; k <= i__1; ++k) {
1993         suma = zero;
1994         sumb = zero;
1995         sum = zero;
1996         i__2 = *n;
1997         for (j = 1; j <= i__2; ++j) {
1998             suma += xpt[k + j * xpt_dim1] * d__[j];
1999             sumb += xpt[k + j * xpt_dim1] * xopt[j];
2000 /* L220: */
2001             sum += bmat[k + j * bmat_dim1] * d__[j];
2002         }
2003         w[k] = suma * (half * suma + sumb);
2004 /* L230: */
2005         vlag[k] = sum;
2006     }
2007     beta = zero;
2008     i__1 = nptm;
2009     for (k = 1; k <= i__1; ++k) {
2010         sum = zero;
2011         i__2 = *npt;
2012         for (i__ = 1; i__ <= i__2; ++i__) {
2013 /* L240: */
2014             sum += zmat[i__ + k * zmat_dim1] * w[i__];
2015         }
2016         if (k < idz) {
2017             beta += sum * sum;
2018             sum = -sum;
2019         } else {
2020             beta -= sum * sum;
2021         }
2022         i__2 = *npt;
2023         for (i__ = 1; i__ <= i__2; ++i__) {
2024 /* L250: */
2025             vlag[i__] += sum * zmat[i__ + k * zmat_dim1];
2026         }
2027     }
2028     bsum = zero;
2029     dx = zero;
2030     i__2 = *n;
2031     for (j = 1; j <= i__2; ++j) {
2032         sum = zero;
2033         i__1 = *npt;
2034         for (i__ = 1; i__ <= i__1; ++i__) {
2035 /* L260: */
2036             sum += w[i__] * bmat[i__ + j * bmat_dim1];
2037         }
2038         bsum += sum * d__[j];
2039         jp = *npt + j;
2040         i__1 = *n;
2041         for (k = 1; k <= i__1; ++k) {
2042 /* L270: */
2043             sum += bmat[jp + k * bmat_dim1] * d__[k];
2044         }
2045         vlag[jp] = sum;
2046         bsum += sum * d__[j];
2047 /* L280: */
2048         dx += d__[j] * xopt[j];
2049     }
2050     beta = dx * dx + dsq * (xoptsq + dx + dx + half * dsq) + beta - bsum;
2051     vlag[kopt] += one;
2052
2053 /* If KNEW is positive and if the cancellation in DENOM is unacceptable, */
2054 /* then BIGDEN calculates an alternative model step, XNEW being used for */
2055 /* working space. */
2056
2057     if (knew > 0) {
2058 /* Computing 2nd power */
2059         d__1 = vlag[knew];
2060         temp = one + alpha * beta / (d__1 * d__1);
2061         if (fabs(temp) <= .8) {
2062             bigden_(n, npt, &xopt[1], &xpt[xpt_offset], &bmat[bmat_offset], &
2063                     zmat[zmat_offset], &idz, ndim, &kopt, &knew, &d__[1], &w[
2064                     1], &vlag[1], &beta, &xnew[1], &w[*ndim + 1], &w[*ndim *
2065                     6 + 1], &xbase[1], lb, ub);
2066         }
2067     }
2068
2069 /* Calculate the next value of the objective function. */
2070
2071 L290:
2072     i__2 = *n;
2073     for (i__ = 1; i__ <= i__2; ++i__) {
2074         xnew[i__] = xopt[i__] + d__[i__];
2075 /* L300: */
2076         x[i__] = xbase[i__] + xnew[i__];
2077     }
2078     ++nf;
2079 L310:
2080     if (nlopt_stop_evals(stop)) rc = NLOPT_MAXEVAL_REACHED;
2081     else if (nlopt_stop_time(stop)) rc = NLOPT_MAXTIME_REACHED;
2082     if (rc != NLOPT_SUCCESS) goto L530;
2083
2084     stop->nevals++;
2085     f = calfun(*n, &x[1], calfun_data);
2086     if (f < stop->minf_max) {
2087        rc = NLOPT_MINF_MAX_REACHED;
2088        goto L530;
2089     }
2090
2091     if (nf <= *npt) {
2092         goto L70;
2093     }
2094     if (knew == -1) {
2095         goto L530;
2096     }
2097
2098 /* Use the quadratic model to predict the change in F due to the step D, */
2099 /* and set DIFF to the error of this prediction. */
2100
2101     vquad = zero;
2102     ih = 0;
2103     i__2 = *n;
2104     for (j = 1; j <= i__2; ++j) {
2105         vquad += d__[j] * gq[j];
2106         i__1 = j;
2107         for (i__ = 1; i__ <= i__1; ++i__) {
2108             ++ih;
2109             temp = d__[i__] * xnew[j] + d__[j] * xopt[i__];
2110             if (i__ == j) {
2111                 temp = half * temp;
2112             }
2113 /* L340: */
2114             vquad += temp * hq[ih];
2115         }
2116     }
2117     i__1 = *npt;
2118     for (k = 1; k <= i__1; ++k) {
2119 /* L350: */
2120         vquad += pq[k] * w[k];
2121     }
2122     diff = f - fopt - vquad;
2123     diffc = diffb;
2124     diffb = diffa;
2125     diffa = fabs(diff);
2126     if (dnorm > rho) {
2127         nfsav = nf;
2128     }
2129
2130 /* Update FOPT and XOPT if the new F is the least value of the objective */
2131 /* function so far. The branch when KNEW is positive occurs if D is not */
2132 /* a trust region step. */
2133
2134     fsave = fopt;
2135     if (f < fopt) {
2136         fopt = f;
2137         xoptsq = zero;
2138         i__1 = *n;
2139         for (i__ = 1; i__ <= i__1; ++i__) {
2140             xopt[i__] = xnew[i__];
2141 /* L360: */
2142 /* Computing 2nd power */
2143             d__1 = xopt[i__];
2144             xoptsq += d__1 * d__1;
2145         }
2146         if (nlopt_stop_ftol(stop, fopt, fsave)) {
2147              rc = NLOPT_FTOL_REACHED;
2148              goto L530;
2149         }
2150
2151     }
2152     ksave = knew;
2153     if (knew > 0) {
2154         goto L410;
2155     }
2156
2157 /* Pick the next value of DELTA after a trust region step. */
2158
2159     if (vquad >= zero) {
2160         goto L530;
2161     }
2162     ratio = (f - fsave) / vquad;
2163     if (ratio <= tenth) {
2164         delta = half * dnorm;
2165     } else if (ratio <= .7) {
2166 /* Computing MAX */
2167         d__1 = half * delta;
2168         delta = max(d__1,dnorm);
2169     } else {
2170 /* Computing MAX */
2171         d__1 = half * delta, d__2 = dnorm + dnorm;
2172         delta = max(d__1,d__2);
2173     }
2174     if (delta <= rho * 1.5) {
2175         delta = rho;
2176     }
2177
2178 /* Set KNEW to the index of the next interpolation point to be deleted. */
2179
2180 /* Computing MAX */
2181     d__2 = tenth * delta;
2182 /* Computing 2nd power */
2183     d__1 = max(d__2,rho);
2184     rhosq = d__1 * d__1;
2185     ktemp = 0;
2186     detrat = zero;
2187     if (f >= fsave) {
2188         ktemp = kopt;
2189         detrat = one;
2190     }
2191     i__1 = *npt;
2192     for (k = 1; k <= i__1; ++k) {
2193         hdiag = zero;
2194         i__2 = nptm;
2195         for (j = 1; j <= i__2; ++j) {
2196             temp = one;
2197             if (j < idz) {
2198                 temp = -one;
2199             }
2200 /* L380: */
2201 /* Computing 2nd power */
2202             d__1 = zmat[k + j * zmat_dim1];
2203             hdiag += temp * (d__1 * d__1);
2204         }
2205 /* Computing 2nd power */
2206         d__2 = vlag[k];
2207         temp = (d__1 = beta * hdiag + d__2 * d__2, fabs(d__1));
2208         distsq = zero;
2209         i__2 = *n;
2210         for (j = 1; j <= i__2; ++j) {
2211 /* L390: */
2212 /* Computing 2nd power */
2213             d__1 = xpt[k + j * xpt_dim1] - xopt[j];
2214             distsq += d__1 * d__1;
2215         }
2216         if (distsq > rhosq) {
2217 /* Computing 3rd power */
2218             d__1 = distsq / rhosq;
2219             temp *= d__1 * (d__1 * d__1);
2220         }
2221         if (temp > detrat && k != ktemp) {
2222             detrat = temp;
2223             knew = k;
2224         }
2225 /* L400: */
2226     }
2227     if (knew == 0) {
2228         goto L460;
2229     }
2230
2231 /* Update BMAT, ZMAT and IDZ, so that the KNEW-th interpolation point */
2232 /* can be moved. Begin the updating of the quadratic model, starting */
2233 /* with the explicit second derivative term. */
2234
2235 L410:
2236     update_(n, npt, &bmat[bmat_offset], &zmat[zmat_offset], &idz, ndim, &vlag[
2237             1], &beta, &knew, &w[1]);
2238     fval[knew] = f;
2239     ih = 0;
2240     i__1 = *n;
2241     for (i__ = 1; i__ <= i__1; ++i__) {
2242         temp = pq[knew] * xpt[knew + i__ * xpt_dim1];
2243         i__2 = i__;
2244         for (j = 1; j <= i__2; ++j) {
2245             ++ih;
2246 /* L420: */
2247             hq[ih] += temp * xpt[knew + j * xpt_dim1];
2248         }
2249     }
2250     pq[knew] = zero;
2251
2252 /* Update the other second derivative parameters, and then the gradient */
2253 /* vector of the model. Also include the new interpolation point. */
2254
2255     i__2 = nptm;
2256     for (j = 1; j <= i__2; ++j) {
2257         temp = diff * zmat[knew + j * zmat_dim1];
2258         if (j < idz) {
2259             temp = -temp;
2260         }
2261         i__1 = *npt;
2262         for (k = 1; k <= i__1; ++k) {
2263 /* L440: */
2264             pq[k] += temp * zmat[k + j * zmat_dim1];
2265         }
2266     }
2267     gqsq = zero;
2268     i__1 = *n;
2269     for (i__ = 1; i__ <= i__1; ++i__) {
2270         gq[i__] += diff * bmat[knew + i__ * bmat_dim1];
2271 /* Computing 2nd power */
2272         d__1 = gq[i__];
2273         gqsq += d__1 * d__1;
2274 /* L450: */
2275         xpt[knew + i__ * xpt_dim1] = xnew[i__];
2276     }
2277
2278 /* If a trust region step makes a small change to the objective function, */
2279 /* then calculate the gradient of the least Frobenius norm interpolant at */
2280 /* XBASE, and store it in W, using VLAG for a vector of right hand sides. */
2281
2282     if (ksave == 0 && delta == rho) {
2283         if (fabs(ratio) > .01) {
2284             itest = 0;
2285         } else {
2286             i__1 = *npt;
2287             for (k = 1; k <= i__1; ++k) {
2288 /* L700: */
2289                 vlag[k] = fval[k] - fval[kopt];
2290             }
2291             gisq = zero;
2292             i__1 = *n;
2293             for (i__ = 1; i__ <= i__1; ++i__) {
2294                 sum = zero;
2295                 i__2 = *npt;
2296                 for (k = 1; k <= i__2; ++k) {
2297 /* L710: */
2298                     sum += bmat[k + i__ * bmat_dim1] * vlag[k];
2299                 }
2300                 gisq += sum * sum;
2301 /* L720: */
2302                 w[i__] = sum;
2303             }
2304
2305 /* Test whether to replace the new quadratic model by the least Frobenius */
2306 /* norm interpolant, making the replacement if the test is satisfied. */
2307
2308             ++itest;
2309             if (gqsq < gisq * 100.) {
2310                 itest = 0;
2311             }
2312             if (itest >= 3) {
2313                 i__1 = *n;
2314                 for (i__ = 1; i__ <= i__1; ++i__) {
2315 /* L730: */
2316                     gq[i__] = w[i__];
2317                 }
2318                 i__1 = nh;
2319                 for (ih = 1; ih <= i__1; ++ih) {
2320 /* L740: */
2321                     hq[ih] = zero;
2322                 }
2323                 i__1 = nptm;
2324                 for (j = 1; j <= i__1; ++j) {
2325                     w[j] = zero;
2326                     i__2 = *npt;
2327                     for (k = 1; k <= i__2; ++k) {
2328 /* L750: */
2329                         w[j] += vlag[k] * zmat[k + j * zmat_dim1];
2330                     }
2331 /* L760: */
2332                     if (j < idz) {
2333                         w[j] = -w[j];
2334                     }
2335                 }
2336                 i__1 = *npt;
2337                 for (k = 1; k <= i__1; ++k) {
2338                     pq[k] = zero;
2339                     i__2 = nptm;
2340                     for (j = 1; j <= i__2; ++j) {
2341 /* L770: */
2342                         pq[k] += zmat[k + j * zmat_dim1] * w[j];
2343                     }
2344                 }
2345                 itest = 0;
2346             }
2347         }
2348     }
2349     if (f < fsave) {
2350         kopt = knew;
2351     }
2352
2353 /* If a trust region step has provided a sufficient decrease in F, then */
2354 /* branch for another trust region calculation. The case KSAVE>0 occurs */
2355 /* when the new function value was calculated by a model step. */
2356
2357     if (f <= fsave + tenth * vquad) {
2358         goto L100;
2359     }
2360     if (ksave > 0) {
2361         goto L100;
2362     }
2363
2364 /* Alternatively, find out if the interpolation points are close enough */
2365 /* to the best point so far. */
2366
2367     knew = 0;
2368 L460:
2369     distsq = delta * 4. * delta;
2370     i__2 = *npt;
2371     for (k = 1; k <= i__2; ++k) {
2372         sum = zero;
2373         i__1 = *n;
2374         for (j = 1; j <= i__1; ++j) {
2375 /* L470: */
2376 /* Computing 2nd power */
2377             d__1 = xpt[k + j * xpt_dim1] - xopt[j];
2378             sum += d__1 * d__1;
2379         }
2380         if (sum > distsq) {
2381             knew = k;
2382             distsq = sum;
2383         }
2384 /* L480: */
2385     }
2386
2387 /* If KNEW is positive, then set DSTEP, and branch back for the next */
2388 /* iteration, which will generate a "model step". */
2389
2390     if (knew > 0) {
2391 /* Computing MAX */
2392 /* Computing MIN */
2393         d__2 = tenth * sqrt(distsq), d__3 = half * delta;
2394         d__1 = min(d__2,d__3);
2395         dstep = max(d__1,rho);
2396         dsq = dstep * dstep;
2397         goto L120;
2398     }
2399     if (ratio > zero) {
2400         goto L100;
2401     }
2402     if (max(delta,dnorm) > rho) {
2403         goto L100;
2404     }
2405
2406 /* The calculations with the current value of RHO are complete. Pick the */
2407 /* next values of RHO and DELTA. */
2408
2409 L490:
2410     if (rho > rhoend) {
2411         delta = half * rho;
2412         ratio = rho / rhoend;
2413         if (ratio <= 16.) {
2414             rho = rhoend;
2415         } else if (ratio <= 250.) {
2416             rho = sqrt(ratio) * rhoend;
2417         } else {
2418             rho = tenth * rho;
2419         }
2420         delta = max(delta,rho);
2421         goto L90;
2422     }
2423
2424 /* Return from the calculation, after another Newton-Raphson step, if */
2425 /* it is too short to have been tried before. */
2426
2427     if (knew == -1) {
2428         goto L290;
2429     }
2430     rc = NLOPT_XTOL_REACHED;
2431 L530:
2432     if (fopt <= f) {
2433         i__2 = *n;
2434         for (i__ = 1; i__ <= i__2; ++i__) {
2435 /* L540: */
2436             x[i__] = xbase[i__] + xopt[i__];
2437         }
2438         f = fopt;
2439     }
2440     *minf = f;
2441     return rc;
2442 } /* newuob_ */
2443
2444 /*************************************************************************/
2445 /* newuoa.f */
2446
2447 nlopt_result newuoa(int n, int npt, double *x,
2448                     const double *lb, const double *ub,
2449                     double rhobeg, nlopt_stopping *stop, double *minf,
2450                     newuoa_func calfun, void *calfun_data)
2451 {
2452     /* Local variables */
2453     int id, np, iw, igq, ihq, ixb, ifv, ipq, ivl, ixn, ixo, ixp, ndim,
2454             nptm, ibmat, izmat;
2455     nlopt_result ret;
2456     double *w;
2457
2458 /* This subroutine seeks the least value of a function of many variables, */
2459 /* by a trust region method that forms quadratic models by interpolation. */
2460 /* There can be some freedom in the interpolation conditions, which is */
2461 /* taken up by minimizing the Frobenius norm of the change to the second */
2462 /* derivative of the quadratic model, beginning with a zero matrix. The */
2463 /* arguments of the subroutine are as follows. */
2464
2465 /* N must be set to the number of variables and must be at least two. */
2466 /* NPT is the number of interpolation conditions. Its value must be in the */
2467 /*   interval [N+2,(N+1)(N+2)/2]. */
2468 /* Initial values of the variables must be set in X(1),X(2),...,X(N). They */
2469 /*   will be changed to the values that give the least calculated F. */
2470 /* RHOBEG and RHOEND must be set to the initial and final values of a trust */
2471 /*   region radius, so both must be positive with RHOEND<=RHOBEG. Typically */
2472 /*   RHOBEG should be about one tenth of the greatest expected change to a */
2473 /*   variable, and RHOEND should indicate the accuracy that is required in */
2474 /*   the final values of the variables. */
2475 /* The value of IPRINT should be set to 0, 1, 2 or 3, which controls the */
2476 /*   amount of printing. Specifically, there is no output if IPRINT=0 and */
2477 /*   there is output only at the return if IPRINT=1. Otherwise, each new */
2478 /*   value of RHO is printed, with the best vector of variables so far and */
2479 /*   the corresponding value of the objective function. Further, each new */
2480 /*   value of F with its variables are output if IPRINT=3. */
2481 /* MAXFUN must be set to an upper bound on the number of calls of CALFUN. */
2482 /* The array W will be used for working space. Its length must be at least */
2483 /* (NPT+13)*(NPT+N)+3*N*(N+3)/2. */
2484
2485 /* SUBROUTINE CALFUN (N,X,F) must be provided by the user. It must set F to */
2486 /* the value of the objective function for the variables X(1),X(2),...,X(N). */
2487
2488 /* Partition the working space array, so that different parts of it can be */
2489 /* treated separately by the subroutine that performs the main calculation. */
2490
2491     /* Parameter adjustments */
2492     --x;
2493
2494     /* Function Body */
2495     np = n + 1;
2496     nptm = npt - np;
2497     if (n < 2 || npt < n + 2 || npt > (n + 2) * np / 2) {
2498          return NLOPT_INVALID_ARGS;
2499     }
2500     ndim = npt + n;
2501     ixb = 1;
2502     ixo = ixb + n;
2503     ixn = ixo + n;
2504     ixp = ixn + n;
2505     ifv = ixp + n * npt;
2506     igq = ifv + npt;
2507     ihq = igq + n;
2508     ipq = ihq + n * np / 2;
2509     ibmat = ipq + npt;
2510     izmat = ibmat + ndim * n;
2511     id = izmat + npt * nptm;
2512     ivl = id + n;
2513     iw = ivl + ndim;
2514
2515     w = (double *) malloc(sizeof(double) * ((npt+13)*(npt+n) + 3*(n*(n+3))/2));
2516     if (!w) return NLOPT_OUT_OF_MEMORY;
2517     --w;
2518
2519 /* The above settings provide a partition of W for subroutine NEWUOB. */
2520 /* The partition requires the first NPT*(NPT+N)+5*N*(N+3)/2 elements of */
2521 /* W plus the space that is needed by the last array of NEWUOB. */
2522
2523     ret = newuob_(&n, &npt, &x[1], &rhobeg,
2524                   lb, ub, stop, minf, calfun, calfun_data,
2525                   &w[ixb], &w[ixo], &w[ixn], &w[ixp], &w[ifv],
2526                   &w[igq], &w[ihq], &w[ipq], &w[ibmat], &w[izmat],
2527                   &ndim, &w[id], &w[ivl], &w[iw]);
2528
2529     ++w;
2530     free(w);
2531     return ret;
2532 } /* newuoa_ */
2533
2534 /*************************************************************************/
2535 /*************************************************************************/