lemon-project-template-glpk

annotate deps/glpk/src/glpios08.c @ 11:4fc6ad2fb8a6

find changesets by author, revision, files, or words in the commit message
Test GLPK in src/main.cc
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 21:43:29 +0100
parents
children
rev   line source
alpar@9 1 /* glpios08.c (clique cut generator) */
alpar@9 2
alpar@9 3 /***********************************************************************
alpar@9 4 * This code is part of GLPK (GNU Linear Programming Kit).
alpar@9 5 *
alpar@9 6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@9 7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
alpar@9 8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@9 9 * E-mail: <mao@gnu.org>.
alpar@9 10 *
alpar@9 11 * GLPK is free software: you can redistribute it and/or modify it
alpar@9 12 * under the terms of the GNU General Public License as published by
alpar@9 13 * the Free Software Foundation, either version 3 of the License, or
alpar@9 14 * (at your option) any later version.
alpar@9 15 *
alpar@9 16 * GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@9 17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@9 18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@9 19 * License for more details.
alpar@9 20 *
alpar@9 21 * You should have received a copy of the GNU General Public License
alpar@9 22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@9 23 ***********************************************************************/
alpar@9 24
alpar@9 25 #include "glpios.h"
alpar@9 26
alpar@9 27 static double get_row_lb(LPX *lp, int i)
alpar@9 28 { /* this routine returns lower bound of row i or -DBL_MAX if the
alpar@9 29 row has no lower bound */
alpar@9 30 double lb;
alpar@9 31 switch (lpx_get_row_type(lp, i))
alpar@9 32 { case LPX_FR:
alpar@9 33 case LPX_UP:
alpar@9 34 lb = -DBL_MAX;
alpar@9 35 break;
alpar@9 36 case LPX_LO:
alpar@9 37 case LPX_DB:
alpar@9 38 case LPX_FX:
alpar@9 39 lb = lpx_get_row_lb(lp, i);
alpar@9 40 break;
alpar@9 41 default:
alpar@9 42 xassert(lp != lp);
alpar@9 43 }
alpar@9 44 return lb;
alpar@9 45 }
alpar@9 46
alpar@9 47 static double get_row_ub(LPX *lp, int i)
alpar@9 48 { /* this routine returns upper bound of row i or +DBL_MAX if the
alpar@9 49 row has no upper bound */
alpar@9 50 double ub;
alpar@9 51 switch (lpx_get_row_type(lp, i))
alpar@9 52 { case LPX_FR:
alpar@9 53 case LPX_LO:
alpar@9 54 ub = +DBL_MAX;
alpar@9 55 break;
alpar@9 56 case LPX_UP:
alpar@9 57 case LPX_DB:
alpar@9 58 case LPX_FX:
alpar@9 59 ub = lpx_get_row_ub(lp, i);
alpar@9 60 break;
alpar@9 61 default:
alpar@9 62 xassert(lp != lp);
alpar@9 63 }
alpar@9 64 return ub;
alpar@9 65 }
alpar@9 66
alpar@9 67 static double get_col_lb(LPX *lp, int j)
alpar@9 68 { /* this routine returns lower bound of column j or -DBL_MAX if
alpar@9 69 the column has no lower bound */
alpar@9 70 double lb;
alpar@9 71 switch (lpx_get_col_type(lp, j))
alpar@9 72 { case LPX_FR:
alpar@9 73 case LPX_UP:
alpar@9 74 lb = -DBL_MAX;
alpar@9 75 break;
alpar@9 76 case LPX_LO:
alpar@9 77 case LPX_DB:
alpar@9 78 case LPX_FX:
alpar@9 79 lb = lpx_get_col_lb(lp, j);
alpar@9 80 break;
alpar@9 81 default:
alpar@9 82 xassert(lp != lp);
alpar@9 83 }
alpar@9 84 return lb;
alpar@9 85 }
alpar@9 86
alpar@9 87 static double get_col_ub(LPX *lp, int j)
alpar@9 88 { /* this routine returns upper bound of column j or +DBL_MAX if
alpar@9 89 the column has no upper bound */
alpar@9 90 double ub;
alpar@9 91 switch (lpx_get_col_type(lp, j))
alpar@9 92 { case LPX_FR:
alpar@9 93 case LPX_LO:
alpar@9 94 ub = +DBL_MAX;
alpar@9 95 break;
alpar@9 96 case LPX_UP:
alpar@9 97 case LPX_DB:
alpar@9 98 case LPX_FX:
alpar@9 99 ub = lpx_get_col_ub(lp, j);
alpar@9 100 break;
alpar@9 101 default:
alpar@9 102 xassert(lp != lp);
alpar@9 103 }
alpar@9 104 return ub;
alpar@9 105 }
alpar@9 106
alpar@9 107 static int is_binary(LPX *lp, int j)
alpar@9 108 { /* this routine checks if variable x[j] is binary */
alpar@9 109 return
alpar@9 110 lpx_get_col_kind(lp, j) == LPX_IV &&
alpar@9 111 lpx_get_col_type(lp, j) == LPX_DB &&
alpar@9 112 lpx_get_col_lb(lp, j) == 0.0 && lpx_get_col_ub(lp, j) == 1.0;
alpar@9 113 }
alpar@9 114
alpar@9 115 static double eval_lf_min(LPX *lp, int len, int ind[], double val[])
alpar@9 116 { /* this routine computes the minimum of a specified linear form
alpar@9 117
alpar@9 118 sum a[j]*x[j]
alpar@9 119 j
alpar@9 120
alpar@9 121 using the formula:
alpar@9 122
alpar@9 123 min = sum a[j]*lb[j] + sum a[j]*ub[j],
alpar@9 124 j in J+ j in J-
alpar@9 125
alpar@9 126 where J+ = {j: a[j] > 0}, J- = {j: a[j] < 0}, lb[j] and ub[j]
alpar@9 127 are lower and upper bound of variable x[j], resp. */
alpar@9 128 int j, t;
alpar@9 129 double lb, ub, sum;
alpar@9 130 sum = 0.0;
alpar@9 131 for (t = 1; t <= len; t++)
alpar@9 132 { j = ind[t];
alpar@9 133 if (val[t] > 0.0)
alpar@9 134 { lb = get_col_lb(lp, j);
alpar@9 135 if (lb == -DBL_MAX)
alpar@9 136 { sum = -DBL_MAX;
alpar@9 137 break;
alpar@9 138 }
alpar@9 139 sum += val[t] * lb;
alpar@9 140 }
alpar@9 141 else if (val[t] < 0.0)
alpar@9 142 { ub = get_col_ub(lp, j);
alpar@9 143 if (ub == +DBL_MAX)
alpar@9 144 { sum = -DBL_MAX;
alpar@9 145 break;
alpar@9 146 }
alpar@9 147 sum += val[t] * ub;
alpar@9 148 }
alpar@9 149 else
alpar@9 150 xassert(val != val);
alpar@9 151 }
alpar@9 152 return sum;
alpar@9 153 }
alpar@9 154
alpar@9 155 static double eval_lf_max(LPX *lp, int len, int ind[], double val[])
alpar@9 156 { /* this routine computes the maximum of a specified linear form
alpar@9 157
alpar@9 158 sum a[j]*x[j]
alpar@9 159 j
alpar@9 160
alpar@9 161 using the formula:
alpar@9 162
alpar@9 163 max = sum a[j]*ub[j] + sum a[j]*lb[j],
alpar@9 164 j in J+ j in J-
alpar@9 165
alpar@9 166 where J+ = {j: a[j] > 0}, J- = {j: a[j] < 0}, lb[j] and ub[j]
alpar@9 167 are lower and upper bound of variable x[j], resp. */
alpar@9 168 int j, t;
alpar@9 169 double lb, ub, sum;
alpar@9 170 sum = 0.0;
alpar@9 171 for (t = 1; t <= len; t++)
alpar@9 172 { j = ind[t];
alpar@9 173 if (val[t] > 0.0)
alpar@9 174 { ub = get_col_ub(lp, j);
alpar@9 175 if (ub == +DBL_MAX)
alpar@9 176 { sum = +DBL_MAX;
alpar@9 177 break;
alpar@9 178 }
alpar@9 179 sum += val[t] * ub;
alpar@9 180 }
alpar@9 181 else if (val[t] < 0.0)
alpar@9 182 { lb = get_col_lb(lp, j);
alpar@9 183 if (lb == -DBL_MAX)
alpar@9 184 { sum = +DBL_MAX;
alpar@9 185 break;
alpar@9 186 }
alpar@9 187 sum += val[t] * lb;
alpar@9 188 }
alpar@9 189 else
alpar@9 190 xassert(val != val);
alpar@9 191 }
alpar@9 192 return sum;
alpar@9 193 }
alpar@9 194
alpar@9 195 /*----------------------------------------------------------------------
alpar@9 196 -- probing - determine logical relation between binary variables.
alpar@9 197 --
alpar@9 198 -- This routine tentatively sets a binary variable to 0 and then to 1
alpar@9 199 -- and examines whether another binary variable is caused to be fixed.
alpar@9 200 --
alpar@9 201 -- The examination is based only on one row (constraint), which is the
alpar@9 202 -- following:
alpar@9 203 --
alpar@9 204 -- L <= sum a[j]*x[j] <= U. (1)
alpar@9 205 -- j
alpar@9 206 --
alpar@9 207 -- Let x[p] be a probing variable, x[q] be an examined variable. Then
alpar@9 208 -- (1) can be written as:
alpar@9 209 --
alpar@9 210 -- L <= sum a[j]*x[j] + a[p]*x[p] + a[q]*x[q] <= U, (2)
alpar@9 211 -- j in J'
alpar@9 212 --
alpar@9 213 -- where J' = {j: j != p and j != q}.
alpar@9 214 --
alpar@9 215 -- Let
alpar@9 216 --
alpar@9 217 -- L' = L - a[p]*x[p], (3)
alpar@9 218 --
alpar@9 219 -- U' = U - a[p]*x[p], (4)
alpar@9 220 --
alpar@9 221 -- where x[p] is assumed to be fixed at 0 or 1. So (2) can be rewritten
alpar@9 222 -- as follows:
alpar@9 223 --
alpar@9 224 -- L' <= sum a[j]*x[j] + a[q]*x[q] <= U', (5)
alpar@9 225 -- j in J'
alpar@9 226 --
alpar@9 227 -- from where we have:
alpar@9 228 --
alpar@9 229 -- L' - sum a[j]*x[j] <= a[q]*x[q] <= U' - sum a[j]*x[j]. (6)
alpar@9 230 -- j in J' j in J'
alpar@9 231 --
alpar@9 232 -- Thus,
alpar@9 233 --
alpar@9 234 -- min a[q]*x[q] = L' - MAX, (7)
alpar@9 235 --
alpar@9 236 -- max a[q]*x[q] = U' - MIN, (8)
alpar@9 237 --
alpar@9 238 -- where
alpar@9 239 --
alpar@9 240 -- MIN = min sum a[j]*x[j], (9)
alpar@9 241 -- j in J'
alpar@9 242 --
alpar@9 243 -- MAX = max sum a[j]*x[j]. (10)
alpar@9 244 -- j in J'
alpar@9 245 --
alpar@9 246 -- Formulae (7) and (8) allows determining implied lower and upper
alpar@9 247 -- bounds of x[q].
alpar@9 248 --
alpar@9 249 -- Parameters len, val, L and U specify the constraint (1).
alpar@9 250 --
alpar@9 251 -- Parameters lf_min and lf_max specify implied lower and upper bounds
alpar@9 252 -- of the linear form (1). It is assumed that these bounds are computed
alpar@9 253 -- with the routines eval_lf_min and eval_lf_max (see above).
alpar@9 254 --
alpar@9 255 -- Parameter p specifies the probing variable x[p], which is set to 0
alpar@9 256 -- (if set is 0) or to 1 (if set is 1).
alpar@9 257 --
alpar@9 258 -- Parameter q specifies the examined variable x[q].
alpar@9 259 --
alpar@9 260 -- On exit the routine returns one of the following codes:
alpar@9 261 --
alpar@9 262 -- 0 - there is no logical relation between x[p] and x[q];
alpar@9 263 -- 1 - x[q] can take only on value 0;
alpar@9 264 -- 2 - x[q] can take only on value 1. */
alpar@9 265
alpar@9 266 static int probing(int len, double val[], double L, double U,
alpar@9 267 double lf_min, double lf_max, int p, int set, int q)
alpar@9 268 { double temp;
alpar@9 269 xassert(1 <= p && p < q && q <= len);
alpar@9 270 /* compute L' (3) */
alpar@9 271 if (L != -DBL_MAX && set) L -= val[p];
alpar@9 272 /* compute U' (4) */
alpar@9 273 if (U != +DBL_MAX && set) U -= val[p];
alpar@9 274 /* compute MIN (9) */
alpar@9 275 if (lf_min != -DBL_MAX)
alpar@9 276 { if (val[p] < 0.0) lf_min -= val[p];
alpar@9 277 if (val[q] < 0.0) lf_min -= val[q];
alpar@9 278 }
alpar@9 279 /* compute MAX (10) */
alpar@9 280 if (lf_max != +DBL_MAX)
alpar@9 281 { if (val[p] > 0.0) lf_max -= val[p];
alpar@9 282 if (val[q] > 0.0) lf_max -= val[q];
alpar@9 283 }
alpar@9 284 /* compute implied lower bound of x[q]; see (7), (8) */
alpar@9 285 if (val[q] > 0.0)
alpar@9 286 { if (L == -DBL_MAX || lf_max == +DBL_MAX)
alpar@9 287 temp = -DBL_MAX;
alpar@9 288 else
alpar@9 289 temp = (L - lf_max) / val[q];
alpar@9 290 }
alpar@9 291 else
alpar@9 292 { if (U == +DBL_MAX || lf_min == -DBL_MAX)
alpar@9 293 temp = -DBL_MAX;
alpar@9 294 else
alpar@9 295 temp = (U - lf_min) / val[q];
alpar@9 296 }
alpar@9 297 if (temp > 0.001) return 2;
alpar@9 298 /* compute implied upper bound of x[q]; see (7), (8) */
alpar@9 299 if (val[q] > 0.0)
alpar@9 300 { if (U == +DBL_MAX || lf_min == -DBL_MAX)
alpar@9 301 temp = +DBL_MAX;
alpar@9 302 else
alpar@9 303 temp = (U - lf_min) / val[q];
alpar@9 304 }
alpar@9 305 else
alpar@9 306 { if (L == -DBL_MAX || lf_max == +DBL_MAX)
alpar@9 307 temp = +DBL_MAX;
alpar@9 308 else
alpar@9 309 temp = (L - lf_max) / val[q];
alpar@9 310 }
alpar@9 311 if (temp < 0.999) return 1;
alpar@9 312 /* there is no logical relation between x[p] and x[q] */
alpar@9 313 return 0;
alpar@9 314 }
alpar@9 315
alpar@9 316 struct COG
alpar@9 317 { /* conflict graph; it represents logical relations between binary
alpar@9 318 variables and has a vertex for each binary variable and its
alpar@9 319 complement, and an edge between two vertices when at most one
alpar@9 320 of the variables represented by the vertices can equal one in
alpar@9 321 an optimal solution */
alpar@9 322 int n;
alpar@9 323 /* number of variables */
alpar@9 324 int nb;
alpar@9 325 /* number of binary variables represented in the graph (note that
alpar@9 326 not all binary variables can be represented); vertices which
alpar@9 327 correspond to binary variables have numbers 1, ..., nb while
alpar@9 328 vertices which correspond to complements of binary variables
alpar@9 329 have numbers nb+1, ..., nb+nb */
alpar@9 330 int ne;
alpar@9 331 /* number of edges in the graph */
alpar@9 332 int *vert; /* int vert[1+n]; */
alpar@9 333 /* if x[j] is a binary variable represented in the graph, vert[j]
alpar@9 334 is the vertex number corresponding to x[j]; otherwise vert[j]
alpar@9 335 is zero */
alpar@9 336 int *orig; /* int list[1:nb]; */
alpar@9 337 /* if vert[j] = k > 0, then orig[k] = j */
alpar@9 338 unsigned char *a;
alpar@9 339 /* adjacency matrix of the graph having 2*nb rows and columns;
alpar@9 340 only strict lower triangle is stored in dense packed form */
alpar@9 341 };
alpar@9 342
alpar@9 343 /*----------------------------------------------------------------------
alpar@9 344 -- lpx_create_cog - create the conflict graph.
alpar@9 345 --
alpar@9 346 -- SYNOPSIS
alpar@9 347 --
alpar@9 348 -- #include "glplpx.h"
alpar@9 349 -- void *lpx_create_cog(LPX *lp);
alpar@9 350 --
alpar@9 351 -- DESCRIPTION
alpar@9 352 --
alpar@9 353 -- The routine lpx_create_cog creates the conflict graph for a given
alpar@9 354 -- problem instance.
alpar@9 355 --
alpar@9 356 -- RETURNS
alpar@9 357 --
alpar@9 358 -- If the graph has been created, the routine returns a pointer to it.
alpar@9 359 -- Otherwise the routine returns NULL. */
alpar@9 360
alpar@9 361 #define MAX_NB 4000
alpar@9 362 #define MAX_ROW_LEN 500
alpar@9 363
alpar@9 364 static void lpx_add_cog_edge(void *_cog, int i, int j);
alpar@9 365
alpar@9 366 static void *lpx_create_cog(LPX *lp)
alpar@9 367 { struct COG *cog = NULL;
alpar@9 368 int m, n, nb, i, j, p, q, len, *ind, *vert, *orig;
alpar@9 369 double L, U, lf_min, lf_max, *val;
alpar@9 370 xprintf("Creating the conflict graph...\n");
alpar@9 371 m = lpx_get_num_rows(lp);
alpar@9 372 n = lpx_get_num_cols(lp);
alpar@9 373 /* determine which binary variables should be included in the
alpar@9 374 conflict graph */
alpar@9 375 nb = 0;
alpar@9 376 vert = xcalloc(1+n, sizeof(int));
alpar@9 377 for (j = 1; j <= n; j++) vert[j] = 0;
alpar@9 378 orig = xcalloc(1+n, sizeof(int));
alpar@9 379 ind = xcalloc(1+n, sizeof(int));
alpar@9 380 val = xcalloc(1+n, sizeof(double));
alpar@9 381 for (i = 1; i <= m; i++)
alpar@9 382 { L = get_row_lb(lp, i);
alpar@9 383 U = get_row_ub(lp, i);
alpar@9 384 if (L == -DBL_MAX && U == +DBL_MAX) continue;
alpar@9 385 len = lpx_get_mat_row(lp, i, ind, val);
alpar@9 386 if (len > MAX_ROW_LEN) continue;
alpar@9 387 lf_min = eval_lf_min(lp, len, ind, val);
alpar@9 388 lf_max = eval_lf_max(lp, len, ind, val);
alpar@9 389 for (p = 1; p <= len; p++)
alpar@9 390 { if (!is_binary(lp, ind[p])) continue;
alpar@9 391 for (q = p+1; q <= len; q++)
alpar@9 392 { if (!is_binary(lp, ind[q])) continue;
alpar@9 393 if (probing(len, val, L, U, lf_min, lf_max, p, 0, q) ||
alpar@9 394 probing(len, val, L, U, lf_min, lf_max, p, 1, q))
alpar@9 395 { /* there is a logical relation */
alpar@9 396 /* include the first variable in the graph */
alpar@9 397 j = ind[p];
alpar@9 398 if (vert[j] == 0) nb++, vert[j] = nb, orig[nb] = j;
alpar@9 399 /* incude the second variable in the graph */
alpar@9 400 j = ind[q];
alpar@9 401 if (vert[j] == 0) nb++, vert[j] = nb, orig[nb] = j;
alpar@9 402 }
alpar@9 403 }
alpar@9 404 }
alpar@9 405 }
alpar@9 406 /* if the graph is either empty or has too many vertices, do not
alpar@9 407 create it */
alpar@9 408 if (nb == 0 || nb > MAX_NB)
alpar@9 409 { xprintf("The conflict graph is either empty or too big\n");
alpar@9 410 xfree(vert);
alpar@9 411 xfree(orig);
alpar@9 412 goto done;
alpar@9 413 }
alpar@9 414 /* create the conflict graph */
alpar@9 415 cog = xmalloc(sizeof(struct COG));
alpar@9 416 cog->n = n;
alpar@9 417 cog->nb = nb;
alpar@9 418 cog->ne = 0;
alpar@9 419 cog->vert = vert;
alpar@9 420 cog->orig = orig;
alpar@9 421 len = nb + nb; /* number of vertices */
alpar@9 422 len = (len * (len - 1)) / 2; /* number of entries in triangle */
alpar@9 423 len = (len + (CHAR_BIT - 1)) / CHAR_BIT; /* bytes needed */
alpar@9 424 cog->a = xmalloc(len);
alpar@9 425 memset(cog->a, 0, len);
alpar@9 426 for (j = 1; j <= nb; j++)
alpar@9 427 { /* add edge between variable and its complement */
alpar@9 428 lpx_add_cog_edge(cog, +orig[j], -orig[j]);
alpar@9 429 }
alpar@9 430 for (i = 1; i <= m; i++)
alpar@9 431 { L = get_row_lb(lp, i);
alpar@9 432 U = get_row_ub(lp, i);
alpar@9 433 if (L == -DBL_MAX && U == +DBL_MAX) continue;
alpar@9 434 len = lpx_get_mat_row(lp, i, ind, val);
alpar@9 435 if (len > MAX_ROW_LEN) continue;
alpar@9 436 lf_min = eval_lf_min(lp, len, ind, val);
alpar@9 437 lf_max = eval_lf_max(lp, len, ind, val);
alpar@9 438 for (p = 1; p <= len; p++)
alpar@9 439 { if (!is_binary(lp, ind[p])) continue;
alpar@9 440 for (q = p+1; q <= len; q++)
alpar@9 441 { if (!is_binary(lp, ind[q])) continue;
alpar@9 442 /* set x[p] to 0 and examine x[q] */
alpar@9 443 switch (probing(len, val, L, U, lf_min, lf_max, p, 0, q))
alpar@9 444 { case 0:
alpar@9 445 /* no logical relation */
alpar@9 446 break;
alpar@9 447 case 1:
alpar@9 448 /* x[p] = 0 implies x[q] = 0 */
alpar@9 449 lpx_add_cog_edge(cog, -ind[p], +ind[q]);
alpar@9 450 break;
alpar@9 451 case 2:
alpar@9 452 /* x[p] = 0 implies x[q] = 1 */
alpar@9 453 lpx_add_cog_edge(cog, -ind[p], -ind[q]);
alpar@9 454 break;
alpar@9 455 default:
alpar@9 456 xassert(lp != lp);
alpar@9 457 }
alpar@9 458 /* set x[p] to 1 and examine x[q] */
alpar@9 459 switch (probing(len, val, L, U, lf_min, lf_max, p, 1, q))
alpar@9 460 { case 0:
alpar@9 461 /* no logical relation */
alpar@9 462 break;
alpar@9 463 case 1:
alpar@9 464 /* x[p] = 1 implies x[q] = 0 */
alpar@9 465 lpx_add_cog_edge(cog, +ind[p], +ind[q]);
alpar@9 466 break;
alpar@9 467 case 2:
alpar@9 468 /* x[p] = 1 implies x[q] = 1 */
alpar@9 469 lpx_add_cog_edge(cog, +ind[p], -ind[q]);
alpar@9 470 break;
alpar@9 471 default:
alpar@9 472 xassert(lp != lp);
alpar@9 473 }
alpar@9 474 }
alpar@9 475 }
alpar@9 476 }
alpar@9 477 xprintf("The conflict graph has 2*%d vertices and %d edges\n",
alpar@9 478 cog->nb, cog->ne);
alpar@9 479 done: xfree(ind);
alpar@9 480 xfree(val);
alpar@9 481 return cog;
alpar@9 482 }
alpar@9 483
alpar@9 484 /*----------------------------------------------------------------------
alpar@9 485 -- lpx_add_cog_edge - add edge to the conflict graph.
alpar@9 486 --
alpar@9 487 -- SYNOPSIS
alpar@9 488 --
alpar@9 489 -- #include "glplpx.h"
alpar@9 490 -- void lpx_add_cog_edge(void *cog, int i, int j);
alpar@9 491 --
alpar@9 492 -- DESCRIPTION
alpar@9 493 --
alpar@9 494 -- The routine lpx_add_cog_edge adds an edge to the conflict graph.
alpar@9 495 -- The edge connects x[i] (if i > 0) or its complement (if i < 0) and
alpar@9 496 -- x[j] (if j > 0) or its complement (if j < 0), where i and j are
alpar@9 497 -- original ordinal numbers of corresponding variables. */
alpar@9 498
alpar@9 499 static void lpx_add_cog_edge(void *_cog, int i, int j)
alpar@9 500 { struct COG *cog = _cog;
alpar@9 501 int k;
alpar@9 502 xassert(i != j);
alpar@9 503 /* determine indices of corresponding vertices */
alpar@9 504 if (i > 0)
alpar@9 505 { xassert(1 <= i && i <= cog->n);
alpar@9 506 i = cog->vert[i];
alpar@9 507 xassert(i != 0);
alpar@9 508 }
alpar@9 509 else
alpar@9 510 { i = -i;
alpar@9 511 xassert(1 <= i && i <= cog->n);
alpar@9 512 i = cog->vert[i];
alpar@9 513 xassert(i != 0);
alpar@9 514 i += cog->nb;
alpar@9 515 }
alpar@9 516 if (j > 0)
alpar@9 517 { xassert(1 <= j && j <= cog->n);
alpar@9 518 j = cog->vert[j];
alpar@9 519 xassert(j != 0);
alpar@9 520 }
alpar@9 521 else
alpar@9 522 { j = -j;
alpar@9 523 xassert(1 <= j && j <= cog->n);
alpar@9 524 j = cog->vert[j];
alpar@9 525 xassert(j != 0);
alpar@9 526 j += cog->nb;
alpar@9 527 }
alpar@9 528 /* only lower triangle is stored, so we need i > j */
alpar@9 529 if (i < j) k = i, i = j, j = k;
alpar@9 530 k = ((i - 1) * (i - 2)) / 2 + (j - 1);
alpar@9 531 cog->a[k / CHAR_BIT] |=
alpar@9 532 (unsigned char)(1 << ((CHAR_BIT - 1) - k % CHAR_BIT));
alpar@9 533 cog->ne++;
alpar@9 534 return;
alpar@9 535 }
alpar@9 536
alpar@9 537 /*----------------------------------------------------------------------
alpar@9 538 -- MAXIMUM WEIGHT CLIQUE
alpar@9 539 --
alpar@9 540 -- Two subroutines sub() and wclique() below are intended to find a
alpar@9 541 -- maximum weight clique in a given undirected graph. These subroutines
alpar@9 542 -- are slightly modified version of the program WCLIQUE developed by
alpar@9 543 -- Patric Ostergard <http://www.tcs.hut.fi/~pat/wclique.html> and based
alpar@9 544 -- on ideas from the article "P. R. J. Ostergard, A new algorithm for
alpar@9 545 -- the maximum-weight clique problem, submitted for publication", which
alpar@9 546 -- in turn is a generalization of the algorithm for unweighted graphs
alpar@9 547 -- presented in "P. R. J. Ostergard, A fast algorithm for the maximum
alpar@9 548 -- clique problem, submitted for publication".
alpar@9 549 --
alpar@9 550 -- USED WITH PERMISSION OF THE AUTHOR OF THE ORIGINAL CODE. */
alpar@9 551
alpar@9 552 struct dsa
alpar@9 553 { /* dynamic storage area */
alpar@9 554 int n;
alpar@9 555 /* number of vertices */
alpar@9 556 int *wt; /* int wt[0:n-1]; */
alpar@9 557 /* weights */
alpar@9 558 unsigned char *a;
alpar@9 559 /* adjacency matrix (packed lower triangle without main diag.) */
alpar@9 560 int record;
alpar@9 561 /* weight of best clique */
alpar@9 562 int rec_level;
alpar@9 563 /* number of vertices in best clique */
alpar@9 564 int *rec; /* int rec[0:n-1]; */
alpar@9 565 /* best clique so far */
alpar@9 566 int *clique; /* int clique[0:n-1]; */
alpar@9 567 /* table for pruning */
alpar@9 568 int *set; /* int set[0:n-1]; */
alpar@9 569 /* current clique */
alpar@9 570 };
alpar@9 571
alpar@9 572 #define n (dsa->n)
alpar@9 573 #define wt (dsa->wt)
alpar@9 574 #define a (dsa->a)
alpar@9 575 #define record (dsa->record)
alpar@9 576 #define rec_level (dsa->rec_level)
alpar@9 577 #define rec (dsa->rec)
alpar@9 578 #define clique (dsa->clique)
alpar@9 579 #define set (dsa->set)
alpar@9 580
alpar@9 581 #if 0
alpar@9 582 static int is_edge(struct dsa *dsa, int i, int j)
alpar@9 583 { /* if there is arc (i,j), the routine returns true; otherwise
alpar@9 584 false; 0 <= i, j < n */
alpar@9 585 int k;
alpar@9 586 xassert(0 <= i && i < n);
alpar@9 587 xassert(0 <= j && j < n);
alpar@9 588 if (i == j) return 0;
alpar@9 589 if (i < j) k = i, i = j, j = k;
alpar@9 590 k = (i * (i - 1)) / 2 + j;
alpar@9 591 return a[k / CHAR_BIT] &
alpar@9 592 (unsigned char)(1 << ((CHAR_BIT - 1) - k % CHAR_BIT));
alpar@9 593 }
alpar@9 594 #else
alpar@9 595 #define is_edge(dsa, i, j) ((i) == (j) ? 0 : \
alpar@9 596 (i) > (j) ? is_edge1(i, j) : is_edge1(j, i))
alpar@9 597 #define is_edge1(i, j) is_edge2(((i) * ((i) - 1)) / 2 + (j))
alpar@9 598 #define is_edge2(k) (a[(k) / CHAR_BIT] & \
alpar@9 599 (unsigned char)(1 << ((CHAR_BIT - 1) - (k) % CHAR_BIT)))
alpar@9 600 #endif
alpar@9 601
alpar@9 602 static void sub(struct dsa *dsa, int ct, int table[], int level,
alpar@9 603 int weight, int l_weight)
alpar@9 604 { int i, j, k, curr_weight, left_weight, *p1, *p2, *newtable;
alpar@9 605 newtable = xcalloc(n, sizeof(int));
alpar@9 606 if (ct <= 0)
alpar@9 607 { /* 0 or 1 elements left; include these */
alpar@9 608 if (ct == 0)
alpar@9 609 { set[level++] = table[0];
alpar@9 610 weight += l_weight;
alpar@9 611 }
alpar@9 612 if (weight > record)
alpar@9 613 { record = weight;
alpar@9 614 rec_level = level;
alpar@9 615 for (i = 0; i < level; i++) rec[i] = set[i];
alpar@9 616 }
alpar@9 617 goto done;
alpar@9 618 }
alpar@9 619 for (i = ct; i >= 0; i--)
alpar@9 620 { if ((level == 0) && (i < ct)) goto done;
alpar@9 621 k = table[i];
alpar@9 622 if ((level > 0) && (clique[k] <= (record - weight)))
alpar@9 623 goto done; /* prune */
alpar@9 624 set[level] = k;
alpar@9 625 curr_weight = weight + wt[k];
alpar@9 626 l_weight -= wt[k];
alpar@9 627 if (l_weight <= (record - curr_weight))
alpar@9 628 goto done; /* prune */
alpar@9 629 p1 = newtable;
alpar@9 630 p2 = table;
alpar@9 631 left_weight = 0;
alpar@9 632 while (p2 < table + i)
alpar@9 633 { j = *p2++;
alpar@9 634 if (is_edge(dsa, j, k))
alpar@9 635 { *p1++ = j;
alpar@9 636 left_weight += wt[j];
alpar@9 637 }
alpar@9 638 }
alpar@9 639 if (left_weight <= (record - curr_weight)) continue;
alpar@9 640 sub(dsa, p1 - newtable - 1, newtable, level + 1, curr_weight,
alpar@9 641 left_weight);
alpar@9 642 }
alpar@9 643 done: xfree(newtable);
alpar@9 644 return;
alpar@9 645 }
alpar@9 646
alpar@9 647 static int wclique(int _n, int w[], unsigned char _a[], int sol[])
alpar@9 648 { struct dsa _dsa, *dsa = &_dsa;
alpar@9 649 int i, j, p, max_wt, max_nwt, wth, *used, *nwt, *pos;
alpar@9 650 glp_long timer;
alpar@9 651 n = _n;
alpar@9 652 wt = &w[1];
alpar@9 653 a = _a;
alpar@9 654 record = 0;
alpar@9 655 rec_level = 0;
alpar@9 656 rec = &sol[1];
alpar@9 657 clique = xcalloc(n, sizeof(int));
alpar@9 658 set = xcalloc(n, sizeof(int));
alpar@9 659 used = xcalloc(n, sizeof(int));
alpar@9 660 nwt = xcalloc(n, sizeof(int));
alpar@9 661 pos = xcalloc(n, sizeof(int));
alpar@9 662 /* start timer */
alpar@9 663 timer = xtime();
alpar@9 664 /* order vertices */
alpar@9 665 for (i = 0; i < n; i++)
alpar@9 666 { nwt[i] = 0;
alpar@9 667 for (j = 0; j < n; j++)
alpar@9 668 if (is_edge(dsa, i, j)) nwt[i] += wt[j];
alpar@9 669 }
alpar@9 670 for (i = 0; i < n; i++)
alpar@9 671 used[i] = 0;
alpar@9 672 for (i = n-1; i >= 0; i--)
alpar@9 673 { max_wt = -1;
alpar@9 674 max_nwt = -1;
alpar@9 675 for (j = 0; j < n; j++)
alpar@9 676 { if ((!used[j]) && ((wt[j] > max_wt) || (wt[j] == max_wt
alpar@9 677 && nwt[j] > max_nwt)))
alpar@9 678 { max_wt = wt[j];
alpar@9 679 max_nwt = nwt[j];
alpar@9 680 p = j;
alpar@9 681 }
alpar@9 682 }
alpar@9 683 pos[i] = p;
alpar@9 684 used[p] = 1;
alpar@9 685 for (j = 0; j < n; j++)
alpar@9 686 if ((!used[j]) && (j != p) && (is_edge(dsa, p, j)))
alpar@9 687 nwt[j] -= wt[p];
alpar@9 688 }
alpar@9 689 /* main routine */
alpar@9 690 wth = 0;
alpar@9 691 for (i = 0; i < n; i++)
alpar@9 692 { wth += wt[pos[i]];
alpar@9 693 sub(dsa, i, pos, 0, 0, wth);
alpar@9 694 clique[pos[i]] = record;
alpar@9 695 #if 0
alpar@9 696 if (utime() >= timer + 5.0)
alpar@9 697 #else
alpar@9 698 if (xdifftime(xtime(), timer) >= 5.0 - 0.001)
alpar@9 699 #endif
alpar@9 700 { /* print current record and reset timer */
alpar@9 701 xprintf("level = %d (%d); best = %d\n", i+1, n, record);
alpar@9 702 #if 0
alpar@9 703 timer = utime();
alpar@9 704 #else
alpar@9 705 timer = xtime();
alpar@9 706 #endif
alpar@9 707 }
alpar@9 708 }
alpar@9 709 xfree(clique);
alpar@9 710 xfree(set);
alpar@9 711 xfree(used);
alpar@9 712 xfree(nwt);
alpar@9 713 xfree(pos);
alpar@9 714 /* return the solution found */
alpar@9 715 for (i = 1; i <= rec_level; i++) sol[i]++;
alpar@9 716 return rec_level;
alpar@9 717 }
alpar@9 718
alpar@9 719 #undef n
alpar@9 720 #undef wt
alpar@9 721 #undef a
alpar@9 722 #undef record
alpar@9 723 #undef rec_level
alpar@9 724 #undef rec
alpar@9 725 #undef clique
alpar@9 726 #undef set
alpar@9 727
alpar@9 728 /*----------------------------------------------------------------------
alpar@9 729 -- lpx_clique_cut - generate cluque cut.
alpar@9 730 --
alpar@9 731 -- SYNOPSIS
alpar@9 732 --
alpar@9 733 -- #include "glplpx.h"
alpar@9 734 -- int lpx_clique_cut(LPX *lp, void *cog, int ind[], double val[]);
alpar@9 735 --
alpar@9 736 -- DESCRIPTION
alpar@9 737 --
alpar@9 738 -- The routine lpx_clique_cut generates a clique cut using the conflict
alpar@9 739 -- graph specified by the parameter cog.
alpar@9 740 --
alpar@9 741 -- If a violated clique cut has been found, it has the following form:
alpar@9 742 --
alpar@9 743 -- sum{j in J} a[j]*x[j] <= b.
alpar@9 744 --
alpar@9 745 -- Variable indices j in J are stored in elements ind[1], ..., ind[len]
alpar@9 746 -- while corresponding constraint coefficients are stored in elements
alpar@9 747 -- val[1], ..., val[len], where len is returned on exit. The right-hand
alpar@9 748 -- side b is stored in element val[0].
alpar@9 749 --
alpar@9 750 -- RETURNS
alpar@9 751 --
alpar@9 752 -- If the cutting plane has been successfully generated, the routine
alpar@9 753 -- returns 1 <= len <= n, which is the number of non-zero coefficients
alpar@9 754 -- in the inequality constraint. Otherwise, the routine returns zero. */
alpar@9 755
alpar@9 756 static int lpx_clique_cut(LPX *lp, void *_cog, int ind[], double val[])
alpar@9 757 { struct COG *cog = _cog;
alpar@9 758 int n = lpx_get_num_cols(lp);
alpar@9 759 int j, t, v, card, temp, len = 0, *w, *sol;
alpar@9 760 double x, sum, b, *vec;
alpar@9 761 /* allocate working arrays */
alpar@9 762 w = xcalloc(1 + 2 * cog->nb, sizeof(int));
alpar@9 763 sol = xcalloc(1 + 2 * cog->nb, sizeof(int));
alpar@9 764 vec = xcalloc(1+n, sizeof(double));
alpar@9 765 /* assign weights to vertices of the conflict graph */
alpar@9 766 for (t = 1; t <= cog->nb; t++)
alpar@9 767 { j = cog->orig[t];
alpar@9 768 x = lpx_get_col_prim(lp, j);
alpar@9 769 temp = (int)(100.0 * x + 0.5);
alpar@9 770 if (temp < 0) temp = 0;
alpar@9 771 if (temp > 100) temp = 100;
alpar@9 772 w[t] = temp;
alpar@9 773 w[cog->nb + t] = 100 - temp;
alpar@9 774 }
alpar@9 775 /* find a clique of maximum weight */
alpar@9 776 card = wclique(2 * cog->nb, w, cog->a, sol);
alpar@9 777 /* compute the clique weight for unscaled values */
alpar@9 778 sum = 0.0;
alpar@9 779 for ( t = 1; t <= card; t++)
alpar@9 780 { v = sol[t];
alpar@9 781 xassert(1 <= v && v <= 2 * cog->nb);
alpar@9 782 if (v <= cog->nb)
alpar@9 783 { /* vertex v corresponds to binary variable x[j] */
alpar@9 784 j = cog->orig[v];
alpar@9 785 x = lpx_get_col_prim(lp, j);
alpar@9 786 sum += x;
alpar@9 787 }
alpar@9 788 else
alpar@9 789 { /* vertex v corresponds to the complement of x[j] */
alpar@9 790 j = cog->orig[v - cog->nb];
alpar@9 791 x = lpx_get_col_prim(lp, j);
alpar@9 792 sum += 1.0 - x;
alpar@9 793 }
alpar@9 794 }
alpar@9 795 /* if the sum of binary variables and their complements in the
alpar@9 796 clique greater than 1, the clique cut is violated */
alpar@9 797 if (sum >= 1.01)
alpar@9 798 { /* construct the inquality */
alpar@9 799 for (j = 1; j <= n; j++) vec[j] = 0;
alpar@9 800 b = 1.0;
alpar@9 801 for (t = 1; t <= card; t++)
alpar@9 802 { v = sol[t];
alpar@9 803 if (v <= cog->nb)
alpar@9 804 { /* vertex v corresponds to binary variable x[j] */
alpar@9 805 j = cog->orig[v];
alpar@9 806 xassert(1 <= j && j <= n);
alpar@9 807 vec[j] += 1.0;
alpar@9 808 }
alpar@9 809 else
alpar@9 810 { /* vertex v corresponds to the complement of x[j] */
alpar@9 811 j = cog->orig[v - cog->nb];
alpar@9 812 xassert(1 <= j && j <= n);
alpar@9 813 vec[j] -= 1.0;
alpar@9 814 b -= 1.0;
alpar@9 815 }
alpar@9 816 }
alpar@9 817 xassert(len == 0);
alpar@9 818 for (j = 1; j <= n; j++)
alpar@9 819 { if (vec[j] != 0.0)
alpar@9 820 { len++;
alpar@9 821 ind[len] = j, val[len] = vec[j];
alpar@9 822 }
alpar@9 823 }
alpar@9 824 ind[0] = 0, val[0] = b;
alpar@9 825 }
alpar@9 826 /* free working arrays */
alpar@9 827 xfree(w);
alpar@9 828 xfree(sol);
alpar@9 829 xfree(vec);
alpar@9 830 /* return to the calling program */
alpar@9 831 return len;
alpar@9 832 }
alpar@9 833
alpar@9 834 /*----------------------------------------------------------------------
alpar@9 835 -- lpx_delete_cog - delete the conflict graph.
alpar@9 836 --
alpar@9 837 -- SYNOPSIS
alpar@9 838 --
alpar@9 839 -- #include "glplpx.h"
alpar@9 840 -- void lpx_delete_cog(void *cog);
alpar@9 841 --
alpar@9 842 -- DESCRIPTION
alpar@9 843 --
alpar@9 844 -- The routine lpx_delete_cog deletes the conflict graph, which the
alpar@9 845 -- parameter cog points to, freeing all the memory allocated to this
alpar@9 846 -- object. */
alpar@9 847
alpar@9 848 static void lpx_delete_cog(void *_cog)
alpar@9 849 { struct COG *cog = _cog;
alpar@9 850 xfree(cog->vert);
alpar@9 851 xfree(cog->orig);
alpar@9 852 xfree(cog->a);
alpar@9 853 xfree(cog);
alpar@9 854 }
alpar@9 855
alpar@9 856 /**********************************************************************/
alpar@9 857
alpar@9 858 void *ios_clq_init(glp_tree *tree)
alpar@9 859 { /* initialize clique cut generator */
alpar@9 860 glp_prob *mip = tree->mip;
alpar@9 861 xassert(mip != NULL);
alpar@9 862 return lpx_create_cog(mip);
alpar@9 863 }
alpar@9 864
alpar@9 865 /***********************************************************************
alpar@9 866 * NAME
alpar@9 867 *
alpar@9 868 * ios_clq_gen - generate clique cuts
alpar@9 869 *
alpar@9 870 * SYNOPSIS
alpar@9 871 *
alpar@9 872 * #include "glpios.h"
alpar@9 873 * void ios_clq_gen(glp_tree *tree, void *gen);
alpar@9 874 *
alpar@9 875 * DESCRIPTION
alpar@9 876 *
alpar@9 877 * The routine ios_clq_gen generates clique cuts for the current point
alpar@9 878 * and adds them to the clique pool. */
alpar@9 879
alpar@9 880 void ios_clq_gen(glp_tree *tree, void *gen)
alpar@9 881 { int n = lpx_get_num_cols(tree->mip);
alpar@9 882 int len, *ind;
alpar@9 883 double *val;
alpar@9 884 xassert(gen != NULL);
alpar@9 885 ind = xcalloc(1+n, sizeof(int));
alpar@9 886 val = xcalloc(1+n, sizeof(double));
alpar@9 887 len = lpx_clique_cut(tree->mip, gen, ind, val);
alpar@9 888 if (len > 0)
alpar@9 889 { /* xprintf("len = %d\n", len); */
alpar@9 890 glp_ios_add_row(tree, NULL, GLP_RF_CLQ, 0, len, ind, val,
alpar@9 891 GLP_UP, val[0]);
alpar@9 892 }
alpar@9 893 xfree(ind);
alpar@9 894 xfree(val);
alpar@9 895 return;
alpar@9 896 }
alpar@9 897
alpar@9 898 /**********************************************************************/
alpar@9 899
alpar@9 900 void ios_clq_term(void *gen)
alpar@9 901 { /* terminate clique cut generator */
alpar@9 902 xassert(gen != NULL);
alpar@9 903 lpx_delete_cog(gen);
alpar@9 904 return;
alpar@9 905 }
alpar@9 906
alpar@9 907 /* eof */