1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpapi07.c	Sun Nov 06 20:59:10 2011 +0100
     1.3 @@ -0,0 +1,451 @@
     1.4 +/* glpapi07.c (exact simplex solver) */
     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 "glpapi.h"
    1.29 +#include "glpssx.h"
    1.30 +
    1.31 +/***********************************************************************
    1.32 +*  NAME
    1.33 +*
    1.34 +*  glp_exact - solve LP problem in exact arithmetic
    1.35 +*
    1.36 +*  SYNOPSIS
    1.37 +*
    1.38 +*  int glp_exact(glp_prob *lp, const glp_smcp *parm);
    1.39 +*
    1.40 +*  DESCRIPTION
    1.41 +*
    1.42 +*  The routine glp_exact is a tentative implementation of the primal
    1.43 +*  two-phase simplex method based on exact (rational) arithmetic. It is
    1.44 +*  similar to the routine glp_simplex, however, for all internal
    1.45 +*  computations it uses arithmetic of rational numbers, which is exact
    1.46 +*  in mathematical sense, i.e. free of round-off errors unlike floating
    1.47 +*  point arithmetic.
    1.48 +*
    1.49 +*  Note that the routine glp_exact uses inly two control parameters
    1.50 +*  passed in the structure glp_smcp, namely, it_lim and tm_lim.
    1.51 +*
    1.52 +*  RETURNS
    1.53 +*
    1.54 +*  0  The LP problem instance has been successfully solved. This code
    1.55 +*     does not necessarily mean that the solver has found optimal
    1.56 +*     solution. It only means that the solution process was successful.
    1.57 +*
    1.58 +*  GLP_EBADB
    1.59 +*     Unable to start the search, because the initial basis specified
    1.60 +*     in the problem object is invalid--the number of basic (auxiliary
    1.61 +*     and structural) variables is not the same as the number of rows in
    1.62 +*     the problem object.
    1.63 +*
    1.64 +*  GLP_ESING
    1.65 +*     Unable to start the search, because the basis matrix correspodning
    1.66 +*     to the initial basis is exactly singular.
    1.67 +*
    1.68 +*  GLP_EBOUND
    1.69 +*     Unable to start the search, because some double-bounded variables
    1.70 +*     have incorrect bounds.
    1.71 +*
    1.72 +*  GLP_EFAIL
    1.73 +*     The problem has no rows/columns.
    1.74 +*
    1.75 +*  GLP_EITLIM
    1.76 +*     The search was prematurely terminated, because the simplex
    1.77 +*     iteration limit has been exceeded.
    1.78 +*
    1.79 +*  GLP_ETMLIM
    1.80 +*     The search was prematurely terminated, because the time limit has
    1.81 +*     been exceeded. */
    1.82 +
    1.83 +static void set_d_eps(mpq_t x, double val)
    1.84 +{     /* convert double val to rational x obtaining a more adequate
    1.85 +         fraction than provided by mpq_set_d due to allowing a small
    1.86 +         approximation error specified by a given relative tolerance;
    1.87 +         for example, mpq_set_d would give the following
    1.88 +         1/3 ~= 0.333333333333333314829616256247391... ->
    1.89 +             -> 6004799503160661/18014398509481984
    1.90 +         while this routine gives exactly 1/3 */
    1.91 +      int s, n, j;
    1.92 +      double f, p, q, eps = 1e-9;
    1.93 +      mpq_t temp;
    1.94 +      xassert(-DBL_MAX <= val && val <= +DBL_MAX);
    1.95 +#if 1 /* 30/VII-2008 */
    1.96 +      if (val == floor(val))
    1.97 +      {  /* if val is integral, do not approximate */
    1.98 +         mpq_set_d(x, val);
    1.99 +         goto done;
   1.100 +      }
   1.101 +#endif
   1.102 +      if (val > 0.0)
   1.103 +         s = +1;
   1.104 +      else if (val < 0.0)
   1.105 +         s = -1;
   1.106 +      else
   1.107 +      {  mpq_set_si(x, 0, 1);
   1.108 +         goto done;
   1.109 +      }
   1.110 +      f = frexp(fabs(val), &n);
   1.111 +      /* |val| = f * 2^n, where 0.5 <= f < 1.0 */
   1.112 +      fp2rat(f, 0.1 * eps, &p, &q);
   1.113 +      /* f ~= p / q, where p and q are integers */
   1.114 +      mpq_init(temp);
   1.115 +      mpq_set_d(x, p);
   1.116 +      mpq_set_d(temp, q);
   1.117 +      mpq_div(x, x, temp);
   1.118 +      mpq_set_si(temp, 1, 1);
   1.119 +      for (j = 1; j <= abs(n); j++)
   1.120 +         mpq_add(temp, temp, temp);
   1.121 +      if (n > 0)
   1.122 +         mpq_mul(x, x, temp);
   1.123 +      else if (n < 0)
   1.124 +         mpq_div(x, x, temp);
   1.125 +      mpq_clear(temp);
   1.126 +      if (s < 0) mpq_neg(x, x);
   1.127 +      /* check that the desired tolerance has been attained */
   1.128 +      xassert(fabs(val - mpq_get_d(x)) <= eps * (1.0 + fabs(val)));
   1.129 +done: return;
   1.130 +}
   1.131 +
   1.132 +static void load_data(SSX *ssx, LPX *lp)
   1.133 +{     /* load LP problem data into simplex solver workspace */
   1.134 +      int m = ssx->m;
   1.135 +      int n = ssx->n;
   1.136 +      int nnz = ssx->A_ptr[n+1]-1;
   1.137 +      int j, k, type, loc, len, *ind;
   1.138 +      double lb, ub, coef, *val;
   1.139 +      xassert(lpx_get_num_rows(lp) == m);
   1.140 +      xassert(lpx_get_num_cols(lp) == n);
   1.141 +      xassert(lpx_get_num_nz(lp) == nnz);
   1.142 +      /* types and bounds of rows and columns */
   1.143 +      for (k = 1; k <= m+n; k++)
   1.144 +      {  if (k <= m)
   1.145 +         {  type = lpx_get_row_type(lp, k);
   1.146 +            lb = lpx_get_row_lb(lp, k);
   1.147 +            ub = lpx_get_row_ub(lp, k);
   1.148 +         }
   1.149 +         else
   1.150 +         {  type = lpx_get_col_type(lp, k-m);
   1.151 +            lb = lpx_get_col_lb(lp, k-m);
   1.152 +            ub = lpx_get_col_ub(lp, k-m);
   1.153 +         }
   1.154 +         switch (type)
   1.155 +         {  case LPX_FR: type = SSX_FR; break;
   1.156 +            case LPX_LO: type = SSX_LO; break;
   1.157 +            case LPX_UP: type = SSX_UP; break;
   1.158 +            case LPX_DB: type = SSX_DB; break;
   1.159 +            case LPX_FX: type = SSX_FX; break;
   1.160 +            default: xassert(type != type);
   1.161 +         }
   1.162 +         ssx->type[k] = type;
   1.163 +         set_d_eps(ssx->lb[k], lb);
   1.164 +         set_d_eps(ssx->ub[k], ub);
   1.165 +      }
   1.166 +      /* optimization direction */
   1.167 +      switch (lpx_get_obj_dir(lp))
   1.168 +      {  case LPX_MIN: ssx->dir = SSX_MIN; break;
   1.169 +         case LPX_MAX: ssx->dir = SSX_MAX; break;
   1.170 +         default: xassert(lp != lp);
   1.171 +      }
   1.172 +      /* objective coefficients */
   1.173 +      for (k = 0; k <= m+n; k++)
   1.174 +      {  if (k == 0)
   1.175 +            coef = lpx_get_obj_coef(lp, 0);
   1.176 +         else if (k <= m)
   1.177 +            coef = 0.0;
   1.178 +         else
   1.179 +            coef = lpx_get_obj_coef(lp, k-m);
   1.180 +         set_d_eps(ssx->coef[k], coef);
   1.181 +      }
   1.182 +      /* constraint coefficients */
   1.183 +      ind = xcalloc(1+m, sizeof(int));
   1.184 +      val = xcalloc(1+m, sizeof(double));
   1.185 +      loc = 0;
   1.186 +      for (j = 1; j <= n; j++)
   1.187 +      {  ssx->A_ptr[j] = loc+1;
   1.188 +         len = lpx_get_mat_col(lp, j, ind, val);
   1.189 +         for (k = 1; k <= len; k++)
   1.190 +         {  loc++;
   1.191 +            ssx->A_ind[loc] = ind[k];
   1.192 +            set_d_eps(ssx->A_val[loc], val[k]);
   1.193 +         }
   1.194 +      }
   1.195 +      xassert(loc == nnz);
   1.196 +      xfree(ind);
   1.197 +      xfree(val);
   1.198 +      return;
   1.199 +}
   1.200 +
   1.201 +static int load_basis(SSX *ssx, LPX *lp)
   1.202 +{     /* load current LP basis into simplex solver workspace */
   1.203 +      int m = ssx->m;
   1.204 +      int n = ssx->n;
   1.205 +      int *type = ssx->type;
   1.206 +      int *stat = ssx->stat;
   1.207 +      int *Q_row = ssx->Q_row;
   1.208 +      int *Q_col = ssx->Q_col;
   1.209 +      int i, j, k;
   1.210 +      xassert(lpx_get_num_rows(lp) == m);
   1.211 +      xassert(lpx_get_num_cols(lp) == n);
   1.212 +      /* statuses of rows and columns */
   1.213 +      for (k = 1; k <= m+n; k++)
   1.214 +      {  if (k <= m)
   1.215 +            stat[k] = lpx_get_row_stat(lp, k);
   1.216 +         else
   1.217 +            stat[k] = lpx_get_col_stat(lp, k-m);
   1.218 +         switch (stat[k])
   1.219 +         {  case LPX_BS:
   1.220 +               stat[k] = SSX_BS;
   1.221 +               break;
   1.222 +            case LPX_NL:
   1.223 +               stat[k] = SSX_NL;
   1.224 +               xassert(type[k] == SSX_LO || type[k] == SSX_DB);
   1.225 +               break;
   1.226 +            case LPX_NU:
   1.227 +               stat[k] = SSX_NU;
   1.228 +               xassert(type[k] == SSX_UP || type[k] == SSX_DB);
   1.229 +               break;
   1.230 +            case LPX_NF:
   1.231 +               stat[k] = SSX_NF;
   1.232 +               xassert(type[k] == SSX_FR);
   1.233 +               break;
   1.234 +            case LPX_NS:
   1.235 +               stat[k] = SSX_NS;
   1.236 +               xassert(type[k] == SSX_FX);
   1.237 +               break;
   1.238 +            default:
   1.239 +               xassert(stat != stat);
   1.240 +         }
   1.241 +      }
   1.242 +      /* build permutation matix Q */
   1.243 +      i = j = 0;
   1.244 +      for (k = 1; k <= m+n; k++)
   1.245 +      {  if (stat[k] == SSX_BS)
   1.246 +         {  i++;
   1.247 +            if (i > m) return 1;
   1.248 +            Q_row[k] = i, Q_col[i] = k;
   1.249 +         }
   1.250 +         else
   1.251 +         {  j++;
   1.252 +            if (j > n) return 1;
   1.253 +            Q_row[k] = m+j, Q_col[m+j] = k;
   1.254 +         }
   1.255 +      }
   1.256 +      xassert(i == m && j == n);
   1.257 +      return 0;
   1.258 +}
   1.259 +
   1.260 +int glp_exact(glp_prob *lp, const glp_smcp *parm)
   1.261 +{     glp_smcp _parm;
   1.262 +      SSX *ssx;
   1.263 +      int m = lpx_get_num_rows(lp);
   1.264 +      int n = lpx_get_num_cols(lp);
   1.265 +      int nnz = lpx_get_num_nz(lp);
   1.266 +      int i, j, k, type, pst, dst, ret, *stat;
   1.267 +      double lb, ub, *prim, *dual, sum;
   1.268 +      if (parm == NULL)
   1.269 +         parm = &_parm, glp_init_smcp((glp_smcp *)parm);
   1.270 +      /* check control parameters */
   1.271 +      if (parm->it_lim < 0)
   1.272 +         xerror("glp_exact: it_lim = %d; invalid parameter\n",
   1.273 +            parm->it_lim);
   1.274 +      if (parm->tm_lim < 0)
   1.275 +         xerror("glp_exact: tm_lim = %d; invalid parameter\n",
   1.276 +            parm->tm_lim);
   1.277 +      /* the problem must have at least one row and one column */
   1.278 +      if (!(m > 0 && n > 0))
   1.279 +      {  xprintf("glp_exact: problem has no rows/columns\n");
   1.280 +         return GLP_EFAIL;
   1.281 +      }
   1.282 +#if 1
   1.283 +      /* basic solution is currently undefined */
   1.284 +      lp->pbs_stat = lp->dbs_stat = GLP_UNDEF;
   1.285 +      lp->obj_val = 0.0;
   1.286 +      lp->some = 0;
   1.287 +#endif
   1.288 +      /* check that all double-bounded variables have correct bounds */
   1.289 +      for (k = 1; k <= m+n; k++)
   1.290 +      {  if (k <= m)
   1.291 +         {  type = lpx_get_row_type(lp, k);
   1.292 +            lb = lpx_get_row_lb(lp, k);
   1.293 +            ub = lpx_get_row_ub(lp, k);
   1.294 +         }
   1.295 +         else
   1.296 +         {  type = lpx_get_col_type(lp, k-m);
   1.297 +            lb = lpx_get_col_lb(lp, k-m);
   1.298 +            ub = lpx_get_col_ub(lp, k-m);
   1.299 +         }
   1.300 +         if (type == LPX_DB && lb >= ub)
   1.301 +         {  xprintf("glp_exact: %s %d has invalid bounds\n",
   1.302 +               k <= m ? "row" : "column", k <= m ? k : k-m);
   1.303 +            return GLP_EBOUND;
   1.304 +         }
   1.305 +      }
   1.306 +      /* create the simplex solver workspace */
   1.307 +      xprintf("glp_exact: %d rows, %d columns, %d non-zeros\n",
   1.308 +         m, n, nnz);
   1.309 +#ifdef HAVE_GMP
   1.310 +      xprintf("GNU MP bignum library is being used\n");
   1.311 +#else
   1.312 +      xprintf("GLPK bignum module is being used\n");
   1.313 +      xprintf("(Consider installing GNU MP to attain a much better perf"
   1.314 +         "ormance.)\n");
   1.315 +#endif
   1.316 +      ssx = ssx_create(m, n, nnz);
   1.317 +      /* load LP problem data into the workspace */
   1.318 +      load_data(ssx, lp);
   1.319 +      /* load current LP basis into the workspace */
   1.320 +      if (load_basis(ssx, lp))
   1.321 +      {  xprintf("glp_exact: initial LP basis is invalid\n");
   1.322 +         ret = GLP_EBADB;
   1.323 +         goto done;
   1.324 +      }
   1.325 +      /* inherit some control parameters from the LP object */
   1.326 +#if 0
   1.327 +      ssx->it_lim = lpx_get_int_parm(lp, LPX_K_ITLIM);
   1.328 +      ssx->it_cnt = lpx_get_int_parm(lp, LPX_K_ITCNT);
   1.329 +      ssx->tm_lim = lpx_get_real_parm(lp, LPX_K_TMLIM);
   1.330 +#else
   1.331 +      ssx->it_lim = parm->it_lim;
   1.332 +      ssx->it_cnt = lp->it_cnt;
   1.333 +      ssx->tm_lim = (double)parm->tm_lim / 1000.0;
   1.334 +#endif
   1.335 +      ssx->out_frq = 5.0;
   1.336 +      ssx->tm_beg = xtime();
   1.337 +      ssx->tm_lag = xlset(0);
   1.338 +      /* solve LP */
   1.339 +      ret = ssx_driver(ssx);
   1.340 +      /* copy back some statistics to the LP object */
   1.341 +#if 0
   1.342 +      lpx_set_int_parm(lp, LPX_K_ITLIM, ssx->it_lim);
   1.343 +      lpx_set_int_parm(lp, LPX_K_ITCNT, ssx->it_cnt);
   1.344 +      lpx_set_real_parm(lp, LPX_K_TMLIM, ssx->tm_lim);
   1.345 +#else
   1.346 +      lp->it_cnt = ssx->it_cnt;
   1.347 +#endif
   1.348 +      /* analyze the return code */
   1.349 +      switch (ret)
   1.350 +      {  case 0:
   1.351 +            /* optimal solution found */
   1.352 +            ret = 0;
   1.353 +            pst = LPX_P_FEAS, dst = LPX_D_FEAS;
   1.354 +            break;
   1.355 +         case 1:
   1.356 +            /* problem has no feasible solution */
   1.357 +            ret = 0;
   1.358 +            pst = LPX_P_NOFEAS, dst = LPX_D_INFEAS;
   1.359 +            break;
   1.360 +         case 2:
   1.361 +            /* problem has unbounded solution */
   1.362 +            ret = 0;
   1.363 +            pst = LPX_P_FEAS, dst = LPX_D_NOFEAS;
   1.364 +#if 1
   1.365 +            xassert(1 <= ssx->q && ssx->q <= n);
   1.366 +            lp->some = ssx->Q_col[m + ssx->q];
   1.367 +            xassert(1 <= lp->some && lp->some <= m+n);
   1.368 +#endif
   1.369 +            break;
   1.370 +         case 3:
   1.371 +            /* iteration limit exceeded (phase I) */
   1.372 +            ret = GLP_EITLIM;
   1.373 +            pst = LPX_P_INFEAS, dst = LPX_D_INFEAS;
   1.374 +            break;
   1.375 +         case 4:
   1.376 +            /* iteration limit exceeded (phase II) */
   1.377 +            ret = GLP_EITLIM;
   1.378 +            pst = LPX_P_FEAS, dst = LPX_D_INFEAS;
   1.379 +            break;
   1.380 +         case 5:
   1.381 +            /* time limit exceeded (phase I) */
   1.382 +            ret = GLP_ETMLIM;
   1.383 +            pst = LPX_P_INFEAS, dst = LPX_D_INFEAS;
   1.384 +            break;
   1.385 +         case 6:
   1.386 +            /* time limit exceeded (phase II) */
   1.387 +            ret = GLP_ETMLIM;
   1.388 +            pst = LPX_P_FEAS, dst = LPX_D_INFEAS;
   1.389 +            break;
   1.390 +         case 7:
   1.391 +            /* initial basis matrix is singular */
   1.392 +            ret = GLP_ESING;
   1.393 +            goto done;
   1.394 +         default:
   1.395 +            xassert(ret != ret);
   1.396 +      }
   1.397 +      /* obtain final basic solution components */
   1.398 +      stat = xcalloc(1+m+n, sizeof(int));
   1.399 +      prim = xcalloc(1+m+n, sizeof(double));
   1.400 +      dual = xcalloc(1+m+n, sizeof(double));
   1.401 +      for (k = 1; k <= m+n; k++)
   1.402 +      {  if (ssx->stat[k] == SSX_BS)
   1.403 +         {  i = ssx->Q_row[k]; /* x[k] = xB[i] */
   1.404 +            xassert(1 <= i && i <= m);
   1.405 +            stat[k] = LPX_BS;
   1.406 +            prim[k] = mpq_get_d(ssx->bbar[i]);
   1.407 +            dual[k] = 0.0;
   1.408 +         }
   1.409 +         else
   1.410 +         {  j = ssx->Q_row[k] - m; /* x[k] = xN[j] */
   1.411 +            xassert(1 <= j && j <= n);
   1.412 +            switch (ssx->stat[k])
   1.413 +            {  case SSX_NF:
   1.414 +                  stat[k] = LPX_NF;
   1.415 +                  prim[k] = 0.0;
   1.416 +                  break;
   1.417 +               case SSX_NL:
   1.418 +                  stat[k] = LPX_NL;
   1.419 +                  prim[k] = mpq_get_d(ssx->lb[k]);
   1.420 +                  break;
   1.421 +               case SSX_NU:
   1.422 +                  stat[k] = LPX_NU;
   1.423 +                  prim[k] = mpq_get_d(ssx->ub[k]);
   1.424 +                  break;
   1.425 +               case SSX_NS:
   1.426 +                  stat[k] = LPX_NS;
   1.427 +                  prim[k] = mpq_get_d(ssx->lb[k]);
   1.428 +                  break;
   1.429 +               default:
   1.430 +                  xassert(ssx != ssx);
   1.431 +            }
   1.432 +            dual[k] = mpq_get_d(ssx->cbar[j]);
   1.433 +         }
   1.434 +      }
   1.435 +      /* and store them into the LP object */
   1.436 +      pst = pst - LPX_P_UNDEF + GLP_UNDEF;
   1.437 +      dst = dst - LPX_D_UNDEF + GLP_UNDEF;
   1.438 +      for (k = 1; k <= m+n; k++)
   1.439 +         stat[k] = stat[k] - LPX_BS + GLP_BS;
   1.440 +      sum = lpx_get_obj_coef(lp, 0);
   1.441 +      for (j = 1; j <= n; j++)
   1.442 +         sum += lpx_get_obj_coef(lp, j) * prim[m+j];
   1.443 +      lpx_put_solution(lp, 1, &pst, &dst, &sum,
   1.444 +         &stat[0], &prim[0], &dual[0], &stat[m], &prim[m], &dual[m]);
   1.445 +      xfree(stat);
   1.446 +      xfree(prim);
   1.447 +      xfree(dual);
   1.448 +done: /* delete the simplex solver workspace */
   1.449 +      ssx_delete(ssx);
   1.450 +      /* return to the application program */
   1.451 +      return ret;
   1.452 +}
   1.453 +
   1.454 +/* eof */