lemon-project-template-glpk: 4fc6ad2fb8a6 deps/glpk/src/glpapi07.c

lemon-project-template-glpk

view deps/glpk/src/glpapi07.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

line source

1 /* glpapi07.c (exact simplex solver) */

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 "glpapi.h"

26 #include "glpssx.h"

28 /***********************************************************************

29 * NAME

30 *

31 * glp_exact - solve LP problem in exact arithmetic

32 *

33 * SYNOPSIS

34 *

35 * int glp_exact(glp_prob *lp, const glp_smcp *parm);

36 *

37 * DESCRIPTION

38 *

39 * The routine glp_exact is a tentative implementation of the primal

40 * two-phase simplex method based on exact (rational) arithmetic. It is

41 * similar to the routine glp_simplex, however, for all internal

42 * computations it uses arithmetic of rational numbers, which is exact

43 * in mathematical sense, i.e. free of round-off errors unlike floating

44 * point arithmetic.

45 *

46 * Note that the routine glp_exact uses inly two control parameters

47 * passed in the structure glp_smcp, namely, it_lim and tm_lim.

48 *

49 * RETURNS

50 *

51 * 0 The LP problem instance has been successfully solved. This code

52 * does not necessarily mean that the solver has found optimal

53 * solution. It only means that the solution process was successful.

54 *

55 * GLP_EBADB

56 * Unable to start the search, because the initial basis specified

57 * in the problem object is invalid--the number of basic (auxiliary

58 * and structural) variables is not the same as the number of rows in

59 * the problem object.

60 *

61 * GLP_ESING

62 * Unable to start the search, because the basis matrix correspodning

63 * to the initial basis is exactly singular.

64 *

65 * GLP_EBOUND

66 * Unable to start the search, because some double-bounded variables

67 * have incorrect bounds.

68 *

69 * GLP_EFAIL

70 * The problem has no rows/columns.

71 *

72 * GLP_EITLIM

73 * The search was prematurely terminated, because the simplex

74 * iteration limit has been exceeded.

75 *

76 * GLP_ETMLIM

77 * The search was prematurely terminated, because the time limit has

78 * been exceeded. */

80 static void set_d_eps(mpq_t x, double val)

81 { /* convert double val to rational x obtaining a more adequate

82 fraction than provided by mpq_set_d due to allowing a small

83 approximation error specified by a given relative tolerance;

84 for example, mpq_set_d would give the following

85 1/3 ~= 0.333333333333333314829616256247391... ->

86 -> 6004799503160661/18014398509481984

87 while this routine gives exactly 1/3 */

88 int s, n, j;

89 double f, p, q, eps = 1e-9;

90 mpq_t temp;

91 xassert(-DBL_MAX <= val && val <= +DBL_MAX);

92 #if 1 /* 30/VII-2008 */

93 if (val == floor(val))

94 { /* if val is integral, do not approximate */

95 mpq_set_d(x, val);

96 goto done;

97 }

98 #endif

99 if (val > 0.0)

100 s = +1;

101 else if (val < 0.0)

102 s = -1;

103 else

104 { mpq_set_si(x, 0, 1);

105 goto done;

106 }

107 f = frexp(fabs(val), &n);

108 /* |val| = f * 2^n, where 0.5 <= f < 1.0 */

109 fp2rat(f, 0.1 * eps, &p, &q);

110 /* f ~= p / q, where p and q are integers */

111 mpq_init(temp);

112 mpq_set_d(x, p);

113 mpq_set_d(temp, q);

114 mpq_div(x, x, temp);

115 mpq_set_si(temp, 1, 1);

116 for (j = 1; j <= abs(n); j++)

117 mpq_add(temp, temp, temp);

118 if (n > 0)

119 mpq_mul(x, x, temp);

120 else if (n < 0)

121 mpq_div(x, x, temp);

122 mpq_clear(temp);

123 if (s < 0) mpq_neg(x, x);

124 /* check that the desired tolerance has been attained */

125 xassert(fabs(val - mpq_get_d(x)) <= eps * (1.0 + fabs(val)));

126 done: return;

127 }

128

129 static void load_data(SSX *ssx, LPX *lp)

130 { /* load LP problem data into simplex solver workspace */

131 int m = ssx->m;

132 int n = ssx->n;

133 int nnz = ssx->A_ptr[n+1]-1;

134 int j, k, type, loc, len, *ind;

135 double lb, ub, coef, *val;

136 xassert(lpx_get_num_rows(lp) == m);

137 xassert(lpx_get_num_cols(lp) == n);

138 xassert(lpx_get_num_nz(lp) == nnz);

139 /* types and bounds of rows and columns */

140 for (k = 1; k <= m+n; k++)

141 { if (k <= m)

142 { type = lpx_get_row_type(lp, k);

143 lb = lpx_get_row_lb(lp, k);

144 ub = lpx_get_row_ub(lp, k);

145 }

146 else

147 { type = lpx_get_col_type(lp, k-m);

148 lb = lpx_get_col_lb(lp, k-m);

149 ub = lpx_get_col_ub(lp, k-m);

150 }

151 switch (type)

152 { case LPX_FR: type = SSX_FR; break;

153 case LPX_LO: type = SSX_LO; break;

154 case LPX_UP: type = SSX_UP; break;

155 case LPX_DB: type = SSX_DB; break;

156 case LPX_FX: type = SSX_FX; break;

157 default: xassert(type != type);

158 }

159 ssx->type[k] = type;

160 set_d_eps(ssx->lb[k], lb);

161 set_d_eps(ssx->ub[k], ub);

162 }

163 /* optimization direction */

164 switch (lpx_get_obj_dir(lp))

165 { case LPX_MIN: ssx->dir = SSX_MIN; break;

166 case LPX_MAX: ssx->dir = SSX_MAX; break;

167 default: xassert(lp != lp);

168 }

169 /* objective coefficients */

170 for (k = 0; k <= m+n; k++)

171 { if (k == 0)

172 coef = lpx_get_obj_coef(lp, 0);

173 else if (k <= m)

174 coef = 0.0;

175 else

176 coef = lpx_get_obj_coef(lp, k-m);

177 set_d_eps(ssx->coef[k], coef);

178 }

179 /* constraint coefficients */

180 ind = xcalloc(1+m, sizeof(int));

181 val = xcalloc(1+m, sizeof(double));

182 loc = 0;

183 for (j = 1; j <= n; j++)

184 { ssx->A_ptr[j] = loc+1;

185 len = lpx_get_mat_col(lp, j, ind, val);

186 for (k = 1; k <= len; k++)

187 { loc++;

188 ssx->A_ind[loc] = ind[k];

189 set_d_eps(ssx->A_val[loc], val[k]);

190 }

191 }

192 xassert(loc == nnz);

193 xfree(ind);

194 xfree(val);

195 return;

196 }

197

198 static int load_basis(SSX *ssx, LPX *lp)

199 { /* load current LP basis into simplex solver workspace */

200 int m = ssx->m;

201 int n = ssx->n;

202 int *type = ssx->type;

203 int *stat = ssx->stat;

204 int *Q_row = ssx->Q_row;

205 int *Q_col = ssx->Q_col;

206 int i, j, k;

207 xassert(lpx_get_num_rows(lp) == m);

208 xassert(lpx_get_num_cols(lp) == n);

209 /* statuses of rows and columns */

210 for (k = 1; k <= m+n; k++)

211 { if (k <= m)

212 stat[k] = lpx_get_row_stat(lp, k);

213 else

214 stat[k] = lpx_get_col_stat(lp, k-m);

215 switch (stat[k])

216 { case LPX_BS:

217 stat[k] = SSX_BS;

218 break;

219 case LPX_NL:

220 stat[k] = SSX_NL;

221 xassert(type[k] == SSX_LO || type[k] == SSX_DB);

222 break;

223 case LPX_NU:

224 stat[k] = SSX_NU;

225 xassert(type[k] == SSX_UP || type[k] == SSX_DB);

226 break;

227 case LPX_NF:

228 stat[k] = SSX_NF;

229 xassert(type[k] == SSX_FR);

230 break;

231 case LPX_NS:

232 stat[k] = SSX_NS;

233 xassert(type[k] == SSX_FX);

234 break;

235 default:

236 xassert(stat != stat);

237 }

238 }

239 /* build permutation matix Q */

240 i = j = 0;

241 for (k = 1; k <= m+n; k++)

242 { if (stat[k] == SSX_BS)

243 { i++;

244 if (i > m) return 1;

245 Q_row[k] = i, Q_col[i] = k;

246 }

247 else

248 { j++;

249 if (j > n) return 1;

250 Q_row[k] = m+j, Q_col[m+j] = k;

251 }

252 }

253 xassert(i == m && j == n);

254 return 0;

255 }

256

257 int glp_exact(glp_prob *lp, const glp_smcp *parm)

258 { glp_smcp _parm;

259 SSX *ssx;

260 int m = lpx_get_num_rows(lp);

261 int n = lpx_get_num_cols(lp);

262 int nnz = lpx_get_num_nz(lp);

263 int i, j, k, type, pst, dst, ret, *stat;

264 double lb, ub, *prim, *dual, sum;

265 if (parm == NULL)

266 parm = &_parm, glp_init_smcp((glp_smcp *)parm);

267 /* check control parameters */

268 if (parm->it_lim < 0)

269 xerror("glp_exact: it_lim = %d; invalid parameter\n",

270 parm->it_lim);

271 if (parm->tm_lim < 0)

272 xerror("glp_exact: tm_lim = %d; invalid parameter\n",

273 parm->tm_lim);

274 /* the problem must have at least one row and one column */

275 if (!(m > 0 && n > 0))

276 { xprintf("glp_exact: problem has no rows/columns\n");

277 return GLP_EFAIL;

278 }

279 #if 1

280 /* basic solution is currently undefined */

281 lp->pbs_stat = lp->dbs_stat = GLP_UNDEF;

282 lp->obj_val = 0.0;

283 lp->some = 0;

284 #endif

285 /* check that all double-bounded variables have correct bounds */

286 for (k = 1; k <= m+n; k++)

287 { if (k <= m)

288 { type = lpx_get_row_type(lp, k);

289 lb = lpx_get_row_lb(lp, k);

290 ub = lpx_get_row_ub(lp, k);

291 }

292 else

293 { type = lpx_get_col_type(lp, k-m);

294 lb = lpx_get_col_lb(lp, k-m);

295 ub = lpx_get_col_ub(lp, k-m);

296 }

297 if (type == LPX_DB && lb >= ub)

298 { xprintf("glp_exact: %s %d has invalid bounds\n",

299 k <= m ? "row" : "column", k <= m ? k : k-m);

300 return GLP_EBOUND;

301 }

302 }

303 /* create the simplex solver workspace */

304 xprintf("glp_exact: %d rows, %d columns, %d non-zeros\n",

305 m, n, nnz);

306 #ifdef HAVE_GMP

307 xprintf("GNU MP bignum library is being used\n");

308 #else

309 xprintf("GLPK bignum module is being used\n");

310 xprintf("(Consider installing GNU MP to attain a much better perf"

311 "ormance.)\n");

312 #endif

313 ssx = ssx_create(m, n, nnz);

314 /* load LP problem data into the workspace */

315 load_data(ssx, lp);

316 /* load current LP basis into the workspace */

317 if (load_basis(ssx, lp))

318 { xprintf("glp_exact: initial LP basis is invalid\n");

319 ret = GLP_EBADB;

320 goto done;

321 }

322 /* inherit some control parameters from the LP object */

323 #if 0

324 ssx->it_lim = lpx_get_int_parm(lp, LPX_K_ITLIM);

325 ssx->it_cnt = lpx_get_int_parm(lp, LPX_K_ITCNT);

326 ssx->tm_lim = lpx_get_real_parm(lp, LPX_K_TMLIM);

327 #else

328 ssx->it_lim = parm->it_lim;

329 ssx->it_cnt = lp->it_cnt;

330 ssx->tm_lim = (double)parm->tm_lim / 1000.0;

331 #endif

332 ssx->out_frq = 5.0;

333 ssx->tm_beg = xtime();

334 ssx->tm_lag = xlset(0);

335 /* solve LP */

336 ret = ssx_driver(ssx);

337 /* copy back some statistics to the LP object */

338 #if 0

339 lpx_set_int_parm(lp, LPX_K_ITLIM, ssx->it_lim);

340 lpx_set_int_parm(lp, LPX_K_ITCNT, ssx->it_cnt);

341 lpx_set_real_parm(lp, LPX_K_TMLIM, ssx->tm_lim);

342 #else

343 lp->it_cnt = ssx->it_cnt;

344 #endif

345 /* analyze the return code */

346 switch (ret)

347 { case 0:

348 /* optimal solution found */

349 ret = 0;

350 pst = LPX_P_FEAS, dst = LPX_D_FEAS;

351 break;

352 case 1:

353 /* problem has no feasible solution */

354 ret = 0;

355 pst = LPX_P_NOFEAS, dst = LPX_D_INFEAS;

356 break;

357 case 2:

358 /* problem has unbounded solution */

359 ret = 0;

360 pst = LPX_P_FEAS, dst = LPX_D_NOFEAS;

361 #if 1

362 xassert(1 <= ssx->q && ssx->q <= n);

363 lp->some = ssx->Q_col[m + ssx->q];

364 xassert(1 <= lp->some && lp->some <= m+n);

365 #endif

366 break;

367 case 3:

368 /* iteration limit exceeded (phase I) */

369 ret = GLP_EITLIM;

370 pst = LPX_P_INFEAS, dst = LPX_D_INFEAS;

371 break;

372 case 4:

373 /* iteration limit exceeded (phase II) */

374 ret = GLP_EITLIM;

375 pst = LPX_P_FEAS, dst = LPX_D_INFEAS;

376 break;

377 case 5:

378 /* time limit exceeded (phase I) */

379 ret = GLP_ETMLIM;

380 pst = LPX_P_INFEAS, dst = LPX_D_INFEAS;

381 break;

382 case 6:

383 /* time limit exceeded (phase II) */

384 ret = GLP_ETMLIM;

385 pst = LPX_P_FEAS, dst = LPX_D_INFEAS;

386 break;

387 case 7:

388 /* initial basis matrix is singular */

389 ret = GLP_ESING;

390 goto done;

391 default:

392 xassert(ret != ret);

393 }

394 /* obtain final basic solution components */

395 stat = xcalloc(1+m+n, sizeof(int));

396 prim = xcalloc(1+m+n, sizeof(double));

397 dual = xcalloc(1+m+n, sizeof(double));

398 for (k = 1; k <= m+n; k++)

399 { if (ssx->stat[k] == SSX_BS)

400 { i = ssx->Q_row[k]; /* x[k] = xB[i] */

401 xassert(1 <= i && i <= m);

402 stat[k] = LPX_BS;

403 prim[k] = mpq_get_d(ssx->bbar[i]);

404 dual[k] = 0.0;

405 }

406 else

407 { j = ssx->Q_row[k] - m; /* x[k] = xN[j] */

408 xassert(1 <= j && j <= n);

409 switch (ssx->stat[k])

410 { case SSX_NF:

411 stat[k] = LPX_NF;

412 prim[k] = 0.0;

413 break;

414 case SSX_NL:

415 stat[k] = LPX_NL;

416 prim[k] = mpq_get_d(ssx->lb[k]);

417 break;

418 case SSX_NU:

419 stat[k] = LPX_NU;

420 prim[k] = mpq_get_d(ssx->ub[k]);

421 break;

422 case SSX_NS:

423 stat[k] = LPX_NS;

424 prim[k] = mpq_get_d(ssx->lb[k]);

425 break;

426 default:

427 xassert(ssx != ssx);

428 }

429 dual[k] = mpq_get_d(ssx->cbar[j]);

430 }

431 }

432 /* and store them into the LP object */

433 pst = pst - LPX_P_UNDEF + GLP_UNDEF;

434 dst = dst - LPX_D_UNDEF + GLP_UNDEF;

435 for (k = 1; k <= m+n; k++)

436 stat[k] = stat[k] - LPX_BS + GLP_BS;

437 sum = lpx_get_obj_coef(lp, 0);

438 for (j = 1; j <= n; j++)

439 sum += lpx_get_obj_coef(lp, j) * prim[m+j];

440 lpx_put_solution(lp, 1, &pst, &dst, &sum,

441 &stat[0], &prim[0], &dual[0], &stat[m], &prim[m], &dual[m]);

442 xfree(stat);

443 xfree(prim);

444 xfree(dual);

445 done: /* delete the simplex solver workspace */

446 ssx_delete(ssx);

447 /* return to the application program */

448 return ret;

449 }

450

451 /* eof */