glpios08.c

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

/* glpios08.c (clique cut generator) */

/***********************************************************************
*  This code is part of GLPK (GNU Linear Programming Kit).
*
*  Copyright (C) 2000,01,02,03,04,05,06,07,08,2009 Andrew Makhorin,
*  Department for Applied Informatics, Moscow Aviation Institute,
*  Moscow, Russia. All rights reserved. E-mail: <[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 <http://www.gnu.org/licenses/>.
***********************************************************************/

#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: