glplpp02.c

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

            }
            else
            {  /* it is its own upper bound, so the row is inactive */
               lpp->row_stat[info->p] = LPX_BS;
               lpp->row_prim[info->p] =
                  info->apq * lpp->col_prim[info->q];
               lpp->row_dual[info->p] = 0.0;
            }
            break;
         case LPX_NF:
            /* the column cannot be free, since the row is always has
               at least one finite bound; see (5) */
            xassert(lpp != lpp);
            /* no break */
         case LPX_NS:
            /* the column is non-basic and fixed; it cannot be fixed
               before tightening its bounds, because in that case this
               transformation entry is not created; therefore we need
               to consider the column as double-bounded whose lower and
               upper bounds are equal to each other, and choose a bound
               using its dual value lambda[q] */
            if (lpp->col_dual[info->q] >= 0.0)
            {  /* lambda[q] >= 0; set the column on lower bound */
               lpp->col_stat[info->q] = LPX_NL;
               goto nl;
            }
            else
            {  /* lambda[q] < 0; set the column on upper bound */
               lpp->col_stat[info->q] = LPX_NU;
               goto nu;
            }
            /* no break */
         default:
            xassert(0);
      }
      return;
}

/*----------------------------------------------------------------------
-- COLUMN SINGLETON (IMPLIED SLACK VARIABLE)
--
-- *Processing*
--
-- Let column q have the only non-zero constraint coefficient in row p,
-- which is an equality constraint:
--
--    y[p] =  sum a[p,j] * x[j] + a[p,q] * x[q] = b[p],              (1)
--           j!=q
--
--    l[q] <= x[q] <= u[q].                                          (2)
--
-- The term a[p,q] * x[q] can be considered as a slack variable of the
-- row p, that allows removing the column q from the problem.
--
-- From (1) it follows that:
--
--    sum a[p,j] * x[j] = b[p] - a[p,q] * x[q].                      (3)
--   j!=q
--
-- So, (1) can be replaced by the following equivalent constraint:
--
--    l' <= y' =  sum a[p,j] * x[j] <= u',                           (4)
--               j!=q
--
-- where y' is an auxiliary variable of this new constraint, and
--
--    if a[p,q] > 0:   l' = b[p] - a[p,q] * u[q],                    (5)
--                     u' = b[p] - a[p,q] * l[q],
--
--    if a[p,q] < 0:   l' = b[p] - a[p,q] * l[q],                    (6)
--                     u' = b[p] - a[p,q] * u[q].
--
-- On removing x[q] from the problem it also should be substituted in
-- the objective function. From (3) it follows that:
--
--    x[q] = - sum (a[p,j] / a[p,q]) * x[j] + b[p] / a[p,q],         (7)
--            j!=q
--
-- and substituting x[q] in the objective function gives:
--
--    Z = sum c[j] * x[j] + c[0] =
--         j
--
--      = sum c[j] * x[j] + c[q] * x[q] + c[0] =                     (8)
--       j!=q
--
--      = sum c'[j] * x[j] + c'[0],
--       j!=q
--
-- where
--
--    c'[j] = c[j] - c[q] * (a[p,j] / a[p,q]),                       (9)
--
--    c'[0] = c[0] + c[q] * (b[p] / a[p,q]).                        (10)
--
-- *Recovering*
--
-- On entry to the recovering routine the row (4) that corresponds to
-- the row p is already recovered, i.e. its status, primal value y' and
-- dual value pi' are known.
--
-- From (3) and (4) it follows that
--
--    y' = b[p] - a[p,q] * x[q].                                    (11)
--
-- Therefore the status of the column q is determined as follows:
--
-- if y' is basic, x[q] is basic;
--
-- if y' is non-basic on lower bound, x[q] is non-basic on upper (if
-- a[p,q] > 0) or lower (if a[p,q] < 0) bound;
--
-- if y' is non-basic on upper bound, x[q] is non-basic on lower (if
-- a[p,q] > 0) or upper (of a[p,q] < 0) bound.
--
-- The primal and dual values of the column q are computed as follows:
--
--    x[q] = (b[p] - y') / a[p,q],                                  (12)
--
--    lambda[q] = - a[p,q] * pi',                                   (13)
--
-- where y' and pi' are primal and dual values of the row (4), resp.
--
-- Being an equality constraint the row p is active.
--
-- From (1) it follows that
--
--    y[p] = y' + a[p,q] * x[q] = b[p],                             (14)
--
-- which is the formula to compute the primal value of the row p.
--
-- The dual equality constraint for the row p is:
--
--    a[p,q] * pi[p] + lambda[q] = c[q],                            (15)
--
-- therefore
--
--    pi[p] = (c[q] - lambda[q]) / a[p,q],                          (16)
--
-- which is the formula to compute the dual value of the row p. */

typedef struct
{     /* column singleton (implied slack variable) */
      int p;
      /* row reference number */
      int q;
      /* column reference number */
      double rhs;
      /* b[p], right-hand side of the row */
      double c;
      /* c[q], objective coefficient of the column */
      double apq;
      /* a[p,q], constraint coefficient */
} COL_SNGTON1;

static void process_col_sngton1(LPP *lpp, LPPCOL *col)
{     COL_SNGTON1 *info;
      LPPROW *row;
      LPPAIJ *aij;
      double lb, ub, eps;
      /* the column must have the only non-zero constraint coefficient
         in the corresponding row */
      xassert(col->ptr != NULL && col->ptr->c_next == NULL);
      /* the corresponding row must be equality constraint */
      aij = col->ptr;
      row = aij->row;
      xassert(row->lb == row->ub);
      /* if the column is fixed, it can be just substituted and removed
         from the problem */
      if (col->lb == col->ub)
      {  process_fixed_col(lpp, col);
         goto done;
      }
      /* create transformation queue entry */
      info = lpp_append_tqe(lpp, LPP_COL_SNGTON1, sizeof(COL_SNGTON1));
      info->p = row->i;
      info->q = col->j;
      info->rhs = row->lb;
      info->c = col->c;
      info->apq = aij->val;
      /* compute new bounds l' and u' of the row; see (4)--(6) */
      if (info->apq > 0.0)
      {  lb = (col->ub ==
            +DBL_MAX ? -DBL_MAX : info->rhs - info->apq * col->ub);
         ub = (col->lb ==
            -DBL_MAX ? +DBL_MAX : info->rhs - info->apq * col->lb);
      }
      else
      {  lb = (col->lb ==
            -DBL_MAX ? -DBL_MAX : info->rhs - info->apq * col->lb);
         ub = (col->ub ==
            +DBL_MAX ? +DBL_MAX : info->rhs - info->apq * col->ub);
      }
      row->lb = lb;
      row->ub = ub;
      /* if new bounds of the row are close to each other, the row can
         be treated as equality constraint */
      if (row->lb != -DBL_MAX && row->ub != +DBL_MAX)
      {  eps = 1e-7 * (1.0 + fabs(row->lb));
         if (fabs(row->lb - row->ub) <= eps)
         {  if (fabs(row->lb) <= fabs(row->ub))
               row->ub = row->lb;
            else
               row->lb = row->ub;
         }
      }
      /* remove the column from the problem */
      lpp_remove_col(lpp, col);
      /* substitute x[q] into the objective function; see (7)--(10) */
      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
         aij->col->c -= info->c * (aij->val / info->apq);
      lpp->c0 += info->c * (info->rhs / info->apq);
done: return;
}

static void recover_col_sngton1(LPP *lpp, COL_SNGTON1 *info)
{     /* the (modified) row is already recovered */
      xassert(1 <= info->p && info->p <= lpp->nrows);
      xassert(lpp->row_stat[info->p] != 0);
      /* while the column is not recovered yet */
      xassert(1 <= info->q && info->q <= lpp->ncols);
      xassert(lpp->col_stat[info->q] == 0);
      /* recover the column status */
      switch (lpp->row_stat[info->p])
      {  case LPX_BS:
            /* y' is basic */
            lpp->col_stat[info->q] = LPX_BS;
            break;
         case LPX_NL:
nl:         /* y' is non-basic on lower bound */
            lpp->col_stat[info->q] =
               (info->apq > 0.0 ? LPX_NU : LPX_NL);
            break;
         case LPX_NU:
nu:         /* y' is non-basic on upper bound */
            lpp->col_stat[info->q] =
               (info->apq > 0.0 ? LPX_NL : LPX_NU);
            break;
         case LPX_NF:
            /* y' is non-basic free */
            xassert(lpp != lpp /* this case is not tested yet */);
            lpp->col_stat[info->q] = LPX_NF;
            break;
         case LPX_NS:
            /* y' is non-basic fixed; this can only mean the row was
               turned to equality constraint, because if the column is
               of fixed type, this transformation entry is not created;
               thus, we need to consider the row as ranged, whose lower
               and upper bounds are equal to each other, and choose an
               appropriate bound using the dual value pi[p] */
#if 0
            xassert(lpp != lpp /* this case is not tested yet */);
#endif
            if (lpp->row_dual[info->p] >= 0.0) goto nl; else goto nu;
            /* no break */
         default:
            xassert(lpp != lpp);
      }
      /* recover primal value x[q] using the formula (12) */
      lpp->col_prim[info->q] =
         (info->rhs - lpp->row_prim[info->p]) / info->apq;
      /* recover dual value lambda[q] using the formula (13) */
      lpp->col_dual[info->q] = - info->apq * lpp->row_dual[info->p];
      /* the row is active equality constraint */
      lpp->row_stat[info->p] = LPX_NS;
      /* recover primal value y[p] using the formula (14) */
      lpp->row_prim[info->p] = info->rhs;
      /* recover dual value pi[p] using the formula (15) */
      lpp->row_dual[info->p] =
         (info->c - lpp->col_dual[info->q]) / info->apq;
      return;
}

/*----------------------------------------------------------------------
-- COLUMN SINGLETON (IMPLIED FREE VARIABLE)
--
-- *Processing*
--
-- Let column q have the only non-zero constraint coefficient in row p,
-- which is an inequality constraint:
--
--    l[p] <= y[p] =  sum a[p,j] * x[j] + a[p,q] * x[q] <= u[p],     (1)
--                   j!=q
--
--    l[q] <= x[q] <= u[q].                                          (2)
--
-- Resolving (1) through x[q] we have:
--
--    x[q] = a'[p,q] * y[p] + sum a'[p,j] * x[j],                    (3)
--                           j!=q
--
-- where
--
--    a'[p,q] = 1 / a[p,q],                                          (4)
--
--    a'[p,j] = - a[p,j] / a[p,q],                                   (5)
--
-- The formula (3) allows computing implied lower and upper bounds of
-- x[q] using explicit bounds of y[p] and x[j]. Let these impled bounds
-- of x[q] are l' and u'. If l' >= l[q] and u' <= u[q], own bounds of
-- x[q] are never reached, in which case x[q] can be considered as an
-- implied free variable.
--
-- Let x[q] be a free variable. Then it is always basic, and its dual
-- value lambda[q] is zero. If x[q] is basic, y[p] must be non-basic
-- (otherwise there were two linearly dependent columns in the basic
-- matrix). The dual value pi[p] can be found from the dual equality
-- constraint for the column q:
--
--    a[p,q] * pi[p] + lambda[q] = c[q],                             (6)
--
-- that gives:
--
--    pi[p] = (c[q] - lambda[q]) / a[p,q] = c[q] / a[p,q].           (7)
--
-- If pi[p] > 0, the row p must be active on its lower bound l[p], and
-- if the row p has no lower bound, it is dual infeasible. Analogously,
-- if pi[p] < 0, the row p must be active on its upper bound u[p], and
-- if the row p has no upper bound, it is dual infeasible. Thus, in both
-- cases the row p becomes an equality constraint, whose right-hand side
-- is l[p] or u[p] depending on the sign of pi[p], while the column q
-- becomes an implied slack variable. If pi[p] = 0 (i.e. if c[q] = 0),
-- the row p can be fixed at any bound.
--
-- *Recovering*
--
-- The only thing needed on recovering is to properly set the status of
-- the row p, because on entry to the recovering routine this row is an
-- active equality constraint while actually it is an active inequality
-- constraint. */

typedef struct
{     /* column singleton (implied free variable) */
      int p;
      /* row reference number */
      int q;
      /* column reference number */
      int stat;
      /* row status:
         LPX_NL - active constraint on lower bound
         LPX_NU - active constraint on upper bound */
} COL_SNGTON2;