COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glpios05.c @ 1:c445c931472f

Last change on this file since 1:c445c931472f was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 14 years ago

Import glpk-4.45

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