COIN-OR::LEMON - Graph Library

source: lemon-project-template-glpk/deps/glpk/src/glpios06.c @ 10:5545663ca997

subpack-glpk
Last change on this file since 10:5545663ca997 was 9:33de93886c88, checked in by Alpar Juttner <alpar@…>, 13 years ago

Import GLPK 4.47

File size: 47.8 KB
Line 
1/* glpios06.c (MIR cut generator) */
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#define _MIR_DEBUG 0
28
29#define MAXAGGR 5
30/* maximal number of rows which can be aggregated */
31
32struct MIR
33{     /* MIR cut generator working area */
34      /*--------------------------------------------------------------*/
35      /* global information valid for the root subproblem */
36      int m;
37      /* number of rows (in the root subproblem) */
38      int n;
39      /* number of columns */
40      char *skip; /* char skip[1+m]; */
41      /* skip[i], 1 <= i <= m, is a flag that means that row i should
42         not be used because (1) it is not suitable, or (2) because it
43         has been used in the aggregated constraint */
44      char *isint; /* char isint[1+m+n]; */
45      /* isint[k], 1 <= k <= m+n, is a flag that means that variable
46         x[k] is integer (otherwise, continuous) */
47      double *lb; /* double lb[1+m+n]; */
48      /* lb[k], 1 <= k <= m+n, is lower bound of x[k]; -DBL_MAX means
49         that x[k] has no lower bound */
50      int *vlb; /* int vlb[1+m+n]; */
51      /* vlb[k] = k', 1 <= k <= m+n, is the number of integer variable,
52         which defines variable lower bound x[k] >= lb[k] * x[k']; zero
53         means that x[k] has simple lower bound */
54      double *ub; /* double ub[1+m+n]; */
55      /* ub[k], 1 <= k <= m+n, is upper bound of x[k]; +DBL_MAX means
56         that x[k] has no upper bound */
57      int *vub; /* int vub[1+m+n]; */
58      /* vub[k] = k', 1 <= k <= m+n, is the number of integer variable,
59         which defines variable upper bound x[k] <= ub[k] * x[k']; zero
60         means that x[k] has simple upper bound */
61      /*--------------------------------------------------------------*/
62      /* current (fractional) point to be separated */
63      double *x; /* double x[1+m+n]; */
64      /* x[k] is current value of auxiliary (1 <= k <= m) or structural
65         (m+1 <= k <= m+n) variable */
66      /*--------------------------------------------------------------*/
67      /* aggregated constraint sum a[k] * x[k] = b, which is a linear
68         combination of original constraints transformed to equalities
69         by introducing auxiliary variables */
70      int agg_cnt;
71      /* number of rows (original constraints) used to build aggregated
72         constraint, 1 <= agg_cnt <= MAXAGGR */
73      int *agg_row; /* int agg_row[1+MAXAGGR]; */
74      /* agg_row[k], 1 <= k <= agg_cnt, is the row number used to build
75         aggregated constraint */
76      IOSVEC *agg_vec; /* IOSVEC agg_vec[1:m+n]; */
77      /* sparse vector of aggregated constraint coefficients, a[k] */
78      double agg_rhs;
79      /* right-hand side of the aggregated constraint, b */
80      /*--------------------------------------------------------------*/
81      /* bound substitution flags for modified constraint */
82      char *subst; /* char subst[1+m+n]; */
83      /* subst[k], 1 <= k <= m+n, is a bound substitution flag used for
84         variable x[k]:
85         '?' - x[k] is missing in modified constraint
86         'L' - x[k] = (lower bound) + x'[k]
87         'U' - x[k] = (upper bound) - x'[k] */
88      /*--------------------------------------------------------------*/
89      /* modified constraint sum a'[k] * x'[k] = b', where x'[k] >= 0,
90         derived from aggregated constraint by substituting bounds;
91         note that due to substitution of variable bounds there may be
92         additional terms in the modified constraint */
93      IOSVEC *mod_vec; /* IOSVEC mod_vec[1:m+n]; */
94      /* sparse vector of modified constraint coefficients, a'[k] */
95      double mod_rhs;
96      /* right-hand side of the modified constraint, b' */
97      /*--------------------------------------------------------------*/
98      /* cutting plane sum alpha[k] * x[k] <= beta */
99      IOSVEC *cut_vec; /* IOSVEC cut_vec[1:m+n]; */
100      /* sparse vector of cutting plane coefficients, alpha[k] */
101      double cut_rhs;
102      /* right-hand size of the cutting plane, beta */
103};
104
105/***********************************************************************
106*  NAME
107*
108*  ios_mir_init - initialize MIR cut generator
109*
110*  SYNOPSIS
111*
112*  #include "glpios.h"
113*  void *ios_mir_init(glp_tree *tree);
114*
115*  DESCRIPTION
116*
117*  The routine ios_mir_init initializes the MIR cut generator assuming
118*  that the current subproblem is the root subproblem.
119*
120*  RETURNS
121*
122*  The routine ios_mir_init returns a pointer to the MIR cut generator
123*  working area. */
124
125static void set_row_attrib(glp_tree *tree, struct MIR *mir)
126{     /* set global row attributes */
127      glp_prob *mip = tree->mip;
128      int m = mir->m;
129      int k;
130      for (k = 1; k <= m; k++)
131      {  GLPROW *row = mip->row[k];
132         mir->skip[k] = 0;
133         mir->isint[k] = 0;
134         switch (row->type)
135         {  case GLP_FR:
136               mir->lb[k] = -DBL_MAX, mir->ub[k] = +DBL_MAX; break;
137            case GLP_LO:
138               mir->lb[k] = row->lb, mir->ub[k] = +DBL_MAX; break;
139            case GLP_UP:
140               mir->lb[k] = -DBL_MAX, mir->ub[k] = row->ub; break;
141            case GLP_DB:
142               mir->lb[k] = row->lb, mir->ub[k] = row->ub; break;
143            case GLP_FX:
144               mir->lb[k] = mir->ub[k] = row->lb; break;
145            default:
146               xassert(row != row);
147         }
148         mir->vlb[k] = mir->vub[k] = 0;
149      }
150      return;
151}
152
153static void set_col_attrib(glp_tree *tree, struct MIR *mir)
154{     /* set global column attributes */
155      glp_prob *mip = tree->mip;
156      int m = mir->m;
157      int n = mir->n;
158      int k;
159      for (k = m+1; k <= m+n; k++)
160      {  GLPCOL *col = mip->col[k-m];
161         switch (col->kind)
162         {  case GLP_CV:
163               mir->isint[k] = 0; break;
164            case GLP_IV:
165               mir->isint[k] = 1; break;
166            default:
167               xassert(col != col);
168         }
169         switch (col->type)
170         {  case GLP_FR:
171               mir->lb[k] = -DBL_MAX, mir->ub[k] = +DBL_MAX; break;
172            case GLP_LO:
173               mir->lb[k] = col->lb, mir->ub[k] = +DBL_MAX; break;
174            case GLP_UP:
175               mir->lb[k] = -DBL_MAX, mir->ub[k] = col->ub; break;
176            case GLP_DB:
177               mir->lb[k] = col->lb, mir->ub[k] = col->ub; break;
178            case GLP_FX:
179               mir->lb[k] = mir->ub[k] = col->lb; break;
180            default:
181               xassert(col != col);
182         }
183         mir->vlb[k] = mir->vub[k] = 0;
184      }
185      return;
186}
187
188static void set_var_bounds(glp_tree *tree, struct MIR *mir)
189{     /* set variable bounds */
190      glp_prob *mip = tree->mip;
191      int m = mir->m;
192      GLPAIJ *aij;
193      int i, k1, k2;
194      double a1, a2;
195      for (i = 1; i <= m; i++)
196      {  /* we need the row to be '>= 0' or '<= 0' */
197         if (!(mir->lb[i] == 0.0 && mir->ub[i] == +DBL_MAX ||
198               mir->lb[i] == -DBL_MAX && mir->ub[i] == 0.0)) continue;
199         /* take first term */
200         aij = mip->row[i]->ptr;
201         if (aij == NULL) continue;
202         k1 = m + aij->col->j, a1 = aij->val;
203         /* take second term */
204         aij = aij->r_next;
205         if (aij == NULL) continue;
206         k2 = m + aij->col->j, a2 = aij->val;
207         /* there must be only two terms */
208         if (aij->r_next != NULL) continue;
209         /* interchange terms, if needed */
210         if (!mir->isint[k1] && mir->isint[k2])
211            ;
212         else if (mir->isint[k1] && !mir->isint[k2])
213         {  k2 = k1, a2 = a1;
214            k1 = m + aij->col->j, a1 = aij->val;
215         }
216         else
217         {  /* both terms are either continuous or integer */
218            continue;
219         }
220         /* x[k2] should be double-bounded */
221         if (mir->lb[k2] == -DBL_MAX || mir->ub[k2] == +DBL_MAX ||
222             mir->lb[k2] == mir->ub[k2]) continue;
223         /* change signs, if necessary */
224         if (mir->ub[i] == 0.0) a1 = - a1, a2 = - a2;
225         /* now the row has the form a1 * x1 + a2 * x2 >= 0, where x1
226            is continuous, x2 is integer */
227         if (a1 > 0.0)
228         {  /* x1 >= - (a2 / a1) * x2 */
229            if (mir->vlb[k1] == 0)
230            {  /* set variable lower bound for x1 */
231               mir->lb[k1] = - a2 / a1;
232               mir->vlb[k1] = k2;
233               /* the row should not be used */
234               mir->skip[i] = 1;
235            }
236         }
237         else /* a1 < 0.0 */
238         {  /* x1 <= - (a2 / a1) * x2 */
239            if (mir->vub[k1] == 0)
240            {  /* set variable upper bound for x1 */
241               mir->ub[k1] = - a2 / a1;
242               mir->vub[k1] = k2;
243               /* the row should not be used */
244               mir->skip[i] = 1;
245            }
246         }
247      }
248      return;
249}
250
251static void mark_useless_rows(glp_tree *tree, struct MIR *mir)
252{     /* mark rows which should not be used */
253      glp_prob *mip = tree->mip;
254      int m = mir->m;
255      GLPAIJ *aij;
256      int i, k, nv;
257      for (i = 1; i <= m; i++)
258      {  /* free rows should not be used */
259         if (mir->lb[i] == -DBL_MAX && mir->ub[i] == +DBL_MAX)
260         {  mir->skip[i] = 1;
261            continue;
262         }
263         nv = 0;
264         for (aij = mip->row[i]->ptr; aij != NULL; aij = aij->r_next)
265         {  k = m + aij->col->j;
266            /* rows with free variables should not be used */
267            if (mir->lb[k] == -DBL_MAX && mir->ub[k] == +DBL_MAX)
268            {  mir->skip[i] = 1;
269               break;
270            }
271            /* rows with integer variables having infinite (lower or
272               upper) bound should not be used */
273            if (mir->isint[k] && mir->lb[k] == -DBL_MAX ||
274                mir->isint[k] && mir->ub[k] == +DBL_MAX)
275            {  mir->skip[i] = 1;
276               break;
277            }
278            /* count non-fixed variables */
279            if (!(mir->vlb[k] == 0 && mir->vub[k] == 0 &&
280                  mir->lb[k] == mir->ub[k])) nv++;
281         }
282         /* rows with all variables fixed should not be used */
283         if (nv == 0)
284         {  mir->skip[i] = 1;
285            continue;
286         }
287      }
288      return;
289}
290
291void *ios_mir_init(glp_tree *tree)
292{     /* initialize MIR cut generator */
293      glp_prob *mip = tree->mip;
294      int m = mip->m;
295      int n = mip->n;
296      struct MIR *mir;
297#if _MIR_DEBUG
298      xprintf("ios_mir_init: warning: debug mode enabled\n");
299#endif
300      /* allocate working area */
301      mir = xmalloc(sizeof(struct MIR));
302      mir->m = m;
303      mir->n = n;
304      mir->skip = xcalloc(1+m, sizeof(char));
305      mir->isint = xcalloc(1+m+n, sizeof(char));
306      mir->lb = xcalloc(1+m+n, sizeof(double));
307      mir->vlb = xcalloc(1+m+n, sizeof(int));
308      mir->ub = xcalloc(1+m+n, sizeof(double));
309      mir->vub = xcalloc(1+m+n, sizeof(int));
310      mir->x = xcalloc(1+m+n, sizeof(double));
311      mir->agg_row = xcalloc(1+MAXAGGR, sizeof(int));
312      mir->agg_vec = ios_create_vec(m+n);
313      mir->subst = xcalloc(1+m+n, sizeof(char));
314      mir->mod_vec = ios_create_vec(m+n);
315      mir->cut_vec = ios_create_vec(m+n);
316      /* set global row attributes */
317      set_row_attrib(tree, mir);
318      /* set global column attributes */
319      set_col_attrib(tree, mir);
320      /* set variable bounds */
321      set_var_bounds(tree, mir);
322      /* mark rows which should not be used */
323      mark_useless_rows(tree, mir);
324      return mir;
325}
326
327/***********************************************************************
328*  NAME
329*
330*  ios_mir_gen - generate MIR cuts
331*
332*  SYNOPSIS
333*
334*  #include "glpios.h"
335*  void ios_mir_gen(glp_tree *tree, void *gen, IOSPOOL *pool);
336*
337*  DESCRIPTION
338*
339*  The routine ios_mir_gen generates MIR cuts for the current point and
340*  adds them to the cut pool. */
341
342static void get_current_point(glp_tree *tree, struct MIR *mir)
343{     /* obtain current point */
344      glp_prob *mip = tree->mip;
345      int m = mir->m;
346      int n = mir->n;
347      int k;
348      for (k = 1; k <= m; k++)
349         mir->x[k] = mip->row[k]->prim;
350      for (k = m+1; k <= m+n; k++)
351         mir->x[k] = mip->col[k-m]->prim;
352      return;
353}
354
355#if _MIR_DEBUG
356static void check_current_point(struct MIR *mir)
357{     /* check current point */
358      int m = mir->m;
359      int n = mir->n;
360      int k, kk;
361      double lb, ub, eps;
362      for (k = 1; k <= m+n; k++)
363      {  /* determine lower bound */
364         lb = mir->lb[k];
365         kk = mir->vlb[k];
366         if (kk != 0)
367         {  xassert(lb != -DBL_MAX);
368            xassert(!mir->isint[k]);
369            xassert(mir->isint[kk]);
370            lb *= mir->x[kk];
371         }
372         /* check lower bound */
373         if (lb != -DBL_MAX)
374         {  eps = 1e-6 * (1.0 + fabs(lb));
375            xassert(mir->x[k] >= lb - eps);
376         }
377         /* determine upper bound */
378         ub = mir->ub[k];
379         kk = mir->vub[k];
380         if (kk != 0)
381         {  xassert(ub != +DBL_MAX);
382            xassert(!mir->isint[k]);
383            xassert(mir->isint[kk]);
384            ub *= mir->x[kk];
385         }
386         /* check upper bound */
387         if (ub != +DBL_MAX)
388         {  eps = 1e-6 * (1.0 + fabs(ub));
389            xassert(mir->x[k] <= ub + eps);
390         }
391      }
392      return;
393}
394#endif
395
396static void initial_agg_row(glp_tree *tree, struct MIR *mir, int i)
397{     /* use original i-th row as initial aggregated constraint */
398      glp_prob *mip = tree->mip;
399      int m = mir->m;
400      GLPAIJ *aij;
401      xassert(1 <= i && i <= m);
402      xassert(!mir->skip[i]);
403      /* mark i-th row in order not to use it in the same aggregated
404         constraint */
405      mir->skip[i] = 2;
406      mir->agg_cnt = 1;
407      mir->agg_row[1] = i;
408      /* use x[i] - sum a[i,j] * x[m+j] = 0, where x[i] is auxiliary
409         variable of row i, x[m+j] are structural variables */
410      ios_clear_vec(mir->agg_vec);
411      ios_set_vj(mir->agg_vec, i, 1.0);
412      for (aij = mip->row[i]->ptr; aij != NULL; aij = aij->r_next)
413         ios_set_vj(mir->agg_vec, m + aij->col->j, - aij->val);
414      mir->agg_rhs = 0.0;
415#if _MIR_DEBUG
416      ios_check_vec(mir->agg_vec);
417#endif
418      return;
419}
420
421#if _MIR_DEBUG
422static void check_agg_row(struct MIR *mir)
423{     /* check aggregated constraint */
424      int m = mir->m;
425      int n = mir->n;
426      int j, k;
427      double r, big;
428      /* compute the residual r = sum a[k] * x[k] - b and determine
429         big = max(1, |a[k]|, |b|) */
430      r = 0.0, big = 1.0;
431      for (j = 1; j <= mir->agg_vec->nnz; j++)
432      {  k = mir->agg_vec->ind[j];
433         xassert(1 <= k && k <= m+n);
434         r += mir->agg_vec->val[j] * mir->x[k];
435         if (big < fabs(mir->agg_vec->val[j]))
436            big = fabs(mir->agg_vec->val[j]);
437      }
438      r -= mir->agg_rhs;
439      if (big < fabs(mir->agg_rhs))
440         big = fabs(mir->agg_rhs);
441      /* the residual must be close to zero */
442      xassert(fabs(r) <= 1e-6 * big);
443      return;
444}
445#endif
446
447static void subst_fixed_vars(struct MIR *mir)
448{     /* substitute fixed variables into aggregated constraint */
449      int m = mir->m;
450      int n = mir->n;
451      int j, k;
452      for (j = 1; j <= mir->agg_vec->nnz; j++)
453      {  k = mir->agg_vec->ind[j];
454         xassert(1 <= k && k <= m+n);
455         if (mir->vlb[k] == 0 && mir->vub[k] == 0 &&
456             mir->lb[k] == mir->ub[k])
457         {  /* x[k] is fixed */
458            mir->agg_rhs -= mir->agg_vec->val[j] * mir->lb[k];
459            mir->agg_vec->val[j] = 0.0;
460         }
461      }
462      /* remove terms corresponding to fixed variables */
463      ios_clean_vec(mir->agg_vec, DBL_EPSILON);
464#if _MIR_DEBUG
465      ios_check_vec(mir->agg_vec);
466#endif
467      return;
468}
469
470static void bound_subst_heur(struct MIR *mir)
471{     /* bound substitution heuristic */
472      int m = mir->m;
473      int n = mir->n;
474      int j, k, kk;
475      double d1, d2;
476      for (j = 1; j <= mir->agg_vec->nnz; j++)
477      {  k = mir->agg_vec->ind[j];
478         xassert(1 <= k && k <= m+n);
479         if (mir->isint[k]) continue; /* skip integer variable */
480         /* compute distance from x[k] to its lower bound */
481         kk = mir->vlb[k];
482         if (kk == 0)
483         {  if (mir->lb[k] == -DBL_MAX)
484               d1 = DBL_MAX;
485            else
486               d1 = mir->x[k] - mir->lb[k];
487         }
488         else
489         {  xassert(1 <= kk && kk <= m+n);
490            xassert(mir->isint[kk]);
491            xassert(mir->lb[k] != -DBL_MAX);
492            d1 = mir->x[k] - mir->lb[k] * mir->x[kk];
493         }
494         /* compute distance from x[k] to its upper bound */
495         kk = mir->vub[k];
496         if (kk == 0)
497         {  if (mir->vub[k] == +DBL_MAX)
498               d2 = DBL_MAX;
499            else
500               d2 = mir->ub[k] - mir->x[k];
501         }
502         else
503         {  xassert(1 <= kk && kk <= m+n);
504            xassert(mir->isint[kk]);
505            xassert(mir->ub[k] != +DBL_MAX);
506            d2 = mir->ub[k] * mir->x[kk] - mir->x[k];
507         }
508         /* x[k] cannot be free */
509         xassert(d1 != DBL_MAX || d2 != DBL_MAX);
510         /* choose the bound which is closer to x[k] */
511         xassert(mir->subst[k] == '?');
512         if (d1 <= d2)
513            mir->subst[k] = 'L';
514         else
515            mir->subst[k] = 'U';
516      }
517      return;
518}
519
520static void build_mod_row(struct MIR *mir)
521{     /* substitute bounds and build modified constraint */
522      int m = mir->m;
523      int n = mir->n;
524      int j, jj, k, kk;
525      /* initially modified constraint is aggregated constraint */
526      ios_copy_vec(mir->mod_vec, mir->agg_vec);
527      mir->mod_rhs = mir->agg_rhs;
528#if _MIR_DEBUG
529      ios_check_vec(mir->mod_vec);
530#endif
531      /* substitute bounds for continuous variables; note that due to
532         substitution of variable bounds additional terms may appear in
533         modified constraint */
534      for (j = mir->mod_vec->nnz; j >= 1; j--)
535      {  k = mir->mod_vec->ind[j];
536         xassert(1 <= k && k <= m+n);
537         if (mir->isint[k]) continue; /* skip integer variable */
538         if (mir->subst[k] == 'L')
539         {  /* x[k] = (lower bound) + x'[k] */
540            xassert(mir->lb[k] != -DBL_MAX);
541            kk = mir->vlb[k];
542            if (kk == 0)
543            {  /* x[k] = lb[k] + x'[k] */
544               mir->mod_rhs -= mir->mod_vec->val[j] * mir->lb[k];
545            }
546            else
547            {  /* x[k] = lb[k] * x[kk] + x'[k] */
548               xassert(mir->isint[kk]);
549               jj = mir->mod_vec->pos[kk];
550               if (jj == 0)
551               {  ios_set_vj(mir->mod_vec, kk, 1.0);
552                  jj = mir->mod_vec->pos[kk];
553                  mir->mod_vec->val[jj] = 0.0;
554               }
555               mir->mod_vec->val[jj] +=
556                  mir->mod_vec->val[j] * mir->lb[k];
557            }
558         }
559         else if (mir->subst[k] == 'U')
560         {  /* x[k] = (upper bound) - x'[k] */
561            xassert(mir->ub[k] != +DBL_MAX);
562            kk = mir->vub[k];
563            if (kk == 0)
564            {  /* x[k] = ub[k] - x'[k] */
565               mir->mod_rhs -= mir->mod_vec->val[j] * mir->ub[k];
566            }
567            else
568            {  /* x[k] = ub[k] * x[kk] - x'[k] */
569               xassert(mir->isint[kk]);
570               jj = mir->mod_vec->pos[kk];
571               if (jj == 0)
572               {  ios_set_vj(mir->mod_vec, kk, 1.0);
573                  jj = mir->mod_vec->pos[kk];
574                  mir->mod_vec->val[jj] = 0.0;
575               }
576               mir->mod_vec->val[jj] +=
577                  mir->mod_vec->val[j] * mir->ub[k];
578            }
579            mir->mod_vec->val[j] = - mir->mod_vec->val[j];
580         }
581         else
582            xassert(k != k);
583      }
584#if _MIR_DEBUG
585      ios_check_vec(mir->mod_vec);
586#endif
587      /* substitute bounds for integer variables */
588      for (j = 1; j <= mir->mod_vec->nnz; j++)
589      {  k = mir->mod_vec->ind[j];
590         xassert(1 <= k && k <= m+n);
591         if (!mir->isint[k]) continue; /* skip continuous variable */
592         xassert(mir->subst[k] == '?');
593         xassert(mir->vlb[k] == 0 && mir->vub[k] == 0);
594         xassert(mir->lb[k] != -DBL_MAX && mir->ub[k] != +DBL_MAX);
595         if (fabs(mir->lb[k]) <= fabs(mir->ub[k]))
596         {  /* x[k] = lb[k] + x'[k] */
597            mir->subst[k] = 'L';
598            mir->mod_rhs -= mir->mod_vec->val[j] * mir->lb[k];
599         }
600         else
601         {  /* x[k] = ub[k] - x'[k] */
602            mir->subst[k] = 'U';
603            mir->mod_rhs -= mir->mod_vec->val[j] * mir->ub[k];
604            mir->mod_vec->val[j] = - mir->mod_vec->val[j];
605         }
606      }
607#if _MIR_DEBUG
608      ios_check_vec(mir->mod_vec);
609#endif
610      return;
611}
612
613#if _MIR_DEBUG
614static void check_mod_row(struct MIR *mir)
615{     /* check modified constraint */
616      int m = mir->m;
617      int n = mir->n;
618      int j, k, kk;
619      double r, big, x;
620      /* compute the residual r = sum a'[k] * x'[k] - b' and determine
621         big = max(1, |a[k]|, |b|) */
622      r = 0.0, big = 1.0;
623      for (j = 1; j <= mir->mod_vec->nnz; j++)
624      {  k = mir->mod_vec->ind[j];
625         xassert(1 <= k && k <= m+n);
626         if (mir->subst[k] == 'L')
627         {  /* x'[k] = x[k] - (lower bound) */
628            xassert(mir->lb[k] != -DBL_MAX);
629            kk = mir->vlb[k];
630            if (kk == 0)
631               x = mir->x[k] - mir->lb[k];
632            else
633               x = mir->x[k] - mir->lb[k] * mir->x[kk];
634         }
635         else if (mir->subst[k] == 'U')
636         {  /* x'[k] = (upper bound) - x[k] */
637            xassert(mir->ub[k] != +DBL_MAX);
638            kk = mir->vub[k];
639            if (kk == 0)
640               x = mir->ub[k] - mir->x[k];
641            else
642               x = mir->ub[k] * mir->x[kk] - mir->x[k];
643         }
644         else
645            xassert(k != k);
646         r += mir->mod_vec->val[j] * x;
647         if (big < fabs(mir->mod_vec->val[j]))
648            big = fabs(mir->mod_vec->val[j]);
649      }
650      r -= mir->mod_rhs;
651      if (big < fabs(mir->mod_rhs))
652         big = fabs(mir->mod_rhs);
653      /* the residual must be close to zero */
654      xassert(fabs(r) <= 1e-6 * big);
655      return;
656}
657#endif
658
659/***********************************************************************
660*  mir_ineq - construct MIR inequality
661*
662*  Given the single constraint mixed integer set
663*
664*                    |N|
665*     X = {(x,s) in Z    x R  : sum   a[j] * x[j] <= b + s},
666*                    +      +  j in N
667*
668*  this routine constructs the mixed integer rounding (MIR) inequality
669*
670*     sum   alpha[j] * x[j] <= beta + gamma * s,
671*    j in N
672*
673*  which is valid for X.
674*
675*  If the MIR inequality has been successfully constructed, the routine
676*  returns zero. Otherwise, if b is close to nearest integer, there may
677*  be numeric difficulties due to big coefficients; so in this case the
678*  routine returns non-zero. */
679
680static int mir_ineq(const int n, const double a[], const double b,
681      double alpha[], double *beta, double *gamma)
682{     int j;
683      double f, t;
684      if (fabs(b - floor(b + .5)) < 0.01)
685         return 1;
686      f = b - floor(b);
687      for (j = 1; j <= n; j++)
688      {  t = (a[j] - floor(a[j])) - f;
689         if (t <= 0.0)
690            alpha[j] = floor(a[j]);
691         else
692            alpha[j] = floor(a[j]) + t / (1.0 - f);
693      }
694      *beta = floor(b);
695      *gamma = 1.0 / (1.0 - f);
696      return 0;
697}
698
699/***********************************************************************
700*  cmir_ineq - construct c-MIR inequality
701*
702*  Given the mixed knapsack set
703*
704*      MK              |N|
705*     X   = {(x,s) in Z    x R  : sum   a[j] * x[j] <= b + s,
706*                      +      +  j in N
707*
708*             x[j] <= u[j]},
709*
710*  a subset C of variables to be complemented, and a divisor delta > 0,
711*  this routine constructs the complemented MIR (c-MIR) inequality
712*
713*     sum   alpha[j] * x[j] <= beta + gamma * s,
714*    j in N
715*                      MK
716*  which is valid for X  .
717*
718*  If the c-MIR inequality has been successfully constructed, the
719*  routine returns zero. Otherwise, if there is a risk of numerical
720*  difficulties due to big coefficients (see comments to the routine
721*  mir_ineq), the routine cmir_ineq returns non-zero. */
722
723static int cmir_ineq(const int n, const double a[], const double b,
724      const double u[], const char cset[], const double delta,
725      double alpha[], double *beta, double *gamma)
726{     int j;
727      double *aa, bb;
728      aa = alpha, bb = b;
729      for (j = 1; j <= n; j++)
730      {  aa[j] = a[j] / delta;
731         if (cset[j])
732            aa[j] = - aa[j], bb -= a[j] * u[j];
733      }
734      bb /= delta;
735      if (mir_ineq(n, aa, bb, alpha, beta, gamma)) return 1;
736      for (j = 1; j <= n; j++)
737      {  if (cset[j])
738            alpha[j] = - alpha[j], *beta += alpha[j] * u[j];
739      }
740      *gamma /= delta;
741      return 0;
742}
743
744/***********************************************************************
745*  cmir_sep - c-MIR separation heuristic
746*
747*  Given the mixed knapsack set
748*
749*      MK              |N|
750*     X   = {(x,s) in Z    x R  : sum   a[j] * x[j] <= b + s,
751*                      +      +  j in N
752*
753*             x[j] <= u[j]}
754*
755*                           *   *
756*  and a fractional point (x , s ), this routine tries to construct
757*  c-MIR inequality
758*
759*     sum   alpha[j] * x[j] <= beta + gamma * s,
760*    j in N
761*                      MK
762*  which is valid for X   and has (desirably maximal) violation at the
763*  fractional point given. This is attained by choosing an appropriate
764*  set C of variables to be complemented and a divisor delta > 0, which
765*  together define corresponding c-MIR inequality.
766*
767*  If a violated c-MIR inequality has been successfully constructed,
768*  the routine returns its violation:
769*
770*                       *                      *
771*     sum   alpha[j] * x [j] - beta - gamma * s ,
772*    j in N
773*
774*  which is positive. In case of failure the routine returns zero. */
775
776struct vset { int j; double v; };
777
778static int cmir_cmp(const void *p1, const void *p2)
779{     const struct vset *v1 = p1, *v2 = p2;
780      if (v1->v < v2->v) return -1;
781      if (v1->v > v2->v) return +1;
782      return 0;
783}
784
785static double cmir_sep(const int n, const double a[], const double b,
786      const double u[], const double x[], const double s,
787      double alpha[], double *beta, double *gamma)
788{     int fail, j, k, nv, v;
789      double delta, eps, d_try[1+3], r, r_best;
790      char *cset;
791      struct vset *vset;
792      /* allocate working arrays */
793      cset = xcalloc(1+n, sizeof(char));
794      vset = xcalloc(1+n, sizeof(struct vset));
795      /* choose initial C */
796      for (j = 1; j <= n; j++)
797         cset[j] = (char)(x[j] >= 0.5 * u[j]);
798      /* choose initial delta */
799      r_best = delta = 0.0;
800      for (j = 1; j <= n; j++)
801      {  xassert(a[j] != 0.0);
802         /* if x[j] is close to its bounds, skip it */
803         eps = 1e-9 * (1.0 + fabs(u[j]));
804         if (x[j] < eps || x[j] > u[j] - eps) continue;
805         /* try delta = |a[j]| to construct c-MIR inequality */
806         fail = cmir_ineq(n, a, b, u, cset, fabs(a[j]), alpha, beta,
807            gamma);
808         if (fail) continue;
809         /* compute violation */
810         r = - (*beta) - (*gamma) * s;
811         for (k = 1; k <= n; k++) r += alpha[k] * x[k];
812         if (r_best < r) r_best = r, delta = fabs(a[j]);
813      }
814      if (r_best < 0.001) r_best = 0.0;
815      if (r_best == 0.0) goto done;
816      xassert(delta > 0.0);
817      /* try to increase violation by dividing delta by 2, 4, and 8,
818         respectively */
819      d_try[1] = delta / 2.0;
820      d_try[2] = delta / 4.0;
821      d_try[3] = delta / 8.0;
822      for (j = 1; j <= 3; j++)
823      {  /* construct c-MIR inequality */
824         fail = cmir_ineq(n, a, b, u, cset, d_try[j], alpha, beta,
825            gamma);
826         if (fail) continue;
827         /* compute violation */
828         r = - (*beta) - (*gamma) * s;
829         for (k = 1; k <= n; k++) r += alpha[k] * x[k];
830         if (r_best < r) r_best = r, delta = d_try[j];
831      }
832      /* build subset of variables lying strictly between their bounds
833         and order it by nondecreasing values of |x[j] - u[j]/2| */
834      nv = 0;
835      for (j = 1; j <= n; j++)
836      {  /* if x[j] is close to its bounds, skip it */
837         eps = 1e-9 * (1.0 + fabs(u[j]));
838         if (x[j] < eps || x[j] > u[j] - eps) continue;
839         /* add x[j] to the subset */
840         nv++;
841         vset[nv].j = j;
842         vset[nv].v = fabs(x[j] - 0.5 * u[j]);
843      }
844      qsort(&vset[1], nv, sizeof(struct vset), cmir_cmp);
845      /* try to increase violation by successively complementing each
846         variable in the subset */
847      for (v = 1; v <= nv; v++)
848      {  j = vset[v].j;
849         /* replace x[j] by its complement or vice versa */
850         cset[j] = (char)!cset[j];
851         /* construct c-MIR inequality */
852         fail = cmir_ineq(n, a, b, u, cset, delta, alpha, beta, gamma);
853         /* restore the variable */
854         cset[j] = (char)!cset[j];
855         /* do not replace the variable in case of failure */
856         if (fail) continue;
857         /* compute violation */
858         r = - (*beta) - (*gamma) * s;
859         for (k = 1; k <= n; k++) r += alpha[k] * x[k];
860         if (r_best < r) r_best = r, cset[j] = (char)!cset[j];
861      }
862      /* construct the best c-MIR inequality chosen */
863      fail = cmir_ineq(n, a, b, u, cset, delta, alpha, beta, gamma);
864      xassert(!fail);
865done: /* free working arrays */
866      xfree(cset);
867      xfree(vset);
868      /* return to the calling routine */
869      return r_best;
870}
871
872static double generate(struct MIR *mir)
873{     /* try to generate violated c-MIR cut for modified constraint */
874      int m = mir->m;
875      int n = mir->n;
876      int j, k, kk, nint;
877      double s, *u, *x, *alpha, r_best = 0.0, b, beta, gamma;
878      ios_copy_vec(mir->cut_vec, mir->mod_vec);
879      mir->cut_rhs = mir->mod_rhs;
880      /* remove small terms, which can appear due to substitution of
881         variable bounds */
882      ios_clean_vec(mir->cut_vec, DBL_EPSILON);
883#if _MIR_DEBUG
884      ios_check_vec(mir->cut_vec);
885#endif
886      /* remove positive continuous terms to obtain MK relaxation */
887      for (j = 1; j <= mir->cut_vec->nnz; j++)
888      {  k = mir->cut_vec->ind[j];
889         xassert(1 <= k && k <= m+n);
890         if (!mir->isint[k] && mir->cut_vec->val[j] > 0.0)
891            mir->cut_vec->val[j] = 0.0;
892      }
893      ios_clean_vec(mir->cut_vec, 0.0);
894#if _MIR_DEBUG
895      ios_check_vec(mir->cut_vec);
896#endif
897      /* move integer terms to the beginning of the sparse vector and
898         determine the number of integer variables */
899      nint = 0;
900      for (j = 1; j <= mir->cut_vec->nnz; j++)
901      {  k = mir->cut_vec->ind[j];
902         xassert(1 <= k && k <= m+n);
903         if (mir->isint[k])
904         {  double temp;
905            nint++;
906            /* interchange elements [nint] and [j] */
907            kk = mir->cut_vec->ind[nint];
908            mir->cut_vec->pos[k] = nint;
909            mir->cut_vec->pos[kk] = j;
910            mir->cut_vec->ind[nint] = k;
911            mir->cut_vec->ind[j] = kk;
912            temp = mir->cut_vec->val[nint];
913            mir->cut_vec->val[nint] = mir->cut_vec->val[j];
914            mir->cut_vec->val[j] = temp;
915         }
916      }
917#if _MIR_DEBUG
918      ios_check_vec(mir->cut_vec);
919#endif
920      /* if there is no integer variable, nothing to generate */
921      if (nint == 0) goto done;
922      /* allocate working arrays */
923      u = xcalloc(1+nint, sizeof(double));
924      x = xcalloc(1+nint, sizeof(double));
925      alpha = xcalloc(1+nint, sizeof(double));
926      /* determine u and x */
927      for (j = 1; j <= nint; j++)
928      {  k = mir->cut_vec->ind[j];
929         xassert(m+1 <= k && k <= m+n);
930         xassert(mir->isint[k]);
931         u[j] = mir->ub[k] - mir->lb[k];
932         xassert(u[j] >= 1.0);
933         if (mir->subst[k] == 'L')
934            x[j] = mir->x[k] - mir->lb[k];
935         else if (mir->subst[k] == 'U')
936            x[j] = mir->ub[k] - mir->x[k];
937         else
938            xassert(k != k);
939         xassert(x[j] >= -0.001);
940         if (x[j] < 0.0) x[j] = 0.0;
941      }
942      /* compute s = - sum of continuous terms */
943      s = 0.0;
944      for (j = nint+1; j <= mir->cut_vec->nnz; j++)
945      {  double x;
946         k = mir->cut_vec->ind[j];
947         xassert(1 <= k && k <= m+n);
948         /* must be continuous */
949         xassert(!mir->isint[k]);
950         if (mir->subst[k] == 'L')
951         {  xassert(mir->lb[k] != -DBL_MAX);
952            kk = mir->vlb[k];
953            if (kk == 0)
954               x = mir->x[k] - mir->lb[k];
955            else
956               x = mir->x[k] - mir->lb[k] * mir->x[kk];
957         }
958         else if (mir->subst[k] == 'U')
959         {  xassert(mir->ub[k] != +DBL_MAX);
960            kk = mir->vub[k];
961            if (kk == 0)
962               x = mir->ub[k] - mir->x[k];
963            else
964               x = mir->ub[k] * mir->x[kk] - mir->x[k];
965         }
966         else
967            xassert(k != k);
968         xassert(x >= -0.001);
969         if (x < 0.0) x = 0.0;
970         s -= mir->cut_vec->val[j] * x;
971      }
972      xassert(s >= 0.0);
973      /* apply heuristic to obtain most violated c-MIR inequality */
974      b = mir->cut_rhs;
975      r_best = cmir_sep(nint, mir->cut_vec->val, b, u, x, s, alpha,
976         &beta, &gamma);
977      if (r_best == 0.0) goto skip;
978      xassert(r_best > 0.0);
979      /* convert to raw cut */
980      /* sum alpha[j] * x[j] <= beta + gamma * s */
981      for (j = 1; j <= nint; j++)
982         mir->cut_vec->val[j] = alpha[j];
983      for (j = nint+1; j <= mir->cut_vec->nnz; j++)
984      {  k = mir->cut_vec->ind[j];
985         if (k <= m+n) mir->cut_vec->val[j] *= gamma;
986      }
987      mir->cut_rhs = beta;
988#if _MIR_DEBUG
989      ios_check_vec(mir->cut_vec);
990#endif
991skip: /* free working arrays */
992      xfree(u);
993      xfree(x);
994      xfree(alpha);
995done: return r_best;
996}
997
998#if _MIR_DEBUG
999static void check_raw_cut(struct MIR *mir, double r_best)
1000{     /* check raw cut before back bound substitution */
1001      int m = mir->m;
1002      int n = mir->n;
1003      int j, k, kk;
1004      double r, big, x;
1005      /* compute the residual r = sum a[k] * x[k] - b and determine
1006         big = max(1, |a[k]|, |b|) */
1007      r = 0.0, big = 1.0;
1008      for (j = 1; j <= mir->cut_vec->nnz; j++)
1009      {  k = mir->cut_vec->ind[j];
1010         xassert(1 <= k && k <= m+n);
1011         if (mir->subst[k] == 'L')
1012         {  xassert(mir->lb[k] != -DBL_MAX);
1013            kk = mir->vlb[k];
1014            if (kk == 0)
1015               x = mir->x[k] - mir->lb[k];
1016            else
1017               x = mir->x[k] - mir->lb[k] * mir->x[kk];
1018         }
1019         else if (mir->subst[k] == 'U')
1020         {  xassert(mir->ub[k] != +DBL_MAX);
1021            kk = mir->vub[k];
1022            if (kk == 0)
1023               x = mir->ub[k] - mir->x[k];
1024            else
1025               x = mir->ub[k] * mir->x[kk] - mir->x[k];
1026         }
1027         else
1028            xassert(k != k);
1029         r += mir->cut_vec->val[j] * x;
1030         if (big < fabs(mir->cut_vec->val[j]))
1031            big = fabs(mir->cut_vec->val[j]);
1032      }
1033      r -= mir->cut_rhs;
1034      if (big < fabs(mir->cut_rhs))
1035         big = fabs(mir->cut_rhs);
1036      /* the residual must be close to r_best */
1037      xassert(fabs(r - r_best) <= 1e-6 * big);
1038      return;
1039}
1040#endif
1041
1042static void back_subst(struct MIR *mir)
1043{     /* back substitution of original bounds */
1044      int m = mir->m;
1045      int n = mir->n;
1046      int j, jj, k, kk;
1047      /* at first, restore bounds of integer variables (because on
1048         restoring variable bounds of continuous variables we need
1049         original, not shifted, bounds of integer variables) */
1050      for (j = 1; j <= mir->cut_vec->nnz; j++)
1051      {  k = mir->cut_vec->ind[j];
1052         xassert(1 <= k && k <= m+n);
1053         if (!mir->isint[k]) continue; /* skip continuous */
1054         if (mir->subst[k] == 'L')
1055         {  /* x'[k] = x[k] - lb[k] */
1056            xassert(mir->lb[k] != -DBL_MAX);
1057            xassert(mir->vlb[k] == 0);
1058            mir->cut_rhs += mir->cut_vec->val[j] * mir->lb[k];
1059         }
1060         else if (mir->subst[k] == 'U')
1061         {  /* x'[k] = ub[k] - x[k] */
1062            xassert(mir->ub[k] != +DBL_MAX);
1063            xassert(mir->vub[k] == 0);
1064            mir->cut_rhs -= mir->cut_vec->val[j] * mir->ub[k];
1065            mir->cut_vec->val[j] = - mir->cut_vec->val[j];
1066         }
1067         else
1068            xassert(k != k);
1069      }
1070      /* now restore bounds of continuous variables */
1071      for (j = 1; j <= mir->cut_vec->nnz; j++)
1072      {  k = mir->cut_vec->ind[j];
1073         xassert(1 <= k && k <= m+n);
1074         if (mir->isint[k]) continue; /* skip integer */
1075         if (mir->subst[k] == 'L')
1076         {  /* x'[k] = x[k] - (lower bound) */
1077            xassert(mir->lb[k] != -DBL_MAX);
1078            kk = mir->vlb[k];
1079            if (kk == 0)
1080            {  /* x'[k] = x[k] - lb[k] */
1081               mir->cut_rhs += mir->cut_vec->val[j] * mir->lb[k];
1082            }
1083            else
1084            {  /* x'[k] = x[k] - lb[k] * x[kk] */
1085               jj = mir->cut_vec->pos[kk];
1086#if 0
1087               xassert(jj != 0);
1088#else
1089               if (jj == 0)
1090               {  ios_set_vj(mir->cut_vec, kk, 1.0);
1091                  jj = mir->cut_vec->pos[kk];
1092                  xassert(jj != 0);
1093                  mir->cut_vec->val[jj] = 0.0;
1094               }
1095#endif
1096               mir->cut_vec->val[jj] -= mir->cut_vec->val[j] *
1097                  mir->lb[k];
1098            }
1099         }
1100         else if (mir->subst[k] == 'U')
1101         {  /* x'[k] = (upper bound) - x[k] */
1102            xassert(mir->ub[k] != +DBL_MAX);
1103            kk = mir->vub[k];
1104            if (kk == 0)
1105            {  /* x'[k] = ub[k] - x[k] */
1106               mir->cut_rhs -= mir->cut_vec->val[j] * mir->ub[k];
1107            }
1108            else
1109            {  /* x'[k] = ub[k] * x[kk] - x[k] */
1110               jj = mir->cut_vec->pos[kk];
1111               if (jj == 0)
1112               {  ios_set_vj(mir->cut_vec, kk, 1.0);
1113                  jj = mir->cut_vec->pos[kk];
1114                  xassert(jj != 0);
1115                  mir->cut_vec->val[jj] = 0.0;
1116               }
1117               mir->cut_vec->val[jj] += mir->cut_vec->val[j] *
1118                  mir->ub[k];
1119            }
1120            mir->cut_vec->val[j] = - mir->cut_vec->val[j];
1121         }
1122         else
1123            xassert(k != k);
1124      }
1125#if _MIR_DEBUG
1126      ios_check_vec(mir->cut_vec);
1127#endif
1128      return;
1129}
1130
1131#if _MIR_DEBUG
1132static void check_cut_row(struct MIR *mir, double r_best)
1133{     /* check the cut after back bound substitution or elimination of
1134         auxiliary variables */
1135      int m = mir->m;
1136      int n = mir->n;
1137      int j, k;
1138      double r, big;
1139      /* compute the residual r = sum a[k] * x[k] - b and determine
1140         big = max(1, |a[k]|, |b|) */
1141      r = 0.0, big = 1.0;
1142      for (j = 1; j <= mir->cut_vec->nnz; j++)
1143      {  k = mir->cut_vec->ind[j];
1144         xassert(1 <= k && k <= m+n);
1145         r += mir->cut_vec->val[j] * mir->x[k];
1146         if (big < fabs(mir->cut_vec->val[j]))
1147            big = fabs(mir->cut_vec->val[j]);
1148      }
1149      r -= mir->cut_rhs;
1150      if (big < fabs(mir->cut_rhs))
1151         big = fabs(mir->cut_rhs);
1152      /* the residual must be close to r_best */
1153      xassert(fabs(r - r_best) <= 1e-6 * big);
1154      return;
1155}
1156#endif
1157
1158static void subst_aux_vars(glp_tree *tree, struct MIR *mir)
1159{     /* final substitution to eliminate auxiliary variables */
1160      glp_prob *mip = tree->mip;
1161      int m = mir->m;
1162      int n = mir->n;
1163      GLPAIJ *aij;
1164      int j, k, kk, jj;
1165      for (j = mir->cut_vec->nnz; j >= 1; j--)
1166      {  k = mir->cut_vec->ind[j];
1167         xassert(1 <= k && k <= m+n);
1168         if (k > m) continue; /* skip structurals */
1169         for (aij = mip->row[k]->ptr; aij != NULL; aij = aij->r_next)
1170         {  kk = m + aij->col->j; /* structural */
1171            jj = mir->cut_vec->pos[kk];
1172            if (jj == 0)
1173            {  ios_set_vj(mir->cut_vec, kk, 1.0);
1174               jj = mir->cut_vec->pos[kk];
1175               mir->cut_vec->val[jj] = 0.0;
1176            }
1177            mir->cut_vec->val[jj] += mir->cut_vec->val[j] * aij->val;
1178         }
1179         mir->cut_vec->val[j] = 0.0;
1180      }
1181      ios_clean_vec(mir->cut_vec, 0.0);
1182      return;
1183}
1184
1185static void add_cut(glp_tree *tree, struct MIR *mir)
1186{     /* add constructed cut inequality to the cut pool */
1187      int m = mir->m;
1188      int n = mir->n;
1189      int j, k, len;
1190      int *ind = xcalloc(1+n, sizeof(int));
1191      double *val = xcalloc(1+n, sizeof(double));
1192      len = 0;
1193      for (j = mir->cut_vec->nnz; j >= 1; j--)
1194      {  k = mir->cut_vec->ind[j];
1195         xassert(m+1 <= k && k <= m+n);
1196         len++, ind[len] = k - m, val[len] = mir->cut_vec->val[j];
1197      }
1198#if 0
1199      ios_add_cut_row(tree, pool, GLP_RF_MIR, len, ind, val, GLP_UP,
1200         mir->cut_rhs);
1201#else
1202      glp_ios_add_row(tree, NULL, GLP_RF_MIR, 0, len, ind, val, GLP_UP,
1203         mir->cut_rhs);
1204#endif
1205      xfree(ind);
1206      xfree(val);
1207      return;
1208}
1209
1210static int aggregate_row(glp_tree *tree, struct MIR *mir)
1211{     /* try to aggregate another row */
1212      glp_prob *mip = tree->mip;
1213      int m = mir->m;
1214      int n = mir->n;
1215      GLPAIJ *aij;
1216      IOSVEC *v;
1217      int ii, j, jj, k, kk, kappa = 0, ret = 0;
1218      double d1, d2, d, d_max = 0.0;
1219      /* choose appropriate structural variable in the aggregated row
1220         to be substituted */
1221      for (j = 1; j <= mir->agg_vec->nnz; j++)
1222      {  k = mir->agg_vec->ind[j];
1223         xassert(1 <= k && k <= m+n);
1224         if (k <= m) continue; /* skip auxiliary var */
1225         if (mir->isint[k]) continue; /* skip integer var */
1226         if (fabs(mir->agg_vec->val[j]) < 0.001) continue;
1227         /* compute distance from x[k] to its lower bound */
1228         kk = mir->vlb[k];
1229         if (kk == 0)
1230         {  if (mir->lb[k] == -DBL_MAX)
1231               d1 = DBL_MAX;
1232            else
1233               d1 = mir->x[k] - mir->lb[k];
1234         }
1235         else
1236         {  xassert(1 <= kk && kk <= m+n);
1237            xassert(mir->isint[kk]);
1238            xassert(mir->lb[k] != -DBL_MAX);
1239            d1 = mir->x[k] - mir->lb[k] * mir->x[kk];
1240         }
1241         /* compute distance from x[k] to its upper bound */
1242         kk = mir->vub[k];
1243         if (kk == 0)
1244         {  if (mir->vub[k] == +DBL_MAX)
1245               d2 = DBL_MAX;
1246            else
1247               d2 = mir->ub[k] - mir->x[k];
1248         }
1249         else
1250         {  xassert(1 <= kk && kk <= m+n);
1251            xassert(mir->isint[kk]);
1252            xassert(mir->ub[k] != +DBL_MAX);
1253            d2 = mir->ub[k] * mir->x[kk] - mir->x[k];
1254         }
1255         /* x[k] cannot be free */
1256         xassert(d1 != DBL_MAX || d2 != DBL_MAX);
1257         /* d = min(d1, d2) */
1258         d = (d1 <= d2 ? d1 : d2);
1259         xassert(d != DBL_MAX);
1260         /* should not be close to corresponding bound */
1261         if (d < 0.001) continue;
1262         if (d_max < d) d_max = d, kappa = k;
1263      }
1264      if (kappa == 0)
1265      {  /* nothing chosen */
1266         ret = 1;
1267         goto done;
1268      }
1269      /* x[kappa] has been chosen */
1270      xassert(m+1 <= kappa && kappa <= m+n);
1271      xassert(!mir->isint[kappa]);
1272      /* find another row, which have not been used yet, to eliminate
1273         x[kappa] from the aggregated row */
1274      for (ii = 1; ii <= m; ii++)
1275      {  if (mir->skip[ii]) continue;
1276         for (aij = mip->row[ii]->ptr; aij != NULL; aij = aij->r_next)
1277            if (aij->col->j == kappa - m) break;
1278         if (aij != NULL && fabs(aij->val) >= 0.001) break;
1279      }
1280      if (ii > m)
1281      {  /* nothing found */
1282         ret = 2;
1283         goto done;
1284      }
1285      /* row ii has been found; include it in the aggregated list */
1286      mir->agg_cnt++;
1287      xassert(mir->agg_cnt <= MAXAGGR);
1288      mir->agg_row[mir->agg_cnt] = ii;
1289      mir->skip[ii] = 2;
1290      /* v := new row */
1291      v = ios_create_vec(m+n);
1292      ios_set_vj(v, ii, 1.0);
1293      for (aij = mip->row[ii]->ptr; aij != NULL; aij = aij->r_next)
1294         ios_set_vj(v, m + aij->col->j, - aij->val);
1295#if _MIR_DEBUG
1296      ios_check_vec(v);
1297#endif
1298      /* perform gaussian elimination to remove x[kappa] */
1299      j = mir->agg_vec->pos[kappa];
1300      xassert(j != 0);
1301      jj = v->pos[kappa];
1302      xassert(jj != 0);
1303      ios_linear_comb(mir->agg_vec,
1304         - mir->agg_vec->val[j] / v->val[jj], v);
1305      ios_delete_vec(v);
1306      ios_set_vj(mir->agg_vec, kappa, 0.0);
1307#if _MIR_DEBUG
1308      ios_check_vec(mir->agg_vec);
1309#endif
1310done: return ret;
1311}
1312
1313void ios_mir_gen(glp_tree *tree, void *gen)
1314{     /* main routine to generate MIR cuts */
1315      glp_prob *mip = tree->mip;
1316      struct MIR *mir = gen;
1317      int m = mir->m;
1318      int n = mir->n;
1319      int i;
1320      double r_best;
1321      xassert(mip->m >= m);
1322      xassert(mip->n == n);
1323      /* obtain current point */
1324      get_current_point(tree, mir);
1325#if _MIR_DEBUG
1326      /* check current point */
1327      check_current_point(mir);
1328#endif
1329      /* reset bound substitution flags */
1330      memset(&mir->subst[1], '?', m+n);
1331      /* try to generate a set of violated MIR cuts */
1332      for (i = 1; i <= m; i++)
1333      {  if (mir->skip[i]) continue;
1334         /* use original i-th row as initial aggregated constraint */
1335         initial_agg_row(tree, mir, i);
1336loop:    ;
1337#if _MIR_DEBUG
1338         /* check aggregated row */
1339         check_agg_row(mir);
1340#endif
1341         /* substitute fixed variables into aggregated constraint */
1342         subst_fixed_vars(mir);
1343#if _MIR_DEBUG
1344         /* check aggregated row */
1345         check_agg_row(mir);
1346#endif
1347#if _MIR_DEBUG
1348         /* check bound substitution flags */
1349         {  int k;
1350            for (k = 1; k <= m+n; k++)
1351               xassert(mir->subst[k] == '?');
1352         }
1353#endif
1354         /* apply bound substitution heuristic */
1355         bound_subst_heur(mir);
1356         /* substitute bounds and build modified constraint */
1357         build_mod_row(mir);
1358#if _MIR_DEBUG
1359         /* check modified row */
1360         check_mod_row(mir);
1361#endif
1362         /* try to generate violated c-MIR cut for modified row */
1363         r_best = generate(mir);
1364         if (r_best > 0.0)
1365         {  /* success */
1366#if _MIR_DEBUG
1367            /* check raw cut before back bound substitution */
1368            check_raw_cut(mir, r_best);
1369#endif
1370            /* back substitution of original bounds */
1371            back_subst(mir);
1372#if _MIR_DEBUG
1373            /* check the cut after back bound substitution */
1374            check_cut_row(mir, r_best);
1375#endif
1376            /* final substitution to eliminate auxiliary variables */
1377            subst_aux_vars(tree, mir);
1378#if _MIR_DEBUG
1379            /* check the cut after elimination of auxiliaries */
1380            check_cut_row(mir, r_best);
1381#endif
1382            /* add constructed cut inequality to the cut pool */
1383            add_cut(tree, mir);
1384         }
1385         /* reset bound substitution flags */
1386         {  int j, k;
1387            for (j = 1; j <= mir->mod_vec->nnz; j++)
1388            {  k = mir->mod_vec->ind[j];
1389               xassert(1 <= k && k <= m+n);
1390               xassert(mir->subst[k] != '?');
1391               mir->subst[k] = '?';
1392            }
1393         }
1394         if (r_best == 0.0)
1395         {  /* failure */
1396            if (mir->agg_cnt < MAXAGGR)
1397            {  /* try to aggregate another row */
1398               if (aggregate_row(tree, mir) == 0) goto loop;
1399            }
1400         }
1401         /* unmark rows used in the aggregated constraint */
1402         {  int k, ii;
1403            for (k = 1; k <= mir->agg_cnt; k++)
1404            {  ii = mir->agg_row[k];
1405               xassert(1 <= ii && ii <= m);
1406               xassert(mir->skip[ii] == 2);
1407               mir->skip[ii] = 0;
1408            }
1409         }
1410      }
1411      return;
1412}
1413
1414/***********************************************************************
1415*  NAME
1416*
1417*  ios_mir_term - terminate MIR cut generator
1418*
1419*  SYNOPSIS
1420*
1421*  #include "glpios.h"
1422*  void ios_mir_term(void *gen);
1423*
1424*  DESCRIPTION
1425*
1426*  The routine ios_mir_term deletes the MIR cut generator working area
1427*  freeing all the memory allocated to it. */
1428
1429void ios_mir_term(void *gen)
1430{     struct MIR *mir = gen;
1431      xfree(mir->skip);
1432      xfree(mir->isint);
1433      xfree(mir->lb);
1434      xfree(mir->vlb);
1435      xfree(mir->ub);
1436      xfree(mir->vub);
1437      xfree(mir->x);
1438      xfree(mir->agg_row);
1439      ios_delete_vec(mir->agg_vec);
1440      xfree(mir->subst);
1441      ios_delete_vec(mir->mod_vec);
1442      ios_delete_vec(mir->cut_vec);
1443      xfree(mir);
1444      return;
1445}
1446
1447/* eof */
Note: See TracBrowser for help on using the repository browser.