/* glpios06.c (MIR cut generator) */

 /***********************************************************************

 *  This code is part of GLPK (GNU Linear Programming Kit).

 *

 *  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,

 *  2009, 2010 Andrew Makhorin, Department for Applied Informatics,

 *  Moscow Aviation Institute, Moscow, Russia. All rights reserved.

 *  E-mail: <mao@gnu.org>.

 *

 *  GLPK is free software: you can redistribute it and/or modify it

 *  under the terms of the GNU General Public License as published by

 *  the Free Software Foundation, either version 3 of the License, or

 *  (at your option) any later version.

 *

 *  GLPK is distributed in the hope that it will be useful, but WITHOUT

 *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY

 *  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public

 *  License for more details.

 *

 *  You should have received a copy of the GNU General Public License

 *  along with GLPK. If not, see <http://www.gnu.org/licenses/>.

 ***********************************************************************/

 #include "glpios.h"

 #define _MIR_DEBUG 0

 #define MAXAGGR 5

 /* maximal number of rows which can be aggregated */

 struct MIR

 {     /* MIR cut generator working area */

       /*--------------------------------------------------------------*/

       /* global information valid for the root subproblem */

       int m;

       /* number of rows (in the root subproblem) */

       int n;

       /* number of columns */

       char *skip; /* char skip[1+m]; */

       /* skip[i], 1 <= i <= m, is a flag that means that row i should

          not be used because (1) it is not suitable, or (2) because it

          has been used in the aggregated constraint */

       char *isint; /* char isint[1+m+n]; */

       /* isint[k], 1 <= k <= m+n, is a flag that means that variable

          x[k] is integer (otherwise, continuous) */

       double *lb; /* double lb[1+m+n]; */

       /* lb[k], 1 <= k <= m+n, is lower bound of x[k]; -DBL_MAX means

          that x[k] has no lower bound */

       int *vlb; /* int vlb[1+m+n]; */

       /* vlb[k] = k', 1 <= k <= m+n, is the number of integer variable,

          which defines variable lower bound x[k] >= lb[k] * x[k']; zero

          means that x[k] has simple lower bound */

       double *ub; /* double ub[1+m+n]; */

       /* ub[k], 1 <= k <= m+n, is upper bound of x[k]; +DBL_MAX means

          that x[k] has no upper bound */

       int *vub; /* int vub[1+m+n]; */

       /* vub[k] = k', 1 <= k <= m+n, is the number of integer variable,

          which defines variable upper bound x[k] <= ub[k] * x[k']; zero

          means that x[k] has simple upper bound */

       /*--------------------------------------------------------------*/

       /* current (fractional) point to be separated */

       double *x; /* double x[1+m+n]; */

       /* x[k] is current value of auxiliary (1 <= k <= m) or structural

          (m+1 <= k <= m+n) variable */

       /*--------------------------------------------------------------*/

       /* aggregated constraint sum a[k] * x[k] = b, which is a linear

          combination of original constraints transformed to equalities

          by introducing auxiliary variables */

       int agg_cnt;

       /* number of rows (original constraints) used to build aggregated

          constraint, 1 <= agg_cnt <= MAXAGGR */

       int *agg_row; /* int agg_row[1+MAXAGGR]; */

       /* agg_row[k], 1 <= k <= agg_cnt, is the row number used to build

          aggregated constraint */

       IOSVEC *agg_vec; /* IOSVEC agg_vec[1:m+n]; */

       /* sparse vector of aggregated constraint coefficients, a[k] */

       double agg_rhs;

       /* right-hand side of the aggregated constraint, b */

       /*--------------------------------------------------------------*/

       /* bound substitution flags for modified constraint */

       char *subst; /* char subst[1+m+n]; */

       /* subst[k], 1 <= k <= m+n, is a bound substitution flag used for

          variable x[k]:

          '?' - x[k] is missing in modified constraint

          'L' - x[k] = (lower bound) + x'[k]

          'U' - x[k] = (upper bound) - x'[k] */

       /*--------------------------------------------------------------*/

       /* modified constraint sum a'[k] * x'[k] = b', where x'[k] >= 0,

          derived from aggregated constraint by substituting bounds;

          note that due to substitution of variable bounds there may be

          additional terms in the modified constraint */

       IOSVEC *mod_vec; /* IOSVEC mod_vec[1:m+n]; */

       /* sparse vector of modified constraint coefficients, a'[k] */

       double mod_rhs;

       /* right-hand side of the modified constraint, b' */

       /*--------------------------------------------------------------*/

       /* cutting plane sum alpha[k] * x[k] <= beta */

       IOSVEC *cut_vec; /* IOSVEC cut_vec[1:m+n]; */

       /* sparse vector of cutting plane coefficients, alpha[k] */

       double cut_rhs;

       /* right-hand size of the cutting plane, beta */

 };

 /***********************************************************************

 *  NAME

 *

 *  ios_mir_init - initialize MIR cut generator

 *

 *  SYNOPSIS

 *

 *  #include "glpios.h"

 *  void *ios_mir_init(glp_tree *tree);

 *

 *  DESCRIPTION

 *

 *  The routine ios_mir_init initializes the MIR cut generator assuming

 *  that the current subproblem is the root subproblem.

 *

 *  RETURNS

 *

 *  The routine ios_mir_init returns a pointer to the MIR cut generator

 *  working area. */

 static void set_row_attrib(glp_tree *tree, struct MIR *mir)

 {     /* set global row attributes */

       glp_prob *mip = tree->mip;

       int m = mir->m;

       int k;

       for (k = 1; k <= m; k++)

       {  GLPROW *row = mip->row[k];

          mir->skip[k] = 0;

          mir->isint[k] = 0;

          switch (row->type)

          {  case GLP_FR:

                mir->lb[k] = -DBL_MAX, mir->ub[k] = +DBL_MAX; break;

             case GLP_LO:

                mir->lb[k] = row->lb, mir->ub[k] = +DBL_MAX; break;

             case GLP_UP:

                mir->lb[k] = -DBL_MAX, mir->ub[k] = row->ub; break;

             case GLP_DB:

                mir->lb[k] = row->lb, mir->ub[k] = row->ub; break;

             case GLP_FX:

                mir->lb[k] = mir->ub[k] = row->lb; break;

             default:

                xassert(row != row);

          }

          mir->vlb[k] = mir->vub[k] = 0;

       }

       return;

 }

 static void set_col_attrib(glp_tree *tree, struct MIR *mir)

 {     /* set global column attributes */

       glp_prob *mip = tree->mip;

       int m = mir->m;

       int n = mir->n;

       int k;

       for (k = m+1; k <= m+n; k++)

       {  GLPCOL *col = mip->col[k-m];

          switch (col->kind)

          {  case GLP_CV:

                mir->isint[k] = 0; break;

             case GLP_IV:

                mir->isint[k] = 1; break;

             default:

                xassert(col != col);

          }

          switch (col->type)

          {  case GLP_FR:

                mir->lb[k] = -DBL_MAX, mir->ub[k] = +DBL_MAX; break;

             case GLP_LO:

                mir->lb[k] = col->lb, mir->ub[k] = +DBL_MAX; break;

             case GLP_UP:

                mir->lb[k] = -DBL_MAX, mir->ub[k] = col->ub; break;

             case GLP_DB:

                mir->lb[k] = col->lb, mir->ub[k] = col->ub; break;

             case GLP_FX:

                mir->lb[k] = mir->ub[k] = col->lb; break;

             default:

                xassert(col != col);

          }

          mir->vlb[k] = mir->vub[k] = 0;

       }

       return;

 }

 static void set_var_bounds(glp_tree *tree, struct MIR *mir)

 {     /* set variable bounds */

       glp_prob *mip = tree->mip;

       int m = mir->m;

       GLPAIJ *aij;

       int i, k1, k2;

       double a1, a2;

       for (i = 1; i <= m; i++)

       {  /* we need the row to be '>= 0' or '<= 0' */

          if (!(mir->lb[i] == 0.0 && mir->ub[i] == +DBL_MAX ||

                mir->lb[i] == -DBL_MAX && mir->ub[i] == 0.0)) continue;

          /* take first term */

          aij = mip->row[i]->ptr;

          if (aij == NULL) continue;

          k1 = m + aij->col->j, a1 = aij->val;

          /* take second term */

          aij = aij->r_next;

          if (aij == NULL) continue;

          k2 = m + aij->col->j, a2 = aij->val;

          /* there must be only two terms */

          if (aij->r_next != NULL) continue;

          /* interchange terms, if needed */

          if (!mir->isint[k1] && mir->isint[k2])

             ;

          else if (mir->isint[k1] && !mir->isint[k2])

          {  k2 = k1, a2 = a1;

             k1 = m + aij->col->j, a1 = aij->val;

          }

          else

          {  /* both terms are either continuous or integer */

             continue;

          }

          /* x[k2] should be double-bounded */

          if (mir->lb[k2] == -DBL_MAX || mir->ub[k2] == +DBL_MAX ||

              mir->lb[k2] == mir->ub[k2]) continue;

          /* change signs, if necessary */

          if (mir->ub[i] == 0.0) a1 = - a1, a2 = - a2;

          /* now the row has the form a1 * x1 + a2 * x2 >= 0, where x1

             is continuous, x2 is integer */

          if (a1 > 0.0)

          {  /* x1 >= - (a2 / a1) * x2 */

             if (mir->vlb[k1] == 0)

             {  /* set variable lower bound for x1 */

                mir->lb[k1] = - a2 / a1;

                mir->vlb[k1] = k2;

                /* the row should not be used */

                mir->skip[i] = 1;

             }

          }

          else /* a1 < 0.0 */

          {  /* x1 <= - (a2 / a1) * x2 */

             if (mir->vub[k1] == 0)

             {  /* set variable upper bound for x1 */

                mir->ub[k1] = - a2 / a1;

                mir->vub[k1] = k2;

                /* the row should not be used */

                mir->skip[i] = 1;

             }

          }

       }

       return;

 }

 static void mark_useless_rows(glp_tree *tree, struct MIR *mir)

 {     /* mark rows which should not be used */

       glp_prob *mip = tree->mip;

       int m = mir->m;

       GLPAIJ *aij;

       int i, k, nv;

       for (i = 1; i <= m; i++)

       {  /* free rows should not be used */

          if (mir->lb[i] == -DBL_MAX && mir->ub[i] == +DBL_MAX)

          {  mir->skip[i] = 1;

             continue;

          }

          nv = 0;

          for (aij = mip->row[i]->ptr; aij != NULL; aij = aij->r_next)

          {  k = m + aij->col->j;

             /* rows with free variables should not be used */

             if (mir->lb[k] == -DBL_MAX && mir->ub[k] == +DBL_MAX)

             {  mir->skip[i] = 1;

                break;

             }

             /* rows with integer variables having infinite (lower or

                upper) bound should not be used */

             if (mir->isint[k] && mir->lb[k] == -DBL_MAX ||

                 mir->isint[k] && mir->ub[k] == +DBL_MAX)

             {  mir->skip[i] = 1;

                break;

             }

             /* count non-fixed variables */

             if (!(mir->vlb[k] == 0 && mir->vub[k] == 0 &&

                   mir->lb[k] == mir->ub[k])) nv++;

          }

          /* rows with all variables fixed should not be used */

          if (nv == 0)

          {  mir->skip[i] = 1;

             continue;

          }

       }

       return;

 }

 void *ios_mir_init(glp_tree *tree)

 {     /* initialize MIR cut generator */

       glp_prob *mip = tree->mip;

       int m = mip->m;

       int n = mip->n;

       struct MIR *mir;

 #if _MIR_DEBUG

       xprintf("ios_mir_init: warning: debug mode enabled\n");

 #endif

       /* allocate working area */

       mir = xmalloc(sizeof(struct MIR));

       mir->m = m;

       mir->n = n;

       mir->skip = xcalloc(1+m, sizeof(char));

       mir->isint = xcalloc(1+m+n, sizeof(char));

       mir->lb = xcalloc(1+m+n, sizeof(double));

       mir->vlb = xcalloc(1+m+n, sizeof(int));

       mir->ub = xcalloc(1+m+n, sizeof(double));

       mir->vub = xcalloc(1+m+n, sizeof(int));

       mir->x = xcalloc(1+m+n, sizeof(double));

       mir->agg_row = xcalloc(1+MAXAGGR, sizeof(int));

       mir->agg_vec = ios_create_vec(m+n);

       mir->subst = xcalloc(1+m+n, sizeof(char));

       mir->mod_vec = ios_create_vec(m+n);

       mir->cut_vec = ios_create_vec(m+n);

       /* set global row attributes */

       set_row_attrib(tree, mir);

       /* set global column attributes */

       set_col_attrib(tree, mir);

       /* set variable bounds */

       set_var_bounds(tree, mir);

       /* mark rows which should not be used */

       mark_useless_rows(tree, mir);

       return mir;

 }

 /***********************************************************************

 *  NAME

 *

 *  ios_mir_gen - generate MIR cuts

 *

 *  SYNOPSIS

 *

 *  #include "glpios.h"

 *  void ios_mir_gen(glp_tree *tree, void *gen, IOSPOOL *pool);

 *

 *  DESCRIPTION

 *

 *  The routine ios_mir_gen generates MIR cuts for the current point and

 *  adds them to the cut pool. */

 static void get_current_point(glp_tree *tree, struct MIR *mir)

 {     /* obtain current point */

       glp_prob *mip = tree->mip;

       int m = mir->m;

       int n = mir->n;

       int k;

       for (k = 1; k <= m; k++)

          mir->x[k] = mip->row[k]->prim;

       for (k = m+1; k <= m+n; k++)

          mir->x[k] = mip->col[k-m]->prim;

       return;

 }

 #if _MIR_DEBUG

 static void check_current_point(struct MIR *mir)

 {     /* check current point */

       int m = mir->m;

       int n = mir->n;

       int k, kk;

       double lb, ub, eps;

       for (k = 1; k <= m+n; k++)

       {  /* determine lower bound */

          lb = mir->lb[k];

          kk = mir->vlb[k];

          if (kk != 0)

          {  xassert(lb != -DBL_MAX);

             xassert(!mir->isint[k]);

             xassert(mir->isint[kk]);

             lb *= mir->x[kk];

          }

          /* check lower bound */

          if (lb != -DBL_MAX)

          {  eps = 1e-6 * (1.0 + fabs(lb));

             xassert(mir->x[k] >= lb - eps);

          }

          /* determine upper bound */

          ub = mir->ub[k];

          kk = mir->vub[k];

          if (kk != 0)

          {  xassert(ub != +DBL_MAX);

             xassert(!mir->isint[k]);

             xassert(mir->isint[kk]);

             ub *= mir->x[kk];

          }

          /* check upper bound */

          if (ub != +DBL_MAX)

          {  eps = 1e-6 * (1.0 + fabs(ub));

             xassert(mir->x[k] <= ub + eps);

          }

       }

       return;

 }

 #endif

 static void initial_agg_row(glp_tree *tree, struct MIR *mir, int i)

 {     /* use original i-th row as initial aggregated constraint */

       glp_prob *mip = tree->mip;

       int m = mir->m;

       GLPAIJ *aij;

       xassert(1 <= i && i <= m);

       xassert(!mir->skip[i]);

       /* mark i-th row in order not to use it in the same aggregated

          constraint */

       mir->skip[i] = 2;

       mir->agg_cnt = 1;

       mir->agg_row[1] = i;

       /* use x[i] - sum a[i,j] * x[m+j] = 0, where x[i] is auxiliary

          variable of row i, x[m+j] are structural variables */

       ios_clear_vec(mir->agg_vec);

       ios_set_vj(mir->agg_vec, i, 1.0);

       for (aij = mip->row[i]->ptr; aij != NULL; aij = aij->r_next)

          ios_set_vj(mir->agg_vec, m + aij->col->j, - aij->val);

       mir->agg_rhs = 0.0;

 #if _MIR_DEBUG

       ios_check_vec(mir->agg_vec);

 #endif

       return;

 }

 #if _MIR_DEBUG

 static void check_agg_row(struct MIR *mir)

 {     /* check aggregated constraint */

       int m = mir->m;

       int n = mir->n;

       int j, k;

       double r, big;

       /* compute the residual r = sum a[k] * x[k] - b and determine

          big = max(1, |a[k]|, |b|) */

       r = 0.0, big = 1.0;

       for (j = 1; j <= mir->agg_vec->nnz; j++)

       {  k = mir->agg_vec->ind[j];

          xassert(1 <= k && k <= m+n);

          r += mir->agg_vec->val[j] * mir->x[k];

          if (big < fabs(mir->agg_vec->val[j]))

             big = fabs(mir->agg_vec->val[j]);

       }

       r -= mir->agg_rhs;

       if (big < fabs(mir->agg_rhs))

          big = fabs(mir->agg_rhs);

       /* the residual must be close to zero */

       xassert(fabs(r) <= 1e-6 * big);

       return;

 }

 #endif

 static void subst_fixed_vars(struct MIR *mir)

 {     /* substitute fixed variables into aggregated constraint */

       int m = mir->m;

       int n = mir->n;

       int j, k;

       for (j = 1; j <= mir->agg_vec->nnz; j++)

       {  k = mir->agg_vec->ind[j];

          xassert(1 <= k && k <= m+n);

          if (mir->vlb[k] == 0 && mir->vub[k] == 0 &&

              mir->lb[k] == mir->ub[k])

          {  /* x[k] is fixed */

             mir->agg_rhs -= mir->agg_vec->val[j] * mir->lb[k];

             mir->agg_vec->val[j] = 0.0;

          }

       }

       /* remove terms corresponding to fixed variables */

       ios_clean_vec(mir->agg_vec, DBL_EPSILON);

 #if _MIR_DEBUG

       ios_check_vec(mir->agg_vec);

 #endif

       return;

 }

 static void bound_subst_heur(struct MIR *mir)

 {     /* bound substitution heuristic */

       int m = mir->m;

       int n = mir->n;

       int j, k, kk;

       double d1, d2;

       for (j = 1; j <= mir->agg_vec->nnz; j++)

       {  k = mir->agg_vec->ind[j];

          xassert(1 <= k && k <= m+n);

          if (mir->isint[k]) continue; /* skip integer variable */

          /* compute distance from x[k] to its lower bound */

          kk = mir->vlb[k];

          if (kk == 0)

          {  if (mir->lb[k] == -DBL_MAX)

                d1 = DBL_MAX;

             else

                d1 = mir->x[k] - mir->lb[k];

          }

          else

          {  xassert(1 <= kk && kk <= m+n);

             xassert(mir->isint[kk]);

             xassert(mir->lb[k] != -DBL_MAX);

             d1 = mir->x[k] - mir->lb[k] * mir->x[kk];

          }

          /* compute distance from x[k] to its upper bound */

          kk = mir->vub[k];

          if (kk == 0)

          {  if (mir->vub[k] == +DBL_MAX)

                d2 = DBL_MAX;

             else

                d2 = mir->ub[k] - mir->x[k];

          }

          else

          {  xassert(1 <= kk && kk <= m+n);

             xassert(mir->isint[kk]);

             xassert(mir->ub[k] != +DBL_MAX);

             d2 = mir->ub[k] * mir->x[kk] - mir->x[k];

          }

          /* x[k] cannot be free */

          xassert(d1 != DBL_MAX || d2 != DBL_MAX);

          /* choose the bound which is closer to x[k] */

          xassert(mir->subst[k] == '?');

          if (d1 <= d2)

             mir->subst[k] = 'L';

          else

             mir->subst[k] = 'U';

       }

       return;

 }

 static void build_mod_row(struct MIR *mir)

 {     /* substitute bounds and build modified constraint */

       int m = mir->m;

       int n = mir->n;

       int j, jj, k, kk;

       /* initially modified constraint is aggregated constraint */

       ios_copy_vec(mir->mod_vec, mir->agg_vec);

       mir->mod_rhs = mir->agg_rhs;

 #if _MIR_DEBUG

       ios_check_vec(mir->mod_vec);

 #endif

       /* substitute bounds for continuous variables; note that due to

          substitution of variable bounds additional terms may appear in

          modified constraint */

       for (j = mir->mod_vec->nnz; j >= 1; j--)

       {  k = mir->mod_vec->ind[j];

          xassert(1 <= k && k <= m+n);

          if (mir->isint[k]) continue; /* skip integer variable */

          if (mir->subst[k] == 'L')

          {  /* x[k] = (lower bound) + x'[k] */

             xassert(mir->lb[k] != -DBL_MAX);

             kk = mir->vlb[k];

             if (kk == 0)

             {  /* x[k] = lb[k] + x'[k] */

                mir->mod_rhs -= mir->mod_vec->val[j] * mir->lb[k];

             }

             else

             {  /* x[k] = lb[k] * x[kk] + x'[k] */

                xassert(mir->isint[kk]);

                jj = mir->mod_vec->pos[kk];

                if (jj == 0)

                {  ios_set_vj(mir->mod_vec, kk, 1.0);

                   jj = mir->mod_vec->pos[kk];

                   mir->mod_vec->val[jj] = 0.0;

                }

                mir->mod_vec->val[jj] +=

                   mir->mod_vec->val[j] * mir->lb[k];

             }

          }

          else if (mir->subst[k] == 'U')

          {  /* x[k] = (upper bound) - x'[k] */

             xassert(mir->ub[k] != +DBL_MAX);

             kk = mir->vub[k];

             if (kk == 0)

             {  /* x[k] = ub[k] - x'[k] */

                mir->mod_rhs -= mir->mod_vec->val[j] * mir->ub[k];

             }

             else

             {  /* x[k] = ub[k] * x[kk] - x'[k] */

                xassert(mir->isint[kk]);

                jj = mir->mod_vec->pos[kk];

                if (jj == 0)

                {  ios_set_vj(mir->mod_vec, kk, 1.0);

                   jj = mir->mod_vec->pos[kk];

                   mir->mod_vec->val[jj] = 0.0;

                }

                mir->mod_vec->val[jj] +=

                   mir->mod_vec->val[j] * mir->ub[k];

             }

             mir->mod_vec->val[j] = - mir->mod_vec->val[j];

          }

          else

             xassert(k != k);

       }

 #if _MIR_DEBUG

       ios_check_vec(mir->mod_vec);

 #endif

       /* substitute bounds for integer variables */

       for (j = 1; j <= mir->mod_vec->nnz; j++)

       {  k = mir->mod_vec->ind[j];

          xassert(1 <= k && k <= m+n);

          if (!mir->isint[k]) continue; /* skip continuous variable */

          xassert(mir->subst[k] == '?');

          xassert(mir->vlb[k] == 0 && mir->vub[k] == 0);

          xassert(mir->lb[k] != -DBL_MAX && mir->ub[k] != +DBL_MAX);

          if (fabs(mir->lb[k]) <= fabs(mir->ub[k]))

          {  /* x[k] = lb[k] + x'[k] */

             mir->subst[k] = 'L';

             mir->mod_rhs -= mir->mod_vec->val[j] * mir->lb[k];

          }

          else

          {  /* x[k] = ub[k] - x'[k] */

             mir->subst[k] = 'U';

             mir->mod_rhs -= mir->mod_vec->val[j] * mir->ub[k];

             mir->mod_vec->val[j] = - mir->mod_vec->val[j];

          }

       }

 #if _MIR_DEBUG

       ios_check_vec(mir->mod_vec);

 #endif

       return;

 }

 #if _MIR_DEBUG

 static void check_mod_row(struct MIR *mir)

 {     /* check modified constraint */

       int m = mir->m;

       int n = mir->n;

       int j, k, kk;

       double r, big, x;

       /* compute the residual r = sum a'[k] * x'[k] - b' and determine

          big = max(1, |a[k]|, |b|) */

       r = 0.0, big = 1.0;

       for (j = 1; j <= mir->mod_vec->nnz; j++)

       {  k = mir->mod_vec->ind[j];

          xassert(1 <= k && k <= m+n);

          if (mir->subst[k] == 'L')

          {  /* x'[k] = x[k] - (lower bound) */

             xassert(mir->lb[k] != -DBL_MAX);

             kk = mir->vlb[k];

             if (kk == 0)

                x = mir->x[k] - mir->lb[k];

             else

                x = mir->x[k] - mir->lb[k] * mir->x[kk];

          }

          else if (mir->subst[k] == 'U')

          {  /* x'[k] = (upper bound) - x[k] */

             xassert(mir->ub[k] != +DBL_MAX);

             kk = mir->vub[k];

             if (kk == 0)

                x = mir->ub[k] - mir->x[k];

             else

                x = mir->ub[k] * mir->x[kk] - mir->x[k];

          }

          else

             xassert(k != k);

          r += mir->mod_vec->val[j] * x;

          if (big < fabs(mir->mod_vec->val[j]))

             big = fabs(mir->mod_vec->val[j]);

       }

       r -= mir->mod_rhs;

       if (big < fabs(mir->mod_rhs))

          big = fabs(mir->mod_rhs);

       /* the residual must be close to zero */

       xassert(fabs(r) <= 1e-6 * big);

       return;

 }

 #endif

 /***********************************************************************

 *  mir_ineq - construct MIR inequality

 *

 *  Given the single constraint mixed integer set

 *

 *                    |N|

 *     X = {(x,s) in Z    x R  : sum   a[j] * x[j] <= b + s},

 *                    +      +  j in N

 *

 *  this routine constructs the mixed integer rounding (MIR) inequality

 *

 *     sum   alpha[j] * x[j] <= beta + gamma * s,

 *    j in N

 *

 *  which is valid for X.

 *

 *  If the MIR inequality has been successfully constructed, the routine

 *  returns zero. Otherwise, if b is close to nearest integer, there may

 *  be numeric difficulties due to big coefficients; so in this case the

 *  routine returns non-zero. */

 static int mir_ineq(const int n, const double a[], const double b,

       double alpha[], double *beta, double *gamma)

 {     int j;

       double f, t;

       if (fabs(b - floor(b + .5)) < 0.01)

          return 1;

       f = b - floor(b);

       for (j = 1; j <= n; j++)

       {  t = (a[j] - floor(a[j])) - f;

          if (t <= 0.0)

             alpha[j] = floor(a[j]);

          else

             alpha[j] = floor(a[j]) + t / (1.0 - f);

       }

       *beta = floor(b);

       *gamma = 1.0 / (1.0 - f);

       return 0;

 }

 /***********************************************************************

 *  cmir_ineq - construct c-MIR inequality

 *

 *  Given the mixed knapsack set

 *

 *      MK              |N|

 *     X   = {(x,s) in Z    x R  : sum   a[j] * x[j] <= b + s,

 *                      +      +  j in N

 *

 *             x[j] <= u[j]},

 *

 *  a subset C of variables to be complemented, and a divisor delta > 0,

 *  this routine constructs the complemented MIR (c-MIR) inequality

 *

 *     sum   alpha[j] * x[j] <= beta + gamma * s,

 *    j in N

 *                      MK

 *  which is valid for X  .

 *

 *  If the c-MIR inequality has been successfully constructed, the

 *  routine returns zero. Otherwise, if there is a risk of numerical

 *  difficulties due to big coefficients (see comments to the routine

 *  mir_ineq), the routine cmir_ineq returns non-zero. */

 static int cmir_ineq(const int n, const double a[], const double b,

       const double u[], const char cset[], const double delta,

       double alpha[], double *beta, double *gamma)

 {     int j;

       double *aa, bb;

       aa = alpha, bb = b;

       for (j = 1; j <= n; j++)

       {  aa[j] = a[j] / delta;

          if (cset[j])

             aa[j] = - aa[j], bb -= a[j] * u[j];

       }

       bb /= delta;

       if (mir_ineq(n, aa, bb, alpha, beta, gamma)) return 1;

       for (j = 1; j <= n; j++)

       {  if (cset[j])

             alpha[j] = - alpha[j], *beta += alpha[j] * u[j];

       }

       *gamma /= delta;

       return 0;

 }

 /***********************************************************************

 *  cmir_sep - c-MIR separation heuristic

 *

 *  Given the mixed knapsack set

 *

 *      MK              |N|

 *     X   = {(x,s) in Z    x R  : sum   a[j] * x[j] <= b + s,

 *                      +      +  j in N

 *

 *             x[j] <= u[j]}

 *

 *                           *   *

 *  and a fractional point (x , s ), this routine tries to construct

 *  c-MIR inequality

 *

 *     sum   alpha[j] * x[j] <= beta + gamma * s,

 *    j in N

 *                      MK

 *  which is valid for X   and has (desirably maximal) violation at the

 *  fractional point given. This is attained by choosing an appropriate

 *  set C of variables to be complemented and a divisor delta > 0, which

 *  together define corresponding c-MIR inequality.

 *

 *  If a violated c-MIR inequality has been successfully constructed,

 *  the routine returns its violation:

 *

 *                       *                      *

 *     sum   alpha[j] * x [j] - beta - gamma * s ,

 *    j in N

 *

 *  which is positive. In case of failure the routine returns zero. */

 struct vset { int j; double v; };

 static int cmir_cmp(const void *p1, const void *p2)

 {     const struct vset *v1 = p1, *v2 = p2;

       if (v1->v < v2->v) return -1;

       if (v1->v > v2->v) return +1;

       return 0;

 }

 static double cmir_sep(const int n, const double a[], const double b,

       const double u[], const double x[], const double s,

       double alpha[], double *beta, double *gamma)

 {     int fail, j, k, nv, v;

       double delta, eps, d_try[1+3], r, r_best;

       char *cset;

       struct vset *vset;

       /* allocate working arrays */

       cset = xcalloc(1+n, sizeof(char));

       vset = xcalloc(1+n, sizeof(struct vset));

       /* choose initial C */

       for (j = 1; j <= n; j++)

          cset[j] = (char)(x[j] >= 0.5 * u[j]);

       /* choose initial delta */

       r_best = delta = 0.0;

       for (j = 1; j <= n; j++)

       {  xassert(a[j] != 0.0);

          /* if x[j] is close to its bounds, skip it */

          eps = 1e-9 * (1.0 + fabs(u[j]));

          if (x[j] < eps || x[j] > u[j] - eps) continue;

          /* try delta = |a[j]| to construct c-MIR inequality */

          fail = cmir_ineq(n, a, b, u, cset, fabs(a[j]), alpha, beta,

             gamma);

          if (fail) continue;

          /* compute violation */

          r = - (*beta) - (*gamma) * s;

          for (k = 1; k <= n; k++) r += alpha[k] * x[k];

          if (r_best < r) r_best = r, delta = fabs(a[j]);

       }

       if (r_best < 0.001) r_best = 0.0;

       if (r_best == 0.0) goto done;

       xassert(delta > 0.0);

       /* try to increase violation by dividing delta by 2, 4, and 8,

          respectively */

       d_try[1] = delta / 2.0;

       d_try[2] = delta / 4.0;

       d_try[3] = delta / 8.0;

       for (j = 1; j <= 3; j++)

       {  /* construct c-MIR inequality */

          fail = cmir_ineq(n, a, b, u, cset, d_try[j], alpha, beta,

             gamma);

          if (fail) continue;

          /* compute violation */

          r = - (*beta) - (*gamma) * s;

          for (k = 1; k <= n; k++) r += alpha[k] * x[k];

          if (r_best < r) r_best = r, delta = d_try[j];

       }

       /* build subset of variables lying strictly between their bounds

          and order it by nondecreasing values of |x[j] - u[j]/2| */

       nv = 0;

       for (j = 1; j <= n; j++)

       {  /* if x[j] is close to its bounds, skip it */

          eps = 1e-9 * (1.0 + fabs(u[j]));

          if (x[j] < eps || x[j] > u[j] - eps) continue;

          /* add x[j] to the subset */

          nv++;

          vset[nv].j = j;

          vset[nv].v = fabs(x[j] - 0.5 * u[j]);

       }

       qsort(&vset[1], nv, sizeof(struct vset), cmir_cmp);

       /* try to increase violation by successively complementing each

          variable in the subset */

       for (v = 1; v <= nv; v++)

       {  j = vset[v].j;

          /* replace x[j] by its complement or vice versa */

          cset[j] = (char)!cset[j];

          /* construct c-MIR inequality */

          fail = cmir_ineq(n, a, b, u, cset, delta, alpha, beta, gamma);

          /* restore the variable */

          cset[j] = (char)!cset[j];

          /* do not replace the variable in case of failure */

          if (fail) continue;

          /* compute violation */

          r = - (*beta) - (*gamma) * s;

          for (k = 1; k <= n; k++) r += alpha[k] * x[k];

          if (r_best < r) r_best = r, cset[j] = (char)!cset[j];

       }

       /* construct the best c-MIR inequality chosen */

       fail = cmir_ineq(n, a, b, u, cset, delta, alpha, beta, gamma);

       xassert(!fail);

 done: /* free working arrays */

       xfree(cset);

       xfree(vset);

       /* return to the calling routine */

       return r_best;

 }

 static double generate(struct MIR *mir)

 {     /* try to generate violated c-MIR cut for modified constraint */

       int m = mir->m;

       int n = mir->n;

       int j, k, kk, nint;

       double s, *u, *x, *alpha, r_best = 0.0, b, beta, gamma;

       ios_copy_vec(mir->cut_vec, mir->mod_vec);

       mir->cut_rhs = mir->mod_rhs;

       /* remove small terms, which can appear due to substitution of

          variable bounds */

       ios_clean_vec(mir->cut_vec, DBL_EPSILON);

 #if _MIR_DEBUG

       ios_check_vec(mir->cut_vec);

 #endif

       /* remove positive continuous terms to obtain MK relaxation */

       for (j = 1; j <= mir->cut_vec->nnz; j++)

       {  k = mir->cut_vec->ind[j];

          xassert(1 <= k && k <= m+n);

          if (!mir->isint[k] && mir->cut_vec->val[j] > 0.0)

             mir->cut_vec->val[j] = 0.0;

       }

       ios_clean_vec(mir->cut_vec, 0.0);

 #if _MIR_DEBUG

       ios_check_vec(mir->cut_vec);

 #endif

       /* move integer terms to the beginning of the sparse vector and

          determine the number of integer variables */

       nint = 0;

       for (j = 1; j <= mir->cut_vec->nnz; j++)

       {  k = mir->cut_vec->ind[j];

          xassert(1 <= k && k <= m+n);

          if (mir->isint[k])

          {  double temp;

             nint++;

             /* interchange elements [nint] and [j] */

             kk = mir->cut_vec->ind[nint];

             mir->cut_vec->pos[k] = nint;

             mir->cut_vec->pos[kk] = j;

             mir->cut_vec->ind[nint] = k;

             mir->cut_vec->ind[j] = kk;

             temp = mir->cut_vec->val[nint];

             mir->cut_vec->val[nint] = mir->cut_vec->val[j];

             mir->cut_vec->val[j] = temp;

          }

       }

 #if _MIR_DEBUG

       ios_check_vec(mir->cut_vec);

 #endif

       /* if there is no integer variable, nothing to generate */

       if (nint == 0) goto done;

       /* allocate working arrays */

       u = xcalloc(1+nint, sizeof(double));

       x = xcalloc(1+nint, sizeof(double));

       alpha = xcalloc(1+nint, sizeof(double));

       /* determine u and x */

       for (j = 1; j <= nint; j++)

       {  k = mir->cut_vec->ind[j];

          xassert(m+1 <= k && k <= m+n);

          xassert(mir->isint[k]);

          u[j] = mir->ub[k] - mir->lb[k];

          xassert(u[j] >= 1.0);

          if (mir->subst[k] == 'L')

             x[j] = mir->x[k] - mir->lb[k];

          else if (mir->subst[k] == 'U')

             x[j] = mir->ub[k] - mir->x[k];

          else

             xassert(k != k);

          xassert(x[j] >= -0.001);

          if (x[j] < 0.0) x[j] = 0.0;

       }

       /* compute s = - sum of continuous terms */

       s = 0.0;

       for (j = nint+1; j <= mir->cut_vec->nnz; j++)

       {  double x;

          k = mir->cut_vec->ind[j];

          xassert(1 <= k && k <= m+n);

          /* must be continuous */

          xassert(!mir->isint[k]);

          if (mir->subst[k] == 'L')

          {  xassert(mir->lb[k] != -DBL_MAX);

             kk = mir->vlb[k];

             if (kk == 0)

                x = mir->x[k] - mir->lb[k];

             else

                x = mir->x[k] - mir->lb[k] * mir->x[kk];

          }

          else if (mir->subst[k] == 'U')

          {  xassert(mir->ub[k] != +DBL_MAX);

             kk = mir->vub[k];

             if (kk == 0)

                x = mir->ub[k] - mir->x[k];

             else

                x = mir->ub[k] * mir->x[kk] - mir->x[k];

          }

          else

             xassert(k != k);

          xassert(x >= -0.001);

          if (x < 0.0) x = 0.0;

          s -= mir->cut_vec->val[j] * x;

       }

       xassert(s >= 0.0);

       /* apply heuristic to obtain most violated c-MIR inequality */

       b = mir->cut_rhs;

       r_best = cmir_sep(nint, mir->cut_vec->val, b, u, x, s, alpha,

          &beta, &gamma);

       if (r_best == 0.0) goto skip;

       xassert(r_best > 0.0);

       /* convert to raw cut */

       /* sum alpha[j] * x[j] <= beta + gamma * s */

       for (j = 1; j <= nint; j++)

          mir->cut_vec->val[j] = alpha[j];

       for (j = nint+1; j <= mir->cut_vec->nnz; j++)

       {  k = mir->cut_vec->ind[j];

          if (k <= m+n) mir->cut_vec->val[j] *= gamma;

       }

       mir->cut_rhs = beta;

 #if _MIR_DEBUG

       ios_check_vec(mir->cut_vec);

 #endif

 skip: /* free working arrays */

       xfree(u);

       xfree(x);

       xfree(alpha);

 done: return r_best;

 }

 #if _MIR_DEBUG

 static void check_raw_cut(struct MIR *mir, double r_best)

 {     /* check raw cut before back bound substitution */

       int m = mir->m;

       int n = mir->n;

       int j, k, kk;

       double r, big, x;

       /* compute the residual r = sum a[k] * x[k] - b and determine

          big = max(1, |a[k]|, |b|) */

       r = 0.0, big = 1.0;

       for (j = 1; j <= mir->cut_vec->nnz; j++)

       {  k = mir->cut_vec->ind[j];

          xassert(1 <= k && k <= m+n);

          if (mir->subst[k] == 'L')

          {  xassert(mir->lb[k] != -DBL_MAX);

             kk = mir->vlb[k];

             if (kk == 0)

                x = mir->x[k] - mir->lb[k];

             else

                x = mir->x[k] - mir->lb[k] * mir->x[kk];

          }

          else if (mir->subst[k] == 'U')

          {  xassert(mir->ub[k] != +DBL_MAX);

             kk = mir->vub[k];

             if (kk == 0)

                x = mir->ub[k] - mir->x[k];

             else

                x = mir->ub[k] * mir->x[kk] - mir->x[k];

          }

          else

             xassert(k != k);

          r += mir->cut_vec->val[j] * x;

          if (big < fabs(mir->cut_vec->val[j]))

             big = fabs(mir->cut_vec->val[j]);

       }

       r -= mir->cut_rhs;

       if (big < fabs(mir->cut_rhs))

          big = fabs(mir->cut_rhs);

       /* the residual must be close to r_best */

       xassert(fabs(r - r_best) <= 1e-6 * big);

       return;

 }

 #endif

 static void back_subst(struct MIR *mir)

 {     /* back substitution of original bounds */

       int m = mir->m;

       int n = mir->n;

       int j, jj, k, kk;

       /* at first, restore bounds of integer variables (because on

          restoring variable bounds of continuous variables we need

          original, not shifted, bounds of integer variables) */

       for (j = 1; j <= mir->cut_vec->nnz; j++)

       {  k = mir->cut_vec->ind[j];

          xassert(1 <= k && k <= m+n);

          if (!mir->isint[k]) continue; /* skip continuous */

          if (mir->subst[k] == 'L')

          {  /* x'[k] = x[k] - lb[k] */

             xassert(mir->lb[k] != -DBL_MAX);

             xassert(mir->vlb[k] == 0);

             mir->cut_rhs += mir->cut_vec->val[j] * mir->lb[k];

          }

          else if (mir->subst[k] == 'U')

          {  /* x'[k] = ub[k] - x[k] */

             xassert(mir->ub[k] != +DBL_MAX);

             xassert(mir->vub[k] == 0);

             mir->cut_rhs -= mir->cut_vec->val[j] * mir->ub[k];

             mir->cut_vec->val[j] = - mir->cut_vec->val[j];

          }

          else

             xassert(k != k);

       }

       /* now restore bounds of continuous variables */

       for (j = 1; j <= mir->cut_vec->nnz; j++)

       {  k = mir->cut_vec->ind[j];

          xassert(1 <= k && k <= m+n);

          if (mir->isint[k]) continue; /* skip integer */

          if (mir->subst[k] == 'L')

          {  /* x'[k] = x[k] - (lower bound) */

             xassert(mir->lb[k] != -DBL_MAX);

             kk = mir->vlb[k];

             if (kk == 0)

             {  /* x'[k] = x[k] - lb[k] */

                mir->cut_rhs += mir->cut_vec->val[j] * mir->lb[k];

             }

             else

             {  /* x'[k] = x[k] - lb[k] * x[kk] */

                jj = mir->cut_vec->pos[kk];

 #if 0

                xassert(jj != 0);

 #else

                if (jj == 0)

                {  ios_set_vj(mir->cut_vec, kk, 1.0);

                   jj = mir->cut_vec->pos[kk];

                   xassert(jj != 0);

                   mir->cut_vec->val[jj] = 0.0;

                }

 #endif

                mir->cut_vec->val[jj] -= mir->cut_vec->val[j] *

                   mir->lb[k];

             }

          }

          else if (mir->subst[k] == 'U')

          {  /* x'[k] = (upper bound) - x[k] */

             xassert(mir->ub[k] != +DBL_MAX);

             kk = mir->vub[k];

             if (kk == 0)

             {  /* x'[k] = ub[k] - x[k] */

                mir->cut_rhs -= mir->cut_vec->val[j] * mir->ub[k];

             }

             else

             {  /* x'[k] = ub[k] * x[kk] - x[k] */

                jj = mir->cut_vec->pos[kk];

                if (jj == 0)

                {  ios_set_vj(mir->cut_vec, kk, 1.0);

                   jj = mir->cut_vec->pos[kk];

                   xassert(jj != 0);

                   mir->cut_vec->val[jj] = 0.0;

                }

                mir->cut_vec->val[jj] += mir->cut_vec->val[j] *

                   mir->ub[k];

             }

             mir->cut_vec->val[j] = - mir->cut_vec->val[j];

          }

          else

             xassert(k != k);

       }

 #if _MIR_DEBUG

       ios_check_vec(mir->cut_vec);

 #endif

       return;

 }

 #if _MIR_DEBUG

 static void check_cut_row(struct MIR *mir, double r_best)

 {     /* check the cut after back bound substitution or elimination of

          auxiliary variables */

       int m = mir->m;

       int n = mir->n;

       int j, k;

       double r, big;

       /* compute the residual r = sum a[k] * x[k] - b and determine

          big = max(1, |a[k]|, |b|) */

       r = 0.0, big = 1.0;

       for (j = 1; j <= mir->cut_vec->nnz; j++)

       {  k = mir->cut_vec->ind[j];

          xassert(1 <= k && k <= m+n);

          r += mir->cut_vec->val[j] * mir->x[k];

          if (big < fabs(mir->cut_vec->val[j]))

             big = fabs(mir->cut_vec->val[j]);

       }

       r -= mir->cut_rhs;

       if (big < fabs(mir->cut_rhs))

          big = fabs(mir->cut_rhs);

       /* the residual must be close to r_best */

       xassert(fabs(r - r_best) <= 1e-6 * big);

       return;

 }

 #endif

 static void subst_aux_vars(glp_tree *tree, struct MIR *mir)

 {     /* final substitution to eliminate auxiliary variables */

       glp_prob *mip = tree->mip;

       int m = mir->m;

       int n = mir->n;

       GLPAIJ *aij;

       int j, k, kk, jj;

       for (j = mir->cut_vec->nnz; j >= 1; j--)

       {  k = mir->cut_vec->ind[j];

          xassert(1 <= k && k <= m+n);

          if (k > m) continue; /* skip structurals */

          for (aij = mip->row[k]->ptr; aij != NULL; aij = aij->r_next)

          {  kk = m + aij->col->j; /* structural */

             jj = mir->cut_vec->pos[kk];

             if (jj == 0)

             {  ios_set_vj(mir->cut_vec, kk, 1.0);

                jj = mir->cut_vec->pos[kk];

                mir->cut_vec->val[jj] = 0.0;

             }

             mir->cut_vec->val[jj] += mir->cut_vec->val[j] * aij->val;

          }

          mir->cut_vec->val[j] = 0.0;

       }

       ios_clean_vec(mir->cut_vec, 0.0);

       return;

 }

 static void add_cut(glp_tree *tree, struct MIR *mir)

 {     /* add constructed cut inequality to the cut pool */

       int m = mir->m;

       int n = mir->n;

       int j, k, len;

       int *ind = xcalloc(1+n, sizeof(int));

       double *val = xcalloc(1+n, sizeof(double));

       len = 0;

       for (j = mir->cut_vec->nnz; j >= 1; j--)

       {  k = mir->cut_vec->ind[j];

          xassert(m+1 <= k && k <= m+n);

          len++, ind[len] = k - m, val[len] = mir->cut_vec->val[j];

       }

 #if 0

       ios_add_cut_row(tree, pool, GLP_RF_MIR, len, ind, val, GLP_UP,

          mir->cut_rhs);

 #else

       glp_ios_add_row(tree, NULL, GLP_RF_MIR, 0, len, ind, val, GLP_UP,

          mir->cut_rhs);

 #endif

       xfree(ind);

       xfree(val);

       return;

 }

 static int aggregate_row(glp_tree *tree, struct MIR *mir)

 {     /* try to aggregate another row */

       glp_prob *mip = tree->mip;

       int m = mir->m;

       int n = mir->n;

       GLPAIJ *aij;

       IOSVEC *v;

       int ii, j, jj, k, kk, kappa = 0, ret = 0;

       double d1, d2, d, d_max = 0.0;

       /* choose appropriate structural variable in the aggregated row

          to be substituted */

       for (j = 1; j <= mir->agg_vec->nnz; j++)

       {  k = mir->agg_vec->ind[j];

          xassert(1 <= k && k <= m+n);

          if (k <= m) continue; /* skip auxiliary var */

          if (mir->isint[k]) continue; /* skip integer var */

          if (fabs(mir->agg_vec->val[j]) < 0.001) continue;

          /* compute distance from x[k] to its lower bound */

          kk = mir->vlb[k];

          if (kk == 0)

          {  if (mir->lb[k] == -DBL_MAX)

                d1 = DBL_MAX;

             else

                d1 = mir->x[k] - mir->lb[k];

          }

          else

          {  xassert(1 <= kk && kk <= m+n);

             xassert(mir->isint[kk]);

             xassert(mir->lb[k] != -DBL_MAX);

             d1 = mir->x[k] - mir->lb[k] * mir->x[kk];

          }

          /* compute distance from x[k] to its upper bound */

          kk = mir->vub[k];

          if (kk == 0)

          {  if (mir->vub[k] == +DBL_MAX)

                d2 = DBL_MAX;

             else

                d2 = mir->ub[k] - mir->x[k];

          }

          else

          {  xassert(1 <= kk && kk <= m+n);

             xassert(mir->isint[kk]);

             xassert(mir->ub[k] != +DBL_MAX);

             d2 = mir->ub[k] * mir->x[kk] - mir->x[k];

          }

          /* x[k] cannot be free */

          xassert(d1 != DBL_MAX || d2 != DBL_MAX);

          /* d = min(d1, d2) */

          d = (d1 <= d2 ? d1 : d2);

          xassert(d != DBL_MAX);

          /* should not be close to corresponding bound */

          if (d < 0.001) continue;

          if (d_max < d) d_max = d, kappa = k;

       }

       if (kappa == 0)

       {  /* nothing chosen */

          ret = 1;

          goto done;

       }

       /* x[kappa] has been chosen */

       xassert(m+1 <= kappa && kappa <= m+n);

       xassert(!mir->isint[kappa]);

       /* find another row, which have not been used yet, to eliminate

          x[kappa] from the aggregated row */

       for (ii = 1; ii <= m; ii++)

       {  if (mir->skip[ii]) continue;

          for (aij = mip->row[ii]->ptr; aij != NULL; aij = aij->r_next)

             if (aij->col->j == kappa - m) break;

          if (aij != NULL && fabs(aij->val) >= 0.001) break;

       }

       if (ii > m)

       {  /* nothing found */

          ret = 2;

          goto done;

       }

       /* row ii has been found; include it in the aggregated list */

       mir->agg_cnt++;

       xassert(mir->agg_cnt <= MAXAGGR);

       mir->agg_row[mir->agg_cnt] = ii;

       mir->skip[ii] = 2;

       /* v := new row */

       v = ios_create_vec(m+n);

       ios_set_vj(v, ii, 1.0);

       for (aij = mip->row[ii]->ptr; aij != NULL; aij = aij->r_next)

          ios_set_vj(v, m + aij->col->j, - aij->val);

 #if _MIR_DEBUG

       ios_check_vec(v);

 #endif

       /* perform gaussian elimination to remove x[kappa] */

       j = mir->agg_vec->pos[kappa];

       xassert(j != 0);

       jj = v->pos[kappa];

       xassert(jj != 0);

       ios_linear_comb(mir->agg_vec,

          - mir->agg_vec->val[j] / v->val[jj], v);

       ios_delete_vec(v);

       ios_set_vj(mir->agg_vec, kappa, 0.0);

 #if _MIR_DEBUG

       ios_check_vec(mir->agg_vec);

 #endif

 done: return ret;

 }

 void ios_mir_gen(glp_tree *tree, void *gen)

 {     /* main routine to generate MIR cuts */

       glp_prob *mip = tree->mip;

       struct MIR *mir = gen;

       int m = mir->m;

       int n = mir->n;

       int i;

       double r_best;

       xassert(mip->m >= m);

       xassert(mip->n == n);

       /* obtain current point */

       get_current_point(tree, mir);

 #if _MIR_DEBUG

       /* check current point */

       check_current_point(mir);

 #endif

       /* reset bound substitution flags */

       memset(&mir->subst[1], '?', m+n);

       /* try to generate a set of violated MIR cuts */

       for (i = 1; i <= m; i++)

       {  if (mir->skip[i]) continue;

          /* use original i-th row as initial aggregated constraint */

          initial_agg_row(tree, mir, i);

 loop:    ;

 #if _MIR_DEBUG

          /* check aggregated row */

          check_agg_row(mir);

 #endif

          /* substitute fixed variables into aggregated constraint */

          subst_fixed_vars(mir);

 #if _MIR_DEBUG

          /* check aggregated row */

          check_agg_row(mir);

 #endif

 #if _MIR_DEBUG

          /* check bound substitution flags */

          {  int k;

             for (k = 1; k <= m+n; k++)

                xassert(mir->subst[k] == '?');

          }

 #endif

          /* apply bound substitution heuristic */

          bound_subst_heur(mir);

          /* substitute bounds and build modified constraint */

          build_mod_row(mir);

 #if _MIR_DEBUG

          /* check modified row */

          check_mod_row(mir);

 #endif

          /* try to generate violated c-MIR cut for modified row */

          r_best = generate(mir);

          if (r_best > 0.0)

          {  /* success */

 #if _MIR_DEBUG

             /* check raw cut before back bound substitution */

             check_raw_cut(mir, r_best);

 #endif

             /* back substitution of original bounds */

             back_subst(mir);

 #if _MIR_DEBUG

             /* check the cut after back bound substitution */

             check_cut_row(mir, r_best);

 #endif

             /* final substitution to eliminate auxiliary variables */

             subst_aux_vars(tree, mir);

 #if _MIR_DEBUG

             /* check the cut after elimination of auxiliaries */

             check_cut_row(mir, r_best);

 #endif

             /* add constructed cut inequality to the cut pool */

             add_cut(tree, mir);

          }

          /* reset bound substitution flags */

          {  int j, k;

             for (j = 1; j <= mir->mod_vec->nnz; j++)

             {  k = mir->mod_vec->ind[j];

                xassert(1 <= k && k <= m+n);

                xassert(mir->subst[k] != '?');

                mir->subst[k] = '?';

             }

          }

          if (r_best == 0.0)

          {  /* failure */

             if (mir->agg_cnt < MAXAGGR)

             {  /* try to aggregate another row */

                if (aggregate_row(tree, mir) == 0) goto loop;

             }

          }

          /* unmark rows used in the aggregated constraint */

          {  int k, ii;

             for (k = 1; k <= mir->agg_cnt; k++)

             {  ii = mir->agg_row[k];

                xassert(1 <= ii && ii <= m);

                xassert(mir->skip[ii] == 2);

                mir->skip[ii] = 0;

             }

          }

       }

       return;

 }

 /***********************************************************************

 *  NAME

 *

 *  ios_mir_term - terminate MIR cut generator

 *

 *  SYNOPSIS

 *

 *  #include "glpios.h"

 *  void ios_mir_term(void *gen);

 *

 *  DESCRIPTION

 *

 *  The routine ios_mir_term deletes the MIR cut generator working area

 *  freeing all the memory allocated to it. */

 void ios_mir_term(void *gen)

 {     struct MIR *mir = gen;

       xfree(mir->skip);

       xfree(mir->isint);

       xfree(mir->lb);

       xfree(mir->vlb);

       xfree(mir->ub);

       xfree(mir->vub);

       xfree(mir->x);

       xfree(mir->agg_row);

       ios_delete_vec(mir->agg_vec);

       xfree(mir->subst);

       ios_delete_vec(mir->mod_vec);

       ios_delete_vec(mir->cut_vec);

       xfree(mir);

       return;

 }

 /* eof */