1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpnpp02.c	Sun Nov 06 20:59:10 2011 +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, 2011 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 */