1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpssx02.c	Sun Nov 06 20:59:10 2011 +0100
     1.3 @@ -0,0 +1,479 @@
     1.4 +/* glpssx02.c */
     1.5 +
     1.6 +/***********************************************************************
     1.7 +*  This code is part of GLPK (GNU Linear Programming Kit).
     1.8 +*
     1.9 +*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
    1.10 +*  2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
    1.11 +*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
    1.12 +*  E-mail: <mao@gnu.org>.
    1.13 +*
    1.14 +*  GLPK is free software: you can redistribute it and/or modify it
    1.15 +*  under the terms of the GNU General Public License as published by
    1.16 +*  the Free Software Foundation, either version 3 of the License, or
    1.17 +*  (at your option) any later version.
    1.18 +*
    1.19 +*  GLPK is distributed in the hope that it will be useful, but WITHOUT
    1.20 +*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    1.21 +*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
    1.22 +*  License for more details.
    1.23 +*
    1.24 +*  You should have received a copy of the GNU General Public License
    1.25 +*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
    1.26 +***********************************************************************/
    1.27 +
    1.28 +#include "glpenv.h"
    1.29 +#include "glpssx.h"
    1.30 +
    1.31 +static void show_progress(SSX *ssx, int phase)
    1.32 +{     /* this auxiliary routine displays information about progress of
    1.33 +         the search */
    1.34 +      int i, def = 0;
    1.35 +      for (i = 1; i <= ssx->m; i++)
    1.36 +         if (ssx->type[ssx->Q_col[i]] == SSX_FX) def++;
    1.37 +      xprintf("%s%6d:   %s = %22.15g   (%d)\n", phase == 1 ? " " : "*",
    1.38 +         ssx->it_cnt, phase == 1 ? "infsum" : "objval",
    1.39 +         mpq_get_d(ssx->bbar[0]), def);
    1.40 +#if 0
    1.41 +      ssx->tm_lag = utime();
    1.42 +#else
    1.43 +      ssx->tm_lag = xtime();
    1.44 +#endif
    1.45 +      return;
    1.46 +}
    1.47 +
    1.48 +/*----------------------------------------------------------------------
    1.49 +// ssx_phase_I - find primal feasible solution.
    1.50 +//
    1.51 +// This routine implements phase I of the primal simplex method.
    1.52 +//
    1.53 +// On exit the routine returns one of the following codes:
    1.54 +//
    1.55 +// 0 - feasible solution found;
    1.56 +// 1 - problem has no feasible solution;
    1.57 +// 2 - iterations limit exceeded;
    1.58 +// 3 - time limit exceeded.
    1.59 +----------------------------------------------------------------------*/
    1.60 +
    1.61 +int ssx_phase_I(SSX *ssx)
    1.62 +{     int m = ssx->m;
    1.63 +      int n = ssx->n;
    1.64 +      int *type = ssx->type;
    1.65 +      mpq_t *lb = ssx->lb;
    1.66 +      mpq_t *ub = ssx->ub;
    1.67 +      mpq_t *coef = ssx->coef;
    1.68 +      int *A_ptr = ssx->A_ptr;
    1.69 +      int *A_ind = ssx->A_ind;
    1.70 +      mpq_t *A_val = ssx->A_val;
    1.71 +      int *Q_col = ssx->Q_col;
    1.72 +      mpq_t *bbar = ssx->bbar;
    1.73 +      mpq_t *pi = ssx->pi;
    1.74 +      mpq_t *cbar = ssx->cbar;
    1.75 +      int *orig_type, orig_dir;
    1.76 +      mpq_t *orig_lb, *orig_ub, *orig_coef;
    1.77 +      int i, k, ret;
    1.78 +      /* save components of the original LP problem, which are changed
    1.79 +         by the routine */
    1.80 +      orig_type = xcalloc(1+m+n, sizeof(int));
    1.81 +      orig_lb = xcalloc(1+m+n, sizeof(mpq_t));
    1.82 +      orig_ub = xcalloc(1+m+n, sizeof(mpq_t));
    1.83 +      orig_coef = xcalloc(1+m+n, sizeof(mpq_t));
    1.84 +      for (k = 1; k <= m+n; k++)
    1.85 +      {  orig_type[k] = type[k];
    1.86 +         mpq_init(orig_lb[k]);
    1.87 +         mpq_set(orig_lb[k], lb[k]);
    1.88 +         mpq_init(orig_ub[k]);
    1.89 +         mpq_set(orig_ub[k], ub[k]);
    1.90 +      }
    1.91 +      orig_dir = ssx->dir;
    1.92 +      for (k = 0; k <= m+n; k++)
    1.93 +      {  mpq_init(orig_coef[k]);
    1.94 +         mpq_set(orig_coef[k], coef[k]);
    1.95 +      }
    1.96 +      /* build an artificial basic solution, which is primal feasible,
    1.97 +         and also build an auxiliary objective function to minimize the
    1.98 +         sum of infeasibilities for the original problem */
    1.99 +      ssx->dir = SSX_MIN;
   1.100 +      for (k = 0; k <= m+n; k++) mpq_set_si(coef[k], 0, 1);
   1.101 +      mpq_set_si(bbar[0], 0, 1);
   1.102 +      for (i = 1; i <= m; i++)
   1.103 +      {  int t;
   1.104 +         k = Q_col[i]; /* x[k] = xB[i] */
   1.105 +         t = type[k];
   1.106 +         if (t == SSX_LO || t == SSX_DB || t == SSX_FX)
   1.107 +         {  /* in the original problem x[k] has lower bound */
   1.108 +            if (mpq_cmp(bbar[i], lb[k]) < 0)
   1.109 +            {  /* which is violated */
   1.110 +               type[k] = SSX_UP;
   1.111 +               mpq_set(ub[k], lb[k]);
   1.112 +               mpq_set_si(lb[k], 0, 1);
   1.113 +               mpq_set_si(coef[k], -1, 1);
   1.114 +               mpq_add(bbar[0], bbar[0], ub[k]);
   1.115 +               mpq_sub(bbar[0], bbar[0], bbar[i]);
   1.116 +            }
   1.117 +         }
   1.118 +         if (t == SSX_UP || t == SSX_DB || t == SSX_FX)
   1.119 +         {  /* in the original problem x[k] has upper bound */
   1.120 +            if (mpq_cmp(bbar[i], ub[k]) > 0)
   1.121 +            {  /* which is violated */
   1.122 +               type[k] = SSX_LO;
   1.123 +               mpq_set(lb[k], ub[k]);
   1.124 +               mpq_set_si(ub[k], 0, 1);
   1.125 +               mpq_set_si(coef[k], +1, 1);
   1.126 +               mpq_add(bbar[0], bbar[0], bbar[i]);
   1.127 +               mpq_sub(bbar[0], bbar[0], lb[k]);
   1.128 +            }
   1.129 +         }
   1.130 +      }
   1.131 +      /* now the initial basic solution should be primal feasible due
   1.132 +         to changes of bounds of some basic variables, which turned to
   1.133 +         implicit artifical variables */
   1.134 +      /* compute simplex multipliers and reduced costs */
   1.135 +      ssx_eval_pi(ssx);
   1.136 +      ssx_eval_cbar(ssx);
   1.137 +      /* display initial progress of the search */
   1.138 +      show_progress(ssx, 1);
   1.139 +      /* main loop starts here */
   1.140 +      for (;;)
   1.141 +      {  /* display current progress of the search */
   1.142 +#if 0
   1.143 +         if (utime() - ssx->tm_lag >= ssx->out_frq - 0.001)
   1.144 +#else
   1.145 +         if (xdifftime(xtime(), ssx->tm_lag) >= ssx->out_frq - 0.001)
   1.146 +#endif
   1.147 +            show_progress(ssx, 1);
   1.148 +         /* we do not need to wait until all artificial variables have
   1.149 +            left the basis */
   1.150 +         if (mpq_sgn(bbar[0]) == 0)
   1.151 +         {  /* the sum of infeasibilities is zero, therefore the current
   1.152 +               solution is primal feasible for the original problem */
   1.153 +            ret = 0;
   1.154 +            break;
   1.155 +         }
   1.156 +         /* check if the iterations limit has been exhausted */
   1.157 +         if (ssx->it_lim == 0)
   1.158 +         {  ret = 2;
   1.159 +            break;
   1.160 +         }
   1.161 +         /* check if the time limit has been exhausted */
   1.162 +#if 0
   1.163 +         if (ssx->tm_lim >= 0.0 && ssx->tm_lim <= utime() - ssx->tm_beg)
   1.164 +#else
   1.165 +         if (ssx->tm_lim >= 0.0 &&
   1.166 +             ssx->tm_lim <= xdifftime(xtime(), ssx->tm_beg))
   1.167 +#endif
   1.168 +         {  ret = 3;
   1.169 +            break;
   1.170 +         }
   1.171 +         /* choose non-basic variable xN[q] */
   1.172 +         ssx_chuzc(ssx);
   1.173 +         /* if xN[q] cannot be chosen, the sum of infeasibilities is
   1.174 +            minimal but non-zero; therefore the original problem has no
   1.175 +            primal feasible solution */
   1.176 +         if (ssx->q == 0)
   1.177 +         {  ret = 1;
   1.178 +            break;
   1.179 +         }
   1.180 +         /* compute q-th column of the simplex table */
   1.181 +         ssx_eval_col(ssx);
   1.182 +         /* choose basic variable xB[p] */
   1.183 +         ssx_chuzr(ssx);
   1.184 +         /* the sum of infeasibilities cannot be negative, therefore
   1.185 +            the auxiliary lp problem cannot have unbounded solution */
   1.186 +         xassert(ssx->p != 0);
   1.187 +         /* update values of basic variables */
   1.188 +         ssx_update_bbar(ssx);
   1.189 +         if (ssx->p > 0)
   1.190 +         {  /* compute p-th row of the inverse inv(B) */
   1.191 +            ssx_eval_rho(ssx);
   1.192 +            /* compute p-th row of the simplex table */
   1.193 +            ssx_eval_row(ssx);
   1.194 +            xassert(mpq_cmp(ssx->aq[ssx->p], ssx->ap[ssx->q]) == 0);
   1.195 +            /* update simplex multipliers */
   1.196 +            ssx_update_pi(ssx);
   1.197 +            /* update reduced costs of non-basic variables */
   1.198 +            ssx_update_cbar(ssx);
   1.199 +         }
   1.200 +         /* xB[p] is leaving the basis; if it is implicit artificial
   1.201 +            variable, the corresponding residual vanishes; therefore
   1.202 +            bounds of this variable should be restored to the original
   1.203 +            values */
   1.204 +         if (ssx->p > 0)
   1.205 +         {  k = Q_col[ssx->p]; /* x[k] = xB[p] */
   1.206 +            if (type[k] != orig_type[k])
   1.207 +            {  /* x[k] is implicit artificial variable */
   1.208 +               type[k] = orig_type[k];
   1.209 +               mpq_set(lb[k], orig_lb[k]);
   1.210 +               mpq_set(ub[k], orig_ub[k]);
   1.211 +               xassert(ssx->p_stat == SSX_NL || ssx->p_stat == SSX_NU);
   1.212 +               ssx->p_stat = (ssx->p_stat == SSX_NL ? SSX_NU : SSX_NL);
   1.213 +               if (type[k] == SSX_FX) ssx->p_stat = SSX_NS;
   1.214 +               /* nullify the objective coefficient at x[k] */
   1.215 +               mpq_set_si(coef[k], 0, 1);
   1.216 +               /* since coef[k] has been changed, we need to compute
   1.217 +                  new reduced cost of x[k], which it will have in the
   1.218 +                  adjacent basis */
   1.219 +               /* the formula d[j] = cN[j] - pi' * N[j] is used (note
   1.220 +                  that the vector pi is not changed, because it depends
   1.221 +                  on objective coefficients at basic variables, but in
   1.222 +                  the adjacent basis, for which the vector pi has been
   1.223 +                  just recomputed, x[k] is non-basic) */
   1.224 +               if (k <= m)
   1.225 +               {  /* x[k] is auxiliary variable */
   1.226 +                  mpq_neg(cbar[ssx->q], pi[k]);
   1.227 +               }
   1.228 +               else
   1.229 +               {  /* x[k] is structural variable */
   1.230 +                  int ptr;
   1.231 +                  mpq_t temp;
   1.232 +                  mpq_init(temp);
   1.233 +                  mpq_set_si(cbar[ssx->q], 0, 1);
   1.234 +                  for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++)
   1.235 +                  {  mpq_mul(temp, pi[A_ind[ptr]], A_val[ptr]);
   1.236 +                     mpq_add(cbar[ssx->q], cbar[ssx->q], temp);
   1.237 +                  }
   1.238 +                  mpq_clear(temp);
   1.239 +               }
   1.240 +            }
   1.241 +         }
   1.242 +         /* jump to the adjacent vertex of the polyhedron */
   1.243 +         ssx_change_basis(ssx);
   1.244 +         /* one simplex iteration has been performed */
   1.245 +         if (ssx->it_lim > 0) ssx->it_lim--;
   1.246 +         ssx->it_cnt++;
   1.247 +      }
   1.248 +      /* display final progress of the search */
   1.249 +      show_progress(ssx, 1);
   1.250 +      /* restore components of the original problem, which were changed
   1.251 +         by the routine */
   1.252 +      for (k = 1; k <= m+n; k++)
   1.253 +      {  type[k] = orig_type[k];
   1.254 +         mpq_set(lb[k], orig_lb[k]);
   1.255 +         mpq_clear(orig_lb[k]);
   1.256 +         mpq_set(ub[k], orig_ub[k]);
   1.257 +         mpq_clear(orig_ub[k]);
   1.258 +      }
   1.259 +      ssx->dir = orig_dir;
   1.260 +      for (k = 0; k <= m+n; k++)
   1.261 +      {  mpq_set(coef[k], orig_coef[k]);
   1.262 +         mpq_clear(orig_coef[k]);
   1.263 +      }
   1.264 +      xfree(orig_type);
   1.265 +      xfree(orig_lb);
   1.266 +      xfree(orig_ub);
   1.267 +      xfree(orig_coef);
   1.268 +      /* return to the calling program */
   1.269 +      return ret;
   1.270 +}
   1.271 +
   1.272 +/*----------------------------------------------------------------------
   1.273 +// ssx_phase_II - find optimal solution.
   1.274 +//
   1.275 +// This routine implements phase II of the primal simplex method.
   1.276 +//
   1.277 +// On exit the routine returns one of the following codes:
   1.278 +//
   1.279 +// 0 - optimal solution found;
   1.280 +// 1 - problem has unbounded solution;
   1.281 +// 2 - iterations limit exceeded;
   1.282 +// 3 - time limit exceeded.
   1.283 +----------------------------------------------------------------------*/
   1.284 +
   1.285 +int ssx_phase_II(SSX *ssx)
   1.286 +{     int ret;
   1.287 +      /* display initial progress of the search */
   1.288 +      show_progress(ssx, 2);
   1.289 +      /* main loop starts here */
   1.290 +      for (;;)
   1.291 +      {  /* display current progress of the search */
   1.292 +#if 0
   1.293 +         if (utime() - ssx->tm_lag >= ssx->out_frq - 0.001)
   1.294 +#else
   1.295 +         if (xdifftime(xtime(), ssx->tm_lag) >= ssx->out_frq - 0.001)
   1.296 +#endif
   1.297 +            show_progress(ssx, 2);
   1.298 +         /* check if the iterations limit has been exhausted */
   1.299 +         if (ssx->it_lim == 0)
   1.300 +         {  ret = 2;
   1.301 +            break;
   1.302 +         }
   1.303 +         /* check if the time limit has been exhausted */
   1.304 +#if 0
   1.305 +         if (ssx->tm_lim >= 0.0 && ssx->tm_lim <= utime() - ssx->tm_beg)
   1.306 +#else
   1.307 +         if (ssx->tm_lim >= 0.0 &&
   1.308 +             ssx->tm_lim <= xdifftime(xtime(), ssx->tm_beg))
   1.309 +#endif
   1.310 +         {  ret = 3;
   1.311 +            break;
   1.312 +         }
   1.313 +         /* choose non-basic variable xN[q] */
   1.314 +         ssx_chuzc(ssx);
   1.315 +         /* if xN[q] cannot be chosen, the current basic solution is
   1.316 +            dual feasible and therefore optimal */
   1.317 +         if (ssx->q == 0)
   1.318 +         {  ret = 0;
   1.319 +            break;
   1.320 +         }
   1.321 +         /* compute q-th column of the simplex table */
   1.322 +         ssx_eval_col(ssx);
   1.323 +         /* choose basic variable xB[p] */
   1.324 +         ssx_chuzr(ssx);
   1.325 +         /* if xB[p] cannot be chosen, the problem has no dual feasible
   1.326 +            solution (i.e. unbounded) */
   1.327 +         if (ssx->p == 0)
   1.328 +         {  ret = 1;
   1.329 +            break;
   1.330 +         }
   1.331 +         /* update values of basic variables */
   1.332 +         ssx_update_bbar(ssx);
   1.333 +         if (ssx->p > 0)
   1.334 +         {  /* compute p-th row of the inverse inv(B) */
   1.335 +            ssx_eval_rho(ssx);
   1.336 +            /* compute p-th row of the simplex table */
   1.337 +            ssx_eval_row(ssx);
   1.338 +            xassert(mpq_cmp(ssx->aq[ssx->p], ssx->ap[ssx->q]) == 0);
   1.339 +#if 0
   1.340 +            /* update simplex multipliers */
   1.341 +            ssx_update_pi(ssx);
   1.342 +#endif
   1.343 +            /* update reduced costs of non-basic variables */
   1.344 +            ssx_update_cbar(ssx);
   1.345 +         }
   1.346 +         /* jump to the adjacent vertex of the polyhedron */
   1.347 +         ssx_change_basis(ssx);
   1.348 +         /* one simplex iteration has been performed */
   1.349 +         if (ssx->it_lim > 0) ssx->it_lim--;
   1.350 +         ssx->it_cnt++;
   1.351 +      }
   1.352 +      /* display final progress of the search */
   1.353 +      show_progress(ssx, 2);
   1.354 +      /* return to the calling program */
   1.355 +      return ret;
   1.356 +}
   1.357 +
   1.358 +/*----------------------------------------------------------------------
   1.359 +// ssx_driver - base driver to exact simplex method.
   1.360 +//
   1.361 +// This routine is a base driver to a version of the primal simplex
   1.362 +// method using exact (bignum) arithmetic.
   1.363 +//
   1.364 +// On exit the routine returns one of the following codes:
   1.365 +//
   1.366 +// 0 - optimal solution found;
   1.367 +// 1 - problem has no feasible solution;
   1.368 +// 2 - problem has unbounded solution;
   1.369 +// 3 - iterations limit exceeded (phase I);
   1.370 +// 4 - iterations limit exceeded (phase II);
   1.371 +// 5 - time limit exceeded (phase I);
   1.372 +// 6 - time limit exceeded (phase II);
   1.373 +// 7 - initial basis matrix is exactly singular.
   1.374 +----------------------------------------------------------------------*/
   1.375 +
   1.376 +int ssx_driver(SSX *ssx)
   1.377 +{     int m = ssx->m;
   1.378 +      int *type = ssx->type;
   1.379 +      mpq_t *lb = ssx->lb;
   1.380 +      mpq_t *ub = ssx->ub;
   1.381 +      int *Q_col = ssx->Q_col;
   1.382 +      mpq_t *bbar = ssx->bbar;
   1.383 +      int i, k, ret;
   1.384 +      ssx->tm_beg = xtime();
   1.385 +      /* factorize the initial basis matrix */
   1.386 +      if (ssx_factorize(ssx))
   1.387 +      {  xprintf("Initial basis matrix is singular\n");
   1.388 +         ret = 7;
   1.389 +         goto done;
   1.390 +      }
   1.391 +      /* compute values of basic variables */
   1.392 +      ssx_eval_bbar(ssx);
   1.393 +      /* check if the initial basic solution is primal feasible */
   1.394 +      for (i = 1; i <= m; i++)
   1.395 +      {  int t;
   1.396 +         k = Q_col[i]; /* x[k] = xB[i] */
   1.397 +         t = type[k];
   1.398 +         if (t == SSX_LO || t == SSX_DB || t == SSX_FX)
   1.399 +         {  /* x[k] has lower bound */
   1.400 +            if (mpq_cmp(bbar[i], lb[k]) < 0)
   1.401 +            {  /* which is violated */
   1.402 +               break;
   1.403 +            }
   1.404 +         }
   1.405 +         if (t == SSX_UP || t == SSX_DB || t == SSX_FX)
   1.406 +         {  /* x[k] has upper bound */
   1.407 +            if (mpq_cmp(bbar[i], ub[k]) > 0)
   1.408 +            {  /* which is violated */
   1.409 +               break;
   1.410 +            }
   1.411 +         }
   1.412 +      }
   1.413 +      if (i > m)
   1.414 +      {  /* no basic variable violates its bounds */
   1.415 +         ret = 0;
   1.416 +         goto skip;
   1.417 +      }
   1.418 +      /* phase I: find primal feasible solution */
   1.419 +      ret = ssx_phase_I(ssx);
   1.420 +      switch (ret)
   1.421 +      {  case 0:
   1.422 +            ret = 0;
   1.423 +            break;
   1.424 +         case 1:
   1.425 +            xprintf("PROBLEM HAS NO FEASIBLE SOLUTION\n");
   1.426 +            ret = 1;
   1.427 +            break;
   1.428 +         case 2:
   1.429 +            xprintf("ITERATIONS LIMIT EXCEEDED; SEARCH TERMINATED\n");
   1.430 +            ret = 3;
   1.431 +            break;
   1.432 +         case 3:
   1.433 +            xprintf("TIME LIMIT EXCEEDED; SEARCH TERMINATED\n");
   1.434 +            ret = 5;
   1.435 +            break;
   1.436 +         default:
   1.437 +            xassert(ret != ret);
   1.438 +      }
   1.439 +      /* compute values of basic variables (actually only the objective
   1.440 +         value needs to be computed) */
   1.441 +      ssx_eval_bbar(ssx);
   1.442 +skip: /* compute simplex multipliers */
   1.443 +      ssx_eval_pi(ssx);
   1.444 +      /* compute reduced costs of non-basic variables */
   1.445 +      ssx_eval_cbar(ssx);
   1.446 +      /* if phase I failed, do not start phase II */
   1.447 +      if (ret != 0) goto done;
   1.448 +      /* phase II: find optimal solution */
   1.449 +      ret = ssx_phase_II(ssx);
   1.450 +      switch (ret)
   1.451 +      {  case 0:
   1.452 +            xprintf("OPTIMAL SOLUTION FOUND\n");
   1.453 +            ret = 0;
   1.454 +            break;
   1.455 +         case 1:
   1.456 +            xprintf("PROBLEM HAS UNBOUNDED SOLUTION\n");
   1.457 +            ret = 2;
   1.458 +            break;
   1.459 +         case 2:
   1.460 +            xprintf("ITERATIONS LIMIT EXCEEDED; SEARCH TERMINATED\n");
   1.461 +            ret = 4;
   1.462 +            break;
   1.463 +         case 3:
   1.464 +            xprintf("TIME LIMIT EXCEEDED; SEARCH TERMINATED\n");
   1.465 +            ret = 6;
   1.466 +            break;
   1.467 +         default:
   1.468 +            xassert(ret != ret);
   1.469 +      }
   1.470 +done: /* decrease the time limit by the spent amount of time */
   1.471 +      if (ssx->tm_lim >= 0.0)
   1.472 +#if 0
   1.473 +      {  ssx->tm_lim -= utime() - ssx->tm_beg;
   1.474 +#else
   1.475 +      {  ssx->tm_lim -= xdifftime(xtime(), ssx->tm_beg);
   1.476 +#endif
   1.477 +         if (ssx->tm_lim < 0.0) ssx->tm_lim = 0.0;
   1.478 +      }
   1.479 +      return ret;
   1.480 +}
   1.481 +
   1.482 +/* eof */