alpar@1: /* glpios02.c (preprocess current subproblem) */ alpar@1: alpar@1: /*********************************************************************** alpar@1: * This code is part of GLPK (GNU Linear Programming Kit). alpar@1: * alpar@1: * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, alpar@1: * 2009, 2010 Andrew Makhorin, Department for Applied Informatics, alpar@1: * Moscow Aviation Institute, Moscow, Russia. All rights reserved. alpar@1: * E-mail: . alpar@1: * alpar@1: * GLPK is free software: you can redistribute it and/or modify it alpar@1: * under the terms of the GNU General Public License as published by alpar@1: * the Free Software Foundation, either version 3 of the License, or alpar@1: * (at your option) any later version. alpar@1: * alpar@1: * GLPK is distributed in the hope that it will be useful, but WITHOUT alpar@1: * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY alpar@1: * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public alpar@1: * License for more details. alpar@1: * alpar@1: * You should have received a copy of the GNU General Public License alpar@1: * along with GLPK. If not, see . alpar@1: ***********************************************************************/ alpar@1: alpar@1: #include "glpios.h" alpar@1: alpar@1: /*********************************************************************** alpar@1: * prepare_row_info - prepare row info to determine implied bounds alpar@1: * alpar@1: * Given a row (linear form) alpar@1: * alpar@1: * n alpar@1: * sum a[j] * x[j] (1) alpar@1: * j=1 alpar@1: * alpar@1: * and bounds of columns (variables) alpar@1: * alpar@1: * l[j] <= x[j] <= u[j] (2) alpar@1: * alpar@1: * this routine computes f_min, j_min, f_max, j_max needed to determine alpar@1: * implied bounds. alpar@1: * alpar@1: * ALGORITHM alpar@1: * alpar@1: * Let J+ = {j : a[j] > 0} and J- = {j : a[j] < 0}. alpar@1: * alpar@1: * Parameters f_min and j_min are computed as follows: alpar@1: * alpar@1: * 1) if there is no x[k] such that k in J+ and l[k] = -inf or k in J- alpar@1: * and u[k] = +inf, then alpar@1: * alpar@1: * f_min := sum a[j] * l[j] + sum a[j] * u[j] alpar@1: * j in J+ j in J- alpar@1: * (3) alpar@1: * j_min := 0 alpar@1: * alpar@1: * 2) if there is exactly one x[k] such that k in J+ and l[k] = -inf alpar@1: * or k in J- and u[k] = +inf, then alpar@1: * alpar@1: * f_min := sum a[j] * l[j] + sum a[j] * u[j] alpar@1: * j in J+\{k} j in J-\{k} alpar@1: * (4) alpar@1: * j_min := k alpar@1: * alpar@1: * 3) if there are two or more x[k] such that k in J+ and l[k] = -inf alpar@1: * or k in J- and u[k] = +inf, then alpar@1: * alpar@1: * f_min := -inf alpar@1: * (5) alpar@1: * j_min := 0 alpar@1: * alpar@1: * Parameters f_max and j_max are computed in a similar way as follows: alpar@1: * alpar@1: * 1) if there is no x[k] such that k in J+ and u[k] = +inf or k in J- alpar@1: * and l[k] = -inf, then alpar@1: * alpar@1: * f_max := sum a[j] * u[j] + sum a[j] * l[j] alpar@1: * j in J+ j in J- alpar@1: * (6) alpar@1: * j_max := 0 alpar@1: * alpar@1: * 2) if there is exactly one x[k] such that k in J+ and u[k] = +inf alpar@1: * or k in J- and l[k] = -inf, then alpar@1: * alpar@1: * f_max := sum a[j] * u[j] + sum a[j] * l[j] alpar@1: * j in J+\{k} j in J-\{k} alpar@1: * (7) alpar@1: * j_max := k alpar@1: * alpar@1: * 3) if there are two or more x[k] such that k in J+ and u[k] = +inf alpar@1: * or k in J- and l[k] = -inf, then alpar@1: * alpar@1: * f_max := +inf alpar@1: * (8) alpar@1: * j_max := 0 */ alpar@1: alpar@1: struct f_info alpar@1: { int j_min, j_max; alpar@1: double f_min, f_max; alpar@1: }; alpar@1: alpar@1: static void prepare_row_info(int n, const double a[], const double l[], alpar@1: const double u[], struct f_info *f) alpar@1: { int j, j_min, j_max; alpar@1: double f_min, f_max; alpar@1: xassert(n >= 0); alpar@1: /* determine f_min and j_min */ alpar@1: f_min = 0.0, j_min = 0; alpar@1: for (j = 1; j <= n; j++) alpar@1: { if (a[j] > 0.0) alpar@1: { if (l[j] == -DBL_MAX) alpar@1: { if (j_min == 0) alpar@1: j_min = j; alpar@1: else alpar@1: { f_min = -DBL_MAX, j_min = 0; alpar@1: break; alpar@1: } alpar@1: } alpar@1: else alpar@1: f_min += a[j] * l[j]; alpar@1: } alpar@1: else if (a[j] < 0.0) alpar@1: { if (u[j] == +DBL_MAX) alpar@1: { if (j_min == 0) alpar@1: j_min = j; alpar@1: else alpar@1: { f_min = -DBL_MAX, j_min = 0; alpar@1: break; alpar@1: } alpar@1: } alpar@1: else alpar@1: f_min += a[j] * u[j]; alpar@1: } alpar@1: else alpar@1: xassert(a != a); alpar@1: } alpar@1: f->f_min = f_min, f->j_min = j_min; alpar@1: /* determine f_max and j_max */ alpar@1: f_max = 0.0, j_max = 0; alpar@1: for (j = 1; j <= n; j++) alpar@1: { if (a[j] > 0.0) alpar@1: { if (u[j] == +DBL_MAX) alpar@1: { if (j_max == 0) alpar@1: j_max = j; alpar@1: else alpar@1: { f_max = +DBL_MAX, j_max = 0; alpar@1: break; alpar@1: } alpar@1: } alpar@1: else alpar@1: f_max += a[j] * u[j]; alpar@1: } alpar@1: else if (a[j] < 0.0) alpar@1: { if (l[j] == -DBL_MAX) alpar@1: { if (j_max == 0) alpar@1: j_max = j; alpar@1: else alpar@1: { f_max = +DBL_MAX, j_max = 0; alpar@1: break; alpar@1: } alpar@1: } alpar@1: else alpar@1: f_max += a[j] * l[j]; alpar@1: } alpar@1: else alpar@1: xassert(a != a); alpar@1: } alpar@1: f->f_max = f_max, f->j_max = j_max; alpar@1: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * row_implied_bounds - determine row implied bounds alpar@1: * alpar@1: * Given a row (linear form) alpar@1: * alpar@1: * n alpar@1: * sum a[j] * x[j] alpar@1: * j=1 alpar@1: * alpar@1: * and bounds of columns (variables) alpar@1: * alpar@1: * l[j] <= x[j] <= u[j] alpar@1: * alpar@1: * this routine determines implied bounds of the row. alpar@1: * alpar@1: * ALGORITHM alpar@1: * alpar@1: * Let J+ = {j : a[j] > 0} and J- = {j : a[j] < 0}. alpar@1: * alpar@1: * The implied lower bound of the row is computed as follows: alpar@1: * alpar@1: * L' := sum a[j] * l[j] + sum a[j] * u[j] (9) alpar@1: * j in J+ j in J- alpar@1: * alpar@1: * and as it follows from (3), (4), and (5): alpar@1: * alpar@1: * L' := if j_min = 0 then f_min else -inf (10) alpar@1: * alpar@1: * The implied upper bound of the row is computed as follows: alpar@1: * alpar@1: * U' := sum a[j] * u[j] + sum a[j] * l[j] (11) alpar@1: * j in J+ j in J- alpar@1: * alpar@1: * and as it follows from (6), (7), and (8): alpar@1: * alpar@1: * U' := if j_max = 0 then f_max else +inf (12) alpar@1: * alpar@1: * The implied bounds are stored in locations LL and UU. */ alpar@1: alpar@1: static void row_implied_bounds(const struct f_info *f, double *LL, alpar@1: double *UU) alpar@1: { *LL = (f->j_min == 0 ? f->f_min : -DBL_MAX); alpar@1: *UU = (f->j_max == 0 ? f->f_max : +DBL_MAX); alpar@1: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * col_implied_bounds - determine column implied bounds alpar@1: * alpar@1: * Given a row (constraint) alpar@1: * alpar@1: * n alpar@1: * L <= sum a[j] * x[j] <= U (13) alpar@1: * j=1 alpar@1: * alpar@1: * and bounds of columns (variables) alpar@1: * alpar@1: * l[j] <= x[j] <= u[j] alpar@1: * alpar@1: * this routine determines implied bounds of variable x[k]. alpar@1: * alpar@1: * It is assumed that if L != -inf, the lower bound of the row can be alpar@1: * active, and if U != +inf, the upper bound of the row can be active. alpar@1: * alpar@1: * ALGORITHM alpar@1: * alpar@1: * From (13) it follows that alpar@1: * alpar@1: * L <= sum a[j] * x[j] + a[k] * x[k] <= U alpar@1: * j!=k alpar@1: * or alpar@1: * alpar@1: * L - sum a[j] * x[j] <= a[k] * x[k] <= U - sum a[j] * x[j] alpar@1: * j!=k j!=k alpar@1: * alpar@1: * Thus, if the row lower bound L can be active, implied lower bound of alpar@1: * term a[k] * x[k] can be determined as follows: alpar@1: * alpar@1: * ilb(a[k] * x[k]) = min(L - sum a[j] * x[j]) = alpar@1: * j!=k alpar@1: * (14) alpar@1: * = L - max sum a[j] * x[j] alpar@1: * j!=k alpar@1: * alpar@1: * where, as it follows from (6), (7), and (8) alpar@1: * alpar@1: * / f_max - a[k] * u[k], j_max = 0, a[k] > 0 alpar@1: * | alpar@1: * | f_max - a[k] * l[k], j_max = 0, a[k] < 0 alpar@1: * max sum a[j] * x[j] = { alpar@1: * j!=k | f_max, j_max = k alpar@1: * | alpar@1: * \ +inf, j_max != 0 alpar@1: * alpar@1: * and if the upper bound U can be active, implied upper bound of term alpar@1: * a[k] * x[k] can be determined as follows: alpar@1: * alpar@1: * iub(a[k] * x[k]) = max(U - sum a[j] * x[j]) = alpar@1: * j!=k alpar@1: * (15) alpar@1: * = U - min sum a[j] * x[j] alpar@1: * j!=k alpar@1: * alpar@1: * where, as it follows from (3), (4), and (5) alpar@1: * alpar@1: * / f_min - a[k] * l[k], j_min = 0, a[k] > 0 alpar@1: * | alpar@1: * | f_min - a[k] * u[k], j_min = 0, a[k] < 0 alpar@1: * min sum a[j] * x[j] = { alpar@1: * j!=k | f_min, j_min = k alpar@1: * | alpar@1: * \ -inf, j_min != 0 alpar@1: * alpar@1: * Since alpar@1: * alpar@1: * ilb(a[k] * x[k]) <= a[k] * x[k] <= iub(a[k] * x[k]) alpar@1: * alpar@1: * implied lower and upper bounds of x[k] are determined as follows: alpar@1: * alpar@1: * l'[k] := if a[k] > 0 then ilb / a[k] else ulb / a[k] (16) alpar@1: * alpar@1: * u'[k] := if a[k] > 0 then ulb / a[k] else ilb / a[k] (17) alpar@1: * alpar@1: * The implied bounds are stored in locations ll and uu. */ alpar@1: alpar@1: static void col_implied_bounds(const struct f_info *f, int n, alpar@1: const double a[], double L, double U, const double l[], alpar@1: const double u[], int k, double *ll, double *uu) alpar@1: { double ilb, iub; alpar@1: xassert(n >= 0); alpar@1: xassert(1 <= k && k <= n); alpar@1: /* determine implied lower bound of term a[k] * x[k] (14) */ alpar@1: if (L == -DBL_MAX || f->f_max == +DBL_MAX) alpar@1: ilb = -DBL_MAX; alpar@1: else if (f->j_max == 0) alpar@1: { if (a[k] > 0.0) alpar@1: { xassert(u[k] != +DBL_MAX); alpar@1: ilb = L - (f->f_max - a[k] * u[k]); alpar@1: } alpar@1: else if (a[k] < 0.0) alpar@1: { xassert(l[k] != -DBL_MAX); alpar@1: ilb = L - (f->f_max - a[k] * l[k]); alpar@1: } alpar@1: else alpar@1: xassert(a != a); alpar@1: } alpar@1: else if (f->j_max == k) alpar@1: ilb = L - f->f_max; alpar@1: else alpar@1: ilb = -DBL_MAX; alpar@1: /* determine implied upper bound of term a[k] * x[k] (15) */ alpar@1: if (U == +DBL_MAX || f->f_min == -DBL_MAX) alpar@1: iub = +DBL_MAX; alpar@1: else if (f->j_min == 0) alpar@1: { if (a[k] > 0.0) alpar@1: { xassert(l[k] != -DBL_MAX); alpar@1: iub = U - (f->f_min - a[k] * l[k]); alpar@1: } alpar@1: else if (a[k] < 0.0) alpar@1: { xassert(u[k] != +DBL_MAX); alpar@1: iub = U - (f->f_min - a[k] * u[k]); alpar@1: } alpar@1: else alpar@1: xassert(a != a); alpar@1: } alpar@1: else if (f->j_min == k) alpar@1: iub = U - f->f_min; alpar@1: else alpar@1: iub = +DBL_MAX; alpar@1: /* determine implied bounds of x[k] (16) and (17) */ alpar@1: #if 1 alpar@1: /* do not use a[k] if it has small magnitude to prevent wrong alpar@1: implied bounds; for example, 1e-15 * x1 >= x2 + x3, where alpar@1: x1 >= -10, x2, x3 >= 0, would lead to wrong conclusion that alpar@1: x1 >= 0 */ alpar@1: if (fabs(a[k]) < 1e-6) alpar@1: *ll = -DBL_MAX, *uu = +DBL_MAX; else alpar@1: #endif alpar@1: if (a[k] > 0.0) alpar@1: { *ll = (ilb == -DBL_MAX ? -DBL_MAX : ilb / a[k]); alpar@1: *uu = (iub == +DBL_MAX ? +DBL_MAX : iub / a[k]); alpar@1: } alpar@1: else if (a[k] < 0.0) alpar@1: { *ll = (iub == +DBL_MAX ? -DBL_MAX : iub / a[k]); alpar@1: *uu = (ilb == -DBL_MAX ? +DBL_MAX : ilb / a[k]); alpar@1: } alpar@1: else alpar@1: xassert(a != a); alpar@1: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * check_row_bounds - check and relax original row bounds alpar@1: * alpar@1: * Given a row (constraint) alpar@1: * alpar@1: * n alpar@1: * L <= sum a[j] * x[j] <= U alpar@1: * j=1 alpar@1: * alpar@1: * and bounds of columns (variables) alpar@1: * alpar@1: * l[j] <= x[j] <= u[j] alpar@1: * alpar@1: * this routine checks the original row bounds L and U for feasibility alpar@1: * and redundancy. If the original lower bound L or/and upper bound U alpar@1: * cannot be active due to bounds of variables, the routine remove them alpar@1: * replacing by -inf or/and +inf, respectively. alpar@1: * alpar@1: * If no primal infeasibility is detected, the routine returns zero, alpar@1: * otherwise non-zero. */ alpar@1: alpar@1: static int check_row_bounds(const struct f_info *f, double *L_, alpar@1: double *U_) alpar@1: { int ret = 0; alpar@1: double L = *L_, U = *U_, LL, UU; alpar@1: /* determine implied bounds of the row */ alpar@1: row_implied_bounds(f, &LL, &UU); alpar@1: /* check if the original lower bound is infeasible */ alpar@1: if (L != -DBL_MAX) alpar@1: { double eps = 1e-3 * (1.0 + fabs(L)); alpar@1: if (UU < L - eps) alpar@1: { ret = 1; alpar@1: goto done; alpar@1: } alpar@1: } alpar@1: /* check if the original upper bound is infeasible */ alpar@1: if (U != +DBL_MAX) alpar@1: { double eps = 1e-3 * (1.0 + fabs(U)); alpar@1: if (LL > U + eps) alpar@1: { ret = 1; alpar@1: goto done; alpar@1: } alpar@1: } alpar@1: /* check if the original lower bound is redundant */ alpar@1: if (L != -DBL_MAX) alpar@1: { double eps = 1e-12 * (1.0 + fabs(L)); alpar@1: if (LL > L - eps) alpar@1: { /* it cannot be active, so remove it */ alpar@1: *L_ = -DBL_MAX; alpar@1: } alpar@1: } alpar@1: /* check if the original upper bound is redundant */ alpar@1: if (U != +DBL_MAX) alpar@1: { double eps = 1e-12 * (1.0 + fabs(U)); alpar@1: if (UU < U + eps) alpar@1: { /* it cannot be active, so remove it */ alpar@1: *U_ = +DBL_MAX; alpar@1: } alpar@1: } alpar@1: done: return ret; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * check_col_bounds - check and tighten original column bounds alpar@1: * alpar@1: * Given a row (constraint) alpar@1: * alpar@1: * n alpar@1: * L <= sum a[j] * x[j] <= U alpar@1: * j=1 alpar@1: * alpar@1: * and bounds of columns (variables) alpar@1: * alpar@1: * l[j] <= x[j] <= u[j] alpar@1: * alpar@1: * for column (variable) x[j] this routine checks the original column alpar@1: * bounds l[j] and u[j] for feasibility and redundancy. If the original alpar@1: * lower bound l[j] or/and upper bound u[j] cannot be active due to alpar@1: * bounds of the constraint and other variables, the routine tighten alpar@1: * them replacing by corresponding implied bounds, if possible. alpar@1: * alpar@1: * NOTE: It is assumed that if L != -inf, the row lower bound can be alpar@1: * active, and if U != +inf, the row upper bound can be active. alpar@1: * alpar@1: * The flag means that variable x[j] is required to be integer. alpar@1: * alpar@1: * New actual bounds for x[j] are stored in locations lj and uj. alpar@1: * alpar@1: * If no primal infeasibility is detected, the routine returns zero, alpar@1: * otherwise non-zero. */ alpar@1: alpar@1: static int check_col_bounds(const struct f_info *f, int n, alpar@1: const double a[], double L, double U, const double l[], alpar@1: const double u[], int flag, int j, double *_lj, double *_uj) alpar@1: { int ret = 0; alpar@1: double lj, uj, ll, uu; alpar@1: xassert(n >= 0); alpar@1: xassert(1 <= j && j <= n); alpar@1: lj = l[j], uj = u[j]; alpar@1: /* determine implied bounds of the column */ alpar@1: col_implied_bounds(f, n, a, L, U, l, u, j, &ll, &uu); alpar@1: /* if x[j] is integral, round its implied bounds */ alpar@1: if (flag) alpar@1: { if (ll != -DBL_MAX) alpar@1: ll = (ll - floor(ll) < 1e-3 ? floor(ll) : ceil(ll)); alpar@1: if (uu != +DBL_MAX) alpar@1: uu = (ceil(uu) - uu < 1e-3 ? ceil(uu) : floor(uu)); alpar@1: } alpar@1: /* check if the original lower bound is infeasible */ alpar@1: if (lj != -DBL_MAX) alpar@1: { double eps = 1e-3 * (1.0 + fabs(lj)); alpar@1: if (uu < lj - eps) alpar@1: { ret = 1; alpar@1: goto done; alpar@1: } alpar@1: } alpar@1: /* check if the original upper bound is infeasible */ alpar@1: if (uj != +DBL_MAX) alpar@1: { double eps = 1e-3 * (1.0 + fabs(uj)); alpar@1: if (ll > uj + eps) alpar@1: { ret = 1; alpar@1: goto done; alpar@1: } alpar@1: } alpar@1: /* check if the original lower bound is redundant */ alpar@1: if (ll != -DBL_MAX) alpar@1: { double eps = 1e-3 * (1.0 + fabs(ll)); alpar@1: if (lj < ll - eps) alpar@1: { /* it cannot be active, so tighten it */ alpar@1: lj = ll; alpar@1: } alpar@1: } alpar@1: /* check if the original upper bound is redundant */ alpar@1: if (uu != +DBL_MAX) alpar@1: { double eps = 1e-3 * (1.0 + fabs(uu)); alpar@1: if (uj > uu + eps) alpar@1: { /* it cannot be active, so tighten it */ alpar@1: uj = uu; alpar@1: } alpar@1: } alpar@1: /* due to round-off errors it may happen that lj > uj (although alpar@1: lj < uj + eps, since no primal infeasibility is detected), so alpar@1: adjuct the new actual bounds to provide lj <= uj */ alpar@1: if (!(lj == -DBL_MAX || uj == +DBL_MAX)) alpar@1: { double t1 = fabs(lj), t2 = fabs(uj); alpar@1: double eps = 1e-10 * (1.0 + (t1 <= t2 ? t1 : t2)); alpar@1: if (lj > uj - eps) alpar@1: { if (lj == l[j]) alpar@1: uj = lj; alpar@1: else if (uj == u[j]) alpar@1: lj = uj; alpar@1: else if (t1 <= t2) alpar@1: uj = lj; alpar@1: else alpar@1: lj = uj; alpar@1: } alpar@1: } alpar@1: *_lj = lj, *_uj = uj; alpar@1: done: return ret; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * check_efficiency - check if change in column bounds is efficient alpar@1: * alpar@1: * Given the original bounds of a column l and u and its new actual alpar@1: * bounds l' and u' (possibly tighten by the routine check_col_bounds) alpar@1: * this routine checks if the change in the column bounds is efficient alpar@1: * enough. If so, the routine returns non-zero, otherwise zero. alpar@1: * alpar@1: * The flag means that the variable is required to be integer. */ alpar@1: alpar@1: static int check_efficiency(int flag, double l, double u, double ll, alpar@1: double uu) alpar@1: { int eff = 0; alpar@1: /* check efficiency for lower bound */ alpar@1: if (l < ll) alpar@1: { if (flag || l == -DBL_MAX) alpar@1: eff++; alpar@1: else alpar@1: { double r; alpar@1: if (u == +DBL_MAX) alpar@1: r = 1.0 + fabs(l); alpar@1: else alpar@1: r = 1.0 + (u - l); alpar@1: if (ll - l >= 0.25 * r) alpar@1: eff++; alpar@1: } alpar@1: } alpar@1: /* check efficiency for upper bound */ alpar@1: if (u > uu) alpar@1: { if (flag || u == +DBL_MAX) alpar@1: eff++; alpar@1: else alpar@1: { double r; alpar@1: if (l == -DBL_MAX) alpar@1: r = 1.0 + fabs(u); alpar@1: else alpar@1: r = 1.0 + (u - l); alpar@1: if (u - uu >= 0.25 * r) alpar@1: eff++; alpar@1: } alpar@1: } alpar@1: return eff; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * basic_preprocessing - perform basic preprocessing alpar@1: * alpar@1: * This routine performs basic preprocessing of the specified MIP that alpar@1: * includes relaxing some row bounds and tightening some column bounds. alpar@1: * alpar@1: * On entry the arrays L and U contains original row bounds, and the alpar@1: * arrays l and u contains original column bounds: alpar@1: * alpar@1: * L[0] is the lower bound of the objective row; alpar@1: * L[i], i = 1,...,m, is the lower bound of i-th row; alpar@1: * U[0] is the upper bound of the objective row; alpar@1: * U[i], i = 1,...,m, is the upper bound of i-th row; alpar@1: * l[0] is not used; alpar@1: * l[j], j = 1,...,n, is the lower bound of j-th column; alpar@1: * u[0] is not used; alpar@1: * u[j], j = 1,...,n, is the upper bound of j-th column. alpar@1: * alpar@1: * On exit the arrays L, U, l, and u contain new actual bounds of rows alpar@1: * and column in the same locations. alpar@1: * alpar@1: * The parameters nrs and num specify an initial list of rows to be alpar@1: * processed: alpar@1: * alpar@1: * nrs is the number of rows in the initial list, 0 <= nrs <= m+1; alpar@1: * num[0] is not used; alpar@1: * num[1,...,nrs] are row numbers (0 means the objective row). alpar@1: * alpar@1: * The parameter max_pass specifies the maximal number of times that alpar@1: * each row can be processed, max_pass > 0. alpar@1: * alpar@1: * If no primal infeasibility is detected, the routine returns zero, alpar@1: * otherwise non-zero. */ alpar@1: alpar@1: static int basic_preprocessing(glp_prob *mip, double L[], double U[], alpar@1: double l[], double u[], int nrs, const int num[], int max_pass) alpar@1: { int m = mip->m; alpar@1: int n = mip->n; alpar@1: struct f_info f; alpar@1: int i, j, k, len, size, ret = 0; alpar@1: int *ind, *list, *mark, *pass; alpar@1: double *val, *lb, *ub; alpar@1: xassert(0 <= nrs && nrs <= m+1); alpar@1: xassert(max_pass > 0); alpar@1: /* allocate working arrays */ alpar@1: ind = xcalloc(1+n, sizeof(int)); alpar@1: list = xcalloc(1+m+1, sizeof(int)); alpar@1: mark = xcalloc(1+m+1, sizeof(int)); alpar@1: memset(&mark[0], 0, (m+1) * sizeof(int)); alpar@1: pass = xcalloc(1+m+1, sizeof(int)); alpar@1: memset(&pass[0], 0, (m+1) * sizeof(int)); alpar@1: val = xcalloc(1+n, sizeof(double)); alpar@1: lb = xcalloc(1+n, sizeof(double)); alpar@1: ub = xcalloc(1+n, sizeof(double)); alpar@1: /* initialize the list of rows to be processed */ alpar@1: size = 0; alpar@1: for (k = 1; k <= nrs; k++) alpar@1: { i = num[k]; alpar@1: xassert(0 <= i && i <= m); alpar@1: /* duplicate row numbers are not allowed */ alpar@1: xassert(!mark[i]); alpar@1: list[++size] = i, mark[i] = 1; alpar@1: } alpar@1: xassert(size == nrs); alpar@1: /* process rows in the list until it becomes empty */ alpar@1: while (size > 0) alpar@1: { /* get a next row from the list */ alpar@1: i = list[size--], mark[i] = 0; alpar@1: /* increase the row processing count */ alpar@1: pass[i]++; alpar@1: /* if the row is free, skip it */ alpar@1: if (L[i] == -DBL_MAX && U[i] == +DBL_MAX) continue; alpar@1: /* obtain coefficients of the row */ alpar@1: len = 0; alpar@1: if (i == 0) alpar@1: { for (j = 1; j <= n; j++) alpar@1: { GLPCOL *col = mip->col[j]; alpar@1: if (col->coef != 0.0) alpar@1: len++, ind[len] = j, val[len] = col->coef; alpar@1: } alpar@1: } alpar@1: else alpar@1: { GLPROW *row = mip->row[i]; alpar@1: GLPAIJ *aij; alpar@1: for (aij = row->ptr; aij != NULL; aij = aij->r_next) alpar@1: len++, ind[len] = aij->col->j, val[len] = aij->val; alpar@1: } alpar@1: /* determine lower and upper bounds of columns corresponding alpar@1: to non-zero row coefficients */ alpar@1: for (k = 1; k <= len; k++) alpar@1: j = ind[k], lb[k] = l[j], ub[k] = u[j]; alpar@1: /* prepare the row info to determine implied bounds */ alpar@1: prepare_row_info(len, val, lb, ub, &f); alpar@1: /* check and relax bounds of the row */ alpar@1: if (check_row_bounds(&f, &L[i], &U[i])) alpar@1: { /* the feasible region is empty */ alpar@1: ret = 1; alpar@1: goto done; alpar@1: } alpar@1: /* if the row became free, drop it */ alpar@1: if (L[i] == -DBL_MAX && U[i] == +DBL_MAX) continue; alpar@1: /* process columns having non-zero coefficients in the row */ alpar@1: for (k = 1; k <= len; k++) alpar@1: { GLPCOL *col; alpar@1: int flag, eff; alpar@1: double ll, uu; alpar@1: /* take a next column in the row */ alpar@1: j = ind[k], col = mip->col[j]; alpar@1: flag = col->kind != GLP_CV; alpar@1: /* check and tighten bounds of the column */ alpar@1: if (check_col_bounds(&f, len, val, L[i], U[i], lb, ub, alpar@1: flag, k, &ll, &uu)) alpar@1: { /* the feasible region is empty */ alpar@1: ret = 1; alpar@1: goto done; alpar@1: } alpar@1: /* check if change in the column bounds is efficient */ alpar@1: eff = check_efficiency(flag, l[j], u[j], ll, uu); alpar@1: /* set new actual bounds of the column */ alpar@1: l[j] = ll, u[j] = uu; alpar@1: /* if the change is efficient, add all rows affected by the alpar@1: corresponding column, to the list */ alpar@1: if (eff > 0) alpar@1: { GLPAIJ *aij; alpar@1: for (aij = col->ptr; aij != NULL; aij = aij->c_next) alpar@1: { int ii = aij->row->i; alpar@1: /* if the row was processed maximal number of times, alpar@1: skip it */ alpar@1: if (pass[ii] >= max_pass) continue; alpar@1: /* if the row is free, skip it */ alpar@1: if (L[ii] == -DBL_MAX && U[ii] == +DBL_MAX) continue; alpar@1: /* put the row into the list */ alpar@1: if (mark[ii] == 0) alpar@1: { xassert(size <= m); alpar@1: list[++size] = ii, mark[ii] = 1; alpar@1: } alpar@1: } alpar@1: } alpar@1: } alpar@1: } alpar@1: done: /* free working arrays */ alpar@1: xfree(ind); alpar@1: xfree(list); alpar@1: xfree(mark); alpar@1: xfree(pass); alpar@1: xfree(val); alpar@1: xfree(lb); alpar@1: xfree(ub); alpar@1: return ret; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * ios_preprocess_node - preprocess current subproblem alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * #include "glpios.h" alpar@1: * int ios_preprocess_node(glp_tree *tree, int max_pass); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * The routine ios_preprocess_node performs basic preprocessing of the alpar@1: * current subproblem. alpar@1: * alpar@1: * RETURNS alpar@1: * alpar@1: * If no primal infeasibility is detected, the routine returns zero, alpar@1: * otherwise non-zero. */ alpar@1: alpar@1: int ios_preprocess_node(glp_tree *tree, int max_pass) alpar@1: { glp_prob *mip = tree->mip; alpar@1: int m = mip->m; alpar@1: int n = mip->n; alpar@1: int i, j, nrs, *num, ret = 0; alpar@1: double *L, *U, *l, *u; alpar@1: /* the current subproblem must exist */ alpar@1: xassert(tree->curr != NULL); alpar@1: /* determine original row bounds */ alpar@1: L = xcalloc(1+m, sizeof(double)); alpar@1: U = xcalloc(1+m, sizeof(double)); alpar@1: switch (mip->mip_stat) alpar@1: { case GLP_UNDEF: alpar@1: L[0] = -DBL_MAX, U[0] = +DBL_MAX; alpar@1: break; alpar@1: case GLP_FEAS: alpar@1: switch (mip->dir) alpar@1: { case GLP_MIN: alpar@1: L[0] = -DBL_MAX, U[0] = mip->mip_obj - mip->c0; alpar@1: break; alpar@1: case GLP_MAX: alpar@1: L[0] = mip->mip_obj - mip->c0, U[0] = +DBL_MAX; alpar@1: break; alpar@1: default: alpar@1: xassert(mip != mip); alpar@1: } alpar@1: break; alpar@1: default: alpar@1: xassert(mip != mip); alpar@1: } alpar@1: for (i = 1; i <= m; i++) alpar@1: { L[i] = glp_get_row_lb(mip, i); alpar@1: U[i] = glp_get_row_ub(mip, i); alpar@1: } alpar@1: /* determine original column bounds */ alpar@1: l = xcalloc(1+n, sizeof(double)); alpar@1: u = xcalloc(1+n, sizeof(double)); alpar@1: for (j = 1; j <= n; j++) alpar@1: { l[j] = glp_get_col_lb(mip, j); alpar@1: u[j] = glp_get_col_ub(mip, j); alpar@1: } alpar@1: /* build the initial list of rows to be analyzed */ alpar@1: nrs = m + 1; alpar@1: num = xcalloc(1+nrs, sizeof(int)); alpar@1: for (i = 1; i <= nrs; i++) num[i] = i - 1; alpar@1: /* perform basic preprocessing */ alpar@1: if (basic_preprocessing(mip , L, U, l, u, nrs, num, max_pass)) alpar@1: { ret = 1; alpar@1: goto done; alpar@1: } alpar@1: /* set new actual (relaxed) row bounds */ alpar@1: for (i = 1; i <= m; i++) alpar@1: { /* consider only non-active rows to keep dual feasibility */ alpar@1: if (glp_get_row_stat(mip, i) == GLP_BS) alpar@1: { if (L[i] == -DBL_MAX && U[i] == +DBL_MAX) alpar@1: glp_set_row_bnds(mip, i, GLP_FR, 0.0, 0.0); alpar@1: else if (U[i] == +DBL_MAX) alpar@1: glp_set_row_bnds(mip, i, GLP_LO, L[i], 0.0); alpar@1: else if (L[i] == -DBL_MAX) alpar@1: glp_set_row_bnds(mip, i, GLP_UP, 0.0, U[i]); alpar@1: } alpar@1: } alpar@1: /* set new actual (tightened) column bounds */ alpar@1: for (j = 1; j <= n; j++) alpar@1: { int type; alpar@1: if (l[j] == -DBL_MAX && u[j] == +DBL_MAX) alpar@1: type = GLP_FR; alpar@1: else if (u[j] == +DBL_MAX) alpar@1: type = GLP_LO; alpar@1: else if (l[j] == -DBL_MAX) alpar@1: type = GLP_UP; alpar@1: else if (l[j] != u[j]) alpar@1: type = GLP_DB; alpar@1: else alpar@1: type = GLP_FX; alpar@1: glp_set_col_bnds(mip, j, type, l[j], u[j]); alpar@1: } alpar@1: done: /* free working arrays and return */ alpar@1: xfree(L); alpar@1: xfree(U); alpar@1: xfree(l); alpar@1: xfree(u); alpar@1: xfree(num); alpar@1: return ret; alpar@1: } alpar@1: alpar@1: /* eof */