1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpios08.c	Sun Nov 06 20:59:10 2011 +0100
     1.3 @@ -0,0 +1,907 @@
     1.4 +/* glpios08.c (clique 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, 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 "glpios.h"
    1.29 +
    1.30 +static double get_row_lb(LPX *lp, int i)
    1.31 +{     /* this routine returns lower bound of row i or -DBL_MAX if the
    1.32 +         row has no lower bound */
    1.33 +      double lb;
    1.34 +      switch (lpx_get_row_type(lp, i))
    1.35 +      {  case LPX_FR:
    1.36 +         case LPX_UP:
    1.37 +            lb = -DBL_MAX;
    1.38 +            break;
    1.39 +         case LPX_LO:
    1.40 +         case LPX_DB:
    1.41 +         case LPX_FX:
    1.42 +            lb = lpx_get_row_lb(lp, i);
    1.43 +            break;
    1.44 +         default:
    1.45 +            xassert(lp != lp);
    1.46 +      }
    1.47 +      return lb;
    1.48 +}
    1.49 +
    1.50 +static double get_row_ub(LPX *lp, int i)
    1.51 +{     /* this routine returns upper bound of row i or +DBL_MAX if the
    1.52 +         row has no upper bound */
    1.53 +      double ub;
    1.54 +      switch (lpx_get_row_type(lp, i))
    1.55 +      {  case LPX_FR:
    1.56 +         case LPX_LO:
    1.57 +            ub = +DBL_MAX;
    1.58 +            break;
    1.59 +         case LPX_UP:
    1.60 +         case LPX_DB:
    1.61 +         case LPX_FX:
    1.62 +            ub = lpx_get_row_ub(lp, i);
    1.63 +            break;
    1.64 +         default:
    1.65 +            xassert(lp != lp);
    1.66 +      }
    1.67 +      return ub;
    1.68 +}
    1.69 +
    1.70 +static double get_col_lb(LPX *lp, int j)
    1.71 +{     /* this routine returns lower bound of column j or -DBL_MAX if
    1.72 +         the column has no lower bound */
    1.73 +      double lb;
    1.74 +      switch (lpx_get_col_type(lp, j))
    1.75 +      {  case LPX_FR:
    1.76 +         case LPX_UP:
    1.77 +            lb = -DBL_MAX;
    1.78 +            break;
    1.79 +         case LPX_LO:
    1.80 +         case LPX_DB:
    1.81 +         case LPX_FX:
    1.82 +            lb = lpx_get_col_lb(lp, j);
    1.83 +            break;
    1.84 +         default:
    1.85 +            xassert(lp != lp);
    1.86 +      }
    1.87 +      return lb;
    1.88 +}
    1.89 +
    1.90 +static double get_col_ub(LPX *lp, int j)
    1.91 +{     /* this routine returns upper bound of column j or +DBL_MAX if
    1.92 +         the column has no upper bound */
    1.93 +      double ub;
    1.94 +      switch (lpx_get_col_type(lp, j))
    1.95 +      {  case LPX_FR:
    1.96 +         case LPX_LO:
    1.97 +            ub = +DBL_MAX;
    1.98 +            break;
    1.99 +         case LPX_UP:
   1.100 +         case LPX_DB:
   1.101 +         case LPX_FX:
   1.102 +            ub = lpx_get_col_ub(lp, j);
   1.103 +            break;
   1.104 +         default:
   1.105 +            xassert(lp != lp);
   1.106 +      }
   1.107 +      return ub;
   1.108 +}
   1.109 +
   1.110 +static int is_binary(LPX *lp, int j)
   1.111 +{     /* this routine checks if variable x[j] is binary */
   1.112 +      return
   1.113 +         lpx_get_col_kind(lp, j) == LPX_IV &&
   1.114 +         lpx_get_col_type(lp, j) == LPX_DB &&
   1.115 +         lpx_get_col_lb(lp, j) == 0.0 && lpx_get_col_ub(lp, j) == 1.0;
   1.116 +}
   1.117 +
   1.118 +static double eval_lf_min(LPX *lp, int len, int ind[], double val[])
   1.119 +{     /* this routine computes the minimum of a specified linear form
   1.120 +
   1.121 +            sum a[j]*x[j]
   1.122 +             j
   1.123 +
   1.124 +         using the formula:
   1.125 +
   1.126 +            min =   sum   a[j]*lb[j] +   sum   a[j]*ub[j],
   1.127 +                  j in J+              j in J-
   1.128 +
   1.129 +         where J+ = {j: a[j] > 0}, J- = {j: a[j] < 0}, lb[j] and ub[j]
   1.130 +         are lower and upper bound of variable x[j], resp. */
   1.131 +      int j, t;
   1.132 +      double lb, ub, sum;
   1.133 +      sum = 0.0;
   1.134 +      for (t = 1; t <= len; t++)
   1.135 +      {  j = ind[t];
   1.136 +         if (val[t] > 0.0)
   1.137 +         {  lb = get_col_lb(lp, j);
   1.138 +            if (lb == -DBL_MAX)
   1.139 +            {  sum = -DBL_MAX;
   1.140 +               break;
   1.141 +            }
   1.142 +            sum += val[t] * lb;
   1.143 +         }
   1.144 +         else if (val[t] < 0.0)
   1.145 +         {  ub = get_col_ub(lp, j);
   1.146 +            if (ub == +DBL_MAX)
   1.147 +            {  sum = -DBL_MAX;
   1.148 +               break;
   1.149 +            }
   1.150 +            sum += val[t] * ub;
   1.151 +         }
   1.152 +         else
   1.153 +            xassert(val != val);
   1.154 +      }
   1.155 +      return sum;
   1.156 +}
   1.157 +
   1.158 +static double eval_lf_max(LPX *lp, int len, int ind[], double val[])
   1.159 +{     /* this routine computes the maximum of a specified linear form
   1.160 +
   1.161 +            sum a[j]*x[j]
   1.162 +             j
   1.163 +
   1.164 +         using the formula:
   1.165 +
   1.166 +            max =   sum   a[j]*ub[j] +   sum   a[j]*lb[j],
   1.167 +                  j in J+              j in J-
   1.168 +
   1.169 +         where J+ = {j: a[j] > 0}, J- = {j: a[j] < 0}, lb[j] and ub[j]
   1.170 +         are lower and upper bound of variable x[j], resp. */
   1.171 +      int j, t;
   1.172 +      double lb, ub, sum;
   1.173 +      sum = 0.0;
   1.174 +      for (t = 1; t <= len; t++)
   1.175 +      {  j = ind[t];
   1.176 +         if (val[t] > 0.0)
   1.177 +         {  ub = get_col_ub(lp, j);
   1.178 +            if (ub == +DBL_MAX)
   1.179 +            {  sum = +DBL_MAX;
   1.180 +               break;
   1.181 +            }
   1.182 +            sum += val[t] * ub;
   1.183 +         }
   1.184 +         else if (val[t] < 0.0)
   1.185 +         {  lb = get_col_lb(lp, j);
   1.186 +            if (lb == -DBL_MAX)
   1.187 +            {  sum = +DBL_MAX;
   1.188 +               break;
   1.189 +            }
   1.190 +            sum += val[t] * lb;
   1.191 +         }
   1.192 +         else
   1.193 +            xassert(val != val);
   1.194 +      }
   1.195 +      return sum;
   1.196 +}
   1.197 +
   1.198 +/*----------------------------------------------------------------------
   1.199 +-- probing - determine logical relation between binary variables.
   1.200 +--
   1.201 +-- This routine tentatively sets a binary variable to 0 and then to 1
   1.202 +-- and examines whether another binary variable is caused to be fixed.
   1.203 +--
   1.204 +-- The examination is based only on one row (constraint), which is the
   1.205 +-- following:
   1.206 +--
   1.207 +--    L <= sum a[j]*x[j] <= U.                                       (1)
   1.208 +--          j
   1.209 +--
   1.210 +-- Let x[p] be a probing variable, x[q] be an examined variable. Then
   1.211 +-- (1) can be written as:
   1.212 +--
   1.213 +--    L <=   sum  a[j]*x[j] + a[p]*x[p] + a[q]*x[q] <= U,            (2)
   1.214 +--         j in J'
   1.215 +--
   1.216 +-- where J' = {j: j != p and j != q}.
   1.217 +--
   1.218 +-- Let
   1.219 +--
   1.220 +--    L' = L - a[p]*x[p],                                            (3)
   1.221 +--
   1.222 +--    U' = U - a[p]*x[p],                                            (4)
   1.223 +--
   1.224 +-- where x[p] is assumed to be fixed at 0 or 1. So (2) can be rewritten
   1.225 +-- as follows:
   1.226 +--
   1.227 +--    L' <=   sum  a[j]*x[j] + a[q]*x[q] <= U',                      (5)
   1.228 +--          j in J'
   1.229 +--
   1.230 +-- from where we have:
   1.231 +--
   1.232 +--    L' -  sum  a[j]*x[j] <= a[q]*x[q] <= U' -  sum  a[j]*x[j].     (6)
   1.233 +--        j in J'                              j in J'
   1.234 +--
   1.235 +-- Thus,
   1.236 +--
   1.237 +--    min a[q]*x[q] = L' - MAX,                                      (7)
   1.238 +--
   1.239 +--    max a[q]*x[q] = U' - MIN,                                      (8)
   1.240 +--
   1.241 +-- where
   1.242 +--
   1.243 +--    MIN = min  sum  a[j]*x[j],                                     (9)
   1.244 +--             j in J'
   1.245 +--
   1.246 +--    MAX = max  sum  a[j]*x[j].                                    (10)
   1.247 +--             j in J'
   1.248 +--
   1.249 +-- Formulae (7) and (8) allows determining implied lower and upper
   1.250 +-- bounds of x[q].
   1.251 +--
   1.252 +-- Parameters len, val, L and U specify the constraint (1).
   1.253 +--
   1.254 +-- Parameters lf_min and lf_max specify implied lower and upper bounds
   1.255 +-- of the linear form (1). It is assumed that these bounds are computed
   1.256 +-- with the routines eval_lf_min and eval_lf_max (see above).
   1.257 +--
   1.258 +-- Parameter p specifies the probing variable x[p], which is set to 0
   1.259 +-- (if set is 0) or to 1 (if set is 1).
   1.260 +--
   1.261 +-- Parameter q specifies the examined variable x[q].
   1.262 +--
   1.263 +-- On exit the routine returns one of the following codes:
   1.264 +--
   1.265 +-- 0 - there is no logical relation between x[p] and x[q];
   1.266 +-- 1 - x[q] can take only on value 0;
   1.267 +-- 2 - x[q] can take only on value 1. */
   1.268 +
   1.269 +static int probing(int len, double val[], double L, double U,
   1.270 +      double lf_min, double lf_max, int p, int set, int q)
   1.271 +{     double temp;
   1.272 +      xassert(1 <= p && p < q && q <= len);
   1.273 +      /* compute L' (3) */
   1.274 +      if (L != -DBL_MAX && set) L -= val[p];
   1.275 +      /* compute U' (4) */
   1.276 +      if (U != +DBL_MAX && set) U -= val[p];
   1.277 +      /* compute MIN (9) */
   1.278 +      if (lf_min != -DBL_MAX)
   1.279 +      {  if (val[p] < 0.0) lf_min -= val[p];
   1.280 +         if (val[q] < 0.0) lf_min -= val[q];
   1.281 +      }
   1.282 +      /* compute MAX (10) */
   1.283 +      if (lf_max != +DBL_MAX)
   1.284 +      {  if (val[p] > 0.0) lf_max -= val[p];
   1.285 +         if (val[q] > 0.0) lf_max -= val[q];
   1.286 +      }
   1.287 +      /* compute implied lower bound of x[q]; see (7), (8) */
   1.288 +      if (val[q] > 0.0)
   1.289 +      {  if (L == -DBL_MAX || lf_max == +DBL_MAX)
   1.290 +            temp = -DBL_MAX;
   1.291 +         else
   1.292 +            temp = (L - lf_max) / val[q];
   1.293 +      }
   1.294 +      else
   1.295 +      {  if (U == +DBL_MAX || lf_min == -DBL_MAX)
   1.296 +            temp = -DBL_MAX;
   1.297 +         else
   1.298 +            temp = (U - lf_min) / val[q];
   1.299 +      }
   1.300 +      if (temp > 0.001) return 2;
   1.301 +      /* compute implied upper bound of x[q]; see (7), (8) */
   1.302 +      if (val[q] > 0.0)
   1.303 +      {  if (U == +DBL_MAX || lf_min == -DBL_MAX)
   1.304 +            temp = +DBL_MAX;
   1.305 +         else
   1.306 +            temp = (U - lf_min) / val[q];
   1.307 +      }
   1.308 +      else
   1.309 +      {  if (L == -DBL_MAX || lf_max == +DBL_MAX)
   1.310 +            temp = +DBL_MAX;
   1.311 +         else
   1.312 +            temp = (L - lf_max) / val[q];
   1.313 +      }
   1.314 +      if (temp < 0.999) return 1;
   1.315 +      /* there is no logical relation between x[p] and x[q] */
   1.316 +      return 0;
   1.317 +}
   1.318 +
   1.319 +struct COG
   1.320 +{     /* conflict graph; it represents logical relations between binary
   1.321 +         variables and has a vertex for each binary variable and its
   1.322 +         complement, and an edge between two vertices when at most one
   1.323 +         of the variables represented by the vertices can equal one in
   1.324 +         an optimal solution */
   1.325 +      int n;
   1.326 +      /* number of variables */
   1.327 +      int nb;
   1.328 +      /* number of binary variables represented in the graph (note that
   1.329 +         not all binary variables can be represented); vertices which
   1.330 +         correspond to binary variables have numbers 1, ..., nb while
   1.331 +         vertices which correspond to complements of binary variables
   1.332 +         have numbers nb+1, ..., nb+nb */
   1.333 +      int ne;
   1.334 +      /* number of edges in the graph */
   1.335 +      int *vert; /* int vert[1+n]; */
   1.336 +      /* if x[j] is a binary variable represented in the graph, vert[j]
   1.337 +         is the vertex number corresponding to x[j]; otherwise vert[j]
   1.338 +         is zero */
   1.339 +      int *orig; /* int list[1:nb]; */
   1.340 +      /* if vert[j] = k > 0, then orig[k] = j */
   1.341 +      unsigned char *a;
   1.342 +      /* adjacency matrix of the graph having 2*nb rows and columns;
   1.343 +         only strict lower triangle is stored in dense packed form */
   1.344 +};
   1.345 +
   1.346 +/*----------------------------------------------------------------------
   1.347 +-- lpx_create_cog - create the conflict graph.
   1.348 +--
   1.349 +-- SYNOPSIS
   1.350 +--
   1.351 +-- #include "glplpx.h"
   1.352 +-- void *lpx_create_cog(LPX *lp);
   1.353 +--
   1.354 +-- DESCRIPTION
   1.355 +--
   1.356 +-- The routine lpx_create_cog creates the conflict graph for a given
   1.357 +-- problem instance.
   1.358 +--
   1.359 +-- RETURNS
   1.360 +--
   1.361 +-- If the graph has been created, the routine returns a pointer to it.
   1.362 +-- Otherwise the routine returns NULL. */
   1.363 +
   1.364 +#define MAX_NB 4000
   1.365 +#define MAX_ROW_LEN 500
   1.366 +
   1.367 +static void lpx_add_cog_edge(void *_cog, int i, int j);
   1.368 +
   1.369 +static void *lpx_create_cog(LPX *lp)
   1.370 +{     struct COG *cog = NULL;
   1.371 +      int m, n, nb, i, j, p, q, len, *ind, *vert, *orig;
   1.372 +      double L, U, lf_min, lf_max, *val;
   1.373 +      xprintf("Creating the conflict graph...\n");
   1.374 +      m = lpx_get_num_rows(lp);
   1.375 +      n = lpx_get_num_cols(lp);
   1.376 +      /* determine which binary variables should be included in the
   1.377 +         conflict graph */
   1.378 +      nb = 0;
   1.379 +      vert = xcalloc(1+n, sizeof(int));
   1.380 +      for (j = 1; j <= n; j++) vert[j] = 0;
   1.381 +      orig = xcalloc(1+n, sizeof(int));
   1.382 +      ind = xcalloc(1+n, sizeof(int));
   1.383 +      val = xcalloc(1+n, sizeof(double));
   1.384 +      for (i = 1; i <= m; i++)
   1.385 +      {  L = get_row_lb(lp, i);
   1.386 +         U = get_row_ub(lp, i);
   1.387 +         if (L == -DBL_MAX && U == +DBL_MAX) continue;
   1.388 +         len = lpx_get_mat_row(lp, i, ind, val);
   1.389 +         if (len > MAX_ROW_LEN) continue;
   1.390 +         lf_min = eval_lf_min(lp, len, ind, val);
   1.391 +         lf_max = eval_lf_max(lp, len, ind, val);
   1.392 +         for (p = 1; p <= len; p++)
   1.393 +         {  if (!is_binary(lp, ind[p])) continue;
   1.394 +            for (q = p+1; q <= len; q++)
   1.395 +            {  if (!is_binary(lp, ind[q])) continue;
   1.396 +               if (probing(len, val, L, U, lf_min, lf_max, p, 0, q) ||
   1.397 +                   probing(len, val, L, U, lf_min, lf_max, p, 1, q))
   1.398 +               {  /* there is a logical relation */
   1.399 +                  /* include the first variable in the graph */
   1.400 +                  j = ind[p];
   1.401 +                  if (vert[j] == 0) nb++, vert[j] = nb, orig[nb] = j;
   1.402 +                  /* incude the second variable in the graph */
   1.403 +                  j = ind[q];
   1.404 +                  if (vert[j] == 0) nb++, vert[j] = nb, orig[nb] = j;
   1.405 +               }
   1.406 +            }
   1.407 +         }
   1.408 +      }
   1.409 +      /* if the graph is either empty or has too many vertices, do not
   1.410 +         create it */
   1.411 +      if (nb == 0 || nb > MAX_NB)
   1.412 +      {  xprintf("The conflict graph is either empty or too big\n");
   1.413 +         xfree(vert);
   1.414 +         xfree(orig);
   1.415 +         goto done;
   1.416 +      }
   1.417 +      /* create the conflict graph */
   1.418 +      cog = xmalloc(sizeof(struct COG));
   1.419 +      cog->n = n;
   1.420 +      cog->nb = nb;
   1.421 +      cog->ne = 0;
   1.422 +      cog->vert = vert;
   1.423 +      cog->orig = orig;
   1.424 +      len = nb + nb; /* number of vertices */
   1.425 +      len = (len * (len - 1)) / 2; /* number of entries in triangle */
   1.426 +      len = (len + (CHAR_BIT - 1)) / CHAR_BIT; /* bytes needed */
   1.427 +      cog->a = xmalloc(len);
   1.428 +      memset(cog->a, 0, len);
   1.429 +      for (j = 1; j <= nb; j++)
   1.430 +      {  /* add edge between variable and its complement */
   1.431 +         lpx_add_cog_edge(cog, +orig[j], -orig[j]);
   1.432 +      }
   1.433 +      for (i = 1; i <= m; i++)
   1.434 +      {  L = get_row_lb(lp, i);
   1.435 +         U = get_row_ub(lp, i);
   1.436 +         if (L == -DBL_MAX && U == +DBL_MAX) continue;
   1.437 +         len = lpx_get_mat_row(lp, i, ind, val);
   1.438 +         if (len > MAX_ROW_LEN) continue;
   1.439 +         lf_min = eval_lf_min(lp, len, ind, val);
   1.440 +         lf_max = eval_lf_max(lp, len, ind, val);
   1.441 +         for (p = 1; p <= len; p++)
   1.442 +         {  if (!is_binary(lp, ind[p])) continue;
   1.443 +            for (q = p+1; q <= len; q++)
   1.444 +            {  if (!is_binary(lp, ind[q])) continue;
   1.445 +               /* set x[p] to 0 and examine x[q] */
   1.446 +               switch (probing(len, val, L, U, lf_min, lf_max, p, 0, q))
   1.447 +               {  case 0:
   1.448 +                     /* no logical relation */
   1.449 +                     break;
   1.450 +                  case 1:
   1.451 +                     /* x[p] = 0 implies x[q] = 0 */
   1.452 +                     lpx_add_cog_edge(cog, -ind[p], +ind[q]);
   1.453 +                     break;
   1.454 +                  case 2:
   1.455 +                     /* x[p] = 0 implies x[q] = 1 */
   1.456 +                     lpx_add_cog_edge(cog, -ind[p], -ind[q]);
   1.457 +                     break;
   1.458 +                  default:
   1.459 +                     xassert(lp != lp);
   1.460 +               }
   1.461 +               /* set x[p] to 1 and examine x[q] */
   1.462 +               switch (probing(len, val, L, U, lf_min, lf_max, p, 1, q))
   1.463 +               {  case 0:
   1.464 +                     /* no logical relation */
   1.465 +                     break;
   1.466 +                  case 1:
   1.467 +                     /* x[p] = 1 implies x[q] = 0 */
   1.468 +                     lpx_add_cog_edge(cog, +ind[p], +ind[q]);
   1.469 +                     break;
   1.470 +                  case 2:
   1.471 +                     /* x[p] = 1 implies x[q] = 1 */
   1.472 +                     lpx_add_cog_edge(cog, +ind[p], -ind[q]);
   1.473 +                     break;
   1.474 +                  default:
   1.475 +                     xassert(lp != lp);
   1.476 +               }
   1.477 +            }
   1.478 +         }
   1.479 +      }
   1.480 +      xprintf("The conflict graph has 2*%d vertices and %d edges\n",
   1.481 +         cog->nb, cog->ne);
   1.482 +done: xfree(ind);
   1.483 +      xfree(val);
   1.484 +      return cog;
   1.485 +}
   1.486 +
   1.487 +/*----------------------------------------------------------------------
   1.488 +-- lpx_add_cog_edge - add edge to the conflict graph.
   1.489 +--
   1.490 +-- SYNOPSIS
   1.491 +--
   1.492 +-- #include "glplpx.h"
   1.493 +-- void lpx_add_cog_edge(void *cog, int i, int j);
   1.494 +--
   1.495 +-- DESCRIPTION
   1.496 +--
   1.497 +-- The routine lpx_add_cog_edge adds an edge to the conflict graph.
   1.498 +-- The edge connects x[i] (if i > 0) or its complement (if i < 0) and
   1.499 +-- x[j] (if j > 0) or its complement (if j < 0), where i and j are
   1.500 +-- original ordinal numbers of corresponding variables. */
   1.501 +
   1.502 +static void lpx_add_cog_edge(void *_cog, int i, int j)
   1.503 +{     struct COG *cog = _cog;
   1.504 +      int k;
   1.505 +      xassert(i != j);
   1.506 +      /* determine indices of corresponding vertices */
   1.507 +      if (i > 0)
   1.508 +      {  xassert(1 <= i && i <= cog->n);
   1.509 +         i = cog->vert[i];
   1.510 +         xassert(i != 0);
   1.511 +      }
   1.512 +      else
   1.513 +      {  i = -i;
   1.514 +         xassert(1 <= i && i <= cog->n);
   1.515 +         i = cog->vert[i];
   1.516 +         xassert(i != 0);
   1.517 +         i += cog->nb;
   1.518 +      }
   1.519 +      if (j > 0)
   1.520 +      {  xassert(1 <= j && j <= cog->n);
   1.521 +         j = cog->vert[j];
   1.522 +         xassert(j != 0);
   1.523 +      }
   1.524 +      else
   1.525 +      {  j = -j;
   1.526 +         xassert(1 <= j && j <= cog->n);
   1.527 +         j = cog->vert[j];
   1.528 +         xassert(j != 0);
   1.529 +         j += cog->nb;
   1.530 +      }
   1.531 +      /* only lower triangle is stored, so we need i > j */
   1.532 +      if (i < j) k = i, i = j, j = k;
   1.533 +      k = ((i - 1) * (i - 2)) / 2 + (j - 1);
   1.534 +      cog->a[k / CHAR_BIT] |=
   1.535 +         (unsigned char)(1 << ((CHAR_BIT - 1) - k % CHAR_BIT));
   1.536 +      cog->ne++;
   1.537 +      return;
   1.538 +}
   1.539 +
   1.540 +/*----------------------------------------------------------------------
   1.541 +-- MAXIMUM WEIGHT CLIQUE
   1.542 +--
   1.543 +-- Two subroutines sub() and wclique() below are intended to find a
   1.544 +-- maximum weight clique in a given undirected graph. These subroutines
   1.545 +-- are slightly modified version of the program WCLIQUE developed by
   1.546 +-- Patric Ostergard <http://www.tcs.hut.fi/~pat/wclique.html> and based
   1.547 +-- on ideas from the article "P. R. J. Ostergard, A new algorithm for
   1.548 +-- the maximum-weight clique problem, submitted for publication", which
   1.549 +-- in turn is a generalization of the algorithm for unweighted graphs
   1.550 +-- presented in "P. R. J. Ostergard, A fast algorithm for the maximum
   1.551 +-- clique problem, submitted for publication".
   1.552 +--
   1.553 +-- USED WITH PERMISSION OF THE AUTHOR OF THE ORIGINAL CODE. */
   1.554 +
   1.555 +struct dsa
   1.556 +{     /* dynamic storage area */
   1.557 +      int n;
   1.558 +      /* number of vertices */
   1.559 +      int *wt; /* int wt[0:n-1]; */
   1.560 +      /* weights */
   1.561 +      unsigned char *a;
   1.562 +      /* adjacency matrix (packed lower triangle without main diag.) */
   1.563 +      int record;
   1.564 +      /* weight of best clique */
   1.565 +      int rec_level;
   1.566 +      /* number of vertices in best clique */
   1.567 +      int *rec; /* int rec[0:n-1]; */
   1.568 +      /* best clique so far */
   1.569 +      int *clique; /* int clique[0:n-1]; */
   1.570 +      /* table for pruning */
   1.571 +      int *set; /* int set[0:n-1]; */
   1.572 +      /* current clique */
   1.573 +};
   1.574 +
   1.575 +#define n         (dsa->n)
   1.576 +#define wt        (dsa->wt)
   1.577 +#define a         (dsa->a)
   1.578 +#define record    (dsa->record)
   1.579 +#define rec_level (dsa->rec_level)
   1.580 +#define rec       (dsa->rec)
   1.581 +#define clique    (dsa->clique)
   1.582 +#define set       (dsa->set)
   1.583 +
   1.584 +#if 0
   1.585 +static int is_edge(struct dsa *dsa, int i, int j)
   1.586 +{     /* if there is arc (i,j), the routine returns true; otherwise
   1.587 +         false; 0 <= i, j < n */
   1.588 +      int k;
   1.589 +      xassert(0 <= i && i < n);
   1.590 +      xassert(0 <= j && j < n);
   1.591 +      if (i == j) return 0;
   1.592 +      if (i < j) k = i, i = j, j = k;
   1.593 +      k = (i * (i - 1)) / 2 + j;
   1.594 +      return a[k / CHAR_BIT] &
   1.595 +         (unsigned char)(1 << ((CHAR_BIT - 1) - k % CHAR_BIT));
   1.596 +}
   1.597 +#else
   1.598 +#define is_edge(dsa, i, j) ((i) == (j) ? 0 : \
   1.599 +      (i) > (j) ? is_edge1(i, j) : is_edge1(j, i))
   1.600 +#define is_edge1(i, j) is_edge2(((i) * ((i) - 1)) / 2 + (j))
   1.601 +#define is_edge2(k) (a[(k) / CHAR_BIT] & \
   1.602 +      (unsigned char)(1 << ((CHAR_BIT - 1) - (k) % CHAR_BIT)))
   1.603 +#endif
   1.604 +
   1.605 +static void sub(struct dsa *dsa, int ct, int table[], int level,
   1.606 +      int weight, int l_weight)
   1.607 +{     int i, j, k, curr_weight, left_weight, *p1, *p2, *newtable;
   1.608 +      newtable = xcalloc(n, sizeof(int));
   1.609 +      if (ct <= 0)
   1.610 +      {  /* 0 or 1 elements left; include these */
   1.611 +         if (ct == 0)
   1.612 +         {  set[level++] = table[0];
   1.613 +            weight += l_weight;
   1.614 +         }
   1.615 +         if (weight > record)
   1.616 +         {  record = weight;
   1.617 +            rec_level = level;
   1.618 +            for (i = 0; i < level; i++) rec[i] = set[i];
   1.619 +         }
   1.620 +         goto done;
   1.621 +      }
   1.622 +      for (i = ct; i >= 0; i--)
   1.623 +      {  if ((level == 0) && (i < ct)) goto done;
   1.624 +         k = table[i];
   1.625 +         if ((level > 0) && (clique[k] <= (record - weight)))
   1.626 +            goto done; /* prune */
   1.627 +         set[level] = k;
   1.628 +         curr_weight = weight + wt[k];
   1.629 +         l_weight -= wt[k];
   1.630 +         if (l_weight <= (record - curr_weight))
   1.631 +            goto done; /* prune */
   1.632 +         p1 = newtable;
   1.633 +         p2 = table;
   1.634 +         left_weight = 0;
   1.635 +         while (p2 < table + i)
   1.636 +         {  j = *p2++;
   1.637 +            if (is_edge(dsa, j, k))
   1.638 +            {  *p1++ = j;
   1.639 +               left_weight += wt[j];
   1.640 +            }
   1.641 +         }
   1.642 +         if (left_weight <= (record - curr_weight)) continue;
   1.643 +         sub(dsa, p1 - newtable - 1, newtable, level + 1, curr_weight,
   1.644 +            left_weight);
   1.645 +      }
   1.646 +done: xfree(newtable);
   1.647 +      return;
   1.648 +}
   1.649 +
   1.650 +static int wclique(int _n, int w[], unsigned char _a[], int sol[])
   1.651 +{     struct dsa _dsa, *dsa = &_dsa;
   1.652 +      int i, j, p, max_wt, max_nwt, wth, *used, *nwt, *pos;
   1.653 +      glp_long timer;
   1.654 +      n = _n;
   1.655 +      wt = &w[1];
   1.656 +      a = _a;
   1.657 +      record = 0;
   1.658 +      rec_level = 0;
   1.659 +      rec = &sol[1];
   1.660 +      clique = xcalloc(n, sizeof(int));
   1.661 +      set = xcalloc(n, sizeof(int));
   1.662 +      used = xcalloc(n, sizeof(int));
   1.663 +      nwt = xcalloc(n, sizeof(int));
   1.664 +      pos = xcalloc(n, sizeof(int));
   1.665 +      /* start timer */
   1.666 +      timer = xtime();
   1.667 +      /* order vertices */
   1.668 +      for (i = 0; i < n; i++)
   1.669 +      {  nwt[i] = 0;
   1.670 +         for (j = 0; j < n; j++)
   1.671 +            if (is_edge(dsa, i, j)) nwt[i] += wt[j];
   1.672 +      }
   1.673 +      for (i = 0; i < n; i++)
   1.674 +         used[i] = 0;
   1.675 +      for (i = n-1; i >= 0; i--)
   1.676 +      {  max_wt = -1;
   1.677 +         max_nwt = -1;
   1.678 +         for (j = 0; j < n; j++)
   1.679 +         {  if ((!used[j]) && ((wt[j] > max_wt) || (wt[j] == max_wt
   1.680 +               && nwt[j] > max_nwt)))
   1.681 +            {  max_wt = wt[j];
   1.682 +               max_nwt = nwt[j];
   1.683 +               p = j;
   1.684 +            }
   1.685 +         }
   1.686 +         pos[i] = p;
   1.687 +         used[p] = 1;
   1.688 +         for (j = 0; j < n; j++)
   1.689 +            if ((!used[j]) && (j != p) && (is_edge(dsa, p, j)))
   1.690 +               nwt[j] -= wt[p];
   1.691 +      }
   1.692 +      /* main routine */
   1.693 +      wth = 0;
   1.694 +      for (i = 0; i < n; i++)
   1.695 +      {  wth += wt[pos[i]];
   1.696 +         sub(dsa, i, pos, 0, 0, wth);
   1.697 +         clique[pos[i]] = record;
   1.698 +#if 0
   1.699 +         if (utime() >= timer + 5.0)
   1.700 +#else
   1.701 +         if (xdifftime(xtime(), timer) >= 5.0 - 0.001)
   1.702 +#endif
   1.703 +         {  /* print current record and reset timer */
   1.704 +            xprintf("level = %d (%d); best = %d\n", i+1, n, record);
   1.705 +#if 0
   1.706 +            timer = utime();
   1.707 +#else
   1.708 +            timer = xtime();
   1.709 +#endif
   1.710 +         }
   1.711 +      }
   1.712 +      xfree(clique);
   1.713 +      xfree(set);
   1.714 +      xfree(used);
   1.715 +      xfree(nwt);
   1.716 +      xfree(pos);
   1.717 +      /* return the solution found */
   1.718 +      for (i = 1; i <= rec_level; i++) sol[i]++;
   1.719 +      return rec_level;
   1.720 +}
   1.721 +
   1.722 +#undef n
   1.723 +#undef wt
   1.724 +#undef a
   1.725 +#undef record
   1.726 +#undef rec_level
   1.727 +#undef rec
   1.728 +#undef clique
   1.729 +#undef set
   1.730 +
   1.731 +/*----------------------------------------------------------------------
   1.732 +-- lpx_clique_cut - generate cluque cut.
   1.733 +--
   1.734 +-- SYNOPSIS
   1.735 +--
   1.736 +-- #include "glplpx.h"
   1.737 +-- int lpx_clique_cut(LPX *lp, void *cog, int ind[], double val[]);
   1.738 +--
   1.739 +-- DESCRIPTION
   1.740 +--
   1.741 +-- The routine lpx_clique_cut generates a clique cut using the conflict
   1.742 +-- graph specified by the parameter cog.
   1.743 +--
   1.744 +-- If a violated clique cut has been found, it has the following form:
   1.745 +--
   1.746 +--    sum{j in J} a[j]*x[j] <= b.
   1.747 +--
   1.748 +-- Variable indices j in J are stored in elements ind[1], ..., ind[len]
   1.749 +-- while corresponding constraint coefficients are stored in elements
   1.750 +-- val[1], ..., val[len], where len is returned on exit. The right-hand
   1.751 +-- side b is stored in element val[0].
   1.752 +--
   1.753 +-- RETURNS
   1.754 +--
   1.755 +-- If the cutting plane has been successfully generated, the routine
   1.756 +-- returns 1 <= len <= n, which is the number of non-zero coefficients
   1.757 +-- in the inequality constraint. Otherwise, the routine returns zero. */
   1.758 +
   1.759 +static int lpx_clique_cut(LPX *lp, void *_cog, int ind[], double val[])
   1.760 +{     struct COG *cog = _cog;
   1.761 +      int n = lpx_get_num_cols(lp);
   1.762 +      int j, t, v, card, temp, len = 0, *w, *sol;
   1.763 +      double x, sum, b, *vec;
   1.764 +      /* allocate working arrays */
   1.765 +      w = xcalloc(1 + 2 * cog->nb, sizeof(int));
   1.766 +      sol = xcalloc(1 + 2 * cog->nb, sizeof(int));
   1.767 +      vec = xcalloc(1+n, sizeof(double));
   1.768 +      /* assign weights to vertices of the conflict graph */
   1.769 +      for (t = 1; t <= cog->nb; t++)
   1.770 +      {  j = cog->orig[t];
   1.771 +         x = lpx_get_col_prim(lp, j);
   1.772 +         temp = (int)(100.0 * x + 0.5);
   1.773 +         if (temp < 0) temp = 0;
   1.774 +         if (temp > 100) temp = 100;
   1.775 +         w[t] = temp;
   1.776 +         w[cog->nb + t] = 100 - temp;
   1.777 +      }
   1.778 +      /* find a clique of maximum weight */
   1.779 +      card = wclique(2 * cog->nb, w, cog->a, sol);
   1.780 +      /* compute the clique weight for unscaled values */
   1.781 +      sum = 0.0;
   1.782 +      for ( t = 1; t <= card; t++)
   1.783 +      {  v = sol[t];
   1.784 +         xassert(1 <= v && v <= 2 * cog->nb);
   1.785 +         if (v <= cog->nb)
   1.786 +         {  /* vertex v corresponds to binary variable x[j] */
   1.787 +            j = cog->orig[v];
   1.788 +            x = lpx_get_col_prim(lp, j);
   1.789 +            sum += x;
   1.790 +         }
   1.791 +         else
   1.792 +         {  /* vertex v corresponds to the complement of x[j] */
   1.793 +            j = cog->orig[v - cog->nb];
   1.794 +            x = lpx_get_col_prim(lp, j);
   1.795 +            sum += 1.0 - x;
   1.796 +         }
   1.797 +      }
   1.798 +      /* if the sum of binary variables and their complements in the
   1.799 +         clique greater than 1, the clique cut is violated */
   1.800 +      if (sum >= 1.01)
   1.801 +      {  /* construct the inquality */
   1.802 +         for (j = 1; j <= n; j++) vec[j] = 0;
   1.803 +         b = 1.0;
   1.804 +         for (t = 1; t <= card; t++)
   1.805 +         {  v = sol[t];
   1.806 +            if (v <= cog->nb)
   1.807 +            {  /* vertex v corresponds to binary variable x[j] */
   1.808 +               j = cog->orig[v];
   1.809 +               xassert(1 <= j && j <= n);
   1.810 +               vec[j] += 1.0;
   1.811 +            }
   1.812 +            else
   1.813 +            {  /* vertex v corresponds to the complement of x[j] */
   1.814 +               j = cog->orig[v - cog->nb];
   1.815 +               xassert(1 <= j && j <= n);
   1.816 +               vec[j] -= 1.0;
   1.817 +               b -= 1.0;
   1.818 +            }
   1.819 +         }
   1.820 +         xassert(len == 0);
   1.821 +         for (j = 1; j <= n; j++)
   1.822 +         {  if (vec[j] != 0.0)
   1.823 +            {  len++;
   1.824 +               ind[len] = j, val[len] = vec[j];
   1.825 +            }
   1.826 +         }
   1.827 +         ind[0] = 0, val[0] = b;
   1.828 +      }
   1.829 +      /* free working arrays */
   1.830 +      xfree(w);
   1.831 +      xfree(sol);
   1.832 +      xfree(vec);
   1.833 +      /* return to the calling program */
   1.834 +      return len;
   1.835 +}
   1.836 +
   1.837 +/*----------------------------------------------------------------------
   1.838 +-- lpx_delete_cog - delete the conflict graph.
   1.839 +--
   1.840 +-- SYNOPSIS
   1.841 +--
   1.842 +-- #include "glplpx.h"
   1.843 +-- void lpx_delete_cog(void *cog);
   1.844 +--
   1.845 +-- DESCRIPTION
   1.846 +--
   1.847 +-- The routine lpx_delete_cog deletes the conflict graph, which the
   1.848 +-- parameter cog points to, freeing all the memory allocated to this
   1.849 +-- object. */
   1.850 +
   1.851 +static void lpx_delete_cog(void *_cog)
   1.852 +{     struct COG *cog = _cog;
   1.853 +      xfree(cog->vert);
   1.854 +      xfree(cog->orig);
   1.855 +      xfree(cog->a);
   1.856 +      xfree(cog);
   1.857 +}
   1.858 +
   1.859 +/**********************************************************************/
   1.860 +
   1.861 +void *ios_clq_init(glp_tree *tree)
   1.862 +{     /* initialize clique cut generator */
   1.863 +      glp_prob *mip = tree->mip;
   1.864 +      xassert(mip != NULL);
   1.865 +      return lpx_create_cog(mip);
   1.866 +}
   1.867 +
   1.868 +/***********************************************************************
   1.869 +*  NAME
   1.870 +*
   1.871 +*  ios_clq_gen - generate clique cuts
   1.872 +*
   1.873 +*  SYNOPSIS
   1.874 +*
   1.875 +*  #include "glpios.h"
   1.876 +*  void ios_clq_gen(glp_tree *tree, void *gen);
   1.877 +*
   1.878 +*  DESCRIPTION
   1.879 +*
   1.880 +*  The routine ios_clq_gen generates clique cuts for the current point
   1.881 +*  and adds them to the clique pool. */
   1.882 +
   1.883 +void ios_clq_gen(glp_tree *tree, void *gen)
   1.884 +{     int n = lpx_get_num_cols(tree->mip);
   1.885 +      int len, *ind;
   1.886 +      double *val;
   1.887 +      xassert(gen != NULL);
   1.888 +      ind = xcalloc(1+n, sizeof(int));
   1.889 +      val = xcalloc(1+n, sizeof(double));
   1.890 +      len = lpx_clique_cut(tree->mip, gen, ind, val);
   1.891 +      if (len > 0)
   1.892 +      {  /* xprintf("len = %d\n", len); */
   1.893 +         glp_ios_add_row(tree, NULL, GLP_RF_CLQ, 0, len, ind, val,
   1.894 +            GLP_UP, val[0]);
   1.895 +      }
   1.896 +      xfree(ind);
   1.897 +      xfree(val);
   1.898 +      return;
   1.899 +}
   1.900 +
   1.901 +/**********************************************************************/
   1.902 +
   1.903 +void ios_clq_term(void *gen)
   1.904 +{     /* terminate clique cut generator */
   1.905 +      xassert(gen != NULL);
   1.906 +      lpx_delete_cog(gen);
   1.907 +      return;
   1.908 +}
   1.909 +
   1.910 +/* eof */