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     if (lb && ub) { /* SGJ, 2008: make sure we are within bounds */
1748          for (j = 1; j <= i__1; ++j) {
1749               if (x[j] < lb[j-1]) x[j] = lb[j-1];
1750               else if (x[j] > ub[j-1]) x[j] = ub[j-1];
1751          }
1752     }
1753     goto L310;
1754 L70:
1755     fval[nf] = f;
1756     if (nf == 1) {
1757         fbeg = f;
1758         fopt = f;
1759         kopt = 1;
1760     } else if (f < fopt) {
1761         fopt = f;
1762         kopt = nf;
1763     }
1764
1765 /* Set the nonzero initial elements of BMAT and the quadratic model in */
1766 /* the cases when NF is at most 2*N+1. */
1767
1768     if (nfm <= *n << 1) {
1769         if (nfm >= 1 && nfm <= *n) {
1770             gq[nfm] = (f - fbeg) / *rhobeg;
1771             if (*npt < nf + *n) {
1772                 bmat[nfm * bmat_dim1 + 1] = -one / *rhobeg;
1773                 bmat[nf + nfm * bmat_dim1] = one / *rhobeg;
1774                 bmat[*npt + nfm + nfm * bmat_dim1] = -half * rhosq;
1775             }
1776         } else if (nfm > *n) {
1777             bmat[nf - *n + nfmm * bmat_dim1] = half / *rhobeg;
1778             bmat[nf + nfmm * bmat_dim1] = -half / *rhobeg;
1779             zmat[nfmm * zmat_dim1 + 1] = -reciq - reciq;
1780             zmat[nf - *n + nfmm * zmat_dim1] = reciq;
1781             zmat[nf + nfmm * zmat_dim1] = reciq;
1782             ih = nfmm * (nfmm + 1) / 2;
1783             temp = (fbeg - f) / *rhobeg;
1784             hq[ih] = (gq[nfmm] - temp) / *rhobeg;
1785             gq[nfmm] = half * (gq[nfmm] + temp);
1786         }
1787
1788 /* Set the off-diagonal second derivatives of the Lagrange functions and */
1789 /* the initial quadratic model. */
1790
1791     } else {
1792         ih = ipt * (ipt - 1) / 2 + jpt;
1793         if (xipt < zero) {
1794             ipt += *n;
1795         }
1796         if (xjpt < zero) {
1797             jpt += *n;
1798         }
1799         zmat[nfmm * zmat_dim1 + 1] = recip;
1800         zmat[nf + nfmm * zmat_dim1] = recip;
1801         zmat[ipt + 1 + nfmm * zmat_dim1] = -recip;
1802         zmat[jpt + 1 + nfmm * zmat_dim1] = -recip;
1803         hq[ih] = (fbeg - fval[ipt + 1] - fval[jpt + 1] + f) / (xipt * xjpt);
1804     }
1805     if (nf < *npt) {
1806         goto L50;
1807     }
1808
1809 /* Begin the iterative procedure, because the initial model is complete. */
1810
1811     rho = *rhobeg;
1812     delta = rho;
1813     idz = 1;
1814     diffa = zero;
1815     diffb = zero;
1816     itest = 0;
1817     xoptsq = zero;
1818     i__1 = *n;
1819     for (i__ = 1; i__ <= i__1; ++i__) {
1820         xopt[i__] = xpt[kopt + i__ * xpt_dim1];
1821 /* L80: */
1822 /* Computing 2nd power */
1823         d__1 = xopt[i__];
1824         xoptsq += d__1 * d__1;
1825     }
1826 L90:
1827     nfsav = nf;
1828
1829 /* Generate the next trust region step and test its length. Set KNEW */
1830 /* to -1 if the purpose of the next F will be to improve the model. */
1831
1832 L100:
1833     knew = 0;
1834     rc2 = trsapp_(n, npt, &xopt[1], &xpt[xpt_offset], &gq[1], &hq[1], &pq[1], &
1835             delta, &d__[1], &w[1], &w[np], &w[np + *n], &w[np + (*n << 1)], &
1836             crvmin, &xbase[1], lb, ub);
1837     if (rc2 < 0) { rc = rc2; goto L530; }
1838     dsq = zero;
1839     i__1 = *n;
1840     for (i__ = 1; i__ <= i__1; ++i__) {
1841 /* L110: */
1842 /* Computing 2nd power */
1843         d__1 = d__[i__];
1844         dsq += d__1 * d__1;
1845     }
1846 /* Computing MIN */
1847     d__1 = delta, d__2 = sqrt(dsq);
1848     dnorm = min(d__1,d__2);
1849     if (dnorm < half * rho) {
1850         knew = -1;
1851         delta = tenth * delta;
1852         ratio = -1.;
1853         if (delta <= rho * 1.5) {
1854             delta = rho;
1855         }
1856         if (nf <= nfsav + 2) {
1857             goto L460;
1858         }
1859         temp = crvmin * .125 * rho * rho;
1860 /* Computing MAX */
1861         d__1 = max(diffa,diffb);
1862         if (temp <= max(d__1,diffc)) {
1863             goto L460;
1864         }
1865         goto L490;
1866     }
1867
1868 /* Shift XBASE if XOPT may be too far from XBASE. First make the changes */
1869 /* to BMAT that do not depend on ZMAT. */
1870
1871 L120:
1872     if (dsq <= xoptsq * .001) {
1873         tempq = xoptsq * .25;
1874         i__1 = *npt;
1875         for (k = 1; k <= i__1; ++k) {
1876             sum = zero;
1877             i__2 = *n;
1878             for (i__ = 1; i__ <= i__2; ++i__) {
1879 /* L130: */
1880                 sum += xpt[k + i__ * xpt_dim1] * xopt[i__];
1881             }
1882             temp = pq[k] * sum;
1883             sum -= half * xoptsq;
1884             w[*npt + k] = sum;
1885             i__2 = *n;
1886             for (i__ = 1; i__ <= i__2; ++i__) {
1887                 gq[i__] += temp * xpt[k + i__ * xpt_dim1];
1888                 xpt[k + i__ * xpt_dim1] -= half * xopt[i__];
1889                 vlag[i__] = bmat[k + i__ * bmat_dim1];
1890                 w[i__] = sum * xpt[k + i__ * xpt_dim1] + tempq * xopt[i__];
1891                 ip = *npt + i__;
1892                 i__3 = i__;
1893                 for (j = 1; j <= i__3; ++j) {
1894 /* L140: */
1895                     bmat[ip + j * bmat_dim1] = bmat[ip + j * bmat_dim1] +
1896                             vlag[i__] * w[j] + w[i__] * vlag[j];
1897                 }
1898             }
1899         }
1900
1901 /* Then the revisions of BMAT that depend on ZMAT are calculated. */
1902
1903         i__3 = nptm;
1904         for (k = 1; k <= i__3; ++k) {
1905             sumz = zero;
1906             i__2 = *npt;
1907             for (i__ = 1; i__ <= i__2; ++i__) {
1908                 sumz += zmat[i__ + k * zmat_dim1];
1909 /* L150: */
1910                 w[i__] = w[*npt + i__] * zmat[i__ + k * zmat_dim1];
1911             }
1912             i__2 = *n;
1913             for (j = 1; j <= i__2; ++j) {
1914                 sum = tempq * sumz * xopt[j];
1915                 i__1 = *npt;
1916                 for (i__ = 1; i__ <= i__1; ++i__) {
1917 /* L160: */
1918                     sum += w[i__] * xpt[i__ + j * xpt_dim1];
1919                 }
1920                 vlag[j] = sum;
1921                 if (k < idz) {
1922                     sum = -sum;
1923                 }
1924                 i__1 = *npt;
1925                 for (i__ = 1; i__ <= i__1; ++i__) {
1926 /* L170: */
1927                     bmat[i__ + j * bmat_dim1] += sum * zmat[i__ + k *
1928                             zmat_dim1];
1929                 }
1930             }
1931             i__1 = *n;
1932             for (i__ = 1; i__ <= i__1; ++i__) {
1933                 ip = i__ + *npt;
1934                 temp = vlag[i__];
1935                 if (k < idz) {
1936                     temp = -temp;
1937                 }
1938                 i__2 = i__;
1939                 for (j = 1; j <= i__2; ++j) {
1940 /* L180: */
1941                     bmat[ip + j * bmat_dim1] += temp * vlag[j];
1942                 }
1943             }
1944         }
1945
1946 /* The following instructions complete the shift of XBASE, including */
1947 /* the changes to the parameters of the quadratic model. */
1948
1949         ih = 0;
1950         i__2 = *n;
1951         for (j = 1; j <= i__2; ++j) {
1952             w[j] = zero;
1953             i__1 = *npt;
1954             for (k = 1; k <= i__1; ++k) {
1955                 w[j] += pq[k] * xpt[k + j * xpt_dim1];
1956 /* L190: */
1957                 xpt[k + j * xpt_dim1] -= half * xopt[j];
1958             }
1959             i__1 = j;
1960             for (i__ = 1; i__ <= i__1; ++i__) {
1961                 ++ih;
1962                 if (i__ < j) {
1963                     gq[j] += hq[ih] * xopt[i__];
1964                 }
1965                 gq[i__] += hq[ih] * xopt[j];
1966                 hq[ih] = hq[ih] + w[i__] * xopt[j] + xopt[i__] * w[j];
1967 /* L200: */
1968                 bmat[*npt + i__ + j * bmat_dim1] = bmat[*npt + j + i__ *
1969                         bmat_dim1];
1970             }
1971         }
1972         i__1 = *n;
1973         for (j = 1; j <= i__1; ++j) {
1974             xbase[j] += xopt[j];
1975 /* L210: */
1976             xopt[j] = zero;
1977         }
1978         xoptsq = zero;
1979     }
1980
1981 /* Pick the model step if KNEW is positive. A different choice of D */
1982 /* may be made later, if the choice of D by BIGLAG causes substantial */
1983 /* cancellation in DENOM. */
1984
1985     if (knew > 0) {
1986        rc2 =
1987          biglag_(n, npt, &xopt[1], &xpt[xpt_offset], &bmat[bmat_offset], &zmat[
1988                 zmat_offset], &idz, ndim, &knew, &dstep, &d__[1], &alpha, &
1989                 vlag[1], &vlag[*npt + 1], &w[1], &w[np], &w[np + *n],
1990                 &xbase[1], lb, ub);
1991        if (rc2 < 0) { rc = rc2; goto L530; }
1992     }
1993
1994 /* Calculate VLAG and BETA for the current choice of D. The first NPT */
1995 /* components of W_check will be held in W. */
1996
1997     i__1 = *npt;
1998     for (k = 1; k <= i__1; ++k) {
1999         suma = zero;
2000         sumb = zero;
2001         sum = zero;
2002         i__2 = *n;
2003         for (j = 1; j <= i__2; ++j) {
2004             suma += xpt[k + j * xpt_dim1] * d__[j];
2005             sumb += xpt[k + j * xpt_dim1] * xopt[j];
2006 /* L220: */
2007             sum += bmat[k + j * bmat_dim1] * d__[j];
2008         }
2009         w[k] = suma * (half * suma + sumb);
2010 /* L230: */
2011         vlag[k] = sum;
2012     }
2013     beta = zero;
2014     i__1 = nptm;
2015     for (k = 1; k <= i__1; ++k) {
2016         sum = zero;
2017         i__2 = *npt;
2018         for (i__ = 1; i__ <= i__2; ++i__) {
2019 /* L240: */
2020             sum += zmat[i__ + k * zmat_dim1] * w[i__];
2021         }
2022         if (k < idz) {
2023             beta += sum * sum;
2024             sum = -sum;
2025         } else {
2026             beta -= sum * sum;
2027         }
2028         i__2 = *npt;
2029         for (i__ = 1; i__ <= i__2; ++i__) {
2030 /* L250: */
2031             vlag[i__] += sum * zmat[i__ + k * zmat_dim1];
2032         }
2033     }
2034     bsum = zero;
2035     dx = zero;
2036     i__2 = *n;
2037     for (j = 1; j <= i__2; ++j) {
2038         sum = zero;
2039         i__1 = *npt;
2040         for (i__ = 1; i__ <= i__1; ++i__) {
2041 /* L260: */
2042             sum += w[i__] * bmat[i__ + j * bmat_dim1];
2043         }
2044         bsum += sum * d__[j];
2045         jp = *npt + j;
2046         i__1 = *n;
2047         for (k = 1; k <= i__1; ++k) {
2048 /* L270: */
2049             sum += bmat[jp + k * bmat_dim1] * d__[k];
2050         }
2051         vlag[jp] = sum;
2052         bsum += sum * d__[j];
2053 /* L280: */
2054         dx += d__[j] * xopt[j];
2055     }
2056     beta = dx * dx + dsq * (xoptsq + dx + dx + half * dsq) + beta - bsum;
2057     vlag[kopt] += one;
2058
2059 /* If KNEW is positive and if the cancellation in DENOM is unacceptable, */
2060 /* then BIGDEN calculates an alternative model step, XNEW being used for */
2061 /* working space. */
2062
2063     if (knew > 0) {
2064 /* Computing 2nd power */
2065         d__1 = vlag[knew];
2066         temp = one + alpha * beta / (d__1 * d__1);
2067         if (fabs(temp) <= .8) {
2068             bigden_(n, npt, &xopt[1], &xpt[xpt_offset], &bmat[bmat_offset], &
2069                     zmat[zmat_offset], &idz, ndim, &kopt, &knew, &d__[1], &w[
2070                     1], &vlag[1], &beta, &xnew[1], &w[*ndim + 1], &w[*ndim *
2071                     6 + 1], &xbase[1], lb, ub);
2072         }
2073     }
2074
2075 /* Calculate the next value of the objective function. */
2076
2077 L290:
2078     i__2 = *n;
2079     for (i__ = 1; i__ <= i__2; ++i__) {
2080         xnew[i__] = xopt[i__] + d__[i__];
2081 /* L300: */
2082         x[i__] = xbase[i__] + xnew[i__];
2083     }
2084     ++nf;
2085 L310:
2086     if (nlopt_stop_evals(stop)) rc = NLOPT_MAXEVAL_REACHED;
2087     else if (nlopt_stop_time(stop)) rc = NLOPT_MAXTIME_REACHED;
2088     if (rc != NLOPT_SUCCESS) goto L530;
2089
2090     stop->nevals++;
2091     f = calfun(*n, &x[1], calfun_data);
2092     if (f < stop->minf_max) {
2093        rc = NLOPT_MINF_MAX_REACHED;
2094        goto L530;
2095     }
2096
2097     if (nf <= *npt) {
2098         goto L70;
2099     }
2100     if (knew == -1) {
2101         goto L530;
2102     }
2103
2104 /* Use the quadratic model to predict the change in F due to the step D, */
2105 /* and set DIFF to the error of this prediction. */
2106
2107     vquad = zero;
2108     ih = 0;
2109     i__2 = *n;
2110     for (j = 1; j <= i__2; ++j) {
2111         vquad += d__[j] * gq[j];
2112         i__1 = j;
2113         for (i__ = 1; i__ <= i__1; ++i__) {
2114             ++ih;
2115             temp = d__[i__] * xnew[j] + d__[j] * xopt[i__];
2116             if (i__ == j) {
2117                 temp = half * temp;
2118             }
2119 /* L340: */
2120             vquad += temp * hq[ih];
2121         }
2122     }
2123     i__1 = *npt;
2124     for (k = 1; k <= i__1; ++k) {
2125 /* L350: */
2126         vquad += pq[k] * w[k];
2127     }
2128     diff = f - fopt - vquad;
2129     diffc = diffb;
2130     diffb = diffa;
2131     diffa = fabs(diff);
2132     if (dnorm > rho) {
2133         nfsav = nf;
2134     }
2135
2136 /* Update FOPT and XOPT if the new F is the least value of the objective */
2137 /* function so far. The branch when KNEW is positive occurs if D is not */
2138 /* a trust region step. */
2139
2140     fsave = fopt;
2141     if (f < fopt) {
2142         fopt = f;
2143         xoptsq = zero;
2144         i__1 = *n;
2145         for (i__ = 1; i__ <= i__1; ++i__) {
2146             xopt[i__] = xnew[i__];
2147 /* L360: */
2148 /* Computing 2nd power */
2149             d__1 = xopt[i__];
2150             xoptsq += d__1 * d__1;
2151         }
2152         if (nlopt_stop_ftol(stop, fopt, fsave)) {
2153              rc = NLOPT_FTOL_REACHED;
2154              goto L530;
2155         }
2156
2157     }
2158     ksave = knew;
2159     if (knew > 0) {
2160         goto L410;
2161     }
2162
2163 /* Pick the next value of DELTA after a trust region step. */
2164
2165     if (vquad >= zero) {
2166         goto L530;
2167     }
2168     ratio = (f - fsave) / vquad;
2169     if (ratio <= tenth) {
2170         delta = half * dnorm;
2171     } else if (ratio <= .7) {
2172 /* Computing MAX */
2173         d__1 = half * delta;
2174         delta = max(d__1,dnorm);
2175     } else {
2176 /* Computing MAX */
2177         d__1 = half * delta, d__2 = dnorm + dnorm;
2178         delta = max(d__1,d__2);
2179     }
2180     if (delta <= rho * 1.5) {
2181         delta = rho;
2182     }
2183
2184 /* Set KNEW to the index of the next interpolation point to be deleted. */
2185
2186 /* Computing MAX */
2187     d__2 = tenth * delta;
2188 /* Computing 2nd power */
2189     d__1 = max(d__2,rho);
2190     rhosq = d__1 * d__1;
2191     ktemp = 0;
2192     detrat = zero;
2193     if (f >= fsave) {
2194         ktemp = kopt;
2195         detrat = one;
2196     }
2197     i__1 = *npt;
2198     for (k = 1; k <= i__1; ++k) {
2199         hdiag = zero;
2200         i__2 = nptm;
2201         for (j = 1; j <= i__2; ++j) {
2202             temp = one;
2203             if (j < idz) {
2204                 temp = -one;
2205             }
2206 /* L380: */
2207 /* Computing 2nd power */
2208             d__1 = zmat[k + j * zmat_dim1];
2209             hdiag += temp * (d__1 * d__1);
2210         }
2211 /* Computing 2nd power */
2212         d__2 = vlag[k];
2213         temp = (d__1 = beta * hdiag + d__2 * d__2, fabs(d__1));
2214         distsq = zero;
2215         i__2 = *n;
2216         for (j = 1; j <= i__2; ++j) {
2217 /* L390: */
2218 /* Computing 2nd power */
2219             d__1 = xpt[k + j * xpt_dim1] - xopt[j];
2220             distsq += d__1 * d__1;
2221         }
2222         if (distsq > rhosq) {
2223 /* Computing 3rd power */
2224             d__1 = distsq / rhosq;
2225             temp *= d__1 * (d__1 * d__1);
2226         }
2227         if (temp > detrat && k != ktemp) {
2228             detrat = temp;
2229             knew = k;
2230         }
2231 /* L400: */
2232     }
2233     if (knew == 0) {
2234         goto L460;
2235     }
2236
2237 /* Update BMAT, ZMAT and IDZ, so that the KNEW-th interpolation point */
2238 /* can be moved. Begin the updating of the quadratic model, starting */
2239 /* with the explicit second derivative term. */
2240
2241 L410:
2242     update_(n, npt, &bmat[bmat_offset], &zmat[zmat_offset], &idz, ndim, &vlag[
2243             1], &beta, &knew, &w[1]);
2244     fval[knew] = f;
2245     ih = 0;
2246     i__1 = *n;
2247     for (i__ = 1; i__ <= i__1; ++i__) {
2248         temp = pq[knew] * xpt[knew + i__ * xpt_dim1];
2249         i__2 = i__;
2250         for (j = 1; j <= i__2; ++j) {
2251             ++ih;
2252 /* L420: */
2253             hq[ih] += temp * xpt[knew + j * xpt_dim1];
2254         }
2255     }
2256     pq[knew] = zero;
2257
2258 /* Update the other second derivative parameters, and then the gradient */
2259 /* vector of the model. Also include the new interpolation point. */
2260
2261     i__2 = nptm;
2262     for (j = 1; j <= i__2; ++j) {
2263         temp = diff * zmat[knew + j * zmat_dim1];
2264         if (j < idz) {
2265             temp = -temp;
2266         }
2267         i__1 = *npt;
2268         for (k = 1; k <= i__1; ++k) {
2269 /* L440: */
2270             pq[k] += temp * zmat[k + j * zmat_dim1];
2271         }
2272     }
2273     gqsq = zero;
2274     i__1 = *n;
2275     for (i__ = 1; i__ <= i__1; ++i__) {
2276         gq[i__] += diff * bmat[knew + i__ * bmat_dim1];
2277 /* Computing 2nd power */
2278         d__1 = gq[i__];
2279         gqsq += d__1 * d__1;
2280 /* L450: */
2281         xpt[knew + i__ * xpt_dim1] = xnew[i__];
2282     }
2283
2284 /* If a trust region step makes a small change to the objective function, */
2285 /* then calculate the gradient of the least Frobenius norm interpolant at */
2286 /* XBASE, and store it in W, using VLAG for a vector of right hand sides. */
2287
2288     if (ksave == 0 && delta == rho) {
2289         if (fabs(ratio) > .01) {
2290             itest = 0;
2291         } else {
2292             i__1 = *npt;
2293             for (k = 1; k <= i__1; ++k) {
2294 /* L700: */
2295                 vlag[k] = fval[k] - fval[kopt];
2296             }
2297             gisq = zero;
2298             i__1 = *n;
2299             for (i__ = 1; i__ <= i__1; ++i__) {
2300                 sum = zero;
2301                 i__2 = *npt;
2302                 for (k = 1; k <= i__2; ++k) {
2303 /* L710: */
2304                     sum += bmat[k + i__ * bmat_dim1] * vlag[k];
2305                 }
2306                 gisq += sum * sum;
2307 /* L720: */
2308                 w[i__] = sum;
2309             }
2310
2311 /* Test whether to replace the new quadratic model by the least Frobenius */
2312 /* norm interpolant, making the replacement if the test is satisfied. */
2313
2314             ++itest;
2315             if (gqsq < gisq * 100.) {
2316                 itest = 0;
2317             }
2318             if (itest >= 3) {
2319                 i__1 = *n;
2320                 for (i__ = 1; i__ <= i__1; ++i__) {
2321 /* L730: */
2322                     gq[i__] = w[i__];
2323                 }
2324                 i__1 = nh;
2325                 for (ih = 1; ih <= i__1; ++ih) {
2326 /* L740: */
2327                     hq[ih] = zero;
2328                 }
2329                 i__1 = nptm;
2330                 for (j = 1; j <= i__1; ++j) {
2331                     w[j] = zero;
2332                     i__2 = *npt;
2333                     for (k = 1; k <= i__2; ++k) {
2334 /* L750: */
2335                         w[j] += vlag[k] * zmat[k + j * zmat_dim1];
2336                     }
2337 /* L760: */
2338                     if (j < idz) {
2339                         w[j] = -w[j];
2340                     }
2341                 }
2342                 i__1 = *npt;
2343                 for (k = 1; k <= i__1; ++k) {
2344                     pq[k] = zero;
2345                     i__2 = nptm;
2346                     for (j = 1; j <= i__2; ++j) {
2347 /* L770: */
2348                         pq[k] += zmat[k + j * zmat_dim1] * w[j];
2349                     }
2350                 }
2351                 itest = 0;
2352             }
2353         }
2354     }
2355     if (f < fsave) {
2356         kopt = knew;
2357     }
2358
2359 /* If a trust region step has provided a sufficient decrease in F, then */
2360 /* branch for another trust region calculation. The case KSAVE>0 occurs */
2361 /* when the new function value was calculated by a model step. */
2362
2363     if (f <= fsave + tenth * vquad) {
2364         goto L100;
2365     }
2366     if (ksave > 0) {
2367         goto L100;
2368     }
2369
2370 /* Alternatively, find out if the interpolation points are close enough */
2371 /* to the best point so far. */
2372
2373     knew = 0;
2374 L460:
2375     distsq = delta * 4. * delta;
2376     i__2 = *npt;
2377     for (k = 1; k <= i__2; ++k) {
2378         sum = zero;
2379         i__1 = *n;
2380         for (j = 1; j <= i__1; ++j) {
2381 /* L470: */
2382 /* Computing 2nd power */
2383             d__1 = xpt[k + j * xpt_dim1] - xopt[j];
2384             sum += d__1 * d__1;
2385         }
2386         if (sum > distsq) {
2387             knew = k;
2388             distsq = sum;
2389         }
2390 /* L480: */
2391     }
2392
2393 /* If KNEW is positive, then set DSTEP, and branch back for the next */
2394 /* iteration, which will generate a "model step". */
2395
2396     if (knew > 0) {
2397 /* Computing MAX */
2398 /* Computing MIN */
2399         d__2 = tenth * sqrt(distsq), d__3 = half * delta;
2400         d__1 = min(d__2,d__3);
2401         dstep = max(d__1,rho);
2402         dsq = dstep * dstep;
2403         goto L120;
2404     }
2405     if (ratio > zero) {
2406         goto L100;
2407     }
2408     if (max(delta,dnorm) > rho) {
2409         goto L100;
2410     }
2411
2412 /* The calculations with the current value of RHO are complete. Pick the */
2413 /* next values of RHO and DELTA. */
2414
2415 L490:
2416     if (rho > rhoend) {
2417         delta = half * rho;
2418         ratio = rho / rhoend;
2419         if (ratio <= 16.) {
2420             rho = rhoend;
2421         } else if (ratio <= 250.) {
2422             rho = sqrt(ratio) * rhoend;
2423         } else {
2424             rho = tenth * rho;
2425         }
2426         delta = max(delta,rho);
2427         goto L90;
2428     }
2429
2430 /* Return from the calculation, after another Newton-Raphson step, if */
2431 /* it is too short to have been tried before. */
2432
2433     if (knew == -1) {
2434         goto L290;
2435     }
2436     rc = NLOPT_XTOL_REACHED;
2437 L530:
2438     if (fopt <= f) {
2439         i__2 = *n;
2440         for (i__ = 1; i__ <= i__2; ++i__) {
2441 /* L540: */
2442             x[i__] = xbase[i__] + xopt[i__];
2443         }
2444         f = fopt;
2445     }
2446     *minf = f;
2447     return rc;
2448 } /* newuob_ */
2449
2450 /*************************************************************************/
2451 /* newuoa.f */
2452
2453 nlopt_result newuoa(int n, int npt, double *x,
2454                     const double *lb, const double *ub,
2455                     double rhobeg, nlopt_stopping *stop, double *minf,
2456                     newuoa_func calfun, void *calfun_data)
2457 {
2458     /* Local variables */
2459     int id, np, iw, igq, ihq, ixb, ifv, ipq, ivl, ixn, ixo, ixp, ndim,
2460             nptm, ibmat, izmat;
2461     nlopt_result ret;
2462     double *w;
2463
2464 /* This subroutine seeks the least value of a function of many variables, */
2465 /* by a trust region method that forms quadratic models by interpolation. */
2466 /* There can be some freedom in the interpolation conditions, which is */
2467 /* taken up by minimizing the Frobenius norm of the change to the second */
2468 /* derivative of the quadratic model, beginning with a zero matrix. The */
2469 /* arguments of the subroutine are as follows. */
2470
2471 /* N must be set to the number of variables and must be at least two. */
2472 /* NPT is the number of interpolation conditions. Its value must be in the */
2473 /*   interval [N+2,(N+1)(N+2)/2]. */
2474 /* Initial values of the variables must be set in X(1),X(2),...,X(N). They */
2475 /*   will be changed to the values that give the least calculated F. */
2476 /* RHOBEG and RHOEND must be set to the initial and final values of a trust */
2477 /*   region radius, so both must be positive with RHOEND<=RHOBEG. Typically */
2478 /*   RHOBEG should be about one tenth of the greatest expected change to a */
2479 /*   variable, and RHOEND should indicate the accuracy that is required in */
2480 /*   the final values of the variables. */
2481 /* The value of IPRINT should be set to 0, 1, 2 or 3, which controls the */
2482 /*   amount of printing. Specifically, there is no output if IPRINT=0 and */
2483 /*   there is output only at the return if IPRINT=1. Otherwise, each new */
2484 /*   value of RHO is printed, with the best vector of variables so far and */
2485 /*   the corresponding value of the objective function. Further, each new */
2486 /*   value of F with its variables are output if IPRINT=3. */
2487 /* MAXFUN must be set to an upper bound on the number of calls of CALFUN. */
2488 /* The array W will be used for working space. Its length must be at least */
2489 /* (NPT+13)*(NPT+N)+3*N*(N+3)/2. */
2490
2491 /* SUBROUTINE CALFUN (N,X,F) must be provided by the user. It must set F to */
2492 /* the value of the objective function for the variables X(1),X(2),...,X(N). */
2493
2494 /* Partition the working space array, so that different parts of it can be */
2495 /* treated separately by the subroutine that performs the main calculation. */
2496
2497     /* Parameter adjustments */
2498     --x;
2499
2500     /* Function Body */
2501     np = n + 1;
2502     nptm = npt - np;
2503     if (n < 2 || npt < n + 2 || npt > (n + 2) * np / 2) {
2504          return NLOPT_INVALID_ARGS;
2505     }
2506     ndim = npt + n;
2507     ixb = 1;
2508     ixo = ixb + n;
2509     ixn = ixo + n;
2510     ixp = ixn + n;
2511     ifv = ixp + n * npt;
2512     igq = ifv + npt;
2513     ihq = igq + n;
2514     ipq = ihq + n * np / 2;
2515     ibmat = ipq + npt;
2516     izmat = ibmat + ndim * n;
2517     id = izmat + npt * nptm;
2518     ivl = id + n;
2519     iw = ivl + ndim;
2520
2521     w = (double *) malloc(sizeof(double) * ((npt+13)*(npt+n) + 3*(n*(n+3))/2));
2522     if (!w) return NLOPT_OUT_OF_MEMORY;
2523     --w;
2524
2525 /* The above settings provide a partition of W for subroutine NEWUOB. */
2526 /* The partition requires the first NPT*(NPT+N)+5*N*(N+3)/2 elements of */
2527 /* W plus the space that is needed by the last array of NEWUOB. */
2528
2529     ret = newuob_(&n, &npt, &x[1], &rhobeg,
2530                   lb, ub, stop, minf, calfun, calfun_data,
2531                   &w[ixb], &w[ixo], &w[ixn], &w[ixp], &w[ifv],
2532                   &w[igq], &w[ihq], &w[ipq], &w[ibmat], &w[izmat],
2533                   &ndim, &w[id], &w[ivl], &w[iw]);
2534
2535     ++w;
2536     free(w);
2537     return ret;
2538 } /* newuoa_ */
2539
2540 /*************************************************************************/
2541 /*************************************************************************/