glplpx06.c

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

      for (k = 1; k <= m+n; k++)
      {  if (k <= m)
         {  type = lpx_get_row_type(lp, k);
            lb = lpx_get_row_lb(lp, k);
            ub = lpx_get_row_ub(lp, k);
            prim = row_prim[k];
         }
         else
         {  type = lpx_get_col_type(lp, k-m);
            lb = lpx_get_col_lb(lp, k-m);
            ub = lpx_get_col_ub(lp, k-m);
            prim = col_prim[k-m];
         }
         if (type == LPX_LO || type == LPX_DB || type == LPX_FX)
         {  /* variable x[k] has lower bound */
            if (prim < lb - tol_bnd * (1.0 + fabs(lb)))
            {  p_stat = LPX_P_INFEAS;
               break;
            }
         }
         if (type == LPX_UP || type == LPX_DB || type == LPX_FX)
         {  /* variable x[k] has upper bound */
            if (prim > ub + tol_bnd * (1.0 + fabs(ub)))
            {  p_stat = LPX_P_INFEAS;
               break;
            }
         }
      }
      /* compute dual basic solution components */
      lpx_eval_b_dual(lp, row_dual, col_dual);
      /* determine dual status of basic solution */
      tol_dj = 3.0 * lpx_get_real_parm(lp, LPX_K_TOLDJ);
      dir = (lpx_get_obj_dir(lp) == LPX_MIN ? +1.0 : -1.0);
      d_stat = LPX_D_FEAS;
      for (k = 1; k <= m+n; k++)
      {  if (k <= m)
         {  stat = lpx_get_row_stat(lp, k);
            dual = row_dual[k];
         }
         else
         {  stat = lpx_get_col_stat(lp, k-m);
            dual = col_dual[k-m];
         }
         if (stat == LPX_BS || stat == LPX_NL || stat == LPX_NF)
         {  /* reduced cost of x[k] must be non-negative (minimization)
               or non-positive (maximization) */
            if (dir * dual < - tol_dj)
            {  d_stat = LPX_D_INFEAS;
               break;
            }
         }
         if (stat == LPX_BS || stat == LPX_NU || stat == LPX_NF)
         {  /* reduced cost of x[k] must be non-positive (minimization)
               or non-negative (maximization) */
            if (dir * dual > + tol_dj)
            {  d_stat = LPX_D_INFEAS;
               break;
            }
         }
      }
      /* store basic solution components */
      p_stat = p_stat - LPX_P_UNDEF + GLP_UNDEF;
      d_stat = d_stat - LPX_D_UNDEF + GLP_UNDEF;
      sum = lpx_get_obj_coef(lp, 0);
      for (j = 1; j <= n; j++)
         sum += lpx_get_obj_coef(lp, j) * col_prim[j];
      lpx_put_solution(lp, 0, &p_stat, &d_stat, &sum,
         NULL, row_prim, row_dual, NULL, col_prim, col_dual);
      xassert(lpx_is_b_avail(lp));
      /* free working arrays */
      xfree(row_prim);
      xfree(row_dual);
      xfree(col_prim);
      xfree(col_dual);
done: /* return to the calling program */
      return ret;
}

/*----------------------------------------------------------------------
-- lpx_transform_row - transform explicitly specified row.
--
-- *Synopsis*
--
-- #include "glplpx.h"
-- int lpx_transform_row(LPX *lp, int len, int ind[], double val[]);
--
-- *Description*
--
-- The routine lpx_transform_row performs the same operation as the
-- routine lpx_eval_tab_row with exception that the transformed row is
-- specified explicitly as a sparse vector.
--
-- The explicitly specified row may be thought as a linear form:
--
--    x = a[1]*x[m+1] + a[2]*x[m+2] + ... + a[n]*x[m+n],             (1)
--
-- where x is an auxiliary variable for this row, a[j] are coefficients
-- of the linear form, x[m+j] are structural variables.
--
-- On entry column indices and numerical values of non-zero elements of
-- the row should be stored in locations ind[1], ..., ind[len] and
-- val[1], ..., val[len], where len is the number of non-zero elements.
--
-- This routine uses the system of equality constraints and the current
-- basis in order to express the auxiliary variable x in (1) through the
-- current non-basic variables (as if the transformed row were added to
-- the problem object and its auxiliary variable were basic), i.e. the
-- resultant row has the form:
--
--    x = alfa[1]*xN[1] + alfa[2]*xN[2] + ... + alfa[n]*xN[n],       (2)
--
-- where xN[j] are non-basic (auxiliary or structural) variables, n is
-- the number of columns in the LP problem object.
--
-- On exit the routine stores indices and numerical values of non-zero
-- elements of the resultant row (2) in locations ind[1], ..., ind[len']
-- and val[1], ..., val[len'], where 0 <= len' <= n is the number of
-- non-zero elements in the resultant row returned by the routine. Note
-- that indices (numbers) of non-basic variables stored in the array ind
-- correspond to original ordinal numbers of variables: indices 1 to m
-- mean auxiliary variables and indices m+1 to m+n mean structural ones.
--
-- *Returns*
--
-- The routine returns len', which is the number of non-zero elements in
-- the resultant row stored in the arrays ind and val.
--
-- *Background*
--
-- The explicitly specified row (1) is transformed in the same way as
-- it were the objective function row.
--
-- From (1) it follows that:
--
--    x = aB * xB + aN * xN,                                         (3)
--
-- where xB is the vector of basic variables, xN is the vector of
-- non-basic variables.
--
-- The simplex table, which corresponds to the current basis, is:
--
--    xB = [-inv(B) * N] * xN.                                       (4)
--
-- Therefore substituting xB from (4) to (3) we have:
--
--    x = aB * [-inv(B) * N] * xN + aN * xN =
--                                                                   (5)
--      = rho * (-N) * xN + aN * xN = alfa * xN,
--
-- where:
--
--    rho = inv(B') * aB,                                            (6)
--
-- and
--
--    alfa = aN + rho * (-N)                                         (7)
--
-- is the resultant row computed by the routine. */

int lpx_transform_row(LPX *lp, int len, int ind[], double val[])
{     int i, j, k, m, n, t, lll, *iii;
      double alfa, *a, *aB, *rho, *vvv;
      if (!lpx_is_b_avail(lp))
         xfault("lpx_transform_row: LP basis is not available\n");
      m = lpx_get_num_rows(lp);
      n = lpx_get_num_cols(lp);
      /* unpack the row to be transformed to the array a */
      a = xcalloc(1+n, sizeof(double));
      for (j = 1; j <= n; j++) a[j] = 0.0;
      if (!(0 <= len && len <= n))
         xfault("lpx_transform_row: len = %d; invalid row length\n",
            len);
      for (t = 1; t <= len; t++)
      {  j = ind[t];
         if (!(1 <= j && j <= n))
            xfault("lpx_transform_row: ind[%d] = %d; column index out o"
               "f range\n", t, j);
         if (val[t] == 0.0)
            xfault("lpx_transform_row: val[%d] = 0; zero coefficient no"
               "t allowed\n", t);
         if (a[j] != 0.0)
            xfault("lpx_transform_row: ind[%d] = %d; duplicate column i"
               "ndices not allowed\n", t, j);
         a[j] = val[t];
      }
      /* construct the vector aB */
      aB = 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);
         aB[i] = (k <= m ? 0.0 : a[k-m]);
      }
      /* solve the system B'*rho = aB to compute the vector rho */
      rho = aB, lpx_btran(lp, rho);
      /* compute coefficients at non-basic auxiliary variables */
      len = 0;
      for (i = 1; i <= m; i++)
      {  if (lpx_get_row_stat(lp, i) != LPX_BS)
         {  alfa = - rho[i];
            if (alfa != 0.0)
            {  len++;
               ind[len] = i;
               val[len] = alfa;
            }
         }
      }
      /* compute coefficients at non-basic structural variables */
      iii = xcalloc(1+m, sizeof(int));
      vvv = xcalloc(1+m, sizeof(double));
      for (j = 1; j <= n; j++)
      {  if (lpx_get_col_stat(lp, j) != LPX_BS)
         {  alfa = a[j];
            lll = lpx_get_mat_col(lp, j, iii, vvv);
            for (t = 1; t <= lll; t++) alfa += vvv[t] * rho[iii[t]];
            if (alfa != 0.0)
            {  len++;
               ind[len] = m+j;
               val[len] = alfa;
            }
         }
      }
      xassert(len <= n);
      xfree(iii);
      xfree(vvv);
      xfree(aB);
      xfree(a);
      return len;
}

/*----------------------------------------------------------------------
-- lpx_transform_col - transform explicitly specified column.
--
-- *Synopsis*
--
-- #include "glplpx.h"
-- int lpx_transform_col(LPX *lp, int len, int ind[], double val[]);
--
-- *Description*
--
-- The routine lpx_transform_col performs the same operation as the
-- routine lpx_eval_tab_col with exception that the transformed column
-- is specified explicitly as a sparse vector.
--
-- The explicitly specified column may be thought as if it were added
-- to the original system of equality constraints:
--
--    x[1] = a[1,1]*x[m+1] + ... + a[1,n]*x[m+n] + a[1]*x
--    x[2] = a[2,1]*x[m+1] + ... + a[2,n]*x[m+n] + a[2]*x            (1)
--       .  .  .  .  .  .  .  .  .  .  .  .  .  .  .
--    x[m] = a[m,1]*x[m+1] + ... + a[m,n]*x[m+n] + a[m]*x
--
-- where x[i] are auxiliary variables, x[m+j] are structural variables,
-- x is a structural variable for the explicitly specified column, a[i]
-- are constraint coefficients for x.
--
-- On entry row indices and numerical values of non-zero elements of
-- the column should be stored in locations ind[1], ..., ind[len] and
-- val[1], ..., val[len], where len is the number of non-zero elements.
--
-- This routine uses the system of equality constraints and the current
-- basis in order to express the current basic variables through the
-- structural variable x in (1) (as if the transformed column were added
-- to the problem object and the variable x were non-basic), i.e. the
-- resultant column has the form:
--
--    xB[1] = ... + alfa[1]*x
--    xB[2] = ... + alfa[2]*x                                        (2)
--       .  .  .  .  .  .
--    xB[m] = ... + alfa[m]*x
--
-- where xB are basic (auxiliary and structural) variables, m is the
-- number of rows in the problem object.
--
-- On exit the routine stores indices and numerical values of non-zero
-- elements of the resultant column (2) in locations ind[1], ...,
-- ind[len'] and val[1], ..., val[len'], where 0 <= len' <= m is the
-- number of non-zero element in the resultant column returned by the
-- routine. Note that indices (numbers) of basic variables stored in
-- the array ind correspond to original ordinal numbers of variables:
-- indices 1 to m mean auxiliary variables and indices m+1 to m+n mean
-- structural ones.
--
-- *Returns*
--
-- The routine returns len', which is the number of non-zero elements
-- in the resultant column stored in the arrays ind and val.
--
-- *Background*
--
-- The explicitly specified column (1) is transformed in the same way
-- as any other column of the constraint matrix using the formula:
--
--    alfa = inv(B) * a,                                             (3)
--
-- where alfa is the resultant column computed by the routine. */

int lpx_transform_col(LPX *lp, int len, int ind[], double val[])
{     int i, m, t;
      double *a, *alfa;
      if (!lpx_is_b_avail(lp))
         xfault("lpx_transform_col: LP basis is not available\n");
      m = lpx_get_num_rows(lp);
      /* unpack the column to be transformed to the array a */
      a = xcalloc(1+m, sizeof(double));
      for (i = 1; i <= m; i++) a[i] = 0.0;
      if (!(0 <= len && len <= m))
         xfault("lpx_transform_col: len = %d; invalid column length\n",
            len);
      for (t = 1; t <= len; t++)
      {  i = ind[t];
         if (!(1 <= i && i <= m))
            xfault("lpx_transform_col: ind[%d] = %d; row index out of r"
               "ange\n", t, i);
         if (val[t] == 0.0)
            xfault("lpx_transform_col: val[%d] = 0; zero coefficient no"
               "t allowed\n", t);
         if (a[i] != 0.0)
            xfault("lpx_transform_col: ind[%d] = %d; duplicate row indi"
               "ces not allowed\n", t, i);
         a[i] = val[t];
      }
      /* solve the system B*a = alfa to compute the vector alfa */
      alfa = a, lpx_ftran(lp, alfa);
      /* store resultant coefficients */
      len = 0;
      for (i = 1; i <= m; i++)
      {  if (alfa[i] != 0.0)
         {  len++;
            ind[len] = lpx_get_b_info(lp, i);
            val[len] = alfa[i];
         }
      }