vendor/scs/src/scs.c in scs-0.4.0 vs vendor/scs/src/scs.c in scs-0.4.1
- old
+ new
@@ -74,11 +74,11 @@
static void print_init_header(const ScsData *d, const ScsCone *k,
const ScsSettings *stgs) {
scs_int i;
char *cone_str = SCS(get_cone_header)(k);
- const char *lin_sys_method = SCS(get_lin_sys_method)();
+ const char *lin_sys_method = scs_get_lin_sys_method();
#ifdef USE_LAPACK
scs_int acceleration_lookback = stgs->acceleration_lookback;
scs_int acceleration_interval = stgs->acceleration_interval;
#else
scs_int acceleration_lookback = 0;
@@ -193,21 +193,28 @@
SCS(un_normalize_sol)(w->scal, sol);
}
}
static void compute_residuals(ScsResiduals *r, scs_int m, scs_int n) {
- r->res_pri = SAFEDIV_POS(NORM(r->ax_s_btau, m), r->tau);
- r->res_dual = SAFEDIV_POS(NORM(r->px_aty_ctau, n), r->tau);
+ scs_float nm_ax_s, nm_px, nm_aty;
+ scs_float nm_ax_s_btau = NORM(r->ax_s_btau, m);
+ scs_float nm_px_aty_ctau = NORM(r->px_aty_ctau, n);
+
+ r->res_pri = SAFEDIV_POS(nm_ax_s_btau, r->tau);
+ r->res_dual = SAFEDIV_POS(nm_px_aty_ctau, r->tau);
r->res_unbdd_a = NAN;
r->res_unbdd_p = NAN;
r->res_infeas = NAN;
if (r->ctx_tau < 0) {
- r->res_unbdd_a = SAFEDIV_POS(NORM(r->ax_s, m), -r->ctx_tau);
- r->res_unbdd_p = SAFEDIV_POS(NORM(r->px, n), -r->ctx_tau);
+ nm_ax_s = NORM(r->ax_s, m);
+ nm_px = NORM(r->px, n);
+ r->res_unbdd_a = SAFEDIV_POS(nm_ax_s, -r->ctx_tau);
+ r->res_unbdd_p = SAFEDIV_POS(nm_px, -r->ctx_tau);
}
if (r->bty_tau < 0) {
- r->res_infeas = SAFEDIV_POS(NORM(r->aty, n), -r->bty_tau);
+ nm_aty = NORM(r->aty, n);
+ r->res_infeas = SAFEDIV_POS(nm_aty, -r->bty_tau);
}
}
static void unnormalize_residuals(ScsWork *w) {
ScsResiduals *r_n = w->r_normalized; /* normalized residuals */
@@ -345,43 +352,46 @@
return ip;
}
static scs_float root_plus(ScsWork *w, scs_float *p, scs_float *mu,
scs_float eta) {
- scs_float a, b, c, tau_scale = w->diag_r[w->d->n + w->d->m];
+ scs_float a, b, c, rad, tau_scale = w->diag_r[w->d->n + w->d->m];
a = tau_scale + dot_r(w, w->g, w->g);
b = dot_r(w, mu, w->g) - 2 * dot_r(w, p, w->g) - eta * tau_scale;
c = dot_r(w, p, p) - dot_r(w, p, mu);
- return (-b + SQRTF(MAX(b * b - 4 * a * c, 0.))) / (2 * a);
+ rad = b * b - 4 * a * c;
+ return (-b + SQRTF(MAX(rad, 0.))) / (2 * a);
}
-/* status < 0 indicates failure */
+/* status != 0 indicates failure */
static scs_int project_lin_sys(ScsWork *w, scs_int iter) {
scs_int n = w->d->n, m = w->d->m, l = n + m + 1, status, i;
scs_float *warm_start = SCS_NULL;
scs_float tol = -1.0; /* only used for indirect methods, overridden later */
memcpy(w->u_t, w->v, l * sizeof(scs_float));
for (i = 0; i < l - 1; ++i) {
w->u_t[i] *= (i < n ? 1 : -1) * w->diag_r[i];
}
#if INDIRECT > 0
+ scs_float nm_ax_s_btau, nm_px_aty_ctau, nm_ws;
/* compute warm start using the cone projection output */
+ nm_ax_s_btau = CG_NORM(w->r_normalized->ax_s_btau, m);
+ nm_px_aty_ctau = CG_NORM(w->r_normalized->px_aty_ctau, n);
warm_start = w->lin_sys_warm_start;
- memcpy(warm_start, w->u, (l - 1) * sizeof(scs_float));
/* warm_start = u[:n] + tau * g[:n] */
- SCS(add_scaled_array)(warm_start, w->g, l - 1, w->u[l - 1]);
+ memcpy(warm_start, w->u, n * sizeof(scs_float));
+ SCS(add_scaled_array)(warm_start, w->g, n, w->u[l - 1]);
/* use normalized residuals to compute tolerance */
- tol = MIN(CG_NORM(w->r_normalized->ax_s_btau, w->d->m),
- CG_NORM(w->r_normalized->px_aty_ctau, w->d->n));
+ tol = MIN(nm_ax_s_btau, nm_px_aty_ctau);
/* tol ~ O(1/k^(1+eps)) guarantees convergence */
/* use warm-start to calculate tolerance rather than w->u_t, since warm_start
* should be approximately equal to the true solution */
- tol = CG_TOL_FACTOR * MIN(tol, CG_NORM(warm_start, w->d->n) /
- POWF((scs_float)iter + 1, CG_RATE));
+ nm_ws = CG_NORM(warm_start, n) / POWF((scs_float)iter + 1, CG_RATE);
+ tol = CG_TOL_FACTOR * MIN(tol, nm_ws);
tol = MAX(CG_BEST_TOL, tol);
#endif
- status = SCS(solve_lin_sys)(w->p, w->u_t, warm_start, tol);
+ status = scs_solve_lin_sys(w->p, w->u_t, warm_start, tol);
if (iter < FEASIBLE_ITERS) {
w->u_t[l - 1] = 1.;
} else {
w->u_t[l - 1] = root_plus(w, w->u_t, w->v, w->v[l - 1]);
}
@@ -512,29 +522,34 @@
}
/* sets solutions, re-scales by inner prods if infeasible or unbounded */
static void finalize(ScsWork *w, ScsSolution *sol, ScsInfo *info,
scs_int iter) {
+ scs_float nm_s, nm_y, sty;
setx(w, sol);
sety(w, sol);
sets(w, sol);
if (w->stgs->normalize) {
SCS(un_normalize_sol)(w->scal, sol);
}
populate_residual_struct(w, iter);
+
+ nm_s = SCS(norm_inf)(sol->s, w->d->m);
+ nm_y = SCS(norm_inf)(sol->y, w->d->m);
+ sty = SCS(dot)(sol->s, sol->y, w->d->m);
+
info->setup_time = w->setup_time;
info->iter = iter;
info->res_infeas = w->r_orig->res_infeas;
info->res_unbdd_a = w->r_orig->res_unbdd_a;
info->res_unbdd_p = w->r_orig->res_unbdd_p;
info->scale = w->stgs->scale;
info->scale_updates = w->scale_updates;
info->rejected_accel_steps = w->rejected_accel_steps;
info->accepted_accel_steps = w->accepted_accel_steps;
- info->comp_slack = ABS(SCS(dot)(sol->s, sol->y, w->d->m));
- if (info->comp_slack > 1e-5 * MAX(SCS(norm_inf)(sol->s, w->d->m),
- SCS(norm_inf)(sol->y, w->d->m))) {
+ info->comp_slack = ABS(sty);
+ if (info->comp_slack > 1e-5 * MAX(nm_s, nm_y)) {
scs_printf("WARNING - large complementary slackness residual: %f\n",
info->comp_slack);
}
switch (info->status_val) {
case SCS_SOLVED:
@@ -561,11 +576,13 @@
scs_printf("%*.2e ", (int)HSPACE, r->res_dual);
scs_printf("%*.2e ", (int)HSPACE, r->gap);
/* report mid point of primal and dual objective values */
scs_printf("%*.2e ", (int)HSPACE, 0.5 * (r->pobj + r->dobj));
scs_printf("%*.2e ", (int)HSPACE, w->stgs->scale);
- scs_printf("%*.2e ", (int)HSPACE, SCS(tocq)(solve_timer) / 1e3);
+ /* Report TOTAL time, including setup */
+ scs_printf("%*.2e ", (int)HSPACE,
+ (SCS(tocq)(solve_timer) + w->setup_time) / 1e3);
scs_printf("\n");
#if VERBOSITY > 0
scs_printf("Norm u = %4f, ", SCS(norm_2)(w->u, w->d->n + w->d->m + 1));
scs_printf("Norm u_t = %4f, ", SCS(norm_2)(w->u_t, w->d->n + w->d->m + 1));
@@ -647,30 +664,33 @@
mexEvalString("drawnow;");
#endif
}
static scs_int has_converged(ScsWork *w, scs_int iter) {
+ scs_float abs_xt_p_x, abs_ctx, abs_bty;
+ scs_float nm_s, nm_px, nm_aty, nm_ax;
+ scs_float grl, prl, drl;
scs_float eps_abs = w->stgs->eps_abs;
scs_float eps_rel = w->stgs->eps_rel;
scs_float eps_infeas = w->stgs->eps_infeas;
- scs_float grl, prl, drl;
ScsResiduals *r = w->r_orig;
- scs_float *b = w->b_orig;
- scs_float *c = w->c_orig;
- scs_float *s = w->xys_orig->s;
if (r->tau > 0.) {
+ abs_xt_p_x = ABS(r->xt_p_x);
+ abs_ctx = ABS(r->ctx);
+ abs_bty = ABS(r->bty);
+
+ nm_s = NORM(w->xys_orig->s, w->d->m);
+ nm_px = NORM(r->px, w->d->n);
+ nm_aty = NORM(r->aty, w->d->n);
+ nm_ax = NORM(r->ax, w->d->m);
/* xt_p_x, ctx, bty already have tau divided out */
- grl = MAX(MAX(ABS(r->xt_p_x), ABS(r->ctx)), ABS(r->bty));
+ grl = MAX(MAX(abs_xt_p_x, abs_ctx), abs_bty);
/* s, ax, px, aty do *not* have tau divided out, so need to divide */
- prl = MAX(MAX(NORM(b, w->d->m) * r->tau, NORM(s, w->d->m)),
- NORM(r->ax, w->d->m)) /
- r->tau;
- drl = MAX(MAX(NORM(c, w->d->n) * r->tau, NORM(r->px, w->d->n)),
- NORM(r->aty, w->d->n)) /
- r->tau;
+ prl = MAX(MAX(w->nm_b_orig * r->tau, nm_s), nm_ax) / r->tau;
+ drl = MAX(MAX(w->nm_c_orig * r->tau, nm_px), nm_aty) / r->tau;
if (isless(r->res_pri, eps_abs + eps_rel * prl) &&
isless(r->res_dual, eps_abs + eps_rel * drl) &&
isless(r->gap, eps_abs + eps_rel * grl)) {
return SCS_SOLVED;
}
@@ -683,11 +703,11 @@
return SCS_INFEASIBLE;
}
return 0;
}
-#if NOVALIDATE == 0
+#if NO_VALIDATE == 0
static scs_int validate(const ScsData *d, const ScsCone *k,
const ScsSettings *stgs) {
if (d->m <= 0 || d->n <= 0) {
scs_printf("m and n must both be greater than 0; m = %li, n = %li\n",
(long)d->m, (long)d->n);
@@ -760,16 +780,19 @@
memcpy(w->b_orig, b, w->d->m * sizeof(scs_float));
memcpy(w->d->b, b, w->d->m * sizeof(scs_float));
} else {
memcpy(w->d->b, w->b_orig, w->d->m * sizeof(scs_float));
}
+ w->nm_b_orig = NORM(w->b_orig, w->d->m);
+
if (c) {
memcpy(w->c_orig, c, w->d->n * sizeof(scs_float));
memcpy(w->d->c, c, w->d->n * sizeof(scs_float));
} else {
memcpy(w->d->c, w->c_orig, w->d->n * sizeof(scs_float));
}
+ w->nm_c_orig = NORM(w->c_orig, w->d->n);
/* normalize */
if (w->scal) {
SCS(normalize_b_c)(w->scal, w->d->b, w->d->c);
}
@@ -824,11 +847,11 @@
w->v = (scs_float *)scs_calloc(l, sizeof(scs_float));
w->v_prev = (scs_float *)scs_calloc(l, sizeof(scs_float));
w->rsk = (scs_float *)scs_calloc(l, sizeof(scs_float));
w->h = (scs_float *)scs_calloc((l - 1), sizeof(scs_float));
w->g = (scs_float *)scs_calloc((l - 1), sizeof(scs_float));
- w->lin_sys_warm_start = (scs_float *)scs_calloc((l - 1), sizeof(scs_float));
+ w->lin_sys_warm_start = (scs_float *)scs_calloc(w->d->n, sizeof(scs_float));
w->diag_r = (scs_float *)scs_calloc(l, sizeof(scs_float));
/* x,y,s struct */
w->xys_orig = (ScsSolution *)scs_calloc(1, sizeof(ScsSolution));
w->xys_orig->x = (scs_float *)scs_calloc(w->d->n, sizeof(scs_float));
w->xys_orig->s = (scs_float *)scs_calloc(w->d->m, sizeof(scs_float));
@@ -862,11 +885,11 @@
w->scal = SCS_NULL;
}
/* set w->*_orig and performs normalization if appropriate */
scs_update(w, w->d->b, w->d->c);
- if (!(w->p = SCS(init_lin_sys_work)(w->d->A, w->d->P, w->diag_r))) {
+ if (!(w->p = scs_init_lin_sys_work(w->d->A, w->d->P, w->diag_r))) {
scs_printf("ERROR: init_lin_sys_work failure\n");
return SCS_NULL;
}
if (w->stgs->acceleration_lookback) {
/* TODO(HACK!) negative acceleration_lookback interpreted as type-II */
@@ -889,11 +912,11 @@
static void update_work_cache(ScsWork *w) {
/* g = (I + M)^{-1} h */
memcpy(w->g, w->h, (w->d->n + w->d->m) * sizeof(scs_float));
SCS(scale_array)(&(w->g[w->d->n]), -1., w->d->m);
- SCS(solve_lin_sys)(w->p, w->g, SCS_NULL, CG_BEST_TOL);
+ scs_solve_lin_sys(w->p, w->g, SCS_NULL, CG_BEST_TOL);
return;
}
/* Reset quantities specific to current solve */
static void reset_tracking(ScsWork *w) {
@@ -932,28 +955,31 @@
return (factor > SQRTF(10.) || factor < 1. / SQRTF(10.));
}
static void update_scale(ScsWork *w, const ScsCone *k, scs_int iter) {
scs_int i;
- scs_float factor, new_scale;
+ scs_float factor, new_scale, relative_res_pri, relative_res_dual;
+ scs_float denom_pri, denom_dual;
ScsResiduals *r = w->r_orig;
- ScsSolution *xys = w->xys_orig;
- scs_float *b = w->b_orig;
- scs_float *c = w->c_orig;
+ scs_float nm_ax = SCALE_NORM(r->ax, w->d->m);
+ scs_float nm_s = SCALE_NORM(w->xys_orig->s, w->d->m);
+ scs_float nm_px_aty_ctau = SCALE_NORM(r->px_aty_ctau, w->d->n);
+ scs_float nm_px = SCALE_NORM(r->px, w->d->n);
+ scs_float nm_aty = SCALE_NORM(r->aty, w->d->n);
+ scs_float nm_ax_s_btau = SCALE_NORM(r->ax_s_btau, w->d->m);
+
scs_int iters_since_last_update = iter - w->last_scale_update_iter;
/* ||Ax + s - b * tau|| */
- scs_float relative_res_pri = SAFEDIV_POS(
- SCALE_NORM(r->ax_s_btau, w->d->m),
- MAX(MAX(SCALE_NORM(r->ax, w->d->m), SCALE_NORM(xys->s, w->d->m)),
- SCALE_NORM(b, w->d->m) * r->tau));
+ denom_pri = MAX(nm_ax, nm_s);
+ denom_pri = MAX(denom_pri, w->nm_b_orig * r->tau);
+ relative_res_pri = SAFEDIV_POS(nm_ax_s_btau, denom_pri);
/* ||Px + A'y + c * tau|| */
- scs_float relative_res_dual = SAFEDIV_POS(
- SCALE_NORM(r->px_aty_ctau, w->d->n),
- MAX(MAX(SCALE_NORM(r->px, w->d->n), SCALE_NORM(r->aty, w->d->n)),
- SCALE_NORM(c, w->d->n) * r->tau));
+ denom_dual = MAX(nm_px, nm_aty);
+ denom_dual = MAX(denom_dual, w->nm_c_orig * r->tau);
+ relative_res_dual = SAFEDIV_POS(nm_px_aty_ctau, denom_dual);
/* higher scale makes res_pri go down faster, so increase if res_pri larger */
w->sum_log_scale_factor += log(relative_res_pri) - log(relative_res_dual);
w->n_log_scale_factor++;
@@ -979,11 +1005,11 @@
/* update diag r vector */
set_diag_r(w);
/* update linear systems */
- SCS(update_lin_sys_diag_r)(w->p, w->diag_r);
+ scs_update_lin_sys_diag_r(w->p, w->diag_r);
/* update pre-solved quantities */
update_work_cache(w);
/* reset acceleration so that old iterates aren't affecting new values */
@@ -1023,11 +1049,11 @@
stgs->warm_start = warm_start;
/* initialize ctrl-c support */
scs_start_interrupt_listener();
SCS(tic)(&solve_timer);
- strcpy(info->lin_sys_solver, SCS(get_lin_sys_method)());
+ strcpy(info->lin_sys_solver, scs_get_lin_sys_method());
info->status_val = SCS_UNFINISHED; /* not yet converged */
update_work(w, sol);
if (w->stgs->verbose) {
print_header(w, k);
@@ -1055,11 +1081,11 @@
/* store v_prev = v, *after* normalizing */
memcpy(w->v_prev, w->v, l * sizeof(scs_float));
/******************* linear system solve ********************/
SCS(tic)(&lin_sys_timer);
- if (project_lin_sys(w, i) < 0) {
+ if (project_lin_sys(w, i) != 0) {
return failure(w, w->d->m, w->d->n, sol, info, SCS_FAILED,
"error in project_lin_sys", "failure");
}
total_lin_sys_time += SCS(tocq)(&lin_sys_timer);
@@ -1160,11 +1186,11 @@
void scs_finish(ScsWork *w) {
if (w) {
SCS(finish_cone)(w->cone_work);
if (w->p) {
- SCS(free_lin_sys_work)(w->p);
+ scs_free_lin_sys_work(w->p);
}
if (w->accel) {
aa_finish(w->accel);
}
free_work(w);
@@ -1177,11 +1203,11 @@
scs_start_interrupt_listener();
if (!d || !k) {
scs_printf("ERROR: Missing ScsData or ScsCone input\n");
return SCS_NULL;
}
-#if NOVALIDATE == 0
+#if NO_VALIDATE == 0
if (validate(d, k, stgs) < 0) {
scs_printf("ERROR: Validation returned failure\n");
return SCS_NULL;
}
#endif
@@ -1189,10 +1215,15 @@
scs_printf("size of scs_int = %lu, size of scs_float = %lu\n",
sizeof(scs_int), sizeof(scs_float));
#endif
SCS(tic)(&init_timer);
if (stgs->write_data_filename) {
+ scs_printf("Writing raw problem data to %s\n", stgs->write_data_filename);
SCS(write_data)(d, k, stgs);
+ }
+ if (stgs->log_csv_filename) {
+ scs_printf("Logging run data to %s\n", stgs->log_csv_filename);
+ /* logging done every iteration */
}
w = init_work(d, k, stgs);
if (w) {
w->setup_time = SCS(tocq)(&init_timer);
}