diff -r d59bea55db9b -r c445c931472f src/glpios08.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/glpios08.c Mon Dec 06 13:09:21 2010 +0100 @@ -0,0 +1,907 @@ +/* glpios08.c (clique cut generator) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, +* 2009, 2010 Andrew Makhorin, Department for Applied Informatics, +* Moscow Aviation Institute, Moscow, Russia. All rights reserved. +* E-mail: . +* +* 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 . +***********************************************************************/ + +#include "glpios.h" + +static double get_row_lb(LPX *lp, int i) +{ /* this routine returns lower bound of row i or -DBL_MAX if the + row has no lower bound */ + double lb; + switch (lpx_get_row_type(lp, i)) + { case LPX_FR: + case LPX_UP: + lb = -DBL_MAX; + break; + case LPX_LO: + case LPX_DB: + case LPX_FX: + lb = lpx_get_row_lb(lp, i); + break; + default: + xassert(lp != lp); + } + return lb; +} + +static double get_row_ub(LPX *lp, int i) +{ /* this routine returns upper bound of row i or +DBL_MAX if the + row has no upper bound */ + double ub; + switch (lpx_get_row_type(lp, i)) + { case LPX_FR: + case LPX_LO: + ub = +DBL_MAX; + break; + case LPX_UP: + case LPX_DB: + case LPX_FX: + ub = lpx_get_row_ub(lp, i); + break; + default: + xassert(lp != lp); + } + return ub; +} + +static double get_col_lb(LPX *lp, int j) +{ /* this routine returns lower bound of column j or -DBL_MAX if + the column has no lower bound */ + double lb; + switch (lpx_get_col_type(lp, j)) + { case LPX_FR: + case LPX_UP: + lb = -DBL_MAX; + break; + case LPX_LO: + case LPX_DB: + case LPX_FX: + lb = lpx_get_col_lb(lp, j); + break; + default: + xassert(lp != lp); + } + return lb; +} + +static double get_col_ub(LPX *lp, int j) +{ /* this routine returns upper bound of column j or +DBL_MAX if + the column has no upper bound */ + double ub; + switch (lpx_get_col_type(lp, j)) + { case LPX_FR: + case LPX_LO: + ub = +DBL_MAX; + break; + case LPX_UP: + case LPX_DB: + case LPX_FX: + ub = lpx_get_col_ub(lp, j); + break; + default: + xassert(lp != lp); + } + return ub; +} + +static int is_binary(LPX *lp, int j) +{ /* this routine checks if variable x[j] is binary */ + return + lpx_get_col_kind(lp, j) == LPX_IV && + lpx_get_col_type(lp, j) == LPX_DB && + lpx_get_col_lb(lp, j) == 0.0 && lpx_get_col_ub(lp, j) == 1.0; +} + +static double eval_lf_min(LPX *lp, int len, int ind[], double val[]) +{ /* this routine computes the minimum of a specified linear form + + sum a[j]*x[j] + j + + using the formula: + + min = sum a[j]*lb[j] + sum a[j]*ub[j], + j in J+ j in J- + + where J+ = {j: a[j] > 0}, J- = {j: a[j] < 0}, lb[j] and ub[j] + are lower and upper bound of variable x[j], resp. */ + int j, t; + double lb, ub, sum; + sum = 0.0; + for (t = 1; t <= len; t++) + { j = ind[t]; + if (val[t] > 0.0) + { lb = get_col_lb(lp, j); + if (lb == -DBL_MAX) + { sum = -DBL_MAX; + break; + } + sum += val[t] * lb; + } + else if (val[t] < 0.0) + { ub = get_col_ub(lp, j); + if (ub == +DBL_MAX) + { sum = -DBL_MAX; + break; + } + sum += val[t] * ub; + } + else + xassert(val != val); + } + return sum; +} + +static double eval_lf_max(LPX *lp, int len, int ind[], double val[]) +{ /* this routine computes the maximum of a specified linear form + + sum a[j]*x[j] + j + + using the formula: + + max = sum a[j]*ub[j] + sum a[j]*lb[j], + j in J+ j in J- + + where J+ = {j: a[j] > 0}, J- = {j: a[j] < 0}, lb[j] and ub[j] + are lower and upper bound of variable x[j], resp. */ + int j, t; + double lb, ub, sum; + sum = 0.0; + for (t = 1; t <= len; t++) + { j = ind[t]; + if (val[t] > 0.0) + { ub = get_col_ub(lp, j); + if (ub == +DBL_MAX) + { sum = +DBL_MAX; + break; + } + sum += val[t] * ub; + } + else if (val[t] < 0.0) + { lb = get_col_lb(lp, j); + if (lb == -DBL_MAX) + { sum = +DBL_MAX; + break; + } + sum += val[t] * lb; + } + else + xassert(val != val); + } + return sum; +} + +/*---------------------------------------------------------------------- +-- probing - determine logical relation between binary variables. +-- +-- This routine tentatively sets a binary variable to 0 and then to 1 +-- and examines whether another binary variable is caused to be fixed. +-- +-- The examination is based only on one row (constraint), which is the +-- following: +-- +-- L <= sum a[j]*x[j] <= U. (1) +-- j +-- +-- Let x[p] be a probing variable, x[q] be an examined variable. Then +-- (1) can be written as: +-- +-- L <= sum a[j]*x[j] + a[p]*x[p] + a[q]*x[q] <= U, (2) +-- j in J' +-- +-- where J' = {j: j != p and j != q}. +-- +-- Let +-- +-- L' = L - a[p]*x[p], (3) +-- +-- U' = U - a[p]*x[p], (4) +-- +-- where x[p] is assumed to be fixed at 0 or 1. So (2) can be rewritten +-- as follows: +-- +-- L' <= sum a[j]*x[j] + a[q]*x[q] <= U', (5) +-- j in J' +-- +-- from where we have: +-- +-- L' - sum a[j]*x[j] <= a[q]*x[q] <= U' - sum a[j]*x[j]. (6) +-- j in J' j in J' +-- +-- Thus, +-- +-- min a[q]*x[q] = L' - MAX, (7) +-- +-- max a[q]*x[q] = U' - MIN, (8) +-- +-- where +-- +-- MIN = min sum a[j]*x[j], (9) +-- j in J' +-- +-- MAX = max sum a[j]*x[j]. (10) +-- j in J' +-- +-- Formulae (7) and (8) allows determining implied lower and upper +-- bounds of x[q]. +-- +-- Parameters len, val, L and U specify the constraint (1). +-- +-- Parameters lf_min and lf_max specify implied lower and upper bounds +-- of the linear form (1). It is assumed that these bounds are computed +-- with the routines eval_lf_min and eval_lf_max (see above). +-- +-- Parameter p specifies the probing variable x[p], which is set to 0 +-- (if set is 0) or to 1 (if set is 1). +-- +-- Parameter q specifies the examined variable x[q]. +-- +-- On exit the routine returns one of the following codes: +-- +-- 0 - there is no logical relation between x[p] and x[q]; +-- 1 - x[q] can take only on value 0; +-- 2 - x[q] can take only on value 1. */ + +static int probing(int len, double val[], double L, double U, + double lf_min, double lf_max, int p, int set, int q) +{ double temp; + xassert(1 <= p && p < q && q <= len); + /* compute L' (3) */ + if (L != -DBL_MAX && set) L -= val[p]; + /* compute U' (4) */ + if (U != +DBL_MAX && set) U -= val[p]; + /* compute MIN (9) */ + if (lf_min != -DBL_MAX) + { if (val[p] < 0.0) lf_min -= val[p]; + if (val[q] < 0.0) lf_min -= val[q]; + } + /* compute MAX (10) */ + if (lf_max != +DBL_MAX) + { if (val[p] > 0.0) lf_max -= val[p]; + if (val[q] > 0.0) lf_max -= val[q]; + } + /* compute implied lower bound of x[q]; see (7), (8) */ + if (val[q] > 0.0) + { if (L == -DBL_MAX || lf_max == +DBL_MAX) + temp = -DBL_MAX; + else + temp = (L - lf_max) / val[q]; + } + else + { if (U == +DBL_MAX || lf_min == -DBL_MAX) + temp = -DBL_MAX; + else + temp = (U - lf_min) / val[q]; + } + if (temp > 0.001) return 2; + /* compute implied upper bound of x[q]; see (7), (8) */ + if (val[q] > 0.0) + { if (U == +DBL_MAX || lf_min == -DBL_MAX) + temp = +DBL_MAX; + else + temp = (U - lf_min) / val[q]; + } + else + { if (L == -DBL_MAX || lf_max == +DBL_MAX) + temp = +DBL_MAX; + else + temp = (L - lf_max) / val[q]; + } + if (temp < 0.999) return 1; + /* there is no logical relation between x[p] and x[q] */ + return 0; +} + +struct COG +{ /* conflict graph; it represents logical relations between binary + variables and has a vertex for each binary variable and its + complement, and an edge between two vertices when at most one + of the variables represented by the vertices can equal one in + an optimal solution */ + int n; + /* number of variables */ + int nb; + /* number of binary variables represented in the graph (note that + not all binary variables can be represented); vertices which + correspond to binary variables have numbers 1, ..., nb while + vertices which correspond to complements of binary variables + have numbers nb+1, ..., nb+nb */ + int ne; + /* number of edges in the graph */ + int *vert; /* int vert[1+n]; */ + /* if x[j] is a binary variable represented in the graph, vert[j] + is the vertex number corresponding to x[j]; otherwise vert[j] + is zero */ + int *orig; /* int list[1:nb]; */ + /* if vert[j] = k > 0, then orig[k] = j */ + unsigned char *a; + /* adjacency matrix of the graph having 2*nb rows and columns; + only strict lower triangle is stored in dense packed form */ +}; + +/*---------------------------------------------------------------------- +-- lpx_create_cog - create the conflict graph. +-- +-- SYNOPSIS +-- +-- #include "glplpx.h" +-- void *lpx_create_cog(LPX *lp); +-- +-- DESCRIPTION +-- +-- The routine lpx_create_cog creates the conflict graph for a given +-- problem instance. +-- +-- RETURNS +-- +-- If the graph has been created, the routine returns a pointer to it. +-- Otherwise the routine returns NULL. */ + +#define MAX_NB 4000 +#define MAX_ROW_LEN 500 + +static void lpx_add_cog_edge(void *_cog, int i, int j); + +static void *lpx_create_cog(LPX *lp) +{ struct COG *cog = NULL; + int m, n, nb, i, j, p, q, len, *ind, *vert, *orig; + double L, U, lf_min, lf_max, *val; + xprintf("Creating the conflict graph...\n"); + m = lpx_get_num_rows(lp); + n = lpx_get_num_cols(lp); + /* determine which binary variables should be included in the + conflict graph */ + nb = 0; + vert = xcalloc(1+n, sizeof(int)); + for (j = 1; j <= n; j++) vert[j] = 0; + orig = xcalloc(1+n, sizeof(int)); + ind = xcalloc(1+n, sizeof(int)); + val = xcalloc(1+n, sizeof(double)); + for (i = 1; i <= m; i++) + { L = get_row_lb(lp, i); + U = get_row_ub(lp, i); + if (L == -DBL_MAX && U == +DBL_MAX) continue; + len = lpx_get_mat_row(lp, i, ind, val); + if (len > MAX_ROW_LEN) continue; + lf_min = eval_lf_min(lp, len, ind, val); + lf_max = eval_lf_max(lp, len, ind, val); + for (p = 1; p <= len; p++) + { if (!is_binary(lp, ind[p])) continue; + for (q = p+1; q <= len; q++) + { if (!is_binary(lp, ind[q])) continue; + if (probing(len, val, L, U, lf_min, lf_max, p, 0, q) || + probing(len, val, L, U, lf_min, lf_max, p, 1, q)) + { /* there is a logical relation */ + /* include the first variable in the graph */ + j = ind[p]; + if (vert[j] == 0) nb++, vert[j] = nb, orig[nb] = j; + /* incude the second variable in the graph */ + j = ind[q]; + if (vert[j] == 0) nb++, vert[j] = nb, orig[nb] = j; + } + } + } + } + /* if the graph is either empty or has too many vertices, do not + create it */ + if (nb == 0 || nb > MAX_NB) + { xprintf("The conflict graph is either empty or too big\n"); + xfree(vert); + xfree(orig); + goto done; + } + /* create the conflict graph */ + cog = xmalloc(sizeof(struct COG)); + cog->n = n; + cog->nb = nb; + cog->ne = 0; + cog->vert = vert; + cog->orig = orig; + len = nb + nb; /* number of vertices */ + len = (len * (len - 1)) / 2; /* number of entries in triangle */ + len = (len + (CHAR_BIT - 1)) / CHAR_BIT; /* bytes needed */ + cog->a = xmalloc(len); + memset(cog->a, 0, len); + for (j = 1; j <= nb; j++) + { /* add edge between variable and its complement */ + lpx_add_cog_edge(cog, +orig[j], -orig[j]); + } + for (i = 1; i <= m; i++) + { L = get_row_lb(lp, i); + U = get_row_ub(lp, i); + if (L == -DBL_MAX && U == +DBL_MAX) continue; + len = lpx_get_mat_row(lp, i, ind, val); + if (len > MAX_ROW_LEN) continue; + lf_min = eval_lf_min(lp, len, ind, val); + lf_max = eval_lf_max(lp, len, ind, val); + for (p = 1; p <= len; p++) + { if (!is_binary(lp, ind[p])) continue; + for (q = p+1; q <= len; q++) + { if (!is_binary(lp, ind[q])) continue; + /* set x[p] to 0 and examine x[q] */ + switch (probing(len, val, L, U, lf_min, lf_max, p, 0, q)) + { case 0: + /* no logical relation */ + break; + case 1: + /* x[p] = 0 implies x[q] = 0 */ + lpx_add_cog_edge(cog, -ind[p], +ind[q]); + break; + case 2: + /* x[p] = 0 implies x[q] = 1 */ + lpx_add_cog_edge(cog, -ind[p], -ind[q]); + break; + default: + xassert(lp != lp); + } + /* set x[p] to 1 and examine x[q] */ + switch (probing(len, val, L, U, lf_min, lf_max, p, 1, q)) + { case 0: + /* no logical relation */ + break; + case 1: + /* x[p] = 1 implies x[q] = 0 */ + lpx_add_cog_edge(cog, +ind[p], +ind[q]); + break; + case 2: + /* x[p] = 1 implies x[q] = 1 */ + lpx_add_cog_edge(cog, +ind[p], -ind[q]); + break; + default: + xassert(lp != lp); + } + } + } + } + xprintf("The conflict graph has 2*%d vertices and %d edges\n", + cog->nb, cog->ne); +done: xfree(ind); + xfree(val); + return cog; +} + +/*---------------------------------------------------------------------- +-- lpx_add_cog_edge - add edge to the conflict graph. +-- +-- SYNOPSIS +-- +-- #include "glplpx.h" +-- void lpx_add_cog_edge(void *cog, int i, int j); +-- +-- DESCRIPTION +-- +-- The routine lpx_add_cog_edge adds an edge to the conflict graph. +-- The edge connects x[i] (if i > 0) or its complement (if i < 0) and +-- x[j] (if j > 0) or its complement (if j < 0), where i and j are +-- original ordinal numbers of corresponding variables. */ + +static void lpx_add_cog_edge(void *_cog, int i, int j) +{ struct COG *cog = _cog; + int k; + xassert(i != j); + /* determine indices of corresponding vertices */ + if (i > 0) + { xassert(1 <= i && i <= cog->n); + i = cog->vert[i]; + xassert(i != 0); + } + else + { i = -i; + xassert(1 <= i && i <= cog->n); + i = cog->vert[i]; + xassert(i != 0); + i += cog->nb; + } + if (j > 0) + { xassert(1 <= j && j <= cog->n); + j = cog->vert[j]; + xassert(j != 0); + } + else + { j = -j; + xassert(1 <= j && j <= cog->n); + j = cog->vert[j]; + xassert(j != 0); + j += cog->nb; + } + /* only lower triangle is stored, so we need i > j */ + if (i < j) k = i, i = j, j = k; + k = ((i - 1) * (i - 2)) / 2 + (j - 1); + cog->a[k / CHAR_BIT] |= + (unsigned char)(1 << ((CHAR_BIT - 1) - k % CHAR_BIT)); + cog->ne++; + return; +} + +/*---------------------------------------------------------------------- +-- MAXIMUM WEIGHT CLIQUE +-- +-- 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 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. */ + +struct dsa +{ /* dynamic storage area */ + int n; + /* number of vertices */ + int *wt; /* int wt[0:n-1]; */ + /* weights */ + 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 (dsa->n) +#define wt (dsa->wt) +#define a (dsa->a) +#define record (dsa->record) +#define rec_level (dsa->rec_level) +#define rec (dsa->rec) +#define clique (dsa->clique) +#define set (dsa->set) + +#if 0 +static int is_edge(struct dsa *dsa, 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(dsa, 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 dsa *dsa, 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(dsa, j, k)) + { *p1++ = j; + left_weight += wt[j]; + } + } + if (left_weight <= (record - curr_weight)) continue; + sub(dsa, p1 - newtable - 1, newtable, level + 1, curr_weight, + left_weight); + } +done: xfree(newtable); + return; +} + +static int wclique(int _n, int w[], unsigned char _a[], int sol[]) +{ struct dsa _dsa, *dsa = &_dsa; + int i, j, p, max_wt, max_nwt, wth, *used, *nwt, *pos; + glp_long timer; + n = _n; + wt = &w[1]; + a = _a; + record = 0; + rec_level = 0; + rec = &sol[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(dsa, 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(dsa, p, j))) + nwt[j] -= wt[p]; + } + /* main routine */ + wth = 0; + for (i = 0; i < n; i++) + { wth += wt[pos[i]]; + sub(dsa, i, pos, 0, 0, wth); + clique[pos[i]] = record; +#if 0 + if (utime() >= timer + 5.0) +#else + if (xdifftime(xtime(), timer) >= 5.0 - 0.001) +#endif + { /* print current record and reset timer */ + xprintf("level = %d (%d); best = %d\n", i+1, n, record); +#if 0 + timer = utime(); +#else + timer = xtime(); +#endif + } + } + xfree(clique); + xfree(set); + xfree(used); + xfree(nwt); + xfree(pos); + /* return the solution found */ + for (i = 1; i <= rec_level; i++) sol[i]++; + return rec_level; +} + +#undef n +#undef wt +#undef a +#undef record +#undef rec_level +#undef rec +#undef clique +#undef set + +/*---------------------------------------------------------------------- +-- lpx_clique_cut - generate cluque cut. +-- +-- SYNOPSIS +-- +-- #include "glplpx.h" +-- int lpx_clique_cut(LPX *lp, void *cog, int ind[], double val[]); +-- +-- DESCRIPTION +-- +-- The routine lpx_clique_cut generates a clique cut using the conflict +-- graph specified by the parameter cog. +-- +-- If a violated clique cut has been found, it has the following form: +-- +-- sum{j in J} a[j]*x[j] <= b. +-- +-- Variable indices j in J are stored in elements ind[1], ..., ind[len] +-- while corresponding constraint coefficients are stored in elements +-- val[1], ..., val[len], where len is returned on exit. The right-hand +-- side b is stored in element val[0]. +-- +-- RETURNS +-- +-- If the cutting plane has been successfully generated, the routine +-- returns 1 <= len <= n, which is the number of non-zero coefficients +-- in the inequality constraint. Otherwise, the routine returns zero. */ + +static int lpx_clique_cut(LPX *lp, void *_cog, int ind[], double val[]) +{ struct COG *cog = _cog; + int n = lpx_get_num_cols(lp); + int j, t, v, card, temp, len = 0, *w, *sol; + double x, sum, b, *vec; + /* allocate working arrays */ + w = xcalloc(1 + 2 * cog->nb, sizeof(int)); + sol = xcalloc(1 + 2 * cog->nb, sizeof(int)); + vec = xcalloc(1+n, sizeof(double)); + /* assign weights to vertices of the conflict graph */ + for (t = 1; t <= cog->nb; t++) + { j = cog->orig[t]; + x = lpx_get_col_prim(lp, j); + temp = (int)(100.0 * x + 0.5); + if (temp < 0) temp = 0; + if (temp > 100) temp = 100; + w[t] = temp; + w[cog->nb + t] = 100 - temp; + } + /* find a clique of maximum weight */ + card = wclique(2 * cog->nb, w, cog->a, sol); + /* compute the clique weight for unscaled values */ + sum = 0.0; + for ( t = 1; t <= card; t++) + { v = sol[t]; + xassert(1 <= v && v <= 2 * cog->nb); + if (v <= cog->nb) + { /* vertex v corresponds to binary variable x[j] */ + j = cog->orig[v]; + x = lpx_get_col_prim(lp, j); + sum += x; + } + else + { /* vertex v corresponds to the complement of x[j] */ + j = cog->orig[v - cog->nb]; + x = lpx_get_col_prim(lp, j); + sum += 1.0 - x; + } + } + /* if the sum of binary variables and their complements in the + clique greater than 1, the clique cut is violated */ + if (sum >= 1.01) + { /* construct the inquality */ + for (j = 1; j <= n; j++) vec[j] = 0; + b = 1.0; + for (t = 1; t <= card; t++) + { v = sol[t]; + if (v <= cog->nb) + { /* vertex v corresponds to binary variable x[j] */ + j = cog->orig[v]; + xassert(1 <= j && j <= n); + vec[j] += 1.0; + } + else + { /* vertex v corresponds to the complement of x[j] */ + j = cog->orig[v - cog->nb]; + xassert(1 <= j && j <= n); + vec[j] -= 1.0; + b -= 1.0; + } + } + xassert(len == 0); + for (j = 1; j <= n; j++) + { if (vec[j] != 0.0) + { len++; + ind[len] = j, val[len] = vec[j]; + } + } + ind[0] = 0, val[0] = b; + } + /* free working arrays */ + xfree(w); + xfree(sol); + xfree(vec); + /* return to the calling program */ + return len; +} + +/*---------------------------------------------------------------------- +-- lpx_delete_cog - delete the conflict graph. +-- +-- SYNOPSIS +-- +-- #include "glplpx.h" +-- void lpx_delete_cog(void *cog); +-- +-- DESCRIPTION +-- +-- The routine lpx_delete_cog deletes the conflict graph, which the +-- parameter cog points to, freeing all the memory allocated to this +-- object. */ + +static void lpx_delete_cog(void *_cog) +{ struct COG *cog = _cog; + xfree(cog->vert); + xfree(cog->orig); + xfree(cog->a); + xfree(cog); +} + +/**********************************************************************/ + +void *ios_clq_init(glp_tree *tree) +{ /* initialize clique cut generator */ + glp_prob *mip = tree->mip; + xassert(mip != NULL); + return lpx_create_cog(mip); +} + +/*********************************************************************** +* NAME +* +* ios_clq_gen - generate clique cuts +* +* SYNOPSIS +* +* #include "glpios.h" +* void ios_clq_gen(glp_tree *tree, void *gen); +* +* DESCRIPTION +* +* The routine ios_clq_gen generates clique cuts for the current point +* and adds them to the clique pool. */ + +void ios_clq_gen(glp_tree *tree, void *gen) +{ int n = lpx_get_num_cols(tree->mip); + int len, *ind; + double *val; + xassert(gen != NULL); + ind = xcalloc(1+n, sizeof(int)); + val = xcalloc(1+n, sizeof(double)); + len = lpx_clique_cut(tree->mip, gen, ind, val); + if (len > 0) + { /* xprintf("len = %d\n", len); */ + glp_ios_add_row(tree, NULL, GLP_RF_CLQ, 0, len, ind, val, + GLP_UP, val[0]); + } + xfree(ind); + xfree(val); + return; +} + +/**********************************************************************/ + +void ios_clq_term(void *gen) +{ /* terminate clique cut generator */ + xassert(gen != NULL); + lpx_delete_cog(gen); + return; +} + +/* eof */