lemon-project-template-glpk

comparison deps/glpk/src/glpios08.c @ 9:33de93886c88

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