COIN-OR::LEMON - Graph Library

source: lemon-project-template-glpk/deps/glpk/src/glpnet08.c

subpack-glpk
Last change on this file was 9:33de93886c88, checked in by Alpar Juttner <alpar@…>, 13 years ago

Import GLPK 4.47

File size: 7.3 KB
Line 
1/* glpnet08.c */
2
3/***********************************************************************
4*  This code is part of GLPK (GNU Linear Programming Kit).
5*
6*  Two subroutines sub() and wclique() below are intended to find a
7*  maximum weight clique in a given undirected graph. These subroutines
8*  are slightly modified version of the program WCLIQUE developed by
9*  Patric Ostergard <http://www.tcs.hut.fi/~pat/wclique.html> and based
10*  on ideas from the article "P. R. J. Ostergard, A new algorithm for
11*  the maximum-weight clique problem, submitted for publication", which
12*  in turn is a generalization of the algorithm for unweighted graphs
13*  presented in "P. R. J. Ostergard, A fast algorithm for the maximum
14*  clique problem, submitted for publication".
15*
16*  USED WITH PERMISSION OF THE AUTHOR OF THE ORIGINAL CODE.
17*
18*  Changes were made by Andrew Makhorin <mao@gnu.org>.
19*
20*  GLPK is free software: you can redistribute it and/or modify it
21*  under the terms of the GNU General Public License as published by
22*  the Free Software Foundation, either version 3 of the License, or
23*  (at your option) any later version.
24*
25*  GLPK is distributed in the hope that it will be useful, but WITHOUT
26*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
27*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
28*  License for more details.
29*
30*  You should have received a copy of the GNU General Public License
31*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
32***********************************************************************/
33
34#include "glpenv.h"
35#include "glpnet.h"
36
37/***********************************************************************
38*  NAME
39*
40*  wclique - find maximum weight clique with Ostergard's algorithm
41*
42*  SYNOPSIS
43*
44*  int wclique(int n, const int w[], const unsigned char a[],
45*     int ind[]);
46*
47*  DESCRIPTION
48*
49*  The routine wclique finds a maximum weight clique in an undirected
50*  graph with Ostergard's algorithm.
51*
52*  INPUT PARAMETERS
53*
54*  n is the number of vertices, n > 0.
55*
56*  w[i], i = 1,...,n, is a weight of vertex i.
57*
58*  a[*] is the strict (without main diagonal) lower triangle of the
59*  graph adjacency matrix in packed format.
60*
61*  OUTPUT PARAMETER
62*
63*  ind[k], k = 1,...,size, is the number of a vertex included in the
64*  clique found, 1 <= ind[k] <= n, where size is the number of vertices
65*  in the clique returned on exit.
66*
67*  RETURNS
68*
69*  The routine returns the clique size, i.e. the number of vertices in
70*  the clique. */
71
72struct csa
73{     /* common storage area */
74      int n;
75      /* number of vertices */
76      const int *wt; /* int wt[0:n-1]; */
77      /* weights */
78      const unsigned char *a;
79      /* adjacency matrix (packed lower triangle without main diag.) */
80      int record;
81      /* weight of best clique */
82      int rec_level;
83      /* number of vertices in best clique */
84      int *rec; /* int rec[0:n-1]; */
85      /* best clique so far */
86      int *clique; /* int clique[0:n-1]; */
87      /* table for pruning */
88      int *set; /* int set[0:n-1]; */
89      /* current clique */
90};
91
92#define n         (csa->n)
93#define wt        (csa->wt)
94#define a         (csa->a)
95#define record    (csa->record)
96#define rec_level (csa->rec_level)
97#define rec       (csa->rec)
98#define clique    (csa->clique)
99#define set       (csa->set)
100
101#if 0
102static int is_edge(struct csa *csa, int i, int j)
103{     /* if there is arc (i,j), the routine returns true; otherwise
104         false; 0 <= i, j < n */
105      int k;
106      xassert(0 <= i && i < n);
107      xassert(0 <= j && j < n);
108      if (i == j) return 0;
109      if (i < j) k = i, i = j, j = k;
110      k = (i * (i - 1)) / 2 + j;
111      return a[k / CHAR_BIT] &
112         (unsigned char)(1 << ((CHAR_BIT - 1) - k % CHAR_BIT));
113}
114#else
115#define is_edge(csa, i, j) ((i) == (j) ? 0 : \
116      (i) > (j) ? is_edge1(i, j) : is_edge1(j, i))
117#define is_edge1(i, j) is_edge2(((i) * ((i) - 1)) / 2 + (j))
118#define is_edge2(k) (a[(k) / CHAR_BIT] & \
119      (unsigned char)(1 << ((CHAR_BIT - 1) - (k) % CHAR_BIT)))
120#endif
121
122static void sub(struct csa *csa, int ct, int table[], int level,
123      int weight, int l_weight)
124{     int i, j, k, curr_weight, left_weight, *p1, *p2, *newtable;
125      newtable = xcalloc(n, sizeof(int));
126      if (ct <= 0)
127      {  /* 0 or 1 elements left; include these */
128         if (ct == 0)
129         {  set[level++] = table[0];
130            weight += l_weight;
131         }
132         if (weight > record)
133         {  record = weight;
134            rec_level = level;
135            for (i = 0; i < level; i++) rec[i] = set[i];
136         }
137         goto done;
138      }
139      for (i = ct; i >= 0; i--)
140      {  if ((level == 0) && (i < ct)) goto done;
141         k = table[i];
142         if ((level > 0) && (clique[k] <= (record - weight)))
143            goto done; /* prune */
144         set[level] = k;
145         curr_weight = weight + wt[k];
146         l_weight -= wt[k];
147         if (l_weight <= (record - curr_weight))
148            goto done; /* prune */
149         p1 = newtable;
150         p2 = table;
151         left_weight = 0;
152         while (p2 < table + i)
153         {  j = *p2++;
154            if (is_edge(csa, j, k))
155            {  *p1++ = j;
156               left_weight += wt[j];
157            }
158         }
159         if (left_weight <= (record - curr_weight)) continue;
160         sub(csa, p1 - newtable - 1, newtable, level + 1, curr_weight,
161            left_weight);
162      }
163done: xfree(newtable);
164      return;
165}
166
167int wclique(int _n, const int w[], const unsigned char _a[], int ind[])
168{     struct csa _csa, *csa = &_csa;
169      int i, j, p, max_wt, max_nwt, wth, *used, *nwt, *pos;
170      glp_long timer;
171      n = _n;
172      xassert(n > 0);
173      wt = &w[1];
174      a = _a;
175      record = 0;
176      rec_level = 0;
177      rec = &ind[1];
178      clique = xcalloc(n, sizeof(int));
179      set = xcalloc(n, sizeof(int));
180      used = xcalloc(n, sizeof(int));
181      nwt = xcalloc(n, sizeof(int));
182      pos = xcalloc(n, sizeof(int));
183      /* start timer */
184      timer = xtime();
185      /* order vertices */
186      for (i = 0; i < n; i++)
187      {  nwt[i] = 0;
188         for (j = 0; j < n; j++)
189            if (is_edge(csa, i, j)) nwt[i] += wt[j];
190      }
191      for (i = 0; i < n; i++)
192         used[i] = 0;
193      for (i = n-1; i >= 0; i--)
194      {  max_wt = -1;
195         max_nwt = -1;
196         for (j = 0; j < n; j++)
197         {  if ((!used[j]) && ((wt[j] > max_wt) || (wt[j] == max_wt
198               && nwt[j] > max_nwt)))
199            {  max_wt = wt[j];
200               max_nwt = nwt[j];
201               p = j;
202            }
203         }
204         pos[i] = p;
205         used[p] = 1;
206         for (j = 0; j < n; j++)
207            if ((!used[j]) && (j != p) && (is_edge(csa, p, j)))
208               nwt[j] -= wt[p];
209      }
210      /* main routine */
211      wth = 0;
212      for (i = 0; i < n; i++)
213      {  wth += wt[pos[i]];
214         sub(csa, i, pos, 0, 0, wth);
215         clique[pos[i]] = record;
216         if (xdifftime(xtime(), timer) >= 5.0 - 0.001)
217         {  /* print current record and reset timer */
218            xprintf("level = %d (%d); best = %d\n", i+1, n, record);
219            timer = xtime();
220         }
221      }
222      xfree(clique);
223      xfree(set);
224      xfree(used);
225      xfree(nwt);
226      xfree(pos);
227      /* return the solution found */
228      for (i = 1; i <= rec_level; i++) ind[i]++;
229      return rec_level;
230}
231
232#undef n
233#undef wt
234#undef a
235#undef record
236#undef rec_level
237#undef rec
238#undef clique
239#undef set
240
241/* eof */
Note: See TracBrowser for help on using the repository browser.