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