1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpnpp02.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,1433 @@
1.4 +/* glpnpp02.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_free_row - process free (unbounded) row
1.34 +*
1.35 +* SYNOPSIS
1.36 +*
1.37 +* #include "glpnpp.h"
1.38 +* void npp_free_row(NPP *npp, NPPROW *p);
1.39 +*
1.40 +* DESCRIPTION
1.41 +*
1.42 +* The routine npp_free_row processes row p, which is free (i.e. has
1.43 +* no finite bounds):
1.44 +*
1.45 +* -inf < sum a[p,j] x[j] < +inf. (1)
1.46 +* j
1.47 +*
1.48 +* PROBLEM TRANSFORMATION
1.49 +*
1.50 +* Constraint (1) cannot be active, so it is redundant and can be
1.51 +* removed from the original problem.
1.52 +*
1.53 +* Removing row p leads to removing a column of multiplier pi[p] for
1.54 +* this row in the dual system. Since row p has no bounds, pi[p] = 0,
1.55 +* so removing the column does not affect the dual solution.
1.56 +*
1.57 +* RECOVERING BASIC SOLUTION
1.58 +*
1.59 +* In solution to the original problem row p is inactive constraint,
1.60 +* so it is assigned status GLP_BS, and multiplier pi[p] is assigned
1.61 +* zero value.
1.62 +*
1.63 +* RECOVERING INTERIOR-POINT SOLUTION
1.64 +*
1.65 +* In solution to the original problem row p is inactive constraint,
1.66 +* so its multiplier pi[p] is assigned zero value.
1.67 +*
1.68 +* RECOVERING MIP SOLUTION
1.69 +*
1.70 +* None needed. */
1.71 +
1.72 +struct free_row
1.73 +{ /* free (unbounded) row */
1.74 + int p;
1.75 + /* row reference number */
1.76 +};
1.77 +
1.78 +static int rcv_free_row(NPP *npp, void *info);
1.79 +
1.80 +void npp_free_row(NPP *npp, NPPROW *p)
1.81 +{ /* process free (unbounded) row */
1.82 + struct free_row *info;
1.83 + /* the row must be free */
1.84 + xassert(p->lb == -DBL_MAX && p->ub == +DBL_MAX);
1.85 + /* create transformation stack entry */
1.86 + info = npp_push_tse(npp,
1.87 + rcv_free_row, sizeof(struct free_row));
1.88 + info->p = p->i;
1.89 + /* remove the row from the problem */
1.90 + npp_del_row(npp, p);
1.91 + return;
1.92 +}
1.93 +
1.94 +static int rcv_free_row(NPP *npp, void *_info)
1.95 +{ /* recover free (unbounded) row */
1.96 + struct free_row *info = _info;
1.97 + if (npp->sol == GLP_SOL)
1.98 + npp->r_stat[info->p] = GLP_BS;
1.99 + if (npp->sol != GLP_MIP)
1.100 + npp->r_pi[info->p] = 0.0;
1.101 + return 0;
1.102 +}
1.103 +
1.104 +/***********************************************************************
1.105 +* NAME
1.106 +*
1.107 +* npp_geq_row - process row of 'not less than' type
1.108 +*
1.109 +* SYNOPSIS
1.110 +*
1.111 +* #include "glpnpp.h"
1.112 +* void npp_geq_row(NPP *npp, NPPROW *p);
1.113 +*
1.114 +* DESCRIPTION
1.115 +*
1.116 +* The routine npp_geq_row processes row p, which is 'not less than'
1.117 +* inequality constraint:
1.118 +*
1.119 +* L[p] <= sum a[p,j] x[j] (<= U[p]), (1)
1.120 +* j
1.121 +*
1.122 +* where L[p] < U[p], and upper bound may not exist (U[p] = +oo).
1.123 +*
1.124 +* PROBLEM TRANSFORMATION
1.125 +*
1.126 +* Constraint (1) can be replaced by equality constraint:
1.127 +*
1.128 +* sum a[p,j] x[j] - s = L[p], (2)
1.129 +* j
1.130 +*
1.131 +* where
1.132 +*
1.133 +* 0 <= s (<= U[p] - L[p]) (3)
1.134 +*
1.135 +* is a non-negative surplus variable.
1.136 +*
1.137 +* Since in the primal system there appears column s having the only
1.138 +* non-zero coefficient in row p, in the dual system there appears a
1.139 +* new row:
1.140 +*
1.141 +* (-1) pi[p] + lambda = 0, (4)
1.142 +*
1.143 +* where (-1) is coefficient of column s in row p, pi[p] is multiplier
1.144 +* of row p, lambda is multiplier of column q, 0 is coefficient of
1.145 +* column s in the objective row.
1.146 +*
1.147 +* RECOVERING BASIC SOLUTION
1.148 +*
1.149 +* Status of row p in solution to the original problem is determined
1.150 +* by its status and status of column q in solution to the transformed
1.151 +* problem as follows:
1.152 +*
1.153 +* +--------------------------------------+------------------+
1.154 +* | Transformed problem | Original problem |
1.155 +* +-----------------+--------------------+------------------+
1.156 +* | Status of row p | Status of column s | Status of row p |
1.157 +* +-----------------+--------------------+------------------+
1.158 +* | GLP_BS | GLP_BS | N/A |
1.159 +* | GLP_BS | GLP_NL | GLP_BS |
1.160 +* | GLP_BS | GLP_NU | GLP_BS |
1.161 +* | GLP_NS | GLP_BS | GLP_BS |
1.162 +* | GLP_NS | GLP_NL | GLP_NL |
1.163 +* | GLP_NS | GLP_NU | GLP_NU |
1.164 +* +-----------------+--------------------+------------------+
1.165 +*
1.166 +* Value of row multiplier pi[p] in solution to the original problem
1.167 +* is the same as in solution to the transformed problem.
1.168 +*
1.169 +* 1. In solution to the transformed problem row p and column q cannot
1.170 +* be basic at the same time; otherwise the basis matrix would have
1.171 +* two linear dependent columns: unity column of auxiliary variable
1.172 +* of row p and unity column of variable s.
1.173 +*
1.174 +* 2. Though in the transformed problem row p is equality constraint,
1.175 +* it may be basic due to primal degenerate solution.
1.176 +*
1.177 +* RECOVERING INTERIOR-POINT SOLUTION
1.178 +*
1.179 +* Value of row multiplier pi[p] in solution to the original problem
1.180 +* is the same as in solution to the transformed problem.
1.181 +*
1.182 +* RECOVERING MIP SOLUTION
1.183 +*
1.184 +* None needed. */
1.185 +
1.186 +struct ineq_row
1.187 +{ /* inequality constraint row */
1.188 + int p;
1.189 + /* row reference number */
1.190 + int s;
1.191 + /* column reference number for slack/surplus variable */
1.192 +};
1.193 +
1.194 +static int rcv_geq_row(NPP *npp, void *info);
1.195 +
1.196 +void npp_geq_row(NPP *npp, NPPROW *p)
1.197 +{ /* process row of 'not less than' type */
1.198 + struct ineq_row *info;
1.199 + NPPCOL *s;
1.200 + /* the row must have lower bound */
1.201 + xassert(p->lb != -DBL_MAX);
1.202 + xassert(p->lb < p->ub);
1.203 + /* create column for surplus variable */
1.204 + s = npp_add_col(npp);
1.205 + s->lb = 0.0;
1.206 + s->ub = (p->ub == +DBL_MAX ? +DBL_MAX : p->ub - p->lb);
1.207 + /* and add it to the transformed problem */
1.208 + npp_add_aij(npp, p, s, -1.0);
1.209 + /* create transformation stack entry */
1.210 + info = npp_push_tse(npp,
1.211 + rcv_geq_row, sizeof(struct ineq_row));
1.212 + info->p = p->i;
1.213 + info->s = s->j;
1.214 + /* replace the row by equality constraint */
1.215 + p->ub = p->lb;
1.216 + return;
1.217 +}
1.218 +
1.219 +static int rcv_geq_row(NPP *npp, void *_info)
1.220 +{ /* recover row of 'not less than' type */
1.221 + struct ineq_row *info = _info;
1.222 + if (npp->sol == GLP_SOL)
1.223 + { if (npp->r_stat[info->p] == GLP_BS)
1.224 + { if (npp->c_stat[info->s] == GLP_BS)
1.225 + { npp_error();
1.226 + return 1;
1.227 + }
1.228 + else if (npp->c_stat[info->s] == GLP_NL ||
1.229 + npp->c_stat[info->s] == GLP_NU)
1.230 + npp->r_stat[info->p] = GLP_BS;
1.231 + else
1.232 + { npp_error();
1.233 + return 1;
1.234 + }
1.235 + }
1.236 + else if (npp->r_stat[info->p] == GLP_NS)
1.237 + { if (npp->c_stat[info->s] == GLP_BS)
1.238 + npp->r_stat[info->p] = GLP_BS;
1.239 + else if (npp->c_stat[info->s] == GLP_NL)
1.240 + npp->r_stat[info->p] = GLP_NL;
1.241 + else if (npp->c_stat[info->s] == GLP_NU)
1.242 + npp->r_stat[info->p] = GLP_NU;
1.243 + else
1.244 + { npp_error();
1.245 + return 1;
1.246 + }
1.247 + }
1.248 + else
1.249 + { npp_error();
1.250 + return 1;
1.251 + }
1.252 + }
1.253 + return 0;
1.254 +}
1.255 +
1.256 +/***********************************************************************
1.257 +* NAME
1.258 +*
1.259 +* npp_leq_row - process row of 'not greater than' type
1.260 +*
1.261 +* SYNOPSIS
1.262 +*
1.263 +* #include "glpnpp.h"
1.264 +* void npp_leq_row(NPP *npp, NPPROW *p);
1.265 +*
1.266 +* DESCRIPTION
1.267 +*
1.268 +* The routine npp_leq_row processes row p, which is 'not greater than'
1.269 +* inequality constraint:
1.270 +*
1.271 +* (L[p] <=) sum a[p,j] x[j] <= U[p], (1)
1.272 +* j
1.273 +*
1.274 +* where L[p] < U[p], and lower bound may not exist (L[p] = +oo).
1.275 +*
1.276 +* PROBLEM TRANSFORMATION
1.277 +*
1.278 +* Constraint (1) can be replaced by equality constraint:
1.279 +*
1.280 +* sum a[p,j] x[j] + s = L[p], (2)
1.281 +* j
1.282 +*
1.283 +* where
1.284 +*
1.285 +* 0 <= s (<= U[p] - L[p]) (3)
1.286 +*
1.287 +* is a non-negative slack variable.
1.288 +*
1.289 +* Since in the primal system there appears column s having the only
1.290 +* non-zero coefficient in row p, in the dual system there appears a
1.291 +* new row:
1.292 +*
1.293 +* (+1) pi[p] + lambda = 0, (4)
1.294 +*
1.295 +* where (+1) is coefficient of column s in row p, pi[p] is multiplier
1.296 +* of row p, lambda is multiplier of column q, 0 is coefficient of
1.297 +* column s in the objective row.
1.298 +*
1.299 +* RECOVERING BASIC SOLUTION
1.300 +*
1.301 +* Status of row p in solution to the original problem is determined
1.302 +* by its status and status of column q in solution to the transformed
1.303 +* problem as follows:
1.304 +*
1.305 +* +--------------------------------------+------------------+
1.306 +* | Transformed problem | Original problem |
1.307 +* +-----------------+--------------------+------------------+
1.308 +* | Status of row p | Status of column s | Status of row p |
1.309 +* +-----------------+--------------------+------------------+
1.310 +* | GLP_BS | GLP_BS | N/A |
1.311 +* | GLP_BS | GLP_NL | GLP_BS |
1.312 +* | GLP_BS | GLP_NU | GLP_BS |
1.313 +* | GLP_NS | GLP_BS | GLP_BS |
1.314 +* | GLP_NS | GLP_NL | GLP_NU |
1.315 +* | GLP_NS | GLP_NU | GLP_NL |
1.316 +* +-----------------+--------------------+------------------+
1.317 +*
1.318 +* Value of row multiplier pi[p] in solution to the original problem
1.319 +* is the same as in solution to the transformed problem.
1.320 +*
1.321 +* 1. In solution to the transformed problem row p and column q cannot
1.322 +* be basic at the same time; otherwise the basis matrix would have
1.323 +* two linear dependent columns: unity column of auxiliary variable
1.324 +* of row p and unity column of variable s.
1.325 +*
1.326 +* 2. Though in the transformed problem row p is equality constraint,
1.327 +* it may be basic due to primal degeneracy.
1.328 +*
1.329 +* RECOVERING INTERIOR-POINT SOLUTION
1.330 +*
1.331 +* Value of row multiplier pi[p] in solution to the original problem
1.332 +* is the same as in solution to the transformed problem.
1.333 +*
1.334 +* RECOVERING MIP SOLUTION
1.335 +*
1.336 +* None needed. */
1.337 +
1.338 +static int rcv_leq_row(NPP *npp, void *info);
1.339 +
1.340 +void npp_leq_row(NPP *npp, NPPROW *p)
1.341 +{ /* process row of 'not greater than' type */
1.342 + struct ineq_row *info;
1.343 + NPPCOL *s;
1.344 + /* the row must have upper bound */
1.345 + xassert(p->ub != +DBL_MAX);
1.346 + xassert(p->lb < p->ub);
1.347 + /* create column for slack variable */
1.348 + s = npp_add_col(npp);
1.349 + s->lb = 0.0;
1.350 + s->ub = (p->lb == -DBL_MAX ? +DBL_MAX : p->ub - p->lb);
1.351 + /* and add it to the transformed problem */
1.352 + npp_add_aij(npp, p, s, +1.0);
1.353 + /* create transformation stack entry */
1.354 + info = npp_push_tse(npp,
1.355 + rcv_leq_row, sizeof(struct ineq_row));
1.356 + info->p = p->i;
1.357 + info->s = s->j;
1.358 + /* replace the row by equality constraint */
1.359 + p->lb = p->ub;
1.360 + return;
1.361 +}
1.362 +
1.363 +static int rcv_leq_row(NPP *npp, void *_info)
1.364 +{ /* recover row of 'not greater than' type */
1.365 + struct ineq_row *info = _info;
1.366 + if (npp->sol == GLP_SOL)
1.367 + { if (npp->r_stat[info->p] == GLP_BS)
1.368 + { if (npp->c_stat[info->s] == GLP_BS)
1.369 + { npp_error();
1.370 + return 1;
1.371 + }
1.372 + else if (npp->c_stat[info->s] == GLP_NL ||
1.373 + npp->c_stat[info->s] == GLP_NU)
1.374 + npp->r_stat[info->p] = GLP_BS;
1.375 + else
1.376 + { npp_error();
1.377 + return 1;
1.378 + }
1.379 + }
1.380 + else if (npp->r_stat[info->p] == GLP_NS)
1.381 + { if (npp->c_stat[info->s] == GLP_BS)
1.382 + npp->r_stat[info->p] = GLP_BS;
1.383 + else if (npp->c_stat[info->s] == GLP_NL)
1.384 + npp->r_stat[info->p] = GLP_NU;
1.385 + else if (npp->c_stat[info->s] == GLP_NU)
1.386 + npp->r_stat[info->p] = GLP_NL;
1.387 + else
1.388 + { npp_error();
1.389 + return 1;
1.390 + }
1.391 + }
1.392 + else
1.393 + { npp_error();
1.394 + return 1;
1.395 + }
1.396 + }
1.397 + return 0;
1.398 +}
1.399 +
1.400 +/***********************************************************************
1.401 +* NAME
1.402 +*
1.403 +* npp_free_col - process free (unbounded) column
1.404 +*
1.405 +* SYNOPSIS
1.406 +*
1.407 +* #include "glpnpp.h"
1.408 +* void npp_free_col(NPP *npp, NPPCOL *q);
1.409 +*
1.410 +* DESCRIPTION
1.411 +*
1.412 +* The routine npp_free_col processes column q, which is free (i.e. has
1.413 +* no finite bounds):
1.414 +*
1.415 +* -oo < x[q] < +oo. (1)
1.416 +*
1.417 +* PROBLEM TRANSFORMATION
1.418 +*
1.419 +* Free (unbounded) variable can be replaced by the difference of two
1.420 +* non-negative variables:
1.421 +*
1.422 +* x[q] = s' - s'', s', s'' >= 0. (2)
1.423 +*
1.424 +* Assuming that in the transformed problem x[q] becomes s',
1.425 +* transformation (2) causes new column s'' to appear, which differs
1.426 +* from column s' only in the sign of coefficients in constraint and
1.427 +* objective rows. Thus, if in the dual system the following row
1.428 +* corresponds to column s':
1.429 +*
1.430 +* sum a[i,q] pi[i] + lambda' = c[q], (3)
1.431 +* i
1.432 +*
1.433 +* the row which corresponds to column s'' is the following:
1.434 +*
1.435 +* sum (-a[i,q]) pi[i] + lambda'' = -c[q]. (4)
1.436 +* i
1.437 +*
1.438 +* Then from (3) and (4) it follows that:
1.439 +*
1.440 +* lambda' + lambda'' = 0 => lambda' = lmabda'' = 0, (5)
1.441 +*
1.442 +* where lambda' and lambda'' are multipliers for columns s' and s'',
1.443 +* resp.
1.444 +*
1.445 +* RECOVERING BASIC SOLUTION
1.446 +*
1.447 +* With respect to (5) status of column q in solution to the original
1.448 +* problem is determined by statuses of columns s' and s'' in solution
1.449 +* to the transformed problem as follows:
1.450 +*
1.451 +* +--------------------------------------+------------------+
1.452 +* | Transformed problem | Original problem |
1.453 +* +------------------+-------------------+------------------+
1.454 +* | Status of col s' | Status of col s'' | Status of col q |
1.455 +* +------------------+-------------------+------------------+
1.456 +* | GLP_BS | GLP_BS | N/A |
1.457 +* | GLP_BS | GLP_NL | GLP_BS |
1.458 +* | GLP_NL | GLP_BS | GLP_BS |
1.459 +* | GLP_NL | GLP_NL | GLP_NF |
1.460 +* +------------------+-------------------+------------------+
1.461 +*
1.462 +* Value of column q is computed with formula (2).
1.463 +*
1.464 +* 1. In solution to the transformed problem columns s' and s'' cannot
1.465 +* be basic at the same time, because they differ only in the sign,
1.466 +* hence, are linear dependent.
1.467 +*
1.468 +* 2. Though column q is free, it can be non-basic due to dual
1.469 +* degeneracy.
1.470 +*
1.471 +* 3. If column q is integral, columns s' and s'' are also integral.
1.472 +*
1.473 +* RECOVERING INTERIOR-POINT SOLUTION
1.474 +*
1.475 +* Value of column q is computed with formula (2).
1.476 +*
1.477 +* RECOVERING MIP SOLUTION
1.478 +*
1.479 +* Value of column q is computed with formula (2). */
1.480 +
1.481 +struct free_col
1.482 +{ /* free (unbounded) column */
1.483 + int q;
1.484 + /* column reference number for variables x[q] and s' */
1.485 + int s;
1.486 + /* column reference number for variable s'' */
1.487 +};
1.488 +
1.489 +static int rcv_free_col(NPP *npp, void *info);
1.490 +
1.491 +void npp_free_col(NPP *npp, NPPCOL *q)
1.492 +{ /* process free (unbounded) column */
1.493 + struct free_col *info;
1.494 + NPPCOL *s;
1.495 + NPPAIJ *aij;
1.496 + /* the column must be free */
1.497 + xassert(q->lb == -DBL_MAX && q->ub == +DBL_MAX);
1.498 + /* variable x[q] becomes s' */
1.499 + q->lb = 0.0, q->ub = +DBL_MAX;
1.500 + /* create variable s'' */
1.501 + s = npp_add_col(npp);
1.502 + s->is_int = q->is_int;
1.503 + s->lb = 0.0, s->ub = +DBL_MAX;
1.504 + /* duplicate objective coefficient */
1.505 + s->coef = -q->coef;
1.506 + /* duplicate column of the constraint matrix */
1.507 + for (aij = q->ptr; aij != NULL; aij = aij->c_next)
1.508 + npp_add_aij(npp, aij->row, s, -aij->val);
1.509 + /* create transformation stack entry */
1.510 + info = npp_push_tse(npp,
1.511 + rcv_free_col, sizeof(struct free_col));
1.512 + info->q = q->j;
1.513 + info->s = s->j;
1.514 + return;
1.515 +}
1.516 +
1.517 +static int rcv_free_col(NPP *npp, void *_info)
1.518 +{ /* recover free (unbounded) column */
1.519 + struct free_col *info = _info;
1.520 + if (npp->sol == GLP_SOL)
1.521 + { if (npp->c_stat[info->q] == GLP_BS)
1.522 + { if (npp->c_stat[info->s] == GLP_BS)
1.523 + { npp_error();
1.524 + return 1;
1.525 + }
1.526 + else if (npp->c_stat[info->s] == GLP_NL)
1.527 + npp->c_stat[info->q] = GLP_BS;
1.528 + else
1.529 + { npp_error();
1.530 + return -1;
1.531 + }
1.532 + }
1.533 + else if (npp->c_stat[info->q] == GLP_NL)
1.534 + { if (npp->c_stat[info->s] == GLP_BS)
1.535 + npp->c_stat[info->q] = GLP_BS;
1.536 + else if (npp->c_stat[info->s] == GLP_NL)
1.537 + npp->c_stat[info->q] = GLP_NF;
1.538 + else
1.539 + { npp_error();
1.540 + return -1;
1.541 + }
1.542 + }
1.543 + else
1.544 + { npp_error();
1.545 + return -1;
1.546 + }
1.547 + }
1.548 + /* compute value of x[q] with formula (2) */
1.549 + npp->c_value[info->q] -= npp->c_value[info->s];
1.550 + return 0;
1.551 +}
1.552 +
1.553 +/***********************************************************************
1.554 +* NAME
1.555 +*
1.556 +* npp_lbnd_col - process column with (non-zero) lower bound
1.557 +*
1.558 +* SYNOPSIS
1.559 +*
1.560 +* #include "glpnpp.h"
1.561 +* void npp_lbnd_col(NPP *npp, NPPCOL *q);
1.562 +*
1.563 +* DESCRIPTION
1.564 +*
1.565 +* The routine npp_lbnd_col processes column q, which has (non-zero)
1.566 +* lower bound:
1.567 +*
1.568 +* l[q] <= x[q] (<= u[q]), (1)
1.569 +*
1.570 +* where l[q] < u[q], and upper bound may not exist (u[q] = +oo).
1.571 +*
1.572 +* PROBLEM TRANSFORMATION
1.573 +*
1.574 +* Column q can be replaced as follows:
1.575 +*
1.576 +* x[q] = l[q] + s, (2)
1.577 +*
1.578 +* where
1.579 +*
1.580 +* 0 <= s (<= u[q] - l[q]) (3)
1.581 +*
1.582 +* is a non-negative variable.
1.583 +*
1.584 +* Substituting x[q] from (2) into the objective row, we have:
1.585 +*
1.586 +* z = sum c[j] x[j] + c0 =
1.587 +* j
1.588 +*
1.589 +* = sum c[j] x[j] + c[q] x[q] + c0 =
1.590 +* j!=q
1.591 +*
1.592 +* = sum c[j] x[j] + c[q] (l[q] + s) + c0 =
1.593 +* j!=q
1.594 +*
1.595 +* = sum c[j] x[j] + c[q] s + c~0,
1.596 +*
1.597 +* where
1.598 +*
1.599 +* c~0 = c0 + c[q] l[q] (4)
1.600 +*
1.601 +* is the constant term of the objective in the transformed problem.
1.602 +* Similarly, substituting x[q] into constraint row i, we have:
1.603 +*
1.604 +* L[i] <= sum a[i,j] x[j] <= U[i] ==>
1.605 +* j
1.606 +*
1.607 +* L[i] <= sum a[i,j] x[j] + a[i,q] x[q] <= U[i] ==>
1.608 +* j!=q
1.609 +*
1.610 +* L[i] <= sum a[i,j] x[j] + a[i,q] (l[q] + s) <= U[i] ==>
1.611 +* j!=q
1.612 +*
1.613 +* L~[i] <= sum a[i,j] x[j] + a[i,q] s <= U~[i],
1.614 +* j!=q
1.615 +*
1.616 +* where
1.617 +*
1.618 +* L~[i] = L[i] - a[i,q] l[q], U~[i] = U[i] - a[i,q] l[q] (5)
1.619 +*
1.620 +* are lower and upper bounds of row i in the transformed problem,
1.621 +* resp.
1.622 +*
1.623 +* Transformation (2) does not affect the dual system.
1.624 +*
1.625 +* RECOVERING BASIC SOLUTION
1.626 +*
1.627 +* Status of column q in solution to the original problem is the same
1.628 +* as in solution to the transformed problem (GLP_BS, GLP_NL or GLP_NU).
1.629 +* Value of column q is computed with formula (2).
1.630 +*
1.631 +* RECOVERING INTERIOR-POINT SOLUTION
1.632 +*
1.633 +* Value of column q is computed with formula (2).
1.634 +*
1.635 +* RECOVERING MIP SOLUTION
1.636 +*
1.637 +* Value of column q is computed with formula (2). */
1.638 +
1.639 +struct bnd_col
1.640 +{ /* bounded column */
1.641 + int q;
1.642 + /* column reference number for variables x[q] and s */
1.643 + double bnd;
1.644 + /* lower/upper bound l[q] or u[q] */
1.645 +};
1.646 +
1.647 +static int rcv_lbnd_col(NPP *npp, void *info);
1.648 +
1.649 +void npp_lbnd_col(NPP *npp, NPPCOL *q)
1.650 +{ /* process column with (non-zero) lower bound */
1.651 + struct bnd_col *info;
1.652 + NPPROW *i;
1.653 + NPPAIJ *aij;
1.654 + /* the column must have non-zero lower bound */
1.655 + xassert(q->lb != 0.0);
1.656 + xassert(q->lb != -DBL_MAX);
1.657 + xassert(q->lb < q->ub);
1.658 + /* create transformation stack entry */
1.659 + info = npp_push_tse(npp,
1.660 + rcv_lbnd_col, sizeof(struct bnd_col));
1.661 + info->q = q->j;
1.662 + info->bnd = q->lb;
1.663 + /* substitute x[q] into objective row */
1.664 + npp->c0 += q->coef * q->lb;
1.665 + /* substitute x[q] into constraint rows */
1.666 + for (aij = q->ptr; aij != NULL; aij = aij->c_next)
1.667 + { i = aij->row;
1.668 + if (i->lb == i->ub)
1.669 + i->ub = (i->lb -= aij->val * q->lb);
1.670 + else
1.671 + { if (i->lb != -DBL_MAX)
1.672 + i->lb -= aij->val * q->lb;
1.673 + if (i->ub != +DBL_MAX)
1.674 + i->ub -= aij->val * q->lb;
1.675 + }
1.676 + }
1.677 + /* column x[q] becomes column s */
1.678 + if (q->ub != +DBL_MAX)
1.679 + q->ub -= q->lb;
1.680 + q->lb = 0.0;
1.681 + return;
1.682 +}
1.683 +
1.684 +static int rcv_lbnd_col(NPP *npp, void *_info)
1.685 +{ /* recover column with (non-zero) lower bound */
1.686 + struct bnd_col *info = _info;
1.687 + if (npp->sol == GLP_SOL)
1.688 + { if (npp->c_stat[info->q] == GLP_BS ||
1.689 + npp->c_stat[info->q] == GLP_NL ||
1.690 + npp->c_stat[info->q] == GLP_NU)
1.691 + npp->c_stat[info->q] = npp->c_stat[info->q];
1.692 + else
1.693 + { npp_error();
1.694 + return 1;
1.695 + }
1.696 + }
1.697 + /* compute value of x[q] with formula (2) */
1.698 + npp->c_value[info->q] = info->bnd + npp->c_value[info->q];
1.699 + return 0;
1.700 +}
1.701 +
1.702 +/***********************************************************************
1.703 +* NAME
1.704 +*
1.705 +* npp_ubnd_col - process column with upper bound
1.706 +*
1.707 +* SYNOPSIS
1.708 +*
1.709 +* #include "glpnpp.h"
1.710 +* void npp_ubnd_col(NPP *npp, NPPCOL *q);
1.711 +*
1.712 +* DESCRIPTION
1.713 +*
1.714 +* The routine npp_ubnd_col processes column q, which has upper bound:
1.715 +*
1.716 +* (l[q] <=) x[q] <= u[q], (1)
1.717 +*
1.718 +* where l[q] < u[q], and lower bound may not exist (l[q] = -oo).
1.719 +*
1.720 +* PROBLEM TRANSFORMATION
1.721 +*
1.722 +* Column q can be replaced as follows:
1.723 +*
1.724 +* x[q] = u[q] - s, (2)
1.725 +*
1.726 +* where
1.727 +*
1.728 +* 0 <= s (<= u[q] - l[q]) (3)
1.729 +*
1.730 +* is a non-negative variable.
1.731 +*
1.732 +* Substituting x[q] from (2) into the objective row, we have:
1.733 +*
1.734 +* z = sum c[j] x[j] + c0 =
1.735 +* j
1.736 +*
1.737 +* = sum c[j] x[j] + c[q] x[q] + c0 =
1.738 +* j!=q
1.739 +*
1.740 +* = sum c[j] x[j] + c[q] (u[q] - s) + c0 =
1.741 +* j!=q
1.742 +*
1.743 +* = sum c[j] x[j] - c[q] s + c~0,
1.744 +*
1.745 +* where
1.746 +*
1.747 +* c~0 = c0 + c[q] u[q] (4)
1.748 +*
1.749 +* is the constant term of the objective in the transformed problem.
1.750 +* Similarly, substituting x[q] into constraint row i, we have:
1.751 +*
1.752 +* L[i] <= sum a[i,j] x[j] <= U[i] ==>
1.753 +* j
1.754 +*
1.755 +* L[i] <= sum a[i,j] x[j] + a[i,q] x[q] <= U[i] ==>
1.756 +* j!=q
1.757 +*
1.758 +* L[i] <= sum a[i,j] x[j] + a[i,q] (u[q] - s) <= U[i] ==>
1.759 +* j!=q
1.760 +*
1.761 +* L~[i] <= sum a[i,j] x[j] - a[i,q] s <= U~[i],
1.762 +* j!=q
1.763 +*
1.764 +* where
1.765 +*
1.766 +* L~[i] = L[i] - a[i,q] u[q], U~[i] = U[i] - a[i,q] u[q] (5)
1.767 +*
1.768 +* are lower and upper bounds of row i in the transformed problem,
1.769 +* resp.
1.770 +*
1.771 +* Note that in the transformed problem coefficients c[q] and a[i,q]
1.772 +* change their sign. Thus, the row of the dual system corresponding to
1.773 +* column q:
1.774 +*
1.775 +* sum a[i,q] pi[i] + lambda[q] = c[q] (6)
1.776 +* i
1.777 +*
1.778 +* in the transformed problem becomes the following:
1.779 +*
1.780 +* sum (-a[i,q]) pi[i] + lambda[s] = -c[q]. (7)
1.781 +* i
1.782 +*
1.783 +* Therefore:
1.784 +*
1.785 +* lambda[q] = - lambda[s], (8)
1.786 +*
1.787 +* where lambda[q] is multiplier for column q, lambda[s] is multiplier
1.788 +* for column s.
1.789 +*
1.790 +* RECOVERING BASIC SOLUTION
1.791 +*
1.792 +* With respect to (8) status of column q in solution to the original
1.793 +* problem is determined by status of column s in solution to the
1.794 +* transformed problem as follows:
1.795 +*
1.796 +* +-----------------------+--------------------+
1.797 +* | Status of column s | Status of column q |
1.798 +* | (transformed problem) | (original problem) |
1.799 +* +-----------------------+--------------------+
1.800 +* | GLP_BS | GLP_BS |
1.801 +* | GLP_NL | GLP_NU |
1.802 +* | GLP_NU | GLP_NL |
1.803 +* +-----------------------+--------------------+
1.804 +*
1.805 +* Value of column q is computed with formula (2).
1.806 +*
1.807 +* RECOVERING INTERIOR-POINT SOLUTION
1.808 +*
1.809 +* Value of column q is computed with formula (2).
1.810 +*
1.811 +* RECOVERING MIP SOLUTION
1.812 +*
1.813 +* Value of column q is computed with formula (2). */
1.814 +
1.815 +static int rcv_ubnd_col(NPP *npp, void *info);
1.816 +
1.817 +void npp_ubnd_col(NPP *npp, NPPCOL *q)
1.818 +{ /* process column with upper bound */
1.819 + struct bnd_col *info;
1.820 + NPPROW *i;
1.821 + NPPAIJ *aij;
1.822 + /* the column must have upper bound */
1.823 + xassert(q->ub != +DBL_MAX);
1.824 + xassert(q->lb < q->ub);
1.825 + /* create transformation stack entry */
1.826 + info = npp_push_tse(npp,
1.827 + rcv_ubnd_col, sizeof(struct bnd_col));
1.828 + info->q = q->j;
1.829 + info->bnd = q->ub;
1.830 + /* substitute x[q] into objective row */
1.831 + npp->c0 += q->coef * q->ub;
1.832 + q->coef = -q->coef;
1.833 + /* substitute x[q] into constraint rows */
1.834 + for (aij = q->ptr; aij != NULL; aij = aij->c_next)
1.835 + { i = aij->row;
1.836 + if (i->lb == i->ub)
1.837 + i->ub = (i->lb -= aij->val * q->ub);
1.838 + else
1.839 + { if (i->lb != -DBL_MAX)
1.840 + i->lb -= aij->val * q->ub;
1.841 + if (i->ub != +DBL_MAX)
1.842 + i->ub -= aij->val * q->ub;
1.843 + }
1.844 + aij->val = -aij->val;
1.845 + }
1.846 + /* column x[q] becomes column s */
1.847 + if (q->lb != -DBL_MAX)
1.848 + q->ub -= q->lb;
1.849 + else
1.850 + q->ub = +DBL_MAX;
1.851 + q->lb = 0.0;
1.852 + return;
1.853 +}
1.854 +
1.855 +static int rcv_ubnd_col(NPP *npp, void *_info)
1.856 +{ /* recover column with upper bound */
1.857 + struct bnd_col *info = _info;
1.858 + if (npp->sol == GLP_BS)
1.859 + { if (npp->c_stat[info->q] == GLP_BS)
1.860 + npp->c_stat[info->q] = GLP_BS;
1.861 + else if (npp->c_stat[info->q] == GLP_NL)
1.862 + npp->c_stat[info->q] = GLP_NU;
1.863 + else if (npp->c_stat[info->q] == GLP_NU)
1.864 + npp->c_stat[info->q] = GLP_NL;
1.865 + else
1.866 + { npp_error();
1.867 + return 1;
1.868 + }
1.869 + }
1.870 + /* compute value of x[q] with formula (2) */
1.871 + npp->c_value[info->q] = info->bnd - npp->c_value[info->q];
1.872 + return 0;
1.873 +}
1.874 +
1.875 +/***********************************************************************
1.876 +* NAME
1.877 +*
1.878 +* npp_dbnd_col - process non-negative column with upper bound
1.879 +*
1.880 +* SYNOPSIS
1.881 +*
1.882 +* #include "glpnpp.h"
1.883 +* void npp_dbnd_col(NPP *npp, NPPCOL *q);
1.884 +*
1.885 +* DESCRIPTION
1.886 +*
1.887 +* The routine npp_dbnd_col processes column q, which is non-negative
1.888 +* and has upper bound:
1.889 +*
1.890 +* 0 <= x[q] <= u[q], (1)
1.891 +*
1.892 +* where u[q] > 0.
1.893 +*
1.894 +* PROBLEM TRANSFORMATION
1.895 +*
1.896 +* Upper bound of column q can be replaced by the following equality
1.897 +* constraint:
1.898 +*
1.899 +* x[q] + s = u[q], (2)
1.900 +*
1.901 +* where s >= 0 is a non-negative complement variable.
1.902 +*
1.903 +* Since in the primal system along with new row (2) there appears a
1.904 +* new column s having the only non-zero coefficient in this row, in
1.905 +* the dual system there appears a new row:
1.906 +*
1.907 +* (+1)pi + lambda[s] = 0, (3)
1.908 +*
1.909 +* where (+1) is coefficient at column s in row (2), pi is multiplier
1.910 +* for row (2), lambda[s] is multiplier for column s, 0 is coefficient
1.911 +* at column s in the objective row.
1.912 +*
1.913 +* RECOVERING BASIC SOLUTION
1.914 +*
1.915 +* Status of column q in solution to the original problem is determined
1.916 +* by its status and status of column s in solution to the transformed
1.917 +* problem as follows:
1.918 +*
1.919 +* +-----------------------------------+------------------+
1.920 +* | Transformed problem | Original problem |
1.921 +* +-----------------+-----------------+------------------+
1.922 +* | Status of col q | Status of col s | Status of col q |
1.923 +* +-----------------+-----------------+------------------+
1.924 +* | GLP_BS | GLP_BS | GLP_BS |
1.925 +* | GLP_BS | GLP_NL | GLP_NU |
1.926 +* | GLP_NL | GLP_BS | GLP_NL |
1.927 +* | GLP_NL | GLP_NL | GLP_NL (*) |
1.928 +* +-----------------+-----------------+------------------+
1.929 +*
1.930 +* Value of column q in solution to the original problem is the same as
1.931 +* in solution to the transformed problem.
1.932 +*
1.933 +* 1. Formally, in solution to the transformed problem columns q and s
1.934 +* cannot be non-basic at the same time, since the constraint (2)
1.935 +* would be violated. However, if u[q] is close to zero, violation
1.936 +* may be less than a working precision even if both columns q and s
1.937 +* are non-basic. In this degenerate case row (2) can be only basic,
1.938 +* i.e. non-active constraint (otherwise corresponding row of the
1.939 +* basis matrix would be zero). This allows to pivot out auxiliary
1.940 +* variable and pivot in column s, in which case the row becomes
1.941 +* active while column s becomes basic.
1.942 +*
1.943 +* 2. If column q is integral, column s is also integral.
1.944 +*
1.945 +* RECOVERING INTERIOR-POINT SOLUTION
1.946 +*
1.947 +* Value of column q in solution to the original problem is the same as
1.948 +* in solution to the transformed problem.
1.949 +*
1.950 +* RECOVERING MIP SOLUTION
1.951 +*
1.952 +* Value of column q in solution to the original problem is the same as
1.953 +* in solution to the transformed problem. */
1.954 +
1.955 +struct dbnd_col
1.956 +{ /* double-bounded column */
1.957 + int q;
1.958 + /* column reference number for variable x[q] */
1.959 + int s;
1.960 + /* column reference number for complement variable s */
1.961 +};
1.962 +
1.963 +static int rcv_dbnd_col(NPP *npp, void *info);
1.964 +
1.965 +void npp_dbnd_col(NPP *npp, NPPCOL *q)
1.966 +{ /* process non-negative column with upper bound */
1.967 + struct dbnd_col *info;
1.968 + NPPROW *p;
1.969 + NPPCOL *s;
1.970 + /* the column must be non-negative with upper bound */
1.971 + xassert(q->lb == 0.0);
1.972 + xassert(q->ub > 0.0);
1.973 + xassert(q->ub != +DBL_MAX);
1.974 + /* create variable s */
1.975 + s = npp_add_col(npp);
1.976 + s->is_int = q->is_int;
1.977 + s->lb = 0.0, s->ub = +DBL_MAX;
1.978 + /* create equality constraint (2) */
1.979 + p = npp_add_row(npp);
1.980 + p->lb = p->ub = q->ub;
1.981 + npp_add_aij(npp, p, q, +1.0);
1.982 + npp_add_aij(npp, p, s, +1.0);
1.983 + /* create transformation stack entry */
1.984 + info = npp_push_tse(npp,
1.985 + rcv_dbnd_col, sizeof(struct dbnd_col));
1.986 + info->q = q->j;
1.987 + info->s = s->j;
1.988 + /* remove upper bound of x[q] */
1.989 + q->ub = +DBL_MAX;
1.990 + return;
1.991 +}
1.992 +
1.993 +static int rcv_dbnd_col(NPP *npp, void *_info)
1.994 +{ /* recover non-negative column with upper bound */
1.995 + struct dbnd_col *info = _info;
1.996 + if (npp->sol == GLP_BS)
1.997 + { if (npp->c_stat[info->q] == GLP_BS)
1.998 + { if (npp->c_stat[info->s] == GLP_BS)
1.999 + npp->c_stat[info->q] = GLP_BS;
1.1000 + else if (npp->c_stat[info->s] == GLP_NL)
1.1001 + npp->c_stat[info->q] = GLP_NU;
1.1002 + else
1.1003 + { npp_error();
1.1004 + return 1;
1.1005 + }
1.1006 + }
1.1007 + else if (npp->c_stat[info->q] == GLP_NL)
1.1008 + { if (npp->c_stat[info->s] == GLP_BS ||
1.1009 + npp->c_stat[info->s] == GLP_NL)
1.1010 + npp->c_stat[info->q] = GLP_NL;
1.1011 + else
1.1012 + { npp_error();
1.1013 + return 1;
1.1014 + }
1.1015 + }
1.1016 + else
1.1017 + { npp_error();
1.1018 + return 1;
1.1019 + }
1.1020 + }
1.1021 + return 0;
1.1022 +}
1.1023 +
1.1024 +/***********************************************************************
1.1025 +* NAME
1.1026 +*
1.1027 +* npp_fixed_col - process fixed column
1.1028 +*
1.1029 +* SYNOPSIS
1.1030 +*
1.1031 +* #include "glpnpp.h"
1.1032 +* void npp_fixed_col(NPP *npp, NPPCOL *q);
1.1033 +*
1.1034 +* DESCRIPTION
1.1035 +*
1.1036 +* The routine npp_fixed_col processes column q, which is fixed:
1.1037 +*
1.1038 +* x[q] = s[q], (1)
1.1039 +*
1.1040 +* where s[q] is a fixed column value.
1.1041 +*
1.1042 +* PROBLEM TRANSFORMATION
1.1043 +*
1.1044 +* The value of a fixed column can be substituted into the objective
1.1045 +* and constraint rows that allows removing the column from the problem.
1.1046 +*
1.1047 +* Substituting x[q] = s[q] into the objective row, we have:
1.1048 +*
1.1049 +* z = sum c[j] x[j] + c0 =
1.1050 +* j
1.1051 +*
1.1052 +* = sum c[j] x[j] + c[q] x[q] + c0 =
1.1053 +* j!=q
1.1054 +*
1.1055 +* = sum c[j] x[j] + c[q] s[q] + c0 =
1.1056 +* j!=q
1.1057 +*
1.1058 +* = sum c[j] x[j] + c~0,
1.1059 +* j!=q
1.1060 +*
1.1061 +* where
1.1062 +*
1.1063 +* c~0 = c0 + c[q] s[q] (2)
1.1064 +*
1.1065 +* is the constant term of the objective in the transformed problem.
1.1066 +* Similarly, substituting x[q] = s[q] into constraint row i, we have:
1.1067 +*
1.1068 +* L[i] <= sum a[i,j] x[j] <= U[i] ==>
1.1069 +* j
1.1070 +*
1.1071 +* L[i] <= sum a[i,j] x[j] + a[i,q] x[q] <= U[i] ==>
1.1072 +* j!=q
1.1073 +*
1.1074 +* L[i] <= sum a[i,j] x[j] + a[i,q] s[q] <= U[i] ==>
1.1075 +* j!=q
1.1076 +*
1.1077 +* L~[i] <= sum a[i,j] x[j] + a[i,q] s <= U~[i],
1.1078 +* j!=q
1.1079 +*
1.1080 +* where
1.1081 +*
1.1082 +* L~[i] = L[i] - a[i,q] s[q], U~[i] = U[i] - a[i,q] s[q] (3)
1.1083 +*
1.1084 +* are lower and upper bounds of row i in the transformed problem,
1.1085 +* resp.
1.1086 +*
1.1087 +* RECOVERING BASIC SOLUTION
1.1088 +*
1.1089 +* Column q is assigned status GLP_NS and its value is assigned s[q].
1.1090 +*
1.1091 +* RECOVERING INTERIOR-POINT SOLUTION
1.1092 +*
1.1093 +* Value of column q is assigned s[q].
1.1094 +*
1.1095 +* RECOVERING MIP SOLUTION
1.1096 +*
1.1097 +* Value of column q is assigned s[q]. */
1.1098 +
1.1099 +struct fixed_col
1.1100 +{ /* fixed column */
1.1101 + int q;
1.1102 + /* column reference number for variable x[q] */
1.1103 + double s;
1.1104 + /* value, at which x[q] is fixed */
1.1105 +};
1.1106 +
1.1107 +static int rcv_fixed_col(NPP *npp, void *info);
1.1108 +
1.1109 +void npp_fixed_col(NPP *npp, NPPCOL *q)
1.1110 +{ /* process fixed column */
1.1111 + struct fixed_col *info;
1.1112 + NPPROW *i;
1.1113 + NPPAIJ *aij;
1.1114 + /* the column must be fixed */
1.1115 + xassert(q->lb == q->ub);
1.1116 + /* create transformation stack entry */
1.1117 + info = npp_push_tse(npp,
1.1118 + rcv_fixed_col, sizeof(struct fixed_col));
1.1119 + info->q = q->j;
1.1120 + info->s = q->lb;
1.1121 + /* substitute x[q] = s[q] into objective row */
1.1122 + npp->c0 += q->coef * q->lb;
1.1123 + /* substitute x[q] = s[q] into constraint rows */
1.1124 + for (aij = q->ptr; aij != NULL; aij = aij->c_next)
1.1125 + { i = aij->row;
1.1126 + if (i->lb == i->ub)
1.1127 + i->ub = (i->lb -= aij->val * q->lb);
1.1128 + else
1.1129 + { if (i->lb != -DBL_MAX)
1.1130 + i->lb -= aij->val * q->lb;
1.1131 + if (i->ub != +DBL_MAX)
1.1132 + i->ub -= aij->val * q->lb;
1.1133 + }
1.1134 + }
1.1135 + /* remove the column from the problem */
1.1136 + npp_del_col(npp, q);
1.1137 + return;
1.1138 +}
1.1139 +
1.1140 +static int rcv_fixed_col(NPP *npp, void *_info)
1.1141 +{ /* recover fixed column */
1.1142 + struct fixed_col *info = _info;
1.1143 + if (npp->sol == GLP_SOL)
1.1144 + npp->c_stat[info->q] = GLP_NS;
1.1145 + npp->c_value[info->q] = info->s;
1.1146 + return 0;
1.1147 +}
1.1148 +
1.1149 +/***********************************************************************
1.1150 +* NAME
1.1151 +*
1.1152 +* npp_make_equality - process row with almost identical bounds
1.1153 +*
1.1154 +* SYNOPSIS
1.1155 +*
1.1156 +* #include "glpnpp.h"
1.1157 +* int npp_make_equality(NPP *npp, NPPROW *p);
1.1158 +*
1.1159 +* DESCRIPTION
1.1160 +*
1.1161 +* The routine npp_make_equality processes row p:
1.1162 +*
1.1163 +* L[p] <= sum a[p,j] x[j] <= U[p], (1)
1.1164 +* j
1.1165 +*
1.1166 +* where -oo < L[p] < U[p] < +oo, i.e. which is double-sided inequality
1.1167 +* constraint.
1.1168 +*
1.1169 +* RETURNS
1.1170 +*
1.1171 +* 0 - row bounds have not been changed;
1.1172 +*
1.1173 +* 1 - row has been replaced by equality constraint.
1.1174 +*
1.1175 +* PROBLEM TRANSFORMATION
1.1176 +*
1.1177 +* If bounds of row (1) are very close to each other:
1.1178 +*
1.1179 +* U[p] - L[p] <= eps, (2)
1.1180 +*
1.1181 +* where eps is an absolute tolerance for row value, the row can be
1.1182 +* replaced by the following almost equivalent equiality constraint:
1.1183 +*
1.1184 +* sum a[p,j] x[j] = b, (3)
1.1185 +* j
1.1186 +*
1.1187 +* where b = (L[p] + U[p]) / 2. If the right-hand side in (3) happens
1.1188 +* to be very close to its nearest integer:
1.1189 +*
1.1190 +* |b - floor(b + 0.5)| <= eps, (4)
1.1191 +*
1.1192 +* it is reasonable to use this nearest integer as the right-hand side.
1.1193 +*
1.1194 +* RECOVERING BASIC SOLUTION
1.1195 +*
1.1196 +* Status of row p in solution to the original problem is determined
1.1197 +* by its status and the sign of its multiplier pi[p] in solution to
1.1198 +* the transformed problem as follows:
1.1199 +*
1.1200 +* +-----------------------+---------+--------------------+
1.1201 +* | Status of row p | Sign of | Status of row p |
1.1202 +* | (transformed problem) | pi[p] | (original problem) |
1.1203 +* +-----------------------+---------+--------------------+
1.1204 +* | GLP_BS | + / - | GLP_BS |
1.1205 +* | GLP_NS | + | GLP_NL |
1.1206 +* | GLP_NS | - | GLP_NU |
1.1207 +* +-----------------------+---------+--------------------+
1.1208 +*
1.1209 +* Value of row multiplier pi[p] in solution to the original problem is
1.1210 +* the same as in solution to the transformed problem.
1.1211 +*
1.1212 +* RECOVERING INTERIOR POINT SOLUTION
1.1213 +*
1.1214 +* Value of row multiplier pi[p] in solution to the original problem is
1.1215 +* the same as in solution to the transformed problem.
1.1216 +*
1.1217 +* RECOVERING MIP SOLUTION
1.1218 +*
1.1219 +* None needed. */
1.1220 +
1.1221 +struct make_equality
1.1222 +{ /* row with almost identical bounds */
1.1223 + int p;
1.1224 + /* row reference number */
1.1225 +};
1.1226 +
1.1227 +static int rcv_make_equality(NPP *npp, void *info);
1.1228 +
1.1229 +int npp_make_equality(NPP *npp, NPPROW *p)
1.1230 +{ /* process row with almost identical bounds */
1.1231 + struct make_equality *info;
1.1232 + double b, eps, nint;
1.1233 + /* the row must be double-sided inequality */
1.1234 + xassert(p->lb != -DBL_MAX);
1.1235 + xassert(p->ub != +DBL_MAX);
1.1236 + xassert(p->lb < p->ub);
1.1237 + /* check row bounds */
1.1238 + eps = 1e-9 + 1e-12 * fabs(p->lb);
1.1239 + if (p->ub - p->lb > eps) return 0;
1.1240 + /* row bounds are very close to each other */
1.1241 + /* create transformation stack entry */
1.1242 + info = npp_push_tse(npp,
1.1243 + rcv_make_equality, sizeof(struct make_equality));
1.1244 + info->p = p->i;
1.1245 + /* compute right-hand side */
1.1246 + b = 0.5 * (p->ub + p->lb);
1.1247 + nint = floor(b + 0.5);
1.1248 + if (fabs(b - nint) <= eps) b = nint;
1.1249 + /* replace row p by almost equivalent equality constraint */
1.1250 + p->lb = p->ub = b;
1.1251 + return 1;
1.1252 +}
1.1253 +
1.1254 +int rcv_make_equality(NPP *npp, void *_info)
1.1255 +{ /* recover row with almost identical bounds */
1.1256 + struct make_equality *info = _info;
1.1257 + if (npp->sol == GLP_SOL)
1.1258 + { if (npp->r_stat[info->p] == GLP_BS)
1.1259 + npp->r_stat[info->p] = GLP_BS;
1.1260 + else if (npp->r_stat[info->p] == GLP_NS)
1.1261 + { if (npp->r_pi[info->p] >= 0.0)
1.1262 + npp->r_stat[info->p] = GLP_NL;
1.1263 + else
1.1264 + npp->r_stat[info->p] = GLP_NU;
1.1265 + }
1.1266 + else
1.1267 + { npp_error();
1.1268 + return 1;
1.1269 + }
1.1270 + }
1.1271 + return 0;
1.1272 +}
1.1273 +
1.1274 +/***********************************************************************
1.1275 +* NAME
1.1276 +*
1.1277 +* npp_make_fixed - process column with almost identical bounds
1.1278 +*
1.1279 +* SYNOPSIS
1.1280 +*
1.1281 +* #include "glpnpp.h"
1.1282 +* int npp_make_fixed(NPP *npp, NPPCOL *q);
1.1283 +*
1.1284 +* DESCRIPTION
1.1285 +*
1.1286 +* The routine npp_make_fixed processes column q:
1.1287 +*
1.1288 +* l[q] <= x[q] <= u[q], (1)
1.1289 +*
1.1290 +* where -oo < l[q] < u[q] < +oo, i.e. which has both lower and upper
1.1291 +* bounds.
1.1292 +*
1.1293 +* RETURNS
1.1294 +*
1.1295 +* 0 - column bounds have not been changed;
1.1296 +*
1.1297 +* 1 - column has been fixed.
1.1298 +*
1.1299 +* PROBLEM TRANSFORMATION
1.1300 +*
1.1301 +* If bounds of column (1) are very close to each other:
1.1302 +*
1.1303 +* u[q] - l[q] <= eps, (2)
1.1304 +*
1.1305 +* where eps is an absolute tolerance for column value, the column can
1.1306 +* be fixed:
1.1307 +*
1.1308 +* x[q] = s[q], (3)
1.1309 +*
1.1310 +* where s[q] = (l[q] + u[q]) / 2. And if the fixed column value s[q]
1.1311 +* happens to be very close to its nearest integer:
1.1312 +*
1.1313 +* |s[q] - floor(s[q] + 0.5)| <= eps, (4)
1.1314 +*
1.1315 +* it is reasonable to use this nearest integer as the fixed value.
1.1316 +*
1.1317 +* RECOVERING BASIC SOLUTION
1.1318 +*
1.1319 +* In the dual system of the original (as well as transformed) problem
1.1320 +* column q corresponds to the following row:
1.1321 +*
1.1322 +* sum a[i,q] pi[i] + lambda[q] = c[q]. (5)
1.1323 +* i
1.1324 +*
1.1325 +* Since multipliers pi[i] are known for all rows from solution to the
1.1326 +* transformed problem, formula (5) allows computing value of multiplier
1.1327 +* (reduced cost) for column q:
1.1328 +*
1.1329 +* lambda[q] = c[q] - sum a[i,q] pi[i]. (6)
1.1330 +* i
1.1331 +*
1.1332 +* Status of column q in solution to the original problem is determined
1.1333 +* by its status and the sign of its multiplier lambda[q] in solution to
1.1334 +* the transformed problem as follows:
1.1335 +*
1.1336 +* +-----------------------+-----------+--------------------+
1.1337 +* | Status of column q | Sign of | Status of column q |
1.1338 +* | (transformed problem) | lambda[q] | (original problem) |
1.1339 +* +-----------------------+-----------+--------------------+
1.1340 +* | GLP_BS | + / - | GLP_BS |
1.1341 +* | GLP_NS | + | GLP_NL |
1.1342 +* | GLP_NS | - | GLP_NU |
1.1343 +* +-----------------------+-----------+--------------------+
1.1344 +*
1.1345 +* Value of column q in solution to the original problem is the same as
1.1346 +* in solution to the transformed problem.
1.1347 +*
1.1348 +* RECOVERING INTERIOR POINT SOLUTION
1.1349 +*
1.1350 +* Value of column q in solution to the original problem is the same as
1.1351 +* in solution to the transformed problem.
1.1352 +*
1.1353 +* RECOVERING MIP SOLUTION
1.1354 +*
1.1355 +* None needed. */
1.1356 +
1.1357 +struct make_fixed
1.1358 +{ /* column with almost identical bounds */
1.1359 + int q;
1.1360 + /* column reference number */
1.1361 + double c;
1.1362 + /* objective coefficient at x[q] */
1.1363 + NPPLFE *ptr;
1.1364 + /* list of non-zero coefficients a[i,q] */
1.1365 +};
1.1366 +
1.1367 +static int rcv_make_fixed(NPP *npp, void *info);
1.1368 +
1.1369 +int npp_make_fixed(NPP *npp, NPPCOL *q)
1.1370 +{ /* process column with almost identical bounds */
1.1371 + struct make_fixed *info;
1.1372 + NPPAIJ *aij;
1.1373 + NPPLFE *lfe;
1.1374 + double s, eps, nint;
1.1375 + /* the column must be double-bounded */
1.1376 + xassert(q->lb != -DBL_MAX);
1.1377 + xassert(q->ub != +DBL_MAX);
1.1378 + xassert(q->lb < q->ub);
1.1379 + /* check column bounds */
1.1380 + eps = 1e-9 + 1e-12 * fabs(q->lb);
1.1381 + if (q->ub - q->lb > eps) return 0;
1.1382 + /* column bounds are very close to each other */
1.1383 + /* create transformation stack entry */
1.1384 + info = npp_push_tse(npp,
1.1385 + rcv_make_fixed, sizeof(struct make_fixed));
1.1386 + info->q = q->j;
1.1387 + info->c = q->coef;
1.1388 + info->ptr = NULL;
1.1389 + /* save column coefficients a[i,q] (needed for basic solution
1.1390 + only) */
1.1391 + if (npp->sol == GLP_SOL)
1.1392 + { for (aij = q->ptr; aij != NULL; aij = aij->c_next)
1.1393 + { lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE));
1.1394 + lfe->ref = aij->row->i;
1.1395 + lfe->val = aij->val;
1.1396 + lfe->next = info->ptr;
1.1397 + info->ptr = lfe;
1.1398 + }
1.1399 + }
1.1400 + /* compute column fixed value */
1.1401 + s = 0.5 * (q->ub + q->lb);
1.1402 + nint = floor(s + 0.5);
1.1403 + if (fabs(s - nint) <= eps) s = nint;
1.1404 + /* make column q fixed */
1.1405 + q->lb = q->ub = s;
1.1406 + return 1;
1.1407 +}
1.1408 +
1.1409 +static int rcv_make_fixed(NPP *npp, void *_info)
1.1410 +{ /* recover column with almost identical bounds */
1.1411 + struct make_fixed *info = _info;
1.1412 + NPPLFE *lfe;
1.1413 + double lambda;
1.1414 + if (npp->sol == GLP_SOL)
1.1415 + { if (npp->c_stat[info->q] == GLP_BS)
1.1416 + npp->c_stat[info->q] = GLP_BS;
1.1417 + else if (npp->c_stat[info->q] == GLP_NS)
1.1418 + { /* compute multiplier for column q with formula (6) */
1.1419 + lambda = info->c;
1.1420 + for (lfe = info->ptr; lfe != NULL; lfe = lfe->next)
1.1421 + lambda -= lfe->val * npp->r_pi[lfe->ref];
1.1422 + /* assign status to non-basic column */
1.1423 + if (lambda >= 0.0)
1.1424 + npp->c_stat[info->q] = GLP_NL;
1.1425 + else
1.1426 + npp->c_stat[info->q] = GLP_NU;
1.1427 + }
1.1428 + else
1.1429 + { npp_error();
1.1430 + return 1;
1.1431 + }
1.1432 + }
1.1433 + return 0;
1.1434 +}
1.1435 +
1.1436 +/* eof */