lemon-project-template-glpk: 33de93886c88 deps/glpk/src/glpios05.c

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 *

7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,

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 }

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 */