1 /* glpios02.c (preprocess current subproblem) */
3 /***********************************************************************
4 * This code is part of GLPK (GNU Linear Programming Kit).
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>.
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.
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.
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 ***********************************************************************/
27 /***********************************************************************
28 * prepare_row_info - prepare row info to determine implied bounds
30 * Given a row (linear form)
36 * and bounds of columns (variables)
38 * l[j] <= x[j] <= u[j] (2)
40 * this routine computes f_min, j_min, f_max, j_max needed to determine
45 * Let J+ = {j : a[j] > 0} and J- = {j : a[j] < 0}.
47 * Parameters f_min and j_min are computed as follows:
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
52 * f_min := sum a[j] * l[j] + sum a[j] * u[j]
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
60 * f_min := sum a[j] * l[j] + sum a[j] * u[j]
61 * j in J+\{k} j in J-\{k}
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
72 * Parameters f_max and j_max are computed in a similar way as follows:
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
77 * f_max := sum a[j] * u[j] + sum a[j] * l[j]
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
85 * f_max := sum a[j] * u[j] + sum a[j] * l[j]
86 * j in J+\{k} j in J-\{k}
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
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;
107 /* determine f_min and j_min */
108 f_min = 0.0, j_min = 0;
109 for (j = 1; j <= n; j++)
111 { if (l[j] == -DBL_MAX)
115 { f_min = -DBL_MAX, j_min = 0;
120 f_min += a[j] * l[j];
123 { if (u[j] == +DBL_MAX)
127 { f_min = -DBL_MAX, j_min = 0;
132 f_min += a[j] * u[j];
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++)
142 { if (u[j] == +DBL_MAX)
146 { f_max = +DBL_MAX, j_max = 0;
151 f_max += a[j] * u[j];
154 { if (l[j] == -DBL_MAX)
158 { f_max = +DBL_MAX, j_max = 0;
163 f_max += a[j] * l[j];
168 f->f_max = f_max, f->j_max = j_max;
172 /***********************************************************************
173 * row_implied_bounds - determine row implied bounds
175 * Given a row (linear form)
181 * and bounds of columns (variables)
183 * l[j] <= x[j] <= u[j]
185 * this routine determines implied bounds of the row.
189 * Let J+ = {j : a[j] > 0} and J- = {j : a[j] < 0}.
191 * The implied lower bound of the row is computed as follows:
193 * L' := sum a[j] * l[j] + sum a[j] * u[j] (9)
196 * and as it follows from (3), (4), and (5):
198 * L' := if j_min = 0 then f_min else -inf (10)
200 * The implied upper bound of the row is computed as follows:
202 * U' := sum a[j] * u[j] + sum a[j] * l[j] (11)
205 * and as it follows from (6), (7), and (8):
207 * U' := if j_max = 0 then f_max else +inf (12)
209 * The implied bounds are stored in locations LL and UU. */
211 static void row_implied_bounds(const struct f_info *f, double *LL,
213 { *LL = (f->j_min == 0 ? f->f_min : -DBL_MAX);
214 *UU = (f->j_max == 0 ? f->f_max : +DBL_MAX);
218 /***********************************************************************
219 * col_implied_bounds - determine column implied bounds
221 * Given a row (constraint)
224 * L <= sum a[j] * x[j] <= U (13)
227 * and bounds of columns (variables)
229 * l[j] <= x[j] <= u[j]
231 * this routine determines implied bounds of variable x[k].
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.
238 * From (13) it follows that
240 * L <= sum a[j] * x[j] + a[k] * x[k] <= U
244 * L - sum a[j] * x[j] <= a[k] * x[k] <= U - sum a[j] * x[j]
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:
250 * ilb(a[k] * x[k]) = min(L - sum a[j] * x[j]) =
253 * = L - max sum a[j] * x[j]
256 * where, as it follows from (6), (7), and (8)
258 * / f_max - a[k] * u[k], j_max = 0, a[k] > 0
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
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:
269 * iub(a[k] * x[k]) = max(U - sum a[j] * x[j]) =
272 * = U - min sum a[j] * x[j]
275 * where, as it follows from (3), (4), and (5)
277 * / f_min - a[k] * l[k], j_min = 0, a[k] > 0
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
287 * ilb(a[k] * x[k]) <= a[k] * x[k] <= iub(a[k] * x[k])
289 * implied lower and upper bounds of x[k] are determined as follows:
291 * l'[k] := if a[k] > 0 then ilb / a[k] else ulb / a[k] (16)
293 * u'[k] := if a[k] > 0 then ulb / a[k] else ilb / a[k] (17)
295 * The implied bounds are stored in locations ll and uu. */
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)
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)
306 else if (f->j_max == 0)
308 { xassert(u[k] != +DBL_MAX);
309 ilb = L - (f->f_max - a[k] * u[k]);
312 { xassert(l[k] != -DBL_MAX);
313 ilb = L - (f->f_max - a[k] * l[k]);
318 else if (f->j_max == k)
322 /* determine implied upper bound of term a[k] * x[k] (15) */
323 if (U == +DBL_MAX || f->f_min == -DBL_MAX)
325 else if (f->j_min == 0)
327 { xassert(l[k] != -DBL_MAX);
328 iub = U - (f->f_min - a[k] * l[k]);
331 { xassert(u[k] != +DBL_MAX);
332 iub = U - (f->f_min - a[k] * u[k]);
337 else if (f->j_min == k)
341 /* determine implied bounds of x[k] (16) and (17) */
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
347 if (fabs(a[k]) < 1e-6)
348 *ll = -DBL_MAX, *uu = +DBL_MAX; else
351 { *ll = (ilb == -DBL_MAX ? -DBL_MAX : ilb / a[k]);
352 *uu = (iub == +DBL_MAX ? +DBL_MAX : iub / a[k]);
355 { *ll = (iub == +DBL_MAX ? -DBL_MAX : iub / a[k]);
356 *uu = (ilb == -DBL_MAX ? +DBL_MAX : ilb / a[k]);
363 /***********************************************************************
364 * check_row_bounds - check and relax original row bounds
366 * Given a row (constraint)
369 * L <= sum a[j] * x[j] <= U
372 * and bounds of columns (variables)
374 * l[j] <= x[j] <= u[j]
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.
381 * If no primal infeasibility is detected, the routine returns zero,
382 * otherwise non-zero. */
384 static int check_row_bounds(const struct f_info *f, double *L_,
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 */
392 { double eps = 1e-3 * (1.0 + fabs(L));
398 /* check if the original upper bound is infeasible */
400 { double eps = 1e-3 * (1.0 + fabs(U));
406 /* check if the original lower bound is redundant */
408 { double eps = 1e-12 * (1.0 + fabs(L));
410 { /* it cannot be active, so remove it */
414 /* check if the original upper bound is redundant */
416 { double eps = 1e-12 * (1.0 + fabs(U));
418 { /* it cannot be active, so remove it */
425 /***********************************************************************
426 * check_col_bounds - check and tighten original column bounds
428 * Given a row (constraint)
431 * L <= sum a[j] * x[j] <= U
434 * and bounds of columns (variables)
436 * l[j] <= x[j] <= u[j]
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.
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.
447 * The flag means that variable x[j] is required to be integer.
449 * New actual bounds for x[j] are stored in locations lj and uj.
451 * If no primal infeasibility is detected, the routine returns zero,
452 * otherwise non-zero. */
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)
458 double lj, uj, ll, uu;
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 */
466 { if (ll != -DBL_MAX)
467 ll = (ll - floor(ll) < 1e-3 ? floor(ll) : ceil(ll));
469 uu = (ceil(uu) - uu < 1e-3 ? ceil(uu) : floor(uu));
471 /* check if the original lower bound is infeasible */
473 { double eps = 1e-3 * (1.0 + fabs(lj));
479 /* check if the original upper bound is infeasible */
481 { double eps = 1e-3 * (1.0 + fabs(uj));
487 /* check if the original lower bound is redundant */
489 { double eps = 1e-3 * (1.0 + fabs(ll));
491 { /* it cannot be active, so tighten it */
495 /* check if the original upper bound is redundant */
497 { double eps = 1e-3 * (1.0 + fabs(uu));
499 { /* it cannot be active, so tighten it */
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));
520 *_lj = lj, *_uj = uj;
524 /***********************************************************************
525 * check_efficiency - check if change in column bounds is efficient
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.
532 * The flag means that the variable is required to be integer. */
534 static int check_efficiency(int flag, double l, double u, double ll,
537 /* check efficiency for lower bound */
539 { if (flag || l == -DBL_MAX)
547 if (ll - l >= 0.25 * r)
551 /* check efficiency for upper bound */
553 { if (flag || u == +DBL_MAX)
561 if (u - uu >= 0.25 * r)
568 /***********************************************************************
569 * basic_preprocessing - perform basic preprocessing
571 * This routine performs basic preprocessing of the specified MIP that
572 * includes relaxing some row bounds and tightening some column bounds.
574 * On entry the arrays L and U contains original row bounds, and the
575 * arrays l and u contains original column bounds:
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;
582 * l[j], j = 1,...,n, is the lower bound of j-th column;
584 * u[j], j = 1,...,n, is the upper bound of j-th column.
586 * On exit the arrays L, U, l, and u contain new actual bounds of rows
587 * and column in the same locations.
589 * The parameters nrs and num specify an initial list of rows to be
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).
596 * The parameter max_pass specifies the maximal number of times that
597 * each row can be processed, max_pass > 0.
599 * If no primal infeasibility is detected, the routine returns zero,
600 * otherwise non-zero. */
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)
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 */
624 for (k = 1; k <= nrs; k++)
626 xassert(0 <= i && i <= m);
627 /* duplicate row numbers are not allowed */
629 list[++size] = i, mark[i] = 1;
631 xassert(size == nrs);
632 /* process rows in the list until it becomes empty */
634 { /* get a next row from the list */
635 i = list[size--], mark[i] = 0;
636 /* increase the row processing count */
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 */
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;
650 { GLPROW *row = mip->row[i];
652 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
653 len++, ind[len] = aij->col->j, val[len] = aij->val;
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 */
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++)
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,
680 { /* the feasible region is empty */
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 */
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,
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 */
701 { xassert(size <= m);
702 list[++size] = ii, mark[ii] = 1;
708 done: /* free working arrays */
719 /***********************************************************************
722 * ios_preprocess_node - preprocess current subproblem
726 * #include "glpios.h"
727 * int ios_preprocess_node(glp_tree *tree, int max_pass);
731 * The routine ios_preprocess_node performs basic preprocessing of the
732 * current subproblem.
736 * If no primal infeasibility is detected, the routine returns zero,
737 * otherwise non-zero. */
739 int ios_preprocess_node(glp_tree *tree, int max_pass)
740 { glp_prob *mip = tree->mip;
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)
752 L[0] = -DBL_MAX, U[0] = +DBL_MAX;
757 L[0] = -DBL_MAX, U[0] = mip->mip_obj - mip->c0;
760 L[0] = mip->mip_obj - mip->c0, U[0] = +DBL_MAX;
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);
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);
780 /* build the initial list of rows to be analyzed */
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))
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]);
801 /* set new actual (tightened) column bounds */
802 for (j = 1; j <= n; j++)
804 if (l[j] == -DBL_MAX && u[j] == +DBL_MAX)
806 else if (u[j] == +DBL_MAX)
808 else if (l[j] == -DBL_MAX)
810 else if (l[j] != u[j])
814 glp_set_col_bnds(mip, j, type, l[j], u[j]);
816 done: /* free working arrays and return */