COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glpios08.c @ 1:c445c931472f

Last change on this file since 1:c445c931472f was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 14 years ago

Import glpk-4.45

  • Generated files and doc/notes are removed
File size: 27.6 KB
Line 
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 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
27static 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
47static 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
67static 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
87static 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
107static 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
115static 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
155static 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
266static 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
316struct 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
364static void lpx_add_cog_edge(void *_cog, int i, int j);
365
366static 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);
479done: 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
499static 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
552struct 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
582static 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
602static 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      }
643done: xfree(newtable);
644      return;
645}
646
647static 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
756static 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
848static 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
858void *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
880void 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
900void 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 */
Note: See TracBrowser for help on using the repository browser.