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