cobyla/cobyla.c

   1 /* cobyla : contrained optimization by linear approximation */
   2
   3 /*
   4  * Copyright (c) 1992, Michael J. D. Powell (M.J.D.Powell@damtp.cam.ac.uk)
   5  * Copyright (c) 2004, Jean-Sebastien Roy (js@jeannot.org)
   6  *
   7  * Permission is hereby granted, free of charge, to any person obtaining a
   8  * copy of this software and associated documentation files (the
   9  * "Software"), to deal in the Software without restriction, including
  10  * without limitation the rights to use, copy, modify, merge, publish,
  11  * distribute, sublicense, and/or sell copies of the Software, and to
  12  * permit persons to whom the Software is furnished to do so, subject to
  13  * the following conditions:
  14  *
  15  * The above copyright notice and this permission notice shall be included
  16  * in all copies or substantial portions of the Software.
  17  *
  18  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  19  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  20  * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
  21  * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
  22  * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
  23  * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
  24  * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  25  */
  26
  27 /*
  28  * This software is a C version of COBYLA2, a contrained optimization by linear
  29  * approximation package developed by Michael J. D. Powell in Fortran.
  30  *
  31  * The original source code can be found at :
  32  * http://plato.la.asu.edu/topics/problems/nlores.html
  33  */
  34
  35 static char const rcsid[] =
  36   "@(#) $Jeannot: cobyla.c,v 1.11 2004/04/18 09:51:36 js Exp $";
  37
  38 #include <stdlib.h>
  39 #include <stdio.h>
  40 #include <math.h>
  41
  42 #include "cobyla.h"
  43
  44 #define min(a,b) ((a) <= (b) ? (a) : (b))
  45 #define max(a,b) ((a) >= (b) ? (a) : (b))
  46 #define abs(x) ((x) >= 0 ? (x) : -(x))
  47
  48 /**************************************************************************/
  49 /* SGJ, 2008: NLopt-style interface function: */
  50
  51 typedef struct {
  52      nlopt_func f;
  53      void *f_data;
  54      int m_orig;
  55      nlopt_func fc;
  56      char *fc_data;
  57      ptrdiff_t fc_datum_size;
  58      double *xtmp;
  59      const double *lb, *ub;
  60 } func_wrap_state;
  61
  62 static int func_wrap(int n, int m, double *x, double *f, double *con,
  63                      void *state)
  64 {
  65      int i, j;
  66      func_wrap_state *s = (func_wrap_state *) state;
  67      double *xtmp = s->xtmp;
  68      const double *lb = s->lb, *ub = s->ub;
  69
  70      /* in nlopt, we guarante that the function is never evaluated outside
  71         the lb and ub bounds, so we need force this with xtmp */
  72      for (j = 0; j < n; ++j) {
  73           if (x[j] < lb[j]) xtmp[j] = lb[j];
  74           else if (x[j] > ub[j]) xtmp[j] = ub[j];
  75           else xtmp[j] = x[j];
  76      }
  77
  78      *f = s->f(n, xtmp, NULL, s->f_data);
  79      for (i = 0; i < s->m_orig; ++i)
  80           con[i] = -s->fc(n, xtmp, NULL, s->fc_data + i * s->fc_datum_size);
  81      for (j = 0; j < n; ++j) {
  82           if (!nlopt_isinf(lb[j]))
  83                con[i++] = x[j] - lb[j];
  84           if (!nlopt_isinf(ub[j]))
  85                con[i++] = ub[j] - x[j];
  86      }
  87      if (m != i) return 1; /* ... bug?? */
  88      return 0;
  89 }
  90
  91 nlopt_result cobyla_minimize(int n, nlopt_func f, void *f_data,
  92                              int m, nlopt_func fc,
  93                              void *fc_data_, ptrdiff_t fc_datum_size,
  94                              const double *lb, const double *ub, /* bounds */
  95                              double *x, /* in: initial guess, out: minimizer */
  96                              double *minf,
  97                              nlopt_stopping *stop,
  98                              double rhobegin)
  99 {
 100      int j;
 101      func_wrap_state s;
 102      nlopt_result ret;
 103
 104      s.f = f; s.f_data = f_data;
 105      s.m_orig = m;
 106      s.fc = fc; s.fc_data = (char*) fc_data_; s.fc_datum_size = fc_datum_size;
 107      s.lb = lb; s.ub = ub;
 108      s.xtmp = (double *) malloc(sizeof(double) * n);
 109      if (!s.xtmp) return NLOPT_OUT_OF_MEMORY;
 110
 111      /* add constraints for lower/upper bounds (if any) */
 112      for (j = 0; j < n; ++j) {
 113           if (!nlopt_isinf(lb[j]))
 114                ++m;
 115           if (!nlopt_isinf(ub[j]))
 116                ++m;
 117      }
 118
 119      ret = cobyla(n, m, x, minf, rhobegin, stop, COBYLA_MSG_NONE,
 120                   func_wrap, &s);
 121
 122      free(s.xtmp);
 123      return ret;
 124 }
 125
 126 /**************************************************************************/
 127
 128 static nlopt_result cobylb(int *n, int *m, int *mpp, double *x, double *minf, double *rhobeg,
 129   nlopt_stopping *stop, int *iprint, double *con, double *sim,
 130   double *simi, double *datmat, double *a, double *vsig, double *veta,
 131   double *sigbar, double *dx, double *w, int *iact, cobyla_function *calcfc,
 132   void *state);
 133 static void trstlp(int *n, int *m, double *a, double *b, double *rho,
 134   double *dx, int *ifull, int *iact, double *z__, double *zdota, double *vmultc,
 135   double *sdirn, double *dxnew, double *vmultd);
 136
 137 /* ------------------------------------------------------------------------ */
 138
 139 nlopt_result cobyla(int n, int m, double *x, double *minf, double rhobeg, nlopt_stopping *stop, int iprint,
 140   cobyla_function *calcfc, void *state)
 141 {
 142   int icon, isim, isigb, idatm, iveta, isimi, ivsig, iwork, ia, idx, mpp;
 143   int *iact;
 144   double *w;
 145   nlopt_result rc;
 146
 147 /*
 148  * This subroutine minimizes an objective function F(X) subject to M
 149  * inequality constraints on X, where X is a vector of variables that has
 150  * N components. The algorithm employs linear approximations to the
 151  * objective and constraint functions, the approximations being formed by
 152  * linear interpolation at N+1 points in the space of the variables.
 153  * We regard these interpolation points as vertices of a simplex. The
 154  * parameter RHO controls the size of the simplex and it is reduced
 155  * automatically from RHOBEG to RHOEND. For each RHO the subroutine tries
 156  * to achieve a good vector of variables for the current size, and then
 157  * RHO is reduced until the value RHOEND is reached. Therefore RHOBEG and
 158  * RHOEND should be set to reasonable initial changes to and the required
 159  * accuracy in the variables respectively, but this accuracy should be
 160  * viewed as a subject for experimentation because it is not guaranteed.
 161  * The subroutine has an advantage over many of its competitors, however,
 162  * which is that it treats each constraint individually when calculating
 163  * a change to the variables, instead of lumping the constraints together
 164  * into a single penalty function. The name of the subroutine is derived
 165  * from the phrase Constrained Optimization BY Linear Approximations.
 166  *
 167  * The user must set the values of N, M, RHOBEG and RHOEND, and must
 168  * provide an initial vector of variables in X. Further, the value of
 169  * IPRINT should be set to 0, 1, 2 or 3, which controls the amount of
 170  * printing during the calculation. Specifically, there is no output if
 171  * IPRINT=0 and there is output only at the end of the calculation if
 172  * IPRINT=1. Otherwise each new value of RHO and SIGMA is printed.
 173  * Further, the vector of variables and some function information are
 174  * given either when RHO is reduced or when each new value of F(X) is
 175  * computed in the cases IPRINT=2 or IPRINT=3 respectively. Here SIGMA
 176  * is a penalty parameter, it being assumed that a change to X is an
 177  * improvement if it reduces the merit function
 178  *      F(X)+SIGMA*MAX(0.0,-C1(X),-C2(X),...,-CM(X)),
 179  * where C1,C2,...,CM denote the constraint functions that should become
 180  * nonnegative eventually, at least to the precision of RHOEND. In the
 181  * printed output the displayed term that is multiplied by SIGMA is
 182  * called MAXCV, which stands for 'MAXimum Constraint Violation'. The
 183  * argument MAXFUN is an int variable that must be set by the user to a
 184  * limit on the number of calls of CALCFC, the purpose of this routine being
 185  * given below. The value of MAXFUN will be altered to the number of calls
 186  * of CALCFC that are made. The arguments W and IACT provide real and
 187  * int arrays that are used as working space. Their lengths must be at
 188  * least N*(3*N+2*M+11)+4*M+6 and M+1 respectively.
 189  *
 190  * In order to define the objective and constraint functions, we require
 191  * a subroutine that has the name and arguments
 192  *      SUBROUTINE CALCFC (N,M,X,F,CON)
 193  *      DIMENSION X(*),CON(*)  .
 194  * The values of N and M are fixed and have been defined already, while
 195  * X is now the current vector of variables. The subroutine should return
 196  * the objective and constraint functions at X in F and CON(1),CON(2),
 197  * ...,CON(M). Note that we are trying to adjust X so that F(X) is as
 198  * small as possible subject to the constraint functions being nonnegative.
 199  *
 200  * Partition the working space array W to provide the storage that is needed
 201  * for the main calculation.
 202  */
 203
 204   stop->nevals = 0;
 205
 206   if (n == 0)
 207   {
 208     if (iprint>=1) fprintf(stderr, "cobyla: N==0.\n");
 209     return NLOPT_SUCCESS;
 210   }
 211
 212   if (n < 0 || m < 0)
 213   {
 214     if (iprint>=1) fprintf(stderr, "cobyla: N<0 or M<0.\n");
 215     return NLOPT_INVALID_ARGS;
 216   }
 217
 218   /* workspace allocation */
 219   w = malloc((n*(3*n+2*m+11)+4*m+6)*sizeof(*w));
 220   if (w == NULL)
 221   {
 222     if (iprint>=1) fprintf(stderr, "cobyla: memory allocation error.\n");
 223     return NLOPT_OUT_OF_MEMORY;
 224   }
 225   iact = malloc((m+1)*sizeof(*iact));
 226   if (iact == NULL)
 227   {
 228     if (iprint>=1) fprintf(stderr, "cobyla: memory allocation error.\n");
 229     free(w);
 230     return NLOPT_OUT_OF_MEMORY;
 231   }
 232
 233   /* Parameter adjustments */
 234   --iact;
 235   --w;
 236   --x;
 237
 238   /* Function Body */
 239   mpp = m + 2;
 240   icon = 1;
 241   isim = icon + mpp;
 242   isimi = isim + n * n + n;
 243   idatm = isimi + n * n;
 244   ia = idatm + n * mpp + mpp;
 245   ivsig = ia + m * n + n;
 246   iveta = ivsig + n;
 247   isigb = iveta + n;
 248   idx = isigb + n;
 249   iwork = idx + n;
 250   rc = cobylb(&n, &m, &mpp, &x[1], minf, &rhobeg, stop, &iprint,
 251       &w[icon], &w[isim], &w[isimi], &w[idatm], &w[ia], &w[ivsig], &w[iveta],
 252       &w[isigb], &w[idx], &w[iwork], &iact[1], calcfc, state);
 253
 254   /* Parameter adjustments (reverse) */
 255   ++iact;
 256   ++w;
 257
 258   free(w);
 259   free(iact);
 260
 261   return rc;
 262 } /* cobyla */
 263
 264 /* ------------------------------------------------------------------------- */
 265 static nlopt_result cobylb(int *n, int *m, int *mpp, double
 266     *x, double *minf, double *rhobeg, nlopt_stopping *stop, int *iprint,
 267     double *con, double *sim, double *simi,
 268     double *datmat, double *a, double *vsig, double *veta,
 269      double *sigbar, double *dx, double *w, int *iact, cobyla_function *calcfc,
 270      void *state)
 271 {
 272   /* System generated locals */
 273   int sim_dim1, sim_offset, simi_dim1, simi_offset, datmat_dim1,
 274       datmat_offset, a_dim1, a_offset, i__1, i__2, i__3;
 275   double d__1, d__2;
 276
 277   /* Local variables */
 278   double alpha, delta, denom, tempa, barmu;
 279   double beta, cmin = 0.0, cmax = 0.0;
 280   double cvmaxm, dxsign, prerem = 0.0;
 281   double edgmax, pareta, prerec = 0.0, phimin, parsig = 0.0;
 282   double gamma_;
 283   double phi, rho, sum = 0.0;
 284   double ratio, vmold, parmu, error, vmnew;
 285   double resmax, cvmaxp;
 286   double resnew, trured;
 287   double temp, wsig, f;
 288   double weta;
 289   int i__, j, k, l;
 290   int idxnew;
 291   int iflag = 0;
 292   int iptemp;
 293   int isdirn, izdota;
 294   int ivmc;
 295   int ivmd;
 296   int mp, np, iz, ibrnch;
 297   int nbest, ifull, iptem, jdrop;
 298   nlopt_result rc = NLOPT_SUCCESS;
 299   double rhoend;
 300
 301   /* SGJ, 2008: compute rhoend from NLopt stop info */
 302   rhoend = stop->xtol_rel * (*rhobeg);
 303   for (j = 0; j < *n; ++j)
 304        if (rhoend < stop->xtol_abs[j])
 305             rhoend = stop->xtol_abs[j];
 306
 307   /* SGJ, 2008: added code to keep track of minimum feasible function val */
 308   *minf = HUGE_VAL;
 309
 310 /* Set the initial values of some parameters. The last column of SIM holds */
 311 /* the optimal vertex of the current simplex, and the preceding N columns */
 312 /* hold the displacements from the optimal vertex to the other vertices. */
 313 /* Further, SIMI holds the inverse of the matrix that is contained in the */
 314 /* first N columns of SIM. */
 315
 316   /* Parameter adjustments */
 317   a_dim1 = *n;
 318   a_offset = 1 + a_dim1 * 1;
 319   a -= a_offset;
 320   simi_dim1 = *n;
 321   simi_offset = 1 + simi_dim1 * 1;
 322   simi -= simi_offset;
 323   sim_dim1 = *n;
 324   sim_offset = 1 + sim_dim1 * 1;
 325   sim -= sim_offset;
 326   datmat_dim1 = *mpp;
 327   datmat_offset = 1 + datmat_dim1 * 1;
 328   datmat -= datmat_offset;
 329   --x;
 330   --con;
 331   --vsig;
 332   --veta;
 333   --sigbar;
 334   --dx;
 335   --w;
 336   --iact;
 337
 338   /* Function Body */
 339   iptem = min(*n,4);
 340   iptemp = iptem + 1;
 341   np = *n + 1;
 342   mp = *m + 1;
 343   alpha = .25;
 344   beta = 2.1;
 345   gamma_ = .5;
 346   delta = 1.1;
 347   rho = *rhobeg;
 348   parmu = 0.;
 349   if (*iprint >= 2) {
 350     fprintf(stderr,
 351       "cobyla: the initial value of RHO is %12.6E and PARMU is set to zero.\n",
 352       rho);
 353   }
 354   temp = 1. / rho;
 355   i__1 = *n;
 356   for (i__ = 1; i__ <= i__1; ++i__) {
 357     sim[i__ + np * sim_dim1] = x[i__];
 358     i__2 = *n;
 359     for (j = 1; j <= i__2; ++j) {
 360       sim[i__ + j * sim_dim1] = 0.;
 361       simi[i__ + j * simi_dim1] = 0.;
 362     }
 363     sim[i__ + i__ * sim_dim1] = rho;
 364     simi[i__ + i__ * simi_dim1] = temp;
 365   }
 366   jdrop = np;
 367   ibrnch = 0;
 368
 369 /* Make the next call of the user-supplied subroutine CALCFC. These */
 370 /* instructions are also used for calling CALCFC during the iterations of */
 371 /* the algorithm. */
 372
 373 L40:
 374   if (nlopt_stop_evals(stop)) rc = NLOPT_MAXEVAL_REACHED;
 375   else if (nlopt_stop_time(stop)) rc = NLOPT_MAXTIME_REACHED;
 376   if (rc != NLOPT_SUCCESS) goto L600;
 377
 378   stop->nevals++;
 379   if (calcfc(*n, *m, &x[1], &f, &con[1], state))
 380   {
 381     if (*iprint >= 1) {
 382       fprintf(stderr, "cobyla: user requested end of minimization.\n");
 383     }
 384     rc = NLOPT_FAILURE;
 385     goto L600;
 386   }
 387
 388   resmax = 0.;
 389   if (*m > 0) {
 390     i__1 = *m;
 391     for (k = 1; k <= i__1; ++k) {
 392       d__1 = resmax, d__2 = -con[k];
 393       resmax = max(d__1,d__2);
 394     }
 395   }
 396
 397   /* SGJ, 2008: check for minf_max reached by feasible point */
 398   if (f < stop->minf_max && resmax <= 0) {
 399        rc = NLOPT_MINF_MAX_REACHED;
 400        goto L620; /* not L600 because we want to use current x, f, resmax */
 401   }
 402
 403   if (stop->nevals == *iprint - 1 || *iprint == 3) {
 404     fprintf(stderr, "cobyla: NFVALS = %4d, F =%13.6E, MAXCV =%13.6E\n",
 405             stop->nevals, f, resmax);
 406     i__1 = iptem;
 407     fprintf(stderr, "cobyla: X =");
 408     for (i__ = 1; i__ <= i__1; ++i__) {
 409       if (i__>1) fprintf(stderr, "  ");
 410       fprintf(stderr, "%13.6E", x[i__]);
 411     }
 412     if (iptem < *n) {
 413       i__1 = *n;
 414       for (i__ = iptemp; i__ <= i__1; ++i__) {
 415         if (!((i__-1) % 4)) fprintf(stderr, "\ncobyla:  ");
 416         fprintf(stderr, "%15.6E", x[i__]);
 417       }
 418     }
 419     fprintf(stderr, "\n");
 420   }
 421   con[mp] = f;
 422   con[*mpp] = resmax;
 423   if (ibrnch == 1) {
 424     goto L440;
 425   }
 426
 427 /* Set the recently calculated function values in a column of DATMAT. This */
 428 /* array has a column for each vertex of the current simplex, the entries of */
 429 /* each column being the values of the constraint functions (if any) */
 430 /* followed by the objective function and the greatest constraint violation */
 431 /* at the vertex. */
 432
 433   i__1 = *mpp;
 434   for (k = 1; k <= i__1; ++k) {
 435     datmat[k + jdrop * datmat_dim1] = con[k];
 436   }
 437   if (stop->nevals > np) {
 438     goto L130;
 439   }
 440
 441 /* Exchange the new vertex of the initial simplex with the optimal vertex if */
 442 /* necessary. Then, if the initial simplex is not complete, pick its next */
 443 /* vertex and calculate the function values there. */
 444
 445   if (jdrop <= *n) {
 446     if (datmat[mp + np * datmat_dim1] <= f) {
 447       x[jdrop] = sim[jdrop + np * sim_dim1];
 448     } else { /* improvement in function val */
 449       sim[jdrop + np * sim_dim1] = x[jdrop];
 450       i__1 = *mpp;
 451       for (k = 1; k <= i__1; ++k) {
 452         datmat[k + jdrop * datmat_dim1] = datmat[k + np * datmat_dim1]
 453             ;
 454         datmat[k + np * datmat_dim1] = con[k];
 455       }
 456       i__1 = jdrop;
 457       for (k = 1; k <= i__1; ++k) {
 458         sim[jdrop + k * sim_dim1] = -rho;
 459         temp = 0.f;
 460         i__2 = jdrop;
 461         for (i__ = k; i__ <= i__2; ++i__) {
 462           temp -= simi[i__ + k * simi_dim1];
 463         }
 464         simi[jdrop + k * simi_dim1] = temp;
 465       }
 466     }
 467   }
 468   if (stop->nevals <= *n) {
 469     jdrop = stop->nevals;
 470     x[jdrop] += rho;
 471     goto L40;
 472   }
 473 L130:
 474   ibrnch = 1;
 475
 476 /* Identify the optimal vertex of the current simplex. */
 477
 478 L140:
 479   phimin = datmat[mp + np * datmat_dim1] + parmu * datmat[*mpp + np *
 480       datmat_dim1];
 481   nbest = np;
 482   i__1 = *n;
 483   for (j = 1; j <= i__1; ++j) {
 484     temp = datmat[mp + j * datmat_dim1] + parmu * datmat[*mpp + j *
 485         datmat_dim1];
 486     if (temp < phimin) {
 487       nbest = j;
 488       phimin = temp;
 489     } else if (temp == phimin && parmu == 0.) {
 490       if (datmat[*mpp + j * datmat_dim1] < datmat[*mpp + nbest *
 491           datmat_dim1]) {
 492         nbest = j;
 493       }
 494     }
 495   }
 496
 497 /* Switch the best vertex into pole position if it is not there already, */
 498 /* and also update SIM, SIMI and DATMAT. */
 499
 500   if (nbest <= *n) {
 501     i__1 = *mpp;
 502     for (i__ = 1; i__ <= i__1; ++i__) {
 503       temp = datmat[i__ + np * datmat_dim1];
 504       datmat[i__ + np * datmat_dim1] = datmat[i__ + nbest * datmat_dim1]
 505           ;
 506       datmat[i__ + nbest * datmat_dim1] = temp;
 507     }
 508     i__1 = *n;
 509     for (i__ = 1; i__ <= i__1; ++i__) {
 510       temp = sim[i__ + nbest * sim_dim1];
 511       sim[i__ + nbest * sim_dim1] = 0.;
 512       sim[i__ + np * sim_dim1] += temp;
 513       tempa = 0.;
 514       i__2 = *n;
 515       for (k = 1; k <= i__2; ++k) {
 516         sim[i__ + k * sim_dim1] -= temp;
 517         tempa -= simi[k + i__ * simi_dim1];
 518       }
 519       simi[nbest + i__ * simi_dim1] = tempa;
 520     }
 521   }
 522
 523 /* Make an error return if SIGI is a poor approximation to the inverse of */
 524 /* the leading N by N submatrix of SIG. */
 525
 526   error = 0.;
 527   i__1 = *n;
 528   for (i__ = 1; i__ <= i__1; ++i__) {
 529     i__2 = *n;
 530     for (j = 1; j <= i__2; ++j) {
 531       temp = 0.;
 532       if (i__ == j) {
 533         temp += -1.;
 534       }
 535       i__3 = *n;
 536       for (k = 1; k <= i__3; ++k) if (sim[k + j * sim_dim1] != 0) {
 537         temp += simi[i__ + k * simi_dim1] * sim[k + j * sim_dim1];
 538       }
 539       d__1 = error, d__2 = abs(temp);
 540       error = max(d__1,d__2);
 541     }
 542   }
 543   if (error > .1) {
 544     if (*iprint >= 1) {
 545       fprintf(stderr, "cobyla: rounding errors are becoming damaging.\n");
 546     }
 547     rc = NLOPT_FAILURE;
 548     goto L600;
 549   }
 550
 551 /* Calculate the coefficients of the linear approximations to the objective */
 552 /* and constraint functions, placing minus the objective function gradient */
 553 /* after the constraint gradients in the array A. The vector W is used for */
 554 /* working space. */
 555
 556   i__2 = mp;
 557   for (k = 1; k <= i__2; ++k) {
 558     con[k] = -datmat[k + np * datmat_dim1];
 559     i__1 = *n;
 560     for (j = 1; j <= i__1; ++j) {
 561       w[j] = datmat[k + j * datmat_dim1] + con[k];
 562     }
 563     i__1 = *n;
 564     for (i__ = 1; i__ <= i__1; ++i__) {
 565       temp = 0.;
 566       i__3 = *n;
 567       for (j = 1; j <= i__3; ++j) {
 568         temp += w[j] * simi[j + i__ * simi_dim1];
 569       }
 570       if (k == mp) {
 571         temp = -temp;
 572       }
 573       a[i__ + k * a_dim1] = temp;
 574     }
 575   }
 576
 577 /* Calculate the values of sigma and eta, and set IFLAG=0 if the current */
 578 /* simplex is not acceptable. */
 579
 580   iflag = 1;
 581   parsig = alpha * rho;
 582   pareta = beta * rho;
 583   i__1 = *n;
 584   for (j = 1; j <= i__1; ++j) {
 585     wsig = 0.;
 586     weta = 0.;
 587     i__2 = *n;
 588     for (i__ = 1; i__ <= i__2; ++i__) {
 589       d__1 = simi[j + i__ * simi_dim1];
 590       wsig += d__1 * d__1;
 591       d__1 = sim[i__ + j * sim_dim1];
 592       weta += d__1 * d__1;
 593     }
 594     vsig[j] = 1. / sqrt(wsig);
 595     veta[j] = sqrt(weta);
 596     if (vsig[j] < parsig || veta[j] > pareta) {
 597       iflag = 0;
 598     }
 599   }
 600
 601 /* If a new vertex is needed to improve acceptability, then decide which */
 602 /* vertex to drop from the simplex. */
 603
 604   if (ibrnch == 1 || iflag == 1) {
 605     goto L370;
 606   }
 607   jdrop = 0;
 608   temp = pareta;
 609   i__1 = *n;
 610   for (j = 1; j <= i__1; ++j) {
 611     if (veta[j] > temp) {
 612       jdrop = j;
 613       temp = veta[j];
 614     }
 615   }
 616   if (jdrop == 0) {
 617     i__1 = *n;
 618     for (j = 1; j <= i__1; ++j) {
 619       if (vsig[j] < temp) {
 620         jdrop = j;
 621         temp = vsig[j];
 622       }
 623     }
 624   }
 625
 626 /* Calculate the step to the new vertex and its sign. */
 627
 628   temp = gamma_ * rho * vsig[jdrop];
 629   i__1 = *n;
 630   for (i__ = 1; i__ <= i__1; ++i__) {
 631     dx[i__] = temp * simi[jdrop + i__ * simi_dim1];
 632   }
 633   cvmaxp = 0.;
 634   cvmaxm = 0.;
 635   i__1 = mp;
 636   for (k = 1; k <= i__1; ++k) {
 637     sum = 0.;
 638     i__2 = *n;
 639     for (i__ = 1; i__ <= i__2; ++i__) {
 640       sum += a[i__ + k * a_dim1] * dx[i__];
 641     }
 642     if (k < mp) {
 643       temp = datmat[k + np * datmat_dim1];
 644       d__1 = cvmaxp, d__2 = -sum - temp;
 645       cvmaxp = max(d__1,d__2);
 646       d__1 = cvmaxm, d__2 = sum - temp;
 647       cvmaxm = max(d__1,d__2);
 648     }
 649   }
 650   dxsign = 1.;
 651   if (parmu * (cvmaxp - cvmaxm) > sum + sum) {
 652     dxsign = -1.;
 653   }
 654
 655 /* Update the elements of SIM and SIMI, and set the next X. */
 656
 657   temp = 0.;
 658   i__1 = *n;
 659   for (i__ = 1; i__ <= i__1; ++i__) {
 660     dx[i__] = dxsign * dx[i__];
 661     sim[i__ + jdrop * sim_dim1] = dx[i__];
 662     temp += simi[jdrop + i__ * simi_dim1] * dx[i__];
 663   }
 664   i__1 = *n;
 665   for (i__ = 1; i__ <= i__1; ++i__) {
 666     simi[jdrop + i__ * simi_dim1] /= temp;
 667   }
 668   i__1 = *n;
 669   for (j = 1; j <= i__1; ++j) {
 670     if (j != jdrop) {
 671       temp = 0.;
 672       i__2 = *n;
 673       for (i__ = 1; i__ <= i__2; ++i__) {
 674         temp += simi[j + i__ * simi_dim1] * dx[i__];
 675       }
 676       i__2 = *n;
 677       for (i__ = 1; i__ <= i__2; ++i__) {
 678         simi[j + i__ * simi_dim1] -= temp * simi[jdrop + i__ *
 679             simi_dim1];
 680       }
 681     }
 682     x[j] = sim[j + np * sim_dim1] + dx[j];
 683   }
 684   goto L40;
 685
 686 /* Calculate DX=x(*)-x(0). Branch if the length of DX is less than 0.5*RHO. */
 687
 688 L370:
 689   iz = 1;
 690   izdota = iz + *n * *n;
 691   ivmc = izdota + *n;
 692   isdirn = ivmc + mp;
 693   idxnew = isdirn + *n;
 694   ivmd = idxnew + *n;
 695   trstlp(n, m, &a[a_offset], &con[1], &rho, &dx[1], &ifull, &iact[1], &w[
 696       iz], &w[izdota], &w[ivmc], &w[isdirn], &w[idxnew], &w[ivmd]);
 697   if (ifull == 0) {
 698     temp = 0.;
 699     i__1 = *n;
 700     for (i__ = 1; i__ <= i__1; ++i__) {
 701       d__1 = dx[i__];
 702       temp += d__1 * d__1;
 703     }
 704     if (temp < rho * .25 * rho) {
 705       ibrnch = 1;
 706       goto L550;
 707     }
 708   }
 709
 710 /* Predict the change to F and the new maximum constraint violation if the */
 711 /* variables are altered from x(0) to x(0)+DX. */
 712
 713   resnew = 0.;
 714   con[mp] = 0.;
 715   i__1 = mp;
 716   for (k = 1; k <= i__1; ++k) {
 717     sum = con[k];
 718     i__2 = *n;
 719     for (i__ = 1; i__ <= i__2; ++i__) {
 720       sum -= a[i__ + k * a_dim1] * dx[i__];
 721     }
 722     if (k < mp) {
 723       resnew = max(resnew,sum);
 724     }
 725   }
 726
 727 /* Increase PARMU if necessary and branch back if this change alters the */
 728 /* optimal vertex. Otherwise PREREM and PREREC will be set to the predicted */
 729 /* reductions in the merit function and the maximum constraint violation */
 730 /* respectively. */
 731
 732   barmu = 0.;
 733   prerec = datmat[*mpp + np * datmat_dim1] - resnew;
 734   if (prerec > 0.) {
 735     barmu = sum / prerec;
 736   }
 737   if (parmu < barmu * 1.5) {
 738     parmu = barmu * 2.;
 739     if (*iprint >= 2) {
 740       fprintf(stderr, "cobyla: increase in PARMU to %12.6E\n", parmu);
 741     }
 742     phi = datmat[mp + np * datmat_dim1] + parmu * datmat[*mpp + np *
 743         datmat_dim1];
 744     i__1 = *n;
 745     for (j = 1; j <= i__1; ++j) {
 746       temp = datmat[mp + j * datmat_dim1] + parmu * datmat[*mpp + j *
 747           datmat_dim1];
 748       if (temp < phi) {
 749         goto L140;
 750       }
 751       if (temp == phi && parmu == 0.f) {
 752         if (datmat[*mpp + j * datmat_dim1] < datmat[*mpp + np *
 753             datmat_dim1]) {
 754           goto L140;
 755         }
 756       }
 757     }
 758   }
 759   prerem = parmu * prerec - sum;
 760
 761 /* Calculate the constraint and objective functions at x(*). Then find the */
 762 /* actual reduction in the merit function. */
 763
 764   i__1 = *n;
 765   for (i__ = 1; i__ <= i__1; ++i__) {
 766     x[i__] = sim[i__ + np * sim_dim1] + dx[i__];
 767   }
 768   ibrnch = 1;
 769   goto L40;
 770 L440:
 771   vmold = datmat[mp + np * datmat_dim1] + parmu * datmat[*mpp + np *
 772       datmat_dim1];
 773   vmnew = f + parmu * resmax;
 774   trured = vmold - vmnew;
 775   if (parmu == 0. && f == datmat[mp + np * datmat_dim1]) {
 776     prerem = prerec;
 777     trured = datmat[*mpp + np * datmat_dim1] - resmax;
 778   }
 779
 780 /* Begin the operations that decide whether x(*) should replace one of the */
 781 /* vertices of the current simplex, the change being mandatory if TRURED is */
 782 /* positive. Firstly, JDROP is set to the index of the vertex that is to be */
 783 /* replaced. */
 784
 785   ratio = 0.;
 786   if (trured <= 0.f) {
 787     ratio = 1.f;
 788   }
 789   jdrop = 0;
 790   i__1 = *n;
 791   for (j = 1; j <= i__1; ++j) {
 792     temp = 0.;
 793     i__2 = *n;
 794     for (i__ = 1; i__ <= i__2; ++i__) {
 795       temp += simi[j + i__ * simi_dim1] * dx[i__];
 796     }
 797     temp = abs(temp);
 798     if (temp > ratio) {
 799       jdrop = j;
 800       ratio = temp;
 801     }
 802     sigbar[j] = temp * vsig[j];
 803   }
 804
 805 /* Calculate the value of ell. */
 806
 807   edgmax = delta * rho;
 808   l = 0;
 809   i__1 = *n;
 810   for (j = 1; j <= i__1; ++j) {
 811     if (sigbar[j] >= parsig || sigbar[j] >= vsig[j]) {
 812       temp = veta[j];
 813       if (trured > 0.) {
 814         temp = 0.;
 815         i__2 = *n;
 816         for (i__ = 1; i__ <= i__2; ++i__) {
 817           d__1 = dx[i__] - sim[i__ + j * sim_dim1];
 818           temp += d__1 * d__1;
 819         }
 820         temp = sqrt(temp);
 821       }
 822       if (temp > edgmax) {
 823         l = j;
 824         edgmax = temp;
 825       }
 826     }
 827   }
 828   if (l > 0) {
 829     jdrop = l;
 830   }
 831   if (jdrop == 0) {
 832     goto L550;
 833   }
 834
 835 /* Revise the simplex by updating the elements of SIM, SIMI and DATMAT. */
 836
 837   temp = 0.;
 838   i__1 = *n;
 839   for (i__ = 1; i__ <= i__1; ++i__) {
 840     sim[i__ + jdrop * sim_dim1] = dx[i__];
 841     temp += simi[jdrop + i__ * simi_dim1] * dx[i__];
 842   }
 843   i__1 = *n;
 844   for (i__ = 1; i__ <= i__1; ++i__) {
 845     simi[jdrop + i__ * simi_dim1] /= temp;
 846   }
 847   i__1 = *n;
 848   for (j = 1; j <= i__1; ++j) {
 849     if (j != jdrop) {
 850       temp = 0.;
 851       i__2 = *n;
 852       for (i__ = 1; i__ <= i__2; ++i__) {
 853         temp += simi[j + i__ * simi_dim1] * dx[i__];
 854       }
 855       i__2 = *n;
 856       for (i__ = 1; i__ <= i__2; ++i__) {
 857         simi[j + i__ * simi_dim1] -= temp * simi[jdrop + i__ *
 858             simi_dim1];
 859       }
 860     }
 861   }
 862   i__1 = *mpp;
 863   for (k = 1; k <= i__1; ++k) {
 864     datmat[k + jdrop * datmat_dim1] = con[k];
 865   }
 866
 867 /* Branch back for further iterations with the current RHO. */
 868
 869   if (trured > 0. && trured >= prerem * .1) {
 870     goto L140;
 871   }
 872 L550:
 873   if (iflag == 0) {
 874     ibrnch = 0;
 875     goto L140;
 876   }
 877
 878   /* SGJ, 2008: convergence tests for function vals; note that current
 879      best val is stored in datmat[mp + np * datmat_dim1], or in f if
 880      ifull == 1, and prev.  best is in *minf.  This seems like a
 881      sensible place to put the convergence test, as it is the same
 882      place that Powell checks the x tolerance (rho > rhoend). */
 883   {
 884        double fbest = ifull == 1 ? f : datmat[mp + np * datmat_dim1];
 885        if (fbest < *minf && nlopt_stop_ftol(stop, fbest, *minf)) {
 886             rc = NLOPT_FTOL_REACHED;
 887             goto L600;
 888        }
 889        *minf = fbest;
 890   }
 891
 892 /* Otherwise reduce RHO if it is not at its least value and reset PARMU. */
 893
 894   if (rho > rhoend) {
 895     rho *= .5;
 896     if (rho <= rhoend * 1.5) {
 897       rho = rhoend;
 898     }
 899     if (parmu > 0.) {
 900       denom = 0.;
 901       i__1 = mp;
 902       for (k = 1; k <= i__1; ++k) {
 903         cmin = datmat[k + np * datmat_dim1];
 904         cmax = cmin;
 905         i__2 = *n;
 906         for (i__ = 1; i__ <= i__2; ++i__) {
 907           d__1 = cmin, d__2 = datmat[k + i__ * datmat_dim1];
 908           cmin = min(d__1,d__2);
 909           d__1 = cmax, d__2 = datmat[k + i__ * datmat_dim1];
 910           cmax = max(d__1,d__2);
 911         }
 912         if (k <= *m && cmin < cmax * .5) {
 913           temp = max(cmax,0.) - cmin;
 914           if (denom <= 0.) {
 915             denom = temp;
 916           } else {
 917             denom = min(denom,temp);
 918           }
 919         }
 920       }
 921       if (denom == 0.) {
 922         parmu = 0.;
 923       } else if (cmax - cmin < parmu * denom) {
 924         parmu = (cmax - cmin) / denom;
 925       }
 926     }
 927     if (*iprint >= 2) {
 928       fprintf(stderr, "cobyla: reduction in RHO to %12.6E and PARMU =%13.6E\n",
 929         rho, parmu);
 930     }
 931     if (*iprint == 2) {
 932       fprintf(stderr, "cobyla: NFVALS = %4d, F =%13.6E, MAXCV =%13.6E\n",
 933         stop->nevals, datmat[mp + np * datmat_dim1], datmat[*mpp + np * datmat_dim1]);
 934
 935       fprintf(stderr, "cobyla: X =");
 936       i__1 = iptem;
 937       for (i__ = 1; i__ <= i__1; ++i__) {
 938         if (i__>1) fprintf(stderr, "  ");
 939         fprintf(stderr, "%13.6E", sim[i__ + np * sim_dim1]);
 940       }
 941       if (iptem < *n) {
 942         i__1 = *n;
 943         for (i__ = iptemp; i__ <= i__1; ++i__) {
 944           if (!((i__-1) % 4)) fprintf(stderr, "\ncobyla:  ");
 945           fprintf(stderr, "%15.6E", x[i__]);
 946         }
 947       }
 948       fprintf(stderr, "\n");
 949     }
 950     goto L140;
 951   }
 952   else
 953        rc = NLOPT_XTOL_REACHED;
 954
 955 /* Return the best calculated values of the variables. */
 956
 957   if (*iprint >= 1) {
 958     fprintf(stderr, "cobyla: normal return.\n");
 959   }
 960   if (ifull == 1) {
 961     goto L620;
 962   }
 963 L600:
 964   i__1 = *n;
 965   for (i__ = 1; i__ <= i__1; ++i__) {
 966     x[i__] = sim[i__ + np * sim_dim1];
 967   }
 968   f = datmat[mp + np * datmat_dim1];
 969   resmax = datmat[*mpp + np * datmat_dim1];
 970 L620:
 971   *minf = f;
 972   if (*iprint >= 1) {
 973     fprintf(stderr, "cobyla: NFVALS = %4d, F =%13.6E, MAXCV =%13.6E\n",
 974             stop->nevals, f, resmax);
 975     i__1 = iptem;
 976     fprintf(stderr, "cobyla: X =");
 977     for (i__ = 1; i__ <= i__1; ++i__) {
 978       if (i__>1) fprintf(stderr, "  ");
 979       fprintf(stderr, "%13.6E", x[i__]);
 980     }
 981     if (iptem < *n) {
 982       i__1 = *n;
 983       for (i__ = iptemp; i__ <= i__1; ++i__) {
 984         if (!((i__-1) % 4)) fprintf(stderr, "\ncobyla:  ");
 985         fprintf(stderr, "%15.6E", x[i__]);
 986       }
 987     }
 988     fprintf(stderr, "\n");
 989   }
 990   return rc;
 991 } /* cobylb */
 992
 993 /* ------------------------------------------------------------------------- */
 994 static void trstlp(int *n, int *m, double *a,
 995     double *b, double *rho, double *dx, int *ifull,
 996     int *iact, double *z__, double *zdota, double *vmultc,
 997      double *sdirn, double *dxnew, double *vmultd)
 998 {
 999   /* System generated locals */
1000   int a_dim1, a_offset, z_dim1, z_offset, i__1, i__2;
1001   double d__1, d__2;
1002
1003   /* Local variables */
1004   double alpha, tempa;
1005   double beta;
1006   double optnew, stpful, sum, tot, acca, accb;
1007   double ratio, vsave, zdotv, zdotw, dd;
1008   double sd;
1009   double sp, ss, resold = 0.0, zdvabs, zdwabs, sumabs, resmax, optold;
1010   double spabs;
1011   double temp, step;
1012   int icount;
1013   int iout, i__, j, k;
1014   int isave;
1015   int kk;
1016   int kl, kp, kw;
1017   int nact, icon = 0, mcon;
1018   int nactx = 0;
1019
1020
1021 /* This subroutine calculates an N-component vector DX by applying the */
1022 /* following two stages. In the first stage, DX is set to the shortest */
1023 /* vector that minimizes the greatest violation of the constraints */
1024 /*   A(1,K)*DX(1)+A(2,K)*DX(2)+...+A(N,K)*DX(N) .GE. B(K), K=2,3,...,M, */
1025 /* subject to the Euclidean length of DX being at most RHO. If its length is */
1026 /* strictly less than RHO, then we use the resultant freedom in DX to */
1027 /* minimize the objective function */
1028 /*      -A(1,M+1)*DX(1)-A(2,M+1)*DX(2)-...-A(N,M+1)*DX(N) */
1029 /* subject to no increase in any greatest constraint violation. This */
1030 /* notation allows the gradient of the objective function to be regarded as */
1031 /* the gradient of a constraint. Therefore the two stages are distinguished */
1032 /* by MCON .EQ. M and MCON .GT. M respectively. It is possible that a */
1033 /* degeneracy may prevent DX from attaining the target length RHO. Then the */
1034 /* value IFULL=0 would be set, but usually IFULL=1 on return. */
1035
1036 /* In general NACT is the number of constraints in the active set and */
1037 /* IACT(1),...,IACT(NACT) are their indices, while the remainder of IACT */
1038 /* contains a permutation of the remaining constraint indices. Further, Z is */
1039 /* an orthogonal matrix whose first NACT columns can be regarded as the */
1040 /* result of Gram-Schmidt applied to the active constraint gradients. For */
1041 /* J=1,2,...,NACT, the number ZDOTA(J) is the scalar product of the J-th */
1042 /* column of Z with the gradient of the J-th active constraint. DX is the */
1043 /* current vector of variables and here the residuals of the active */
1044 /* constraints should be zero. Further, the active constraints have */
1045 /* nonnegative Lagrange multipliers that are held at the beginning of */
1046 /* VMULTC. The remainder of this vector holds the residuals of the inactive */
1047 /* constraints at DX, the ordering of the components of VMULTC being in */
1048 /* agreement with the permutation of the indices of the constraints that is */
1049 /* in IACT. All these residuals are nonnegative, which is achieved by the */
1050 /* shift RESMAX that makes the least residual zero. */
1051
1052 /* Initialize Z and some other variables. The value of RESMAX will be */
1053 /* appropriate to DX=0, while ICON will be the index of a most violated */
1054 /* constraint if RESMAX is positive. Usually during the first stage the */
1055 /* vector SDIRN gives a search direction that reduces all the active */
1056 /* constraint violations by one simultaneously. */
1057
1058   /* Parameter adjustments */
1059   z_dim1 = *n;
1060   z_offset = 1 + z_dim1 * 1;
1061   z__ -= z_offset;
1062   a_dim1 = *n;
1063   a_offset = 1 + a_dim1 * 1;
1064   a -= a_offset;
1065   --b;
1066   --dx;
1067   --iact;
1068   --zdota;
1069   --vmultc;
1070   --sdirn;
1071   --dxnew;
1072   --vmultd;
1073
1074   /* Function Body */
1075   *ifull = 1;
1076   mcon = *m;
1077   nact = 0;
1078   resmax = 0.;
1079   i__1 = *n;
1080   for (i__ = 1; i__ <= i__1; ++i__) {
1081     i__2 = *n;
1082     for (j = 1; j <= i__2; ++j) {
1083       z__[i__ + j * z_dim1] = 0.;
1084     }
1085     z__[i__ + i__ * z_dim1] = 1.;
1086     dx[i__] = 0.;
1087   }
1088   if (*m >= 1) {
1089     i__1 = *m;
1090     for (k = 1; k <= i__1; ++k) {
1091       if (b[k] > resmax) {
1092         resmax = b[k];
1093         icon = k;
1094       }
1095     }
1096     i__1 = *m;
1097     for (k = 1; k <= i__1; ++k) {
1098       iact[k] = k;
1099       vmultc[k] = resmax - b[k];
1100     }
1101   }
1102   if (resmax == 0.) {
1103     goto L480;
1104   }
1105   i__1 = *n;
1106   for (i__ = 1; i__ <= i__1; ++i__) {
1107     sdirn[i__] = 0.;
1108   }
1109
1110 /* End the current stage of the calculation if 3 consecutive iterations */
1111 /* have either failed to reduce the best calculated value of the objective */
1112 /* function or to increase the number of active constraints since the best */
1113 /* value was calculated. This strategy prevents cycling, but there is a */
1114 /* remote possibility that it will cause premature termination. */
1115
1116 L60:
1117   optold = 0.;
1118   icount = 0;
1119 L70:
1120   if (mcon == *m) {
1121     optnew = resmax;
1122   } else {
1123     optnew = 0.;
1124     i__1 = *n;
1125     for (i__ = 1; i__ <= i__1; ++i__) {
1126       optnew -= dx[i__] * a[i__ + mcon * a_dim1];
1127     }
1128   }
1129   if (icount == 0 || optnew < optold) {
1130     optold = optnew;
1131     nactx = nact;
1132     icount = 3;
1133   } else if (nact > nactx) {
1134     nactx = nact;
1135     icount = 3;
1136   } else {
1137     --icount;
1138     if (icount == 0) {
1139       goto L490;
1140     }
1141   }
1142
1143 /* If ICON exceeds NACT, then we add the constraint with index IACT(ICON) to */
1144 /* the active set. Apply Givens rotations so that the last N-NACT-1 columns */
1145 /* of Z are orthogonal to the gradient of the new constraint, a scalar */
1146 /* product being set to zero if its nonzero value could be due to computer */
1147 /* rounding errors. The array DXNEW is used for working space. */
1148
1149   if (icon <= nact) {
1150     goto L260;
1151   }
1152   kk = iact[icon];
1153   i__1 = *n;
1154   for (i__ = 1; i__ <= i__1; ++i__) {
1155     dxnew[i__] = a[i__ + kk * a_dim1];
1156   }
1157   tot = 0.;
1158   k = *n;
1159 L100:
1160   if (k > nact) {
1161     sp = 0.;
1162     spabs = 0.;
1163     i__1 = *n;
1164     for (i__ = 1; i__ <= i__1; ++i__) {
1165       temp = z__[i__ + k * z_dim1] * dxnew[i__];
1166       sp += temp;
1167       spabs += abs(temp);
1168     }
1169     acca = spabs + abs(sp) * .1;
1170     accb = spabs + abs(sp) * .2;
1171     if (spabs >= acca || acca >= accb) {
1172       sp = 0.;
1173     }
1174     if (tot == 0.) {
1175       tot = sp;
1176     } else {
1177       kp = k + 1;
1178       temp = sqrt(sp * sp + tot * tot);
1179       alpha = sp / temp;
1180       beta = tot / temp;
1181       tot = temp;
1182       i__1 = *n;
1183       for (i__ = 1; i__ <= i__1; ++i__) {
1184         temp = alpha * z__[i__ + k * z_dim1] + beta * z__[i__ + kp *
1185             z_dim1];
1186         z__[i__ + kp * z_dim1] = alpha * z__[i__ + kp * z_dim1] -
1187             beta * z__[i__ + k * z_dim1];
1188         z__[i__ + k * z_dim1] = temp;
1189       }
1190     }
1191     --k;
1192     goto L100;
1193   }
1194
1195 /* Add the new constraint if this can be done without a deletion from the */
1196 /* active set. */
1197
1198   if (tot != 0.) {
1199     ++nact;
1200     zdota[nact] = tot;
1201     vmultc[icon] = vmultc[nact];
1202     vmultc[nact] = 0.;
1203     goto L210;
1204   }
1205
1206 /* The next instruction is reached if a deletion has to be made from the */
1207 /* active set in order to make room for the new active constraint, because */
1208 /* the new constraint gradient is a linear combination of the gradients of */
1209 /* the old active constraints. Set the elements of VMULTD to the multipliers */
1210 /* of the linear combination. Further, set IOUT to the index of the */
1211 /* constraint to be deleted, but branch if no suitable index can be found. */
1212
1213   ratio = -1.;
1214   k = nact;
1215 L130:
1216   zdotv = 0.;
1217   zdvabs = 0.;
1218   i__1 = *n;
1219   for (i__ = 1; i__ <= i__1; ++i__) {
1220     temp = z__[i__ + k * z_dim1] * dxnew[i__];
1221     zdotv += temp;
1222     zdvabs += abs(temp);
1223   }
1224   acca = zdvabs + abs(zdotv) * .1;
1225   accb = zdvabs + abs(zdotv) * .2;
1226   if (zdvabs < acca && acca < accb) {
1227     temp = zdotv / zdota[k];
1228     if (temp > 0. && iact[k] <= *m) {
1229       tempa = vmultc[k] / temp;
1230       if (ratio < 0. || tempa < ratio) {
1231         ratio = tempa;
1232         iout = k;
1233       }
1234     }
1235     if (k >= 2) {
1236       kw = iact[k];
1237       i__1 = *n;
1238       for (i__ = 1; i__ <= i__1; ++i__) {
1239         dxnew[i__] -= temp * a[i__ + kw * a_dim1];
1240       }
1241     }
1242     vmultd[k] = temp;
1243   } else {
1244     vmultd[k] = 0.;
1245   }
1246   --k;
1247   if (k > 0) {
1248     goto L130;
1249   }
1250   if (ratio < 0.) {
1251     goto L490;
1252   }
1253
1254 /* Revise the Lagrange multipliers and reorder the active constraints so */
1255 /* that the one to be replaced is at the end of the list. Also calculate the */
1256 /* new value of ZDOTA(NACT) and branch if it is not acceptable. */
1257
1258   i__1 = nact;
1259   for (k = 1; k <= i__1; ++k) {
1260     d__1 = 0., d__2 = vmultc[k] - ratio * vmultd[k];
1261     vmultc[k] = max(d__1,d__2);
1262   }
1263   if (icon < nact) {
1264     isave = iact[icon];
1265     vsave = vmultc[icon];
1266     k = icon;
1267 L170:
1268     kp = k + 1;
1269     kw = iact[kp];
1270     sp = 0.;
1271     i__1 = *n;
1272     for (i__ = 1; i__ <= i__1; ++i__) {
1273       sp += z__[i__ + k * z_dim1] * a[i__ + kw * a_dim1];
1274     }
1275     d__1 = zdota[kp];
1276     temp = sqrt(sp * sp + d__1 * d__1);
1277     alpha = zdota[kp] / temp;
1278     beta = sp / temp;
1279     zdota[kp] = alpha * zdota[k];
1280     zdota[k] = temp;
1281     i__1 = *n;
1282     for (i__ = 1; i__ <= i__1; ++i__) {
1283       temp = alpha * z__[i__ + kp * z_dim1] + beta * z__[i__ + k *
1284           z_dim1];
1285       z__[i__ + kp * z_dim1] = alpha * z__[i__ + k * z_dim1] - beta *
1286           z__[i__ + kp * z_dim1];
1287       z__[i__ + k * z_dim1] = temp;
1288     }
1289     iact[k] = kw;
1290     vmultc[k] = vmultc[kp];
1291     k = kp;
1292     if (k < nact) {
1293       goto L170;
1294     }
1295     iact[k] = isave;
1296     vmultc[k] = vsave;
1297   }
1298   temp = 0.;
1299   i__1 = *n;
1300   for (i__ = 1; i__ <= i__1; ++i__) {
1301     temp += z__[i__ + nact * z_dim1] * a[i__ + kk * a_dim1];
1302   }
1303   if (temp == 0.) {
1304     goto L490;
1305   }
1306   zdota[nact] = temp;
1307   vmultc[icon] = 0.;
1308   vmultc[nact] = ratio;
1309
1310 /* Update IACT and ensure that the objective function continues to be */
1311 /* treated as the last active constraint when MCON>M. */
1312
1313 L210:
1314   iact[icon] = iact[nact];
1315   iact[nact] = kk;
1316   if (mcon > *m && kk != mcon) {
1317     k = nact - 1;
1318     sp = 0.;
1319     i__1 = *n;
1320     for (i__ = 1; i__ <= i__1; ++i__) {
1321       sp += z__[i__ + k * z_dim1] * a[i__ + kk * a_dim1];
1322     }
1323     d__1 = zdota[nact];
1324     temp = sqrt(sp * sp + d__1 * d__1);
1325     alpha = zdota[nact] / temp;
1326     beta = sp / temp;
1327     zdota[nact] = alpha * zdota[k];
1328     zdota[k] = temp;
1329     i__1 = *n;
1330     for (i__ = 1; i__ <= i__1; ++i__) {
1331       temp = alpha * z__[i__ + nact * z_dim1] + beta * z__[i__ + k *
1332           z_dim1];
1333       z__[i__ + nact * z_dim1] = alpha * z__[i__ + k * z_dim1] - beta *
1334           z__[i__ + nact * z_dim1];
1335       z__[i__ + k * z_dim1] = temp;
1336     }
1337     iact[nact] = iact[k];
1338     iact[k] = kk;
1339     temp = vmultc[k];
1340     vmultc[k] = vmultc[nact];
1341     vmultc[nact] = temp;
1342   }
1343
1344 /* If stage one is in progress, then set SDIRN to the direction of the next */
1345 /* change to the current vector of variables. */
1346
1347   if (mcon > *m) {
1348     goto L320;
1349   }
1350   kk = iact[nact];
1351   temp = 0.;
1352   i__1 = *n;
1353   for (i__ = 1; i__ <= i__1; ++i__) {
1354     temp += sdirn[i__] * a[i__ + kk * a_dim1];
1355   }
1356   temp += -1.;
1357   temp /= zdota[nact];
1358   i__1 = *n;
1359   for (i__ = 1; i__ <= i__1; ++i__) {
1360     sdirn[i__] -= temp * z__[i__ + nact * z_dim1];
1361   }
1362   goto L340;
1363
1364 /* Delete the constraint that has the index IACT(ICON) from the active set. */
1365
1366 L260:
1367   if (icon < nact) {
1368     isave = iact[icon];
1369     vsave = vmultc[icon];
1370     k = icon;
1371 L270:
1372     kp = k + 1;
1373     kk = iact[kp];
1374     sp = 0.;
1375     i__1 = *n;
1376     for (i__ = 1; i__ <= i__1; ++i__) {
1377       sp += z__[i__ + k * z_dim1] * a[i__ + kk * a_dim1];
1378     }
1379     d__1 = zdota[kp];
1380     temp = sqrt(sp * sp + d__1 * d__1);
1381     alpha = zdota[kp] / temp;
1382     beta = sp / temp;
1383     zdota[kp] = alpha * zdota[k];
1384     zdota[k] = temp;
1385     i__1 = *n;
1386     for (i__ = 1; i__ <= i__1; ++i__) {
1387       temp = alpha * z__[i__ + kp * z_dim1] + beta * z__[i__ + k *
1388           z_dim1];
1389       z__[i__ + kp * z_dim1] = alpha * z__[i__ + k * z_dim1] - beta *
1390           z__[i__ + kp * z_dim1];
1391       z__[i__ + k * z_dim1] = temp;
1392     }
1393     iact[k] = kk;
1394     vmultc[k] = vmultc[kp];
1395     k = kp;
1396     if (k < nact) {
1397       goto L270;
1398     }
1399     iact[k] = isave;
1400     vmultc[k] = vsave;
1401   }
1402   --nact;
1403
1404 /* If stage one is in progress, then set SDIRN to the direction of the next */
1405 /* change to the current vector of variables. */
1406
1407   if (mcon > *m) {
1408     goto L320;
1409   }
1410   temp = 0.;
1411   i__1 = *n;
1412   for (i__ = 1; i__ <= i__1; ++i__) {
1413     temp += sdirn[i__] * z__[i__ + (nact + 1) * z_dim1];
1414   }
1415   i__1 = *n;
1416   for (i__ = 1; i__ <= i__1; ++i__) {
1417     sdirn[i__] -= temp * z__[i__ + (nact + 1) * z_dim1];
1418   }
1419   goto L340;
1420
1421 /* Pick the next search direction of stage two. */
1422
1423 L320:
1424   temp = 1. / zdota[nact];
1425   i__1 = *n;
1426   for (i__ = 1; i__ <= i__1; ++i__) {
1427     sdirn[i__] = temp * z__[i__ + nact * z_dim1];
1428   }
1429
1430 /* Calculate the step to the boundary of the trust region or take the step */
1431 /* that reduces RESMAX to zero. The two statements below that include the */
1432 /* factor 1.0E-6 prevent some harmless underflows that occurred in a test */
1433 /* calculation. Further, we skip the step if it could be zero within a */
1434 /* reasonable tolerance for computer rounding errors. */
1435
1436 L340:
1437   dd = *rho * *rho;
1438   sd = 0.;
1439   ss = 0.;
1440   i__1 = *n;
1441   for (i__ = 1; i__ <= i__1; ++i__) {
1442     if ((d__1 = dx[i__], abs(d__1)) >= *rho * 1e-6f) {
1443       d__2 = dx[i__];
1444       dd -= d__2 * d__2;
1445     }
1446     sd += dx[i__] * sdirn[i__];
1447     d__1 = sdirn[i__];
1448     ss += d__1 * d__1;
1449   }
1450   if (dd <= 0.) {
1451     goto L490;
1452   }
1453   temp = sqrt(ss * dd);
1454   if (abs(sd) >= temp * 1e-6f) {
1455     temp = sqrt(ss * dd + sd * sd);
1456   }
1457   stpful = dd / (temp + sd);
1458   step = stpful;
1459   if (mcon == *m) {
1460     acca = step + resmax * .1;
1461     accb = step + resmax * .2;
1462     if (step >= acca || acca >= accb) {
1463       goto L480;
1464     }
1465     step = min(step,resmax);
1466   }
1467
1468 /* Set DXNEW to the new variables if STEP is the steplength, and reduce */
1469 /* RESMAX to the corresponding maximum residual if stage one is being done. */
1470 /* Because DXNEW will be changed during the calculation of some Lagrange */
1471 /* multipliers, it will be restored to the following value later. */
1472
1473   i__1 = *n;
1474   for (i__ = 1; i__ <= i__1; ++i__) {
1475     dxnew[i__] = dx[i__] + step * sdirn[i__];
1476   }
1477   if (mcon == *m) {
1478     resold = resmax;
1479     resmax = 0.;
1480     i__1 = nact;
1481     for (k = 1; k <= i__1; ++k) {
1482       kk = iact[k];
1483       temp = b[kk];
1484       i__2 = *n;
1485       for (i__ = 1; i__ <= i__2; ++i__) {
1486         temp -= a[i__ + kk * a_dim1] * dxnew[i__];
1487       }
1488       resmax = max(resmax,temp);
1489     }
1490   }
1491
1492 /* Set VMULTD to the VMULTC vector that would occur if DX became DXNEW. A */
1493 /* device is included to force VMULTD(K)=0.0 if deviations from this value */
1494 /* can be attributed to computer rounding errors. First calculate the new */
1495 /* Lagrange multipliers. */
1496
1497   k = nact;
1498 L390:
1499   zdotw = 0.;
1500   zdwabs = 0.;
1501   i__1 = *n;
1502   for (i__ = 1; i__ <= i__1; ++i__) {
1503     temp = z__[i__ + k * z_dim1] * dxnew[i__];
1504     zdotw += temp;
1505     zdwabs += abs(temp);
1506   }
1507   acca = zdwabs + abs(zdotw) * .1;
1508   accb = zdwabs + abs(zdotw) * .2;
1509   if (zdwabs >= acca || acca >= accb) {
1510     zdotw = 0.;
1511   }
1512   vmultd[k] = zdotw / zdota[k];
1513   if (k >= 2) {
1514     kk = iact[k];
1515     i__1 = *n;
1516     for (i__ = 1; i__ <= i__1; ++i__) {
1517       dxnew[i__] -= vmultd[k] * a[i__ + kk * a_dim1];
1518     }
1519     --k;
1520     goto L390;
1521   }
1522   if (mcon > *m) {
1523     d__1 = 0., d__2 = vmultd[nact];
1524     vmultd[nact] = max(d__1,d__2);
1525   }
1526
1527 /* Complete VMULTC by finding the new constraint residuals. */
1528
1529   i__1 = *n;
1530   for (i__ = 1; i__ <= i__1; ++i__) {
1531     dxnew[i__] = dx[i__] + step * sdirn[i__];
1532   }
1533   if (mcon > nact) {
1534     kl = nact + 1;
1535     i__1 = mcon;
1536     for (k = kl; k <= i__1; ++k) {
1537       kk = iact[k];
1538       sum = resmax - b[kk];
1539       sumabs = resmax + (d__1 = b[kk], abs(d__1));
1540       i__2 = *n;
1541       for (i__ = 1; i__ <= i__2; ++i__) {
1542         temp = a[i__ + kk * a_dim1] * dxnew[i__];
1543         sum += temp;
1544         sumabs += abs(temp);
1545       }
1546       acca = sumabs + abs(sum) * .1f;
1547       accb = sumabs + abs(sum) * .2f;
1548       if (sumabs >= acca || acca >= accb) {
1549         sum = 0.f;
1550       }
1551       vmultd[k] = sum;
1552     }
1553   }
1554
1555 /* Calculate the fraction of the step from DX to DXNEW that will be taken. */
1556
1557   ratio = 1.;
1558   icon = 0;
1559   i__1 = mcon;
1560   for (k = 1; k <= i__1; ++k) {
1561     if (vmultd[k] < 0.) {
1562       temp = vmultc[k] / (vmultc[k] - vmultd[k]);
1563       if (temp < ratio) {
1564         ratio = temp;
1565         icon = k;
1566       }
1567     }
1568   }
1569
1570 /* Update DX, VMULTC and RESMAX. */
1571
1572   temp = 1. - ratio;
1573   i__1 = *n;
1574   for (i__ = 1; i__ <= i__1; ++i__) {
1575     dx[i__] = temp * dx[i__] + ratio * dxnew[i__];
1576   }
1577   i__1 = mcon;
1578   for (k = 1; k <= i__1; ++k) {
1579     d__1 = 0., d__2 = temp * vmultc[k] + ratio * vmultd[k];
1580     vmultc[k] = max(d__1,d__2);
1581   }
1582   if (mcon == *m) {
1583     resmax = resold + ratio * (resmax - resold);
1584   }
1585
1586 /* If the full step is not acceptable then begin another iteration. */
1587 /* Otherwise switch to stage two or end the calculation. */
1588
1589   if (icon > 0) {
1590     goto L70;
1591   }
1592   if (step == stpful) {
1593     goto L500;
1594   }
1595 L480:
1596   mcon = *m + 1;
1597   icon = mcon;
1598   iact[mcon] = mcon;
1599   vmultc[mcon] = 0.;
1600   goto L60;
1601
1602 /* We employ any freedom that may be available to reduce the objective */
1603 /* function before returning a DX whose length is less than RHO. */
1604
1605 L490:
1606   if (mcon == *m) {
1607     goto L480;
1608   }
1609   *ifull = 0;
1610 L500:
1611   return;
1612 } /* trstlp */