glpios07.c

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

/*----------------------------------------------------------------------
-- lpx_cover_cut - generate mixed cover cut.
--
-- SYNOPSIS
--
-- #include "glplpx.h"
-- int lpx_cover_cut(LPX *lp, int len, int ind[], double val[],
--    double work[]);
--
-- DESCRIPTION
--
-- The routine lpx_cover_cut generates a mixed cover cut for a given
-- row of the MIP problem.
--
-- The given row of the MIP problem should be explicitly specified in
-- the form:
--
--    sum{j in J} a[j]*x[j] <= b.                                    (1)
--
-- On entry indices (ordinal numbers) of structural variables, which
-- have non-zero constraint coefficients, should be placed in locations
-- ind[1], ..., ind[len], and corresponding constraint coefficients
-- should be placed in locations val[1], ..., val[len]. The right-hand
-- side b should be stored in location val[0].
--
-- The working array work should have at least nb locations, where nb
-- is the number of binary variables in (1).
--
-- The routine generates a mixed cover cut in the same form as (1) and
-- stores the cut coefficients and right-hand side in the same way as
-- just described above.
--
-- 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_cover_cut(LPX *lp, int len, int ind[], double val[],
      double work[])
{     int cov[1+4], j, k, nb, newlen, r;
      double f_min, f_max, alfa, beta, u, *x = work, y;
      /* substitute and remove fixed variables */
      newlen = 0;
      for (k = 1; k <= len; k++)
      {  j = ind[k];
         if (lpx_get_col_type(lp, j) == LPX_FX)
            val[0] -= val[k] * lpx_get_col_lb(lp, j);
         else
         {  newlen++;
            ind[newlen] = ind[k];
            val[newlen] = val[k];
         }
      }
      len = newlen;
      /* move binary variables to the beginning of the list so that
         elements 1, 2, ..., nb correspond to binary variables, and
         elements nb+1, nb+2, ..., len correspond to rest variables */
      nb = 0;
      for (k = 1; k <= len; k++)
      {  j = ind[k];
         if (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)
         {  /* binary variable */
            int ind_k;
            double val_k;
            nb++;
            ind_k = ind[nb], val_k = val[nb];
            ind[nb] = ind[k], val[nb] = val[k];
            ind[k] = ind_k, val[k] = val_k;
         }
      }
      /* now the specified row has the form:
         sum a[j]*x[j] + sum a[j]*y[j] <= b,
         where x[j] are binary variables, y[j] are rest variables */
      /* at least two binary variables are needed */
      if (nb < 2) return 0;
      /* compute implied lower and upper bounds for sum a[j]*y[j] */
      f_min = f_max = 0.0;
      for (k = nb+1; k <= len; k++)
      {  j = ind[k];
         /* both bounds must be finite */
         if (lpx_get_col_type(lp, j) != LPX_DB) return 0;
         if (val[k] > 0.0)
         {  f_min += val[k] * lpx_get_col_lb(lp, j);
            f_max += val[k] * lpx_get_col_ub(lp, j);
         }
         else
         {  f_min += val[k] * lpx_get_col_ub(lp, j);
            f_max += val[k] * lpx_get_col_lb(lp, j);
         }
      }
      /* sum a[j]*x[j] + sum a[j]*y[j] <= b ===>
         sum a[j]*x[j] + (sum a[j]*y[j] - f_min) <= b - f_min ===>
         sum a[j]*x[j] + y <= b - f_min,
         where y = sum a[j]*y[j] - f_min;
         note that 0 <= y <= u, u = f_max - f_min */
      /* determine upper bound of y */
      u = f_max - f_min;
      /* determine value of y at the current point */
      y = 0.0;
      for (k = nb+1; k <= len; k++)
      {  j = ind[k];
         y += val[k] * lpx_get_col_prim(lp, j);
      }
      y -= f_min;
      if (y < 0.0) y = 0.0;
      if (y > u) y = u;
      /* modify the right-hand side b */
      val[0] -= f_min;
      /* now the transformed row has the form:
         sum a[j]*x[j] + y <= b, where 0 <= y <= u */
      /* determine values of x[j] at the current point */
      for (k = 1; k <= nb; k++)
      {  j = ind[k];
         x[k] = lpx_get_col_prim(lp, j);
         if (x[k] < 0.0) x[k] = 0.0;
         if (x[k] > 1.0) x[k] = 1.0;
      }
      /* if a[j] < 0, replace x[j] by its complement 1 - x'[j] */
      for (k = 1; k <= nb; k++)
      {  if (val[k] < 0.0)
         {  ind[k] = - ind[k];
            val[k] = - val[k];
            val[0] += val[k];
            x[k] = 1.0 - x[k];
         }
      }
      /* try to generate a mixed cover cut for the transformed row */
      r = cover(nb, val, val[0], u, x, y, cov, &alfa, &beta);
      if (r == 0) return 0;
      xassert(2 <= r && r <= 4);
      /* now the cut is in the form:
         sum{j in C} x[j] + alfa * y <= beta */
      /* store the right-hand side beta */
      ind[0] = 0, val[0] = beta;
      /* restore the original ordinal numbers of x[j] */
      for (j = 1; j <= r; j++) cov[j] = ind[cov[j]];
      /* store cut coefficients at binary variables complementing back
         the variables having negative row coefficients */
      xassert(r <= nb);
      for (k = 1; k <= r; k++)
      {  if (cov[k] > 0)
         {  ind[k] = +cov[k];
            val[k] = +1.0;
         }
         else
         {  ind[k] = -cov[k];
            val[k] = -1.0;
            val[0] -= 1.0;
         }
      }
      /* substitute y = sum a[j]*y[j] - f_min */
      for (k = nb+1; k <= len; k++)
      {  r++;
         ind[r] = ind[k];
         val[r] = alfa * val[k];
      }
      val[0] += alfa * f_min;
      xassert(r <= len);
      len = r;
      return len;
}

/*----------------------------------------------------------------------
-- lpx_eval_row - compute explictily specified row.
--
-- SYNOPSIS
--
-- #include "glplpx.h"
-- double lpx_eval_row(LPX *lp, int len, int ind[], double val[]);
--
-- DESCRIPTION
--
-- The routine lpx_eval_row computes the primal value of an explicitly
-- specified row using current values of structural variables.
--
-- The explicitly specified row may be thought as a linear form:
--
--    y = a[1]*x[m+1] + a[2]*x[m+2] + ... + a[n]*x[m+n],
--
-- where y is an auxiliary variable for this row, a[j] are coefficients
-- of the linear form, x[m+j] are structural variables.
--
-- On entry column indices and numerical values of non-zero elements of
-- the row should be stored in locations ind[1], ..., ind[len] and
-- val[1], ..., val[len], where len is the number of non-zero elements.
-- The array ind and val are not changed on exit.
--
-- RETURNS
--
-- The routine returns a computed value of y, the auxiliary variable of
-- the specified row. */

static double lpx_eval_row(LPX *lp, int len, int ind[], double val[])
{     int n = lpx_get_num_cols(lp);
      int j, k;
      double sum = 0.0;
      if (len < 0)
         xerror("lpx_eval_row: len = %d; invalid row length\n", len);
      for (k = 1; k <= len; k++)
      {  j = ind[k];
         if (!(1 <= j && j <= n))
            xerror("lpx_eval_row: j = %d; column number out of range\n",
               j);
         sum += val[k] * lpx_get_col_prim(lp, j);
      }
      return sum;
}

/***********************************************************************
*  NAME
*
*  ios_cov_gen - generate mixed cover cuts
*
*  SYNOPSIS
*
*  #include "glpios.h"
*  void ios_cov_gen(glp_tree *tree);
*
*  DESCRIPTION
*
*  The routine ios_cov_gen generates mixed cover cuts for the current
*  point and adds them to the cut pool. */

void ios_cov_gen(glp_tree *tree)
{     glp_prob *prob = tree->mip;
      int m = lpx_get_num_rows(prob);
      int n = lpx_get_num_cols(prob);
      int i, k, type, kase, len, *ind;
      double r, *val, *work;
      xassert(lpx_get_status(prob) == LPX_OPT);
      /* allocate working arrays */
      ind = xcalloc(1+n, sizeof(int));
      val = xcalloc(1+n, sizeof(double));
      work = xcalloc(1+n, sizeof(double));
      /* look through all rows */
      for (i = 1; i <= m; i++)
      for (kase = 1; kase <= 2; kase++)
      {  type = lpx_get_row_type(prob, i);
         if (kase == 1)
         {  /* consider rows of '<=' type */
            if (!(type == LPX_UP || type == LPX_DB)) continue;
            len = lpx_get_mat_row(prob, i, ind, val);
            val[0] = lpx_get_row_ub(prob, i);
         }
         else
         {  /* consider rows of '>=' type */
            if (!(type == LPX_LO || type == LPX_DB)) continue;
            len = lpx_get_mat_row(prob, i, ind, val);
            for (k = 1; k <= len; k++) val[k] = - val[k];
            val[0] = - lpx_get_row_lb(prob, i);
         }
         /* generate mixed cover cut:
            sum{j in J} a[j] * x[j] <= b */
         len = lpx_cover_cut(prob, len, ind, val, work);
         if (len == 0) continue;
         /* at the current point the cut inequality is violated, i.e.
            sum{j in J} a[j] * x[j] - b > 0 */
         r = lpx_eval_row(prob, len, ind, val) - val[0];
         if (r < 1e-3) continue;
         /* add the cut to the cut pool */
         glp_ios_add_row(tree, NULL, GLP_RF_COV, 0, len, ind, val,
            GLP_UP, val[0]);
      }
      /* free working arrays */
      xfree(ind);
      xfree(val);
      xfree(work);
      return;
}

/* eof */