lemon-project-template-glpk

annotate deps/glpk/src/glpnpp04.c @ 11:4fc6ad2fb8a6

find changesets by author, revision, files, or words in the commit message
Test GLPK in src/main.cc
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 21:43:29 +0100
parents
children
rev   line source
alpar@9 1 /* glpnpp04.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_binarize_prob - binarize MIP problem
alpar@9 31 *
alpar@9 32 * SYNOPSIS
alpar@9 33 *
alpar@9 34 * #include "glpnpp.h"
alpar@9 35 * int npp_binarize_prob(NPP *npp);
alpar@9 36 *
alpar@9 37 * DESCRIPTION
alpar@9 38 *
alpar@9 39 * The routine npp_binarize_prob replaces in the original MIP problem
alpar@9 40 * every integer variable:
alpar@9 41 *
alpar@9 42 * l[q] <= x[q] <= u[q], (1)
alpar@9 43 *
alpar@9 44 * where l[q] < u[q], by an equivalent sum of binary variables.
alpar@9 45 *
alpar@9 46 * RETURNS
alpar@9 47 *
alpar@9 48 * The routine returns the number of integer variables for which the
alpar@9 49 * transformation failed, because u[q] - l[q] > d_max.
alpar@9 50 *
alpar@9 51 * PROBLEM TRANSFORMATION
alpar@9 52 *
alpar@9 53 * If variable x[q] has non-zero lower bound, it is first processed
alpar@9 54 * with the routine npp_lbnd_col. Thus, we can assume that:
alpar@9 55 *
alpar@9 56 * 0 <= x[q] <= u[q]. (2)
alpar@9 57 *
alpar@9 58 * If u[q] = 1, variable x[q] is already binary, so further processing
alpar@9 59 * is not needed. Let, therefore, that 2 <= u[q] <= d_max, and n be a
alpar@9 60 * smallest integer such that u[q] <= 2^n - 1 (n >= 2, since u[q] >= 2).
alpar@9 61 * Then variable x[q] can be replaced by the following sum:
alpar@9 62 *
alpar@9 63 * n-1
alpar@9 64 * x[q] = sum 2^k x[k], (3)
alpar@9 65 * k=0
alpar@9 66 *
alpar@9 67 * where x[k] are binary columns (variables). If u[q] < 2^n - 1, the
alpar@9 68 * following additional inequality constraint must be also included in
alpar@9 69 * the transformed problem:
alpar@9 70 *
alpar@9 71 * n-1
alpar@9 72 * sum 2^k x[k] <= u[q]. (4)
alpar@9 73 * k=0
alpar@9 74 *
alpar@9 75 * Note: Assuming that in the transformed problem x[q] becomes binary
alpar@9 76 * variable x[0], this transformation causes new n-1 binary variables
alpar@9 77 * to appear.
alpar@9 78 *
alpar@9 79 * Substituting x[q] from (3) to the objective row gives:
alpar@9 80 *
alpar@9 81 * z = sum c[j] x[j] + c[0] =
alpar@9 82 * j
alpar@9 83 *
alpar@9 84 * = sum c[j] x[j] + c[q] x[q] + c[0] =
alpar@9 85 * j!=q
alpar@9 86 * n-1
alpar@9 87 * = sum c[j] x[j] + c[q] sum 2^k x[k] + c[0] =
alpar@9 88 * j!=q k=0
alpar@9 89 * n-1
alpar@9 90 * = sum c[j] x[j] + sum c[k] x[k] + c[0],
alpar@9 91 * j!=q k=0
alpar@9 92 *
alpar@9 93 * where:
alpar@9 94 *
alpar@9 95 * c[k] = 2^k c[q], k = 0, ..., n-1. (5)
alpar@9 96 *
alpar@9 97 * And substituting x[q] from (3) to i-th constraint row i gives:
alpar@9 98 *
alpar@9 99 * L[i] <= sum a[i,j] x[j] <= U[i] ==>
alpar@9 100 * j
alpar@9 101 *
alpar@9 102 * L[i] <= sum a[i,j] x[j] + a[i,q] x[q] <= U[i] ==>
alpar@9 103 * j!=q
alpar@9 104 * n-1
alpar@9 105 * L[i] <= sum a[i,j] x[j] + a[i,q] sum 2^k x[k] <= U[i] ==>
alpar@9 106 * j!=q k=0
alpar@9 107 * n-1
alpar@9 108 * L[i] <= sum a[i,j] x[j] + sum a[i,k] x[k] <= U[i],
alpar@9 109 * j!=q k=0
alpar@9 110 *
alpar@9 111 * where:
alpar@9 112 *
alpar@9 113 * a[i,k] = 2^k a[i,q], k = 0, ..., n-1. (6)
alpar@9 114 *
alpar@9 115 * RECOVERING SOLUTION
alpar@9 116 *
alpar@9 117 * Value of variable x[q] is computed with formula (3). */
alpar@9 118
alpar@9 119 struct binarize
alpar@9 120 { int q;
alpar@9 121 /* column reference number for x[q] = x[0] */
alpar@9 122 int j;
alpar@9 123 /* column reference number for x[1]; x[2] has reference number
alpar@9 124 j+1, x[3] - j+2, etc. */
alpar@9 125 int n;
alpar@9 126 /* total number of binary variables, n >= 2 */
alpar@9 127 };
alpar@9 128
alpar@9 129 static int rcv_binarize_prob(NPP *npp, void *info);
alpar@9 130
alpar@9 131 int npp_binarize_prob(NPP *npp)
alpar@9 132 { /* binarize MIP problem */
alpar@9 133 struct binarize *info;
alpar@9 134 NPPROW *row;
alpar@9 135 NPPCOL *col, *bin;
alpar@9 136 NPPAIJ *aij;
alpar@9 137 int u, n, k, temp, nfails, nvars, nbins, nrows;
alpar@9 138 /* new variables will be added to the end of the column list, so
alpar@9 139 we go from the end to beginning of the column list */
alpar@9 140 nfails = nvars = nbins = nrows = 0;
alpar@9 141 for (col = npp->c_tail; col != NULL; col = col->prev)
alpar@9 142 { /* skip continuous variable */
alpar@9 143 if (!col->is_int) continue;
alpar@9 144 /* skip fixed variable */
alpar@9 145 if (col->lb == col->ub) continue;
alpar@9 146 /* skip binary variable */
alpar@9 147 if (col->lb == 0.0 && col->ub == 1.0) continue;
alpar@9 148 /* check if the transformation is applicable */
alpar@9 149 if (col->lb < -1e6 || col->ub > +1e6 ||
alpar@9 150 col->ub - col->lb > 4095.0)
alpar@9 151 { /* unfortunately, not */
alpar@9 152 nfails++;
alpar@9 153 continue;
alpar@9 154 }
alpar@9 155 /* process integer non-binary variable x[q] */
alpar@9 156 nvars++;
alpar@9 157 /* make x[q] non-negative, if its lower bound is non-zero */
alpar@9 158 if (col->lb != 0.0)
alpar@9 159 npp_lbnd_col(npp, col);
alpar@9 160 /* now 0 <= x[q] <= u[q] */
alpar@9 161 xassert(col->lb == 0.0);
alpar@9 162 u = (int)col->ub;
alpar@9 163 xassert(col->ub == (double)u);
alpar@9 164 /* if x[q] is binary, further processing is not needed */
alpar@9 165 if (u == 1) continue;
alpar@9 166 /* determine smallest n such that u <= 2^n - 1 (thus, n is the
alpar@9 167 number of binary variables needed) */
alpar@9 168 n = 2, temp = 4;
alpar@9 169 while (u >= temp)
alpar@9 170 n++, temp += temp;
alpar@9 171 nbins += n;
alpar@9 172 /* create transformation stack entry */
alpar@9 173 info = npp_push_tse(npp,
alpar@9 174 rcv_binarize_prob, sizeof(struct binarize));
alpar@9 175 info->q = col->j;
alpar@9 176 info->j = 0; /* will be set below */
alpar@9 177 info->n = n;
alpar@9 178 /* if u < 2^n - 1, we need one additional row for (4) */
alpar@9 179 if (u < temp - 1)
alpar@9 180 { row = npp_add_row(npp), nrows++;
alpar@9 181 row->lb = -DBL_MAX, row->ub = u;
alpar@9 182 }
alpar@9 183 else
alpar@9 184 row = NULL;
alpar@9 185 /* in the transformed problem variable x[q] becomes binary
alpar@9 186 variable x[0], so its objective and constraint coefficients
alpar@9 187 are not changed */
alpar@9 188 col->ub = 1.0;
alpar@9 189 /* include x[0] into constraint (4) */
alpar@9 190 if (row != NULL)
alpar@9 191 npp_add_aij(npp, row, col, 1.0);
alpar@9 192 /* add other binary variables x[1], ..., x[n-1] */
alpar@9 193 for (k = 1, temp = 2; k < n; k++, temp += temp)
alpar@9 194 { /* add new binary variable x[k] */
alpar@9 195 bin = npp_add_col(npp);
alpar@9 196 bin->is_int = 1;
alpar@9 197 bin->lb = 0.0, bin->ub = 1.0;
alpar@9 198 bin->coef = (double)temp * col->coef;
alpar@9 199 /* store column reference number for x[1] */
alpar@9 200 if (info->j == 0)
alpar@9 201 info->j = bin->j;
alpar@9 202 else
alpar@9 203 xassert(info->j + (k-1) == bin->j);
alpar@9 204 /* duplicate constraint coefficients for x[k]; this also
alpar@9 205 automatically includes x[k] into constraint (4) */
alpar@9 206 for (aij = col->ptr; aij != NULL; aij = aij->c_next)
alpar@9 207 npp_add_aij(npp, aij->row, bin, (double)temp * aij->val);
alpar@9 208 }
alpar@9 209 }
alpar@9 210 if (nvars > 0)
alpar@9 211 xprintf("%d integer variable(s) were replaced by %d binary one"
alpar@9 212 "s\n", nvars, nbins);
alpar@9 213 if (nrows > 0)
alpar@9 214 xprintf("%d row(s) were added due to binarization\n", nrows);
alpar@9 215 if (nfails > 0)
alpar@9 216 xprintf("Binarization failed for %d integer variable(s)\n",
alpar@9 217 nfails);
alpar@9 218 return nfails;
alpar@9 219 }
alpar@9 220
alpar@9 221 static int rcv_binarize_prob(NPP *npp, void *_info)
alpar@9 222 { /* recovery binarized variable */
alpar@9 223 struct binarize *info = _info;
alpar@9 224 int k, temp;
alpar@9 225 double sum;
alpar@9 226 /* compute value of x[q]; see formula (3) */
alpar@9 227 sum = npp->c_value[info->q];
alpar@9 228 for (k = 1, temp = 2; k < info->n; k++, temp += temp)
alpar@9 229 sum += (double)temp * npp->c_value[info->j + (k-1)];
alpar@9 230 npp->c_value[info->q] = sum;
alpar@9 231 return 0;
alpar@9 232 }
alpar@9 233
alpar@9 234 /**********************************************************************/
alpar@9 235
alpar@9 236 struct elem
alpar@9 237 { /* linear form element a[j] x[j] */
alpar@9 238 double aj;
alpar@9 239 /* non-zero coefficient value */
alpar@9 240 NPPCOL *xj;
alpar@9 241 /* pointer to variable (column) */
alpar@9 242 struct elem *next;
alpar@9 243 /* pointer to another term */
alpar@9 244 };
alpar@9 245
alpar@9 246 static struct elem *copy_form(NPP *npp, NPPROW *row, double s)
alpar@9 247 { /* copy linear form */
alpar@9 248 NPPAIJ *aij;
alpar@9 249 struct elem *ptr, *e;
alpar@9 250 ptr = NULL;
alpar@9 251 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
alpar@9 252 { e = dmp_get_atom(npp->pool, sizeof(struct elem));
alpar@9 253 e->aj = s * aij->val;
alpar@9 254 e->xj = aij->col;
alpar@9 255 e->next = ptr;
alpar@9 256 ptr = e;
alpar@9 257 }
alpar@9 258 return ptr;
alpar@9 259 }
alpar@9 260
alpar@9 261 static void drop_form(NPP *npp, struct elem *ptr)
alpar@9 262 { /* drop linear form */
alpar@9 263 struct elem *e;
alpar@9 264 while (ptr != NULL)
alpar@9 265 { e = ptr;
alpar@9 266 ptr = e->next;
alpar@9 267 dmp_free_atom(npp->pool, e, sizeof(struct elem));
alpar@9 268 }
alpar@9 269 return;
alpar@9 270 }
alpar@9 271
alpar@9 272 /***********************************************************************
alpar@9 273 * NAME
alpar@9 274 *
alpar@9 275 * npp_is_packing - test if constraint is packing inequality
alpar@9 276 *
alpar@9 277 * SYNOPSIS
alpar@9 278 *
alpar@9 279 * #include "glpnpp.h"
alpar@9 280 * int npp_is_packing(NPP *npp, NPPROW *row);
alpar@9 281 *
alpar@9 282 * RETURNS
alpar@9 283 *
alpar@9 284 * If the specified row (constraint) is packing inequality (see below),
alpar@9 285 * the routine npp_is_packing returns non-zero. Otherwise, it returns
alpar@9 286 * zero.
alpar@9 287 *
alpar@9 288 * PACKING INEQUALITIES
alpar@9 289 *
alpar@9 290 * In canonical format the packing inequality is the following:
alpar@9 291 *
alpar@9 292 * sum x[j] <= 1, (1)
alpar@9 293 * j in J
alpar@9 294 *
alpar@9 295 * where all variables x[j] are binary. This inequality expresses the
alpar@9 296 * condition that in any integer feasible solution at most one variable
alpar@9 297 * from set J can take non-zero (unity) value while other variables
alpar@9 298 * must be equal to zero. W.l.o.g. it is assumed that |J| >= 2, because
alpar@9 299 * if J is empty or |J| = 1, the inequality (1) is redundant.
alpar@9 300 *
alpar@9 301 * In general case the packing inequality may include original variables
alpar@9 302 * x[j] as well as their complements x~[j]:
alpar@9 303 *
alpar@9 304 * sum x[j] + sum x~[j] <= 1, (2)
alpar@9 305 * j in Jp j in Jn
alpar@9 306 *
alpar@9 307 * where Jp and Jn are not intersected. Therefore, using substitution
alpar@9 308 * x~[j] = 1 - x[j] gives the packing inequality in generalized format:
alpar@9 309 *
alpar@9 310 * sum x[j] - sum x[j] <= 1 - |Jn|. (3)
alpar@9 311 * j in Jp j in Jn */
alpar@9 312
alpar@9 313 int npp_is_packing(NPP *npp, NPPROW *row)
alpar@9 314 { /* test if constraint is packing inequality */
alpar@9 315 NPPCOL *col;
alpar@9 316 NPPAIJ *aij;
alpar@9 317 int b;
alpar@9 318 xassert(npp == npp);
alpar@9 319 if (!(row->lb == -DBL_MAX && row->ub != +DBL_MAX))
alpar@9 320 return 0;
alpar@9 321 b = 1;
alpar@9 322 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
alpar@9 323 { col = aij->col;
alpar@9 324 if (!(col->is_int && col->lb == 0.0 && col->ub == 1.0))
alpar@9 325 return 0;
alpar@9 326 if (aij->val == +1.0)
alpar@9 327 ;
alpar@9 328 else if (aij->val == -1.0)
alpar@9 329 b--;
alpar@9 330 else
alpar@9 331 return 0;
alpar@9 332 }
alpar@9 333 if (row->ub != (double)b) return 0;
alpar@9 334 return 1;
alpar@9 335 }
alpar@9 336
alpar@9 337 /***********************************************************************
alpar@9 338 * NAME
alpar@9 339 *
alpar@9 340 * npp_hidden_packing - identify hidden packing inequality
alpar@9 341 *
alpar@9 342 * SYNOPSIS
alpar@9 343 *
alpar@9 344 * #include "glpnpp.h"
alpar@9 345 * int npp_hidden_packing(NPP *npp, NPPROW *row);
alpar@9 346 *
alpar@9 347 * DESCRIPTION
alpar@9 348 *
alpar@9 349 * The routine npp_hidden_packing processes specified inequality
alpar@9 350 * constraint, which includes only binary variables, and the number of
alpar@9 351 * the variables is not less than two. If the original inequality is
alpar@9 352 * equivalent to a packing inequality, the routine replaces it by this
alpar@9 353 * equivalent inequality. If the original constraint is double-sided
alpar@9 354 * inequality, it is replaced by a pair of single-sided inequalities,
alpar@9 355 * if necessary.
alpar@9 356 *
alpar@9 357 * RETURNS
alpar@9 358 *
alpar@9 359 * If the original inequality constraint was replaced by equivalent
alpar@9 360 * packing inequality, the routine npp_hidden_packing returns non-zero.
alpar@9 361 * Otherwise, it returns zero.
alpar@9 362 *
alpar@9 363 * PROBLEM TRANSFORMATION
alpar@9 364 *
alpar@9 365 * Consider an inequality constraint:
alpar@9 366 *
alpar@9 367 * sum a[j] x[j] <= b, (1)
alpar@9 368 * j in J
alpar@9 369 *
alpar@9 370 * where all variables x[j] are binary, and |J| >= 2. (In case of '>='
alpar@9 371 * inequality it can be transformed to '<=' format by multiplying both
alpar@9 372 * its sides by -1.)
alpar@9 373 *
alpar@9 374 * Let Jp = {j: a[j] > 0}, Jn = {j: a[j] < 0}. Performing substitution
alpar@9 375 * x[j] = 1 - x~[j] for all j in Jn, we have:
alpar@9 376 *
alpar@9 377 * sum a[j] x[j] <= b ==>
alpar@9 378 * j in J
alpar@9 379 *
alpar@9 380 * sum a[j] x[j] + sum a[j] x[j] <= b ==>
alpar@9 381 * j in Jp j in Jn
alpar@9 382 *
alpar@9 383 * sum a[j] x[j] + sum a[j] (1 - x~[j]) <= b ==>
alpar@9 384 * j in Jp j in Jn
alpar@9 385 *
alpar@9 386 * sum a[j] x[j] - sum a[j] x~[j] <= b - sum a[j].
alpar@9 387 * j in Jp j in Jn j in Jn
alpar@9 388 *
alpar@9 389 * Thus, meaning the transformation above, we can assume that in
alpar@9 390 * inequality (1) all coefficients a[j] are positive. Moreover, we can
alpar@9 391 * assume that a[j] <= b. In fact, let a[j] > b; then the following
alpar@9 392 * three cases are possible:
alpar@9 393 *
alpar@9 394 * 1) b < 0. In this case inequality (1) is infeasible, so the problem
alpar@9 395 * has no feasible solution (see the routine npp_analyze_row);
alpar@9 396 *
alpar@9 397 * 2) b = 0. In this case inequality (1) is a forcing inequality on its
alpar@9 398 * upper bound (see the routine npp_forcing row), from which it
alpar@9 399 * follows that all variables x[j] should be fixed at zero;
alpar@9 400 *
alpar@9 401 * 3) b > 0. In this case inequality (1) defines an implied zero upper
alpar@9 402 * bound for variable x[j] (see the routine npp_implied_bounds), from
alpar@9 403 * which it follows that x[j] should be fixed at zero.
alpar@9 404 *
alpar@9 405 * It is assumed that all three cases listed above have been recognized
alpar@9 406 * by the routine npp_process_prob, which performs basic MIP processing
alpar@9 407 * prior to a call the routine npp_hidden_packing. So, if one of these
alpar@9 408 * cases occurs, we should just skip processing such constraint.
alpar@9 409 *
alpar@9 410 * Thus, let 0 < a[j] <= b. Then it is obvious that constraint (1) is
alpar@9 411 * equivalent to packing inquality only if:
alpar@9 412 *
alpar@9 413 * a[j] + a[k] > b + eps (2)
alpar@9 414 *
alpar@9 415 * for all j, k in J, j != k, where eps is an absolute tolerance for
alpar@9 416 * row (linear form) value. Checking the condition (2) for all j and k,
alpar@9 417 * j != k, requires time O(|J|^2). However, this time can be reduced to
alpar@9 418 * O(|J|), if use minimal a[j] and a[k], in which case it is sufficient
alpar@9 419 * to check the condition (2) only once.
alpar@9 420 *
alpar@9 421 * Once the original inequality (1) is replaced by equivalent packing
alpar@9 422 * inequality, we need to perform back substitution x~[j] = 1 - x[j] for
alpar@9 423 * all j in Jn (see above).
alpar@9 424 *
alpar@9 425 * RECOVERING SOLUTION
alpar@9 426 *
alpar@9 427 * None needed. */
alpar@9 428
alpar@9 429 static int hidden_packing(NPP *npp, struct elem *ptr, double *_b)
alpar@9 430 { /* process inequality constraint: sum a[j] x[j] <= b;
alpar@9 431 0 - specified row is NOT hidden packing inequality;
alpar@9 432 1 - specified row is packing inequality;
alpar@9 433 2 - specified row is hidden packing inequality. */
alpar@9 434 struct elem *e, *ej, *ek;
alpar@9 435 int neg;
alpar@9 436 double b = *_b, eps;
alpar@9 437 xassert(npp == npp);
alpar@9 438 /* a[j] must be non-zero, x[j] must be binary, for all j in J */
alpar@9 439 for (e = ptr; e != NULL; e = e->next)
alpar@9 440 { xassert(e->aj != 0.0);
alpar@9 441 xassert(e->xj->is_int);
alpar@9 442 xassert(e->xj->lb == 0.0 && e->xj->ub == 1.0);
alpar@9 443 }
alpar@9 444 /* check if the specified inequality constraint already has the
alpar@9 445 form of packing inequality */
alpar@9 446 neg = 0; /* neg is |Jn| */
alpar@9 447 for (e = ptr; e != NULL; e = e->next)
alpar@9 448 { if (e->aj == +1.0)
alpar@9 449 ;
alpar@9 450 else if (e->aj == -1.0)
alpar@9 451 neg++;
alpar@9 452 else
alpar@9 453 break;
alpar@9 454 }
alpar@9 455 if (e == NULL)
alpar@9 456 { /* all coefficients a[j] are +1 or -1; check rhs b */
alpar@9 457 if (b == (double)(1 - neg))
alpar@9 458 { /* it is packing inequality; no processing is needed */
alpar@9 459 return 1;
alpar@9 460 }
alpar@9 461 }
alpar@9 462 /* substitute x[j] = 1 - x~[j] for all j in Jn to make all a[j]
alpar@9 463 positive; the result is a~[j] = |a[j]| and new rhs b */
alpar@9 464 for (e = ptr; e != NULL; e = e->next)
alpar@9 465 if (e->aj < 0) b -= e->aj;
alpar@9 466 /* now a[j] > 0 for all j in J (actually |a[j]| are used) */
alpar@9 467 /* if a[j] > b, skip processing--this case must not appear */
alpar@9 468 for (e = ptr; e != NULL; e = e->next)
alpar@9 469 if (fabs(e->aj) > b) return 0;
alpar@9 470 /* now 0 < a[j] <= b for all j in J */
alpar@9 471 /* find two minimal coefficients a[j] and a[k], j != k */
alpar@9 472 ej = NULL;
alpar@9 473 for (e = ptr; e != NULL; e = e->next)
alpar@9 474 if (ej == NULL || fabs(ej->aj) > fabs(e->aj)) ej = e;
alpar@9 475 xassert(ej != NULL);
alpar@9 476 ek = NULL;
alpar@9 477 for (e = ptr; e != NULL; e = e->next)
alpar@9 478 if (e != ej)
alpar@9 479 if (ek == NULL || fabs(ek->aj) > fabs(e->aj)) ek = e;
alpar@9 480 xassert(ek != NULL);
alpar@9 481 /* the specified constraint is equivalent to packing inequality
alpar@9 482 iff a[j] + a[k] > b + eps */
alpar@9 483 eps = 1e-3 + 1e-6 * fabs(b);
alpar@9 484 if (fabs(ej->aj) + fabs(ek->aj) <= b + eps) return 0;
alpar@9 485 /* perform back substitution x~[j] = 1 - x[j] and construct the
alpar@9 486 final equivalent packing inequality in generalized format */
alpar@9 487 b = 1.0;
alpar@9 488 for (e = ptr; e != NULL; e = e->next)
alpar@9 489 { if (e->aj > 0.0)
alpar@9 490 e->aj = +1.0;
alpar@9 491 else /* e->aj < 0.0 */
alpar@9 492 e->aj = -1.0, b -= 1.0;
alpar@9 493 }
alpar@9 494 *_b = b;
alpar@9 495 return 2;
alpar@9 496 }
alpar@9 497
alpar@9 498 int npp_hidden_packing(NPP *npp, NPPROW *row)
alpar@9 499 { /* identify hidden packing inequality */
alpar@9 500 NPPROW *copy;
alpar@9 501 NPPAIJ *aij;
alpar@9 502 struct elem *ptr, *e;
alpar@9 503 int kase, ret, count = 0;
alpar@9 504 double b;
alpar@9 505 /* the row must be inequality constraint */
alpar@9 506 xassert(row->lb < row->ub);
alpar@9 507 for (kase = 0; kase <= 1; kase++)
alpar@9 508 { if (kase == 0)
alpar@9 509 { /* process row upper bound */
alpar@9 510 if (row->ub == +DBL_MAX) continue;
alpar@9 511 ptr = copy_form(npp, row, +1.0);
alpar@9 512 b = + row->ub;
alpar@9 513 }
alpar@9 514 else
alpar@9 515 { /* process row lower bound */
alpar@9 516 if (row->lb == -DBL_MAX) continue;
alpar@9 517 ptr = copy_form(npp, row, -1.0);
alpar@9 518 b = - row->lb;
alpar@9 519 }
alpar@9 520 /* now the inequality has the form "sum a[j] x[j] <= b" */
alpar@9 521 ret = hidden_packing(npp, ptr, &b);
alpar@9 522 xassert(0 <= ret && ret <= 2);
alpar@9 523 if (kase == 1 && ret == 1 || ret == 2)
alpar@9 524 { /* the original inequality has been identified as hidden
alpar@9 525 packing inequality */
alpar@9 526 count++;
alpar@9 527 #ifdef GLP_DEBUG
alpar@9 528 xprintf("Original constraint:\n");
alpar@9 529 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
alpar@9 530 xprintf(" %+g x%d", aij->val, aij->col->j);
alpar@9 531 if (row->lb != -DBL_MAX) xprintf(", >= %g", row->lb);
alpar@9 532 if (row->ub != +DBL_MAX) xprintf(", <= %g", row->ub);
alpar@9 533 xprintf("\n");
alpar@9 534 xprintf("Equivalent packing inequality:\n");
alpar@9 535 for (e = ptr; e != NULL; e = e->next)
alpar@9 536 xprintf(" %sx%d", e->aj > 0.0 ? "+" : "-", e->xj->j);
alpar@9 537 xprintf(", <= %g\n", b);
alpar@9 538 #endif
alpar@9 539 if (row->lb == -DBL_MAX || row->ub == +DBL_MAX)
alpar@9 540 { /* the original row is single-sided inequality; no copy
alpar@9 541 is needed */
alpar@9 542 copy = NULL;
alpar@9 543 }
alpar@9 544 else
alpar@9 545 { /* the original row is double-sided inequality; we need
alpar@9 546 to create its copy for other bound before replacing it
alpar@9 547 with the equivalent inequality */
alpar@9 548 copy = npp_add_row(npp);
alpar@9 549 if (kase == 0)
alpar@9 550 { /* the copy is for lower bound */
alpar@9 551 copy->lb = row->lb, copy->ub = +DBL_MAX;
alpar@9 552 }
alpar@9 553 else
alpar@9 554 { /* the copy is for upper bound */
alpar@9 555 copy->lb = -DBL_MAX, copy->ub = row->ub;
alpar@9 556 }
alpar@9 557 /* copy original row coefficients */
alpar@9 558 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
alpar@9 559 npp_add_aij(npp, copy, aij->col, aij->val);
alpar@9 560 }
alpar@9 561 /* replace the original inequality by equivalent one */
alpar@9 562 npp_erase_row(npp, row);
alpar@9 563 row->lb = -DBL_MAX, row->ub = b;
alpar@9 564 for (e = ptr; e != NULL; e = e->next)
alpar@9 565 npp_add_aij(npp, row, e->xj, e->aj);
alpar@9 566 /* continue processing lower bound for the copy */
alpar@9 567 if (copy != NULL) row = copy;
alpar@9 568 }
alpar@9 569 drop_form(npp, ptr);
alpar@9 570 }
alpar@9 571 return count;
alpar@9 572 }
alpar@9 573
alpar@9 574 /***********************************************************************
alpar@9 575 * NAME
alpar@9 576 *
alpar@9 577 * npp_implied_packing - identify implied packing inequality
alpar@9 578 *
alpar@9 579 * SYNOPSIS
alpar@9 580 *
alpar@9 581 * #include "glpnpp.h"
alpar@9 582 * int npp_implied_packing(NPP *npp, NPPROW *row, int which,
alpar@9 583 * NPPCOL *var[], char set[]);
alpar@9 584 *
alpar@9 585 * DESCRIPTION
alpar@9 586 *
alpar@9 587 * The routine npp_implied_packing processes specified row (constraint)
alpar@9 588 * of general format:
alpar@9 589 *
alpar@9 590 * L <= sum a[j] x[j] <= U. (1)
alpar@9 591 * j
alpar@9 592 *
alpar@9 593 * If which = 0, only lower bound L, which must exist, is considered,
alpar@9 594 * while upper bound U is ignored. Similarly, if which = 1, only upper
alpar@9 595 * bound U, which must exist, is considered, while lower bound L is
alpar@9 596 * ignored. Thus, if the specified row is a double-sided inequality or
alpar@9 597 * equality constraint, this routine should be called twice for both
alpar@9 598 * lower and upper bounds.
alpar@9 599 *
alpar@9 600 * The routine npp_implied_packing attempts to find a non-trivial (i.e.
alpar@9 601 * having not less than two binary variables) packing inequality:
alpar@9 602 *
alpar@9 603 * sum x[j] - sum x[j] <= 1 - |Jn|, (2)
alpar@9 604 * j in Jp j in Jn
alpar@9 605 *
alpar@9 606 * which is relaxation of the constraint (1) in the sense that any
alpar@9 607 * solution satisfying to that constraint also satisfies to the packing
alpar@9 608 * inequality (2). If such relaxation exists, the routine stores
alpar@9 609 * pointers to descriptors of corresponding binary variables and their
alpar@9 610 * flags, resp., to locations var[1], var[2], ..., var[len] and set[1],
alpar@9 611 * set[2], ..., set[len], where set[j] = 0 means that j in Jp and
alpar@9 612 * set[j] = 1 means that j in Jn.
alpar@9 613 *
alpar@9 614 * RETURNS
alpar@9 615 *
alpar@9 616 * The routine npp_implied_packing returns len, which is the total
alpar@9 617 * number of binary variables in the packing inequality found, len >= 2.
alpar@9 618 * However, if the relaxation does not exist, the routine returns zero.
alpar@9 619 *
alpar@9 620 * ALGORITHM
alpar@9 621 *
alpar@9 622 * If which = 0, the constraint coefficients (1) are multiplied by -1
alpar@9 623 * and b is assigned -L; if which = 1, the constraint coefficients (1)
alpar@9 624 * are not changed and b is assigned +U. In both cases the specified
alpar@9 625 * constraint gets the following format:
alpar@9 626 *
alpar@9 627 * sum a[j] x[j] <= b. (3)
alpar@9 628 * j
alpar@9 629 *
alpar@9 630 * (Note that (3) is a relaxation of (1), because one of bounds L or U
alpar@9 631 * is ignored.)
alpar@9 632 *
alpar@9 633 * Let J be set of binary variables, Kp be set of non-binary (integer
alpar@9 634 * or continuous) variables with a[j] > 0, and Kn be set of non-binary
alpar@9 635 * variables with a[j] < 0. Then the inequality (3) can be written as
alpar@9 636 * follows:
alpar@9 637 *
alpar@9 638 * sum a[j] x[j] <= b - sum a[j] x[j] - sum a[j] x[j]. (4)
alpar@9 639 * j in J j in Kp j in Kn
alpar@9 640 *
alpar@9 641 * To get rid of non-binary variables we can replace the inequality (4)
alpar@9 642 * by the following relaxed inequality:
alpar@9 643 *
alpar@9 644 * sum a[j] x[j] <= b~, (5)
alpar@9 645 * j in J
alpar@9 646 *
alpar@9 647 * where:
alpar@9 648 *
alpar@9 649 * b~ = sup(b - sum a[j] x[j] - sum a[j] x[j]) =
alpar@9 650 * j in Kp j in Kn
alpar@9 651 *
alpar@9 652 * = b - inf sum a[j] x[j] - inf sum a[j] x[j] = (6)
alpar@9 653 * j in Kp j in Kn
alpar@9 654 *
alpar@9 655 * = b - sum a[j] l[j] - sum a[j] u[j].
alpar@9 656 * j in Kp j in Kn
alpar@9 657 *
alpar@9 658 * Note that if lower bound l[j] (if j in Kp) or upper bound u[j]
alpar@9 659 * (if j in Kn) of some non-binary variable x[j] does not exist, then
alpar@9 660 * formally b = +oo, in which case further analysis is not performed.
alpar@9 661 *
alpar@9 662 * Let Bp = {j in J: a[j] > 0}, Bn = {j in J: a[j] < 0}. To make all
alpar@9 663 * the inequality coefficients in (5) positive, we replace all x[j] in
alpar@9 664 * Bn by their complementaries, substituting x[j] = 1 - x~[j] for all
alpar@9 665 * j in Bn, that gives:
alpar@9 666 *
alpar@9 667 * sum a[j] x[j] - sum a[j] x~[j] <= b~ - sum a[j]. (7)
alpar@9 668 * j in Bp j in Bn j in Bn
alpar@9 669 *
alpar@9 670 * This inequality is a relaxation of the original constraint (1), and
alpar@9 671 * it is a binary knapsack inequality. Writing it in the standard format
alpar@9 672 * we have:
alpar@9 673 *
alpar@9 674 * sum alfa[j] z[j] <= beta, (8)
alpar@9 675 * j in J
alpar@9 676 *
alpar@9 677 * where:
alpar@9 678 * ( + a[j], if j in Bp,
alpar@9 679 * alfa[j] = < (9)
alpar@9 680 * ( - a[j], if j in Bn,
alpar@9 681 *
alpar@9 682 * ( x[j], if j in Bp,
alpar@9 683 * z[j] = < (10)
alpar@9 684 * ( 1 - x[j], if j in Bn,
alpar@9 685 *
alpar@9 686 * beta = b~ - sum a[j]. (11)
alpar@9 687 * j in Bn
alpar@9 688 *
alpar@9 689 * In the inequality (8) all coefficients are positive, therefore, the
alpar@9 690 * packing relaxation to be found for this inequality is the following:
alpar@9 691 *
alpar@9 692 * sum z[j] <= 1. (12)
alpar@9 693 * j in P
alpar@9 694 *
alpar@9 695 * It is obvious that set P within J, which we would like to find, must
alpar@9 696 * satisfy to the following condition:
alpar@9 697 *
alpar@9 698 * alfa[j] + alfa[k] > beta + eps for all j, k in P, j != k, (13)
alpar@9 699 *
alpar@9 700 * where eps is an absolute tolerance for value of the linear form.
alpar@9 701 * Thus, it is natural to take P = {j: alpha[j] > (beta + eps) / 2}.
alpar@9 702 * Moreover, if in the equality (8) there exist coefficients alfa[k],
alpar@9 703 * for which alfa[k] <= (beta + eps) / 2, but which, nevertheless,
alpar@9 704 * satisfies to the condition (13) for all j in P, *one* corresponding
alpar@9 705 * variable z[k] (having, for example, maximal coefficient alfa[k]) can
alpar@9 706 * be included in set P, that allows increasing the number of binary
alpar@9 707 * variables in (12) by one.
alpar@9 708 *
alpar@9 709 * Once the set P has been built, for the inequality (12) we need to
alpar@9 710 * perform back substitution according to (10) in order to express it
alpar@9 711 * through the original binary variables. As the result of such back
alpar@9 712 * substitution the relaxed packing inequality get its final format (2),
alpar@9 713 * where Jp = J intersect Bp, and Jn = J intersect Bn. */
alpar@9 714
alpar@9 715 int npp_implied_packing(NPP *npp, NPPROW *row, int which,
alpar@9 716 NPPCOL *var[], char set[])
alpar@9 717 { struct elem *ptr, *e, *i, *k;
alpar@9 718 int len = 0;
alpar@9 719 double b, eps;
alpar@9 720 /* build inequality (3) */
alpar@9 721 if (which == 0)
alpar@9 722 { ptr = copy_form(npp, row, -1.0);
alpar@9 723 xassert(row->lb != -DBL_MAX);
alpar@9 724 b = - row->lb;
alpar@9 725 }
alpar@9 726 else if (which == 1)
alpar@9 727 { ptr = copy_form(npp, row, +1.0);
alpar@9 728 xassert(row->ub != +DBL_MAX);
alpar@9 729 b = + row->ub;
alpar@9 730 }
alpar@9 731 /* remove non-binary variables to build relaxed inequality (5);
alpar@9 732 compute its right-hand side b~ with formula (6) */
alpar@9 733 for (e = ptr; e != NULL; e = e->next)
alpar@9 734 { if (!(e->xj->is_int && e->xj->lb == 0.0 && e->xj->ub == 1.0))
alpar@9 735 { /* x[j] is non-binary variable */
alpar@9 736 if (e->aj > 0.0)
alpar@9 737 { if (e->xj->lb == -DBL_MAX) goto done;
alpar@9 738 b -= e->aj * e->xj->lb;
alpar@9 739 }
alpar@9 740 else /* e->aj < 0.0 */
alpar@9 741 { if (e->xj->ub == +DBL_MAX) goto done;
alpar@9 742 b -= e->aj * e->xj->ub;
alpar@9 743 }
alpar@9 744 /* a[j] = 0 means that variable x[j] is removed */
alpar@9 745 e->aj = 0.0;
alpar@9 746 }
alpar@9 747 }
alpar@9 748 /* substitute x[j] = 1 - x~[j] to build knapsack inequality (8);
alpar@9 749 compute its right-hand side beta with formula (11) */
alpar@9 750 for (e = ptr; e != NULL; e = e->next)
alpar@9 751 if (e->aj < 0.0) b -= e->aj;
alpar@9 752 /* if beta is close to zero, the knapsack inequality is either
alpar@9 753 infeasible or forcing inequality; this must never happen, so
alpar@9 754 we skip further analysis */
alpar@9 755 if (b < 1e-3) goto done;
alpar@9 756 /* build set P as well as sets Jp and Jn, and determine x[k] as
alpar@9 757 explained above in comments to the routine */
alpar@9 758 eps = 1e-3 + 1e-6 * b;
alpar@9 759 i = k = NULL;
alpar@9 760 for (e = ptr; e != NULL; e = e->next)
alpar@9 761 { /* note that alfa[j] = |a[j]| */
alpar@9 762 if (fabs(e->aj) > 0.5 * (b + eps))
alpar@9 763 { /* alfa[j] > (b + eps) / 2; include x[j] in set P, i.e. in
alpar@9 764 set Jp or Jn */
alpar@9 765 var[++len] = e->xj;
alpar@9 766 set[len] = (char)(e->aj > 0.0 ? 0 : 1);
alpar@9 767 /* alfa[i] = min alfa[j] over all j included in set P */
alpar@9 768 if (i == NULL || fabs(i->aj) > fabs(e->aj)) i = e;
alpar@9 769 }
alpar@9 770 else if (fabs(e->aj) >= 1e-3)
alpar@9 771 { /* alfa[k] = max alfa[j] over all j not included in set P;
alpar@9 772 we skip coefficient a[j] if it is close to zero to avoid
alpar@9 773 numerically unreliable results */
alpar@9 774 if (k == NULL || fabs(k->aj) < fabs(e->aj)) k = e;
alpar@9 775 }
alpar@9 776 }
alpar@9 777 /* if alfa[k] satisfies to condition (13) for all j in P, include
alpar@9 778 x[k] in P */
alpar@9 779 if (i != NULL && k != NULL && fabs(i->aj) + fabs(k->aj) > b + eps)
alpar@9 780 { var[++len] = k->xj;
alpar@9 781 set[len] = (char)(k->aj > 0.0 ? 0 : 1);
alpar@9 782 }
alpar@9 783 /* trivial packing inequality being redundant must never appear,
alpar@9 784 so we just ignore it */
alpar@9 785 if (len < 2) len = 0;
alpar@9 786 done: drop_form(npp, ptr);
alpar@9 787 return len;
alpar@9 788 }
alpar@9 789
alpar@9 790 /***********************************************************************
alpar@9 791 * NAME
alpar@9 792 *
alpar@9 793 * npp_is_covering - test if constraint is covering inequality
alpar@9 794 *
alpar@9 795 * SYNOPSIS
alpar@9 796 *
alpar@9 797 * #include "glpnpp.h"
alpar@9 798 * int npp_is_covering(NPP *npp, NPPROW *row);
alpar@9 799 *
alpar@9 800 * RETURNS
alpar@9 801 *
alpar@9 802 * If the specified row (constraint) is covering inequality (see below),
alpar@9 803 * the routine npp_is_covering returns non-zero. Otherwise, it returns
alpar@9 804 * zero.
alpar@9 805 *
alpar@9 806 * COVERING INEQUALITIES
alpar@9 807 *
alpar@9 808 * In canonical format the covering inequality is the following:
alpar@9 809 *
alpar@9 810 * sum x[j] >= 1, (1)
alpar@9 811 * j in J
alpar@9 812 *
alpar@9 813 * where all variables x[j] are binary. This inequality expresses the
alpar@9 814 * condition that in any integer feasible solution variables in set J
alpar@9 815 * cannot be all equal to zero at the same time, i.e. at least one
alpar@9 816 * variable must take non-zero (unity) value. W.l.o.g. it is assumed
alpar@9 817 * that |J| >= 2, because if J is empty, the inequality (1) is
alpar@9 818 * infeasible, and if |J| = 1, the inequality (1) is a forcing row.
alpar@9 819 *
alpar@9 820 * In general case the covering inequality may include original
alpar@9 821 * variables x[j] as well as their complements x~[j]:
alpar@9 822 *
alpar@9 823 * sum x[j] + sum x~[j] >= 1, (2)
alpar@9 824 * j in Jp j in Jn
alpar@9 825 *
alpar@9 826 * where Jp and Jn are not intersected. Therefore, using substitution
alpar@9 827 * x~[j] = 1 - x[j] gives the packing inequality in generalized format:
alpar@9 828 *
alpar@9 829 * sum x[j] - sum x[j] >= 1 - |Jn|. (3)
alpar@9 830 * j in Jp j in Jn
alpar@9 831 *
alpar@9 832 * (May note that the inequality (3) cuts off infeasible solutions,
alpar@9 833 * where x[j] = 0 for all j in Jp and x[j] = 1 for all j in Jn.)
alpar@9 834 *
alpar@9 835 * NOTE: If |J| = 2, the inequality (3) is equivalent to packing
alpar@9 836 * inequality (see the routine npp_is_packing). */
alpar@9 837
alpar@9 838 int npp_is_covering(NPP *npp, NPPROW *row)
alpar@9 839 { /* test if constraint is covering inequality */
alpar@9 840 NPPCOL *col;
alpar@9 841 NPPAIJ *aij;
alpar@9 842 int b;
alpar@9 843 xassert(npp == npp);
alpar@9 844 if (!(row->lb != -DBL_MAX && row->ub == +DBL_MAX))
alpar@9 845 return 0;
alpar@9 846 b = 1;
alpar@9 847 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
alpar@9 848 { col = aij->col;
alpar@9 849 if (!(col->is_int && col->lb == 0.0 && col->ub == 1.0))
alpar@9 850 return 0;
alpar@9 851 if (aij->val == +1.0)
alpar@9 852 ;
alpar@9 853 else if (aij->val == -1.0)
alpar@9 854 b--;
alpar@9 855 else
alpar@9 856 return 0;
alpar@9 857 }
alpar@9 858 if (row->lb != (double)b) return 0;
alpar@9 859 return 1;
alpar@9 860 }
alpar@9 861
alpar@9 862 /***********************************************************************
alpar@9 863 * NAME
alpar@9 864 *
alpar@9 865 * npp_hidden_covering - identify hidden covering inequality
alpar@9 866 *
alpar@9 867 * SYNOPSIS
alpar@9 868 *
alpar@9 869 * #include "glpnpp.h"
alpar@9 870 * int npp_hidden_covering(NPP *npp, NPPROW *row);
alpar@9 871 *
alpar@9 872 * DESCRIPTION
alpar@9 873 *
alpar@9 874 * The routine npp_hidden_covering processes specified inequality
alpar@9 875 * constraint, which includes only binary variables, and the number of
alpar@9 876 * the variables is not less than three. If the original inequality is
alpar@9 877 * equivalent to a covering inequality (see below), the routine
alpar@9 878 * replaces it by the equivalent inequality. If the original constraint
alpar@9 879 * is double-sided inequality, it is replaced by a pair of single-sided
alpar@9 880 * inequalities, if necessary.
alpar@9 881 *
alpar@9 882 * RETURNS
alpar@9 883 *
alpar@9 884 * If the original inequality constraint was replaced by equivalent
alpar@9 885 * covering inequality, the routine npp_hidden_covering returns
alpar@9 886 * non-zero. Otherwise, it returns zero.
alpar@9 887 *
alpar@9 888 * PROBLEM TRANSFORMATION
alpar@9 889 *
alpar@9 890 * Consider an inequality constraint:
alpar@9 891 *
alpar@9 892 * sum a[j] x[j] >= b, (1)
alpar@9 893 * j in J
alpar@9 894 *
alpar@9 895 * where all variables x[j] are binary, and |J| >= 3. (In case of '<='
alpar@9 896 * inequality it can be transformed to '>=' format by multiplying both
alpar@9 897 * its sides by -1.)
alpar@9 898 *
alpar@9 899 * Let Jp = {j: a[j] > 0}, Jn = {j: a[j] < 0}. Performing substitution
alpar@9 900 * x[j] = 1 - x~[j] for all j in Jn, we have:
alpar@9 901 *
alpar@9 902 * sum a[j] x[j] >= b ==>
alpar@9 903 * j in J
alpar@9 904 *
alpar@9 905 * sum a[j] x[j] + sum a[j] x[j] >= b ==>
alpar@9 906 * j in Jp j in Jn
alpar@9 907 *
alpar@9 908 * sum a[j] x[j] + sum a[j] (1 - x~[j]) >= b ==>
alpar@9 909 * j in Jp j in Jn
alpar@9 910 *
alpar@9 911 * sum m a[j] x[j] - sum a[j] x~[j] >= b - sum a[j].
alpar@9 912 * j in Jp j in Jn j in Jn
alpar@9 913 *
alpar@9 914 * Thus, meaning the transformation above, we can assume that in
alpar@9 915 * inequality (1) all coefficients a[j] are positive. Moreover, we can
alpar@9 916 * assume that b > 0, because otherwise the inequality (1) would be
alpar@9 917 * redundant (see the routine npp_analyze_row). It is then obvious that
alpar@9 918 * constraint (1) is equivalent to covering inequality only if:
alpar@9 919 *
alpar@9 920 * a[j] >= b, (2)
alpar@9 921 *
alpar@9 922 * for all j in J.
alpar@9 923 *
alpar@9 924 * Once the original inequality (1) is replaced by equivalent covering
alpar@9 925 * inequality, we need to perform back substitution x~[j] = 1 - x[j] for
alpar@9 926 * all j in Jn (see above).
alpar@9 927 *
alpar@9 928 * RECOVERING SOLUTION
alpar@9 929 *
alpar@9 930 * None needed. */
alpar@9 931
alpar@9 932 static int hidden_covering(NPP *npp, struct elem *ptr, double *_b)
alpar@9 933 { /* process inequality constraint: sum a[j] x[j] >= b;
alpar@9 934 0 - specified row is NOT hidden covering inequality;
alpar@9 935 1 - specified row is covering inequality;
alpar@9 936 2 - specified row is hidden covering inequality. */
alpar@9 937 struct elem *e;
alpar@9 938 int neg;
alpar@9 939 double b = *_b, eps;
alpar@9 940 xassert(npp == npp);
alpar@9 941 /* a[j] must be non-zero, x[j] must be binary, for all j in J */
alpar@9 942 for (e = ptr; e != NULL; e = e->next)
alpar@9 943 { xassert(e->aj != 0.0);
alpar@9 944 xassert(e->xj->is_int);
alpar@9 945 xassert(e->xj->lb == 0.0 && e->xj->ub == 1.0);
alpar@9 946 }
alpar@9 947 /* check if the specified inequality constraint already has the
alpar@9 948 form of covering inequality */
alpar@9 949 neg = 0; /* neg is |Jn| */
alpar@9 950 for (e = ptr; e != NULL; e = e->next)
alpar@9 951 { if (e->aj == +1.0)
alpar@9 952 ;
alpar@9 953 else if (e->aj == -1.0)
alpar@9 954 neg++;
alpar@9 955 else
alpar@9 956 break;
alpar@9 957 }
alpar@9 958 if (e == NULL)
alpar@9 959 { /* all coefficients a[j] are +1 or -1; check rhs b */
alpar@9 960 if (b == (double)(1 - neg))
alpar@9 961 { /* it is covering inequality; no processing is needed */
alpar@9 962 return 1;
alpar@9 963 }
alpar@9 964 }
alpar@9 965 /* substitute x[j] = 1 - x~[j] for all j in Jn to make all a[j]
alpar@9 966 positive; the result is a~[j] = |a[j]| and new rhs b */
alpar@9 967 for (e = ptr; e != NULL; e = e->next)
alpar@9 968 if (e->aj < 0) b -= e->aj;
alpar@9 969 /* now a[j] > 0 for all j in J (actually |a[j]| are used) */
alpar@9 970 /* if b <= 0, skip processing--this case must not appear */
alpar@9 971 if (b < 1e-3) return 0;
alpar@9 972 /* now a[j] > 0 for all j in J, and b > 0 */
alpar@9 973 /* the specified constraint is equivalent to covering inequality
alpar@9 974 iff a[j] >= b for all j in J */
alpar@9 975 eps = 1e-9 + 1e-12 * fabs(b);
alpar@9 976 for (e = ptr; e != NULL; e = e->next)
alpar@9 977 if (fabs(e->aj) < b - eps) return 0;
alpar@9 978 /* perform back substitution x~[j] = 1 - x[j] and construct the
alpar@9 979 final equivalent covering inequality in generalized format */
alpar@9 980 b = 1.0;
alpar@9 981 for (e = ptr; e != NULL; e = e->next)
alpar@9 982 { if (e->aj > 0.0)
alpar@9 983 e->aj = +1.0;
alpar@9 984 else /* e->aj < 0.0 */
alpar@9 985 e->aj = -1.0, b -= 1.0;
alpar@9 986 }
alpar@9 987 *_b = b;
alpar@9 988 return 2;
alpar@9 989 }
alpar@9 990
alpar@9 991 int npp_hidden_covering(NPP *npp, NPPROW *row)
alpar@9 992 { /* identify hidden covering inequality */
alpar@9 993 NPPROW *copy;
alpar@9 994 NPPAIJ *aij;
alpar@9 995 struct elem *ptr, *e;
alpar@9 996 int kase, ret, count = 0;
alpar@9 997 double b;
alpar@9 998 /* the row must be inequality constraint */
alpar@9 999 xassert(row->lb < row->ub);
alpar@9 1000 for (kase = 0; kase <= 1; kase++)
alpar@9 1001 { if (kase == 0)
alpar@9 1002 { /* process row lower bound */
alpar@9 1003 if (row->lb == -DBL_MAX) continue;
alpar@9 1004 ptr = copy_form(npp, row, +1.0);
alpar@9 1005 b = + row->lb;
alpar@9 1006 }
alpar@9 1007 else
alpar@9 1008 { /* process row upper bound */
alpar@9 1009 if (row->ub == +DBL_MAX) continue;
alpar@9 1010 ptr = copy_form(npp, row, -1.0);
alpar@9 1011 b = - row->ub;
alpar@9 1012 }
alpar@9 1013 /* now the inequality has the form "sum a[j] x[j] >= b" */
alpar@9 1014 ret = hidden_covering(npp, ptr, &b);
alpar@9 1015 xassert(0 <= ret && ret <= 2);
alpar@9 1016 if (kase == 1 && ret == 1 || ret == 2)
alpar@9 1017 { /* the original inequality has been identified as hidden
alpar@9 1018 covering inequality */
alpar@9 1019 count++;
alpar@9 1020 #ifdef GLP_DEBUG
alpar@9 1021 xprintf("Original constraint:\n");
alpar@9 1022 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
alpar@9 1023 xprintf(" %+g x%d", aij->val, aij->col->j);
alpar@9 1024 if (row->lb != -DBL_MAX) xprintf(", >= %g", row->lb);
alpar@9 1025 if (row->ub != +DBL_MAX) xprintf(", <= %g", row->ub);
alpar@9 1026 xprintf("\n");
alpar@9 1027 xprintf("Equivalent covering inequality:\n");
alpar@9 1028 for (e = ptr; e != NULL; e = e->next)
alpar@9 1029 xprintf(" %sx%d", e->aj > 0.0 ? "+" : "-", e->xj->j);
alpar@9 1030 xprintf(", >= %g\n", b);
alpar@9 1031 #endif
alpar@9 1032 if (row->lb == -DBL_MAX || row->ub == +DBL_MAX)
alpar@9 1033 { /* the original row is single-sided inequality; no copy
alpar@9 1034 is needed */
alpar@9 1035 copy = NULL;
alpar@9 1036 }
alpar@9 1037 else
alpar@9 1038 { /* the original row is double-sided inequality; we need
alpar@9 1039 to create its copy for other bound before replacing it
alpar@9 1040 with the equivalent inequality */
alpar@9 1041 copy = npp_add_row(npp);
alpar@9 1042 if (kase == 0)
alpar@9 1043 { /* the copy is for upper bound */
alpar@9 1044 copy->lb = -DBL_MAX, copy->ub = row->ub;
alpar@9 1045 }
alpar@9 1046 else
alpar@9 1047 { /* the copy is for lower bound */
alpar@9 1048 copy->lb = row->lb, copy->ub = +DBL_MAX;
alpar@9 1049 }
alpar@9 1050 /* copy original row coefficients */
alpar@9 1051 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
alpar@9 1052 npp_add_aij(npp, copy, aij->col, aij->val);
alpar@9 1053 }
alpar@9 1054 /* replace the original inequality by equivalent one */
alpar@9 1055 npp_erase_row(npp, row);
alpar@9 1056 row->lb = b, row->ub = +DBL_MAX;
alpar@9 1057 for (e = ptr; e != NULL; e = e->next)
alpar@9 1058 npp_add_aij(npp, row, e->xj, e->aj);
alpar@9 1059 /* continue processing upper bound for the copy */
alpar@9 1060 if (copy != NULL) row = copy;
alpar@9 1061 }
alpar@9 1062 drop_form(npp, ptr);
alpar@9 1063 }
alpar@9 1064 return count;
alpar@9 1065 }
alpar@9 1066
alpar@9 1067 /***********************************************************************
alpar@9 1068 * NAME
alpar@9 1069 *
alpar@9 1070 * npp_is_partitioning - test if constraint is partitioning equality
alpar@9 1071 *
alpar@9 1072 * SYNOPSIS
alpar@9 1073 *
alpar@9 1074 * #include "glpnpp.h"
alpar@9 1075 * int npp_is_partitioning(NPP *npp, NPPROW *row);
alpar@9 1076 *
alpar@9 1077 * RETURNS
alpar@9 1078 *
alpar@9 1079 * If the specified row (constraint) is partitioning equality (see
alpar@9 1080 * below), the routine npp_is_partitioning returns non-zero. Otherwise,
alpar@9 1081 * it returns zero.
alpar@9 1082 *
alpar@9 1083 * PARTITIONING EQUALITIES
alpar@9 1084 *
alpar@9 1085 * In canonical format the partitioning equality is the following:
alpar@9 1086 *
alpar@9 1087 * sum x[j] = 1, (1)
alpar@9 1088 * j in J
alpar@9 1089 *
alpar@9 1090 * where all variables x[j] are binary. This equality expresses the
alpar@9 1091 * condition that in any integer feasible solution exactly one variable
alpar@9 1092 * in set J must take non-zero (unity) value while other variables must
alpar@9 1093 * be equal to zero. W.l.o.g. it is assumed that |J| >= 2, because if
alpar@9 1094 * J is empty, the inequality (1) is infeasible, and if |J| = 1, the
alpar@9 1095 * inequality (1) is a fixing row.
alpar@9 1096 *
alpar@9 1097 * In general case the partitioning equality may include original
alpar@9 1098 * variables x[j] as well as their complements x~[j]:
alpar@9 1099 *
alpar@9 1100 * sum x[j] + sum x~[j] = 1, (2)
alpar@9 1101 * j in Jp j in Jn
alpar@9 1102 *
alpar@9 1103 * where Jp and Jn are not intersected. Therefore, using substitution
alpar@9 1104 * x~[j] = 1 - x[j] leads to the partitioning equality in generalized
alpar@9 1105 * format:
alpar@9 1106 *
alpar@9 1107 * sum x[j] - sum x[j] = 1 - |Jn|. (3)
alpar@9 1108 * j in Jp j in Jn */
alpar@9 1109
alpar@9 1110 int npp_is_partitioning(NPP *npp, NPPROW *row)
alpar@9 1111 { /* test if constraint is partitioning equality */
alpar@9 1112 NPPCOL *col;
alpar@9 1113 NPPAIJ *aij;
alpar@9 1114 int b;
alpar@9 1115 xassert(npp == npp);
alpar@9 1116 if (row->lb != row->ub) return 0;
alpar@9 1117 b = 1;
alpar@9 1118 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
alpar@9 1119 { col = aij->col;
alpar@9 1120 if (!(col->is_int && col->lb == 0.0 && col->ub == 1.0))
alpar@9 1121 return 0;
alpar@9 1122 if (aij->val == +1.0)
alpar@9 1123 ;
alpar@9 1124 else if (aij->val == -1.0)
alpar@9 1125 b--;
alpar@9 1126 else
alpar@9 1127 return 0;
alpar@9 1128 }
alpar@9 1129 if (row->lb != (double)b) return 0;
alpar@9 1130 return 1;
alpar@9 1131 }
alpar@9 1132
alpar@9 1133 /***********************************************************************
alpar@9 1134 * NAME
alpar@9 1135 *
alpar@9 1136 * npp_reduce_ineq_coef - reduce inequality constraint coefficients
alpar@9 1137 *
alpar@9 1138 * SYNOPSIS
alpar@9 1139 *
alpar@9 1140 * #include "glpnpp.h"
alpar@9 1141 * int npp_reduce_ineq_coef(NPP *npp, NPPROW *row);
alpar@9 1142 *
alpar@9 1143 * DESCRIPTION
alpar@9 1144 *
alpar@9 1145 * The routine npp_reduce_ineq_coef processes specified inequality
alpar@9 1146 * constraint attempting to replace it by an equivalent constraint,
alpar@9 1147 * where magnitude of coefficients at binary variables is smaller than
alpar@9 1148 * in the original constraint. If the inequality is double-sided, it is
alpar@9 1149 * replaced by a pair of single-sided inequalities, if necessary.
alpar@9 1150 *
alpar@9 1151 * RETURNS
alpar@9 1152 *
alpar@9 1153 * The routine npp_reduce_ineq_coef returns the number of coefficients
alpar@9 1154 * reduced.
alpar@9 1155 *
alpar@9 1156 * BACKGROUND
alpar@9 1157 *
alpar@9 1158 * Consider an inequality constraint:
alpar@9 1159 *
alpar@9 1160 * sum a[j] x[j] >= b. (1)
alpar@9 1161 * j in J
alpar@9 1162 *
alpar@9 1163 * (In case of '<=' inequality it can be transformed to '>=' format by
alpar@9 1164 * multiplying both its sides by -1.) Let x[k] be a binary variable;
alpar@9 1165 * other variables can be integer as well as continuous. We can write
alpar@9 1166 * constraint (1) as follows:
alpar@9 1167 *
alpar@9 1168 * a[k] x[k] + t[k] >= b, (2)
alpar@9 1169 *
alpar@9 1170 * where:
alpar@9 1171 *
alpar@9 1172 * t[k] = sum a[j] x[j]. (3)
alpar@9 1173 * j in J\{k}
alpar@9 1174 *
alpar@9 1175 * Since x[k] is binary, constraint (2) is equivalent to disjunction of
alpar@9 1176 * the following two constraints:
alpar@9 1177 *
alpar@9 1178 * x[k] = 0, t[k] >= b (4)
alpar@9 1179 *
alpar@9 1180 * OR
alpar@9 1181 *
alpar@9 1182 * x[k] = 1, t[k] >= b - a[k]. (5)
alpar@9 1183 *
alpar@9 1184 * Let also that for the partial sum t[k] be known some its implied
alpar@9 1185 * lower bound inf t[k].
alpar@9 1186 *
alpar@9 1187 * Case a[k] > 0. Let inf t[k] < b, since otherwise both constraints
alpar@9 1188 * (4) and (5) and therefore constraint (2) are redundant.
alpar@9 1189 * If inf t[k] > b - a[k], only constraint (5) is redundant, in which
alpar@9 1190 * case it can be replaced with the following redundant and therefore
alpar@9 1191 * equivalent constraint:
alpar@9 1192 *
alpar@9 1193 * t[k] >= b - a'[k] = inf t[k], (6)
alpar@9 1194 *
alpar@9 1195 * where:
alpar@9 1196 *
alpar@9 1197 * a'[k] = b - inf t[k]. (7)
alpar@9 1198 *
alpar@9 1199 * Thus, the original constraint (2) is equivalent to the following
alpar@9 1200 * constraint with coefficient at variable x[k] changed:
alpar@9 1201 *
alpar@9 1202 * a'[k] x[k] + t[k] >= b. (8)
alpar@9 1203 *
alpar@9 1204 * From inf t[k] < b it follows that a'[k] > 0, i.e. the coefficient
alpar@9 1205 * at x[k] keeps its sign. And from inf t[k] > b - a[k] it follows that
alpar@9 1206 * a'[k] < a[k], i.e. the coefficient reduces in magnitude.
alpar@9 1207 *
alpar@9 1208 * Case a[k] < 0. Let inf t[k] < b - a[k], since otherwise both
alpar@9 1209 * constraints (4) and (5) and therefore constraint (2) are redundant.
alpar@9 1210 * If inf t[k] > b, only constraint (4) is redundant, in which case it
alpar@9 1211 * can be replaced with the following redundant and therefore equivalent
alpar@9 1212 * constraint:
alpar@9 1213 *
alpar@9 1214 * t[k] >= b' = inf t[k]. (9)
alpar@9 1215 *
alpar@9 1216 * Rewriting constraint (5) as follows:
alpar@9 1217 *
alpar@9 1218 * t[k] >= b - a[k] = b' - a'[k], (10)
alpar@9 1219 *
alpar@9 1220 * where:
alpar@9 1221 *
alpar@9 1222 * a'[k] = a[k] + b' - b = a[k] + inf t[k] - b, (11)
alpar@9 1223 *
alpar@9 1224 * we can see that disjunction of constraint (9) and (10) is equivalent
alpar@9 1225 * to disjunction of constraint (4) and (5), from which it follows that
alpar@9 1226 * the original constraint (2) is equivalent to the following constraint
alpar@9 1227 * with both coefficient at variable x[k] and right-hand side changed:
alpar@9 1228 *
alpar@9 1229 * a'[k] x[k] + t[k] >= b'. (12)
alpar@9 1230 *
alpar@9 1231 * From inf t[k] < b - a[k] it follows that a'[k] < 0, i.e. the
alpar@9 1232 * coefficient at x[k] keeps its sign. And from inf t[k] > b it follows
alpar@9 1233 * that a'[k] > a[k], i.e. the coefficient reduces in magnitude.
alpar@9 1234 *
alpar@9 1235 * PROBLEM TRANSFORMATION
alpar@9 1236 *
alpar@9 1237 * In the routine npp_reduce_ineq_coef the following implied lower
alpar@9 1238 * bound of the partial sum (3) is used:
alpar@9 1239 *
alpar@9 1240 * inf t[k] = sum a[j] l[j] + sum a[j] u[j], (13)
alpar@9 1241 * j in Jp\{k} k in Jn\{k}
alpar@9 1242 *
alpar@9 1243 * where Jp = {j : a[j] > 0}, Jn = {j : a[j] < 0}, l[j] and u[j] are
alpar@9 1244 * lower and upper bounds, resp., of variable x[j].
alpar@9 1245 *
alpar@9 1246 * In order to compute inf t[k] more efficiently, the following formula,
alpar@9 1247 * which is equivalent to (13), is actually used:
alpar@9 1248 *
alpar@9 1249 * ( h - a[k] l[k] = h, if a[k] > 0,
alpar@9 1250 * inf t[k] = < (14)
alpar@9 1251 * ( h - a[k] u[k] = h - a[k], if a[k] < 0,
alpar@9 1252 *
alpar@9 1253 * where:
alpar@9 1254 *
alpar@9 1255 * h = sum a[j] l[j] + sum a[j] u[j] (15)
alpar@9 1256 * j in Jp j in Jn
alpar@9 1257 *
alpar@9 1258 * is the implied lower bound of row (1).
alpar@9 1259 *
alpar@9 1260 * Reduction of positive coefficient (a[k] > 0) does not change value
alpar@9 1261 * of h, since l[k] = 0. In case of reduction of negative coefficient
alpar@9 1262 * (a[k] < 0) from (11) it follows that:
alpar@9 1263 *
alpar@9 1264 * delta a[k] = a'[k] - a[k] = inf t[k] - b (> 0), (16)
alpar@9 1265 *
alpar@9 1266 * so new value of h (accounting that u[k] = 1) can be computed as
alpar@9 1267 * follows:
alpar@9 1268 *
alpar@9 1269 * h := h + delta a[k] = h + (inf t[k] - b). (17)
alpar@9 1270 *
alpar@9 1271 * RECOVERING SOLUTION
alpar@9 1272 *
alpar@9 1273 * None needed. */
alpar@9 1274
alpar@9 1275 static int reduce_ineq_coef(NPP *npp, struct elem *ptr, double *_b)
alpar@9 1276 { /* process inequality constraint: sum a[j] x[j] >= b */
alpar@9 1277 /* returns: the number of coefficients reduced */
alpar@9 1278 struct elem *e;
alpar@9 1279 int count = 0;
alpar@9 1280 double h, inf_t, new_a, b = *_b;
alpar@9 1281 xassert(npp == npp);
alpar@9 1282 /* compute h; see (15) */
alpar@9 1283 h = 0.0;
alpar@9 1284 for (e = ptr; e != NULL; e = e->next)
alpar@9 1285 { if (e->aj > 0.0)
alpar@9 1286 { if (e->xj->lb == -DBL_MAX) goto done;
alpar@9 1287 h += e->aj * e->xj->lb;
alpar@9 1288 }
alpar@9 1289 else /* e->aj < 0.0 */
alpar@9 1290 { if (e->xj->ub == +DBL_MAX) goto done;
alpar@9 1291 h += e->aj * e->xj->ub;
alpar@9 1292 }
alpar@9 1293 }
alpar@9 1294 /* perform reduction of coefficients at binary variables */
alpar@9 1295 for (e = ptr; e != NULL; e = e->next)
alpar@9 1296 { /* skip non-binary variable */
alpar@9 1297 if (!(e->xj->is_int && e->xj->lb == 0.0 && e->xj->ub == 1.0))
alpar@9 1298 continue;
alpar@9 1299 if (e->aj > 0.0)
alpar@9 1300 { /* compute inf t[k]; see (14) */
alpar@9 1301 inf_t = h;
alpar@9 1302 if (b - e->aj < inf_t && inf_t < b)
alpar@9 1303 { /* compute reduced coefficient a'[k]; see (7) */
alpar@9 1304 new_a = b - inf_t;
alpar@9 1305 if (new_a >= +1e-3 &&
alpar@9 1306 e->aj - new_a >= 0.01 * (1.0 + e->aj))
alpar@9 1307 { /* accept a'[k] */
alpar@9 1308 #ifdef GLP_DEBUG
alpar@9 1309 xprintf("+");
alpar@9 1310 #endif
alpar@9 1311 e->aj = new_a;
alpar@9 1312 count++;
alpar@9 1313 }
alpar@9 1314 }
alpar@9 1315 }
alpar@9 1316 else /* e->aj < 0.0 */
alpar@9 1317 { /* compute inf t[k]; see (14) */
alpar@9 1318 inf_t = h - e->aj;
alpar@9 1319 if (b < inf_t && inf_t < b - e->aj)
alpar@9 1320 { /* compute reduced coefficient a'[k]; see (11) */
alpar@9 1321 new_a = e->aj + (inf_t - b);
alpar@9 1322 if (new_a <= -1e-3 &&
alpar@9 1323 new_a - e->aj >= 0.01 * (1.0 - e->aj))
alpar@9 1324 { /* accept a'[k] */
alpar@9 1325 #ifdef GLP_DEBUG
alpar@9 1326 xprintf("-");
alpar@9 1327 #endif
alpar@9 1328 e->aj = new_a;
alpar@9 1329 /* update h; see (17) */
alpar@9 1330 h += (inf_t - b);
alpar@9 1331 /* compute b'; see (9) */
alpar@9 1332 b = inf_t;
alpar@9 1333 count++;
alpar@9 1334 }
alpar@9 1335 }
alpar@9 1336 }
alpar@9 1337 }
alpar@9 1338 *_b = b;
alpar@9 1339 done: return count;
alpar@9 1340 }
alpar@9 1341
alpar@9 1342 int npp_reduce_ineq_coef(NPP *npp, NPPROW *row)
alpar@9 1343 { /* reduce inequality constraint coefficients */
alpar@9 1344 NPPROW *copy;
alpar@9 1345 NPPAIJ *aij;
alpar@9 1346 struct elem *ptr, *e;
alpar@9 1347 int kase, count[2];
alpar@9 1348 double b;
alpar@9 1349 /* the row must be inequality constraint */
alpar@9 1350 xassert(row->lb < row->ub);
alpar@9 1351 count[0] = count[1] = 0;
alpar@9 1352 for (kase = 0; kase <= 1; kase++)
alpar@9 1353 { if (kase == 0)
alpar@9 1354 { /* process row lower bound */
alpar@9 1355 if (row->lb == -DBL_MAX) continue;
alpar@9 1356 #ifdef GLP_DEBUG
alpar@9 1357 xprintf("L");
alpar@9 1358 #endif
alpar@9 1359 ptr = copy_form(npp, row, +1.0);
alpar@9 1360 b = + row->lb;
alpar@9 1361 }
alpar@9 1362 else
alpar@9 1363 { /* process row upper bound */
alpar@9 1364 if (row->ub == +DBL_MAX) continue;
alpar@9 1365 #ifdef GLP_DEBUG
alpar@9 1366 xprintf("U");
alpar@9 1367 #endif
alpar@9 1368 ptr = copy_form(npp, row, -1.0);
alpar@9 1369 b = - row->ub;
alpar@9 1370 }
alpar@9 1371 /* now the inequality has the form "sum a[j] x[j] >= b" */
alpar@9 1372 count[kase] = reduce_ineq_coef(npp, ptr, &b);
alpar@9 1373 if (count[kase] > 0)
alpar@9 1374 { /* the original inequality has been replaced by equivalent
alpar@9 1375 one with coefficients reduced */
alpar@9 1376 if (row->lb == -DBL_MAX || row->ub == +DBL_MAX)
alpar@9 1377 { /* the original row is single-sided inequality; no copy
alpar@9 1378 is needed */
alpar@9 1379 copy = NULL;
alpar@9 1380 }
alpar@9 1381 else
alpar@9 1382 { /* the original row is double-sided inequality; we need
alpar@9 1383 to create its copy for other bound before replacing it
alpar@9 1384 with the equivalent inequality */
alpar@9 1385 #ifdef GLP_DEBUG
alpar@9 1386 xprintf("*");
alpar@9 1387 #endif
alpar@9 1388 copy = npp_add_row(npp);
alpar@9 1389 if (kase == 0)
alpar@9 1390 { /* the copy is for upper bound */
alpar@9 1391 copy->lb = -DBL_MAX, copy->ub = row->ub;
alpar@9 1392 }
alpar@9 1393 else
alpar@9 1394 { /* the copy is for lower bound */
alpar@9 1395 copy->lb = row->lb, copy->ub = +DBL_MAX;
alpar@9 1396 }
alpar@9 1397 /* copy original row coefficients */
alpar@9 1398 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
alpar@9 1399 npp_add_aij(npp, copy, aij->col, aij->val);
alpar@9 1400 }
alpar@9 1401 /* replace the original inequality by equivalent one */
alpar@9 1402 npp_erase_row(npp, row);
alpar@9 1403 row->lb = b, row->ub = +DBL_MAX;
alpar@9 1404 for (e = ptr; e != NULL; e = e->next)
alpar@9 1405 npp_add_aij(npp, row, e->xj, e->aj);
alpar@9 1406 /* continue processing upper bound for the copy */
alpar@9 1407 if (copy != NULL) row = copy;
alpar@9 1408 }
alpar@9 1409 drop_form(npp, ptr);
alpar@9 1410 }
alpar@9 1411 return count[0] + count[1];
alpar@9 1412 }
alpar@9 1413
alpar@9 1414 /* eof */