/farmR/src/glpios08.c

https://code.google.com/p/javawfm/ · C · 906 lines · 589 code · 36 blank · 281 comment · 174 complexity · 0fc19d229e90379b35dd17e751095a56 MD5 · raw file

  1. /* glpios08.c (clique cut generator) */
  2. /***********************************************************************
  3. * This code is part of GLPK (GNU Linear Programming Kit).
  4. *
  5. * Copyright (C) 2000,01,02,03,04,05,06,07,08,2009 Andrew Makhorin,
  6. * Department for Applied Informatics, Moscow Aviation Institute,
  7. * Moscow, Russia. All rights reserved. E-mail: <mao@mai2.rcnet.ru>.
  8. *
  9. * GLPK is free software: you can redistribute it and/or modify it
  10. * under the terms of the GNU General Public License as published by
  11. * the Free Software Foundation, either version 3 of the License, or
  12. * (at your option) any later version.
  13. *
  14. * GLPK is distributed in the hope that it will be useful, but WITHOUT
  15. * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  16. * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
  17. * License for more details.
  18. *
  19. * You should have received a copy of the GNU General Public License
  20. * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
  21. ***********************************************************************/
  22. #include "glpios.h"
  23. static double get_row_lb(LPX *lp, int i)
  24. { /* this routine returns lower bound of row i or -DBL_MAX if the
  25. row has no lower bound */
  26. double lb;
  27. switch (lpx_get_row_type(lp, i))
  28. { case LPX_FR:
  29. case LPX_UP:
  30. lb = -DBL_MAX;
  31. break;
  32. case LPX_LO:
  33. case LPX_DB:
  34. case LPX_FX:
  35. lb = lpx_get_row_lb(lp, i);
  36. break;
  37. default:
  38. xassert(lp != lp);
  39. }
  40. return lb;
  41. }
  42. static double get_row_ub(LPX *lp, int i)
  43. { /* this routine returns upper bound of row i or +DBL_MAX if the
  44. row has no upper bound */
  45. double ub;
  46. switch (lpx_get_row_type(lp, i))
  47. { case LPX_FR:
  48. case LPX_LO:
  49. ub = +DBL_MAX;
  50. break;
  51. case LPX_UP:
  52. case LPX_DB:
  53. case LPX_FX:
  54. ub = lpx_get_row_ub(lp, i);
  55. break;
  56. default:
  57. xassert(lp != lp);
  58. }
  59. return ub;
  60. }
  61. static double get_col_lb(LPX *lp, int j)
  62. { /* this routine returns lower bound of column j or -DBL_MAX if
  63. the column has no lower bound */
  64. double lb;
  65. switch (lpx_get_col_type(lp, j))
  66. { case LPX_FR:
  67. case LPX_UP:
  68. lb = -DBL_MAX;
  69. break;
  70. case LPX_LO:
  71. case LPX_DB:
  72. case LPX_FX:
  73. lb = lpx_get_col_lb(lp, j);
  74. break;
  75. default:
  76. xassert(lp != lp);
  77. }
  78. return lb;
  79. }
  80. static double get_col_ub(LPX *lp, int j)
  81. { /* this routine returns upper bound of column j or +DBL_MAX if
  82. the column has no upper bound */
  83. double ub;
  84. switch (lpx_get_col_type(lp, j))
  85. { case LPX_FR:
  86. case LPX_LO:
  87. ub = +DBL_MAX;
  88. break;
  89. case LPX_UP:
  90. case LPX_DB:
  91. case LPX_FX:
  92. ub = lpx_get_col_ub(lp, j);
  93. break;
  94. default:
  95. xassert(lp != lp);
  96. }
  97. return ub;
  98. }
  99. static int is_binary(LPX *lp, int j)
  100. { /* this routine checks if variable x[j] is binary */
  101. return
  102. lpx_get_col_kind(lp, j) == LPX_IV &&
  103. lpx_get_col_type(lp, j) == LPX_DB &&
  104. lpx_get_col_lb(lp, j) == 0.0 && lpx_get_col_ub(lp, j) == 1.0;
  105. }
  106. static double eval_lf_min(LPX *lp, int len, int ind[], double val[])
  107. { /* this routine computes the minimum of a specified linear form
  108. sum a[j]*x[j]
  109. j
  110. using the formula:
  111. min = sum a[j]*lb[j] + sum a[j]*ub[j],
  112. j in J+ j in J-
  113. where J+ = {j: a[j] > 0}, J- = {j: a[j] < 0}, lb[j] and ub[j]
  114. are lower and upper bound of variable x[j], resp. */
  115. int j, t;
  116. double lb, ub, sum;
  117. sum = 0.0;
  118. for (t = 1; t <= len; t++)
  119. { j = ind[t];
  120. if (val[t] > 0.0)
  121. { lb = get_col_lb(lp, j);
  122. if (lb == -DBL_MAX)
  123. { sum = -DBL_MAX;
  124. break;
  125. }
  126. sum += val[t] * lb;
  127. }
  128. else if (val[t] < 0.0)
  129. { ub = get_col_ub(lp, j);
  130. if (ub == +DBL_MAX)
  131. { sum = -DBL_MAX;
  132. break;
  133. }
  134. sum += val[t] * ub;
  135. }
  136. else
  137. xassert(val != val);
  138. }
  139. return sum;
  140. }
  141. static double eval_lf_max(LPX *lp, int len, int ind[], double val[])
  142. { /* this routine computes the maximum of a specified linear form
  143. sum a[j]*x[j]
  144. j
  145. using the formula:
  146. max = sum a[j]*ub[j] + sum a[j]*lb[j],
  147. j in J+ j in J-
  148. where J+ = {j: a[j] > 0}, J- = {j: a[j] < 0}, lb[j] and ub[j]
  149. are lower and upper bound of variable x[j], resp. */
  150. int j, t;
  151. double lb, ub, sum;
  152. sum = 0.0;
  153. for (t = 1; t <= len; t++)
  154. { j = ind[t];
  155. if (val[t] > 0.0)
  156. { ub = get_col_ub(lp, j);
  157. if (ub == +DBL_MAX)
  158. { sum = +DBL_MAX;
  159. break;
  160. }
  161. sum += val[t] * ub;
  162. }
  163. else if (val[t] < 0.0)
  164. { lb = get_col_lb(lp, j);
  165. if (lb == -DBL_MAX)
  166. { sum = +DBL_MAX;
  167. break;
  168. }
  169. sum += val[t] * lb;
  170. }
  171. else
  172. xassert(val != val);
  173. }
  174. return sum;
  175. }
  176. /*----------------------------------------------------------------------
  177. -- probing - determine logical relation between binary variables.
  178. --
  179. -- This routine tentatively sets a binary variable to 0 and then to 1
  180. -- and examines whether another binary variable is caused to be fixed.
  181. --
  182. -- The examination is based only on one row (constraint), which is the
  183. -- following:
  184. --
  185. -- L <= sum a[j]*x[j] <= U. (1)
  186. -- j
  187. --
  188. -- Let x[p] be a probing variable, x[q] be an examined variable. Then
  189. -- (1) can be written as:
  190. --
  191. -- L <= sum a[j]*x[j] + a[p]*x[p] + a[q]*x[q] <= U, (2)
  192. -- j in J'
  193. --
  194. -- where J' = {j: j != p and j != q}.
  195. --
  196. -- Let
  197. --
  198. -- L' = L - a[p]*x[p], (3)
  199. --
  200. -- U' = U - a[p]*x[p], (4)
  201. --
  202. -- where x[p] is assumed to be fixed at 0 or 1. So (2) can be rewritten
  203. -- as follows:
  204. --
  205. -- L' <= sum a[j]*x[j] + a[q]*x[q] <= U', (5)
  206. -- j in J'
  207. --
  208. -- from where we have:
  209. --
  210. -- L' - sum a[j]*x[j] <= a[q]*x[q] <= U' - sum a[j]*x[j]. (6)
  211. -- j in J' j in J'
  212. --
  213. -- Thus,
  214. --
  215. -- min a[q]*x[q] = L' - MAX, (7)
  216. --
  217. -- max a[q]*x[q] = U' - MIN, (8)
  218. --
  219. -- where
  220. --
  221. -- MIN = min sum a[j]*x[j], (9)
  222. -- j in J'
  223. --
  224. -- MAX = max sum a[j]*x[j]. (10)
  225. -- j in J'
  226. --
  227. -- Formulae (7) and (8) allows determining implied lower and upper
  228. -- bounds of x[q].
  229. --
  230. -- Parameters len, val, L and U specify the constraint (1).
  231. --
  232. -- Parameters lf_min and lf_max specify implied lower and upper bounds
  233. -- of the linear form (1). It is assumed that these bounds are computed
  234. -- with the routines eval_lf_min and eval_lf_max (see above).
  235. --
  236. -- Parameter p specifies the probing variable x[p], which is set to 0
  237. -- (if set is 0) or to 1 (if set is 1).
  238. --
  239. -- Parameter q specifies the examined variable x[q].
  240. --
  241. -- On exit the routine returns one of the following codes:
  242. --
  243. -- 0 - there is no logical relation between x[p] and x[q];
  244. -- 1 - x[q] can take only on value 0;
  245. -- 2 - x[q] can take only on value 1. */
  246. static int probing(int len, double val[], double L, double U,
  247. double lf_min, double lf_max, int p, int set, int q)
  248. { double temp;
  249. xassert(1 <= p && p < q && q <= len);
  250. /* compute L' (3) */
  251. if (L != -DBL_MAX && set) L -= val[p];
  252. /* compute U' (4) */
  253. if (U != +DBL_MAX && set) U -= val[p];
  254. /* compute MIN (9) */
  255. if (lf_min != -DBL_MAX)
  256. { if (val[p] < 0.0) lf_min -= val[p];
  257. if (val[q] < 0.0) lf_min -= val[q];
  258. }
  259. /* compute MAX (10) */
  260. if (lf_max != +DBL_MAX)
  261. { if (val[p] > 0.0) lf_max -= val[p];
  262. if (val[q] > 0.0) lf_max -= val[q];
  263. }
  264. /* compute implied lower bound of x[q]; see (7), (8) */
  265. if (val[q] > 0.0)
  266. { if (L == -DBL_MAX || lf_max == +DBL_MAX)
  267. temp = -DBL_MAX;
  268. else
  269. temp = (L - lf_max) / val[q];
  270. }
  271. else
  272. { if (U == +DBL_MAX || lf_min == -DBL_MAX)
  273. temp = -DBL_MAX;
  274. else
  275. temp = (U - lf_min) / val[q];
  276. }
  277. if (temp > 0.001) return 2;
  278. /* compute implied upper bound of x[q]; see (7), (8) */
  279. if (val[q] > 0.0)
  280. { if (U == +DBL_MAX || lf_min == -DBL_MAX)
  281. temp = +DBL_MAX;
  282. else
  283. temp = (U - lf_min) / val[q];
  284. }
  285. else
  286. { if (L == -DBL_MAX || lf_max == +DBL_MAX)
  287. temp = +DBL_MAX;
  288. else
  289. temp = (L - lf_max) / val[q];
  290. }
  291. if (temp < 0.999) return 1;
  292. /* there is no logical relation between x[p] and x[q] */
  293. return 0;
  294. }
  295. struct COG
  296. { /* conflict graph; it represents logical relations between binary
  297. variables and has a vertex for each binary variable and its
  298. complement, and an edge between two vertices when at most one
  299. of the variables represented by the vertices can equal one in
  300. an optimal solution */
  301. int n;
  302. /* number of variables */
  303. int nb;
  304. /* number of binary variables represented in the graph (note that
  305. not all binary variables can be represented); vertices which
  306. correspond to binary variables have numbers 1, ..., nb while
  307. vertices which correspond to complements of binary variables
  308. have numbers nb+1, ..., nb+nb */
  309. int ne;
  310. /* number of edges in the graph */
  311. int *vert; /* int vert[1+n]; */
  312. /* if x[j] is a binary variable represented in the graph, vert[j]
  313. is the vertex number corresponding to x[j]; otherwise vert[j]
  314. is zero */
  315. int *orig; /* int list[1:nb]; */
  316. /* if vert[j] = k > 0, then orig[k] = j */
  317. unsigned char *a;
  318. /* adjacency matrix of the graph having 2*nb rows and columns;
  319. only strict lower triangle is stored in dense packed form */
  320. };
  321. /*----------------------------------------------------------------------
  322. -- lpx_create_cog - create the conflict graph.
  323. --
  324. -- SYNOPSIS
  325. --
  326. -- #include "glplpx.h"
  327. -- void *lpx_create_cog(LPX *lp);
  328. --
  329. -- DESCRIPTION
  330. --
  331. -- The routine lpx_create_cog creates the conflict graph for a given
  332. -- problem instance.
  333. --
  334. -- RETURNS
  335. --
  336. -- If the graph has been created, the routine returns a pointer to it.
  337. -- Otherwise the routine returns NULL. */
  338. #define MAX_NB 4000
  339. #define MAX_ROW_LEN 500
  340. static void lpx_add_cog_edge(void *_cog, int i, int j);
  341. static void *lpx_create_cog(LPX *lp)
  342. { struct COG *cog = NULL;
  343. int m, n, nb, i, j, p, q, len, *ind, *vert, *orig;
  344. double L, U, lf_min, lf_max, *val;
  345. xprintf("Creating the conflict graph...\n");
  346. m = lpx_get_num_rows(lp);
  347. n = lpx_get_num_cols(lp);
  348. /* determine which binary variables should be included in the
  349. conflict graph */
  350. nb = 0;
  351. vert = xcalloc(1+n, sizeof(int));
  352. for (j = 1; j <= n; j++) vert[j] = 0;
  353. orig = xcalloc(1+n, sizeof(int));
  354. ind = xcalloc(1+n, sizeof(int));
  355. val = xcalloc(1+n, sizeof(double));
  356. for (i = 1; i <= m; i++)
  357. { L = get_row_lb(lp, i);
  358. U = get_row_ub(lp, i);
  359. if (L == -DBL_MAX && U == +DBL_MAX) continue;
  360. len = lpx_get_mat_row(lp, i, ind, val);
  361. if (len > MAX_ROW_LEN) continue;
  362. lf_min = eval_lf_min(lp, len, ind, val);
  363. lf_max = eval_lf_max(lp, len, ind, val);
  364. for (p = 1; p <= len; p++)
  365. { if (!is_binary(lp, ind[p])) continue;
  366. for (q = p+1; q <= len; q++)
  367. { if (!is_binary(lp, ind[q])) continue;
  368. if (probing(len, val, L, U, lf_min, lf_max, p, 0, q) ||
  369. probing(len, val, L, U, lf_min, lf_max, p, 1, q))
  370. { /* there is a logical relation */
  371. /* include the first variable in the graph */
  372. j = ind[p];
  373. if (vert[j] == 0) nb++, vert[j] = nb, orig[nb] = j;
  374. /* incude the second variable in the graph */
  375. j = ind[q];
  376. if (vert[j] == 0) nb++, vert[j] = nb, orig[nb] = j;
  377. }
  378. }
  379. }
  380. }
  381. /* if the graph is either empty or has too many vertices, do not
  382. create it */
  383. if (nb == 0 || nb > MAX_NB)
  384. { xprintf("The conflict graph is either empty or too big\n");
  385. xfree(vert);
  386. xfree(orig);
  387. goto done;
  388. }
  389. /* create the conflict graph */
  390. cog = xmalloc(sizeof(struct COG));
  391. cog->n = n;
  392. cog->nb = nb;
  393. cog->ne = 0;
  394. cog->vert = vert;
  395. cog->orig = orig;
  396. len = nb + nb; /* number of vertices */
  397. len = (len * (len - 1)) / 2; /* number of entries in triangle */
  398. len = (len + (CHAR_BIT - 1)) / CHAR_BIT; /* bytes needed */
  399. cog->a = xmalloc(len);
  400. memset(cog->a, 0, len);
  401. for (j = 1; j <= nb; j++)
  402. { /* add edge between variable and its complement */
  403. lpx_add_cog_edge(cog, +orig[j], -orig[j]);
  404. }
  405. for (i = 1; i <= m; i++)
  406. { L = get_row_lb(lp, i);
  407. U = get_row_ub(lp, i);
  408. if (L == -DBL_MAX && U == +DBL_MAX) continue;
  409. len = lpx_get_mat_row(lp, i, ind, val);
  410. if (len > MAX_ROW_LEN) continue;
  411. lf_min = eval_lf_min(lp, len, ind, val);
  412. lf_max = eval_lf_max(lp, len, ind, val);
  413. for (p = 1; p <= len; p++)
  414. { if (!is_binary(lp, ind[p])) continue;
  415. for (q = p+1; q <= len; q++)
  416. { if (!is_binary(lp, ind[q])) continue;
  417. /* set x[p] to 0 and examine x[q] */
  418. switch (probing(len, val, L, U, lf_min, lf_max, p, 0, q))
  419. { case 0:
  420. /* no logical relation */
  421. break;
  422. case 1:
  423. /* x[p] = 0 implies x[q] = 0 */
  424. lpx_add_cog_edge(cog, -ind[p], +ind[q]);
  425. break;
  426. case 2:
  427. /* x[p] = 0 implies x[q] = 1 */
  428. lpx_add_cog_edge(cog, -ind[p], -ind[q]);
  429. break;
  430. default:
  431. xassert(lp != lp);
  432. }
  433. /* set x[p] to 1 and examine x[q] */
  434. switch (probing(len, val, L, U, lf_min, lf_max, p, 1, q))
  435. { case 0:
  436. /* no logical relation */
  437. break;
  438. case 1:
  439. /* x[p] = 1 implies x[q] = 0 */
  440. lpx_add_cog_edge(cog, +ind[p], +ind[q]);
  441. break;
  442. case 2:
  443. /* x[p] = 1 implies x[q] = 1 */
  444. lpx_add_cog_edge(cog, +ind[p], -ind[q]);
  445. break;
  446. default:
  447. xassert(lp != lp);
  448. }
  449. }
  450. }
  451. }
  452. xprintf("The conflict graph has 2*%d vertices and %d edges\n",
  453. cog->nb, cog->ne);
  454. done: xfree(ind);
  455. xfree(val);
  456. return cog;
  457. }
  458. /*----------------------------------------------------------------------
  459. -- lpx_add_cog_edge - add edge to the conflict graph.
  460. --
  461. -- SYNOPSIS
  462. --
  463. -- #include "glplpx.h"
  464. -- void lpx_add_cog_edge(void *cog, int i, int j);
  465. --
  466. -- DESCRIPTION
  467. --
  468. -- The routine lpx_add_cog_edge adds an edge to the conflict graph.
  469. -- The edge connects x[i] (if i > 0) or its complement (if i < 0) and
  470. -- x[j] (if j > 0) or its complement (if j < 0), where i and j are
  471. -- original ordinal numbers of corresponding variables. */
  472. static void lpx_add_cog_edge(void *_cog, int i, int j)
  473. { struct COG *cog = _cog;
  474. int k;
  475. xassert(i != j);
  476. /* determine indices of corresponding vertices */
  477. if (i > 0)
  478. { xassert(1 <= i && i <= cog->n);
  479. i = cog->vert[i];
  480. xassert(i != 0);
  481. }
  482. else
  483. { i = -i;
  484. xassert(1 <= i && i <= cog->n);
  485. i = cog->vert[i];
  486. xassert(i != 0);
  487. i += cog->nb;
  488. }
  489. if (j > 0)
  490. { xassert(1 <= j && j <= cog->n);
  491. j = cog->vert[j];
  492. xassert(j != 0);
  493. }
  494. else
  495. { j = -j;
  496. xassert(1 <= j && j <= cog->n);
  497. j = cog->vert[j];
  498. xassert(j != 0);
  499. j += cog->nb;
  500. }
  501. /* only lower triangle is stored, so we need i > j */
  502. if (i < j) k = i, i = j, j = k;
  503. k = ((i - 1) * (i - 2)) / 2 + (j - 1);
  504. cog->a[k / CHAR_BIT] |=
  505. (unsigned char)(1 << ((CHAR_BIT - 1) - k % CHAR_BIT));
  506. cog->ne++;
  507. return;
  508. }
  509. /*----------------------------------------------------------------------
  510. -- MAXIMUM WEIGHT CLIQUE
  511. --
  512. -- Two subroutines sub() and wclique() below are intended to find a
  513. -- maximum weight clique in a given undirected graph. These subroutines
  514. -- are slightly modified version of the program WCLIQUE developed by
  515. -- Patric Ostergard <http://www.tcs.hut.fi/~pat/wclique.html> and based
  516. -- on ideas from the article "P. R. J. Ostergard, A new algorithm for
  517. -- the maximum-weight clique problem, submitted for publication", which
  518. -- in turn is a generalization of the algorithm for unweighted graphs
  519. -- presented in "P. R. J. Ostergard, A fast algorithm for the maximum
  520. -- clique problem, submitted for publication".
  521. --
  522. -- USED WITH PERMISSION OF THE AUTHOR OF THE ORIGINAL CODE. */
  523. struct dsa
  524. { /* dynamic storage area */
  525. int n;
  526. /* number of vertices */
  527. int *wt; /* int wt[0:n-1]; */
  528. /* weights */
  529. unsigned char *a;
  530. /* adjacency matrix (packed lower triangle without main diag.) */
  531. int record;
  532. /* weight of best clique */
  533. int rec_level;
  534. /* number of vertices in best clique */
  535. int *rec; /* int rec[0:n-1]; */
  536. /* best clique so far */
  537. int *clique; /* int clique[0:n-1]; */
  538. /* table for pruning */
  539. int *set; /* int set[0:n-1]; */
  540. /* current clique */
  541. };
  542. #define n (dsa->n)
  543. #define wt (dsa->wt)
  544. #define a (dsa->a)
  545. #define record (dsa->record)
  546. #define rec_level (dsa->rec_level)
  547. #define rec (dsa->rec)
  548. #define clique (dsa->clique)
  549. #define set (dsa->set)
  550. #if 0
  551. static int is_edge(struct dsa *dsa, int i, int j)
  552. { /* if there is arc (i,j), the routine returns true; otherwise
  553. false; 0 <= i, j < n */
  554. int k;
  555. xassert(0 <= i && i < n);
  556. xassert(0 <= j && j < n);
  557. if (i == j) return 0;
  558. if (i < j) k = i, i = j, j = k;
  559. k = (i * (i - 1)) / 2 + j;
  560. return a[k / CHAR_BIT] &
  561. (unsigned char)(1 << ((CHAR_BIT - 1) - k % CHAR_BIT));
  562. }
  563. #else
  564. #define is_edge(dsa, i, j) ((i) == (j) ? 0 : \
  565. (i) > (j) ? is_edge1(i, j) : is_edge1(j, i))
  566. #define is_edge1(i, j) is_edge2(((i) * ((i) - 1)) / 2 + (j))
  567. #define is_edge2(k) (a[(k) / CHAR_BIT] & \
  568. (unsigned char)(1 << ((CHAR_BIT - 1) - (k) % CHAR_BIT)))
  569. #endif
  570. static void sub(struct dsa *dsa, int ct, int table[], int level,
  571. int weight, int l_weight)
  572. { int i, j, k, curr_weight, left_weight, *p1, *p2, *newtable;
  573. newtable = xcalloc(n, sizeof(int));
  574. if (ct <= 0)
  575. { /* 0 or 1 elements left; include these */
  576. if (ct == 0)
  577. { set[level++] = table[0];
  578. weight += l_weight;
  579. }
  580. if (weight > record)
  581. { record = weight;
  582. rec_level = level;
  583. for (i = 0; i < level; i++) rec[i] = set[i];
  584. }
  585. goto done;
  586. }
  587. for (i = ct; i >= 0; i--)
  588. { if ((level == 0) && (i < ct)) goto done;
  589. k = table[i];
  590. if ((level > 0) && (clique[k] <= (record - weight)))
  591. goto done; /* prune */
  592. set[level] = k;
  593. curr_weight = weight + wt[k];
  594. l_weight -= wt[k];
  595. if (l_weight <= (record - curr_weight))
  596. goto done; /* prune */
  597. p1 = newtable;
  598. p2 = table;
  599. left_weight = 0;
  600. while (p2 < table + i)
  601. { j = *p2++;
  602. if (is_edge(dsa, j, k))
  603. { *p1++ = j;
  604. left_weight += wt[j];
  605. }
  606. }
  607. if (left_weight <= (record - curr_weight)) continue;
  608. sub(dsa, p1 - newtable - 1, newtable, level + 1, curr_weight,
  609. left_weight);
  610. }
  611. done: xfree(newtable);
  612. return;
  613. }
  614. static int wclique(int _n, int w[], unsigned char _a[], int sol[])
  615. { struct dsa _dsa, *dsa = &_dsa;
  616. int i, j, p, max_wt, max_nwt, wth, *used, *nwt, *pos;
  617. xlong_t timer;
  618. n = _n;
  619. wt = &w[1];
  620. a = _a;
  621. record = 0;
  622. rec_level = 0;
  623. rec = &sol[1];
  624. clique = xcalloc(n, sizeof(int));
  625. set = xcalloc(n, sizeof(int));
  626. used = xcalloc(n, sizeof(int));
  627. nwt = xcalloc(n, sizeof(int));
  628. pos = xcalloc(n, sizeof(int));
  629. /* start timer */
  630. timer = xtime();
  631. /* order vertices */
  632. for (i = 0; i < n; i++)
  633. { nwt[i] = 0;
  634. for (j = 0; j < n; j++)
  635. if (is_edge(dsa, i, j)) nwt[i] += wt[j];
  636. }
  637. for (i = 0; i < n; i++)
  638. used[i] = 0;
  639. for (i = n-1; i >= 0; i--)
  640. { max_wt = -1;
  641. max_nwt = -1;
  642. for (j = 0; j < n; j++)
  643. { if ((!used[j]) && ((wt[j] > max_wt) || (wt[j] == max_wt
  644. && nwt[j] > max_nwt)))
  645. { max_wt = wt[j];
  646. max_nwt = nwt[j];
  647. p = j;
  648. }
  649. }
  650. pos[i] = p;
  651. used[p] = 1;
  652. for (j = 0; j < n; j++)
  653. if ((!used[j]) && (j != p) && (is_edge(dsa, p, j)))
  654. nwt[j] -= wt[p];
  655. }
  656. /* main routine */
  657. wth = 0;
  658. for (i = 0; i < n; i++)
  659. { wth += wt[pos[i]];
  660. sub(dsa, i, pos, 0, 0, wth);
  661. clique[pos[i]] = record;
  662. #if 0
  663. if (utime() >= timer + 5.0)
  664. #else
  665. if (xdifftime(xtime(), timer) >= 5.0 - 0.001)
  666. #endif
  667. { /* print current record and reset timer */
  668. xprintf("level = %d (%d); best = %d\n", i+1, n, record);
  669. #if 0
  670. timer = utime();
  671. #else
  672. timer = xtime();
  673. #endif
  674. }
  675. }
  676. xfree(clique);
  677. xfree(set);
  678. xfree(used);
  679. xfree(nwt);
  680. xfree(pos);
  681. /* return the solution found */
  682. for (i = 1; i <= rec_level; i++) sol[i]++;
  683. return rec_level;
  684. }
  685. #undef n
  686. #undef wt
  687. #undef a
  688. #undef record
  689. #undef rec_level
  690. #undef rec
  691. #undef clique
  692. #undef set
  693. /*----------------------------------------------------------------------
  694. -- lpx_clique_cut - generate cluque cut.
  695. --
  696. -- SYNOPSIS
  697. --
  698. -- #include "glplpx.h"
  699. -- int lpx_clique_cut(LPX *lp, void *cog, int ind[], double val[]);
  700. --
  701. -- DESCRIPTION
  702. --
  703. -- The routine lpx_clique_cut generates a clique cut using the conflict
  704. -- graph specified by the parameter cog.
  705. --
  706. -- If a violated clique cut has been found, it has the following form:
  707. --
  708. -- sum{j in J} a[j]*x[j] <= b.
  709. --
  710. -- Variable indices j in J are stored in elements ind[1], ..., ind[len]
  711. -- while corresponding constraint coefficients are stored in elements
  712. -- val[1], ..., val[len], where len is returned on exit. The right-hand
  713. -- side b is stored in element val[0].
  714. --
  715. -- RETURNS
  716. --
  717. -- If the cutting plane has been successfully generated, the routine
  718. -- returns 1 <= len <= n, which is the number of non-zero coefficients
  719. -- in the inequality constraint. Otherwise, the routine returns zero. */
  720. static int lpx_clique_cut(LPX *lp, void *_cog, int ind[], double val[])
  721. { struct COG *cog = _cog;
  722. int n = lpx_get_num_cols(lp);
  723. int j, t, v, card, temp, len = 0, *w, *sol;
  724. double x, sum, b, *vec;
  725. /* allocate working arrays */
  726. w = xcalloc(1 + 2 * cog->nb, sizeof(int));
  727. sol = xcalloc(1 + 2 * cog->nb, sizeof(int));
  728. vec = xcalloc(1+n, sizeof(double));
  729. /* assign weights to vertices of the conflict graph */
  730. for (t = 1; t <= cog->nb; t++)
  731. { j = cog->orig[t];
  732. x = lpx_get_col_prim(lp, j);
  733. temp = (int)(100.0 * x + 0.5);
  734. if (temp < 0) temp = 0;
  735. if (temp > 100) temp = 100;
  736. w[t] = temp;
  737. w[cog->nb + t] = 100 - temp;
  738. }
  739. /* find a clique of maximum weight */
  740. card = wclique(2 * cog->nb, w, cog->a, sol);
  741. /* compute the clique weight for unscaled values */
  742. sum = 0.0;
  743. for ( t = 1; t <= card; t++)
  744. { v = sol[t];
  745. xassert(1 <= v && v <= 2 * cog->nb);
  746. if (v <= cog->nb)
  747. { /* vertex v corresponds to binary variable x[j] */
  748. j = cog->orig[v];
  749. x = lpx_get_col_prim(lp, j);
  750. sum += x;
  751. }
  752. else
  753. { /* vertex v corresponds to the complement of x[j] */
  754. j = cog->orig[v - cog->nb];
  755. x = lpx_get_col_prim(lp, j);
  756. sum += 1.0 - x;
  757. }
  758. }
  759. /* if the sum of binary variables and their complements in the
  760. clique greater than 1, the clique cut is violated */
  761. if (sum >= 1.01)
  762. { /* construct the inquality */
  763. for (j = 1; j <= n; j++) vec[j] = 0;
  764. b = 1.0;
  765. for (t = 1; t <= card; t++)
  766. { v = sol[t];
  767. if (v <= cog->nb)
  768. { /* vertex v corresponds to binary variable x[j] */
  769. j = cog->orig[v];
  770. xassert(1 <= j && j <= n);
  771. vec[j] += 1.0;
  772. }
  773. else
  774. { /* vertex v corresponds to the complement of x[j] */
  775. j = cog->orig[v - cog->nb];
  776. xassert(1 <= j && j <= n);
  777. vec[j] -= 1.0;
  778. b -= 1.0;
  779. }
  780. }
  781. xassert(len == 0);
  782. for (j = 1; j <= n; j++)
  783. { if (vec[j] != 0.0)
  784. { len++;
  785. ind[len] = j, val[len] = vec[j];
  786. }
  787. }
  788. ind[0] = 0, val[0] = b;
  789. }
  790. /* free working arrays */
  791. xfree(w);
  792. xfree(sol);
  793. xfree(vec);
  794. /* return to the calling program */
  795. return len;
  796. }
  797. /*----------------------------------------------------------------------
  798. -- lpx_delete_cog - delete the conflict graph.
  799. --
  800. -- SYNOPSIS
  801. --
  802. -- #include "glplpx.h"
  803. -- void lpx_delete_cog(void *cog);
  804. --
  805. -- DESCRIPTION
  806. --
  807. -- The routine lpx_delete_cog deletes the conflict graph, which the
  808. -- parameter cog points to, freeing all the memory allocated to this
  809. -- object. */
  810. static void lpx_delete_cog(void *_cog)
  811. { struct COG *cog = _cog;
  812. xfree(cog->vert);
  813. xfree(cog->orig);
  814. xfree(cog->a);
  815. xfree(cog);
  816. }
  817. /**********************************************************************/
  818. void *ios_clq_init(glp_tree *tree)
  819. { /* initialize clique cut generator */
  820. glp_prob *mip = tree->mip;
  821. xassert(mip != NULL);
  822. return lpx_create_cog(mip);
  823. }
  824. /***********************************************************************
  825. * NAME
  826. *
  827. * ios_clq_gen - generate clique cuts
  828. *
  829. * SYNOPSIS
  830. *
  831. * #include "glpios.h"
  832. * void ios_clq_gen(glp_tree *tree, void *gen);
  833. *
  834. * DESCRIPTION
  835. *
  836. * The routine ios_clq_gen generates clique cuts for the current point
  837. * and adds them to the clique pool. */
  838. void ios_clq_gen(glp_tree *tree, void *gen)
  839. { int n = lpx_get_num_cols(tree->mip);
  840. int len, *ind;
  841. double *val;
  842. xassert(gen != NULL);
  843. ind = xcalloc(1+n, sizeof(int));
  844. val = xcalloc(1+n, sizeof(double));
  845. len = lpx_clique_cut(tree->mip, gen, ind, val);
  846. if (len > 0)
  847. { /* xprintf("len = %d\n", len); */
  848. glp_ios_add_row(tree, NULL, GLP_RF_CLQ, 0, len, ind, val,
  849. GLP_UP, val[0]);
  850. }
  851. xfree(ind);
  852. xfree(val);
  853. return;
  854. }
  855. /**********************************************************************/
  856. void ios_clq_term(void *gen)
  857. { /* terminate clique cut generator */
  858. xassert(gen != NULL);
  859. lpx_delete_cog(gen);
  860. return;
  861. }
  862. /* eof */