glpipp02.c

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

/* glpipp02.c */

/***********************************************************************
*  This code is part of GLPK (GNU Linear Programming Kit).
*
*  Copyright (C) 2000,01,02,03,04,05,06,07,08,2009 Andrew Makhorin,
*  Department for Applied Informatics, Moscow Aviation Institute,
*  Moscow, Russia. All rights reserved. E-mail: <[e-mail]>.
*
*  GLPK is free software: you can redistribute it and/or modify it
*  under the terms of the GNU General Public License as published by
*  the Free Software Foundation, either version 3 of the License, or
*  (at your option) any later version.
*
*  GLPK is distributed in the hope that it will be useful, but WITHOUT
*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
*  License for more details.
*
*  You should have received a copy of the GNU General Public License
*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/

#include "glpipp.h"

/*----------------------------------------------------------------------
-- FREE ROW
--
-- Let row p be free, i.e. it has no bounds:
--
--    -inf < y[p] = sum a[p,j] * x[j] < +inf.                        (1)
--
-- The corresponding "constraint" can never be active, so the free row
-- is redundant and can be removed from the problem.
--
-- RETURNS
--
-- None. */

void ipp_free_row(IPP *ipp, IPPROW *row)
{     /* process free row */
      IPPAIJ *aij;
      /* the row must be free */
      xassert(row->lb == -DBL_MAX && row->ub == +DBL_MAX);
      /* activate corresponding columns for further processing */
      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
         ipp_enque_col(ipp, aij->col);
      /* remove the row from the problem */
      ipp_remove_row(ipp, row);
      return;
}

/*----------------------------------------------------------------------
-- FIXED COLUMN
--
-- Let column q be fixed:
--
--    x[q] = s[q],                                                   (1)
--
-- where s[q] is a given value. Then it can be substituted to all rows
-- and the objective function and removed from the problem.
--
-- Consider a row i:
--
--    lb[i] <= y[i] =  sum a[i,j] * x[j] + a[i,q] * x[q] <= ub[i].   (2)
--                   j != q
--
-- Substituting x[q] from (1) to (2) transforms the row i to:
--
--    lb'[i] <= y'[i] =  sum a[i,j] * x[j] <= ub'[i],                (3)
--                     j != q
--
-- with new lower and upper bounds
--
--    lb'[i] = lb[i] - a[i,q] * s[q]                                 (4)
--
--    ub'[i] = ub[i] - a[i,q] * s[q]                                 (5)
--
-- of a new auxiliary variable y'[i] for the row i.
--
-- RECOVERING
--
-- Value of x[q] is determined using the formula (1).
--
-- RETURNS
--
-- None. */

struct fixed_col
{     /* fixed column */
      int q;
      /* reference number of a fixed column */
      double s;
      /* value at which the column is fixed */
};

void ipp_fixed_col(IPP *ipp, IPPCOL *col)
{     /* process fixed column */
      struct fixed_col *info;
      IPPROW *row;
      IPPAIJ *aij;
      double temp;
      /* the column must be fixed */
      xassert(col->lb == col->ub);
      /* create transformation queue entry */
      info = ipp_append_tqe(ipp, IPP_FIXED_COL, sizeof(*info));
      info->q = col->j;
      info->s = col->lb;
      /* update bounds of corresponding rows and activate them for
         further processing */
      for (aij = col->ptr; aij != NULL; aij = aij->c_next)
      {  row = aij->row;
         /* see (4) and (5) */
         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;
         }
         ipp_enque_row(ipp, row);
      }
      /* update the constant term of the objective function */
      ipp->c0 += col->c * info->s;
      /* remove the column from the problem */
      ipp_remove_col(ipp, col);
      return;
}

void ipp_fixed_col_r(IPP *ipp, void *_info)
{     /* recover fixed column */
      struct fixed_col *info = _info;
      xassert(1 <= info->q && info->q <= ipp->ncols);
      xassert(ipp->col_stat[info->q] == 0);
      ipp->col_stat[info->q] = 1;
      ipp->col_mipx[info->q] = info->s;
      return;
}

/*----------------------------------------------------------------------
-- EMPTY ROW
--
-- Let row p be empty, i.e. all its coefficients are zero:
--
--    l[p] <= y[p] = 0 <= u[p].                                      (1)
--
-- If l[p] <= 0 and u[p] >= 0, the row is redundant and therefore can
-- be removed from the problem. Otherwise, if l[p] > 0 or u[p] < 0, the
-- row is primal infeasible.
--
-- RETURNS
--
-- 0 - the row is primal feasible and has been removed
-- 1 - the row is primal infeasible */

int ipp_empty_row(IPP *ipp, IPPROW *row)
{     /* process empty row */
      double eps;
      /* the row must be empty */
      xassert(row->ptr == NULL);
      /* check for primal infeasibility */
      eps = 1e-5;
      if (row->lb > +eps || row->ub < -eps) return 1;
      /* make the row free */
      row->lb = -DBL_MAX, row->ub = +DBL_MAX;
      /* and activate it for further processing (removing) */
      ipp_enque_row(ipp, row);
      return 0;
}

/*----------------------------------------------------------------------
-- EMPTY COLUMN
--
-- Let column q be empty, i.e. it has zero coefficients in all rows.
-- Then its dual value is equal to its objective coefficient:
--
--    lambda[q] = c[q].                                              (1)
--
-- So, if c[q] > 0 (c[q] < 0), the variable x[q] should be fixed at its
-- lower (upper) bound, and if c[q] = 0, the variable x[q] can be fixed
-- at any value. And being fixed the variable x[q] can be removed from
-- the problem. If the variable x[q] cannot be properly fixed, i.e. if
-- c[q] > 0 (c[q] < 0) and x[q] has no lower (upper) bound, the column
-- is dual infeasible.
--
-- RETURNS
--
-- 0 - the column is dual feasible and has been fixed and removed
-- 1 - the column is dual infeasible */

int ipp_empty_col(IPP *ipp, IPPCOL *col)
{     /* process empty column */
      double eps;
      /* the column must be empty */
      xassert(col->ptr == NULL);
      /* check for dual infeasibility */
      eps = 1e-5;
      if (col->c > +eps && col->lb == -DBL_MAX ||
          col->c < -eps && col->ub == +DBL_MAX) return 1;
      /* fix the column at proper bound */
      if (col->lb == -DBL_MAX && col->ub == +DBL_MAX)
      {  /* free variable */
         col->lb = col->ub = 0.0;
      }
      else if (col->ub == +DBL_MAX)
lo:   {  /* variable with lower bound */
         col->ub = col->lb;
      }
      else if (col->lb == -DBL_MAX)
up:   {  /* variable with upper bound */
         col->lb = col->ub;
      }
      else if (col->lb != col->ub)
      {  /* double bounded variable */
         if (col->c > 0.0) goto lo;
         if (col->c < 0.0) goto up;
         if (fabs(col->lb) <= fabs(col->ub)) goto lo; else goto up;
      }
      else
      {  /* fixed variable */
         /* nop */;
      }
      /* and activate it for further processing (removing) */
      ipp_enque_col(ipp, col);
      return 0;
}

/*----------------------------------------------------------------------
-- ROW SINGLETON
--
-- Let row p be a constraing having only one column:
--
--    l[p] <= y[p] = a[p,q] * x[q] <= u[p].                          (1)
--
-- 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.
--
-- RETURNS
--
-- 0 - the row is primal feasible
-- 1 - the row is primal infeasible */

int ipp_row_sing(IPP *ipp, IPPROW *row)
{     /* process row singleton */
      IPPCOL *col;
      IPPAIJ *aij;
      double lb, ub;
      /* the row must have only one constraint coefficient */
      xassert(row->ptr != NULL && row->ptr->r_next == NULL);
      /* 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);
      }
      /* tight current column bounds */
      col = aij->col;
      switch (ipp_tight_bnds(ipp, col, lb, ub))
      {  case 0:
            /* bounds remain unchanged */
            break;
         case 1:
            /* bounds have been changed, therefore activate the column
               for further processing */
            ipp_enque_col(ipp, col);
            break;
         case 2:
            /* implied bounds are in conflict */
            return 1;
         default:
            xassert(ipp != ipp);
      }
      /* make the row free */
      row->lb = -DBL_MAX, row->ub = +DBL_MAX;
      /* and activate it for further processing (removing) */
      ipp_enque_row(ipp, row);
      return 0;
}

/*----------------------------------------------------------------------
-- GENERAL ROW ANALYSIS
--
-- 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)
--
-- The analysis is based on implied lower and upper bounds l' and u' of
-- the primal value y[p]:
--
--                      ( l[j], if a[p,j] > 0 )
--    l' = sum a[p,j] * <                     > ,                    (3)
--          j           ( u[j], if a[p,j] < 0 )
--
--                      ( u[j], if a[p,j] > 0 )
--    u' = sum a[p,j] * <                     > .                    (4)
--          j           ( l[j], if a[p,j] < 0 )
--
-- If l' > u[p] or u' < l[p], the row is primal infeasible.
--
-- If l' = u[p], all the variables x[j] with non-zero coefficients in
-- the row p can be fixed on their bounds as follows:
--
--    if a[p,j] > 0, x[j] is fixed on l[j];                          (5)
--
--    if a[p,j] < 0, x[j] is fixed on u[j].                          (6)
--
-- If u' = l[p], all the variables x[j] with non-zero coefficients in
-- the row p can be fixed on their bounds as follows:
--
--    if a[p,j] > 0, x[j] is fixed on u[j];                          (7)
--
--    if a[p,j] < 0, x[j] is fixed on l[j].                          (8)
--
-- In both cases l' = u[p] and u' = l[p] after fixing variables the row
-- becomes redundant and therefore can be removed from the problem.
--
-- If l' > l[p], the row cannot be active on its lower bound, so the
-- lower bound can be removed. Analogously, If u' < u[p], the row cannot
-- be active on its upper bound, so the upper bound can be removed. If
-- both lower and upper bounds have been removed, the row becomes free