vendor/scs/src/scs.c in scs-0.3.1 vs vendor/scs/src/scs.c in scs-0.3.2

- old
+ new

@@ -5,10 +5,11 @@ #include "linalg.h" #include "linsys.h" #include "normalize.h" #include "rw.h" #include "scs_matrix.h" +#include "scs_work.h" #include "util.h" /* printing header */ static const char *HEADER[] = { " iter ", " pri res ", " dua res ", " gap ", @@ -39,23 +40,20 @@ scs_free(w->rsk); scs_free(w->h); scs_free(w->g); scs_free(w->b_normalized); scs_free(w->c_normalized); - scs_free(w->rho_y_vec); scs_free(w->lin_sys_warm_start); - if (w->cone_boundaries) { - scs_free(w->cone_boundaries); - } + scs_free(w->diag_r); if (w->scal) { scs_free(w->scal->D); scs_free(w->scal->E); scs_free(w->scal); } SCS(free_sol)(w->xys_orig); free_residuals(w->r_orig); - if (w->stgs->normalize) { + if (w->stgs && w->stgs->normalize) { SCS(free_sol)(w->xys_normalized); free_residuals(w->r_normalized); } scs_free(w); } @@ -76,11 +74,11 @@ for (i = 0; i < LINE_LEN; ++i) { scs_printf("-"); } scs_printf("\n\t SCS v%s - Splitting Conic Solver\n\t(c) Brendan " "O'Donoghue, Stanford University, 2012\n", - SCS(version)()); + scs_version()); for (i = 0; i < LINE_LEN; ++i) { scs_printf("-"); } scs_printf("\n"); scs_printf("problem: variables n: %i, constraints m: %i\n", (int)d->n, @@ -155,27 +153,26 @@ scs_end_interrupt_listener(); return status; } /* given x,y,s warm start, set v = [x; s / R + y; 1] - * where R = diag(w->rho_y_vec). */ static void warm_start_vars(ScsWork *w, ScsSolution *sol) { scs_int n = w->n, m = w->m, i; scs_float *v = w->v; /* normalize the warm-start */ if (w->stgs->normalize) { - SCS(normalize_sol)(w, sol); + SCS(normalize_sol)(w->scal, sol); } memcpy(v, sol->x, n * sizeof(scs_float)); for (i = 0; i < m; ++i) { - v[i + n] = sol->y[i] + sol->s[i] / w->rho_y_vec[i]; + v[i + n] = sol->y[i] + sol->s[i] / w->diag_r[i + n]; } v[n + m] = 1.0; /* tau = 1 */ /* un-normalize so sol unchanged */ if (w->stgs->normalize) { - SCS(un_normalize_sol)(w, sol); + 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); @@ -219,16 +216,16 @@ r->bty = r_n->bty / pd; r->pobj = r_n->pobj / pd; r->dobj = r_n->dobj / pd; r->gap = r_n->gap / pd; - SCS(un_normalize_primal)(w, r->ax); - SCS(un_normalize_primal)(w, r->ax_s); - SCS(un_normalize_primal)(w, r->ax_s_btau); - SCS(un_normalize_dual)(w, r->aty); - SCS(un_normalize_dual)(w, r->px); - SCS(un_normalize_dual)(w, r->px_aty_ctau); + SCS(un_normalize_primal)(w->scal, r->ax); + SCS(un_normalize_primal)(w->scal, r->ax_s); + SCS(un_normalize_primal)(w->scal, r->ax_s_btau); + SCS(un_normalize_dual)(w->scal, r->aty); + SCS(un_normalize_dual)(w->scal, r->px); + SCS(un_normalize_dual)(w->scal, r->px_aty_ctau); compute_residuals(r, w->m, w->n); } /* calculates un-normalized residual quantities */ @@ -304,62 +301,49 @@ if (w->stgs->normalize) { memcpy(w->xys_orig->x, w->xys_normalized->x, n * sizeof(scs_float)); memcpy(w->xys_orig->y, w->xys_normalized->y, m * sizeof(scs_float)); memcpy(w->xys_orig->s, w->xys_normalized->s, m * sizeof(scs_float)); - SCS(un_normalize_sol)(w, w->xys_orig); + SCS(un_normalize_sol)(w->scal, w->xys_orig); unnormalize_residuals(w); } } static void cold_start_vars(ScsWork *w) { scs_int l = w->n + w->m + 1; memset(w->v, 0, l * sizeof(scs_float)); w->v[l - 1] = 1.; } -/* utility function that scales first n entries in inner prod by rho_x */ -/* and last m entries by 1 / rho_y_vec, assumes length of array is n + m */ -/* See .note_on_scale in repo for explanation */ -static scs_float dot_with_diag_scaling(ScsWork *w, const scs_float *x, - const scs_float *y) { - scs_int i, n = w->n, len = w->n + w->m; +/* utility function that computes x'Ry */ +static inline scs_float dot_r(ScsWork *w, const scs_float *x, + const scs_float *y) { + scs_int i; scs_float ip = 0.0; - for (i = 0; i < n; ++i) { - ip += w->stgs->rho_x * x[i] * y[i]; + for (i = 0; i < w->n + w->m; ++i) { + ip += x[i] * y[i] * w->diag_r[i]; } - for (i = n; i < len; ++i) { - ip += x[i] * y[i] * w->rho_y_vec[i - n]; - } return ip; } -static inline scs_float get_tau_scale(ScsWork *w) { - return TAU_FACTOR; /* TAU_FACTOR * w->scale; */ -} - static scs_float root_plus(ScsWork *w, scs_float *p, scs_float *mu, scs_float eta) { - scs_float b, c, tau, a, tau_scale; - tau_scale = get_tau_scale(w); - a = tau_scale + dot_with_diag_scaling(w, w->g, w->g); - b = (dot_with_diag_scaling(w, mu, w->g) - - 2 * dot_with_diag_scaling(w, p, w->g) - eta * tau_scale); - c = dot_with_diag_scaling(w, p, p) - dot_with_diag_scaling(w, p, mu); - tau = (-b + SQRTF(MAX(b * b - 4 * a * c, 0.))) / (2 * a); - return tau; + scs_float a, b, c, tau_scale = w->diag_r[w->n + w->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); } /* status < 0 indicates failure */ static scs_int project_lin_sys(ScsWork *w, scs_int iter) { scs_int n = w->n, m = w->m, l = n + m + 1, status, i; scs_float *warm_start = SCS_NULL; - scs_float tol = -1.0; /* only used for indirect methods, overriden later */ + scs_float tol = -1.0; /* only used for indirect methods, overridden later */ memcpy(w->u_t, w->v, l * sizeof(scs_float)); - SCS(scale_array)(w->u_t, w->stgs->rho_x, n); - for (i = n; i < l - 1; ++i) { - w->u_t[i] *= -w->rho_y_vec[i - n]; + for (i = 0; i < l - 1; ++i) { + w->u_t[i] *= (i < n ? 1 : -1) * w->diag_r[i]; } #if INDIRECT > 0 /* compute warm start using the cone projection output */ warm_start = w->lin_sys_warm_start; memcpy(warm_start, w->u, (l - 1) * sizeof(scs_float)); @@ -393,43 +377,31 @@ * since it depends on v^k MUST be called before update_dual_vars is done * (no effect of w->stgs->alpha here). */ static void compute_rsk(ScsWork *w) { scs_int i, l = w->m + w->n + 1; - /* r, should = 0 so skip */ - /* - for (i = 0; i < w->n; ++i) { - w->rsk[i] = w->stgs->rho_x * (w->v[i] + w->u[i] - 2 * w->u_t[i]); + for (i = 0; i < l; ++i) { + w->rsk[i] = (w->v[i] + w->u[i] - 2 * w->u_t[i]) * w->diag_r[i]; } - */ - /* s */ - for (i = w->n; i < l - 1; ++i) { - w->rsk[i] = (w->v[i] + w->u[i] - 2 * w->u_t[i]) * w->rho_y_vec[i - w->n]; - } - /* kappa, incorporates tau scaling parameter */ - w->rsk[l - 1] = - get_tau_scale(w) * (w->v[l - 1] + w->u[l - 1] - 2 * w->u_t[l - 1]); } static void update_dual_vars(ScsWork *w) { scs_int i, l = w->n + w->m + 1; - scs_float a = w->stgs->alpha; - /* compute and store [r;s;kappa] */ - compute_rsk(w); for (i = 0; i < l; ++i) { - w->v[i] += a * (w->u[i] - w->u_t[i]); + w->v[i] += w->stgs->alpha * (w->u[i] - w->u_t[i]); } } /* status < 0 indicates failure */ static scs_int project_cones(ScsWork *w, const ScsCone *k, scs_int iter) { scs_int i, n = w->n, l = w->n + w->m + 1, status; for (i = 0; i < l; ++i) { w->u[i] = 2 * w->u_t[i] - w->v[i]; } /* u = [x;y;tau] */ - status = SCS(proj_dual_cone)(&(w->u[n]), k, w->cone_work, w->stgs->normalize); + status = + SCS(proj_dual_cone)(&(w->u[n]), w->cone_work, w->scal, &(w->diag_r[n])); if (iter < FEASIBLE_ITERS) { w->u[l - 1] = 1.0; } else { w->u[l - 1] = MAX(w->u[l - 1], 0.); } @@ -524,11 +496,11 @@ scs_int iter) { setx(w, sol); sety(w, sol); sets(w, sol); if (w->stgs->normalize) { - SCS(un_normalize_sol)(w, sol); + SCS(un_normalize_sol)(w->scal, sol); } populate_residual_struct(w, iter); info->setup_time = w->setup_time; info->iter = iter; info->res_infeas = w->r_orig->res_infeas; @@ -576,10 +548,11 @@ #if VERBOSITY > 0 scs_printf("Norm u = %4f, ", SCS(norm_2)(w->u, w->n + w->m + 1)); scs_printf("Norm u_t = %4f, ", SCS(norm_2)(w->u_t, w->n + w->m + 1)); scs_printf("Norm v = %4f, ", SCS(norm_2)(w->v, w->n + w->m + 1)); + scs_printf("Norm rsk = %4f, ", SCS(norm_2)(w->rsk, w->n + w->m + 1)); scs_printf("Norm x = %4f, ", SCS(norm_2)(w->xys_orig->x, w->n)); scs_printf("Norm y = %4f, ", SCS(norm_2)(w->xys_orig->y, w->m)); scs_printf("Norm s = %4f, ", SCS(norm_2)(w->xys_orig->s, w->m)); scs_printf("Norm |Ax + s| = %1.2e, ", SCS(norm_2)(r->ax_s, w->m)); scs_printf("tau = %4f, ", w->u[w->n + w->m]); @@ -737,10 +710,14 @@ } if (stgs->scale <= 0) { scs_printf("scale must be positive (1 works well).\n"); return -1; } + if (stgs->acceleration_interval <= 0) { + scs_printf("acceleration_interval must be positive (10 works well).\n"); + return -1; + } return 0; } #endif static ScsResiduals *init_residuals(const ScsData *d) { @@ -753,10 +730,22 @@ r->aty = (scs_float *)scs_calloc(d->n, sizeof(scs_float)); r->px_aty_ctau = (scs_float *)scs_calloc(d->n, sizeof(scs_float)); return r; } +/* Sets the diag_r vector, given the scale parameters in work */ +static void set_diag_r(ScsWork *w) { + scs_int i; + for (i = 0; i < w->n; ++i) { + w->diag_r[i] = w->stgs->rho_x; + } + /* use cone information to set R_y */ + SCS(set_r_y)(w->cone_work, w->scale, &(w->diag_r[w->n])); + /* if modified need to SCS(enforce_cone_boundaries)(...) */ + w->diag_r[w->n + w->m] = TAU_FACTOR; /* TODO: is this the best choice? */ +} + static ScsWork *init_work(const ScsData *d, const ScsCone *k, const ScsSettings *stgs) { ScsWork *w = (ScsWork *)scs_calloc(1, sizeof(ScsWork)); scs_int l = d->n + d->m + 1; if (stgs->verbose) { @@ -785,11 +774,11 @@ 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->rho_y_vec = (scs_float *)scs_calloc(d->m, 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(d->n, sizeof(scs_float)); w->xys_orig->s = (scs_float *)scs_calloc(d->m, sizeof(scs_float)); w->xys_orig->y = (scs_float *)scs_calloc(d->m, sizeof(scs_float)); @@ -802,12 +791,17 @@ w->c_orig = d->c; w->b_normalized = (scs_float *)scs_calloc(d->m, sizeof(scs_float)); w->c_normalized = (scs_float *)scs_calloc(d->n, sizeof(scs_float)); memcpy(w->b_normalized, w->b_orig, w->m * sizeof(scs_float)); memcpy(w->c_normalized, w->c_orig, w->n * sizeof(scs_float)); - SCS(set_rho_y_vec)(k, w->scale, w->rho_y_vec, w->m); + if (!(w->cone_work = SCS(init_cone)(k, w->m))) { + scs_printf("ERROR: init_cone failure\n"); + return SCS_NULL; + } + set_diag_r(w); + if (!w->c_normalized) { scs_printf("ERROR: work memory allocation failure\n"); return SCS_NULL; } @@ -827,28 +821,18 @@ scs_printf("ERROR: copy P matrix failed\n"); return SCS_NULL; } #endif /* this allocates memory that must be freed */ - w->cone_boundaries_len = SCS(set_cone_boundaries)(k, &w->cone_boundaries); - w->scal = (ScsScaling *)scs_calloc(1, sizeof(ScsScaling)); - SCS(normalize) - (w->P, w->A, w->b_normalized, w->c_normalized, w->scal, w->cone_boundaries, - w->cone_boundaries_len); + w->scal = SCS(normalize_a_p)(w->P, w->A, w->b_normalized, w->c_normalized, + w->cone_work); } else { w->xys_normalized = w->xys_orig; w->r_normalized = w->r_orig; - w->cone_boundaries_len = 0; - w->cone_boundaries = SCS_NULL; w->scal = SCS_NULL; } - if (!(w->cone_work = SCS(init_cone)(k, w->scal, w->m))) { - scs_printf("ERROR: init_cone failure\n"); - return SCS_NULL; - } - if (!(w->p = - SCS(init_lin_sys_work)(w->A, w->P, w->rho_y_vec, w->stgs->rho_x))) { + if (!(w->p = SCS(init_lin_sys_work)(w->A, w->P, w->diag_r))) { scs_printf("ERROR: init_lin_sys_work failure\n"); return SCS_NULL; } /* Acceleration */ w->rejected_accel_steps = 0; @@ -896,15 +880,15 @@ update_work_cache(w); return 0; } /* will update if the factor is outside of range */ -scs_int should_update_rho_y_vec(scs_float factor, scs_int iter) { +scs_int should_update_r(scs_float factor, scs_int iter) { return (factor > SQRTF(10.) || factor < 1. / SQRTF(10.)); } -static void maybe_update_scale(ScsWork *w, const ScsCone *k, scs_int iter) { +static void update_scale(ScsWork *w, const ScsCone *k, scs_int iter) { scs_int i; scs_float factor, new_scale; ScsResiduals *r = w->r_orig; ScsSolution *xys = w->xys_orig; @@ -937,42 +921,47 @@ } new_scale = MIN(MAX(w->scale * factor, MIN_SCALE_VALUE), MAX_SCALE_VALUE); if (new_scale == w->scale) { return; } - if (should_update_rho_y_vec(factor, iters_since_last_update)) { + if (should_update_r(factor, iters_since_last_update)) { w->scale_updates++; w->sum_log_scale_factor = 0; w->n_log_scale_factor = 0; w->last_scale_update_iter = iter; w->scale = new_scale; - SCS(set_rho_y_vec)(k, w->scale, w->rho_y_vec, w->m); - SCS(update_lin_sys_rho_y_vec)(w->p, w->rho_y_vec); + /* update diag r vector */ + set_diag_r(w); + + /* update linear systems */ + 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 */ if (w->accel) { aa_reset(w->accel); } /* update v, using fact that rsk, u, u_t vectors should be the same */ - /* solve: R (v^+ + u - 2u_t) = rsk => v^+ = R^-1 rsk + 2u_t - u */ - /* only elements with scale on diag */ - for (i = w->n; i < w->n + w->m; i++) { - w->v[i] = w->rsk[i] / w->rho_y_vec[i - w->n] + 2 * w->u_t[i] - w->u[i]; + /* solve: R^+ (v^+ + u - 2u_t) = rsk = R(v + u - 2u_t) + * => v^+ = R+^-1 rsk + 2u_t - u + */ + for (i = 0; i < w->n + w->m + 1; i++) { + w->v[i] = w->rsk[i] / w->diag_r[i] + 2 * w->u_t[i] - w->u[i]; } } } /* scs is homogeneous so scale the iterate to keep norm reasonable */ static inline void normalize_v(scs_float *v, scs_int len) { scs_float v_norm = SCS(norm_2)(v, len); /* always l2 norm */ SCS(scale_array)(v, SQRTF((scs_float)len) * ITERATE_NORM / v_norm, len); } -scs_int SCS(solve)(ScsWork *w, ScsSolution *sol, ScsInfo *info) { +scs_int scs_solve(ScsWork *w, ScsSolution *sol, ScsInfo *info) { scs_int i; SCS(timer) solve_timer, lin_sys_timer, cone_timer, accel_timer; scs_float total_accel_time = 0.0, total_cone_time = 0.0, total_lin_sys_time = 0.0; if (!sol || !w || !info) { @@ -984,10 +973,11 @@ const ScsCone *k = w->k; const ScsSettings *stgs = w->stgs; /* initialize ctrl-c support */ scs_start_interrupt_listener(); SCS(tic)(&solve_timer); + strcpy(info->lin_sys_solver, SCS(get_lin_sys_method)()); info->status_val = SCS_UNFINISHED; /* not yet converged */ update_work(d, w, sol); if (w->stgs->verbose) { print_header(w, k); @@ -1013,28 +1003,29 @@ } /* store v_prev = v, *after* normalizing */ memcpy(w->v_prev, w->v, l * sizeof(scs_float)); - /* linear system solve */ + /******************* linear system solve ********************/ SCS(tic)(&lin_sys_timer); if (project_lin_sys(w, i) < 0) { return failure(w, w->m, w->n, sol, info, SCS_FAILED, "error in project_lin_sys", "failure"); } total_lin_sys_time += SCS(tocq)(&lin_sys_timer); - /* project onto the cones */ + /****************** project onto the cones ******************/ SCS(tic)(&cone_timer); if (project_cones(w, k, i) < 0) { return failure(w, w->m, w->n, sol, info, SCS_FAILED, "error in project_cones", "failure"); } total_cone_time += SCS(tocq)(&cone_timer); - /* dual variable step */ - update_dual_vars(w); + /* compute [r;s;kappa], must be before dual var update */ + /* since Moreau decomp logic relies on v at start */ + compute_rsk(w); if (i % CONVERGED_INTERVAL == 0) { if (scs_is_interrupted()) { return failure(w, w->m, w->n, sol, info, SCS_SIGINT, "interrupted", "interrupted"); @@ -1049,10 +1040,25 @@ break; } } } + /* Compute residuals. */ + if (w->stgs->verbose && i % PRINT_INTERVAL == 0) { + populate_residual_struct(w, i); + print_summary(w, i, &solve_timer); + } + + /* If residuals are fresh then maybe compute new scale. */ + if (w->stgs->adaptive_scale && i == w->r_orig->last_iter) { + update_scale(w, k, i); + } + + /****************** dual variable step **********************/ + /* do this after update_scale due to remapping that happens there */ + update_dual_vars(w); + /* AA safeguard check. * Perform safeguarding *after* convergence check to prevent safeguard * overwriting converged iterate, since safeguard is on `v` and convergence * is on `u`. */ @@ -1064,21 +1070,10 @@ } else { w->accepted_accel_steps++; } } - /* Compute residuals. */ - if (w->stgs->verbose && i % PRINT_INTERVAL == 0) { - populate_residual_struct(w, i); - print_summary(w, i, &solve_timer); - } - - /* If residuals are fresh then maybe compute new scale. */ - if (w->stgs->adaptive_scale && i == w->r_orig->last_iter) { - maybe_update_scale(w, k, i); - } - /* Log *after* updating scale so residual recalc does not affect alg */ if (w->stgs->log_csv_filename) { /* calc residuals every iter if logging to csv */ populate_residual_struct(w, i); SCS(log_data_to_csv)(d, k, stgs, w, i, &solve_timer); @@ -1111,16 +1106,16 @@ scs_end_interrupt_listener(); return info->status_val; } -void SCS(finish)(ScsWork *w) { +void scs_finish(ScsWork *w) { if (w) { SCS(finish_cone)(w->cone_work); if (w->stgs && w->stgs->normalize) { #ifndef COPYAMATRIX - SCS(un_normalize)(w->A, w->P, w->scal); + SCS(un_normalize_a_p)(w->A, w->P, w->scal); #else SCS(free_scs_matrix)(w->A); SCS(free_scs_matrix)(w->P); #endif } @@ -1132,12 +1127,11 @@ } free_work(w); } } -ScsWork *SCS(init)(const ScsData *d, const ScsCone *k, - const ScsSettings *stgs) { +ScsWork *scs_init(const ScsData *d, const ScsCone *k, const ScsSettings *stgs) { ScsWork *w; SCS(timer) init_timer; scs_start_interrupt_listener(); if (!d || !k) { scs_printf("ERROR: Missing ScsData or ScsCone input\n"); @@ -1159,24 +1153,24 @@ } scs_end_interrupt_listener(); return w; } -/* this just calls SCS(init), SCS(solve), and SCS(finish) */ +/* this just calls scs_init, scs_solve, and scs_finish */ scs_int scs(const ScsData *d, const ScsCone *k, const ScsSettings *stgs, ScsSolution *sol, ScsInfo *info) { scs_int status; - ScsWork *w = SCS(init)(d, k, stgs); + ScsWork *w = scs_init(d, k, stgs); #if VERBOSITY > 0 scs_printf("size of scs_int = %lu, size of scs_float = %lu\n", sizeof(scs_int), sizeof(scs_float)); #endif if (w) { - SCS(solve)(w, sol, info); + scs_solve(w, sol, info); status = info->status_val; } else { status = failure(SCS_NULL, d ? d->m : -1, d ? d->n : -1, sol, info, SCS_FAILED, "could not initialize work", "failure"); } - SCS(finish)(w); + scs_finish(w); return status; }