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