glpipp02.c

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

-- and therefore itself can be removed from the problem.
--
-- RETURNS
--
-- 0 - the row is primal feasible
-- 1 - the row is primal infeasible */

int ipp_analyze_row(IPP *ipp, IPPROW *row)
{     IPPCOL *col;
      IPPAIJ *aij;
      double lb, ub, eps, s;
      /* compute the implied lower bound l'; see (3) */
      lb = 0.0;
      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
      {  if (lb == -DBL_MAX) break;
         col = aij->col;
         if (aij->val > 0.0)
         {  if (col->lb == -DBL_MAX)
               lb = -DBL_MAX;
            else
               lb += aij->val * col->lb;
         }
         else
         {  if (col->ub == +DBL_MAX)
               lb = -DBL_MAX;
            else
               lb += aij->val * col->ub;
         }
      }
      /* compute the implied upper bound u'; see (4) */
      ub = 0.0;
      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
      {  if (ub == +DBL_MAX) break;
         col = aij->col;
         if (aij->val > 0.0)
         {  if (col->ub == +DBL_MAX)
               ub = +DBL_MAX;
            else
               ub += aij->val * col->ub;
         }
         else
         {  if (col->lb == -DBL_MAX)
               ub = +DBL_MAX;
            else
               ub += aij->val * col->lb;
         }
      }
      /* check for primal infeasibility */
      if (row->lb != -DBL_MAX)
      {  eps = 1e-5 * (1.0 + fabs(row->lb));
         if (ub < row->lb - eps) return 1;
      }
      if (row->ub != +DBL_MAX)
      {  eps = 1e-5 * (1.0 + fabs(row->ub));
         if (lb > row->ub + eps) return 1;
      }
      /* check if the row is implicitly fixed on its lower bound */
      if (row->lb != -DBL_MAX)
      {  eps = 1e-7 * (1.0 + fabs(row->lb));
         if (ub <= row->lb + eps)
         {  /* fix all the columns; see (7)-(8) */
            for (aij = row->ptr; aij != NULL; aij = aij->r_next)
            {  col = aij->col;
               s = (aij->val > 0.0 ? col->ub : col->lb);
               switch (ipp_tight_bnds(ipp, col, s, s))
               {  case 0:
                     /* bounds remain unchanged */
                     break;
                  case 1:
                     /* bounds have been changed, therefore activate
                        the column for further processing (removing) */
                     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);
            goto done;
         }
      }
      /* check if the row is implicitly fixed on its upper bound */
      if (row->ub != +DBL_MAX)
      {  eps = 1e-7 * (1.0 + fabs(row->ub));
         if (lb >= row->ub - eps)
         {  /* fix all the columns; see (5)-(6) */
            for (aij = row->ptr; aij != NULL; aij = aij->r_next)
            {  col = aij->col;
               s = (aij->val > 0.0 ? col->lb : col->ub);
               switch (ipp_tight_bnds(ipp, col, s, s))
               {  case 0:
                     /* bounds remain unchanged */
                     break;
                  case 1:
                     /* bounds have been changed, therefore activate
                        the column for further processing (removing) */
                     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);
            goto done;
         }
      }
      /* check whether the row can reach its lower bound */
      if (row->lb != -DBL_MAX)
      {  eps = 1.001e-7 * (1.0 + fabs(row->lb));
         if (lb >= row->lb - eps)
         {  /* the row cannot reach its lower bound, so the lower bound
               is redundant and therefore can be removed (note that the
               row cannot be an equality constraint at this point, since
               the case l' > u[p] = l[p] would be processed above */
            xassert(row->lb != row->ub);
            row->lb = -DBL_MAX;
            ipp_enque_row(ipp, row);
         }
      }
      /* check whether the row can reach its upper bound */
      if (row->ub != +DBL_MAX)
      {  eps = 1.001e-7 * (1.0 + fabs(row->ub));
         if (ub <= row->ub + eps)
         {  /* the row cannot reach its upper bound, so the upper bound
               is redundant and therefore can be removed (note that the
               row cannot be an equality constraint at this point, since
               the case u' < l[p] = u[p] would be processed above */
            xassert(row->lb != row->ub);
            row->ub = +DBL_MAX;
            ipp_enque_row(ipp, row);
         }
      }
done: return 0;
}

/*----------------------------------------------------------------------
-- GENERAL COLUMN ANALYSIS
--
-- Let column q be of general kind. Then its dual equality constraint
-- is the following:
--
--    sum a[i,q] * pi[i] + lambda[q] = c[q],                         (1)
--     i
--
-- from which it follows that:
--
--    lambda[q] = c[q] - sum a[i,q] * pi[i],                         (2)
--                        i
--
-- where c[q] is an objective coefficient of the column, a[i,q] are
-- constraints coefficents, pi[i] are dual values of the corresponding
-- rows.
--
-- If row i has no lower bound, pi[i] <= 0, and if row i has no upper
-- bound, pi[i] >= 0. This allows determining implied lower and upper
-- bounds of lambda[q].
--
-- If lambda[q] > 0, column q can be fixed on its lower bound, and if
-- lambda[q] < 0, column q can be fixed on its upper bound.
--
-- RETURNS
--
-- 0 - the column is dual feasible
-- 1 - the column is dual infeasible */

int ipp_analyze_col(IPP *ipp, IPPCOL *col)
{     IPPROW *row;
      IPPAIJ *aij;
      if (col->c > +1e-5)
      {  /* check if lambda[q] > 0 */
         for (aij = col->ptr; aij != NULL; aij = aij->c_next)
         {  row = aij->row;
            if (aij->val > 0.0)
            {  /* pi[i] <= 0 iff row i has no lower bound */
               if (row->lb != -DBL_MAX) goto done;
            }
            else
            {  /* pi[i] >= 0 iff row i has no upper bound */
               if (row->ub != +DBL_MAX) goto done;
            }
         }
         /* lambda[q] > 0; fix the column on its lower bound */
         if (col->lb == -DBL_MAX) return 1;
         ipp_tight_bnds(ipp, col, col->lb, col->lb);
         /* activate the column for further processing (removing) */
         ipp_enque_col(ipp, col);
      }
      else if (col->c < -1e-5)
      {  /* check if lambda[q] < 0 */
         for (aij = col->ptr; aij != NULL; aij = aij->c_next)
         {  row = aij->row;
            if (aij->val > 0.0)
            {  /* pi[i] >= 0 iff row i has no upper bound */
               if (row->ub != +DBL_MAX) goto done;
            }
            else
            {  /* pi[i] <= 0 iff row i has no lower bound */
               if (row->lb != -DBL_MAX) goto done;
            }
         }
         /* lambda[q] < 0; fix the column on its upper bound */
         if (col->ub == +DBL_MAX) return 1;
         ipp_tight_bnds(ipp, col, col->ub, col->ub);
         /* activate the column for further processing (removing) */
         ipp_enque_col(ipp, col);
      }
done: return 0;
}

/*----------------------------------------------------------------------
-- ipp_basic_tech - basic MIP presolve analysis.
--
-- SYNOPSIS
--
-- #include "glpipp.h"
-- int ipp_basic_tech(IPP *ipp);
--
-- DESCRIPTION
--
-- The routine ipp_basic_tech performs a basic MIP presolve analysis.
--
-- RETURNS
--
-- The routine returns one of the following codes:
--
-- 0 - neither primal nor dual infeasibility have been detected;
-- 1 - the original problem is primal infeasible;
-- 2 - the original problem is dual infeasible.
--
-- If the return code is non-zero, the resultant problem is invalid. */

int ipp_basic_tech(IPP *ipp)
{     IPPROW *row;
      IPPCOL *col;
      int nrows, ncols;
      /* the presolver uses two queues (organized in the form of doubly
         liked lists): one for active rows and other for active columns;
         initially these queues contain all rows and columns of the
         original problem; once a next active row or column has been
         selected for processing, it is removed from the corresponding
         queue; if a transformation applied to the selected row/column
         may affect other rows/columns, the latters are placed in the
         queue for the repeated processing */
      /* activate all rows and columns */
      nrows = 0;
      for (row = ipp->row_ptr; row != NULL; row = row->next)
         ipp_enque_row(ipp, row), nrows++;
      ncols = 0;
      for (col = ipp->col_ptr; col != NULL; col = col->next)
         ipp_enque_col(ipp, col), ncols++;
      /* process active rows and columns until both queues are empty */
      while (!(ipp->row_que == NULL && ipp->col_que == NULL))
      {  /* process active rows */
         while (ipp->row_que != NULL)
         {  /* select an active row */
            row = ipp->row_que;
            /* and remove it from the queue */
            ipp_deque_row(ipp, row);
            /* process the row */
            if (row->lb == -DBL_MAX && row->ub == +DBL_MAX)
            {  /* process free row */
               ipp_free_row(ipp, row);
            }
            else if (row->ptr == NULL)
            {  /* process empty row */
               if (ipp_empty_row(ipp, row)) return 1;
            }
            else if (row->ptr != NULL && row->ptr->r_next == NULL)
            {  /* process row singleton */
               if (ipp_row_sing(ipp, row)) return 1;
            }
            else
            {  /* general row analysis */
               if (ipp_analyze_row(ipp, row)) return 1;
            }
         }
         /* process active columns */
         while (ipp->col_que != NULL)
         {  /* select an active column */
            col = ipp->col_que;
            /* and remove it from the queue */
            ipp_deque_col(ipp, col);
            /* process the column */
            if (col->lb == col->ub)
            {  /* process fixed column */
               ipp_fixed_col(ipp, col);
            }
            else if (col->ptr == NULL)
            {  /* process empty column */
               if (ipp_empty_col(ipp, col)) return 2;
            }
            else
            {  /* general column analysis */
               if (ipp_analyze_col(ipp, col)) return 2;
            }
         }
      }
      for (row = ipp->row_ptr; row != NULL; row = row->next) nrows--;
      for (col = ipp->col_ptr; col != NULL; col = col->next) ncols--;
      xprintf("ipp_basic_tech:  %d row(s) and %d column(s) removed\n",
         nrows, ncols);
      return 0;
}

/*----------------------------------------------------------------------
-- REDUCE COLUMN BOUNDS
--
-- Given a row (constraint) the routine reduce_bounds tries to reduce
-- column bounds replacing them by implied bounds.
--
-- To obtain an implied lower/upper bound of a column x[j] the routine
-- solves the following LP relaxation:
--
--    minimize/maximize x[j]                                         (1)
--
--    subject to        L <= sum a[j] * x[j] <= U                    (2)
--                            j
--
--                      l[j] <= x[j] <= u[j]                         (3)
--
-- where (2) is a given row (constraint).
--
-- If a bound of a column has been changed, the column is placed in the
-- active queue for further processing. */

static int reduce_bounds(IPP *ipp, IPPROW *row)
{     IPPCOL *col, *c_min, *c_max;
      IPPAIJ *aij;
      int flag;
      double f_min, f_max, ff_min, ff_max, lb, ub, delta;
      /* compute implied lower bound of the row */
      c_min = NULL, f_min = 0.0;
      for (aij = row->ptr; aij != NULL; aij = aij->r_next)