PageRenderTime 64ms CodeModel.GetById 43ms app.highlight 16ms RepoModel.GetById 2ms app.codeStats 0ms

/farmR/src/glpios05.c

https://code.google.com/p/javawfm/
C | 280 lines | 175 code | 11 blank | 94 comment | 43 complexity | 797a8bab4cca3ffb9510ac0755d78313 MD5 | raw file
  1/* glpios05.c (Gomory's mixed integer cut generator) */
  2
  3/***********************************************************************
  4*  This code is part of GLPK (GNU Linear Programming Kit).
  5*
  6*  Copyright (C) 2000,01,02,03,04,05,06,07,08,2009 Andrew Makhorin,
  7*  Department for Applied Informatics, Moscow Aviation Institute,
  8*  Moscow, Russia. All rights reserved. E-mail: <mao@mai2.rcnet.ru>.
  9*
 10*  GLPK is free software: you can redistribute it and/or modify it
 11*  under the terms of the GNU General Public License as published by
 12*  the Free Software Foundation, either version 3 of the License, or
 13*  (at your option) any later version.
 14*
 15*  GLPK is distributed in the hope that it will be useful, but WITHOUT
 16*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 17*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
 18*  License for more details.
 19*
 20*  You should have received a copy of the GNU General Public License
 21*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
 22***********************************************************************/
 23
 24#include "glpios.h"
 25
 26/***********************************************************************
 27*  NAME
 28*
 29*  ios_gmi_gen - generate Gomory's mixed integer cuts.
 30*
 31*  SYNOPSIS
 32*
 33*  #include "glpios.h"
 34*  void ios_gmi_gen(glp_tree *tree, IOSPOOL *pool);
 35*
 36*  DESCRIPTION
 37*
 38*  The routine ios_gmi_gen generates Gomory's mixed integer cuts for
 39*  the current point and adds them to the cut pool. */
 40
 41#define MAXCUTS 50
 42/* maximal number of cuts to be generated for one round */
 43
 44struct worka
 45{     /* Gomory's cut generator working area */
 46      int *ind; /* int ind[1+n]; */
 47      double *val; /* double val[1+n]; */
 48      double *phi; /* double phi[1+m+n]; */
 49};
 50
 51#define f(x) ((x) - floor(x))
 52/* compute fractional part of x */
 53
 54static void gen_cut(glp_tree *tree, struct worka *worka, int j)
 55{     /* this routine tries to generate Gomory's mixed integer cut for
 56         specified structural variable x[m+j] of integer kind, which is
 57         basic and has fractional value in optimal solution to current
 58         LP relaxation */
 59      glp_prob *mip = tree->mip;
 60      int m = mip->m;
 61      int n = mip->n;
 62      int *ind = worka->ind;
 63      double *val = worka->val;
 64      double *phi = worka->phi;
 65      int i, k, len, kind, stat;
 66      double lb, ub, alfa, beta, ksi, phi1, rhs;
 67      /* compute row of the simplex tableau, which (row) corresponds
 68         to specified basic variable xB[i] = x[m+j]; see (23) */
 69      len = glp_eval_tab_row(mip, m+j, ind, val);
 70      /* determine beta[i], which a value of xB[i] in optimal solution
 71         to current LP relaxation; note that this value is the same as
 72         if it would be computed with formula (27); it is assumed that
 73         beta[i] is fractional enough */
 74      beta = mip->col[j]->prim;
 75      /* compute cut coefficients phi and right-hand side rho, which
 76         correspond to formula (30); dense format is used, because rows
 77         of the simplex tableau is usually dense */
 78      for (k = 1; k <= m+n; k++) phi[k] = 0.0;
 79      rhs = f(beta); /* initial value of rho; see (28), (32) */
 80      for (j = 1; j <= len; j++)
 81      {  /* determine original number of non-basic variable xN[j] */
 82         k = ind[j];
 83         xassert(1 <= k && k <= m+n);
 84         /* determine the kind, bounds and current status of xN[j] in
 85            optimal solution to LP relaxation */
 86         if (k <= m)
 87         {  /* auxiliary variable */
 88            GLPROW *row = mip->row[k];
 89            kind = GLP_CV;
 90            lb = row->lb;
 91            ub = row->ub;
 92            stat = row->stat;
 93         }
 94         else
 95         {  /* structural variable */
 96            GLPCOL *col = mip->col[k-m];
 97            kind = col->kind;
 98            lb = col->lb;
 99            ub = col->ub;
100            stat = col->stat;
101         }
102         /* xN[j] cannot be basic */
103         xassert(stat != GLP_BS);
104         /* determine row coefficient ksi[i,j] at xN[j]; see (23) */
105         ksi = val[j];
106         /* if ksi[i,j] is too large in the magnitude, do not generate
107            the cut */
108         if (fabs(ksi) > 1e+05) goto fini;
109         /* if ksi[i,j] is too small in the magnitude, skip it */
110         if (fabs(ksi) < 1e-10) goto skip;
111         /* compute row coefficient alfa[i,j] at y[j]; see (26) */
112         switch (stat)
113         {  case GLP_NF:
114               /* xN[j] is free (unbounded) having non-zero ksi[i,j];
115                  do not generate the cut */
116               goto fini;
117            case GLP_NL:
118               /* xN[j] has active lower bound */
119               alfa = - ksi;
120               break;
121            case GLP_NU:
122               /* xN[j] has active upper bound */
123               alfa = + ksi;
124               break;
125            case GLP_NS:
126               /* xN[j] is fixed; skip it */
127               goto skip;
128            default:
129               xassert(stat != stat);
130         }
131         /* compute cut coefficient phi'[j] at y[j]; see (21), (28) */
132         switch (kind)
133         {  case GLP_IV:
134               /* y[j] is integer */
135               if (fabs(alfa - floor(alfa + 0.5)) < 1e-10)
136               {  /* alfa[i,j] is close to nearest integer; skip it */
137                  goto skip;
138               }
139               else if (f(alfa) <= f(beta))
140                  phi1 = f(alfa);
141               else
142                  phi1 = (f(beta) / (1.0 - f(beta))) * (1.0 - f(alfa));
143               break;
144            case GLP_CV:
145               /* y[j] is continuous */
146               if (alfa >= 0.0)
147                  phi1 = + alfa;
148               else
149                  phi1 = (f(beta) / (1.0 - f(beta))) * (- alfa);
150               break;
151            default:
152               xassert(kind != kind);
153         }
154         /* compute cut coefficient phi[j] at xN[j] and update right-
155            hand side rho; see (31), (32) */
156         switch (stat)
157         {  case GLP_NL:
158               /* xN[j] has active lower bound */
159               phi[k] = + phi1;
160               rhs += phi1 * lb;
161               break;
162            case GLP_NU:
163               /* xN[j] has active upper bound */
164               phi[k] = - phi1;
165               rhs -= phi1 * ub;
166               break;
167            default:
168               xassert(stat != stat);
169         }
170skip:    ;
171      }
172      /* now the cut has the form sum_k phi[k] * x[k] >= rho, where cut
173         coefficients are stored in the array phi in dense format;
174         x[1,...,m] are auxiliary variables, x[m+1,...,m+n] are struc-
175         tural variables; see (30) */
176      /* eliminate auxiliary variables in order to express the cut only
177         through structural variables; see (33) */
178      for (i = 1; i <= m; i++)
179      {  GLPROW *row;
180         GLPAIJ *aij;
181         if (fabs(phi[i]) < 1e-10) continue;
182         /* auxiliary variable x[i] has non-zero cut coefficient */
183         row = mip->row[i];
184         /* x[i] cannot be fixed */
185         xassert(row->type != GLP_FX);
186         /* substitute x[i] = sum_j a[i,j] * x[m+j] */
187         for (aij = row->ptr; aij != NULL; aij = aij->r_next)
188            phi[m+aij->col->j] += phi[i] * aij->val;
189      }
190      /* convert the final cut to sparse format and substitute fixed
191         (structural) variables */
192      len = 0;
193      for (j = 1; j <= n; j++)
194      {  GLPCOL *col;
195         if (fabs(phi[m+j]) < 1e-10) continue;
196         /* structural variable x[m+j] has non-zero cut coefficient */
197         col = mip->col[j];
198         if (col->type == GLP_FX)
199         {  /* eliminate x[m+j] */
200            rhs -= phi[m+j] * col->lb;
201         }
202         else
203         {  len++;
204            ind[len] = j;
205            val[len] = phi[m+j];
206         }
207      }
208      if (fabs(rhs) < 1e-12) rhs = 0.0;
209      /* if the cut inequality seems to be badly scaled, reject it to
210         avoid numeric difficulties */
211      for (k = 1; k <= len; k++)
212      {  if (fabs(val[k]) < 1e-03) goto fini;
213         if (fabs(val[k]) > 1e+03) goto fini;
214      }
215      /* add the cut to the cut pool for further consideration */
216#if 0
217      ios_add_cut_row(tree, pool, GLP_RF_GMI, len, ind, val, GLP_LO,
218         rhs);
219#else
220      glp_ios_add_row(tree, NULL, GLP_RF_GMI, 0, len, ind, val, GLP_LO,
221         rhs);
222#endif
223fini: return;
224}
225
226struct var { int j; double f; };
227
228static int fcmp(const void *p1, const void *p2)
229{     const struct var *v1 = p1, *v2 = p2;
230      if (v1->f > v2->f) return -1;
231      if (v1->f < v2->f) return +1;
232      return 0;
233}
234
235void ios_gmi_gen(glp_tree *tree)
236{     /* main routine to generate Gomory's cuts */
237      glp_prob *mip = tree->mip;
238      int m = mip->m;
239      int n = mip->n;
240      struct var *var;
241      int k, nv, j, size;
242      struct worka _worka, *worka = &_worka;
243      /* allocate working arrays */
244      var = xcalloc(1+n, sizeof(struct var));
245      worka->ind = xcalloc(1+n, sizeof(int));
246      worka->val = xcalloc(1+n, sizeof(double));
247      worka->phi = xcalloc(1+m+n, sizeof(double));
248      /* build the list of integer structural variables, which are
249         basic and have fractional value in optimal solution to current
250         LP relaxation */
251      nv = 0;
252      for (j = 1; j <= n; j++)
253      {  GLPCOL *col = mip->col[j];
254         double frac;
255         if (col->kind != GLP_IV) continue;
256         if (col->type == GLP_FX) continue;
257         if (col->stat != GLP_BS) continue;
258         frac = f(col->prim);
259         if (!(0.05 <= frac && frac <= 0.95)) continue;
260         /* add variable to the list */
261         nv++, var[nv].j = j, var[nv].f = frac;
262      }
263      /* order the list by descending fractionality */
264      qsort(&var[1], nv, sizeof(struct var), fcmp);
265      /* try to generate cuts by one for each variable in the list, but
266         not more than MAXCUTS cuts */
267      size = glp_ios_pool_size(tree);
268      for (k = 1; k <= nv; k++)
269      {  if (glp_ios_pool_size(tree) - size >= MAXCUTS) break;
270         gen_cut(tree, worka, var[k].j);
271      }
272      /* free working arrays */
273      xfree(var);
274      xfree(worka->ind);
275      xfree(worka->val);
276      xfree(worka->phi);
277      return;
278}
279
280/* eof */