From d43ea2e3d16cb6bce051f6c8a752c27b29b738e9 Mon Sep 17 00:00:00 2001 From: Chris Matrakidis Date: Mon, 4 Apr 2016 19:38:28 +0300 Subject: [PATCH 1/2] keep track of perturbed coefficients in dual simplex (conservative) --- src/simplex/spydual.c | 51 +++++++++++++++++++++++++++++++-------------------- 1 file changed, 31 insertions(+), 20 deletions(-) diff --git a/src/simplex/spydual.c b/src/simplex/spydual.c index e41baf2..0b7ad07 100644 --- a/src/simplex/spydual.c +++ b/src/simplex/spydual.c @@ -902,9 +902,15 @@ done: /* the choice has been made */ * objective coefficient cN[j] to provide d[j] = 0, i.e. to make all * d[j] strongly feasible. Otherwise, if d[j] has a feasible value, the * routine attempts to reduce (or remove) perturbation in cN[j] by -* shifting d[j] to its zero bound keeping strong feasibility. */ +* shifting d[j] to its zero bound keeping strong feasibility. +* +* If no objective coefficient are modified, the routine returns zero. +* Otherwise the routine returns an upper bound of the number of +* modifications away from the original value minus the number of +* coefficients restored to their original bound. The sum of the +* returned values will always be positive or zero. */ -static void play_coef(struct csa *csa, int all) +static int play_coef(struct csa *csa, int all) { SPXLP *lp = csa->lp; int m = lp->m; int n = lp->n; @@ -922,6 +928,7 @@ static void play_coef(struct csa *csa, int all) #endif int j, k; double temp; + int ret = 0; /* reduced costs d = (d[j]) should be valid */ xassert(csa->d_st); /* walk thru the list of non-basic variables xN = (xN[j]) */ @@ -937,17 +944,19 @@ static void play_coef(struct csa *csa, int all) { /* xN[j] is free variable */ /* strong feasibility means d[j] = 0 */ c[k] -= d[j], d[j] = 0.0; + if (c[k] != orig_c[k]) + ret++; } else if (!flag[j]) { /* xN[j] has its lower bound active */ /* strong feasibility means d[j] >= 0 */ if (d[j] < 0.0) - c[k] -= d[j], d[j] = 0.0; + c[k] -= d[j], d[j] = 0.0, ret++; else if (c[k] > orig_c[k]) { /* remove/reduce perturbation of cN[j] */ temp = c[k] - orig_c[k]; /* > 0 */ if (temp < d[j]) - c[k] = orig_c[k], d[j] -= temp; + c[k] = orig_c[k], d[j] -= temp, ret--; else c[k] -= d[j], d[j] = 0.0; } @@ -956,19 +965,19 @@ static void play_coef(struct csa *csa, int all) { /* xN[j] has its upper bound active */ /* strong feasibility means d[j] <= 0 */ if (d[j] > 0.0) - c[k] -= d[j], d[j] = 0.0; + c[k] -= d[j], d[j] = 0.0, ret++; else if (c[k] < orig_c[k]) { /* remove/reduce perturbation of cN[j] */ temp = c[k] - orig_c[k]; /* < 0 */ if (temp > d[j]) - c[k] = orig_c[k], d[j] -= temp; + c[k] = orig_c[k], d[j] -= temp, ret--; else c[k] -= d[j], d[j] = 0.0; } } } } - return; + return ret; } #endif @@ -1291,7 +1300,7 @@ loop: /* main loop starts here */ #if 0 /* 29/III-2016 */ j = check_feas(csa, tol_dj, tol_dj1, 0); #else - if (perturb) + if (perturb > 1) j = check_feas(csa, 3.14 * tol_dj, 3.14 * tol_dj1, 0); else j = check_feas(csa, tol_dj, tol_dj1, 0); @@ -1336,7 +1345,7 @@ loop: /* main loop starts here */ /* check_feas guarantees that all d[j] are feasible within * a tolerance */ if (perturb) - play_coef(csa, 1); + perturb += play_coef(csa, 1); #endif } /* at this point the search phase is determined */ @@ -1376,7 +1385,7 @@ loop: /* main loop starts here */ /* check if the objective limit has been reached */ #if PERTURB /* FIXME */ - if (!perturb) + if (perturb <= 1) #endif if (csa->phase == 2 && csa->obj_lim != DBL_MAX && spx_eval_obj(lp, beta) >= csa->obj_lim) @@ -1405,12 +1414,12 @@ loop: /* main loop starts here */ if (csa->it_cnt - csa->it_beg >= csa->it_lim) { #if PERTURB - if (perturb) + if (perturb > 1) { /* remove perturbation */ - perturb = 0; memcpy(csa->lp->c, csa->orig_c, (1+n) * sizeof(double)); csa->phase = csa->d_st = 0; } + perturb = 0; #endif if (csa->beta_st != 1) csa->beta_st = 0; @@ -1441,12 +1450,12 @@ loop: /* main loop starts here */ if (1000.0 * xdifftime(xtime(), csa->tm_beg) >= csa->tm_lim) { #if PERTURB - if (perturb) + if (perturb > 1) { /* remove perturbation */ - perturb = 0; memcpy(csa->lp->c, csa->orig_c, (1+n) * sizeof(double)); csa->phase = csa->d_st = 0; } + perturb = 0; #endif if (csa->beta_st != 1) csa->beta_st = 0; @@ -1505,10 +1514,12 @@ loop: /* main loop starts here */ { #if PERTURB if (perturb && csa->phase == 2) - { /* remove perturbation */ + { if (perturb > 1) + { /* remove perturbation */ + memcpy(csa->lp->c, csa->orig_c, (1+n) * sizeof(double)); + csa->phase = csa->d_st = 0; + } perturb = 0; - memcpy(csa->lp->c, csa->orig_c, (1+n) * sizeof(double)); - csa->phase = csa->d_st = 0; } #endif if (csa->beta_st != 1) @@ -1582,12 +1593,12 @@ loop: /* main loop starts here */ if (csa->q == 0) { #if PERTURB - if (perturb) + if (perturb > 1) { /* remove perturbation */ - perturb = 0; memcpy(csa->lp->c, csa->orig_c, (1+n) * sizeof(double)); csa->phase = csa->d_st = 0; } + perturb = 0; #endif if (csa->beta_st != 1) csa->beta_st = 0; @@ -1747,7 +1758,7 @@ loop: /* main loop starts here */ #endif #if PERTURB if (perturb && csa->d_st) - play_coef(csa, 0); + perturb += play_coef(csa, 0); #endif /* dual simplex iteration complete */ csa->it_cnt++; -- 1.8.1.msysgit.1