lemon-project-template-glpk
view deps/glpk/src/glpios05.c @ 9:33de93886c88
Import GLPK 4.47
author | Alpar Juttner <alpar@cs.elte.hu> |
---|---|
date | Sun, 06 Nov 2011 20:59:10 +0100 |
parents | |
children |
line source
1 /* glpios05.c (Gomory's mixed integer cut generator) */
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, 2011 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 ***********************************************************************/
25 #include "glpios.h"
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. */
42 #define MAXCUTS 50
43 /* maximal number of cuts to be generated for one round */
45 struct 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 };
52 #define f(x) ((x) - floor(x))
53 /* compute fractional part of x */
55 static 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 }
171 skip: ;
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
224 fini: return;
225 }
227 struct var { int j; double f; };
229 static 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 }
236 void 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 }
281 /* eof */