glplpp02.c

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

/* glplpp02.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 "glpapi.h"
#include "glplpp.h"
#define dmp_get_atomv dmp_get_atom

/*----------------------------------------------------------------------
-- EMPTY ROW
--
-- *Processing*
--
-- 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.
--
-- *Recovering*
--
-- Being redundant the empty row is non-active constraint.
--
-- Primal value y[p] is zero as it follows from (1).
--
-- Dual value pi[p] is zero, because the constraint is non-active. */

typedef struct
{     /* empty row */
      int p;
      /* reference number of empty row */
} EMPTY_ROW;

static int process_empty_row(LPP *lpp, LPPROW *row)
{     EMPTY_ROW *info;
      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;
      /* create transformation queue entry */
      info = lpp_append_tqe(lpp, LPP_EMPTY_ROW, sizeof(EMPTY_ROW));
      info->p = row->i;
      /* remove the row from the problem */
      lpp_remove_row(lpp, row);
      return 0;
}

static void recover_empty_row(LPP *lpp, EMPTY_ROW *info)
{     xassert(1 <= info->p && info->p <= lpp->nrows);
      xassert(lpp->row_stat[info->p] == 0);
      /* the empty row is non-active */
      lpp->row_stat[info->p] = LPX_BS;
      /* both primal and dual values are zero */
      lpp->row_prim[info->p] = 0.0;
      lpp->row_dual[info->p] = 0.0;
      return;
}

/*----------------------------------------------------------------------
-- EMPTY COLUMN
--
-- *Processing*
--
-- 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.
--
-- *Recovering*
--
-- Being fixed on its bound the empty column is non-basic.
--
-- Primal value x[q] is a value, at which the column is fixed.
--
-- Dual value lambda[q] is the objective coefficient c[q] as it follows
-- from (1). */

typedef struct
{     /* empty column */
      int q;
      /* reference number of empty column */
      int stat;
      /* column status:
         LPX_NL - non-basic variable on lower bound
         LPX_NU - non-basic variable on upper bound
         LPX_NF - non-basic free variable
         LPX_NS - non-basic fixed variable */
      double prim;
      /* primal value */
      double dual;
      /* dual value */
} EMPTY_COL;

static int process_empty_col(LPP *lpp, LPPCOL *col)
{     EMPTY_COL *info;
      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;
      /* create transformation queue entry */
      info = lpp_append_tqe(lpp, LPP_EMPTY_COL, sizeof(EMPTY_COL));
      info->q = col->j;
      if (col->lb == -DBL_MAX && col->ub == +DBL_MAX)
      {  /* free variable */
         info->stat = LPX_NF;
         info->prim = 0.0;
      }
      else if (col->ub == +DBL_MAX)
lo:   {  /* variable with lower bound */
         info->stat = LPX_NL;
         info->prim = col->lb;
      }
      else if (col->lb == -DBL_MAX)
up:   {  /* variable with upper bound */
         info->stat = LPX_NU;
         info->prim = 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 */
         info->stat = LPX_NS;
         info->prim = col->lb;
      }
      info->dual = col->c;
      /* update the constant term of the objective function */
      lpp->c0 += col->c * info->prim;
      /* remove the column from the problem */
      lpp_remove_col(lpp, col);
      return 0;
}

static void recover_empty_col(LPP *lpp, EMPTY_COL *info)
{     xassert(1 <= info->q && info->q <= lpp->ncols);
      xassert(lpp->col_stat[info->q] == 0);
      lpp->col_stat[info->q] = info->stat;
      lpp->col_prim[info->q] = info->prim;
      lpp->col_dual[info->q] = info->dual;
      return;
}

/*----------------------------------------------------------------------
-- FREE ROW
--
-- *Processing*
--
-- 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.
--
-- *Recovering*
--
-- Being redundant the free row is non-active constraint.
--
-- Primal value y[p] is computed using the formula (1).
--
-- Dual value pi[p] is zero, because the constraint is non-active. */

typedef struct
{     /* free row */
      int p;
      /* reference number of a free row */
      LPPLFE *ptr;
      /* list of non-zero coefficients a[p,j] */
} FREE_ROW;

static void process_free_row(LPP *lpp, LPPROW *row)
{     FREE_ROW *info;
      LPPCOL *col;
      LPPAIJ *aij;
      LPPLFE *lfe;
      /* the row must be free */
      xassert(row->lb == -DBL_MAX && row->ub == +DBL_MAX);
      /* create transformation queue entry */
      info = lpp_append_tqe(lpp, LPP_FREE_ROW, sizeof(FREE_ROW));
      info->p = row->i;
      info->ptr = NULL;
      /* save coefficients of the row */
      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
      {  col = aij->col;
         lfe = dmp_get_atomv(lpp->tqe_pool, sizeof(LPPLFE));
         lfe->ref = col->j;
         lfe->val = aij->val;
         lfe->next = info->ptr;
         info->ptr = lfe;
      }
      /* remove the row from the problem */
      lpp_remove_row(lpp, row);
      return;
}

static void recover_free_row(LPP *lpp, FREE_ROW *info)
{     LPPLFE *lfe;
      double sum;
      xassert(1 <= info->p && info->p <= lpp->nrows);
      xassert(lpp->row_stat[info->p] == 0);
      /* the free row is non-active */
      lpp->row_stat[info->p] = LPX_BS;
      /* compute primal value y[p] using the formula (1) */
      sum = 0.0;
      for (lfe = info->ptr; lfe != NULL; lfe = lfe->next)
      {  xassert(1 <= lfe->ref && lfe->ref <= lpp->ncols);
         xassert(lpp->col_stat[lfe->ref] != 0);
         sum += lfe->val * lpp->col_prim[lfe->ref];
      }
      lpp->row_prim[info->p] = sum;
      /* set dual value pi[p] to zero */
      lpp->row_dual[info->p] = 0.0;
      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*
--
-- The fixed column is non-basic.
--
-- Primal value x[q] is a value s[q], at which the column is fixed.
--
-- Dual value can be computed using the dual equality constraint:
--
--    sum a[i,q] * pi[i] + lambda[q] = c[q],                         (6)
--     i
--
-- from which it follows that:
--
--    lambda[q] = c[q] - sum a[i,q] * pi[i],                         (7)
--                        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.
--
-- After recovering the fixed column its primal value s[q] should be
-- also substututed to all rows in order to recover the corresponding
-- primal values (see (2) and (3)):
--
--    y[i] = y'[i] + a[i,q] * s[q].                                  (8)
--
-- Note that the fixed column doesn't affect dual values pi[i], so they
-- are not changed on recovering. */

typedef struct
{     /* fixed column */
      int q;
      /* reference number of a fixed column */
      double val;
      /* value, at which the column is fixed */
      double c;
      /* objective coefficient of the column */
      LPPLFE *ptr;
      /* list of non-zero coefficients a[i,q] */
} FIXED_COL;

static void process_fixed_col(LPP *lpp, LPPCOL *col)
{     FIXED_COL *info;
      LPPROW *row;
      LPPAIJ *aij;
      LPPLFE *lfe;
      /* the column must be fixed */
      xassert(col->lb == col->ub);
      /* create transformation queue entry */
      info = lpp_append_tqe(lpp, LPP_FIXED_COL, sizeof(FIXED_COL));
      info->q = col->j;
      info->val = col->lb;
      info->c = col->c;
      info->ptr = NULL;
      /* save coefficients of the column and update bounds of rows */
      for (aij = col->ptr; aij != NULL; aij = aij->c_next)
      {  row = aij->row;
         lfe = dmp_get_atomv(lpp->tqe_pool, sizeof(LPPLFE));
         lfe->ref = row->i;
         lfe->val = aij->val;
         lfe->next = info->ptr;
         info->ptr = lfe;
         /* see (4) and (5) */
         if (row->lb == row->ub)