1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpios05.c	Sun Nov 06 20:59:10 2011 +0100
     1.3 @@ -0,0 +1,281 @@
     1.4 +/* glpios05.c (Gomory's mixed integer 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 +/***********************************************************************
    1.31 +*  NAME
    1.32 +*
    1.33 +*  ios_gmi_gen - generate Gomory's mixed integer cuts.
    1.34 +*
    1.35 +*  SYNOPSIS
    1.36 +*
    1.37 +*  #include "glpios.h"
    1.38 +*  void ios_gmi_gen(glp_tree *tree, IOSPOOL *pool);
    1.39 +*
    1.40 +*  DESCRIPTION
    1.41 +*
    1.42 +*  The routine ios_gmi_gen generates Gomory's mixed integer cuts for
    1.43 +*  the current point and adds them to the cut pool. */
    1.44 +
    1.45 +#define MAXCUTS 50
    1.46 +/* maximal number of cuts to be generated for one round */
    1.47 +
    1.48 +struct worka
    1.49 +{     /* Gomory's cut generator working area */
    1.50 +      int *ind; /* int ind[1+n]; */
    1.51 +      double *val; /* double val[1+n]; */
    1.52 +      double *phi; /* double phi[1+m+n]; */
    1.53 +};
    1.54 +
    1.55 +#define f(x) ((x) - floor(x))
    1.56 +/* compute fractional part of x */
    1.57 +
    1.58 +static void gen_cut(glp_tree *tree, struct worka *worka, int j)
    1.59 +{     /* this routine tries to generate Gomory's mixed integer cut for
    1.60 +         specified structural variable x[m+j] of integer kind, which is
    1.61 +         basic and has fractional value in optimal solution to current
    1.62 +         LP relaxation */
    1.63 +      glp_prob *mip = tree->mip;
    1.64 +      int m = mip->m;
    1.65 +      int n = mip->n;
    1.66 +      int *ind = worka->ind;
    1.67 +      double *val = worka->val;
    1.68 +      double *phi = worka->phi;
    1.69 +      int i, k, len, kind, stat;
    1.70 +      double lb, ub, alfa, beta, ksi, phi1, rhs;
    1.71 +      /* compute row of the simplex tableau, which (row) corresponds
    1.72 +         to specified basic variable xB[i] = x[m+j]; see (23) */
    1.73 +      len = glp_eval_tab_row(mip, m+j, ind, val);
    1.74 +      /* determine beta[i], which a value of xB[i] in optimal solution
    1.75 +         to current LP relaxation; note that this value is the same as
    1.76 +         if it would be computed with formula (27); it is assumed that
    1.77 +         beta[i] is fractional enough */
    1.78 +      beta = mip->col[j]->prim;
    1.79 +      /* compute cut coefficients phi and right-hand side rho, which
    1.80 +         correspond to formula (30); dense format is used, because rows
    1.81 +         of the simplex tableau is usually dense */
    1.82 +      for (k = 1; k <= m+n; k++) phi[k] = 0.0;
    1.83 +      rhs = f(beta); /* initial value of rho; see (28), (32) */
    1.84 +      for (j = 1; j <= len; j++)
    1.85 +      {  /* determine original number of non-basic variable xN[j] */
    1.86 +         k = ind[j];
    1.87 +         xassert(1 <= k && k <= m+n);
    1.88 +         /* determine the kind, bounds and current status of xN[j] in
    1.89 +            optimal solution to LP relaxation */
    1.90 +         if (k <= m)
    1.91 +         {  /* auxiliary variable */
    1.92 +            GLPROW *row = mip->row[k];
    1.93 +            kind = GLP_CV;
    1.94 +            lb = row->lb;
    1.95 +            ub = row->ub;
    1.96 +            stat = row->stat;
    1.97 +         }
    1.98 +         else
    1.99 +         {  /* structural variable */
   1.100 +            GLPCOL *col = mip->col[k-m];
   1.101 +            kind = col->kind;
   1.102 +            lb = col->lb;
   1.103 +            ub = col->ub;
   1.104 +            stat = col->stat;
   1.105 +         }
   1.106 +         /* xN[j] cannot be basic */
   1.107 +         xassert(stat != GLP_BS);
   1.108 +         /* determine row coefficient ksi[i,j] at xN[j]; see (23) */
   1.109 +         ksi = val[j];
   1.110 +         /* if ksi[i,j] is too large in the magnitude, do not generate
   1.111 +            the cut */
   1.112 +         if (fabs(ksi) > 1e+05) goto fini;
   1.113 +         /* if ksi[i,j] is too small in the magnitude, skip it */
   1.114 +         if (fabs(ksi) < 1e-10) goto skip;
   1.115 +         /* compute row coefficient alfa[i,j] at y[j]; see (26) */
   1.116 +         switch (stat)
   1.117 +         {  case GLP_NF:
   1.118 +               /* xN[j] is free (unbounded) having non-zero ksi[i,j];
   1.119 +                  do not generate the cut */
   1.120 +               goto fini;
   1.121 +            case GLP_NL:
   1.122 +               /* xN[j] has active lower bound */
   1.123 +               alfa = - ksi;
   1.124 +               break;
   1.125 +            case GLP_NU:
   1.126 +               /* xN[j] has active upper bound */
   1.127 +               alfa = + ksi;
   1.128 +               break;
   1.129 +            case GLP_NS:
   1.130 +               /* xN[j] is fixed; skip it */
   1.131 +               goto skip;
   1.132 +            default:
   1.133 +               xassert(stat != stat);
   1.134 +         }
   1.135 +         /* compute cut coefficient phi'[j] at y[j]; see (21), (28) */
   1.136 +         switch (kind)
   1.137 +         {  case GLP_IV:
   1.138 +               /* y[j] is integer */
   1.139 +               if (fabs(alfa - floor(alfa + 0.5)) < 1e-10)
   1.140 +               {  /* alfa[i,j] is close to nearest integer; skip it */
   1.141 +                  goto skip;
   1.142 +               }
   1.143 +               else if (f(alfa) <= f(beta))
   1.144 +                  phi1 = f(alfa);
   1.145 +               else
   1.146 +                  phi1 = (f(beta) / (1.0 - f(beta))) * (1.0 - f(alfa));
   1.147 +               break;
   1.148 +            case GLP_CV:
   1.149 +               /* y[j] is continuous */
   1.150 +               if (alfa >= 0.0)
   1.151 +                  phi1 = + alfa;
   1.152 +               else
   1.153 +                  phi1 = (f(beta) / (1.0 - f(beta))) * (- alfa);
   1.154 +               break;
   1.155 +            default:
   1.156 +               xassert(kind != kind);
   1.157 +         }
   1.158 +         /* compute cut coefficient phi[j] at xN[j] and update right-
   1.159 +            hand side rho; see (31), (32) */
   1.160 +         switch (stat)
   1.161 +         {  case GLP_NL:
   1.162 +               /* xN[j] has active lower bound */
   1.163 +               phi[k] = + phi1;
   1.164 +               rhs += phi1 * lb;
   1.165 +               break;
   1.166 +            case GLP_NU:
   1.167 +               /* xN[j] has active upper bound */
   1.168 +               phi[k] = - phi1;
   1.169 +               rhs -= phi1 * ub;
   1.170 +               break;
   1.171 +            default:
   1.172 +               xassert(stat != stat);
   1.173 +         }
   1.174 +skip:    ;
   1.175 +      }
   1.176 +      /* now the cut has the form sum_k phi[k] * x[k] >= rho, where cut
   1.177 +         coefficients are stored in the array phi in dense format;
   1.178 +         x[1,...,m] are auxiliary variables, x[m+1,...,m+n] are struc-
   1.179 +         tural variables; see (30) */
   1.180 +      /* eliminate auxiliary variables in order to express the cut only
   1.181 +         through structural variables; see (33) */
   1.182 +      for (i = 1; i <= m; i++)
   1.183 +      {  GLPROW *row;
   1.184 +         GLPAIJ *aij;
   1.185 +         if (fabs(phi[i]) < 1e-10) continue;
   1.186 +         /* auxiliary variable x[i] has non-zero cut coefficient */
   1.187 +         row = mip->row[i];
   1.188 +         /* x[i] cannot be fixed */
   1.189 +         xassert(row->type != GLP_FX);
   1.190 +         /* substitute x[i] = sum_j a[i,j] * x[m+j] */
   1.191 +         for (aij = row->ptr; aij != NULL; aij = aij->r_next)
   1.192 +            phi[m+aij->col->j] += phi[i] * aij->val;
   1.193 +      }
   1.194 +      /* convert the final cut to sparse format and substitute fixed
   1.195 +         (structural) variables */
   1.196 +      len = 0;
   1.197 +      for (j = 1; j <= n; j++)
   1.198 +      {  GLPCOL *col;
   1.199 +         if (fabs(phi[m+j]) < 1e-10) continue;
   1.200 +         /* structural variable x[m+j] has non-zero cut coefficient */
   1.201 +         col = mip->col[j];
   1.202 +         if (col->type == GLP_FX)
   1.203 +         {  /* eliminate x[m+j] */
   1.204 +            rhs -= phi[m+j] * col->lb;
   1.205 +         }
   1.206 +         else
   1.207 +         {  len++;
   1.208 +            ind[len] = j;
   1.209 +            val[len] = phi[m+j];
   1.210 +         }
   1.211 +      }
   1.212 +      if (fabs(rhs) < 1e-12) rhs = 0.0;
   1.213 +      /* if the cut inequality seems to be badly scaled, reject it to
   1.214 +         avoid numeric difficulties */
   1.215 +      for (k = 1; k <= len; k++)
   1.216 +      {  if (fabs(val[k]) < 1e-03) goto fini;
   1.217 +         if (fabs(val[k]) > 1e+03) goto fini;
   1.218 +      }
   1.219 +      /* add the cut to the cut pool for further consideration */
   1.220 +#if 0
   1.221 +      ios_add_cut_row(tree, pool, GLP_RF_GMI, len, ind, val, GLP_LO,
   1.222 +         rhs);
   1.223 +#else
   1.224 +      glp_ios_add_row(tree, NULL, GLP_RF_GMI, 0, len, ind, val, GLP_LO,
   1.225 +         rhs);
   1.226 +#endif
   1.227 +fini: return;
   1.228 +}
   1.229 +
   1.230 +struct var { int j; double f; };
   1.231 +
   1.232 +static int fcmp(const void *p1, const void *p2)
   1.233 +{     const struct var *v1 = p1, *v2 = p2;
   1.234 +      if (v1->f > v2->f) return -1;
   1.235 +      if (v1->f < v2->f) return +1;
   1.236 +      return 0;
   1.237 +}
   1.238 +
   1.239 +void ios_gmi_gen(glp_tree *tree)
   1.240 +{     /* main routine to generate Gomory's cuts */
   1.241 +      glp_prob *mip = tree->mip;
   1.242 +      int m = mip->m;
   1.243 +      int n = mip->n;
   1.244 +      struct var *var;
   1.245 +      int k, nv, j, size;
   1.246 +      struct worka _worka, *worka = &_worka;
   1.247 +      /* allocate working arrays */
   1.248 +      var = xcalloc(1+n, sizeof(struct var));
   1.249 +      worka->ind = xcalloc(1+n, sizeof(int));
   1.250 +      worka->val = xcalloc(1+n, sizeof(double));
   1.251 +      worka->phi = xcalloc(1+m+n, sizeof(double));
   1.252 +      /* build the list of integer structural variables, which are
   1.253 +         basic and have fractional value in optimal solution to current
   1.254 +         LP relaxation */
   1.255 +      nv = 0;
   1.256 +      for (j = 1; j <= n; j++)
   1.257 +      {  GLPCOL *col = mip->col[j];
   1.258 +         double frac;
   1.259 +         if (col->kind != GLP_IV) continue;
   1.260 +         if (col->type == GLP_FX) continue;
   1.261 +         if (col->stat != GLP_BS) continue;
   1.262 +         frac = f(col->prim);
   1.263 +         if (!(0.05 <= frac && frac <= 0.95)) continue;
   1.264 +         /* add variable to the list */
   1.265 +         nv++, var[nv].j = j, var[nv].f = frac;
   1.266 +      }
   1.267 +      /* order the list by descending fractionality */
   1.268 +      qsort(&var[1], nv, sizeof(struct var), fcmp);
   1.269 +      /* try to generate cuts by one for each variable in the list, but
   1.270 +         not more than MAXCUTS cuts */
   1.271 +      size = glp_ios_pool_size(tree);
   1.272 +      for (k = 1; k <= nv; k++)
   1.273 +      {  if (glp_ios_pool_size(tree) - size >= MAXCUTS) break;
   1.274 +         gen_cut(tree, worka, var[k].j);
   1.275 +      }
   1.276 +      /* free working arrays */
   1.277 +      xfree(var);
   1.278 +      xfree(worka->ind);
   1.279 +      xfree(worka->val);
   1.280 +      xfree(worka->phi);
   1.281 +      return;
   1.282 +}
   1.283 +
   1.284 +/* eof */