diff --git a/include/glbopts.h b/include/glbopts.h index 4d7b7fdf..ed2a15b5 100644 --- a/include/glbopts.h +++ b/include/glbopts.h @@ -39,7 +39,8 @@ extern "C" { #define WRITE_DATA_FILENAME (0) #define LOG_CSV_FILENAME (0) #define TIME_LIMIT_SECS (0.) - +/* Tolerance to check negativity condition for infeasibility */ +#define INFEAS_NEGATIVITY_TOL (1e-9) /* redefine printfs as needed */ #if NO_PRINTING > 0 /* Disable all printing */ #define scs_printf(...) /* No-op */ diff --git a/src/scs.c b/src/scs.c index 9cb1ee58..7dfd092a 100644 --- a/src/scs.c +++ b/src/scs.c @@ -207,13 +207,13 @@ static void compute_residuals(ScsResiduals *r, scs_int m, scs_int n) { r->res_unbdd_a = NAN; r->res_unbdd_p = NAN; r->res_infeas = NAN; - if (r->ctx_tau < 0) { + if (r->ctx_tau < -INFEAS_NEGATIVITY_TOL) { 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) { + if (r->bty_tau < -INFEAS_NEGATIVITY_TOL) { nm_aty = NORM(r->aty, n); r->res_infeas = SAFEDIV_POS(nm_aty, -r->bty_tau); } @@ -587,22 +587,22 @@ static void print_summary(ScsWork *w, scs_int i, SCS(timer) * solve_timer) { 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)); - scs_printf("Norm v = %4f, ", SCS(norm_2)(w->v, w->d->n + w->d->m + 1)); - scs_printf("Norm rsk = %4f, ", SCS(norm_2)(w->rsk, w->d->n + w->d->m + 1)); - scs_printf("Norm x = %4f, ", SCS(norm_2)(w->xys_orig->x, w->d->n)); - scs_printf("Norm y = %4f, ", SCS(norm_2)(w->xys_orig->y, w->d->m)); - scs_printf("Norm s = %4f, ", SCS(norm_2)(w->xys_orig->s, w->d->m)); - scs_printf("Norm |Ax + s| = %1.2e, ", SCS(norm_2)(r->ax_s, w->d->m)); - scs_printf("tau = %4f, ", w->u[w->d->n + w->d->m]); - scs_printf("kappa = %4f, ", w->rsk[w->d->n + w->d->m]); - scs_printf("|u - u_t| = %1.2e, ", + scs_printf("Norm u = %1.6e, ", SCS(norm_2)(w->u, w->d->n + w->d->m + 1)); + scs_printf("Norm u_t = %1.6e, ", SCS(norm_2)(w->u_t, w->d->n + w->d->m + 1)); + scs_printf("Norm v = %1.6e, ", SCS(norm_2)(w->v, w->d->n + w->d->m + 1)); + scs_printf("Norm rsk = %1.6e, ", SCS(norm_2)(w->rsk, w->d->n + w->d->m + 1)); + scs_printf("Norm x = %1.6e, ", SCS(norm_2)(w->xys_orig->x, w->d->n)); + scs_printf("Norm y = %1.6e, ", SCS(norm_2)(w->xys_orig->y, w->d->m)); + scs_printf("Norm s = %1.6e, ", SCS(norm_2)(w->xys_orig->s, w->d->m)); + scs_printf("Norm |Ax + s| = %1.6e, ", SCS(norm_2)(r->ax_s, w->d->m)); + scs_printf("tau = %1.6e, ", w->u[w->d->n + w->d->m]); + scs_printf("kappa = %1.6e, ", w->rsk[w->d->n + w->d->m]); + scs_printf("|u - u_t| = %1.6e, ", SCS(norm_diff)(w->u, w->u_t, w->d->n + w->d->m + 1)); - scs_printf("res_infeas = %1.2e, ", r->res_infeas); - scs_printf("res_unbdd_a = %1.2e, ", r->res_unbdd_a); - scs_printf("res_unbdd_p = %1.2e, ", r->res_unbdd_p); - scs_printf("ctx_tau = %1.2e, ", r->ctx_tau); + scs_printf("res_infeas = %1.6e, ", r->res_infeas); + scs_printf("res_unbdd_a = %1.6e, ", r->res_unbdd_a); + scs_printf("res_unbdd_p = %1.6e, ", r->res_unbdd_p); + scs_printf("ctx_tau = %1.6e, ", r->ctx_tau); scs_printf("bty_tau = %1.2e\n", r->bty_tau); #endif diff --git a/test/problems/mpc_bug.h b/test/problems/mpc_bug.h new file mode 100644 index 00000000..e8fc2e03 --- /dev/null +++ b/test/problems/mpc_bug.h @@ -0,0 +1,19 @@ +#include "glbopts.h" +#include "problems/test_prob_from_data_file.h" +#include "scs.h" + +static const char *mpc_bug(void) { + const char *fail; + scs_float OPT1 = -0.473957794500; /* from scs */ + scs_float OPT2 = -0.029336830816; /* from scs */ + scs_float OPT3 = -0.002215217478; /* from scs */ + fail = _test_prob_from_data("test/problems/mpc_bug1", OPT1); + if (fail) { + return fail; + } + fail = _test_prob_from_data("test/problems/mpc_bug2", OPT2); + if (fail) { + return fail; + } + return _test_prob_from_data("test/problems/mpc_bug3", OPT3); +} diff --git a/test/problems/mpc_bug1 b/test/problems/mpc_bug1 new file mode 100644 index 00000000..989fbb3e Binary files /dev/null and b/test/problems/mpc_bug1 differ diff --git a/test/problems/mpc_bug2 b/test/problems/mpc_bug2 new file mode 100644 index 00000000..4fc1ca6d Binary files /dev/null and b/test/problems/mpc_bug2 differ diff --git a/test/problems/mpc_bug3 b/test/problems/mpc_bug3 new file mode 100644 index 00000000..bfed77f2 Binary files /dev/null and b/test/problems/mpc_bug3 differ diff --git a/test/problems/test_prob_from_data_file.h b/test/problems/test_prob_from_data_file.h index 5e972bdb..3b993039 100644 --- a/test/problems/test_prob_from_data_file.h +++ b/test/problems/test_prob_from_data_file.h @@ -19,6 +19,8 @@ static const char *_test_prob_from_data(const char *file, scs_float OPT) { scs_float perr, derr; scs_int success; const char *fail; + scs_float xt_p_x; + scs_float *px = SCS_NULL; read_status = SCS(read_data)(file, &d, &k, &stgs); @@ -28,12 +30,24 @@ static const char *_test_prob_from_data(const char *file, scs_float OPT) { stgs->eps_abs = 1e-6; stgs->eps_rel = 1e-6; + /* Force verbosity for the test */ + stgs->verbose = 1; sol = (ScsSolution *)scs_calloc(1, sizeof(ScsSolution)); exitflag = scs(d, k, stgs, sol, &info); - perr = SCS(dot)(d->c, sol->x, d->n) - OPT; - derr = -SCS(dot)(d->b, sol->y, d->m) - OPT; + if (d->P) { + /* px = Px */ + px = (scs_float *)scs_calloc(d->n, sizeof(scs_float)); + memset(px, 0, d->n * sizeof(scs_float)); + SCS(accum_by_p)(d->P, sol->x, px); + xt_p_x = SCS(dot)(px, sol->x, d->n); + } else { + xt_p_x = 0.; + } + + perr = 0.5 * xt_p_x + SCS(dot)(d->c, sol->x, d->n) - OPT; + derr = -0.5 * xt_p_x - SCS(dot)(d->b, sol->y, d->m) - OPT; scs_printf("primal obj error %4e\n", perr); scs_printf("dual obj error %4e\n", derr); @@ -47,7 +61,9 @@ static const char *_test_prob_from_data(const char *file, scs_float OPT) { SCS(free_cone)(k); SCS(free_sol)(sol); scs_free(stgs); - + if (px) { + scs_free(px); + } if (fail) { scs_printf("%s: FAILED\n", file); } diff --git a/test/run_tests.c b/test/run_tests.c index 337b5d7a..15094ce7 100644 --- a/test/run_tests.c +++ b/test/run_tests.c @@ -43,9 +43,11 @@ _SKIP(random_prob) #if NO_READ_WRITE == 0 /* reads / writes */ #include "problems/hs21_tiny_qp_rw.h" #include "problems/max_ent.h" +#include "problems/mpc_bug.h" #else _SKIP(hs21_tiny_qp_rw) _SKIP(max_ent) +_SKIP(mpc_bug) #endif static const char *all_tests(void) { @@ -61,6 +63,7 @@ static const char *all_tests(void) { mu_run_test(unbounded_tiny_qp); mu_run_test(random_prob); mu_run_test(max_ent); + mu_run_test(mpc_bug); mu_run_test(test_exp_cone); return 0; }