1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpios06.c	Sun Nov 06 20:59:10 2011 +0100
     1.3 @@ -0,0 +1,1447 @@
     1.4 +/* glpios06.c (MIR cut generator) */
     1.5 +
     1.6 +/***********************************************************************
     1.7 +*  This code is part of GLPK (GNU Linear Programming Kit).
     1.8 +*
     1.9 +*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
    1.10 +*  2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
    1.11 +*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
    1.12 +*  E-mail: <mao@gnu.org>.
    1.13 +*
    1.14 +*  GLPK is free software: you can redistribute it and/or modify it
    1.15 +*  under the terms of the GNU General Public License as published by
    1.16 +*  the Free Software Foundation, either version 3 of the License, or
    1.17 +*  (at your option) any later version.
    1.18 +*
    1.19 +*  GLPK is distributed in the hope that it will be useful, but WITHOUT
    1.20 +*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    1.21 +*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
    1.22 +*  License for more details.
    1.23 +*
    1.24 +*  You should have received a copy of the GNU General Public License
    1.25 +*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
    1.26 +***********************************************************************/
    1.27 +
    1.28 +#include "glpios.h"
    1.29 +
    1.30 +#define _MIR_DEBUG 0
    1.31 +
    1.32 +#define MAXAGGR 5
    1.33 +/* maximal number of rows which can be aggregated */
    1.34 +
    1.35 +struct MIR
    1.36 +{     /* MIR cut generator working area */
    1.37 +      /*--------------------------------------------------------------*/
    1.38 +      /* global information valid for the root subproblem */
    1.39 +      int m;
    1.40 +      /* number of rows (in the root subproblem) */
    1.41 +      int n;
    1.42 +      /* number of columns */
    1.43 +      char *skip; /* char skip[1+m]; */
    1.44 +      /* skip[i], 1 <= i <= m, is a flag that means that row i should
    1.45 +         not be used because (1) it is not suitable, or (2) because it
    1.46 +         has been used in the aggregated constraint */
    1.47 +      char *isint; /* char isint[1+m+n]; */
    1.48 +      /* isint[k], 1 <= k <= m+n, is a flag that means that variable
    1.49 +         x[k] is integer (otherwise, continuous) */
    1.50 +      double *lb; /* double lb[1+m+n]; */
    1.51 +      /* lb[k], 1 <= k <= m+n, is lower bound of x[k]; -DBL_MAX means
    1.52 +         that x[k] has no lower bound */
    1.53 +      int *vlb; /* int vlb[1+m+n]; */
    1.54 +      /* vlb[k] = k', 1 <= k <= m+n, is the number of integer variable,
    1.55 +         which defines variable lower bound x[k] >= lb[k] * x[k']; zero
    1.56 +         means that x[k] has simple lower bound */
    1.57 +      double *ub; /* double ub[1+m+n]; */
    1.58 +      /* ub[k], 1 <= k <= m+n, is upper bound of x[k]; +DBL_MAX means
    1.59 +         that x[k] has no upper bound */
    1.60 +      int *vub; /* int vub[1+m+n]; */
    1.61 +      /* vub[k] = k', 1 <= k <= m+n, is the number of integer variable,
    1.62 +         which defines variable upper bound x[k] <= ub[k] * x[k']; zero
    1.63 +         means that x[k] has simple upper bound */
    1.64 +      /*--------------------------------------------------------------*/
    1.65 +      /* current (fractional) point to be separated */
    1.66 +      double *x; /* double x[1+m+n]; */
    1.67 +      /* x[k] is current value of auxiliary (1 <= k <= m) or structural
    1.68 +         (m+1 <= k <= m+n) variable */
    1.69 +      /*--------------------------------------------------------------*/
    1.70 +      /* aggregated constraint sum a[k] * x[k] = b, which is a linear
    1.71 +         combination of original constraints transformed to equalities
    1.72 +         by introducing auxiliary variables */
    1.73 +      int agg_cnt;
    1.74 +      /* number of rows (original constraints) used to build aggregated
    1.75 +         constraint, 1 <= agg_cnt <= MAXAGGR */
    1.76 +      int *agg_row; /* int agg_row[1+MAXAGGR]; */
    1.77 +      /* agg_row[k], 1 <= k <= agg_cnt, is the row number used to build
    1.78 +         aggregated constraint */
    1.79 +      IOSVEC *agg_vec; /* IOSVEC agg_vec[1:m+n]; */
    1.80 +      /* sparse vector of aggregated constraint coefficients, a[k] */
    1.81 +      double agg_rhs;
    1.82 +      /* right-hand side of the aggregated constraint, b */
    1.83 +      /*--------------------------------------------------------------*/
    1.84 +      /* bound substitution flags for modified constraint */
    1.85 +      char *subst; /* char subst[1+m+n]; */
    1.86 +      /* subst[k], 1 <= k <= m+n, is a bound substitution flag used for
    1.87 +         variable x[k]:
    1.88 +         '?' - x[k] is missing in modified constraint
    1.89 +         'L' - x[k] = (lower bound) + x'[k]
    1.90 +         'U' - x[k] = (upper bound) - x'[k] */
    1.91 +      /*--------------------------------------------------------------*/
    1.92 +      /* modified constraint sum a'[k] * x'[k] = b', where x'[k] >= 0,
    1.93 +         derived from aggregated constraint by substituting bounds;
    1.94 +         note that due to substitution of variable bounds there may be
    1.95 +         additional terms in the modified constraint */
    1.96 +      IOSVEC *mod_vec; /* IOSVEC mod_vec[1:m+n]; */
    1.97 +      /* sparse vector of modified constraint coefficients, a'[k] */
    1.98 +      double mod_rhs;
    1.99 +      /* right-hand side of the modified constraint, b' */
   1.100 +      /*--------------------------------------------------------------*/
   1.101 +      /* cutting plane sum alpha[k] * x[k] <= beta */
   1.102 +      IOSVEC *cut_vec; /* IOSVEC cut_vec[1:m+n]; */
   1.103 +      /* sparse vector of cutting plane coefficients, alpha[k] */
   1.104 +      double cut_rhs;
   1.105 +      /* right-hand size of the cutting plane, beta */
   1.106 +};
   1.107 +
   1.108 +/***********************************************************************
   1.109 +*  NAME
   1.110 +*
   1.111 +*  ios_mir_init - initialize MIR cut generator
   1.112 +*
   1.113 +*  SYNOPSIS
   1.114 +*
   1.115 +*  #include "glpios.h"
   1.116 +*  void *ios_mir_init(glp_tree *tree);
   1.117 +*
   1.118 +*  DESCRIPTION
   1.119 +*
   1.120 +*  The routine ios_mir_init initializes the MIR cut generator assuming
   1.121 +*  that the current subproblem is the root subproblem.
   1.122 +*
   1.123 +*  RETURNS
   1.124 +*
   1.125 +*  The routine ios_mir_init returns a pointer to the MIR cut generator
   1.126 +*  working area. */
   1.127 +
   1.128 +static void set_row_attrib(glp_tree *tree, struct MIR *mir)
   1.129 +{     /* set global row attributes */
   1.130 +      glp_prob *mip = tree->mip;
   1.131 +      int m = mir->m;
   1.132 +      int k;
   1.133 +      for (k = 1; k <= m; k++)
   1.134 +      {  GLPROW *row = mip->row[k];
   1.135 +         mir->skip[k] = 0;
   1.136 +         mir->isint[k] = 0;
   1.137 +         switch (row->type)
   1.138 +         {  case GLP_FR:
   1.139 +               mir->lb[k] = -DBL_MAX, mir->ub[k] = +DBL_MAX; break;
   1.140 +            case GLP_LO:
   1.141 +               mir->lb[k] = row->lb, mir->ub[k] = +DBL_MAX; break;
   1.142 +            case GLP_UP:
   1.143 +               mir->lb[k] = -DBL_MAX, mir->ub[k] = row->ub; break;
   1.144 +            case GLP_DB:
   1.145 +               mir->lb[k] = row->lb, mir->ub[k] = row->ub; break;
   1.146 +            case GLP_FX:
   1.147 +               mir->lb[k] = mir->ub[k] = row->lb; break;
   1.148 +            default:
   1.149 +               xassert(row != row);
   1.150 +         }
   1.151 +         mir->vlb[k] = mir->vub[k] = 0;
   1.152 +      }
   1.153 +      return;
   1.154 +}
   1.155 +
   1.156 +static void set_col_attrib(glp_tree *tree, struct MIR *mir)
   1.157 +{     /* set global column attributes */
   1.158 +      glp_prob *mip = tree->mip;
   1.159 +      int m = mir->m;
   1.160 +      int n = mir->n;
   1.161 +      int k;
   1.162 +      for (k = m+1; k <= m+n; k++)
   1.163 +      {  GLPCOL *col = mip->col[k-m];
   1.164 +         switch (col->kind)
   1.165 +         {  case GLP_CV:
   1.166 +               mir->isint[k] = 0; break;
   1.167 +            case GLP_IV:
   1.168 +               mir->isint[k] = 1; break;
   1.169 +            default:
   1.170 +               xassert(col != col);
   1.171 +         }
   1.172 +         switch (col->type)
   1.173 +         {  case GLP_FR:
   1.174 +               mir->lb[k] = -DBL_MAX, mir->ub[k] = +DBL_MAX; break;
   1.175 +            case GLP_LO:
   1.176 +               mir->lb[k] = col->lb, mir->ub[k] = +DBL_MAX; break;
   1.177 +            case GLP_UP:
   1.178 +               mir->lb[k] = -DBL_MAX, mir->ub[k] = col->ub; break;
   1.179 +            case GLP_DB:
   1.180 +               mir->lb[k] = col->lb, mir->ub[k] = col->ub; break;
   1.181 +            case GLP_FX:
   1.182 +               mir->lb[k] = mir->ub[k] = col->lb; break;
   1.183 +            default:
   1.184 +               xassert(col != col);
   1.185 +         }
   1.186 +         mir->vlb[k] = mir->vub[k] = 0;
   1.187 +      }
   1.188 +      return;
   1.189 +}
   1.190 +
   1.191 +static void set_var_bounds(glp_tree *tree, struct MIR *mir)
   1.192 +{     /* set variable bounds */
   1.193 +      glp_prob *mip = tree->mip;
   1.194 +      int m = mir->m;
   1.195 +      GLPAIJ *aij;
   1.196 +      int i, k1, k2;
   1.197 +      double a1, a2;
   1.198 +      for (i = 1; i <= m; i++)
   1.199 +      {  /* we need the row to be '>= 0' or '<= 0' */
   1.200 +         if (!(mir->lb[i] == 0.0 && mir->ub[i] == +DBL_MAX ||
   1.201 +               mir->lb[i] == -DBL_MAX && mir->ub[i] == 0.0)) continue;
   1.202 +         /* take first term */
   1.203 +         aij = mip->row[i]->ptr;
   1.204 +         if (aij == NULL) continue;
   1.205 +         k1 = m + aij->col->j, a1 = aij->val;
   1.206 +         /* take second term */
   1.207 +         aij = aij->r_next;
   1.208 +         if (aij == NULL) continue;
   1.209 +         k2 = m + aij->col->j, a2 = aij->val;
   1.210 +         /* there must be only two terms */
   1.211 +         if (aij->r_next != NULL) continue;
   1.212 +         /* interchange terms, if needed */
   1.213 +         if (!mir->isint[k1] && mir->isint[k2])
   1.214 +            ;
   1.215 +         else if (mir->isint[k1] && !mir->isint[k2])
   1.216 +         {  k2 = k1, a2 = a1;
   1.217 +            k1 = m + aij->col->j, a1 = aij->val;
   1.218 +         }
   1.219 +         else
   1.220 +         {  /* both terms are either continuous or integer */
   1.221 +            continue;
   1.222 +         }
   1.223 +         /* x[k2] should be double-bounded */
   1.224 +         if (mir->lb[k2] == -DBL_MAX || mir->ub[k2] == +DBL_MAX ||
   1.225 +             mir->lb[k2] == mir->ub[k2]) continue;
   1.226 +         /* change signs, if necessary */
   1.227 +         if (mir->ub[i] == 0.0) a1 = - a1, a2 = - a2;
   1.228 +         /* now the row has the form a1 * x1 + a2 * x2 >= 0, where x1
   1.229 +            is continuous, x2 is integer */
   1.230 +         if (a1 > 0.0)
   1.231 +         {  /* x1 >= - (a2 / a1) * x2 */
   1.232 +            if (mir->vlb[k1] == 0)
   1.233 +            {  /* set variable lower bound for x1 */
   1.234 +               mir->lb[k1] = - a2 / a1;
   1.235 +               mir->vlb[k1] = k2;
   1.236 +               /* the row should not be used */
   1.237 +               mir->skip[i] = 1;
   1.238 +            }
   1.239 +         }
   1.240 +         else /* a1 < 0.0 */
   1.241 +         {  /* x1 <= - (a2 / a1) * x2 */
   1.242 +            if (mir->vub[k1] == 0)
   1.243 +            {  /* set variable upper bound for x1 */
   1.244 +               mir->ub[k1] = - a2 / a1;
   1.245 +               mir->vub[k1] = k2;
   1.246 +               /* the row should not be used */
   1.247 +               mir->skip[i] = 1;
   1.248 +            }
   1.249 +         }
   1.250 +      }
   1.251 +      return;
   1.252 +}
   1.253 +
   1.254 +static void mark_useless_rows(glp_tree *tree, struct MIR *mir)
   1.255 +{     /* mark rows which should not be used */
   1.256 +      glp_prob *mip = tree->mip;
   1.257 +      int m = mir->m;
   1.258 +      GLPAIJ *aij;
   1.259 +      int i, k, nv;
   1.260 +      for (i = 1; i <= m; i++)
   1.261 +      {  /* free rows should not be used */
   1.262 +         if (mir->lb[i] == -DBL_MAX && mir->ub[i] == +DBL_MAX)
   1.263 +         {  mir->skip[i] = 1;
   1.264 +            continue;
   1.265 +         }
   1.266 +         nv = 0;
   1.267 +         for (aij = mip->row[i]->ptr; aij != NULL; aij = aij->r_next)
   1.268 +         {  k = m + aij->col->j;
   1.269 +            /* rows with free variables should not be used */
   1.270 +            if (mir->lb[k] == -DBL_MAX && mir->ub[k] == +DBL_MAX)
   1.271 +            {  mir->skip[i] = 1;
   1.272 +               break;
   1.273 +            }
   1.274 +            /* rows with integer variables having infinite (lower or
   1.275 +               upper) bound should not be used */
   1.276 +            if (mir->isint[k] && mir->lb[k] == -DBL_MAX ||
   1.277 +                mir->isint[k] && mir->ub[k] == +DBL_MAX)
   1.278 +            {  mir->skip[i] = 1;
   1.279 +               break;
   1.280 +            }
   1.281 +            /* count non-fixed variables */
   1.282 +            if (!(mir->vlb[k] == 0 && mir->vub[k] == 0 &&
   1.283 +                  mir->lb[k] == mir->ub[k])) nv++;
   1.284 +         }
   1.285 +         /* rows with all variables fixed should not be used */
   1.286 +         if (nv == 0)
   1.287 +         {  mir->skip[i] = 1;
   1.288 +            continue;
   1.289 +         }
   1.290 +      }
   1.291 +      return;
   1.292 +}
   1.293 +
   1.294 +void *ios_mir_init(glp_tree *tree)
   1.295 +{     /* initialize MIR cut generator */
   1.296 +      glp_prob *mip = tree->mip;
   1.297 +      int m = mip->m;
   1.298 +      int n = mip->n;
   1.299 +      struct MIR *mir;
   1.300 +#if _MIR_DEBUG
   1.301 +      xprintf("ios_mir_init: warning: debug mode enabled\n");
   1.302 +#endif
   1.303 +      /* allocate working area */
   1.304 +      mir = xmalloc(sizeof(struct MIR));
   1.305 +      mir->m = m;
   1.306 +      mir->n = n;
   1.307 +      mir->skip = xcalloc(1+m, sizeof(char));
   1.308 +      mir->isint = xcalloc(1+m+n, sizeof(char));
   1.309 +      mir->lb = xcalloc(1+m+n, sizeof(double));
   1.310 +      mir->vlb = xcalloc(1+m+n, sizeof(int));
   1.311 +      mir->ub = xcalloc(1+m+n, sizeof(double));
   1.312 +      mir->vub = xcalloc(1+m+n, sizeof(int));
   1.313 +      mir->x = xcalloc(1+m+n, sizeof(double));
   1.314 +      mir->agg_row = xcalloc(1+MAXAGGR, sizeof(int));
   1.315 +      mir->agg_vec = ios_create_vec(m+n);
   1.316 +      mir->subst = xcalloc(1+m+n, sizeof(char));
   1.317 +      mir->mod_vec = ios_create_vec(m+n);
   1.318 +      mir->cut_vec = ios_create_vec(m+n);
   1.319 +      /* set global row attributes */
   1.320 +      set_row_attrib(tree, mir);
   1.321 +      /* set global column attributes */
   1.322 +      set_col_attrib(tree, mir);
   1.323 +      /* set variable bounds */
   1.324 +      set_var_bounds(tree, mir);
   1.325 +      /* mark rows which should not be used */
   1.326 +      mark_useless_rows(tree, mir);
   1.327 +      return mir;
   1.328 +}
   1.329 +
   1.330 +/***********************************************************************
   1.331 +*  NAME
   1.332 +*
   1.333 +*  ios_mir_gen - generate MIR cuts
   1.334 +*
   1.335 +*  SYNOPSIS
   1.336 +*
   1.337 +*  #include "glpios.h"
   1.338 +*  void ios_mir_gen(glp_tree *tree, void *gen, IOSPOOL *pool);
   1.339 +*
   1.340 +*  DESCRIPTION
   1.341 +*
   1.342 +*  The routine ios_mir_gen generates MIR cuts for the current point and
   1.343 +*  adds them to the cut pool. */
   1.344 +
   1.345 +static void get_current_point(glp_tree *tree, struct MIR *mir)
   1.346 +{     /* obtain current point */
   1.347 +      glp_prob *mip = tree->mip;
   1.348 +      int m = mir->m;
   1.349 +      int n = mir->n;
   1.350 +      int k;
   1.351 +      for (k = 1; k <= m; k++)
   1.352 +         mir->x[k] = mip->row[k]->prim;
   1.353 +      for (k = m+1; k <= m+n; k++)
   1.354 +         mir->x[k] = mip->col[k-m]->prim;
   1.355 +      return;
   1.356 +}
   1.357 +
   1.358 +#if _MIR_DEBUG
   1.359 +static void check_current_point(struct MIR *mir)
   1.360 +{     /* check current point */
   1.361 +      int m = mir->m;
   1.362 +      int n = mir->n;
   1.363 +      int k, kk;
   1.364 +      double lb, ub, eps;
   1.365 +      for (k = 1; k <= m+n; k++)
   1.366 +      {  /* determine lower bound */
   1.367 +         lb = mir->lb[k];
   1.368 +         kk = mir->vlb[k];
   1.369 +         if (kk != 0)
   1.370 +         {  xassert(lb != -DBL_MAX);
   1.371 +            xassert(!mir->isint[k]);
   1.372 +            xassert(mir->isint[kk]);
   1.373 +            lb *= mir->x[kk];
   1.374 +         }
   1.375 +         /* check lower bound */
   1.376 +         if (lb != -DBL_MAX)
   1.377 +         {  eps = 1e-6 * (1.0 + fabs(lb));
   1.378 +            xassert(mir->x[k] >= lb - eps);
   1.379 +         }
   1.380 +         /* determine upper bound */
   1.381 +         ub = mir->ub[k];
   1.382 +         kk = mir->vub[k];
   1.383 +         if (kk != 0)
   1.384 +         {  xassert(ub != +DBL_MAX);
   1.385 +            xassert(!mir->isint[k]);
   1.386 +            xassert(mir->isint[kk]);
   1.387 +            ub *= mir->x[kk];
   1.388 +         }
   1.389 +         /* check upper bound */
   1.390 +         if (ub != +DBL_MAX)
   1.391 +         {  eps = 1e-6 * (1.0 + fabs(ub));
   1.392 +            xassert(mir->x[k] <= ub + eps);
   1.393 +         }
   1.394 +      }
   1.395 +      return;
   1.396 +}
   1.397 +#endif
   1.398 +
   1.399 +static void initial_agg_row(glp_tree *tree, struct MIR *mir, int i)
   1.400 +{     /* use original i-th row as initial aggregated constraint */
   1.401 +      glp_prob *mip = tree->mip;
   1.402 +      int m = mir->m;
   1.403 +      GLPAIJ *aij;
   1.404 +      xassert(1 <= i && i <= m);
   1.405 +      xassert(!mir->skip[i]);
   1.406 +      /* mark i-th row in order not to use it in the same aggregated
   1.407 +         constraint */
   1.408 +      mir->skip[i] = 2;
   1.409 +      mir->agg_cnt = 1;
   1.410 +      mir->agg_row[1] = i;
   1.411 +      /* use x[i] - sum a[i,j] * x[m+j] = 0, where x[i] is auxiliary
   1.412 +         variable of row i, x[m+j] are structural variables */
   1.413 +      ios_clear_vec(mir->agg_vec);
   1.414 +      ios_set_vj(mir->agg_vec, i, 1.0);
   1.415 +      for (aij = mip->row[i]->ptr; aij != NULL; aij = aij->r_next)
   1.416 +         ios_set_vj(mir->agg_vec, m + aij->col->j, - aij->val);
   1.417 +      mir->agg_rhs = 0.0;
   1.418 +#if _MIR_DEBUG
   1.419 +      ios_check_vec(mir->agg_vec);
   1.420 +#endif
   1.421 +      return;
   1.422 +}
   1.423 +
   1.424 +#if _MIR_DEBUG
   1.425 +static void check_agg_row(struct MIR *mir)
   1.426 +{     /* check aggregated constraint */
   1.427 +      int m = mir->m;
   1.428 +      int n = mir->n;
   1.429 +      int j, k;
   1.430 +      double r, big;
   1.431 +      /* compute the residual r = sum a[k] * x[k] - b and determine
   1.432 +         big = max(1, |a[k]|, |b|) */
   1.433 +      r = 0.0, big = 1.0;
   1.434 +      for (j = 1; j <= mir->agg_vec->nnz; j++)
   1.435 +      {  k = mir->agg_vec->ind[j];
   1.436 +         xassert(1 <= k && k <= m+n);
   1.437 +         r += mir->agg_vec->val[j] * mir->x[k];
   1.438 +         if (big < fabs(mir->agg_vec->val[j]))
   1.439 +            big = fabs(mir->agg_vec->val[j]);
   1.440 +      }
   1.441 +      r -= mir->agg_rhs;
   1.442 +      if (big < fabs(mir->agg_rhs))
   1.443 +         big = fabs(mir->agg_rhs);
   1.444 +      /* the residual must be close to zero */
   1.445 +      xassert(fabs(r) <= 1e-6 * big);
   1.446 +      return;
   1.447 +}
   1.448 +#endif
   1.449 +
   1.450 +static void subst_fixed_vars(struct MIR *mir)
   1.451 +{     /* substitute fixed variables into aggregated constraint */
   1.452 +      int m = mir->m;
   1.453 +      int n = mir->n;
   1.454 +      int j, k;
   1.455 +      for (j = 1; j <= mir->agg_vec->nnz; j++)
   1.456 +      {  k = mir->agg_vec->ind[j];
   1.457 +         xassert(1 <= k && k <= m+n);
   1.458 +         if (mir->vlb[k] == 0 && mir->vub[k] == 0 &&
   1.459 +             mir->lb[k] == mir->ub[k])
   1.460 +         {  /* x[k] is fixed */
   1.461 +            mir->agg_rhs -= mir->agg_vec->val[j] * mir->lb[k];
   1.462 +            mir->agg_vec->val[j] = 0.0;
   1.463 +         }
   1.464 +      }
   1.465 +      /* remove terms corresponding to fixed variables */
   1.466 +      ios_clean_vec(mir->agg_vec, DBL_EPSILON);
   1.467 +#if _MIR_DEBUG
   1.468 +      ios_check_vec(mir->agg_vec);
   1.469 +#endif
   1.470 +      return;
   1.471 +}
   1.472 +
   1.473 +static void bound_subst_heur(struct MIR *mir)
   1.474 +{     /* bound substitution heuristic */
   1.475 +      int m = mir->m;
   1.476 +      int n = mir->n;
   1.477 +      int j, k, kk;
   1.478 +      double d1, d2;
   1.479 +      for (j = 1; j <= mir->agg_vec->nnz; j++)
   1.480 +      {  k = mir->agg_vec->ind[j];
   1.481 +         xassert(1 <= k && k <= m+n);
   1.482 +         if (mir->isint[k]) continue; /* skip integer variable */
   1.483 +         /* compute distance from x[k] to its lower bound */
   1.484 +         kk = mir->vlb[k];
   1.485 +         if (kk == 0)
   1.486 +         {  if (mir->lb[k] == -DBL_MAX)
   1.487 +               d1 = DBL_MAX;
   1.488 +            else
   1.489 +               d1 = mir->x[k] - mir->lb[k];
   1.490 +         }
   1.491 +         else
   1.492 +         {  xassert(1 <= kk && kk <= m+n);
   1.493 +            xassert(mir->isint[kk]);
   1.494 +            xassert(mir->lb[k] != -DBL_MAX);
   1.495 +            d1 = mir->x[k] - mir->lb[k] * mir->x[kk];
   1.496 +         }
   1.497 +         /* compute distance from x[k] to its upper bound */
   1.498 +         kk = mir->vub[k];
   1.499 +         if (kk == 0)
   1.500 +         {  if (mir->vub[k] == +DBL_MAX)
   1.501 +               d2 = DBL_MAX;
   1.502 +            else
   1.503 +               d2 = mir->ub[k] - mir->x[k];
   1.504 +         }
   1.505 +         else
   1.506 +         {  xassert(1 <= kk && kk <= m+n);
   1.507 +            xassert(mir->isint[kk]);
   1.508 +            xassert(mir->ub[k] != +DBL_MAX);
   1.509 +            d2 = mir->ub[k] * mir->x[kk] - mir->x[k];
   1.510 +         }
   1.511 +         /* x[k] cannot be free */
   1.512 +         xassert(d1 != DBL_MAX || d2 != DBL_MAX);
   1.513 +         /* choose the bound which is closer to x[k] */
   1.514 +         xassert(mir->subst[k] == '?');
   1.515 +         if (d1 <= d2)
   1.516 +            mir->subst[k] = 'L';
   1.517 +         else
   1.518 +            mir->subst[k] = 'U';
   1.519 +      }
   1.520 +      return;
   1.521 +}
   1.522 +
   1.523 +static void build_mod_row(struct MIR *mir)
   1.524 +{     /* substitute bounds and build modified constraint */
   1.525 +      int m = mir->m;
   1.526 +      int n = mir->n;
   1.527 +      int j, jj, k, kk;
   1.528 +      /* initially modified constraint is aggregated constraint */
   1.529 +      ios_copy_vec(mir->mod_vec, mir->agg_vec);
   1.530 +      mir->mod_rhs = mir->agg_rhs;
   1.531 +#if _MIR_DEBUG
   1.532 +      ios_check_vec(mir->mod_vec);
   1.533 +#endif
   1.534 +      /* substitute bounds for continuous variables; note that due to
   1.535 +         substitution of variable bounds additional terms may appear in
   1.536 +         modified constraint */
   1.537 +      for (j = mir->mod_vec->nnz; j >= 1; j--)
   1.538 +      {  k = mir->mod_vec->ind[j];
   1.539 +         xassert(1 <= k && k <= m+n);
   1.540 +         if (mir->isint[k]) continue; /* skip integer variable */
   1.541 +         if (mir->subst[k] == 'L')
   1.542 +         {  /* x[k] = (lower bound) + x'[k] */
   1.543 +            xassert(mir->lb[k] != -DBL_MAX);
   1.544 +            kk = mir->vlb[k];
   1.545 +            if (kk == 0)
   1.546 +            {  /* x[k] = lb[k] + x'[k] */
   1.547 +               mir->mod_rhs -= mir->mod_vec->val[j] * mir->lb[k];
   1.548 +            }
   1.549 +            else
   1.550 +            {  /* x[k] = lb[k] * x[kk] + x'[k] */
   1.551 +               xassert(mir->isint[kk]);
   1.552 +               jj = mir->mod_vec->pos[kk];
   1.553 +               if (jj == 0)
   1.554 +               {  ios_set_vj(mir->mod_vec, kk, 1.0);
   1.555 +                  jj = mir->mod_vec->pos[kk];
   1.556 +                  mir->mod_vec->val[jj] = 0.0;
   1.557 +               }
   1.558 +               mir->mod_vec->val[jj] +=
   1.559 +                  mir->mod_vec->val[j] * mir->lb[k];
   1.560 +            }
   1.561 +         }
   1.562 +         else if (mir->subst[k] == 'U')
   1.563 +         {  /* x[k] = (upper bound) - x'[k] */
   1.564 +            xassert(mir->ub[k] != +DBL_MAX);
   1.565 +            kk = mir->vub[k];
   1.566 +            if (kk == 0)
   1.567 +            {  /* x[k] = ub[k] - x'[k] */
   1.568 +               mir->mod_rhs -= mir->mod_vec->val[j] * mir->ub[k];
   1.569 +            }
   1.570 +            else
   1.571 +            {  /* x[k] = ub[k] * x[kk] - x'[k] */
   1.572 +               xassert(mir->isint[kk]);
   1.573 +               jj = mir->mod_vec->pos[kk];
   1.574 +               if (jj == 0)
   1.575 +               {  ios_set_vj(mir->mod_vec, kk, 1.0);
   1.576 +                  jj = mir->mod_vec->pos[kk];
   1.577 +                  mir->mod_vec->val[jj] = 0.0;
   1.578 +               }
   1.579 +               mir->mod_vec->val[jj] +=
   1.580 +                  mir->mod_vec->val[j] * mir->ub[k];
   1.581 +            }
   1.582 +            mir->mod_vec->val[j] = - mir->mod_vec->val[j];
   1.583 +         }
   1.584 +         else
   1.585 +            xassert(k != k);
   1.586 +      }
   1.587 +#if _MIR_DEBUG
   1.588 +      ios_check_vec(mir->mod_vec);
   1.589 +#endif
   1.590 +      /* substitute bounds for integer variables */
   1.591 +      for (j = 1; j <= mir->mod_vec->nnz; j++)
   1.592 +      {  k = mir->mod_vec->ind[j];
   1.593 +         xassert(1 <= k && k <= m+n);
   1.594 +         if (!mir->isint[k]) continue; /* skip continuous variable */
   1.595 +         xassert(mir->subst[k] == '?');
   1.596 +         xassert(mir->vlb[k] == 0 && mir->vub[k] == 0);
   1.597 +         xassert(mir->lb[k] != -DBL_MAX && mir->ub[k] != +DBL_MAX);
   1.598 +         if (fabs(mir->lb[k]) <= fabs(mir->ub[k]))
   1.599 +         {  /* x[k] = lb[k] + x'[k] */
   1.600 +            mir->subst[k] = 'L';
   1.601 +            mir->mod_rhs -= mir->mod_vec->val[j] * mir->lb[k];
   1.602 +         }
   1.603 +         else
   1.604 +         {  /* x[k] = ub[k] - x'[k] */
   1.605 +            mir->subst[k] = 'U';
   1.606 +            mir->mod_rhs -= mir->mod_vec->val[j] * mir->ub[k];
   1.607 +            mir->mod_vec->val[j] = - mir->mod_vec->val[j];
   1.608 +         }
   1.609 +      }
   1.610 +#if _MIR_DEBUG
   1.611 +      ios_check_vec(mir->mod_vec);
   1.612 +#endif
   1.613 +      return;
   1.614 +}
   1.615 +
   1.616 +#if _MIR_DEBUG
   1.617 +static void check_mod_row(struct MIR *mir)
   1.618 +{     /* check modified constraint */
   1.619 +      int m = mir->m;
   1.620 +      int n = mir->n;
   1.621 +      int j, k, kk;
   1.622 +      double r, big, x;
   1.623 +      /* compute the residual r = sum a'[k] * x'[k] - b' and determine
   1.624 +         big = max(1, |a[k]|, |b|) */
   1.625 +      r = 0.0, big = 1.0;
   1.626 +      for (j = 1; j <= mir->mod_vec->nnz; j++)
   1.627 +      {  k = mir->mod_vec->ind[j];
   1.628 +         xassert(1 <= k && k <= m+n);
   1.629 +         if (mir->subst[k] == 'L')
   1.630 +         {  /* x'[k] = x[k] - (lower bound) */
   1.631 +            xassert(mir->lb[k] != -DBL_MAX);
   1.632 +            kk = mir->vlb[k];
   1.633 +            if (kk == 0)
   1.634 +               x = mir->x[k] - mir->lb[k];
   1.635 +            else
   1.636 +               x = mir->x[k] - mir->lb[k] * mir->x[kk];
   1.637 +         }
   1.638 +         else if (mir->subst[k] == 'U')
   1.639 +         {  /* x'[k] = (upper bound) - x[k] */
   1.640 +            xassert(mir->ub[k] != +DBL_MAX);
   1.641 +            kk = mir->vub[k];
   1.642 +            if (kk == 0)
   1.643 +               x = mir->ub[k] - mir->x[k];
   1.644 +            else
   1.645 +               x = mir->ub[k] * mir->x[kk] - mir->x[k];
   1.646 +         }
   1.647 +         else
   1.648 +            xassert(k != k);
   1.649 +         r += mir->mod_vec->val[j] * x;
   1.650 +         if (big < fabs(mir->mod_vec->val[j]))
   1.651 +            big = fabs(mir->mod_vec->val[j]);
   1.652 +      }
   1.653 +      r -= mir->mod_rhs;
   1.654 +      if (big < fabs(mir->mod_rhs))
   1.655 +         big = fabs(mir->mod_rhs);
   1.656 +      /* the residual must be close to zero */
   1.657 +      xassert(fabs(r) <= 1e-6 * big);
   1.658 +      return;
   1.659 +}
   1.660 +#endif
   1.661 +
   1.662 +/***********************************************************************
   1.663 +*  mir_ineq - construct MIR inequality
   1.664 +*
   1.665 +*  Given the single constraint mixed integer set
   1.666 +*
   1.667 +*                    |N|
   1.668 +*     X = {(x,s) in Z    x R  : sum   a[j] * x[j] <= b + s},
   1.669 +*                    +      +  j in N
   1.670 +*
   1.671 +*  this routine constructs the mixed integer rounding (MIR) inequality
   1.672 +*
   1.673 +*     sum   alpha[j] * x[j] <= beta + gamma * s,
   1.674 +*    j in N
   1.675 +*
   1.676 +*  which is valid for X.
   1.677 +*
   1.678 +*  If the MIR inequality has been successfully constructed, the routine
   1.679 +*  returns zero. Otherwise, if b is close to nearest integer, there may
   1.680 +*  be numeric difficulties due to big coefficients; so in this case the
   1.681 +*  routine returns non-zero. */
   1.682 +
   1.683 +static int mir_ineq(const int n, const double a[], const double b,
   1.684 +      double alpha[], double *beta, double *gamma)
   1.685 +{     int j;
   1.686 +      double f, t;
   1.687 +      if (fabs(b - floor(b + .5)) < 0.01)
   1.688 +         return 1;
   1.689 +      f = b - floor(b);
   1.690 +      for (j = 1; j <= n; j++)
   1.691 +      {  t = (a[j] - floor(a[j])) - f;
   1.692 +         if (t <= 0.0)
   1.693 +            alpha[j] = floor(a[j]);
   1.694 +         else
   1.695 +            alpha[j] = floor(a[j]) + t / (1.0 - f);
   1.696 +      }
   1.697 +      *beta = floor(b);
   1.698 +      *gamma = 1.0 / (1.0 - f);
   1.699 +      return 0;
   1.700 +}
   1.701 +
   1.702 +/***********************************************************************
   1.703 +*  cmir_ineq - construct c-MIR inequality
   1.704 +*
   1.705 +*  Given the mixed knapsack set
   1.706 +*
   1.707 +*      MK              |N|
   1.708 +*     X   = {(x,s) in Z    x R  : sum   a[j] * x[j] <= b + s,
   1.709 +*                      +      +  j in N
   1.710 +*
   1.711 +*             x[j] <= u[j]},
   1.712 +*
   1.713 +*  a subset C of variables to be complemented, and a divisor delta > 0,
   1.714 +*  this routine constructs the complemented MIR (c-MIR) inequality
   1.715 +*
   1.716 +*     sum   alpha[j] * x[j] <= beta + gamma * s,
   1.717 +*    j in N
   1.718 +*                      MK
   1.719 +*  which is valid for X  .
   1.720 +*
   1.721 +*  If the c-MIR inequality has been successfully constructed, the
   1.722 +*  routine returns zero. Otherwise, if there is a risk of numerical
   1.723 +*  difficulties due to big coefficients (see comments to the routine
   1.724 +*  mir_ineq), the routine cmir_ineq returns non-zero. */
   1.725 +
   1.726 +static int cmir_ineq(const int n, const double a[], const double b,
   1.727 +      const double u[], const char cset[], const double delta,
   1.728 +      double alpha[], double *beta, double *gamma)
   1.729 +{     int j;
   1.730 +      double *aa, bb;
   1.731 +      aa = alpha, bb = b;
   1.732 +      for (j = 1; j <= n; j++)
   1.733 +      {  aa[j] = a[j] / delta;
   1.734 +         if (cset[j])
   1.735 +            aa[j] = - aa[j], bb -= a[j] * u[j];
   1.736 +      }
   1.737 +      bb /= delta;
   1.738 +      if (mir_ineq(n, aa, bb, alpha, beta, gamma)) return 1;
   1.739 +      for (j = 1; j <= n; j++)
   1.740 +      {  if (cset[j])
   1.741 +            alpha[j] = - alpha[j], *beta += alpha[j] * u[j];
   1.742 +      }
   1.743 +      *gamma /= delta;
   1.744 +      return 0;
   1.745 +}
   1.746 +
   1.747 +/***********************************************************************
   1.748 +*  cmir_sep - c-MIR separation heuristic
   1.749 +*
   1.750 +*  Given the mixed knapsack set
   1.751 +*
   1.752 +*      MK              |N|
   1.753 +*     X   = {(x,s) in Z    x R  : sum   a[j] * x[j] <= b + s,
   1.754 +*                      +      +  j in N
   1.755 +*
   1.756 +*             x[j] <= u[j]}
   1.757 +*
   1.758 +*                           *   *
   1.759 +*  and a fractional point (x , s ), this routine tries to construct
   1.760 +*  c-MIR inequality
   1.761 +*
   1.762 +*     sum   alpha[j] * x[j] <= beta + gamma * s,
   1.763 +*    j in N
   1.764 +*                      MK
   1.765 +*  which is valid for X   and has (desirably maximal) violation at the
   1.766 +*  fractional point given. This is attained by choosing an appropriate
   1.767 +*  set C of variables to be complemented and a divisor delta > 0, which
   1.768 +*  together define corresponding c-MIR inequality.
   1.769 +*
   1.770 +*  If a violated c-MIR inequality has been successfully constructed,
   1.771 +*  the routine returns its violation:
   1.772 +*
   1.773 +*                       *                      *
   1.774 +*     sum   alpha[j] * x [j] - beta - gamma * s ,
   1.775 +*    j in N
   1.776 +*
   1.777 +*  which is positive. In case of failure the routine returns zero. */
   1.778 +
   1.779 +struct vset { int j; double v; };
   1.780 +
   1.781 +static int cmir_cmp(const void *p1, const void *p2)
   1.782 +{     const struct vset *v1 = p1, *v2 = p2;
   1.783 +      if (v1->v < v2->v) return -1;
   1.784 +      if (v1->v > v2->v) return +1;
   1.785 +      return 0;
   1.786 +}
   1.787 +
   1.788 +static double cmir_sep(const int n, const double a[], const double b,
   1.789 +      const double u[], const double x[], const double s,
   1.790 +      double alpha[], double *beta, double *gamma)
   1.791 +{     int fail, j, k, nv, v;
   1.792 +      double delta, eps, d_try[1+3], r, r_best;
   1.793 +      char *cset;
   1.794 +      struct vset *vset;
   1.795 +      /* allocate working arrays */
   1.796 +      cset = xcalloc(1+n, sizeof(char));
   1.797 +      vset = xcalloc(1+n, sizeof(struct vset));
   1.798 +      /* choose initial C */
   1.799 +      for (j = 1; j <= n; j++)
   1.800 +         cset[j] = (char)(x[j] >= 0.5 * u[j]);
   1.801 +      /* choose initial delta */
   1.802 +      r_best = delta = 0.0;
   1.803 +      for (j = 1; j <= n; j++)
   1.804 +      {  xassert(a[j] != 0.0);
   1.805 +         /* if x[j] is close to its bounds, skip it */
   1.806 +         eps = 1e-9 * (1.0 + fabs(u[j]));
   1.807 +         if (x[j] < eps || x[j] > u[j] - eps) continue;
   1.808 +         /* try delta = |a[j]| to construct c-MIR inequality */
   1.809 +         fail = cmir_ineq(n, a, b, u, cset, fabs(a[j]), alpha, beta,
   1.810 +            gamma);
   1.811 +         if (fail) continue;
   1.812 +         /* compute violation */
   1.813 +         r = - (*beta) - (*gamma) * s;
   1.814 +         for (k = 1; k <= n; k++) r += alpha[k] * x[k];
   1.815 +         if (r_best < r) r_best = r, delta = fabs(a[j]);
   1.816 +      }
   1.817 +      if (r_best < 0.001) r_best = 0.0;
   1.818 +      if (r_best == 0.0) goto done;
   1.819 +      xassert(delta > 0.0);
   1.820 +      /* try to increase violation by dividing delta by 2, 4, and 8,
   1.821 +         respectively */
   1.822 +      d_try[1] = delta / 2.0;
   1.823 +      d_try[2] = delta / 4.0;
   1.824 +      d_try[3] = delta / 8.0;
   1.825 +      for (j = 1; j <= 3; j++)
   1.826 +      {  /* construct c-MIR inequality */
   1.827 +         fail = cmir_ineq(n, a, b, u, cset, d_try[j], alpha, beta,
   1.828 +            gamma);
   1.829 +         if (fail) continue;
   1.830 +         /* compute violation */
   1.831 +         r = - (*beta) - (*gamma) * s;
   1.832 +         for (k = 1; k <= n; k++) r += alpha[k] * x[k];
   1.833 +         if (r_best < r) r_best = r, delta = d_try[j];
   1.834 +      }
   1.835 +      /* build subset of variables lying strictly between their bounds
   1.836 +         and order it by nondecreasing values of |x[j] - u[j]/2| */
   1.837 +      nv = 0;
   1.838 +      for (j = 1; j <= n; j++)
   1.839 +      {  /* if x[j] is close to its bounds, skip it */
   1.840 +         eps = 1e-9 * (1.0 + fabs(u[j]));
   1.841 +         if (x[j] < eps || x[j] > u[j] - eps) continue;
   1.842 +         /* add x[j] to the subset */
   1.843 +         nv++;
   1.844 +         vset[nv].j = j;
   1.845 +         vset[nv].v = fabs(x[j] - 0.5 * u[j]);
   1.846 +      }
   1.847 +      qsort(&vset[1], nv, sizeof(struct vset), cmir_cmp);
   1.848 +      /* try to increase violation by successively complementing each
   1.849 +         variable in the subset */
   1.850 +      for (v = 1; v <= nv; v++)
   1.851 +      {  j = vset[v].j;
   1.852 +         /* replace x[j] by its complement or vice versa */
   1.853 +         cset[j] = (char)!cset[j];
   1.854 +         /* construct c-MIR inequality */
   1.855 +         fail = cmir_ineq(n, a, b, u, cset, delta, alpha, beta, gamma);
   1.856 +         /* restore the variable */
   1.857 +         cset[j] = (char)!cset[j];
   1.858 +         /* do not replace the variable in case of failure */
   1.859 +         if (fail) continue;
   1.860 +         /* compute violation */
   1.861 +         r = - (*beta) - (*gamma) * s;
   1.862 +         for (k = 1; k <= n; k++) r += alpha[k] * x[k];
   1.863 +         if (r_best < r) r_best = r, cset[j] = (char)!cset[j];
   1.864 +      }
   1.865 +      /* construct the best c-MIR inequality chosen */
   1.866 +      fail = cmir_ineq(n, a, b, u, cset, delta, alpha, beta, gamma);
   1.867 +      xassert(!fail);
   1.868 +done: /* free working arrays */
   1.869 +      xfree(cset);
   1.870 +      xfree(vset);
   1.871 +      /* return to the calling routine */
   1.872 +      return r_best;
   1.873 +}
   1.874 +
   1.875 +static double generate(struct MIR *mir)
   1.876 +{     /* try to generate violated c-MIR cut for modified constraint */
   1.877 +      int m = mir->m;
   1.878 +      int n = mir->n;
   1.879 +      int j, k, kk, nint;
   1.880 +      double s, *u, *x, *alpha, r_best = 0.0, b, beta, gamma;
   1.881 +      ios_copy_vec(mir->cut_vec, mir->mod_vec);
   1.882 +      mir->cut_rhs = mir->mod_rhs;
   1.883 +      /* remove small terms, which can appear due to substitution of
   1.884 +         variable bounds */
   1.885 +      ios_clean_vec(mir->cut_vec, DBL_EPSILON);
   1.886 +#if _MIR_DEBUG
   1.887 +      ios_check_vec(mir->cut_vec);
   1.888 +#endif
   1.889 +      /* remove positive continuous terms to obtain MK relaxation */
   1.890 +      for (j = 1; j <= mir->cut_vec->nnz; j++)
   1.891 +      {  k = mir->cut_vec->ind[j];
   1.892 +         xassert(1 <= k && k <= m+n);
   1.893 +         if (!mir->isint[k] && mir->cut_vec->val[j] > 0.0)
   1.894 +            mir->cut_vec->val[j] = 0.0;
   1.895 +      }
   1.896 +      ios_clean_vec(mir->cut_vec, 0.0);
   1.897 +#if _MIR_DEBUG
   1.898 +      ios_check_vec(mir->cut_vec);
   1.899 +#endif
   1.900 +      /* move integer terms to the beginning of the sparse vector and
   1.901 +         determine the number of integer variables */
   1.902 +      nint = 0;
   1.903 +      for (j = 1; j <= mir->cut_vec->nnz; j++)
   1.904 +      {  k = mir->cut_vec->ind[j];
   1.905 +         xassert(1 <= k && k <= m+n);
   1.906 +         if (mir->isint[k])
   1.907 +         {  double temp;
   1.908 +            nint++;
   1.909 +            /* interchange elements [nint] and [j] */
   1.910 +            kk = mir->cut_vec->ind[nint];
   1.911 +            mir->cut_vec->pos[k] = nint;
   1.912 +            mir->cut_vec->pos[kk] = j;
   1.913 +            mir->cut_vec->ind[nint] = k;
   1.914 +            mir->cut_vec->ind[j] = kk;
   1.915 +            temp = mir->cut_vec->val[nint];
   1.916 +            mir->cut_vec->val[nint] = mir->cut_vec->val[j];
   1.917 +            mir->cut_vec->val[j] = temp;
   1.918 +         }
   1.919 +      }
   1.920 +#if _MIR_DEBUG
   1.921 +      ios_check_vec(mir->cut_vec);
   1.922 +#endif
   1.923 +      /* if there is no integer variable, nothing to generate */
   1.924 +      if (nint == 0) goto done;
   1.925 +      /* allocate working arrays */
   1.926 +      u = xcalloc(1+nint, sizeof(double));
   1.927 +      x = xcalloc(1+nint, sizeof(double));
   1.928 +      alpha = xcalloc(1+nint, sizeof(double));
   1.929 +      /* determine u and x */
   1.930 +      for (j = 1; j <= nint; j++)
   1.931 +      {  k = mir->cut_vec->ind[j];
   1.932 +         xassert(m+1 <= k && k <= m+n);
   1.933 +         xassert(mir->isint[k]);
   1.934 +         u[j] = mir->ub[k] - mir->lb[k];
   1.935 +         xassert(u[j] >= 1.0);
   1.936 +         if (mir->subst[k] == 'L')
   1.937 +            x[j] = mir->x[k] - mir->lb[k];
   1.938 +         else if (mir->subst[k] == 'U')
   1.939 +            x[j] = mir->ub[k] - mir->x[k];
   1.940 +         else
   1.941 +            xassert(k != k);
   1.942 +         xassert(x[j] >= -0.001);
   1.943 +         if (x[j] < 0.0) x[j] = 0.0;
   1.944 +      }
   1.945 +      /* compute s = - sum of continuous terms */
   1.946 +      s = 0.0;
   1.947 +      for (j = nint+1; j <= mir->cut_vec->nnz; j++)
   1.948 +      {  double x;
   1.949 +         k = mir->cut_vec->ind[j];
   1.950 +         xassert(1 <= k && k <= m+n);
   1.951 +         /* must be continuous */
   1.952 +         xassert(!mir->isint[k]);
   1.953 +         if (mir->subst[k] == 'L')
   1.954 +         {  xassert(mir->lb[k] != -DBL_MAX);
   1.955 +            kk = mir->vlb[k];
   1.956 +            if (kk == 0)
   1.957 +               x = mir->x[k] - mir->lb[k];
   1.958 +            else
   1.959 +               x = mir->x[k] - mir->lb[k] * mir->x[kk];
   1.960 +         }
   1.961 +         else if (mir->subst[k] == 'U')
   1.962 +         {  xassert(mir->ub[k] != +DBL_MAX);
   1.963 +            kk = mir->vub[k];
   1.964 +            if (kk == 0)
   1.965 +               x = mir->ub[k] - mir->x[k];
   1.966 +            else
   1.967 +               x = mir->ub[k] * mir->x[kk] - mir->x[k];
   1.968 +         }
   1.969 +         else
   1.970 +            xassert(k != k);
   1.971 +         xassert(x >= -0.001);
   1.972 +         if (x < 0.0) x = 0.0;
   1.973 +         s -= mir->cut_vec->val[j] * x;
   1.974 +      }
   1.975 +      xassert(s >= 0.0);
   1.976 +      /* apply heuristic to obtain most violated c-MIR inequality */
   1.977 +      b = mir->cut_rhs;
   1.978 +      r_best = cmir_sep(nint, mir->cut_vec->val, b, u, x, s, alpha,
   1.979 +         &beta, &gamma);
   1.980 +      if (r_best == 0.0) goto skip;
   1.981 +      xassert(r_best > 0.0);
   1.982 +      /* convert to raw cut */
   1.983 +      /* sum alpha[j] * x[j] <= beta + gamma * s */
   1.984 +      for (j = 1; j <= nint; j++)
   1.985 +         mir->cut_vec->val[j] = alpha[j];
   1.986 +      for (j = nint+1; j <= mir->cut_vec->nnz; j++)
   1.987 +      {  k = mir->cut_vec->ind[j];
   1.988 +         if (k <= m+n) mir->cut_vec->val[j] *= gamma;
   1.989 +      }
   1.990 +      mir->cut_rhs = beta;
   1.991 +#if _MIR_DEBUG
   1.992 +      ios_check_vec(mir->cut_vec);
   1.993 +#endif
   1.994 +skip: /* free working arrays */
   1.995 +      xfree(u);
   1.996 +      xfree(x);
   1.997 +      xfree(alpha);
   1.998 +done: return r_best;
   1.999 +}
  1.1000 +
  1.1001 +#if _MIR_DEBUG
  1.1002 +static void check_raw_cut(struct MIR *mir, double r_best)
  1.1003 +{     /* check raw cut before back bound substitution */
  1.1004 +      int m = mir->m;
  1.1005 +      int n = mir->n;
  1.1006 +      int j, k, kk;
  1.1007 +      double r, big, x;
  1.1008 +      /* compute the residual r = sum a[k] * x[k] - b and determine
  1.1009 +         big = max(1, |a[k]|, |b|) */
  1.1010 +      r = 0.0, big = 1.0;
  1.1011 +      for (j = 1; j <= mir->cut_vec->nnz; j++)
  1.1012 +      {  k = mir->cut_vec->ind[j];
  1.1013 +         xassert(1 <= k && k <= m+n);
  1.1014 +         if (mir->subst[k] == 'L')
  1.1015 +         {  xassert(mir->lb[k] != -DBL_MAX);
  1.1016 +            kk = mir->vlb[k];
  1.1017 +            if (kk == 0)
  1.1018 +               x = mir->x[k] - mir->lb[k];
  1.1019 +            else
  1.1020 +               x = mir->x[k] - mir->lb[k] * mir->x[kk];
  1.1021 +         }
  1.1022 +         else if (mir->subst[k] == 'U')
  1.1023 +         {  xassert(mir->ub[k] != +DBL_MAX);
  1.1024 +            kk = mir->vub[k];
  1.1025 +            if (kk == 0)
  1.1026 +               x = mir->ub[k] - mir->x[k];
  1.1027 +            else
  1.1028 +               x = mir->ub[k] * mir->x[kk] - mir->x[k];
  1.1029 +         }
  1.1030 +         else
  1.1031 +            xassert(k != k);
  1.1032 +         r += mir->cut_vec->val[j] * x;
  1.1033 +         if (big < fabs(mir->cut_vec->val[j]))
  1.1034 +            big = fabs(mir->cut_vec->val[j]);
  1.1035 +      }
  1.1036 +      r -= mir->cut_rhs;
  1.1037 +      if (big < fabs(mir->cut_rhs))
  1.1038 +         big = fabs(mir->cut_rhs);
  1.1039 +      /* the residual must be close to r_best */
  1.1040 +      xassert(fabs(r - r_best) <= 1e-6 * big);
  1.1041 +      return;
  1.1042 +}
  1.1043 +#endif
  1.1044 +
  1.1045 +static void back_subst(struct MIR *mir)
  1.1046 +{     /* back substitution of original bounds */
  1.1047 +      int m = mir->m;
  1.1048 +      int n = mir->n;
  1.1049 +      int j, jj, k, kk;
  1.1050 +      /* at first, restore bounds of integer variables (because on
  1.1051 +         restoring variable bounds of continuous variables we need
  1.1052 +         original, not shifted, bounds of integer variables) */
  1.1053 +      for (j = 1; j <= mir->cut_vec->nnz; j++)
  1.1054 +      {  k = mir->cut_vec->ind[j];
  1.1055 +         xassert(1 <= k && k <= m+n);
  1.1056 +         if (!mir->isint[k]) continue; /* skip continuous */
  1.1057 +         if (mir->subst[k] == 'L')
  1.1058 +         {  /* x'[k] = x[k] - lb[k] */
  1.1059 +            xassert(mir->lb[k] != -DBL_MAX);
  1.1060 +            xassert(mir->vlb[k] == 0);
  1.1061 +            mir->cut_rhs += mir->cut_vec->val[j] * mir->lb[k];
  1.1062 +         }
  1.1063 +         else if (mir->subst[k] == 'U')
  1.1064 +         {  /* x'[k] = ub[k] - x[k] */
  1.1065 +            xassert(mir->ub[k] != +DBL_MAX);
  1.1066 +            xassert(mir->vub[k] == 0);
  1.1067 +            mir->cut_rhs -= mir->cut_vec->val[j] * mir->ub[k];
  1.1068 +            mir->cut_vec->val[j] = - mir->cut_vec->val[j];
  1.1069 +         }
  1.1070 +         else
  1.1071 +            xassert(k != k);
  1.1072 +      }
  1.1073 +      /* now restore bounds of continuous variables */
  1.1074 +      for (j = 1; j <= mir->cut_vec->nnz; j++)
  1.1075 +      {  k = mir->cut_vec->ind[j];
  1.1076 +         xassert(1 <= k && k <= m+n);
  1.1077 +         if (mir->isint[k]) continue; /* skip integer */
  1.1078 +         if (mir->subst[k] == 'L')
  1.1079 +         {  /* x'[k] = x[k] - (lower bound) */
  1.1080 +            xassert(mir->lb[k] != -DBL_MAX);
  1.1081 +            kk = mir->vlb[k];
  1.1082 +            if (kk == 0)
  1.1083 +            {  /* x'[k] = x[k] - lb[k] */
  1.1084 +               mir->cut_rhs += mir->cut_vec->val[j] * mir->lb[k];
  1.1085 +            }
  1.1086 +            else
  1.1087 +            {  /* x'[k] = x[k] - lb[k] * x[kk] */
  1.1088 +               jj = mir->cut_vec->pos[kk];
  1.1089 +#if 0
  1.1090 +               xassert(jj != 0);
  1.1091 +#else
  1.1092 +               if (jj == 0)
  1.1093 +               {  ios_set_vj(mir->cut_vec, kk, 1.0);
  1.1094 +                  jj = mir->cut_vec->pos[kk];
  1.1095 +                  xassert(jj != 0);
  1.1096 +                  mir->cut_vec->val[jj] = 0.0;
  1.1097 +               }
  1.1098 +#endif
  1.1099 +               mir->cut_vec->val[jj] -= mir->cut_vec->val[j] *
  1.1100 +                  mir->lb[k];
  1.1101 +            }
  1.1102 +         }
  1.1103 +         else if (mir->subst[k] == 'U')
  1.1104 +         {  /* x'[k] = (upper bound) - x[k] */
  1.1105 +            xassert(mir->ub[k] != +DBL_MAX);
  1.1106 +            kk = mir->vub[k];
  1.1107 +            if (kk == 0)
  1.1108 +            {  /* x'[k] = ub[k] - x[k] */
  1.1109 +               mir->cut_rhs -= mir->cut_vec->val[j] * mir->ub[k];
  1.1110 +            }
  1.1111 +            else
  1.1112 +            {  /* x'[k] = ub[k] * x[kk] - x[k] */
  1.1113 +               jj = mir->cut_vec->pos[kk];
  1.1114 +               if (jj == 0)
  1.1115 +               {  ios_set_vj(mir->cut_vec, kk, 1.0);
  1.1116 +                  jj = mir->cut_vec->pos[kk];
  1.1117 +                  xassert(jj != 0);
  1.1118 +                  mir->cut_vec->val[jj] = 0.0;
  1.1119 +               }
  1.1120 +               mir->cut_vec->val[jj] += mir->cut_vec->val[j] *
  1.1121 +                  mir->ub[k];
  1.1122 +            }
  1.1123 +            mir->cut_vec->val[j] = - mir->cut_vec->val[j];
  1.1124 +         }
  1.1125 +         else
  1.1126 +            xassert(k != k);
  1.1127 +      }
  1.1128 +#if _MIR_DEBUG
  1.1129 +      ios_check_vec(mir->cut_vec);
  1.1130 +#endif
  1.1131 +      return;
  1.1132 +}
  1.1133 +
  1.1134 +#if _MIR_DEBUG
  1.1135 +static void check_cut_row(struct MIR *mir, double r_best)
  1.1136 +{     /* check the cut after back bound substitution or elimination of
  1.1137 +         auxiliary variables */
  1.1138 +      int m = mir->m;
  1.1139 +      int n = mir->n;
  1.1140 +      int j, k;
  1.1141 +      double r, big;
  1.1142 +      /* compute the residual r = sum a[k] * x[k] - b and determine
  1.1143 +         big = max(1, |a[k]|, |b|) */
  1.1144 +      r = 0.0, big = 1.0;
  1.1145 +      for (j = 1; j <= mir->cut_vec->nnz; j++)
  1.1146 +      {  k = mir->cut_vec->ind[j];
  1.1147 +         xassert(1 <= k && k <= m+n);
  1.1148 +         r += mir->cut_vec->val[j] * mir->x[k];
  1.1149 +         if (big < fabs(mir->cut_vec->val[j]))
  1.1150 +            big = fabs(mir->cut_vec->val[j]);
  1.1151 +      }
  1.1152 +      r -= mir->cut_rhs;
  1.1153 +      if (big < fabs(mir->cut_rhs))
  1.1154 +         big = fabs(mir->cut_rhs);
  1.1155 +      /* the residual must be close to r_best */
  1.1156 +      xassert(fabs(r - r_best) <= 1e-6 * big);
  1.1157 +      return;
  1.1158 +}
  1.1159 +#endif
  1.1160 +
  1.1161 +static void subst_aux_vars(glp_tree *tree, struct MIR *mir)
  1.1162 +{     /* final substitution to eliminate auxiliary variables */
  1.1163 +      glp_prob *mip = tree->mip;
  1.1164 +      int m = mir->m;
  1.1165 +      int n = mir->n;
  1.1166 +      GLPAIJ *aij;
  1.1167 +      int j, k, kk, jj;
  1.1168 +      for (j = mir->cut_vec->nnz; j >= 1; j--)
  1.1169 +      {  k = mir->cut_vec->ind[j];
  1.1170 +         xassert(1 <= k && k <= m+n);
  1.1171 +         if (k > m) continue; /* skip structurals */
  1.1172 +         for (aij = mip->row[k]->ptr; aij != NULL; aij = aij->r_next)
  1.1173 +         {  kk = m + aij->col->j; /* structural */
  1.1174 +            jj = mir->cut_vec->pos[kk];
  1.1175 +            if (jj == 0)
  1.1176 +            {  ios_set_vj(mir->cut_vec, kk, 1.0);
  1.1177 +               jj = mir->cut_vec->pos[kk];
  1.1178 +               mir->cut_vec->val[jj] = 0.0;
  1.1179 +            }
  1.1180 +            mir->cut_vec->val[jj] += mir->cut_vec->val[j] * aij->val;
  1.1181 +         }
  1.1182 +         mir->cut_vec->val[j] = 0.0;
  1.1183 +      }
  1.1184 +      ios_clean_vec(mir->cut_vec, 0.0);
  1.1185 +      return;
  1.1186 +}
  1.1187 +
  1.1188 +static void add_cut(glp_tree *tree, struct MIR *mir)
  1.1189 +{     /* add constructed cut inequality to the cut pool */
  1.1190 +      int m = mir->m;
  1.1191 +      int n = mir->n;
  1.1192 +      int j, k, len;
  1.1193 +      int *ind = xcalloc(1+n, sizeof(int));
  1.1194 +      double *val = xcalloc(1+n, sizeof(double));
  1.1195 +      len = 0;
  1.1196 +      for (j = mir->cut_vec->nnz; j >= 1; j--)
  1.1197 +      {  k = mir->cut_vec->ind[j];
  1.1198 +         xassert(m+1 <= k && k <= m+n);
  1.1199 +         len++, ind[len] = k - m, val[len] = mir->cut_vec->val[j];
  1.1200 +      }
  1.1201 +#if 0
  1.1202 +      ios_add_cut_row(tree, pool, GLP_RF_MIR, len, ind, val, GLP_UP,
  1.1203 +         mir->cut_rhs);
  1.1204 +#else
  1.1205 +      glp_ios_add_row(tree, NULL, GLP_RF_MIR, 0, len, ind, val, GLP_UP,
  1.1206 +         mir->cut_rhs);
  1.1207 +#endif
  1.1208 +      xfree(ind);
  1.1209 +      xfree(val);
  1.1210 +      return;
  1.1211 +}
  1.1212 +
  1.1213 +static int aggregate_row(glp_tree *tree, struct MIR *mir)
  1.1214 +{     /* try to aggregate another row */
  1.1215 +      glp_prob *mip = tree->mip;
  1.1216 +      int m = mir->m;
  1.1217 +      int n = mir->n;
  1.1218 +      GLPAIJ *aij;
  1.1219 +      IOSVEC *v;
  1.1220 +      int ii, j, jj, k, kk, kappa = 0, ret = 0;
  1.1221 +      double d1, d2, d, d_max = 0.0;
  1.1222 +      /* choose appropriate structural variable in the aggregated row
  1.1223 +         to be substituted */
  1.1224 +      for (j = 1; j <= mir->agg_vec->nnz; j++)
  1.1225 +      {  k = mir->agg_vec->ind[j];
  1.1226 +         xassert(1 <= k && k <= m+n);
  1.1227 +         if (k <= m) continue; /* skip auxiliary var */
  1.1228 +         if (mir->isint[k]) continue; /* skip integer var */
  1.1229 +         if (fabs(mir->agg_vec->val[j]) < 0.001) continue;
  1.1230 +         /* compute distance from x[k] to its lower bound */
  1.1231 +         kk = mir->vlb[k];
  1.1232 +         if (kk == 0)
  1.1233 +         {  if (mir->lb[k] == -DBL_MAX)
  1.1234 +               d1 = DBL_MAX;
  1.1235 +            else
  1.1236 +               d1 = mir->x[k] - mir->lb[k];
  1.1237 +         }
  1.1238 +         else
  1.1239 +         {  xassert(1 <= kk && kk <= m+n);
  1.1240 +            xassert(mir->isint[kk]);
  1.1241 +            xassert(mir->lb[k] != -DBL_MAX);
  1.1242 +            d1 = mir->x[k] - mir->lb[k] * mir->x[kk];
  1.1243 +         }
  1.1244 +         /* compute distance from x[k] to its upper bound */
  1.1245 +         kk = mir->vub[k];
  1.1246 +         if (kk == 0)
  1.1247 +         {  if (mir->vub[k] == +DBL_MAX)
  1.1248 +               d2 = DBL_MAX;
  1.1249 +            else
  1.1250 +               d2 = mir->ub[k] - mir->x[k];
  1.1251 +         }
  1.1252 +         else
  1.1253 +         {  xassert(1 <= kk && kk <= m+n);
  1.1254 +            xassert(mir->isint[kk]);
  1.1255 +            xassert(mir->ub[k] != +DBL_MAX);
  1.1256 +            d2 = mir->ub[k] * mir->x[kk] - mir->x[k];
  1.1257 +         }
  1.1258 +         /* x[k] cannot be free */
  1.1259 +         xassert(d1 != DBL_MAX || d2 != DBL_MAX);
  1.1260 +         /* d = min(d1, d2) */
  1.1261 +         d = (d1 <= d2 ? d1 : d2);
  1.1262 +         xassert(d != DBL_MAX);
  1.1263 +         /* should not be close to corresponding bound */
  1.1264 +         if (d < 0.001) continue;
  1.1265 +         if (d_max < d) d_max = d, kappa = k;
  1.1266 +      }
  1.1267 +      if (kappa == 0)
  1.1268 +      {  /* nothing chosen */
  1.1269 +         ret = 1;
  1.1270 +         goto done;
  1.1271 +      }
  1.1272 +      /* x[kappa] has been chosen */
  1.1273 +      xassert(m+1 <= kappa && kappa <= m+n);
  1.1274 +      xassert(!mir->isint[kappa]);
  1.1275 +      /* find another row, which have not been used yet, to eliminate
  1.1276 +         x[kappa] from the aggregated row */
  1.1277 +      for (ii = 1; ii <= m; ii++)
  1.1278 +      {  if (mir->skip[ii]) continue;
  1.1279 +         for (aij = mip->row[ii]->ptr; aij != NULL; aij = aij->r_next)
  1.1280 +            if (aij->col->j == kappa - m) break;
  1.1281 +         if (aij != NULL && fabs(aij->val) >= 0.001) break;
  1.1282 +      }
  1.1283 +      if (ii > m)
  1.1284 +      {  /* nothing found */
  1.1285 +         ret = 2;
  1.1286 +         goto done;
  1.1287 +      }
  1.1288 +      /* row ii has been found; include it in the aggregated list */
  1.1289 +      mir->agg_cnt++;
  1.1290 +      xassert(mir->agg_cnt <= MAXAGGR);
  1.1291 +      mir->agg_row[mir->agg_cnt] = ii;
  1.1292 +      mir->skip[ii] = 2;
  1.1293 +      /* v := new row */
  1.1294 +      v = ios_create_vec(m+n);
  1.1295 +      ios_set_vj(v, ii, 1.0);
  1.1296 +      for (aij = mip->row[ii]->ptr; aij != NULL; aij = aij->r_next)
  1.1297 +         ios_set_vj(v, m + aij->col->j, - aij->val);
  1.1298 +#if _MIR_DEBUG
  1.1299 +      ios_check_vec(v);
  1.1300 +#endif
  1.1301 +      /* perform gaussian elimination to remove x[kappa] */
  1.1302 +      j = mir->agg_vec->pos[kappa];
  1.1303 +      xassert(j != 0);
  1.1304 +      jj = v->pos[kappa];
  1.1305 +      xassert(jj != 0);
  1.1306 +      ios_linear_comb(mir->agg_vec,
  1.1307 +         - mir->agg_vec->val[j] / v->val[jj], v);
  1.1308 +      ios_delete_vec(v);
  1.1309 +      ios_set_vj(mir->agg_vec, kappa, 0.0);
  1.1310 +#if _MIR_DEBUG
  1.1311 +      ios_check_vec(mir->agg_vec);
  1.1312 +#endif
  1.1313 +done: return ret;
  1.1314 +}
  1.1315 +
  1.1316 +void ios_mir_gen(glp_tree *tree, void *gen)
  1.1317 +{     /* main routine to generate MIR cuts */
  1.1318 +      glp_prob *mip = tree->mip;
  1.1319 +      struct MIR *mir = gen;
  1.1320 +      int m = mir->m;
  1.1321 +      int n = mir->n;
  1.1322 +      int i;
  1.1323 +      double r_best;
  1.1324 +      xassert(mip->m >= m);
  1.1325 +      xassert(mip->n == n);
  1.1326 +      /* obtain current point */
  1.1327 +      get_current_point(tree, mir);
  1.1328 +#if _MIR_DEBUG
  1.1329 +      /* check current point */
  1.1330 +      check_current_point(mir);
  1.1331 +#endif
  1.1332 +      /* reset bound substitution flags */
  1.1333 +      memset(&mir->subst[1], '?', m+n);
  1.1334 +      /* try to generate a set of violated MIR cuts */
  1.1335 +      for (i = 1; i <= m; i++)
  1.1336 +      {  if (mir->skip[i]) continue;
  1.1337 +         /* use original i-th row as initial aggregated constraint */
  1.1338 +         initial_agg_row(tree, mir, i);
  1.1339 +loop:    ;
  1.1340 +#if _MIR_DEBUG
  1.1341 +         /* check aggregated row */
  1.1342 +         check_agg_row(mir);
  1.1343 +#endif
  1.1344 +         /* substitute fixed variables into aggregated constraint */
  1.1345 +         subst_fixed_vars(mir);
  1.1346 +#if _MIR_DEBUG
  1.1347 +         /* check aggregated row */
  1.1348 +         check_agg_row(mir);
  1.1349 +#endif
  1.1350 +#if _MIR_DEBUG
  1.1351 +         /* check bound substitution flags */
  1.1352 +         {  int k;
  1.1353 +            for (k = 1; k <= m+n; k++)
  1.1354 +               xassert(mir->subst[k] == '?');
  1.1355 +         }
  1.1356 +#endif
  1.1357 +         /* apply bound substitution heuristic */
  1.1358 +         bound_subst_heur(mir);
  1.1359 +         /* substitute bounds and build modified constraint */
  1.1360 +         build_mod_row(mir);
  1.1361 +#if _MIR_DEBUG
  1.1362 +         /* check modified row */
  1.1363 +         check_mod_row(mir);
  1.1364 +#endif
  1.1365 +         /* try to generate violated c-MIR cut for modified row */
  1.1366 +         r_best = generate(mir);
  1.1367 +         if (r_best > 0.0)
  1.1368 +         {  /* success */
  1.1369 +#if _MIR_DEBUG
  1.1370 +            /* check raw cut before back bound substitution */
  1.1371 +            check_raw_cut(mir, r_best);
  1.1372 +#endif
  1.1373 +            /* back substitution of original bounds */
  1.1374 +            back_subst(mir);
  1.1375 +#if _MIR_DEBUG
  1.1376 +            /* check the cut after back bound substitution */
  1.1377 +            check_cut_row(mir, r_best);
  1.1378 +#endif
  1.1379 +            /* final substitution to eliminate auxiliary variables */
  1.1380 +            subst_aux_vars(tree, mir);
  1.1381 +#if _MIR_DEBUG
  1.1382 +            /* check the cut after elimination of auxiliaries */
  1.1383 +            check_cut_row(mir, r_best);
  1.1384 +#endif
  1.1385 +            /* add constructed cut inequality to the cut pool */
  1.1386 +            add_cut(tree, mir);
  1.1387 +         }
  1.1388 +         /* reset bound substitution flags */
  1.1389 +         {  int j, k;
  1.1390 +            for (j = 1; j <= mir->mod_vec->nnz; j++)
  1.1391 +            {  k = mir->mod_vec->ind[j];
  1.1392 +               xassert(1 <= k && k <= m+n);
  1.1393 +               xassert(mir->subst[k] != '?');
  1.1394 +               mir->subst[k] = '?';
  1.1395 +            }
  1.1396 +         }
  1.1397 +         if (r_best == 0.0)
  1.1398 +         {  /* failure */
  1.1399 +            if (mir->agg_cnt < MAXAGGR)
  1.1400 +            {  /* try to aggregate another row */
  1.1401 +               if (aggregate_row(tree, mir) == 0) goto loop;
  1.1402 +            }
  1.1403 +         }
  1.1404 +         /* unmark rows used in the aggregated constraint */
  1.1405 +         {  int k, ii;
  1.1406 +            for (k = 1; k <= mir->agg_cnt; k++)
  1.1407 +            {  ii = mir->agg_row[k];
  1.1408 +               xassert(1 <= ii && ii <= m);
  1.1409 +               xassert(mir->skip[ii] == 2);
  1.1410 +               mir->skip[ii] = 0;
  1.1411 +            }
  1.1412 +         }
  1.1413 +      }
  1.1414 +      return;
  1.1415 +}
  1.1416 +
  1.1417 +/***********************************************************************
  1.1418 +*  NAME
  1.1419 +*
  1.1420 +*  ios_mir_term - terminate MIR cut generator
  1.1421 +*
  1.1422 +*  SYNOPSIS
  1.1423 +*
  1.1424 +*  #include "glpios.h"
  1.1425 +*  void ios_mir_term(void *gen);
  1.1426 +*
  1.1427 +*  DESCRIPTION
  1.1428 +*
  1.1429 +*  The routine ios_mir_term deletes the MIR cut generator working area
  1.1430 +*  freeing all the memory allocated to it. */
  1.1431 +
  1.1432 +void ios_mir_term(void *gen)
  1.1433 +{     struct MIR *mir = gen;
  1.1434 +      xfree(mir->skip);
  1.1435 +      xfree(mir->isint);
  1.1436 +      xfree(mir->lb);
  1.1437 +      xfree(mir->vlb);
  1.1438 +      xfree(mir->ub);
  1.1439 +      xfree(mir->vub);
  1.1440 +      xfree(mir->x);
  1.1441 +      xfree(mir->agg_row);
  1.1442 +      ios_delete_vec(mir->agg_vec);
  1.1443 +      xfree(mir->subst);
  1.1444 +      ios_delete_vec(mir->mod_vec);
  1.1445 +      ios_delete_vec(mir->cut_vec);
  1.1446 +      xfree(mir);
  1.1447 +      return;
  1.1448 +}
  1.1449 +
  1.1450 +/* eof */