|
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 */ |