lemon-project-template-glpk

comparison deps/glpk/src/glpios02.c @ 9:33de93886c88

find changesets by author, revision, files, or words in the commit message
Import GLPK 4.47
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 20:59:10 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:7f71bb379e88
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, 2011 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 */