src/algs/mma/ccsa_quadratic.c

   1 /* Copyright (c) 2007-2014 Massachusetts Institute of Technology
   2  *
   3  * Permission is hereby granted, free of charge, to any person obtaining
   4  * a copy of this software and associated documentation files (the
   5  * "Software"), to deal in the Software without restriction, including
   6  * without limitation the rights to use, copy, modify, merge, publish,
   7  * distribute, sublicense, and/or sell copies of the Software, and to
   8  * permit persons to whom the Software is furnished to do so, subject to
   9  * the following conditions:
  10  *
  11  * The above copyright notice and this permission notice shall be
  12  * included in all copies or substantial portions of the Software.
  13  *
  14  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  15  * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  16  * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  17  * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
  18  * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
  19  * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
  20  * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  21  */
  22
  23 /* In this file we implement Svanberg's CCSA algorithm with the
  24    simple linear approximation + quadratic penalty function.
  25
  26    We also allow the user to specify an optional "preconditioner": an
  27    approximate Hessian (which must be symmetric & positive
  28    semidefinite) that can be added into the approximation.  [X. Liang
  29    and I went through the convergence proof in Svanberg's paper
  30    and it does not seem to be modified by such preconditioning, as
  31    long as the preconditioner eigenvalues are bounded above for all x.]
  32
  33    For the non-preconditioned case the trust-region subproblem is
  34    separable and can be solved by a dual method.  For the preconditioned
  35    case the subproblem is still convex but in general is non-separable
  36    so we solve by calling the same algorithm recursively, under the
  37    assumption that the subproblem objective is cheap to evaluate.
  38 */
  39
  40 #include <stdlib.h>
  41 #include <math.h>
  42 #include <string.h>
  43 #include <stdio.h>
  44
  45 #include "mma.h"
  46 #include "nlopt-util.h"
  47
  48 unsigned ccsa_verbose = 0; /* > 0 for verbose output */
  49
  50 #define MIN(a,b) ((a) < (b) ? (a) : (b))
  51 #define MAX(a,b) ((a) > (b) ? (a) : (b))
  52
  53 /* magic minimum value for rho in CCSA ... the 2002 paper says it should
  54    be a "fixed, strictly positive `small' number, e.g. 1e-5"
  55    ... grrr, I hate these magic numbers, which seem like they
  56    should depend on the objective function in some way ... in particular,
  57    note that rho is dimensionful (= dimensions of objective function) */
  58 #define CCSA_RHOMIN 1e-5
  59
  60 /***********************************************************************/
  61 /* function for CCSA's dual solution of the approximate problem */
  62
  63 typedef struct {
  64      int count; /* evaluation count, incremented each call */
  65      unsigned n; /* must be set on input to dimension of x */
  66      const double *x, *lb, *ub, *sigma, *dfdx; /* arrays of length n */
  67      const double *dfcdx; /* m-by-n array of fc gradients */
  68      double fval, rho; /* must be set on input */
  69      const double *fcval, *rhoc; /* arrays of length m */
  70      double *xcur; /* array of length n, output each time */
  71      double gval, wval, *gcval; /* output each time (array length m) */
  72      nlopt_precond pre; void *pre_data;
  73      nlopt_precond *prec; void **prec_data; /* length = # constraints */
  74      double *scratch; /* length = 2*n */
  75 } dual_data;
  76
  77 static double sqr(double x) { return x * x; }
  78
  79 static double dual_func(unsigned m, const double *y, double *grad, void *d_)
  80 {
  81      dual_data *d = (dual_data *) d_;
  82      unsigned n = d->n;
  83      const double *x = d->x, *lb = d->lb, *ub = d->ub, *sigma = d->sigma,
  84           *dfdx = d->dfdx;
  85      const double *dfcdx = d->dfcdx;
  86      double rho = d->rho, fval = d->fval;
  87      const double *rhoc = d->rhoc, *fcval = d->fcval;
  88      double *xcur = d->xcur;
  89      double *gcval = d->gcval;
  90      unsigned i, j;
  91      double val;
  92
  93      d->count++;
  94
  95      val = d->gval = fval;
  96      d->wval = 0;
  97      for (i = 0; i < m; ++i)
  98           val += y[i] * (gcval[i] = fcval[i]);
  99
 100      for (j = 0; j < n; ++j) {
 101           double u, v, dx, sigma2, dx2, dx2sig;
 102
 103           /* first, compute xcur[j] = x+dx for y.  Because this objective is
 104              separable, we can minimize over x analytically, and the minimum
 105              dx is given by the solution of a linear equation
 106                      u dx + v sigma^2 = 0, i.e. dx = -sigma^2 v/u
 107              where u and v are defined by the sums below.  However,
 108              we also have to check that |dx| <= sigma and that
 109              lb <= x+dx <= ub. */
 110
 111           if (sigma[j] == 0) { /* special case for lb[i] == ub[i] dims, dx=0 */
 112                xcur[j] = x[j];
 113                continue;
 114           }
 115
 116           u = rho;
 117           v = dfdx[j];
 118           for (i = 0; i < m; ++i) {
 119                u += rhoc[i] * y[i];
 120                v += dfcdx[i*n + j] * y[i];
 121           }
 122           dx = -(sigma2 = sqr(sigma[j])) * v/u;
 123
 124           /* if dx is out of bounds, we are guaranteed by convexity
 125              that the minimum is at the bound on the side of dx */
 126           if (fabs(dx) > sigma[j]) dx = copysign(sigma[j], dx);
 127           xcur[j] = x[j] + dx;
 128           if (xcur[j] > ub[j]) xcur[j] = ub[j];
 129           else if (xcur[j] < lb[j]) xcur[j] = lb[j];
 130           dx = xcur[j] - x[j];
 131
 132           /* function value: */
 133           dx2 = dx * dx;
 134           val += v * dx + 0.5 * u * dx2 / sigma2;
 135
 136           /* update gval, wval, gcval (approximant functions) */
 137           d->gval += dfdx[j] * dx + rho * (dx2sig = 0.5*dx2/sigma2);
 138           d->wval += dx2sig;
 139           for (i = 0; i < m; ++i)
 140                gcval[i] += dfcdx[i*n+j] * dx + rhoc[i] * dx2sig;
 141      }
 142
 143      /* gradient is easy to compute: since we are at a minimum x (dval/dx=0),
 144         we only need the partial derivative with respect to y, and
 145         we negate because we are maximizing: */
 146      if (grad) for (i = 0; i < m; ++i) grad[i] = -gcval[i];
 147      return -val;
 148 }
 149
 150 /***********************************************************************/
 151
 152 /* compute g(x - x0) and its gradient */
 153 static double gfunc(unsigned n, double f, const double *dfdx,
 154                     double rho, const double *sigma,
 155                     const double *x0,
 156                     nlopt_precond pre, void *pre_data, double *scratch,
 157                     const double *x, double *grad)
 158 {
 159      double *dx = scratch, *Hdx = scratch + n;
 160      double val = f;
 161      unsigned j;
 162
 163      for (j = 0; j < n; ++j) {
 164           double sigma2inv = 1.0 / sqr(sigma[j]);
 165           dx[j] = x[j] - x0[j];
 166           val += dfdx[j] * dx[j] + (0.5*rho) * sqr(dx[j]) * sigma2inv;
 167           if (grad) grad[j] = dfdx[j] + rho * dx[j] * sigma2inv;
 168      }
 169
 170      if (pre) {
 171           pre(n, x0, dx, Hdx, pre_data);
 172           for (j = 0; j < n; ++j)
 173                val += 0.5 * dx[j] * Hdx[j];
 174           if (grad)
 175                for (j = 0; j < n; ++j)
 176                     grad[j] += Hdx[j];
 177      }
 178
 179      return val;
 180 }
 181
 182 static double g0(unsigned n, const double *x, double *grad, void *d_)
 183 {
 184      dual_data *d = (dual_data *) d_;
 185      d->count++;
 186      return gfunc(n, d->fval, d->dfdx, d->rho, d->sigma,
 187                   d->x,
 188                   d->pre, d->pre_data, d->scratch,
 189                   x, grad);
 190 }
 191
 192
 193 static void gi(unsigned m, double *result,
 194                unsigned n, const double *x, double *grad, void *d_)
 195 {
 196      dual_data *d = (dual_data *) d_;
 197      unsigned i;
 198      for (i = 0; i < m; ++i)
 199           result[i] = gfunc(n, d->fcval[i], d->dfcdx + i*n, d->rhoc[i],
 200                             d->sigma,
 201                             d->x,
 202                             d->prec ? d->prec[i] : NULL,
 203                             d->prec_data ? d->prec_data[i] : NULL,
 204                             d->scratch,
 205                             x, grad);
 206 }
 207
 208
 209 /***********************************************************************/
 210
 211 nlopt_result ccsa_quadratic_minimize(
 212      unsigned n, nlopt_func f, void *f_data,
 213      unsigned m, nlopt_constraint *fc,
 214
 215      nlopt_precond pre,
 216
 217      const double *lb, const double *ub, /* bounds */
 218      double *x, /* in: initial guess, out: minimizer */
 219      double *minf,
 220      nlopt_stopping *stop,
 221      nlopt_opt dual_opt)
 222 {
 223      nlopt_result ret = NLOPT_SUCCESS;
 224      double *xcur, rho, *sigma, *dfdx, *dfdx_cur, *xprev, *xprevprev, fcur;
 225      double *dfcdx, *dfcdx_cur;
 226      double *fcval, *fcval_cur, *rhoc, *gcval, *y, *dual_lb, *dual_ub;
 227      double *pre_lb = NULL, *pre_ub = NULL;
 228      unsigned i, ifc, j, k = 0;
 229      dual_data dd;
 230      int feasible;
 231      double infeasibility;
 232      unsigned mfc;
 233      unsigned no_precond;
 234      nlopt_opt pre_opt = NULL;
 235
 236      m = nlopt_count_constraints(mfc = m, fc);
 237      if (nlopt_get_dimension(dual_opt) != m) {
 238          nlopt_stop_msg(stop, "dual optimizer has wrong dimension %d != %d",
 239                         nlopt_get_dimension(dual_opt), m);
 240          return NLOPT_INVALID_ARGS;
 241      }
 242      sigma = (double *) malloc(sizeof(double) * (6*n + 2*m*n + m*7));
 243      if (!sigma) return NLOPT_OUT_OF_MEMORY;
 244      dfdx = sigma + n;
 245      dfdx_cur = dfdx + n;
 246      xcur = dfdx_cur + n;
 247      xprev = xcur + n;
 248      xprevprev = xprev + n;
 249      fcval = xprevprev + n;
 250      fcval_cur = fcval + m;
 251      rhoc = fcval_cur + m;
 252      gcval = rhoc + m;
 253      dual_lb = gcval + m;
 254      dual_ub = dual_lb + m;
 255      y = dual_ub + m;
 256      dfcdx = y + m;
 257      dfcdx_cur = dfcdx + m*n;
 258
 259      dd.n = n;
 260      dd.x = x;
 261      dd.lb = lb;
 262      dd.ub = ub;
 263      dd.sigma = sigma;
 264      dd.dfdx = dfdx;
 265      dd.dfcdx = dfcdx;
 266      dd.fcval = fcval;
 267      dd.rhoc = rhoc;
 268      dd.xcur = xcur;
 269      dd.gcval = gcval;
 270      dd.pre = pre; dd.pre_data = f_data;
 271      dd.prec = NULL; dd.prec_data = NULL;
 272      dd.scratch = NULL;
 273
 274      if (m) {
 275           dd.prec = (nlopt_precond *) malloc(sizeof(nlopt_precond) * m);
 276           dd.prec_data = (void **) malloc(sizeof(void *) * m);
 277           if (!dd.prec || !dd.prec_data) {
 278                ret = NLOPT_OUT_OF_MEMORY;
 279                goto done;
 280           }
 281           for (i = ifc = 0; ifc < mfc; ++ifc) {
 282                unsigned inext = i + fc[ifc].m;
 283                for (; i < inext; ++i) {
 284                     dd.prec[i] = fc[ifc].pre;
 285                     dd.prec_data[i] = fc[ifc].f_data;
 286                }
 287           }
 288      }
 289
 290      no_precond = pre == NULL;
 291      if (dd.prec)
 292           for (i = 0; i < m; ++i)
 293                no_precond = no_precond && dd.prec[i] == NULL;
 294
 295      if (!no_precond) {
 296           dd.scratch = (double*) malloc(sizeof(double) * (4*n));
 297           if (!dd.scratch) {
 298                free(sigma);
 299                return NLOPT_OUT_OF_MEMORY;
 300           }
 301           pre_lb = dd.scratch + 2*n;
 302           pre_ub = pre_lb + n;
 303
 304           pre_opt = nlopt_create(nlopt_get_algorithm(dual_opt), n);
 305           if (!pre_opt) {
 306               nlopt_stop_msg(stop, "failure creating precond. optimizer");
 307               ret = NLOPT_FAILURE;
 308               goto done;
 309           }
 310           ret = nlopt_set_min_objective(pre_opt, g0, &dd);
 311           if (ret < 0) goto done;
 312           ret = nlopt_add_inequality_mconstraint(pre_opt, m, gi, &dd, NULL);
 313           if (ret < 0) goto done;
 314           ret = nlopt_set_ftol_rel(pre_opt, nlopt_get_ftol_rel(dual_opt));
 315           if (ret < 0) goto done;
 316           ret = nlopt_set_ftol_abs(pre_opt, nlopt_get_ftol_abs(dual_opt));
 317           if (ret < 0) goto done;
 318           ret = nlopt_set_maxeval(pre_opt, nlopt_get_maxeval(dual_opt));
 319           if (ret < 0) goto done;
 320      }
 321
 322      for (j = 0; j < n; ++j) {
 323           if (nlopt_isinf(ub[j]) || nlopt_isinf(lb[j]))
 324                sigma[j] = 1.0; /* arbitrary default */
 325           else
 326                sigma[j] = 0.5 * (ub[j] - lb[j]);
 327      }
 328      rho = 1.0;
 329      for (i = 0; i < m; ++i) {
 330           rhoc[i] = 1.0;
 331           dual_lb[i] = y[i] = 0.0;
 332           dual_ub[i] = HUGE_VAL;
 333      }
 334
 335      dd.fval = fcur = *minf = f(n, x, dfdx, f_data);
 336      ++ *(stop->nevals_p);
 337      memcpy(xcur, x, sizeof(double) * n);
 338      if (nlopt_stop_forced(stop)) { ret = NLOPT_FORCED_STOP; goto done; }
 339
 340      feasible = 1; infeasibility = 0;
 341      for (i = ifc = 0; ifc < mfc; ++ifc) {
 342           nlopt_eval_constraint(fcval + i, dfcdx + i*n,
 343                                 fc + ifc, n, x);
 344           i += fc[ifc].m;
 345           if (nlopt_stop_forced(stop)) { ret = NLOPT_FORCED_STOP; goto done; }
 346      }
 347      for (i = 0; i < m; ++i) {
 348           feasible = feasible && fcval[i] <= 0;
 349           if (fcval[i] > infeasibility) infeasibility = fcval[i];
 350      }
 351      /* For non-feasible initial points, set a finite (large)
 352         upper-bound on the dual variables.  What this means is that,
 353         if no feasible solution is found from the dual problem, it
 354         will minimize the dual objective with the unfeasible
 355         constraint weighted by 1e40 -- basically, minimizing the
 356         unfeasible constraint until it becomes feasible or until we at
 357         least obtain a step towards a feasible point.
 358
 359         Svanberg suggested a different approach in his 1987 paper, basically
 360         introducing additional penalty variables for unfeasible constraints,
 361         but this is easier to implement and at least as efficient. */
 362      if (!feasible)
 363           for (i = 0; i < m; ++i) dual_ub[i] = 1e40;
 364
 365      nlopt_set_min_objective(dual_opt, dual_func, &dd);
 366      nlopt_set_lower_bounds(dual_opt, dual_lb);
 367      nlopt_set_upper_bounds(dual_opt, dual_ub);
 368      nlopt_set_stopval(dual_opt, -HUGE_VAL);
 369      nlopt_remove_inequality_constraints(dual_opt);
 370      nlopt_remove_equality_constraints(dual_opt);
 371
 372      while (1) { /* outer iterations */
 373           double fprev = fcur;
 374           if (nlopt_stop_forced(stop)) ret = NLOPT_FORCED_STOP;
 375           else if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
 376           else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
 377           else if (feasible && *minf < stop->minf_max)
 378                ret = NLOPT_MINF_MAX_REACHED;
 379           if (ret != NLOPT_SUCCESS) goto done;
 380           if (++k > 1) memcpy(xprevprev, xprev, sizeof(double) * n);
 381           memcpy(xprev, xcur, sizeof(double) * n);
 382
 383           while (1) { /* inner iterations */
 384                double min_dual, infeasibility_cur;
 385                int feasible_cur, inner_done;
 386                unsigned save_verbose;
 387                nlopt_result reti;
 388
 389                if (no_precond) {
 390                     /* solve dual problem */
 391                     dd.rho = rho; dd.count = 0;
 392                     save_verbose = ccsa_verbose;
 393                     ccsa_verbose = 0; /* no recursive verbosity */
 394                     reti = nlopt_optimize_limited(dual_opt, y, &min_dual,
 395                                                   0,
 396                                                   stop->maxtime
 397                                                   - (nlopt_seconds()
 398                                                      - stop->start));
 399                     ccsa_verbose = save_verbose;
 400                     if (reti < 0 || reti == NLOPT_MAXTIME_REACHED) {
 401                          ret = reti;
 402                          goto done;
 403                     }
 404
 405                     dual_func(m, y, NULL, &dd); /* evaluate final xcur etc. */
 406                }
 407                else {
 408                     double pre_min;
 409                     for (j = 0; j < n; ++j) {
 410                          pre_lb[j] = MAX(lb[j], x[j] - sigma[j]);
 411                          pre_ub[j] = MIN(ub[j], x[j] + sigma[j]);
 412                          xcur[j] = x[j];
 413                     }
 414                     nlopt_set_lower_bounds(pre_opt, pre_lb);
 415                     nlopt_set_upper_bounds(pre_opt, pre_ub);
 416
 417                     dd.rho = rho; dd.count = 0;
 418                     save_verbose = ccsa_verbose;
 419                     ccsa_verbose = 0; /* no recursive verbosity */
 420                     reti = nlopt_optimize_limited(pre_opt, xcur, &pre_min,
 421                                                   0, stop->maxtime
 422                                                   - (nlopt_seconds()
 423                                                      - stop->start));
 424                     ccsa_verbose = save_verbose;
 425                     if (reti < 0 || reti == NLOPT_MAXTIME_REACHED) {
 426                          ret = reti;
 427                          goto done;
 428                     }
 429
 430                     /* evaluate final xcur etc */
 431                     dd.gval = g0(n, xcur, NULL, &dd);
 432                     gi(m, dd.gcval, n, xcur, NULL, &dd);
 433                }
 434
 435                if (ccsa_verbose) {
 436                     printf("CCSA dual converged in %d iters to g=%g:\n",
 437                            dd.count, dd.gval);
 438                     for (i = 0; i < MIN(ccsa_verbose, m); ++i)
 439                          printf("    CCSA y[%u]=%g, gc[%u]=%g\n",
 440                                 i, y[i], i, dd.gcval[i]);
 441                }
 442
 443                fcur = f(n, xcur, dfdx_cur, f_data);
 444                ++ *(stop->nevals_p);
 445                if (nlopt_stop_forced(stop)) {
 446                     ret = NLOPT_FORCED_STOP; goto done; }
 447                feasible_cur = 1; infeasibility_cur = 0;
 448                inner_done = dd.gval >= fcur;
 449                for (i = ifc = 0; ifc < mfc; ++ifc) {
 450                     nlopt_eval_constraint(fcval_cur + i, dfcdx_cur + i*n,
 451                                           fc + ifc, n, xcur);
 452                     i += fc[ifc].m;
 453                     if (nlopt_stop_forced(stop)) {
 454                          ret = NLOPT_FORCED_STOP; goto done; }
 455                }
 456                for (i = ifc = 0; ifc < mfc; ++ifc) {
 457                     unsigned i0 = i, inext = i + fc[ifc].m;
 458                     for (; i < inext; ++i) {
 459                          feasible_cur = feasible_cur
 460                               && fcval_cur[i] <= fc[ifc].tol[i-i0];
 461                          inner_done = inner_done &&
 462                               (dd.gcval[i] >= fcval_cur[i]);
 463                          if (fcval_cur[i] > infeasibility_cur)
 464                               infeasibility_cur = fcval_cur[i];
 465                     }
 466                }
 467
 468                if ((fcur < *minf && (inner_done || feasible_cur || !feasible))
 469                     || (!feasible && infeasibility_cur < infeasibility)) {
 470                     if (ccsa_verbose && !feasible_cur)
 471                          printf("CCSA - using infeasible point?\n");
 472                     dd.fval = *minf = fcur;
 473                     infeasibility = infeasibility_cur;
 474                     memcpy(fcval, fcval_cur, sizeof(double)*m);
 475                     memcpy(x, xcur, sizeof(double)*n);
 476                     memcpy(dfdx, dfdx_cur, sizeof(double)*n);
 477                     memcpy(dfcdx, dfcdx_cur, sizeof(double)*n*m);
 478
 479                     /* once we have reached a feasible solution, the
 480                        algorithm should never make the solution infeasible
 481                        again (if inner_done), although the constraints may
 482                        be violated slightly by rounding errors etc. so we
 483                        must be a little careful about checking feasibility */
 484                     if (infeasibility_cur == 0) {
 485                          if (!feasible) { /* reset upper bounds to infin. */
 486                               for (i = 0; i < m; ++i) dual_ub[i] = HUGE_VAL;
 487                               nlopt_set_upper_bounds(dual_opt, dual_ub);
 488                          }
 489                          feasible = 1;
 490                     }
 491
 492                }
 493                if (nlopt_stop_forced(stop)) ret = NLOPT_FORCED_STOP;
 494                else if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
 495                else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
 496                else if (feasible && *minf < stop->minf_max)
 497                     ret = NLOPT_MINF_MAX_REACHED;
 498                if (ret != NLOPT_SUCCESS) goto done;
 499
 500                if (inner_done) break;
 501
 502                if (fcur > dd.gval)
 503                     rho = MIN(10*rho, 1.1 * (rho + (fcur-dd.gval) / dd.wval));
 504                for (i = 0; i < m; ++i)
 505                     if (fcval_cur[i] > dd.gcval[i])
 506                          rhoc[i] =
 507                               MIN(10*rhoc[i],
 508                                   1.1 * (rhoc[i] + (fcval_cur[i]-dd.gcval[i])
 509                                          / dd.wval));
 510
 511                if (ccsa_verbose)
 512                     printf("CCSA inner iteration: rho -> %g\n", rho);
 513                for (i = 0; i < MIN(ccsa_verbose, m); ++i)
 514                     printf("                CCSA rhoc[%u] -> %g\n", i,rhoc[i]);
 515           }
 516
 517           if (nlopt_stop_ftol(stop, fcur, fprev))
 518                ret = NLOPT_FTOL_REACHED;
 519           if (nlopt_stop_x(stop, xcur, xprev))
 520                ret = NLOPT_XTOL_REACHED;
 521           if (ret != NLOPT_SUCCESS) goto done;
 522
 523           /* update rho and sigma for iteration k+1 */
 524           rho = MAX(0.1 * rho, CCSA_RHOMIN);
 525           if (ccsa_verbose)
 526                printf("CCSA outer iteration: rho -> %g\n", rho);
 527           for (i = 0; i < m; ++i)
 528                rhoc[i] = MAX(0.1 * rhoc[i], CCSA_RHOMIN);
 529           for (i = 0; i < MIN(ccsa_verbose, m); ++i)
 530                printf("                 CCSA rhoc[%u] -> %g\n", i, rhoc[i]);
 531           if (k > 1) {
 532                for (j = 0; j < n; ++j) {
 533                     double dx2 = (xcur[j]-xprev[j]) * (xprev[j]-xprevprev[j]);
 534                     double gam = dx2 < 0 ? 0.7 : (dx2 > 0 ? 1.2 : 1);
 535                     sigma[j] *= gam;
 536                     if (!nlopt_isinf(ub[j]) && !nlopt_isinf(lb[j])) {
 537                          sigma[j] = MIN(sigma[j], 10*(ub[j]-lb[j]));
 538                          /* use a smaller lower bound than Svanberg's
 539                             0.01*(ub-lb), which seems unnecessarily large */
 540                          sigma[j] = MAX(sigma[j], 1e-8*(ub[j]-lb[j]));
 541                     }
 542                }
 543                for (j = 0; j < MIN(ccsa_verbose, n); ++j)
 544                     printf("                 CCSA sigma[%u] -> %g\n",
 545                            j, sigma[j]);
 546           }
 547      }
 548
 549  done:
 550      nlopt_destroy(pre_opt);
 551      if (dd.scratch) free(dd.scratch);
 552      if (m) {
 553           free(dd.prec_data);
 554           free(dd.prec);
 555      }
 556      free(sigma);
 557      return ret;
 558 }