src/glpios02.c
author Alpar Juttner <alpar@cs.elte.hu>
Sun, 05 Dec 2010 17:35:23 +0100
changeset 2 4c8956a7bdf4
permissions -rw-r--r--
Set up CMAKE build environment
     1 /* glpios02.c (preprocess current subproblem) */
     2 
     3 /***********************************************************************
     4 *  This code is part of GLPK (GNU Linear Programming Kit).
     5 *
     6 *  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
     7 *  2009, 2010 Andrew Makhorin, Department for Applied Informatics,
     8 *  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
     9 *  E-mail: <mao@gnu.org>.
    10 *
    11 *  GLPK is free software: you can redistribute it and/or modify it
    12 *  under the terms of the GNU General Public License as published by
    13 *  the Free Software Foundation, either version 3 of the License, or
    14 *  (at your option) any later version.
    15 *
    16 *  GLPK is distributed in the hope that it will be useful, but WITHOUT
    17 *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    18 *  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
    19 *  License for more details.
    20 *
    21 *  You should have received a copy of the GNU General Public License
    22 *  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
    23 ***********************************************************************/
    24 
    25 #include "glpios.h"
    26 
    27 /***********************************************************************
    28 *  prepare_row_info - prepare row info to determine implied bounds
    29 *
    30 *  Given a row (linear form)
    31 *
    32 *      n
    33 *     sum a[j] * x[j]                                                (1)
    34 *     j=1
    35 *
    36 *  and bounds of columns (variables)
    37 *
    38 *     l[j] <= x[j] <= u[j]                                           (2)
    39 *
    40 *  this routine computes f_min, j_min, f_max, j_max needed to determine
    41 *  implied bounds.
    42 *
    43 *  ALGORITHM
    44 *
    45 *  Let J+ = {j : a[j] > 0} and J- = {j : a[j] < 0}.
    46 *
    47 *  Parameters f_min and j_min are computed as follows:
    48 *
    49 *  1) if there is no x[k] such that k in J+ and l[k] = -inf or k in J-
    50 *     and u[k] = +inf, then
    51 *
    52 *     f_min :=   sum   a[j] * l[j] +   sum   a[j] * u[j]
    53 *              j in J+               j in J-
    54 *                                                                    (3)
    55 *     j_min := 0
    56 *
    57 *  2) if there is exactly one x[k] such that k in J+ and l[k] = -inf
    58 *     or k in J- and u[k] = +inf, then
    59 *
    60 *     f_min :=   sum       a[j] * l[j] +   sum       a[j] * u[j]
    61 *              j in J+\{k}               j in J-\{k}
    62 *                                                                    (4)
    63 *     j_min := k
    64 *
    65 *  3) if there are two or more x[k] such that k in J+ and l[k] = -inf
    66 *     or k in J- and u[k] = +inf, then
    67 *
    68 *     f_min := -inf
    69 *                                                                    (5)
    70 *     j_min := 0
    71 *
    72 *  Parameters f_max and j_max are computed in a similar way as follows:
    73 *
    74 *  1) if there is no x[k] such that k in J+ and u[k] = +inf or k in J-
    75 *     and l[k] = -inf, then
    76 *
    77 *     f_max :=   sum   a[j] * u[j] +   sum   a[j] * l[j]
    78 *              j in J+               j in J-
    79 *                                                                    (6)
    80 *     j_max := 0
    81 *
    82 *  2) if there is exactly one x[k] such that k in J+ and u[k] = +inf
    83 *     or k in J- and l[k] = -inf, then
    84 *
    85 *     f_max :=   sum       a[j] * u[j] +   sum       a[j] * l[j]
    86 *              j in J+\{k}               j in J-\{k}
    87 *                                                                    (7)
    88 *     j_max := k
    89 *
    90 *  3) if there are two or more x[k] such that k in J+ and u[k] = +inf
    91 *     or k in J- and l[k] = -inf, then
    92 *
    93 *     f_max := +inf
    94 *                                                                    (8)
    95 *     j_max := 0                                                      */
    96 
    97 struct f_info
    98 {     int j_min, j_max;
    99       double f_min, f_max;
   100 };
   101 
   102 static void prepare_row_info(int n, const double a[], const double l[],
   103       const double u[], struct f_info *f)
   104 {     int j, j_min, j_max;
   105       double f_min, f_max;
   106       xassert(n >= 0);
   107       /* determine f_min and j_min */
   108       f_min = 0.0, j_min = 0;
   109       for (j = 1; j <= n; j++)
   110       {  if (a[j] > 0.0)
   111          {  if (l[j] == -DBL_MAX)
   112             {  if (j_min == 0)
   113                   j_min = j;
   114                else
   115                {  f_min = -DBL_MAX, j_min = 0;
   116                   break;
   117                }
   118             }
   119             else
   120                f_min += a[j] * l[j];
   121          }
   122          else if (a[j] < 0.0)
   123          {  if (u[j] == +DBL_MAX)
   124             {  if (j_min == 0)
   125                   j_min = j;
   126                else
   127                {  f_min = -DBL_MAX, j_min = 0;
   128                   break;
   129                }
   130             }
   131             else
   132                f_min += a[j] * u[j];
   133          }
   134          else
   135             xassert(a != a);
   136       }
   137       f->f_min = f_min, f->j_min = j_min;
   138       /* determine f_max and j_max */
   139       f_max = 0.0, j_max = 0;
   140       for (j = 1; j <= n; j++)
   141       {  if (a[j] > 0.0)
   142          {  if (u[j] == +DBL_MAX)
   143             {  if (j_max == 0)
   144                   j_max = j;
   145                else
   146                {  f_max = +DBL_MAX, j_max = 0;
   147                   break;
   148                }
   149             }
   150             else
   151                f_max += a[j] * u[j];
   152          }
   153          else if (a[j] < 0.0)
   154          {  if (l[j] == -DBL_MAX)
   155             {  if (j_max == 0)
   156                   j_max = j;
   157                else
   158                {  f_max = +DBL_MAX, j_max = 0;
   159                   break;
   160                }
   161             }
   162             else
   163                f_max += a[j] * l[j];
   164          }
   165          else
   166             xassert(a != a);
   167       }
   168       f->f_max = f_max, f->j_max = j_max;
   169       return;
   170 }
   171 
   172 /***********************************************************************
   173 *  row_implied_bounds - determine row implied bounds
   174 *
   175 *  Given a row (linear form)
   176 *
   177 *      n
   178 *     sum a[j] * x[j]
   179 *     j=1
   180 *
   181 *  and bounds of columns (variables)
   182 *
   183 *     l[j] <= x[j] <= u[j]
   184 *
   185 *  this routine determines implied bounds of the row.
   186 *
   187 *  ALGORITHM
   188 *
   189 *  Let J+ = {j : a[j] > 0} and J- = {j : a[j] < 0}.
   190 *
   191 *  The implied lower bound of the row is computed as follows:
   192 *
   193 *     L' :=   sum   a[j] * l[j] +   sum   a[j] * u[j]                (9)
   194 *           j in J+               j in J-
   195 *
   196 *  and as it follows from (3), (4), and (5):
   197 *
   198 *     L' := if j_min = 0 then f_min else -inf                       (10)
   199 *
   200 *  The implied upper bound of the row is computed as follows:
   201 *
   202 *     U' :=   sum   a[j] * u[j] +   sum   a[j] * l[j]               (11)
   203 *           j in J+               j in J-
   204 *
   205 *  and as it follows from (6), (7), and (8):
   206 *
   207 *     U' := if j_max = 0 then f_max else +inf                       (12)
   208 *
   209 *  The implied bounds are stored in locations LL and UU. */
   210 
   211 static void row_implied_bounds(const struct f_info *f, double *LL,
   212       double *UU)
   213 {     *LL = (f->j_min == 0 ? f->f_min : -DBL_MAX);
   214       *UU = (f->j_max == 0 ? f->f_max : +DBL_MAX);
   215       return;
   216 }
   217 
   218 /***********************************************************************
   219 *  col_implied_bounds - determine column implied bounds
   220 *
   221 *  Given a row (constraint)
   222 *
   223 *           n
   224 *     L <= sum a[j] * x[j] <= U                                     (13)
   225 *          j=1
   226 *
   227 *  and bounds of columns (variables)
   228 *
   229 *     l[j] <= x[j] <= u[j]
   230 *
   231 *  this routine determines implied bounds of variable x[k].
   232 *
   233 *  It is assumed that if L != -inf, the lower bound of the row can be
   234 *  active, and if U != +inf, the upper bound of the row can be active.
   235 *
   236 *  ALGORITHM
   237 *
   238 *  From (13) it follows that
   239 *
   240 *     L <= sum a[j] * x[j] + a[k] * x[k] <= U
   241 *          j!=k
   242 *  or
   243 *
   244 *     L - sum a[j] * x[j] <= a[k] * x[k] <= U - sum a[j] * x[j]
   245 *         j!=k                                  j!=k
   246 *
   247 *  Thus, if the row lower bound L can be active, implied lower bound of
   248 *  term a[k] * x[k] can be determined as follows:
   249 *
   250 *     ilb(a[k] * x[k]) = min(L - sum a[j] * x[j]) =
   251 *                                j!=k
   252 *                                                                   (14)
   253 *                      = L - max sum a[j] * x[j]
   254 *                            j!=k
   255 *
   256 *  where, as it follows from (6), (7), and (8)
   257 *
   258 *                           / f_max - a[k] * u[k], j_max = 0, a[k] > 0
   259 *                           |
   260 *                           | f_max - a[k] * l[k], j_max = 0, a[k] < 0
   261 *     max sum a[j] * x[j] = {
   262 *         j!=k              | f_max,               j_max = k
   263 *                           |
   264 *                           \ +inf,                j_max != 0
   265 *
   266 *  and if the upper bound U can be active, implied upper bound of term
   267 *  a[k] * x[k] can be determined as follows:
   268 *
   269 *     iub(a[k] * x[k]) = max(U - sum a[j] * x[j]) =
   270 *                                j!=k
   271 *                                                                   (15)
   272 *                      = U - min sum a[j] * x[j]
   273 *                            j!=k
   274 *
   275 *  where, as it follows from (3), (4), and (5)
   276 *
   277 *                           / f_min - a[k] * l[k], j_min = 0, a[k] > 0
   278 *                           |
   279 *                           | f_min - a[k] * u[k], j_min = 0, a[k] < 0
   280 *     min sum a[j] * x[j] = {
   281 *         j!=k              | f_min,               j_min = k
   282 *                           |
   283 *                           \ -inf,                j_min != 0
   284 *
   285 *  Since
   286 *
   287 *     ilb(a[k] * x[k]) <= a[k] * x[k] <= iub(a[k] * x[k])
   288 *
   289 *  implied lower and upper bounds of x[k] are determined as follows:
   290 *
   291 *     l'[k] := if a[k] > 0 then ilb / a[k] else ulb / a[k]          (16)
   292 *
   293 *     u'[k] := if a[k] > 0 then ulb / a[k] else ilb / a[k]          (17)
   294 *
   295 *  The implied bounds are stored in locations ll and uu. */
   296 
   297 static void col_implied_bounds(const struct f_info *f, int n,
   298       const double a[], double L, double U, const double l[],
   299       const double u[], int k, double *ll, double *uu)
   300 {     double ilb, iub;
   301       xassert(n >= 0);
   302       xassert(1 <= k && k <= n);
   303       /* determine implied lower bound of term a[k] * x[k] (14) */
   304       if (L == -DBL_MAX || f->f_max == +DBL_MAX)
   305          ilb = -DBL_MAX;
   306       else if (f->j_max == 0)
   307       {  if (a[k] > 0.0)
   308          {  xassert(u[k] != +DBL_MAX);
   309             ilb = L - (f->f_max - a[k] * u[k]);
   310          }
   311          else if (a[k] < 0.0)
   312          {  xassert(l[k] != -DBL_MAX);
   313             ilb = L - (f->f_max - a[k] * l[k]);
   314          }
   315          else
   316             xassert(a != a);
   317       }
   318       else if (f->j_max == k)
   319          ilb = L - f->f_max;
   320       else
   321          ilb = -DBL_MAX;
   322       /* determine implied upper bound of term a[k] * x[k] (15) */
   323       if (U == +DBL_MAX || f->f_min == -DBL_MAX)
   324          iub = +DBL_MAX;
   325       else if (f->j_min == 0)
   326       {  if (a[k] > 0.0)
   327          {  xassert(l[k] != -DBL_MAX);
   328             iub = U - (f->f_min - a[k] * l[k]);
   329          }
   330          else if (a[k] < 0.0)
   331          {  xassert(u[k] != +DBL_MAX);
   332             iub = U - (f->f_min - a[k] * u[k]);
   333          }
   334          else
   335             xassert(a != a);
   336       }
   337       else if (f->j_min == k)
   338          iub = U - f->f_min;
   339       else
   340          iub = +DBL_MAX;
   341       /* determine implied bounds of x[k] (16) and (17) */
   342 #if 1
   343       /* do not use a[k] if it has small magnitude to prevent wrong
   344          implied bounds; for example, 1e-15 * x1 >= x2 + x3, where
   345          x1 >= -10, x2, x3 >= 0, would lead to wrong conclusion that
   346          x1 >= 0 */
   347       if (fabs(a[k]) < 1e-6)
   348          *ll = -DBL_MAX, *uu = +DBL_MAX; else
   349 #endif
   350       if (a[k] > 0.0)
   351       {  *ll = (ilb == -DBL_MAX ? -DBL_MAX : ilb / a[k]);
   352          *uu = (iub == +DBL_MAX ? +DBL_MAX : iub / a[k]);
   353       }
   354       else if (a[k] < 0.0)
   355       {  *ll = (iub == +DBL_MAX ? -DBL_MAX : iub / a[k]);
   356          *uu = (ilb == -DBL_MAX ? +DBL_MAX : ilb / a[k]);
   357       }
   358       else
   359          xassert(a != a);
   360       return;
   361 }
   362 
   363 /***********************************************************************
   364 *  check_row_bounds - check and relax original row bounds
   365 *
   366 *  Given a row (constraint)
   367 *
   368 *           n
   369 *     L <= sum a[j] * x[j] <= U
   370 *          j=1
   371 *
   372 *  and bounds of columns (variables)
   373 *
   374 *     l[j] <= x[j] <= u[j]
   375 *
   376 *  this routine checks the original row bounds L and U for feasibility
   377 *  and redundancy. If the original lower bound L or/and upper bound U
   378 *  cannot be active due to bounds of variables, the routine remove them
   379 *  replacing by -inf or/and +inf, respectively.
   380 *
   381 *  If no primal infeasibility is detected, the routine returns zero,
   382 *  otherwise non-zero. */
   383 
   384 static int check_row_bounds(const struct f_info *f, double *L_,
   385       double *U_)
   386 {     int ret = 0;
   387       double L = *L_, U = *U_, LL, UU;
   388       /* determine implied bounds of the row */
   389       row_implied_bounds(f, &LL, &UU);
   390       /* check if the original lower bound is infeasible */
   391       if (L != -DBL_MAX)
   392       {  double eps = 1e-3 * (1.0 + fabs(L));
   393          if (UU < L - eps)
   394          {  ret = 1;
   395             goto done;
   396          }
   397       }
   398       /* check if the original upper bound is infeasible */
   399       if (U != +DBL_MAX)
   400       {  double eps = 1e-3 * (1.0 + fabs(U));
   401          if (LL > U + eps)
   402          {  ret = 1;
   403             goto done;
   404          }
   405       }
   406       /* check if the original lower bound is redundant */
   407       if (L != -DBL_MAX)
   408       {  double eps = 1e-12 * (1.0 + fabs(L));
   409          if (LL > L - eps)
   410          {  /* it cannot be active, so remove it */
   411             *L_ = -DBL_MAX;
   412          }
   413       }
   414       /* check if the original upper bound is redundant */
   415       if (U != +DBL_MAX)
   416       {  double eps = 1e-12 * (1.0 + fabs(U));
   417          if (UU < U + eps)
   418          {  /* it cannot be active, so remove it */
   419             *U_ = +DBL_MAX;
   420          }
   421       }
   422 done: return ret;
   423 }
   424 
   425 /***********************************************************************
   426 *  check_col_bounds - check and tighten original column bounds
   427 *
   428 *  Given a row (constraint)
   429 *
   430 *           n
   431 *     L <= sum a[j] * x[j] <= U
   432 *          j=1
   433 *
   434 *  and bounds of columns (variables)
   435 *
   436 *     l[j] <= x[j] <= u[j]
   437 *
   438 *  for column (variable) x[j] this routine checks the original column
   439 *  bounds l[j] and u[j] for feasibility and redundancy. If the original
   440 *  lower bound l[j] or/and upper bound u[j] cannot be active due to
   441 *  bounds of the constraint and other variables, the routine tighten
   442 *  them replacing by corresponding implied bounds, if possible.
   443 *
   444 *  NOTE: It is assumed that if L != -inf, the row lower bound can be
   445 *        active, and if U != +inf, the row upper bound can be active.
   446 *
   447 *  The flag means that variable x[j] is required to be integer.
   448 *
   449 *  New actual bounds for x[j] are stored in locations lj and uj.
   450 *
   451 *  If no primal infeasibility is detected, the routine returns zero,
   452 *  otherwise non-zero. */
   453 
   454 static int check_col_bounds(const struct f_info *f, int n,
   455       const double a[], double L, double U, const double l[],
   456       const double u[], int flag, int j, double *_lj, double *_uj)
   457 {     int ret = 0;
   458       double lj, uj, ll, uu;
   459       xassert(n >= 0);
   460       xassert(1 <= j && j <= n);
   461       lj = l[j], uj = u[j];
   462       /* determine implied bounds of the column */
   463       col_implied_bounds(f, n, a, L, U, l, u, j, &ll, &uu);
   464       /* if x[j] is integral, round its implied bounds */
   465       if (flag)
   466       {  if (ll != -DBL_MAX)
   467             ll = (ll - floor(ll) < 1e-3 ? floor(ll) : ceil(ll));
   468          if (uu != +DBL_MAX)
   469             uu = (ceil(uu) - uu < 1e-3 ? ceil(uu) : floor(uu));
   470       }
   471       /* check if the original lower bound is infeasible */
   472       if (lj != -DBL_MAX)
   473       {  double eps = 1e-3 * (1.0 + fabs(lj));
   474          if (uu < lj - eps)
   475          {  ret = 1;
   476             goto done;
   477          }
   478       }
   479       /* check if the original upper bound is infeasible */
   480       if (uj != +DBL_MAX)
   481       {  double eps = 1e-3 * (1.0 + fabs(uj));
   482          if (ll > uj + eps)
   483          {  ret = 1;
   484             goto done;
   485          }
   486       }
   487       /* check if the original lower bound is redundant */
   488       if (ll != -DBL_MAX)
   489       {  double eps = 1e-3 * (1.0 + fabs(ll));
   490          if (lj < ll - eps)
   491          {  /* it cannot be active, so tighten it */
   492             lj = ll;
   493          }
   494       }
   495       /* check if the original upper bound is redundant */
   496       if (uu != +DBL_MAX)
   497       {  double eps = 1e-3 * (1.0 + fabs(uu));
   498          if (uj > uu + eps)
   499          {  /* it cannot be active, so tighten it */
   500             uj = uu;
   501          }
   502       }
   503       /* due to round-off errors it may happen that lj > uj (although
   504          lj < uj + eps, since no primal infeasibility is detected), so
   505          adjuct the new actual bounds to provide lj <= uj */
   506       if (!(lj == -DBL_MAX || uj == +DBL_MAX))
   507       {  double t1 = fabs(lj), t2 = fabs(uj);
   508          double eps = 1e-10 * (1.0 + (t1 <= t2 ? t1 : t2));
   509          if (lj > uj - eps)
   510          {  if (lj == l[j])
   511                uj = lj;
   512             else if (uj == u[j])
   513                lj = uj;
   514             else if (t1 <= t2)
   515                uj = lj;
   516             else
   517                lj = uj;
   518          }
   519       }
   520       *_lj = lj, *_uj = uj;
   521 done: return ret;
   522 }
   523 
   524 /***********************************************************************
   525 *  check_efficiency - check if change in column bounds is efficient
   526 *
   527 *  Given the original bounds of a column l and u and its new actual
   528 *  bounds l' and u' (possibly tighten by the routine check_col_bounds)
   529 *  this routine checks if the change in the column bounds is efficient
   530 *  enough. If so, the routine returns non-zero, otherwise zero.
   531 *
   532 *  The flag means that the variable is required to be integer. */
   533 
   534 static int check_efficiency(int flag, double l, double u, double ll,
   535       double uu)
   536 {     int eff = 0;
   537       /* check efficiency for lower bound */
   538       if (l < ll)
   539       {  if (flag || l == -DBL_MAX)
   540             eff++;
   541          else
   542          {  double r;
   543             if (u == +DBL_MAX)
   544                r = 1.0 + fabs(l);
   545             else
   546                r = 1.0 + (u - l);
   547             if (ll - l >= 0.25 * r)
   548                eff++;
   549          }
   550       }
   551       /* check efficiency for upper bound */
   552       if (u > uu)
   553       {  if (flag || u == +DBL_MAX)
   554             eff++;
   555          else
   556          {  double r;
   557             if (l == -DBL_MAX)
   558                r = 1.0 + fabs(u);
   559             else
   560                r = 1.0 + (u - l);
   561             if (u - uu >= 0.25 * r)
   562                eff++;
   563          }
   564       }
   565       return eff;
   566 }
   567 
   568 /***********************************************************************
   569 *  basic_preprocessing - perform basic preprocessing
   570 *
   571 *  This routine performs basic preprocessing of the specified MIP that
   572 *  includes relaxing some row bounds and tightening some column bounds.
   573 *
   574 *  On entry the arrays L and U contains original row bounds, and the
   575 *  arrays l and u contains original column bounds:
   576 *
   577 *  L[0] is the lower bound of the objective row;
   578 *  L[i], i = 1,...,m, is the lower bound of i-th row;
   579 *  U[0] is the upper bound of the objective row;
   580 *  U[i], i = 1,...,m, is the upper bound of i-th row;
   581 *  l[0] is not used;
   582 *  l[j], j = 1,...,n, is the lower bound of j-th column;
   583 *  u[0] is not used;
   584 *  u[j], j = 1,...,n, is the upper bound of j-th column.
   585 *
   586 *  On exit the arrays L, U, l, and u contain new actual bounds of rows
   587 *  and column in the same locations.
   588 *
   589 *  The parameters nrs and num specify an initial list of rows to be
   590 *  processed:
   591 *
   592 *  nrs is the number of rows in the initial list, 0 <= nrs <= m+1;
   593 *  num[0] is not used;
   594 *  num[1,...,nrs] are row numbers (0 means the objective row).
   595 *
   596 *  The parameter max_pass specifies the maximal number of times that
   597 *  each row can be processed, max_pass > 0.
   598 *
   599 *  If no primal infeasibility is detected, the routine returns zero,
   600 *  otherwise non-zero. */
   601 
   602 static int basic_preprocessing(glp_prob *mip, double L[], double U[],
   603       double l[], double u[], int nrs, const int num[], int max_pass)
   604 {     int m = mip->m;
   605       int n = mip->n;
   606       struct f_info f;
   607       int i, j, k, len, size, ret = 0;
   608       int *ind, *list, *mark, *pass;
   609       double *val, *lb, *ub;
   610       xassert(0 <= nrs && nrs <= m+1);
   611       xassert(max_pass > 0);
   612       /* allocate working arrays */
   613       ind = xcalloc(1+n, sizeof(int));
   614       list = xcalloc(1+m+1, sizeof(int));
   615       mark = xcalloc(1+m+1, sizeof(int));
   616       memset(&mark[0], 0, (m+1) * sizeof(int));
   617       pass = xcalloc(1+m+1, sizeof(int));
   618       memset(&pass[0], 0, (m+1) * sizeof(int));
   619       val = xcalloc(1+n, sizeof(double));
   620       lb = xcalloc(1+n, sizeof(double));
   621       ub = xcalloc(1+n, sizeof(double));
   622       /* initialize the list of rows to be processed */
   623       size = 0;
   624       for (k = 1; k <= nrs; k++)
   625       {  i = num[k];
   626          xassert(0 <= i && i <= m);
   627          /* duplicate row numbers are not allowed */
   628          xassert(!mark[i]);
   629          list[++size] = i, mark[i] = 1;
   630       }
   631       xassert(size == nrs);
   632       /* process rows in the list until it becomes empty */
   633       while (size > 0)
   634       {  /* get a next row from the list */
   635          i = list[size--], mark[i] = 0;
   636          /* increase the row processing count */
   637          pass[i]++;
   638          /* if the row is free, skip it */
   639          if (L[i] == -DBL_MAX && U[i] == +DBL_MAX) continue;
   640          /* obtain coefficients of the row */
   641          len = 0;
   642          if (i == 0)
   643          {  for (j = 1; j <= n; j++)
   644             {  GLPCOL *col = mip->col[j];
   645                if (col->coef != 0.0)
   646                   len++, ind[len] = j, val[len] = col->coef;
   647             }
   648          }
   649          else
   650          {  GLPROW *row = mip->row[i];
   651             GLPAIJ *aij;
   652             for (aij = row->ptr; aij != NULL; aij = aij->r_next)
   653                len++, ind[len] = aij->col->j, val[len] = aij->val;
   654          }
   655          /* determine lower and upper bounds of columns corresponding
   656             to non-zero row coefficients */
   657          for (k = 1; k <= len; k++)
   658             j = ind[k], lb[k] = l[j], ub[k] = u[j];
   659          /* prepare the row info to determine implied bounds */
   660          prepare_row_info(len, val, lb, ub, &f);
   661          /* check and relax bounds of the row */
   662          if (check_row_bounds(&f, &L[i], &U[i]))
   663          {  /* the feasible region is empty */
   664             ret = 1;
   665             goto done;
   666          }
   667          /* if the row became free, drop it */
   668          if (L[i] == -DBL_MAX && U[i] == +DBL_MAX) continue;
   669          /* process columns having non-zero coefficients in the row */
   670          for (k = 1; k <= len; k++)
   671          {  GLPCOL *col;
   672             int flag, eff;
   673             double ll, uu;
   674             /* take a next column in the row */
   675             j = ind[k], col = mip->col[j];
   676             flag = col->kind != GLP_CV;
   677             /* check and tighten bounds of the column */
   678             if (check_col_bounds(&f, len, val, L[i], U[i], lb, ub,
   679                 flag, k, &ll, &uu))
   680             {  /* the feasible region is empty */
   681                ret = 1;
   682                goto done;
   683             }
   684             /* check if change in the column bounds is efficient */
   685             eff = check_efficiency(flag, l[j], u[j], ll, uu);
   686             /* set new actual bounds of the column */
   687             l[j] = ll, u[j] = uu;
   688             /* if the change is efficient, add all rows affected by the
   689                corresponding column, to the list */
   690             if (eff > 0)
   691             {  GLPAIJ *aij;
   692                for (aij = col->ptr; aij != NULL; aij = aij->c_next)
   693                {  int ii = aij->row->i;
   694                   /* if the row was processed maximal number of times,
   695                      skip it */
   696                   if (pass[ii] >= max_pass) continue;
   697                   /* if the row is free, skip it */
   698                   if (L[ii] == -DBL_MAX && U[ii] == +DBL_MAX) continue;
   699                   /* put the row into the list */
   700                   if (mark[ii] == 0)
   701                   {  xassert(size <= m);
   702                      list[++size] = ii, mark[ii] = 1;
   703                   }
   704                }
   705             }
   706          }
   707       }
   708 done: /* free working arrays */
   709       xfree(ind);
   710       xfree(list);
   711       xfree(mark);
   712       xfree(pass);
   713       xfree(val);
   714       xfree(lb);
   715       xfree(ub);
   716       return ret;
   717 }
   718 
   719 /***********************************************************************
   720 *  NAME
   721 *
   722 *  ios_preprocess_node - preprocess current subproblem
   723 *
   724 *  SYNOPSIS
   725 *
   726 *  #include "glpios.h"
   727 *  int ios_preprocess_node(glp_tree *tree, int max_pass);
   728 *
   729 *  DESCRIPTION
   730 *
   731 *  The routine ios_preprocess_node performs basic preprocessing of the
   732 *  current subproblem.
   733 *
   734 *  RETURNS
   735 *
   736 *  If no primal infeasibility is detected, the routine returns zero,
   737 *  otherwise non-zero. */
   738 
   739 int ios_preprocess_node(glp_tree *tree, int max_pass)
   740 {     glp_prob *mip = tree->mip;
   741       int m = mip->m;
   742       int n = mip->n;
   743       int i, j, nrs, *num, ret = 0;
   744       double *L, *U, *l, *u;
   745       /* the current subproblem must exist */
   746       xassert(tree->curr != NULL);
   747       /* determine original row bounds */
   748       L = xcalloc(1+m, sizeof(double));
   749       U = xcalloc(1+m, sizeof(double));
   750       switch (mip->mip_stat)
   751       {  case GLP_UNDEF:
   752             L[0] = -DBL_MAX, U[0] = +DBL_MAX;
   753             break;
   754          case GLP_FEAS:
   755             switch (mip->dir)
   756             {  case GLP_MIN:
   757                   L[0] = -DBL_MAX, U[0] = mip->mip_obj - mip->c0;
   758                   break;
   759                case GLP_MAX:
   760                   L[0] = mip->mip_obj - mip->c0, U[0] = +DBL_MAX;
   761                   break;
   762                default:
   763                   xassert(mip != mip);
   764             }
   765             break;
   766          default:
   767             xassert(mip != mip);
   768       }
   769       for (i = 1; i <= m; i++)
   770       {  L[i] = glp_get_row_lb(mip, i);
   771          U[i] = glp_get_row_ub(mip, i);
   772       }
   773       /* determine original column bounds */
   774       l = xcalloc(1+n, sizeof(double));
   775       u = xcalloc(1+n, sizeof(double));
   776       for (j = 1; j <= n; j++)
   777       {  l[j] = glp_get_col_lb(mip, j);
   778          u[j] = glp_get_col_ub(mip, j);
   779       }
   780       /* build the initial list of rows to be analyzed */
   781       nrs = m + 1;
   782       num = xcalloc(1+nrs, sizeof(int));
   783       for (i = 1; i <= nrs; i++) num[i] = i - 1;
   784       /* perform basic preprocessing */
   785       if (basic_preprocessing(mip , L, U, l, u, nrs, num, max_pass))
   786       {  ret = 1;
   787          goto done;
   788       }
   789       /* set new actual (relaxed) row bounds */
   790       for (i = 1; i <= m; i++)
   791       {  /* consider only non-active rows to keep dual feasibility */
   792          if (glp_get_row_stat(mip, i) == GLP_BS)
   793          {  if (L[i] == -DBL_MAX && U[i] == +DBL_MAX)
   794                glp_set_row_bnds(mip, i, GLP_FR, 0.0, 0.0);
   795             else if (U[i] == +DBL_MAX)
   796                glp_set_row_bnds(mip, i, GLP_LO, L[i], 0.0);
   797             else if (L[i] == -DBL_MAX)
   798                glp_set_row_bnds(mip, i, GLP_UP, 0.0, U[i]);
   799          }
   800       }
   801       /* set new actual (tightened) column bounds */
   802       for (j = 1; j <= n; j++)
   803       {  int type;
   804          if (l[j] == -DBL_MAX && u[j] == +DBL_MAX)
   805             type = GLP_FR;
   806          else if (u[j] == +DBL_MAX)
   807             type = GLP_LO;
   808          else if (l[j] == -DBL_MAX)
   809             type = GLP_UP;
   810          else if (l[j] != u[j])
   811             type = GLP_DB;
   812          else
   813             type = GLP_FX;
   814          glp_set_col_bnds(mip, j, type, l[j], u[j]);
   815       }
   816 done: /* free working arrays and return */
   817       xfree(L);
   818       xfree(U);
   819       xfree(l);
   820       xfree(u);
   821       xfree(num);
   822       return ret;
   823 }
   824 
   825 /* eof */