COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glpios02.c @ 1:c445c931472f

Last change on this file since 1:c445c931472f was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 10 years ago

Import glpk-4.45

  • Generated files and doc/notes are removed
File size: 26.3 KB
Line 
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
97struct f_info
98{     int j_min, j_max;
99      double f_min, f_max;
100};
101
102static 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
211static 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
297static 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
384static 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      }
422done: 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
454static 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;
521done: 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
534static 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
602static 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      }
708done: /* 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
739int 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      }
816done: /* 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 */
Note: See TracBrowser for help on using the repository browser.