PageRenderTime 61ms CodeModel.GetById 17ms app.highlight 36ms RepoModel.GetById 2ms app.codeStats 0ms

/farmR/src/glpios08.c

https://code.google.com/p/javawfm/
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 */