glplpp02.c

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

static int process_col_sngton2(LPP *lpp, LPPCOL *col)
{     COL_SNGTON2 *info;
      LPPROW *row;
      LPPAIJ *apq, *aij;
      double lb, ub, eps, temp;
      /* 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 inequality constraint */
      apq = col->ptr;
      row = apq->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;
      }
      /* if the row is free, it can be just removed from the problem */
      if (row->lb == -DBL_MAX && row->ub == +DBL_MAX)
      {  process_free_row(lpp, row);
         goto done;
      }
      /* compute l', an implied lower bound of x[q]; see (3)--(5) */
      temp = 1.0 / apq->val;
      if (temp > 0.0)
         lb = (row->lb == -DBL_MAX ? -DBL_MAX : temp * row->lb);
      else
         lb = (row->ub == +DBL_MAX ? -DBL_MAX : temp * row->ub);
      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
      {  if (lb == -DBL_MAX) break;
         if (aij == apq) continue;
         temp = - aij->val / apq->val;
         if (temp > 0.0)
         {  if (aij->col->lb == -DBL_MAX)
               lb = -DBL_MAX;
            else
               lb += temp * aij->col->lb;
         }
         else
         {  if (aij->col->ub == +DBL_MAX)
               lb = -DBL_MAX;
            else
               lb += temp * aij->col->ub;
         }
      }
      /* compute u', an implied upper bound of x[q]; see (3)--(5) */
      temp = 1.0 / apq->val;
      if (temp > 0.0)
         ub = (row->ub == +DBL_MAX ? +DBL_MAX : temp * row->ub);
      else
         ub = (row->lb == -DBL_MAX ? +DBL_MAX : temp * row->lb);
      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
      {  if (ub == +DBL_MAX) break;
         if (aij == apq) continue;
         temp = - aij->val / apq->val;
         if (temp > 0.0)
         {  if (aij->col->ub == +DBL_MAX)
               ub = +DBL_MAX;
            else
               ub += temp * aij->col->ub;
         }
         else
         {  if (aij->col->lb == -DBL_MAX)
               ub = +DBL_MAX;
            else
               ub += temp * aij->col->lb;
         }
      }
      /* check if x[q] can reach its own lower bound */
      if (col->lb != -DBL_MAX)
      {  eps = 1e-7 * (1.0 + fabs(col->lb));
         if (lb < col->lb - eps) goto done; /* yes, it can */
      }
      /* check if x[q] can reach its own upper bound */
      if (col->ub != +DBL_MAX)
      {  eps = 1e-7 * (1.0 * fabs(col->ub));
         if (ub > col->ub + eps) goto done; /* yes, it can */
      }
      /* create transformation queue entry */
      info = lpp_append_tqe(lpp, LPP_COL_SNGTON2, sizeof(COL_SNGTON2));
      info->p = row->i;
      info->q = col->j;
      info->stat = 0;
      /* x[q] is implied free variable */
      col->lb = -DBL_MAX, col->ub = +DBL_MAX;
      /* since x[q] is always basic, the row p must be active; compute
         its dual value pi[p]; see (7) */
      temp = col->c / apq->val;
      /* analyze the dual value pi[p] */
      if (temp > 0.0)
      {  /* pi[p] > 0; the row p must be active on its lower bound */
         if (row->lb == -DBL_MAX) return 1; /* dual infeasibility */
         info->stat = LPX_NL;
         row->ub = row->lb;
      }
      else if (temp < 0.0)
      {  /* pi[p] < 0; the row p must be active on its upper bound */
         if (row->ub == +DBL_MAX) return 1; /* dual infeasibility */
         info->stat = LPX_NU;
         row->lb = row->ub;
      }
      else
      {  /* pi[p] = 0; the row must be active on any bound */
         if (row->lb != -DBL_MAX)
         {  info->stat = LPX_NL;
            row->ub = row->lb;
         }
         else
         {  xassert(row->ub != +DBL_MAX);
            info->stat = LPX_NU;
            row->lb = row->ub;
         }
      }
      /* now x[q] can be considered as an implied slack variable */
      process_col_sngton1(lpp, col);
done: return 0;
}

static void recover_col_sngton2(LPP *lpp, COL_SNGTON2 *info)
{     /* the (modified) row is already recovered and must be an active
         equality constraint */
      xassert(1 <= info->p && info->p <= lpp->nrows);
      xassert(lpp->row_stat[info->p] == LPX_NS);
      /* the (modified) column is already recovered and must be a basic
         variable */
      xassert(1 <= info->q && info->q <= lpp->ncols);
      xassert(lpp->col_stat[info->q] == LPX_BS);
      /* set proper status of the row */
      lpp->row_stat[info->p] = info->stat;
      return;
}

/*----------------------------------------------------------------------
-- FORCING ROW
--
-- *Processing*
--
-- Let row p be of general kind:
--
--    l[p] <= y[p] = sum a[p,j] * x[j] <= u[p],                      (1)
--                    j
--
--    l[j] <= x[j] <= u[j].                                          (2)
--
-- If l[p] = u' or u[p] = l', where l' and u' are implied lower and
-- upper bounds of the row (see (3) and (4) in comments to the routine
-- analyze_row), the constraint forces the corresponding variables x[j]
-- to be set on their bounds in order to provide primal feasibility:
--
-- if l[p] = u' and a[p,j] > 0 or if u[p] = l' and a[p,j] < 0, x[j] can
-- only be set on its upper bound u[j], and
--
-- if l[p] = u' and a[p,j] < 0 or if u[p] = l' and a[p,j] > 0, x[j] can
-- only be set on its lower bound l[j].
--
-- Being set on their bounds the variables x[j] become fixed and can be
-- removed from the problem. At the same time the row becomes redundant
-- and therefore can also be removed from the problem.
--
-- *Recovering*
--
-- On entry to the recovering routine all x[j] are non-basic and fixed,
-- while the row p is still removed.
--
-- Since the corresponding basic solution is assumed to be optimal, the
-- following dual equality constraints are satisfied:
--
--    sum a[i,j] * pi[i] + lambda'[j] = c[j],                        (3)
--   i!=p
--
-- where lambda'[j] is a dual value of the column j on entry to the
-- recovering routine.
--
-- On recovering the row p an additional term that corresponds to this
-- row appears in the dual constraints (3):
--
--    sum a[i,j] * pi[i] + a[p,j] * pi[p] + lambda[j] = c[j],        (4)
--   i!=p
--
-- where pi[p] is recovered dual value of the row p, and lambda[j] is
-- recovered dual value of the column j.
--
-- From (3) and (4) it follows that:
--
--    lambda[j] = lambda'[j] - a[p,j] * pi[p].                       (5)
--
-- Now note that actually x[j] are non-fixed. Therefore, the following
-- dual feasibility condition must be satisfied for all x[j]:
--
--    lambda[j] >= 0, if x[j] was forced on its lower bound l[j],    (6)
--
--    lambda[j] <= 0, if x[j] was forced on its upper bound u[j].    (7)
--
-- Let the row p is non-active (i.e. y[p] is basic). Then pi[p] = 0 and
-- therefore lambda[j] = lambda'[j]. If the conditions (6)--(7) are
-- satisfied for all x[j], recovering is finished. Otherwise, we should
-- change (increase or decrease) pi[p] while at least one lambda[j] in
-- (5) has wrong sign. (Note that it is *always* possible to attain dual
-- feasibility by changing pi[p].) Once signs of all lambda[j] become
-- correct, there is always some lambda[q], which reaches zero last. In
-- this case the row p is active (i.e. y[p] is non-basic) with non-zero
-- dual value pi[p], all columns (except the column q) are non-basic
-- with dual values lambda[j], which are computed using the formula (5),
-- and the column q is basic with zero dual value lambda[q]. Should also
-- note that due to primal degeneracy changing pi[p] doesn't affect any
-- primal values y[p] and x[j]. */

typedef struct
{     /* forcing row */
      int p;
      /* reference number of a forcing row */
      int stat;
      /* status assigned to the row if it becomes active constraint:
         LPX_NS - equality constraint
         LPX_NL - inequality constraint on lower bound
         LPX_NU - inequality constraint on upper bound */
      double bnd;
      /* primal value of the row (one of its bounds) */
      LPPLFX *ptr;
      /* list of non-zero coefficients a[p,j] with additional flags of
         the corresponding variables x[j]:
         LPX_NL - x[j] is forced on its lower bound
         LPX_NU - x[j] is forced on its upper bound */
} FORCING_ROW;

static void process_forcing_row(LPP *lpp, LPPROW *row, int at)
{     FORCING_ROW *info;
      LPPCOL *col;
      LPPAIJ *aij, *next_aij;
      LPPLFX *lfx;
      /* if there are fixed columns in the row, they can be substituted
         and thus removed from the problem */
      for (aij = row->ptr; aij != NULL; aij = next_aij)
      {  col = aij->col;
         next_aij = aij->r_next;
         if (col->lb == col->ub) process_fixed_col(lpp, col);
      }
      /* the row may be empty; in this case it being redundant can just
         be removed from the problem */
      if (row->ptr == NULL)
      {  row->lb = -DBL_MAX, row->ub = +DBL_MAX;
         xassert(process_empty_row(lpp, row) == 0);
         goto done;
      }
      /* create transformation queue entry */
      info = lpp_append_tqe(lpp, LPP_FORCING_ROW, sizeof(FORCING_ROW));
      info->p = row->i;
      if (row->lb == row->ub)
      {  /* equality constraint */
         info->stat = LPX_NS;
         info->bnd = row->lb;
      }
      else if (at == 0)
      {  /* inequality constraint; the case l[p] = u' */
         info->stat = LPX_NL;
         xassert(row->lb != -DBL_MAX);
         info->bnd = row->lb;
      }
      else
      {  /* inequality constraint; the case u[p] = l' */
         info->stat = LPX_NU;
         xassert(row->ub != +DBL_MAX);
         info->bnd = row->ub;
      }
      info->ptr = NULL;
      /* formally the row p is removed from the problem at this point;
         however, its constraint coefficients are still needed, so its
         physical removal will be performed later */
      /* walk through the list of constraint coefficients of the row,
         save them, and fix the corresponding columns at appropriate
         bounds */
      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
      {  col = aij->col;
         /* save the constraint coefficient a[p,j] */
         lfx = dmp_get_atomv(lpp->tqe_pool, sizeof(LPPLFX));
         lfx->ref = col->j;
         lfx->flag = 0; /* will be set a bit later */
         lfx->val = aij->val;
         lfx->next = info->ptr;
         info->ptr = lfx;
         /* there must be no fixed columns */
         xassert(col->lb != col->ub);
         /* fix the column j at its lower or upper bound */
         if (at == 0 && aij->val < 0.0 || at != 0 && aij->val > 0.0)
         {  /* fix at lower bound */
            lfx->flag = LPX_NL;
            xassert(col->lb != -DBL_MAX);
            col->ub = col->lb;
         }
         else
         {  /* fix at upper bound */
            lfx->flag = LPX_NU;
            xassert(col->ub != +DBL_MAX);
            col->lb = col->ub;
         }
         /* before removing the column j from the problem we need to
            remove a[p,j] from the column list, because the row p is
            formally removed */
         if (aij->c_prev == NULL)
            aij->col->ptr = aij->c_next;
         else
            aij->c_prev->c_next = aij->c_next;
         if (aij->c_next == NULL)
            ;
         else
            aij->c_next->c_prev = aij->c_prev;
         /* remove the column j from the problem */
         process_fixed_col(lpp, col);
      }
      /* remove all a[p,j] from the row p */
      while (row->ptr != NULL)
      {  /* get a next element in the row */
         aij = row->ptr;
         /* remove it from the row list */
         row->ptr = aij->r_next;
         /* and return it to its pool */
         dmp_free_atom(lpp->aij_pool, aij, sizeof(LPPAIJ));
      }
      /* physically remove the row p (now it is empty) */
      lpp_remove_row(lpp, row);
done: return;
}

static void recover_forcing_row(LPP *lpp, FORCING_ROW *info)
{     LPPLFX *lfx, *that;
      double big, lambda, pi, temp;
      /* the row is not recovered yet */
      xassert(1 <= info->p && info->p <= lpp->nrows);
      xassert(lpp->row_stat[info->p] == 0);
      /* while all the corresponding columns are already recovered;
         they were fixed, so currently all they must be non-basic */
      for (lfx = info->ptr; lfx != NULL; lfx = lfx->next)
      {  xassert(1 <= lfx->ref && lfx->ref <= lpp->ncols);
         xassert(lpp->col_stat[lfx->ref] == LPX_NS);
      }
      /* choose a column q, whose dual value lambda[q] has wrong sign
         and reaches zero value last on changing pi[p] */
      that = NULL, big = 0.0;
      for (lfx = info->ptr; lfx != NULL; lfx = lfx->next)
      {  lambda = lpp->col_dual[lfx->ref];
         temp = fabs(lambda / lfx->val);
         switch (lfx->flag)