glpipp02.c

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

      {  col = aij->col;
         if (aij->val > 0.0 && col->lb == -DBL_MAX ||
             aij->val < 0.0 && col->ub == +DBL_MAX)
         {  if (c_min == NULL)
               c_min = col;
            else
            {  f_min = -DBL_MAX;
               break;
            }
         }
         else
            f_min += aij->val * (aij->val > 0.0 ? col->lb : col->ub);
      }
      /* compute implied upper bound of the row */
      c_max = NULL, f_max = 0.0;
      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
      {  col = aij->col;
         if (aij->val > 0.0 && col->ub == +DBL_MAX ||
             aij->val < 0.0 && col->lb == -DBL_MAX)
         {  if (c_max == NULL)
               c_max = col;
            else
            {  f_max = +DBL_MAX;
               break;
            }
         }
         else
            f_max += aij->val * (aij->val > 0.0 ? col->ub : col->lb);
      }
      /* process all columns in the row */
      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
      {  col = aij->col;
         /* compute implied lower bound of the row minus a[j] * x[j] */
         if (f_min == -DBL_MAX)
            ff_min = -DBL_MAX;
         else if (c_min == NULL)
            ff_min = f_min -
               aij->val * (aij->val > 0.0 ? col->lb : col->ub);
         else if (c_min == col)
            ff_min = f_min;
         else
            ff_min = -DBL_MAX;
         /* compute implied upper bound of the row minus a[j] * x[j] */
         if (f_max == +DBL_MAX)
            ff_max = +DBL_MAX;
         else if (c_max == NULL)
            ff_max = f_max -
               aij->val * (aij->val > 0.0 ? col->ub : col->lb);
         else if (c_max == col)
            ff_max = f_max;
         else
            ff_max = +DBL_MAX;
         /* compute implied lower and upper bounds of x[j] */
#if 1
         /* do not use aij if it has small magnitude to prevent wrong
            implied bounds; for example, 1e-15 * x1 >= x2 + x3, where
            x1 >= -10, x2, x3 >= 0, would lead to wrong conclusion that
            x1 >= 0 */
         if (fabs(aij->val) < 1e-6)
            lb = -DBL_MAX, ub = +DBL_MAX; else
#endif
         if (aij->val > 0.0)
         {  if (row->lb == -DBL_MAX || ff_max == +DBL_MAX)
               lb = -DBL_MAX;
            else
               lb = (row->lb - ff_max) / aij->val;
            if (row->ub == +DBL_MAX || ff_min == -DBL_MAX)
               ub = +DBL_MAX;
            else
               ub = (row->ub - ff_min) / aij->val;
         }
         else
         {  if (row->ub == +DBL_MAX || ff_min == -DBL_MAX)
               lb = -DBL_MAX;
            else
               lb = (row->ub - ff_min) / aij->val;
            if (row->lb == -DBL_MAX || ff_max == +DBL_MAX)
               ub = +DBL_MAX;
            else
               ub = (row->lb - ff_max) / aij->val;
         }
         /* to prevent infinite reducing we replace current bounds of
            x[j] by its implied bounds only if there is a significant
            change (not less than 10%) in the bounds */
         flag = 0;
         /* estimate significance in lower bound */
         if (lb != -DBL_MAX)
         {  if (col->i_flag)
               delta = 0.001; /* even tiny change is significant */
            else
               delta = 0.10 * (1.0 + fabs(lb));
            if (lb - delta >= col->lb) flag = 1;
         }
         /* estimate significance in upper bound */
         if (ub != +DBL_MAX)
         {  if (col->i_flag)
               delta = 0.001; /* even tiny change is significant */
            else
               delta = 0.10 * (1.0 + fabs(ub));
            if (ub + delta <= col->ub) flag = 1;
         }
         /* if the change is significant, perform replacing */
#if 0
         if (flag)
#else
         if (flag && lb < +1e15 && ub > -1e15)
#endif
         {  switch(ipp_tight_bnds(ipp, col, lb, ub))
            {  case 0:
#if 0
                  /* bounds remain unchanged; can never be */
                  xassert(ipp != ipp);
#else
                  /* can be if lb or ub is large and x[j] is integer */
                  break;
#endif
               case 1:
                  /* bounds have been changed */
                  /* activate x[j] for further processing */
                  ipp_enque_col(ipp, col);
                  break;
               case 2:
                  /* new bounds are primal infeasible */
                  return 1;
               default:
                  xassert(ipp != ipp);
            }
         }
      }
      return 0;
}

/*----------------------------------------------------------------------
-- ipp_reduce_bnds - reduce column bounds.
--
-- SYNOPSIS
--
-- #include "glpipp.h"
-- int ipp_reduce_bnds(IPP *ipp);
--
-- DESCRIPTION
--
-- The routine ipp_reduce_bnds tries to reduce column bounds replacing
-- them by implied bounds.
--
-- Since changes in bounds of one column may involve changes in bounds
-- of other columns, the routine repeats reducing passes while there is
-- a significant change in bounds of at least one column.
--
-- RETURNS
--
-- 0 - no primal infeasibility is detected
-- 1 - implied bounds of some column are primal infeasible */

int ipp_reduce_bnds(IPP *ipp)
{     IPPROW *row;
      IPPCOL *col;
      IPPAIJ *aij;
      int pass = 0, total = 0, count;
      /* activate all rows */
      for (row = ipp->row_ptr; row != NULL; row = row->next)
         ipp_enque_row(ipp, row);
      /* deactivate all columns */
      for (col = ipp->col_ptr; col != NULL; col = col->next)
         ipp_deque_col(ipp, col);
loop: /* start a next pass for all active rows */
      pass++;
      while (ipp->row_que != NULL)
      {  /* select an active row */
         row = ipp->row_que;
         /* and remove it from the queue */
         ipp_deque_row(ipp, row);
         /* use the row to reduce column bounds */
         if (reduce_bounds(ipp, row)) return 1;
      }
      /* now all rows are inactive while all columns whose bounds were
         reduced are active */
      count = 0;
      while (ipp->col_que != NULL)
      {  count++;
         /* select an active column */
         col = ipp->col_que;
         /* and remove it from the queue */
         ipp_deque_col(ipp, col);
         /* activate corresponding rows for next pass */
         for (aij = col->ptr; aij != NULL; aij = aij->c_next)
         {  row = aij->row;
            ipp_enque_row(ipp, row);
         }
      }
      total += count;
      if (count > 0) goto loop;
      xprintf(
         "ipp_reduce_bnds: %d pass(es) made, %d bound(s) reduced\n",
         pass, total);
      return 0;
}

/*----------------------------------------------------------------------
-- SHIFTED COLUMN
--
-- Let column q has a finite non-zero lower bound:
--
--    l[q] <= x[q] <= u[q].                                          (1)
--
-- Then it can be replaced by:
--
--    x[q] = x'[q] + l[q],                                           (2)
--
-- where:
--
--    0 <= x'[q] <= u[q] - l[q]                                      (3)
--
-- is a non-negative column.
--
-- RECOVERING
--
-- Value of x[q] is recovered with formula (2) assuming that value of
-- x'[q] is already recovered. */

struct shift_col
{     /* shifted column */
      int q;
      /* reference number of a shifted column */
      double s;
      /* shifting value = original non-zero lower bound */
};

void ipp_shift_col(IPP *ipp, IPPCOL *col)
{     /* process shifted column */
      struct shift_col *info;
      IPPROW *row;
      IPPAIJ *aij;
      double temp;
      /* the column must have finite non-zero lower bound */
      xassert(col->lb != -DBL_MAX && col->lb != 0.0);
      /* create transformation queue entry */
      info = ipp_append_tqe(ipp, IPP_SHIFT_COL, sizeof(*info));
      info->q = col->j;
      info->s = col->lb;
      /* update bounds of corresponding rows */
      for (aij = col->ptr; aij != NULL; aij = aij->c_next)
      {  row = aij->row;
         temp = aij->val * info->s;
         if (row->lb == row->ub)
         {  row->lb -= temp;
            row->ub = row->lb;
         }
         else
         {  if (row->lb != -DBL_MAX) row->lb -= temp;
            if (row->ub != +DBL_MAX) row->ub -= temp;
         }
      }
      /* update the constant term of the objective function */
      ipp->c0 += col->c * info->s;
      /* set shifted bounds of the column */
      col->lb = 0.0;
      if (col->ub != +DBL_MAX) col->ub -= info->s;
      return;
}

void ipp_shift_col_r(IPP *ipp, void *_info)
{     /* recover shifted column */
      struct shift_col *info = _info;
      xassert(1 <= info->q && info->q <= ipp->ncols);
      xassert(ipp->col_stat[info->q] == 1);
      ipp->col_mipx[info->q] += info->s;
      return;
}

/*----------------------------------------------------------------------
-- NON-BINARY COLUMN
--
-- Let column q is integer non-negative:
--
--    0 <= x[q] <= u[q]                                              (1)
--
-- where u[q] > 1, i.e. the column is non-binary. In this case it can
-- be replaced by:
--
--    x[q] = 2^0 * z[0] + 2^1 * z[1] + ... + 2^(t-1) * z[t-1],       (2)
--
-- where z[0], z[1], ..., z[t-1] are binary variables, and t is minimal
-- number of bits sufficient to represent u[q].
--
-- To keep the upper bound of x[q] the following additional constraint
-- should be added to the problem:
--
--    2^0 * z[0] + 2^1 * z[1] + ... + 2^(t-1) * z[t-1] <= u[q].      (3)
--
-- However, if u[q] = 2^t - 1, the constraint (3) is redundant and may
-- be omitted.
--
-- RECOVERING
--
-- Value of x[q] is recovered with formula (2) assuming that values of
-- z[0], z[1], ..., z[t-1] are already recovered. */

struct nonbin_col
{     /* non-binary column */
      int q;
      /* reference number of a non-binary column */
      IPPLFE *ptr;
      /* 2^0 * z[0] + 2^1 * z[1] + ... + 2^(t-1) * z[t-1] */
};

int ipp_nonbin_col(IPP *ipp, IPPCOL *col)
{     /* process non-binary column */
      struct nonbin_col *info;
      IPPROW *row;
      IPPCOL *bin;
      IPPAIJ *aij;
      IPPLFE *lfe;
      int u, t, two_t, k, two_k;
      /* the column must be integral */
      xassert(col->i_flag);
      /* its lower bound must be zero */
      xassert(col->lb == 0.0);
      /* its upper bound must be greater than one */
      xassert(col->ub >= 2.0);
      /* and must be not greater than 2^15-1 = 32767 (implementation
         restriction) */
      xassert(col->ub <= 32767.0);
      /* create transformation queue entry */
      info = ipp_append_tqe(ipp, IPP_NONBIN_COL, sizeof(*info));
      info->q = col->j;
      info->ptr = NULL;
      /* determine t, minimal number of bits sufficient to represent
         the upper bound */
      u = (int)col->ub;
      xassert(col->ub == (double)u);
      for (t = 2, two_t = 4; t <= 15; t++, two_t += two_t)
         if (u <= two_t - 1) break;
      xassert(t <= 15);
      /* create additional constraint (3), if necessary */
      if (u <= two_t - 2)
         row = ipp_add_row(ipp, -DBL_MAX, (double)u);
      /* create binary columns z[0], z[1], ..., z[t-1] */
      for (k = 0, two_k = 1; k < t; k++, two_k += two_k)
      {  bin = ipp_add_col(ipp, 1, 0.0, 1.0, 0.0);
         lfe = dmp_get_atom(ipp->tqe_pool, sizeof(IPPLFE));
         lfe->ref = bin->j;
         lfe->val = (double)two_k;
         lfe->next = info->ptr;
         info->ptr = lfe;