lemon-project-template-glpk

annotate deps/glpk/src/glpios07.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 /* glpios07.c (mixed cover cut generator) */
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 "glpios.h"
alpar@9 26
alpar@9 27 /*----------------------------------------------------------------------
alpar@9 28 -- COVER INEQUALITIES
alpar@9 29 --
alpar@9 30 -- Consider the set of feasible solutions to 0-1 knapsack problem:
alpar@9 31 --
alpar@9 32 -- sum a[j]*x[j] <= b, (1)
alpar@9 33 -- j in J
alpar@9 34 --
alpar@9 35 -- x[j] is binary, (2)
alpar@9 36 --
alpar@9 37 -- where, wlog, we assume that a[j] > 0 (since 0-1 variables can be
alpar@9 38 -- complemented) and a[j] <= b (since a[j] > b implies x[j] = 0).
alpar@9 39 --
alpar@9 40 -- A set C within J is called a cover if
alpar@9 41 --
alpar@9 42 -- sum a[j] > b. (3)
alpar@9 43 -- j in C
alpar@9 44 --
alpar@9 45 -- For any cover C the inequality
alpar@9 46 --
alpar@9 47 -- sum x[j] <= |C| - 1 (4)
alpar@9 48 -- j in C
alpar@9 49 --
alpar@9 50 -- is called a cover inequality and is valid for (1)-(2).
alpar@9 51 --
alpar@9 52 -- MIXED COVER INEQUALITIES
alpar@9 53 --
alpar@9 54 -- Consider the set of feasible solutions to mixed knapsack problem:
alpar@9 55 --
alpar@9 56 -- sum a[j]*x[j] + y <= b, (5)
alpar@9 57 -- j in J
alpar@9 58 --
alpar@9 59 -- x[j] is binary, (6)
alpar@9 60 --
alpar@9 61 -- 0 <= y <= u is continuous, (7)
alpar@9 62 --
alpar@9 63 -- where again we assume that a[j] > 0.
alpar@9 64 --
alpar@9 65 -- Let C within J be some set. From (1)-(4) it follows that
alpar@9 66 --
alpar@9 67 -- sum a[j] > b - y (8)
alpar@9 68 -- j in C
alpar@9 69 --
alpar@9 70 -- implies
alpar@9 71 --
alpar@9 72 -- sum x[j] <= |C| - 1. (9)
alpar@9 73 -- j in C
alpar@9 74 --
alpar@9 75 -- Thus, we need to modify the inequality (9) in such a way that it be
alpar@9 76 -- a constraint only if the condition (8) is satisfied.
alpar@9 77 --
alpar@9 78 -- Consider the following inequality:
alpar@9 79 --
alpar@9 80 -- sum x[j] <= |C| - t. (10)
alpar@9 81 -- j in C
alpar@9 82 --
alpar@9 83 -- If 0 < t <= 1, then (10) is equivalent to (9), because all x[j] are
alpar@9 84 -- binary variables. On the other hand, if t <= 0, (10) being satisfied
alpar@9 85 -- for any values of x[j] is not a constraint.
alpar@9 86 --
alpar@9 87 -- Let
alpar@9 88 --
alpar@9 89 -- t' = sum a[j] + y - b. (11)
alpar@9 90 -- j in C
alpar@9 91 --
alpar@9 92 -- It is understood that the condition t' > 0 is equivalent to (8).
alpar@9 93 -- Besides, from (6)-(7) it follows that t' has an implied upper bound:
alpar@9 94 --
alpar@9 95 -- t'max = sum a[j] + u - b. (12)
alpar@9 96 -- j in C
alpar@9 97 --
alpar@9 98 -- This allows to express the parameter t having desired properties:
alpar@9 99 --
alpar@9 100 -- t = t' / t'max. (13)
alpar@9 101 --
alpar@9 102 -- In fact, t <= 1 by definition, and t > 0 being equivalent to t' > 0
alpar@9 103 -- is equivalent to (8).
alpar@9 104 --
alpar@9 105 -- Thus, the inequality (10), where t is given by formula (13) is valid
alpar@9 106 -- for (5)-(7).
alpar@9 107 --
alpar@9 108 -- Note that if u = 0, then y = 0, so t = 1, and the conditions (8) and
alpar@9 109 -- (10) is transformed to the conditions (3) and (4).
alpar@9 110 --
alpar@9 111 -- GENERATING MIXED COVER CUTS
alpar@9 112 --
alpar@9 113 -- To generate a mixed cover cut in the form (10) we need to find such
alpar@9 114 -- set C which satisfies to the inequality (8) and for which, in turn,
alpar@9 115 -- the inequality (10) is violated in the current point.
alpar@9 116 --
alpar@9 117 -- Substituting t from (13) to (10) gives:
alpar@9 118 --
alpar@9 119 -- 1
alpar@9 120 -- sum x[j] <= |C| - ----- (sum a[j] + y - b), (14)
alpar@9 121 -- j in C t'max j in C
alpar@9 122 --
alpar@9 123 -- and finally we have the cut inequality in the standard form:
alpar@9 124 --
alpar@9 125 -- sum x[j] + alfa * y <= beta, (15)
alpar@9 126 -- j in C
alpar@9 127 --
alpar@9 128 -- where:
alpar@9 129 --
alpar@9 130 -- alfa = 1 / t'max, (16)
alpar@9 131 --
alpar@9 132 -- beta = |C| - alfa * (sum a[j] - b). (17)
alpar@9 133 -- j in C */
alpar@9 134
alpar@9 135 #if 1
alpar@9 136 #define MAXTRY 1000
alpar@9 137 #else
alpar@9 138 #define MAXTRY 10000
alpar@9 139 #endif
alpar@9 140
alpar@9 141 static int cover2(int n, double a[], double b, double u, double x[],
alpar@9 142 double y, int cov[], double *_alfa, double *_beta)
alpar@9 143 { /* try to generate mixed cover cut using two-element cover */
alpar@9 144 int i, j, try = 0, ret = 0;
alpar@9 145 double eps, alfa, beta, temp, rmax = 0.001;
alpar@9 146 eps = 0.001 * (1.0 + fabs(b));
alpar@9 147 for (i = 0+1; i <= n; i++)
alpar@9 148 for (j = i+1; j <= n; j++)
alpar@9 149 { /* C = {i, j} */
alpar@9 150 try++;
alpar@9 151 if (try > MAXTRY) goto done;
alpar@9 152 /* check if condition (8) is satisfied */
alpar@9 153 if (a[i] + a[j] + y > b + eps)
alpar@9 154 { /* compute parameters for inequality (15) */
alpar@9 155 temp = a[i] + a[j] - b;
alpar@9 156 alfa = 1.0 / (temp + u);
alpar@9 157 beta = 2.0 - alfa * temp;
alpar@9 158 /* compute violation of inequality (15) */
alpar@9 159 temp = x[i] + x[j] + alfa * y - beta;
alpar@9 160 /* choose C providing maximum violation */
alpar@9 161 if (rmax < temp)
alpar@9 162 { rmax = temp;
alpar@9 163 cov[1] = i;
alpar@9 164 cov[2] = j;
alpar@9 165 *_alfa = alfa;
alpar@9 166 *_beta = beta;
alpar@9 167 ret = 1;
alpar@9 168 }
alpar@9 169 }
alpar@9 170 }
alpar@9 171 done: return ret;
alpar@9 172 }
alpar@9 173
alpar@9 174 static int cover3(int n, double a[], double b, double u, double x[],
alpar@9 175 double y, int cov[], double *_alfa, double *_beta)
alpar@9 176 { /* try to generate mixed cover cut using three-element cover */
alpar@9 177 int i, j, k, try = 0, ret = 0;
alpar@9 178 double eps, alfa, beta, temp, rmax = 0.001;
alpar@9 179 eps = 0.001 * (1.0 + fabs(b));
alpar@9 180 for (i = 0+1; i <= n; i++)
alpar@9 181 for (j = i+1; j <= n; j++)
alpar@9 182 for (k = j+1; k <= n; k++)
alpar@9 183 { /* C = {i, j, k} */
alpar@9 184 try++;
alpar@9 185 if (try > MAXTRY) goto done;
alpar@9 186 /* check if condition (8) is satisfied */
alpar@9 187 if (a[i] + a[j] + a[k] + y > b + eps)
alpar@9 188 { /* compute parameters for inequality (15) */
alpar@9 189 temp = a[i] + a[j] + a[k] - b;
alpar@9 190 alfa = 1.0 / (temp + u);
alpar@9 191 beta = 3.0 - alfa * temp;
alpar@9 192 /* compute violation of inequality (15) */
alpar@9 193 temp = x[i] + x[j] + x[k] + alfa * y - beta;
alpar@9 194 /* choose C providing maximum violation */
alpar@9 195 if (rmax < temp)
alpar@9 196 { rmax = temp;
alpar@9 197 cov[1] = i;
alpar@9 198 cov[2] = j;
alpar@9 199 cov[3] = k;
alpar@9 200 *_alfa = alfa;
alpar@9 201 *_beta = beta;
alpar@9 202 ret = 1;
alpar@9 203 }
alpar@9 204 }
alpar@9 205 }
alpar@9 206 done: return ret;
alpar@9 207 }
alpar@9 208
alpar@9 209 static int cover4(int n, double a[], double b, double u, double x[],
alpar@9 210 double y, int cov[], double *_alfa, double *_beta)
alpar@9 211 { /* try to generate mixed cover cut using four-element cover */
alpar@9 212 int i, j, k, l, try = 0, ret = 0;
alpar@9 213 double eps, alfa, beta, temp, rmax = 0.001;
alpar@9 214 eps = 0.001 * (1.0 + fabs(b));
alpar@9 215 for (i = 0+1; i <= n; i++)
alpar@9 216 for (j = i+1; j <= n; j++)
alpar@9 217 for (k = j+1; k <= n; k++)
alpar@9 218 for (l = k+1; l <= n; l++)
alpar@9 219 { /* C = {i, j, k, l} */
alpar@9 220 try++;
alpar@9 221 if (try > MAXTRY) goto done;
alpar@9 222 /* check if condition (8) is satisfied */
alpar@9 223 if (a[i] + a[j] + a[k] + a[l] + y > b + eps)
alpar@9 224 { /* compute parameters for inequality (15) */
alpar@9 225 temp = a[i] + a[j] + a[k] + a[l] - b;
alpar@9 226 alfa = 1.0 / (temp + u);
alpar@9 227 beta = 4.0 - alfa * temp;
alpar@9 228 /* compute violation of inequality (15) */
alpar@9 229 temp = x[i] + x[j] + x[k] + x[l] + alfa * y - beta;
alpar@9 230 /* choose C providing maximum violation */
alpar@9 231 if (rmax < temp)
alpar@9 232 { rmax = temp;
alpar@9 233 cov[1] = i;
alpar@9 234 cov[2] = j;
alpar@9 235 cov[3] = k;
alpar@9 236 cov[4] = l;
alpar@9 237 *_alfa = alfa;
alpar@9 238 *_beta = beta;
alpar@9 239 ret = 1;
alpar@9 240 }
alpar@9 241 }
alpar@9 242 }
alpar@9 243 done: return ret;
alpar@9 244 }
alpar@9 245
alpar@9 246 static int cover(int n, double a[], double b, double u, double x[],
alpar@9 247 double y, int cov[], double *alfa, double *beta)
alpar@9 248 { /* try to generate mixed cover cut;
alpar@9 249 input (see (5)):
alpar@9 250 n is the number of binary variables;
alpar@9 251 a[1:n] are coefficients at binary variables;
alpar@9 252 b is the right-hand side;
alpar@9 253 u is upper bound of continuous variable;
alpar@9 254 x[1:n] are values of binary variables at current point;
alpar@9 255 y is value of continuous variable at current point;
alpar@9 256 output (see (15), (16), (17)):
alpar@9 257 cov[1:r] are indices of binary variables included in cover C,
alpar@9 258 where r is the set cardinality returned on exit;
alpar@9 259 alfa coefficient at continuous variable;
alpar@9 260 beta is the right-hand side; */
alpar@9 261 int j;
alpar@9 262 /* perform some sanity checks */
alpar@9 263 xassert(n >= 2);
alpar@9 264 for (j = 1; j <= n; j++) xassert(a[j] > 0.0);
alpar@9 265 #if 1 /* ??? */
alpar@9 266 xassert(b > -1e-5);
alpar@9 267 #else
alpar@9 268 xassert(b > 0.0);
alpar@9 269 #endif
alpar@9 270 xassert(u >= 0.0);
alpar@9 271 for (j = 1; j <= n; j++) xassert(0.0 <= x[j] && x[j] <= 1.0);
alpar@9 272 xassert(0.0 <= y && y <= u);
alpar@9 273 /* try to generate mixed cover cut */
alpar@9 274 if (cover2(n, a, b, u, x, y, cov, alfa, beta)) return 2;
alpar@9 275 if (cover3(n, a, b, u, x, y, cov, alfa, beta)) return 3;
alpar@9 276 if (cover4(n, a, b, u, x, y, cov, alfa, beta)) return 4;
alpar@9 277 return 0;
alpar@9 278 }
alpar@9 279
alpar@9 280 /*----------------------------------------------------------------------
alpar@9 281 -- lpx_cover_cut - generate mixed cover cut.
alpar@9 282 --
alpar@9 283 -- SYNOPSIS
alpar@9 284 --
alpar@9 285 -- #include "glplpx.h"
alpar@9 286 -- int lpx_cover_cut(LPX *lp, int len, int ind[], double val[],
alpar@9 287 -- double work[]);
alpar@9 288 --
alpar@9 289 -- DESCRIPTION
alpar@9 290 --
alpar@9 291 -- The routine lpx_cover_cut generates a mixed cover cut for a given
alpar@9 292 -- row of the MIP problem.
alpar@9 293 --
alpar@9 294 -- The given row of the MIP problem should be explicitly specified in
alpar@9 295 -- the form:
alpar@9 296 --
alpar@9 297 -- sum{j in J} a[j]*x[j] <= b. (1)
alpar@9 298 --
alpar@9 299 -- On entry indices (ordinal numbers) of structural variables, which
alpar@9 300 -- have non-zero constraint coefficients, should be placed in locations
alpar@9 301 -- ind[1], ..., ind[len], and corresponding constraint coefficients
alpar@9 302 -- should be placed in locations val[1], ..., val[len]. The right-hand
alpar@9 303 -- side b should be stored in location val[0].
alpar@9 304 --
alpar@9 305 -- The working array work should have at least nb locations, where nb
alpar@9 306 -- is the number of binary variables in (1).
alpar@9 307 --
alpar@9 308 -- The routine generates a mixed cover cut in the same form as (1) and
alpar@9 309 -- stores the cut coefficients and right-hand side in the same way as
alpar@9 310 -- just described above.
alpar@9 311 --
alpar@9 312 -- RETURNS
alpar@9 313 --
alpar@9 314 -- If the cutting plane has been successfully generated, the routine
alpar@9 315 -- returns 1 <= len' <= n, which is the number of non-zero coefficients
alpar@9 316 -- in the inequality constraint. Otherwise, the routine returns zero. */
alpar@9 317
alpar@9 318 static int lpx_cover_cut(LPX *lp, int len, int ind[], double val[],
alpar@9 319 double work[])
alpar@9 320 { int cov[1+4], j, k, nb, newlen, r;
alpar@9 321 double f_min, f_max, alfa, beta, u, *x = work, y;
alpar@9 322 /* substitute and remove fixed variables */
alpar@9 323 newlen = 0;
alpar@9 324 for (k = 1; k <= len; k++)
alpar@9 325 { j = ind[k];
alpar@9 326 if (lpx_get_col_type(lp, j) == LPX_FX)
alpar@9 327 val[0] -= val[k] * lpx_get_col_lb(lp, j);
alpar@9 328 else
alpar@9 329 { newlen++;
alpar@9 330 ind[newlen] = ind[k];
alpar@9 331 val[newlen] = val[k];
alpar@9 332 }
alpar@9 333 }
alpar@9 334 len = newlen;
alpar@9 335 /* move binary variables to the beginning of the list so that
alpar@9 336 elements 1, 2, ..., nb correspond to binary variables, and
alpar@9 337 elements nb+1, nb+2, ..., len correspond to rest variables */
alpar@9 338 nb = 0;
alpar@9 339 for (k = 1; k <= len; k++)
alpar@9 340 { j = ind[k];
alpar@9 341 if (lpx_get_col_kind(lp, j) == LPX_IV &&
alpar@9 342 lpx_get_col_type(lp, j) == LPX_DB &&
alpar@9 343 lpx_get_col_lb(lp, j) == 0.0 &&
alpar@9 344 lpx_get_col_ub(lp, j) == 1.0)
alpar@9 345 { /* binary variable */
alpar@9 346 int ind_k;
alpar@9 347 double val_k;
alpar@9 348 nb++;
alpar@9 349 ind_k = ind[nb], val_k = val[nb];
alpar@9 350 ind[nb] = ind[k], val[nb] = val[k];
alpar@9 351 ind[k] = ind_k, val[k] = val_k;
alpar@9 352 }
alpar@9 353 }
alpar@9 354 /* now the specified row has the form:
alpar@9 355 sum a[j]*x[j] + sum a[j]*y[j] <= b,
alpar@9 356 where x[j] are binary variables, y[j] are rest variables */
alpar@9 357 /* at least two binary variables are needed */
alpar@9 358 if (nb < 2) return 0;
alpar@9 359 /* compute implied lower and upper bounds for sum a[j]*y[j] */
alpar@9 360 f_min = f_max = 0.0;
alpar@9 361 for (k = nb+1; k <= len; k++)
alpar@9 362 { j = ind[k];
alpar@9 363 /* both bounds must be finite */
alpar@9 364 if (lpx_get_col_type(lp, j) != LPX_DB) return 0;
alpar@9 365 if (val[k] > 0.0)
alpar@9 366 { f_min += val[k] * lpx_get_col_lb(lp, j);
alpar@9 367 f_max += val[k] * lpx_get_col_ub(lp, j);
alpar@9 368 }
alpar@9 369 else
alpar@9 370 { f_min += val[k] * lpx_get_col_ub(lp, j);
alpar@9 371 f_max += val[k] * lpx_get_col_lb(lp, j);
alpar@9 372 }
alpar@9 373 }
alpar@9 374 /* sum a[j]*x[j] + sum a[j]*y[j] <= b ===>
alpar@9 375 sum a[j]*x[j] + (sum a[j]*y[j] - f_min) <= b - f_min ===>
alpar@9 376 sum a[j]*x[j] + y <= b - f_min,
alpar@9 377 where y = sum a[j]*y[j] - f_min;
alpar@9 378 note that 0 <= y <= u, u = f_max - f_min */
alpar@9 379 /* determine upper bound of y */
alpar@9 380 u = f_max - f_min;
alpar@9 381 /* determine value of y at the current point */
alpar@9 382 y = 0.0;
alpar@9 383 for (k = nb+1; k <= len; k++)
alpar@9 384 { j = ind[k];
alpar@9 385 y += val[k] * lpx_get_col_prim(lp, j);
alpar@9 386 }
alpar@9 387 y -= f_min;
alpar@9 388 if (y < 0.0) y = 0.0;
alpar@9 389 if (y > u) y = u;
alpar@9 390 /* modify the right-hand side b */
alpar@9 391 val[0] -= f_min;
alpar@9 392 /* now the transformed row has the form:
alpar@9 393 sum a[j]*x[j] + y <= b, where 0 <= y <= u */
alpar@9 394 /* determine values of x[j] at the current point */
alpar@9 395 for (k = 1; k <= nb; k++)
alpar@9 396 { j = ind[k];
alpar@9 397 x[k] = lpx_get_col_prim(lp, j);
alpar@9 398 if (x[k] < 0.0) x[k] = 0.0;
alpar@9 399 if (x[k] > 1.0) x[k] = 1.0;
alpar@9 400 }
alpar@9 401 /* if a[j] < 0, replace x[j] by its complement 1 - x'[j] */
alpar@9 402 for (k = 1; k <= nb; k++)
alpar@9 403 { if (val[k] < 0.0)
alpar@9 404 { ind[k] = - ind[k];
alpar@9 405 val[k] = - val[k];
alpar@9 406 val[0] += val[k];
alpar@9 407 x[k] = 1.0 - x[k];
alpar@9 408 }
alpar@9 409 }
alpar@9 410 /* try to generate a mixed cover cut for the transformed row */
alpar@9 411 r = cover(nb, val, val[0], u, x, y, cov, &alfa, &beta);
alpar@9 412 if (r == 0) return 0;
alpar@9 413 xassert(2 <= r && r <= 4);
alpar@9 414 /* now the cut is in the form:
alpar@9 415 sum{j in C} x[j] + alfa * y <= beta */
alpar@9 416 /* store the right-hand side beta */
alpar@9 417 ind[0] = 0, val[0] = beta;
alpar@9 418 /* restore the original ordinal numbers of x[j] */
alpar@9 419 for (j = 1; j <= r; j++) cov[j] = ind[cov[j]];
alpar@9 420 /* store cut coefficients at binary variables complementing back
alpar@9 421 the variables having negative row coefficients */
alpar@9 422 xassert(r <= nb);
alpar@9 423 for (k = 1; k <= r; k++)
alpar@9 424 { if (cov[k] > 0)
alpar@9 425 { ind[k] = +cov[k];
alpar@9 426 val[k] = +1.0;
alpar@9 427 }
alpar@9 428 else
alpar@9 429 { ind[k] = -cov[k];
alpar@9 430 val[k] = -1.0;
alpar@9 431 val[0] -= 1.0;
alpar@9 432 }
alpar@9 433 }
alpar@9 434 /* substitute y = sum a[j]*y[j] - f_min */
alpar@9 435 for (k = nb+1; k <= len; k++)
alpar@9 436 { r++;
alpar@9 437 ind[r] = ind[k];
alpar@9 438 val[r] = alfa * val[k];
alpar@9 439 }
alpar@9 440 val[0] += alfa * f_min;
alpar@9 441 xassert(r <= len);
alpar@9 442 len = r;
alpar@9 443 return len;
alpar@9 444 }
alpar@9 445
alpar@9 446 /*----------------------------------------------------------------------
alpar@9 447 -- lpx_eval_row - compute explictily specified row.
alpar@9 448 --
alpar@9 449 -- SYNOPSIS
alpar@9 450 --
alpar@9 451 -- #include "glplpx.h"
alpar@9 452 -- double lpx_eval_row(LPX *lp, int len, int ind[], double val[]);
alpar@9 453 --
alpar@9 454 -- DESCRIPTION
alpar@9 455 --
alpar@9 456 -- The routine lpx_eval_row computes the primal value of an explicitly
alpar@9 457 -- specified row using current values of structural variables.
alpar@9 458 --
alpar@9 459 -- The explicitly specified row may be thought as a linear form:
alpar@9 460 --
alpar@9 461 -- y = a[1]*x[m+1] + a[2]*x[m+2] + ... + a[n]*x[m+n],
alpar@9 462 --
alpar@9 463 -- where y is an auxiliary variable for this row, a[j] are coefficients
alpar@9 464 -- of the linear form, x[m+j] are structural variables.
alpar@9 465 --
alpar@9 466 -- On entry column indices and numerical values of non-zero elements of
alpar@9 467 -- the row should be stored in locations ind[1], ..., ind[len] and
alpar@9 468 -- val[1], ..., val[len], where len is the number of non-zero elements.
alpar@9 469 -- The array ind and val are not changed on exit.
alpar@9 470 --
alpar@9 471 -- RETURNS
alpar@9 472 --
alpar@9 473 -- The routine returns a computed value of y, the auxiliary variable of
alpar@9 474 -- the specified row. */
alpar@9 475
alpar@9 476 static double lpx_eval_row(LPX *lp, int len, int ind[], double val[])
alpar@9 477 { int n = lpx_get_num_cols(lp);
alpar@9 478 int j, k;
alpar@9 479 double sum = 0.0;
alpar@9 480 if (len < 0)
alpar@9 481 xerror("lpx_eval_row: len = %d; invalid row length\n", len);
alpar@9 482 for (k = 1; k <= len; k++)
alpar@9 483 { j = ind[k];
alpar@9 484 if (!(1 <= j && j <= n))
alpar@9 485 xerror("lpx_eval_row: j = %d; column number out of range\n",
alpar@9 486 j);
alpar@9 487 sum += val[k] * lpx_get_col_prim(lp, j);
alpar@9 488 }
alpar@9 489 return sum;
alpar@9 490 }
alpar@9 491
alpar@9 492 /***********************************************************************
alpar@9 493 * NAME
alpar@9 494 *
alpar@9 495 * ios_cov_gen - generate mixed cover cuts
alpar@9 496 *
alpar@9 497 * SYNOPSIS
alpar@9 498 *
alpar@9 499 * #include "glpios.h"
alpar@9 500 * void ios_cov_gen(glp_tree *tree);
alpar@9 501 *
alpar@9 502 * DESCRIPTION
alpar@9 503 *
alpar@9 504 * The routine ios_cov_gen generates mixed cover cuts for the current
alpar@9 505 * point and adds them to the cut pool. */
alpar@9 506
alpar@9 507 void ios_cov_gen(glp_tree *tree)
alpar@9 508 { glp_prob *prob = tree->mip;
alpar@9 509 int m = lpx_get_num_rows(prob);
alpar@9 510 int n = lpx_get_num_cols(prob);
alpar@9 511 int i, k, type, kase, len, *ind;
alpar@9 512 double r, *val, *work;
alpar@9 513 xassert(lpx_get_status(prob) == LPX_OPT);
alpar@9 514 /* allocate working arrays */
alpar@9 515 ind = xcalloc(1+n, sizeof(int));
alpar@9 516 val = xcalloc(1+n, sizeof(double));
alpar@9 517 work = xcalloc(1+n, sizeof(double));
alpar@9 518 /* look through all rows */
alpar@9 519 for (i = 1; i <= m; i++)
alpar@9 520 for (kase = 1; kase <= 2; kase++)
alpar@9 521 { type = lpx_get_row_type(prob, i);
alpar@9 522 if (kase == 1)
alpar@9 523 { /* consider rows of '<=' type */
alpar@9 524 if (!(type == LPX_UP || type == LPX_DB)) continue;
alpar@9 525 len = lpx_get_mat_row(prob, i, ind, val);
alpar@9 526 val[0] = lpx_get_row_ub(prob, i);
alpar@9 527 }
alpar@9 528 else
alpar@9 529 { /* consider rows of '>=' type */
alpar@9 530 if (!(type == LPX_LO || type == LPX_DB)) continue;
alpar@9 531 len = lpx_get_mat_row(prob, i, ind, val);
alpar@9 532 for (k = 1; k <= len; k++) val[k] = - val[k];
alpar@9 533 val[0] = - lpx_get_row_lb(prob, i);
alpar@9 534 }
alpar@9 535 /* generate mixed cover cut:
alpar@9 536 sum{j in J} a[j] * x[j] <= b */
alpar@9 537 len = lpx_cover_cut(prob, len, ind, val, work);
alpar@9 538 if (len == 0) continue;
alpar@9 539 /* at the current point the cut inequality is violated, i.e.
alpar@9 540 sum{j in J} a[j] * x[j] - b > 0 */
alpar@9 541 r = lpx_eval_row(prob, len, ind, val) - val[0];
alpar@9 542 if (r < 1e-3) continue;
alpar@9 543 /* add the cut to the cut pool */
alpar@9 544 glp_ios_add_row(tree, NULL, GLP_RF_COV, 0, len, ind, val,
alpar@9 545 GLP_UP, val[0]);
alpar@9 546 }
alpar@9 547 /* free working arrays */
alpar@9 548 xfree(ind);
alpar@9 549 xfree(val);
alpar@9 550 xfree(work);
alpar@9 551 return;
alpar@9 552 }
alpar@9 553
alpar@9 554 /* eof */