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 /* SGJ, 2008: modified COBYLA code to take explicit account of bound
  45    constraints.  Since bound constraints are linear, these should
  46    already be handled exactly when COBYLA optimizes it linear model.
  47    However, they were not handled properly when COBYLA created its
  48    initial simplex, or when COBYLA updated unacceptable simplices.
  49    Since NLopt guarantees that the objective will not be evaluated
  50    outside the bound constraints, this required us to handle such
  51    points by putting a slope discontinuity into the objective &
  52    constraints (below), which slows convergence considerably for
  53    smooth functions.  Instead, handling them explicitly prevents this
  54    problem */
  55 #define ENFORCE_BOUNDS 1
  56
  57 #define MIN2(a,b) ((a) <= (b) ? (a) : (b))
  58 #define MAX2(a,b) ((a) >= (b) ? (a) : (b))
  59
  60 #define U(n) ((unsigned) (n))
  61
  62 /**************************************************************************/
  63 /* SGJ, 2008: NLopt-style interface function: */
  64
  65 typedef struct {
  66      nlopt_func f;
  67      void *f_data;
  68      unsigned m_orig;
  69      nlopt_constraint *fc;
  70      unsigned p;
  71      nlopt_constraint *h;
  72      double *xtmp;
  73      double *lb, *ub;
  74      double *con_tol, *scale;
  75      nlopt_stopping *stop;
  76 } func_wrap_state;
  77
  78 static int func_wrap(int ni, int mi, double *x, double *f, double *con,
  79                      func_wrap_state *s)
  80 {
  81      unsigned n = U(ni);
  82      unsigned i, j, k;
  83      double *xtmp = s->xtmp;
  84      const double *lb = s->lb, *ub = s->ub;
  85
  86      (void) mi; /* unused */
  87
  88      /* in nlopt, we guarante that the function is never evaluated outside
  89         the lb and ub bounds, so we need force this with xtmp ... note
  90         that this leads to discontinuity in the first derivative, which
  91         slows convergence if we don't enable the ENFORCE_BOUNDS feature
  92         above. */
  93      for (j = 0; j < n; ++j) {
  94           if (x[j] < lb[j]) xtmp[j] = lb[j];
  95           else if (x[j] > ub[j]) xtmp[j] = ub[j];
  96           else xtmp[j] = x[j];
  97      }
  98      nlopt_unscale(n, s->scale, xtmp, xtmp);
  99
 100      *f = s->f(n, xtmp, NULL, s->f_data);
 101      if (nlopt_stop_forced(s->stop)) return 1;
 102      i = 0;
 103      for (j = 0; j < s->m_orig; ++j) {
 104           nlopt_eval_constraint(con + i, NULL, s->fc+j, n, xtmp);
 105           if (nlopt_stop_forced(s->stop)) return 1;
 106           for (k = 0; k < s->fc[j].m; ++k)
 107                con[i + k] = -con[i + k];
 108           i += s->fc[j].m;
 109      }
 110      for (j = 0; j < s->p; ++j) {
 111           nlopt_eval_constraint(con + i, NULL, s->h+j, n, xtmp);
 112           if (nlopt_stop_forced(s->stop)) return 1;
 113           for (k = 0; k < s->h[j].m; ++k)
 114                con[(i + s->h[j].m) + k] = -con[i + k];
 115           i += 2 * s->h[j].m;
 116      }
 117      for (j = 0; j < n; ++j) {
 118           if (!nlopt_isinf(lb[j]))
 119                con[i++] = x[j] - lb[j];
 120           if (!nlopt_isinf(ub[j]))
 121                con[i++] = ub[j] - x[j];
 122      }
 123      return 0;
 124 }
 125
 126 /*
 127  * Verbosity level
 128  */
 129 typedef enum {
 130   COBYLA_MSG_NONE = 0, /* No messages */
 131   COBYLA_MSG_EXIT = 1, /* Exit reasons */
 132   COBYLA_MSG_ITER = 2, /* Rho and Sigma changes */
 133   COBYLA_MSG_INFO = 3 /* Informational messages */
 134 } cobyla_message;
 135
 136 /*
 137  * A function as required by cobyla
 138  * state is a void pointer provided to the function at each call
 139  *
 140  * n     : the number of variables
 141  * m     : the number of constraints
 142  * x     : on input, then vector of variables (should not be modified)
 143  * f     : on output, the value of the function
 144  * con   : on output, the value of the constraints (vector of size m)
 145  * state : on input, the value of the state variable as provided to cobyla
 146  *
 147  * COBYLA will try to make all the values of the constraints positive.
 148  * So if you want to input a constraint j such as x[i] <= MAX, set:
 149  *   con[j] = MAX - x[i]
 150  * The function must returns 0 if no error occurs or 1 to immediately end the
 151  * minimization.
 152  *
 153  */
 154 typedef int cobyla_function(int n, int m, double *x, double *f, double *con,
 155   func_wrap_state *state);
 156
 157 /*
 158  * cobyla : minimize a function subject to constraints
 159  *
 160  * n         : number of variables (>=0)
 161  * m         : number of constraints (>=0)
 162  * x         : on input, initial estimate ; on output, the solution
 163  * minf      : on output, minimum objective function
 164  * rhobeg    : a reasonable initial change to the variables
 165  * stop      : the NLopt stopping criteria
 166  * lb, ub    : lower and upper bounds on x
 167  * message   : see the cobyla_message enum
 168  * calcfc    : the function to minimize (see cobyla_function)
 169  * state     : used by function (see cobyla_function)
 170  *
 171  * The cobyla function returns the usual nlopt_result codes.
 172  *
 173  */
 174 extern nlopt_result cobyla(int n, int m, double *x, double *minf, double rhobeg, double rhoend, nlopt_stopping *stop, const double *lb, const double *ub,
 175   int message, cobyla_function *calcfc, func_wrap_state *state);
 176
 177 nlopt_result cobyla_minimize(unsigned n, nlopt_func f, void *f_data,
 178                              unsigned m, nlopt_constraint *fc,
 179                              unsigned p, nlopt_constraint *h,
 180                              const double *lb, const double *ub, /* bounds */
 181                              double *x, /* in: initial guess, out: minimizer */
 182                              double *minf,
 183                              nlopt_stopping *stop,
 184                              const double *dx)
 185 {
 186      unsigned i, j;
 187      func_wrap_state s;
 188      nlopt_result ret;
 189      double rhobeg, rhoend;
 190
 191      s.f = f; s.f_data = f_data;
 192      s.m_orig = m;
 193      s.fc = fc;
 194      s.p = p;
 195      s.h = h;
 196      s.stop = stop;
 197      s.lb = s.ub = s.xtmp = s.con_tol = s.scale = NULL;
 198
 199      s.scale = nlopt_compute_rescaling(n, dx);
 200      if (!s.scale) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
 201
 202      s.lb = nlopt_new_rescaled(n, s.scale, lb);
 203      if (!s.lb) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
 204      s.ub = nlopt_new_rescaled(n, s.scale, ub);
 205      if (!s.ub) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
 206
 207      s.xtmp = (double *) malloc(sizeof(double) * n);
 208      if (!s.xtmp) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
 209
 210      /* SGJ, 2008: compute rhoend from NLopt stop info */
 211      rhobeg = dx[0] / s.scale[0];
 212      rhoend = stop->xtol_rel * (rhobeg);
 213      for (j = 0; j < n; ++j)
 214           if (rhoend < stop->xtol_abs[j] / s.scale[j])
 215                rhoend = stop->xtol_abs[j] / s.scale[j];
 216
 217      /* each equality constraint gives two inequality constraints */
 218      m = nlopt_count_constraints(m, fc) + 2 * nlopt_count_constraints(p, h);
 219
 220      /* add constraints for lower/upper bounds (if any) */
 221      for (j = 0; j < n; ++j) {
 222           if (!nlopt_isinf(lb[j]))
 223                ++m;
 224           if (!nlopt_isinf(ub[j]))
 225                ++m;
 226      }
 227
 228      s.con_tol = (double *) malloc(sizeof(double) * m);
 229      if (m && !s.con_tol) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
 230
 231      for (j = 0; j < m; ++j) s.con_tol[j] = 0;
 232      for (j = i = 0; i < s.m_orig; ++i) {
 233           unsigned ji = j, jnext = j + fc[i].m;
 234           for (; j < jnext; ++j) s.con_tol[j] = fc[i].tol[j - ji];
 235      }
 236      for (i = 0; i < s.p; ++i) {
 237           unsigned ji = j, jnext = j + h[i].m;
 238           for (; j < jnext; ++j) s.con_tol[j] = h[i].tol[j - ji];
 239           ji = j; jnext = j + h[i].m;
 240           for (; j < jnext; ++j) s.con_tol[j] = h[i].tol[j - ji];
 241      }
 242
 243      nlopt_rescale(n, s.scale, x, x);
 244      ret = cobyla((int) n, (int) m, x, minf, rhobeg, rhoend,
 245                   stop, s.lb, s.ub, COBYLA_MSG_NONE,
 246                   func_wrap, &s);
 247      nlopt_unscale(n, s.scale, x, x);
 248
 249      /* make sure e.g. rounding errors didn't push us slightly out of bounds */
 250      for (j = 0; j < n; ++j) {
 251           if (x[j] < lb[j]) x[j] = lb[j];
 252           if (x[j] > ub[j]) x[j] = ub[j];
 253      }
 254
 255 done:
 256      free(s.con_tol);
 257      free(s.xtmp);
 258      free(s.ub);
 259      free(s.lb);
 260      free(s.scale);
 261      return ret;
 262 }
 263
 264 /**************************************************************************/
 265
 266 /* SGJ, 2010: I find that it seems to increase robustness of the algorithm
 267    if, on "simplex" steps (which are intended to improve the linear
 268    independence of the simplex, not improve the objective), we multiply
 269    the steps by pseudorandom numbers in [0,1).  Intuitively, pseudorandom
 270    steps are likely to avoid any accidental dependency in the simplex
 271    vertices.  However, since I still want COBYLA to be a deterministic
 272    algorithm, and I don't care all that much about the quality of the
 273    randomness here, I implement a simple linear congruential generator
 274    rather than use nlopt_urand, and set the initial seed deterministically. */
 275
 276 #if defined(HAVE_STDINT_H)
 277 #  include <stdint.h>
 278 #endif
 279
 280 #ifndef HAVE_UINT32_T
 281 #  if SIZEOF_UNSIGNED_LONG == 4
 282 typedef unsigned long uint32_t;
 283 #  elif SIZEOF_UNSIGNED_INT == 4
 284 typedef unsigned int uint32_t;
 285 #  else
 286 #    error No 32-bit unsigned integer type
 287 #  endif
 288 #endif
 289
 290 /* a simple linear congruential generator */
 291
 292 static uint32_t lcg_rand(uint32_t *seed)
 293 {
 294      return (*seed = *seed * 1103515245 + 12345);
 295 }
 296
 297 static double lcg_urand(uint32_t *seed, double a, double b)
 298 {
 299      return a + lcg_rand(seed) * (b - a) / ((uint32_t) -1);
 300 }
 301
 302 /**************************************************************************/
 303
 304 static nlopt_result cobylb(int *n, int *m, int *mpp, double *x, double *minf, double *rhobeg, double rhoend,
 305   nlopt_stopping *stop, const double *lb, const double *ub, int *iprint, double *con, double *sim,
 306   double *simi, double *datmat, double *a, double *vsig, double *veta,
 307   double *sigbar, double *dx, double *w, int *iact, cobyla_function *calcfc,
 308   func_wrap_state *state);
 309 static nlopt_result trstlp(int *n, int *m, double *a, double *b, double *rho,
 310   double *dx, int *ifull, int *iact, double *z__, double *zdota, double *vmultc,
 311   double *sdirn, double *dxnew, double *vmultd);
 312
 313 /* ------------------------------------------------------------------------ */
 314
 315 nlopt_result cobyla(int n, int m, double *x, double *minf, double rhobeg, double rhoend, nlopt_stopping *stop, const double *lb, const double *ub, int iprint,
 316   cobyla_function *calcfc, func_wrap_state *state)
 317 {
 318   int icon, isim, isigb, idatm, iveta, isimi, ivsig, iwork, ia, idx, mpp;
 319   int *iact;
 320   double *w;
 321   nlopt_result rc;
 322
 323 /*
 324  * This subroutine minimizes an objective function F(X) subject to M
 325  * inequality constraints on X, where X is a vector of variables that has
 326  * N components. The algorithm employs linear approximations to the
 327  * objective and constraint functions, the approximations being formed by
 328  * linear interpolation at N+1 points in the space of the variables.
 329  * We regard these interpolation points as vertices of a simplex. The
 330  * parameter RHO controls the size of the simplex and it is reduced
 331  * automatically from RHOBEG to RHOEND. For each RHO the subroutine tries
 332  * to achieve a good vector of variables for the current size, and then
 333  * RHO is reduced until the value RHOEND is reached. Therefore RHOBEG and
 334  * RHOEND should be set to reasonable initial changes to and the required
 335  * accuracy in the variables respectively, but this accuracy should be
 336  * viewed as a subject for experimentation because it is not guaranteed.
 337  * The subroutine has an advantage over many of its competitors, however,
 338  * which is that it treats each constraint individually when calculating
 339  * a change to the variables, instead of lumping the constraints together
 340  * into a single penalty function. The name of the subroutine is derived
 341  * from the phrase Constrained Optimization BY Linear Approximations.
 342  *
 343  * The user must set the values of N, M, RHOBEG and RHOEND, and must
 344  * provide an initial vector of variables in X. Further, the value of
 345  * IPRINT should be set to 0, 1, 2 or 3, which controls the amount of
 346  * printing during the calculation. Specifically, there is no output if
 347  * IPRINT=0 and there is output only at the end of the calculation if
 348  * IPRINT=1. Otherwise each new value of RHO and SIGMA is printed.
 349  * Further, the vector of variables and some function information are
 350  * given either when RHO is reduced or when each new value of F(X) is
 351  * computed in the cases IPRINT=2 or IPRINT=3 respectively. Here SIGMA
 352  * is a penalty parameter, it being assumed that a change to X is an
 353  * improvement if it reduces the merit function
 354  *      F(X)+SIGMA*MAX(0.0,-C1(X),-C2(X),...,-CM(X)),
 355  * where C1,C2,...,CM denote the constraint functions that should become
 356  * nonnegative eventually, at least to the precision of RHOEND. In the
 357  * printed output the displayed term that is multiplied by SIGMA is
 358  * called MAXCV, which stands for 'MAXimum Constraint Violation'. The
 359  * argument MAXFUN is an int variable that must be set by the user to a
 360  * limit on the number of calls of CALCFC, the purpose of this routine being
 361  * given below. The value of MAXFUN will be altered to the number of calls
 362  * of CALCFC that are made. The arguments W and IACT provide real and
 363  * int arrays that are used as working space. Their lengths must be at
 364  * least N*(3*N+2*M+11)+4*M+6 and M+1 respectively.
 365  *
 366  * In order to define the objective and constraint functions, we require
 367  * a subroutine that has the name and arguments
 368  *      SUBROUTINE CALCFC (N,M,X,F,CON)
 369  *      DIMENSION X(*),CON(*)  .
 370  * The values of N and M are fixed and have been defined already, while
 371  * X is now the current vector of variables. The subroutine should return
 372  * the objective and constraint functions at X in F and CON(1),CON(2),
 373  * ...,CON(M). Note that we are trying to adjust X so that F(X) is as
 374  * small as possible subject to the constraint functions being nonnegative.
 375  *
 376  * Partition the working space array W to provide the storage that is needed
 377  * for the main calculation.
 378  */
 379
 380   stop->nevals = 0;
 381
 382   if (n == 0)
 383   {
 384     if (iprint>=1) fprintf(stderr, "cobyla: N==0.\n");
 385     return NLOPT_SUCCESS;
 386   }
 387
 388   if (n < 0 || m < 0)
 389   {
 390     if (iprint>=1) fprintf(stderr, "cobyla: N<0 or M<0.\n");
 391     return NLOPT_INVALID_ARGS;
 392   }
 393
 394   /* workspace allocation */
 395   w = (double*) malloc(U(n*(3*n+2*m+11)+4*m+6)*sizeof(*w));
 396   if (w == NULL)
 397   {
 398     if (iprint>=1) fprintf(stderr, "cobyla: memory allocation error.\n");
 399     return NLOPT_OUT_OF_MEMORY;
 400   }
 401   iact = (int*)malloc(U(m+1)*sizeof(*iact));
 402   if (iact == NULL)
 403   {
 404     if (iprint>=1) fprintf(stderr, "cobyla: memory allocation error.\n");
 405     free(w);
 406     return NLOPT_OUT_OF_MEMORY;
 407   }
 408
 409   /* Parameter adjustments */
 410   --iact;
 411   --w;
 412   --x;
 413   --lb; --ub;
 414
 415   /* Function Body */
 416   mpp = m + 2;
 417   icon = 1;
 418   isim = icon + mpp;
 419   isimi = isim + n * n + n;
 420   idatm = isimi + n * n;
 421   ia = idatm + n * mpp + mpp;
 422   ivsig = ia + m * n + n;
 423   iveta = ivsig + n;
 424   isigb = iveta + n;
 425   idx = isigb + n;
 426   iwork = idx + n;
 427   rc = cobylb(&n, &m, &mpp, &x[1], minf, &rhobeg, rhoend, stop, &lb[1], &ub[1], &iprint,
 428       &w[icon], &w[isim], &w[isimi], &w[idatm], &w[ia], &w[ivsig], &w[iveta],
 429       &w[isigb], &w[idx], &w[iwork], &iact[1], calcfc, state);
 430
 431   /* Parameter adjustments (reverse) */
 432   ++iact;
 433   ++w;
 434
 435   free(w);
 436   free(iact);
 437
 438   return rc;
 439 } /* cobyla */
 440
 441 /* ------------------------------------------------------------------------- */
 442 static nlopt_result cobylb(int *n, int *m, int *mpp,
 443    double *x, double *minf, double *rhobeg, double rhoend,
 444     nlopt_stopping *stop, const double *lb, const double *ub,
 445     int *iprint, double *con, double *sim, double *simi,
 446     double *datmat, double *a, double *vsig, double *veta,
 447      double *sigbar, double *dx, double *w, int *iact, cobyla_function *calcfc,
 448      func_wrap_state *state)
 449 {
 450   /* System generated locals */
 451   int sim_dim1, sim_offset, simi_dim1, simi_offset, datmat_dim1,
 452       datmat_offset, a_dim1, a_offset, i__1, i__2, i__3;
 453   double d__1, d__2;
 454
 455   /* Local variables */
 456   double alpha, delta, denom, tempa, barmu;
 457   double beta, cmin = 0.0, cmax = 0.0;
 458   double cvmaxm, dxsign, prerem = 0.0;
 459   double edgmax, pareta, prerec = 0.0, phimin, parsig = 0.0;
 460   double gamma_;
 461   double phi, rho, sum = 0.0;
 462   double ratio, vmold, parmu, error, vmnew;
 463   double resmax, cvmaxp;
 464   double resnew, trured;
 465   double temp, wsig, f;
 466   double weta;
 467   int i__, j, k, l;
 468   int idxnew;
 469   int iflag = 0;
 470   int iptemp;
 471   int isdirn, izdota;
 472   int ivmc;
 473   int ivmd;
 474   int mp, np, iz, ibrnch;
 475   int nbest, ifull, iptem, jdrop;
 476   nlopt_result rc = NLOPT_SUCCESS;
 477   uint32_t seed = (uint32_t) (*n + *m); /* arbitrary deterministic LCG seed */
 478   int feasible;
 479
 480   /* SGJ, 2008: added code to keep track of minimum feasible function val */
 481   *minf = HUGE_VAL;
 482
 483 /* Set the initial values of some parameters. The last column of SIM holds */
 484 /* the optimal vertex of the current simplex, and the preceding N columns */
 485 /* hold the displacements from the optimal vertex to the other vertices. */
 486 /* Further, SIMI holds the inverse of the matrix that is contained in the */
 487 /* first N columns of SIM. */
 488
 489   /* Parameter adjustments */
 490   a_dim1 = *n;
 491   a_offset = 1 + a_dim1 * 1;
 492   a -= a_offset;
 493   simi_dim1 = *n;
 494   simi_offset = 1 + simi_dim1 * 1;
 495   simi -= simi_offset;
 496   sim_dim1 = *n;
 497   sim_offset = 1 + sim_dim1 * 1;
 498   sim -= sim_offset;
 499   datmat_dim1 = *mpp;
 500   datmat_offset = 1 + datmat_dim1 * 1;
 501   datmat -= datmat_offset;
 502   --x;
 503   --con;
 504   --vsig;
 505   --veta;
 506   --sigbar;
 507   --dx;
 508   --w;
 509   --iact;
 510   --lb; --ub;
 511
 512   /* Function Body */
 513   iptem = MIN2(*n,4);
 514   iptemp = iptem + 1;
 515   np = *n + 1;
 516   mp = *m + 1;
 517   alpha = .25;
 518   beta = 2.1;
 519   gamma_ = .5;
 520   delta = 1.1;
 521   rho = *rhobeg;
 522   parmu = 0.;
 523   if (*iprint >= 2) {
 524     fprintf(stderr,
 525       "cobyla: the initial value of RHO is %12.6E and PARMU is set to zero.\n",
 526       rho);
 527   }
 528   temp = 1. / rho;
 529   i__1 = *n;
 530   for (i__ = 1; i__ <= i__1; ++i__) {
 531     double rhocur;
 532     sim[i__ + np * sim_dim1] = x[i__];
 533     i__2 = *n;
 534     for (j = 1; j <= i__2; ++j) {
 535       sim[i__ + j * sim_dim1] = 0.;
 536       simi[i__ + j * simi_dim1] = 0.;
 537     }
 538     rhocur = rho;
 539 #if ENFORCE_BOUNDS
 540     /* SGJ: make sure step rhocur stays inside [lb,ub] */
 541     if (x[i__] + rhocur > ub[i__]) {
 542          if (x[i__] - rhocur >= lb[i__])
 543               rhocur = -rhocur;
 544          else if (ub[i__] - x[i__] > x[i__] - lb[i__])
 545               rhocur = 0.5 * (ub[i__] - x[i__]);
 546          else
 547               rhocur = 0.5 * (x[i__] - lb[i__]);
 548     }
 549 #endif
 550     sim[i__ + i__ * sim_dim1] = rhocur;
 551     simi[i__ + i__ * simi_dim1] = 1.0 / rhocur;
 552   }
 553   jdrop = np;
 554   ibrnch = 0;
 555
 556 /* Make the next call of the user-supplied subroutine CALCFC. These */
 557 /* instructions are also used for calling CALCFC during the iterations of */
 558 /* the algorithm. */
 559
 560   /* SGJ comment: why the hell couldn't he just use a damn subroutine?
 561      #*&!%*@ Fortran-66 spaghetti code */
 562
 563 L40:
 564   if (nlopt_stop_forced(stop)) rc = NLOPT_FORCED_STOP;
 565   else if (nlopt_stop_evals(stop)) rc = NLOPT_MAXEVAL_REACHED;
 566   else if (nlopt_stop_time(stop)) rc = NLOPT_MAXTIME_REACHED;
 567   if (rc != NLOPT_SUCCESS) goto L600;
 568
 569   stop->nevals++;
 570   if (calcfc(*n, *m, &x[1], &f, &con[1], state))
 571   {
 572     if (*iprint >= 1) {
 573       fprintf(stderr, "cobyla: user requested end of minimization.\n");
 574     }
 575     rc = NLOPT_FORCED_STOP;
 576     goto L600;
 577   }
 578
 579   resmax = 0.;
 580   feasible = 1; /* SGJ, 2010 */
 581   if (*m > 0) {
 582     i__1 = *m;
 583     for (k = 1; k <= i__1; ++k) {
 584       d__1 = resmax, d__2 = -con[k];
 585       resmax = MAX2(d__1,d__2);
 586       if (d__2 > state->con_tol[k-1])
 587            feasible = 0; /* SGJ, 2010 */
 588     }
 589   }
 590
 591   /* SGJ, 2008: check for minf_max reached by feasible point */
 592   if (f < stop->minf_max && feasible) {
 593        rc = NLOPT_MINF_MAX_REACHED;
 594        goto L620; /* not L600 because we want to use current x, f, resmax */
 595   }
 596
 597   if (stop->nevals == *iprint - 1 || *iprint == 3) {
 598     fprintf(stderr, "cobyla: NFVALS = %4d, F =%13.6E, MAXCV =%13.6E\n",
 599             stop->nevals, f, resmax);
 600     i__1 = iptem;
 601     fprintf(stderr, "cobyla: X =");
 602     for (i__ = 1; i__ <= i__1; ++i__) {
 603       if (i__>1) fprintf(stderr, "  ");
 604       fprintf(stderr, "%13.6E", x[i__]);
 605     }
 606     if (iptem < *n) {
 607       i__1 = *n;
 608       for (i__ = iptemp; i__ <= i__1; ++i__) {
 609         if (!((i__-1) % 4)) fprintf(stderr, "\ncobyla:  ");
 610         fprintf(stderr, "%15.6E", x[i__]);
 611       }
 612     }
 613     fprintf(stderr, "\n");
 614   }
 615   con[mp] = f;
 616   con[*mpp] = resmax;
 617   if (ibrnch == 1) {
 618     goto L440;
 619   }
 620
 621 /* Set the recently calculated function values in a column of DATMAT. This */
 622 /* array has a column for each vertex of the current simplex, the entries of */
 623 /* each column being the values of the constraint functions (if any) */
 624 /* followed by the objective function and the greatest constraint violation */
 625 /* at the vertex. */
 626
 627   i__1 = *mpp;
 628   for (k = 1; k <= i__1; ++k) {
 629     datmat[k + jdrop * datmat_dim1] = con[k];
 630   }
 631   if (stop->nevals > np) {
 632     goto L130;
 633   }
 634
 635 /* Exchange the new vertex of the initial simplex with the optimal vertex if */
 636 /* necessary. Then, if the initial simplex is not complete, pick its next */
 637 /* vertex and calculate the function values there. */
 638
 639   if (jdrop <= *n) {
 640     if (datmat[mp + np * datmat_dim1] <= f) {
 641       x[jdrop] = sim[jdrop + np * sim_dim1];
 642     } else { /* improvement in function val */
 643       double rhocur = x[jdrop] - sim[jdrop + np * sim_dim1];
 644       /* SGJ: use rhocur instead of rho.  In original code, rhocur == rho
 645               always, but I want to change this to ensure that simplex points
 646               fall within [lb,ub]. */
 647       sim[jdrop + np * sim_dim1] = x[jdrop];
 648       i__1 = *mpp;
 649       for (k = 1; k <= i__1; ++k) {
 650         datmat[k + jdrop * datmat_dim1] = datmat[k + np * datmat_dim1]
 651             ;
 652         datmat[k + np * datmat_dim1] = con[k];
 653       }
 654       i__1 = jdrop;
 655       for (k = 1; k <= i__1; ++k) {
 656         sim[jdrop + k * sim_dim1] = -rhocur;
 657         temp = 0.f;
 658         i__2 = jdrop;
 659         for (i__ = k; i__ <= i__2; ++i__) {
 660           temp -= simi[i__ + k * simi_dim1];
 661         }
 662         simi[jdrop + k * simi_dim1] = temp;
 663       }
 664     }
 665   }
 666   if (stop->nevals <= *n) { /* evaluating initial simplex */
 667     jdrop = stop->nevals;
 668     /* SGJ: was += rho, but using sim[jdrop,jdrop] enforces consistency
 669             if we change the stepsize above to stay in [lb,ub]. */
 670     x[jdrop] += sim[jdrop + jdrop * sim_dim1];
 671     goto L40;
 672   }
 673 L130:
 674   ibrnch = 1;
 675
 676 /* Identify the optimal vertex of the current simplex. */
 677
 678 L140:
 679   phimin = datmat[mp + np * datmat_dim1] + parmu * datmat[*mpp + np *
 680       datmat_dim1];
 681   nbest = np;
 682   i__1 = *n;
 683   for (j = 1; j <= i__1; ++j) {
 684     temp = datmat[mp + j * datmat_dim1] + parmu * datmat[*mpp + j *
 685         datmat_dim1];
 686     if (temp < phimin) {
 687       nbest = j;
 688       phimin = temp;
 689     } else if (temp == phimin && parmu == 0.) {
 690       if (datmat[*mpp + j * datmat_dim1] < datmat[*mpp + nbest *
 691           datmat_dim1]) {
 692         nbest = j;
 693       }
 694     }
 695   }
 696
 697 /* Switch the best vertex into pole position if it is not there already, */
 698 /* and also update SIM, SIMI and DATMAT. */
 699
 700   if (nbest <= *n) {
 701     i__1 = *mpp;
 702     for (i__ = 1; i__ <= i__1; ++i__) {
 703       temp = datmat[i__ + np * datmat_dim1];
 704       datmat[i__ + np * datmat_dim1] = datmat[i__ + nbest * datmat_dim1]
 705           ;
 706       datmat[i__ + nbest * datmat_dim1] = temp;
 707     }
 708     i__1 = *n;
 709     for (i__ = 1; i__ <= i__1; ++i__) {
 710       temp = sim[i__ + nbest * sim_dim1];
 711       sim[i__ + nbest * sim_dim1] = 0.;
 712       sim[i__ + np * sim_dim1] += temp;
 713       tempa = 0.;
 714       i__2 = *n;
 715       for (k = 1; k <= i__2; ++k) {
 716         sim[i__ + k * sim_dim1] -= temp;
 717         tempa -= simi[k + i__ * simi_dim1];
 718       }
 719       simi[nbest + i__ * simi_dim1] = tempa;
 720     }
 721   }
 722
 723 /* Make an error return if SIGI is a poor approximation to the inverse of */
 724 /* the leading N by N submatrix of SIG. */
 725
 726   error = 0.;
 727   i__1 = *n;
 728   for (i__ = 1; i__ <= i__1; ++i__) {
 729     i__2 = *n;
 730     for (j = 1; j <= i__2; ++j) {
 731       temp = 0.;
 732       if (i__ == j) {
 733         temp += -1.;
 734       }
 735       i__3 = *n;
 736       for (k = 1; k <= i__3; ++k) if (sim[k + j * sim_dim1] != 0) {
 737         temp += simi[i__ + k * simi_dim1] * sim[k + j * sim_dim1];
 738       }
 739       d__1 = error, d__2 = fabs(temp);
 740       error = MAX2(d__1,d__2);
 741     }
 742   }
 743   if (error > .1) {
 744     if (*iprint >= 1) {
 745       fprintf(stderr, "cobyla: rounding errors are becoming damaging.\n");
 746     }
 747     rc = NLOPT_FAILURE;
 748     goto L600;
 749   }
 750
 751 /* Calculate the coefficients of the linear approximations to the objective */
 752 /* and constraint functions, placing minus the objective function gradient */
 753 /* after the constraint gradients in the array A. The vector W is used for */
 754 /* working space. */
 755
 756   i__2 = mp;
 757   for (k = 1; k <= i__2; ++k) {
 758     con[k] = -datmat[k + np * datmat_dim1];
 759     i__1 = *n;
 760     for (j = 1; j <= i__1; ++j) {
 761       w[j] = datmat[k + j * datmat_dim1] + con[k];
 762     }
 763     i__1 = *n;
 764     for (i__ = 1; i__ <= i__1; ++i__) {
 765       temp = 0.;
 766       i__3 = *n;
 767       for (j = 1; j <= i__3; ++j) {
 768         temp += w[j] * simi[j + i__ * simi_dim1];
 769       }
 770       if (k == mp) {
 771         temp = -temp;
 772       }
 773       a[i__ + k * a_dim1] = temp;
 774     }
 775   }
 776
 777 /* Calculate the values of sigma and eta, and set IFLAG=0 if the current */
 778 /* simplex is not acceptable. */
 779
 780   iflag = 1;
 781   parsig = alpha * rho;
 782   pareta = beta * rho;
 783   i__1 = *n;
 784   for (j = 1; j <= i__1; ++j) {
 785     wsig = 0.;
 786     weta = 0.;
 787     i__2 = *n;
 788     for (i__ = 1; i__ <= i__2; ++i__) {
 789       d__1 = simi[j + i__ * simi_dim1];
 790       wsig += d__1 * d__1;
 791       d__1 = sim[i__ + j * sim_dim1];
 792       weta += d__1 * d__1;
 793     }
 794     vsig[j] = 1. / sqrt(wsig);
 795     veta[j] = sqrt(weta);
 796     if (vsig[j] < parsig || veta[j] > pareta) {
 797       iflag = 0;
 798     }
 799   }
 800
 801 /* If a new vertex is needed to improve acceptability, then decide which */
 802 /* vertex to drop from the simplex. */
 803
 804   if (ibrnch == 1 || iflag == 1) {
 805     goto L370;
 806   }
 807   jdrop = 0;
 808   temp = pareta;
 809   i__1 = *n;
 810   for (j = 1; j <= i__1; ++j) {
 811     if (veta[j] > temp) {
 812       jdrop = j;
 813       temp = veta[j];
 814     }
 815   }
 816   if (jdrop == 0) {
 817     i__1 = *n;
 818     for (j = 1; j <= i__1; ++j) {
 819       if (vsig[j] < temp) {
 820         jdrop = j;
 821         temp = vsig[j];
 822       }
 823     }
 824   }
 825
 826 /* Calculate the step to the new vertex and its sign. */
 827
 828   temp = gamma_ * rho * vsig[jdrop];
 829   i__1 = *n;
 830   for (i__ = 1; i__ <= i__1; ++i__) {
 831     dx[i__] = temp * simi[jdrop + i__ * simi_dim1];
 832   }
 833   cvmaxp = 0.;
 834   cvmaxm = 0.;
 835   i__1 = mp;
 836   for (k = 1; k <= i__1; ++k) {
 837     sum = 0.;
 838     i__2 = *n;
 839     for (i__ = 1; i__ <= i__2; ++i__) {
 840       sum += a[i__ + k * a_dim1] * dx[i__];
 841     }
 842     if (k < mp) {
 843       temp = datmat[k + np * datmat_dim1];
 844       d__1 = cvmaxp, d__2 = -sum - temp;
 845       cvmaxp = MAX2(d__1,d__2);
 846       d__1 = cvmaxm, d__2 = sum - temp;
 847       cvmaxm = MAX2(d__1,d__2);
 848     }
 849   }
 850   dxsign = 1.;
 851   if (parmu * (cvmaxp - cvmaxm) > sum + sum) {
 852     dxsign = -1.;
 853   }
 854
 855 /* Update the elements of SIM and SIMI, and set the next X. */
 856
 857   temp = 0.;
 858   i__1 = *n;
 859   for (i__ = 1; i__ <= i__1; ++i__) {
 860     /* SGJ, 2010: pseudo-randomize simplex steps (see LCG comments above) */
 861     dx[i__] = dxsign * dx[i__] * lcg_urand(&seed, 0.01, 1);
 862     /* SGJ: make sure dx step says in [lb,ub] */
 863 #if ENFORCE_BOUNDS
 864     {
 865          double xi = sim[i__ + np * sim_dim1];
 866     fixdx:
 867          if (xi + dx[i__] > ub[i__])
 868               dx[i__] = -dx[i__];
 869          if (xi + dx[i__] < lb[i__]) {
 870               if (xi - dx[i__] <= ub[i__])
 871                    dx[i__] = -dx[i__];
 872               else { /* try again with halved step */
 873                    dx[i__] *= 0.5;
 874                    goto fixdx;
 875               }
 876          }
 877     }
 878 #endif
 879     sim[i__ + jdrop * sim_dim1] = dx[i__];
 880     temp += simi[jdrop + i__ * simi_dim1] * dx[i__];
 881   }
 882   i__1 = *n;
 883   for (i__ = 1; i__ <= i__1; ++i__) {
 884     simi[jdrop + i__ * simi_dim1] /= temp;
 885   }
 886   i__1 = *n;
 887   for (j = 1; j <= i__1; ++j) {
 888     if (j != jdrop) {
 889       temp = 0.;
 890       i__2 = *n;
 891       for (i__ = 1; i__ <= i__2; ++i__) {
 892         temp += simi[j + i__ * simi_dim1] * dx[i__];
 893       }
 894       i__2 = *n;
 895       for (i__ = 1; i__ <= i__2; ++i__) {
 896         simi[j + i__ * simi_dim1] -= temp * simi[jdrop + i__ *
 897             simi_dim1];
 898       }
 899     }
 900     x[j] = sim[j + np * sim_dim1] + dx[j];
 901   }
 902   goto L40;
 903
 904 /* Calculate DX=x(*)-x(0). Branch if the length of DX is less than 0.5*RHO. */
 905
 906 L370:
 907   iz = 1;
 908   izdota = iz + *n * *n;
 909   ivmc = izdota + *n;
 910   isdirn = ivmc + mp;
 911   idxnew = isdirn + *n;
 912   ivmd = idxnew + *n;
 913   rc = trstlp(n, m, &a[a_offset], &con[1], &rho, &dx[1], &ifull, &iact[1], &w[
 914       iz], &w[izdota], &w[ivmc], &w[isdirn], &w[idxnew], &w[ivmd]);
 915   if (rc != NLOPT_SUCCESS) goto L600;
 916 #if ENFORCE_BOUNDS
 917   /* SGJ: since the bound constraints are linear, we should never get
 918      a dx that lies outside the [lb,ub] constraints here, but we'll
 919      enforce this anyway just to be paranoid */
 920   i__1 = *n;
 921   for (i__ = 1; i__ <= i__1; ++i__) {
 922        double xi = sim[i__ + np * sim_dim1];
 923        if (xi + dx[i__] > ub[i__]) dx[i__] = ub[i__] - xi;
 924        if (xi + dx[i__] < lb[i__]) dx[i__] = xi - lb[i__];
 925   }
 926 #endif
 927   if (ifull == 0) {
 928     temp = 0.;
 929     i__1 = *n;
 930     for (i__ = 1; i__ <= i__1; ++i__) {
 931       d__1 = dx[i__];
 932       temp += d__1 * d__1;
 933     }
 934     if (temp < rho * .25 * rho) {
 935       ibrnch = 1;
 936       goto L550;
 937     }
 938   }
 939
 940 /* Predict the change to F and the new maximum constraint violation if the */
 941 /* variables are altered from x(0) to x(0)+DX. */
 942
 943   resnew = 0.;
 944   con[mp] = 0.;
 945   i__1 = mp;
 946   for (k = 1; k <= i__1; ++k) {
 947     sum = con[k];
 948     i__2 = *n;
 949     for (i__ = 1; i__ <= i__2; ++i__) {
 950       sum -= a[i__ + k * a_dim1] * dx[i__];
 951     }
 952     if (k < mp) {
 953       resnew = MAX2(resnew,sum);
 954     }
 955   }
 956
 957 /* Increase PARMU if necessary and branch back if this change alters the */
 958 /* optimal vertex. Otherwise PREREM and PREREC will be set to the predicted */
 959 /* reductions in the merit function and the maximum constraint violation */
 960 /* respectively. */
 961
 962   barmu = 0.;
 963   prerec = datmat[*mpp + np * datmat_dim1] - resnew;
 964   if (prerec > 0.) {
 965     barmu = sum / prerec;
 966   }
 967   if (parmu < barmu * 1.5) {
 968     parmu = barmu * 2.;
 969     if (*iprint >= 2) {
 970       fprintf(stderr, "cobyla: increase in PARMU to %12.6E\n", parmu);
 971     }
 972     phi = datmat[mp + np * datmat_dim1] + parmu * datmat[*mpp + np *
 973         datmat_dim1];
 974     i__1 = *n;
 975     for (j = 1; j <= i__1; ++j) {
 976       temp = datmat[mp + j * datmat_dim1] + parmu * datmat[*mpp + j *
 977           datmat_dim1];
 978       if (temp < phi) {
 979         goto L140;
 980       }
 981       if (temp == phi && parmu == 0.f) {
 982         if (datmat[*mpp + j * datmat_dim1] < datmat[*mpp + np *
 983             datmat_dim1]) {
 984           goto L140;
 985         }
 986       }
 987     }
 988   }
 989   prerem = parmu * prerec - sum;
 990
 991 /* Calculate the constraint and objective functions at x(*). Then find the */
 992 /* actual reduction in the merit function. */
 993
 994   i__1 = *n;
 995   for (i__ = 1; i__ <= i__1; ++i__) {
 996     x[i__] = sim[i__ + np * sim_dim1] + dx[i__];
 997   }
 998   ibrnch = 1;
 999   goto L40;
1000 L440:
1001   vmold = datmat[mp + np * datmat_dim1] + parmu * datmat[*mpp + np *
1002       datmat_dim1];
1003   vmnew = f + parmu * resmax;
1004   trured = vmold - vmnew;
1005   if (parmu == 0. && f == datmat[mp + np * datmat_dim1]) {
1006     prerem = prerec;
1007     trured = datmat[*mpp + np * datmat_dim1] - resmax;
1008   }
1009
1010 /* Begin the operations that decide whether x(*) should replace one of the */
1011 /* vertices of the current simplex, the change being mandatory if TRURED is */
1012 /* positive. Firstly, JDROP is set to the index of the vertex that is to be */
1013 /* replaced. */
1014
1015   ratio = 0.;
1016   if (trured <= 0.f) {
1017     ratio = 1.f;
1018   }
1019   jdrop = 0;
1020   i__1 = *n;
1021   for (j = 1; j <= i__1; ++j) {
1022     temp = 0.;
1023     i__2 = *n;
1024     for (i__ = 1; i__ <= i__2; ++i__) {
1025       temp += simi[j + i__ * simi_dim1] * dx[i__];
1026     }
1027     temp = fabs(temp);
1028     if (temp > ratio) {
1029       jdrop = j;
1030       ratio = temp;
1031     }
1032     sigbar[j] = temp * vsig[j];
1033   }
1034
1035 /* Calculate the value of ell. */
1036
1037   edgmax = delta * rho;
1038   l = 0;
1039   i__1 = *n;
1040   for (j = 1; j <= i__1; ++j) {
1041     if (sigbar[j] >= parsig || sigbar[j] >= vsig[j]) {
1042       temp = veta[j];
1043       if (trured > 0.) {
1044         temp = 0.;
1045         i__2 = *n;
1046         for (i__ = 1; i__ <= i__2; ++i__) {
1047           d__1 = dx[i__] - sim[i__ + j * sim_dim1];
1048           temp += d__1 * d__1;
1049         }
1050         temp = sqrt(temp);
1051       }
1052       if (temp > edgmax) {
1053         l = j;
1054         edgmax = temp;
1055       }
1056     }
1057   }
1058   if (l > 0) {
1059     jdrop = l;
1060   }
1061   if (jdrop == 0) {
1062     goto L550;
1063   }
1064
1065 /* Revise the simplex by updating the elements of SIM, SIMI and DATMAT. */
1066
1067   temp = 0.;
1068   i__1 = *n;
1069   for (i__ = 1; i__ <= i__1; ++i__) {
1070     sim[i__ + jdrop * sim_dim1] = dx[i__];
1071     temp += simi[jdrop + i__ * simi_dim1] * dx[i__];
1072   }
1073   i__1 = *n;
1074   for (i__ = 1; i__ <= i__1; ++i__) {
1075     simi[jdrop + i__ * simi_dim1] /= temp;
1076   }
1077   i__1 = *n;
1078   for (j = 1; j <= i__1; ++j) {
1079     if (j != jdrop) {
1080       temp = 0.;
1081       i__2 = *n;
1082       for (i__ = 1; i__ <= i__2; ++i__) {
1083         temp += simi[j + i__ * simi_dim1] * dx[i__];
1084       }
1085       i__2 = *n;
1086       for (i__ = 1; i__ <= i__2; ++i__) {
1087         simi[j + i__ * simi_dim1] -= temp * simi[jdrop + i__ *
1088             simi_dim1];
1089       }
1090     }
1091   }
1092   i__1 = *mpp;
1093   for (k = 1; k <= i__1; ++k) {
1094     datmat[k + jdrop * datmat_dim1] = con[k];
1095   }
1096
1097 /* Branch back for further iterations with the current RHO. */
1098
1099   if (trured > 0. && trured >= prerem * .1) {
1100     /* SGJ, 2010: following a suggestion in the SAS manual (which
1101        mentions a similar modification to COBYLA, although they didn't
1102        publish their source code), increase rho if predicted reduction
1103        is sufficiently close to the actual reduction.  Otherwise,
1104        COBLYA seems to easily get stuck making very small steps.
1105        Also require iflag != 0 (i.e., acceptable simplex), again
1106        following SAS suggestion (otherwise I get convergence failure
1107        in some cases.) */
1108     if (trured >= prerem * 0.9 && trured <= prerem * 1.1 && iflag) {
1109          rho *= 2.0;
1110     }
1111     goto L140;
1112   }
1113 L550:
1114   if (iflag == 0) {
1115     ibrnch = 0;
1116     goto L140;
1117   }
1118
1119   /* SGJ, 2008: convergence tests for function vals; note that current
1120      best val is stored in datmat[mp + np * datmat_dim1], or in f if
1121      ifull == 1, and previous best is in *minf.  This seems like a
1122      sensible place to put the convergence test, as it is the same
1123      place that Powell checks the x tolerance (rho > rhoend). */
1124   {
1125        double fbest = ifull == 1 ? f : datmat[mp + np * datmat_dim1];
1126        if (fbest < *minf && nlopt_stop_ftol(stop, fbest, *minf)) {
1127             rc = NLOPT_FTOL_REACHED;
1128             goto L600;
1129        }
1130        *minf = fbest;
1131   }
1132
1133 /* Otherwise reduce RHO if it is not at its least value and reset PARMU. */
1134
1135   if (rho > rhoend) {
1136     rho *= .5;
1137     if (rho <= rhoend * 1.5) {
1138       rho = rhoend;
1139     }
1140     if (parmu > 0.) {
1141       denom = 0.;
1142       i__1 = mp;
1143       for (k = 1; k <= i__1; ++k) {
1144         cmin = datmat[k + np * datmat_dim1];
1145         cmax = cmin;
1146         i__2 = *n;
1147         for (i__ = 1; i__ <= i__2; ++i__) {
1148           d__1 = cmin, d__2 = datmat[k + i__ * datmat_dim1];
1149           cmin = MIN2(d__1,d__2);
1150           d__1 = cmax, d__2 = datmat[k + i__ * datmat_dim1];
1151           cmax = MAX2(d__1,d__2);
1152         }
1153         if (k <= *m && cmin < cmax * .5) {
1154           temp = MAX2(cmax,0.) - cmin;
1155           if (denom <= 0.) {
1156             denom = temp;
1157           } else {
1158             denom = MIN2(denom,temp);
1159           }
1160         }
1161       }
1162       if (denom == 0.) {
1163         parmu = 0.;
1164       } else if (cmax - cmin < parmu * denom) {
1165         parmu = (cmax - cmin) / denom;
1166       }
1167     }
1168     if (*iprint >= 2) {
1169       fprintf(stderr, "cobyla: reduction in RHO to %12.6E and PARMU =%13.6E\n",
1170         rho, parmu);
1171     }
1172     if (*iprint == 2) {
1173       fprintf(stderr, "cobyla: NFVALS = %4d, F =%13.6E, MAXCV =%13.6E\n",
1174         stop->nevals, datmat[mp + np * datmat_dim1], datmat[*mpp + np * datmat_dim1]);
1175
1176       fprintf(stderr, "cobyla: X =");
1177       i__1 = iptem;
1178       for (i__ = 1; i__ <= i__1; ++i__) {
1179         if (i__>1) fprintf(stderr, "  ");
1180         fprintf(stderr, "%13.6E", sim[i__ + np * sim_dim1]);
1181       }
1182       if (iptem < *n) {
1183         i__1 = *n;
1184         for (i__ = iptemp; i__ <= i__1; ++i__) {
1185           if (!((i__-1) % 4)) fprintf(stderr, "\ncobyla:  ");
1186           fprintf(stderr, "%15.6E", x[i__]);
1187         }
1188       }
1189       fprintf(stderr, "\n");
1190     }
1191     goto L140;
1192   }
1193   else
1194        rc = NLOPT_XTOL_REACHED;
1195
1196 /* Return the best calculated values of the variables. */
1197
1198   if (*iprint >= 1) {
1199     fprintf(stderr, "cobyla: normal return.\n");
1200   }
1201   if (ifull == 1) {
1202     goto L620;
1203   }
1204 L600:
1205   i__1 = *n;
1206   for (i__ = 1; i__ <= i__1; ++i__) {
1207     x[i__] = sim[i__ + np * sim_dim1];
1208   }
1209   f = datmat[mp + np * datmat_dim1];
1210   resmax = datmat[*mpp + np * datmat_dim1];
1211 L620:
1212   *minf = f;
1213   if (*iprint >= 1) {
1214     fprintf(stderr, "cobyla: NFVALS = %4d, F =%13.6E, MAXCV =%13.6E\n",
1215             stop->nevals, f, resmax);
1216     i__1 = iptem;
1217     fprintf(stderr, "cobyla: X =");
1218     for (i__ = 1; i__ <= i__1; ++i__) {
1219       if (i__>1) fprintf(stderr, "  ");
1220       fprintf(stderr, "%13.6E", x[i__]);
1221     }
1222     if (iptem < *n) {
1223       i__1 = *n;
1224       for (i__ = iptemp; i__ <= i__1; ++i__) {
1225         if (!((i__-1) % 4)) fprintf(stderr, "\ncobyla:  ");
1226         fprintf(stderr, "%15.6E", x[i__]);
1227       }
1228     }
1229     fprintf(stderr, "\n");
1230   }
1231   return rc;
1232 } /* cobylb */
1233
1234 /* ------------------------------------------------------------------------- */
1235 static nlopt_result trstlp(int *n, int *m, double *a,
1236     double *b, double *rho, double *dx, int *ifull,
1237     int *iact, double *z__, double *zdota, double *vmultc,
1238      double *sdirn, double *dxnew, double *vmultd)
1239 {
1240   /* System generated locals */
1241   int a_dim1, a_offset, z_dim1, z_offset, i__1, i__2;
1242   double d__1, d__2;
1243
1244   /* Local variables */
1245   double alpha, tempa;
1246   double beta;
1247   double optnew, stpful, sum, tot, acca, accb;
1248   double ratio, vsave, zdotv, zdotw, dd;
1249   double sd;
1250   double sp, ss, resold = 0.0, zdvabs, zdwabs, sumabs, resmax, optold;
1251   double spabs;
1252   double temp, step;
1253   int icount;
1254   int iout, i__, j, k;
1255   int isave;
1256   int kk;
1257   int kl, kp, kw;
1258   int nact, icon = 0, mcon;
1259   int nactx = 0;
1260
1261
1262 /* This subroutine calculates an N-component vector DX by applying the */
1263 /* following two stages. In the first stage, DX is set to the shortest */
1264 /* vector that minimizes the greatest violation of the constraints */
1265 /*   A(1,K)*DX(1)+A(2,K)*DX(2)+...+A(N,K)*DX(N) .GE. B(K), K=2,3,...,M, */
1266 /* subject to the Euclidean length of DX being at most RHO. If its length is */
1267 /* strictly less than RHO, then we use the resultant freedom in DX to */
1268 /* minimize the objective function */
1269 /*      -A(1,M+1)*DX(1)-A(2,M+1)*DX(2)-...-A(N,M+1)*DX(N) */
1270 /* subject to no increase in any greatest constraint violation. This */
1271 /* notation allows the gradient of the objective function to be regarded as */
1272 /* the gradient of a constraint. Therefore the two stages are distinguished */
1273 /* by MCON .EQ. M and MCON .GT. M respectively. It is possible that a */
1274 /* degeneracy may prevent DX from attaining the target length RHO. Then the */
1275 /* value IFULL=0 would be set, but usually IFULL=1 on return. */
1276
1277 /* In general NACT is the number of constraints in the active set and */
1278 /* IACT(1),...,IACT(NACT) are their indices, while the remainder of IACT */
1279 /* contains a permutation of the remaining constraint indices. Further, Z is */
1280 /* an orthogonal matrix whose first NACT columns can be regarded as the */
1281 /* result of Gram-Schmidt applied to the active constraint gradients. For */
1282 /* J=1,2,...,NACT, the number ZDOTA(J) is the scalar product of the J-th */
1283 /* column of Z with the gradient of the J-th active constraint. DX is the */
1284 /* current vector of variables and here the residuals of the active */
1285 /* constraints should be zero. Further, the active constraints have */
1286 /* nonnegative Lagrange multipliers that are held at the beginning of */
1287 /* VMULTC. The remainder of this vector holds the residuals of the inactive */
1288 /* constraints at DX, the ordering of the components of VMULTC being in */
1289 /* agreement with the permutation of the indices of the constraints that is */
1290 /* in IACT. All these residuals are nonnegative, which is achieved by the */
1291 /* shift RESMAX that makes the least residual zero. */
1292
1293 /* Initialize Z and some other variables. The value of RESMAX will be */
1294 /* appropriate to DX=0, while ICON will be the index of a most violated */
1295 /* constraint if RESMAX is positive. Usually during the first stage the */
1296 /* vector SDIRN gives a search direction that reduces all the active */
1297 /* constraint violations by one simultaneously. */
1298
1299   /* Parameter adjustments */
1300   z_dim1 = *n;
1301   z_offset = 1 + z_dim1 * 1;
1302   z__ -= z_offset;
1303   a_dim1 = *n;
1304   a_offset = 1 + a_dim1 * 1;
1305   a -= a_offset;
1306   --b;
1307   --dx;
1308   --iact;
1309   --zdota;
1310   --vmultc;
1311   --sdirn;
1312   --dxnew;
1313   --vmultd;
1314
1315   /* Function Body */
1316   *ifull = 1;
1317   mcon = *m;
1318   nact = 0;
1319   resmax = 0.;
1320   i__1 = *n;
1321   for (i__ = 1; i__ <= i__1; ++i__) {
1322     i__2 = *n;
1323     for (j = 1; j <= i__2; ++j) {
1324       z__[i__ + j * z_dim1] = 0.;
1325     }
1326     z__[i__ + i__ * z_dim1] = 1.;
1327     dx[i__] = 0.;
1328   }
1329   if (*m >= 1) {
1330     i__1 = *m;
1331     for (k = 1; k <= i__1; ++k) {
1332       if (b[k] > resmax) {
1333         resmax = b[k];
1334         icon = k;
1335       }
1336     }
1337     i__1 = *m;
1338     for (k = 1; k <= i__1; ++k) {
1339       iact[k] = k;
1340       vmultc[k] = resmax - b[k];
1341     }
1342   }
1343   if (resmax == 0.) {
1344     goto L480;
1345   }
1346   i__1 = *n;
1347   for (i__ = 1; i__ <= i__1; ++i__) {
1348     sdirn[i__] = 0.;
1349   }
1350
1351 /* End the current stage of the calculation if 3 consecutive iterations */
1352 /* have either failed to reduce the best calculated value of the objective */
1353 /* function or to increase the number of active constraints since the best */
1354 /* value was calculated. This strategy prevents cycling, but there is a */
1355 /* remote possibility that it will cause premature termination. */
1356
1357 L60:
1358   optold = 0.;
1359   icount = 0;
1360 L70:
1361   if (mcon == *m) {
1362     optnew = resmax;
1363   } else {
1364     optnew = 0.;
1365     i__1 = *n;
1366     for (i__ = 1; i__ <= i__1; ++i__) {
1367       optnew -= dx[i__] * a[i__ + mcon * a_dim1];
1368     }
1369   }
1370   if (icount == 0 || optnew < optold) {
1371     optold = optnew;
1372     nactx = nact;
1373     icount = 3;
1374   } else if (nact > nactx) {
1375     nactx = nact;
1376     icount = 3;
1377   } else {
1378     --icount;
1379     if (icount == 0) {
1380       goto L490;
1381     }
1382   }
1383
1384 /* If ICON exceeds NACT, then we add the constraint with index IACT(ICON) to */
1385 /* the active set. Apply Givens rotations so that the last N-NACT-1 columns */
1386 /* of Z are orthogonal to the gradient of the new constraint, a scalar */
1387 /* product being set to zero if its nonzero value could be due to computer */
1388 /* rounding errors. The array DXNEW is used for working space. */
1389
1390   if (icon <= nact) {
1391     goto L260;
1392   }
1393   kk = iact[icon];
1394   i__1 = *n;
1395   for (i__ = 1; i__ <= i__1; ++i__) {
1396     dxnew[i__] = a[i__ + kk * a_dim1];
1397   }
1398   tot = 0.;
1399   k = *n;
1400 L100:
1401   if (k > nact) {
1402     sp = 0.;
1403     spabs = 0.;
1404     i__1 = *n;
1405     for (i__ = 1; i__ <= i__1; ++i__) {
1406       temp = z__[i__ + k * z_dim1] * dxnew[i__];
1407       sp += temp;
1408       spabs += fabs(temp);
1409     }
1410     acca = spabs + fabs(sp) * .1;
1411     accb = spabs + fabs(sp) * .2;
1412     if (spabs >= acca || acca >= accb) {
1413       sp = 0.;
1414     }
1415     if (tot == 0.) {
1416       tot = sp;
1417     } else {
1418       kp = k + 1;
1419       temp = sqrt(sp * sp + tot * tot);
1420       alpha = sp / temp;
1421       beta = tot / temp;
1422       tot = temp;
1423       i__1 = *n;
1424       for (i__ = 1; i__ <= i__1; ++i__) {
1425         temp = alpha * z__[i__ + k * z_dim1] + beta * z__[i__ + kp *
1426             z_dim1];
1427         z__[i__ + kp * z_dim1] = alpha * z__[i__ + kp * z_dim1] -
1428             beta * z__[i__ + k * z_dim1];
1429         z__[i__ + k * z_dim1] = temp;
1430       }
1431     }
1432     --k;
1433     goto L100;
1434   }
1435
1436 /* Add the new constraint if this can be done without a deletion from the */
1437 /* active set. */
1438
1439   if (tot != 0.) {
1440     ++nact;
1441     zdota[nact] = tot;
1442     vmultc[icon] = vmultc[nact];
1443     vmultc[nact] = 0.;
1444     goto L210;
1445   }
1446
1447 /* The next instruction is reached if a deletion has to be made from the */
1448 /* active set in order to make room for the new active constraint, because */
1449 /* the new constraint gradient is a linear combination of the gradients of */
1450 /* the old active constraints. Set the elements of VMULTD to the multipliers */
1451 /* of the linear combination. Further, set IOUT to the index of the */
1452 /* constraint to be deleted, but branch if no suitable index can be found. */
1453
1454   ratio = -1.;
1455   k = nact;
1456 L130:
1457   zdotv = 0.;
1458   zdvabs = 0.;
1459   i__1 = *n;
1460   for (i__ = 1; i__ <= i__1; ++i__) {
1461     temp = z__[i__ + k * z_dim1] * dxnew[i__];
1462     zdotv += temp;
1463     zdvabs += fabs(temp);
1464   }
1465   acca = zdvabs + fabs(zdotv) * .1;
1466   accb = zdvabs + fabs(zdotv) * .2;
1467   if (zdvabs < acca && acca < accb) {
1468     temp = zdotv / zdota[k];
1469     if (temp > 0. && iact[k] <= *m) {
1470       tempa = vmultc[k] / temp;
1471       if (ratio < 0. || tempa < ratio) {
1472         ratio = tempa;
1473         iout = k;
1474       }
1475     }
1476     if (k >= 2) {
1477       kw = iact[k];
1478       i__1 = *n;
1479       for (i__ = 1; i__ <= i__1; ++i__) {
1480         dxnew[i__] -= temp * a[i__ + kw * a_dim1];
1481       }
1482     }
1483     vmultd[k] = temp;
1484   } else {
1485     vmultd[k] = 0.;
1486   }
1487   --k;
1488   if (k > 0) {
1489     goto L130;
1490   }
1491   if (ratio < 0.) {
1492     goto L490;
1493   }
1494
1495 /* Revise the Lagrange multipliers and reorder the active constraints so */
1496 /* that the one to be replaced is at the end of the list. Also calculate the */
1497 /* new value of ZDOTA(NACT) and branch if it is not acceptable. */
1498
1499   i__1 = nact;
1500   for (k = 1; k <= i__1; ++k) {
1501     d__1 = 0., d__2 = vmultc[k] - ratio * vmultd[k];
1502     vmultc[k] = MAX2(d__1,d__2);
1503   }
1504   if (icon < nact) {
1505     isave = iact[icon];
1506     vsave = vmultc[icon];
1507     k = icon;
1508 L170:
1509     kp = k + 1;
1510     kw = iact[kp];
1511     sp = 0.;
1512     i__1 = *n;
1513     for (i__ = 1; i__ <= i__1; ++i__) {
1514       sp += z__[i__ + k * z_dim1] * a[i__ + kw * a_dim1];
1515     }
1516     d__1 = zdota[kp];
1517     temp = sqrt(sp * sp + d__1 * d__1);
1518     alpha = zdota[kp] / temp;
1519     beta = sp / temp;
1520     zdota[kp] = alpha * zdota[k];
1521     zdota[k] = temp;
1522     i__1 = *n;
1523     for (i__ = 1; i__ <= i__1; ++i__) {
1524       temp = alpha * z__[i__ + kp * z_dim1] + beta * z__[i__ + k *
1525           z_dim1];
1526       z__[i__ + kp * z_dim1] = alpha * z__[i__ + k * z_dim1] - beta *
1527           z__[i__ + kp * z_dim1];
1528       z__[i__ + k * z_dim1] = temp;
1529     }
1530     iact[k] = kw;
1531     vmultc[k] = vmultc[kp];
1532     k = kp;
1533     if (k < nact) {
1534       goto L170;
1535     }
1536     iact[k] = isave;
1537     vmultc[k] = vsave;
1538   }
1539   temp = 0.;
1540   i__1 = *n;
1541   for (i__ = 1; i__ <= i__1; ++i__) {
1542     temp += z__[i__ + nact * z_dim1] * a[i__ + kk * a_dim1];
1543   }
1544   if (temp == 0.) {
1545     goto L490;
1546   }
1547   zdota[nact] = temp;
1548   vmultc[icon] = 0.;
1549   vmultc[nact] = ratio;
1550
1551 /* Update IACT and ensure that the objective function continues to be */
1552 /* treated as the last active constraint when MCON>M. */
1553
1554 L210:
1555   iact[icon] = iact[nact];
1556   iact[nact] = kk;
1557   if (mcon > *m && kk != mcon) {
1558     k = nact - 1;
1559     sp = 0.;
1560     i__1 = *n;
1561     for (i__ = 1; i__ <= i__1; ++i__) {
1562       sp += z__[i__ + k * z_dim1] * a[i__ + kk * a_dim1];
1563     }
1564     d__1 = zdota[nact];
1565     temp = sqrt(sp * sp + d__1 * d__1);
1566     alpha = zdota[nact] / temp;
1567     beta = sp / temp;
1568     zdota[nact] = alpha * zdota[k];
1569     zdota[k] = temp;
1570     i__1 = *n;
1571     for (i__ = 1; i__ <= i__1; ++i__) {
1572       temp = alpha * z__[i__ + nact * z_dim1] + beta * z__[i__ + k *
1573           z_dim1];
1574       z__[i__ + nact * z_dim1] = alpha * z__[i__ + k * z_dim1] - beta *
1575           z__[i__ + nact * z_dim1];
1576       z__[i__ + k * z_dim1] = temp;
1577     }
1578     iact[nact] = iact[k];
1579     iact[k] = kk;
1580     temp = vmultc[k];
1581     vmultc[k] = vmultc[nact];
1582     vmultc[nact] = temp;
1583   }
1584
1585 /* If stage one is in progress, then set SDIRN to the direction of the next */
1586 /* change to the current vector of variables. */
1587
1588   if (mcon > *m) {
1589     goto L320;
1590   }
1591   kk = iact[nact];
1592   temp = 0.;
1593   i__1 = *n;
1594   for (i__ = 1; i__ <= i__1; ++i__) {
1595     temp += sdirn[i__] * a[i__ + kk * a_dim1];
1596   }
1597   temp += -1.;
1598   temp /= zdota[nact];
1599   i__1 = *n;
1600   for (i__ = 1; i__ <= i__1; ++i__) {
1601     sdirn[i__] -= temp * z__[i__ + nact * z_dim1];
1602   }
1603   goto L340;
1604
1605 /* Delete the constraint that has the index IACT(ICON) from the active set. */
1606
1607 L260:
1608   if (icon < nact) {
1609     isave = iact[icon];
1610     vsave = vmultc[icon];
1611     k = icon;
1612 L270:
1613     kp = k + 1;
1614     kk = iact[kp];
1615     sp = 0.;
1616     i__1 = *n;
1617     for (i__ = 1; i__ <= i__1; ++i__) {
1618       sp += z__[i__ + k * z_dim1] * a[i__ + kk * a_dim1];
1619     }
1620     d__1 = zdota[kp];
1621     temp = sqrt(sp * sp + d__1 * d__1);
1622     alpha = zdota[kp] / temp;
1623     beta = sp / temp;
1624     zdota[kp] = alpha * zdota[k];
1625     zdota[k] = temp;
1626     i__1 = *n;
1627     for (i__ = 1; i__ <= i__1; ++i__) {
1628       temp = alpha * z__[i__ + kp * z_dim1] + beta * z__[i__ + k *
1629           z_dim1];
1630       z__[i__ + kp * z_dim1] = alpha * z__[i__ + k * z_dim1] - beta *
1631           z__[i__ + kp * z_dim1];
1632       z__[i__ + k * z_dim1] = temp;
1633     }
1634     iact[k] = kk;
1635     vmultc[k] = vmultc[kp];
1636     k = kp;
1637     if (k < nact) {
1638       goto L270;
1639     }
1640     iact[k] = isave;
1641     vmultc[k] = vsave;
1642   }
1643   --nact;
1644
1645 /* If stage one is in progress, then set SDIRN to the direction of the next */
1646 /* change to the current vector of variables. */
1647
1648   if (mcon > *m) {
1649     goto L320;
1650   }
1651   temp = 0.;
1652   i__1 = *n;
1653   for (i__ = 1; i__ <= i__1; ++i__) {
1654     temp += sdirn[i__] * z__[i__ + (nact + 1) * z_dim1];
1655   }
1656   i__1 = *n;
1657   for (i__ = 1; i__ <= i__1; ++i__) {
1658     sdirn[i__] -= temp * z__[i__ + (nact + 1) * z_dim1];
1659   }
1660   goto L340;
1661
1662 /* Pick the next search direction of stage two. */
1663
1664 L320:
1665   temp = 1. / zdota[nact];
1666   i__1 = *n;
1667   for (i__ = 1; i__ <= i__1; ++i__) {
1668     sdirn[i__] = temp * z__[i__ + nact * z_dim1];
1669   }
1670
1671 /* Calculate the step to the boundary of the trust region or take the step */
1672 /* that reduces RESMAX to zero. The two statements below that include the */
1673 /* factor 1.0E-6 prevent some harmless underflows that occurred in a test */
1674 /* calculation. Further, we skip the step if it could be zero within a */
1675 /* reasonable tolerance for computer rounding errors. */
1676
1677 L340:
1678   dd = *rho * *rho;
1679   sd = 0.;
1680   ss = 0.;
1681   i__1 = *n;
1682   for (i__ = 1; i__ <= i__1; ++i__) {
1683     if ((d__1 = dx[i__], fabs(d__1)) >= *rho * 1e-6f) {
1684       d__2 = dx[i__];
1685       dd -= d__2 * d__2;
1686     }
1687     sd += dx[i__] * sdirn[i__];
1688     d__1 = sdirn[i__];
1689     ss += d__1 * d__1;
1690   }
1691   if (dd <= 0.) {
1692     goto L490;
1693   }
1694   temp = sqrt(ss * dd);
1695   if (fabs(sd) >= temp * 1e-6f) {
1696     temp = sqrt(ss * dd + sd * sd);
1697   }
1698   stpful = dd / (temp + sd);
1699   step = stpful;
1700   if (mcon == *m) {
1701     acca = step + resmax * .1;
1702     accb = step + resmax * .2;
1703     if (step >= acca || acca >= accb) {
1704       goto L480;
1705     }
1706     step = MIN2(step,resmax);
1707   }
1708
1709   /* SGJ, 2010: check for error here */
1710   if (nlopt_isinf(step)) return NLOPT_ROUNDOFF_LIMITED;
1711
1712 /* Set DXNEW to the new variables if STEP is the steplength, and reduce */
1713 /* RESMAX to the corresponding maximum residual if stage one is being done. */
1714 /* Because DXNEW will be changed during the calculation of some Lagrange */
1715 /* multipliers, it will be restored to the following value later. */
1716
1717   i__1 = *n;
1718   for (i__ = 1; i__ <= i__1; ++i__) {
1719     dxnew[i__] = dx[i__] + step * sdirn[i__];
1720   }
1721   if (mcon == *m) {
1722     resold = resmax;
1723     resmax = 0.;
1724     i__1 = nact;
1725     for (k = 1; k <= i__1; ++k) {
1726       kk = iact[k];
1727       temp = b[kk];
1728       i__2 = *n;
1729       for (i__ = 1; i__ <= i__2; ++i__) {
1730         temp -= a[i__ + kk * a_dim1] * dxnew[i__];
1731       }
1732       resmax = MAX2(resmax,temp);
1733     }
1734   }
1735
1736 /* Set VMULTD to the VMULTC vector that would occur if DX became DXNEW. A */
1737 /* device is included to force VMULTD(K)=0.0 if deviations from this value */
1738 /* can be attributed to computer rounding errors. First calculate the new */
1739 /* Lagrange multipliers. */
1740
1741   k = nact;
1742 L390:
1743   zdotw = 0.;
1744   zdwabs = 0.;
1745   i__1 = *n;
1746   for (i__ = 1; i__ <= i__1; ++i__) {
1747     temp = z__[i__ + k * z_dim1] * dxnew[i__];
1748     zdotw += temp;
1749     zdwabs += fabs(temp);
1750   }
1751   acca = zdwabs + fabs(zdotw) * .1;
1752   accb = zdwabs + fabs(zdotw) * .2;
1753   if (zdwabs >= acca || acca >= accb) {
1754     zdotw = 0.;
1755   }
1756   vmultd[k] = zdotw / zdota[k];
1757   if (k >= 2) {
1758     kk = iact[k];
1759     i__1 = *n;
1760     for (i__ = 1; i__ <= i__1; ++i__) {
1761       dxnew[i__] -= vmultd[k] * a[i__ + kk * a_dim1];
1762     }
1763     --k;
1764     goto L390;
1765   }
1766   if (mcon > *m) {
1767     d__1 = 0., d__2 = vmultd[nact];
1768     vmultd[nact] = MAX2(d__1,d__2);
1769   }
1770
1771 /* Complete VMULTC by finding the new constraint residuals. */
1772
1773   i__1 = *n;
1774   for (i__ = 1; i__ <= i__1; ++i__) {
1775     dxnew[i__] = dx[i__] + step * sdirn[i__];
1776   }
1777   if (mcon > nact) {
1778     kl = nact + 1;
1779     i__1 = mcon;
1780     for (k = kl; k <= i__1; ++k) {
1781       kk = iact[k];
1782       sum = resmax - b[kk];
1783       sumabs = resmax + (d__1 = b[kk], fabs(d__1));
1784       i__2 = *n;
1785       for (i__ = 1; i__ <= i__2; ++i__) {
1786         temp = a[i__ + kk * a_dim1] * dxnew[i__];
1787         sum += temp;
1788         sumabs += fabs(temp);
1789       }
1790       acca = sumabs + fabs(sum) * .1f;
1791       accb = sumabs + fabs(sum) * .2f;
1792       if (sumabs >= acca || acca >= accb) {
1793         sum = 0.f;
1794       }
1795       vmultd[k] = sum;
1796     }
1797   }
1798
1799 /* Calculate the fraction of the step from DX to DXNEW that will be taken. */
1800
1801   ratio = 1.;
1802   icon = 0;
1803   i__1 = mcon;
1804   for (k = 1; k <= i__1; ++k) {
1805     if (vmultd[k] < 0.) {
1806       temp = vmultc[k] / (vmultc[k] - vmultd[k]);
1807       if (temp < ratio) {
1808         ratio = temp;
1809         icon = k;
1810       }
1811     }
1812   }
1813
1814 /* Update DX, VMULTC and RESMAX. */
1815
1816   temp = 1. - ratio;
1817   i__1 = *n;
1818   for (i__ = 1; i__ <= i__1; ++i__) {
1819     dx[i__] = temp * dx[i__] + ratio * dxnew[i__];
1820   }
1821   i__1 = mcon;
1822   for (k = 1; k <= i__1; ++k) {
1823     d__1 = 0., d__2 = temp * vmultc[k] + ratio * vmultd[k];
1824     vmultc[k] = MAX2(d__1,d__2);
1825   }
1826   if (mcon == *m) {
1827     resmax = resold + ratio * (resmax - resold);
1828   }
1829
1830 /* If the full step is not acceptable then begin another iteration. */
1831 /* Otherwise switch to stage two or end the calculation. */
1832
1833   if (icon > 0) {
1834     goto L70;
1835   }
1836   if (step == stpful) {
1837     goto L500;
1838   }
1839 L480:
1840   mcon = *m + 1;
1841   icon = mcon;
1842   iact[mcon] = mcon;
1843   vmultc[mcon] = 0.;
1844   goto L60;
1845
1846 /* We employ any freedom that may be available to reduce the objective */
1847 /* function before returning a DX whose length is less than RHO. */
1848
1849 L490:
1850   if (mcon == *m) {
1851     goto L480;
1852   }
1853   *ifull = 0;
1854 L500:
1855   return NLOPT_SUCCESS;
1856 } /* trstlp */