glplpp02.c

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

            row->ub = (row->lb -= aij->val * col->lb);
         else
         {  if (row->lb != -DBL_MAX) row->lb -= aij->val * col->lb;
            if (row->ub != +DBL_MAX) row->ub -= aij->val * col->lb;
         }
      }
      /* update the constant term of the objective function */
      lpp->c0 += col->c * col->lb;
      /* remove the column from the problem */
      lpp_remove_col(lpp, col);
      return;
}

static void recover_fixed_col(LPP *lpp, FIXED_COL *info)
{     LPPLFE *lfe;
      double dual;
      xassert(1 <= info->q && info->q <= lpp->ncols);
      xassert(lpp->col_stat[info->q] == 0);
      /* the fixed column is non-basic */
      lpp->col_stat[info->q] = LPX_NS;
      /* set primal value x[q] */
      lpp->col_prim[info->q] = info->val;
      /* compute dual value lambda[q] using the formula (7) and update
         primal values of rows using the formula (8) */
      dual = info->c;
      for (lfe = info->ptr; lfe != NULL; lfe = lfe->next)
      {  xassert(1 <= lfe->ref && lfe->ref <= lpp->nrows);
         xassert(lpp->row_stat[lfe->ref] != 0);
         dual -= lfe->val * lpp->row_dual[lfe->ref];
         lpp->row_prim[lfe->ref] += lfe->val * info->val;
      }
      lpp->col_dual[info->q] = dual;
      return;
}

/*----------------------------------------------------------------------
-- ROW SINGLETON (EQUALITY CONSTRAINT)
--
-- *Processing*
--
-- Let row p be an equality constraint that has the only column:
--
--    y[p] = a[p,q] * x[q] = b[p],                                   (1)
--
-- in which case it implies fixing x[q]:
--
--    x[q] = b[p] / a[p,q].                                          (2)
--
-- So, if (2) doesn't conflict with bounds of x[q], the column q can be
-- fixed and therefore removed from the problem. Then the row p becomes
-- redundant and therefore can be removed from the problem. Otherwise,
-- if (2) conflicts with bounds of x[q], the row is primal infeasible.
--
-- Note that at first the row p is removed. Then the column q is fixed
-- and processed as a fixed column (see section "Fixed Column").
--
-- *Recovering*
--
-- On entry to the recovering routine the column q is already recovered
-- as if it were a fixed column, i.e. it is non-basic with primal value
-- x[q] and dual value lambda[q]. The routine makes it basic with the
-- same primal value and zero dual value.
--
-- Then the recovering routine makes the row p non-basic, and computes
-- its primal value y[p] using the formula (1) and its dual value pi[p]
-- using the formula:
--
--    pi[p] = lambda[q] / a[p,q],                                    (3)
--
-- where lambda[q] is a dual value, which the column q has on entry to
-- the recovering routine. */

typedef struct
{     /* row singleton (equality constraint) */
      int p;
      /* row reference number */
      int q;
      /* column reference number */
      double apq;
      /* constraint coefficient */
} ROW_SNGTON1;

static int process_row_sngton1(LPP *lpp, LPPROW *row)
{     ROW_SNGTON1 *info;
      LPPCOL *col;
      LPPAIJ *aij;
      double val, eps;
      /* the row must have the only non-zero constraint coefficient and
         be an equality constraint */
      xassert(row->ptr != NULL && row->ptr->r_next == NULL);
      xassert(row->lb == row->ub);
      /* compute the implied value; see (2) */
      aij = row->ptr;
      val = row->lb / aij->val;
      /* check for primal infeasibility */
      col = aij->col;
      if (col->lb != -DBL_MAX)
      {  eps = 1e-5 * (1.0 + fabs(col->lb));
         if (val < col->lb - eps) return 1;
      }
      if (col->ub != +DBL_MAX)
      {  eps = 1e-5 * (1.0 + fabs(col->ub));
         if (val > col->ub + eps) return 1;
      }
      /* create transformation queue entry */
      info = lpp_append_tqe(lpp, LPP_ROW_SNGTON1, sizeof(ROW_SNGTON1));
      info->p = row->i;
      info->q = col->j;
      info->apq = aij->val;
      /* remove the row from the problem */
      lpp_remove_row(lpp, row);
      /* fix the column */
      col->lb = col->ub = val;
      /* and also remove it from the problem */
      process_fixed_col(lpp, col);
      return 0;
}

static void recover_row_sngton1(LPP *lpp, ROW_SNGTON1 *info)
{     /* the row is not recovered yet */
      xassert(1 <= info->p && info->p <= lpp->nrows);
      xassert(lpp->row_stat[info->p] == 0);
      /* while the column is already recovered; it was removed as if it
         were fixed, so currently it must be non-basic */
      xassert(1 <= info->q && info->q <= lpp->ncols);
      xassert(lpp->col_stat[info->q] == LPX_NS);
      /* the row is actually active */
      lpp->row_stat[info->p] = LPX_NS;
      /* compute primal value y[p] using the formula (1) */
      lpp->row_prim[info->p] = info->apq * lpp->col_prim[info->q];
      /* compute dual value pi[p] using the formula (3) */
      lpp->row_dual[info->p] = lpp->col_dual[info->q] / info->apq;
      /* the column is actually basic */
      lpp->col_stat[info->q] = LPX_BS;
      lpp->col_dual[info->q] = 0.0;
      return;
}

/*----------------------------------------------------------------------
-- ROW SINGLETON (INEQUALITY CONSTRAINT)
--
-- *Processing*
--
-- Let row p be an inequality constraint that has the only column:
--
--    l[p] <= y[p] = a[p,q] * x[q] <= u[p],                          (1)
--
-- where l[p] != u[p] and at least one of the bounds l[p] and u[p] is
-- finite. In this case the row implies bounds of x[q]:
--
--    l' <= x[q] <= u',                                              (2)
--
-- where:
--
--    if a[p,q] > 0:   l' = l[p] / a[p,q],   u' = u[p] / a[p,q];     (3)
--
--    if a[p,q] < 0:   l' = u[p] / a[p,q],   u' = l[p] / a[p,q].     (4)
--
-- If the bounds l' and u' conflict with own bounds of x[q], i.e. if
-- l' > u[q] or u' < l[q], the row is primal infeasible. Otherwise the
-- range of x[q] can be tightened as follows:
--
--    max(l[q], l') <= x[q] <= min(u[q], u'),                        (5)
--
-- in which case the row becomes redundant and therefore can be removed
-- from the problem.
--
-- *Recovering*
--
-- On entry to the recovering routine the column q is already recovered
-- as if it had the modified bounds (5). Then the row p is recovered as
-- follows.
--
-- Let the column q be basic. In this case the row p is inactive, so
-- its primal value y[p] is computed with the formula (1) and its dual
-- value pi[p] is set to zero.
--
-- Let the column q be non-basic on lower bound. If l' <= l[q], the row
-- p is inactive and recovered in the same way as above. Otherwise, if
-- l' > l[q], the column q is actually basic while the row p is active
-- on its lower bound (a[p,q] > 0) or on its upper bound (a[p,q] < 0).
-- In the latter case the row's primal value y[p] is computed using the
-- formula (1) and the dual value pi[p] is computed using the formula:
--
--    pi[p] = lambda[q] / a[p,q],                                    (6)
--
-- where lambda[q] is a dual value, which the column q has on entry to
-- the recovering routine.
--
-- Let the column q be non-basic on upper bound. If u' >= u[q], the row
-- p is inactive and recovered in the same way as above. Otherwise, if
-- u' < u[q], the column q is actually basic while the row p is active
-- on its lower bound (a[p,q] < 0) or on its upper bound (a[p,q] > 0).
-- Then the row's primal value y[p] is computed using the formula (1)
-- and the dual value pi[p] is computed using the formula (6). */

typedef struct
{     /* row singleton (inequality constraint) */
      int p;
      /* row reference number */
      int q;
      /* column reference number */
      double apq;
      /* constraint coefficient */
      int lb_changed;
      /* this flag is set if the lower bound of the column was changed,
         i.e. if l' > l[q]; see (2)--(5) */
      int ub_changed;
      /* this flag is set if the upper bound of the column was changed,
         i.e. if u' < u[q]; see (2)--(5) */
} ROW_SNGTON2;

static int process_row_sngton2(LPP *lpp, LPPROW *row)
{     ROW_SNGTON2 *info;
      LPPCOL *col;
      LPPAIJ *aij;
      double lb, ub, eps;
      /* the row must have the only non-zero constraint coefficient and
         be an inequality constraint */
      xassert(row->ptr != NULL && row->ptr->r_next == NULL);
      xassert(row->lb != row->ub);
      /* if the row is free, just remove it from the problem */
      if (row->lb == -DBL_MAX && row->ub == +DBL_MAX)
      {  process_free_row(lpp, row);
         goto done;
      }
      /* compute the impled bounds l' and u'; see (2)--(4) */
      aij = row->ptr;
      if (aij->val > 0.0)
      {  lb = (row->lb == -DBL_MAX ? -DBL_MAX : row->lb / aij->val);
         ub = (row->ub == +DBL_MAX ? +DBL_MAX : row->ub / aij->val);
      }
      else
      {  lb = (row->ub == +DBL_MAX ? -DBL_MAX : row->ub / aij->val);
         ub = (row->lb == -DBL_MAX ? +DBL_MAX : row->lb / aij->val);
      }
      /* check for primal infeasibility */
      col = aij->col;
      if (col->lb != -DBL_MAX)
      {  eps = 1e-5 * (1.0 + fabs(col->lb));
         if (ub < col->lb - eps) return 1;
      }
      if (col->ub != +DBL_MAX)
      {  eps = 1e-5 * (1.0 + fabs(col->ub));
         if (lb > col->ub + eps) return 1;
      }
      /* if the column q is fixed, it can be substituted and removed,
         in which case the row p becomes empty and being feasible (as
         already checked above) can be also removed */
      if (col->lb == col->ub)
      {  /* substitute and remove the column q */
         process_fixed_col(lpp, col);
         /* now the row p became empty and redundant */
         xassert(row->ptr == NULL);
         /* remove the row p */
         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_ROW_SNGTON2, sizeof(ROW_SNGTON2));
      info->p = row->i;
      info->q = col->j;
      info->apq = aij->val;
      info->lb_changed = (lb != -DBL_MAX && lb > col->lb);
      info->ub_changed = (ub != +DBL_MAX && ub < col->ub);
      /* tighten bounds of the column, if necessary */
      if (info->lb_changed) col->lb = lb;
      if (info->ub_changed) col->ub = ub;
      /* the row is now redundant; remove it from the problem */
      lpp_remove_row(lpp, row);
      /* if modified bounds of the column are close to each other, the
         column can be fixed, substituted, and therefore removed */
      if (col->lb != -DBL_MAX && col->ub != +DBL_MAX)
      {  eps = 1e-7 * (1.0 + fabs(col->lb));
         if (fabs(col->lb - col->ub) <= eps)
         {  if (fabs(col->lb) <= fabs(col->ub))
               col->ub = col->lb;
            else
               col->lb = col->ub;
            process_fixed_col(lpp, col);
         }
      }
done: return 0;
}

static void recover_row_sngton2(LPP *lpp, ROW_SNGTON2 *info)
{     /* the row is not recovered yet */
      xassert(1 <= info->p && info->p <= lpp->nrows);
      xassert(lpp->row_stat[info->p] == 0);
      /* while the column is already recovered */
      xassert(1 <= info->q && info->q <= lpp->ncols);
      xassert(lpp->col_stat[info->q] != 0);
      /* analyze the column status */
      switch (lpp->col_stat[info->q])
      {  case LPX_BS:
            /* the column is basic, 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_NL:
nl:         /* the column is non-basic on lower bound */
            if (info->lb_changed)
            {  /* it is not its own lower bound, so actually the row is
                  active */
               if (info->apq > 0.0)
                  lpp->row_stat[info->p] = LPX_NL;
               else
                  lpp->row_stat[info->p] = LPX_NU;
               lpp->row_prim[info->p] =
                  info->apq * lpp->col_prim[info->q];
               lpp->row_dual[info->p] =
                  lpp->col_dual[info->q] / info->apq;
               /* and the column is basic */
               lpp->col_stat[info->q] = LPX_BS;
               lpp->col_dual[info->q] = 0.0;
            }
            else
            {  /* it is its own lower bound, so the row is incative */
               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_NU:
nu:         /* the column is non-basic on upper bound */
            if (info->ub_changed)
            {  /* it is not its own upper bound, so actually the row is
                  active */
               if (info->apq > 0.0)
                  lpp->row_stat[info->p] = LPX_NU;
               else
                  lpp->row_stat[info->p] = LPX_NL;
               lpp->row_prim[info->p] =
                  info->apq * lpp->col_prim[info->q];
               lpp->row_dual[info->p] =
                  lpp->col_dual[info->q] / info->apq;
               /* and the column is basic */
               lpp->col_stat[info->q] = LPX_BS;
               lpp->col_dual[info->q] = 0.0;