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

lemon-project-template-glpk

view deps/glpk/src/glpssx01.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 /* glpssx01.c */

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

26 #include "glpssx.h"

27 #define xfault xerror

29 /*----------------------------------------------------------------------

30 // ssx_create - create simplex solver workspace.

31 //

32 // This routine creates the workspace used by simplex solver routines,

33 // and returns a pointer to it.

34 //

35 // Parameters m, n, and nnz specify, respectively, the number of rows,

36 // columns, and non-zero constraint coefficients.

37 //

38 // This routine only allocates the memory for the workspace components,

39 // so the workspace needs to be saturated by data. */

41 SSX *ssx_create(int m, int n, int nnz)

42 { SSX *ssx;

43 int i, j, k;

44 if (m < 1)

45 xfault("ssx_create: m = %d; invalid number of rows\n", m);

46 if (n < 1)

47 xfault("ssx_create: n = %d; invalid number of columns\n", n);

48 if (nnz < 0)

49 xfault("ssx_create: nnz = %d; invalid number of non-zero const"

50 "raint coefficients\n", nnz);

51 ssx = xmalloc(sizeof(SSX));

52 ssx->m = m;

53 ssx->n = n;

54 ssx->type = xcalloc(1+m+n, sizeof(int));

55 ssx->lb = xcalloc(1+m+n, sizeof(mpq_t));

56 for (k = 1; k <= m+n; k++) mpq_init(ssx->lb[k]);

57 ssx->ub = xcalloc(1+m+n, sizeof(mpq_t));

58 for (k = 1; k <= m+n; k++) mpq_init(ssx->ub[k]);

59 ssx->coef = xcalloc(1+m+n, sizeof(mpq_t));

60 for (k = 0; k <= m+n; k++) mpq_init(ssx->coef[k]);

61 ssx->A_ptr = xcalloc(1+n+1, sizeof(int));

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

63 ssx->A_ind = xcalloc(1+nnz, sizeof(int));

64 ssx->A_val = xcalloc(1+nnz, sizeof(mpq_t));

65 for (k = 1; k <= nnz; k++) mpq_init(ssx->A_val[k]);

66 ssx->stat = xcalloc(1+m+n, sizeof(int));

67 ssx->Q_row = xcalloc(1+m+n, sizeof(int));

68 ssx->Q_col = xcalloc(1+m+n, sizeof(int));

69 ssx->binv = bfx_create_binv();

70 ssx->bbar = xcalloc(1+m, sizeof(mpq_t));

71 for (i = 0; i <= m; i++) mpq_init(ssx->bbar[i]);

72 ssx->pi = xcalloc(1+m, sizeof(mpq_t));

73 for (i = 1; i <= m; i++) mpq_init(ssx->pi[i]);

74 ssx->cbar = xcalloc(1+n, sizeof(mpq_t));

75 for (j = 1; j <= n; j++) mpq_init(ssx->cbar[j]);

76 ssx->rho = xcalloc(1+m, sizeof(mpq_t));

77 for (i = 1; i <= m; i++) mpq_init(ssx->rho[i]);

78 ssx->ap = xcalloc(1+n, sizeof(mpq_t));

79 for (j = 1; j <= n; j++) mpq_init(ssx->ap[j]);

80 ssx->aq = xcalloc(1+m, sizeof(mpq_t));

81 for (i = 1; i <= m; i++) mpq_init(ssx->aq[i]);

82 mpq_init(ssx->delta);

83 return ssx;

84 }

86 /*----------------------------------------------------------------------

87 // ssx_factorize - factorize the current basis matrix.

88 //

89 // This routine computes factorization of the current basis matrix B

90 // and returns the singularity flag. If the matrix B is non-singular,

91 // the flag is zero, otherwise non-zero. */

93 static int basis_col(void *info, int j, int ind[], mpq_t val[])

94 { /* this auxiliary routine provides row indices and numeric values

95 of non-zero elements in j-th column of the matrix B */

96 SSX *ssx = info;

97 int m = ssx->m;

98 int n = ssx->n;

99 int *A_ptr = ssx->A_ptr;

100 int *A_ind = ssx->A_ind;

101 mpq_t *A_val = ssx->A_val;

102 int *Q_col = ssx->Q_col;

103 int k, len, ptr;

104 xassert(1 <= j && j <= m);

105 k = Q_col[j]; /* x[k] = xB[j] */

106 xassert(1 <= k && k <= m+n);

107 /* j-th column of the matrix B is k-th column of the augmented

108 constraint matrix (I | -A) */

109 if (k <= m)

110 { /* it is a column of the unity matrix I */

111 len = 1, ind[1] = k, mpq_set_si(val[1], 1, 1);

112 }

113 else

114 { /* it is a column of the original constraint matrix -A */

115 len = 0;

116 for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++)

117 { len++;

118 ind[len] = A_ind[ptr];

119 mpq_neg(val[len], A_val[ptr]);

120 }

121 }

122 return len;

123 }

124

125 int ssx_factorize(SSX *ssx)

126 { int ret;

127 ret = bfx_factorize(ssx->binv, ssx->m, basis_col, ssx);

128 return ret;

129 }

130

131 /*----------------------------------------------------------------------

132 // ssx_get_xNj - determine value of non-basic variable.

133 //

134 // This routine determines the value of non-basic variable xN[j] in the

135 // current basic solution defined as follows:

136 //

137 // 0, if xN[j] is free variable

138 // lN[j], if xN[j] is on its lower bound

139 // uN[j], if xN[j] is on its upper bound

140 // lN[j] = uN[j], if xN[j] is fixed variable

141 //

142 // where lN[j] and uN[j] are lower and upper bounds of xN[j]. */

143

144 void ssx_get_xNj(SSX *ssx, int j, mpq_t x)

145 { int m = ssx->m;

146 int n = ssx->n;

147 mpq_t *lb = ssx->lb;

148 mpq_t *ub = ssx->ub;

149 int *stat = ssx->stat;

150 int *Q_col = ssx->Q_col;

151 int k;

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

153 k = Q_col[m+j]; /* x[k] = xN[j] */

154 xassert(1 <= k && k <= m+n);

155 switch (stat[k])

156 { case SSX_NL:

157 /* xN[j] is on its lower bound */

158 mpq_set(x, lb[k]); break;

159 case SSX_NU:

160 /* xN[j] is on its upper bound */

161 mpq_set(x, ub[k]); break;

162 case SSX_NF:

163 /* xN[j] is free variable */

164 mpq_set_si(x, 0, 1); break;

165 case SSX_NS:

166 /* xN[j] is fixed variable */

167 mpq_set(x, lb[k]); break;

168 default:

169 xassert(stat != stat);

170 }

171 return;

172 }

173

174 /*----------------------------------------------------------------------

175 // ssx_eval_bbar - compute values of basic variables.

176 //

177 // This routine computes values of basic variables xB in the current

178 // basic solution as follows:

179 //

180 // beta = - inv(B) * N * xN,

181 //

182 // where B is the basis matrix, N is the matrix of non-basic columns,

183 // xN is a vector of current values of non-basic variables. */

184

185 void ssx_eval_bbar(SSX *ssx)

186 { int m = ssx->m;

187 int n = ssx->n;

188 mpq_t *coef = ssx->coef;

189 int *A_ptr = ssx->A_ptr;

190 int *A_ind = ssx->A_ind;

191 mpq_t *A_val = ssx->A_val;

192 int *Q_col = ssx->Q_col;

193 mpq_t *bbar = ssx->bbar;

194 int i, j, k, ptr;

195 mpq_t x, temp;

196 mpq_init(x);

197 mpq_init(temp);

198 /* bbar := 0 */

199 for (i = 1; i <= m; i++)

200 mpq_set_si(bbar[i], 0, 1);

201 /* bbar := - N * xN = - N[1] * xN[1] - ... - N[n] * xN[n] */

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

203 { ssx_get_xNj(ssx, j, x);

204 if (mpq_sgn(x) == 0) continue;

205 k = Q_col[m+j]; /* x[k] = xN[j] */

206 if (k <= m)

207 { /* N[j] is a column of the unity matrix I */

208 mpq_sub(bbar[k], bbar[k], x);

209 }

210 else

211 { /* N[j] is a column of the original constraint matrix -A */

212 for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++)

213 { mpq_mul(temp, A_val[ptr], x);

214 mpq_add(bbar[A_ind[ptr]], bbar[A_ind[ptr]], temp);

215 }

216 }

217 }

218 /* bbar := inv(B) * bbar */

219 bfx_ftran(ssx->binv, bbar, 0);

220 #if 1

221 /* compute value of the objective function */

222 /* bbar[0] := c[0] */

223 mpq_set(bbar[0], coef[0]);

224 /* bbar[0] := bbar[0] + sum{i in B} cB[i] * xB[i] */

225 for (i = 1; i <= m; i++)

226 { k = Q_col[i]; /* x[k] = xB[i] */

227 if (mpq_sgn(coef[k]) == 0) continue;

228 mpq_mul(temp, coef[k], bbar[i]);

229 mpq_add(bbar[0], bbar[0], temp);

230 }

231 /* bbar[0] := bbar[0] + sum{j in N} cN[j] * xN[j] */

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

233 { k = Q_col[m+j]; /* x[k] = xN[j] */

234 if (mpq_sgn(coef[k]) == 0) continue;

235 ssx_get_xNj(ssx, j, x);

236 mpq_mul(temp, coef[k], x);

237 mpq_add(bbar[0], bbar[0], temp);

238 }

239 #endif

240 mpq_clear(x);

241 mpq_clear(temp);

242 return;

243 }

244

245 /*----------------------------------------------------------------------

246 // ssx_eval_pi - compute values of simplex multipliers.

247 //

248 // This routine computes values of simplex multipliers (shadow prices)

249 // pi in the current basic solution as follows:

250 //

251 // pi = inv(B') * cB,

252 //

253 // where B' is a matrix transposed to the basis matrix B, cB is a vector

254 // of objective coefficients at basic variables xB. */

255

256 void ssx_eval_pi(SSX *ssx)

257 { int m = ssx->m;

258 mpq_t *coef = ssx->coef;

259 int *Q_col = ssx->Q_col;

260 mpq_t *pi = ssx->pi;

261 int i;

262 /* pi := cB */

263 for (i = 1; i <= m; i++) mpq_set(pi[i], coef[Q_col[i]]);

264 /* pi := inv(B') * cB */

265 bfx_btran(ssx->binv, pi);

266 return;

267 }

268

269 /*----------------------------------------------------------------------

270 // ssx_eval_dj - compute reduced cost of non-basic variable.

271 //

272 // This routine computes reduced cost d[j] of non-basic variable xN[j]

273 // in the current basic solution as follows:

274 //

275 // d[j] = cN[j] - N[j] * pi,

276 //

277 // where cN[j] is an objective coefficient at xN[j], N[j] is a column

278 // of the augmented constraint matrix (I | -A) corresponding to xN[j],

279 // pi is the vector of simplex multipliers (shadow prices). */

280

281 void ssx_eval_dj(SSX *ssx, int j, mpq_t dj)

282 { int m = ssx->m;

283 int n = ssx->n;

284 mpq_t *coef = ssx->coef;

285 int *A_ptr = ssx->A_ptr;

286 int *A_ind = ssx->A_ind;

287 mpq_t *A_val = ssx->A_val;

288 int *Q_col = ssx->Q_col;

289 mpq_t *pi = ssx->pi;

290 int k, ptr, end;

291 mpq_t temp;

292 mpq_init(temp);

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

294 k = Q_col[m+j]; /* x[k] = xN[j] */

295 xassert(1 <= k && k <= m+n);

296 /* j-th column of the matrix N is k-th column of the augmented

297 constraint matrix (I | -A) */

298 if (k <= m)

299 { /* it is a column of the unity matrix I */

300 mpq_sub(dj, coef[k], pi[k]);

301 }

302 else

303 { /* it is a column of the original constraint matrix -A */

304 mpq_set(dj, coef[k]);

305 for (ptr = A_ptr[k-m], end = A_ptr[k-m+1]; ptr < end; ptr++)

306 { mpq_mul(temp, A_val[ptr], pi[A_ind[ptr]]);

307 mpq_add(dj, dj, temp);

308 }

309 }

310 mpq_clear(temp);

311 return;

312 }

313

314 /*----------------------------------------------------------------------

315 // ssx_eval_cbar - compute reduced costs of all non-basic variables.

316 //

317 // This routine computes the vector of reduced costs pi in the current

318 // basic solution for all non-basic variables, including fixed ones. */

319

320 void ssx_eval_cbar(SSX *ssx)

321 { int n = ssx->n;

322 mpq_t *cbar = ssx->cbar;

323 int j;

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

325 ssx_eval_dj(ssx, j, cbar[j]);

326 return;

327 }

328

329 /*----------------------------------------------------------------------

330 // ssx_eval_rho - compute p-th row of the inverse.

331 //

332 // This routine computes p-th row of the matrix inv(B), where B is the

333 // current basis matrix.

334 //

335 // p-th row of the inverse is computed using the following formula:

336 //

337 // rho = inv(B') * e[p],

338 //

339 // where B' is a matrix transposed to B, e[p] is a unity vector, which

340 // contains one in p-th position. */

341

342 void ssx_eval_rho(SSX *ssx)

343 { int m = ssx->m;

344 int p = ssx->p;

345 mpq_t *rho = ssx->rho;

346 int i;

347 xassert(1 <= p && p <= m);

348 /* rho := 0 */

349 for (i = 1; i <= m; i++) mpq_set_si(rho[i], 0, 1);

350 /* rho := e[p] */

351 mpq_set_si(rho[p], 1, 1);

352 /* rho := inv(B') * rho */

353 bfx_btran(ssx->binv, rho);

354 return;

355 }

356

357 /*----------------------------------------------------------------------

358 // ssx_eval_row - compute pivot row of the simplex table.

359 //

360 // This routine computes p-th (pivot) row of the current simplex table

361 // A~ = - inv(B) * N using the following formula:

362 //

363 // A~[p] = - N' * inv(B') * e[p] = - N' * rho[p],

364 //

365 // where N' is a matrix transposed to the matrix N, rho[p] is p-th row

366 // of the inverse inv(B). */

367

368 void ssx_eval_row(SSX *ssx)

369 { int m = ssx->m;

370 int n = ssx->n;

371 int *A_ptr = ssx->A_ptr;

372 int *A_ind = ssx->A_ind;

373 mpq_t *A_val = ssx->A_val;

374 int *Q_col = ssx->Q_col;

375 mpq_t *rho = ssx->rho;

376 mpq_t *ap = ssx->ap;

377 int j, k, ptr;

378 mpq_t temp;

379 mpq_init(temp);

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

381 { /* ap[j] := - N'[j] * rho (inner product) */

382 k = Q_col[m+j]; /* x[k] = xN[j] */

383 if (k <= m)

384 mpq_neg(ap[j], rho[k]);

385 else

386 { mpq_set_si(ap[j], 0, 1);

387 for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++)

388 { mpq_mul(temp, A_val[ptr], rho[A_ind[ptr]]);

389 mpq_add(ap[j], ap[j], temp);

390 }

391 }

392 }

393 mpq_clear(temp);

394 return;

395 }

396

397 /*----------------------------------------------------------------------

398 // ssx_eval_col - compute pivot column of the simplex table.

399 //

400 // This routine computes q-th (pivot) column of the current simplex

401 // table A~ = - inv(B) * N using the following formula:

402 //

403 // A~[q] = - inv(B) * N[q],

404 //

405 // where N[q] is q-th column of the matrix N corresponding to chosen

406 // non-basic variable xN[q]. */

407

408 void ssx_eval_col(SSX *ssx)

409 { int m = ssx->m;

410 int n = ssx->n;

411 int *A_ptr = ssx->A_ptr;

412 int *A_ind = ssx->A_ind;

413 mpq_t *A_val = ssx->A_val;

414 int *Q_col = ssx->Q_col;

415 int q = ssx->q;

416 mpq_t *aq = ssx->aq;

417 int i, k, ptr;

418 xassert(1 <= q && q <= n);

419 /* aq := 0 */

420 for (i = 1; i <= m; i++) mpq_set_si(aq[i], 0, 1);

421 /* aq := N[q] */

422 k = Q_col[m+q]; /* x[k] = xN[q] */

423 if (k <= m)

424 { /* N[q] is a column of the unity matrix I */

425 mpq_set_si(aq[k], 1, 1);

426 }

427 else

428 { /* N[q] is a column of the original constraint matrix -A */

429 for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++)

430 mpq_neg(aq[A_ind[ptr]], A_val[ptr]);

431 }

432 /* aq := inv(B) * aq */

433 bfx_ftran(ssx->binv, aq, 1);

434 /* aq := - aq */

435 for (i = 1; i <= m; i++) mpq_neg(aq[i], aq[i]);

436 return;

437 }

438

439 /*----------------------------------------------------------------------

440 // ssx_chuzc - choose pivot column.

441 //

442 // This routine chooses non-basic variable xN[q] whose reduced cost

443 // indicates possible improving of the objective function to enter it

444 // in the basis.

445 //

446 // Currently the standard (textbook) pricing is used, i.e. that

447 // non-basic variable is preferred which has greatest reduced cost (in

448 // magnitude).

449 //

450 // If xN[q] has been chosen, the routine stores its number q and also

451 // sets the flag q_dir that indicates direction in which xN[q] has to

452 // change (+1 means increasing, -1 means decreasing).

453 //

454 // If the choice cannot be made, because the current basic solution is

455 // dual feasible, the routine sets the number q to 0. */

456

457 void ssx_chuzc(SSX *ssx)

458 { int m = ssx->m;

459 int n = ssx->n;

460 int dir = (ssx->dir == SSX_MIN ? +1 : -1);

461 int *Q_col = ssx->Q_col;

462 int *stat = ssx->stat;

463 mpq_t *cbar = ssx->cbar;

464 int j, k, s, q, q_dir;

465 double best, temp;

466 /* nothing is chosen so far */

467 q = 0, q_dir = 0, best = 0.0;

468 /* look through the list of non-basic variables */

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

470 { k = Q_col[m+j]; /* x[k] = xN[j] */

471 s = dir * mpq_sgn(cbar[j]);

472 if ((stat[k] == SSX_NF || stat[k] == SSX_NL) && s < 0 ||

473 (stat[k] == SSX_NF || stat[k] == SSX_NU) && s > 0)

474 { /* reduced cost of xN[j] indicates possible improving of

475 the objective function */

476 temp = fabs(mpq_get_d(cbar[j]));

477 xassert(temp != 0.0);

478 if (q == 0 || best < temp)

479 q = j, q_dir = - s, best = temp;

480 }

481 }

482 ssx->q = q, ssx->q_dir = q_dir;

483 return;

484 }

485

486 /*----------------------------------------------------------------------

487 // ssx_chuzr - choose pivot row.

488 //

489 // This routine looks through elements of q-th column of the simplex

490 // table and chooses basic variable xB[p] which should leave the basis.

491 //

492 // The choice is based on the standard (textbook) ratio test.

493 //

494 // If xB[p] has been chosen, the routine stores its number p and also

495 // sets its non-basic status p_stat which should be assigned to xB[p]

496 // when it has left the basis and become xN[q].

497 //

498 // Special case p < 0 means that xN[q] is double-bounded variable and

499 // it reaches its opposite bound before any basic variable does that,

500 // so the current basis remains unchanged.

501 //

502 // If the choice cannot be made, because xN[q] can infinitely change in

503 // the feasible direction, the routine sets the number p to 0. */

504

505 void ssx_chuzr(SSX *ssx)

506 { int m = ssx->m;

507 int n = ssx->n;

508 int *type = ssx->type;

509 mpq_t *lb = ssx->lb;

510 mpq_t *ub = ssx->ub;

511 int *Q_col = ssx->Q_col;

512 mpq_t *bbar = ssx->bbar;

513 int q = ssx->q;

514 mpq_t *aq = ssx->aq;

515 int q_dir = ssx->q_dir;

516 int i, k, s, t, p, p_stat;

517 mpq_t teta, temp;

518 mpq_init(teta);

519 mpq_init(temp);

520 xassert(1 <= q && q <= n);

521 xassert(q_dir == +1 || q_dir == -1);

522 /* nothing is chosen so far */

523 p = 0, p_stat = 0;

524 /* look through the list of basic variables */

525 for (i = 1; i <= m; i++)

526 { s = q_dir * mpq_sgn(aq[i]);

527 if (s < 0)

528 { /* xB[i] decreases */

529 k = Q_col[i]; /* x[k] = xB[i] */

530 t = type[k];

531 if (t == SSX_LO || t == SSX_DB || t == SSX_FX)

532 { /* xB[i] has finite lower bound */

533 mpq_sub(temp, bbar[i], lb[k]);

534 mpq_div(temp, temp, aq[i]);

535 mpq_abs(temp, temp);

536 if (p == 0 || mpq_cmp(teta, temp) > 0)

537 { p = i;

538 p_stat = (t == SSX_FX ? SSX_NS : SSX_NL);

539 mpq_set(teta, temp);

540 }

541 }

542 }

543 else if (s > 0)

544 { /* xB[i] increases */

545 k = Q_col[i]; /* x[k] = xB[i] */

546 t = type[k];

547 if (t == SSX_UP || t == SSX_DB || t == SSX_FX)

548 { /* xB[i] has finite upper bound */

549 mpq_sub(temp, bbar[i], ub[k]);

550 mpq_div(temp, temp, aq[i]);

551 mpq_abs(temp, temp);

552 if (p == 0 || mpq_cmp(teta, temp) > 0)

553 { p = i;

554 p_stat = (t == SSX_FX ? SSX_NS : SSX_NU);

555 mpq_set(teta, temp);

556 }

557 }

558 }

559 /* if something has been chosen and the ratio test indicates

560 exact degeneracy, the search can be finished */

561 if (p != 0 && mpq_sgn(teta) == 0) break;

562 }

563 /* if xN[q] is double-bounded, check if it can reach its opposite

564 bound before any basic variable */

565 k = Q_col[m+q]; /* x[k] = xN[q] */

566 if (type[k] == SSX_DB)

567 { mpq_sub(temp, ub[k], lb[k]);

568 if (p == 0 || mpq_cmp(teta, temp) > 0)

569 { p = -1;

570 p_stat = -1;

571 mpq_set(teta, temp);

572 }

573 }

574 ssx->p = p;

575 ssx->p_stat = p_stat;

576 /* if xB[p] has been chosen, determine its actual change in the

577 adjacent basis (it has the same sign as q_dir) */

578 if (p != 0)

579 { xassert(mpq_sgn(teta) >= 0);

580 if (q_dir > 0)

581 mpq_set(ssx->delta, teta);

582 else

583 mpq_neg(ssx->delta, teta);

584 }

585 mpq_clear(teta);

586 mpq_clear(temp);

587 return;

588 }

589

590 /*----------------------------------------------------------------------

591 // ssx_update_bbar - update values of basic variables.

592 //

593 // This routine recomputes the current values of basic variables for

594 // the adjacent basis.

595 //

596 // The simplex table for the current basis is the following:

597 //

598 // xB[i] = sum{j in 1..n} alfa[i,j] * xN[q], i = 1,...,m

599 //

600 // therefore

601 //

602 // delta xB[i] = alfa[i,q] * delta xN[q], i = 1,...,m

603 //

604 // where delta xN[q] = xN.new[q] - xN[q] is the change of xN[q] in the

605 // adjacent basis, and delta xB[i] = xB.new[i] - xB[i] is the change of

606 // xB[i]. This gives formulae for recomputing values of xB[i]:

607 //

608 // xB.new[p] = xN[q] + delta xN[q]

609 //

610 // (because xN[q] becomes xB[p] in the adjacent basis), and

611 //

612 // xB.new[i] = xB[i] + alfa[i,q] * delta xN[q], i != p

613 //

614 // for other basic variables. */

615

616 void ssx_update_bbar(SSX *ssx)

617 { int m = ssx->m;

618 int n = ssx->n;

619 mpq_t *bbar = ssx->bbar;

620 mpq_t *cbar = ssx->cbar;

621 int p = ssx->p;

622 int q = ssx->q;

623 mpq_t *aq = ssx->aq;

624 int i;

625 mpq_t temp;

626 mpq_init(temp);

627 xassert(1 <= q && q <= n);

628 if (p < 0)

629 { /* xN[q] is double-bounded and goes to its opposite bound */

630 /* nop */;

631 }

632 else

633 { /* xN[q] becomes xB[p] in the adjacent basis */

634 /* xB.new[p] = xN[q] + delta xN[q] */

635 xassert(1 <= p && p <= m);

636 ssx_get_xNj(ssx, q, temp);

637 mpq_add(bbar[p], temp, ssx->delta);

638 }

639 /* update values of other basic variables depending on xN[q] */

640 for (i = 1; i <= m; i++)

641 { if (i == p) continue;

642 /* xB.new[i] = xB[i] + alfa[i,q] * delta xN[q] */

643 if (mpq_sgn(aq[i]) == 0) continue;

644 mpq_mul(temp, aq[i], ssx->delta);

645 mpq_add(bbar[i], bbar[i], temp);

646 }

647 #if 1

648 /* update value of the objective function */

649 /* z.new = z + d[q] * delta xN[q] */

650 mpq_mul(temp, cbar[q], ssx->delta);

651 mpq_add(bbar[0], bbar[0], temp);

652 #endif

653 mpq_clear(temp);

654 return;

655 }

656

657 /*----------------------------------------------------------------------

658 -- ssx_update_pi - update simplex multipliers.

659 --

660 -- This routine recomputes the vector of simplex multipliers for the

661 -- adjacent basis. */

662

663 void ssx_update_pi(SSX *ssx)

664 { int m = ssx->m;

665 int n = ssx->n;

666 mpq_t *pi = ssx->pi;

667 mpq_t *cbar = ssx->cbar;

668 int p = ssx->p;

669 int q = ssx->q;

670 mpq_t *aq = ssx->aq;

671 mpq_t *rho = ssx->rho;

672 int i;

673 mpq_t new_dq, temp;

674 mpq_init(new_dq);

675 mpq_init(temp);

676 xassert(1 <= p && p <= m);

677 xassert(1 <= q && q <= n);

678 /* compute d[q] in the adjacent basis */

679 mpq_div(new_dq, cbar[q], aq[p]);

680 /* update the vector of simplex multipliers */

681 for (i = 1; i <= m; i++)

682 { if (mpq_sgn(rho[i]) == 0) continue;

683 mpq_mul(temp, new_dq, rho[i]);

684 mpq_sub(pi[i], pi[i], temp);

685 }

686 mpq_clear(new_dq);

687 mpq_clear(temp);

688 return;

689 }

690

691 /*----------------------------------------------------------------------

692 // ssx_update_cbar - update reduced costs of non-basic variables.

693 //

694 // This routine recomputes the vector of reduced costs of non-basic

695 // variables for the adjacent basis. */

696

697 void ssx_update_cbar(SSX *ssx)

698 { int m = ssx->m;

699 int n = ssx->n;

700 mpq_t *cbar = ssx->cbar;

701 int p = ssx->p;

702 int q = ssx->q;

703 mpq_t *ap = ssx->ap;

704 int j;

705 mpq_t temp;

706 mpq_init(temp);

707 xassert(1 <= p && p <= m);

708 xassert(1 <= q && q <= n);

709 /* compute d[q] in the adjacent basis */

710 /* d.new[q] = d[q] / alfa[p,q] */

711 mpq_div(cbar[q], cbar[q], ap[q]);

712 /* update reduced costs of other non-basic variables */

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

714 { if (j == q) continue;

715 /* d.new[j] = d[j] - (alfa[p,j] / alfa[p,q]) * d[q] */

716 if (mpq_sgn(ap[j]) == 0) continue;

717 mpq_mul(temp, ap[j], cbar[q]);

718 mpq_sub(cbar[j], cbar[j], temp);

719 }

720 mpq_clear(temp);

721 return;

722 }

723

724 /*----------------------------------------------------------------------

725 // ssx_change_basis - change current basis to adjacent one.

726 //

727 // This routine changes the current basis to the adjacent one swapping

728 // basic variable xB[p] and non-basic variable xN[q]. */

729

730 void ssx_change_basis(SSX *ssx)

731 { int m = ssx->m;

732 int n = ssx->n;

733 int *type = ssx->type;

734 int *stat = ssx->stat;

735 int *Q_row = ssx->Q_row;

736 int *Q_col = ssx->Q_col;

737 int p = ssx->p;

738 int q = ssx->q;

739 int p_stat = ssx->p_stat;

740 int k, kp, kq;

741 if (p < 0)

742 { /* special case: xN[q] goes to its opposite bound */

743 xassert(1 <= q && q <= n);

744 k = Q_col[m+q]; /* x[k] = xN[q] */

745 xassert(type[k] == SSX_DB);

746 switch (stat[k])

747 { case SSX_NL:

748 stat[k] = SSX_NU;

749 break;

750 case SSX_NU:

751 stat[k] = SSX_NL;

752 break;

753 default:

754 xassert(stat != stat);

755 }

756 }

757 else

758 { /* xB[p] leaves the basis, xN[q] enters the basis */

759 xassert(1 <= p && p <= m);

760 xassert(1 <= q && q <= n);

761 kp = Q_col[p]; /* x[kp] = xB[p] */

762 kq = Q_col[m+q]; /* x[kq] = xN[q] */

763 /* check non-basic status of xB[p] which becomes xN[q] */

764 switch (type[kp])

765 { case SSX_FR:

766 xassert(p_stat == SSX_NF);

767 break;

768 case SSX_LO:

769 xassert(p_stat == SSX_NL);

770 break;

771 case SSX_UP:

772 xassert(p_stat == SSX_NU);

773 break;

774 case SSX_DB:

775 xassert(p_stat == SSX_NL || p_stat == SSX_NU);

776 break;

777 case SSX_FX:

778 xassert(p_stat == SSX_NS);

779 break;

780 default:

781 xassert(type != type);

782 }

783 /* swap xB[p] and xN[q] */

784 stat[kp] = (char)p_stat, stat[kq] = SSX_BS;

785 Q_row[kp] = m+q, Q_row[kq] = p;

786 Q_col[p] = kq, Q_col[m+q] = kp;

787 /* update factorization of the basis matrix */

788 if (bfx_update(ssx->binv, p))

789 { if (ssx_factorize(ssx))

790 xassert(("Internal error: basis matrix is singular", 0));

791 }

792 }

793 return;

794 }

795

796 /*----------------------------------------------------------------------

797 // ssx_delete - delete simplex solver workspace.

798 //

799 // This routine deletes the simplex solver workspace freeing all the

800 // memory allocated to this object. */

801

802 void ssx_delete(SSX *ssx)

803 { int m = ssx->m;

804 int n = ssx->n;

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

806 int i, j, k;

807 xfree(ssx->type);

808 for (k = 1; k <= m+n; k++) mpq_clear(ssx->lb[k]);

809 xfree(ssx->lb);

810 for (k = 1; k <= m+n; k++) mpq_clear(ssx->ub[k]);

811 xfree(ssx->ub);

812 for (k = 0; k <= m+n; k++) mpq_clear(ssx->coef[k]);

813 xfree(ssx->coef);

814 xfree(ssx->A_ptr);

815 xfree(ssx->A_ind);

816 for (k = 1; k <= nnz; k++) mpq_clear(ssx->A_val[k]);

817 xfree(ssx->A_val);

818 xfree(ssx->stat);

819 xfree(ssx->Q_row);

820 xfree(ssx->Q_col);

821 bfx_delete_binv(ssx->binv);

822 for (i = 0; i <= m; i++) mpq_clear(ssx->bbar[i]);

823 xfree(ssx->bbar);

824 for (i = 1; i <= m; i++) mpq_clear(ssx->pi[i]);

825 xfree(ssx->pi);

826 for (j = 1; j <= n; j++) mpq_clear(ssx->cbar[j]);

827 xfree(ssx->cbar);

828 for (i = 1; i <= m; i++) mpq_clear(ssx->rho[i]);

829 xfree(ssx->rho);

830 for (j = 1; j <= n; j++) mpq_clear(ssx->ap[j]);

831 xfree(ssx->ap);

832 for (i = 1; i <= m; i++) mpq_clear(ssx->aq[i]);

833 xfree(ssx->aq);

834 mpq_clear(ssx->delta);

835 xfree(ssx);

836 return;

837 }

838

839 /* eof */