/* glpios02.c (preprocess current subproblem) */

/***********************************************************************

*  This code is part of GLPK (GNU Linear Programming Kit).

*

*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,

*  2009, 2010 Andrew Makhorin, Department for Applied Informatics,

*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.

*  E-mail: <mao@gnu.org>.

*

*  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 */

      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 */