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, *pre_ub;
 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) return NLOPT_INVALID_ARGS;
 238      sigma = (double *) malloc(sizeof(double) * (6*n + 2*m*n + m*7));
 239      if (!sigma) return NLOPT_OUT_OF_MEMORY;
 240      dfdx = sigma + n;
 241      dfdx_cur = dfdx + n;
 242      xcur = dfdx_cur + n;
 243      xprev = xcur + n;
 244      xprevprev = xprev + n;
 245      fcval = xprevprev + n;
 246      fcval_cur = fcval + m;
 247      rhoc = fcval_cur + m;
 248      gcval = rhoc + m;
 249      dual_lb = gcval + m;
 250      dual_ub = dual_lb + m;
 251      y = dual_ub + m;
 252      dfcdx = y + m;
 253      dfcdx_cur = dfcdx + m*n;
 254
 255      dd.n = n;
 256      dd.x = x;
 257      dd.lb = lb;
 258      dd.ub = ub;
 259      dd.sigma = sigma;
 260      dd.dfdx = dfdx;
 261      dd.dfcdx = dfcdx;
 262      dd.fcval = fcval;
 263      dd.rhoc = rhoc;
 264      dd.xcur = xcur;
 265      dd.gcval = gcval;
 266      dd.pre = pre; dd.pre_data = f_data;
 267      dd.prec = NULL; dd.prec_data = NULL;
 268      dd.scratch = NULL;
 269
 270      if (m) {
 271           dd.prec = (nlopt_precond *) malloc(sizeof(nlopt_precond) * m);
 272           dd.prec_data = (void **) malloc(sizeof(void *) * m);
 273           if (!dd.prec || !dd.prec_data) {
 274                ret = NLOPT_OUT_OF_MEMORY;
 275                goto done;
 276           }
 277           for (i = ifc = 0; ifc < mfc; ++ifc) {
 278                unsigned inext = i + fc[ifc].m;
 279                for (; i < inext; ++i) {
 280                     dd.prec[i] = fc[ifc].pre;
 281                     dd.prec_data[i] = fc[ifc].f_data;
 282                }
 283           }
 284      }
 285
 286      no_precond = pre == NULL;
 287      if (dd.prec)
 288           for (i = 0; i < m; ++i)
 289                no_precond = no_precond && dd.prec[i] == NULL;
 290
 291      if (!no_precond) {
 292           dd.scratch = (double*) malloc(sizeof(double) * (4*n));
 293           if (!dd.scratch) {
 294                free(sigma);
 295                return NLOPT_OUT_OF_MEMORY;
 296           }
 297           pre_lb = dd.scratch + 2*n;
 298           pre_ub = pre_lb + n;
 299
 300           pre_opt = nlopt_create(nlopt_get_algorithm(dual_opt), n);
 301           if (!pre_opt) { ret = NLOPT_FAILURE; goto done; }
 302           ret = nlopt_set_min_objective(pre_opt, g0, &dd);
 303           if (ret < 0) goto done;
 304           ret = nlopt_add_inequality_mconstraint(pre_opt, m, gi, &dd, NULL);
 305           if (ret < 0) goto done;
 306           ret = nlopt_set_ftol_rel(pre_opt, nlopt_get_ftol_rel(dual_opt));
 307           if (ret < 0) goto done;
 308           ret = nlopt_set_ftol_abs(pre_opt, nlopt_get_ftol_abs(dual_opt));
 309           if (ret < 0) goto done;
 310           ret = nlopt_set_maxeval(pre_opt, nlopt_get_maxeval(dual_opt));
 311           if (ret < 0) goto done;
 312      }
 313
 314      for (j = 0; j < n; ++j) {
 315           if (nlopt_isinf(ub[j]) || nlopt_isinf(lb[j]))
 316                sigma[j] = 1.0; /* arbitrary default */
 317           else
 318                sigma[j] = 0.5 * (ub[j] - lb[j]);
 319      }
 320      rho = 1.0;
 321      for (i = 0; i < m; ++i) {
 322           rhoc[i] = 1.0;
 323           dual_lb[i] = y[i] = 0.0;
 324           dual_ub[i] = HUGE_VAL;
 325      }
 326
 327      dd.fval = fcur = *minf = f(n, x, dfdx, f_data);
 328      stop->nevals++;
 329      memcpy(xcur, x, sizeof(double) * n);
 330      if (nlopt_stop_forced(stop)) { ret = NLOPT_FORCED_STOP; goto done; }
 331
 332      feasible = 1; infeasibility = 0;
 333      for (i = ifc = 0; ifc < mfc; ++ifc) {
 334           nlopt_eval_constraint(fcval + i, dfcdx + i*n,
 335                                 fc + ifc, n, x);
 336           i += fc[ifc].m;
 337           if (nlopt_stop_forced(stop)) { ret = NLOPT_FORCED_STOP; goto done; }
 338      }
 339      for (i = 0; i < m; ++i) {
 340           feasible = feasible && fcval[i] <= 0;
 341           if (fcval[i] > infeasibility) infeasibility = fcval[i];
 342      }
 343      /* For non-feasible initial points, set a finite (large)
 344         upper-bound on the dual variables.  What this means is that,
 345         if no feasible solution is found from the dual problem, it
 346         will minimize the dual objective with the unfeasible
 347         constraint weighted by 1e40 -- basically, minimizing the
 348         unfeasible constraint until it becomes feasible or until we at
 349         least obtain a step towards a feasible point.
 350
 351         Svanberg suggested a different approach in his 1987 paper, basically
 352         introducing additional penalty variables for unfeasible constraints,
 353         but this is easier to implement and at least as efficient. */
 354      if (!feasible)
 355           for (i = 0; i < m; ++i) dual_ub[i] = 1e40;
 356
 357      nlopt_set_min_objective(dual_opt, dual_func, &dd);
 358      nlopt_set_lower_bounds(dual_opt, dual_lb);
 359      nlopt_set_upper_bounds(dual_opt, dual_ub);
 360      nlopt_set_stopval(dual_opt, -HUGE_VAL);
 361      nlopt_remove_inequality_constraints(dual_opt);
 362      nlopt_remove_equality_constraints(dual_opt);
 363
 364      while (1) { /* outer iterations */
 365           double fprev = fcur;
 366           if (nlopt_stop_forced(stop)) ret = NLOPT_FORCED_STOP;
 367           else if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
 368           else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
 369           else if (feasible && *minf < stop->minf_max)
 370                ret = NLOPT_MINF_MAX_REACHED;
 371           if (ret != NLOPT_SUCCESS) goto done;
 372           if (++k > 1) memcpy(xprevprev, xprev, sizeof(double) * n);
 373           memcpy(xprev, xcur, sizeof(double) * n);
 374
 375           while (1) { /* inner iterations */
 376                double min_dual, infeasibility_cur;
 377                int feasible_cur, inner_done;
 378                unsigned save_verbose;
 379                nlopt_result reti;
 380
 381                if (no_precond) {
 382                     /* solve dual problem */
 383                     dd.rho = rho; dd.count = 0;
 384                     save_verbose = ccsa_verbose;
 385                     ccsa_verbose = 0; /* no recursive verbosity */
 386                     reti = nlopt_optimize_limited(dual_opt, y, &min_dual,
 387                                                   0,
 388                                                   stop->maxtime
 389                                                   - (nlopt_seconds()
 390                                                      - stop->start));
 391                     ccsa_verbose = save_verbose;
 392                     if (reti < 0 || reti == NLOPT_MAXTIME_REACHED) {
 393                          ret = reti;
 394                          goto done;
 395                     }
 396
 397                     dual_func(m, y, NULL, &dd); /* evaluate final xcur etc. */
 398                }
 399                else {
 400                     double pre_min;
 401                     for (j = 0; j < n; ++j) {
 402                          pre_lb[j] = MAX(lb[j], x[j] - sigma[j]);
 403                          pre_ub[j] = MIN(ub[j], x[j] + sigma[j]);
 404                          xcur[j] = x[j];
 405                     }
 406                     nlopt_set_lower_bounds(pre_opt, pre_lb);
 407                     nlopt_set_upper_bounds(pre_opt, pre_ub);
 408
 409                     dd.rho = rho; dd.count = 0;
 410                     save_verbose = ccsa_verbose;
 411                     ccsa_verbose = 0; /* no recursive verbosity */
 412                     reti = nlopt_optimize_limited(pre_opt, xcur, &pre_min,
 413                                                   0, stop->maxtime
 414                                                   - (nlopt_seconds()
 415                                                      - stop->start));
 416                     ccsa_verbose = save_verbose;
 417                     if (reti < 0 || reti == NLOPT_MAXTIME_REACHED) {
 418                          ret = reti;
 419                          goto done;
 420                     }
 421
 422                     /* evaluate final xcur etc */
 423                     dd.gval = g0(n, xcur, NULL, &dd);
 424                     gi(m, dd.gcval, n, xcur, NULL, &dd);
 425                }
 426
 427                if (ccsa_verbose) {
 428                     printf("CCSA dual converged in %d iters to g=%g:\n",
 429                            dd.count, dd.gval);
 430                     for (i = 0; i < MIN(ccsa_verbose, m); ++i)
 431                          printf("    CCSA y[%d]=%g, gc[%d]=%g\n",
 432                                 i, y[i], i, dd.gcval[i]);
 433                }
 434
 435                fcur = f(n, xcur, dfdx_cur, f_data);
 436                stop->nevals++;
 437                if (nlopt_stop_forced(stop)) {
 438                     ret = NLOPT_FORCED_STOP; goto done; }
 439                feasible_cur = 1; infeasibility_cur = 0;
 440                inner_done = dd.gval >= fcur;
 441                for (i = ifc = 0; ifc < mfc; ++ifc) {
 442                     nlopt_eval_constraint(fcval_cur + i, dfcdx_cur + i*n,
 443                                           fc + ifc, n, xcur);
 444                     i += fc[ifc].m;
 445                     if (nlopt_stop_forced(stop)) {
 446                          ret = NLOPT_FORCED_STOP; goto done; }
 447                }
 448                for (i = ifc = 0; ifc < mfc; ++ifc) {
 449                     unsigned i0 = i, inext = i + fc[ifc].m;
 450                     for (; i < inext; ++i) {
 451                          feasible_cur = feasible_cur
 452                               && fcval_cur[i] <= fc[ifc].tol[i-i0];
 453                          inner_done = inner_done &&
 454                               (dd.gcval[i] >= fcval_cur[i]);
 455                          if (fcval_cur[i] > infeasibility_cur)
 456                               infeasibility_cur = fcval_cur[i];
 457                     }
 458                }
 459
 460                if ((fcur < *minf && (inner_done || feasible_cur || !feasible))
 461                     || (!feasible && infeasibility_cur < infeasibility)) {
 462                     if (ccsa_verbose && !feasible_cur)
 463                          printf("CCSA - using infeasible point?\n");
 464                     dd.fval = *minf = fcur;
 465                     infeasibility = infeasibility_cur;
 466                     memcpy(fcval, fcval_cur, sizeof(double)*m);
 467                     memcpy(x, xcur, sizeof(double)*n);
 468                     memcpy(dfdx, dfdx_cur, sizeof(double)*n);
 469                     memcpy(dfcdx, dfcdx_cur, sizeof(double)*n*m);
 470
 471                     /* once we have reached a feasible solution, the
 472                        algorithm should never make the solution infeasible
 473                        again (if inner_done), although the constraints may
 474                        be violated slightly by rounding errors etc. so we
 475                        must be a little careful about checking feasibility */
 476                     if (infeasibility_cur == 0) {
 477                          if (!feasible) { /* reset upper bounds to infin. */
 478                               for (i = 0; i < m; ++i) dual_ub[i] = HUGE_VAL;
 479                               nlopt_set_upper_bounds(dual_opt, dual_ub);
 480                          }
 481                          feasible = 1;
 482                     }
 483
 484                }
 485                if (nlopt_stop_forced(stop)) ret = NLOPT_FORCED_STOP;
 486                else if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
 487                else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
 488                else if (feasible && *minf < stop->minf_max)
 489                     ret = NLOPT_MINF_MAX_REACHED;
 490                if (ret != NLOPT_SUCCESS) goto done;
 491
 492                if (inner_done) break;
 493
 494                if (fcur > dd.gval)
 495                     rho = MIN(10*rho, 1.1 * (rho + (fcur-dd.gval) / dd.wval));
 496                for (i = 0; i < m; ++i)
 497                     if (fcval_cur[i] > dd.gcval[i])
 498                          rhoc[i] =
 499                               MIN(10*rhoc[i],
 500                                   1.1 * (rhoc[i] + (fcval_cur[i]-dd.gcval[i])
 501                                          / dd.wval));
 502
 503                if (ccsa_verbose)
 504                     printf("CCSA inner iteration: rho -> %g\n", rho);
 505                for (i = 0; i < MIN(ccsa_verbose, m); ++i)
 506                     printf("                CCSA rhoc[%d] -> %g\n", i,rhoc[i]);
 507           }
 508
 509           if (nlopt_stop_ftol(stop, fcur, fprev))
 510                ret = NLOPT_FTOL_REACHED;
 511           if (nlopt_stop_x(stop, xcur, xprev))
 512                ret = NLOPT_XTOL_REACHED;
 513           if (ret != NLOPT_SUCCESS) goto done;
 514
 515           /* update rho and sigma for iteration k+1 */
 516           rho = MAX(0.1 * rho, CCSA_RHOMIN);
 517           if (ccsa_verbose)
 518                printf("CCSA outer iteration: rho -> %g\n", rho);
 519           for (i = 0; i < m; ++i)
 520                rhoc[i] = MAX(0.1 * rhoc[i], CCSA_RHOMIN);
 521           for (i = 0; i < MIN(ccsa_verbose, m); ++i)
 522                printf("                 CCSA rhoc[%d] -> %g\n", i, rhoc[i]);
 523           if (k > 1) {
 524                for (j = 0; j < n; ++j) {
 525                     double dx2 = (xcur[j]-xprev[j]) * (xprev[j]-xprevprev[j]);
 526                     double gam = dx2 < 0 ? 0.7 : (dx2 > 0 ? 1.2 : 1);
 527                     sigma[j] *= gam;
 528                     if (!nlopt_isinf(ub[j]) && !nlopt_isinf(lb[j])) {
 529                          sigma[j] = MIN(sigma[j], 10*(ub[j]-lb[j]));
 530                          /* use a smaller lower bound than Svanberg's
 531                             0.01*(ub-lb), which seems unnecessarily large */
 532                          sigma[j] = MAX(sigma[j], 1e-8*(ub[j]-lb[j]));
 533                     }
 534                }
 535                for (j = 0; j < MIN(ccsa_verbose, n); ++j)
 536                     printf("                 CCSA sigma[%d] -> %g\n",
 537                            j, sigma[j]);
 538           }
 539      }
 540
 541  done:
 542      nlopt_destroy(pre_opt);
 543      if (dd.scratch) free(dd.scratch);
 544      if (m) {
 545           free(dd.prec_data);
 546           free(dd.prec);
 547      }
 548      free(sigma);
 549      return ret;
 550 }