lemon-project-template-glpk

annotate deps/glpk/src/glpnpp02.c @ 9:33de93886c88

find changesets by author, revision, files, or words in the commit message
Import GLPK 4.47
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 20:59:10 +0100
parents
children
rev   line source
alpar@9 1 /* glpnpp02.c */
alpar@9 2
alpar@9 3 /***********************************************************************
alpar@9 4 * This code is part of GLPK (GNU Linear Programming Kit).
alpar@9 5 *
alpar@9 6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@9 7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
alpar@9 8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@9 9 * E-mail: <mao@gnu.org>.
alpar@9 10 *
alpar@9 11 * GLPK is free software: you can redistribute it and/or modify it
alpar@9 12 * under the terms of the GNU General Public License as published by
alpar@9 13 * the Free Software Foundation, either version 3 of the License, or
alpar@9 14 * (at your option) any later version.
alpar@9 15 *
alpar@9 16 * GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@9 17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@9 18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@9 19 * License for more details.
alpar@9 20 *
alpar@9 21 * You should have received a copy of the GNU General Public License
alpar@9 22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@9 23 ***********************************************************************/
alpar@9 24
alpar@9 25 #include "glpnpp.h"
alpar@9 26
alpar@9 27 /***********************************************************************
alpar@9 28 * NAME
alpar@9 29 *
alpar@9 30 * npp_free_row - process free (unbounded) row
alpar@9 31 *
alpar@9 32 * SYNOPSIS
alpar@9 33 *
alpar@9 34 * #include "glpnpp.h"
alpar@9 35 * void npp_free_row(NPP *npp, NPPROW *p);
alpar@9 36 *
alpar@9 37 * DESCRIPTION
alpar@9 38 *
alpar@9 39 * The routine npp_free_row processes row p, which is free (i.e. has
alpar@9 40 * no finite bounds):
alpar@9 41 *
alpar@9 42 * -inf < sum a[p,j] x[j] < +inf. (1)
alpar@9 43 * j
alpar@9 44 *
alpar@9 45 * PROBLEM TRANSFORMATION
alpar@9 46 *
alpar@9 47 * Constraint (1) cannot be active, so it is redundant and can be
alpar@9 48 * removed from the original problem.
alpar@9 49 *
alpar@9 50 * Removing row p leads to removing a column of multiplier pi[p] for
alpar@9 51 * this row in the dual system. Since row p has no bounds, pi[p] = 0,
alpar@9 52 * so removing the column does not affect the dual solution.
alpar@9 53 *
alpar@9 54 * RECOVERING BASIC SOLUTION
alpar@9 55 *
alpar@9 56 * In solution to the original problem row p is inactive constraint,
alpar@9 57 * so it is assigned status GLP_BS, and multiplier pi[p] is assigned
alpar@9 58 * zero value.
alpar@9 59 *
alpar@9 60 * RECOVERING INTERIOR-POINT SOLUTION
alpar@9 61 *
alpar@9 62 * In solution to the original problem row p is inactive constraint,
alpar@9 63 * so its multiplier pi[p] is assigned zero value.
alpar@9 64 *
alpar@9 65 * RECOVERING MIP SOLUTION
alpar@9 66 *
alpar@9 67 * None needed. */
alpar@9 68
alpar@9 69 struct free_row
alpar@9 70 { /* free (unbounded) row */
alpar@9 71 int p;
alpar@9 72 /* row reference number */
alpar@9 73 };
alpar@9 74
alpar@9 75 static int rcv_free_row(NPP *npp, void *info);
alpar@9 76
alpar@9 77 void npp_free_row(NPP *npp, NPPROW *p)
alpar@9 78 { /* process free (unbounded) row */
alpar@9 79 struct free_row *info;
alpar@9 80 /* the row must be free */
alpar@9 81 xassert(p->lb == -DBL_MAX && p->ub == +DBL_MAX);
alpar@9 82 /* create transformation stack entry */
alpar@9 83 info = npp_push_tse(npp,
alpar@9 84 rcv_free_row, sizeof(struct free_row));
alpar@9 85 info->p = p->i;
alpar@9 86 /* remove the row from the problem */
alpar@9 87 npp_del_row(npp, p);
alpar@9 88 return;
alpar@9 89 }
alpar@9 90
alpar@9 91 static int rcv_free_row(NPP *npp, void *_info)
alpar@9 92 { /* recover free (unbounded) row */
alpar@9 93 struct free_row *info = _info;
alpar@9 94 if (npp->sol == GLP_SOL)
alpar@9 95 npp->r_stat[info->p] = GLP_BS;
alpar@9 96 if (npp->sol != GLP_MIP)
alpar@9 97 npp->r_pi[info->p] = 0.0;
alpar@9 98 return 0;
alpar@9 99 }
alpar@9 100
alpar@9 101 /***********************************************************************
alpar@9 102 * NAME
alpar@9 103 *
alpar@9 104 * npp_geq_row - process row of 'not less than' type
alpar@9 105 *
alpar@9 106 * SYNOPSIS
alpar@9 107 *
alpar@9 108 * #include "glpnpp.h"
alpar@9 109 * void npp_geq_row(NPP *npp, NPPROW *p);
alpar@9 110 *
alpar@9 111 * DESCRIPTION
alpar@9 112 *
alpar@9 113 * The routine npp_geq_row processes row p, which is 'not less than'
alpar@9 114 * inequality constraint:
alpar@9 115 *
alpar@9 116 * L[p] <= sum a[p,j] x[j] (<= U[p]), (1)
alpar@9 117 * j
alpar@9 118 *
alpar@9 119 * where L[p] < U[p], and upper bound may not exist (U[p] = +oo).
alpar@9 120 *
alpar@9 121 * PROBLEM TRANSFORMATION
alpar@9 122 *
alpar@9 123 * Constraint (1) can be replaced by equality constraint:
alpar@9 124 *
alpar@9 125 * sum a[p,j] x[j] - s = L[p], (2)
alpar@9 126 * j
alpar@9 127 *
alpar@9 128 * where
alpar@9 129 *
alpar@9 130 * 0 <= s (<= U[p] - L[p]) (3)
alpar@9 131 *
alpar@9 132 * is a non-negative surplus variable.
alpar@9 133 *
alpar@9 134 * Since in the primal system there appears column s having the only
alpar@9 135 * non-zero coefficient in row p, in the dual system there appears a
alpar@9 136 * new row:
alpar@9 137 *
alpar@9 138 * (-1) pi[p] + lambda = 0, (4)
alpar@9 139 *
alpar@9 140 * where (-1) is coefficient of column s in row p, pi[p] is multiplier
alpar@9 141 * of row p, lambda is multiplier of column q, 0 is coefficient of
alpar@9 142 * column s in the objective row.
alpar@9 143 *
alpar@9 144 * RECOVERING BASIC SOLUTION
alpar@9 145 *
alpar@9 146 * Status of row p in solution to the original problem is determined
alpar@9 147 * by its status and status of column q in solution to the transformed
alpar@9 148 * problem as follows:
alpar@9 149 *
alpar@9 150 * +--------------------------------------+------------------+
alpar@9 151 * | Transformed problem | Original problem |
alpar@9 152 * +-----------------+--------------------+------------------+
alpar@9 153 * | Status of row p | Status of column s | Status of row p |
alpar@9 154 * +-----------------+--------------------+------------------+
alpar@9 155 * | GLP_BS | GLP_BS | N/A |
alpar@9 156 * | GLP_BS | GLP_NL | GLP_BS |
alpar@9 157 * | GLP_BS | GLP_NU | GLP_BS |
alpar@9 158 * | GLP_NS | GLP_BS | GLP_BS |
alpar@9 159 * | GLP_NS | GLP_NL | GLP_NL |
alpar@9 160 * | GLP_NS | GLP_NU | GLP_NU |
alpar@9 161 * +-----------------+--------------------+------------------+
alpar@9 162 *
alpar@9 163 * Value of row multiplier pi[p] in solution to the original problem
alpar@9 164 * is the same as in solution to the transformed problem.
alpar@9 165 *
alpar@9 166 * 1. In solution to the transformed problem row p and column q cannot
alpar@9 167 * be basic at the same time; otherwise the basis matrix would have
alpar@9 168 * two linear dependent columns: unity column of auxiliary variable
alpar@9 169 * of row p and unity column of variable s.
alpar@9 170 *
alpar@9 171 * 2. Though in the transformed problem row p is equality constraint,
alpar@9 172 * it may be basic due to primal degenerate solution.
alpar@9 173 *
alpar@9 174 * RECOVERING INTERIOR-POINT SOLUTION
alpar@9 175 *
alpar@9 176 * Value of row multiplier pi[p] in solution to the original problem
alpar@9 177 * is the same as in solution to the transformed problem.
alpar@9 178 *
alpar@9 179 * RECOVERING MIP SOLUTION
alpar@9 180 *
alpar@9 181 * None needed. */
alpar@9 182
alpar@9 183 struct ineq_row
alpar@9 184 { /* inequality constraint row */
alpar@9 185 int p;
alpar@9 186 /* row reference number */
alpar@9 187 int s;
alpar@9 188 /* column reference number for slack/surplus variable */
alpar@9 189 };
alpar@9 190
alpar@9 191 static int rcv_geq_row(NPP *npp, void *info);
alpar@9 192
alpar@9 193 void npp_geq_row(NPP *npp, NPPROW *p)
alpar@9 194 { /* process row of 'not less than' type */
alpar@9 195 struct ineq_row *info;
alpar@9 196 NPPCOL *s;
alpar@9 197 /* the row must have lower bound */
alpar@9 198 xassert(p->lb != -DBL_MAX);
alpar@9 199 xassert(p->lb < p->ub);
alpar@9 200 /* create column for surplus variable */
alpar@9 201 s = npp_add_col(npp);
alpar@9 202 s->lb = 0.0;
alpar@9 203 s->ub = (p->ub == +DBL_MAX ? +DBL_MAX : p->ub - p->lb);
alpar@9 204 /* and add it to the transformed problem */
alpar@9 205 npp_add_aij(npp, p, s, -1.0);
alpar@9 206 /* create transformation stack entry */
alpar@9 207 info = npp_push_tse(npp,
alpar@9 208 rcv_geq_row, sizeof(struct ineq_row));
alpar@9 209 info->p = p->i;
alpar@9 210 info->s = s->j;
alpar@9 211 /* replace the row by equality constraint */
alpar@9 212 p->ub = p->lb;
alpar@9 213 return;
alpar@9 214 }
alpar@9 215
alpar@9 216 static int rcv_geq_row(NPP *npp, void *_info)
alpar@9 217 { /* recover row of 'not less than' type */
alpar@9 218 struct ineq_row *info = _info;
alpar@9 219 if (npp->sol == GLP_SOL)
alpar@9 220 { if (npp->r_stat[info->p] == GLP_BS)
alpar@9 221 { if (npp->c_stat[info->s] == GLP_BS)
alpar@9 222 { npp_error();
alpar@9 223 return 1;
alpar@9 224 }
alpar@9 225 else if (npp->c_stat[info->s] == GLP_NL ||
alpar@9 226 npp->c_stat[info->s] == GLP_NU)
alpar@9 227 npp->r_stat[info->p] = GLP_BS;
alpar@9 228 else
alpar@9 229 { npp_error();
alpar@9 230 return 1;
alpar@9 231 }
alpar@9 232 }
alpar@9 233 else if (npp->r_stat[info->p] == GLP_NS)
alpar@9 234 { if (npp->c_stat[info->s] == GLP_BS)
alpar@9 235 npp->r_stat[info->p] = GLP_BS;
alpar@9 236 else if (npp->c_stat[info->s] == GLP_NL)
alpar@9 237 npp->r_stat[info->p] = GLP_NL;
alpar@9 238 else if (npp->c_stat[info->s] == GLP_NU)
alpar@9 239 npp->r_stat[info->p] = GLP_NU;
alpar@9 240 else
alpar@9 241 { npp_error();
alpar@9 242 return 1;
alpar@9 243 }
alpar@9 244 }
alpar@9 245 else
alpar@9 246 { npp_error();
alpar@9 247 return 1;
alpar@9 248 }
alpar@9 249 }
alpar@9 250 return 0;
alpar@9 251 }
alpar@9 252
alpar@9 253 /***********************************************************************
alpar@9 254 * NAME
alpar@9 255 *
alpar@9 256 * npp_leq_row - process row of 'not greater than' type
alpar@9 257 *
alpar@9 258 * SYNOPSIS
alpar@9 259 *
alpar@9 260 * #include "glpnpp.h"
alpar@9 261 * void npp_leq_row(NPP *npp, NPPROW *p);
alpar@9 262 *
alpar@9 263 * DESCRIPTION
alpar@9 264 *
alpar@9 265 * The routine npp_leq_row processes row p, which is 'not greater than'
alpar@9 266 * inequality constraint:
alpar@9 267 *
alpar@9 268 * (L[p] <=) sum a[p,j] x[j] <= U[p], (1)
alpar@9 269 * j
alpar@9 270 *
alpar@9 271 * where L[p] < U[p], and lower bound may not exist (L[p] = +oo).
alpar@9 272 *
alpar@9 273 * PROBLEM TRANSFORMATION
alpar@9 274 *
alpar@9 275 * Constraint (1) can be replaced by equality constraint:
alpar@9 276 *
alpar@9 277 * sum a[p,j] x[j] + s = L[p], (2)
alpar@9 278 * j
alpar@9 279 *
alpar@9 280 * where
alpar@9 281 *
alpar@9 282 * 0 <= s (<= U[p] - L[p]) (3)
alpar@9 283 *
alpar@9 284 * is a non-negative slack variable.
alpar@9 285 *
alpar@9 286 * Since in the primal system there appears column s having the only
alpar@9 287 * non-zero coefficient in row p, in the dual system there appears a
alpar@9 288 * new row:
alpar@9 289 *
alpar@9 290 * (+1) pi[p] + lambda = 0, (4)
alpar@9 291 *
alpar@9 292 * where (+1) is coefficient of column s in row p, pi[p] is multiplier
alpar@9 293 * of row p, lambda is multiplier of column q, 0 is coefficient of
alpar@9 294 * column s in the objective row.
alpar@9 295 *
alpar@9 296 * RECOVERING BASIC SOLUTION
alpar@9 297 *
alpar@9 298 * Status of row p in solution to the original problem is determined
alpar@9 299 * by its status and status of column q in solution to the transformed
alpar@9 300 * problem as follows:
alpar@9 301 *
alpar@9 302 * +--------------------------------------+------------------+
alpar@9 303 * | Transformed problem | Original problem |
alpar@9 304 * +-----------------+--------------------+------------------+
alpar@9 305 * | Status of row p | Status of column s | Status of row p |
alpar@9 306 * +-----------------+--------------------+------------------+
alpar@9 307 * | GLP_BS | GLP_BS | N/A |
alpar@9 308 * | GLP_BS | GLP_NL | GLP_BS |
alpar@9 309 * | GLP_BS | GLP_NU | GLP_BS |
alpar@9 310 * | GLP_NS | GLP_BS | GLP_BS |
alpar@9 311 * | GLP_NS | GLP_NL | GLP_NU |
alpar@9 312 * | GLP_NS | GLP_NU | GLP_NL |
alpar@9 313 * +-----------------+--------------------+------------------+
alpar@9 314 *
alpar@9 315 * Value of row multiplier pi[p] in solution to the original problem
alpar@9 316 * is the same as in solution to the transformed problem.
alpar@9 317 *
alpar@9 318 * 1. In solution to the transformed problem row p and column q cannot
alpar@9 319 * be basic at the same time; otherwise the basis matrix would have
alpar@9 320 * two linear dependent columns: unity column of auxiliary variable
alpar@9 321 * of row p and unity column of variable s.
alpar@9 322 *
alpar@9 323 * 2. Though in the transformed problem row p is equality constraint,
alpar@9 324 * it may be basic due to primal degeneracy.
alpar@9 325 *
alpar@9 326 * RECOVERING INTERIOR-POINT SOLUTION
alpar@9 327 *
alpar@9 328 * Value of row multiplier pi[p] in solution to the original problem
alpar@9 329 * is the same as in solution to the transformed problem.
alpar@9 330 *
alpar@9 331 * RECOVERING MIP SOLUTION
alpar@9 332 *
alpar@9 333 * None needed. */
alpar@9 334
alpar@9 335 static int rcv_leq_row(NPP *npp, void *info);
alpar@9 336
alpar@9 337 void npp_leq_row(NPP *npp, NPPROW *p)
alpar@9 338 { /* process row of 'not greater than' type */
alpar@9 339 struct ineq_row *info;
alpar@9 340 NPPCOL *s;
alpar@9 341 /* the row must have upper bound */
alpar@9 342 xassert(p->ub != +DBL_MAX);
alpar@9 343 xassert(p->lb < p->ub);
alpar@9 344 /* create column for slack variable */
alpar@9 345 s = npp_add_col(npp);
alpar@9 346 s->lb = 0.0;
alpar@9 347 s->ub = (p->lb == -DBL_MAX ? +DBL_MAX : p->ub - p->lb);
alpar@9 348 /* and add it to the transformed problem */
alpar@9 349 npp_add_aij(npp, p, s, +1.0);
alpar@9 350 /* create transformation stack entry */
alpar@9 351 info = npp_push_tse(npp,
alpar@9 352 rcv_leq_row, sizeof(struct ineq_row));
alpar@9 353 info->p = p->i;
alpar@9 354 info->s = s->j;
alpar@9 355 /* replace the row by equality constraint */
alpar@9 356 p->lb = p->ub;
alpar@9 357 return;
alpar@9 358 }
alpar@9 359
alpar@9 360 static int rcv_leq_row(NPP *npp, void *_info)
alpar@9 361 { /* recover row of 'not greater than' type */
alpar@9 362 struct ineq_row *info = _info;
alpar@9 363 if (npp->sol == GLP_SOL)
alpar@9 364 { if (npp->r_stat[info->p] == GLP_BS)
alpar@9 365 { if (npp->c_stat[info->s] == GLP_BS)
alpar@9 366 { npp_error();
alpar@9 367 return 1;
alpar@9 368 }
alpar@9 369 else if (npp->c_stat[info->s] == GLP_NL ||
alpar@9 370 npp->c_stat[info->s] == GLP_NU)
alpar@9 371 npp->r_stat[info->p] = GLP_BS;
alpar@9 372 else
alpar@9 373 { npp_error();
alpar@9 374 return 1;
alpar@9 375 }
alpar@9 376 }
alpar@9 377 else if (npp->r_stat[info->p] == GLP_NS)
alpar@9 378 { if (npp->c_stat[info->s] == GLP_BS)
alpar@9 379 npp->r_stat[info->p] = GLP_BS;
alpar@9 380 else if (npp->c_stat[info->s] == GLP_NL)
alpar@9 381 npp->r_stat[info->p] = GLP_NU;
alpar@9 382 else if (npp->c_stat[info->s] == GLP_NU)
alpar@9 383 npp->r_stat[info->p] = GLP_NL;
alpar@9 384 else
alpar@9 385 { npp_error();
alpar@9 386 return 1;
alpar@9 387 }
alpar@9 388 }
alpar@9 389 else
alpar@9 390 { npp_error();
alpar@9 391 return 1;
alpar@9 392 }
alpar@9 393 }
alpar@9 394 return 0;
alpar@9 395 }
alpar@9 396
alpar@9 397 /***********************************************************************
alpar@9 398 * NAME
alpar@9 399 *
alpar@9 400 * npp_free_col - process free (unbounded) column
alpar@9 401 *
alpar@9 402 * SYNOPSIS
alpar@9 403 *
alpar@9 404 * #include "glpnpp.h"
alpar@9 405 * void npp_free_col(NPP *npp, NPPCOL *q);
alpar@9 406 *
alpar@9 407 * DESCRIPTION
alpar@9 408 *
alpar@9 409 * The routine npp_free_col processes column q, which is free (i.e. has
alpar@9 410 * no finite bounds):
alpar@9 411 *
alpar@9 412 * -oo < x[q] < +oo. (1)
alpar@9 413 *
alpar@9 414 * PROBLEM TRANSFORMATION
alpar@9 415 *
alpar@9 416 * Free (unbounded) variable can be replaced by the difference of two
alpar@9 417 * non-negative variables:
alpar@9 418 *
alpar@9 419 * x[q] = s' - s'', s', s'' >= 0. (2)
alpar@9 420 *
alpar@9 421 * Assuming that in the transformed problem x[q] becomes s',
alpar@9 422 * transformation (2) causes new column s'' to appear, which differs
alpar@9 423 * from column s' only in the sign of coefficients in constraint and
alpar@9 424 * objective rows. Thus, if in the dual system the following row
alpar@9 425 * corresponds to column s':
alpar@9 426 *
alpar@9 427 * sum a[i,q] pi[i] + lambda' = c[q], (3)
alpar@9 428 * i
alpar@9 429 *
alpar@9 430 * the row which corresponds to column s'' is the following:
alpar@9 431 *
alpar@9 432 * sum (-a[i,q]) pi[i] + lambda'' = -c[q]. (4)
alpar@9 433 * i
alpar@9 434 *
alpar@9 435 * Then from (3) and (4) it follows that:
alpar@9 436 *
alpar@9 437 * lambda' + lambda'' = 0 => lambda' = lmabda'' = 0, (5)
alpar@9 438 *
alpar@9 439 * where lambda' and lambda'' are multipliers for columns s' and s'',
alpar@9 440 * resp.
alpar@9 441 *
alpar@9 442 * RECOVERING BASIC SOLUTION
alpar@9 443 *
alpar@9 444 * With respect to (5) status of column q in solution to the original
alpar@9 445 * problem is determined by statuses of columns s' and s'' in solution
alpar@9 446 * to the transformed problem as follows:
alpar@9 447 *
alpar@9 448 * +--------------------------------------+------------------+
alpar@9 449 * | Transformed problem | Original problem |
alpar@9 450 * +------------------+-------------------+------------------+
alpar@9 451 * | Status of col s' | Status of col s'' | Status of col q |
alpar@9 452 * +------------------+-------------------+------------------+
alpar@9 453 * | GLP_BS | GLP_BS | N/A |
alpar@9 454 * | GLP_BS | GLP_NL | GLP_BS |
alpar@9 455 * | GLP_NL | GLP_BS | GLP_BS |
alpar@9 456 * | GLP_NL | GLP_NL | GLP_NF |
alpar@9 457 * +------------------+-------------------+------------------+
alpar@9 458 *
alpar@9 459 * Value of column q is computed with formula (2).
alpar@9 460 *
alpar@9 461 * 1. In solution to the transformed problem columns s' and s'' cannot
alpar@9 462 * be basic at the same time, because they differ only in the sign,
alpar@9 463 * hence, are linear dependent.
alpar@9 464 *
alpar@9 465 * 2. Though column q is free, it can be non-basic due to dual
alpar@9 466 * degeneracy.
alpar@9 467 *
alpar@9 468 * 3. If column q is integral, columns s' and s'' are also integral.
alpar@9 469 *
alpar@9 470 * RECOVERING INTERIOR-POINT SOLUTION
alpar@9 471 *
alpar@9 472 * Value of column q is computed with formula (2).
alpar@9 473 *
alpar@9 474 * RECOVERING MIP SOLUTION
alpar@9 475 *
alpar@9 476 * Value of column q is computed with formula (2). */
alpar@9 477
alpar@9 478 struct free_col
alpar@9 479 { /* free (unbounded) column */
alpar@9 480 int q;
alpar@9 481 /* column reference number for variables x[q] and s' */
alpar@9 482 int s;
alpar@9 483 /* column reference number for variable s'' */
alpar@9 484 };
alpar@9 485
alpar@9 486 static int rcv_free_col(NPP *npp, void *info);
alpar@9 487
alpar@9 488 void npp_free_col(NPP *npp, NPPCOL *q)
alpar@9 489 { /* process free (unbounded) column */
alpar@9 490 struct free_col *info;
alpar@9 491 NPPCOL *s;
alpar@9 492 NPPAIJ *aij;
alpar@9 493 /* the column must be free */
alpar@9 494 xassert(q->lb == -DBL_MAX && q->ub == +DBL_MAX);
alpar@9 495 /* variable x[q] becomes s' */
alpar@9 496 q->lb = 0.0, q->ub = +DBL_MAX;
alpar@9 497 /* create variable s'' */
alpar@9 498 s = npp_add_col(npp);
alpar@9 499 s->is_int = q->is_int;
alpar@9 500 s->lb = 0.0, s->ub = +DBL_MAX;
alpar@9 501 /* duplicate objective coefficient */
alpar@9 502 s->coef = -q->coef;
alpar@9 503 /* duplicate column of the constraint matrix */
alpar@9 504 for (aij = q->ptr; aij != NULL; aij = aij->c_next)
alpar@9 505 npp_add_aij(npp, aij->row, s, -aij->val);
alpar@9 506 /* create transformation stack entry */
alpar@9 507 info = npp_push_tse(npp,
alpar@9 508 rcv_free_col, sizeof(struct free_col));
alpar@9 509 info->q = q->j;
alpar@9 510 info->s = s->j;
alpar@9 511 return;
alpar@9 512 }
alpar@9 513
alpar@9 514 static int rcv_free_col(NPP *npp, void *_info)
alpar@9 515 { /* recover free (unbounded) column */
alpar@9 516 struct free_col *info = _info;
alpar@9 517 if (npp->sol == GLP_SOL)
alpar@9 518 { if (npp->c_stat[info->q] == GLP_BS)
alpar@9 519 { if (npp->c_stat[info->s] == GLP_BS)
alpar@9 520 { npp_error();
alpar@9 521 return 1;
alpar@9 522 }
alpar@9 523 else if (npp->c_stat[info->s] == GLP_NL)
alpar@9 524 npp->c_stat[info->q] = GLP_BS;
alpar@9 525 else
alpar@9 526 { npp_error();
alpar@9 527 return -1;
alpar@9 528 }
alpar@9 529 }
alpar@9 530 else if (npp->c_stat[info->q] == GLP_NL)
alpar@9 531 { if (npp->c_stat[info->s] == GLP_BS)
alpar@9 532 npp->c_stat[info->q] = GLP_BS;
alpar@9 533 else if (npp->c_stat[info->s] == GLP_NL)
alpar@9 534 npp->c_stat[info->q] = GLP_NF;
alpar@9 535 else
alpar@9 536 { npp_error();
alpar@9 537 return -1;
alpar@9 538 }
alpar@9 539 }
alpar@9 540 else
alpar@9 541 { npp_error();
alpar@9 542 return -1;
alpar@9 543 }
alpar@9 544 }
alpar@9 545 /* compute value of x[q] with formula (2) */
alpar@9 546 npp->c_value[info->q] -= npp->c_value[info->s];
alpar@9 547 return 0;
alpar@9 548 }
alpar@9 549
alpar@9 550 /***********************************************************************
alpar@9 551 * NAME
alpar@9 552 *
alpar@9 553 * npp_lbnd_col - process column with (non-zero) lower bound
alpar@9 554 *
alpar@9 555 * SYNOPSIS
alpar@9 556 *
alpar@9 557 * #include "glpnpp.h"
alpar@9 558 * void npp_lbnd_col(NPP *npp, NPPCOL *q);
alpar@9 559 *
alpar@9 560 * DESCRIPTION
alpar@9 561 *
alpar@9 562 * The routine npp_lbnd_col processes column q, which has (non-zero)
alpar@9 563 * lower bound:
alpar@9 564 *
alpar@9 565 * l[q] <= x[q] (<= u[q]), (1)
alpar@9 566 *
alpar@9 567 * where l[q] < u[q], and upper bound may not exist (u[q] = +oo).
alpar@9 568 *
alpar@9 569 * PROBLEM TRANSFORMATION
alpar@9 570 *
alpar@9 571 * Column q can be replaced as follows:
alpar@9 572 *
alpar@9 573 * x[q] = l[q] + s, (2)
alpar@9 574 *
alpar@9 575 * where
alpar@9 576 *
alpar@9 577 * 0 <= s (<= u[q] - l[q]) (3)
alpar@9 578 *
alpar@9 579 * is a non-negative variable.
alpar@9 580 *
alpar@9 581 * Substituting x[q] from (2) into the objective row, we have:
alpar@9 582 *
alpar@9 583 * z = sum c[j] x[j] + c0 =
alpar@9 584 * j
alpar@9 585 *
alpar@9 586 * = sum c[j] x[j] + c[q] x[q] + c0 =
alpar@9 587 * j!=q
alpar@9 588 *
alpar@9 589 * = sum c[j] x[j] + c[q] (l[q] + s) + c0 =
alpar@9 590 * j!=q
alpar@9 591 *
alpar@9 592 * = sum c[j] x[j] + c[q] s + c~0,
alpar@9 593 *
alpar@9 594 * where
alpar@9 595 *
alpar@9 596 * c~0 = c0 + c[q] l[q] (4)
alpar@9 597 *
alpar@9 598 * is the constant term of the objective in the transformed problem.
alpar@9 599 * Similarly, substituting x[q] into constraint row i, we have:
alpar@9 600 *
alpar@9 601 * L[i] <= sum a[i,j] x[j] <= U[i] ==>
alpar@9 602 * j
alpar@9 603 *
alpar@9 604 * L[i] <= sum a[i,j] x[j] + a[i,q] x[q] <= U[i] ==>
alpar@9 605 * j!=q
alpar@9 606 *
alpar@9 607 * L[i] <= sum a[i,j] x[j] + a[i,q] (l[q] + s) <= U[i] ==>
alpar@9 608 * j!=q
alpar@9 609 *
alpar@9 610 * L~[i] <= sum a[i,j] x[j] + a[i,q] s <= U~[i],
alpar@9 611 * j!=q
alpar@9 612 *
alpar@9 613 * where
alpar@9 614 *
alpar@9 615 * L~[i] = L[i] - a[i,q] l[q], U~[i] = U[i] - a[i,q] l[q] (5)
alpar@9 616 *
alpar@9 617 * are lower and upper bounds of row i in the transformed problem,
alpar@9 618 * resp.
alpar@9 619 *
alpar@9 620 * Transformation (2) does not affect the dual system.
alpar@9 621 *
alpar@9 622 * RECOVERING BASIC SOLUTION
alpar@9 623 *
alpar@9 624 * Status of column q in solution to the original problem is the same
alpar@9 625 * as in solution to the transformed problem (GLP_BS, GLP_NL or GLP_NU).
alpar@9 626 * Value of column q is computed with formula (2).
alpar@9 627 *
alpar@9 628 * RECOVERING INTERIOR-POINT SOLUTION
alpar@9 629 *
alpar@9 630 * Value of column q is computed with formula (2).
alpar@9 631 *
alpar@9 632 * RECOVERING MIP SOLUTION
alpar@9 633 *
alpar@9 634 * Value of column q is computed with formula (2). */
alpar@9 635
alpar@9 636 struct bnd_col
alpar@9 637 { /* bounded column */
alpar@9 638 int q;
alpar@9 639 /* column reference number for variables x[q] and s */
alpar@9 640 double bnd;
alpar@9 641 /* lower/upper bound l[q] or u[q] */
alpar@9 642 };
alpar@9 643
alpar@9 644 static int rcv_lbnd_col(NPP *npp, void *info);
alpar@9 645
alpar@9 646 void npp_lbnd_col(NPP *npp, NPPCOL *q)
alpar@9 647 { /* process column with (non-zero) lower bound */
alpar@9 648 struct bnd_col *info;
alpar@9 649 NPPROW *i;
alpar@9 650 NPPAIJ *aij;
alpar@9 651 /* the column must have non-zero lower bound */
alpar@9 652 xassert(q->lb != 0.0);
alpar@9 653 xassert(q->lb != -DBL_MAX);
alpar@9 654 xassert(q->lb < q->ub);
alpar@9 655 /* create transformation stack entry */
alpar@9 656 info = npp_push_tse(npp,
alpar@9 657 rcv_lbnd_col, sizeof(struct bnd_col));
alpar@9 658 info->q = q->j;
alpar@9 659 info->bnd = q->lb;
alpar@9 660 /* substitute x[q] into objective row */
alpar@9 661 npp->c0 += q->coef * q->lb;
alpar@9 662 /* substitute x[q] into constraint rows */
alpar@9 663 for (aij = q->ptr; aij != NULL; aij = aij->c_next)
alpar@9 664 { i = aij->row;
alpar@9 665 if (i->lb == i->ub)
alpar@9 666 i->ub = (i->lb -= aij->val * q->lb);
alpar@9 667 else
alpar@9 668 { if (i->lb != -DBL_MAX)
alpar@9 669 i->lb -= aij->val * q->lb;
alpar@9 670 if (i->ub != +DBL_MAX)
alpar@9 671 i->ub -= aij->val * q->lb;
alpar@9 672 }
alpar@9 673 }
alpar@9 674 /* column x[q] becomes column s */
alpar@9 675 if (q->ub != +DBL_MAX)
alpar@9 676 q->ub -= q->lb;
alpar@9 677 q->lb = 0.0;
alpar@9 678 return;
alpar@9 679 }
alpar@9 680
alpar@9 681 static int rcv_lbnd_col(NPP *npp, void *_info)
alpar@9 682 { /* recover column with (non-zero) lower bound */
alpar@9 683 struct bnd_col *info = _info;
alpar@9 684 if (npp->sol == GLP_SOL)
alpar@9 685 { if (npp->c_stat[info->q] == GLP_BS ||
alpar@9 686 npp->c_stat[info->q] == GLP_NL ||
alpar@9 687 npp->c_stat[info->q] == GLP_NU)
alpar@9 688 npp->c_stat[info->q] = npp->c_stat[info->q];
alpar@9 689 else
alpar@9 690 { npp_error();
alpar@9 691 return 1;
alpar@9 692 }
alpar@9 693 }
alpar@9 694 /* compute value of x[q] with formula (2) */
alpar@9 695 npp->c_value[info->q] = info->bnd + npp->c_value[info->q];
alpar@9 696 return 0;
alpar@9 697 }
alpar@9 698
alpar@9 699 /***********************************************************************
alpar@9 700 * NAME
alpar@9 701 *
alpar@9 702 * npp_ubnd_col - process column with upper bound
alpar@9 703 *
alpar@9 704 * SYNOPSIS
alpar@9 705 *
alpar@9 706 * #include "glpnpp.h"
alpar@9 707 * void npp_ubnd_col(NPP *npp, NPPCOL *q);
alpar@9 708 *
alpar@9 709 * DESCRIPTION
alpar@9 710 *
alpar@9 711 * The routine npp_ubnd_col processes column q, which has upper bound:
alpar@9 712 *
alpar@9 713 * (l[q] <=) x[q] <= u[q], (1)
alpar@9 714 *
alpar@9 715 * where l[q] < u[q], and lower bound may not exist (l[q] = -oo).
alpar@9 716 *
alpar@9 717 * PROBLEM TRANSFORMATION
alpar@9 718 *
alpar@9 719 * Column q can be replaced as follows:
alpar@9 720 *
alpar@9 721 * x[q] = u[q] - s, (2)
alpar@9 722 *
alpar@9 723 * where
alpar@9 724 *
alpar@9 725 * 0 <= s (<= u[q] - l[q]) (3)
alpar@9 726 *
alpar@9 727 * is a non-negative variable.
alpar@9 728 *
alpar@9 729 * Substituting x[q] from (2) into the objective row, we have:
alpar@9 730 *
alpar@9 731 * z = sum c[j] x[j] + c0 =
alpar@9 732 * j
alpar@9 733 *
alpar@9 734 * = sum c[j] x[j] + c[q] x[q] + c0 =
alpar@9 735 * j!=q
alpar@9 736 *
alpar@9 737 * = sum c[j] x[j] + c[q] (u[q] - s) + c0 =
alpar@9 738 * j!=q
alpar@9 739 *
alpar@9 740 * = sum c[j] x[j] - c[q] s + c~0,
alpar@9 741 *
alpar@9 742 * where
alpar@9 743 *
alpar@9 744 * c~0 = c0 + c[q] u[q] (4)
alpar@9 745 *
alpar@9 746 * is the constant term of the objective in the transformed problem.
alpar@9 747 * Similarly, substituting x[q] into constraint row i, we have:
alpar@9 748 *
alpar@9 749 * L[i] <= sum a[i,j] x[j] <= U[i] ==>
alpar@9 750 * j
alpar@9 751 *
alpar@9 752 * L[i] <= sum a[i,j] x[j] + a[i,q] x[q] <= U[i] ==>
alpar@9 753 * j!=q
alpar@9 754 *
alpar@9 755 * L[i] <= sum a[i,j] x[j] + a[i,q] (u[q] - s) <= U[i] ==>
alpar@9 756 * j!=q
alpar@9 757 *
alpar@9 758 * L~[i] <= sum a[i,j] x[j] - a[i,q] s <= U~[i],
alpar@9 759 * j!=q
alpar@9 760 *
alpar@9 761 * where
alpar@9 762 *
alpar@9 763 * L~[i] = L[i] - a[i,q] u[q], U~[i] = U[i] - a[i,q] u[q] (5)
alpar@9 764 *
alpar@9 765 * are lower and upper bounds of row i in the transformed problem,
alpar@9 766 * resp.
alpar@9 767 *
alpar@9 768 * Note that in the transformed problem coefficients c[q] and a[i,q]
alpar@9 769 * change their sign. Thus, the row of the dual system corresponding to
alpar@9 770 * column q:
alpar@9 771 *
alpar@9 772 * sum a[i,q] pi[i] + lambda[q] = c[q] (6)
alpar@9 773 * i
alpar@9 774 *
alpar@9 775 * in the transformed problem becomes the following:
alpar@9 776 *
alpar@9 777 * sum (-a[i,q]) pi[i] + lambda[s] = -c[q]. (7)
alpar@9 778 * i
alpar@9 779 *
alpar@9 780 * Therefore:
alpar@9 781 *
alpar@9 782 * lambda[q] = - lambda[s], (8)
alpar@9 783 *
alpar@9 784 * where lambda[q] is multiplier for column q, lambda[s] is multiplier
alpar@9 785 * for column s.
alpar@9 786 *
alpar@9 787 * RECOVERING BASIC SOLUTION
alpar@9 788 *
alpar@9 789 * With respect to (8) status of column q in solution to the original
alpar@9 790 * problem is determined by status of column s in solution to the
alpar@9 791 * transformed problem as follows:
alpar@9 792 *
alpar@9 793 * +-----------------------+--------------------+
alpar@9 794 * | Status of column s | Status of column q |
alpar@9 795 * | (transformed problem) | (original problem) |
alpar@9 796 * +-----------------------+--------------------+
alpar@9 797 * | GLP_BS | GLP_BS |
alpar@9 798 * | GLP_NL | GLP_NU |
alpar@9 799 * | GLP_NU | GLP_NL |
alpar@9 800 * +-----------------------+--------------------+
alpar@9 801 *
alpar@9 802 * Value of column q is computed with formula (2).
alpar@9 803 *
alpar@9 804 * RECOVERING INTERIOR-POINT SOLUTION
alpar@9 805 *
alpar@9 806 * Value of column q is computed with formula (2).
alpar@9 807 *
alpar@9 808 * RECOVERING MIP SOLUTION
alpar@9 809 *
alpar@9 810 * Value of column q is computed with formula (2). */
alpar@9 811
alpar@9 812 static int rcv_ubnd_col(NPP *npp, void *info);
alpar@9 813
alpar@9 814 void npp_ubnd_col(NPP *npp, NPPCOL *q)
alpar@9 815 { /* process column with upper bound */
alpar@9 816 struct bnd_col *info;
alpar@9 817 NPPROW *i;
alpar@9 818 NPPAIJ *aij;
alpar@9 819 /* the column must have upper bound */
alpar@9 820 xassert(q->ub != +DBL_MAX);
alpar@9 821 xassert(q->lb < q->ub);
alpar@9 822 /* create transformation stack entry */
alpar@9 823 info = npp_push_tse(npp,
alpar@9 824 rcv_ubnd_col, sizeof(struct bnd_col));
alpar@9 825 info->q = q->j;
alpar@9 826 info->bnd = q->ub;
alpar@9 827 /* substitute x[q] into objective row */
alpar@9 828 npp->c0 += q->coef * q->ub;
alpar@9 829 q->coef = -q->coef;
alpar@9 830 /* substitute x[q] into constraint rows */
alpar@9 831 for (aij = q->ptr; aij != NULL; aij = aij->c_next)
alpar@9 832 { i = aij->row;
alpar@9 833 if (i->lb == i->ub)
alpar@9 834 i->ub = (i->lb -= aij->val * q->ub);
alpar@9 835 else
alpar@9 836 { if (i->lb != -DBL_MAX)
alpar@9 837 i->lb -= aij->val * q->ub;
alpar@9 838 if (i->ub != +DBL_MAX)
alpar@9 839 i->ub -= aij->val * q->ub;
alpar@9 840 }
alpar@9 841 aij->val = -aij->val;
alpar@9 842 }
alpar@9 843 /* column x[q] becomes column s */
alpar@9 844 if (q->lb != -DBL_MAX)
alpar@9 845 q->ub -= q->lb;
alpar@9 846 else
alpar@9 847 q->ub = +DBL_MAX;
alpar@9 848 q->lb = 0.0;
alpar@9 849 return;
alpar@9 850 }
alpar@9 851
alpar@9 852 static int rcv_ubnd_col(NPP *npp, void *_info)
alpar@9 853 { /* recover column with upper bound */
alpar@9 854 struct bnd_col *info = _info;
alpar@9 855 if (npp->sol == GLP_BS)
alpar@9 856 { if (npp->c_stat[info->q] == GLP_BS)
alpar@9 857 npp->c_stat[info->q] = GLP_BS;
alpar@9 858 else if (npp->c_stat[info->q] == GLP_NL)
alpar@9 859 npp->c_stat[info->q] = GLP_NU;
alpar@9 860 else if (npp->c_stat[info->q] == GLP_NU)
alpar@9 861 npp->c_stat[info->q] = GLP_NL;
alpar@9 862 else
alpar@9 863 { npp_error();
alpar@9 864 return 1;
alpar@9 865 }
alpar@9 866 }
alpar@9 867 /* compute value of x[q] with formula (2) */
alpar@9 868 npp->c_value[info->q] = info->bnd - npp->c_value[info->q];
alpar@9 869 return 0;
alpar@9 870 }
alpar@9 871
alpar@9 872 /***********************************************************************
alpar@9 873 * NAME
alpar@9 874 *
alpar@9 875 * npp_dbnd_col - process non-negative column with upper bound
alpar@9 876 *
alpar@9 877 * SYNOPSIS
alpar@9 878 *
alpar@9 879 * #include "glpnpp.h"
alpar@9 880 * void npp_dbnd_col(NPP *npp, NPPCOL *q);
alpar@9 881 *
alpar@9 882 * DESCRIPTION
alpar@9 883 *
alpar@9 884 * The routine npp_dbnd_col processes column q, which is non-negative
alpar@9 885 * and has upper bound:
alpar@9 886 *
alpar@9 887 * 0 <= x[q] <= u[q], (1)
alpar@9 888 *
alpar@9 889 * where u[q] > 0.
alpar@9 890 *
alpar@9 891 * PROBLEM TRANSFORMATION
alpar@9 892 *
alpar@9 893 * Upper bound of column q can be replaced by the following equality
alpar@9 894 * constraint:
alpar@9 895 *
alpar@9 896 * x[q] + s = u[q], (2)
alpar@9 897 *
alpar@9 898 * where s >= 0 is a non-negative complement variable.
alpar@9 899 *
alpar@9 900 * Since in the primal system along with new row (2) there appears a
alpar@9 901 * new column s having the only non-zero coefficient in this row, in
alpar@9 902 * the dual system there appears a new row:
alpar@9 903 *
alpar@9 904 * (+1)pi + lambda[s] = 0, (3)
alpar@9 905 *
alpar@9 906 * where (+1) is coefficient at column s in row (2), pi is multiplier
alpar@9 907 * for row (2), lambda[s] is multiplier for column s, 0 is coefficient
alpar@9 908 * at column s in the objective row.
alpar@9 909 *
alpar@9 910 * RECOVERING BASIC SOLUTION
alpar@9 911 *
alpar@9 912 * Status of column q in solution to the original problem is determined
alpar@9 913 * by its status and status of column s in solution to the transformed
alpar@9 914 * problem as follows:
alpar@9 915 *
alpar@9 916 * +-----------------------------------+------------------+
alpar@9 917 * | Transformed problem | Original problem |
alpar@9 918 * +-----------------+-----------------+------------------+
alpar@9 919 * | Status of col q | Status of col s | Status of col q |
alpar@9 920 * +-----------------+-----------------+------------------+
alpar@9 921 * | GLP_BS | GLP_BS | GLP_BS |
alpar@9 922 * | GLP_BS | GLP_NL | GLP_NU |
alpar@9 923 * | GLP_NL | GLP_BS | GLP_NL |
alpar@9 924 * | GLP_NL | GLP_NL | GLP_NL (*) |
alpar@9 925 * +-----------------+-----------------+------------------+
alpar@9 926 *
alpar@9 927 * Value of column q in solution to the original problem is the same as
alpar@9 928 * in solution to the transformed problem.
alpar@9 929 *
alpar@9 930 * 1. Formally, in solution to the transformed problem columns q and s
alpar@9 931 * cannot be non-basic at the same time, since the constraint (2)
alpar@9 932 * would be violated. However, if u[q] is close to zero, violation
alpar@9 933 * may be less than a working precision even if both columns q and s
alpar@9 934 * are non-basic. In this degenerate case row (2) can be only basic,
alpar@9 935 * i.e. non-active constraint (otherwise corresponding row of the
alpar@9 936 * basis matrix would be zero). This allows to pivot out auxiliary
alpar@9 937 * variable and pivot in column s, in which case the row becomes
alpar@9 938 * active while column s becomes basic.
alpar@9 939 *
alpar@9 940 * 2. If column q is integral, column s is also integral.
alpar@9 941 *
alpar@9 942 * RECOVERING INTERIOR-POINT SOLUTION
alpar@9 943 *
alpar@9 944 * Value of column q in solution to the original problem is the same as
alpar@9 945 * in solution to the transformed problem.
alpar@9 946 *
alpar@9 947 * RECOVERING MIP SOLUTION
alpar@9 948 *
alpar@9 949 * Value of column q in solution to the original problem is the same as
alpar@9 950 * in solution to the transformed problem. */
alpar@9 951
alpar@9 952 struct dbnd_col
alpar@9 953 { /* double-bounded column */
alpar@9 954 int q;
alpar@9 955 /* column reference number for variable x[q] */
alpar@9 956 int s;
alpar@9 957 /* column reference number for complement variable s */
alpar@9 958 };
alpar@9 959
alpar@9 960 static int rcv_dbnd_col(NPP *npp, void *info);
alpar@9 961
alpar@9 962 void npp_dbnd_col(NPP *npp, NPPCOL *q)
alpar@9 963 { /* process non-negative column with upper bound */
alpar@9 964 struct dbnd_col *info;
alpar@9 965 NPPROW *p;
alpar@9 966 NPPCOL *s;
alpar@9 967 /* the column must be non-negative with upper bound */
alpar@9 968 xassert(q->lb == 0.0);
alpar@9 969 xassert(q->ub > 0.0);
alpar@9 970 xassert(q->ub != +DBL_MAX);
alpar@9 971 /* create variable s */
alpar@9 972 s = npp_add_col(npp);
alpar@9 973 s->is_int = q->is_int;
alpar@9 974 s->lb = 0.0, s->ub = +DBL_MAX;
alpar@9 975 /* create equality constraint (2) */
alpar@9 976 p = npp_add_row(npp);
alpar@9 977 p->lb = p->ub = q->ub;
alpar@9 978 npp_add_aij(npp, p, q, +1.0);
alpar@9 979 npp_add_aij(npp, p, s, +1.0);
alpar@9 980 /* create transformation stack entry */
alpar@9 981 info = npp_push_tse(npp,
alpar@9 982 rcv_dbnd_col, sizeof(struct dbnd_col));
alpar@9 983 info->q = q->j;
alpar@9 984 info->s = s->j;
alpar@9 985 /* remove upper bound of x[q] */
alpar@9 986 q->ub = +DBL_MAX;
alpar@9 987 return;
alpar@9 988 }
alpar@9 989
alpar@9 990 static int rcv_dbnd_col(NPP *npp, void *_info)
alpar@9 991 { /* recover non-negative column with upper bound */
alpar@9 992 struct dbnd_col *info = _info;
alpar@9 993 if (npp->sol == GLP_BS)
alpar@9 994 { if (npp->c_stat[info->q] == GLP_BS)
alpar@9 995 { if (npp->c_stat[info->s] == GLP_BS)
alpar@9 996 npp->c_stat[info->q] = GLP_BS;
alpar@9 997 else if (npp->c_stat[info->s] == GLP_NL)
alpar@9 998 npp->c_stat[info->q] = GLP_NU;
alpar@9 999 else
alpar@9 1000 { npp_error();
alpar@9 1001 return 1;
alpar@9 1002 }
alpar@9 1003 }
alpar@9 1004 else if (npp->c_stat[info->q] == GLP_NL)
alpar@9 1005 { if (npp->c_stat[info->s] == GLP_BS ||
alpar@9 1006 npp->c_stat[info->s] == GLP_NL)
alpar@9 1007 npp->c_stat[info->q] = GLP_NL;
alpar@9 1008 else
alpar@9 1009 { npp_error();
alpar@9 1010 return 1;
alpar@9 1011 }
alpar@9 1012 }
alpar@9 1013 else
alpar@9 1014 { npp_error();
alpar@9 1015 return 1;
alpar@9 1016 }
alpar@9 1017 }
alpar@9 1018 return 0;
alpar@9 1019 }
alpar@9 1020
alpar@9 1021 /***********************************************************************
alpar@9 1022 * NAME
alpar@9 1023 *
alpar@9 1024 * npp_fixed_col - process fixed column
alpar@9 1025 *
alpar@9 1026 * SYNOPSIS
alpar@9 1027 *
alpar@9 1028 * #include "glpnpp.h"
alpar@9 1029 * void npp_fixed_col(NPP *npp, NPPCOL *q);
alpar@9 1030 *
alpar@9 1031 * DESCRIPTION
alpar@9 1032 *
alpar@9 1033 * The routine npp_fixed_col processes column q, which is fixed:
alpar@9 1034 *
alpar@9 1035 * x[q] = s[q], (1)
alpar@9 1036 *
alpar@9 1037 * where s[q] is a fixed column value.
alpar@9 1038 *
alpar@9 1039 * PROBLEM TRANSFORMATION
alpar@9 1040 *
alpar@9 1041 * The value of a fixed column can be substituted into the objective
alpar@9 1042 * and constraint rows that allows removing the column from the problem.
alpar@9 1043 *
alpar@9 1044 * Substituting x[q] = s[q] into the objective row, we have:
alpar@9 1045 *
alpar@9 1046 * z = sum c[j] x[j] + c0 =
alpar@9 1047 * j
alpar@9 1048 *
alpar@9 1049 * = sum c[j] x[j] + c[q] x[q] + c0 =
alpar@9 1050 * j!=q
alpar@9 1051 *
alpar@9 1052 * = sum c[j] x[j] + c[q] s[q] + c0 =
alpar@9 1053 * j!=q
alpar@9 1054 *
alpar@9 1055 * = sum c[j] x[j] + c~0,
alpar@9 1056 * j!=q
alpar@9 1057 *
alpar@9 1058 * where
alpar@9 1059 *
alpar@9 1060 * c~0 = c0 + c[q] s[q] (2)
alpar@9 1061 *
alpar@9 1062 * is the constant term of the objective in the transformed problem.
alpar@9 1063 * Similarly, substituting x[q] = s[q] into constraint row i, we have:
alpar@9 1064 *
alpar@9 1065 * L[i] <= sum a[i,j] x[j] <= U[i] ==>
alpar@9 1066 * j
alpar@9 1067 *
alpar@9 1068 * L[i] <= sum a[i,j] x[j] + a[i,q] x[q] <= U[i] ==>
alpar@9 1069 * j!=q
alpar@9 1070 *
alpar@9 1071 * L[i] <= sum a[i,j] x[j] + a[i,q] s[q] <= U[i] ==>
alpar@9 1072 * j!=q
alpar@9 1073 *
alpar@9 1074 * L~[i] <= sum a[i,j] x[j] + a[i,q] s <= U~[i],
alpar@9 1075 * j!=q
alpar@9 1076 *
alpar@9 1077 * where
alpar@9 1078 *
alpar@9 1079 * L~[i] = L[i] - a[i,q] s[q], U~[i] = U[i] - a[i,q] s[q] (3)
alpar@9 1080 *
alpar@9 1081 * are lower and upper bounds of row i in the transformed problem,
alpar@9 1082 * resp.
alpar@9 1083 *
alpar@9 1084 * RECOVERING BASIC SOLUTION
alpar@9 1085 *
alpar@9 1086 * Column q is assigned status GLP_NS and its value is assigned s[q].
alpar@9 1087 *
alpar@9 1088 * RECOVERING INTERIOR-POINT SOLUTION
alpar@9 1089 *
alpar@9 1090 * Value of column q is assigned s[q].
alpar@9 1091 *
alpar@9 1092 * RECOVERING MIP SOLUTION
alpar@9 1093 *
alpar@9 1094 * Value of column q is assigned s[q]. */
alpar@9 1095
alpar@9 1096 struct fixed_col
alpar@9 1097 { /* fixed column */
alpar@9 1098 int q;
alpar@9 1099 /* column reference number for variable x[q] */
alpar@9 1100 double s;
alpar@9 1101 /* value, at which x[q] is fixed */
alpar@9 1102 };
alpar@9 1103
alpar@9 1104 static int rcv_fixed_col(NPP *npp, void *info);
alpar@9 1105
alpar@9 1106 void npp_fixed_col(NPP *npp, NPPCOL *q)
alpar@9 1107 { /* process fixed column */
alpar@9 1108 struct fixed_col *info;
alpar@9 1109 NPPROW *i;
alpar@9 1110 NPPAIJ *aij;
alpar@9 1111 /* the column must be fixed */
alpar@9 1112 xassert(q->lb == q->ub);
alpar@9 1113 /* create transformation stack entry */
alpar@9 1114 info = npp_push_tse(npp,
alpar@9 1115 rcv_fixed_col, sizeof(struct fixed_col));
alpar@9 1116 info->q = q->j;
alpar@9 1117 info->s = q->lb;
alpar@9 1118 /* substitute x[q] = s[q] into objective row */
alpar@9 1119 npp->c0 += q->coef * q->lb;
alpar@9 1120 /* substitute x[q] = s[q] into constraint rows */
alpar@9 1121 for (aij = q->ptr; aij != NULL; aij = aij->c_next)
alpar@9 1122 { i = aij->row;
alpar@9 1123 if (i->lb == i->ub)
alpar@9 1124 i->ub = (i->lb -= aij->val * q->lb);
alpar@9 1125 else
alpar@9 1126 { if (i->lb != -DBL_MAX)
alpar@9 1127 i->lb -= aij->val * q->lb;
alpar@9 1128 if (i->ub != +DBL_MAX)
alpar@9 1129 i->ub -= aij->val * q->lb;
alpar@9 1130 }
alpar@9 1131 }
alpar@9 1132 /* remove the column from the problem */
alpar@9 1133 npp_del_col(npp, q);
alpar@9 1134 return;
alpar@9 1135 }
alpar@9 1136
alpar@9 1137 static int rcv_fixed_col(NPP *npp, void *_info)
alpar@9 1138 { /* recover fixed column */
alpar@9 1139 struct fixed_col *info = _info;
alpar@9 1140 if (npp->sol == GLP_SOL)
alpar@9 1141 npp->c_stat[info->q] = GLP_NS;
alpar@9 1142 npp->c_value[info->q] = info->s;
alpar@9 1143 return 0;
alpar@9 1144 }
alpar@9 1145
alpar@9 1146 /***********************************************************************
alpar@9 1147 * NAME
alpar@9 1148 *
alpar@9 1149 * npp_make_equality - process row with almost identical bounds
alpar@9 1150 *
alpar@9 1151 * SYNOPSIS
alpar@9 1152 *
alpar@9 1153 * #include "glpnpp.h"
alpar@9 1154 * int npp_make_equality(NPP *npp, NPPROW *p);
alpar@9 1155 *
alpar@9 1156 * DESCRIPTION
alpar@9 1157 *
alpar@9 1158 * The routine npp_make_equality processes row p:
alpar@9 1159 *
alpar@9 1160 * L[p] <= sum a[p,j] x[j] <= U[p], (1)
alpar@9 1161 * j
alpar@9 1162 *
alpar@9 1163 * where -oo < L[p] < U[p] < +oo, i.e. which is double-sided inequality
alpar@9 1164 * constraint.
alpar@9 1165 *
alpar@9 1166 * RETURNS
alpar@9 1167 *
alpar@9 1168 * 0 - row bounds have not been changed;
alpar@9 1169 *
alpar@9 1170 * 1 - row has been replaced by equality constraint.
alpar@9 1171 *
alpar@9 1172 * PROBLEM TRANSFORMATION
alpar@9 1173 *
alpar@9 1174 * If bounds of row (1) are very close to each other:
alpar@9 1175 *
alpar@9 1176 * U[p] - L[p] <= eps, (2)
alpar@9 1177 *
alpar@9 1178 * where eps is an absolute tolerance for row value, the row can be
alpar@9 1179 * replaced by the following almost equivalent equiality constraint:
alpar@9 1180 *
alpar@9 1181 * sum a[p,j] x[j] = b, (3)
alpar@9 1182 * j
alpar@9 1183 *
alpar@9 1184 * where b = (L[p] + U[p]) / 2. If the right-hand side in (3) happens
alpar@9 1185 * to be very close to its nearest integer:
alpar@9 1186 *
alpar@9 1187 * |b - floor(b + 0.5)| <= eps, (4)
alpar@9 1188 *
alpar@9 1189 * it is reasonable to use this nearest integer as the right-hand side.
alpar@9 1190 *
alpar@9 1191 * RECOVERING BASIC SOLUTION
alpar@9 1192 *
alpar@9 1193 * Status of row p in solution to the original problem is determined
alpar@9 1194 * by its status and the sign of its multiplier pi[p] in solution to
alpar@9 1195 * the transformed problem as follows:
alpar@9 1196 *
alpar@9 1197 * +-----------------------+---------+--------------------+
alpar@9 1198 * | Status of row p | Sign of | Status of row p |
alpar@9 1199 * | (transformed problem) | pi[p] | (original problem) |
alpar@9 1200 * +-----------------------+---------+--------------------+
alpar@9 1201 * | GLP_BS | + / - | GLP_BS |
alpar@9 1202 * | GLP_NS | + | GLP_NL |
alpar@9 1203 * | GLP_NS | - | GLP_NU |
alpar@9 1204 * +-----------------------+---------+--------------------+
alpar@9 1205 *
alpar@9 1206 * Value of row multiplier pi[p] in solution to the original problem is
alpar@9 1207 * the same as in solution to the transformed problem.
alpar@9 1208 *
alpar@9 1209 * RECOVERING INTERIOR POINT SOLUTION
alpar@9 1210 *
alpar@9 1211 * Value of row multiplier pi[p] in solution to the original problem is
alpar@9 1212 * the same as in solution to the transformed problem.
alpar@9 1213 *
alpar@9 1214 * RECOVERING MIP SOLUTION
alpar@9 1215 *
alpar@9 1216 * None needed. */
alpar@9 1217
alpar@9 1218 struct make_equality
alpar@9 1219 { /* row with almost identical bounds */
alpar@9 1220 int p;
alpar@9 1221 /* row reference number */
alpar@9 1222 };
alpar@9 1223
alpar@9 1224 static int rcv_make_equality(NPP *npp, void *info);
alpar@9 1225
alpar@9 1226 int npp_make_equality(NPP *npp, NPPROW *p)
alpar@9 1227 { /* process row with almost identical bounds */
alpar@9 1228 struct make_equality *info;
alpar@9 1229 double b, eps, nint;
alpar@9 1230 /* the row must be double-sided inequality */
alpar@9 1231 xassert(p->lb != -DBL_MAX);
alpar@9 1232 xassert(p->ub != +DBL_MAX);
alpar@9 1233 xassert(p->lb < p->ub);
alpar@9 1234 /* check row bounds */
alpar@9 1235 eps = 1e-9 + 1e-12 * fabs(p->lb);
alpar@9 1236 if (p->ub - p->lb > eps) return 0;
alpar@9 1237 /* row bounds are very close to each other */
alpar@9 1238 /* create transformation stack entry */
alpar@9 1239 info = npp_push_tse(npp,
alpar@9 1240 rcv_make_equality, sizeof(struct make_equality));
alpar@9 1241 info->p = p->i;
alpar@9 1242 /* compute right-hand side */
alpar@9 1243 b = 0.5 * (p->ub + p->lb);
alpar@9 1244 nint = floor(b + 0.5);
alpar@9 1245 if (fabs(b - nint) <= eps) b = nint;
alpar@9 1246 /* replace row p by almost equivalent equality constraint */
alpar@9 1247 p->lb = p->ub = b;
alpar@9 1248 return 1;
alpar@9 1249 }
alpar@9 1250
alpar@9 1251 int rcv_make_equality(NPP *npp, void *_info)
alpar@9 1252 { /* recover row with almost identical bounds */
alpar@9 1253 struct make_equality *info = _info;
alpar@9 1254 if (npp->sol == GLP_SOL)
alpar@9 1255 { if (npp->r_stat[info->p] == GLP_BS)
alpar@9 1256 npp->r_stat[info->p] = GLP_BS;
alpar@9 1257 else if (npp->r_stat[info->p] == GLP_NS)
alpar@9 1258 { if (npp->r_pi[info->p] >= 0.0)
alpar@9 1259 npp->r_stat[info->p] = GLP_NL;
alpar@9 1260 else
alpar@9 1261 npp->r_stat[info->p] = GLP_NU;
alpar@9 1262 }
alpar@9 1263 else
alpar@9 1264 { npp_error();
alpar@9 1265 return 1;
alpar@9 1266 }
alpar@9 1267 }
alpar@9 1268 return 0;
alpar@9 1269 }
alpar@9 1270
alpar@9 1271 /***********************************************************************
alpar@9 1272 * NAME
alpar@9 1273 *
alpar@9 1274 * npp_make_fixed - process column with almost identical bounds
alpar@9 1275 *
alpar@9 1276 * SYNOPSIS
alpar@9 1277 *
alpar@9 1278 * #include "glpnpp.h"
alpar@9 1279 * int npp_make_fixed(NPP *npp, NPPCOL *q);
alpar@9 1280 *
alpar@9 1281 * DESCRIPTION
alpar@9 1282 *
alpar@9 1283 * The routine npp_make_fixed processes column q:
alpar@9 1284 *
alpar@9 1285 * l[q] <= x[q] <= u[q], (1)
alpar@9 1286 *
alpar@9 1287 * where -oo < l[q] < u[q] < +oo, i.e. which has both lower and upper
alpar@9 1288 * bounds.
alpar@9 1289 *
alpar@9 1290 * RETURNS
alpar@9 1291 *
alpar@9 1292 * 0 - column bounds have not been changed;
alpar@9 1293 *
alpar@9 1294 * 1 - column has been fixed.
alpar@9 1295 *
alpar@9 1296 * PROBLEM TRANSFORMATION
alpar@9 1297 *
alpar@9 1298 * If bounds of column (1) are very close to each other:
alpar@9 1299 *
alpar@9 1300 * u[q] - l[q] <= eps, (2)
alpar@9 1301 *
alpar@9 1302 * where eps is an absolute tolerance for column value, the column can
alpar@9 1303 * be fixed:
alpar@9 1304 *
alpar@9 1305 * x[q] = s[q], (3)
alpar@9 1306 *
alpar@9 1307 * where s[q] = (l[q] + u[q]) / 2. And if the fixed column value s[q]
alpar@9 1308 * happens to be very close to its nearest integer:
alpar@9 1309 *
alpar@9 1310 * |s[q] - floor(s[q] + 0.5)| <= eps, (4)
alpar@9 1311 *
alpar@9 1312 * it is reasonable to use this nearest integer as the fixed value.
alpar@9 1313 *
alpar@9 1314 * RECOVERING BASIC SOLUTION
alpar@9 1315 *
alpar@9 1316 * In the dual system of the original (as well as transformed) problem
alpar@9 1317 * column q corresponds to the following row:
alpar@9 1318 *
alpar@9 1319 * sum a[i,q] pi[i] + lambda[q] = c[q]. (5)
alpar@9 1320 * i
alpar@9 1321 *
alpar@9 1322 * Since multipliers pi[i] are known for all rows from solution to the
alpar@9 1323 * transformed problem, formula (5) allows computing value of multiplier
alpar@9 1324 * (reduced cost) for column q:
alpar@9 1325 *
alpar@9 1326 * lambda[q] = c[q] - sum a[i,q] pi[i]. (6)
alpar@9 1327 * i
alpar@9 1328 *
alpar@9 1329 * Status of column q in solution to the original problem is determined
alpar@9 1330 * by its status and the sign of its multiplier lambda[q] in solution to
alpar@9 1331 * the transformed problem as follows:
alpar@9 1332 *
alpar@9 1333 * +-----------------------+-----------+--------------------+
alpar@9 1334 * | Status of column q | Sign of | Status of column q |
alpar@9 1335 * | (transformed problem) | lambda[q] | (original problem) |
alpar@9 1336 * +-----------------------+-----------+--------------------+
alpar@9 1337 * | GLP_BS | + / - | GLP_BS |
alpar@9 1338 * | GLP_NS | + | GLP_NL |
alpar@9 1339 * | GLP_NS | - | GLP_NU |
alpar@9 1340 * +-----------------------+-----------+--------------------+
alpar@9 1341 *
alpar@9 1342 * Value of column q in solution to the original problem is the same as
alpar@9 1343 * in solution to the transformed problem.
alpar@9 1344 *
alpar@9 1345 * RECOVERING INTERIOR POINT SOLUTION
alpar@9 1346 *
alpar@9 1347 * Value of column q in solution to the original problem is the same as
alpar@9 1348 * in solution to the transformed problem.
alpar@9 1349 *
alpar@9 1350 * RECOVERING MIP SOLUTION
alpar@9 1351 *
alpar@9 1352 * None needed. */
alpar@9 1353
alpar@9 1354 struct make_fixed
alpar@9 1355 { /* column with almost identical bounds */
alpar@9 1356 int q;
alpar@9 1357 /* column reference number */
alpar@9 1358 double c;
alpar@9 1359 /* objective coefficient at x[q] */
alpar@9 1360 NPPLFE *ptr;
alpar@9 1361 /* list of non-zero coefficients a[i,q] */
alpar@9 1362 };
alpar@9 1363
alpar@9 1364 static int rcv_make_fixed(NPP *npp, void *info);
alpar@9 1365
alpar@9 1366 int npp_make_fixed(NPP *npp, NPPCOL *q)
alpar@9 1367 { /* process column with almost identical bounds */
alpar@9 1368 struct make_fixed *info;
alpar@9 1369 NPPAIJ *aij;
alpar@9 1370 NPPLFE *lfe;
alpar@9 1371 double s, eps, nint;
alpar@9 1372 /* the column must be double-bounded */
alpar@9 1373 xassert(q->lb != -DBL_MAX);
alpar@9 1374 xassert(q->ub != +DBL_MAX);
alpar@9 1375 xassert(q->lb < q->ub);
alpar@9 1376 /* check column bounds */
alpar@9 1377 eps = 1e-9 + 1e-12 * fabs(q->lb);
alpar@9 1378 if (q->ub - q->lb > eps) return 0;
alpar@9 1379 /* column bounds are very close to each other */
alpar@9 1380 /* create transformation stack entry */
alpar@9 1381 info = npp_push_tse(npp,
alpar@9 1382 rcv_make_fixed, sizeof(struct make_fixed));
alpar@9 1383 info->q = q->j;
alpar@9 1384 info->c = q->coef;
alpar@9 1385 info->ptr = NULL;
alpar@9 1386 /* save column coefficients a[i,q] (needed for basic solution
alpar@9 1387 only) */
alpar@9 1388 if (npp->sol == GLP_SOL)
alpar@9 1389 { for (aij = q->ptr; aij != NULL; aij = aij->c_next)
alpar@9 1390 { lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE));
alpar@9 1391 lfe->ref = aij->row->i;
alpar@9 1392 lfe->val = aij->val;
alpar@9 1393 lfe->next = info->ptr;
alpar@9 1394 info->ptr = lfe;
alpar@9 1395 }
alpar@9 1396 }
alpar@9 1397 /* compute column fixed value */
alpar@9 1398 s = 0.5 * (q->ub + q->lb);
alpar@9 1399 nint = floor(s + 0.5);
alpar@9 1400 if (fabs(s - nint) <= eps) s = nint;
alpar@9 1401 /* make column q fixed */
alpar@9 1402 q->lb = q->ub = s;
alpar@9 1403 return 1;
alpar@9 1404 }
alpar@9 1405
alpar@9 1406 static int rcv_make_fixed(NPP *npp, void *_info)
alpar@9 1407 { /* recover column with almost identical bounds */
alpar@9 1408 struct make_fixed *info = _info;
alpar@9 1409 NPPLFE *lfe;
alpar@9 1410 double lambda;
alpar@9 1411 if (npp->sol == GLP_SOL)
alpar@9 1412 { if (npp->c_stat[info->q] == GLP_BS)
alpar@9 1413 npp->c_stat[info->q] = GLP_BS;
alpar@9 1414 else if (npp->c_stat[info->q] == GLP_NS)
alpar@9 1415 { /* compute multiplier for column q with formula (6) */
alpar@9 1416 lambda = info->c;
alpar@9 1417 for (lfe = info->ptr; lfe != NULL; lfe = lfe->next)
alpar@9 1418 lambda -= lfe->val * npp->r_pi[lfe->ref];
alpar@9 1419 /* assign status to non-basic column */
alpar@9 1420 if (lambda >= 0.0)
alpar@9 1421 npp->c_stat[info->q] = GLP_NL;
alpar@9 1422 else
alpar@9 1423 npp->c_stat[info->q] = GLP_NU;
alpar@9 1424 }
alpar@9 1425 else
alpar@9 1426 { npp_error();
alpar@9 1427 return 1;
alpar@9 1428 }
alpar@9 1429 }
alpar@9 1430 return 0;
alpar@9 1431 }
alpar@9 1432
alpar@9 1433 /* eof */