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

lemon-project-template-glpk

view deps/glpk/src/glpnpp03.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 /* glpnpp03.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 "glpnpp.h"

27 /***********************************************************************

28 * NAME

29 *

30 * npp_empty_row - process empty row

31 *

32 * SYNOPSIS

33 *

34 * #include "glpnpp.h"

35 * int npp_empty_row(NPP *npp, NPPROW *p);

36 *

37 * DESCRIPTION

38 *

39 * The routine npp_empty_row processes row p, which is empty, i.e.

40 * coefficients at all columns in this row are zero:

41 *

42 * L[p] <= sum 0 x[j] <= U[p], (1)

43 *

44 * where L[p] <= U[p].

45 *

46 * RETURNS

47 *

48 * 0 - success;

49 *

50 * 1 - problem has no primal feasible solution.

51 *

52 * PROBLEM TRANSFORMATION

53 *

54 * If the following conditions hold:

55 *

56 * L[p] <= +eps, U[p] >= -eps, (2)

57 *

58 * where eps is an absolute tolerance for row value, the row p is

59 * redundant. In this case it can be replaced by equivalent redundant

60 * row, which is free (unbounded), and then removed from the problem.

61 * Otherwise, the row p is infeasible and, thus, the problem has no

62 * primal feasible solution.

63 *

64 * RECOVERING BASIC SOLUTION

65 *

66 * See the routine npp_free_row.

67 *

68 * RECOVERING INTERIOR-POINT SOLUTION

69 *

70 * See the routine npp_free_row.

71 *

72 * RECOVERING MIP SOLUTION

73 *

74 * None needed. */

76 int npp_empty_row(NPP *npp, NPPROW *p)

77 { /* process empty row */

78 double eps = 1e-3;

79 /* the row must be empty */

80 xassert(p->ptr == NULL);

81 /* check primal feasibility */

82 if (p->lb > +eps || p->ub < -eps)

83 return 1;

84 /* replace the row by equivalent free (unbounded) row */

85 p->lb = -DBL_MAX, p->ub = +DBL_MAX;

86 /* and process it */

87 npp_free_row(npp, p);

88 return 0;

89 }

91 /***********************************************************************

92 * NAME

93 *

94 * npp_empty_col - process empty column

95 *

96 * SYNOPSIS

97 *

98 * #include "glpnpp.h"

99 * int npp_empty_col(NPP *npp, NPPCOL *q);

100 *

101 * DESCRIPTION

102 *

103 * The routine npp_empty_col processes column q:

104 *

105 * l[q] <= x[q] <= u[q], (1)

106 *

107 * where l[q] <= u[q], which is empty, i.e. has zero coefficients in

108 * all constraint rows.

109 *

110 * RETURNS

111 *

112 * 0 - success;

113 *

114 * 1 - problem has no dual feasible solution.

115 *

116 * PROBLEM TRANSFORMATION

117 *

118 * The row of the dual system corresponding to the empty column is the

119 * following:

120 *

121 * sum 0 pi[i] + lambda[q] = c[q], (2)

122 * i

123 *

124 * from which it follows that:

125 *

126 * lambda[q] = c[q]. (3)

127 *

128 * If the following condition holds:

129 *

130 * c[q] < - eps, (4)

131 *

132 * where eps is an absolute tolerance for column multiplier, the lower

133 * column bound l[q] must be active to provide dual feasibility (note

134 * that being preprocessed the problem is always minimization). In this

135 * case the column can be fixed on its lower bound and removed from the

136 * problem (if the column is integral, its bounds are also assumed to

137 * be integral). And if the column has no lower bound (l[q] = -oo), the

138 * problem has no dual feasible solution.

139 *

140 * If the following condition holds:

141 *

142 * c[q] > + eps, (5)

143 *

144 * the upper column bound u[q] must be active to provide dual

145 * feasibility. In this case the column can be fixed on its upper bound

146 * and removed from the problem. And if the column has no upper bound

147 * (u[q] = +oo), the problem has no dual feasible solution.

148 *

149 * Finally, if the following condition holds:

150 *

151 * - eps <= c[q] <= +eps, (6)

152 *

153 * dual feasibility does not depend on a particular value of column q.

154 * In this case the column can be fixed either on its lower bound (if

155 * l[q] > -oo) or on its upper bound (if u[q] < +oo) or at zero (if the

156 * column is unbounded) and then removed from the problem.

157 *

158 * RECOVERING BASIC SOLUTION

159 *

160 * See the routine npp_fixed_col. Having been recovered the column

161 * is assigned status GLP_NS. However, if actually it is not fixed

162 * (l[q] < u[q]), its status should be changed to GLP_NL, GLP_NU, or

163 * GLP_NF depending on which bound it was fixed on transformation stage.

164 *

165 * RECOVERING INTERIOR-POINT SOLUTION

166 *

167 * See the routine npp_fixed_col.

168 *

169 * RECOVERING MIP SOLUTION

170 *

171 * See the routine npp_fixed_col. */

172

173 struct empty_col

174 { /* empty column */

175 int q;

176 /* column reference number */

177 char stat;

178 /* status in basic solution */

179 };

180

181 static int rcv_empty_col(NPP *npp, void *info);

182

183 int npp_empty_col(NPP *npp, NPPCOL *q)

184 { /* process empty column */

185 struct empty_col *info;

186 double eps = 1e-3;

187 /* the column must be empty */

188 xassert(q->ptr == NULL);

189 /* check dual feasibility */

190 if (q->coef > +eps && q->lb == -DBL_MAX)

191 return 1;

192 if (q->coef < -eps && q->ub == +DBL_MAX)

193 return 1;

194 /* create transformation stack entry */

195 info = npp_push_tse(npp,

196 rcv_empty_col, sizeof(struct empty_col));

197 info->q = q->j;

198 /* fix the column */

199 if (q->lb == -DBL_MAX && q->ub == +DBL_MAX)

200 { /* free column */

201 info->stat = GLP_NF;

202 q->lb = q->ub = 0.0;

203 }

204 else if (q->ub == +DBL_MAX)

205 lo: { /* column with lower bound */

206 info->stat = GLP_NL;

207 q->ub = q->lb;

208 }

209 else if (q->lb == -DBL_MAX)

210 up: { /* column with upper bound */

211 info->stat = GLP_NU;

212 q->lb = q->ub;

213 }

214 else if (q->lb != q->ub)

215 { /* double-bounded column */

216 if (q->coef >= +DBL_EPSILON) goto lo;

217 if (q->coef <= -DBL_EPSILON) goto up;

218 if (fabs(q->lb) <= fabs(q->ub)) goto lo; else goto up;

219 }

220 else

221 { /* fixed column */

222 info->stat = GLP_NS;

223 }

224 /* process fixed column */

225 npp_fixed_col(npp, q);

226 return 0;

227 }

228

229 static int rcv_empty_col(NPP *npp, void *_info)

230 { /* recover empty column */

231 struct empty_col *info = _info;

232 if (npp->sol == GLP_SOL)

233 npp->c_stat[info->q] = info->stat;

234 return 0;

235 }

236

237 /***********************************************************************

238 * NAME

239 *

240 * npp_implied_value - process implied column value

241 *

242 * SYNOPSIS

243 *

244 * #include "glpnpp.h"

245 * int npp_implied_value(NPP *npp, NPPCOL *q, double s);

246 *

247 * DESCRIPTION

248 *

249 * For column q:

250 *

251 * l[q] <= x[q] <= u[q], (1)

252 *

253 * where l[q] < u[q], the routine npp_implied_value processes its

254 * implied value s[q]. If this implied value satisfies to the current

255 * column bounds and integrality condition, the routine fixes column q

256 * at the given point. Note that the column is kept in the problem in

257 * any case.

258 *

259 * RETURNS

260 *

261 * 0 - column has been fixed;

262 *

263 * 1 - implied value violates to current column bounds;

264 *

265 * 2 - implied value violates integrality condition.

266 *

267 * ALGORITHM

268 *

269 * Implied column value s[q] satisfies to the current column bounds if

270 * the following condition holds:

271 *

272 * l[q] - eps <= s[q] <= u[q] + eps, (2)

273 *

274 * where eps is an absolute tolerance for column value. If the column

275 * is integral, the following condition also must hold:

276 *

277 * |s[q] - floor(s[q]+0.5)| <= eps, (3)

278 *

279 * where floor(s[q]+0.5) is the nearest integer to s[q].

280 *

281 * If both condition (2) and (3) are satisfied, the column can be fixed

282 * at the value s[q], or, if it is integral, at floor(s[q]+0.5).

283 * Otherwise, if s[q] violates (2) or (3), the problem has no feasible

284 * solution.

285 *

286 * Note: If s[q] is close to l[q] or u[q], it seems to be reasonable to

287 * fix the column at its lower or upper bound, resp. rather than at the

288 * implied value. */

289

290 int npp_implied_value(NPP *npp, NPPCOL *q, double s)

291 { /* process implied column value */

292 double eps, nint;

293 xassert(npp == npp);

294 /* column must not be fixed */

295 xassert(q->lb < q->ub);

296 /* check integrality */

297 if (q->is_int)

298 { nint = floor(s + 0.5);

299 if (fabs(s - nint) <= 1e-5)

300 s = nint;

301 else

302 return 2;

303 }

304 /* check current column lower bound */

305 if (q->lb != -DBL_MAX)

306 { eps = (q->is_int ? 1e-5 : 1e-5 + 1e-8 * fabs(q->lb));

307 if (s < q->lb - eps) return 1;

308 /* if s[q] is close to l[q], fix column at its lower bound

309 rather than at the implied value */

310 if (s < q->lb + 1e-3 * eps)

311 { q->ub = q->lb;

312 return 0;

313 }

314 }

315 /* check current column upper bound */

316 if (q->ub != +DBL_MAX)

317 { eps = (q->is_int ? 1e-5 : 1e-5 + 1e-8 * fabs(q->ub));

318 if (s > q->ub + eps) return 1;

319 /* if s[q] is close to u[q], fix column at its upper bound

320 rather than at the implied value */

321 if (s > q->ub - 1e-3 * eps)

322 { q->lb = q->ub;

323 return 0;

324 }

325 }

326 /* fix column at the implied value */

327 q->lb = q->ub = s;

328 return 0;

329 }

330

331 /***********************************************************************

332 * NAME

333 *

334 * npp_eq_singlet - process row singleton (equality constraint)

335 *

336 * SYNOPSIS

337 *

338 * #include "glpnpp.h"

339 * int npp_eq_singlet(NPP *npp, NPPROW *p);

340 *

341 * DESCRIPTION

342 *

343 * The routine npp_eq_singlet processes row p, which is equiality

344 * constraint having the only non-zero coefficient:

345 *

346 * a[p,q] x[q] = b. (1)

347 *

348 * RETURNS

349 *

350 * 0 - success;

351 *

352 * 1 - problem has no primal feasible solution;

353 *

354 * 2 - problem has no integer feasible solution.

355 *

356 * PROBLEM TRANSFORMATION

357 *

358 * The equality constraint defines implied value of column q:

359 *

360 * x[q] = s[q] = b / a[p,q]. (2)

361 *

362 * If the implied value s[q] satisfies to the column bounds (see the

363 * routine npp_implied_value), the column can be fixed at s[q] and

364 * removed from the problem. In this case row p becomes redundant, so

365 * it can be replaced by equivalent free row and also removed from the

366 * problem.

367 *

368 * Note that the routine removes from the problem only row p. Column q

369 * becomes fixed, however, it is kept in the problem.

370 *

371 * RECOVERING BASIC SOLUTION

372 *

373 * In solution to the original problem row p is assigned status GLP_NS

374 * (active equality constraint), and column q is assigned status GLP_BS

375 * (basic column).

376 *

377 * Multiplier for row p can be computed as follows. In the dual system

378 * of the original problem column q corresponds to the following row:

379 *

380 * sum a[i,q] pi[i] + lambda[q] = c[q] ==>

381 * i

382 *

383 * sum a[i,q] pi[i] + a[p,q] pi[p] + lambda[q] = c[q].

384 * i!=p

385 *

386 * Therefore:

387 *

388 * 1

389 * pi[p] = ------ (c[q] - lambda[q] - sum a[i,q] pi[i]), (3)

390 * a[p,q] i!=q

391 *

392 * where lambda[q] = 0 (since column[q] is basic), and pi[i] for all

393 * i != p are known in solution to the transformed problem.

394 *

395 * Value of column q in solution to the original problem is assigned

396 * its implied value s[q].

397 *

398 * RECOVERING INTERIOR-POINT SOLUTION

399 *

400 * Multiplier for row p is computed with formula (3). Value of column

401 * q is assigned its implied value s[q].

402 *

403 * RECOVERING MIP SOLUTION

404 *

405 * Value of column q is assigned its implied value s[q]. */

406

407 struct eq_singlet

408 { /* row singleton (equality constraint) */

409 int p;

410 /* row reference number */

411 int q;

412 /* column reference number */

413 double apq;

414 /* constraint coefficient a[p,q] */

415 double c;

416 /* objective coefficient at x[q] */

417 NPPLFE *ptr;

418 /* list of non-zero coefficients a[i,q], i != p */

419 };

420

421 static int rcv_eq_singlet(NPP *npp, void *info);

422

423 int npp_eq_singlet(NPP *npp, NPPROW *p)

424 { /* process row singleton (equality constraint) */

425 struct eq_singlet *info;

426 NPPCOL *q;

427 NPPAIJ *aij;

428 NPPLFE *lfe;

429 int ret;

430 double s;

431 /* the row must be singleton equality constraint */

432 xassert(p->lb == p->ub);

433 xassert(p->ptr != NULL && p->ptr->r_next == NULL);

434 /* compute and process implied column value */

435 aij = p->ptr;

436 q = aij->col;

437 s = p->lb / aij->val;

438 ret = npp_implied_value(npp, q, s);

439 xassert(0 <= ret && ret <= 2);

440 if (ret != 0) return ret;

441 /* create transformation stack entry */

442 info = npp_push_tse(npp,

443 rcv_eq_singlet, sizeof(struct eq_singlet));

444 info->p = p->i;

445 info->q = q->j;

446 info->apq = aij->val;

447 info->c = q->coef;

448 info->ptr = NULL;

449 /* save column coefficients a[i,q], i != p (not needed for MIP

450 solution) */

451 if (npp->sol != GLP_MIP)

452 { for (aij = q->ptr; aij != NULL; aij = aij->c_next)

453 { if (aij->row == p) continue; /* skip a[p,q] */

454 lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE));

455 lfe->ref = aij->row->i;

456 lfe->val = aij->val;

457 lfe->next = info->ptr;

458 info->ptr = lfe;

459 }

460 }

461 /* remove the row from the problem */

462 npp_del_row(npp, p);

463 return 0;

464 }

465

466 static int rcv_eq_singlet(NPP *npp, void *_info)

467 { /* recover row singleton (equality constraint) */

468 struct eq_singlet *info = _info;

469 NPPLFE *lfe;

470 double temp;

471 if (npp->sol == GLP_SOL)

472 { /* column q must be already recovered as GLP_NS */

473 if (npp->c_stat[info->q] != GLP_NS)

474 { npp_error();

475 return 1;

476 }

477 npp->r_stat[info->p] = GLP_NS;

478 npp->c_stat[info->q] = GLP_BS;

479 }

480 if (npp->sol != GLP_MIP)

481 { /* compute multiplier for row p with formula (3) */

482 temp = info->c;

483 for (lfe = info->ptr; lfe != NULL; lfe = lfe->next)

484 temp -= lfe->val * npp->r_pi[lfe->ref];

485 npp->r_pi[info->p] = temp / info->apq;

486 }

487 return 0;

488 }

489

490 /***********************************************************************

491 * NAME

492 *

493 * npp_implied_lower - process implied column lower bound

494 *

495 * SYNOPSIS

496 *

497 * #include "glpnpp.h"

498 * int npp_implied_lower(NPP *npp, NPPCOL *q, double l);

499 *

500 * DESCRIPTION

501 *

502 * For column q:

503 *

504 * l[q] <= x[q] <= u[q], (1)

505 *

506 * where l[q] < u[q], the routine npp_implied_lower processes its

507 * implied lower bound l'[q]. As the result the current column lower

508 * bound may increase. Note that the column is kept in the problem in

509 * any case.

510 *

511 * RETURNS

512 *

513 * 0 - current column lower bound has not changed;

514 *

515 * 1 - current column lower bound has changed, but not significantly;

516 *

517 * 2 - current column lower bound has significantly changed;

518 *

519 * 3 - column has been fixed on its upper bound;

520 *

521 * 4 - implied lower bound violates current column upper bound.

522 *

523 * ALGORITHM

524 *

525 * If column q is integral, before processing its implied lower bound

526 * should be rounded up:

527 *

528 * ( floor(l'[q]+0.5), if |l'[q] - floor(l'[q]+0.5)| <= eps

529 * l'[q] := < (2)

530 * ( ceil(l'[q]), otherwise

531 *

532 * where floor(l'[q]+0.5) is the nearest integer to l'[q], ceil(l'[q])

533 * is smallest integer not less than l'[q], and eps is an absolute

534 * tolerance for column value.

535 *

536 * Processing implied column lower bound l'[q] includes the following

537 * cases:

538 *

539 * 1) if l'[q] < l[q] + eps, implied lower bound is redundant;

540 *

541 * 2) if l[q] + eps <= l[q] <= u[q] + eps, current column lower bound

542 * l[q] can be strengthened by replacing it with l'[q]. If in this

543 * case new column lower bound becomes close to current column upper

544 * bound u[q], the column can be fixed on its upper bound;

545 *

546 * 3) if l'[q] > u[q] + eps, implied lower bound violates current

547 * column upper bound u[q], in which case the problem has no primal

548 * feasible solution. */

549

550 int npp_implied_lower(NPP *npp, NPPCOL *q, double l)

551 { /* process implied column lower bound */

552 int ret;

553 double eps, nint;

554 xassert(npp == npp);

555 /* column must not be fixed */

556 xassert(q->lb < q->ub);

557 /* implied lower bound must be finite */

558 xassert(l != -DBL_MAX);

559 /* if column is integral, round up l'[q] */

560 if (q->is_int)

561 { nint = floor(l + 0.5);

562 if (fabs(l - nint) <= 1e-5)

563 l = nint;

564 else

565 l = ceil(l);

566 }

567 /* check current column lower bound */

568 if (q->lb != -DBL_MAX)

569 { eps = (q->is_int ? 1e-3 : 1e-3 + 1e-6 * fabs(q->lb));

570 if (l < q->lb + eps)

571 { ret = 0; /* redundant */

572 goto done;

573 }

574 }

575 /* check current column upper bound */

576 if (q->ub != +DBL_MAX)

577 { eps = (q->is_int ? 1e-5 : 1e-5 + 1e-8 * fabs(q->ub));

578 if (l > q->ub + eps)

579 { ret = 4; /* infeasible */

580 goto done;

581 }

582 /* if l'[q] is close to u[q], fix column at its upper bound */

583 if (l > q->ub - 1e-3 * eps)

584 { q->lb = q->ub;

585 ret = 3; /* fixed */

586 goto done;

587 }

588 }

589 /* check if column lower bound changes significantly */

590 if (q->lb == -DBL_MAX)

591 ret = 2; /* significantly */

592 else if (q->is_int && l > q->lb + 0.5)

593 ret = 2; /* significantly */

594 else if (l > q->lb + 0.30 * (1.0 + fabs(q->lb)))

595 ret = 2; /* significantly */

596 else

597 ret = 1; /* not significantly */

598 /* set new column lower bound */

599 q->lb = l;

600 done: return ret;

601 }

602

603 /***********************************************************************

604 * NAME

605 *

606 * npp_implied_upper - process implied column upper bound

607 *

608 * SYNOPSIS

609 *

610 * #include "glpnpp.h"

611 * int npp_implied_upper(NPP *npp, NPPCOL *q, double u);

612 *

613 * DESCRIPTION

614 *

615 * For column q:

616 *

617 * l[q] <= x[q] <= u[q], (1)

618 *

619 * where l[q] < u[q], the routine npp_implied_upper processes its

620 * implied upper bound u'[q]. As the result the current column upper

621 * bound may decrease. Note that the column is kept in the problem in

622 * any case.

623 *

624 * RETURNS

625 *

626 * 0 - current column upper bound has not changed;

627 *

628 * 1 - current column upper bound has changed, but not significantly;

629 *

630 * 2 - current column upper bound has significantly changed;

631 *

632 * 3 - column has been fixed on its lower bound;

633 *

634 * 4 - implied upper bound violates current column lower bound.

635 *

636 * ALGORITHM

637 *

638 * If column q is integral, before processing its implied upper bound

639 * should be rounded down:

640 *

641 * ( floor(u'[q]+0.5), if |u'[q] - floor(l'[q]+0.5)| <= eps

642 * u'[q] := < (2)

643 * ( floor(l'[q]), otherwise

644 *

645 * where floor(u'[q]+0.5) is the nearest integer to u'[q],

646 * floor(u'[q]) is largest integer not greater than u'[q], and eps is

647 * an absolute tolerance for column value.

648 *

649 * Processing implied column upper bound u'[q] includes the following

650 * cases:

651 *

652 * 1) if u'[q] > u[q] - eps, implied upper bound is redundant;

653 *

654 * 2) if l[q] - eps <= u[q] <= u[q] - eps, current column upper bound

655 * u[q] can be strengthened by replacing it with u'[q]. If in this

656 * case new column upper bound becomes close to current column lower

657 * bound, the column can be fixed on its lower bound;

658 *

659 * 3) if u'[q] < l[q] - eps, implied upper bound violates current

660 * column lower bound l[q], in which case the problem has no primal

661 * feasible solution. */

662

663 int npp_implied_upper(NPP *npp, NPPCOL *q, double u)

664 { int ret;

665 double eps, nint;

666 xassert(npp == npp);

667 /* column must not be fixed */

668 xassert(q->lb < q->ub);

669 /* implied upper bound must be finite */

670 xassert(u != +DBL_MAX);

671 /* if column is integral, round down u'[q] */

672 if (q->is_int)

673 { nint = floor(u + 0.5);

674 if (fabs(u - nint) <= 1e-5)

675 u = nint;

676 else

677 u = floor(u);

678 }

679 /* check current column upper bound */

680 if (q->ub != +DBL_MAX)

681 { eps = (q->is_int ? 1e-3 : 1e-3 + 1e-6 * fabs(q->ub));

682 if (u > q->ub - eps)

683 { ret = 0; /* redundant */

684 goto done;

685 }

686 }

687 /* check current column lower bound */

688 if (q->lb != -DBL_MAX)

689 { eps = (q->is_int ? 1e-5 : 1e-5 + 1e-8 * fabs(q->lb));

690 if (u < q->lb - eps)

691 { ret = 4; /* infeasible */

692 goto done;

693 }

694 /* if u'[q] is close to l[q], fix column at its lower bound */

695 if (u < q->lb + 1e-3 * eps)

696 { q->ub = q->lb;

697 ret = 3; /* fixed */

698 goto done;

699 }

700 }

701 /* check if column upper bound changes significantly */

702 if (q->ub == +DBL_MAX)

703 ret = 2; /* significantly */

704 else if (q->is_int && u < q->ub - 0.5)

705 ret = 2; /* significantly */

706 else if (u < q->ub - 0.30 * (1.0 + fabs(q->ub)))

707 ret = 2; /* significantly */

708 else

709 ret = 1; /* not significantly */

710 /* set new column upper bound */

711 q->ub = u;

712 done: return ret;

713 }

714

715 /***********************************************************************

716 * NAME

717 *

718 * npp_ineq_singlet - process row singleton (inequality constraint)

719 *

720 * SYNOPSIS

721 *

722 * #include "glpnpp.h"

723 * int npp_ineq_singlet(NPP *npp, NPPROW *p);

724 *

725 * DESCRIPTION

726 *

727 * The routine npp_ineq_singlet processes row p, which is inequality

728 * constraint having the only non-zero coefficient:

729 *

730 * L[p] <= a[p,q] * x[q] <= U[p], (1)

731 *

732 * where L[p] < U[p], L[p] > -oo and/or U[p] < +oo.

733 *

734 * RETURNS

735 *

736 * 0 - current column bounds have not changed;

737 *

738 * 1 - current column bounds have changed, but not significantly;

739 *

740 * 2 - current column bounds have significantly changed;

741 *

742 * 3 - column has been fixed on its lower or upper bound;

743 *

744 * 4 - problem has no primal feasible solution.

745 *

746 * PROBLEM TRANSFORMATION

747 *

748 * Inequality constraint (1) defines implied bounds of column q:

749 *

750 * ( L[p] / a[p,q], if a[p,q] > 0

751 * l'[q] = < (2)

752 * ( U[p] / a[p,q], if a[p,q] < 0

753 *

754 * ( U[p] / a[p,q], if a[p,q] > 0

755 * u'[q] = < (3)

756 * ( L[p] / a[p,q], if a[p,q] < 0

757 *

758 * If these implied bounds do not violate current bounds of column q:

759 *

760 * l[q] <= x[q] <= u[q], (4)

761 *

762 * they can be used to strengthen the current column bounds:

763 *

764 * l[q] := max(l[q], l'[q]), (5)

765 *

766 * u[q] := min(u[q], u'[q]). (6)

767 *

768 * (See the routines npp_implied_lower and npp_implied_upper.)

769 *

770 * Once bounds of row p (1) have been carried over column q, the row

771 * becomes redundant, so it can be replaced by equivalent free row and

772 * removed from the problem.

773 *

774 * Note that the routine removes from the problem only row p. Column q,

775 * even it has been fixed, is kept in the problem.

776 *

777 * RECOVERING BASIC SOLUTION

778 *

779 * Note that the row in the dual system corresponding to column q is

780 * the following:

781 *

782 * sum a[i,q] pi[i] + lambda[q] = c[q] ==>

783 * i

784 * (7)

785 * sum a[i,q] pi[i] + a[p,q] pi[p] + lambda[q] = c[q],

786 * i!=p

787 *

788 * where pi[i] for all i != p are known in solution to the transformed

789 * problem. Row p does not exist in the transformed problem, so it has

790 * zero multiplier there. This allows computing multiplier for column q

791 * in solution to the transformed problem:

792 *

793 * lambda~[q] = c[q] - sum a[i,q] pi[i]. (8)

794 * i!=p

795 *

796 * Let in solution to the transformed problem column q be non-basic

797 * with lower bound active (GLP_NL, lambda~[q] >= 0), and this lower

798 * bound be implied one l'[q]. From the original problem's standpoint

799 * this then means that actually the original column lower bound l[q]

800 * is inactive, and active is that row bound L[p] or U[p] that defines

801 * the implied bound l'[q] (2). In this case in solution to the

802 * original problem column q is assigned status GLP_BS while row p is

803 * assigned status GLP_NL (if a[p,q] > 0) or GLP_NU (if a[p,q] < 0).

804 * Since now column q is basic, its multiplier lambda[q] is zero. This

805 * allows using (7) and (8) to find multiplier for row p in solution to

806 * the original problem:

807 *

808 * 1

809 * pi[p] = ------ (c[q] - sum a[i,q] pi[i]) = lambda~[q] / a[p,q] (9)

810 * a[p,q] i!=p

811 *

812 * Now let in solution to the transformed problem column q be non-basic

813 * with upper bound active (GLP_NU, lambda~[q] <= 0), and this upper

814 * bound be implied one u'[q]. As in the previous case this then means

815 * that from the original problem's standpoint actually the original

816 * column upper bound u[q] is inactive, and active is that row bound

817 * L[p] or U[p] that defines the implied bound u'[q] (3). In this case

818 * in solution to the original problem column q is assigned status

819 * GLP_BS, row p is assigned status GLP_NU (if a[p,q] > 0) or GLP_NL

820 * (if a[p,q] < 0), and its multiplier is computed with formula (9).

821 *

822 * Strengthening bounds of column q according to (5) and (6) may make

823 * it fixed. Thus, if in solution to the transformed problem column q is

824 * non-basic and fixed (GLP_NS), we can suppose that if lambda~[q] > 0,

825 * column q has active lower bound (GLP_NL), and if lambda~[q] < 0,

826 * column q has active upper bound (GLP_NU), reducing this case to two

827 * previous ones. If, however, lambda~[q] is close to zero or

828 * corresponding bound of row p does not exist (this may happen if

829 * lambda~[q] has wrong sign due to round-off errors, in which case it

830 * is expected to be close to zero, since solution is assumed to be dual

831 * feasible), column q can be assigned status GLP_BS (basic), and row p

832 * can be made active on its existing bound. In the latter case row

833 * multiplier pi[p] computed with formula (9) will be also close to

834 * zero, and dual feasibility will be kept.

835 *

836 * In all other cases, namely, if in solution to the transformed

837 * problem column q is basic (GLP_BS), or non-basic with original lower

838 * bound l[q] active (GLP_NL), or non-basic with original upper bound

839 * u[q] active (GLP_NU), constraint (1) is inactive. So in solution to

840 * the original problem status of column q remains unchanged, row p is

841 * assigned status GLP_BS, and its multiplier pi[p] is assigned zero

842 * value.

843 *

844 * RECOVERING INTERIOR-POINT SOLUTION

845 *

846 * First, value of multiplier for column q in solution to the original

847 * problem is computed with formula (8). If lambda~[q] > 0 and column q

848 * has implied lower bound, or if lambda~[q] < 0 and column q has

849 * implied upper bound, this means that from the original problem's

850 * standpoint actually row p has corresponding active bound, in which

851 * case its multiplier pi[p] is computed with formula (9). In other

852 * cases, when the sign of lambda~[q] corresponds to original bound of

853 * column q, or when lambda~[q] =~ 0, value of row multiplier pi[p] is

854 * assigned zero value.

855 *

856 * RECOVERING MIP SOLUTION

857 *

858 * None needed. */

859

860 struct ineq_singlet

861 { /* row singleton (inequality constraint) */

862 int p;

863 /* row reference number */

864 int q;

865 /* column reference number */

866 double apq;

867 /* constraint coefficient a[p,q] */

868 double c;

869 /* objective coefficient at x[q] */

870 double lb;

871 /* row lower bound */

872 double ub;

873 /* row upper bound */

874 char lb_changed;

875 /* this flag is set if column lower bound was changed */

876 char ub_changed;

877 /* this flag is set if column upper bound was changed */

878 NPPLFE *ptr;

879 /* list of non-zero coefficients a[i,q], i != p */

880 };

881

882 static int rcv_ineq_singlet(NPP *npp, void *info);

883

884 int npp_ineq_singlet(NPP *npp, NPPROW *p)

885 { /* process row singleton (inequality constraint) */

886 struct ineq_singlet *info;

887 NPPCOL *q;

888 NPPAIJ *apq, *aij;

889 NPPLFE *lfe;

890 int lb_changed, ub_changed;

891 double ll, uu;

892 /* the row must be singleton inequality constraint */

893 xassert(p->lb != -DBL_MAX || p->ub != +DBL_MAX);

894 xassert(p->lb < p->ub);

895 xassert(p->ptr != NULL && p->ptr->r_next == NULL);

896 /* compute implied column bounds */

897 apq = p->ptr;

898 q = apq->col;

899 xassert(q->lb < q->ub);

900 if (apq->val > 0.0)

901 { ll = (p->lb == -DBL_MAX ? -DBL_MAX : p->lb / apq->val);

902 uu = (p->ub == +DBL_MAX ? +DBL_MAX : p->ub / apq->val);

903 }

904 else

905 { ll = (p->ub == +DBL_MAX ? -DBL_MAX : p->ub / apq->val);

906 uu = (p->lb == -DBL_MAX ? +DBL_MAX : p->lb / apq->val);

907 }

908 /* process implied column lower bound */

909 if (ll == -DBL_MAX)

910 lb_changed = 0;

911 else

912 { lb_changed = npp_implied_lower(npp, q, ll);

913 xassert(0 <= lb_changed && lb_changed <= 4);

914 if (lb_changed == 4) return 4; /* infeasible */

915 }

916 /* process implied column upper bound */

917 if (uu == +DBL_MAX)

918 ub_changed = 0;

919 else if (lb_changed == 3)

920 { /* column was fixed on its upper bound due to l'[q] = u[q] */

921 /* note that L[p] < U[p], so l'[q] = u[q] < u'[q] */

922 ub_changed = 0;

923 }

924 else

925 { ub_changed = npp_implied_upper(npp, q, uu);

926 xassert(0 <= ub_changed && ub_changed <= 4);

927 if (ub_changed == 4) return 4; /* infeasible */

928 }

929 /* if neither lower nor upper column bound was changed, the row

930 is originally redundant and can be replaced by free row */

931 if (!lb_changed && !ub_changed)

932 { p->lb = -DBL_MAX, p->ub = +DBL_MAX;

933 npp_free_row(npp, p);

934 return 0;

935 }

936 /* create transformation stack entry */

937 info = npp_push_tse(npp,

938 rcv_ineq_singlet, sizeof(struct ineq_singlet));

939 info->p = p->i;

940 info->q = q->j;

941 info->apq = apq->val;

942 info->c = q->coef;

943 info->lb = p->lb;

944 info->ub = p->ub;

945 info->lb_changed = (char)lb_changed;

946 info->ub_changed = (char)ub_changed;

947 info->ptr = NULL;

948 /* save column coefficients a[i,q], i != p (not needed for MIP

949 solution) */

950 if (npp->sol != GLP_MIP)

951 { for (aij = q->ptr; aij != NULL; aij = aij->c_next)

952 { if (aij == apq) continue; /* skip a[p,q] */

953 lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE));

954 lfe->ref = aij->row->i;

955 lfe->val = aij->val;

956 lfe->next = info->ptr;

957 info->ptr = lfe;

958 }

959 }

960 /* remove the row from the problem */

961 npp_del_row(npp, p);

962 return lb_changed >= ub_changed ? lb_changed : ub_changed;

963 }

964

965 static int rcv_ineq_singlet(NPP *npp, void *_info)

966 { /* recover row singleton (inequality constraint) */

967 struct ineq_singlet *info = _info;

968 NPPLFE *lfe;

969 double lambda;

970 if (npp->sol == GLP_MIP) goto done;

971 /* compute lambda~[q] in solution to the transformed problem

972 with formula (8) */

973 lambda = info->c;

974 for (lfe = info->ptr; lfe != NULL; lfe = lfe->next)

975 lambda -= lfe->val * npp->r_pi[lfe->ref];

976 if (npp->sol == GLP_SOL)

977 { /* recover basic solution */

978 if (npp->c_stat[info->q] == GLP_BS)

979 { /* column q is basic, so row p is inactive */

980 npp->r_stat[info->p] = GLP_BS;

981 npp->r_pi[info->p] = 0.0;

982 }

983 else if (npp->c_stat[info->q] == GLP_NL)

984 nl: { /* column q is non-basic with lower bound active */

985 if (info->lb_changed)

986 { /* it is implied bound, so actually row p is active

987 while column q is basic */

988 npp->r_stat[info->p] =

989 (char)(info->apq > 0.0 ? GLP_NL : GLP_NU);

990 npp->c_stat[info->q] = GLP_BS;

991 npp->r_pi[info->p] = lambda / info->apq;

992 }

993 else

994 { /* it is original bound, so row p is inactive */

995 npp->r_stat[info->p] = GLP_BS;

996 npp->r_pi[info->p] = 0.0;

997 }

998 }

999 else if (npp->c_stat[info->q] == GLP_NU)

1000 nu: { /* column q is non-basic with upper bound active */

1001 if (info->ub_changed)

1002 { /* it is implied bound, so actually row p is active

1003 while column q is basic */

1004 npp->r_stat[info->p] =

1005 (char)(info->apq > 0.0 ? GLP_NU : GLP_NL);

1006 npp->c_stat[info->q] = GLP_BS;

1007 npp->r_pi[info->p] = lambda / info->apq;

1008 }

1009 else

1010 { /* it is original bound, so row p is inactive */

1011 npp->r_stat[info->p] = GLP_BS;

1012 npp->r_pi[info->p] = 0.0;

1013 }

1014 }

1015 else if (npp->c_stat[info->q] == GLP_NS)

1016 { /* column q is non-basic and fixed; note, however, that in

1017 in the original problem it is non-fixed */

1018 if (lambda > +1e-7)

1019 { if (info->apq > 0.0 && info->lb != -DBL_MAX ||

1020 info->apq < 0.0 && info->ub != +DBL_MAX ||

1021 !info->lb_changed)

1022 { /* either corresponding bound of row p exists or

1023 column q remains non-basic with its original lower

1024 bound active */

1025 npp->c_stat[info->q] = GLP_NL;

1026 goto nl;

1027 }

1028 }

1029 if (lambda < -1e-7)

1030 { if (info->apq > 0.0 && info->ub != +DBL_MAX ||

1031 info->apq < 0.0 && info->lb != -DBL_MAX ||

1032 !info->ub_changed)

1033 { /* either corresponding bound of row p exists or

1034 column q remains non-basic with its original upper

1035 bound active */

1036 npp->c_stat[info->q] = GLP_NU;

1037 goto nu;

1038 }

1039 }

1040 /* either lambda~[q] is close to zero, or corresponding

1041 bound of row p does not exist, because lambda~[q] has

1042 wrong sign due to round-off errors; in the latter case

1043 lambda~[q] is also assumed to be close to zero; so, we

1044 can make row p active on its existing bound and column q

1045 basic; pi[p] will have wrong sign, but it also will be

1046 close to zero (rarus casus of dual degeneracy) */

1047 if (info->lb != -DBL_MAX && info->ub == +DBL_MAX)

1048 { /* row lower bound exists, but upper bound doesn't */

1049 npp->r_stat[info->p] = GLP_NL;

1050 }

1051 else if (info->lb == -DBL_MAX && info->ub != +DBL_MAX)

1052 { /* row upper bound exists, but lower bound doesn't */

1053 npp->r_stat[info->p] = GLP_NU;

1054 }

1055 else if (info->lb != -DBL_MAX && info->ub != +DBL_MAX)

1056 { /* both row lower and upper bounds exist */

1057 /* to choose proper active row bound we should not use

1058 lambda~[q], because its value being close to zero is

1059 unreliable; so we choose that bound which provides

1060 primal feasibility for original constraint (1) */

1061 if (info->apq * npp->c_value[info->q] <=

1062 0.5 * (info->lb + info->ub))

1063 npp->r_stat[info->p] = GLP_NL;

1064 else

1065 npp->r_stat[info->p] = GLP_NU;

1066 }

1067 else

1068 { npp_error();

1069 return 1;

1070 }

1071 npp->c_stat[info->q] = GLP_BS;

1072 npp->r_pi[info->p] = lambda / info->apq;

1073 }

1074 else

1075 { npp_error();

1076 return 1;

1077 }

1078 }

1079 if (npp->sol == GLP_IPT)

1080 { /* recover interior-point solution */

1081 if (lambda > +DBL_EPSILON && info->lb_changed ||

1082 lambda < -DBL_EPSILON && info->ub_changed)

1083 { /* actually row p has corresponding active bound */

1084 npp->r_pi[info->p] = lambda / info->apq;

1085 }

1086 else

1087 { /* either bounds of column q are both inactive or its

1088 original bound is active */

1089 npp->r_pi[info->p] = 0.0;

1090 }

1091 }

1092 done: return 0;

1093 }

1094

1095 /***********************************************************************

1096 * NAME

1097 *

1098 * npp_implied_slack - process column singleton (implied slack variable)

1099 *

1100 * SYNOPSIS

1101 *

1102 * #include "glpnpp.h"

1103 * void npp_implied_slack(NPP *npp, NPPCOL *q);

1104 *

1105 * DESCRIPTION

1106 *

1107 * The routine npp_implied_slack processes column q:

1108 *

1109 * l[q] <= x[q] <= u[q], (1)

1110 *

1111 * where l[q] < u[q], having the only non-zero coefficient in row p,

1112 * which is equality constraint:

1113 *

1114 * sum a[p,j] x[j] + a[p,q] x[q] = b. (2)

1115 * j!=q

1116 *

1117 * PROBLEM TRANSFORMATION

1118 *

1119 * (If x[q] is integral, this transformation must not be used.)

1120 *

1121 * The term a[p,q] x[q] in constraint (2) can be considered as a slack

1122 * variable that allows to carry bounds of column q over row p and then

1123 * remove column q from the problem.

1124 *

1125 * Constraint (2) can be written as follows:

1126 *

1127 * sum a[p,j] x[j] = b - a[p,q] x[q]. (3)

1128 * j!=q

1129 *

1130 * According to (1) constraint (3) is equivalent to the following

1131 * inequality constraint:

1132 *

1133 * L[p] <= sum a[p,j] x[j] <= U[p], (4)

1134 * j!=q

1135 *

1136 * where

1137 *

1138 * ( b - a[p,q] u[q], if a[p,q] > 0

1139 * L[p] = < (5)

1140 * ( b - a[p,q] l[q], if a[p,q] < 0

1141 *

1142 * ( b - a[p,q] l[q], if a[p,q] > 0

1143 * U[p] = < (6)

1144 * ( b - a[p,q] u[q], if a[p,q] < 0

1145 *

1146 * From (2) it follows that:

1147 *

1148 * 1

1149 * x[q] = ------ (b - sum a[p,j] x[j]). (7)

1150 * a[p,q] j!=q

1151 *

1152 * In order to eliminate x[q] from the objective row we substitute it

1153 * from (6) to that row:

1154 *

1155 * z = sum c[j] x[j] + c[q] x[q] + c[0] =

1156 * j!=q

1157 * 1

1158 * = sum c[j] x[j] + c[q] [------ (b - sum a[p,j] x[j])] + c0 =

1159 * j!=q a[p,q] j!=q

1160 *

1161 * = sum c~[j] x[j] + c~[0],

1162 * j!=q

1163 * a[p,j] b

1164 * c~[j] = c[j] - c[q] ------, c~0 = c0 - c[q] ------ (8)

1165 * a[p,q] a[p,q]

1166 *

1167 * are values of objective coefficients and constant term, resp., in

1168 * the transformed problem.

1169 *

1170 * Note that column q is column singleton, so in the dual system of the

1171 * original problem it corresponds to the following row singleton:

1172 *

1173 * a[p,q] pi[p] + lambda[q] = c[q]. (9)

1174 *

1175 * In the transformed problem row (9) would be the following:

1176 *

1177 * a[p,q] pi~[p] + lambda[q] = c~[q] = 0. (10)

1178 *

1179 * Subtracting (10) from (9) we have:

1180 *

1181 * a[p,q] (pi[p] - pi~[p]) = c[q]

1182 *

1183 * that gives the following formula to compute multiplier for row p in

1184 * solution to the original problem using its value in solution to the

1185 * transformed problem:

1186 *

1187 * pi[p] = pi~[p] + c[q] / a[p,q]. (11)

1188 *

1189 * RECOVERING BASIC SOLUTION

1190 *

1191 * Status of column q in solution to the original problem is defined

1192 * by status of row p in solution to the transformed problem and the

1193 * sign of coefficient a[p,q] in the original inequality constraint (2)

1194 * as follows:

1195 *

1196 * +-----------------------+---------+--------------------+

1197 * | Status of row p | Sign of | Status of column q |

1198 * | (transformed problem) | a[p,q] | (original problem) |

1199 * +-----------------------+---------+--------------------+

1200 * | GLP_BS | + / - | GLP_BS |

1201 * | GLP_NL | + | GLP_NU |

1202 * | GLP_NL | - | GLP_NL |

1203 * | GLP_NU | + | GLP_NL |

1204 * | GLP_NU | - | GLP_NU |

1205 * | GLP_NF | + / - | GLP_NF |

1206 * +-----------------------+---------+--------------------+

1207 *

1208 * Value of column q is computed with formula (7). Since originally row

1209 * p is equality constraint, its status is assigned GLP_NS, and value of

1210 * its multiplier pi[p] is computed with formula (11).

1211 *

1212 * RECOVERING INTERIOR-POINT SOLUTION

1213 *

1214 * Value of column q is computed with formula (7). Row multiplier value

1215 * pi[p] is computed with formula (11).

1216 *

1217 * RECOVERING MIP SOLUTION

1218 *

1219 * Value of column q is computed with formula (7). */

1220

1221 struct implied_slack

1222 { /* column singleton (implied slack variable) */

1223 int p;

1224 /* row reference number */

1225 int q;

1226 /* column reference number */

1227 double apq;

1228 /* constraint coefficient a[p,q] */

1229 double b;

1230 /* right-hand side of original equality constraint */

1231 double c;

1232 /* original objective coefficient at x[q] */

1233 NPPLFE *ptr;

1234 /* list of non-zero coefficients a[p,j], j != q */

1235 };

1236

1237 static int rcv_implied_slack(NPP *npp, void *info);

1238

1239 void npp_implied_slack(NPP *npp, NPPCOL *q)

1240 { /* process column singleton (implied slack variable) */

1241 struct implied_slack *info;

1242 NPPROW *p;

1243 NPPAIJ *aij;

1244 NPPLFE *lfe;

1245 /* the column must be non-integral non-fixed singleton */

1246 xassert(!q->is_int);

1247 xassert(q->lb < q->ub);

1248 xassert(q->ptr != NULL && q->ptr->c_next == NULL);

1249 /* corresponding row must be equality constraint */

1250 aij = q->ptr;

1251 p = aij->row;

1252 xassert(p->lb == p->ub);

1253 /* create transformation stack entry */

1254 info = npp_push_tse(npp,

1255 rcv_implied_slack, sizeof(struct implied_slack));

1256 info->p = p->i;

1257 info->q = q->j;

1258 info->apq = aij->val;

1259 info->b = p->lb;

1260 info->c = q->coef;

1261 info->ptr = NULL;

1262 /* save row coefficients a[p,j], j != q, and substitute x[q]

1263 into the objective row */

1264 for (aij = p->ptr; aij != NULL; aij = aij->r_next)

1265 { if (aij->col == q) continue; /* skip a[p,q] */

1266 lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE));

1267 lfe->ref = aij->col->j;

1268 lfe->val = aij->val;

1269 lfe->next = info->ptr;

1270 info->ptr = lfe;

1271 aij->col->coef -= info->c * (aij->val / info->apq);

1272 }

1273 npp->c0 += info->c * (info->b / info->apq);

1274 /* compute new row bounds */

1275 if (info->apq > 0.0)

1276 { p->lb = (q->ub == +DBL_MAX ?

1277 -DBL_MAX : info->b - info->apq * q->ub);

1278 p->ub = (q->lb == -DBL_MAX ?

1279 +DBL_MAX : info->b - info->apq * q->lb);

1280 }

1281 else

1282 { p->lb = (q->lb == -DBL_MAX ?

1283 -DBL_MAX : info->b - info->apq * q->lb);

1284 p->ub = (q->ub == +DBL_MAX ?

1285 +DBL_MAX : info->b - info->apq * q->ub);

1286 }

1287 /* remove the column from the problem */

1288 npp_del_col(npp, q);

1289 return;

1290 }

1291

1292 static int rcv_implied_slack(NPP *npp, void *_info)

1293 { /* recover column singleton (implied slack variable) */

1294 struct implied_slack *info = _info;

1295 NPPLFE *lfe;

1296 double temp;

1297 if (npp->sol == GLP_SOL)

1298 { /* assign statuses to row p and column q */

1299 if (npp->r_stat[info->p] == GLP_BS ||

1300 npp->r_stat[info->p] == GLP_NF)

1301 npp->c_stat[info->q] = npp->r_stat[info->p];

1302 else if (npp->r_stat[info->p] == GLP_NL)

1303 npp->c_stat[info->q] =

1304 (char)(info->apq > 0.0 ? GLP_NU : GLP_NL);

1305 else if (npp->r_stat[info->p] == GLP_NU)

1306 npp->c_stat[info->q] =

1307 (char)(info->apq > 0.0 ? GLP_NL : GLP_NU);

1308 else

1309 { npp_error();

1310 return 1;

1311 }

1312 npp->r_stat[info->p] = GLP_NS;

1313 }

1314 if (npp->sol != GLP_MIP)

1315 { /* compute multiplier for row p */

1316 npp->r_pi[info->p] += info->c / info->apq;

1317 }

1318 /* compute value of column q */

1319 temp = info->b;

1320 for (lfe = info->ptr; lfe != NULL; lfe = lfe->next)

1321 temp -= lfe->val * npp->c_value[lfe->ref];

1322 npp->c_value[info->q] = temp / info->apq;

1323 return 0;

1324 }

1325

1326 /***********************************************************************

1327 * NAME

1328 *

1329 * npp_implied_free - process column singleton (implied free variable)

1330 *

1331 * SYNOPSIS

1332 *

1333 * #include "glpnpp.h"

1334 * int npp_implied_free(NPP *npp, NPPCOL *q);

1335 *

1336 * DESCRIPTION

1337 *

1338 * The routine npp_implied_free processes column q:

1339 *

1340 * l[q] <= x[q] <= u[q], (1)

1341 *

1342 * having non-zero coefficient in the only row p, which is inequality

1343 * constraint:

1344 *

1345 * L[p] <= sum a[p,j] x[j] + a[p,q] x[q] <= U[p], (2)

1346 * j!=q

1347 *

1348 * where l[q] < u[q], L[p] < U[p], L[p] > -oo and/or U[p] < +oo.

1349 *

1350 * RETURNS

1351 *

1352 * 0 - success;

1353 *

1354 * 1 - column lower and/or upper bound(s) can be active;

1355 *

1356 * 2 - problem has no dual feasible solution.

1357 *

1358 * PROBLEM TRANSFORMATION

1359 *

1360 * Constraint (2) can be written as follows:

1361 *

1362 * L[p] - sum a[p,j] x[j] <= a[p,q] x[q] <= U[p] - sum a[p,j] x[j],

1363 * j!=q j!=q

1364 *

1365 * from which it follows that:

1366 *

1367 * alfa <= a[p,q] x[q] <= beta, (3)

1368 *

1369 * where

1370 *

1371 * alfa = inf(L[p] - sum a[p,j] x[j]) =

1372 * j!=q

1373 *

1374 * = L[p] - sup sum a[p,j] x[j] = (4)

1375 * j!=q

1376 *

1377 * = L[p] - sum a[p,j] u[j] - sum a[p,j] l[j],

1378 * j in Jp j in Jn

1379 *

1380 * beta = sup(L[p] - sum a[p,j] x[j]) =

1381 * j!=q

1382 *

1383 * = L[p] - inf sum a[p,j] x[j] = (5)

1384 * j!=q

1385 *

1386 * = L[p] - sum a[p,j] l[j] - sum a[p,j] u[j],

1387 * j in Jp j in Jn

1388 *

1389 * Jp = {j != q: a[p,j] > 0}, Jn = {j != q: a[p,j] < 0}. (6)

1390 *

1391 * Inequality (3) defines implied bounds of variable x[q]:

1392 *

1393 * l'[q] <= x[q] <= u'[q], (7)

1394 *

1395 * where

1396 *

1397 * ( alfa / a[p,q], if a[p,q] > 0

1398 * l'[q] = < (8a)

1399 * ( beta / a[p,q], if a[p,q] < 0

1400 *

1401 * ( beta / a[p,q], if a[p,q] > 0

1402 * u'[q] = < (8b)

1403 * ( alfa / a[p,q], if a[p,q] < 0

1404 *

1405 * Thus, if l'[q] > l[q] - eps and u'[q] < u[q] + eps, where eps is

1406 * an absolute tolerance for column value, column bounds (1) cannot be

1407 * active, in which case column q can be replaced by equivalent free

1408 * (unbounded) column.

1409 *

1410 * Note that column q is column singleton, so in the dual system of the

1411 * original problem it corresponds to the following row singleton:

1412 *

1413 * a[p,q] pi[p] + lambda[q] = c[q], (9)

1414 *

1415 * from which it follows that:

1416 *

1417 * pi[p] = (c[q] - lambda[q]) / a[p,q]. (10)

1418 *

1419 * Let x[q] be implied free (unbounded) variable. Then column q can be

1420 * only basic, so its multiplier lambda[q] is equal to zero, and from

1421 * (10) we have:

1422 *

1423 * pi[p] = c[q] / a[p,q]. (11)

1424 *

1425 * There are possible three cases:

1426 *

1427 * 1) pi[p] < -eps, where eps is an absolute tolerance for row

1428 * multiplier. In this case, to provide dual feasibility of the

1429 * original problem, row p must be active on its lower bound, and

1430 * if its lower bound does not exist (L[p] = -oo), the problem has

1431 * no dual feasible solution;

1432 *

1433 * 2) pi[p] > +eps. In this case row p must be active on its upper

1434 * bound, and if its upper bound does not exist (U[p] = +oo), the

1435 * problem has no dual feasible solution;

1436 *

1437 * 3) -eps <= pi[p] <= +eps. In this case any (either lower or upper)

1438 * bound of row p can be active, because this does not affect dual

1439 * feasibility.

1440 *

1441 * Thus, in all three cases original inequality constraint (2) can be

1442 * replaced by equality constraint, where the right-hand side is either

1443 * lower or upper bound of row p, and bounds of column q can be removed

1444 * that makes it free (unbounded). (May note that this transformation

1445 * can be followed by transformation "Column singleton (implied slack

1446 * variable)" performed by the routine npp_implied_slack.)

1447 *

1448 * RECOVERING BASIC SOLUTION

1449 *

1450 * Status of row p in solution to the original problem is determined

1451 * by its status in solution to the transformed problem and its bound,

1452 * which was choosen to be active:

1453 *

1454 * +-----------------------+--------+--------------------+

1455 * | Status of row p | Active | Status of row p |

1456 * | (transformed problem) | bound | (original problem) |

1457 * +-----------------------+--------+--------------------+

1458 * | GLP_BS | L[p] | GLP_BS |

1459 * | GLP_BS | U[p] | GLP_BS |

1460 * | GLP_NS | L[p] | GLP_NL |

1461 * | GLP_NS | U[p] | GLP_NU |

1462 * +-----------------------+--------+--------------------+

1463 *

1464 * Value of row multiplier pi[p] (as well as value of column q) in

1465 * solution to the original problem is the same as in solution to the

1466 * transformed problem.

1467 *

1468 * RECOVERING INTERIOR-POINT SOLUTION

1469 *

1470 * Value of row multiplier pi[p] in solution to the original problem is

1471 * the same as in solution to the transformed problem.

1472 *

1473 * RECOVERING MIP SOLUTION

1474 *

1475 * None needed. */

1476

1477 struct implied_free

1478 { /* column singleton (implied free variable) */

1479 int p;

1480 /* row reference number */

1481 char stat;

1482 /* row status:

1483 GLP_NL - active constraint on lower bound

1484 GLP_NU - active constraint on upper bound */

1485 };

1486

1487 static int rcv_implied_free(NPP *npp, void *info);

1488

1489 int npp_implied_free(NPP *npp, NPPCOL *q)

1490 { /* process column singleton (implied free variable) */

1491 struct implied_free *info;

1492 NPPROW *p;

1493 NPPAIJ *apq, *aij;

1494 double alfa, beta, l, u, pi, eps;

1495 /* the column must be non-fixed singleton */

1496 xassert(q->lb < q->ub);

1497 xassert(q->ptr != NULL && q->ptr->c_next == NULL);

1498 /* corresponding row must be inequality constraint */

1499 apq = q->ptr;

1500 p = apq->row;

1501 xassert(p->lb != -DBL_MAX || p->ub != +DBL_MAX);

1502 xassert(p->lb < p->ub);

1503 /* compute alfa */

1504 alfa = p->lb;

1505 if (alfa != -DBL_MAX)

1506 { for (aij = p->ptr; aij != NULL; aij = aij->r_next)

1507 { if (aij == apq) continue; /* skip a[p,q] */

1508 if (aij->val > 0.0)

1509 { if (aij->col->ub == +DBL_MAX)

1510 { alfa = -DBL_MAX;

1511 break;

1512 }

1513 alfa -= aij->val * aij->col->ub;

1514 }

1515 else /* < 0.0 */

1516 { if (aij->col->lb == -DBL_MAX)

1517 { alfa = -DBL_MAX;

1518 break;

1519 }

1520 alfa -= aij->val * aij->col->lb;

1521 }

1522 }

1523 }

1524 /* compute beta */

1525 beta = p->ub;

1526 if (beta != +DBL_MAX)

1527 { for (aij = p->ptr; aij != NULL; aij = aij->r_next)

1528 { if (aij == apq) continue; /* skip a[p,q] */

1529 if (aij->val > 0.0)

1530 { if (aij->col->lb == -DBL_MAX)

1531 { beta = +DBL_MAX;

1532 break;

1533 }

1534 beta -= aij->val * aij->col->lb;

1535 }

1536 else /* < 0.0 */

1537 { if (aij->col->ub == +DBL_MAX)

1538 { beta = +DBL_MAX;

1539 break;

1540 }

1541 beta -= aij->val * aij->col->ub;

1542 }

1543 }

1544 }

1545 /* compute implied column lower bound l'[q] */

1546 if (apq->val > 0.0)

1547 l = (alfa == -DBL_MAX ? -DBL_MAX : alfa / apq->val);

1548 else /* < 0.0 */

1549 l = (beta == +DBL_MAX ? -DBL_MAX : beta / apq->val);

1550 /* compute implied column upper bound u'[q] */

1551 if (apq->val > 0.0)

1552 u = (beta == +DBL_MAX ? +DBL_MAX : beta / apq->val);

1553 else

1554 u = (alfa == -DBL_MAX ? +DBL_MAX : alfa / apq->val);

1555 /* check if column lower bound l[q] can be active */

1556 if (q->lb != -DBL_MAX)

1557 { eps = 1e-9 + 1e-12 * fabs(q->lb);

1558 if (l < q->lb - eps) return 1; /* yes, it can */

1559 }

1560 /* check if column upper bound u[q] can be active */

1561 if (q->ub != +DBL_MAX)

1562 { eps = 1e-9 + 1e-12 * fabs(q->ub);

1563 if (u > q->ub + eps) return 1; /* yes, it can */

1564 }

1565 /* okay; make column q free (unbounded) */

1566 q->lb = -DBL_MAX, q->ub = +DBL_MAX;

1567 /* create transformation stack entry */

1568 info = npp_push_tse(npp,

1569 rcv_implied_free, sizeof(struct implied_free));

1570 info->p = p->i;

1571 info->stat = -1;

1572 /* compute row multiplier pi[p] */

1573 pi = q->coef / apq->val;

1574 /* check dual feasibility for row p */

1575 if (pi > +DBL_EPSILON)

1576 { /* lower bound L[p] must be active */

1577 if (p->lb != -DBL_MAX)

1578 nl: { info->stat = GLP_NL;

1579 p->ub = p->lb;

1580 }

1581 else

1582 { if (pi > +1e-5) return 2; /* dual infeasibility */

1583 /* take a chance on U[p] */

1584 xassert(p->ub != +DBL_MAX);

1585 goto nu;

1586 }

1587 }

1588 else if (pi < -DBL_EPSILON)

1589 { /* upper bound U[p] must be active */

1590 if (p->ub != +DBL_MAX)

1591 nu: { info->stat = GLP_NU;

1592 p->lb = p->ub;

1593 }

1594 else

1595 { if (pi < -1e-5) return 2; /* dual infeasibility */

1596 /* take a chance on L[p] */

1597 xassert(p->lb != -DBL_MAX);

1598 goto nl;

1599 }

1600 }

1601 else

1602 { /* any bound (either L[p] or U[p]) can be made active */

1603 if (p->ub == +DBL_MAX)

1604 { xassert(p->lb != -DBL_MAX);

1605 goto nl;

1606 }

1607 if (p->lb == -DBL_MAX)

1608 { xassert(p->ub != +DBL_MAX);

1609 goto nu;

1610 }

1611 if (fabs(p->lb) <= fabs(p->ub)) goto nl; else goto nu;

1612 }

1613 return 0;

1614 }

1615

1616 static int rcv_implied_free(NPP *npp, void *_info)

1617 { /* recover column singleton (implied free variable) */

1618 struct implied_free *info = _info;

1619 if (npp->sol == GLP_SOL)

1620 { if (npp->r_stat[info->p] == GLP_BS)

1621 npp->r_stat[info->p] = GLP_BS;

1622 else if (npp->r_stat[info->p] == GLP_NS)

1623 { xassert(info->stat == GLP_NL || info->stat == GLP_NU);

1624 npp->r_stat[info->p] = info->stat;

1625 }

1626 else

1627 { npp_error();

1628 return 1;

1629 }

1630 }

1631 return 0;

1632 }

1633

1634 /***********************************************************************

1635 * NAME

1636 *

1637 * npp_eq_doublet - process row doubleton (equality constraint)

1638 *

1639 * SYNOPSIS

1640 *

1641 * #include "glpnpp.h"

1642 * NPPCOL *npp_eq_doublet(NPP *npp, NPPROW *p);

1643 *

1644 * DESCRIPTION

1645 *

1646 * The routine npp_eq_doublet processes row p, which is equality

1647 * constraint having exactly two non-zero coefficients:

1648 *

1649 * a[p,q] x[q] + a[p,r] x[r] = b. (1)

1650 *

1651 * As the result of processing one of columns q or r is eliminated from

1652 * all other rows and, thus, becomes column singleton of type "implied

1653 * slack variable". Row p is not changed and along with column q and r

1654 * remains in the problem.

1655 *

1656 * RETURNS

1657 *

1658 * The routine npp_eq_doublet returns pointer to the descriptor of that

1659 * column q or r which has been eliminated. If, due to some reason, the

1660 * elimination was not performed, the routine returns NULL.

1661 *

1662 * PROBLEM TRANSFORMATION

1663 *

1664 * First, we decide which column q or r will be eliminated. Let it be

1665 * column q. Consider i-th constraint row, where column q has non-zero

1666 * coefficient a[i,q] != 0:

1667 *

1668 * L[i] <= sum a[i,j] x[j] <= U[i]. (2)

1669 * j

1670 *

1671 * In order to eliminate column q from row (2) we subtract from it row

1672 * (1) multiplied by gamma[i] = a[i,q] / a[p,q], i.e. we replace in the

1673 * transformed problem row (2) by its linear combination with row (1).

1674 * This transformation changes only coefficients in columns q and r,

1675 * and bounds of row i as follows:

1676 *

1677 * a~[i,q] = a[i,q] - gamma[i] a[p,q] = 0, (3)

1678 *

1679 * a~[i,r] = a[i,r] - gamma[i] a[p,r], (4)

1680 *

1681 * L~[i] = L[i] - gamma[i] b, (5)

1682 *

1683 * U~[i] = U[i] - gamma[i] b. (6)

1684 *

1685 * RECOVERING BASIC SOLUTION

1686 *

1687 * The transformation of the primal system of the original problem:

1688 *

1689 * L <= A x <= U (7)

1690 *

1691 * is equivalent to multiplying from the left a transformation matrix F

1692 * by components of this primal system, which in the transformed problem

1693 * becomes the following:

1694 *

1695 * F L <= F A x <= F U ==> L~ <= A~x <= U~. (8)

1696 *

1697 * The matrix F has the following structure:

1698 *

1699 * ( 1 -gamma[1] )

1700 * ( )

1701 * ( 1 -gamma[2] )

1702 * ( )

1703 * ( ... ... )

1704 * ( )

1705 * F = ( 1 -gamma[p-1] ) (9)

1706 * ( )

1707 * ( 1 )

1708 * ( )

1709 * ( -gamma[p+1] 1 )

1710 * ( )

1711 * ( ... ... )

1712 *

1713 * where its column containing elements -gamma[i] corresponds to row p

1714 * of the primal system.

1715 *

1716 * From (8) it follows that the dual system of the original problem:

1717 *

1718 * A'pi + lambda = c, (10)

1719 *

1720 * in the transformed problem becomes the following:

1721 *

1722 * A'F'inv(F')pi + lambda = c ==> (A~)'pi~ + lambda = c, (11)

1723 *

1724 * where:

1725 *

1726 * pi~ = inv(F')pi (12)

1727 *

1728 * is the vector of row multipliers in the transformed problem. Thus:

1729 *

1730 * pi = F'pi~. (13)

1731 *

1732 * Therefore, as it follows from (13), value of multiplier for row p in

1733 * solution to the original problem can be computed as follows:

1734 *

1735 * pi[p] = pi~[p] - sum gamma[i] pi~[i], (14)

1736 * i

1737 *

1738 * where pi~[i] = pi[i] is multiplier for row i (i != p).

1739 *

1740 * Note that the statuses of all rows and columns are not changed.

1741 *

1742 * RECOVERING INTERIOR-POINT SOLUTION

1743 *

1744 * Multiplier for row p in solution to the original problem is computed

1745 * with formula (14).

1746 *

1747 * RECOVERING MIP SOLUTION

1748 *

1749 * None needed. */

1750

1751 struct eq_doublet

1752 { /* row doubleton (equality constraint) */

1753 int p;

1754 /* row reference number */

1755 double apq;

1756 /* constraint coefficient a[p,q] */

1757 NPPLFE *ptr;

1758 /* list of non-zero coefficients a[i,q], i != p */

1759 };

1760

1761 static int rcv_eq_doublet(NPP *npp, void *info);

1762

1763 NPPCOL *npp_eq_doublet(NPP *npp, NPPROW *p)

1764 { /* process row doubleton (equality constraint) */

1765 struct eq_doublet *info;

1766 NPPROW *i;

1767 NPPCOL *q, *r;

1768 NPPAIJ *apq, *apr, *aiq, *air, *next;

1769 NPPLFE *lfe;

1770 double gamma;

1771 /* the row must be doubleton equality constraint */

1772 xassert(p->lb == p->ub);

1773 xassert(p->ptr != NULL && p->ptr->r_next != NULL &&

1774 p->ptr->r_next->r_next == NULL);

1775 /* choose column to be eliminated */

1776 { NPPAIJ *a1, *a2;

1777 a1 = p->ptr, a2 = a1->r_next;

1778 if (fabs(a2->val) < 0.001 * fabs(a1->val))

1779 { /* only first column can be eliminated, because second one

1780 has too small constraint coefficient */

1781 apq = a1, apr = a2;

1782 }

1783 else if (fabs(a1->val) < 0.001 * fabs(a2->val))

1784 { /* only second column can be eliminated, because first one

1785 has too small constraint coefficient */

1786 apq = a2, apr = a1;

1787 }

1788 else

1789 { /* both columns are appropriate; choose that one which is

1790 shorter to minimize fill-in */

1791 if (npp_col_nnz(npp, a1->col) <= npp_col_nnz(npp, a2->col))

1792 { /* first column is shorter */

1793 apq = a1, apr = a2;

1794 }

1795 else

1796 { /* second column is shorter */

1797 apq = a2, apr = a1;

1798 }

1799 }

1800 }

1801 /* now columns q and r have been chosen */

1802 q = apq->col, r = apr->col;

1803 /* create transformation stack entry */

1804 info = npp_push_tse(npp,

1805 rcv_eq_doublet, sizeof(struct eq_doublet));

1806 info->p = p->i;

1807 info->apq = apq->val;

1808 info->ptr = NULL;

1809 /* transform each row i (i != p), where a[i,q] != 0, to eliminate

1810 column q */

1811 for (aiq = q->ptr; aiq != NULL; aiq = next)

1812 { next = aiq->c_next;

1813 if (aiq == apq) continue; /* skip row p */

1814 i = aiq->row; /* row i to be transformed */

1815 /* save constraint coefficient a[i,q] */

1816 if (npp->sol != GLP_MIP)

1817 { lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE));

1818 lfe->ref = i->i;

1819 lfe->val = aiq->val;

1820 lfe->next = info->ptr;

1821 info->ptr = lfe;

1822 }

1823 /* find coefficient a[i,r] in row i */

1824 for (air = i->ptr; air != NULL; air = air->r_next)

1825 if (air->col == r) break;

1826 /* if a[i,r] does not exist, create a[i,r] = 0 */

1827 if (air == NULL)

1828 air = npp_add_aij(npp, i, r, 0.0);

1829 /* compute gamma[i] = a[i,q] / a[p,q] */

1830 gamma = aiq->val / apq->val;

1831 /* (row i) := (row i) - gamma[i] * (row p); see (3)-(6) */

1832 /* new a[i,q] is exact zero due to elimnation; remove it from

1833 row i */

1834 npp_del_aij(npp, aiq);

1835 /* compute new a[i,r] */

1836 air->val -= gamma * apr->val;

1837 /* if new a[i,r] is close to zero due to numeric cancelation,

1838 remove it from row i */

1839 if (fabs(air->val) <= 1e-10)

1840 npp_del_aij(npp, air);

1841 /* compute new lower and upper bounds of row i */

1842 if (i->lb == i->ub)

1843 i->lb = i->ub = (i->lb - gamma * p->lb);

1844 else

1845 { if (i->lb != -DBL_MAX)

1846 i->lb -= gamma * p->lb;

1847 if (i->ub != +DBL_MAX)

1848 i->ub -= gamma * p->lb;

1849 }

1850 }

1851 return q;

1852 }

1853

1854 static int rcv_eq_doublet(NPP *npp, void *_info)

1855 { /* recover row doubleton (equality constraint) */

1856 struct eq_doublet *info = _info;

1857 NPPLFE *lfe;

1858 double gamma, temp;

1859 /* we assume that processing row p is followed by processing

1860 column q as singleton of type "implied slack variable", in

1861 which case row p must always be active equality constraint */

1862 if (npp->sol == GLP_SOL)

1863 { if (npp->r_stat[info->p] != GLP_NS)

1864 { npp_error();

1865 return 1;

1866 }

1867 }

1868 if (npp->sol != GLP_MIP)

1869 { /* compute value of multiplier for row p; see (14) */

1870 temp = npp->r_pi[info->p];

1871 for (lfe = info->ptr; lfe != NULL; lfe = lfe->next)

1872 { gamma = lfe->val / info->apq; /* a[i,q] / a[p,q] */

1873 temp -= gamma * npp->r_pi[lfe->ref];

1874 }

1875 npp->r_pi[info->p] = temp;

1876 }

1877 return 0;

1878 }

1879

1880 /***********************************************************************

1881 * NAME

1882 *

1883 * npp_forcing_row - process forcing row

1884 *

1885 * SYNOPSIS

1886 *

1887 * #include "glpnpp.h"

1888 * int npp_forcing_row(NPP *npp, NPPROW *p, int at);

1889 *

1890 * DESCRIPTION

1891 *

1892 * The routine npp_forcing row processes row p of general format:

1893 *

1894 * L[p] <= sum a[p,j] x[j] <= U[p], (1)

1895 * j

1896 *

1897 * l[j] <= x[j] <= u[j], (2)

1898 *

1899 * where L[p] <= U[p] and l[j] < u[j] for all a[p,j] != 0. It is also

1900 * assumed that:

1901 *

1902 * 1) if at = 0 then |L[p] - U'[p]| <= eps, where U'[p] is implied

1903 * row upper bound (see below), eps is an absolute tolerance for row

1904 * value;

1905 *

1906 * 2) if at = 1 then |U[p] - L'[p]| <= eps, where L'[p] is implied

1907 * row lower bound (see below).

1908 *

1909 * RETURNS

1910 *

1911 * 0 - success;

1912 *

1913 * 1 - cannot fix columns due to too small constraint coefficients.

1914 *

1915 * PROBLEM TRANSFORMATION

1916 *

1917 * Implied lower and upper bounds of row (1) are determined by bounds

1918 * of corresponding columns (variables) as follows:

1919 *

1920 * L'[p] = inf sum a[p,j] x[j] =

1921 * j

1922 * (3)

1923 * = sum a[p,j] l[j] + sum a[p,j] u[j],

1924 * j in Jp j in Jn

1925 *

1926 * U'[p] = sup sum a[p,j] x[j] =

1927 * (4)

1928 * = sum a[p,j] u[j] + sum a[p,j] l[j],

1929 * j in Jp j in Jn

1930 *

1931 * Jp = {j: a[p,j] > 0}, Jn = {j: a[p,j] < 0}. (5)

1932 *

1933 * If L[p] =~ U'[p] (at = 0), solution can be primal feasible only when

1934 * all variables take their boundary values as defined by (4):

1935 *

1936 * ( u[j], if j in Jp

1937 * x[j] = < (6)

1938 * ( l[j], if j in Jn

1939 *

1940 * Similarly, if U[p] =~ L'[p] (at = 1), solution can be primal feasible

1941 * only when all variables take their boundary values as defined by (3):

1942 *

1943 * ( l[j], if j in Jp

1944 * x[j] = < (7)

1945 * ( u[j], if j in Jn

1946 *

1947 * Condition (6) or (7) allows fixing all columns (variables x[j])

1948 * in row (1) on their bounds and then removing them from the problem

1949 * (see the routine npp_fixed_col). Due to this row p becomes redundant,

1950 * so it can be replaced by equivalent free (unbounded) row and also

1951 * removed from the problem (see the routine npp_free_row).

1952 *

1953 * 1. To apply this transformation row (1) should not have coefficients

1954 * whose magnitude is too small, i.e. all a[p,j] should satisfy to

1955 * the following condition:

1956 *

1957 * |a[p,j]| >= eps * max(1, |a[p,k]|), (8)

1958 * k

1959 * where eps is a relative tolerance for constraint coefficients.

1960 * Otherwise, fixing columns may be numerically unreliable and may

1961 * lead to wrong solution.

1962 *

1963 * 2. The routine fixes columns and remove bounds of row p, however,

1964 * it does not remove the row and columns from the problem.

1965 *

1966 * RECOVERING BASIC SOLUTION

1967 *

1968 * In the transformed problem row p being inactive constraint is

1969 * assigned status GLP_BS (as the result of transformation of free

1970 * row), and all columns in this row are assigned status GLP_NS (as the

1971 * result of transformation of fixed columns).

1972 *

1973 * Note that in the dual system of the transformed (as well as original)

1974 * problem every column j in row p corresponds to the following row:

1975 *

1976 * sum a[i,j] pi[i] + a[p,j] pi[p] + lambda[j] = c[j], (9)

1977 * i!=p

1978 *

1979 * from which it follows that:

1980 *

1981 * lambda[j] = c[j] - sum a[i,j] pi[i] - a[p,j] pi[p]. (10)

1982 * i!=p

1983 *

1984 * In the transformed problem values of all multipliers pi[i] are known

1985 * (including pi[i], whose value is zero, since row p is inactive).

1986 * Thus, using formula (10) it is possible to compute values of

1987 * multipliers lambda[j] for all columns in row p.

1988 *

1989 * Note also that in the original problem all columns in row p are

1990 * bounded, not fixed. So status GLP_NS assigned to every such column

1991 * must be changed to GLP_NL or GLP_NU depending on which bound the

1992 * corresponding column has been fixed. This status change may lead to

1993 * dual feasibility violation for solution of the original problem,

1994 * because now column multipliers must satisfy to the following

1995 * condition:

1996 *

1997 * ( >= 0, if status of column j is GLP_NL,

1998 * lambda[j] < (11)

1999 * ( <= 0, if status of column j is GLP_NU.

2000 *

2001 * If this condition holds, solution to the original problem is the

2002 * same as to the transformed problem. Otherwise, we have to perform

2003 * one degenerate pivoting step of the primal simplex method to obtain

2004 * dual feasible (hence, optimal) solution to the original problem as

2005 * follows. If, on problem transformation, row p was made active on its

2006 * lower bound (case at = 0), we change its status to GLP_NL (or GLP_NS)

2007 * and start increasing its multiplier pi[p]. Otherwise, if row p was

2008 * made active on its upper bound (case at = 1), we change its status

2009 * to GLP_NU (or GLP_NS) and start decreasing pi[p]. From (10) it

2010 * follows that:

2011 *

2012 * delta lambda[j] = - a[p,j] * delta pi[p] = - a[p,j] pi[p]. (12)

2013 *

2014 * Simple analysis of formulae (3)-(5) shows that changing pi[p] in the

2015 * specified direction causes increasing lambda[j] for every column j

2016 * assigned status GLP_NL (delta lambda[j] > 0) and decreasing lambda[j]

2017 * for every column j assigned status GLP_NU (delta lambda[j] < 0). It

2018 * is understood that once the last lambda[q], which violates condition

2019 * (11), has reached zero, multipliers lambda[j] for all columns get

2020 * valid signs. Such column q can be determined as follows. Let d[j] be

2021 * initial value of lambda[j] (i.e. reduced cost of column j) in the

2022 * transformed problem computed with formula (10) when pi[p] = 0. Then

2023 * lambda[j] = d[j] + delta lambda[j], and from (12) it follows that

2024 * lambda[j] becomes zero if:

2025 *

2026 * delta lambda[j] = - a[p,j] pi[p] = - d[j] ==>

2027 * (13)

2028 * pi[p] = d[j] / a[p,j].

2029 *

2030 * Therefore, the last column q, for which lambda[q] becomes zero, can

2031 * be determined from the following condition:

2032 *

2033 * |d[q] / a[p,q]| = max |pi[p]| = max |d[j] / a[p,j]|, (14)

2034 * j in D j in D

2035 *

2036 * where D is a set of columns j whose, reduced costs d[j] have invalid

2037 * signs, i.e. violate condition (11). (Thus, if D is empty, solution

2038 * to the original problem is the same as solution to the transformed

2039 * problem, and no correction is needed as was noticed above.) In

2040 * solution to the original problem column q is assigned status GLP_BS,

2041 * since it replaces column of auxiliary variable of row p (becoming

2042 * active) in the basis, and multiplier for row p is assigned its new

2043 * value, which is pi[p] = d[q] / a[p,q]. Note that due to primal

2044 * degeneracy values of all columns having non-zero coefficients in row

2045 * p remain unchanged.

2046 *

2047 * RECOVERING INTERIOR-POINT SOLUTION

2048 *

2049 * Value of multiplier pi[p] in solution to the original problem is

2050 * corrected in the same way as for basic solution. Values of all

2051 * columns having non-zero coefficients in row p remain unchanged.

2052 *

2053 * RECOVERING MIP SOLUTION

2054 *

2055 * None needed. */

2056

2057 struct forcing_col

2058 { /* column fixed on its bound by forcing row */

2059 int j;

2060 /* column reference number */

2061 char stat;

2062 /* original column status:

2063 GLP_NL - fixed on lower bound

2064 GLP_NU - fixed on upper bound */

2065 double a;

2066 /* constraint coefficient a[p,j] */

2067 double c;

2068 /* objective coefficient c[j] */

2069 NPPLFE *ptr;

2070 /* list of non-zero coefficients a[i,j], i != p */

2071 struct forcing_col *next;

2072 /* pointer to another column fixed by forcing row */

2073 };

2074

2075 struct forcing_row

2076 { /* forcing row */

2077 int p;

2078 /* row reference number */

2079 char stat;

2080 /* status assigned to the row if it becomes active:

2081 GLP_NS - active equality constraint

2082 GLP_NL - inequality constraint with lower bound active

2083 GLP_NU - inequality constraint with upper bound active */

2084 struct forcing_col *ptr;

2085 /* list of all columns having non-zero constraint coefficient

2086 a[p,j] in the forcing row */

2087 };

2088

2089 static int rcv_forcing_row(NPP *npp, void *info);

2090

2091 int npp_forcing_row(NPP *npp, NPPROW *p, int at)

2092 { /* process forcing row */

2093 struct forcing_row *info;

2094 struct forcing_col *col = NULL;

2095 NPPCOL *j;

2096 NPPAIJ *apj, *aij;

2097 NPPLFE *lfe;

2098 double big;

2099 xassert(at == 0 || at == 1);

2100 /* determine maximal magnitude of the row coefficients */

2101 big = 1.0;

2102 for (apj = p->ptr; apj != NULL; apj = apj->r_next)

2103 if (big < fabs(apj->val)) big = fabs(apj->val);

2104 /* if there are too small coefficients in the row, transformation

2105 should not be applied */

2106 for (apj = p->ptr; apj != NULL; apj = apj->r_next)

2107 if (fabs(apj->val) < 1e-7 * big) return 1;

2108 /* create transformation stack entry */

2109 info = npp_push_tse(npp,

2110 rcv_forcing_row, sizeof(struct forcing_row));

2111 info->p = p->i;

2112 if (p->lb == p->ub)

2113 { /* equality constraint */

2114 info->stat = GLP_NS;

2115 }

2116 else if (at == 0)

2117 { /* inequality constraint; case L[p] = U'[p] */

2118 info->stat = GLP_NL;

2119 xassert(p->lb != -DBL_MAX);

2120 }

2121 else /* at == 1 */

2122 { /* inequality constraint; case U[p] = L'[p] */

2123 info->stat = GLP_NU;

2124 xassert(p->ub != +DBL_MAX);

2125 }

2126 info->ptr = NULL;

2127 /* scan the forcing row, fix columns at corresponding bounds, and

2128 save column information (the latter is not needed for MIP) */

2129 for (apj = p->ptr; apj != NULL; apj = apj->r_next)

2130 { /* column j has non-zero coefficient in the forcing row */

2131 j = apj->col;

2132 /* it must be non-fixed */

2133 xassert(j->lb < j->ub);

2134 /* allocate stack entry to save column information */

2135 if (npp->sol != GLP_MIP)

2136 { col = dmp_get_atom(npp->stack, sizeof(struct forcing_col));

2137 col->j = j->j;

2138 col->stat = -1; /* will be set below */

2139 col->a = apj->val;

2140 col->c = j->coef;

2141 col->ptr = NULL;

2142 col->next = info->ptr;

2143 info->ptr = col;

2144 }

2145 /* fix column j */

2146 if (at == 0 && apj->val < 0.0 || at != 0 && apj->val > 0.0)

2147 { /* at its lower bound */

2148 if (npp->sol != GLP_MIP)

2149 col->stat = GLP_NL;

2150 xassert(j->lb != -DBL_MAX);

2151 j->ub = j->lb;

2152 }

2153 else

2154 { /* at its upper bound */

2155 if (npp->sol != GLP_MIP)

2156 col->stat = GLP_NU;

2157 xassert(j->ub != +DBL_MAX);

2158 j->lb = j->ub;

2159 }

2160 /* save column coefficients a[i,j], i != p */

2161 if (npp->sol != GLP_MIP)

2162 { for (aij = j->ptr; aij != NULL; aij = aij->c_next)

2163 { if (aij == apj) continue; /* skip a[p,j] */

2164 lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE));

2165 lfe->ref = aij->row->i;

2166 lfe->val = aij->val;

2167 lfe->next = col->ptr;

2168 col->ptr = lfe;

2169 }

2170 }

2171 }

2172 /* make the row free (unbounded) */

2173 p->lb = -DBL_MAX, p->ub = +DBL_MAX;

2174 return 0;

2175 }

2176

2177 static int rcv_forcing_row(NPP *npp, void *_info)

2178 { /* recover forcing row */

2179 struct forcing_row *info = _info;

2180 struct forcing_col *col, *piv;

2181 NPPLFE *lfe;

2182 double d, big, temp;

2183 if (npp->sol == GLP_MIP) goto done;

2184 /* initially solution to the original problem is the same as

2185 to the transformed problem, where row p is inactive constraint

2186 with pi[p] = 0, and all columns are non-basic */

2187 if (npp->sol == GLP_SOL)

2188 { if (npp->r_stat[info->p] != GLP_BS)

2189 { npp_error();

2190 return 1;

2191 }

2192 for (col = info->ptr; col != NULL; col = col->next)

2193 { if (npp->c_stat[col->j] != GLP_NS)

2194 { npp_error();

2195 return 1;

2196 }

2197 npp->c_stat[col->j] = col->stat; /* original status */

2198 }

2199 }

2200 /* compute reduced costs d[j] for all columns with formula (10)

2201 and store them in col.c instead objective coefficients */

2202 for (col = info->ptr; col != NULL; col = col->next)

2203 { d = col->c;

2204 for (lfe = col->ptr; lfe != NULL; lfe = lfe->next)

2205 d -= lfe->val * npp->r_pi[lfe->ref];

2206 col->c = d;

2207 }

2208 /* consider columns j, whose multipliers lambda[j] has wrong

2209 sign in solution to the transformed problem (where lambda[j] =

2210 d[j]), and choose column q, whose multipler lambda[q] reaches

2211 zero last on changing row multiplier pi[p]; see (14) */

2212 piv = NULL, big = 0.0;

2213 for (col = info->ptr; col != NULL; col = col->next)

2214 { d = col->c; /* d[j] */

2215 temp = fabs(d / col->a);

2216 if (col->stat == GLP_NL)

2217 { /* column j has active lower bound */

2218 if (d < 0.0 && big < temp)

2219 piv = col, big = temp;

2220 }

2221 else if (col->stat == GLP_NU)

2222 { /* column j has active upper bound */

2223 if (d > 0.0 && big < temp)

2224 piv = col, big = temp;

2225 }

2226 else

2227 { npp_error();

2228 return 1;

2229 }

2230 }

2231 /* if column q does not exist, no correction is needed */

2232 if (piv != NULL)

2233 { /* correct solution; row p becomes active constraint while

2234 column q becomes basic */

2235 if (npp->sol == GLP_SOL)

2236 { npp->r_stat[info->p] = info->stat;

2237 npp->c_stat[piv->j] = GLP_BS;

2238 }

2239 /* assign new value to row multiplier pi[p] = d[p] / a[p,q] */

2240 npp->r_pi[info->p] = piv->c / piv->a;

2241 }

2242 done: return 0;

2243 }

2244

2245 /***********************************************************************

2246 * NAME

2247 *

2248 * npp_analyze_row - perform general row analysis

2249 *

2250 * SYNOPSIS

2251 *

2252 * #include "glpnpp.h"

2253 * int npp_analyze_row(NPP *npp, NPPROW *p);

2254 *

2255 * DESCRIPTION

2256 *

2257 * The routine npp_analyze_row performs analysis of row p of general

2258 * format:

2259 *

2260 * L[p] <= sum a[p,j] x[j] <= U[p], (1)

2261 * j

2262 *

2263 * l[j] <= x[j] <= u[j], (2)

2264 *

2265 * where L[p] <= U[p] and l[j] <= u[j] for all a[p,j] != 0.

2266 *

2267 * RETURNS

2268 *

2269 * 0x?0 - row lower bound does not exist or is redundant;

2270 *

2271 * 0x?1 - row lower bound can be active;

2272 *

2273 * 0x?2 - row lower bound is a forcing bound;

2274 *

2275 * 0x0? - row upper bound does not exist or is redundant;

2276 *

2277 * 0x1? - row upper bound can be active;

2278 *

2279 * 0x2? - row upper bound is a forcing bound;

2280 *

2281 * 0x33 - row bounds are inconsistent with column bounds.

2282 *

2283 * ALGORITHM

2284 *

2285 * Analysis of row (1) is based on analysis of its implied lower and

2286 * upper bounds, which are determined by bounds of corresponding columns

2287 * (variables) as follows:

2288 *

2289 * L'[p] = inf sum a[p,j] x[j] =

2290 * j

2291 * (3)

2292 * = sum a[p,j] l[j] + sum a[p,j] u[j],

2293 * j in Jp j in Jn

2294 *

2295 * U'[p] = sup sum a[p,j] x[j] =

2296 * (4)

2297 * = sum a[p,j] u[j] + sum a[p,j] l[j],

2298 * j in Jp j in Jn

2299 *

2300 * Jp = {j: a[p,j] > 0}, Jn = {j: a[p,j] < 0}. (5)

2301 *

2302 * (Note that bounds of all columns in row p are assumed to be correct,

2303 * so L'[p] <= U'[p].)

2304 *

2305 * Analysis of row lower bound L[p] includes the following cases:

2306 *

2307 * 1) if L[p] > U'[p] + eps, where eps is an absolute tolerance for row

2308 * value, row lower bound L[p] and implied row upper bound U'[p] are

2309 * inconsistent, ergo, the problem has no primal feasible solution;

2310 *

2311 * 2) if U'[p] - eps <= L[p] <= U'[p] + eps, i.e. if L[p] =~ U'[p],

2312 * the row is a forcing row on its lower bound (see description of

2313 * the routine npp_forcing_row);

2314 *

2315 * 3) if L[p] > L'[p] + eps, row lower bound L[p] can be active (this

2316 * conclusion does not account other rows in the problem);

2317 *

2318 * 4) if L[p] <= L'[p] + eps, row lower bound L[p] cannot be active, so

2319 * it is redundant and can be removed (replaced by -oo).

2320 *

2321 * Analysis of row upper bound U[p] is performed in a similar way and

2322 * includes the following cases:

2323 *

2324 * 1) if U[p] < L'[p] - eps, row upper bound U[p] and implied row lower

2325 * bound L'[p] are inconsistent, ergo the problem has no primal

2326 * feasible solution;

2327 *

2328 * 2) if L'[p] - eps <= U[p] <= L'[p] + eps, i.e. if U[p] =~ L'[p],

2329 * the row is a forcing row on its upper bound (see description of

2330 * the routine npp_forcing_row);

2331 *

2332 * 3) if U[p] < U'[p] - eps, row upper bound U[p] can be active (this

2333 * conclusion does not account other rows in the problem);

2334 *

2335 * 4) if U[p] >= U'[p] - eps, row upper bound U[p] cannot be active, so

2336 * it is redundant and can be removed (replaced by +oo). */

2337

2338 int npp_analyze_row(NPP *npp, NPPROW *p)

2339 { /* perform general row analysis */

2340 NPPAIJ *aij;

2341 int ret = 0x00;

2342 double l, u, eps;

2343 xassert(npp == npp);

2344 /* compute implied lower bound L'[p]; see (3) */

2345 l = 0.0;

2346 for (aij = p->ptr; aij != NULL; aij = aij->r_next)

2347 { if (aij->val > 0.0)

2348 { if (aij->col->lb == -DBL_MAX)

2349 { l = -DBL_MAX;

2350 break;

2351 }

2352 l += aij->val * aij->col->lb;

2353 }

2354 else /* aij->val < 0.0 */

2355 { if (aij->col->ub == +DBL_MAX)

2356 { l = -DBL_MAX;

2357 break;

2358 }

2359 l += aij->val * aij->col->ub;

2360 }

2361 }

2362 /* compute implied upper bound U'[p]; see (4) */

2363 u = 0.0;

2364 for (aij = p->ptr; aij != NULL; aij = aij->r_next)

2365 { if (aij->val > 0.0)

2366 { if (aij->col->ub == +DBL_MAX)

2367 { u = +DBL_MAX;

2368 break;

2369 }

2370 u += aij->val * aij->col->ub;

2371 }

2372 else /* aij->val < 0.0 */

2373 { if (aij->col->lb == -DBL_MAX)

2374 { u = +DBL_MAX;

2375 break;

2376 }

2377 u += aij->val * aij->col->lb;

2378 }

2379 }

2380 /* column bounds are assumed correct, so L'[p] <= U'[p] */

2381 /* check if row lower bound is consistent */

2382 if (p->lb != -DBL_MAX)

2383 { eps = 1e-3 + 1e-6 * fabs(p->lb);

2384 if (p->lb - eps > u)

2385 { ret = 0x33;

2386 goto done;

2387 }

2388 }

2389 /* check if row upper bound is consistent */

2390 if (p->ub != +DBL_MAX)

2391 { eps = 1e-3 + 1e-6 * fabs(p->ub);

2392 if (p->ub + eps < l)

2393 { ret = 0x33;

2394 goto done;

2395 }

2396 }

2397 /* check if row lower bound can be active/forcing */

2398 if (p->lb != -DBL_MAX)

2399 { eps = 1e-9 + 1e-12 * fabs(p->lb);

2400 if (p->lb - eps > l)

2401 { if (p->lb + eps <= u)

2402 ret |= 0x01;

2403 else

2404 ret |= 0x02;

2405 }

2406 }

2407 /* check if row upper bound can be active/forcing */

2408 if (p->ub != +DBL_MAX)

2409 { eps = 1e-9 + 1e-12 * fabs(p->ub);

2410 if (p->ub + eps < u)

2411 { /* check if the upper bound is forcing */

2412 if (p->ub - eps >= l)

2413 ret |= 0x10;

2414 else

2415 ret |= 0x20;

2416 }

2417 }

2418 done: return ret;

2419 }

2420

2421 /***********************************************************************

2422 * NAME

2423 *

2424 * npp_inactive_bound - remove row lower/upper inactive bound

2425 *

2426 * SYNOPSIS

2427 *

2428 * #include "glpnpp.h"

2429 * void npp_inactive_bound(NPP *npp, NPPROW *p, int which);

2430 *

2431 * DESCRIPTION

2432 *

2433 * The routine npp_inactive_bound removes lower (if which = 0) or upper

2434 * (if which = 1) bound of row p:

2435 *

2436 * L[p] <= sum a[p,j] x[j] <= U[p],

2437 *

2438 * which (bound) is assumed to be redundant.

2439 *

2440 * PROBLEM TRANSFORMATION

2441 *

2442 * If which = 0, current lower bound L[p] of row p is assigned -oo.

2443 * If which = 1, current upper bound U[p] of row p is assigned +oo.

2444 *

2445 * RECOVERING BASIC SOLUTION

2446 *

2447 * If in solution to the transformed problem row p is inactive

2448 * constraint (GLP_BS), its status is not changed in solution to the

2449 * original problem. Otherwise, status of row p in solution to the

2450 * original problem is defined by its type before transformation and

2451 * its status in solution to the transformed problem as follows:

2452 *

2453 * +---------------------+-------+---------------+---------------+

2456 * +---------------------+-------+---------------+---------------+

2463 * +---------------------+-------+---------------+---------------+

2464 *

2465 * RECOVERING INTERIOR-POINT SOLUTION

2466 *

2467 * None needed.

2468 *

2469 * RECOVERING MIP SOLUTION

2470 *

2471 * None needed. */

2472

2473 struct inactive_bound

2474 { /* row inactive bound */

2475 int p;

2476 /* row reference number */

2477 char stat;

2478 /* row status (if active constraint) */

2479 };

2480

2481 static int rcv_inactive_bound(NPP *npp, void *info);

2482

2483 void npp_inactive_bound(NPP *npp, NPPROW *p, int which)

2484 { /* remove row lower/upper inactive bound */

2485 struct inactive_bound *info;

2486 if (npp->sol == GLP_SOL)

2487 { /* create transformation stack entry */

2488 info = npp_push_tse(npp,

2489 rcv_inactive_bound, sizeof(struct inactive_bound));

2490 info->p = p->i;

2491 if (p->ub == +DBL_MAX)

2492 info->stat = GLP_NL;

2493 else if (p->lb == -DBL_MAX)

2494 info->stat = GLP_NU;

2495 else if (p->lb != p->ub)

2496 info->stat = (char)(which == 0 ? GLP_NU : GLP_NL);

2497 else

2498 info->stat = GLP_NS;

2499 }

2500 /* remove row inactive bound */

2501 if (which == 0)

2502 { xassert(p->lb != -DBL_MAX);

2503 p->lb = -DBL_MAX;

2504 }

2505 else if (which == 1)

2506 { xassert(p->ub != +DBL_MAX);

2507 p->ub = +DBL_MAX;

2508 }

2509 else

2510 xassert(which != which);

2511 return;

2512 }

2513

2514 static int rcv_inactive_bound(NPP *npp, void *_info)

2515 { /* recover row status */

2516 struct inactive_bound *info = _info;

2517 if (npp->sol != GLP_SOL)

2518 { npp_error();

2519 return 1;

2520 }

2521 if (npp->r_stat[info->p] == GLP_BS)

2522 npp->r_stat[info->p] = GLP_BS;

2523 else

2524 npp->r_stat[info->p] = info->stat;

2525 return 0;

2526 }

2527

2528 /***********************************************************************

2529 * NAME

2530 *

2531 * npp_implied_bounds - determine implied column bounds

2532 *

2533 * SYNOPSIS

2534 *

2535 * #include "glpnpp.h"

2536 * void npp_implied_bounds(NPP *npp, NPPROW *p);

2537 *

2538 * DESCRIPTION

2539 *

2540 * The routine npp_implied_bounds inspects general row (constraint) p:

2541 *

2542 * L[p] <= sum a[p,j] x[j] <= U[p], (1)

2543 *

2544 * l[j] <= x[j] <= u[j], (2)

2545 *

2546 * where L[p] <= U[p] and l[j] <= u[j] for all a[p,j] != 0, to compute

2547 * implied bounds of columns (variables x[j]) in this row.

2548 *

2549 * The routine stores implied column bounds l'[j] and u'[j] in column

2550 * descriptors (NPPCOL); it does not change current column bounds l[j]

2551 * and u[j]. (Implied column bounds can be then used to strengthen the

2552 * current column bounds; see the routines npp_implied_lower and

2553 * npp_implied_upper).

2554 *

2555 * ALGORITHM

2556 *

2557 * Current column bounds (2) define implied lower and upper bounds of

2558 * row (1) as follows:

2559 *

2560 * L'[p] = inf sum a[p,j] x[j] =

2561 * j

2562 * (3)

2563 * = sum a[p,j] l[j] + sum a[p,j] u[j],

2564 * j in Jp j in Jn

2565 *

2566 * U'[p] = sup sum a[p,j] x[j] =

2567 * (4)

2568 * = sum a[p,j] u[j] + sum a[p,j] l[j],

2569 * j in Jp j in Jn

2570 *

2571 * Jp = {j: a[p,j] > 0}, Jn = {j: a[p,j] < 0}. (5)

2572 *

2573 * (Note that bounds of all columns in row p are assumed to be correct,

2574 * so L'[p] <= U'[p].)

2575 *

2576 * If L[p] > L'[p] and/or U[p] < U'[p], the lower and/or upper bound of

2577 * row (1) can be active, in which case such row defines implied bounds

2578 * of its variables.

2579 *

2580 * Let x[k] be some variable having in row (1) coefficient a[p,k] != 0.

2581 * Consider a case when row lower bound can be active (L[p] > L'[p]):

2582 *

2583 * sum a[p,j] x[j] >= L[p] ==>

2584 * j

2585 *

2586 * sum a[p,j] x[j] + a[p,k] x[k] >= L[p] ==>

2587 * j!=k

2588 * (6)

2589 * a[p,k] x[k] >= L[p] - sum a[p,j] x[j] ==>

2590 * j!=k

2591 *

2592 * a[p,k] x[k] >= L[p,k],

2593 *

2594 * where

2595 *

2596 * L[p,k] = inf(L[p] - sum a[p,j] x[j]) =

2597 * j!=k

2598 *

2599 * = L[p] - sup sum a[p,j] x[j] = (7)

2600 * j!=k

2601 *

2602 * = L[p] - sum a[p,j] u[j] - sum a[p,j] l[j].

2603 * j in Jp\{k} j in Jn\{k}

2604 *

2605 * Thus:

2606 *

2607 * x[k] >= l'[k] = L[p,k] / a[p,k], if a[p,k] > 0, (8)

2608 *

2609 * x[k] <= u'[k] = L[p,k] / a[p,k], if a[p,k] < 0. (9)

2610 *

2611 * where l'[k] and u'[k] are implied lower and upper bounds of variable

2612 * x[k], resp.

2613 *

2614 * Now consider a similar case when row upper bound can be active

2615 * (U[p] < U'[p]):

2616 *

2617 * sum a[p,j] x[j] <= U[p] ==>

2618 * j

2619 *

2620 * sum a[p,j] x[j] + a[p,k] x[k] <= U[p] ==>

2621 * j!=k

2622 * (10)

2623 * a[p,k] x[k] <= U[p] - sum a[p,j] x[j] ==>

2624 * j!=k

2625 *

2626 * a[p,k] x[k] <= U[p,k],

2627 *

2628 * where:

2629 *

2630 * U[p,k] = sup(U[p] - sum a[p,j] x[j]) =

2631 * j!=k

2632 *

2633 * = U[p] - inf sum a[p,j] x[j] = (11)

2634 * j!=k

2635 *

2636 * = U[p] - sum a[p,j] l[j] - sum a[p,j] u[j].

2637 * j in Jp\{k} j in Jn\{k}

2638 *

2639 * Thus:

2640 *

2641 * x[k] <= u'[k] = U[p,k] / a[p,k], if a[p,k] > 0, (12)

2642 *

2643 * x[k] >= l'[k] = U[p,k] / a[p,k], if a[p,k] < 0. (13)

2644 *

2645 * Note that in formulae (8), (9), (12), and (13) coefficient a[p,k]

2646 * must not be too small in magnitude relatively to other non-zero

2647 * coefficients in row (1), i.e. the following condition must hold:

2648 *

2649 * |a[p,k]| >= eps * max(1, |a[p,j]|), (14)

2650 * j

2651 *

2652 * where eps is a relative tolerance for constraint coefficients.

2653 * Otherwise the implied column bounds can be numerical inreliable. For

2654 * example, using formula (8) for the following inequality constraint:

2655 *

2656 * 1e-12 x1 - x2 - x3 >= 0,

2657 *

2658 * where x1 >= -1, x2, x3, >= 0, may lead to numerically unreliable

2659 * conclusion that x1 >= 0.

2660 *

2661 * Using formulae (8), (9), (12), and (13) to compute implied bounds

2662 * for one variable requires |J| operations, where J = {j: a[p,j] != 0},

2663 * because this needs computing L[p,k] and U[p,k]. Thus, computing

2664 * implied bounds for all variables in row (1) would require |J|^2

2665 * operations, that is not a good technique. However, the total number

2666 * of operations can be reduced to |J| as follows.

2667 *

2668 * Let a[p,k] > 0. Then from (7) and (11) we have:

2669 *

2670 * L[p,k] = L[p] - (U'[p] - a[p,k] u[k]) =

2671 *

2672 * = L[p] - U'[p] + a[p,k] u[k],

2673 *

2674 * U[p,k] = U[p] - (L'[p] - a[p,k] l[k]) =

2675 *

2676 * = U[p] - L'[p] + a[p,k] l[k],

2677 *

2678 * where L'[p] and U'[p] are implied row lower and upper bounds defined

2679 * by formulae (3) and (4). Substituting these expressions into (8) and

2680 * (12) gives:

2681 *

2682 * l'[k] = L[p,k] / a[p,k] = u[k] + (L[p] - U'[p]) / a[p,k], (15)

2683 *

2684 * u'[k] = U[p,k] / a[p,k] = l[k] + (U[p] - L'[p]) / a[p,k]. (16)

2685 *

2686 * Similarly, if a[p,k] < 0, according to (7) and (11) we have:

2687 *

2688 * L[p,k] = L[p] - (U'[p] - a[p,k] l[k]) =

2689 *

2690 * = L[p] - U'[p] + a[p,k] l[k],

2691 *

2692 * U[p,k] = U[p] - (L'[p] - a[p,k] u[k]) =

2693 *

2694 * = U[p] - L'[p] + a[p,k] u[k],

2695 *

2696 * and substituting these expressions into (8) and (12) gives:

2697 *

2698 * l'[k] = U[p,k] / a[p,k] = u[k] + (U[p] - L'[p]) / a[p,k], (17)

2699 *

2700 * u'[k] = L[p,k] / a[p,k] = l[k] + (L[p] - U'[p]) / a[p,k]. (18)

2701 *

2702 * Note that formulae (15)-(18) can be used only if L'[p] and U'[p]

2703 * exist. However, if for some variable x[j] it happens that l[j] = -oo

2704 * and/or u[j] = +oo, values of L'[p] (if a[p,j] > 0) and/or U'[p] (if

2705 * a[p,j] < 0) are undefined. Consider, therefore, the most general

2706 * situation, when some column bounds (2) may not exist.

2707 *

2708 * Let:

2709 *

2710 * J' = {j : (a[p,j] > 0 and l[j] = -oo) or

2711 * (19)

2712 * (a[p,j] < 0 and u[j] = +oo)}.

2713 *

2714 * Then (assuming that row upper bound U[p] can be active) the following

2715 * three cases are possible:

2716 *

2717 * 1) |J'| = 0. In this case L'[p] exists, thus, for all variables x[j]

2718 * in row (1) we can use formulae (16) and (17);

2719 *

2720 * 2) J' = {k}. In this case L'[p] = -oo, however, U[p,k] (11) exists,

2721 * so for variable x[k] we can use formulae (12) and (13). Note that

2722 * for all other variables x[j] (j != k) l'[j] = -oo (if a[p,j] < 0)

2723 * or u'[j] = +oo (if a[p,j] > 0);

2724 *

2725 * 3) |J'| > 1. In this case for all variables x[j] in row [1] we have

2726 * l'[j] = -oo (if a[p,j] < 0) or u'[j] = +oo (if a[p,j] > 0).

2727 *

2728 * Similarly, let:

2729 *

2730 * J'' = {j : (a[p,j] > 0 and u[j] = +oo) or

2731 * (20)

2732 * (a[p,j] < 0 and l[j] = -oo)}.

2733 *

2734 * Then (assuming that row lower bound L[p] can be active) the following

2735 * three cases are possible:

2736 *

2737 * 1) |J''| = 0. In this case U'[p] exists, thus, for all variables x[j]

2738 * in row (1) we can use formulae (15) and (18);

2739 *

2740 * 2) J'' = {k}. In this case U'[p] = +oo, however, L[p,k] (7) exists,

2741 * so for variable x[k] we can use formulae (8) and (9). Note that

2742 * for all other variables x[j] (j != k) l'[j] = -oo (if a[p,j] > 0)

2743 * or u'[j] = +oo (if a[p,j] < 0);

2744 *

2745 * 3) |J''| > 1. In this case for all variables x[j] in row (1) we have

2746 * l'[j] = -oo (if a[p,j] > 0) or u'[j] = +oo (if a[p,j] < 0). */

2747

2748 void npp_implied_bounds(NPP *npp, NPPROW *p)

2749 { NPPAIJ *apj, *apk;

2750 double big, eps, temp;

2751 xassert(npp == npp);

2752 /* initialize implied bounds for all variables and determine

2753 maximal magnitude of row coefficients a[p,j] */

2754 big = 1.0;

2755 for (apj = p->ptr; apj != NULL; apj = apj->r_next)

2756 { apj->col->ll.ll = -DBL_MAX, apj->col->uu.uu = +DBL_MAX;

2757 if (big < fabs(apj->val)) big = fabs(apj->val);

2758 }

2759 eps = 1e-6 * big;

2760 /* process row lower bound (assuming that it can be active) */

2761 if (p->lb != -DBL_MAX)

2762 { apk = NULL;

2763 for (apj = p->ptr; apj != NULL; apj = apj->r_next)

2764 { if (apj->val > 0.0 && apj->col->ub == +DBL_MAX ||

2765 apj->val < 0.0 && apj->col->lb == -DBL_MAX)

2766 { if (apk == NULL)

2767 apk = apj;

2768 else

2769 goto skip1;

2770 }

2771 }

2772 /* if a[p,k] = NULL then |J'| = 0 else J' = { k } */

2773 temp = p->lb;

2774 for (apj = p->ptr; apj != NULL; apj = apj->r_next)

2775 { if (apj == apk)

2776 /* skip a[p,k] */;

2777 else if (apj->val > 0.0)

2778 temp -= apj->val * apj->col->ub;

2779 else /* apj->val < 0.0 */

2780 temp -= apj->val * apj->col->lb;

2781 }

2782 /* compute column implied bounds */

2783 if (apk == NULL)

2784 { /* temp = L[p] - U'[p] */

2785 for (apj = p->ptr; apj != NULL; apj = apj->r_next)

2786 { if (apj->val >= +eps)

2787 { /* l'[j] := u[j] + (L[p] - U'[p]) / a[p,j] */

2788 apj->col->ll.ll = apj->col->ub + temp / apj->val;

2789 }

2790 else if (apj->val <= -eps)

2791 { /* u'[j] := l[j] + (L[p] - U'[p]) / a[p,j] */

2792 apj->col->uu.uu = apj->col->lb + temp / apj->val;

2793 }

2794 }

2795 }

2796 else

2797 { /* temp = L[p,k] */

2798 if (apk->val >= +eps)

2799 { /* l'[k] := L[p,k] / a[p,k] */

2800 apk->col->ll.ll = temp / apk->val;

2801 }

2802 else if (apk->val <= -eps)

2803 { /* u'[k] := L[p,k] / a[p,k] */

2804 apk->col->uu.uu = temp / apk->val;

2805 }

2806 }

2807 skip1: ;

2808 }

2809 /* process row upper bound (assuming that it can be active) */

2810 if (p->ub != +DBL_MAX)

2811 { apk = NULL;

2812 for (apj = p->ptr; apj != NULL; apj = apj->r_next)

2813 { if (apj->val > 0.0 && apj->col->lb == -DBL_MAX ||

2814 apj->val < 0.0 && apj->col->ub == +DBL_MAX)

2815 { if (apk == NULL)

2816 apk = apj;

2817 else

2818 goto skip2;

2819 }

2820 }

2821 /* if a[p,k] = NULL then |J''| = 0 else J'' = { k } */

2822 temp = p->ub;

2823 for (apj = p->ptr; apj != NULL; apj = apj->r_next)

2824 { if (apj == apk)

2825 /* skip a[p,k] */;

2826 else if (apj->val > 0.0)

2827 temp -= apj->val * apj->col->lb;

2828 else /* apj->val < 0.0 */

2829 temp -= apj->val * apj->col->ub;

2830 }

2831 /* compute column implied bounds */

2832 if (apk == NULL)

2833 { /* temp = U[p] - L'[p] */

2834 for (apj = p->ptr; apj != NULL; apj = apj->r_next)

2835 { if (apj->val >= +eps)

2836 { /* u'[j] := l[j] + (U[p] - L'[p]) / a[p,j] */

2837 apj->col->uu.uu = apj->col->lb + temp / apj->val;

2838 }

2839 else if (apj->val <= -eps)

2840 { /* l'[j] := u[j] + (U[p] - L'[p]) / a[p,j] */

2841 apj->col->ll.ll = apj->col->ub + temp / apj->val;

2842 }

2843 }

2844 }

2845 else

2846 { /* temp = U[p,k] */

2847 if (apk->val >= +eps)

2848 { /* u'[k] := U[p,k] / a[p,k] */

2849 apk->col->uu.uu = temp / apk->val;

2850 }

2851 else if (apk->val <= -eps)

2852 { /* l'[k] := U[p,k] / a[p,k] */

2853 apk->col->ll.ll = temp / apk->val;

2854 }

2855 }

2856 skip2: ;

2857 }

2858 return;

2859 }

2860

2861 /* eof */