src/glpios07.c
changeset 1 c445c931472f
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/src/glpios07.c	Mon Dec 06 13:09:21 2010 +0100
     1.3 @@ -0,0 +1,554 @@
     1.4 +/* glpios07.c (mixed cover cut generator) */
     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 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 "glpios.h"
    1.29 +
    1.30 +/*----------------------------------------------------------------------
    1.31 +-- COVER INEQUALITIES
    1.32 +--
    1.33 +-- Consider the set of feasible solutions to 0-1 knapsack problem:
    1.34 +--
    1.35 +--    sum a[j]*x[j] <= b,                                            (1)
    1.36 +--  j in J
    1.37 +--
    1.38 +--    x[j] is binary,                                                (2)
    1.39 +--
    1.40 +-- where, wlog, we assume that a[j] > 0 (since 0-1 variables can be
    1.41 +-- complemented) and a[j] <= b (since a[j] > b implies x[j] = 0).
    1.42 +--
    1.43 +-- A set C within J is called a cover if
    1.44 +--
    1.45 +--    sum a[j] > b.                                                  (3)
    1.46 +--  j in C
    1.47 +--
    1.48 +-- For any cover C the inequality
    1.49 +--
    1.50 +--    sum x[j] <= |C| - 1                                            (4)
    1.51 +--  j in C
    1.52 +--
    1.53 +-- is called a cover inequality and is valid for (1)-(2).
    1.54 +--
    1.55 +-- MIXED COVER INEQUALITIES
    1.56 +--
    1.57 +-- Consider the set of feasible solutions to mixed knapsack problem:
    1.58 +--
    1.59 +--    sum a[j]*x[j] + y <= b,                                        (5)
    1.60 +--  j in J
    1.61 +--
    1.62 +--    x[j] is binary,                                                (6)
    1.63 +--
    1.64 +--    0 <= y <= u is continuous,                                     (7)
    1.65 +--
    1.66 +-- where again we assume that a[j] > 0.
    1.67 +--
    1.68 +-- Let C within J be some set. From (1)-(4) it follows that
    1.69 +--
    1.70 +--    sum a[j] > b - y                                               (8)
    1.71 +--  j in C
    1.72 +--
    1.73 +-- implies
    1.74 +--
    1.75 +--    sum x[j] <= |C| - 1.                                           (9)
    1.76 +--  j in C
    1.77 +--
    1.78 +-- Thus, we need to modify the inequality (9) in such a way that it be
    1.79 +-- a constraint only if the condition (8) is satisfied.
    1.80 +--
    1.81 +-- Consider the following inequality:
    1.82 +--
    1.83 +--    sum x[j] <= |C| - t.                                          (10)
    1.84 +--  j in C
    1.85 +--
    1.86 +-- If 0 < t <= 1, then (10) is equivalent to (9), because all x[j] are
    1.87 +-- binary variables. On the other hand, if t <= 0, (10) being satisfied
    1.88 +-- for any values of x[j] is not a constraint.
    1.89 +--
    1.90 +-- Let
    1.91 +--
    1.92 +--    t' = sum a[j] + y - b.                                        (11)
    1.93 +--       j in C
    1.94 +--
    1.95 +-- It is understood that the condition t' > 0 is equivalent to (8).
    1.96 +-- Besides, from (6)-(7) it follows that t' has an implied upper bound:
    1.97 +--
    1.98 +--    t'max = sum a[j] + u - b.                                     (12)
    1.99 +--          j in C
   1.100 +--
   1.101 +-- This allows to express the parameter t having desired properties:
   1.102 +--
   1.103 +--    t = t' / t'max.                                               (13)
   1.104 +--
   1.105 +-- In fact, t <= 1 by definition, and t > 0 being equivalent to t' > 0
   1.106 +-- is equivalent to (8).
   1.107 +--
   1.108 +-- Thus, the inequality (10), where t is given by formula (13) is valid
   1.109 +-- for (5)-(7).
   1.110 +--
   1.111 +-- Note that if u = 0, then y = 0, so t = 1, and the conditions (8) and
   1.112 +-- (10) is transformed to the conditions (3) and (4).
   1.113 +--
   1.114 +-- GENERATING MIXED COVER CUTS
   1.115 +--
   1.116 +-- To generate a mixed cover cut in the form (10) we need to find such
   1.117 +-- set C which satisfies to the inequality (8) and for which, in turn,
   1.118 +-- the inequality (10) is violated in the current point.
   1.119 +--
   1.120 +-- Substituting t from (13) to (10) gives:
   1.121 +--
   1.122 +--                        1
   1.123 +--    sum x[j] <= |C| - -----  (sum a[j] + y - b),                  (14)
   1.124 +--  j in C              t'max j in C
   1.125 +--
   1.126 +-- and finally we have the cut inequality in the standard form:
   1.127 +--
   1.128 +--    sum x[j] + alfa * y <= beta,                                  (15)
   1.129 +--  j in C
   1.130 +--
   1.131 +-- where:
   1.132 +--
   1.133 +--    alfa = 1 / t'max,                                             (16)
   1.134 +--
   1.135 +--    beta = |C| - alfa *  (sum a[j] - b).                          (17)
   1.136 +--                        j in C                                      */
   1.137 +
   1.138 +#if 1
   1.139 +#define MAXTRY 1000
   1.140 +#else
   1.141 +#define MAXTRY 10000
   1.142 +#endif
   1.143 +
   1.144 +static int cover2(int n, double a[], double b, double u, double x[],
   1.145 +      double y, int cov[], double *_alfa, double *_beta)
   1.146 +{     /* try to generate mixed cover cut using two-element cover */
   1.147 +      int i, j, try = 0, ret = 0;
   1.148 +      double eps, alfa, beta, temp, rmax = 0.001;
   1.149 +      eps = 0.001 * (1.0 + fabs(b));
   1.150 +      for (i = 0+1; i <= n; i++)
   1.151 +      for (j = i+1; j <= n; j++)
   1.152 +      {  /* C = {i, j} */
   1.153 +         try++;
   1.154 +         if (try > MAXTRY) goto done;
   1.155 +         /* check if condition (8) is satisfied */
   1.156 +         if (a[i] + a[j] + y > b + eps)
   1.157 +         {  /* compute parameters for inequality (15) */
   1.158 +            temp = a[i] + a[j] - b;
   1.159 +            alfa = 1.0 / (temp + u);
   1.160 +            beta = 2.0 - alfa * temp;
   1.161 +            /* compute violation of inequality (15) */
   1.162 +            temp = x[i] + x[j] + alfa * y - beta;
   1.163 +            /* choose C providing maximum violation */
   1.164 +            if (rmax < temp)
   1.165 +            {  rmax = temp;
   1.166 +               cov[1] = i;
   1.167 +               cov[2] = j;
   1.168 +               *_alfa = alfa;
   1.169 +               *_beta = beta;
   1.170 +               ret = 1;
   1.171 +            }
   1.172 +         }
   1.173 +      }
   1.174 +done: return ret;
   1.175 +}
   1.176 +
   1.177 +static int cover3(int n, double a[], double b, double u, double x[],
   1.178 +      double y, int cov[], double *_alfa, double *_beta)
   1.179 +{     /* try to generate mixed cover cut using three-element cover */
   1.180 +      int i, j, k, try = 0, ret = 0;
   1.181 +      double eps, alfa, beta, temp, rmax = 0.001;
   1.182 +      eps = 0.001 * (1.0 + fabs(b));
   1.183 +      for (i = 0+1; i <= n; i++)
   1.184 +      for (j = i+1; j <= n; j++)
   1.185 +      for (k = j+1; k <= n; k++)
   1.186 +      {  /* C = {i, j, k} */
   1.187 +         try++;
   1.188 +         if (try > MAXTRY) goto done;
   1.189 +         /* check if condition (8) is satisfied */
   1.190 +         if (a[i] + a[j] + a[k] + y > b + eps)
   1.191 +         {  /* compute parameters for inequality (15) */
   1.192 +            temp = a[i] + a[j] + a[k] - b;
   1.193 +            alfa = 1.0 / (temp + u);
   1.194 +            beta = 3.0 - alfa * temp;
   1.195 +            /* compute violation of inequality (15) */
   1.196 +            temp = x[i] + x[j] + x[k] + alfa * y - beta;
   1.197 +            /* choose C providing maximum violation */
   1.198 +            if (rmax < temp)
   1.199 +            {  rmax = temp;
   1.200 +               cov[1] = i;
   1.201 +               cov[2] = j;
   1.202 +               cov[3] = k;
   1.203 +               *_alfa = alfa;
   1.204 +               *_beta = beta;
   1.205 +               ret = 1;
   1.206 +            }
   1.207 +         }
   1.208 +      }
   1.209 +done: return ret;
   1.210 +}
   1.211 +
   1.212 +static int cover4(int n, double a[], double b, double u, double x[],
   1.213 +      double y, int cov[], double *_alfa, double *_beta)
   1.214 +{     /* try to generate mixed cover cut using four-element cover */
   1.215 +      int i, j, k, l, try = 0, ret = 0;
   1.216 +      double eps, alfa, beta, temp, rmax = 0.001;
   1.217 +      eps = 0.001 * (1.0 + fabs(b));
   1.218 +      for (i = 0+1; i <= n; i++)
   1.219 +      for (j = i+1; j <= n; j++)
   1.220 +      for (k = j+1; k <= n; k++)
   1.221 +      for (l = k+1; l <= n; l++)
   1.222 +      {  /* C = {i, j, k, l} */
   1.223 +         try++;
   1.224 +         if (try > MAXTRY) goto done;
   1.225 +         /* check if condition (8) is satisfied */
   1.226 +         if (a[i] + a[j] + a[k] + a[l] + y > b + eps)
   1.227 +         {  /* compute parameters for inequality (15) */
   1.228 +            temp = a[i] + a[j] + a[k] + a[l] - b;
   1.229 +            alfa = 1.0 / (temp + u);
   1.230 +            beta = 4.0 - alfa * temp;
   1.231 +            /* compute violation of inequality (15) */
   1.232 +            temp = x[i] + x[j] + x[k] + x[l] + alfa * y - beta;
   1.233 +            /* choose C providing maximum violation */
   1.234 +            if (rmax < temp)
   1.235 +            {  rmax = temp;
   1.236 +               cov[1] = i;
   1.237 +               cov[2] = j;
   1.238 +               cov[3] = k;
   1.239 +               cov[4] = l;
   1.240 +               *_alfa = alfa;
   1.241 +               *_beta = beta;
   1.242 +               ret = 1;
   1.243 +            }
   1.244 +         }
   1.245 +      }
   1.246 +done: return ret;
   1.247 +}
   1.248 +
   1.249 +static int cover(int n, double a[], double b, double u, double x[],
   1.250 +      double y, int cov[], double *alfa, double *beta)
   1.251 +{     /* try to generate mixed cover cut;
   1.252 +         input (see (5)):
   1.253 +         n        is the number of binary variables;
   1.254 +         a[1:n]   are coefficients at binary variables;
   1.255 +         b        is the right-hand side;
   1.256 +         u        is upper bound of continuous variable;
   1.257 +         x[1:n]   are values of binary variables at current point;
   1.258 +         y        is value of continuous variable at current point;
   1.259 +         output (see (15), (16), (17)):
   1.260 +         cov[1:r] are indices of binary variables included in cover C,
   1.261 +                  where r is the set cardinality returned on exit;
   1.262 +         alfa     coefficient at continuous variable;
   1.263 +         beta     is the right-hand side; */
   1.264 +      int j;
   1.265 +      /* perform some sanity checks */
   1.266 +      xassert(n >= 2);
   1.267 +      for (j = 1; j <= n; j++) xassert(a[j] > 0.0);
   1.268 +#if 1 /* ??? */
   1.269 +      xassert(b > -1e-5);
   1.270 +#else
   1.271 +      xassert(b > 0.0);
   1.272 +#endif
   1.273 +      xassert(u >= 0.0);
   1.274 +      for (j = 1; j <= n; j++) xassert(0.0 <= x[j] && x[j] <= 1.0);
   1.275 +      xassert(0.0 <= y && y <= u);
   1.276 +      /* try to generate mixed cover cut */
   1.277 +      if (cover2(n, a, b, u, x, y, cov, alfa, beta)) return 2;
   1.278 +      if (cover3(n, a, b, u, x, y, cov, alfa, beta)) return 3;
   1.279 +      if (cover4(n, a, b, u, x, y, cov, alfa, beta)) return 4;
   1.280 +      return 0;
   1.281 +}
   1.282 +
   1.283 +/*----------------------------------------------------------------------
   1.284 +-- lpx_cover_cut - generate mixed cover cut.
   1.285 +--
   1.286 +-- SYNOPSIS
   1.287 +--
   1.288 +-- #include "glplpx.h"
   1.289 +-- int lpx_cover_cut(LPX *lp, int len, int ind[], double val[],
   1.290 +--    double work[]);
   1.291 +--
   1.292 +-- DESCRIPTION
   1.293 +--
   1.294 +-- The routine lpx_cover_cut generates a mixed cover cut for a given
   1.295 +-- row of the MIP problem.
   1.296 +--
   1.297 +-- The given row of the MIP problem should be explicitly specified in
   1.298 +-- the form:
   1.299 +--
   1.300 +--    sum{j in J} a[j]*x[j] <= b.                                    (1)
   1.301 +--
   1.302 +-- On entry indices (ordinal numbers) of structural variables, which
   1.303 +-- have non-zero constraint coefficients, should be placed in locations
   1.304 +-- ind[1], ..., ind[len], and corresponding constraint coefficients
   1.305 +-- should be placed in locations val[1], ..., val[len]. The right-hand
   1.306 +-- side b should be stored in location val[0].
   1.307 +--
   1.308 +-- The working array work should have at least nb locations, where nb
   1.309 +-- is the number of binary variables in (1).
   1.310 +--
   1.311 +-- The routine generates a mixed cover cut in the same form as (1) and
   1.312 +-- stores the cut coefficients and right-hand side in the same way as
   1.313 +-- just described above.
   1.314 +--
   1.315 +-- RETURNS
   1.316 +--
   1.317 +-- If the cutting plane has been successfully generated, the routine
   1.318 +-- returns 1 <= len' <= n, which is the number of non-zero coefficients
   1.319 +-- in the inequality constraint. Otherwise, the routine returns zero. */
   1.320 +
   1.321 +static int lpx_cover_cut(LPX *lp, int len, int ind[], double val[],
   1.322 +      double work[])
   1.323 +{     int cov[1+4], j, k, nb, newlen, r;
   1.324 +      double f_min, f_max, alfa, beta, u, *x = work, y;
   1.325 +      /* substitute and remove fixed variables */
   1.326 +      newlen = 0;
   1.327 +      for (k = 1; k <= len; k++)
   1.328 +      {  j = ind[k];
   1.329 +         if (lpx_get_col_type(lp, j) == LPX_FX)
   1.330 +            val[0] -= val[k] * lpx_get_col_lb(lp, j);
   1.331 +         else
   1.332 +         {  newlen++;
   1.333 +            ind[newlen] = ind[k];
   1.334 +            val[newlen] = val[k];
   1.335 +         }
   1.336 +      }
   1.337 +      len = newlen;
   1.338 +      /* move binary variables to the beginning of the list so that
   1.339 +         elements 1, 2, ..., nb correspond to binary variables, and
   1.340 +         elements nb+1, nb+2, ..., len correspond to rest variables */
   1.341 +      nb = 0;
   1.342 +      for (k = 1; k <= len; k++)
   1.343 +      {  j = ind[k];
   1.344 +         if (lpx_get_col_kind(lp, j) == LPX_IV &&
   1.345 +             lpx_get_col_type(lp, j) == LPX_DB &&
   1.346 +             lpx_get_col_lb(lp, j) == 0.0 &&
   1.347 +             lpx_get_col_ub(lp, j) == 1.0)
   1.348 +         {  /* binary variable */
   1.349 +            int ind_k;
   1.350 +            double val_k;
   1.351 +            nb++;
   1.352 +            ind_k = ind[nb], val_k = val[nb];
   1.353 +            ind[nb] = ind[k], val[nb] = val[k];
   1.354 +            ind[k] = ind_k, val[k] = val_k;
   1.355 +         }
   1.356 +      }
   1.357 +      /* now the specified row has the form:
   1.358 +         sum a[j]*x[j] + sum a[j]*y[j] <= b,
   1.359 +         where x[j] are binary variables, y[j] are rest variables */
   1.360 +      /* at least two binary variables are needed */
   1.361 +      if (nb < 2) return 0;
   1.362 +      /* compute implied lower and upper bounds for sum a[j]*y[j] */
   1.363 +      f_min = f_max = 0.0;
   1.364 +      for (k = nb+1; k <= len; k++)
   1.365 +      {  j = ind[k];
   1.366 +         /* both bounds must be finite */
   1.367 +         if (lpx_get_col_type(lp, j) != LPX_DB) return 0;
   1.368 +         if (val[k] > 0.0)
   1.369 +         {  f_min += val[k] * lpx_get_col_lb(lp, j);
   1.370 +            f_max += val[k] * lpx_get_col_ub(lp, j);
   1.371 +         }
   1.372 +         else
   1.373 +         {  f_min += val[k] * lpx_get_col_ub(lp, j);
   1.374 +            f_max += val[k] * lpx_get_col_lb(lp, j);
   1.375 +         }
   1.376 +      }
   1.377 +      /* sum a[j]*x[j] + sum a[j]*y[j] <= b ===>
   1.378 +         sum a[j]*x[j] + (sum a[j]*y[j] - f_min) <= b - f_min ===>
   1.379 +         sum a[j]*x[j] + y <= b - f_min,
   1.380 +         where y = sum a[j]*y[j] - f_min;
   1.381 +         note that 0 <= y <= u, u = f_max - f_min */
   1.382 +      /* determine upper bound of y */
   1.383 +      u = f_max - f_min;
   1.384 +      /* determine value of y at the current point */
   1.385 +      y = 0.0;
   1.386 +      for (k = nb+1; k <= len; k++)
   1.387 +      {  j = ind[k];
   1.388 +         y += val[k] * lpx_get_col_prim(lp, j);
   1.389 +      }
   1.390 +      y -= f_min;
   1.391 +      if (y < 0.0) y = 0.0;
   1.392 +      if (y > u) y = u;
   1.393 +      /* modify the right-hand side b */
   1.394 +      val[0] -= f_min;
   1.395 +      /* now the transformed row has the form:
   1.396 +         sum a[j]*x[j] + y <= b, where 0 <= y <= u */
   1.397 +      /* determine values of x[j] at the current point */
   1.398 +      for (k = 1; k <= nb; k++)
   1.399 +      {  j = ind[k];
   1.400 +         x[k] = lpx_get_col_prim(lp, j);
   1.401 +         if (x[k] < 0.0) x[k] = 0.0;
   1.402 +         if (x[k] > 1.0) x[k] = 1.0;
   1.403 +      }
   1.404 +      /* if a[j] < 0, replace x[j] by its complement 1 - x'[j] */
   1.405 +      for (k = 1; k <= nb; k++)
   1.406 +      {  if (val[k] < 0.0)
   1.407 +         {  ind[k] = - ind[k];
   1.408 +            val[k] = - val[k];
   1.409 +            val[0] += val[k];
   1.410 +            x[k] = 1.0 - x[k];
   1.411 +         }
   1.412 +      }
   1.413 +      /* try to generate a mixed cover cut for the transformed row */
   1.414 +      r = cover(nb, val, val[0], u, x, y, cov, &alfa, &beta);
   1.415 +      if (r == 0) return 0;
   1.416 +      xassert(2 <= r && r <= 4);
   1.417 +      /* now the cut is in the form:
   1.418 +         sum{j in C} x[j] + alfa * y <= beta */
   1.419 +      /* store the right-hand side beta */
   1.420 +      ind[0] = 0, val[0] = beta;
   1.421 +      /* restore the original ordinal numbers of x[j] */
   1.422 +      for (j = 1; j <= r; j++) cov[j] = ind[cov[j]];
   1.423 +      /* store cut coefficients at binary variables complementing back
   1.424 +         the variables having negative row coefficients */
   1.425 +      xassert(r <= nb);
   1.426 +      for (k = 1; k <= r; k++)
   1.427 +      {  if (cov[k] > 0)
   1.428 +         {  ind[k] = +cov[k];
   1.429 +            val[k] = +1.0;
   1.430 +         }
   1.431 +         else
   1.432 +         {  ind[k] = -cov[k];
   1.433 +            val[k] = -1.0;
   1.434 +            val[0] -= 1.0;
   1.435 +         }
   1.436 +      }
   1.437 +      /* substitute y = sum a[j]*y[j] - f_min */
   1.438 +      for (k = nb+1; k <= len; k++)
   1.439 +      {  r++;
   1.440 +         ind[r] = ind[k];
   1.441 +         val[r] = alfa * val[k];
   1.442 +      }
   1.443 +      val[0] += alfa * f_min;
   1.444 +      xassert(r <= len);
   1.445 +      len = r;
   1.446 +      return len;
   1.447 +}
   1.448 +
   1.449 +/*----------------------------------------------------------------------
   1.450 +-- lpx_eval_row - compute explictily specified row.
   1.451 +--
   1.452 +-- SYNOPSIS
   1.453 +--
   1.454 +-- #include "glplpx.h"
   1.455 +-- double lpx_eval_row(LPX *lp, int len, int ind[], double val[]);
   1.456 +--
   1.457 +-- DESCRIPTION
   1.458 +--
   1.459 +-- The routine lpx_eval_row computes the primal value of an explicitly
   1.460 +-- specified row using current values of structural variables.
   1.461 +--
   1.462 +-- The explicitly specified row may be thought as a linear form:
   1.463 +--
   1.464 +--    y = a[1]*x[m+1] + a[2]*x[m+2] + ... + a[n]*x[m+n],
   1.465 +--
   1.466 +-- where y is an auxiliary variable for this row, a[j] are coefficients
   1.467 +-- of the linear form, x[m+j] are structural variables.
   1.468 +--
   1.469 +-- On entry column indices and numerical values of non-zero elements of
   1.470 +-- the row should be stored in locations ind[1], ..., ind[len] and
   1.471 +-- val[1], ..., val[len], where len is the number of non-zero elements.
   1.472 +-- The array ind and val are not changed on exit.
   1.473 +--
   1.474 +-- RETURNS
   1.475 +--
   1.476 +-- The routine returns a computed value of y, the auxiliary variable of
   1.477 +-- the specified row. */
   1.478 +
   1.479 +static double lpx_eval_row(LPX *lp, int len, int ind[], double val[])
   1.480 +{     int n = lpx_get_num_cols(lp);
   1.481 +      int j, k;
   1.482 +      double sum = 0.0;
   1.483 +      if (len < 0)
   1.484 +         xerror("lpx_eval_row: len = %d; invalid row length\n", len);
   1.485 +      for (k = 1; k <= len; k++)
   1.486 +      {  j = ind[k];
   1.487 +         if (!(1 <= j && j <= n))
   1.488 +            xerror("lpx_eval_row: j = %d; column number out of range\n",
   1.489 +               j);
   1.490 +         sum += val[k] * lpx_get_col_prim(lp, j);
   1.491 +      }
   1.492 +      return sum;
   1.493 +}
   1.494 +
   1.495 +/***********************************************************************
   1.496 +*  NAME
   1.497 +*
   1.498 +*  ios_cov_gen - generate mixed cover cuts
   1.499 +*
   1.500 +*  SYNOPSIS
   1.501 +*
   1.502 +*  #include "glpios.h"
   1.503 +*  void ios_cov_gen(glp_tree *tree);
   1.504 +*
   1.505 +*  DESCRIPTION
   1.506 +*
   1.507 +*  The routine ios_cov_gen generates mixed cover cuts for the current
   1.508 +*  point and adds them to the cut pool. */
   1.509 +
   1.510 +void ios_cov_gen(glp_tree *tree)
   1.511 +{     glp_prob *prob = tree->mip;
   1.512 +      int m = lpx_get_num_rows(prob);
   1.513 +      int n = lpx_get_num_cols(prob);
   1.514 +      int i, k, type, kase, len, *ind;
   1.515 +      double r, *val, *work;
   1.516 +      xassert(lpx_get_status(prob) == LPX_OPT);
   1.517 +      /* allocate working arrays */
   1.518 +      ind = xcalloc(1+n, sizeof(int));
   1.519 +      val = xcalloc(1+n, sizeof(double));
   1.520 +      work = xcalloc(1+n, sizeof(double));
   1.521 +      /* look through all rows */
   1.522 +      for (i = 1; i <= m; i++)
   1.523 +      for (kase = 1; kase <= 2; kase++)
   1.524 +      {  type = lpx_get_row_type(prob, i);
   1.525 +         if (kase == 1)
   1.526 +         {  /* consider rows of '<=' type */
   1.527 +            if (!(type == LPX_UP || type == LPX_DB)) continue;
   1.528 +            len = lpx_get_mat_row(prob, i, ind, val);
   1.529 +            val[0] = lpx_get_row_ub(prob, i);
   1.530 +         }
   1.531 +         else
   1.532 +         {  /* consider rows of '>=' type */
   1.533 +            if (!(type == LPX_LO || type == LPX_DB)) continue;
   1.534 +            len = lpx_get_mat_row(prob, i, ind, val);
   1.535 +            for (k = 1; k <= len; k++) val[k] = - val[k];
   1.536 +            val[0] = - lpx_get_row_lb(prob, i);
   1.537 +         }
   1.538 +         /* generate mixed cover cut:
   1.539 +            sum{j in J} a[j] * x[j] <= b */
   1.540 +         len = lpx_cover_cut(prob, len, ind, val, work);
   1.541 +         if (len == 0) continue;
   1.542 +         /* at the current point the cut inequality is violated, i.e.
   1.543 +            sum{j in J} a[j] * x[j] - b > 0 */
   1.544 +         r = lpx_eval_row(prob, len, ind, val) - val[0];
   1.545 +         if (r < 1e-3) continue;
   1.546 +         /* add the cut to the cut pool */
   1.547 +         glp_ios_add_row(tree, NULL, GLP_RF_COV, 0, len, ind, val,
   1.548 +            GLP_UP, val[0]);
   1.549 +      }
   1.550 +      /* free working arrays */
   1.551 +      xfree(ind);
   1.552 +      xfree(val);
   1.553 +      xfree(work);
   1.554 +      return;
   1.555 +}
   1.556 +
   1.557 +/* eof */