glplpx06.c

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

/* glplpx06.c (LP basis and simplex table routines) */

/***********************************************************************
*  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 "glplib.h"
#define xfault xerror
#define lpx_get_b_info glp_get_bhead
#define lpx_get_row_b_ind glp_get_row_bind
#define lpx_get_col_b_ind glp_get_col_bind
#define lpx_ftran glp_ftran
#define lpx_btran glp_btran

/*----------------------------------------------------------------------
-- lpx_eval_b_prim - compute primal basic solution components.
--
-- *Synopsis*
--
-- #include "glplpx.h"
-- void lpx_eval_b_prim(LPX *lp, double row_prim[], double col_prim[]);
--
-- *Description*
--
-- The routine lpx_eval_b_prim computes primal values of all auxiliary
-- and structural variables for the specified problem object.
--
-- NOTE: This routine is intended for internal use only.
--
-- *Background*
--
-- By definition (see glplpx.h):
--
--    B * xB + N * xN = 0,                                           (1)
--
-- where the (basis) matrix B and the matrix N are built of columns of
-- the augmented constraint matrix A~ = (I | -A).
--
-- The formula (1) allows computing primal values of basic variables xB
-- as follows:
--
--    xB = - inv(B) * N * xN,                                        (2)
--
-- where values of non-basic variables xN are defined by corresponding
-- settings for the current basis. */

void lpx_eval_b_prim(LPX *lp, double row_prim[], double col_prim[])
{     int i, j, k, m, n, stat, len, *ind;
      double xN, *NxN, *xB, *val;
      if (!lpx_is_b_avail(lp))
         xfault("lpx_eval_b_prim: LP basis is not available\n");
      m = lpx_get_num_rows(lp);
      n = lpx_get_num_cols(lp);
      /* store values of non-basic auxiliary and structural variables
         and compute the right-hand side vector (-N*xN) */
      NxN = xcalloc(1+m, sizeof(double));
      for (i = 1; i <= m; i++) NxN[i] = 0.0;
      /* walk through auxiliary variables */
      for (i = 1; i <= m; i++)
      {  /* obtain status of i-th auxiliary variable */
         stat = lpx_get_row_stat(lp, i);
         /* if it is basic, skip it */
         if (stat == LPX_BS) continue;
         /* i-th auxiliary variable is non-basic; get its value */
         switch (stat)
         {  case LPX_NL: xN = lpx_get_row_lb(lp, i); break;
            case LPX_NU: xN = lpx_get_row_ub(lp, i); break;
            case LPX_NF: xN = 0.0; break;
            case LPX_NS: xN = lpx_get_row_lb(lp, i); break;
            default: xassert(lp != lp);
         }
         /* store the value of non-basic auxiliary variable */
         row_prim[i] = xN;
         /* and add corresponding term to the right-hand side vector */
         NxN[i] -= xN;
      }
      /* walk through structural variables */
      ind = xcalloc(1+m, sizeof(int));
      val = xcalloc(1+m, sizeof(double));
      for (j = 1; j <= n; j++)
      {  /* obtain status of j-th structural variable */
         stat = lpx_get_col_stat(lp, j);
         /* if it basic, skip it */
         if (stat == LPX_BS) continue;
         /* j-th structural variable is non-basic; get its value */
         switch (stat)
         {  case LPX_NL: xN = lpx_get_col_lb(lp, j); break;
            case LPX_NU: xN = lpx_get_col_ub(lp, j); break;
            case LPX_NF: xN = 0.0; break;
            case LPX_NS: xN = lpx_get_col_lb(lp, j); break;
            default: xassert(lp != lp);
         }
         /* store the value of non-basic structural variable */
         col_prim[j] = xN;
         /* and add corresponding term to the right-hand side vector */
         if (xN != 0.0)
         {  len = lpx_get_mat_col(lp, j, ind, val);
            for (k = 1; k <= len; k++) NxN[ind[k]] += val[k] * xN;
         }
      }
      xfree(ind);
      xfree(val);
      /* solve the system B*xB = (-N*xN) to compute the vector xB */
      xB = NxN, lpx_ftran(lp, xB);
      /* store values of basic auxiliary and structural variables */
      for (i = 1; i <= m; i++)
      {  k = lpx_get_b_info(lp, i);
         xassert(1 <= k && k <= m+n);
         if (k <= m)
            row_prim[k] = xB[i];
         else
            col_prim[k-m] = xB[i];
      }
      xfree(NxN);
      return;
}

/*----------------------------------------------------------------------
-- lpx_eval_b_dual - compute dual basic solution components.
--
-- *Synopsis*
--
-- #include "glplpx.h"
-- void lpx_eval_b_dual(LPX *lp, double row_dual[], double col_dual[]);
--
-- *Description*
--
-- The routine lpx_eval_b_dual computes dual values (that is reduced
-- costs) of all auxiliary and structural variables for the specified
-- problem object.
--
-- NOTE: This routine is intended for internal use only.
--
-- *Background*
--
-- By definition (see glplpx.h):
--
--    B * xB + N * xN = 0,                                           (1)
--
-- where the (basis) matrix B and the matrix N are built of columns of
-- the augmented constraint matrix A~ = (I | -A).
--
-- The objective function can be written as:
--
--    Z = cB' * xB + cN' * xN + c0,                                  (2)
--
-- where cB and cN are objective coefficients at, respectively, basic
-- and non-basic variables.
--
-- From (1) it follows that:
--
--    xB = - inv(B) * N * xN,                                        (3)
--
-- so substituting xB from (3) to (2) we have:
--
--    Z = - cB' * inv(B) * N * xN + cN' * xN + c0 =
--
--      = (cN' - cB' * inv(B) * N) * xN + c0 =
--                                                                   (4)
--      = (cN - N' * inv(B') * cB)' * xN + c0 =
--
--      = d' * xN + c0,
--
-- where
--
--    d = cN - N' * inv(B') * cB                                     (5)
--
-- is the vector of dual values (reduced costs) of non-basic variables.
--
-- The routine first computes the vector pi:
--
--    pi = inv(B') * cB,                                             (6)
--
-- and then computes the vector d as follows:
--
--    d = cN - N' * pi.                                              (7)
--
-- Note that dual values of basic variables are zero by definition. */

void lpx_eval_b_dual(LPX *lp, double row_dual[], double col_dual[])
{     int i, j, k, m, n, len, *ind;
      double dj, *cB, *pi, *val;
      if (!lpx_is_b_avail(lp))
         xfault("lpx_eval_b_dual: LP basis is not available\n");
      m = lpx_get_num_rows(lp);
      n = lpx_get_num_cols(lp);
      /* store zero reduced costs of basic auxiliary and structural
         variables and build the vector cB of objective coefficients at
         basic variables */
      cB = xcalloc(1+m, sizeof(double));
      for (i = 1; i <= m; i++)
      {  k = lpx_get_b_info(lp, i);
         /* xB[i] is k-th original variable */
         xassert(1 <= k && k <= m+n);
         if (k <= m)
         {  row_dual[k] = 0.0;
            cB[i] = 0.0;
         }
         else
         {  col_dual[k-m] = 0.0;
            cB[i] = lpx_get_obj_coef(lp, k-m);
         }
      }
      /* solve the system B'*pi = cB to compute the vector pi */
      pi = cB, lpx_btran(lp, pi);
      /* compute reduced costs of non-basic auxiliary variables */
      for (i = 1; i <= m; i++)
      {  if (lpx_get_row_stat(lp, i) != LPX_BS)
            row_dual[i] = - pi[i];
      }
      /* compute reduced costs of non-basic structural variables */
      ind = xcalloc(1+m, sizeof(int));
      val = xcalloc(1+m, sizeof(double));
      for (j = 1; j <= n; j++)
      {  if (lpx_get_col_stat(lp, j) != LPX_BS)
         {  dj = lpx_get_obj_coef(lp, j);
            len = lpx_get_mat_col(lp, j, ind, val);
            for (k = 1; k <= len; k++) dj += val[k] * pi[ind[k]];
            col_dual[j] = dj;
         }
      }
      xfree(ind);
      xfree(val);
      xfree(cB);
      return;
}

/*----------------------------------------------------------------------
-- lpx_warm_up - "warm up" LP basis.
--
-- *Synopsis*
--
-- #include "glplpx.h"
-- int lpx_warm_up(LPX *lp);
--
-- *Description*
--
-- The routine lpx_warm_up "warms up" the LP basis for the specified
-- problem object using current statuses assigned to rows and columns
-- (i.e. to auxiliary and structural variables).
--
-- "Warming up" includes reinverting (factorizing) the basis matrix (if
-- neccesary), computing primal and dual components of basic solution,
-- and determining primal and dual statuses of the basic solution.
--
-- *Returns*
--
-- The routine lpx_warm_up returns one of the following exit codes:
--
-- LPX_E_OK       the LP basis has been successfully "warmed up".
--
-- LPX_E_EMPTY    the problem has no rows and/or columns.
--
-- LPX_E_BADB     the LP basis is invalid, because the number of basic
--                variables is not the same as the number of rows.
--
-- LPX_E_SING     the basis matrix is singular or ill-conditioned. */

int lpx_warm_up(LPX *lp)
{     int m, n, j, k, ret, type, stat, p_stat, d_stat;
      double lb, ub, prim, dual, tol_bnd, tol_dj, dir;
      double *row_prim, *row_dual, *col_prim, *col_dual, sum;
      m = lpx_get_num_rows(lp);
      n = lpx_get_num_cols(lp);
      /* reinvert the basis matrix, if necessary */
      if (lpx_is_b_avail(lp))
         ret = LPX_E_OK;
      else
      {  if (m == 0 || n == 0)
         {  ret = LPX_E_EMPTY;
            goto done;
         }
#if 0
         ret = lpx_invert(lp);
         switch (ret)
         {  case 0:
               ret = LPX_E_OK;
               break;
            case 1:
            case 2:
               ret = LPX_E_SING;
               goto done;
            case 3:
               ret = LPX_E_BADB;
               goto done;
            default:
               xassert(ret != ret);
         }
#else
         switch (glp_factorize(lp))
         {  case 0:
               ret = LPX_E_OK;
               break;
            case GLP_EBADB:
               ret = LPX_E_BADB;
               goto done;
            case GLP_ESING:
            case GLP_ECOND:
               ret = LPX_E_SING;
               goto done;
            default:
               xassert(lp != lp);
         }
#endif
      }
      /* allocate working arrays */
      row_prim = xcalloc(1+m, sizeof(double));
      row_dual = xcalloc(1+m, sizeof(double));
      col_prim = xcalloc(1+n, sizeof(double));
      col_dual = xcalloc(1+n, sizeof(double));
      /* compute primal basic solution components */
      lpx_eval_b_prim(lp, row_prim, col_prim);
      /* determine primal status of basic solution */
      tol_bnd = 3.0 * lpx_get_real_parm(lp, LPX_K_TOLBND);
      p_stat = LPX_P_FEAS;