1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpios02.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,825 @@
1.4 +/* glpios02.c (preprocess current subproblem) */
1.5 +
1.6 +/***********************************************************************
1.7 +* This code is part of GLPK (GNU Linear Programming Kit).
1.8 +*
1.9 +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
1.10 +* 2009, 2010 Andrew Makhorin, Department for Applied Informatics,
1.11 +* Moscow Aviation Institute, Moscow, Russia. All rights reserved.
1.12 +* E-mail: <mao@gnu.org>.
1.13 +*
1.14 +* GLPK is free software: you can redistribute it and/or modify it
1.15 +* under the terms of the GNU General Public License as published by
1.16 +* the Free Software Foundation, either version 3 of the License, or
1.17 +* (at your option) any later version.
1.18 +*
1.19 +* GLPK is distributed in the hope that it will be useful, but WITHOUT
1.20 +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.21 +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
1.22 +* License for more details.
1.23 +*
1.24 +* You should have received a copy of the GNU General Public License
1.25 +* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
1.26 +***********************************************************************/
1.27 +
1.28 +#include "glpios.h"
1.29 +
1.30 +/***********************************************************************
1.31 +* prepare_row_info - prepare row info to determine implied bounds
1.32 +*
1.33 +* Given a row (linear form)
1.34 +*
1.35 +* n
1.36 +* sum a[j] * x[j] (1)
1.37 +* j=1
1.38 +*
1.39 +* and bounds of columns (variables)
1.40 +*
1.41 +* l[j] <= x[j] <= u[j] (2)
1.42 +*
1.43 +* this routine computes f_min, j_min, f_max, j_max needed to determine
1.44 +* implied bounds.
1.45 +*
1.46 +* ALGORITHM
1.47 +*
1.48 +* Let J+ = {j : a[j] > 0} and J- = {j : a[j] < 0}.
1.49 +*
1.50 +* Parameters f_min and j_min are computed as follows:
1.51 +*
1.52 +* 1) if there is no x[k] such that k in J+ and l[k] = -inf or k in J-
1.53 +* and u[k] = +inf, then
1.54 +*
1.55 +* f_min := sum a[j] * l[j] + sum a[j] * u[j]
1.56 +* j in J+ j in J-
1.57 +* (3)
1.58 +* j_min := 0
1.59 +*
1.60 +* 2) if there is exactly one x[k] such that k in J+ and l[k] = -inf
1.61 +* or k in J- and u[k] = +inf, then
1.62 +*
1.63 +* f_min := sum a[j] * l[j] + sum a[j] * u[j]
1.64 +* j in J+\{k} j in J-\{k}
1.65 +* (4)
1.66 +* j_min := k
1.67 +*
1.68 +* 3) if there are two or more x[k] such that k in J+ and l[k] = -inf
1.69 +* or k in J- and u[k] = +inf, then
1.70 +*
1.71 +* f_min := -inf
1.72 +* (5)
1.73 +* j_min := 0
1.74 +*
1.75 +* Parameters f_max and j_max are computed in a similar way as follows:
1.76 +*
1.77 +* 1) if there is no x[k] such that k in J+ and u[k] = +inf or k in J-
1.78 +* and l[k] = -inf, then
1.79 +*
1.80 +* f_max := sum a[j] * u[j] + sum a[j] * l[j]
1.81 +* j in J+ j in J-
1.82 +* (6)
1.83 +* j_max := 0
1.84 +*
1.85 +* 2) if there is exactly one x[k] such that k in J+ and u[k] = +inf
1.86 +* or k in J- and l[k] = -inf, then
1.87 +*
1.88 +* f_max := sum a[j] * u[j] + sum a[j] * l[j]
1.89 +* j in J+\{k} j in J-\{k}
1.90 +* (7)
1.91 +* j_max := k
1.92 +*
1.93 +* 3) if there are two or more x[k] such that k in J+ and u[k] = +inf
1.94 +* or k in J- and l[k] = -inf, then
1.95 +*
1.96 +* f_max := +inf
1.97 +* (8)
1.98 +* j_max := 0 */
1.99 +
1.100 +struct f_info
1.101 +{ int j_min, j_max;
1.102 + double f_min, f_max;
1.103 +};
1.104 +
1.105 +static void prepare_row_info(int n, const double a[], const double l[],
1.106 + const double u[], struct f_info *f)
1.107 +{ int j, j_min, j_max;
1.108 + double f_min, f_max;
1.109 + xassert(n >= 0);
1.110 + /* determine f_min and j_min */
1.111 + f_min = 0.0, j_min = 0;
1.112 + for (j = 1; j <= n; j++)
1.113 + { if (a[j] > 0.0)
1.114 + { if (l[j] == -DBL_MAX)
1.115 + { if (j_min == 0)
1.116 + j_min = j;
1.117 + else
1.118 + { f_min = -DBL_MAX, j_min = 0;
1.119 + break;
1.120 + }
1.121 + }
1.122 + else
1.123 + f_min += a[j] * l[j];
1.124 + }
1.125 + else if (a[j] < 0.0)
1.126 + { if (u[j] == +DBL_MAX)
1.127 + { if (j_min == 0)
1.128 + j_min = j;
1.129 + else
1.130 + { f_min = -DBL_MAX, j_min = 0;
1.131 + break;
1.132 + }
1.133 + }
1.134 + else
1.135 + f_min += a[j] * u[j];
1.136 + }
1.137 + else
1.138 + xassert(a != a);
1.139 + }
1.140 + f->f_min = f_min, f->j_min = j_min;
1.141 + /* determine f_max and j_max */
1.142 + f_max = 0.0, j_max = 0;
1.143 + for (j = 1; j <= n; j++)
1.144 + { if (a[j] > 0.0)
1.145 + { if (u[j] == +DBL_MAX)
1.146 + { if (j_max == 0)
1.147 + j_max = j;
1.148 + else
1.149 + { f_max = +DBL_MAX, j_max = 0;
1.150 + break;
1.151 + }
1.152 + }
1.153 + else
1.154 + f_max += a[j] * u[j];
1.155 + }
1.156 + else if (a[j] < 0.0)
1.157 + { if (l[j] == -DBL_MAX)
1.158 + { if (j_max == 0)
1.159 + j_max = j;
1.160 + else
1.161 + { f_max = +DBL_MAX, j_max = 0;
1.162 + break;
1.163 + }
1.164 + }
1.165 + else
1.166 + f_max += a[j] * l[j];
1.167 + }
1.168 + else
1.169 + xassert(a != a);
1.170 + }
1.171 + f->f_max = f_max, f->j_max = j_max;
1.172 + return;
1.173 +}
1.174 +
1.175 +/***********************************************************************
1.176 +* row_implied_bounds - determine row implied bounds
1.177 +*
1.178 +* Given a row (linear form)
1.179 +*
1.180 +* n
1.181 +* sum a[j] * x[j]
1.182 +* j=1
1.183 +*
1.184 +* and bounds of columns (variables)
1.185 +*
1.186 +* l[j] <= x[j] <= u[j]
1.187 +*
1.188 +* this routine determines implied bounds of the row.
1.189 +*
1.190 +* ALGORITHM
1.191 +*
1.192 +* Let J+ = {j : a[j] > 0} and J- = {j : a[j] < 0}.
1.193 +*
1.194 +* The implied lower bound of the row is computed as follows:
1.195 +*
1.196 +* L' := sum a[j] * l[j] + sum a[j] * u[j] (9)
1.197 +* j in J+ j in J-
1.198 +*
1.199 +* and as it follows from (3), (4), and (5):
1.200 +*
1.201 +* L' := if j_min = 0 then f_min else -inf (10)
1.202 +*
1.203 +* The implied upper bound of the row is computed as follows:
1.204 +*
1.205 +* U' := sum a[j] * u[j] + sum a[j] * l[j] (11)
1.206 +* j in J+ j in J-
1.207 +*
1.208 +* and as it follows from (6), (7), and (8):
1.209 +*
1.210 +* U' := if j_max = 0 then f_max else +inf (12)
1.211 +*
1.212 +* The implied bounds are stored in locations LL and UU. */
1.213 +
1.214 +static void row_implied_bounds(const struct f_info *f, double *LL,
1.215 + double *UU)
1.216 +{ *LL = (f->j_min == 0 ? f->f_min : -DBL_MAX);
1.217 + *UU = (f->j_max == 0 ? f->f_max : +DBL_MAX);
1.218 + return;
1.219 +}
1.220 +
1.221 +/***********************************************************************
1.222 +* col_implied_bounds - determine column implied bounds
1.223 +*
1.224 +* Given a row (constraint)
1.225 +*
1.226 +* n
1.227 +* L <= sum a[j] * x[j] <= U (13)
1.228 +* j=1
1.229 +*
1.230 +* and bounds of columns (variables)
1.231 +*
1.232 +* l[j] <= x[j] <= u[j]
1.233 +*
1.234 +* this routine determines implied bounds of variable x[k].
1.235 +*
1.236 +* It is assumed that if L != -inf, the lower bound of the row can be
1.237 +* active, and if U != +inf, the upper bound of the row can be active.
1.238 +*
1.239 +* ALGORITHM
1.240 +*
1.241 +* From (13) it follows that
1.242 +*
1.243 +* L <= sum a[j] * x[j] + a[k] * x[k] <= U
1.244 +* j!=k
1.245 +* or
1.246 +*
1.247 +* L - sum a[j] * x[j] <= a[k] * x[k] <= U - sum a[j] * x[j]
1.248 +* j!=k j!=k
1.249 +*
1.250 +* Thus, if the row lower bound L can be active, implied lower bound of
1.251 +* term a[k] * x[k] can be determined as follows:
1.252 +*
1.253 +* ilb(a[k] * x[k]) = min(L - sum a[j] * x[j]) =
1.254 +* j!=k
1.255 +* (14)
1.256 +* = L - max sum a[j] * x[j]
1.257 +* j!=k
1.258 +*
1.259 +* where, as it follows from (6), (7), and (8)
1.260 +*
1.261 +* / f_max - a[k] * u[k], j_max = 0, a[k] > 0
1.262 +* |
1.263 +* | f_max - a[k] * l[k], j_max = 0, a[k] < 0
1.264 +* max sum a[j] * x[j] = {
1.265 +* j!=k | f_max, j_max = k
1.266 +* |
1.267 +* \ +inf, j_max != 0
1.268 +*
1.269 +* and if the upper bound U can be active, implied upper bound of term
1.270 +* a[k] * x[k] can be determined as follows:
1.271 +*
1.272 +* iub(a[k] * x[k]) = max(U - sum a[j] * x[j]) =
1.273 +* j!=k
1.274 +* (15)
1.275 +* = U - min sum a[j] * x[j]
1.276 +* j!=k
1.277 +*
1.278 +* where, as it follows from (3), (4), and (5)
1.279 +*
1.280 +* / f_min - a[k] * l[k], j_min = 0, a[k] > 0
1.281 +* |
1.282 +* | f_min - a[k] * u[k], j_min = 0, a[k] < 0
1.283 +* min sum a[j] * x[j] = {
1.284 +* j!=k | f_min, j_min = k
1.285 +* |
1.286 +* \ -inf, j_min != 0
1.287 +*
1.288 +* Since
1.289 +*
1.290 +* ilb(a[k] * x[k]) <= a[k] * x[k] <= iub(a[k] * x[k])
1.291 +*
1.292 +* implied lower and upper bounds of x[k] are determined as follows:
1.293 +*
1.294 +* l'[k] := if a[k] > 0 then ilb / a[k] else ulb / a[k] (16)
1.295 +*
1.296 +* u'[k] := if a[k] > 0 then ulb / a[k] else ilb / a[k] (17)
1.297 +*
1.298 +* The implied bounds are stored in locations ll and uu. */
1.299 +
1.300 +static void col_implied_bounds(const struct f_info *f, int n,
1.301 + const double a[], double L, double U, const double l[],
1.302 + const double u[], int k, double *ll, double *uu)
1.303 +{ double ilb, iub;
1.304 + xassert(n >= 0);
1.305 + xassert(1 <= k && k <= n);
1.306 + /* determine implied lower bound of term a[k] * x[k] (14) */
1.307 + if (L == -DBL_MAX || f->f_max == +DBL_MAX)
1.308 + ilb = -DBL_MAX;
1.309 + else if (f->j_max == 0)
1.310 + { if (a[k] > 0.0)
1.311 + { xassert(u[k] != +DBL_MAX);
1.312 + ilb = L - (f->f_max - a[k] * u[k]);
1.313 + }
1.314 + else if (a[k] < 0.0)
1.315 + { xassert(l[k] != -DBL_MAX);
1.316 + ilb = L - (f->f_max - a[k] * l[k]);
1.317 + }
1.318 + else
1.319 + xassert(a != a);
1.320 + }
1.321 + else if (f->j_max == k)
1.322 + ilb = L - f->f_max;
1.323 + else
1.324 + ilb = -DBL_MAX;
1.325 + /* determine implied upper bound of term a[k] * x[k] (15) */
1.326 + if (U == +DBL_MAX || f->f_min == -DBL_MAX)
1.327 + iub = +DBL_MAX;
1.328 + else if (f->j_min == 0)
1.329 + { if (a[k] > 0.0)
1.330 + { xassert(l[k] != -DBL_MAX);
1.331 + iub = U - (f->f_min - a[k] * l[k]);
1.332 + }
1.333 + else if (a[k] < 0.0)
1.334 + { xassert(u[k] != +DBL_MAX);
1.335 + iub = U - (f->f_min - a[k] * u[k]);
1.336 + }
1.337 + else
1.338 + xassert(a != a);
1.339 + }
1.340 + else if (f->j_min == k)
1.341 + iub = U - f->f_min;
1.342 + else
1.343 + iub = +DBL_MAX;
1.344 + /* determine implied bounds of x[k] (16) and (17) */
1.345 +#if 1
1.346 + /* do not use a[k] if it has small magnitude to prevent wrong
1.347 + implied bounds; for example, 1e-15 * x1 >= x2 + x3, where
1.348 + x1 >= -10, x2, x3 >= 0, would lead to wrong conclusion that
1.349 + x1 >= 0 */
1.350 + if (fabs(a[k]) < 1e-6)
1.351 + *ll = -DBL_MAX, *uu = +DBL_MAX; else
1.352 +#endif
1.353 + if (a[k] > 0.0)
1.354 + { *ll = (ilb == -DBL_MAX ? -DBL_MAX : ilb / a[k]);
1.355 + *uu = (iub == +DBL_MAX ? +DBL_MAX : iub / a[k]);
1.356 + }
1.357 + else if (a[k] < 0.0)
1.358 + { *ll = (iub == +DBL_MAX ? -DBL_MAX : iub / a[k]);
1.359 + *uu = (ilb == -DBL_MAX ? +DBL_MAX : ilb / a[k]);
1.360 + }
1.361 + else
1.362 + xassert(a != a);
1.363 + return;
1.364 +}
1.365 +
1.366 +/***********************************************************************
1.367 +* check_row_bounds - check and relax original row bounds
1.368 +*
1.369 +* Given a row (constraint)
1.370 +*
1.371 +* n
1.372 +* L <= sum a[j] * x[j] <= U
1.373 +* j=1
1.374 +*
1.375 +* and bounds of columns (variables)
1.376 +*
1.377 +* l[j] <= x[j] <= u[j]
1.378 +*
1.379 +* this routine checks the original row bounds L and U for feasibility
1.380 +* and redundancy. If the original lower bound L or/and upper bound U
1.381 +* cannot be active due to bounds of variables, the routine remove them
1.382 +* replacing by -inf or/and +inf, respectively.
1.383 +*
1.384 +* If no primal infeasibility is detected, the routine returns zero,
1.385 +* otherwise non-zero. */
1.386 +
1.387 +static int check_row_bounds(const struct f_info *f, double *L_,
1.388 + double *U_)
1.389 +{ int ret = 0;
1.390 + double L = *L_, U = *U_, LL, UU;
1.391 + /* determine implied bounds of the row */
1.392 + row_implied_bounds(f, &LL, &UU);
1.393 + /* check if the original lower bound is infeasible */
1.394 + if (L != -DBL_MAX)
1.395 + { double eps = 1e-3 * (1.0 + fabs(L));
1.396 + if (UU < L - eps)
1.397 + { ret = 1;
1.398 + goto done;
1.399 + }
1.400 + }
1.401 + /* check if the original upper bound is infeasible */
1.402 + if (U != +DBL_MAX)
1.403 + { double eps = 1e-3 * (1.0 + fabs(U));
1.404 + if (LL > U + eps)
1.405 + { ret = 1;
1.406 + goto done;
1.407 + }
1.408 + }
1.409 + /* check if the original lower bound is redundant */
1.410 + if (L != -DBL_MAX)
1.411 + { double eps = 1e-12 * (1.0 + fabs(L));
1.412 + if (LL > L - eps)
1.413 + { /* it cannot be active, so remove it */
1.414 + *L_ = -DBL_MAX;
1.415 + }
1.416 + }
1.417 + /* check if the original upper bound is redundant */
1.418 + if (U != +DBL_MAX)
1.419 + { double eps = 1e-12 * (1.0 + fabs(U));
1.420 + if (UU < U + eps)
1.421 + { /* it cannot be active, so remove it */
1.422 + *U_ = +DBL_MAX;
1.423 + }
1.424 + }
1.425 +done: return ret;
1.426 +}
1.427 +
1.428 +/***********************************************************************
1.429 +* check_col_bounds - check and tighten original column bounds
1.430 +*
1.431 +* Given a row (constraint)
1.432 +*
1.433 +* n
1.434 +* L <= sum a[j] * x[j] <= U
1.435 +* j=1
1.436 +*
1.437 +* and bounds of columns (variables)
1.438 +*
1.439 +* l[j] <= x[j] <= u[j]
1.440 +*
1.441 +* for column (variable) x[j] this routine checks the original column
1.442 +* bounds l[j] and u[j] for feasibility and redundancy. If the original
1.443 +* lower bound l[j] or/and upper bound u[j] cannot be active due to
1.444 +* bounds of the constraint and other variables, the routine tighten
1.445 +* them replacing by corresponding implied bounds, if possible.
1.446 +*
1.447 +* NOTE: It is assumed that if L != -inf, the row lower bound can be
1.448 +* active, and if U != +inf, the row upper bound can be active.
1.449 +*
1.450 +* The flag means that variable x[j] is required to be integer.
1.451 +*
1.452 +* New actual bounds for x[j] are stored in locations lj and uj.
1.453 +*
1.454 +* If no primal infeasibility is detected, the routine returns zero,
1.455 +* otherwise non-zero. */
1.456 +
1.457 +static int check_col_bounds(const struct f_info *f, int n,
1.458 + const double a[], double L, double U, const double l[],
1.459 + const double u[], int flag, int j, double *_lj, double *_uj)
1.460 +{ int ret = 0;
1.461 + double lj, uj, ll, uu;
1.462 + xassert(n >= 0);
1.463 + xassert(1 <= j && j <= n);
1.464 + lj = l[j], uj = u[j];
1.465 + /* determine implied bounds of the column */
1.466 + col_implied_bounds(f, n, a, L, U, l, u, j, &ll, &uu);
1.467 + /* if x[j] is integral, round its implied bounds */
1.468 + if (flag)
1.469 + { if (ll != -DBL_MAX)
1.470 + ll = (ll - floor(ll) < 1e-3 ? floor(ll) : ceil(ll));
1.471 + if (uu != +DBL_MAX)
1.472 + uu = (ceil(uu) - uu < 1e-3 ? ceil(uu) : floor(uu));
1.473 + }
1.474 + /* check if the original lower bound is infeasible */
1.475 + if (lj != -DBL_MAX)
1.476 + { double eps = 1e-3 * (1.0 + fabs(lj));
1.477 + if (uu < lj - eps)
1.478 + { ret = 1;
1.479 + goto done;
1.480 + }
1.481 + }
1.482 + /* check if the original upper bound is infeasible */
1.483 + if (uj != +DBL_MAX)
1.484 + { double eps = 1e-3 * (1.0 + fabs(uj));
1.485 + if (ll > uj + eps)
1.486 + { ret = 1;
1.487 + goto done;
1.488 + }
1.489 + }
1.490 + /* check if the original lower bound is redundant */
1.491 + if (ll != -DBL_MAX)
1.492 + { double eps = 1e-3 * (1.0 + fabs(ll));
1.493 + if (lj < ll - eps)
1.494 + { /* it cannot be active, so tighten it */
1.495 + lj = ll;
1.496 + }
1.497 + }
1.498 + /* check if the original upper bound is redundant */
1.499 + if (uu != +DBL_MAX)
1.500 + { double eps = 1e-3 * (1.0 + fabs(uu));
1.501 + if (uj > uu + eps)
1.502 + { /* it cannot be active, so tighten it */
1.503 + uj = uu;
1.504 + }
1.505 + }
1.506 + /* due to round-off errors it may happen that lj > uj (although
1.507 + lj < uj + eps, since no primal infeasibility is detected), so
1.508 + adjuct the new actual bounds to provide lj <= uj */
1.509 + if (!(lj == -DBL_MAX || uj == +DBL_MAX))
1.510 + { double t1 = fabs(lj), t2 = fabs(uj);
1.511 + double eps = 1e-10 * (1.0 + (t1 <= t2 ? t1 : t2));
1.512 + if (lj > uj - eps)
1.513 + { if (lj == l[j])
1.514 + uj = lj;
1.515 + else if (uj == u[j])
1.516 + lj = uj;
1.517 + else if (t1 <= t2)
1.518 + uj = lj;
1.519 + else
1.520 + lj = uj;
1.521 + }
1.522 + }
1.523 + *_lj = lj, *_uj = uj;
1.524 +done: return ret;
1.525 +}
1.526 +
1.527 +/***********************************************************************
1.528 +* check_efficiency - check if change in column bounds is efficient
1.529 +*
1.530 +* Given the original bounds of a column l and u and its new actual
1.531 +* bounds l' and u' (possibly tighten by the routine check_col_bounds)
1.532 +* this routine checks if the change in the column bounds is efficient
1.533 +* enough. If so, the routine returns non-zero, otherwise zero.
1.534 +*
1.535 +* The flag means that the variable is required to be integer. */
1.536 +
1.537 +static int check_efficiency(int flag, double l, double u, double ll,
1.538 + double uu)
1.539 +{ int eff = 0;
1.540 + /* check efficiency for lower bound */
1.541 + if (l < ll)
1.542 + { if (flag || l == -DBL_MAX)
1.543 + eff++;
1.544 + else
1.545 + { double r;
1.546 + if (u == +DBL_MAX)
1.547 + r = 1.0 + fabs(l);
1.548 + else
1.549 + r = 1.0 + (u - l);
1.550 + if (ll - l >= 0.25 * r)
1.551 + eff++;
1.552 + }
1.553 + }
1.554 + /* check efficiency for upper bound */
1.555 + if (u > uu)
1.556 + { if (flag || u == +DBL_MAX)
1.557 + eff++;
1.558 + else
1.559 + { double r;
1.560 + if (l == -DBL_MAX)
1.561 + r = 1.0 + fabs(u);
1.562 + else
1.563 + r = 1.0 + (u - l);
1.564 + if (u - uu >= 0.25 * r)
1.565 + eff++;
1.566 + }
1.567 + }
1.568 + return eff;
1.569 +}
1.570 +
1.571 +/***********************************************************************
1.572 +* basic_preprocessing - perform basic preprocessing
1.573 +*
1.574 +* This routine performs basic preprocessing of the specified MIP that
1.575 +* includes relaxing some row bounds and tightening some column bounds.
1.576 +*
1.577 +* On entry the arrays L and U contains original row bounds, and the
1.578 +* arrays l and u contains original column bounds:
1.579 +*
1.580 +* L[0] is the lower bound of the objective row;
1.581 +* L[i], i = 1,...,m, is the lower bound of i-th row;
1.582 +* U[0] is the upper bound of the objective row;
1.583 +* U[i], i = 1,...,m, is the upper bound of i-th row;
1.584 +* l[0] is not used;
1.585 +* l[j], j = 1,...,n, is the lower bound of j-th column;
1.586 +* u[0] is not used;
1.587 +* u[j], j = 1,...,n, is the upper bound of j-th column.
1.588 +*
1.589 +* On exit the arrays L, U, l, and u contain new actual bounds of rows
1.590 +* and column in the same locations.
1.591 +*
1.592 +* The parameters nrs and num specify an initial list of rows to be
1.593 +* processed:
1.594 +*
1.595 +* nrs is the number of rows in the initial list, 0 <= nrs <= m+1;
1.596 +* num[0] is not used;
1.597 +* num[1,...,nrs] are row numbers (0 means the objective row).
1.598 +*
1.599 +* The parameter max_pass specifies the maximal number of times that
1.600 +* each row can be processed, max_pass > 0.
1.601 +*
1.602 +* If no primal infeasibility is detected, the routine returns zero,
1.603 +* otherwise non-zero. */
1.604 +
1.605 +static int basic_preprocessing(glp_prob *mip, double L[], double U[],
1.606 + double l[], double u[], int nrs, const int num[], int max_pass)
1.607 +{ int m = mip->m;
1.608 + int n = mip->n;
1.609 + struct f_info f;
1.610 + int i, j, k, len, size, ret = 0;
1.611 + int *ind, *list, *mark, *pass;
1.612 + double *val, *lb, *ub;
1.613 + xassert(0 <= nrs && nrs <= m+1);
1.614 + xassert(max_pass > 0);
1.615 + /* allocate working arrays */
1.616 + ind = xcalloc(1+n, sizeof(int));
1.617 + list = xcalloc(1+m+1, sizeof(int));
1.618 + mark = xcalloc(1+m+1, sizeof(int));
1.619 + memset(&mark[0], 0, (m+1) * sizeof(int));
1.620 + pass = xcalloc(1+m+1, sizeof(int));
1.621 + memset(&pass[0], 0, (m+1) * sizeof(int));
1.622 + val = xcalloc(1+n, sizeof(double));
1.623 + lb = xcalloc(1+n, sizeof(double));
1.624 + ub = xcalloc(1+n, sizeof(double));
1.625 + /* initialize the list of rows to be processed */
1.626 + size = 0;
1.627 + for (k = 1; k <= nrs; k++)
1.628 + { i = num[k];
1.629 + xassert(0 <= i && i <= m);
1.630 + /* duplicate row numbers are not allowed */
1.631 + xassert(!mark[i]);
1.632 + list[++size] = i, mark[i] = 1;
1.633 + }
1.634 + xassert(size == nrs);
1.635 + /* process rows in the list until it becomes empty */
1.636 + while (size > 0)
1.637 + { /* get a next row from the list */
1.638 + i = list[size--], mark[i] = 0;
1.639 + /* increase the row processing count */
1.640 + pass[i]++;
1.641 + /* if the row is free, skip it */
1.642 + if (L[i] == -DBL_MAX && U[i] == +DBL_MAX) continue;
1.643 + /* obtain coefficients of the row */
1.644 + len = 0;
1.645 + if (i == 0)
1.646 + { for (j = 1; j <= n; j++)
1.647 + { GLPCOL *col = mip->col[j];
1.648 + if (col->coef != 0.0)
1.649 + len++, ind[len] = j, val[len] = col->coef;
1.650 + }
1.651 + }
1.652 + else
1.653 + { GLPROW *row = mip->row[i];
1.654 + GLPAIJ *aij;
1.655 + for (aij = row->ptr; aij != NULL; aij = aij->r_next)
1.656 + len++, ind[len] = aij->col->j, val[len] = aij->val;
1.657 + }
1.658 + /* determine lower and upper bounds of columns corresponding
1.659 + to non-zero row coefficients */
1.660 + for (k = 1; k <= len; k++)
1.661 + j = ind[k], lb[k] = l[j], ub[k] = u[j];
1.662 + /* prepare the row info to determine implied bounds */
1.663 + prepare_row_info(len, val, lb, ub, &f);
1.664 + /* check and relax bounds of the row */
1.665 + if (check_row_bounds(&f, &L[i], &U[i]))
1.666 + { /* the feasible region is empty */
1.667 + ret = 1;
1.668 + goto done;
1.669 + }
1.670 + /* if the row became free, drop it */
1.671 + if (L[i] == -DBL_MAX && U[i] == +DBL_MAX) continue;
1.672 + /* process columns having non-zero coefficients in the row */
1.673 + for (k = 1; k <= len; k++)
1.674 + { GLPCOL *col;
1.675 + int flag, eff;
1.676 + double ll, uu;
1.677 + /* take a next column in the row */
1.678 + j = ind[k], col = mip->col[j];
1.679 + flag = col->kind != GLP_CV;
1.680 + /* check and tighten bounds of the column */
1.681 + if (check_col_bounds(&f, len, val, L[i], U[i], lb, ub,
1.682 + flag, k, &ll, &uu))
1.683 + { /* the feasible region is empty */
1.684 + ret = 1;
1.685 + goto done;
1.686 + }
1.687 + /* check if change in the column bounds is efficient */
1.688 + eff = check_efficiency(flag, l[j], u[j], ll, uu);
1.689 + /* set new actual bounds of the column */
1.690 + l[j] = ll, u[j] = uu;
1.691 + /* if the change is efficient, add all rows affected by the
1.692 + corresponding column, to the list */
1.693 + if (eff > 0)
1.694 + { GLPAIJ *aij;
1.695 + for (aij = col->ptr; aij != NULL; aij = aij->c_next)
1.696 + { int ii = aij->row->i;
1.697 + /* if the row was processed maximal number of times,
1.698 + skip it */
1.699 + if (pass[ii] >= max_pass) continue;
1.700 + /* if the row is free, skip it */
1.701 + if (L[ii] == -DBL_MAX && U[ii] == +DBL_MAX) continue;
1.702 + /* put the row into the list */
1.703 + if (mark[ii] == 0)
1.704 + { xassert(size <= m);
1.705 + list[++size] = ii, mark[ii] = 1;
1.706 + }
1.707 + }
1.708 + }
1.709 + }
1.710 + }
1.711 +done: /* free working arrays */
1.712 + xfree(ind);
1.713 + xfree(list);
1.714 + xfree(mark);
1.715 + xfree(pass);
1.716 + xfree(val);
1.717 + xfree(lb);
1.718 + xfree(ub);
1.719 + return ret;
1.720 +}
1.721 +
1.722 +/***********************************************************************
1.723 +* NAME
1.724 +*
1.725 +* ios_preprocess_node - preprocess current subproblem
1.726 +*
1.727 +* SYNOPSIS
1.728 +*
1.729 +* #include "glpios.h"
1.730 +* int ios_preprocess_node(glp_tree *tree, int max_pass);
1.731 +*
1.732 +* DESCRIPTION
1.733 +*
1.734 +* The routine ios_preprocess_node performs basic preprocessing of the
1.735 +* current subproblem.
1.736 +*
1.737 +* RETURNS
1.738 +*
1.739 +* If no primal infeasibility is detected, the routine returns zero,
1.740 +* otherwise non-zero. */
1.741 +
1.742 +int ios_preprocess_node(glp_tree *tree, int max_pass)
1.743 +{ glp_prob *mip = tree->mip;
1.744 + int m = mip->m;
1.745 + int n = mip->n;
1.746 + int i, j, nrs, *num, ret = 0;
1.747 + double *L, *U, *l, *u;
1.748 + /* the current subproblem must exist */
1.749 + xassert(tree->curr != NULL);
1.750 + /* determine original row bounds */
1.751 + L = xcalloc(1+m, sizeof(double));
1.752 + U = xcalloc(1+m, sizeof(double));
1.753 + switch (mip->mip_stat)
1.754 + { case GLP_UNDEF:
1.755 + L[0] = -DBL_MAX, U[0] = +DBL_MAX;
1.756 + break;
1.757 + case GLP_FEAS:
1.758 + switch (mip->dir)
1.759 + { case GLP_MIN:
1.760 + L[0] = -DBL_MAX, U[0] = mip->mip_obj - mip->c0;
1.761 + break;
1.762 + case GLP_MAX:
1.763 + L[0] = mip->mip_obj - mip->c0, U[0] = +DBL_MAX;
1.764 + break;
1.765 + default:
1.766 + xassert(mip != mip);
1.767 + }
1.768 + break;
1.769 + default:
1.770 + xassert(mip != mip);
1.771 + }
1.772 + for (i = 1; i <= m; i++)
1.773 + { L[i] = glp_get_row_lb(mip, i);
1.774 + U[i] = glp_get_row_ub(mip, i);
1.775 + }
1.776 + /* determine original column bounds */
1.777 + l = xcalloc(1+n, sizeof(double));
1.778 + u = xcalloc(1+n, sizeof(double));
1.779 + for (j = 1; j <= n; j++)
1.780 + { l[j] = glp_get_col_lb(mip, j);
1.781 + u[j] = glp_get_col_ub(mip, j);
1.782 + }
1.783 + /* build the initial list of rows to be analyzed */
1.784 + nrs = m + 1;
1.785 + num = xcalloc(1+nrs, sizeof(int));
1.786 + for (i = 1; i <= nrs; i++) num[i] = i - 1;
1.787 + /* perform basic preprocessing */
1.788 + if (basic_preprocessing(mip , L, U, l, u, nrs, num, max_pass))
1.789 + { ret = 1;
1.790 + goto done;
1.791 + }
1.792 + /* set new actual (relaxed) row bounds */
1.793 + for (i = 1; i <= m; i++)
1.794 + { /* consider only non-active rows to keep dual feasibility */
1.795 + if (glp_get_row_stat(mip, i) == GLP_BS)
1.796 + { if (L[i] == -DBL_MAX && U[i] == +DBL_MAX)
1.797 + glp_set_row_bnds(mip, i, GLP_FR, 0.0, 0.0);
1.798 + else if (U[i] == +DBL_MAX)
1.799 + glp_set_row_bnds(mip, i, GLP_LO, L[i], 0.0);
1.800 + else if (L[i] == -DBL_MAX)
1.801 + glp_set_row_bnds(mip, i, GLP_UP, 0.0, U[i]);
1.802 + }
1.803 + }
1.804 + /* set new actual (tightened) column bounds */
1.805 + for (j = 1; j <= n; j++)
1.806 + { int type;
1.807 + if (l[j] == -DBL_MAX && u[j] == +DBL_MAX)
1.808 + type = GLP_FR;
1.809 + else if (u[j] == +DBL_MAX)
1.810 + type = GLP_LO;
1.811 + else if (l[j] == -DBL_MAX)
1.812 + type = GLP_UP;
1.813 + else if (l[j] != u[j])
1.814 + type = GLP_DB;
1.815 + else
1.816 + type = GLP_FX;
1.817 + glp_set_col_bnds(mip, j, type, l[j], u[j]);
1.818 + }
1.819 +done: /* free working arrays and return */
1.820 + xfree(L);
1.821 + xfree(U);
1.822 + xfree(l);
1.823 + xfree(u);
1.824 + xfree(num);
1.825 + return ret;
1.826 +}
1.827 +
1.828 +/* eof */