lemon-project-template-glpk

comparison deps/glpk/src/glpnet08.c @ 11:4fc6ad2fb8a6

find changesets by author, revision, files, or words in the commit message
Test GLPK in src/main.cc
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 21:43:29 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f9221b17a1de
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
72 struct 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
102 static 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
122 static 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 }
163 done: xfree(newtable);
164 return;
165 }
166
167 int 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 */