1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpnpp03.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,2861 @@
1.4 +/* glpnpp03.c */
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 "glpnpp.h"
1.29 +
1.30 +/***********************************************************************
1.31 +* NAME
1.32 +*
1.33 +* npp_empty_row - process empty row
1.34 +*
1.35 +* SYNOPSIS
1.36 +*
1.37 +* #include "glpnpp.h"
1.38 +* int npp_empty_row(NPP *npp, NPPROW *p);
1.39 +*
1.40 +* DESCRIPTION
1.41 +*
1.42 +* The routine npp_empty_row processes row p, which is empty, i.e.
1.43 +* coefficients at all columns in this row are zero:
1.44 +*
1.45 +* L[p] <= sum 0 x[j] <= U[p], (1)
1.46 +*
1.47 +* where L[p] <= U[p].
1.48 +*
1.49 +* RETURNS
1.50 +*
1.51 +* 0 - success;
1.52 +*
1.53 +* 1 - problem has no primal feasible solution.
1.54 +*
1.55 +* PROBLEM TRANSFORMATION
1.56 +*
1.57 +* If the following conditions hold:
1.58 +*
1.59 +* L[p] <= +eps, U[p] >= -eps, (2)
1.60 +*
1.61 +* where eps is an absolute tolerance for row value, the row p is
1.62 +* redundant. In this case it can be replaced by equivalent redundant
1.63 +* row, which is free (unbounded), and then removed from the problem.
1.64 +* Otherwise, the row p is infeasible and, thus, the problem has no
1.65 +* primal feasible solution.
1.66 +*
1.67 +* RECOVERING BASIC SOLUTION
1.68 +*
1.69 +* See the routine npp_free_row.
1.70 +*
1.71 +* RECOVERING INTERIOR-POINT SOLUTION
1.72 +*
1.73 +* See the routine npp_free_row.
1.74 +*
1.75 +* RECOVERING MIP SOLUTION
1.76 +*
1.77 +* None needed. */
1.78 +
1.79 +int npp_empty_row(NPP *npp, NPPROW *p)
1.80 +{ /* process empty row */
1.81 + double eps = 1e-3;
1.82 + /* the row must be empty */
1.83 + xassert(p->ptr == NULL);
1.84 + /* check primal feasibility */
1.85 + if (p->lb > +eps || p->ub < -eps)
1.86 + return 1;
1.87 + /* replace the row by equivalent free (unbounded) row */
1.88 + p->lb = -DBL_MAX, p->ub = +DBL_MAX;
1.89 + /* and process it */
1.90 + npp_free_row(npp, p);
1.91 + return 0;
1.92 +}
1.93 +
1.94 +/***********************************************************************
1.95 +* NAME
1.96 +*
1.97 +* npp_empty_col - process empty column
1.98 +*
1.99 +* SYNOPSIS
1.100 +*
1.101 +* #include "glpnpp.h"
1.102 +* int npp_empty_col(NPP *npp, NPPCOL *q);
1.103 +*
1.104 +* DESCRIPTION
1.105 +*
1.106 +* The routine npp_empty_col processes column q:
1.107 +*
1.108 +* l[q] <= x[q] <= u[q], (1)
1.109 +*
1.110 +* where l[q] <= u[q], which is empty, i.e. has zero coefficients in
1.111 +* all constraint rows.
1.112 +*
1.113 +* RETURNS
1.114 +*
1.115 +* 0 - success;
1.116 +*
1.117 +* 1 - problem has no dual feasible solution.
1.118 +*
1.119 +* PROBLEM TRANSFORMATION
1.120 +*
1.121 +* The row of the dual system corresponding to the empty column is the
1.122 +* following:
1.123 +*
1.124 +* sum 0 pi[i] + lambda[q] = c[q], (2)
1.125 +* i
1.126 +*
1.127 +* from which it follows that:
1.128 +*
1.129 +* lambda[q] = c[q]. (3)
1.130 +*
1.131 +* If the following condition holds:
1.132 +*
1.133 +* c[q] < - eps, (4)
1.134 +*
1.135 +* where eps is an absolute tolerance for column multiplier, the lower
1.136 +* column bound l[q] must be active to provide dual feasibility (note
1.137 +* that being preprocessed the problem is always minimization). In this
1.138 +* case the column can be fixed on its lower bound and removed from the
1.139 +* problem (if the column is integral, its bounds are also assumed to
1.140 +* be integral). And if the column has no lower bound (l[q] = -oo), the
1.141 +* problem has no dual feasible solution.
1.142 +*
1.143 +* If the following condition holds:
1.144 +*
1.145 +* c[q] > + eps, (5)
1.146 +*
1.147 +* the upper column bound u[q] must be active to provide dual
1.148 +* feasibility. In this case the column can be fixed on its upper bound
1.149 +* and removed from the problem. And if the column has no upper bound
1.150 +* (u[q] = +oo), the problem has no dual feasible solution.
1.151 +*
1.152 +* Finally, if the following condition holds:
1.153 +*
1.154 +* - eps <= c[q] <= +eps, (6)
1.155 +*
1.156 +* dual feasibility does not depend on a particular value of column q.
1.157 +* In this case the column can be fixed either on its lower bound (if
1.158 +* l[q] > -oo) or on its upper bound (if u[q] < +oo) or at zero (if the
1.159 +* column is unbounded) and then removed from the problem.
1.160 +*
1.161 +* RECOVERING BASIC SOLUTION
1.162 +*
1.163 +* See the routine npp_fixed_col. Having been recovered the column
1.164 +* is assigned status GLP_NS. However, if actually it is not fixed
1.165 +* (l[q] < u[q]), its status should be changed to GLP_NL, GLP_NU, or
1.166 +* GLP_NF depending on which bound it was fixed on transformation stage.
1.167 +*
1.168 +* RECOVERING INTERIOR-POINT SOLUTION
1.169 +*
1.170 +* See the routine npp_fixed_col.
1.171 +*
1.172 +* RECOVERING MIP SOLUTION
1.173 +*
1.174 +* See the routine npp_fixed_col. */
1.175 +
1.176 +struct empty_col
1.177 +{ /* empty column */
1.178 + int q;
1.179 + /* column reference number */
1.180 + char stat;
1.181 + /* status in basic solution */
1.182 +};
1.183 +
1.184 +static int rcv_empty_col(NPP *npp, void *info);
1.185 +
1.186 +int npp_empty_col(NPP *npp, NPPCOL *q)
1.187 +{ /* process empty column */
1.188 + struct empty_col *info;
1.189 + double eps = 1e-3;
1.190 + /* the column must be empty */
1.191 + xassert(q->ptr == NULL);
1.192 + /* check dual feasibility */
1.193 + if (q->coef > +eps && q->lb == -DBL_MAX)
1.194 + return 1;
1.195 + if (q->coef < -eps && q->ub == +DBL_MAX)
1.196 + return 1;
1.197 + /* create transformation stack entry */
1.198 + info = npp_push_tse(npp,
1.199 + rcv_empty_col, sizeof(struct empty_col));
1.200 + info->q = q->j;
1.201 + /* fix the column */
1.202 + if (q->lb == -DBL_MAX && q->ub == +DBL_MAX)
1.203 + { /* free column */
1.204 + info->stat = GLP_NF;
1.205 + q->lb = q->ub = 0.0;
1.206 + }
1.207 + else if (q->ub == +DBL_MAX)
1.208 +lo: { /* column with lower bound */
1.209 + info->stat = GLP_NL;
1.210 + q->ub = q->lb;
1.211 + }
1.212 + else if (q->lb == -DBL_MAX)
1.213 +up: { /* column with upper bound */
1.214 + info->stat = GLP_NU;
1.215 + q->lb = q->ub;
1.216 + }
1.217 + else if (q->lb != q->ub)
1.218 + { /* double-bounded column */
1.219 + if (q->coef >= +DBL_EPSILON) goto lo;
1.220 + if (q->coef <= -DBL_EPSILON) goto up;
1.221 + if (fabs(q->lb) <= fabs(q->ub)) goto lo; else goto up;
1.222 + }
1.223 + else
1.224 + { /* fixed column */
1.225 + info->stat = GLP_NS;
1.226 + }
1.227 + /* process fixed column */
1.228 + npp_fixed_col(npp, q);
1.229 + return 0;
1.230 +}
1.231 +
1.232 +static int rcv_empty_col(NPP *npp, void *_info)
1.233 +{ /* recover empty column */
1.234 + struct empty_col *info = _info;
1.235 + if (npp->sol == GLP_SOL)
1.236 + npp->c_stat[info->q] = info->stat;
1.237 + return 0;
1.238 +}
1.239 +
1.240 +/***********************************************************************
1.241 +* NAME
1.242 +*
1.243 +* npp_implied_value - process implied column value
1.244 +*
1.245 +* SYNOPSIS
1.246 +*
1.247 +* #include "glpnpp.h"
1.248 +* int npp_implied_value(NPP *npp, NPPCOL *q, double s);
1.249 +*
1.250 +* DESCRIPTION
1.251 +*
1.252 +* For column q:
1.253 +*
1.254 +* l[q] <= x[q] <= u[q], (1)
1.255 +*
1.256 +* where l[q] < u[q], the routine npp_implied_value processes its
1.257 +* implied value s[q]. If this implied value satisfies to the current
1.258 +* column bounds and integrality condition, the routine fixes column q
1.259 +* at the given point. Note that the column is kept in the problem in
1.260 +* any case.
1.261 +*
1.262 +* RETURNS
1.263 +*
1.264 +* 0 - column has been fixed;
1.265 +*
1.266 +* 1 - implied value violates to current column bounds;
1.267 +*
1.268 +* 2 - implied value violates integrality condition.
1.269 +*
1.270 +* ALGORITHM
1.271 +*
1.272 +* Implied column value s[q] satisfies to the current column bounds if
1.273 +* the following condition holds:
1.274 +*
1.275 +* l[q] - eps <= s[q] <= u[q] + eps, (2)
1.276 +*
1.277 +* where eps is an absolute tolerance for column value. If the column
1.278 +* is integral, the following condition also must hold:
1.279 +*
1.280 +* |s[q] - floor(s[q]+0.5)| <= eps, (3)
1.281 +*
1.282 +* where floor(s[q]+0.5) is the nearest integer to s[q].
1.283 +*
1.284 +* If both condition (2) and (3) are satisfied, the column can be fixed
1.285 +* at the value s[q], or, if it is integral, at floor(s[q]+0.5).
1.286 +* Otherwise, if s[q] violates (2) or (3), the problem has no feasible
1.287 +* solution.
1.288 +*
1.289 +* Note: If s[q] is close to l[q] or u[q], it seems to be reasonable to
1.290 +* fix the column at its lower or upper bound, resp. rather than at the
1.291 +* implied value. */
1.292 +
1.293 +int npp_implied_value(NPP *npp, NPPCOL *q, double s)
1.294 +{ /* process implied column value */
1.295 + double eps, nint;
1.296 + xassert(npp == npp);
1.297 + /* column must not be fixed */
1.298 + xassert(q->lb < q->ub);
1.299 + /* check integrality */
1.300 + if (q->is_int)
1.301 + { nint = floor(s + 0.5);
1.302 + if (fabs(s - nint) <= 1e-5)
1.303 + s = nint;
1.304 + else
1.305 + return 2;
1.306 + }
1.307 + /* check current column lower bound */
1.308 + if (q->lb != -DBL_MAX)
1.309 + { eps = (q->is_int ? 1e-5 : 1e-5 + 1e-8 * fabs(q->lb));
1.310 + if (s < q->lb - eps) return 1;
1.311 + /* if s[q] is close to l[q], fix column at its lower bound
1.312 + rather than at the implied value */
1.313 + if (s < q->lb + 1e-3 * eps)
1.314 + { q->ub = q->lb;
1.315 + return 0;
1.316 + }
1.317 + }
1.318 + /* check current column upper bound */
1.319 + if (q->ub != +DBL_MAX)
1.320 + { eps = (q->is_int ? 1e-5 : 1e-5 + 1e-8 * fabs(q->ub));
1.321 + if (s > q->ub + eps) return 1;
1.322 + /* if s[q] is close to u[q], fix column at its upper bound
1.323 + rather than at the implied value */
1.324 + if (s > q->ub - 1e-3 * eps)
1.325 + { q->lb = q->ub;
1.326 + return 0;
1.327 + }
1.328 + }
1.329 + /* fix column at the implied value */
1.330 + q->lb = q->ub = s;
1.331 + return 0;
1.332 +}
1.333 +
1.334 +/***********************************************************************
1.335 +* NAME
1.336 +*
1.337 +* npp_eq_singlet - process row singleton (equality constraint)
1.338 +*
1.339 +* SYNOPSIS
1.340 +*
1.341 +* #include "glpnpp.h"
1.342 +* int npp_eq_singlet(NPP *npp, NPPROW *p);
1.343 +*
1.344 +* DESCRIPTION
1.345 +*
1.346 +* The routine npp_eq_singlet processes row p, which is equiality
1.347 +* constraint having the only non-zero coefficient:
1.348 +*
1.349 +* a[p,q] x[q] = b. (1)
1.350 +*
1.351 +* RETURNS
1.352 +*
1.353 +* 0 - success;
1.354 +*
1.355 +* 1 - problem has no primal feasible solution;
1.356 +*
1.357 +* 2 - problem has no integer feasible solution.
1.358 +*
1.359 +* PROBLEM TRANSFORMATION
1.360 +*
1.361 +* The equality constraint defines implied value of column q:
1.362 +*
1.363 +* x[q] = s[q] = b / a[p,q]. (2)
1.364 +*
1.365 +* If the implied value s[q] satisfies to the column bounds (see the
1.366 +* routine npp_implied_value), the column can be fixed at s[q] and
1.367 +* removed from the problem. In this case row p becomes redundant, so
1.368 +* it can be replaced by equivalent free row and also removed from the
1.369 +* problem.
1.370 +*
1.371 +* Note that the routine removes from the problem only row p. Column q
1.372 +* becomes fixed, however, it is kept in the problem.
1.373 +*
1.374 +* RECOVERING BASIC SOLUTION
1.375 +*
1.376 +* In solution to the original problem row p is assigned status GLP_NS
1.377 +* (active equality constraint), and column q is assigned status GLP_BS
1.378 +* (basic column).
1.379 +*
1.380 +* Multiplier for row p can be computed as follows. In the dual system
1.381 +* of the original problem column q corresponds to the following row:
1.382 +*
1.383 +* sum a[i,q] pi[i] + lambda[q] = c[q] ==>
1.384 +* i
1.385 +*
1.386 +* sum a[i,q] pi[i] + a[p,q] pi[p] + lambda[q] = c[q].
1.387 +* i!=p
1.388 +*
1.389 +* Therefore:
1.390 +*
1.391 +* 1
1.392 +* pi[p] = ------ (c[q] - lambda[q] - sum a[i,q] pi[i]), (3)
1.393 +* a[p,q] i!=q
1.394 +*
1.395 +* where lambda[q] = 0 (since column[q] is basic), and pi[i] for all
1.396 +* i != p are known in solution to the transformed problem.
1.397 +*
1.398 +* Value of column q in solution to the original problem is assigned
1.399 +* its implied value s[q].
1.400 +*
1.401 +* RECOVERING INTERIOR-POINT SOLUTION
1.402 +*
1.403 +* Multiplier for row p is computed with formula (3). Value of column
1.404 +* q is assigned its implied value s[q].
1.405 +*
1.406 +* RECOVERING MIP SOLUTION
1.407 +*
1.408 +* Value of column q is assigned its implied value s[q]. */
1.409 +
1.410 +struct eq_singlet
1.411 +{ /* row singleton (equality constraint) */
1.412 + int p;
1.413 + /* row reference number */
1.414 + int q;
1.415 + /* column reference number */
1.416 + double apq;
1.417 + /* constraint coefficient a[p,q] */
1.418 + double c;
1.419 + /* objective coefficient at x[q] */
1.420 + NPPLFE *ptr;
1.421 + /* list of non-zero coefficients a[i,q], i != p */
1.422 +};
1.423 +
1.424 +static int rcv_eq_singlet(NPP *npp, void *info);
1.425 +
1.426 +int npp_eq_singlet(NPP *npp, NPPROW *p)
1.427 +{ /* process row singleton (equality constraint) */
1.428 + struct eq_singlet *info;
1.429 + NPPCOL *q;
1.430 + NPPAIJ *aij;
1.431 + NPPLFE *lfe;
1.432 + int ret;
1.433 + double s;
1.434 + /* the row must be singleton equality constraint */
1.435 + xassert(p->lb == p->ub);
1.436 + xassert(p->ptr != NULL && p->ptr->r_next == NULL);
1.437 + /* compute and process implied column value */
1.438 + aij = p->ptr;
1.439 + q = aij->col;
1.440 + s = p->lb / aij->val;
1.441 + ret = npp_implied_value(npp, q, s);
1.442 + xassert(0 <= ret && ret <= 2);
1.443 + if (ret != 0) return ret;
1.444 + /* create transformation stack entry */
1.445 + info = npp_push_tse(npp,
1.446 + rcv_eq_singlet, sizeof(struct eq_singlet));
1.447 + info->p = p->i;
1.448 + info->q = q->j;
1.449 + info->apq = aij->val;
1.450 + info->c = q->coef;
1.451 + info->ptr = NULL;
1.452 + /* save column coefficients a[i,q], i != p (not needed for MIP
1.453 + solution) */
1.454 + if (npp->sol != GLP_MIP)
1.455 + { for (aij = q->ptr; aij != NULL; aij = aij->c_next)
1.456 + { if (aij->row == p) continue; /* skip a[p,q] */
1.457 + lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE));
1.458 + lfe->ref = aij->row->i;
1.459 + lfe->val = aij->val;
1.460 + lfe->next = info->ptr;
1.461 + info->ptr = lfe;
1.462 + }
1.463 + }
1.464 + /* remove the row from the problem */
1.465 + npp_del_row(npp, p);
1.466 + return 0;
1.467 +}
1.468 +
1.469 +static int rcv_eq_singlet(NPP *npp, void *_info)
1.470 +{ /* recover row singleton (equality constraint) */
1.471 + struct eq_singlet *info = _info;
1.472 + NPPLFE *lfe;
1.473 + double temp;
1.474 + if (npp->sol == GLP_SOL)
1.475 + { /* column q must be already recovered as GLP_NS */
1.476 + if (npp->c_stat[info->q] != GLP_NS)
1.477 + { npp_error();
1.478 + return 1;
1.479 + }
1.480 + npp->r_stat[info->p] = GLP_NS;
1.481 + npp->c_stat[info->q] = GLP_BS;
1.482 + }
1.483 + if (npp->sol != GLP_MIP)
1.484 + { /* compute multiplier for row p with formula (3) */
1.485 + temp = info->c;
1.486 + for (lfe = info->ptr; lfe != NULL; lfe = lfe->next)
1.487 + temp -= lfe->val * npp->r_pi[lfe->ref];
1.488 + npp->r_pi[info->p] = temp / info->apq;
1.489 + }
1.490 + return 0;
1.491 +}
1.492 +
1.493 +/***********************************************************************
1.494 +* NAME
1.495 +*
1.496 +* npp_implied_lower - process implied column lower bound
1.497 +*
1.498 +* SYNOPSIS
1.499 +*
1.500 +* #include "glpnpp.h"
1.501 +* int npp_implied_lower(NPP *npp, NPPCOL *q, double l);
1.502 +*
1.503 +* DESCRIPTION
1.504 +*
1.505 +* For column q:
1.506 +*
1.507 +* l[q] <= x[q] <= u[q], (1)
1.508 +*
1.509 +* where l[q] < u[q], the routine npp_implied_lower processes its
1.510 +* implied lower bound l'[q]. As the result the current column lower
1.511 +* bound may increase. Note that the column is kept in the problem in
1.512 +* any case.
1.513 +*
1.514 +* RETURNS
1.515 +*
1.516 +* 0 - current column lower bound has not changed;
1.517 +*
1.518 +* 1 - current column lower bound has changed, but not significantly;
1.519 +*
1.520 +* 2 - current column lower bound has significantly changed;
1.521 +*
1.522 +* 3 - column has been fixed on its upper bound;
1.523 +*
1.524 +* 4 - implied lower bound violates current column upper bound.
1.525 +*
1.526 +* ALGORITHM
1.527 +*
1.528 +* If column q is integral, before processing its implied lower bound
1.529 +* should be rounded up:
1.530 +*
1.531 +* ( floor(l'[q]+0.5), if |l'[q] - floor(l'[q]+0.5)| <= eps
1.532 +* l'[q] := < (2)
1.533 +* ( ceil(l'[q]), otherwise
1.534 +*
1.535 +* where floor(l'[q]+0.5) is the nearest integer to l'[q], ceil(l'[q])
1.536 +* is smallest integer not less than l'[q], and eps is an absolute
1.537 +* tolerance for column value.
1.538 +*
1.539 +* Processing implied column lower bound l'[q] includes the following
1.540 +* cases:
1.541 +*
1.542 +* 1) if l'[q] < l[q] + eps, implied lower bound is redundant;
1.543 +*
1.544 +* 2) if l[q] + eps <= l[q] <= u[q] + eps, current column lower bound
1.545 +* l[q] can be strengthened by replacing it with l'[q]. If in this
1.546 +* case new column lower bound becomes close to current column upper
1.547 +* bound u[q], the column can be fixed on its upper bound;
1.548 +*
1.549 +* 3) if l'[q] > u[q] + eps, implied lower bound violates current
1.550 +* column upper bound u[q], in which case the problem has no primal
1.551 +* feasible solution. */
1.552 +
1.553 +int npp_implied_lower(NPP *npp, NPPCOL *q, double l)
1.554 +{ /* process implied column lower bound */
1.555 + int ret;
1.556 + double eps, nint;
1.557 + xassert(npp == npp);
1.558 + /* column must not be fixed */
1.559 + xassert(q->lb < q->ub);
1.560 + /* implied lower bound must be finite */
1.561 + xassert(l != -DBL_MAX);
1.562 + /* if column is integral, round up l'[q] */
1.563 + if (q->is_int)
1.564 + { nint = floor(l + 0.5);
1.565 + if (fabs(l - nint) <= 1e-5)
1.566 + l = nint;
1.567 + else
1.568 + l = ceil(l);
1.569 + }
1.570 + /* check current column lower bound */
1.571 + if (q->lb != -DBL_MAX)
1.572 + { eps = (q->is_int ? 1e-3 : 1e-3 + 1e-6 * fabs(q->lb));
1.573 + if (l < q->lb + eps)
1.574 + { ret = 0; /* redundant */
1.575 + goto done;
1.576 + }
1.577 + }
1.578 + /* check current column upper bound */
1.579 + if (q->ub != +DBL_MAX)
1.580 + { eps = (q->is_int ? 1e-5 : 1e-5 + 1e-8 * fabs(q->ub));
1.581 + if (l > q->ub + eps)
1.582 + { ret = 4; /* infeasible */
1.583 + goto done;
1.584 + }
1.585 + /* if l'[q] is close to u[q], fix column at its upper bound */
1.586 + if (l > q->ub - 1e-3 * eps)
1.587 + { q->lb = q->ub;
1.588 + ret = 3; /* fixed */
1.589 + goto done;
1.590 + }
1.591 + }
1.592 + /* check if column lower bound changes significantly */
1.593 + if (q->lb == -DBL_MAX)
1.594 + ret = 2; /* significantly */
1.595 + else if (q->is_int && l > q->lb + 0.5)
1.596 + ret = 2; /* significantly */
1.597 + else if (l > q->lb + 0.30 * (1.0 + fabs(q->lb)))
1.598 + ret = 2; /* significantly */
1.599 + else
1.600 + ret = 1; /* not significantly */
1.601 + /* set new column lower bound */
1.602 + q->lb = l;
1.603 +done: return ret;
1.604 +}
1.605 +
1.606 +/***********************************************************************
1.607 +* NAME
1.608 +*
1.609 +* npp_implied_upper - process implied column upper bound
1.610 +*
1.611 +* SYNOPSIS
1.612 +*
1.613 +* #include "glpnpp.h"
1.614 +* int npp_implied_upper(NPP *npp, NPPCOL *q, double u);
1.615 +*
1.616 +* DESCRIPTION
1.617 +*
1.618 +* For column q:
1.619 +*
1.620 +* l[q] <= x[q] <= u[q], (1)
1.621 +*
1.622 +* where l[q] < u[q], the routine npp_implied_upper processes its
1.623 +* implied upper bound u'[q]. As the result the current column upper
1.624 +* bound may decrease. Note that the column is kept in the problem in
1.625 +* any case.
1.626 +*
1.627 +* RETURNS
1.628 +*
1.629 +* 0 - current column upper bound has not changed;
1.630 +*
1.631 +* 1 - current column upper bound has changed, but not significantly;
1.632 +*
1.633 +* 2 - current column upper bound has significantly changed;
1.634 +*
1.635 +* 3 - column has been fixed on its lower bound;
1.636 +*
1.637 +* 4 - implied upper bound violates current column lower bound.
1.638 +*
1.639 +* ALGORITHM
1.640 +*
1.641 +* If column q is integral, before processing its implied upper bound
1.642 +* should be rounded down:
1.643 +*
1.644 +* ( floor(u'[q]+0.5), if |u'[q] - floor(l'[q]+0.5)| <= eps
1.645 +* u'[q] := < (2)
1.646 +* ( floor(l'[q]), otherwise
1.647 +*
1.648 +* where floor(u'[q]+0.5) is the nearest integer to u'[q],
1.649 +* floor(u'[q]) is largest integer not greater than u'[q], and eps is
1.650 +* an absolute tolerance for column value.
1.651 +*
1.652 +* Processing implied column upper bound u'[q] includes the following
1.653 +* cases:
1.654 +*
1.655 +* 1) if u'[q] > u[q] - eps, implied upper bound is redundant;
1.656 +*
1.657 +* 2) if l[q] - eps <= u[q] <= u[q] - eps, current column upper bound
1.658 +* u[q] can be strengthened by replacing it with u'[q]. If in this
1.659 +* case new column upper bound becomes close to current column lower
1.660 +* bound, the column can be fixed on its lower bound;
1.661 +*
1.662 +* 3) if u'[q] < l[q] - eps, implied upper bound violates current
1.663 +* column lower bound l[q], in which case the problem has no primal
1.664 +* feasible solution. */
1.665 +
1.666 +int npp_implied_upper(NPP *npp, NPPCOL *q, double u)
1.667 +{ int ret;
1.668 + double eps, nint;
1.669 + xassert(npp == npp);
1.670 + /* column must not be fixed */
1.671 + xassert(q->lb < q->ub);
1.672 + /* implied upper bound must be finite */
1.673 + xassert(u != +DBL_MAX);
1.674 + /* if column is integral, round down u'[q] */
1.675 + if (q->is_int)
1.676 + { nint = floor(u + 0.5);
1.677 + if (fabs(u - nint) <= 1e-5)
1.678 + u = nint;
1.679 + else
1.680 + u = floor(u);
1.681 + }
1.682 + /* check current column upper bound */
1.683 + if (q->ub != +DBL_MAX)
1.684 + { eps = (q->is_int ? 1e-3 : 1e-3 + 1e-6 * fabs(q->ub));
1.685 + if (u > q->ub - eps)
1.686 + { ret = 0; /* redundant */
1.687 + goto done;
1.688 + }
1.689 + }
1.690 + /* check current column lower bound */
1.691 + if (q->lb != -DBL_MAX)
1.692 + { eps = (q->is_int ? 1e-5 : 1e-5 + 1e-8 * fabs(q->lb));
1.693 + if (u < q->lb - eps)
1.694 + { ret = 4; /* infeasible */
1.695 + goto done;
1.696 + }
1.697 + /* if u'[q] is close to l[q], fix column at its lower bound */
1.698 + if (u < q->lb + 1e-3 * eps)
1.699 + { q->ub = q->lb;
1.700 + ret = 3; /* fixed */
1.701 + goto done;
1.702 + }
1.703 + }
1.704 + /* check if column upper bound changes significantly */
1.705 + if (q->ub == +DBL_MAX)
1.706 + ret = 2; /* significantly */
1.707 + else if (q->is_int && u < q->ub - 0.5)
1.708 + ret = 2; /* significantly */
1.709 + else if (u < q->ub - 0.30 * (1.0 + fabs(q->ub)))
1.710 + ret = 2; /* significantly */
1.711 + else
1.712 + ret = 1; /* not significantly */
1.713 + /* set new column upper bound */
1.714 + q->ub = u;
1.715 +done: return ret;
1.716 +}
1.717 +
1.718 +/***********************************************************************
1.719 +* NAME
1.720 +*
1.721 +* npp_ineq_singlet - process row singleton (inequality constraint)
1.722 +*
1.723 +* SYNOPSIS
1.724 +*
1.725 +* #include "glpnpp.h"
1.726 +* int npp_ineq_singlet(NPP *npp, NPPROW *p);
1.727 +*
1.728 +* DESCRIPTION
1.729 +*
1.730 +* The routine npp_ineq_singlet processes row p, which is inequality
1.731 +* constraint having the only non-zero coefficient:
1.732 +*
1.733 +* L[p] <= a[p,q] * x[q] <= U[p], (1)
1.734 +*
1.735 +* where L[p] < U[p], L[p] > -oo and/or U[p] < +oo.
1.736 +*
1.737 +* RETURNS
1.738 +*
1.739 +* 0 - current column bounds have not changed;
1.740 +*
1.741 +* 1 - current column bounds have changed, but not significantly;
1.742 +*
1.743 +* 2 - current column bounds have significantly changed;
1.744 +*
1.745 +* 3 - column has been fixed on its lower or upper bound;
1.746 +*
1.747 +* 4 - problem has no primal feasible solution.
1.748 +*
1.749 +* PROBLEM TRANSFORMATION
1.750 +*
1.751 +* Inequality constraint (1) defines implied bounds of column q:
1.752 +*
1.753 +* ( L[p] / a[p,q], if a[p,q] > 0
1.754 +* l'[q] = < (2)
1.755 +* ( U[p] / a[p,q], if a[p,q] < 0
1.756 +*
1.757 +* ( U[p] / a[p,q], if a[p,q] > 0
1.758 +* u'[q] = < (3)
1.759 +* ( L[p] / a[p,q], if a[p,q] < 0
1.760 +*
1.761 +* If these implied bounds do not violate current bounds of column q:
1.762 +*
1.763 +* l[q] <= x[q] <= u[q], (4)
1.764 +*
1.765 +* they can be used to strengthen the current column bounds:
1.766 +*
1.767 +* l[q] := max(l[q], l'[q]), (5)
1.768 +*
1.769 +* u[q] := min(u[q], u'[q]). (6)
1.770 +*
1.771 +* (See the routines npp_implied_lower and npp_implied_upper.)
1.772 +*
1.773 +* Once bounds of row p (1) have been carried over column q, the row
1.774 +* becomes redundant, so it can be replaced by equivalent free row and
1.775 +* removed from the problem.
1.776 +*
1.777 +* Note that the routine removes from the problem only row p. Column q,
1.778 +* even it has been fixed, is kept in the problem.
1.779 +*
1.780 +* RECOVERING BASIC SOLUTION
1.781 +*
1.782 +* Note that the row in the dual system corresponding to column q is
1.783 +* the following:
1.784 +*
1.785 +* sum a[i,q] pi[i] + lambda[q] = c[q] ==>
1.786 +* i
1.787 +* (7)
1.788 +* sum a[i,q] pi[i] + a[p,q] pi[p] + lambda[q] = c[q],
1.789 +* i!=p
1.790 +*
1.791 +* where pi[i] for all i != p are known in solution to the transformed
1.792 +* problem. Row p does not exist in the transformed problem, so it has
1.793 +* zero multiplier there. This allows computing multiplier for column q
1.794 +* in solution to the transformed problem:
1.795 +*
1.796 +* lambda~[q] = c[q] - sum a[i,q] pi[i]. (8)
1.797 +* i!=p
1.798 +*
1.799 +* Let in solution to the transformed problem column q be non-basic
1.800 +* with lower bound active (GLP_NL, lambda~[q] >= 0), and this lower
1.801 +* bound be implied one l'[q]. From the original problem's standpoint
1.802 +* this then means that actually the original column lower bound l[q]
1.803 +* is inactive, and active is that row bound L[p] or U[p] that defines
1.804 +* the implied bound l'[q] (2). In this case in solution to the
1.805 +* original problem column q is assigned status GLP_BS while row p is
1.806 +* assigned status GLP_NL (if a[p,q] > 0) or GLP_NU (if a[p,q] < 0).
1.807 +* Since now column q is basic, its multiplier lambda[q] is zero. This
1.808 +* allows using (7) and (8) to find multiplier for row p in solution to
1.809 +* the original problem:
1.810 +*
1.811 +* 1
1.812 +* pi[p] = ------ (c[q] - sum a[i,q] pi[i]) = lambda~[q] / a[p,q] (9)
1.813 +* a[p,q] i!=p
1.814 +*
1.815 +* Now let in solution to the transformed problem column q be non-basic
1.816 +* with upper bound active (GLP_NU, lambda~[q] <= 0), and this upper
1.817 +* bound be implied one u'[q]. As in the previous case this then means
1.818 +* that from the original problem's standpoint actually the original
1.819 +* column upper bound u[q] is inactive, and active is that row bound
1.820 +* L[p] or U[p] that defines the implied bound u'[q] (3). In this case
1.821 +* in solution to the original problem column q is assigned status
1.822 +* GLP_BS, row p is assigned status GLP_NU (if a[p,q] > 0) or GLP_NL
1.823 +* (if a[p,q] < 0), and its multiplier is computed with formula (9).
1.824 +*
1.825 +* Strengthening bounds of column q according to (5) and (6) may make
1.826 +* it fixed. Thus, if in solution to the transformed problem column q is
1.827 +* non-basic and fixed (GLP_NS), we can suppose that if lambda~[q] > 0,
1.828 +* column q has active lower bound (GLP_NL), and if lambda~[q] < 0,
1.829 +* column q has active upper bound (GLP_NU), reducing this case to two
1.830 +* previous ones. If, however, lambda~[q] is close to zero or
1.831 +* corresponding bound of row p does not exist (this may happen if
1.832 +* lambda~[q] has wrong sign due to round-off errors, in which case it
1.833 +* is expected to be close to zero, since solution is assumed to be dual
1.834 +* feasible), column q can be assigned status GLP_BS (basic), and row p
1.835 +* can be made active on its existing bound. In the latter case row
1.836 +* multiplier pi[p] computed with formula (9) will be also close to
1.837 +* zero, and dual feasibility will be kept.
1.838 +*
1.839 +* In all other cases, namely, if in solution to the transformed
1.840 +* problem column q is basic (GLP_BS), or non-basic with original lower
1.841 +* bound l[q] active (GLP_NL), or non-basic with original upper bound
1.842 +* u[q] active (GLP_NU), constraint (1) is inactive. So in solution to
1.843 +* the original problem status of column q remains unchanged, row p is
1.844 +* assigned status GLP_BS, and its multiplier pi[p] is assigned zero
1.845 +* value.
1.846 +*
1.847 +* RECOVERING INTERIOR-POINT SOLUTION
1.848 +*
1.849 +* First, value of multiplier for column q in solution to the original
1.850 +* problem is computed with formula (8). If lambda~[q] > 0 and column q
1.851 +* has implied lower bound, or if lambda~[q] < 0 and column q has
1.852 +* implied upper bound, this means that from the original problem's
1.853 +* standpoint actually row p has corresponding active bound, in which
1.854 +* case its multiplier pi[p] is computed with formula (9). In other
1.855 +* cases, when the sign of lambda~[q] corresponds to original bound of
1.856 +* column q, or when lambda~[q] =~ 0, value of row multiplier pi[p] is
1.857 +* assigned zero value.
1.858 +*
1.859 +* RECOVERING MIP SOLUTION
1.860 +*
1.861 +* None needed. */
1.862 +
1.863 +struct ineq_singlet
1.864 +{ /* row singleton (inequality constraint) */
1.865 + int p;
1.866 + /* row reference number */
1.867 + int q;
1.868 + /* column reference number */
1.869 + double apq;
1.870 + /* constraint coefficient a[p,q] */
1.871 + double c;
1.872 + /* objective coefficient at x[q] */
1.873 + double lb;
1.874 + /* row lower bound */
1.875 + double ub;
1.876 + /* row upper bound */
1.877 + char lb_changed;
1.878 + /* this flag is set if column lower bound was changed */
1.879 + char ub_changed;
1.880 + /* this flag is set if column upper bound was changed */
1.881 + NPPLFE *ptr;
1.882 + /* list of non-zero coefficients a[i,q], i != p */
1.883 +};
1.884 +
1.885 +static int rcv_ineq_singlet(NPP *npp, void *info);
1.886 +
1.887 +int npp_ineq_singlet(NPP *npp, NPPROW *p)
1.888 +{ /* process row singleton (inequality constraint) */
1.889 + struct ineq_singlet *info;
1.890 + NPPCOL *q;
1.891 + NPPAIJ *apq, *aij;
1.892 + NPPLFE *lfe;
1.893 + int lb_changed, ub_changed;
1.894 + double ll, uu;
1.895 + /* the row must be singleton inequality constraint */
1.896 + xassert(p->lb != -DBL_MAX || p->ub != +DBL_MAX);
1.897 + xassert(p->lb < p->ub);
1.898 + xassert(p->ptr != NULL && p->ptr->r_next == NULL);
1.899 + /* compute implied column bounds */
1.900 + apq = p->ptr;
1.901 + q = apq->col;
1.902 + xassert(q->lb < q->ub);
1.903 + if (apq->val > 0.0)
1.904 + { ll = (p->lb == -DBL_MAX ? -DBL_MAX : p->lb / apq->val);
1.905 + uu = (p->ub == +DBL_MAX ? +DBL_MAX : p->ub / apq->val);
1.906 + }
1.907 + else
1.908 + { ll = (p->ub == +DBL_MAX ? -DBL_MAX : p->ub / apq->val);
1.909 + uu = (p->lb == -DBL_MAX ? +DBL_MAX : p->lb / apq->val);
1.910 + }
1.911 + /* process implied column lower bound */
1.912 + if (ll == -DBL_MAX)
1.913 + lb_changed = 0;
1.914 + else
1.915 + { lb_changed = npp_implied_lower(npp, q, ll);
1.916 + xassert(0 <= lb_changed && lb_changed <= 4);
1.917 + if (lb_changed == 4) return 4; /* infeasible */
1.918 + }
1.919 + /* process implied column upper bound */
1.920 + if (uu == +DBL_MAX)
1.921 + ub_changed = 0;
1.922 + else if (lb_changed == 3)
1.923 + { /* column was fixed on its upper bound due to l'[q] = u[q] */
1.924 + /* note that L[p] < U[p], so l'[q] = u[q] < u'[q] */
1.925 + ub_changed = 0;
1.926 + }
1.927 + else
1.928 + { ub_changed = npp_implied_upper(npp, q, uu);
1.929 + xassert(0 <= ub_changed && ub_changed <= 4);
1.930 + if (ub_changed == 4) return 4; /* infeasible */
1.931 + }
1.932 + /* if neither lower nor upper column bound was changed, the row
1.933 + is originally redundant and can be replaced by free row */
1.934 + if (!lb_changed && !ub_changed)
1.935 + { p->lb = -DBL_MAX, p->ub = +DBL_MAX;
1.936 + npp_free_row(npp, p);
1.937 + return 0;
1.938 + }
1.939 + /* create transformation stack entry */
1.940 + info = npp_push_tse(npp,
1.941 + rcv_ineq_singlet, sizeof(struct ineq_singlet));
1.942 + info->p = p->i;
1.943 + info->q = q->j;
1.944 + info->apq = apq->val;
1.945 + info->c = q->coef;
1.946 + info->lb = p->lb;
1.947 + info->ub = p->ub;
1.948 + info->lb_changed = (char)lb_changed;
1.949 + info->ub_changed = (char)ub_changed;
1.950 + info->ptr = NULL;
1.951 + /* save column coefficients a[i,q], i != p (not needed for MIP
1.952 + solution) */
1.953 + if (npp->sol != GLP_MIP)
1.954 + { for (aij = q->ptr; aij != NULL; aij = aij->c_next)
1.955 + { if (aij == apq) continue; /* skip a[p,q] */
1.956 + lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE));
1.957 + lfe->ref = aij->row->i;
1.958 + lfe->val = aij->val;
1.959 + lfe->next = info->ptr;
1.960 + info->ptr = lfe;
1.961 + }
1.962 + }
1.963 + /* remove the row from the problem */
1.964 + npp_del_row(npp, p);
1.965 + return lb_changed >= ub_changed ? lb_changed : ub_changed;
1.966 +}
1.967 +
1.968 +static int rcv_ineq_singlet(NPP *npp, void *_info)
1.969 +{ /* recover row singleton (inequality constraint) */
1.970 + struct ineq_singlet *info = _info;
1.971 + NPPLFE *lfe;
1.972 + double lambda;
1.973 + if (npp->sol == GLP_MIP) goto done;
1.974 + /* compute lambda~[q] in solution to the transformed problem
1.975 + with formula (8) */
1.976 + lambda = info->c;
1.977 + for (lfe = info->ptr; lfe != NULL; lfe = lfe->next)
1.978 + lambda -= lfe->val * npp->r_pi[lfe->ref];
1.979 + if (npp->sol == GLP_SOL)
1.980 + { /* recover basic solution */
1.981 + if (npp->c_stat[info->q] == GLP_BS)
1.982 + { /* column q is basic, so row p is inactive */
1.983 + npp->r_stat[info->p] = GLP_BS;
1.984 + npp->r_pi[info->p] = 0.0;
1.985 + }
1.986 + else if (npp->c_stat[info->q] == GLP_NL)
1.987 +nl: { /* column q is non-basic with lower bound active */
1.988 + if (info->lb_changed)
1.989 + { /* it is implied bound, so actually row p is active
1.990 + while column q is basic */
1.991 + npp->r_stat[info->p] =
1.992 + (char)(info->apq > 0.0 ? GLP_NL : GLP_NU);
1.993 + npp->c_stat[info->q] = GLP_BS;
1.994 + npp->r_pi[info->p] = lambda / info->apq;
1.995 + }
1.996 + else
1.997 + { /* it is original bound, so row p is inactive */
1.998 + npp->r_stat[info->p] = GLP_BS;
1.999 + npp->r_pi[info->p] = 0.0;
1.1000 + }
1.1001 + }
1.1002 + else if (npp->c_stat[info->q] == GLP_NU)
1.1003 +nu: { /* column q is non-basic with upper bound active */
1.1004 + if (info->ub_changed)
1.1005 + { /* it is implied bound, so actually row p is active
1.1006 + while column q is basic */
1.1007 + npp->r_stat[info->p] =
1.1008 + (char)(info->apq > 0.0 ? GLP_NU : GLP_NL);
1.1009 + npp->c_stat[info->q] = GLP_BS;
1.1010 + npp->r_pi[info->p] = lambda / info->apq;
1.1011 + }
1.1012 + else
1.1013 + { /* it is original bound, so row p is inactive */
1.1014 + npp->r_stat[info->p] = GLP_BS;
1.1015 + npp->r_pi[info->p] = 0.0;
1.1016 + }
1.1017 + }
1.1018 + else if (npp->c_stat[info->q] == GLP_NS)
1.1019 + { /* column q is non-basic and fixed; note, however, that in
1.1020 + in the original problem it is non-fixed */
1.1021 + if (lambda > +1e-7)
1.1022 + { if (info->apq > 0.0 && info->lb != -DBL_MAX ||
1.1023 + info->apq < 0.0 && info->ub != +DBL_MAX ||
1.1024 + !info->lb_changed)
1.1025 + { /* either corresponding bound of row p exists or
1.1026 + column q remains non-basic with its original lower
1.1027 + bound active */
1.1028 + npp->c_stat[info->q] = GLP_NL;
1.1029 + goto nl;
1.1030 + }
1.1031 + }
1.1032 + if (lambda < -1e-7)
1.1033 + { if (info->apq > 0.0 && info->ub != +DBL_MAX ||
1.1034 + info->apq < 0.0 && info->lb != -DBL_MAX ||
1.1035 + !info->ub_changed)
1.1036 + { /* either corresponding bound of row p exists or
1.1037 + column q remains non-basic with its original upper
1.1038 + bound active */
1.1039 + npp->c_stat[info->q] = GLP_NU;
1.1040 + goto nu;
1.1041 + }
1.1042 + }
1.1043 + /* either lambda~[q] is close to zero, or corresponding
1.1044 + bound of row p does not exist, because lambda~[q] has
1.1045 + wrong sign due to round-off errors; in the latter case
1.1046 + lambda~[q] is also assumed to be close to zero; so, we
1.1047 + can make row p active on its existing bound and column q
1.1048 + basic; pi[p] will have wrong sign, but it also will be
1.1049 + close to zero (rarus casus of dual degeneracy) */
1.1050 + if (info->lb != -DBL_MAX && info->ub == +DBL_MAX)
1.1051 + { /* row lower bound exists, but upper bound doesn't */
1.1052 + npp->r_stat[info->p] = GLP_NL;
1.1053 + }
1.1054 + else if (info->lb == -DBL_MAX && info->ub != +DBL_MAX)
1.1055 + { /* row upper bound exists, but lower bound doesn't */
1.1056 + npp->r_stat[info->p] = GLP_NU;
1.1057 + }
1.1058 + else if (info->lb != -DBL_MAX && info->ub != +DBL_MAX)
1.1059 + { /* both row lower and upper bounds exist */
1.1060 + /* to choose proper active row bound we should not use
1.1061 + lambda~[q], because its value being close to zero is
1.1062 + unreliable; so we choose that bound which provides
1.1063 + primal feasibility for original constraint (1) */
1.1064 + if (info->apq * npp->c_value[info->q] <=
1.1065 + 0.5 * (info->lb + info->ub))
1.1066 + npp->r_stat[info->p] = GLP_NL;
1.1067 + else
1.1068 + npp->r_stat[info->p] = GLP_NU;
1.1069 + }
1.1070 + else
1.1071 + { npp_error();
1.1072 + return 1;
1.1073 + }
1.1074 + npp->c_stat[info->q] = GLP_BS;
1.1075 + npp->r_pi[info->p] = lambda / info->apq;
1.1076 + }
1.1077 + else
1.1078 + { npp_error();
1.1079 + return 1;
1.1080 + }
1.1081 + }
1.1082 + if (npp->sol == GLP_IPT)
1.1083 + { /* recover interior-point solution */
1.1084 + if (lambda > +DBL_EPSILON && info->lb_changed ||
1.1085 + lambda < -DBL_EPSILON && info->ub_changed)
1.1086 + { /* actually row p has corresponding active bound */
1.1087 + npp->r_pi[info->p] = lambda / info->apq;
1.1088 + }
1.1089 + else
1.1090 + { /* either bounds of column q are both inactive or its
1.1091 + original bound is active */
1.1092 + npp->r_pi[info->p] = 0.0;
1.1093 + }
1.1094 + }
1.1095 +done: return 0;
1.1096 +}
1.1097 +
1.1098 +/***********************************************************************
1.1099 +* NAME
1.1100 +*
1.1101 +* npp_implied_slack - process column singleton (implied slack variable)
1.1102 +*
1.1103 +* SYNOPSIS
1.1104 +*
1.1105 +* #include "glpnpp.h"
1.1106 +* void npp_implied_slack(NPP *npp, NPPCOL *q);
1.1107 +*
1.1108 +* DESCRIPTION
1.1109 +*
1.1110 +* The routine npp_implied_slack processes column q:
1.1111 +*
1.1112 +* l[q] <= x[q] <= u[q], (1)
1.1113 +*
1.1114 +* where l[q] < u[q], having the only non-zero coefficient in row p,
1.1115 +* which is equality constraint:
1.1116 +*
1.1117 +* sum a[p,j] x[j] + a[p,q] x[q] = b. (2)
1.1118 +* j!=q
1.1119 +*
1.1120 +* PROBLEM TRANSFORMATION
1.1121 +*
1.1122 +* (If x[q] is integral, this transformation must not be used.)
1.1123 +*
1.1124 +* The term a[p,q] x[q] in constraint (2) can be considered as a slack
1.1125 +* variable that allows to carry bounds of column q over row p and then
1.1126 +* remove column q from the problem.
1.1127 +*
1.1128 +* Constraint (2) can be written as follows:
1.1129 +*
1.1130 +* sum a[p,j] x[j] = b - a[p,q] x[q]. (3)
1.1131 +* j!=q
1.1132 +*
1.1133 +* According to (1) constraint (3) is equivalent to the following
1.1134 +* inequality constraint:
1.1135 +*
1.1136 +* L[p] <= sum a[p,j] x[j] <= U[p], (4)
1.1137 +* j!=q
1.1138 +*
1.1139 +* where
1.1140 +*
1.1141 +* ( b - a[p,q] u[q], if a[p,q] > 0
1.1142 +* L[p] = < (5)
1.1143 +* ( b - a[p,q] l[q], if a[p,q] < 0
1.1144 +*
1.1145 +* ( b - a[p,q] l[q], if a[p,q] > 0
1.1146 +* U[p] = < (6)
1.1147 +* ( b - a[p,q] u[q], if a[p,q] < 0
1.1148 +*
1.1149 +* From (2) it follows that:
1.1150 +*
1.1151 +* 1
1.1152 +* x[q] = ------ (b - sum a[p,j] x[j]). (7)
1.1153 +* a[p,q] j!=q
1.1154 +*
1.1155 +* In order to eliminate x[q] from the objective row we substitute it
1.1156 +* from (6) to that row:
1.1157 +*
1.1158 +* z = sum c[j] x[j] + c[q] x[q] + c[0] =
1.1159 +* j!=q
1.1160 +* 1
1.1161 +* = sum c[j] x[j] + c[q] [------ (b - sum a[p,j] x[j])] + c0 =
1.1162 +* j!=q a[p,q] j!=q
1.1163 +*
1.1164 +* = sum c~[j] x[j] + c~[0],
1.1165 +* j!=q
1.1166 +* a[p,j] b
1.1167 +* c~[j] = c[j] - c[q] ------, c~0 = c0 - c[q] ------ (8)
1.1168 +* a[p,q] a[p,q]
1.1169 +*
1.1170 +* are values of objective coefficients and constant term, resp., in
1.1171 +* the transformed problem.
1.1172 +*
1.1173 +* Note that column q is column singleton, so in the dual system of the
1.1174 +* original problem it corresponds to the following row singleton:
1.1175 +*
1.1176 +* a[p,q] pi[p] + lambda[q] = c[q]. (9)
1.1177 +*
1.1178 +* In the transformed problem row (9) would be the following:
1.1179 +*
1.1180 +* a[p,q] pi~[p] + lambda[q] = c~[q] = 0. (10)
1.1181 +*
1.1182 +* Subtracting (10) from (9) we have:
1.1183 +*
1.1184 +* a[p,q] (pi[p] - pi~[p]) = c[q]
1.1185 +*
1.1186 +* that gives the following formula to compute multiplier for row p in
1.1187 +* solution to the original problem using its value in solution to the
1.1188 +* transformed problem:
1.1189 +*
1.1190 +* pi[p] = pi~[p] + c[q] / a[p,q]. (11)
1.1191 +*
1.1192 +* RECOVERING BASIC SOLUTION
1.1193 +*
1.1194 +* Status of column q in solution to the original problem is defined
1.1195 +* by status of row p in solution to the transformed problem and the
1.1196 +* sign of coefficient a[p,q] in the original inequality constraint (2)
1.1197 +* as follows:
1.1198 +*
1.1199 +* +-----------------------+---------+--------------------+
1.1200 +* | Status of row p | Sign of | Status of column q |
1.1201 +* | (transformed problem) | a[p,q] | (original problem) |
1.1202 +* +-----------------------+---------+--------------------+
1.1203 +* | GLP_BS | + / - | GLP_BS |
1.1204 +* | GLP_NL | + | GLP_NU |
1.1205 +* | GLP_NL | - | GLP_NL |
1.1206 +* | GLP_NU | + | GLP_NL |
1.1207 +* | GLP_NU | - | GLP_NU |
1.1208 +* | GLP_NF | + / - | GLP_NF |
1.1209 +* +-----------------------+---------+--------------------+
1.1210 +*
1.1211 +* Value of column q is computed with formula (7). Since originally row
1.1212 +* p is equality constraint, its status is assigned GLP_NS, and value of
1.1213 +* its multiplier pi[p] is computed with formula (11).
1.1214 +*
1.1215 +* RECOVERING INTERIOR-POINT SOLUTION
1.1216 +*
1.1217 +* Value of column q is computed with formula (7). Row multiplier value
1.1218 +* pi[p] is computed with formula (11).
1.1219 +*
1.1220 +* RECOVERING MIP SOLUTION
1.1221 +*
1.1222 +* Value of column q is computed with formula (7). */
1.1223 +
1.1224 +struct implied_slack
1.1225 +{ /* column singleton (implied slack variable) */
1.1226 + int p;
1.1227 + /* row reference number */
1.1228 + int q;
1.1229 + /* column reference number */
1.1230 + double apq;
1.1231 + /* constraint coefficient a[p,q] */
1.1232 + double b;
1.1233 + /* right-hand side of original equality constraint */
1.1234 + double c;
1.1235 + /* original objective coefficient at x[q] */
1.1236 + NPPLFE *ptr;
1.1237 + /* list of non-zero coefficients a[p,j], j != q */
1.1238 +};
1.1239 +
1.1240 +static int rcv_implied_slack(NPP *npp, void *info);
1.1241 +
1.1242 +void npp_implied_slack(NPP *npp, NPPCOL *q)
1.1243 +{ /* process column singleton (implied slack variable) */
1.1244 + struct implied_slack *info;
1.1245 + NPPROW *p;
1.1246 + NPPAIJ *aij;
1.1247 + NPPLFE *lfe;
1.1248 + /* the column must be non-integral non-fixed singleton */
1.1249 + xassert(!q->is_int);
1.1250 + xassert(q->lb < q->ub);
1.1251 + xassert(q->ptr != NULL && q->ptr->c_next == NULL);
1.1252 + /* corresponding row must be equality constraint */
1.1253 + aij = q->ptr;
1.1254 + p = aij->row;
1.1255 + xassert(p->lb == p->ub);
1.1256 + /* create transformation stack entry */
1.1257 + info = npp_push_tse(npp,
1.1258 + rcv_implied_slack, sizeof(struct implied_slack));
1.1259 + info->p = p->i;
1.1260 + info->q = q->j;
1.1261 + info->apq = aij->val;
1.1262 + info->b = p->lb;
1.1263 + info->c = q->coef;
1.1264 + info->ptr = NULL;
1.1265 + /* save row coefficients a[p,j], j != q, and substitute x[q]
1.1266 + into the objective row */
1.1267 + for (aij = p->ptr; aij != NULL; aij = aij->r_next)
1.1268 + { if (aij->col == q) continue; /* skip a[p,q] */
1.1269 + lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE));
1.1270 + lfe->ref = aij->col->j;
1.1271 + lfe->val = aij->val;
1.1272 + lfe->next = info->ptr;
1.1273 + info->ptr = lfe;
1.1274 + aij->col->coef -= info->c * (aij->val / info->apq);
1.1275 + }
1.1276 + npp->c0 += info->c * (info->b / info->apq);
1.1277 + /* compute new row bounds */
1.1278 + if (info->apq > 0.0)
1.1279 + { p->lb = (q->ub == +DBL_MAX ?
1.1280 + -DBL_MAX : info->b - info->apq * q->ub);
1.1281 + p->ub = (q->lb == -DBL_MAX ?
1.1282 + +DBL_MAX : info->b - info->apq * q->lb);
1.1283 + }
1.1284 + else
1.1285 + { p->lb = (q->lb == -DBL_MAX ?
1.1286 + -DBL_MAX : info->b - info->apq * q->lb);
1.1287 + p->ub = (q->ub == +DBL_MAX ?
1.1288 + +DBL_MAX : info->b - info->apq * q->ub);
1.1289 + }
1.1290 + /* remove the column from the problem */
1.1291 + npp_del_col(npp, q);
1.1292 + return;
1.1293 +}
1.1294 +
1.1295 +static int rcv_implied_slack(NPP *npp, void *_info)
1.1296 +{ /* recover column singleton (implied slack variable) */
1.1297 + struct implied_slack *info = _info;
1.1298 + NPPLFE *lfe;
1.1299 + double temp;
1.1300 + if (npp->sol == GLP_SOL)
1.1301 + { /* assign statuses to row p and column q */
1.1302 + if (npp->r_stat[info->p] == GLP_BS ||
1.1303 + npp->r_stat[info->p] == GLP_NF)
1.1304 + npp->c_stat[info->q] = npp->r_stat[info->p];
1.1305 + else if (npp->r_stat[info->p] == GLP_NL)
1.1306 + npp->c_stat[info->q] =
1.1307 + (char)(info->apq > 0.0 ? GLP_NU : GLP_NL);
1.1308 + else if (npp->r_stat[info->p] == GLP_NU)
1.1309 + npp->c_stat[info->q] =
1.1310 + (char)(info->apq > 0.0 ? GLP_NL : GLP_NU);
1.1311 + else
1.1312 + { npp_error();
1.1313 + return 1;
1.1314 + }
1.1315 + npp->r_stat[info->p] = GLP_NS;
1.1316 + }
1.1317 + if (npp->sol != GLP_MIP)
1.1318 + { /* compute multiplier for row p */
1.1319 + npp->r_pi[info->p] += info->c / info->apq;
1.1320 + }
1.1321 + /* compute value of column q */
1.1322 + temp = info->b;
1.1323 + for (lfe = info->ptr; lfe != NULL; lfe = lfe->next)
1.1324 + temp -= lfe->val * npp->c_value[lfe->ref];
1.1325 + npp->c_value[info->q] = temp / info->apq;
1.1326 + return 0;
1.1327 +}
1.1328 +
1.1329 +/***********************************************************************
1.1330 +* NAME
1.1331 +*
1.1332 +* npp_implied_free - process column singleton (implied free variable)
1.1333 +*
1.1334 +* SYNOPSIS
1.1335 +*
1.1336 +* #include "glpnpp.h"
1.1337 +* int npp_implied_free(NPP *npp, NPPCOL *q);
1.1338 +*
1.1339 +* DESCRIPTION
1.1340 +*
1.1341 +* The routine npp_implied_free processes column q:
1.1342 +*
1.1343 +* l[q] <= x[q] <= u[q], (1)
1.1344 +*
1.1345 +* having non-zero coefficient in the only row p, which is inequality
1.1346 +* constraint:
1.1347 +*
1.1348 +* L[p] <= sum a[p,j] x[j] + a[p,q] x[q] <= U[p], (2)
1.1349 +* j!=q
1.1350 +*
1.1351 +* where l[q] < u[q], L[p] < U[p], L[p] > -oo and/or U[p] < +oo.
1.1352 +*
1.1353 +* RETURNS
1.1354 +*
1.1355 +* 0 - success;
1.1356 +*
1.1357 +* 1 - column lower and/or upper bound(s) can be active;
1.1358 +*
1.1359 +* 2 - problem has no dual feasible solution.
1.1360 +*
1.1361 +* PROBLEM TRANSFORMATION
1.1362 +*
1.1363 +* Constraint (2) can be written as follows:
1.1364 +*
1.1365 +* L[p] - sum a[p,j] x[j] <= a[p,q] x[q] <= U[p] - sum a[p,j] x[j],
1.1366 +* j!=q j!=q
1.1367 +*
1.1368 +* from which it follows that:
1.1369 +*
1.1370 +* alfa <= a[p,q] x[q] <= beta, (3)
1.1371 +*
1.1372 +* where
1.1373 +*
1.1374 +* alfa = inf(L[p] - sum a[p,j] x[j]) =
1.1375 +* j!=q
1.1376 +*
1.1377 +* = L[p] - sup sum a[p,j] x[j] = (4)
1.1378 +* j!=q
1.1379 +*
1.1380 +* = L[p] - sum a[p,j] u[j] - sum a[p,j] l[j],
1.1381 +* j in Jp j in Jn
1.1382 +*
1.1383 +* beta = sup(L[p] - sum a[p,j] x[j]) =
1.1384 +* j!=q
1.1385 +*
1.1386 +* = L[p] - inf sum a[p,j] x[j] = (5)
1.1387 +* j!=q
1.1388 +*
1.1389 +* = L[p] - sum a[p,j] l[j] - sum a[p,j] u[j],
1.1390 +* j in Jp j in Jn
1.1391 +*
1.1392 +* Jp = {j != q: a[p,j] > 0}, Jn = {j != q: a[p,j] < 0}. (6)
1.1393 +*
1.1394 +* Inequality (3) defines implied bounds of variable x[q]:
1.1395 +*
1.1396 +* l'[q] <= x[q] <= u'[q], (7)
1.1397 +*
1.1398 +* where
1.1399 +*
1.1400 +* ( alfa / a[p,q], if a[p,q] > 0
1.1401 +* l'[q] = < (8a)
1.1402 +* ( beta / a[p,q], if a[p,q] < 0
1.1403 +*
1.1404 +* ( beta / a[p,q], if a[p,q] > 0
1.1405 +* u'[q] = < (8b)
1.1406 +* ( alfa / a[p,q], if a[p,q] < 0
1.1407 +*
1.1408 +* Thus, if l'[q] > l[q] - eps and u'[q] < u[q] + eps, where eps is
1.1409 +* an absolute tolerance for column value, column bounds (1) cannot be
1.1410 +* active, in which case column q can be replaced by equivalent free
1.1411 +* (unbounded) column.
1.1412 +*
1.1413 +* Note that column q is column singleton, so in the dual system of the
1.1414 +* original problem it corresponds to the following row singleton:
1.1415 +*
1.1416 +* a[p,q] pi[p] + lambda[q] = c[q], (9)
1.1417 +*
1.1418 +* from which it follows that:
1.1419 +*
1.1420 +* pi[p] = (c[q] - lambda[q]) / a[p,q]. (10)
1.1421 +*
1.1422 +* Let x[q] be implied free (unbounded) variable. Then column q can be
1.1423 +* only basic, so its multiplier lambda[q] is equal to zero, and from
1.1424 +* (10) we have:
1.1425 +*
1.1426 +* pi[p] = c[q] / a[p,q]. (11)
1.1427 +*
1.1428 +* There are possible three cases:
1.1429 +*
1.1430 +* 1) pi[p] < -eps, where eps is an absolute tolerance for row
1.1431 +* multiplier. In this case, to provide dual feasibility of the
1.1432 +* original problem, row p must be active on its lower bound, and
1.1433 +* if its lower bound does not exist (L[p] = -oo), the problem has
1.1434 +* no dual feasible solution;
1.1435 +*
1.1436 +* 2) pi[p] > +eps. In this case row p must be active on its upper
1.1437 +* bound, and if its upper bound does not exist (U[p] = +oo), the
1.1438 +* problem has no dual feasible solution;
1.1439 +*
1.1440 +* 3) -eps <= pi[p] <= +eps. In this case any (either lower or upper)
1.1441 +* bound of row p can be active, because this does not affect dual
1.1442 +* feasibility.
1.1443 +*
1.1444 +* Thus, in all three cases original inequality constraint (2) can be
1.1445 +* replaced by equality constraint, where the right-hand side is either
1.1446 +* lower or upper bound of row p, and bounds of column q can be removed
1.1447 +* that makes it free (unbounded). (May note that this transformation
1.1448 +* can be followed by transformation "Column singleton (implied slack
1.1449 +* variable)" performed by the routine npp_implied_slack.)
1.1450 +*
1.1451 +* RECOVERING BASIC SOLUTION
1.1452 +*
1.1453 +* Status of row p in solution to the original problem is determined
1.1454 +* by its status in solution to the transformed problem and its bound,
1.1455 +* which was choosen to be active:
1.1456 +*
1.1457 +* +-----------------------+--------+--------------------+
1.1458 +* | Status of row p | Active | Status of row p |
1.1459 +* | (transformed problem) | bound | (original problem) |
1.1460 +* +-----------------------+--------+--------------------+
1.1461 +* | GLP_BS | L[p] | GLP_BS |
1.1462 +* | GLP_BS | U[p] | GLP_BS |
1.1463 +* | GLP_NS | L[p] | GLP_NL |
1.1464 +* | GLP_NS | U[p] | GLP_NU |
1.1465 +* +-----------------------+--------+--------------------+
1.1466 +*
1.1467 +* Value of row multiplier pi[p] (as well as value of column q) in
1.1468 +* solution to the original problem is the same as in solution to the
1.1469 +* transformed problem.
1.1470 +*
1.1471 +* RECOVERING INTERIOR-POINT SOLUTION
1.1472 +*
1.1473 +* Value of row multiplier pi[p] in solution to the original problem is
1.1474 +* the same as in solution to the transformed problem.
1.1475 +*
1.1476 +* RECOVERING MIP SOLUTION
1.1477 +*
1.1478 +* None needed. */
1.1479 +
1.1480 +struct implied_free
1.1481 +{ /* column singleton (implied free variable) */
1.1482 + int p;
1.1483 + /* row reference number */
1.1484 + char stat;
1.1485 + /* row status:
1.1486 + GLP_NL - active constraint on lower bound
1.1487 + GLP_NU - active constraint on upper bound */
1.1488 +};
1.1489 +
1.1490 +static int rcv_implied_free(NPP *npp, void *info);
1.1491 +
1.1492 +int npp_implied_free(NPP *npp, NPPCOL *q)
1.1493 +{ /* process column singleton (implied free variable) */
1.1494 + struct implied_free *info;
1.1495 + NPPROW *p;
1.1496 + NPPAIJ *apq, *aij;
1.1497 + double alfa, beta, l, u, pi, eps;
1.1498 + /* the column must be non-fixed singleton */
1.1499 + xassert(q->lb < q->ub);
1.1500 + xassert(q->ptr != NULL && q->ptr->c_next == NULL);
1.1501 + /* corresponding row must be inequality constraint */
1.1502 + apq = q->ptr;
1.1503 + p = apq->row;
1.1504 + xassert(p->lb != -DBL_MAX || p->ub != +DBL_MAX);
1.1505 + xassert(p->lb < p->ub);
1.1506 + /* compute alfa */
1.1507 + alfa = p->lb;
1.1508 + if (alfa != -DBL_MAX)
1.1509 + { for (aij = p->ptr; aij != NULL; aij = aij->r_next)
1.1510 + { if (aij == apq) continue; /* skip a[p,q] */
1.1511 + if (aij->val > 0.0)
1.1512 + { if (aij->col->ub == +DBL_MAX)
1.1513 + { alfa = -DBL_MAX;
1.1514 + break;
1.1515 + }
1.1516 + alfa -= aij->val * aij->col->ub;
1.1517 + }
1.1518 + else /* < 0.0 */
1.1519 + { if (aij->col->lb == -DBL_MAX)
1.1520 + { alfa = -DBL_MAX;
1.1521 + break;
1.1522 + }
1.1523 + alfa -= aij->val * aij->col->lb;
1.1524 + }
1.1525 + }
1.1526 + }
1.1527 + /* compute beta */
1.1528 + beta = p->ub;
1.1529 + if (beta != +DBL_MAX)
1.1530 + { for (aij = p->ptr; aij != NULL; aij = aij->r_next)
1.1531 + { if (aij == apq) continue; /* skip a[p,q] */
1.1532 + if (aij->val > 0.0)
1.1533 + { if (aij->col->lb == -DBL_MAX)
1.1534 + { beta = +DBL_MAX;
1.1535 + break;
1.1536 + }
1.1537 + beta -= aij->val * aij->col->lb;
1.1538 + }
1.1539 + else /* < 0.0 */
1.1540 + { if (aij->col->ub == +DBL_MAX)
1.1541 + { beta = +DBL_MAX;
1.1542 + break;
1.1543 + }
1.1544 + beta -= aij->val * aij->col->ub;
1.1545 + }
1.1546 + }
1.1547 + }
1.1548 + /* compute implied column lower bound l'[q] */
1.1549 + if (apq->val > 0.0)
1.1550 + l = (alfa == -DBL_MAX ? -DBL_MAX : alfa / apq->val);
1.1551 + else /* < 0.0 */
1.1552 + l = (beta == +DBL_MAX ? -DBL_MAX : beta / apq->val);
1.1553 + /* compute implied column upper bound u'[q] */
1.1554 + if (apq->val > 0.0)
1.1555 + u = (beta == +DBL_MAX ? +DBL_MAX : beta / apq->val);
1.1556 + else
1.1557 + u = (alfa == -DBL_MAX ? +DBL_MAX : alfa / apq->val);
1.1558 + /* check if column lower bound l[q] can be active */
1.1559 + if (q->lb != -DBL_MAX)
1.1560 + { eps = 1e-9 + 1e-12 * fabs(q->lb);
1.1561 + if (l < q->lb - eps) return 1; /* yes, it can */
1.1562 + }
1.1563 + /* check if column upper bound u[q] can be active */
1.1564 + if (q->ub != +DBL_MAX)
1.1565 + { eps = 1e-9 + 1e-12 * fabs(q->ub);
1.1566 + if (u > q->ub + eps) return 1; /* yes, it can */
1.1567 + }
1.1568 + /* okay; make column q free (unbounded) */
1.1569 + q->lb = -DBL_MAX, q->ub = +DBL_MAX;
1.1570 + /* create transformation stack entry */
1.1571 + info = npp_push_tse(npp,
1.1572 + rcv_implied_free, sizeof(struct implied_free));
1.1573 + info->p = p->i;
1.1574 + info->stat = -1;
1.1575 + /* compute row multiplier pi[p] */
1.1576 + pi = q->coef / apq->val;
1.1577 + /* check dual feasibility for row p */
1.1578 + if (pi > +DBL_EPSILON)
1.1579 + { /* lower bound L[p] must be active */
1.1580 + if (p->lb != -DBL_MAX)
1.1581 +nl: { info->stat = GLP_NL;
1.1582 + p->ub = p->lb;
1.1583 + }
1.1584 + else
1.1585 + { if (pi > +1e-5) return 2; /* dual infeasibility */
1.1586 + /* take a chance on U[p] */
1.1587 + xassert(p->ub != +DBL_MAX);
1.1588 + goto nu;
1.1589 + }
1.1590 + }
1.1591 + else if (pi < -DBL_EPSILON)
1.1592 + { /* upper bound U[p] must be active */
1.1593 + if (p->ub != +DBL_MAX)
1.1594 +nu: { info->stat = GLP_NU;
1.1595 + p->lb = p->ub;
1.1596 + }
1.1597 + else
1.1598 + { if (pi < -1e-5) return 2; /* dual infeasibility */
1.1599 + /* take a chance on L[p] */
1.1600 + xassert(p->lb != -DBL_MAX);
1.1601 + goto nl;
1.1602 + }
1.1603 + }
1.1604 + else
1.1605 + { /* any bound (either L[p] or U[p]) can be made active */
1.1606 + if (p->ub == +DBL_MAX)
1.1607 + { xassert(p->lb != -DBL_MAX);
1.1608 + goto nl;
1.1609 + }
1.1610 + if (p->lb == -DBL_MAX)
1.1611 + { xassert(p->ub != +DBL_MAX);
1.1612 + goto nu;
1.1613 + }
1.1614 + if (fabs(p->lb) <= fabs(p->ub)) goto nl; else goto nu;
1.1615 + }
1.1616 + return 0;
1.1617 +}
1.1618 +
1.1619 +static int rcv_implied_free(NPP *npp, void *_info)
1.1620 +{ /* recover column singleton (implied free variable) */
1.1621 + struct implied_free *info = _info;
1.1622 + if (npp->sol == GLP_SOL)
1.1623 + { if (npp->r_stat[info->p] == GLP_BS)
1.1624 + npp->r_stat[info->p] = GLP_BS;
1.1625 + else if (npp->r_stat[info->p] == GLP_NS)
1.1626 + { xassert(info->stat == GLP_NL || info->stat == GLP_NU);
1.1627 + npp->r_stat[info->p] = info->stat;
1.1628 + }
1.1629 + else
1.1630 + { npp_error();
1.1631 + return 1;
1.1632 + }
1.1633 + }
1.1634 + return 0;
1.1635 +}
1.1636 +
1.1637 +/***********************************************************************
1.1638 +* NAME
1.1639 +*
1.1640 +* npp_eq_doublet - process row doubleton (equality constraint)
1.1641 +*
1.1642 +* SYNOPSIS
1.1643 +*
1.1644 +* #include "glpnpp.h"
1.1645 +* NPPCOL *npp_eq_doublet(NPP *npp, NPPROW *p);
1.1646 +*
1.1647 +* DESCRIPTION
1.1648 +*
1.1649 +* The routine npp_eq_doublet processes row p, which is equality
1.1650 +* constraint having exactly two non-zero coefficients:
1.1651 +*
1.1652 +* a[p,q] x[q] + a[p,r] x[r] = b. (1)
1.1653 +*
1.1654 +* As the result of processing one of columns q or r is eliminated from
1.1655 +* all other rows and, thus, becomes column singleton of type "implied
1.1656 +* slack variable". Row p is not changed and along with column q and r
1.1657 +* remains in the problem.
1.1658 +*
1.1659 +* RETURNS
1.1660 +*
1.1661 +* The routine npp_eq_doublet returns pointer to the descriptor of that
1.1662 +* column q or r which has been eliminated. If, due to some reason, the
1.1663 +* elimination was not performed, the routine returns NULL.
1.1664 +*
1.1665 +* PROBLEM TRANSFORMATION
1.1666 +*
1.1667 +* First, we decide which column q or r will be eliminated. Let it be
1.1668 +* column q. Consider i-th constraint row, where column q has non-zero
1.1669 +* coefficient a[i,q] != 0:
1.1670 +*
1.1671 +* L[i] <= sum a[i,j] x[j] <= U[i]. (2)
1.1672 +* j
1.1673 +*
1.1674 +* In order to eliminate column q from row (2) we subtract from it row
1.1675 +* (1) multiplied by gamma[i] = a[i,q] / a[p,q], i.e. we replace in the
1.1676 +* transformed problem row (2) by its linear combination with row (1).
1.1677 +* This transformation changes only coefficients in columns q and r,
1.1678 +* and bounds of row i as follows:
1.1679 +*
1.1680 +* a~[i,q] = a[i,q] - gamma[i] a[p,q] = 0, (3)
1.1681 +*
1.1682 +* a~[i,r] = a[i,r] - gamma[i] a[p,r], (4)
1.1683 +*
1.1684 +* L~[i] = L[i] - gamma[i] b, (5)
1.1685 +*
1.1686 +* U~[i] = U[i] - gamma[i] b. (6)
1.1687 +*
1.1688 +* RECOVERING BASIC SOLUTION
1.1689 +*
1.1690 +* The transformation of the primal system of the original problem:
1.1691 +*
1.1692 +* L <= A x <= U (7)
1.1693 +*
1.1694 +* is equivalent to multiplying from the left a transformation matrix F
1.1695 +* by components of this primal system, which in the transformed problem
1.1696 +* becomes the following:
1.1697 +*
1.1698 +* F L <= F A x <= F U ==> L~ <= A~x <= U~. (8)
1.1699 +*
1.1700 +* The matrix F has the following structure:
1.1701 +*
1.1702 +* ( 1 -gamma[1] )
1.1703 +* ( )
1.1704 +* ( 1 -gamma[2] )
1.1705 +* ( )
1.1706 +* ( ... ... )
1.1707 +* ( )
1.1708 +* F = ( 1 -gamma[p-1] ) (9)
1.1709 +* ( )
1.1710 +* ( 1 )
1.1711 +* ( )
1.1712 +* ( -gamma[p+1] 1 )
1.1713 +* ( )
1.1714 +* ( ... ... )
1.1715 +*
1.1716 +* where its column containing elements -gamma[i] corresponds to row p
1.1717 +* of the primal system.
1.1718 +*
1.1719 +* From (8) it follows that the dual system of the original problem:
1.1720 +*
1.1721 +* A'pi + lambda = c, (10)
1.1722 +*
1.1723 +* in the transformed problem becomes the following:
1.1724 +*
1.1725 +* A'F'inv(F')pi + lambda = c ==> (A~)'pi~ + lambda = c, (11)
1.1726 +*
1.1727 +* where:
1.1728 +*
1.1729 +* pi~ = inv(F')pi (12)
1.1730 +*
1.1731 +* is the vector of row multipliers in the transformed problem. Thus:
1.1732 +*
1.1733 +* pi = F'pi~. (13)
1.1734 +*
1.1735 +* Therefore, as it follows from (13), value of multiplier for row p in
1.1736 +* solution to the original problem can be computed as follows:
1.1737 +*
1.1738 +* pi[p] = pi~[p] - sum gamma[i] pi~[i], (14)
1.1739 +* i
1.1740 +*
1.1741 +* where pi~[i] = pi[i] is multiplier for row i (i != p).
1.1742 +*
1.1743 +* Note that the statuses of all rows and columns are not changed.
1.1744 +*
1.1745 +* RECOVERING INTERIOR-POINT SOLUTION
1.1746 +*
1.1747 +* Multiplier for row p in solution to the original problem is computed
1.1748 +* with formula (14).
1.1749 +*
1.1750 +* RECOVERING MIP SOLUTION
1.1751 +*
1.1752 +* None needed. */
1.1753 +
1.1754 +struct eq_doublet
1.1755 +{ /* row doubleton (equality constraint) */
1.1756 + int p;
1.1757 + /* row reference number */
1.1758 + double apq;
1.1759 + /* constraint coefficient a[p,q] */
1.1760 + NPPLFE *ptr;
1.1761 + /* list of non-zero coefficients a[i,q], i != p */
1.1762 +};
1.1763 +
1.1764 +static int rcv_eq_doublet(NPP *npp, void *info);
1.1765 +
1.1766 +NPPCOL *npp_eq_doublet(NPP *npp, NPPROW *p)
1.1767 +{ /* process row doubleton (equality constraint) */
1.1768 + struct eq_doublet *info;
1.1769 + NPPROW *i;
1.1770 + NPPCOL *q, *r;
1.1771 + NPPAIJ *apq, *apr, *aiq, *air, *next;
1.1772 + NPPLFE *lfe;
1.1773 + double gamma;
1.1774 + /* the row must be doubleton equality constraint */
1.1775 + xassert(p->lb == p->ub);
1.1776 + xassert(p->ptr != NULL && p->ptr->r_next != NULL &&
1.1777 + p->ptr->r_next->r_next == NULL);
1.1778 + /* choose column to be eliminated */
1.1779 + { NPPAIJ *a1, *a2;
1.1780 + a1 = p->ptr, a2 = a1->r_next;
1.1781 + if (fabs(a2->val) < 0.001 * fabs(a1->val))
1.1782 + { /* only first column can be eliminated, because second one
1.1783 + has too small constraint coefficient */
1.1784 + apq = a1, apr = a2;
1.1785 + }
1.1786 + else if (fabs(a1->val) < 0.001 * fabs(a2->val))
1.1787 + { /* only second column can be eliminated, because first one
1.1788 + has too small constraint coefficient */
1.1789 + apq = a2, apr = a1;
1.1790 + }
1.1791 + else
1.1792 + { /* both columns are appropriate; choose that one which is
1.1793 + shorter to minimize fill-in */
1.1794 + if (npp_col_nnz(npp, a1->col) <= npp_col_nnz(npp, a2->col))
1.1795 + { /* first column is shorter */
1.1796 + apq = a1, apr = a2;
1.1797 + }
1.1798 + else
1.1799 + { /* second column is shorter */
1.1800 + apq = a2, apr = a1;
1.1801 + }
1.1802 + }
1.1803 + }
1.1804 + /* now columns q and r have been chosen */
1.1805 + q = apq->col, r = apr->col;
1.1806 + /* create transformation stack entry */
1.1807 + info = npp_push_tse(npp,
1.1808 + rcv_eq_doublet, sizeof(struct eq_doublet));
1.1809 + info->p = p->i;
1.1810 + info->apq = apq->val;
1.1811 + info->ptr = NULL;
1.1812 + /* transform each row i (i != p), where a[i,q] != 0, to eliminate
1.1813 + column q */
1.1814 + for (aiq = q->ptr; aiq != NULL; aiq = next)
1.1815 + { next = aiq->c_next;
1.1816 + if (aiq == apq) continue; /* skip row p */
1.1817 + i = aiq->row; /* row i to be transformed */
1.1818 + /* save constraint coefficient a[i,q] */
1.1819 + if (npp->sol != GLP_MIP)
1.1820 + { lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE));
1.1821 + lfe->ref = i->i;
1.1822 + lfe->val = aiq->val;
1.1823 + lfe->next = info->ptr;
1.1824 + info->ptr = lfe;
1.1825 + }
1.1826 + /* find coefficient a[i,r] in row i */
1.1827 + for (air = i->ptr; air != NULL; air = air->r_next)
1.1828 + if (air->col == r) break;
1.1829 + /* if a[i,r] does not exist, create a[i,r] = 0 */
1.1830 + if (air == NULL)
1.1831 + air = npp_add_aij(npp, i, r, 0.0);
1.1832 + /* compute gamma[i] = a[i,q] / a[p,q] */
1.1833 + gamma = aiq->val / apq->val;
1.1834 + /* (row i) := (row i) - gamma[i] * (row p); see (3)-(6) */
1.1835 + /* new a[i,q] is exact zero due to elimnation; remove it from
1.1836 + row i */
1.1837 + npp_del_aij(npp, aiq);
1.1838 + /* compute new a[i,r] */
1.1839 + air->val -= gamma * apr->val;
1.1840 + /* if new a[i,r] is close to zero due to numeric cancelation,
1.1841 + remove it from row i */
1.1842 + if (fabs(air->val) <= 1e-10)
1.1843 + npp_del_aij(npp, air);
1.1844 + /* compute new lower and upper bounds of row i */
1.1845 + if (i->lb == i->ub)
1.1846 + i->lb = i->ub = (i->lb - gamma * p->lb);
1.1847 + else
1.1848 + { if (i->lb != -DBL_MAX)
1.1849 + i->lb -= gamma * p->lb;
1.1850 + if (i->ub != +DBL_MAX)
1.1851 + i->ub -= gamma * p->lb;
1.1852 + }
1.1853 + }
1.1854 + return q;
1.1855 +}
1.1856 +
1.1857 +static int rcv_eq_doublet(NPP *npp, void *_info)
1.1858 +{ /* recover row doubleton (equality constraint) */
1.1859 + struct eq_doublet *info = _info;
1.1860 + NPPLFE *lfe;
1.1861 + double gamma, temp;
1.1862 + /* we assume that processing row p is followed by processing
1.1863 + column q as singleton of type "implied slack variable", in
1.1864 + which case row p must always be active equality constraint */
1.1865 + if (npp->sol == GLP_SOL)
1.1866 + { if (npp->r_stat[info->p] != GLP_NS)
1.1867 + { npp_error();
1.1868 + return 1;
1.1869 + }
1.1870 + }
1.1871 + if (npp->sol != GLP_MIP)
1.1872 + { /* compute value of multiplier for row p; see (14) */
1.1873 + temp = npp->r_pi[info->p];
1.1874 + for (lfe = info->ptr; lfe != NULL; lfe = lfe->next)
1.1875 + { gamma = lfe->val / info->apq; /* a[i,q] / a[p,q] */
1.1876 + temp -= gamma * npp->r_pi[lfe->ref];
1.1877 + }
1.1878 + npp->r_pi[info->p] = temp;
1.1879 + }
1.1880 + return 0;
1.1881 +}
1.1882 +
1.1883 +/***********************************************************************
1.1884 +* NAME
1.1885 +*
1.1886 +* npp_forcing_row - process forcing row
1.1887 +*
1.1888 +* SYNOPSIS
1.1889 +*
1.1890 +* #include "glpnpp.h"
1.1891 +* int npp_forcing_row(NPP *npp, NPPROW *p, int at);
1.1892 +*
1.1893 +* DESCRIPTION
1.1894 +*
1.1895 +* The routine npp_forcing row processes row p of general format:
1.1896 +*
1.1897 +* L[p] <= sum a[p,j] x[j] <= U[p], (1)
1.1898 +* j
1.1899 +*
1.1900 +* l[j] <= x[j] <= u[j], (2)
1.1901 +*
1.1902 +* where L[p] <= U[p] and l[j] < u[j] for all a[p,j] != 0. It is also
1.1903 +* assumed that:
1.1904 +*
1.1905 +* 1) if at = 0 then |L[p] - U'[p]| <= eps, where U'[p] is implied
1.1906 +* row upper bound (see below), eps is an absolute tolerance for row
1.1907 +* value;
1.1908 +*
1.1909 +* 2) if at = 1 then |U[p] - L'[p]| <= eps, where L'[p] is implied
1.1910 +* row lower bound (see below).
1.1911 +*
1.1912 +* RETURNS
1.1913 +*
1.1914 +* 0 - success;
1.1915 +*
1.1916 +* 1 - cannot fix columns due to too small constraint coefficients.
1.1917 +*
1.1918 +* PROBLEM TRANSFORMATION
1.1919 +*
1.1920 +* Implied lower and upper bounds of row (1) are determined by bounds
1.1921 +* of corresponding columns (variables) as follows:
1.1922 +*
1.1923 +* L'[p] = inf sum a[p,j] x[j] =
1.1924 +* j
1.1925 +* (3)
1.1926 +* = sum a[p,j] l[j] + sum a[p,j] u[j],
1.1927 +* j in Jp j in Jn
1.1928 +*
1.1929 +* U'[p] = sup sum a[p,j] x[j] =
1.1930 +* (4)
1.1931 +* = sum a[p,j] u[j] + sum a[p,j] l[j],
1.1932 +* j in Jp j in Jn
1.1933 +*
1.1934 +* Jp = {j: a[p,j] > 0}, Jn = {j: a[p,j] < 0}. (5)
1.1935 +*
1.1936 +* If L[p] =~ U'[p] (at = 0), solution can be primal feasible only when
1.1937 +* all variables take their boundary values as defined by (4):
1.1938 +*
1.1939 +* ( u[j], if j in Jp
1.1940 +* x[j] = < (6)
1.1941 +* ( l[j], if j in Jn
1.1942 +*
1.1943 +* Similarly, if U[p] =~ L'[p] (at = 1), solution can be primal feasible
1.1944 +* only when all variables take their boundary values as defined by (3):
1.1945 +*
1.1946 +* ( l[j], if j in Jp
1.1947 +* x[j] = < (7)
1.1948 +* ( u[j], if j in Jn
1.1949 +*
1.1950 +* Condition (6) or (7) allows fixing all columns (variables x[j])
1.1951 +* in row (1) on their bounds and then removing them from the problem
1.1952 +* (see the routine npp_fixed_col). Due to this row p becomes redundant,
1.1953 +* so it can be replaced by equivalent free (unbounded) row and also
1.1954 +* removed from the problem (see the routine npp_free_row).
1.1955 +*
1.1956 +* 1. To apply this transformation row (1) should not have coefficients
1.1957 +* whose magnitude is too small, i.e. all a[p,j] should satisfy to
1.1958 +* the following condition:
1.1959 +*
1.1960 +* |a[p,j]| >= eps * max(1, |a[p,k]|), (8)
1.1961 +* k
1.1962 +* where eps is a relative tolerance for constraint coefficients.
1.1963 +* Otherwise, fixing columns may be numerically unreliable and may
1.1964 +* lead to wrong solution.
1.1965 +*
1.1966 +* 2. The routine fixes columns and remove bounds of row p, however,
1.1967 +* it does not remove the row and columns from the problem.
1.1968 +*
1.1969 +* RECOVERING BASIC SOLUTION
1.1970 +*
1.1971 +* In the transformed problem row p being inactive constraint is
1.1972 +* assigned status GLP_BS (as the result of transformation of free
1.1973 +* row), and all columns in this row are assigned status GLP_NS (as the
1.1974 +* result of transformation of fixed columns).
1.1975 +*
1.1976 +* Note that in the dual system of the transformed (as well as original)
1.1977 +* problem every column j in row p corresponds to the following row:
1.1978 +*
1.1979 +* sum a[i,j] pi[i] + a[p,j] pi[p] + lambda[j] = c[j], (9)
1.1980 +* i!=p
1.1981 +*
1.1982 +* from which it follows that:
1.1983 +*
1.1984 +* lambda[j] = c[j] - sum a[i,j] pi[i] - a[p,j] pi[p]. (10)
1.1985 +* i!=p
1.1986 +*
1.1987 +* In the transformed problem values of all multipliers pi[i] are known
1.1988 +* (including pi[i], whose value is zero, since row p is inactive).
1.1989 +* Thus, using formula (10) it is possible to compute values of
1.1990 +* multipliers lambda[j] for all columns in row p.
1.1991 +*
1.1992 +* Note also that in the original problem all columns in row p are
1.1993 +* bounded, not fixed. So status GLP_NS assigned to every such column
1.1994 +* must be changed to GLP_NL or GLP_NU depending on which bound the
1.1995 +* corresponding column has been fixed. This status change may lead to
1.1996 +* dual feasibility violation for solution of the original problem,
1.1997 +* because now column multipliers must satisfy to the following
1.1998 +* condition:
1.1999 +*
1.2000 +* ( >= 0, if status of column j is GLP_NL,
1.2001 +* lambda[j] < (11)
1.2002 +* ( <= 0, if status of column j is GLP_NU.
1.2003 +*
1.2004 +* If this condition holds, solution to the original problem is the
1.2005 +* same as to the transformed problem. Otherwise, we have to perform
1.2006 +* one degenerate pivoting step of the primal simplex method to obtain
1.2007 +* dual feasible (hence, optimal) solution to the original problem as
1.2008 +* follows. If, on problem transformation, row p was made active on its
1.2009 +* lower bound (case at = 0), we change its status to GLP_NL (or GLP_NS)
1.2010 +* and start increasing its multiplier pi[p]. Otherwise, if row p was
1.2011 +* made active on its upper bound (case at = 1), we change its status
1.2012 +* to GLP_NU (or GLP_NS) and start decreasing pi[p]. From (10) it
1.2013 +* follows that:
1.2014 +*
1.2015 +* delta lambda[j] = - a[p,j] * delta pi[p] = - a[p,j] pi[p]. (12)
1.2016 +*
1.2017 +* Simple analysis of formulae (3)-(5) shows that changing pi[p] in the
1.2018 +* specified direction causes increasing lambda[j] for every column j
1.2019 +* assigned status GLP_NL (delta lambda[j] > 0) and decreasing lambda[j]
1.2020 +* for every column j assigned status GLP_NU (delta lambda[j] < 0). It
1.2021 +* is understood that once the last lambda[q], which violates condition
1.2022 +* (11), has reached zero, multipliers lambda[j] for all columns get
1.2023 +* valid signs. Such column q can be determined as follows. Let d[j] be
1.2024 +* initial value of lambda[j] (i.e. reduced cost of column j) in the
1.2025 +* transformed problem computed with formula (10) when pi[p] = 0. Then
1.2026 +* lambda[j] = d[j] + delta lambda[j], and from (12) it follows that
1.2027 +* lambda[j] becomes zero if:
1.2028 +*
1.2029 +* delta lambda[j] = - a[p,j] pi[p] = - d[j] ==>
1.2030 +* (13)
1.2031 +* pi[p] = d[j] / a[p,j].
1.2032 +*
1.2033 +* Therefore, the last column q, for which lambda[q] becomes zero, can
1.2034 +* be determined from the following condition:
1.2035 +*
1.2036 +* |d[q] / a[p,q]| = max |pi[p]| = max |d[j] / a[p,j]|, (14)
1.2037 +* j in D j in D
1.2038 +*
1.2039 +* where D is a set of columns j whose, reduced costs d[j] have invalid
1.2040 +* signs, i.e. violate condition (11). (Thus, if D is empty, solution
1.2041 +* to the original problem is the same as solution to the transformed
1.2042 +* problem, and no correction is needed as was noticed above.) In
1.2043 +* solution to the original problem column q is assigned status GLP_BS,
1.2044 +* since it replaces column of auxiliary variable of row p (becoming
1.2045 +* active) in the basis, and multiplier for row p is assigned its new
1.2046 +* value, which is pi[p] = d[q] / a[p,q]. Note that due to primal
1.2047 +* degeneracy values of all columns having non-zero coefficients in row
1.2048 +* p remain unchanged.
1.2049 +*
1.2050 +* RECOVERING INTERIOR-POINT SOLUTION
1.2051 +*
1.2052 +* Value of multiplier pi[p] in solution to the original problem is
1.2053 +* corrected in the same way as for basic solution. Values of all
1.2054 +* columns having non-zero coefficients in row p remain unchanged.
1.2055 +*
1.2056 +* RECOVERING MIP SOLUTION
1.2057 +*
1.2058 +* None needed. */
1.2059 +
1.2060 +struct forcing_col
1.2061 +{ /* column fixed on its bound by forcing row */
1.2062 + int j;
1.2063 + /* column reference number */
1.2064 + char stat;
1.2065 + /* original column status:
1.2066 + GLP_NL - fixed on lower bound
1.2067 + GLP_NU - fixed on upper bound */
1.2068 + double a;
1.2069 + /* constraint coefficient a[p,j] */
1.2070 + double c;
1.2071 + /* objective coefficient c[j] */
1.2072 + NPPLFE *ptr;
1.2073 + /* list of non-zero coefficients a[i,j], i != p */
1.2074 + struct forcing_col *next;
1.2075 + /* pointer to another column fixed by forcing row */
1.2076 +};
1.2077 +
1.2078 +struct forcing_row
1.2079 +{ /* forcing row */
1.2080 + int p;
1.2081 + /* row reference number */
1.2082 + char stat;
1.2083 + /* status assigned to the row if it becomes active:
1.2084 + GLP_NS - active equality constraint
1.2085 + GLP_NL - inequality constraint with lower bound active
1.2086 + GLP_NU - inequality constraint with upper bound active */
1.2087 + struct forcing_col *ptr;
1.2088 + /* list of all columns having non-zero constraint coefficient
1.2089 + a[p,j] in the forcing row */
1.2090 +};
1.2091 +
1.2092 +static int rcv_forcing_row(NPP *npp, void *info);
1.2093 +
1.2094 +int npp_forcing_row(NPP *npp, NPPROW *p, int at)
1.2095 +{ /* process forcing row */
1.2096 + struct forcing_row *info;
1.2097 + struct forcing_col *col = NULL;
1.2098 + NPPCOL *j;
1.2099 + NPPAIJ *apj, *aij;
1.2100 + NPPLFE *lfe;
1.2101 + double big;
1.2102 + xassert(at == 0 || at == 1);
1.2103 + /* determine maximal magnitude of the row coefficients */
1.2104 + big = 1.0;
1.2105 + for (apj = p->ptr; apj != NULL; apj = apj->r_next)
1.2106 + if (big < fabs(apj->val)) big = fabs(apj->val);
1.2107 + /* if there are too small coefficients in the row, transformation
1.2108 + should not be applied */
1.2109 + for (apj = p->ptr; apj != NULL; apj = apj->r_next)
1.2110 + if (fabs(apj->val) < 1e-7 * big) return 1;
1.2111 + /* create transformation stack entry */
1.2112 + info = npp_push_tse(npp,
1.2113 + rcv_forcing_row, sizeof(struct forcing_row));
1.2114 + info->p = p->i;
1.2115 + if (p->lb == p->ub)
1.2116 + { /* equality constraint */
1.2117 + info->stat = GLP_NS;
1.2118 + }
1.2119 + else if (at == 0)
1.2120 + { /* inequality constraint; case L[p] = U'[p] */
1.2121 + info->stat = GLP_NL;
1.2122 + xassert(p->lb != -DBL_MAX);
1.2123 + }
1.2124 + else /* at == 1 */
1.2125 + { /* inequality constraint; case U[p] = L'[p] */
1.2126 + info->stat = GLP_NU;
1.2127 + xassert(p->ub != +DBL_MAX);
1.2128 + }
1.2129 + info->ptr = NULL;
1.2130 + /* scan the forcing row, fix columns at corresponding bounds, and
1.2131 + save column information (the latter is not needed for MIP) */
1.2132 + for (apj = p->ptr; apj != NULL; apj = apj->r_next)
1.2133 + { /* column j has non-zero coefficient in the forcing row */
1.2134 + j = apj->col;
1.2135 + /* it must be non-fixed */
1.2136 + xassert(j->lb < j->ub);
1.2137 + /* allocate stack entry to save column information */
1.2138 + if (npp->sol != GLP_MIP)
1.2139 + { col = dmp_get_atom(npp->stack, sizeof(struct forcing_col));
1.2140 + col->j = j->j;
1.2141 + col->stat = -1; /* will be set below */
1.2142 + col->a = apj->val;
1.2143 + col->c = j->coef;
1.2144 + col->ptr = NULL;
1.2145 + col->next = info->ptr;
1.2146 + info->ptr = col;
1.2147 + }
1.2148 + /* fix column j */
1.2149 + if (at == 0 && apj->val < 0.0 || at != 0 && apj->val > 0.0)
1.2150 + { /* at its lower bound */
1.2151 + if (npp->sol != GLP_MIP)
1.2152 + col->stat = GLP_NL;
1.2153 + xassert(j->lb != -DBL_MAX);
1.2154 + j->ub = j->lb;
1.2155 + }
1.2156 + else
1.2157 + { /* at its upper bound */
1.2158 + if (npp->sol != GLP_MIP)
1.2159 + col->stat = GLP_NU;
1.2160 + xassert(j->ub != +DBL_MAX);
1.2161 + j->lb = j->ub;
1.2162 + }
1.2163 + /* save column coefficients a[i,j], i != p */
1.2164 + if (npp->sol != GLP_MIP)
1.2165 + { for (aij = j->ptr; aij != NULL; aij = aij->c_next)
1.2166 + { if (aij == apj) continue; /* skip a[p,j] */
1.2167 + lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE));
1.2168 + lfe->ref = aij->row->i;
1.2169 + lfe->val = aij->val;
1.2170 + lfe->next = col->ptr;
1.2171 + col->ptr = lfe;
1.2172 + }
1.2173 + }
1.2174 + }
1.2175 + /* make the row free (unbounded) */
1.2176 + p->lb = -DBL_MAX, p->ub = +DBL_MAX;
1.2177 + return 0;
1.2178 +}
1.2179 +
1.2180 +static int rcv_forcing_row(NPP *npp, void *_info)
1.2181 +{ /* recover forcing row */
1.2182 + struct forcing_row *info = _info;
1.2183 + struct forcing_col *col, *piv;
1.2184 + NPPLFE *lfe;
1.2185 + double d, big, temp;
1.2186 + if (npp->sol == GLP_MIP) goto done;
1.2187 + /* initially solution to the original problem is the same as
1.2188 + to the transformed problem, where row p is inactive constraint
1.2189 + with pi[p] = 0, and all columns are non-basic */
1.2190 + if (npp->sol == GLP_SOL)
1.2191 + { if (npp->r_stat[info->p] != GLP_BS)
1.2192 + { npp_error();
1.2193 + return 1;
1.2194 + }
1.2195 + for (col = info->ptr; col != NULL; col = col->next)
1.2196 + { if (npp->c_stat[col->j] != GLP_NS)
1.2197 + { npp_error();
1.2198 + return 1;
1.2199 + }
1.2200 + npp->c_stat[col->j] = col->stat; /* original status */
1.2201 + }
1.2202 + }
1.2203 + /* compute reduced costs d[j] for all columns with formula (10)
1.2204 + and store them in col.c instead objective coefficients */
1.2205 + for (col = info->ptr; col != NULL; col = col->next)
1.2206 + { d = col->c;
1.2207 + for (lfe = col->ptr; lfe != NULL; lfe = lfe->next)
1.2208 + d -= lfe->val * npp->r_pi[lfe->ref];
1.2209 + col->c = d;
1.2210 + }
1.2211 + /* consider columns j, whose multipliers lambda[j] has wrong
1.2212 + sign in solution to the transformed problem (where lambda[j] =
1.2213 + d[j]), and choose column q, whose multipler lambda[q] reaches
1.2214 + zero last on changing row multiplier pi[p]; see (14) */
1.2215 + piv = NULL, big = 0.0;
1.2216 + for (col = info->ptr; col != NULL; col = col->next)
1.2217 + { d = col->c; /* d[j] */
1.2218 + temp = fabs(d / col->a);
1.2219 + if (col->stat == GLP_NL)
1.2220 + { /* column j has active lower bound */
1.2221 + if (d < 0.0 && big < temp)
1.2222 + piv = col, big = temp;
1.2223 + }
1.2224 + else if (col->stat == GLP_NU)
1.2225 + { /* column j has active upper bound */
1.2226 + if (d > 0.0 && big < temp)
1.2227 + piv = col, big = temp;
1.2228 + }
1.2229 + else
1.2230 + { npp_error();
1.2231 + return 1;
1.2232 + }
1.2233 + }
1.2234 + /* if column q does not exist, no correction is needed */
1.2235 + if (piv != NULL)
1.2236 + { /* correct solution; row p becomes active constraint while
1.2237 + column q becomes basic */
1.2238 + if (npp->sol == GLP_SOL)
1.2239 + { npp->r_stat[info->p] = info->stat;
1.2240 + npp->c_stat[piv->j] = GLP_BS;
1.2241 + }
1.2242 + /* assign new value to row multiplier pi[p] = d[p] / a[p,q] */
1.2243 + npp->r_pi[info->p] = piv->c / piv->a;
1.2244 + }
1.2245 +done: return 0;
1.2246 +}
1.2247 +
1.2248 +/***********************************************************************
1.2249 +* NAME
1.2250 +*
1.2251 +* npp_analyze_row - perform general row analysis
1.2252 +*
1.2253 +* SYNOPSIS
1.2254 +*
1.2255 +* #include "glpnpp.h"
1.2256 +* int npp_analyze_row(NPP *npp, NPPROW *p);
1.2257 +*
1.2258 +* DESCRIPTION
1.2259 +*
1.2260 +* The routine npp_analyze_row performs analysis of row p of general
1.2261 +* format:
1.2262 +*
1.2263 +* L[p] <= sum a[p,j] x[j] <= U[p], (1)
1.2264 +* j
1.2265 +*
1.2266 +* l[j] <= x[j] <= u[j], (2)
1.2267 +*
1.2268 +* where L[p] <= U[p] and l[j] <= u[j] for all a[p,j] != 0.
1.2269 +*
1.2270 +* RETURNS
1.2271 +*
1.2272 +* 0x?0 - row lower bound does not exist or is redundant;
1.2273 +*
1.2274 +* 0x?1 - row lower bound can be active;
1.2275 +*
1.2276 +* 0x?2 - row lower bound is a forcing bound;
1.2277 +*
1.2278 +* 0x0? - row upper bound does not exist or is redundant;
1.2279 +*
1.2280 +* 0x1? - row upper bound can be active;
1.2281 +*
1.2282 +* 0x2? - row upper bound is a forcing bound;
1.2283 +*
1.2284 +* 0x33 - row bounds are inconsistent with column bounds.
1.2285 +*
1.2286 +* ALGORITHM
1.2287 +*
1.2288 +* Analysis of row (1) is based on analysis of its implied lower and
1.2289 +* upper bounds, which are determined by bounds of corresponding columns
1.2290 +* (variables) as follows:
1.2291 +*
1.2292 +* L'[p] = inf sum a[p,j] x[j] =
1.2293 +* j
1.2294 +* (3)
1.2295 +* = sum a[p,j] l[j] + sum a[p,j] u[j],
1.2296 +* j in Jp j in Jn
1.2297 +*
1.2298 +* U'[p] = sup sum a[p,j] x[j] =
1.2299 +* (4)
1.2300 +* = sum a[p,j] u[j] + sum a[p,j] l[j],
1.2301 +* j in Jp j in Jn
1.2302 +*
1.2303 +* Jp = {j: a[p,j] > 0}, Jn = {j: a[p,j] < 0}. (5)
1.2304 +*
1.2305 +* (Note that bounds of all columns in row p are assumed to be correct,
1.2306 +* so L'[p] <= U'[p].)
1.2307 +*
1.2308 +* Analysis of row lower bound L[p] includes the following cases:
1.2309 +*
1.2310 +* 1) if L[p] > U'[p] + eps, where eps is an absolute tolerance for row
1.2311 +* value, row lower bound L[p] and implied row upper bound U'[p] are
1.2312 +* inconsistent, ergo, the problem has no primal feasible solution;
1.2313 +*
1.2314 +* 2) if U'[p] - eps <= L[p] <= U'[p] + eps, i.e. if L[p] =~ U'[p],
1.2315 +* the row is a forcing row on its lower bound (see description of
1.2316 +* the routine npp_forcing_row);
1.2317 +*
1.2318 +* 3) if L[p] > L'[p] + eps, row lower bound L[p] can be active (this
1.2319 +* conclusion does not account other rows in the problem);
1.2320 +*
1.2321 +* 4) if L[p] <= L'[p] + eps, row lower bound L[p] cannot be active, so
1.2322 +* it is redundant and can be removed (replaced by -oo).
1.2323 +*
1.2324 +* Analysis of row upper bound U[p] is performed in a similar way and
1.2325 +* includes the following cases:
1.2326 +*
1.2327 +* 1) if U[p] < L'[p] - eps, row upper bound U[p] and implied row lower
1.2328 +* bound L'[p] are inconsistent, ergo the problem has no primal
1.2329 +* feasible solution;
1.2330 +*
1.2331 +* 2) if L'[p] - eps <= U[p] <= L'[p] + eps, i.e. if U[p] =~ L'[p],
1.2332 +* the row is a forcing row on its upper bound (see description of
1.2333 +* the routine npp_forcing_row);
1.2334 +*
1.2335 +* 3) if U[p] < U'[p] - eps, row upper bound U[p] can be active (this
1.2336 +* conclusion does not account other rows in the problem);
1.2337 +*
1.2338 +* 4) if U[p] >= U'[p] - eps, row upper bound U[p] cannot be active, so
1.2339 +* it is redundant and can be removed (replaced by +oo). */
1.2340 +
1.2341 +int npp_analyze_row(NPP *npp, NPPROW *p)
1.2342 +{ /* perform general row analysis */
1.2343 + NPPAIJ *aij;
1.2344 + int ret = 0x00;
1.2345 + double l, u, eps;
1.2346 + xassert(npp == npp);
1.2347 + /* compute implied lower bound L'[p]; see (3) */
1.2348 + l = 0.0;
1.2349 + for (aij = p->ptr; aij != NULL; aij = aij->r_next)
1.2350 + { if (aij->val > 0.0)
1.2351 + { if (aij->col->lb == -DBL_MAX)
1.2352 + { l = -DBL_MAX;
1.2353 + break;
1.2354 + }
1.2355 + l += aij->val * aij->col->lb;
1.2356 + }
1.2357 + else /* aij->val < 0.0 */
1.2358 + { if (aij->col->ub == +DBL_MAX)
1.2359 + { l = -DBL_MAX;
1.2360 + break;
1.2361 + }
1.2362 + l += aij->val * aij->col->ub;
1.2363 + }
1.2364 + }
1.2365 + /* compute implied upper bound U'[p]; see (4) */
1.2366 + u = 0.0;
1.2367 + for (aij = p->ptr; aij != NULL; aij = aij->r_next)
1.2368 + { if (aij->val > 0.0)
1.2369 + { if (aij->col->ub == +DBL_MAX)
1.2370 + { u = +DBL_MAX;
1.2371 + break;
1.2372 + }
1.2373 + u += aij->val * aij->col->ub;
1.2374 + }
1.2375 + else /* aij->val < 0.0 */
1.2376 + { if (aij->col->lb == -DBL_MAX)
1.2377 + { u = +DBL_MAX;
1.2378 + break;
1.2379 + }
1.2380 + u += aij->val * aij->col->lb;
1.2381 + }
1.2382 + }
1.2383 + /* column bounds are assumed correct, so L'[p] <= U'[p] */
1.2384 + /* check if row lower bound is consistent */
1.2385 + if (p->lb != -DBL_MAX)
1.2386 + { eps = 1e-3 + 1e-6 * fabs(p->lb);
1.2387 + if (p->lb - eps > u)
1.2388 + { ret = 0x33;
1.2389 + goto done;
1.2390 + }
1.2391 + }
1.2392 + /* check if row upper bound is consistent */
1.2393 + if (p->ub != +DBL_MAX)
1.2394 + { eps = 1e-3 + 1e-6 * fabs(p->ub);
1.2395 + if (p->ub + eps < l)
1.2396 + { ret = 0x33;
1.2397 + goto done;
1.2398 + }
1.2399 + }
1.2400 + /* check if row lower bound can be active/forcing */
1.2401 + if (p->lb != -DBL_MAX)
1.2402 + { eps = 1e-9 + 1e-12 * fabs(p->lb);
1.2403 + if (p->lb - eps > l)
1.2404 + { if (p->lb + eps <= u)
1.2405 + ret |= 0x01;
1.2406 + else
1.2407 + ret |= 0x02;
1.2408 + }
1.2409 + }
1.2410 + /* check if row upper bound can be active/forcing */
1.2411 + if (p->ub != +DBL_MAX)
1.2412 + { eps = 1e-9 + 1e-12 * fabs(p->ub);
1.2413 + if (p->ub + eps < u)
1.2414 + { /* check if the upper bound is forcing */
1.2415 + if (p->ub - eps >= l)
1.2416 + ret |= 0x10;
1.2417 + else
1.2418 + ret |= 0x20;
1.2419 + }
1.2420 + }
1.2421 +done: return ret;
1.2422 +}
1.2423 +
1.2424 +/***********************************************************************
1.2425 +* NAME
1.2426 +*
1.2427 +* npp_inactive_bound - remove row lower/upper inactive bound
1.2428 +*
1.2429 +* SYNOPSIS
1.2430 +*
1.2431 +* #include "glpnpp.h"
1.2432 +* void npp_inactive_bound(NPP *npp, NPPROW *p, int which);
1.2433 +*
1.2434 +* DESCRIPTION
1.2435 +*
1.2436 +* The routine npp_inactive_bound removes lower (if which = 0) or upper
1.2437 +* (if which = 1) bound of row p:
1.2438 +*
1.2439 +* L[p] <= sum a[p,j] x[j] <= U[p],
1.2440 +*
1.2441 +* which (bound) is assumed to be redundant.
1.2442 +*
1.2443 +* PROBLEM TRANSFORMATION
1.2444 +*
1.2445 +* If which = 0, current lower bound L[p] of row p is assigned -oo.
1.2446 +* If which = 1, current upper bound U[p] of row p is assigned +oo.
1.2447 +*
1.2448 +* RECOVERING BASIC SOLUTION
1.2449 +*
1.2450 +* If in solution to the transformed problem row p is inactive
1.2451 +* constraint (GLP_BS), its status is not changed in solution to the
1.2452 +* original problem. Otherwise, status of row p in solution to the
1.2453 +* original problem is defined by its type before transformation and
1.2454 +* its status in solution to the transformed problem as follows:
1.2455 +*
1.2456 +* +---------------------+-------+---------------+---------------+
1.2457 +* | Row | Flag | Row status in | Row status in |
1.2458 +* | type | which | transfmd soln | original soln |
1.2459 +* +---------------------+-------+---------------+---------------+
1.2460 +* | sum >= L[p] | 0 | GLP_NF | GLP_NL |
1.2461 +* | sum <= U[p] | 1 | GLP_NF | GLP_NU |
1.2462 +* | L[p] <= sum <= U[p] | 0 | GLP_NU | GLP_NU |
1.2463 +* | L[p] <= sum <= U[p] | 1 | GLP_NL | GLP_NL |
1.2464 +* | sum = L[p] = U[p] | 0 | GLP_NU | GLP_NS |
1.2465 +* | sum = L[p] = U[p] | 1 | GLP_NL | GLP_NS |
1.2466 +* +---------------------+-------+---------------+---------------+
1.2467 +*
1.2468 +* RECOVERING INTERIOR-POINT SOLUTION
1.2469 +*
1.2470 +* None needed.
1.2471 +*
1.2472 +* RECOVERING MIP SOLUTION
1.2473 +*
1.2474 +* None needed. */
1.2475 +
1.2476 +struct inactive_bound
1.2477 +{ /* row inactive bound */
1.2478 + int p;
1.2479 + /* row reference number */
1.2480 + char stat;
1.2481 + /* row status (if active constraint) */
1.2482 +};
1.2483 +
1.2484 +static int rcv_inactive_bound(NPP *npp, void *info);
1.2485 +
1.2486 +void npp_inactive_bound(NPP *npp, NPPROW *p, int which)
1.2487 +{ /* remove row lower/upper inactive bound */
1.2488 + struct inactive_bound *info;
1.2489 + if (npp->sol == GLP_SOL)
1.2490 + { /* create transformation stack entry */
1.2491 + info = npp_push_tse(npp,
1.2492 + rcv_inactive_bound, sizeof(struct inactive_bound));
1.2493 + info->p = p->i;
1.2494 + if (p->ub == +DBL_MAX)
1.2495 + info->stat = GLP_NL;
1.2496 + else if (p->lb == -DBL_MAX)
1.2497 + info->stat = GLP_NU;
1.2498 + else if (p->lb != p->ub)
1.2499 + info->stat = (char)(which == 0 ? GLP_NU : GLP_NL);
1.2500 + else
1.2501 + info->stat = GLP_NS;
1.2502 + }
1.2503 + /* remove row inactive bound */
1.2504 + if (which == 0)
1.2505 + { xassert(p->lb != -DBL_MAX);
1.2506 + p->lb = -DBL_MAX;
1.2507 + }
1.2508 + else if (which == 1)
1.2509 + { xassert(p->ub != +DBL_MAX);
1.2510 + p->ub = +DBL_MAX;
1.2511 + }
1.2512 + else
1.2513 + xassert(which != which);
1.2514 + return;
1.2515 +}
1.2516 +
1.2517 +static int rcv_inactive_bound(NPP *npp, void *_info)
1.2518 +{ /* recover row status */
1.2519 + struct inactive_bound *info = _info;
1.2520 + if (npp->sol != GLP_SOL)
1.2521 + { npp_error();
1.2522 + return 1;
1.2523 + }
1.2524 + if (npp->r_stat[info->p] == GLP_BS)
1.2525 + npp->r_stat[info->p] = GLP_BS;
1.2526 + else
1.2527 + npp->r_stat[info->p] = info->stat;
1.2528 + return 0;
1.2529 +}
1.2530 +
1.2531 +/***********************************************************************
1.2532 +* NAME
1.2533 +*
1.2534 +* npp_implied_bounds - determine implied column bounds
1.2535 +*
1.2536 +* SYNOPSIS
1.2537 +*
1.2538 +* #include "glpnpp.h"
1.2539 +* void npp_implied_bounds(NPP *npp, NPPROW *p);
1.2540 +*
1.2541 +* DESCRIPTION
1.2542 +*
1.2543 +* The routine npp_implied_bounds inspects general row (constraint) p:
1.2544 +*
1.2545 +* L[p] <= sum a[p,j] x[j] <= U[p], (1)
1.2546 +*
1.2547 +* l[j] <= x[j] <= u[j], (2)
1.2548 +*
1.2549 +* where L[p] <= U[p] and l[j] <= u[j] for all a[p,j] != 0, to compute
1.2550 +* implied bounds of columns (variables x[j]) in this row.
1.2551 +*
1.2552 +* The routine stores implied column bounds l'[j] and u'[j] in column
1.2553 +* descriptors (NPPCOL); it does not change current column bounds l[j]
1.2554 +* and u[j]. (Implied column bounds can be then used to strengthen the
1.2555 +* current column bounds; see the routines npp_implied_lower and
1.2556 +* npp_implied_upper).
1.2557 +*
1.2558 +* ALGORITHM
1.2559 +*
1.2560 +* Current column bounds (2) define implied lower and upper bounds of
1.2561 +* row (1) as follows:
1.2562 +*
1.2563 +* L'[p] = inf sum a[p,j] x[j] =
1.2564 +* j
1.2565 +* (3)
1.2566 +* = sum a[p,j] l[j] + sum a[p,j] u[j],
1.2567 +* j in Jp j in Jn
1.2568 +*
1.2569 +* U'[p] = sup sum a[p,j] x[j] =
1.2570 +* (4)
1.2571 +* = sum a[p,j] u[j] + sum a[p,j] l[j],
1.2572 +* j in Jp j in Jn
1.2573 +*
1.2574 +* Jp = {j: a[p,j] > 0}, Jn = {j: a[p,j] < 0}. (5)
1.2575 +*
1.2576 +* (Note that bounds of all columns in row p are assumed to be correct,
1.2577 +* so L'[p] <= U'[p].)
1.2578 +*
1.2579 +* If L[p] > L'[p] and/or U[p] < U'[p], the lower and/or upper bound of
1.2580 +* row (1) can be active, in which case such row defines implied bounds
1.2581 +* of its variables.
1.2582 +*
1.2583 +* Let x[k] be some variable having in row (1) coefficient a[p,k] != 0.
1.2584 +* Consider a case when row lower bound can be active (L[p] > L'[p]):
1.2585 +*
1.2586 +* sum a[p,j] x[j] >= L[p] ==>
1.2587 +* j
1.2588 +*
1.2589 +* sum a[p,j] x[j] + a[p,k] x[k] >= L[p] ==>
1.2590 +* j!=k
1.2591 +* (6)
1.2592 +* a[p,k] x[k] >= L[p] - sum a[p,j] x[j] ==>
1.2593 +* j!=k
1.2594 +*
1.2595 +* a[p,k] x[k] >= L[p,k],
1.2596 +*
1.2597 +* where
1.2598 +*
1.2599 +* L[p,k] = inf(L[p] - sum a[p,j] x[j]) =
1.2600 +* j!=k
1.2601 +*
1.2602 +* = L[p] - sup sum a[p,j] x[j] = (7)
1.2603 +* j!=k
1.2604 +*
1.2605 +* = L[p] - sum a[p,j] u[j] - sum a[p,j] l[j].
1.2606 +* j in Jp\{k} j in Jn\{k}
1.2607 +*
1.2608 +* Thus:
1.2609 +*
1.2610 +* x[k] >= l'[k] = L[p,k] / a[p,k], if a[p,k] > 0, (8)
1.2611 +*
1.2612 +* x[k] <= u'[k] = L[p,k] / a[p,k], if a[p,k] < 0. (9)
1.2613 +*
1.2614 +* where l'[k] and u'[k] are implied lower and upper bounds of variable
1.2615 +* x[k], resp.
1.2616 +*
1.2617 +* Now consider a similar case when row upper bound can be active
1.2618 +* (U[p] < U'[p]):
1.2619 +*
1.2620 +* sum a[p,j] x[j] <= U[p] ==>
1.2621 +* j
1.2622 +*
1.2623 +* sum a[p,j] x[j] + a[p,k] x[k] <= U[p] ==>
1.2624 +* j!=k
1.2625 +* (10)
1.2626 +* a[p,k] x[k] <= U[p] - sum a[p,j] x[j] ==>
1.2627 +* j!=k
1.2628 +*
1.2629 +* a[p,k] x[k] <= U[p,k],
1.2630 +*
1.2631 +* where:
1.2632 +*
1.2633 +* U[p,k] = sup(U[p] - sum a[p,j] x[j]) =
1.2634 +* j!=k
1.2635 +*
1.2636 +* = U[p] - inf sum a[p,j] x[j] = (11)
1.2637 +* j!=k
1.2638 +*
1.2639 +* = U[p] - sum a[p,j] l[j] - sum a[p,j] u[j].
1.2640 +* j in Jp\{k} j in Jn\{k}
1.2641 +*
1.2642 +* Thus:
1.2643 +*
1.2644 +* x[k] <= u'[k] = U[p,k] / a[p,k], if a[p,k] > 0, (12)
1.2645 +*
1.2646 +* x[k] >= l'[k] = U[p,k] / a[p,k], if a[p,k] < 0. (13)
1.2647 +*
1.2648 +* Note that in formulae (8), (9), (12), and (13) coefficient a[p,k]
1.2649 +* must not be too small in magnitude relatively to other non-zero
1.2650 +* coefficients in row (1), i.e. the following condition must hold:
1.2651 +*
1.2652 +* |a[p,k]| >= eps * max(1, |a[p,j]|), (14)
1.2653 +* j
1.2654 +*
1.2655 +* where eps is a relative tolerance for constraint coefficients.
1.2656 +* Otherwise the implied column bounds can be numerical inreliable. For
1.2657 +* example, using formula (8) for the following inequality constraint:
1.2658 +*
1.2659 +* 1e-12 x1 - x2 - x3 >= 0,
1.2660 +*
1.2661 +* where x1 >= -1, x2, x3, >= 0, may lead to numerically unreliable
1.2662 +* conclusion that x1 >= 0.
1.2663 +*
1.2664 +* Using formulae (8), (9), (12), and (13) to compute implied bounds
1.2665 +* for one variable requires |J| operations, where J = {j: a[p,j] != 0},
1.2666 +* because this needs computing L[p,k] and U[p,k]. Thus, computing
1.2667 +* implied bounds for all variables in row (1) would require |J|^2
1.2668 +* operations, that is not a good technique. However, the total number
1.2669 +* of operations can be reduced to |J| as follows.
1.2670 +*
1.2671 +* Let a[p,k] > 0. Then from (7) and (11) we have:
1.2672 +*
1.2673 +* L[p,k] = L[p] - (U'[p] - a[p,k] u[k]) =
1.2674 +*
1.2675 +* = L[p] - U'[p] + a[p,k] u[k],
1.2676 +*
1.2677 +* U[p,k] = U[p] - (L'[p] - a[p,k] l[k]) =
1.2678 +*
1.2679 +* = U[p] - L'[p] + a[p,k] l[k],
1.2680 +*
1.2681 +* where L'[p] and U'[p] are implied row lower and upper bounds defined
1.2682 +* by formulae (3) and (4). Substituting these expressions into (8) and
1.2683 +* (12) gives:
1.2684 +*
1.2685 +* l'[k] = L[p,k] / a[p,k] = u[k] + (L[p] - U'[p]) / a[p,k], (15)
1.2686 +*
1.2687 +* u'[k] = U[p,k] / a[p,k] = l[k] + (U[p] - L'[p]) / a[p,k]. (16)
1.2688 +*
1.2689 +* Similarly, if a[p,k] < 0, according to (7) and (11) we have:
1.2690 +*
1.2691 +* L[p,k] = L[p] - (U'[p] - a[p,k] l[k]) =
1.2692 +*
1.2693 +* = L[p] - U'[p] + a[p,k] l[k],
1.2694 +*
1.2695 +* U[p,k] = U[p] - (L'[p] - a[p,k] u[k]) =
1.2696 +*
1.2697 +* = U[p] - L'[p] + a[p,k] u[k],
1.2698 +*
1.2699 +* and substituting these expressions into (8) and (12) gives:
1.2700 +*
1.2701 +* l'[k] = U[p,k] / a[p,k] = u[k] + (U[p] - L'[p]) / a[p,k], (17)
1.2702 +*
1.2703 +* u'[k] = L[p,k] / a[p,k] = l[k] + (L[p] - U'[p]) / a[p,k]. (18)
1.2704 +*
1.2705 +* Note that formulae (15)-(18) can be used only if L'[p] and U'[p]
1.2706 +* exist. However, if for some variable x[j] it happens that l[j] = -oo
1.2707 +* and/or u[j] = +oo, values of L'[p] (if a[p,j] > 0) and/or U'[p] (if
1.2708 +* a[p,j] < 0) are undefined. Consider, therefore, the most general
1.2709 +* situation, when some column bounds (2) may not exist.
1.2710 +*
1.2711 +* Let:
1.2712 +*
1.2713 +* J' = {j : (a[p,j] > 0 and l[j] = -oo) or
1.2714 +* (19)
1.2715 +* (a[p,j] < 0 and u[j] = +oo)}.
1.2716 +*
1.2717 +* Then (assuming that row upper bound U[p] can be active) the following
1.2718 +* three cases are possible:
1.2719 +*
1.2720 +* 1) |J'| = 0. In this case L'[p] exists, thus, for all variables x[j]
1.2721 +* in row (1) we can use formulae (16) and (17);
1.2722 +*
1.2723 +* 2) J' = {k}. In this case L'[p] = -oo, however, U[p,k] (11) exists,
1.2724 +* so for variable x[k] we can use formulae (12) and (13). Note that
1.2725 +* for all other variables x[j] (j != k) l'[j] = -oo (if a[p,j] < 0)
1.2726 +* or u'[j] = +oo (if a[p,j] > 0);
1.2727 +*
1.2728 +* 3) |J'| > 1. In this case for all variables x[j] in row [1] we have
1.2729 +* l'[j] = -oo (if a[p,j] < 0) or u'[j] = +oo (if a[p,j] > 0).
1.2730 +*
1.2731 +* Similarly, let:
1.2732 +*
1.2733 +* J'' = {j : (a[p,j] > 0 and u[j] = +oo) or
1.2734 +* (20)
1.2735 +* (a[p,j] < 0 and l[j] = -oo)}.
1.2736 +*
1.2737 +* Then (assuming that row lower bound L[p] can be active) the following
1.2738 +* three cases are possible:
1.2739 +*
1.2740 +* 1) |J''| = 0. In this case U'[p] exists, thus, for all variables x[j]
1.2741 +* in row (1) we can use formulae (15) and (18);
1.2742 +*
1.2743 +* 2) J'' = {k}. In this case U'[p] = +oo, however, L[p,k] (7) exists,
1.2744 +* so for variable x[k] we can use formulae (8) and (9). Note that
1.2745 +* for all other variables x[j] (j != k) l'[j] = -oo (if a[p,j] > 0)
1.2746 +* or u'[j] = +oo (if a[p,j] < 0);
1.2747 +*
1.2748 +* 3) |J''| > 1. In this case for all variables x[j] in row (1) we have
1.2749 +* l'[j] = -oo (if a[p,j] > 0) or u'[j] = +oo (if a[p,j] < 0). */
1.2750 +
1.2751 +void npp_implied_bounds(NPP *npp, NPPROW *p)
1.2752 +{ NPPAIJ *apj, *apk;
1.2753 + double big, eps, temp;
1.2754 + xassert(npp == npp);
1.2755 + /* initialize implied bounds for all variables and determine
1.2756 + maximal magnitude of row coefficients a[p,j] */
1.2757 + big = 1.0;
1.2758 + for (apj = p->ptr; apj != NULL; apj = apj->r_next)
1.2759 + { apj->col->ll.ll = -DBL_MAX, apj->col->uu.uu = +DBL_MAX;
1.2760 + if (big < fabs(apj->val)) big = fabs(apj->val);
1.2761 + }
1.2762 + eps = 1e-6 * big;
1.2763 + /* process row lower bound (assuming that it can be active) */
1.2764 + if (p->lb != -DBL_MAX)
1.2765 + { apk = NULL;
1.2766 + for (apj = p->ptr; apj != NULL; apj = apj->r_next)
1.2767 + { if (apj->val > 0.0 && apj->col->ub == +DBL_MAX ||
1.2768 + apj->val < 0.0 && apj->col->lb == -DBL_MAX)
1.2769 + { if (apk == NULL)
1.2770 + apk = apj;
1.2771 + else
1.2772 + goto skip1;
1.2773 + }
1.2774 + }
1.2775 + /* if a[p,k] = NULL then |J'| = 0 else J' = { k } */
1.2776 + temp = p->lb;
1.2777 + for (apj = p->ptr; apj != NULL; apj = apj->r_next)
1.2778 + { if (apj == apk)
1.2779 + /* skip a[p,k] */;
1.2780 + else if (apj->val > 0.0)
1.2781 + temp -= apj->val * apj->col->ub;
1.2782 + else /* apj->val < 0.0 */
1.2783 + temp -= apj->val * apj->col->lb;
1.2784 + }
1.2785 + /* compute column implied bounds */
1.2786 + if (apk == NULL)
1.2787 + { /* temp = L[p] - U'[p] */
1.2788 + for (apj = p->ptr; apj != NULL; apj = apj->r_next)
1.2789 + { if (apj->val >= +eps)
1.2790 + { /* l'[j] := u[j] + (L[p] - U'[p]) / a[p,j] */
1.2791 + apj->col->ll.ll = apj->col->ub + temp / apj->val;
1.2792 + }
1.2793 + else if (apj->val <= -eps)
1.2794 + { /* u'[j] := l[j] + (L[p] - U'[p]) / a[p,j] */
1.2795 + apj->col->uu.uu = apj->col->lb + temp / apj->val;
1.2796 + }
1.2797 + }
1.2798 + }
1.2799 + else
1.2800 + { /* temp = L[p,k] */
1.2801 + if (apk->val >= +eps)
1.2802 + { /* l'[k] := L[p,k] / a[p,k] */
1.2803 + apk->col->ll.ll = temp / apk->val;
1.2804 + }
1.2805 + else if (apk->val <= -eps)
1.2806 + { /* u'[k] := L[p,k] / a[p,k] */
1.2807 + apk->col->uu.uu = temp / apk->val;
1.2808 + }
1.2809 + }
1.2810 +skip1: ;
1.2811 + }
1.2812 + /* process row upper bound (assuming that it can be active) */
1.2813 + if (p->ub != +DBL_MAX)
1.2814 + { apk = NULL;
1.2815 + for (apj = p->ptr; apj != NULL; apj = apj->r_next)
1.2816 + { if (apj->val > 0.0 && apj->col->lb == -DBL_MAX ||
1.2817 + apj->val < 0.0 && apj->col->ub == +DBL_MAX)
1.2818 + { if (apk == NULL)
1.2819 + apk = apj;
1.2820 + else
1.2821 + goto skip2;
1.2822 + }
1.2823 + }
1.2824 + /* if a[p,k] = NULL then |J''| = 0 else J'' = { k } */
1.2825 + temp = p->ub;
1.2826 + for (apj = p->ptr; apj != NULL; apj = apj->r_next)
1.2827 + { if (apj == apk)
1.2828 + /* skip a[p,k] */;
1.2829 + else if (apj->val > 0.0)
1.2830 + temp -= apj->val * apj->col->lb;
1.2831 + else /* apj->val < 0.0 */
1.2832 + temp -= apj->val * apj->col->ub;
1.2833 + }
1.2834 + /* compute column implied bounds */
1.2835 + if (apk == NULL)
1.2836 + { /* temp = U[p] - L'[p] */
1.2837 + for (apj = p->ptr; apj != NULL; apj = apj->r_next)
1.2838 + { if (apj->val >= +eps)
1.2839 + { /* u'[j] := l[j] + (U[p] - L'[p]) / a[p,j] */
1.2840 + apj->col->uu.uu = apj->col->lb + temp / apj->val;
1.2841 + }
1.2842 + else if (apj->val <= -eps)
1.2843 + { /* l'[j] := u[j] + (U[p] - L'[p]) / a[p,j] */
1.2844 + apj->col->ll.ll = apj->col->ub + temp / apj->val;
1.2845 + }
1.2846 + }
1.2847 + }
1.2848 + else
1.2849 + { /* temp = U[p,k] */
1.2850 + if (apk->val >= +eps)
1.2851 + { /* u'[k] := U[p,k] / a[p,k] */
1.2852 + apk->col->uu.uu = temp / apk->val;
1.2853 + }
1.2854 + else if (apk->val <= -eps)
1.2855 + { /* l'[k] := U[p,k] / a[p,k] */
1.2856 + apk->col->ll.ll = temp / apk->val;
1.2857 + }
1.2858 + }
1.2859 +skip2: ;
1.2860 + }
1.2861 + return;
1.2862 +}
1.2863 +
1.2864 +/* eof */