1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 2 /* */ 3 /* This file is part of the class library */ 4 /* SoPlex --- the Sequential object-oriented simPlex. */ 5 /* */ 6 /* Copyright (c) 1996-2023 Zuse Institute Berlin (ZIB) */ 7 /* */ 8 /* Licensed under the Apache License, Version 2.0 (the "License"); */ 9 /* you may not use this file except in compliance with the License. */ 10 /* You may obtain a copy of the License at */ 11 /* */ 12 /* http://www.apache.org/licenses/LICENSE-2.0 */ 13 /* */ 14 /* Unless required by applicable law or agreed to in writing, software */ 15 /* distributed under the License is distributed on an "AS IS" BASIS, */ 16 /* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */ 17 /* See the License for the specific language governing permissions and */ 18 /* limitations under the License. */ 19 /* */ 20 /* You should have received a copy of the Apache-2.0 license */ 21 /* along with SoPlex; see the file LICENSE. If not email to soplex@zib.de. */ 22 /* */ 23 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 24 25 namespace soplex 26 { 27 /**@todo suspicious: *max is not set, but it is used 28 * (with the default setting *max=1) in selectLeave and selectEnter 29 * The question might be if max shouldn't be updated with themax? 30 * 31 * numCycle and maxCycle are integers. So degeneps will be 32 * exactly delta until numCycle >= maxCycle. Then it will be 33 * 0 until numCycle >= 2 * maxCycle, after wich it becomes 34 * negative. This does not look ok. 35 */ 36 template <class R> 37 R SPxHarrisRT<R>::degenerateEps() const 38 { 39 return this->solver()->delta() 40 * (1.0 - this->solver()->numCycle() / this->solver()->maxCycle()); 41 } 42 43 template <class R> 44 int SPxHarrisRT<R>::maxDelta( 45 R* /*max*/, /* max abs value in upd */ 46 R* val, /* initial and chosen value */ 47 int num, /* # of indices in idx */ 48 const int* idx, /* nonzero indices in upd */ 49 const R* upd, /* update VectorBase<R> for vec */ 50 const R* vec, /* current VectorBase<R> */ 51 const R* low, /* lower bounds for vec */ 52 const R* up /* upper bounds for vec */ 53 ) const 54 { 55 R x; 56 R theval; 57 /**@todo patch suggests using *max instead of themax */ 58 R themax; 59 int sel; 60 int i; 61 R epsilon = this->tolerances()->epsilon(); 62 63 assert(*val >= 0); 64 65 theval = *val; 66 themax = 0; 67 sel = -1; 68 69 while(num--) 70 { 71 i = idx[num]; 72 x = upd[i]; 73 74 if(x > epsilon) 75 { 76 themax = (x > themax) ? x : themax; 77 x = (up[i] - vec[i] + this->delta) / x; 78 79 if(x < theval && up[i] < R(infinity)) 80 theval = x; 81 } 82 else if(x < -epsilon) 83 { 84 themax = (-x > themax) ? -x : themax; 85 x = (low[i] - vec[i] - this->delta) / x; 86 87 if(x < theval && low[i] > R(-infinity)) 88 theval = x; 89 } 90 } 91 92 *val = theval; 93 return sel; 94 } 95 96 /**@todo suspicious: *max is not set, but it is used 97 (with the default setting *max=1) 98 in selectLeave and selectEnter 99 */ 100 template <class R> 101 int SPxHarrisRT<R>::minDelta( 102 R* /*max*/, /* max abs value in upd */ 103 R* val, /* initial and chosen value */ 104 int num, /* # of indices in idx */ 105 const int* idx, /* nonzero indices in upd */ 106 const R* upd, /* update VectorBase<R> for vec */ 107 const R* vec, /* current VectorBase<R> */ 108 const R* low, /* lower bounds for vec */ 109 const R* up /* upper bounds for vec */ 110 ) const 111 { 112 R x; 113 R theval; 114 /**@todo patch suggests using *max instead of themax */ 115 R themax; 116 int sel; 117 int i; 118 R epsilon = this->tolerances()->epsilon(); 119 120 assert(*val < 0); 121 122 theval = *val; 123 themax = 0; 124 sel = -1; 125 126 while(num--) 127 { 128 i = idx[num]; 129 x = upd[i]; 130 131 if(x > epsilon) 132 { 133 themax = (x > themax) ? x : themax; 134 x = (low[i] - vec[i] - this->delta) / x; 135 136 if(x > theval && low[i] > R(-infinity)) 137 theval = x; 138 } 139 else if(x < -epsilon) 140 { 141 themax = (-x > themax) ? -x : themax; 142 x = (up[i] - vec[i] + this->delta) / x; 143 144 if(x > theval && up[i] < R(infinity)) 145 theval = x; 146 } 147 } 148 149 *val = theval; 150 return sel; 151 } 152 153 /** 154 Here comes our implementation of the Harris procedure improved by shifting 155 bounds. The basic idea is to used the tolerated infeasibility within 156 solver()->entertol() for searching numerically stable pivots. 157 158 The algorithms operates in two phases. In a first phase, the maximum \p val 159 is determined, when infeasibility within solver()->entertol() is allowed. In the second 160 phase, between all variables with values < \p val the one is selected which 161 gives the best step forward in the simplex iteration. However, this may not 162 allways yield an improvement. In that case, we shift the variable toward 163 infeasibility and retry. This avoids cycling in the shifted LP. 164 */ 165 template <class R> 166 int SPxHarrisRT<R>::selectLeave(R& val, R, bool) 167 { 168 int i, j; 169 R stab, x, y; 170 R max; 171 R sel; 172 R lastshift; 173 R useeps; 174 int leave = -1; 175 R maxabs = 1; 176 177 R epsilon = this->solver()->epsilon(); 178 R degeneps = degenerateEps(); 179 180 SSVectorBase<R>& upd = this->solver()->fVec().delta(); 181 VectorBase<R>& vec = this->solver()->fVec(); 182 183 const VectorBase<R>& up = this->solver()->ubBound(); 184 const VectorBase<R>& low = this->solver()->lbBound(); 185 186 assert(this->delta > epsilon); 187 assert(epsilon > 0); 188 assert(this->solver()->maxCycle() > 0); 189 190 max = val; 191 lastshift = this->solver()->shift(); 192 193 this->solver()->fVec().delta().setup(); 194 195 if(max > epsilon) 196 { 197 // phase 1: 198 maxDelta( 199 &maxabs, /* max abs value in upd */ 200 &max, /* initial and chosen value */ 201 upd.size(), /* # of indices in upd */ 202 upd.indexMem(), /* nonzero indices in upd */ 203 upd.values(), /* update VectorBase<R> for vec */ 204 vec.get_const_ptr(), /* current VectorBase<R> */ 205 low.get_const_ptr(), /* lower bounds for vec */ 206 up.get_const_ptr() /* upper bounds for vec */ 207 ); 208 209 if(max == val) 210 return -1; 211 212 213 // phase 2: 214 stab = 0; 215 sel = R(-infinity); 216 useeps = maxabs * epsilon * 0.001; 217 218 if(useeps < epsilon) 219 useeps = epsilon; 220 221 for(j = upd.size() - 1; j >= 0; --j) 222 { 223 i = upd.index(j); 224 x = upd[i]; 225 226 if(x > useeps) 227 { 228 y = up[i] - vec[i]; 229 230 if(y < -degeneps) 231 this->solver()->shiftUBbound(i, vec[i]); // ensure simplex improvement 232 else 233 { 234 y /= x; 235 236 if(y <= max && y > sel - epsilon && x > stab) 237 { 238 sel = y; 239 leave = i; 240 stab = x; 241 } 242 } 243 } 244 else if(x < -useeps) 245 { 246 y = low[i] - vec[i]; 247 248 if(y > degeneps) 249 this->solver()->shiftLBbound(i, vec[i]); // ensure simplex improvement 250 else 251 { 252 y /= x; 253 254 if(y <= max && y > sel - epsilon && -x > stab) 255 { 256 sel = y; 257 leave = i; 258 stab = -x; 259 } 260 } 261 } 262 else 263 upd.clearNum(j); 264 } 265 } 266 267 268 else if(max < -epsilon) 269 { 270 // phase 1: 271 minDelta( 272 &maxabs, /* max abs value in upd */ 273 &max, /* initial and chosen value */ 274 upd.size(), /* # of indices in upd */ 275 upd.indexMem(), /* nonzero indices in upd */ 276 upd.values(), /* update VectorBase<R> for vec */ 277 vec.get_const_ptr(), /* current VectorBase<R> */ 278 low.get_const_ptr(), /* lower bounds for vec */ 279 up.get_const_ptr() /* upper bounds for vec */ 280 ); 281 282 if(max == val) 283 return -1; 284 285 // phase 2: 286 stab = 0; 287 sel = R(infinity); 288 useeps = maxabs * epsilon * 0.001; 289 290 if(useeps < epsilon) 291 useeps = epsilon; 292 293 for(j = upd.size() - 1; j >= 0; --j) 294 { 295 i = upd.index(j); 296 x = upd[i]; 297 298 if(x < -useeps) 299 { 300 y = up[i] - vec[i]; 301 302 if(y < -degeneps) 303 this->solver()->shiftUBbound(i, vec[i]); // ensure simplex improvement 304 else 305 { 306 y /= x; 307 308 if(y >= max && y < sel + epsilon && -x > stab) 309 { 310 sel = y; 311 leave = i; 312 stab = -x; 313 } 314 } 315 } 316 else if(x > useeps) 317 { 318 y = low[i] - vec[i]; 319 320 if(y > degeneps) 321 this->solver()->shiftLBbound(i, vec[i]); // ensure simplex improvement 322 else 323 { 324 y /= x; 325 326 if(y >= max && y < sel + epsilon && x > stab) 327 { 328 sel = y; 329 leave = i; 330 stab = x; 331 } 332 } 333 } 334 else 335 upd.clearNum(j); 336 } 337 } 338 339 else 340 return -1; 341 342 343 if(lastshift != this->solver()->shift()) 344 return selectLeave(val, 0, false); 345 346 assert(leave >= 0); 347 348 val = sel; 349 return leave; 350 } 351 352 template <class R> 353 SPxId SPxHarrisRT<R>::selectEnter(R& val, int, bool) 354 { 355 int i, j; 356 SPxId enterId; 357 R stab, x, y; 358 R max = 0.0; 359 R sel = 0.0; 360 R lastshift; 361 R cuseeps; 362 R ruseeps; 363 R cmaxabs = 1; 364 R rmaxabs = 1; 365 int pnr, cnr; 366 367 R minStability = this->tolerances()->scaleAccordingToEpsilon(1e-4); 368 R epsilon = this->solver()->epsilon(); 369 R degeneps = degenerateEps(); 370 371 VectorBase<R>& pvec = this->solver()->pVec(); 372 SSVectorBase<R>& pupd = this->solver()->pVec().delta(); 373 374 VectorBase<R>& cvec = this->solver()->coPvec(); 375 SSVectorBase<R>& cupd = this->solver()->coPvec().delta(); 376 377 const VectorBase<R>& upb = this->solver()->upBound(); 378 const VectorBase<R>& lpb = this->solver()->lpBound(); 379 const VectorBase<R>& ucb = this->solver()->ucBound(); 380 const VectorBase<R>& lcb = this->solver()->lcBound(); 381 382 assert(this->delta > epsilon); 383 assert(epsilon > 0); 384 assert(this->solver()->maxCycle() > 0); 385 386 this->solver()->coPvec().delta().setup(); 387 this->solver()->pVec().delta().setup(); 388 389 if(val > epsilon) 390 { 391 for(;;) 392 { 393 pnr = -1; 394 cnr = -1; 395 max = val; 396 lastshift = this->solver()->shift(); 397 assert(this->delta > epsilon); 398 399 // phase 1: 400 maxDelta( 401 &rmaxabs, /* max abs value in upd */ 402 &max, /* initial and chosen value */ 403 pupd.size(), /* # of indices in pupd */ 404 pupd.indexMem(), /* nonzero indices in pupd */ 405 pupd.values(), /* update VectorBase<R> for vec */ 406 pvec.get_const_ptr(), /* current VectorBase<R> */ 407 lpb.get_const_ptr(), /* lower bounds for vec */ 408 upb.get_const_ptr() /* upper bounds for vec */ 409 ); 410 411 maxDelta( 412 &cmaxabs, /* max abs value in upd */ 413 &max, /* initial and chosen value */ 414 cupd.size(), /* # of indices in cupd */ 415 cupd.indexMem(), /* nonzero indices in cupd */ 416 cupd.values(), /* update VectorBase<R> for vec */ 417 cvec.get_const_ptr(), /* current VectorBase<R> */ 418 lcb.get_const_ptr(), /* lower bounds for vec */ 419 ucb.get_const_ptr() /* upper bounds for vec */ 420 ); 421 422 if(max == val) 423 return enterId; 424 425 426 // phase 2: 427 stab = 0; 428 sel = R(-infinity); 429 ruseeps = rmaxabs * 0.001 * epsilon; 430 431 if(ruseeps < epsilon) 432 ruseeps = epsilon; 433 434 cuseeps = cmaxabs * 0.001 * epsilon; 435 436 if(cuseeps < epsilon) 437 cuseeps = epsilon; 438 439 for(j = pupd.size() - 1; j >= 0; --j) 440 { 441 i = pupd.index(j); 442 x = pupd[i]; 443 444 if(x > ruseeps) 445 { 446 y = upb[i] - pvec[i]; 447 448 if(y < -degeneps) 449 this->solver()->shiftUPbound(i, pvec[i] - degeneps); 450 else 451 { 452 y /= x; 453 454 if(y <= max && x >= stab) // && y > sel-epsilon 455 { 456 enterId = this->solver()->id(i); 457 sel = y; 458 pnr = i; 459 stab = x; 460 } 461 } 462 } 463 else if(x < -ruseeps) 464 { 465 y = lpb[i] - pvec[i]; 466 467 if(y > degeneps) 468 this->solver()->shiftLPbound(i, pvec[i] + degeneps); 469 else 470 { 471 y /= x; 472 473 if(y <= max && -x >= stab) // && y > sel-epsilon 474 { 475 enterId = this->solver()->id(i); 476 sel = y; 477 pnr = i; 478 stab = -x; 479 } 480 } 481 } 482 else 483 { 484 SPxOut::debug(this, "DHARRI01 removing value {}\n", pupd[i]); 485 pupd.clearNum(j); 486 } 487 } 488 489 for(j = cupd.size() - 1; j >= 0; --j) 490 { 491 i = cupd.index(j); 492 x = cupd[i]; 493 494 if(x > cuseeps) 495 { 496 y = ucb[i] - cvec[i]; 497 498 if(y < -degeneps) 499 this->solver()->shiftUCbound(i, cvec[i] - degeneps); 500 else 501 { 502 y /= x; 503 504 if(y <= max && x >= stab) // && y > sel-epsilon 505 { 506 enterId = this->solver()->coId(i); 507 sel = y; 508 cnr = j; 509 stab = x; 510 } 511 } 512 } 513 else if(x < -cuseeps) 514 { 515 y = lcb[i] - cvec[i]; 516 517 if(y > degeneps) 518 this->solver()->shiftLCbound(i, cvec[i] + degeneps); 519 else 520 { 521 y /= x; 522 523 if(y <= max && -x >= stab) // && y > sel-epsilon 524 { 525 enterId = this->solver()->coId(i); 526 sel = y; 527 cnr = j; 528 stab = -x; 529 } 530 } 531 } 532 else 533 { 534 SPxOut::debug(this, "DHARRI02 removing value {}\n", cupd[i]); 535 cupd.clearNum(j); 536 } 537 } 538 539 if(lastshift == this->solver()->shift()) 540 { 541 if(cnr >= 0) 542 { 543 if(this->solver()->isBasic(enterId)) 544 { 545 cupd.clearNum(cnr); 546 continue; 547 } 548 else 549 break; 550 } 551 else if(pnr >= 0) 552 { 553 pvec[pnr] = this->solver()->vector(pnr) * cvec; 554 555 if(this->solver()->isBasic(enterId)) 556 { 557 pupd.setValue(pnr, 0.0); 558 continue; 559 } 560 else 561 { 562 x = pupd[pnr]; 563 564 if(x > 0) 565 { 566 sel = upb[pnr] - pvec[pnr]; 567 568 if(x < minStability && sel < this->delta) 569 { 570 minStability /= 2.0; 571 this->solver()->shiftUPbound(pnr, pvec[pnr]); 572 continue; 573 } 574 } 575 else 576 { 577 sel = lpb[pnr] - pvec[pnr]; 578 579 if(-x < minStability && -sel < this->delta) 580 { 581 minStability /= 2.0; 582 this->solver()->shiftLPbound(pnr, pvec[pnr]); 583 continue; 584 } 585 } 586 587 sel /= x; 588 } 589 } 590 else 591 { 592 val = 0; 593 enterId.inValidate(); 594 return enterId; 595 } 596 597 if(sel > max) // instability detected => recompute 598 continue; // ratio test with corrected value 599 600 break; 601 } 602 } 603 } 604 else if(val < -epsilon) 605 { 606 for(;;) 607 { 608 pnr = -1; 609 cnr = -1; 610 max = val; 611 lastshift = this->solver()->shift(); 612 assert(this->delta > epsilon); 613 614 615 // phase 1: 616 minDelta 617 ( 618 &rmaxabs, /* max abs value in upd */ 619 &max, /* initial and chosen value */ 620 pupd.size(), /* # of indices in pupd */ 621 pupd.indexMem(), /* nonzero indices in pupd */ 622 pupd.values(), /* update VectorBase<R> for vec */ 623 pvec.get_const_ptr(), /* current VectorBase<R> */ 624 lpb.get_const_ptr(), /* lower bounds for vec */ 625 upb.get_const_ptr() /* upper bounds for vec */ 626 ); 627 628 minDelta 629 ( 630 &cmaxabs, /* max abs value in upd */ 631 &max, /* initial and chosen value */ 632 cupd.size(), /* # of indices in cupd */ 633 cupd.indexMem(), /* nonzero indices in cupd */ 634 cupd.values(), /* update VectorBase<R> for vec */ 635 cvec.get_const_ptr(), /* current VectorBase<R> */ 636 lcb.get_const_ptr(), /* lower bounds for vec */ 637 ucb.get_const_ptr() /* upper bounds for vec */ 638 ); 639 640 if(max == val) 641 return enterId; 642 643 644 // phase 2: 645 stab = 0; 646 sel = R(infinity); 647 ruseeps = rmaxabs * epsilon * 0.001; 648 cuseeps = cmaxabs * epsilon * 0.001; 649 650 for(j = pupd.size() - 1; j >= 0; --j) 651 { 652 i = pupd.index(j); 653 x = pupd[i]; 654 655 if(x > ruseeps) 656 { 657 y = lpb[i] - pvec[i]; 658 659 if(y > degeneps) 660 this->solver()->shiftLPbound(i, pvec[i]); // ensure simplex improvement 661 else 662 { 663 y /= x; 664 665 if(y >= max && x > stab) // && y < sel+epsilon 666 { 667 enterId = this->solver()->id(i); 668 sel = y; 669 pnr = i; 670 stab = x; 671 } 672 } 673 } 674 else if(x < -ruseeps) 675 { 676 y = upb[i] - pvec[i]; 677 678 if(y < -degeneps) 679 this->solver()->shiftUPbound(i, pvec[i]); // ensure simplex improvement 680 else 681 { 682 y /= x; 683 684 if(y >= max && -x > stab) // && y < sel+epsilon 685 { 686 enterId = this->solver()->id(i); 687 sel = y; 688 pnr = i; 689 stab = -x; 690 } 691 } 692 } 693 else 694 { 695 SPxOut::debug(this, "DHARRI03 removing value {}\n", pupd[i]); 696 pupd.clearNum(j); 697 } 698 } 699 700 for(j = cupd.size() - 1; j >= 0; --j) 701 { 702 i = cupd.index(j); 703 x = cupd[i]; 704 705 if(x > cuseeps) 706 { 707 y = lcb[i] - cvec[i]; 708 709 if(y > degeneps) 710 this->solver()->shiftLCbound(i, cvec[i]); // ensure simplex improvement 711 else 712 { 713 y /= x; 714 715 if(y >= max && x > stab) // && y < sel+epsilon 716 { 717 enterId = this->solver()->coId(i); 718 sel = y; 719 cnr = j; 720 stab = x; 721 } 722 } 723 } 724 else if(x < -cuseeps) 725 { 726 y = ucb[i] - cvec[i]; 727 728 if(y < -degeneps) 729 this->solver()->shiftUCbound(i, cvec[i]); // ensure simplex improvement 730 else 731 { 732 y /= x; 733 734 if(y >= max && -x > stab) // && y < sel+epsilon 735 { 736 enterId = this->solver()->coId(i); 737 sel = y; 738 cnr = j; 739 stab = -x; 740 } 741 } 742 } 743 else 744 { 745 SPxOut::debug(this, "DHARRI04 removing value {}\n", x); 746 cupd.clearNum(j); 747 } 748 } 749 750 if(lastshift == this->solver()->shift()) 751 { 752 if(cnr >= 0) 753 { 754 if(this->solver()->isBasic(enterId)) 755 { 756 cupd.clearNum(cnr); 757 continue; 758 } 759 else 760 break; 761 } 762 else if(pnr >= 0) 763 { 764 pvec[pnr] = this->solver()->vector(pnr) * cvec; 765 766 if(this->solver()->isBasic(enterId)) 767 { 768 pupd.setValue(pnr, 0.0); 769 continue; 770 } 771 else 772 { 773 x = pupd[pnr]; 774 775 if(x > 0) 776 { 777 sel = lpb[pnr] - pvec[pnr]; 778 779 if(x < minStability && -sel < this->delta) 780 { 781 minStability /= 2; 782 this->solver()->shiftLPbound(pnr, pvec[pnr]); 783 continue; 784 } 785 } 786 else 787 { 788 sel = upb[pnr] - pvec[pnr]; 789 790 if(-x < minStability && sel < this->delta) 791 { 792 minStability /= 2; 793 this->solver()->shiftUPbound(pnr, pvec[pnr]); 794 continue; 795 } 796 } 797 798 sel /= x; 799 } 800 } 801 else 802 { 803 val = 0; 804 enterId.inValidate(); 805 return enterId; 806 } 807 808 if(sel < max) // instability detected => recompute 809 continue; // ratio test with corrected value 810 811 break; 812 } 813 } 814 } 815 816 assert(max * val >= 0); 817 assert(enterId.type() != SPxId::INVALID); 818 819 val = sel; 820 821 return enterId; 822 } 823 } // namespace soplex 824