glpios08.c

著名的大规模线性规划求解器源码GLPK.C语言版本,可以修剪.内有详细帮助文档.

                     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 <http://www.tcs.hut.fi/~pat/wclique.html> 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;
      xlong_t 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 */