lemon-project-template-glpk
comparison deps/glpk/src/glpios05.c @ 11:4fc6ad2fb8a6
Test GLPK in src/main.cc
author | Alpar Juttner <alpar@cs.elte.hu> |
---|---|
date | Sun, 06 Nov 2011 21:43:29 +0100 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:afa655fe5d7a |
---|---|
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, 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 ***********************************************************************/ | |
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 | |
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 }; | |
51 | |
52 #define f(x) ((x) - floor(x)) | |
53 /* compute fractional part of x */ | |
54 | |
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 } | |
226 | |
227 struct var { int j; double f; }; | |
228 | |
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 } | |
235 | |
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 } | |
280 | |
281 /* eof */ |