diff -r d59bea55db9b -r c445c931472f src/glpnet08.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/glpnet08.c	Mon Dec 06 13:09:21 2010 +0100
@@ -0,0 +1,241 @@
+/* glpnet08.c */
+
+/***********************************************************************
+*  This code is part of GLPK (GNU Linear Programming Kit).
+*
+*  Two subroutines sub() and wclique() below are intended to find a
+*  maximum weight clique in a given undirected graph. These subroutines
+*  are slightly modified version of the program WCLIQUE developed by
+*  Patric Ostergard <http://www.tcs.hut.fi/~pat/wclique.html> and based
+*  on ideas from the article "P. R. J. Ostergard, A new algorithm for
+*  the maximum-weight clique problem, submitted for publication", which
+*  in turn is a generalization of the algorithm for unweighted graphs
+*  presented in "P. R. J. Ostergard, A fast algorithm for the maximum
+*  clique problem, submitted for publication".
+*
+*  USED WITH PERMISSION OF THE AUTHOR OF THE ORIGINAL CODE.
+*
+*  Changes were made by Andrew Makhorin <mao@gnu.org>.
+*
+*  GLPK is free software: you can redistribute it and/or modify it
+*  under the terms of the GNU General Public License as published by
+*  the Free Software Foundation, either version 3 of the License, or
+*  (at your option) any later version.
+*
+*  GLPK is distributed in the hope that it will be useful, but WITHOUT
+*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
+*  License for more details.
+*
+*  You should have received a copy of the GNU General Public License
+*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
+***********************************************************************/
+
+#include "glpenv.h"
+#include "glpnet.h"
+
+/***********************************************************************
+*  NAME
+*
+*  wclique - find maximum weight clique with Ostergard's algorithm
+*
+*  SYNOPSIS
+*
+*  int wclique(int n, const int w[], const unsigned char a[],
+*     int ind[]);
+*
+*  DESCRIPTION
+*
+*  The routine wclique finds a maximum weight clique in an undirected
+*  graph with Ostergard's algorithm.
+*
+*  INPUT PARAMETERS
+*
+*  n is the number of vertices, n > 0.
+*
+*  w[i], i = 1,...,n, is a weight of vertex i.
+*
+*  a[*] is the strict (without main diagonal) lower triangle of the
+*  graph adjacency matrix in packed format.
+*
+*  OUTPUT PARAMETER
+*
+*  ind[k], k = 1,...,size, is the number of a vertex included in the
+*  clique found, 1 <= ind[k] <= n, where size is the number of vertices
+*  in the clique returned on exit.
+*
+*  RETURNS
+*
+*  The routine returns the clique size, i.e. the number of vertices in
+*  the clique. */
+
+struct csa
+{     /* common storage area */
+      int n;
+      /* number of vertices */
+      const int *wt; /* int wt[0:n-1]; */
+      /* weights */
+      const unsigned char *a;
+      /* adjacency matrix (packed lower triangle without main diag.) */
+      int record;
+      /* weight of best clique */
+      int rec_level;
+      /* number of vertices in best clique */
+      int *rec; /* int rec[0:n-1]; */
+      /* best clique so far */
+      int *clique; /* int clique[0:n-1]; */
+      /* table for pruning */
+      int *set; /* int set[0:n-1]; */
+      /* current clique */
+};
+
+#define n         (csa->n)
+#define wt        (csa->wt)
+#define a         (csa->a)
+#define record    (csa->record)
+#define rec_level (csa->rec_level)
+#define rec       (csa->rec)
+#define clique    (csa->clique)
+#define set       (csa->set)
+
+#if 0
+static int is_edge(struct csa *csa, int i, int j)
+{     /* if there is arc (i,j), the routine returns true; otherwise
+         false; 0 <= i, j < n */
+      int k;
+      xassert(0 <= i && i < n);
+      xassert(0 <= j && j < n);
+      if (i == j) return 0;
+      if (i < j) k = i, i = j, j = k;
+      k = (i * (i - 1)) / 2 + j;
+      return a[k / CHAR_BIT] &
+         (unsigned char)(1 << ((CHAR_BIT - 1) - k % CHAR_BIT));
+}
+#else
+#define is_edge(csa, i, j) ((i) == (j) ? 0 : \
+      (i) > (j) ? is_edge1(i, j) : is_edge1(j, i))
+#define is_edge1(i, j) is_edge2(((i) * ((i) - 1)) / 2 + (j))
+#define is_edge2(k) (a[(k) / CHAR_BIT] & \
+      (unsigned char)(1 << ((CHAR_BIT - 1) - (k) % CHAR_BIT)))
+#endif
+
+static void sub(struct csa *csa, int ct, int table[], int level,
+      int weight, int l_weight)
+{     int i, j, k, curr_weight, left_weight, *p1, *p2, *newtable;
+      newtable = xcalloc(n, sizeof(int));
+      if (ct <= 0)
+      {  /* 0 or 1 elements left; include these */
+         if (ct == 0)
+         {  set[level++] = table[0];
+            weight += l_weight;
+         }
+         if (weight > record)
+         {  record = weight;
+            rec_level = level;
+            for (i = 0; i < level; i++) rec[i] = set[i];
+         }
+         goto done;
+      }
+      for (i = ct; i >= 0; i--)
+      {  if ((level == 0) && (i < ct)) goto done;
+         k = table[i];
+         if ((level > 0) && (clique[k] <= (record - weight)))
+            goto done; /* prune */
+         set[level] = k;
+         curr_weight = weight + wt[k];
+         l_weight -= wt[k];
+         if (l_weight <= (record - curr_weight))
+            goto done; /* prune */
+         p1 = newtable;
+         p2 = table;
+         left_weight = 0;
+         while (p2 < table + i)
+         {  j = *p2++;
+            if (is_edge(csa, j, k))
+            {  *p1++ = j;
+               left_weight += wt[j];
+            }
+         }
+         if (left_weight <= (record - curr_weight)) continue;
+         sub(csa, p1 - newtable - 1, newtable, level + 1, curr_weight,
+            left_weight);
+      }
+done: xfree(newtable);
+      return;
+}
+
+int wclique(int _n, const int w[], const unsigned char _a[], int ind[])
+{     struct csa _csa, *csa = &_csa;
+      int i, j, p, max_wt, max_nwt, wth, *used, *nwt, *pos;
+      glp_long timer;
+      n = _n;
+      xassert(n > 0);
+      wt = &w[1];
+      a = _a;
+      record = 0;
+      rec_level = 0;
+      rec = &ind[1];
+      clique = xcalloc(n, sizeof(int));
+      set = xcalloc(n, sizeof(int));
+      used = xcalloc(n, sizeof(int));
+      nwt = xcalloc(n, sizeof(int));
+      pos = xcalloc(n, sizeof(int));
+      /* start timer */
+      timer = xtime();
+      /* order vertices */
+      for (i = 0; i < n; i++)
+      {  nwt[i] = 0;
+         for (j = 0; j < n; j++)
+            if (is_edge(csa, i, j)) nwt[i] += wt[j];
+      }
+      for (i = 0; i < n; i++)
+         used[i] = 0;
+      for (i = n-1; i >= 0; i--)
+      {  max_wt = -1;
+         max_nwt = -1;
+         for (j = 0; j < n; j++)
+         {  if ((!used[j]) && ((wt[j] > max_wt) || (wt[j] == max_wt
+               && nwt[j] > max_nwt)))
+            {  max_wt = wt[j];
+               max_nwt = nwt[j];
+               p = j;
+            }
+         }
+         pos[i] = p;
+         used[p] = 1;
+         for (j = 0; j < n; j++)
+            if ((!used[j]) && (j != p) && (is_edge(csa, p, j)))
+               nwt[j] -= wt[p];
+      }
+      /* main routine */
+      wth = 0;
+      for (i = 0; i < n; i++)
+      {  wth += wt[pos[i]];
+         sub(csa, i, pos, 0, 0, wth);
+         clique[pos[i]] = record;
+         if (xdifftime(xtime(), timer) >= 5.0 - 0.001)
+         {  /* print current record and reset timer */
+            xprintf("level = %d (%d); best = %d\n", i+1, n, record);
+            timer = xtime();
+         }
+      }
+      xfree(clique);
+      xfree(set);
+      xfree(used);
+      xfree(nwt);
+      xfree(pos);
+      /* return the solution found */
+      for (i = 1; i <= rec_level; i++) ind[i]++;
+      return rec_level;
+}
+
+#undef n
+#undef wt
+#undef a
+#undef record
+#undef rec_level
+#undef rec
+#undef clique
+#undef set
+
+/* eof */