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