glplpp02.c

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

         {  case LPX_NL:
               /* x[j] >= l[j], therefore lambda[j] >= 0 */
               if (lambda < 0.0 && big < temp) that = lfx, big = temp;
               break;
            case LPX_NU:
               /* x[j] <= u[j], therefore lambda[j] <= 0 */
               if (lambda > 0.0 && big < temp) that = lfx, big = temp;
               break;
            default:
               xassert(lfx != lfx);
         }
      }
      /* recover the row p and all the corresponding columns */
      if (that == NULL)
      {  /* the row p is inactive; all columns are non-basic */
         lpp->row_stat[info->p] = LPX_BS;
         lpp->row_prim[info->p] = info->bnd;
         lpp->row_dual[info->p] = 0.0;
         for (lfx = info->ptr; lfx != NULL; lfx = lfx->next)
            lpp->col_stat[lfx->ref] = lfx->flag;
      }
      else
      {  /* the row p is active, the column q is basic, and all other
            columns are non-basic */
         /* compute dual value pi[p] = lambda'[q] / a[p,q]; see (5) */
         pi = lpp->col_dual[that->ref] / that->val;
         /* recover the row p */
         lpp->row_stat[info->p] = info->stat;
         lpp->row_prim[info->p] = info->bnd;
         lpp->row_dual[info->p] = pi;
         /* recover the corresponding columns */
         for (lfx = info->ptr; lfx != NULL; lfx = lfx->next)
         {  if (lfx == that)
            {  /* recover the column q, which is basic */
               lpp->col_stat[lfx->ref] = LPX_BS;
               lpp->col_dual[lfx->ref] = 0.0;
            }
            else
            {  /* recover some column j, which is non-basic */
               lpp->col_stat[lfx->ref] = lfx->flag;
               lpp->col_dual[lfx->ref] -= lfx->val * pi;
            }
         }
      }
      return;
}

/*----------------------------------------------------------------------
-- 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
-- and therefore itself can be removed from the problem. */

static int analyze_row(LPP *lpp, LPPROW *row)
{     LPPCOL *col;
      LPPAIJ *aij;
      double lb, ub, eps;
      /* 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)
         {  process_forcing_row(lpp, row, 0);
            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)
         {  process_forcing_row(lpp, row, 1);
            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;
            lpp_enque_row(lpp, 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;
            lpp_enque_row(lpp, row);
         }
      }
done: return 0;
}

/*----------------------------------------------------------------------
-- lpp_presolve - LP presolve analysis.
--
-- *Synopsis*
--
-- #include "glplpp.h"
-- int lpp_presolve(LPP *lpp);
--
-- *Description*
--
-- The routine lpp_presolve performs a presolve analysis and transforms
-- the original problem to a resultant problem.
--
-- *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 lpp_presolve(LPP *lpp)
{     LPPROW *row;
      LPPCOL *col;
      /* 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 */
      while (lpp->row_que != NULL || lpp->col_que != NULL)
      {  /* process active rows */
         while (lpp->row_que != NULL)
         {  /* select a next active row */
            row = lpp->row_que;
            /* and remove it from the queue */
            lpp_deque_row(lpp, row);
            /* process the row */
            if (row->ptr == NULL)
            {  /* remove empty row */
               if (process_empty_row(lpp, row)) return 1;
            }
            else if (row->lb == -DBL_MAX && row->ub == +DBL_MAX)
            {  /* remove free row */
               process_free_row(lpp, row);
            }
            else if (row->ptr != NULL && row->ptr->r_next == NULL &&
                     row->lb == row->ub)
            {  /* remove row singleton (equality constraint) */
               if (process_row_sngton1(lpp, row)) return 1;
            }
            else if (row->ptr != NULL && row->ptr->r_next == NULL &&
                     row->lb != row->ub)
            {  /* remove row singleton (inequality constraint) */
               if (process_row_sngton2(lpp, row)) return 1;
            }
            else
            {  /* general row analysis */
               if (analyze_row(lpp, row)) return 1;
            }
         }
         /* process active columns */
         while (lpp->col_que != NULL)
         {  /* select a next active column */
            col = lpp->col_que;
            /* and remove it from the queue */
            lpp_deque_col(lpp, col);
            /* process the column */
            if (col->ptr == NULL)
            {  /* remove empty column */
               if (process_empty_col(lpp, col)) return 2;
            }
            else if (col->lb == col->ub)
            {  /* remove fixed column */
               process_fixed_col(lpp, col);
            }
            else if (col->ptr != NULL && col->ptr->c_next == NULL &&
                     col->ptr->row->lb == col->ptr->row->ub)
            {  /* remove column singleton (implied slack variable) */
               process_col_sngton1(lpp, col);
            }
            else if (col->ptr != NULL && col->ptr->c_next == NULL &&
                     col->ptr->row->lb != col->ptr->row->ub)
            {  /* remove column singleton (implied free variable) */
               if (process_col_sngton2(lpp, col)) return 2;
            }
            else
            {  /* general column analysis */
               /* (not implemented yet) */
            }
         }
      }
      return 0;
}

/*----------------------------------------------------------------------
-- lpp_postsolve - LP postsolve processing.
--
-- *Synopsis*
--
-- #include "glplpp.h"
-- void lpp_postsolve(LPP *lpp);
--
-- *Description*
--
-- The routine lpp_postsolve performs a postsolve processing to recover
-- a solution of the original problem. It is assumed that a solution of
-- the resultant problem is loaded into the presolver workspace. */

void lpp_postsolve(LPP *lpp)
{     LPPTQE *tqe;
      for (tqe = lpp->tqe_list; tqe != NULL; tqe = tqe->next)
      {  switch (tqe->type)
         {  case LPP_EMPTY_ROW:
               recover_empty_row(lpp, tqe->info);
               break;
            case LPP_EMPTY_COL:
               recover_empty_col(lpp, tqe->info);
               break;
            case LPP_FREE_ROW:
               recover_free_row(lpp, tqe->info);
               break;
            case LPP_FIXED_COL:
               recover_fixed_col(lpp, tqe->info);
               break;
            case LPP_ROW_SNGTON1:
               recover_row_sngton1(lpp, tqe->info);
               break;
            case LPP_ROW_SNGTON2:
               recover_row_sngton2(lpp, tqe->info);
               break;
            case LPP_COL_SNGTON1:
               recover_col_sngton1(lpp, tqe->info);
               break;
            case LPP_COL_SNGTON2:
               recover_col_sngton2(lpp, tqe->info);
               break;
            case LPP_FORCING_ROW:
               recover_forcing_row(lpp, tqe->info);
               break;
            default:
               xassert(tqe->type != tqe->type);
         }
      }
      return;
}

/* eof */