glpios02.c

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

/* glpios02.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 "glpios.h"

/***********************************************************************
*  prepare_row_info - prepare row info to determine implied bounds
*
*  Given a row (linear form)
*
*      n
*     sum a[j] * x[j]                                                (1)
*     j=1
*
*  and bounds of columns (variables)
*
*     l[j] <= x[j] <= u[j]                                           (2)
*
*  this routine computes f_min, j_min, f_max, j_max needed to determine
*  implied bounds.
*
*  ALGORITHM
*
*  Let J+ = {j : a[j] > 0} and J- = {j : a[j] < 0}.
*
*  Parameters f_min and j_min are computed as follows:
*
*  1) if there is no x[k] such that k in J+ and l[k] = -inf or k in J-
*     and u[k] = +inf, then
*
*     f_min :=   sum   a[j] * l[j] +   sum   a[j] * u[j]
*              j in J+               j in J-
*                                                                    (3)
*     j_min := 0
*
*  2) if there is exactly one x[k] such that k in J+ and l[k] = -inf
*     or k in J- and u[k] = +inf, then
*
*     f_min :=   sum       a[j] * l[j] +   sum       a[j] * u[j]
*              j in J+\{k}               j in J-\{k}
*                                                                    (4)
*     j_min := k
*
*  3) if there are two or more x[k] such that k in J+ and l[k] = -inf
*     or k in J- and u[k] = +inf, then
*
*     f_min := -inf
*                                                                    (5)
*     j_min := 0
*
*  Parameters f_max and j_max are computed in a similar way as follows:
*
*  1) if there is no x[k] such that k in J+ and u[k] = +inf or k in J-
*     and l[k] = -inf, then
*
*     f_max :=   sum   a[j] * u[j] +   sum   a[j] * l[j]
*              j in J+               j in J-
*                                                                    (6)
*     j_max := 0
*
*  2) if there is exactly one x[k] such that k in J+ and u[k] = +inf
*     or k in J- and l[k] = -inf, then
*
*     f_max :=   sum       a[j] * u[j] +   sum       a[j] * l[j]
*              j in J+\{k}               j in J-\{k}
*                                                                    (7)
*     j_max := k
*
*  3) if there are two or more x[k] such that k in J+ and u[k] = +inf
*     or k in J- and l[k] = -inf, then
*
*     f_max := +inf
*                                                                    (8)
*     j_max := 0                                                      */

struct f_info
{     int j_min, j_max;
      double f_min, f_max;
};

static void prepare_row_info(int n, const double a[], const double l[],
      const double u[], struct f_info *f)
{     int j, j_min, j_max;
      double f_min, f_max;
      xassert(n >= 0);
      /* determine f_min and j_min */
      f_min = 0.0, j_min = 0;
      for (j = 1; j <= n; j++)
      {  if (a[j] > 0.0)
         {  if (l[j] == -DBL_MAX)
            {  if (j_min == 0)
                  j_min = j;
               else
               {  f_min = -DBL_MAX, j_min = 0;
                  break;
               }
            }
            else
               f_min += a[j] * l[j];
         }
         else if (a[j] < 0.0)
         {  if (u[j] == +DBL_MAX)
            {  if (j_min == 0)
                  j_min = j;
               else
               {  f_min = -DBL_MAX, j_min = 0;
                  break;
               }
            }
            else
               f_min += a[j] * u[j];
         }
         else
            xassert(a != a);
      }
      f->f_min = f_min, f->j_min = j_min;
      /* determine f_max and j_max */
      f_max = 0.0, j_max = 0;
      for (j = 1; j <= n; j++)
      {  if (a[j] > 0.0)
         {  if (u[j] == +DBL_MAX)
            {  if (j_max == 0)
                  j_max = j;
               else
               {  f_max = +DBL_MAX, j_max = 0;
                  break;
               }
            }
            else
               f_max += a[j] * u[j];
         }
         else if (a[j] < 0.0)
         {  if (l[j] == -DBL_MAX)
            {  if (j_max == 0)
                  j_max = j;
               else
               {  f_max = +DBL_MAX, j_max = 0;
                  break;
               }
            }
            else
               f_max += a[j] * l[j];
         }
         else
            xassert(a != a);
      }
      f->f_max = f_max, f->j_max = j_max;
      return;
}

/***********************************************************************
*  row_implied_bounds - determine row implied bounds
*
*  Given a row (linear form)
*
*      n
*     sum a[j] * x[j]
*     j=1
*
*  and bounds of columns (variables)
*
*     l[j] <= x[j] <= u[j]
*
*  this routine determines implied bounds of the row.
*
*  ALGORITHM
*
*  Let J+ = {j : a[j] > 0} and J- = {j : a[j] < 0}.
*
*  The implied lower bound of the row is computed as follows:
*
*     L' :=   sum   a[j] * l[j] +   sum   a[j] * u[j]                (9)
*           j in J+               j in J-
*
*  and as it follows from (3), (4), and (5):
*
*     L' := if j_min = 0 then f_min else -inf                       (10)
*
*  The implied upper bound of the row is computed as follows:
*
*     U' :=   sum   a[j] * u[j] +   sum   a[j] * l[j]               (11)
*           j in J+               j in J-
*
*  and as it follows from (6), (7), and (8):
*
*     U' := if j_max = 0 then f_max else +inf                       (12)
*
*  The implied bounds are stored in locations LL and UU. */

static void row_implied_bounds(const struct f_info *f, double *LL,
      double *UU)
{     *LL = (f->j_min == 0 ? f->f_min : -DBL_MAX);
      *UU = (f->j_max == 0 ? f->f_max : +DBL_MAX);
      return;
}

/***********************************************************************
*  col_implied_bounds - determine column implied bounds
*
*  Given a row (constraint)
*
*           n
*     L <= sum a[j] * x[j] <= U                                     (13)
*          j=1
*
*  and bounds of columns (variables)
*
*     l[j] <= x[j] <= u[j]
*
*  this routine determines implied bounds of variable x[k].
*
*  It is assumed that if L != -inf, the lower bound of the row can be
*  active, and if U != +inf, the upper bound of the row can be active.
*
*  ALGORITHM
*
*  From (13) it follows that
*
*     L <= sum a[j] * x[j] + a[k] * x[k] <= U
*          j!=k
*  or
*
*     L - sum a[j] * x[j] <= a[k] * x[k] <= U - sum a[j] * x[j]
*         j!=k                                  j!=k
*
*  Thus, if the row lower bound L can be active, implied lower bound of
*  term a[k] * x[k] can be determined as follows:
*
*     ilb(a[k] * x[k]) = min(L - sum a[j] * x[j]) =
*                                j!=k
*                                                                   (14)
*                      = L - max sum a[j] * x[j]
*                            j!=k
*
*  where, as it follows from (6), (7), and (8)
*
*                           / f_max - a[k] * u[k], j_max = 0, a[k] > 0
*                           |
*                           | f_max - a[k] * l[k], j_max = 0, a[k] < 0
*     max sum a[j] * x[j] = {
*         j!=k              | f_max,               j_max = k
*                           |
*                           \ +inf,                j_max != 0
*
*  and if the upper bound U can be active, implied upper bound of term
*  a[k] * x[k] can be determined as follows:
*
*     iub(a[k] * x[k]) = max(U - sum a[j] * x[j]) =
*                                j!=k
*                                                                   (15)
*                      = U - min sum a[j] * x[j]
*                            j!=k
*
*  where, as it follows from (3), (4), and (5)
*
*                           / f_min - a[k] * l[k], j_min = 0, a[k] > 0
*                           |
*                           | f_min - a[k] * u[k], j_min = 0, a[k] < 0
*     min sum a[j] * x[j] = {
*         j!=k              | f_min,               j_min = k
*                           |
*                           \ -inf,                j_min != 0
*
*  Since
*
*     ilb(a[k] * x[k]) <= a[k] * x[k] <= iub(a[k] * x[k])
*
*  implied lower and upper bounds of x[k] are determined as follows:
*
*     l'[k] := if a[k] > 0 then ilb / a[k] else ulb / a[k]          (16)
*
*     u'[k] := if a[k] > 0 then ulb / a[k] else ilb / a[k]          (17)
*
*  The implied bounds are stored in locations ll and uu. */

static void col_implied_bounds(const struct f_info *f, int n,
      const double a[], double L, double U, const double l[],
      const double u[], int k, double *ll, double *uu)
{     double ilb, iub;
      xassert(n >= 0);
      xassert(1 <= k && k <= n);
      /* determine implied lower bound of term a[k] * x[k] (14) */
      if (L == -DBL_MAX || f->f_max == +DBL_MAX)
         ilb = -DBL_MAX;
      else if (f->j_max == 0)
      {  if (a[k] > 0.0)
         {  xassert(u[k] != +DBL_MAX);
            ilb = L - (f->f_max - a[k] * u[k]);
         }
         else if (a[k] < 0.0)
         {  xassert(l[k] != -DBL_MAX);
            ilb = L - (f->f_max - a[k] * l[k]);
         }
         else
            xassert(a != a);
      }
      else if (f->j_max == k)
         ilb = L - f->f_max;
      else
         ilb = -DBL_MAX;
      /* determine implied upper bound of term a[k] * x[k] (15) */
      if (U == +DBL_MAX || f->f_min == -DBL_MAX)
         iub = +DBL_MAX;
      else if (f->j_min == 0)
      {  if (a[k] > 0.0)
         {  xassert(l[k] != -DBL_MAX);
            iub = U - (f->f_min - a[k] * l[k]);
         }
         else if (a[k] < 0.0)
         {  xassert(u[k] != +DBL_MAX);
            iub = U - (f->f_min - a[k] * u[k]);
         }
         else
            xassert(a != a);
      }
      else if (f->j_min == k)
         iub = U - f->f_min;
      else
         iub = +DBL_MAX;
      /* determine implied bounds of x[k] (16) and (17) */
#if 1
      /* do not use a[k] if it has small magnitude to prevent wrong
         implied bounds; for example, 1e-15 * x1 >= x2 + x3, where
         x1 >= -10, x2, x3 >= 0, would lead to wrong conclusion that
         x1 >= 0 */
      if (fabs(a[k]) < 1e-6)
         *ll = -DBL_MAX, *uu = +DBL_MAX; else
#endif
      if (a[k] > 0.0)
      {  *ll = (ilb == -DBL_MAX ? -DBL_MAX : ilb / a[k]);
         *uu = (iub == +DBL_MAX ? +DBL_MAX : iub / a[k]);
      }
      else if (a[k] < 0.0)
      {  *ll = (iub == +DBL_MAX ? -DBL_MAX : iub / a[k]);
         *uu = (ilb == -DBL_MAX ? +DBL_MAX : ilb / a[k]);
      }
      else
         xassert(a != a);
      return;
}

/***********************************************************************
*  check_row_bounds - check and relax original row bounds
*
*  Given a row (constraint)
*
*           n
*     L <= sum a[j] * x[j] <= U
*          j=1
*
*  and bounds of columns (variables)
*
*     l[j] <= x[j] <= u[j]
*
*  this routine checks the original row bounds L and U for feasibility
*  and redundancy. If the original lower bound L or/and upper bound U
*  cannot be active due to bounds of variables, the routine remove them
*  replacing by -inf or/and +inf, respectively.
*
*  If no primal infeasibility is detected, the routine returns zero,
*  otherwise non-zero. */

static int check_row_bounds(const struct f_info *f, double *L_,
      double *U_)
{     int ret = 0;
      double L = *L_, U = *U_, LL, UU;
      /* determine implied bounds of the row */
      row_implied_bounds(f, &LL, &UU);
      /* check if the original lower bound is infeasible */
      if (L != -DBL_MAX)
      {  double eps = 1e-3 * (1.0 + fabs(L));
         if (UU < L - eps)
         {  ret = 1;
            goto done;
         }
      }
      /* check if the original upper bound is infeasible */
      if (U != +DBL_MAX)
      {  double eps = 1e-3 * (1.0 + fabs(U));
         if (LL > U + eps)
         {  ret = 1;
            goto done;
         }
      }
      /* check if the original lower bound is redundant */
      if (L != -DBL_MAX)
      {  double eps = 1e-12 * (1.0 + fabs(L));
         if (LL > L - eps)
         {  /* it cannot be active, so remove it */
            *L_ = -DBL_MAX;
         }
      }
      /* check if the original upper bound is redundant */