1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpnet08.c	Sun Nov 06 20:59:10 2011 +0100
     1.3 @@ -0,0 +1,241 @@
     1.4 +/* glpnet08.c */
     1.5 +
     1.6 +/***********************************************************************
     1.7 +*  This code is part of GLPK (GNU Linear Programming Kit).
     1.8 +*
     1.9 +*  Two subroutines sub() and wclique() below are intended to find a
    1.10 +*  maximum weight clique in a given undirected graph. These subroutines
    1.11 +*  are slightly modified version of the program WCLIQUE developed by
    1.12 +*  Patric Ostergard <http://www.tcs.hut.fi/~pat/wclique.html> and based
    1.13 +*  on ideas from the article "P. R. J. Ostergard, A new algorithm for
    1.14 +*  the maximum-weight clique problem, submitted for publication", which
    1.15 +*  in turn is a generalization of the algorithm for unweighted graphs
    1.16 +*  presented in "P. R. J. Ostergard, A fast algorithm for the maximum
    1.17 +*  clique problem, submitted for publication".
    1.18 +*
    1.19 +*  USED WITH PERMISSION OF THE AUTHOR OF THE ORIGINAL CODE.
    1.20 +*
    1.21 +*  Changes were made by Andrew Makhorin <mao@gnu.org>.
    1.22 +*
    1.23 +*  GLPK is free software: you can redistribute it and/or modify it
    1.24 +*  under the terms of the GNU General Public License as published by
    1.25 +*  the Free Software Foundation, either version 3 of the License, or
    1.26 +*  (at your option) any later version.
    1.27 +*
    1.28 +*  GLPK is distributed in the hope that it will be useful, but WITHOUT
    1.29 +*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    1.30 +*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
    1.31 +*  License for more details.
    1.32 +*
    1.33 +*  You should have received a copy of the GNU General Public License
    1.34 +*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
    1.35 +***********************************************************************/
    1.36 +
    1.37 +#include "glpenv.h"
    1.38 +#include "glpnet.h"
    1.39 +
    1.40 +/***********************************************************************
    1.41 +*  NAME
    1.42 +*
    1.43 +*  wclique - find maximum weight clique with Ostergard's algorithm
    1.44 +*
    1.45 +*  SYNOPSIS
    1.46 +*
    1.47 +*  int wclique(int n, const int w[], const unsigned char a[],
    1.48 +*     int ind[]);
    1.49 +*
    1.50 +*  DESCRIPTION
    1.51 +*
    1.52 +*  The routine wclique finds a maximum weight clique in an undirected
    1.53 +*  graph with Ostergard's algorithm.
    1.54 +*
    1.55 +*  INPUT PARAMETERS
    1.56 +*
    1.57 +*  n is the number of vertices, n > 0.
    1.58 +*
    1.59 +*  w[i], i = 1,...,n, is a weight of vertex i.
    1.60 +*
    1.61 +*  a[*] is the strict (without main diagonal) lower triangle of the
    1.62 +*  graph adjacency matrix in packed format.
    1.63 +*
    1.64 +*  OUTPUT PARAMETER
    1.65 +*
    1.66 +*  ind[k], k = 1,...,size, is the number of a vertex included in the
    1.67 +*  clique found, 1 <= ind[k] <= n, where size is the number of vertices
    1.68 +*  in the clique returned on exit.
    1.69 +*
    1.70 +*  RETURNS
    1.71 +*
    1.72 +*  The routine returns the clique size, i.e. the number of vertices in
    1.73 +*  the clique. */
    1.74 +
    1.75 +struct csa
    1.76 +{     /* common storage area */
    1.77 +      int n;
    1.78 +      /* number of vertices */
    1.79 +      const int *wt; /* int wt[0:n-1]; */
    1.80 +      /* weights */
    1.81 +      const unsigned char *a;
    1.82 +      /* adjacency matrix (packed lower triangle without main diag.) */
    1.83 +      int record;
    1.84 +      /* weight of best clique */
    1.85 +      int rec_level;
    1.86 +      /* number of vertices in best clique */
    1.87 +      int *rec; /* int rec[0:n-1]; */
    1.88 +      /* best clique so far */
    1.89 +      int *clique; /* int clique[0:n-1]; */
    1.90 +      /* table for pruning */
    1.91 +      int *set; /* int set[0:n-1]; */
    1.92 +      /* current clique */
    1.93 +};
    1.94 +
    1.95 +#define n         (csa->n)
    1.96 +#define wt        (csa->wt)
    1.97 +#define a         (csa->a)
    1.98 +#define record    (csa->record)
    1.99 +#define rec_level (csa->rec_level)
   1.100 +#define rec       (csa->rec)
   1.101 +#define clique    (csa->clique)
   1.102 +#define set       (csa->set)
   1.103 +
   1.104 +#if 0
   1.105 +static int is_edge(struct csa *csa, int i, int j)
   1.106 +{     /* if there is arc (i,j), the routine returns true; otherwise
   1.107 +         false; 0 <= i, j < n */
   1.108 +      int k;
   1.109 +      xassert(0 <= i && i < n);
   1.110 +      xassert(0 <= j && j < n);
   1.111 +      if (i == j) return 0;
   1.112 +      if (i < j) k = i, i = j, j = k;
   1.113 +      k = (i * (i - 1)) / 2 + j;
   1.114 +      return a[k / CHAR_BIT] &
   1.115 +         (unsigned char)(1 << ((CHAR_BIT - 1) - k % CHAR_BIT));
   1.116 +}
   1.117 +#else
   1.118 +#define is_edge(csa, i, j) ((i) == (j) ? 0 : \
   1.119 +      (i) > (j) ? is_edge1(i, j) : is_edge1(j, i))
   1.120 +#define is_edge1(i, j) is_edge2(((i) * ((i) - 1)) / 2 + (j))
   1.121 +#define is_edge2(k) (a[(k) / CHAR_BIT] & \
   1.122 +      (unsigned char)(1 << ((CHAR_BIT - 1) - (k) % CHAR_BIT)))
   1.123 +#endif
   1.124 +
   1.125 +static void sub(struct csa *csa, int ct, int table[], int level,
   1.126 +      int weight, int l_weight)
   1.127 +{     int i, j, k, curr_weight, left_weight, *p1, *p2, *newtable;
   1.128 +      newtable = xcalloc(n, sizeof(int));
   1.129 +      if (ct <= 0)
   1.130 +      {  /* 0 or 1 elements left; include these */
   1.131 +         if (ct == 0)
   1.132 +         {  set[level++] = table[0];
   1.133 +            weight += l_weight;
   1.134 +         }
   1.135 +         if (weight > record)
   1.136 +         {  record = weight;
   1.137 +            rec_level = level;
   1.138 +            for (i = 0; i < level; i++) rec[i] = set[i];
   1.139 +         }
   1.140 +         goto done;
   1.141 +      }
   1.142 +      for (i = ct; i >= 0; i--)
   1.143 +      {  if ((level == 0) && (i < ct)) goto done;
   1.144 +         k = table[i];
   1.145 +         if ((level > 0) && (clique[k] <= (record - weight)))
   1.146 +            goto done; /* prune */
   1.147 +         set[level] = k;
   1.148 +         curr_weight = weight + wt[k];
   1.149 +         l_weight -= wt[k];
   1.150 +         if (l_weight <= (record - curr_weight))
   1.151 +            goto done; /* prune */
   1.152 +         p1 = newtable;
   1.153 +         p2 = table;
   1.154 +         left_weight = 0;
   1.155 +         while (p2 < table + i)
   1.156 +         {  j = *p2++;
   1.157 +            if (is_edge(csa, j, k))
   1.158 +            {  *p1++ = j;
   1.159 +               left_weight += wt[j];
   1.160 +            }
   1.161 +         }
   1.162 +         if (left_weight <= (record - curr_weight)) continue;
   1.163 +         sub(csa, p1 - newtable - 1, newtable, level + 1, curr_weight,
   1.164 +            left_weight);
   1.165 +      }
   1.166 +done: xfree(newtable);
   1.167 +      return;
   1.168 +}
   1.169 +
   1.170 +int wclique(int _n, const int w[], const unsigned char _a[], int ind[])
   1.171 +{     struct csa _csa, *csa = &_csa;
   1.172 +      int i, j, p, max_wt, max_nwt, wth, *used, *nwt, *pos;
   1.173 +      glp_long timer;
   1.174 +      n = _n;
   1.175 +      xassert(n > 0);
   1.176 +      wt = &w[1];
   1.177 +      a = _a;
   1.178 +      record = 0;
   1.179 +      rec_level = 0;
   1.180 +      rec = &ind[1];
   1.181 +      clique = xcalloc(n, sizeof(int));
   1.182 +      set = xcalloc(n, sizeof(int));
   1.183 +      used = xcalloc(n, sizeof(int));
   1.184 +      nwt = xcalloc(n, sizeof(int));
   1.185 +      pos = xcalloc(n, sizeof(int));
   1.186 +      /* start timer */
   1.187 +      timer = xtime();
   1.188 +      /* order vertices */
   1.189 +      for (i = 0; i < n; i++)
   1.190 +      {  nwt[i] = 0;
   1.191 +         for (j = 0; j < n; j++)
   1.192 +            if (is_edge(csa, i, j)) nwt[i] += wt[j];
   1.193 +      }
   1.194 +      for (i = 0; i < n; i++)
   1.195 +         used[i] = 0;
   1.196 +      for (i = n-1; i >= 0; i--)
   1.197 +      {  max_wt = -1;
   1.198 +         max_nwt = -1;
   1.199 +         for (j = 0; j < n; j++)
   1.200 +         {  if ((!used[j]) && ((wt[j] > max_wt) || (wt[j] == max_wt
   1.201 +               && nwt[j] > max_nwt)))
   1.202 +            {  max_wt = wt[j];
   1.203 +               max_nwt = nwt[j];
   1.204 +               p = j;
   1.205 +            }
   1.206 +         }
   1.207 +         pos[i] = p;
   1.208 +         used[p] = 1;
   1.209 +         for (j = 0; j < n; j++)
   1.210 +            if ((!used[j]) && (j != p) && (is_edge(csa, p, j)))
   1.211 +               nwt[j] -= wt[p];
   1.212 +      }
   1.213 +      /* main routine */
   1.214 +      wth = 0;
   1.215 +      for (i = 0; i < n; i++)
   1.216 +      {  wth += wt[pos[i]];
   1.217 +         sub(csa, i, pos, 0, 0, wth);
   1.218 +         clique[pos[i]] = record;
   1.219 +         if (xdifftime(xtime(), timer) >= 5.0 - 0.001)
   1.220 +         {  /* print current record and reset timer */
   1.221 +            xprintf("level = %d (%d); best = %d\n", i+1, n, record);
   1.222 +            timer = xtime();
   1.223 +         }
   1.224 +      }
   1.225 +      xfree(clique);
   1.226 +      xfree(set);
   1.227 +      xfree(used);
   1.228 +      xfree(nwt);
   1.229 +      xfree(pos);
   1.230 +      /* return the solution found */
   1.231 +      for (i = 1; i <= rec_level; i++) ind[i]++;
   1.232 +      return rec_level;
   1.233 +}
   1.234 +
   1.235 +#undef n
   1.236 +#undef wt
   1.237 +#undef a
   1.238 +#undef record
   1.239 +#undef rec_level
   1.240 +#undef rec
   1.241 +#undef clique
   1.242 +#undef set
   1.243 +
   1.244 +/* eof */