1 | /* glpios02.c (preprocess current subproblem) */ |
---|
2 | |
---|
3 | /*********************************************************************** |
---|
4 | * This code is part of GLPK (GNU Linear Programming Kit). |
---|
5 | * |
---|
6 | * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, |
---|
7 | * 2009, 2010 Andrew Makhorin, Department for Applied Informatics, |
---|
8 | * Moscow Aviation Institute, Moscow, Russia. All rights reserved. |
---|
9 | * E-mail: <mao@gnu.org>. |
---|
10 | * |
---|
11 | * GLPK is free software: you can redistribute it and/or modify it |
---|
12 | * under the terms of the GNU General Public License as published by |
---|
13 | * the Free Software Foundation, either version 3 of the License, or |
---|
14 | * (at your option) any later version. |
---|
15 | * |
---|
16 | * GLPK is distributed in the hope that it will be useful, but WITHOUT |
---|
17 | * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
---|
18 | * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public |
---|
19 | * License for more details. |
---|
20 | * |
---|
21 | * You should have received a copy of the GNU General Public License |
---|
22 | * along with GLPK. If not, see <http://www.gnu.org/licenses/>. |
---|
23 | ***********************************************************************/ |
---|
24 | |
---|
25 | #include "glpios.h" |
---|
26 | |
---|
27 | /*********************************************************************** |
---|
28 | * prepare_row_info - prepare row info to determine implied bounds |
---|
29 | * |
---|
30 | * Given a row (linear form) |
---|
31 | * |
---|
32 | * n |
---|
33 | * sum a[j] * x[j] (1) |
---|
34 | * j=1 |
---|
35 | * |
---|
36 | * and bounds of columns (variables) |
---|
37 | * |
---|
38 | * l[j] <= x[j] <= u[j] (2) |
---|
39 | * |
---|
40 | * this routine computes f_min, j_min, f_max, j_max needed to determine |
---|
41 | * implied bounds. |
---|
42 | * |
---|
43 | * ALGORITHM |
---|
44 | * |
---|
45 | * Let J+ = {j : a[j] > 0} and J- = {j : a[j] < 0}. |
---|
46 | * |
---|
47 | * Parameters f_min and j_min are computed as follows: |
---|
48 | * |
---|
49 | * 1) if there is no x[k] such that k in J+ and l[k] = -inf or k in J- |
---|
50 | * and u[k] = +inf, then |
---|
51 | * |
---|
52 | * f_min := sum a[j] * l[j] + sum a[j] * u[j] |
---|
53 | * j in J+ j in J- |
---|
54 | * (3) |
---|
55 | * j_min := 0 |
---|
56 | * |
---|
57 | * 2) if there is exactly one x[k] such that k in J+ and l[k] = -inf |
---|
58 | * or k in J- and u[k] = +inf, then |
---|
59 | * |
---|
60 | * f_min := sum a[j] * l[j] + sum a[j] * u[j] |
---|
61 | * j in J+\{k} j in J-\{k} |
---|
62 | * (4) |
---|
63 | * j_min := k |
---|
64 | * |
---|
65 | * 3) if there are two or more x[k] such that k in J+ and l[k] = -inf |
---|
66 | * or k in J- and u[k] = +inf, then |
---|
67 | * |
---|
68 | * f_min := -inf |
---|
69 | * (5) |
---|
70 | * j_min := 0 |
---|
71 | * |
---|
72 | * Parameters f_max and j_max are computed in a similar way as follows: |
---|
73 | * |
---|
74 | * 1) if there is no x[k] such that k in J+ and u[k] = +inf or k in J- |
---|
75 | * and l[k] = -inf, then |
---|
76 | * |
---|
77 | * f_max := sum a[j] * u[j] + sum a[j] * l[j] |
---|
78 | * j in J+ j in J- |
---|
79 | * (6) |
---|
80 | * j_max := 0 |
---|
81 | * |
---|
82 | * 2) if there is exactly one x[k] such that k in J+ and u[k] = +inf |
---|
83 | * or k in J- and l[k] = -inf, then |
---|
84 | * |
---|
85 | * f_max := sum a[j] * u[j] + sum a[j] * l[j] |
---|
86 | * j in J+\{k} j in J-\{k} |
---|
87 | * (7) |
---|
88 | * j_max := k |
---|
89 | * |
---|
90 | * 3) if there are two or more x[k] such that k in J+ and u[k] = +inf |
---|
91 | * or k in J- and l[k] = -inf, then |
---|
92 | * |
---|
93 | * f_max := +inf |
---|
94 | * (8) |
---|
95 | * j_max := 0 */ |
---|
96 | |
---|
97 | struct f_info |
---|
98 | { int j_min, j_max; |
---|
99 | double f_min, f_max; |
---|
100 | }; |
---|
101 | |
---|
102 | static void prepare_row_info(int n, const double a[], const double l[], |
---|
103 | const double u[], struct f_info *f) |
---|
104 | { int j, j_min, j_max; |
---|
105 | double f_min, f_max; |
---|
106 | xassert(n >= 0); |
---|
107 | /* determine f_min and j_min */ |
---|
108 | f_min = 0.0, j_min = 0; |
---|
109 | for (j = 1; j <= n; j++) |
---|
110 | { if (a[j] > 0.0) |
---|
111 | { if (l[j] == -DBL_MAX) |
---|
112 | { if (j_min == 0) |
---|
113 | j_min = j; |
---|
114 | else |
---|
115 | { f_min = -DBL_MAX, j_min = 0; |
---|
116 | break; |
---|
117 | } |
---|
118 | } |
---|
119 | else |
---|
120 | f_min += a[j] * l[j]; |
---|
121 | } |
---|
122 | else if (a[j] < 0.0) |
---|
123 | { if (u[j] == +DBL_MAX) |
---|
124 | { if (j_min == 0) |
---|
125 | j_min = j; |
---|
126 | else |
---|
127 | { f_min = -DBL_MAX, j_min = 0; |
---|
128 | break; |
---|
129 | } |
---|
130 | } |
---|
131 | else |
---|
132 | f_min += a[j] * u[j]; |
---|
133 | } |
---|
134 | else |
---|
135 | xassert(a != a); |
---|
136 | } |
---|
137 | f->f_min = f_min, f->j_min = j_min; |
---|
138 | /* determine f_max and j_max */ |
---|
139 | f_max = 0.0, j_max = 0; |
---|
140 | for (j = 1; j <= n; j++) |
---|
141 | { if (a[j] > 0.0) |
---|
142 | { if (u[j] == +DBL_MAX) |
---|
143 | { if (j_max == 0) |
---|
144 | j_max = j; |
---|
145 | else |
---|
146 | { f_max = +DBL_MAX, j_max = 0; |
---|
147 | break; |
---|
148 | } |
---|
149 | } |
---|
150 | else |
---|
151 | f_max += a[j] * u[j]; |
---|
152 | } |
---|
153 | else if (a[j] < 0.0) |
---|
154 | { if (l[j] == -DBL_MAX) |
---|
155 | { if (j_max == 0) |
---|
156 | j_max = j; |
---|
157 | else |
---|
158 | { f_max = +DBL_MAX, j_max = 0; |
---|
159 | break; |
---|
160 | } |
---|
161 | } |
---|
162 | else |
---|
163 | f_max += a[j] * l[j]; |
---|
164 | } |
---|
165 | else |
---|
166 | xassert(a != a); |
---|
167 | } |
---|
168 | f->f_max = f_max, f->j_max = j_max; |
---|
169 | return; |
---|
170 | } |
---|
171 | |
---|
172 | /*********************************************************************** |
---|
173 | * row_implied_bounds - determine row implied bounds |
---|
174 | * |
---|
175 | * Given a row (linear form) |
---|
176 | * |
---|
177 | * n |
---|
178 | * sum a[j] * x[j] |
---|
179 | * j=1 |
---|
180 | * |
---|
181 | * and bounds of columns (variables) |
---|
182 | * |
---|
183 | * l[j] <= x[j] <= u[j] |
---|
184 | * |
---|
185 | * this routine determines implied bounds of the row. |
---|
186 | * |
---|
187 | * ALGORITHM |
---|
188 | * |
---|
189 | * Let J+ = {j : a[j] > 0} and J- = {j : a[j] < 0}. |
---|
190 | * |
---|
191 | * The implied lower bound of the row is computed as follows: |
---|
192 | * |
---|
193 | * L' := sum a[j] * l[j] + sum a[j] * u[j] (9) |
---|
194 | * j in J+ j in J- |
---|
195 | * |
---|
196 | * and as it follows from (3), (4), and (5): |
---|
197 | * |
---|
198 | * L' := if j_min = 0 then f_min else -inf (10) |
---|
199 | * |
---|
200 | * The implied upper bound of the row is computed as follows: |
---|
201 | * |
---|
202 | * U' := sum a[j] * u[j] + sum a[j] * l[j] (11) |
---|
203 | * j in J+ j in J- |
---|
204 | * |
---|
205 | * and as it follows from (6), (7), and (8): |
---|
206 | * |
---|
207 | * U' := if j_max = 0 then f_max else +inf (12) |
---|
208 | * |
---|
209 | * The implied bounds are stored in locations LL and UU. */ |
---|
210 | |
---|
211 | static void row_implied_bounds(const struct f_info *f, double *LL, |
---|
212 | double *UU) |
---|
213 | { *LL = (f->j_min == 0 ? f->f_min : -DBL_MAX); |
---|
214 | *UU = (f->j_max == 0 ? f->f_max : +DBL_MAX); |
---|
215 | return; |
---|
216 | } |
---|
217 | |
---|
218 | /*********************************************************************** |
---|
219 | * col_implied_bounds - determine column implied bounds |
---|
220 | * |
---|
221 | * Given a row (constraint) |
---|
222 | * |
---|
223 | * n |
---|
224 | * L <= sum a[j] * x[j] <= U (13) |
---|
225 | * j=1 |
---|
226 | * |
---|
227 | * and bounds of columns (variables) |
---|
228 | * |
---|
229 | * l[j] <= x[j] <= u[j] |
---|
230 | * |
---|
231 | * this routine determines implied bounds of variable x[k]. |
---|
232 | * |
---|
233 | * It is assumed that if L != -inf, the lower bound of the row can be |
---|
234 | * active, and if U != +inf, the upper bound of the row can be active. |
---|
235 | * |
---|
236 | * ALGORITHM |
---|
237 | * |
---|
238 | * From (13) it follows that |
---|
239 | * |
---|
240 | * L <= sum a[j] * x[j] + a[k] * x[k] <= U |
---|
241 | * j!=k |
---|
242 | * or |
---|
243 | * |
---|
244 | * L - sum a[j] * x[j] <= a[k] * x[k] <= U - sum a[j] * x[j] |
---|
245 | * j!=k j!=k |
---|
246 | * |
---|
247 | * Thus, if the row lower bound L can be active, implied lower bound of |
---|
248 | * term a[k] * x[k] can be determined as follows: |
---|
249 | * |
---|
250 | * ilb(a[k] * x[k]) = min(L - sum a[j] * x[j]) = |
---|
251 | * j!=k |
---|
252 | * (14) |
---|
253 | * = L - max sum a[j] * x[j] |
---|
254 | * j!=k |
---|
255 | * |
---|
256 | * where, as it follows from (6), (7), and (8) |
---|
257 | * |
---|
258 | * / f_max - a[k] * u[k], j_max = 0, a[k] > 0 |
---|
259 | * | |
---|
260 | * | f_max - a[k] * l[k], j_max = 0, a[k] < 0 |
---|
261 | * max sum a[j] * x[j] = { |
---|
262 | * j!=k | f_max, j_max = k |
---|
263 | * | |
---|
264 | * \ +inf, j_max != 0 |
---|
265 | * |
---|
266 | * and if the upper bound U can be active, implied upper bound of term |
---|
267 | * a[k] * x[k] can be determined as follows: |
---|
268 | * |
---|
269 | * iub(a[k] * x[k]) = max(U - sum a[j] * x[j]) = |
---|
270 | * j!=k |
---|
271 | * (15) |
---|
272 | * = U - min sum a[j] * x[j] |
---|
273 | * j!=k |
---|
274 | * |
---|
275 | * where, as it follows from (3), (4), and (5) |
---|
276 | * |
---|
277 | * / f_min - a[k] * l[k], j_min = 0, a[k] > 0 |
---|
278 | * | |
---|
279 | * | f_min - a[k] * u[k], j_min = 0, a[k] < 0 |
---|
280 | * min sum a[j] * x[j] = { |
---|
281 | * j!=k | f_min, j_min = k |
---|
282 | * | |
---|
283 | * \ -inf, j_min != 0 |
---|
284 | * |
---|
285 | * Since |
---|
286 | * |
---|
287 | * ilb(a[k] * x[k]) <= a[k] * x[k] <= iub(a[k] * x[k]) |
---|
288 | * |
---|
289 | * implied lower and upper bounds of x[k] are determined as follows: |
---|
290 | * |
---|
291 | * l'[k] := if a[k] > 0 then ilb / a[k] else ulb / a[k] (16) |
---|
292 | * |
---|
293 | * u'[k] := if a[k] > 0 then ulb / a[k] else ilb / a[k] (17) |
---|
294 | * |
---|
295 | * The implied bounds are stored in locations ll and uu. */ |
---|
296 | |
---|
297 | static void col_implied_bounds(const struct f_info *f, int n, |
---|
298 | const double a[], double L, double U, const double l[], |
---|
299 | const double u[], int k, double *ll, double *uu) |
---|
300 | { double ilb, iub; |
---|
301 | xassert(n >= 0); |
---|
302 | xassert(1 <= k && k <= n); |
---|
303 | /* determine implied lower bound of term a[k] * x[k] (14) */ |
---|
304 | if (L == -DBL_MAX || f->f_max == +DBL_MAX) |
---|
305 | ilb = -DBL_MAX; |
---|
306 | else if (f->j_max == 0) |
---|
307 | { if (a[k] > 0.0) |
---|
308 | { xassert(u[k] != +DBL_MAX); |
---|
309 | ilb = L - (f->f_max - a[k] * u[k]); |
---|
310 | } |
---|
311 | else if (a[k] < 0.0) |
---|
312 | { xassert(l[k] != -DBL_MAX); |
---|
313 | ilb = L - (f->f_max - a[k] * l[k]); |
---|
314 | } |
---|
315 | else |
---|
316 | xassert(a != a); |
---|
317 | } |
---|
318 | else if (f->j_max == k) |
---|
319 | ilb = L - f->f_max; |
---|
320 | else |
---|
321 | ilb = -DBL_MAX; |
---|
322 | /* determine implied upper bound of term a[k] * x[k] (15) */ |
---|
323 | if (U == +DBL_MAX || f->f_min == -DBL_MAX) |
---|
324 | iub = +DBL_MAX; |
---|
325 | else if (f->j_min == 0) |
---|
326 | { if (a[k] > 0.0) |
---|
327 | { xassert(l[k] != -DBL_MAX); |
---|
328 | iub = U - (f->f_min - a[k] * l[k]); |
---|
329 | } |
---|
330 | else if (a[k] < 0.0) |
---|
331 | { xassert(u[k] != +DBL_MAX); |
---|
332 | iub = U - (f->f_min - a[k] * u[k]); |
---|
333 | } |
---|
334 | else |
---|
335 | xassert(a != a); |
---|
336 | } |
---|
337 | else if (f->j_min == k) |
---|
338 | iub = U - f->f_min; |
---|
339 | else |
---|
340 | iub = +DBL_MAX; |
---|
341 | /* determine implied bounds of x[k] (16) and (17) */ |
---|
342 | #if 1 |
---|
343 | /* do not use a[k] if it has small magnitude to prevent wrong |
---|
344 | implied bounds; for example, 1e-15 * x1 >= x2 + x3, where |
---|
345 | x1 >= -10, x2, x3 >= 0, would lead to wrong conclusion that |
---|
346 | x1 >= 0 */ |
---|
347 | if (fabs(a[k]) < 1e-6) |
---|
348 | *ll = -DBL_MAX, *uu = +DBL_MAX; else |
---|
349 | #endif |
---|
350 | if (a[k] > 0.0) |
---|
351 | { *ll = (ilb == -DBL_MAX ? -DBL_MAX : ilb / a[k]); |
---|
352 | *uu = (iub == +DBL_MAX ? +DBL_MAX : iub / a[k]); |
---|
353 | } |
---|
354 | else if (a[k] < 0.0) |
---|
355 | { *ll = (iub == +DBL_MAX ? -DBL_MAX : iub / a[k]); |
---|
356 | *uu = (ilb == -DBL_MAX ? +DBL_MAX : ilb / a[k]); |
---|
357 | } |
---|
358 | else |
---|
359 | xassert(a != a); |
---|
360 | return; |
---|
361 | } |
---|
362 | |
---|
363 | /*********************************************************************** |
---|
364 | * check_row_bounds - check and relax original row bounds |
---|
365 | * |
---|
366 | * Given a row (constraint) |
---|
367 | * |
---|
368 | * n |
---|
369 | * L <= sum a[j] * x[j] <= U |
---|
370 | * j=1 |
---|
371 | * |
---|
372 | * and bounds of columns (variables) |
---|
373 | * |
---|
374 | * l[j] <= x[j] <= u[j] |
---|
375 | * |
---|
376 | * this routine checks the original row bounds L and U for feasibility |
---|
377 | * and redundancy. If the original lower bound L or/and upper bound U |
---|
378 | * cannot be active due to bounds of variables, the routine remove them |
---|
379 | * replacing by -inf or/and +inf, respectively. |
---|
380 | * |
---|
381 | * If no primal infeasibility is detected, the routine returns zero, |
---|
382 | * otherwise non-zero. */ |
---|
383 | |
---|
384 | static int check_row_bounds(const struct f_info *f, double *L_, |
---|
385 | double *U_) |
---|
386 | { int ret = 0; |
---|
387 | double L = *L_, U = *U_, LL, UU; |
---|
388 | /* determine implied bounds of the row */ |
---|
389 | row_implied_bounds(f, &LL, &UU); |
---|
390 | /* check if the original lower bound is infeasible */ |
---|
391 | if (L != -DBL_MAX) |
---|
392 | { double eps = 1e-3 * (1.0 + fabs(L)); |
---|
393 | if (UU < L - eps) |
---|
394 | { ret = 1; |
---|
395 | goto done; |
---|
396 | } |
---|
397 | } |
---|
398 | /* check if the original upper bound is infeasible */ |
---|
399 | if (U != +DBL_MAX) |
---|
400 | { double eps = 1e-3 * (1.0 + fabs(U)); |
---|
401 | if (LL > U + eps) |
---|
402 | { ret = 1; |
---|
403 | goto done; |
---|
404 | } |
---|
405 | } |
---|
406 | /* check if the original lower bound is redundant */ |
---|
407 | if (L != -DBL_MAX) |
---|
408 | { double eps = 1e-12 * (1.0 + fabs(L)); |
---|
409 | if (LL > L - eps) |
---|
410 | { /* it cannot be active, so remove it */ |
---|
411 | *L_ = -DBL_MAX; |
---|
412 | } |
---|
413 | } |
---|
414 | /* check if the original upper bound is redundant */ |
---|
415 | if (U != +DBL_MAX) |
---|
416 | { double eps = 1e-12 * (1.0 + fabs(U)); |
---|
417 | if (UU < U + eps) |
---|
418 | { /* it cannot be active, so remove it */ |
---|
419 | *U_ = +DBL_MAX; |
---|
420 | } |
---|
421 | } |
---|
422 | done: return ret; |
---|
423 | } |
---|
424 | |
---|
425 | /*********************************************************************** |
---|
426 | * check_col_bounds - check and tighten original column bounds |
---|
427 | * |
---|
428 | * Given a row (constraint) |
---|
429 | * |
---|
430 | * n |
---|
431 | * L <= sum a[j] * x[j] <= U |
---|
432 | * j=1 |
---|
433 | * |
---|
434 | * and bounds of columns (variables) |
---|
435 | * |
---|
436 | * l[j] <= x[j] <= u[j] |
---|
437 | * |
---|
438 | * for column (variable) x[j] this routine checks the original column |
---|
439 | * bounds l[j] and u[j] for feasibility and redundancy. If the original |
---|
440 | * lower bound l[j] or/and upper bound u[j] cannot be active due to |
---|
441 | * bounds of the constraint and other variables, the routine tighten |
---|
442 | * them replacing by corresponding implied bounds, if possible. |
---|
443 | * |
---|
444 | * NOTE: It is assumed that if L != -inf, the row lower bound can be |
---|
445 | * active, and if U != +inf, the row upper bound can be active. |
---|
446 | * |
---|
447 | * The flag means that variable x[j] is required to be integer. |
---|
448 | * |
---|
449 | * New actual bounds for x[j] are stored in locations lj and uj. |
---|
450 | * |
---|
451 | * If no primal infeasibility is detected, the routine returns zero, |
---|
452 | * otherwise non-zero. */ |
---|
453 | |
---|
454 | static int check_col_bounds(const struct f_info *f, int n, |
---|
455 | const double a[], double L, double U, const double l[], |
---|
456 | const double u[], int flag, int j, double *_lj, double *_uj) |
---|
457 | { int ret = 0; |
---|
458 | double lj, uj, ll, uu; |
---|
459 | xassert(n >= 0); |
---|
460 | xassert(1 <= j && j <= n); |
---|
461 | lj = l[j], uj = u[j]; |
---|
462 | /* determine implied bounds of the column */ |
---|
463 | col_implied_bounds(f, n, a, L, U, l, u, j, &ll, &uu); |
---|
464 | /* if x[j] is integral, round its implied bounds */ |
---|
465 | if (flag) |
---|
466 | { if (ll != -DBL_MAX) |
---|
467 | ll = (ll - floor(ll) < 1e-3 ? floor(ll) : ceil(ll)); |
---|
468 | if (uu != +DBL_MAX) |
---|
469 | uu = (ceil(uu) - uu < 1e-3 ? ceil(uu) : floor(uu)); |
---|
470 | } |
---|
471 | /* check if the original lower bound is infeasible */ |
---|
472 | if (lj != -DBL_MAX) |
---|
473 | { double eps = 1e-3 * (1.0 + fabs(lj)); |
---|
474 | if (uu < lj - eps) |
---|
475 | { ret = 1; |
---|
476 | goto done; |
---|
477 | } |
---|
478 | } |
---|
479 | /* check if the original upper bound is infeasible */ |
---|
480 | if (uj != +DBL_MAX) |
---|
481 | { double eps = 1e-3 * (1.0 + fabs(uj)); |
---|
482 | if (ll > uj + eps) |
---|
483 | { ret = 1; |
---|
484 | goto done; |
---|
485 | } |
---|
486 | } |
---|
487 | /* check if the original lower bound is redundant */ |
---|
488 | if (ll != -DBL_MAX) |
---|
489 | { double eps = 1e-3 * (1.0 + fabs(ll)); |
---|
490 | if (lj < ll - eps) |
---|
491 | { /* it cannot be active, so tighten it */ |
---|
492 | lj = ll; |
---|
493 | } |
---|
494 | } |
---|
495 | /* check if the original upper bound is redundant */ |
---|
496 | if (uu != +DBL_MAX) |
---|
497 | { double eps = 1e-3 * (1.0 + fabs(uu)); |
---|
498 | if (uj > uu + eps) |
---|
499 | { /* it cannot be active, so tighten it */ |
---|
500 | uj = uu; |
---|
501 | } |
---|
502 | } |
---|
503 | /* due to round-off errors it may happen that lj > uj (although |
---|
504 | lj < uj + eps, since no primal infeasibility is detected), so |
---|
505 | adjuct the new actual bounds to provide lj <= uj */ |
---|
506 | if (!(lj == -DBL_MAX || uj == +DBL_MAX)) |
---|
507 | { double t1 = fabs(lj), t2 = fabs(uj); |
---|
508 | double eps = 1e-10 * (1.0 + (t1 <= t2 ? t1 : t2)); |
---|
509 | if (lj > uj - eps) |
---|
510 | { if (lj == l[j]) |
---|
511 | uj = lj; |
---|
512 | else if (uj == u[j]) |
---|
513 | lj = uj; |
---|
514 | else if (t1 <= t2) |
---|
515 | uj = lj; |
---|
516 | else |
---|
517 | lj = uj; |
---|
518 | } |
---|
519 | } |
---|
520 | *_lj = lj, *_uj = uj; |
---|
521 | done: return ret; |
---|
522 | } |
---|
523 | |
---|
524 | /*********************************************************************** |
---|
525 | * check_efficiency - check if change in column bounds is efficient |
---|
526 | * |
---|
527 | * Given the original bounds of a column l and u and its new actual |
---|
528 | * bounds l' and u' (possibly tighten by the routine check_col_bounds) |
---|
529 | * this routine checks if the change in the column bounds is efficient |
---|
530 | * enough. If so, the routine returns non-zero, otherwise zero. |
---|
531 | * |
---|
532 | * The flag means that the variable is required to be integer. */ |
---|
533 | |
---|
534 | static int check_efficiency(int flag, double l, double u, double ll, |
---|
535 | double uu) |
---|
536 | { int eff = 0; |
---|
537 | /* check efficiency for lower bound */ |
---|
538 | if (l < ll) |
---|
539 | { if (flag || l == -DBL_MAX) |
---|
540 | eff++; |
---|
541 | else |
---|
542 | { double r; |
---|
543 | if (u == +DBL_MAX) |
---|
544 | r = 1.0 + fabs(l); |
---|
545 | else |
---|
546 | r = 1.0 + (u - l); |
---|
547 | if (ll - l >= 0.25 * r) |
---|
548 | eff++; |
---|
549 | } |
---|
550 | } |
---|
551 | /* check efficiency for upper bound */ |
---|
552 | if (u > uu) |
---|
553 | { if (flag || u == +DBL_MAX) |
---|
554 | eff++; |
---|
555 | else |
---|
556 | { double r; |
---|
557 | if (l == -DBL_MAX) |
---|
558 | r = 1.0 + fabs(u); |
---|
559 | else |
---|
560 | r = 1.0 + (u - l); |
---|
561 | if (u - uu >= 0.25 * r) |
---|
562 | eff++; |
---|
563 | } |
---|
564 | } |
---|
565 | return eff; |
---|
566 | } |
---|
567 | |
---|
568 | /*********************************************************************** |
---|
569 | * basic_preprocessing - perform basic preprocessing |
---|
570 | * |
---|
571 | * This routine performs basic preprocessing of the specified MIP that |
---|
572 | * includes relaxing some row bounds and tightening some column bounds. |
---|
573 | * |
---|
574 | * On entry the arrays L and U contains original row bounds, and the |
---|
575 | * arrays l and u contains original column bounds: |
---|
576 | * |
---|
577 | * L[0] is the lower bound of the objective row; |
---|
578 | * L[i], i = 1,...,m, is the lower bound of i-th row; |
---|
579 | * U[0] is the upper bound of the objective row; |
---|
580 | * U[i], i = 1,...,m, is the upper bound of i-th row; |
---|
581 | * l[0] is not used; |
---|
582 | * l[j], j = 1,...,n, is the lower bound of j-th column; |
---|
583 | * u[0] is not used; |
---|
584 | * u[j], j = 1,...,n, is the upper bound of j-th column. |
---|
585 | * |
---|
586 | * On exit the arrays L, U, l, and u contain new actual bounds of rows |
---|
587 | * and column in the same locations. |
---|
588 | * |
---|
589 | * The parameters nrs and num specify an initial list of rows to be |
---|
590 | * processed: |
---|
591 | * |
---|
592 | * nrs is the number of rows in the initial list, 0 <= nrs <= m+1; |
---|
593 | * num[0] is not used; |
---|
594 | * num[1,...,nrs] are row numbers (0 means the objective row). |
---|
595 | * |
---|
596 | * The parameter max_pass specifies the maximal number of times that |
---|
597 | * each row can be processed, max_pass > 0. |
---|
598 | * |
---|
599 | * If no primal infeasibility is detected, the routine returns zero, |
---|
600 | * otherwise non-zero. */ |
---|
601 | |
---|
602 | static int basic_preprocessing(glp_prob *mip, double L[], double U[], |
---|
603 | double l[], double u[], int nrs, const int num[], int max_pass) |
---|
604 | { int m = mip->m; |
---|
605 | int n = mip->n; |
---|
606 | struct f_info f; |
---|
607 | int i, j, k, len, size, ret = 0; |
---|
608 | int *ind, *list, *mark, *pass; |
---|
609 | double *val, *lb, *ub; |
---|
610 | xassert(0 <= nrs && nrs <= m+1); |
---|
611 | xassert(max_pass > 0); |
---|
612 | /* allocate working arrays */ |
---|
613 | ind = xcalloc(1+n, sizeof(int)); |
---|
614 | list = xcalloc(1+m+1, sizeof(int)); |
---|
615 | mark = xcalloc(1+m+1, sizeof(int)); |
---|
616 | memset(&mark[0], 0, (m+1) * sizeof(int)); |
---|
617 | pass = xcalloc(1+m+1, sizeof(int)); |
---|
618 | memset(&pass[0], 0, (m+1) * sizeof(int)); |
---|
619 | val = xcalloc(1+n, sizeof(double)); |
---|
620 | lb = xcalloc(1+n, sizeof(double)); |
---|
621 | ub = xcalloc(1+n, sizeof(double)); |
---|
622 | /* initialize the list of rows to be processed */ |
---|
623 | size = 0; |
---|
624 | for (k = 1; k <= nrs; k++) |
---|
625 | { i = num[k]; |
---|
626 | xassert(0 <= i && i <= m); |
---|
627 | /* duplicate row numbers are not allowed */ |
---|
628 | xassert(!mark[i]); |
---|
629 | list[++size] = i, mark[i] = 1; |
---|
630 | } |
---|
631 | xassert(size == nrs); |
---|
632 | /* process rows in the list until it becomes empty */ |
---|
633 | while (size > 0) |
---|
634 | { /* get a next row from the list */ |
---|
635 | i = list[size--], mark[i] = 0; |
---|
636 | /* increase the row processing count */ |
---|
637 | pass[i]++; |
---|
638 | /* if the row is free, skip it */ |
---|
639 | if (L[i] == -DBL_MAX && U[i] == +DBL_MAX) continue; |
---|
640 | /* obtain coefficients of the row */ |
---|
641 | len = 0; |
---|
642 | if (i == 0) |
---|
643 | { for (j = 1; j <= n; j++) |
---|
644 | { GLPCOL *col = mip->col[j]; |
---|
645 | if (col->coef != 0.0) |
---|
646 | len++, ind[len] = j, val[len] = col->coef; |
---|
647 | } |
---|
648 | } |
---|
649 | else |
---|
650 | { GLPROW *row = mip->row[i]; |
---|
651 | GLPAIJ *aij; |
---|
652 | for (aij = row->ptr; aij != NULL; aij = aij->r_next) |
---|
653 | len++, ind[len] = aij->col->j, val[len] = aij->val; |
---|
654 | } |
---|
655 | /* determine lower and upper bounds of columns corresponding |
---|
656 | to non-zero row coefficients */ |
---|
657 | for (k = 1; k <= len; k++) |
---|
658 | j = ind[k], lb[k] = l[j], ub[k] = u[j]; |
---|
659 | /* prepare the row info to determine implied bounds */ |
---|
660 | prepare_row_info(len, val, lb, ub, &f); |
---|
661 | /* check and relax bounds of the row */ |
---|
662 | if (check_row_bounds(&f, &L[i], &U[i])) |
---|
663 | { /* the feasible region is empty */ |
---|
664 | ret = 1; |
---|
665 | goto done; |
---|
666 | } |
---|
667 | /* if the row became free, drop it */ |
---|
668 | if (L[i] == -DBL_MAX && U[i] == +DBL_MAX) continue; |
---|
669 | /* process columns having non-zero coefficients in the row */ |
---|
670 | for (k = 1; k <= len; k++) |
---|
671 | { GLPCOL *col; |
---|
672 | int flag, eff; |
---|
673 | double ll, uu; |
---|
674 | /* take a next column in the row */ |
---|
675 | j = ind[k], col = mip->col[j]; |
---|
676 | flag = col->kind != GLP_CV; |
---|
677 | /* check and tighten bounds of the column */ |
---|
678 | if (check_col_bounds(&f, len, val, L[i], U[i], lb, ub, |
---|
679 | flag, k, &ll, &uu)) |
---|
680 | { /* the feasible region is empty */ |
---|
681 | ret = 1; |
---|
682 | goto done; |
---|
683 | } |
---|
684 | /* check if change in the column bounds is efficient */ |
---|
685 | eff = check_efficiency(flag, l[j], u[j], ll, uu); |
---|
686 | /* set new actual bounds of the column */ |
---|
687 | l[j] = ll, u[j] = uu; |
---|
688 | /* if the change is efficient, add all rows affected by the |
---|
689 | corresponding column, to the list */ |
---|
690 | if (eff > 0) |
---|
691 | { GLPAIJ *aij; |
---|
692 | for (aij = col->ptr; aij != NULL; aij = aij->c_next) |
---|
693 | { int ii = aij->row->i; |
---|
694 | /* if the row was processed maximal number of times, |
---|
695 | skip it */ |
---|
696 | if (pass[ii] >= max_pass) continue; |
---|
697 | /* if the row is free, skip it */ |
---|
698 | if (L[ii] == -DBL_MAX && U[ii] == +DBL_MAX) continue; |
---|
699 | /* put the row into the list */ |
---|
700 | if (mark[ii] == 0) |
---|
701 | { xassert(size <= m); |
---|
702 | list[++size] = ii, mark[ii] = 1; |
---|
703 | } |
---|
704 | } |
---|
705 | } |
---|
706 | } |
---|
707 | } |
---|
708 | done: /* free working arrays */ |
---|
709 | xfree(ind); |
---|
710 | xfree(list); |
---|
711 | xfree(mark); |
---|
712 | xfree(pass); |
---|
713 | xfree(val); |
---|
714 | xfree(lb); |
---|
715 | xfree(ub); |
---|
716 | return ret; |
---|
717 | } |
---|
718 | |
---|
719 | /*********************************************************************** |
---|
720 | * NAME |
---|
721 | * |
---|
722 | * ios_preprocess_node - preprocess current subproblem |
---|
723 | * |
---|
724 | * SYNOPSIS |
---|
725 | * |
---|
726 | * #include "glpios.h" |
---|
727 | * int ios_preprocess_node(glp_tree *tree, int max_pass); |
---|
728 | * |
---|
729 | * DESCRIPTION |
---|
730 | * |
---|
731 | * The routine ios_preprocess_node performs basic preprocessing of the |
---|
732 | * current subproblem. |
---|
733 | * |
---|
734 | * RETURNS |
---|
735 | * |
---|
736 | * If no primal infeasibility is detected, the routine returns zero, |
---|
737 | * otherwise non-zero. */ |
---|
738 | |
---|
739 | int ios_preprocess_node(glp_tree *tree, int max_pass) |
---|
740 | { glp_prob *mip = tree->mip; |
---|
741 | int m = mip->m; |
---|
742 | int n = mip->n; |
---|
743 | int i, j, nrs, *num, ret = 0; |
---|
744 | double *L, *U, *l, *u; |
---|
745 | /* the current subproblem must exist */ |
---|
746 | xassert(tree->curr != NULL); |
---|
747 | /* determine original row bounds */ |
---|
748 | L = xcalloc(1+m, sizeof(double)); |
---|
749 | U = xcalloc(1+m, sizeof(double)); |
---|
750 | switch (mip->mip_stat) |
---|
751 | { case GLP_UNDEF: |
---|
752 | L[0] = -DBL_MAX, U[0] = +DBL_MAX; |
---|
753 | break; |
---|
754 | case GLP_FEAS: |
---|
755 | switch (mip->dir) |
---|
756 | { case GLP_MIN: |
---|
757 | L[0] = -DBL_MAX, U[0] = mip->mip_obj - mip->c0; |
---|
758 | break; |
---|
759 | case GLP_MAX: |
---|
760 | L[0] = mip->mip_obj - mip->c0, U[0] = +DBL_MAX; |
---|
761 | break; |
---|
762 | default: |
---|
763 | xassert(mip != mip); |
---|
764 | } |
---|
765 | break; |
---|
766 | default: |
---|
767 | xassert(mip != mip); |
---|
768 | } |
---|
769 | for (i = 1; i <= m; i++) |
---|
770 | { L[i] = glp_get_row_lb(mip, i); |
---|
771 | U[i] = glp_get_row_ub(mip, i); |
---|
772 | } |
---|
773 | /* determine original column bounds */ |
---|
774 | l = xcalloc(1+n, sizeof(double)); |
---|
775 | u = xcalloc(1+n, sizeof(double)); |
---|
776 | for (j = 1; j <= n; j++) |
---|
777 | { l[j] = glp_get_col_lb(mip, j); |
---|
778 | u[j] = glp_get_col_ub(mip, j); |
---|
779 | } |
---|
780 | /* build the initial list of rows to be analyzed */ |
---|
781 | nrs = m + 1; |
---|
782 | num = xcalloc(1+nrs, sizeof(int)); |
---|
783 | for (i = 1; i <= nrs; i++) num[i] = i - 1; |
---|
784 | /* perform basic preprocessing */ |
---|
785 | if (basic_preprocessing(mip , L, U, l, u, nrs, num, max_pass)) |
---|
786 | { ret = 1; |
---|
787 | goto done; |
---|
788 | } |
---|
789 | /* set new actual (relaxed) row bounds */ |
---|
790 | for (i = 1; i <= m; i++) |
---|
791 | { /* consider only non-active rows to keep dual feasibility */ |
---|
792 | if (glp_get_row_stat(mip, i) == GLP_BS) |
---|
793 | { if (L[i] == -DBL_MAX && U[i] == +DBL_MAX) |
---|
794 | glp_set_row_bnds(mip, i, GLP_FR, 0.0, 0.0); |
---|
795 | else if (U[i] == +DBL_MAX) |
---|
796 | glp_set_row_bnds(mip, i, GLP_LO, L[i], 0.0); |
---|
797 | else if (L[i] == -DBL_MAX) |
---|
798 | glp_set_row_bnds(mip, i, GLP_UP, 0.0, U[i]); |
---|
799 | } |
---|
800 | } |
---|
801 | /* set new actual (tightened) column bounds */ |
---|
802 | for (j = 1; j <= n; j++) |
---|
803 | { int type; |
---|
804 | if (l[j] == -DBL_MAX && u[j] == +DBL_MAX) |
---|
805 | type = GLP_FR; |
---|
806 | else if (u[j] == +DBL_MAX) |
---|
807 | type = GLP_LO; |
---|
808 | else if (l[j] == -DBL_MAX) |
---|
809 | type = GLP_UP; |
---|
810 | else if (l[j] != u[j]) |
---|
811 | type = GLP_DB; |
---|
812 | else |
---|
813 | type = GLP_FX; |
---|
814 | glp_set_col_bnds(mip, j, type, l[j], u[j]); |
---|
815 | } |
---|
816 | done: /* free working arrays and return */ |
---|
817 | xfree(L); |
---|
818 | xfree(U); |
---|
819 | xfree(l); |
---|
820 | xfree(u); |
---|
821 | xfree(num); |
---|
822 | return ret; |
---|
823 | } |
---|
824 | |
---|
825 | /* eof */ |
---|