lemon-project-template-glpk

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