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