glpios02.c

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

      if (U != +DBL_MAX)
      {  double eps = 1e-12 * (1.0 + fabs(U));
         if (UU < U + eps)
         {  /* it cannot be active, so remove it */
            *U_ = +DBL_MAX;
         }
      }
done: return ret;
}

/***********************************************************************
*  check_col_bounds - check and tighten original column 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]
*
*  for column (variable) x[j] this routine checks the original column
*  bounds l[j] and u[j] for feasibility and redundancy. If the original
*  lower bound l[j] or/and upper bound u[j] cannot be active due to
*  bounds of the constraint and other variables, the routine tighten
*  them replacing by corresponding implied bounds, if possible.
*
*  NOTE: It is assumed that if L != -inf, the row lower bound can be
*        active, and if U != +inf, the row upper bound can be active.
*
*  The flag means that variable x[j] is required to be integer.
*
*  New actual bounds for x[j] are stored in locations lj and uj.
*
*  If no primal infeasibility is detected, the routine returns zero,
*  otherwise non-zero. */

static int check_col_bounds(const struct f_info *f, int n,
      const double a[], double L, double U, const double l[],
      const double u[], int flag, int j, double *_lj, double *_uj)
{     int ret = 0;
      double lj, uj, ll, uu;
      xassert(n >= 0);
      xassert(1 <= j && j <= n);
      lj = l[j], uj = u[j];
      /* determine implied bounds of the column */
      col_implied_bounds(f, n, a, L, U, l, u, j, &ll, &uu);
      /* if x[j] is integral, round its implied bounds */
      if (flag)
      {  if (ll != -DBL_MAX)
            ll = (ll - floor(ll) < 1e-3 ? floor(ll) : ceil(ll));
         if (uu != +DBL_MAX)
            uu = (ceil(uu) - uu < 1e-3 ? ceil(uu) : floor(uu));
      }
      /* check if the original lower bound is infeasible */
      if (lj != -DBL_MAX)
      {  double eps = 1e-3 * (1.0 + fabs(lj));
         if (uu < lj - eps)
         {  ret = 1;
            goto done;
         }
      }
      /* check if the original upper bound is infeasible */
      if (uj != +DBL_MAX)
      {  double eps = 1e-3 * (1.0 + fabs(uj));
         if (ll > uj + eps)
         {  ret = 1;
            goto done;
         }
      }
      /* check if the original lower bound is redundant */
      if (ll != -DBL_MAX)
      {  double eps = 1e-3 * (1.0 + fabs(ll));
         if (lj < ll - eps)
         {  /* it cannot be active, so tighten it */
            lj = ll;
         }
      }
      /* check if the original upper bound is redundant */
      if (uu != +DBL_MAX)
      {  double eps = 1e-3 * (1.0 + fabs(uu));
         if (uj > uu + eps)
         {  /* it cannot be active, so tighten it */
            uj = uu;
         }
      }
      /* due to round-off errors it may happen that lj > uj (although
         lj < uj + eps, since no primal infeasibility is detected), so
         adjuct the new actual bounds to provide lj <= uj */
      if (!(lj == -DBL_MAX || uj == +DBL_MAX))
      {  double t1 = fabs(lj), t2 = fabs(uj);
         double eps = 1e-10 * (1.0 + (t1 <= t2 ? t1 : t2));
         if (lj > uj - eps)
         {  if (lj == l[j])
               uj = lj;
            else if (uj == u[j])
               lj = uj;
            else if (t1 <= t2)
               uj = lj;
            else
               lj = uj;
         }
      }
      *_lj = lj, *_uj = uj;
done: return ret;
}

/***********************************************************************
*  check_efficiency - check if change in column bounds is efficient
*
*  Given the original bounds of a column l and u and its new actual
*  bounds l' and u' (possibly tighten by the routine check_col_bounds)
*  this routine checks if the change in the column bounds is efficient
*  enough. If so, the routine returns non-zero, otherwise zero.
*
*  The flag means that the variable is required to be integer. */

static int check_efficiency(int flag, double l, double u, double ll,
      double uu)
{     int eff = 0;
      /* check efficiency for lower bound */
      if (l < ll)
      {  if (flag || l == -DBL_MAX)
            eff++;
         else
         {  double r;
            if (u == +DBL_MAX)
               r = 1.0 + fabs(l);
            else
               r = 1.0 + (u - l);
            if (ll - l >= 0.25 * r)
               eff++;
         }
      }
      /* check efficiency for upper bound */
      if (u > uu)
      {  if (flag || u == +DBL_MAX)
            eff++;
         else
         {  double r;
            if (l == -DBL_MAX)
               r = 1.0 + fabs(u);
            else
               r = 1.0 + (u - l);
            if (u - uu >= 0.25 * r)
               eff++;
         }
      }
      return eff;
}

/***********************************************************************
*  basic_preprocessing - perform basic preprocessing
*
*  This routine performs basic preprocessing of the specified MIP that
*  includes relaxing some row bounds and tightening some column bounds.
*
*  On entry the arrays L and U contains original row bounds, and the
*  arrays l and u contains original column bounds:
*
*  L[0] is the lower bound of the objective row;
*  L[i], i = 1,...,m, is the lower bound of i-th row;
*  U[0] is the upper bound of the objective row;
*  U[i], i = 1,...,m, is the upper bound of i-th row;
*  l[0] is not used;
*  l[j], j = 1,...,n, is the lower bound of j-th column;
*  u[0] is not used;
*  u[j], j = 1,...,n, is the upper bound of j-th column.
*
*  On exit the arrays L, U, l, and u contain new actual bounds of rows
*  and column in the same locations.
*
*  The parameters nrs and num specify an initial list of rows to be
*  processed:
*
*  nrs is the number of rows in the initial list, 0 <= nrs <= m+1;
*  num[0] is not used;
*  num[1,...,nrs] are row numbers (0 means the objective row).
*
*  The parameter max_pass specifies the maximal number of times that
*  each row can be processed, max_pass > 0.
*
*  If no primal infeasibility is detected, the routine returns zero,
*  otherwise non-zero. */

static int basic_preprocessing(glp_prob *mip, double L[], double U[],
      double l[], double u[], int nrs, const int num[], int max_pass)
{     int m = mip->m;
      int n = mip->n;
      struct f_info f;
      int i, j, k, len, size, ret = 0;
      int *ind, *list, *mark, *pass;
      double *val, *lb, *ub;
      xassert(0 <= nrs && nrs <= m+1);
      xassert(max_pass > 0);
      /* allocate working arrays */
      ind = xcalloc(1+n, sizeof(int));
      list = xcalloc(1+m+1, sizeof(int));
      mark = xcalloc(1+m+1, sizeof(int));
      memset(&mark[0], 0, (m+1) * sizeof(int));
      pass = xcalloc(1+m+1, sizeof(int));
      memset(&pass[0], 0, (m+1) * sizeof(int));
      val = xcalloc(1+n, sizeof(double));
      lb = xcalloc(1+n, sizeof(double));
      ub = xcalloc(1+n, sizeof(double));
      /* initialize the list of rows to be processed */
      size = 0;
      for (k = 1; k <= nrs; k++)
      {  i = num[k];
         xassert(0 <= i && i <= m);
         /* duplicate row numbers are not allowed */
         xassert(!mark[i]);
         list[++size] = i, mark[i] = 1;
      }
      xassert(size == nrs);
      /* process rows in the list until it becomes empty */
      while (size > 0)
      {  /* get a next row from the list */
         i = list[size--], mark[i] = 0;
         /* increase the row processing count */
         pass[i]++;
         /* if the row is free, skip it */
         if (L[i] == -DBL_MAX && U[i] == +DBL_MAX) continue;
         /* obtain coefficients of the row */
         len = 0;
         if (i == 0)
         {  for (j = 1; j <= n; j++)
            {  GLPCOL *col = mip->col[j];
               if (col->coef != 0.0)
                  len++, ind[len] = j, val[len] = col->coef;
            }
         }
         else
         {  GLPROW *row = mip->row[i];
            GLPAIJ *aij;
            for (aij = row->ptr; aij != NULL; aij = aij->r_next)
               len++, ind[len] = aij->col->j, val[len] = aij->val;
         }
         /* determine lower and upper bounds of columns corresponding
            to non-zero row coefficients */
         for (k = 1; k <= len; k++)
            j = ind[k], lb[k] = l[j], ub[k] = u[j];
         /* prepare the row info to determine implied bounds */
         prepare_row_info(len, val, lb, ub, &f);
         /* check and relax bounds of the row */
         if (check_row_bounds(&f, &L[i], &U[i]))
         {  /* the feasible region is empty */
            ret = 1;
            goto done;
         }
         /* if the row became free, drop it */
         if (L[i] == -DBL_MAX && U[i] == +DBL_MAX) continue;
         /* process columns having non-zero coefficients in the row */
         for (k = 1; k <= len; k++)
         {  GLPCOL *col;
            int flag, eff;
            double ll, uu;
            /* take a next column in the row */
            j = ind[k], col = mip->col[j];
            flag = col->kind != GLP_CV;
            /* check and tighten bounds of the column */
            if (check_col_bounds(&f, len, val, L[i], U[i], lb, ub,
                flag, k, &ll, &uu))
            {  /* the feasible region is empty */
               ret = 1;
               goto done;
            }
            /* check if change in the column bounds is efficient */
            eff = check_efficiency(flag, l[j], u[j], ll, uu);
            /* set new actual bounds of the column */
            l[j] = ll, u[j] = uu;
            /* if the change is efficient, add all rows affected by the
               corresponding column, to the list */
            if (eff > 0)
            {  GLPAIJ *aij;
               for (aij = col->ptr; aij != NULL; aij = aij->c_next)
               {  int ii = aij->row->i;
                  /* if the row was processed maximal number of times,
                     skip it */
                  if (pass[ii] >= max_pass) continue;
                  /* if the row is free, skip it */
                  if (L[ii] == -DBL_MAX && U[ii] == +DBL_MAX) continue;
                  /* put the row into the list */
                  if (mark[ii] == 0)
                  {  xassert(size <= m);
                     list[++size] = ii, mark[ii] = 1;
                  }
               }
            }
         }
      }
done: /* free working arrays */
      xfree(ind);
      xfree(list);
      xfree(mark);
      xfree(pass);
      xfree(val);
      xfree(lb);
      xfree(ub);
      return ret;
}

/***********************************************************************
*  NAME
*
*  ios_preprocess_node - preprocess current subproblem
*
*  SYNOPSIS
*
*  #include "glpios.h"
*  int ios_preprocess_node(glp_tree *tree, int max_pass);
*
*  DESCRIPTION
*
*  The routine ios_preprocess_node performs basic preprocessing of the
*  current subproblem.
*
*  RETURNS
*
*  If no primal infeasibility is detected, the routine returns zero,
*  otherwise non-zero. */

int ios_preprocess_node(glp_tree *tree, int max_pass)
{     glp_prob *mip = tree->mip;
      int m = mip->m;
      int n = mip->n;
      int i, j, nrs, *num, ret = 0;
      double *L, *U, *l, *u;
      /* the current subproblem must exist */
      xassert(tree->curr != NULL);
      /* determine original row bounds */
      L = xcalloc(1+m, sizeof(double));
      U = xcalloc(1+m, sizeof(double));
      switch (mip->mip_stat)
      {  case GLP_UNDEF:
            L[0] = -DBL_MAX, U[0] = +DBL_MAX;
            break;
         case GLP_FEAS:
            switch (mip->dir)
            {  case GLP_MIN:
                  L[0] = -DBL_MAX, U[0] = mip->mip_obj - mip->c0;
                  break;
               case GLP_MAX:
                  L[0] = mip->mip_obj - mip->c0, U[0] = +DBL_MAX;
                  break;
               default:
                  xassert(mip != mip);
            }
            break;
         default:
            xassert(mip != mip);
      }
      for (i = 1; i <= m; i++)
      {  L[i] = glp_get_row_lb(mip, i);
         U[i] = glp_get_row_ub(mip, i);
      }
      /* determine original column bounds */
      l = xcalloc(1+n, sizeof(double));
      u = xcalloc(1+n, sizeof(double));
      for (j = 1; j <= n; j++)
      {  l[j] = glp_get_col_lb(mip, j);
         u[j] = glp_get_col_ub(mip, j);
      }
      /* build the initial list of rows to be analyzed */
      nrs = m + 1;
      num = xcalloc(1+nrs, sizeof(int));
      for (i = 1; i <= nrs; i++) num[i] = i - 1;
      /* perform basic preprocessing */
      if (basic_preprocessing(mip , L, U, l, u, nrs, num, max_pass))
      {  ret = 1;
         goto done;
      }
      /* set new actual (relaxed) row bounds */
      for (i = 1; i <= m; i++)
      {  /* consider only non-active rows to keep dual feasibility */
         if (glp_get_row_stat(mip, i) == GLP_BS)
         {  if (L[i] == -DBL_MAX && U[i] == +DBL_MAX)
               glp_set_row_bnds(mip, i, GLP_FR, 0.0, 0.0);
            else if (U[i] == +DBL_MAX)
               glp_set_row_bnds(mip, i, GLP_LO, L[i], 0.0);
            else if (L[i] == -DBL_MAX)
               glp_set_row_bnds(mip, i, GLP_UP, 0.0, U[i]);
         }
      }
      /* set new actual (tightened) column bounds */
      for (j = 1; j <= n; j++)
      {  int type;
         if (l[j] == -DBL_MAX && u[j] == +DBL_MAX)
            type = GLP_FR;
         else if (u[j] == +DBL_MAX)
            type = GLP_LO;
         else if (l[j] == -DBL_MAX)
            type = GLP_UP;
         else if (l[j] != u[j])
            type = GLP_DB;
         else
            type = GLP_FX;
         glp_set_col_bnds(mip, j, type, l[j], u[j]);
      }
done: /* free working arrays and return */
      xfree(L);
      xfree(U);
      xfree(l);
      xfree(u);
      xfree(num);
      return ret;
}

/* eof */