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 #include <iostream> 26 #include <assert.h> 27 28 #include "soplex.h" 29 #include "soplex/statistics.h" 30 #include "soplex/slufactor_rational.h" 31 #include "soplex/ratrecon.h" 32 33 namespace soplex 34 { 35 36 /// solves rational LP 37 template <class R> 38 void SoPlexBase<R>::_optimizeRational(volatile bool* interrupt) 39 { 40 #ifndef SOPLEX_WITH_BOOST 41 SPX_MSG_ERROR(std::cerr << "ERROR: rational solve without Boost not defined!" << std::endl;) 42 return; 43 #else 44 bool hasUnboundedRay = false; 45 bool infeasibilityNotCertified = false; 46 bool unboundednessNotCertified = false; 47 48 // start timing 49 _statistics->solvingTime->start(); 50 _statistics->preprocessingTime->start(); 51 _statistics->initialPrecisionTime->start(); 52 53 // remember that last solve was rational 54 _lastSolveMode = SOLVEMODE_RATIONAL; 55 _resetBoostedPrecision(); 56 57 // ensure that the solver has the original problemo 58 if(!_isRealLPLoaded) 59 { 60 assert(_realLP != &_solver); 61 62 _solver.loadLP(*_realLP); 63 spx_free(_realLP); 64 _realLP = &_solver; 65 _isRealLPLoaded = true; 66 } 67 // during the rational solve, we always store basis information in the basis arrays 68 else if(_hasBasis) 69 { 70 _basisStatusRows.reSize(numRows()); 71 _basisStatusCols.reSize(numCols()); 72 _solver.getBasis(_basisStatusRows.get_ptr(), _basisStatusCols.get_ptr(), _basisStatusRows.size(), 73 _basisStatusCols.size()); 74 } 75 76 // store objective, bounds, and sides of Real LP in case they will be modified during iterative refinement 77 _storeLPReal(); 78 79 // deactivate objective limit in floating-point solver 80 if(realParam(SoPlexBase<R>::OBJLIMIT_LOWER) > -realParam(SoPlexBase<R>::INFTY) 81 || realParam(SoPlexBase<R>::OBJLIMIT_UPPER) < realParam(SoPlexBase<R>::INFTY)) 82 { 83 SPX_MSG_INFO2(spxout, spxout << "Deactivating objective limit.\n"); 84 } 85 86 _solver.setTerminationValue(realParam(SoPlexBase<R>::INFTY)); 87 88 _statistics->preprocessingTime->stop(); 89 90 // apply lifting to reduce range of nonzero matrix coefficients 91 if(boolParam(SoPlexBase<R>::LIFTING)) 92 _lift(); 93 94 // force column representation 95 ///@todo implement row objectives with row representation 96 int oldRepresentation = intParam(SoPlexBase<R>::REPRESENTATION); 97 setIntParam(SoPlexBase<R>::REPRESENTATION, SoPlexBase<R>::REPRESENTATION_COLUMN); 98 99 // force ratio test (avoid bound flipping) 100 int oldRatiotester = intParam(SoPlexBase<R>::RATIOTESTER); 101 102 if(oldRatiotester == SoPlexBase<R>::RATIOTESTER_BOUNDFLIPPING) 103 setIntParam(SoPlexBase<R>::RATIOTESTER, SoPlexBase<R>::RATIOTESTER_FAST); 104 105 ///@todo implement handling of row objectives in Cplex interface 106 #ifdef SOPLEX_WITH_CPX 107 int oldEqtrans = boolParam(SoPlexBase<R>::EQTRANS); 108 setBoolParam(SoPlexBase<R>::EQTRANS, true); 109 #endif 110 111 // introduce slack variables to transform inequality constraints into equations 112 if(boolParam(SoPlexBase<R>::EQTRANS)) 113 _transformEquality(); 114 115 _storedBasis = false; 116 117 bool stoppedTime; 118 bool stoppedIter; 119 120 do 121 { 122 bool primalFeasible = false; 123 bool dualFeasible = false; 124 bool infeasible = false; 125 bool unbounded = false; 126 bool error = false; 127 stoppedTime = false; 128 stoppedIter = false; 129 130 // indicate we are solving the original LP 131 _switchToStandardMode(); 132 133 /// solve problem with iterative refinement and the right precision 134 _performOptIRWrapper(_solRational, !unboundednessNotCertified, !infeasibilityNotCertified, 0, 135 primalFeasible, dualFeasible, infeasible, unbounded, stoppedTime, stoppedIter, error); 136 137 #ifdef SOPLEX_WITH_MPFR 138 139 // case: an unrecoverable error occured, and precision boosting has reached its limit 140 if(error && _boostingLimitReached) 141 { 142 _status = SPxSolverBase<R>::ERROR; 143 break; 144 } 145 // case: an error occured, boost precision 146 else if(error) 147 #else 148 149 // case: an unrecoverable error occured, boost precision 150 if(error) 151 #endif 152 { 153 _status = SPxSolverBase<R>::ERROR; 154 #ifdef SOPLEX_WITH_MPFR 155 156 if(boolParam(SoPlexBase<R>::PRECISION_BOOSTING)) 157 { 158 if(_setupBoostedSolverAfterRecovery()) 159 continue; 160 else // precision boosting has reached its limit 161 break; 162 } 163 else 164 break; 165 166 #else 167 break; 168 #endif 169 } 170 // case: stopped due to some limit 171 else if(stoppedTime) 172 { 173 _status = SPxSolverBase<R>::ABORT_TIME; 174 break; 175 } 176 else if(stoppedIter) 177 { 178 _status = SPxSolverBase<R>::ABORT_ITER; 179 break; 180 } 181 // case: unboundedness detected for the first time 182 else if(unbounded && !unboundednessNotCertified) 183 { 184 SolRational solUnbounded; 185 186 // indicate we are testing unboundedness 187 _switchToUnbdMode(); 188 189 _performUnboundedIRStable(solUnbounded, hasUnboundedRay, stoppedTime, stoppedIter, error); 190 191 assert(!hasUnboundedRay || solUnbounded.hasPrimalRay()); 192 assert(!solUnbounded.hasPrimalRay() || hasUnboundedRay); 193 194 if(error) 195 { 196 SPX_MSG_INFO1(spxout, spxout << "Error while testing for unboundedness.\n"); 197 _status = SPxSolverBase<R>::ERROR; 198 #ifdef SOPLEX_WITH_MPFR 199 200 if(boolParam(SoPlexBase<R>::PRECISION_BOOSTING)) 201 { 202 if(_setupBoostedSolverAfterRecovery()) 203 continue; 204 else // precision boosting has reached its limit 205 break; 206 } 207 else 208 break; 209 210 #else 211 break; 212 #endif 213 } 214 215 if(hasUnboundedRay) 216 { 217 SPX_MSG_INFO1(spxout, spxout << "Dual infeasible. Primal unbounded ray available.\n"); 218 } 219 else 220 { 221 SPX_MSG_INFO1(spxout, spxout << "Dual feasible. Rejecting primal unboundedness.\n"); 222 } 223 224 unboundednessNotCertified = !hasUnboundedRay; 225 226 if(stoppedTime) 227 { 228 _status = SPxSolverBase<R>::ABORT_TIME; 229 break; 230 } 231 else if(stoppedIter) 232 { 233 _status = SPxSolverBase<R>::ABORT_ITER; 234 break; 235 } 236 237 // indicate we are testing feasibility 238 _switchToFeasMode(); 239 240 _performFeasIRStable(_solRational, infeasible, stoppedTime, stoppedIter, error); 241 242 ///@todo this should be stored already earlier, possible switch use solRational above and solFeas here 243 if(hasUnboundedRay) 244 { 245 _solRational._primalRay = solUnbounded._primalRay; 246 _solRational._hasPrimalRay = true; 247 } 248 249 if(error) 250 { 251 SPX_MSG_INFO1(spxout, spxout << "Error while testing for feasibility.\n"); 252 _status = SPxSolverBase<R>::ERROR; 253 #ifdef SOPLEX_WITH_MPFR 254 255 if(boolParam(SoPlexBase<R>::PRECISION_BOOSTING)) 256 { 257 if(_setupBoostedSolverAfterRecovery()) 258 continue; 259 else // precision boosting has reached its limit 260 break; 261 } 262 else 263 break; 264 265 #else 266 break; 267 #endif 268 } 269 else if(stoppedTime) 270 { 271 _status = SPxSolverBase<R>::ABORT_TIME; 272 break; 273 } 274 else if(stoppedIter) 275 { 276 _status = SPxSolverBase<R>::ABORT_ITER; 277 break; 278 } 279 else if(infeasible) 280 { 281 SPX_MSG_INFO1(spxout, spxout << "Primal infeasible. Dual Farkas ray available.\n"); 282 _status = SPxSolverBase<R>::INFEASIBLE; 283 break; 284 } 285 else if(hasUnboundedRay) 286 { 287 SPX_MSG_INFO1(spxout, spxout << "Primal feasible and unbounded.\n"); 288 _status = SPxSolverBase<R>::UNBOUNDED; 289 break; 290 } 291 else 292 { 293 SPX_MSG_INFO1(spxout, spxout << "Primal feasible and bounded.\n"); 294 #ifdef SOPLEX_WITH_MPFR 295 296 if(boolParam(SoPlexBase<R>::PRECISION_BOOSTING)) 297 { 298 if(_setupBoostedSolverAfterRecovery()) 299 continue; 300 else // precision boosting has reached its limit 301 break; 302 } 303 304 #endif 305 continue; 306 } 307 } 308 // case: infeasibility detected 309 else if(infeasible && !infeasibilityNotCertified) 310 { 311 _storeBasis(); 312 313 // indicate we are testing feasibility 314 _switchToFeasMode(); 315 316 _performFeasIRStable(_solRational, infeasible, stoppedTime, stoppedIter, error); 317 318 if(error) 319 { 320 SPX_MSG_INFO1(spxout, spxout << "Error while testing for infeasibility.\n"); 321 _status = SPxSolverBase<R>::ERROR; 322 _restoreBasis(); 323 #ifdef SOPLEX_WITH_MPFR 324 325 if(boolParam(SoPlexBase<R>::PRECISION_BOOSTING)) 326 { 327 if(_setupBoostedSolverAfterRecovery()) 328 continue; 329 else // precision boosting has reached its limit 330 break; 331 } 332 else 333 break; 334 335 #else 336 break; 337 #endif 338 } 339 340 infeasibilityNotCertified = !infeasible; 341 342 if(stoppedTime) 343 { 344 _status = SPxSolverBase<R>::ABORT_TIME; 345 _restoreBasis(); 346 break; 347 } 348 else if(stoppedIter) 349 { 350 _status = SPxSolverBase<R>::ABORT_ITER; 351 _restoreBasis(); 352 break; 353 } 354 355 if(infeasible && boolParam(SoPlexBase<R>::TESTDUALINF)) 356 { 357 SolRational solUnbounded; 358 359 // indicate we are testing unboundedness 360 // not sure about this. 361 _switchToUnbdMode(); 362 363 _performUnboundedIRStable(solUnbounded, hasUnboundedRay, stoppedTime, stoppedIter, error); 364 365 assert(!hasUnboundedRay || solUnbounded.hasPrimalRay()); 366 assert(!solUnbounded.hasPrimalRay() || hasUnboundedRay); 367 368 if(error) 369 { 370 SPX_MSG_INFO1(spxout, spxout << "Error while testing for dual infeasibility.\n"); 371 _status = SPxSolverBase<R>::ERROR; 372 _restoreBasis(); 373 #ifdef SOPLEX_WITH_MPFR 374 375 if(boolParam(SoPlexBase<R>::PRECISION_BOOSTING)) 376 { 377 if(_setupBoostedSolverAfterRecovery()) 378 continue; 379 else // precision boosting has reached its limit 380 break; 381 } 382 else 383 break; 384 385 #else 386 break; 387 #endif 388 } 389 390 if(hasUnboundedRay) 391 { 392 SPX_MSG_INFO1(spxout, spxout << "Dual infeasible. Primal unbounded ray available.\n"); 393 _solRational._primalRay = solUnbounded._primalRay; 394 _solRational._hasPrimalRay = true; 395 } 396 else if(solUnbounded._isDualFeasible) 397 { 398 SPX_MSG_INFO1(spxout, spxout << "Dual feasible. Storing dual multipliers.\n"); 399 _solRational._dual = solUnbounded._dual; 400 _solRational._redCost = solUnbounded._redCost; 401 _solRational._isDualFeasible = true; 402 } 403 else 404 { 405 assert(false); 406 SPX_MSG_INFO1(spxout, spxout << "Not dual infeasible.\n"); 407 } 408 } 409 410 _restoreBasis(); 411 412 if(infeasible) 413 { 414 SPX_MSG_INFO1(spxout, spxout << "Primal infeasible. Dual Farkas ray available.\n"); 415 _status = SPxSolverBase<R>::INFEASIBLE; 416 break; 417 } 418 else if(hasUnboundedRay) 419 { 420 SPX_MSG_INFO1(spxout, spxout << "Primal feasible and unbounded.\n"); 421 _status = SPxSolverBase<R>::UNBOUNDED; 422 break; 423 } 424 else 425 { 426 SPX_MSG_INFO1(spxout, spxout << "Primal feasible. Optimizing again.\n"); 427 #ifdef SOPLEX_WITH_MPFR 428 429 if(boolParam(SoPlexBase<R>::PRECISION_BOOSTING)) 430 { 431 if(_setupBoostedSolverAfterRecovery()) 432 continue; 433 else // precision boosting has reached its limit 434 break; 435 } 436 437 #endif 438 continue; 439 } 440 } 441 else if(primalFeasible && dualFeasible) 442 { 443 SPX_MSG_INFO1(spxout, spxout << "Solved to optimality.\n"); 444 _status = SPxSolverBase<R>::OPTIMAL; 445 break; 446 } 447 else 448 { 449 #ifdef SOPLEX_WITH_MPFR 450 451 if(boolParam(SoPlexBase<R>::PRECISION_BOOSTING)) 452 { 453 if(_setupBoostedSolverAfterRecovery()) 454 continue; 455 else // precision boosting has reached its limit 456 break; 457 } 458 else 459 { 460 SPX_MSG_INFO1(spxout, spxout << "Terminating without success.\n"); 461 break; 462 } 463 464 #else 465 SPX_MSG_INFO1(spxout, spxout << "Terminating without success.\n"); 466 break; 467 #endif 468 } 469 } 470 while(!_isSolveStopped(stoppedTime, stoppedIter)); 471 472 // reset old basis flags for future optimization runs 473 _switchedToBoosted = false; 474 _hasOldBasis = false; 475 _hasOldFeasBasis = false; 476 _hasOldUnbdBasis = false; 477 478 ///@todo set status to ABORT_VALUE if optimal solution exceeds objective limit 479 480 if(_status == SPxSolverBase<R>::OPTIMAL || _status == SPxSolverBase<R>::INFEASIBLE 481 || _status == SPxSolverBase<R>::UNBOUNDED) 482 _hasSolRational = true; 483 484 // restore original problem 485 if(boolParam(SoPlexBase<R>::EQTRANS)) 486 _untransformEquality(_solRational); 487 488 #ifdef SOPLEX_WITH_CPX 489 setBoolParam(SoPlexBase<R>::EQTRANS, oldEqtrans); 490 #endif 491 492 // reset representation and ratio test 493 setIntParam(SoPlexBase<R>::REPRESENTATION, oldRepresentation); 494 setIntParam(SoPlexBase<R>::RATIOTESTER, oldRatiotester); 495 496 // undo lifting 497 if(boolParam(SoPlexBase<R>::LIFTING)) 498 _project(_solRational); 499 500 // restore objective, bounds, and sides of Real LP in case they have been modified during iterative refinement 501 _restoreLPReal(); 502 503 // since the Real LP is loaded in the solver, we need to also pass the basis information to the solver if 504 // available 505 if(_hasBasis) 506 { 507 assert(_isRealLPLoaded); 508 _solver.setBasis(_basisStatusRows.get_const_ptr(), _basisStatusCols.get_const_ptr()); 509 _hasBasis = (_solver.basis().status() > SPxBasisBase<R>::NO_PROBLEM); 510 511 // since setBasis always sets the basis status to regular, we need to set it manually here 512 switch(_status) 513 { 514 case SPxSolverBase<R>::OPTIMAL: 515 _solver.setBasisStatus(SPxBasisBase<R>::OPTIMAL); 516 break; 517 518 case SPxSolverBase<R>::INFEASIBLE: 519 _solver.setBasisStatus(SPxBasisBase<R>::INFEASIBLE); 520 break; 521 522 case SPxSolverBase<R>::UNBOUNDED: 523 _solver.setBasisStatus(SPxBasisBase<R>::UNBOUNDED); 524 break; 525 526 default: 527 break; 528 } 529 530 } 531 532 // stop timing 533 _statistics->initialPrecisionTime->stop(); 534 _statistics->extendedPrecisionTime->stop(); 535 _statistics->solvingTime->stop(); 536 #endif 537 } 538 539 /// perform iterative refinement using the right precision 540 template <class R> 541 void SoPlexBase<R>::_performOptIRWrapper( 542 SolRational& sol, 543 bool acceptUnbounded, 544 bool acceptInfeasible, 545 int minIRRoundsRemaining, 546 bool& primalFeasible, 547 bool& dualFeasible, 548 bool& infeasible, 549 bool& unbounded, 550 bool& stoppedTime, 551 bool& stoppedIter, 552 bool& error 553 ) 554 { 555 _solver.setSolvingForBoosted(boolParam(SoPlexBase<R>::PRECISION_BOOSTING)); 556 _boostedSolver.setSolvingForBoosted(boolParam(SoPlexBase<R>::PRECISION_BOOSTING)); 557 558 #ifdef SOPLEX_WITH_MPFR 559 560 if(boolParam(SoPlexBase<R>::ITERATIVE_REFINEMENT)) 561 { 562 // no precision boosting or initial precision -> solve with plain iterative refinement 563 if(!boolParam(SoPlexBase<R>::PRECISION_BOOSTING) || !_switchedToBoosted) 564 { 565 _performOptIRStable(sol, acceptUnbounded, acceptInfeasible, minIRRoundsRemaining, 566 primalFeasible, dualFeasible, infeasible, unbounded, stoppedTime, stoppedIter, error); 567 } 568 else 569 { 570 do 571 { 572 // setup boosted solver for the new iteration 573 _setupBoostedSolver(); 574 // solve problem with iterative refinement and recovery mechanism, using multiprecision 575 // return false if a new boosted iteration is needed, true otherwise 576 bool needNewBoostedIt; 577 _performOptIRStableBoosted(sol, acceptUnbounded, acceptInfeasible, minIRRoundsRemaining, 578 primalFeasible, dualFeasible, infeasible, unbounded, stoppedTime, stoppedIter, error, 579 needNewBoostedIt); 580 581 // update statistics for precision boosting 582 _updateBoostingStatistics(); 583 584 // boost precision if no success 585 if(needNewBoostedIt) 586 { 587 if(!_boostPrecision()) 588 { 589 // error message already displayed 590 error = true; 591 break; 592 } 593 } 594 else 595 break; 596 } 597 while(true); 598 } 599 } 600 else 601 { 602 assert(boolParam( 603 SoPlexBase<R>::PRECISION_BOOSTING)); // at least one option between IR and precision boosting must be activated 604 605 // initial precision solve 606 if(!_switchedToBoosted) 607 { 608 _solveRealForRationalStable(sol, primalFeasible, dualFeasible, infeasible, unbounded, stoppedTime, 609 stoppedIter, error); 610 } 611 else 612 { 613 do 614 { 615 // setup boosted solver for the new iteration 616 _setupBoostedSolver(); 617 // solve problem with iterative refinement and recovery mechanism, using multiprecision 618 // return false if a new boosted iteration is needed, true otherwise 619 bool needNewBoostedIt; 620 _solveRealForRationalBoostedStable(sol, primalFeasible, dualFeasible, infeasible, unbounded, 621 stoppedTime, stoppedIter, error, needNewBoostedIt); 622 623 // update statistics for precision boosting 624 if(_statistics->precBoosts > 625 1) // if == 1 statistics have already been updated in _solveRealForRationalBoostedStable 626 _updateBoostingStatistics(); 627 628 // boost precision if no success 629 if(needNewBoostedIt) 630 { 631 if(!_boostPrecision()) 632 { 633 // error message already displayed 634 error = true; 635 break; 636 } 637 } 638 else 639 break; 640 } 641 while(true); 642 } 643 } 644 645 #else 646 647 if(boolParam(SoPlexBase<R>::PRECISION_BOOSTING)) 648 { 649 SPX_MSG_ERROR(std::cerr << 650 "ERROR: parameter precision_boosting is set to true but SoPlex was compiled without MPFR support " 651 << std::endl;) 652 error = true; 653 return; 654 } 655 656 // solve problem with iterative refinement and recovery mechanism 657 if(!boolParam(SoPlexBase<R>::ITERATIVE_REFINEMENT)) 658 { 659 SPX_MSG_ERROR(std::cerr << 660 "ERROR: parameter iterative_refinement is set to false but SoPlex was compiled without MPFR support, so boosting is not possible" 661 << std::endl;) 662 error = true; 663 return; 664 } 665 666 _performOptIRStable(sol, acceptUnbounded, acceptInfeasible, minIRRoundsRemaining, 667 primalFeasible, dualFeasible, infeasible, unbounded, stoppedTime, stoppedIter, error); 668 #endif 669 } 670 671 672 673 /// stores floating-point solution of original LP as current rational solution and ensure that solution vectors have right dimension; ensure that solution is aligned with basis 674 /// - updates _rationalLP 675 ///@param sol out 676 ///@param primalReal in 677 ///@param dualReal in 678 ///@param dualSize in & out 679 template <class R> 680 template <typename T> // to support both _solver and _boostedSolver 681 void SoPlexBase<R>::_storeRealSolutionAsRational( 682 SolRational& sol, 683 VectorBase<T>& primalReal, 684 VectorBase<T>& dualReal, 685 int& dualSize) 686 { 687 sol._primal.reDim(numColsRational(), false); 688 sol._slacks.reDim(numRowsRational(), false); 689 sol._dual.reDim(numRowsRational(), false); 690 sol._redCost.reDim(numColsRational(), false); 691 sol._isPrimalFeasible = true; 692 sol._isDualFeasible = true; 693 694 for(int c = numColsRational() - 1; c >= 0; c--) 695 { 696 typename SPxSolverBase<R>::VarStatus& basisStatusCol = _basisStatusCols[c]; 697 698 if(basisStatusCol == SPxSolverBase<R>::ON_LOWER) 699 sol._primal[c] = lowerRational(c); 700 else if(basisStatusCol == SPxSolverBase<R>::ON_UPPER) 701 sol._primal[c] = upperRational(c); 702 else if(basisStatusCol == SPxSolverBase<R>::FIXED) 703 { 704 // it may happen that lower and upper are only equal in the Real LP but different in the rational LP; we do 705 // not check this to avoid rational comparisons, but simply switch the basis status to the lower bound; this 706 // is necessary, because for fixed variables any reduced cost is feasible 707 sol._primal[c] = lowerRational(c); 708 basisStatusCol = SPxSolverBase<R>::ON_LOWER; 709 } 710 else if(basisStatusCol == SPxSolverBase<R>::ZERO) 711 sol._primal[c] = 0; 712 else 713 sol._primal[c].assign(primalReal[c]); 714 } 715 716 _rationalLP->computePrimalActivity(sol._primal, sol._slacks); 717 718 assert(dualSize == 0); 719 720 for(int r = numRowsRational() - 1; r >= 0; r--) 721 { 722 typename SPxSolverBase<R>::VarStatus& basisStatusRow = _basisStatusRows[r]; 723 724 // it may happen that left-hand and right-hand side are different in the rational, but equal in the Real LP, 725 // leading to a fixed basis status; this is critical because rows with fixed basis status are ignored in the 726 // computation of the dual violation; to avoid rational comparisons we do not check this but simply switch to 727 // the left-hand side status 728 if(basisStatusRow == SPxSolverBase<R>::FIXED) 729 basisStatusRow = SPxSolverBase<R>::ON_LOWER; 730 731 { 732 sol._dual[r].assign(dualReal[r]); 733 734 if(dualReal[r] != 0.0) 735 dualSize++; 736 } 737 } 738 739 // we assume that the objective function vector has less nonzeros than the reduced cost vector, and so multiplying 740 // with -1 first and subtracting the dual activity should be faster than adding the dual activity and negating 741 // afterwards 742 _rationalLP->getObj(sol._redCost); 743 _rationalLP->subDualActivity(sol._dual, sol._redCost); 744 } 745 746 747 748 /// computes violation of bounds during the refinement loop 749 template <class R> 750 void SoPlexBase<R>::_computeBoundsViolation(SolRational& sol, Rational& boundsViolation) 751 { 752 boundsViolation = 0; 753 754 for(int c = numColsRational() - 1; c >= 0; c--) 755 { 756 // lower bound 757 assert((lowerRational(c) > _rationalNegInfty) == _lowerFinite(_colTypes[c])); 758 759 if(_lowerFinite(_colTypes[c])) 760 { 761 if(lowerRational(c) == 0) 762 { 763 _modLower[c] = sol._primal[c]; 764 _modLower[c] *= -1; 765 766 if(_modLower[c] > boundsViolation) 767 boundsViolation = _modLower[c]; 768 } 769 else 770 { 771 _modLower[c] = lowerRational(c); 772 _modLower[c] -= sol._primal[c]; 773 774 if(_modLower[c] > boundsViolation) 775 boundsViolation = _modLower[c]; 776 } 777 } 778 779 // upper bound 780 assert((upperRational(c) < _rationalPosInfty) == _upperFinite(_colTypes[c])); 781 782 if(_upperFinite(_colTypes[c])) 783 { 784 if(upperRational(c) == 0) 785 { 786 _modUpper[c] = sol._primal[c]; 787 _modUpper[c] *= -1; 788 789 if(_modUpper[c] < -boundsViolation) 790 boundsViolation = -_modUpper[c]; 791 } 792 else 793 { 794 _modUpper[c] = upperRational(c); 795 _modUpper[c] -= sol._primal[c]; 796 797 if(_modUpper[c] < -boundsViolation) 798 boundsViolation = -_modUpper[c]; 799 } 800 } 801 } 802 } 803 804 805 806 /// computes violation of sides during the refinement loop 807 template <class R> 808 void SoPlexBase<R>::_computeSidesViolation(SolRational& sol, Rational& sideViolation) 809 { 810 sideViolation = 0; 811 812 for(int r = numRowsRational() - 1; r >= 0; r--) 813 { 814 const typename SPxSolverBase<R>::VarStatus& basisStatusRow = _basisStatusRows[r]; 815 816 // left-hand side 817 assert((lhsRational(r) > _rationalNegInfty) == _lowerFinite(_rowTypes[r])); 818 819 if(_lowerFinite(_rowTypes[r])) 820 { 821 if(lhsRational(r) == 0) 822 { 823 _modLhs[r] = sol._slacks[r]; 824 _modLhs[r] *= -1; 825 } 826 else 827 { 828 _modLhs[r] = lhsRational(r); 829 _modLhs[r] -= sol._slacks[r]; 830 } 831 832 if(_modLhs[r] > sideViolation) 833 sideViolation = _modLhs[r]; 834 // if the activity is feasible, but too far from the bound, this violates complementary slackness; we 835 // count it as side violation here 836 else if(basisStatusRow == SPxSolverBase<R>::ON_LOWER && _modLhs[r] < -sideViolation) 837 sideViolation = -_modLhs[r]; 838 } 839 840 // right-hand side 841 assert((rhsRational(r) < _rationalPosInfty) == _upperFinite(_rowTypes[r])); 842 843 if(_upperFinite(_rowTypes[r])) 844 { 845 if(rhsRational(r) == 0) 846 { 847 _modRhs[r] = sol._slacks[r]; 848 _modRhs[r] *= -1; 849 } 850 else 851 { 852 _modRhs[r] = rhsRational(r); 853 _modRhs[r] -= sol._slacks[r]; 854 } 855 856 if(_modRhs[r] < -sideViolation) 857 sideViolation = -_modRhs[r]; 858 // if the activity is feasible, but too far from the bound, this violates complementary slackness; we 859 // count it as side violation here 860 else if(basisStatusRow == SPxSolverBase<R>::ON_UPPER && _modRhs[r] > sideViolation) 861 sideViolation = _modRhs[r]; 862 } 863 } 864 } 865 866 867 868 /// computes reduced cost violation 869 template <class R> 870 void SoPlexBase<R>::_computeReducedCostViolation( 871 SolRational& sol, 872 Rational& redCostViolation, 873 const bool& maximizing) 874 { 875 redCostViolation = 0; 876 877 for(int c = numColsRational() - 1; c >= 0; c--) 878 { 879 if(_colTypes[c] == RANGETYPE_FIXED) 880 continue; 881 882 const typename SPxSolverBase<R>::VarStatus& basisStatusCol = _basisStatusCols[c]; 883 assert(basisStatusCol != SPxSolverBase<R>::FIXED); 884 885 if(((maximizing && basisStatusCol != SPxSolverBase<R>::ON_LOWER) || (!maximizing 886 && basisStatusCol != SPxSolverBase<R>::ON_UPPER)) 887 && sol._redCost[c] < -redCostViolation) 888 { 889 SPxOut::debug(this, "basisStatusCol = {} lower tight = {} upper tight = {} sol._redCost[c] = {}\n", 890 basisStatusCol, bool(sol._primal[c] <= lowerRational(c)), 891 bool(sol._primal[c] >= upperRational(c)), sol._redCost[c].str()); 892 redCostViolation = -sol._redCost[c]; 893 } 894 895 if(((maximizing && basisStatusCol != SPxSolverBase<R>::ON_UPPER) || (!maximizing 896 && basisStatusCol != SPxSolverBase<R>::ON_LOWER)) 897 && sol._redCost[c] > redCostViolation) 898 { 899 SPxOut::debug(this, "basisStatusCol = {} lower tight = {} upper tight = {} sol._redCost[c] = {}\n", 900 basisStatusCol, bool(sol._primal[c] <= lowerRational(c)), 901 bool(sol._primal[c] >= upperRational(c)), sol._redCost[c].str()); 902 redCostViolation = sol._redCost[c]; 903 } 904 } 905 } 906 907 908 909 /// computes dual violation 910 template <class R> 911 void SoPlexBase<R>::_computeDualViolation( 912 SolRational& sol, 913 Rational& dualViolation, 914 const bool& maximizing) 915 { 916 dualViolation = 0; 917 918 for(int r = numRowsRational() - 1; r >= 0; r--) 919 { 920 if(_rowTypes[r] == RANGETYPE_FIXED) 921 continue; 922 923 const typename SPxSolverBase<R>::VarStatus& basisStatusRow = _basisStatusRows[r]; 924 assert(basisStatusRow != SPxSolverBase<R>::FIXED); 925 926 if(((maximizing && basisStatusRow != SPxSolverBase<R>::ON_LOWER) || (!maximizing 927 && basisStatusRow != SPxSolverBase<R>::ON_UPPER)) 928 && sol._dual[r] < -dualViolation) 929 { 930 SPxOut::debug(this, "basisStatusRow = {} lower tight = {} upper tight = {} sol._dual[r] = {}\n", 931 basisStatusRow, bool(sol._slacks[r] <= lhsRational(r)), 932 bool(sol._slacks[r] >= rhsRational(r)), sol._dual[r].str()); 933 dualViolation = -sol._dual[r]; 934 } 935 936 if(((maximizing && basisStatusRow != SPxSolverBase<R>::ON_UPPER) || (!maximizing 937 && basisStatusRow != SPxSolverBase<R>::ON_LOWER)) 938 && sol._dual[r] > dualViolation) 939 { 940 SPxOut::debug(this, "basisStatusRow = {} lower tight = {} upper tight = {} sol._dual[r] = {}\n", 941 basisStatusRow, bool(sol._slacks[r] <= lhsRational(r)), 942 bool(sol._slacks[r] >= rhsRational(r)), sol._dual[r].str()); 943 dualViolation = sol._dual[r]; 944 } 945 } 946 } 947 948 949 950 /// checks termination criteria for refinement loop, return true if termination criteria is met, false otherwise 951 ///@param primalFeasible out 952 ///@param dualFeasible out 953 ///@param ... in 954 template <class R> 955 bool SoPlexBase<R>::_isRefinementOver( 956 bool& primalFeasible, 957 bool& dualFeasible, 958 Rational& boundsViolation, 959 Rational& sideViolation, 960 Rational& redCostViolation, 961 Rational& dualViolation, 962 int minIRRoundsRemaining, 963 bool& stoppedTime, 964 bool& stoppedIter, 965 int numFailedRefinements) 966 { 967 // terminate if tolerances are satisfied 968 primalFeasible = (boundsViolation <= _rationalFeastol && sideViolation <= _rationalFeastol); 969 dualFeasible = (redCostViolation <= _rationalOpttol && dualViolation <= _rationalOpttol); 970 971 if(primalFeasible && dualFeasible) 972 { 973 if(minIRRoundsRemaining < 0) 974 { 975 SPX_MSG_INFO1(spxout, spxout << "Tolerances reached.\n"); 976 return true; 977 } 978 else 979 { 980 SPX_MSG_INFO1(spxout, spxout << 981 "Tolerances reached but minIRRoundsRemaining forcing additional refinement rounds.\n"); 982 } 983 } 984 985 // terminate if some limit is reached 986 if(_isSolveStopped(stoppedTime, stoppedIter) || numFailedRefinements > 2) 987 return true; 988 989 return false; 990 } 991 992 993 994 /// checks refinement loop progress 995 ///@param maxViolation out 996 ///@param bestViolation in & out 997 ///@param violationImprovementFactor in 998 ///@param numFailedRefinements in & out 999 ///@param ... out 1000 template <class R> 1001 void SoPlexBase<R>::_checkRefinementProgress( 1002 Rational& boundsViolation, 1003 Rational& sideViolation, 1004 Rational& redCostViolation, 1005 Rational& dualViolation, 1006 Rational& maxViolation, 1007 Rational& bestViolation, 1008 const Rational& violationImprovementFactor, 1009 int& numFailedRefinements) 1010 { 1011 maxViolation = boundsViolation; 1012 1013 if(sideViolation > maxViolation) 1014 maxViolation = sideViolation; 1015 1016 if(redCostViolation > maxViolation) 1017 maxViolation = redCostViolation; 1018 1019 if(dualViolation > maxViolation) 1020 maxViolation = dualViolation; 1021 1022 bestViolation /= violationImprovementFactor; 1023 1024 if(maxViolation > bestViolation) 1025 { 1026 SPX_MSG_INFO2(spxout, spxout << "Failed to reduce violation significantly.\n"); 1027 bestViolation *= violationImprovementFactor; 1028 numFailedRefinements++; 1029 } 1030 else 1031 bestViolation = maxViolation; 1032 } 1033 1034 1035 1036 /// performs rational reconstruction and/or factorization 1037 ///@param lastStallIterations in 1038 ///@param errorCorrectionFactor in 1039 ///@param dualFeasible out 1040 ///@param stoppedTime out 1041 ///@param stoppedIter out 1042 ///@param breakAfter out 1043 ///@param continueAfter out 1044 ///@param ... in & out 1045 template <class R> 1046 void SoPlexBase<R>::_ratrecAndOrRatfac( 1047 int& minIRRoundsRemaining, 1048 int& lastStallIterations, 1049 int& numberOfIterations, 1050 bool& factorSolNewBasis, 1051 int& nextRatrec, 1052 const Rational& errorCorrectionFactor, 1053 Rational& errorCorrection, 1054 Rational& maxViolation, 1055 SolRational& sol, 1056 bool& primalFeasible, 1057 bool& dualFeasible, 1058 bool& stoppedTime, 1059 bool& stoppedIter, 1060 bool& error, 1061 bool& breakAfter, 1062 bool& continueAfter) 1063 { 1064 1065 ///@todo precision-boosting numberOfIterations is not the best name 1066 // It is supposed to designate either the number of refinements (when we use iterative refinement) or 1067 // the number of precision boosts (when not using iterative refinement) 1068 // numberOfSpecificIterations maybe? 1069 1070 breakAfter = false; 1071 continueAfter = false; 1072 1073 // decide whether to perform rational reconstruction and/or factorization 1074 bool forcebasic = boolParam(SoPlexBase<R>::FORCEBASIC); 1075 bool performRatfac = boolParam(SoPlexBase<R>::RATFAC) 1076 && lastStallIterations >= intParam(SoPlexBase<R>::RATFAC_MINSTALLS) && _hasBasis 1077 && factorSolNewBasis; 1078 bool performRatrec = boolParam(SoPlexBase<R>::RATREC) 1079 && (numberOfIterations >= nextRatrec || performRatfac); 1080 1081 // if we want to force the solution to be basic we need to turn rational factorization on 1082 performRatfac = performRatfac || forcebasic; 1083 1084 // attempt rational reconstruction 1085 errorCorrection *= errorCorrectionFactor; 1086 1087 if(performRatrec && maxViolation > 0) 1088 { 1089 SPX_MSG_INFO1(spxout, spxout << "Performing rational reconstruction . . .\n"); 1090 1091 maxViolation *= errorCorrection; // only used for sign check later 1092 invert(maxViolation); 1093 1094 if(_reconstructSolutionRational(sol, _basisStatusRows, _basisStatusCols, maxViolation)) 1095 { 1096 SPX_MSG_INFO1(spxout, spxout << "Tolerances reached.\n"); 1097 primalFeasible = true; 1098 dualFeasible = true; 1099 1100 if(_hasBasis || !forcebasic) 1101 { 1102 breakAfter = true; 1103 return; 1104 } 1105 } 1106 1107 nextRatrec = int(numberOfIterations * realParam(SoPlexBase<R>::RATREC_FREQ)) + 1; 1108 1109 if(boolParam(SoPlexBase<R>::ITERATIVE_REFINEMENT)) 1110 { 1111 SPxOut::debug(this, "Next rational reconstruction after refinement {}.\n", nextRatrec); 1112 } 1113 else 1114 { 1115 assert(boolParam(SoPlexBase<R>::PRECISION_BOOSTING)); 1116 SPxOut::debug(this, "Next rational reconstruction after precision boost {} \n", nextRatrec); 1117 } 1118 } 1119 1120 // solve basis systems exactly 1121 if((performRatfac && maxViolation > 0) || (!_hasBasis && forcebasic)) 1122 { 1123 SPX_MSG_INFO1(spxout, spxout << "Performing rational factorization . . .\n"); 1124 1125 bool optimal; 1126 _factorizeColumnRational(sol, _basisStatusRows, _basisStatusCols, stoppedTime, stoppedIter, error, 1127 optimal); 1128 factorSolNewBasis = false; 1129 1130 if(stoppedTime) 1131 { 1132 SPX_MSG_INFO1(spxout, spxout << "Stopped rational factorization.\n"); 1133 } 1134 else if(error) 1135 { 1136 // message was already printed; reset error flag and continue without factorization 1137 error = false; 1138 } 1139 else if(optimal) 1140 { 1141 SPX_MSG_INFO1(spxout, spxout << "Tolerances reached.\n"); 1142 primalFeasible = true; 1143 dualFeasible = true; 1144 breakAfter = true; 1145 return; 1146 } 1147 else if(boolParam(SoPlexBase<R>::RATFACJUMP)) 1148 { 1149 SPX_MSG_INFO1(spxout, spxout << "Jumping to exact basic solution.\n"); 1150 minIRRoundsRemaining++; 1151 continueAfter = true; 1152 return; 1153 } 1154 } 1155 } 1156 1157 1158 1159 /// forces value of given nonbasic variable to bound 1160 ///@param sol in & out 1161 ///@param ... in 1162 template <class R> 1163 void SoPlexBase<R>::_forceNonbasicToBound( 1164 SolRational& sol, 1165 int& c, 1166 const int& maxDimRational, 1167 bool toLower) 1168 { 1169 int i = _primalDualDiff.size(); 1170 _ensureDSVectorRationalMemory(_primalDualDiff, maxDimRational); 1171 _primalDualDiff.add(c); 1172 1173 if(toLower) 1174 _primalDualDiff.value(i) = lowerRational(c); 1175 else 1176 _primalDualDiff.value(i) = upperRational(c); 1177 1178 _primalDualDiff.value(i) -= sol._primal[c]; 1179 1180 if(toLower) 1181 sol._primal[c] = lowerRational(c); 1182 else 1183 sol._primal[c] = upperRational(c); 1184 } 1185 1186 1187 1188 /// computes primal scaling factor; limit increase in scaling by tolerance used in floating point solve 1189 ///@param maxScale out 1190 ///@param primalScale in & out 1191 ///@param ... in 1192 template <class R> 1193 void SoPlexBase<R>::_computePrimalScalingFactor( 1194 Rational& maxScale, 1195 Rational& primalScale, 1196 Rational& boundsViolation, 1197 Rational& sideViolation, 1198 Rational& redCostViolation) 1199 { 1200 maxScale = primalScale; 1201 maxScale *= _rationalMaxscaleincr; 1202 1203 primalScale = boundsViolation > sideViolation ? boundsViolation : sideViolation; 1204 1205 if(primalScale < redCostViolation) 1206 primalScale = redCostViolation; 1207 1208 assert(primalScale >= 0); 1209 1210 if(primalScale > 0) 1211 { 1212 invert(primalScale); 1213 1214 if(primalScale > maxScale) 1215 primalScale = maxScale; 1216 } 1217 else 1218 primalScale = maxScale; 1219 1220 if(boolParam(SoPlexBase<R>::POWERSCALING)) 1221 powRound(primalScale); 1222 } 1223 1224 1225 1226 /// computes dual scaling factor; limit increase in scaling by tolerance used in floating point solve 1227 ///@param maxScale out 1228 ///@param dualScale in & out 1229 ///@param ... in 1230 template <class R> 1231 void SoPlexBase<R>::_computeDualScalingFactor( 1232 Rational& maxScale, 1233 Rational& primalScale, 1234 Rational& dualScale, 1235 Rational& redCostViolation, 1236 Rational& dualViolation) 1237 { 1238 maxScale = dualScale; 1239 maxScale *= _rationalMaxscaleincr; 1240 1241 dualScale = redCostViolation > dualViolation ? redCostViolation : dualViolation; 1242 assert(dualScale >= 0); 1243 1244 if(dualScale > 0) 1245 { 1246 invert(dualScale); 1247 1248 if(dualScale > maxScale) 1249 dualScale = maxScale; 1250 } 1251 else 1252 dualScale = maxScale; 1253 1254 if(boolParam(SoPlexBase<R>::POWERSCALING)) 1255 powRound(dualScale); 1256 1257 if(dualScale > primalScale) 1258 dualScale = primalScale; 1259 1260 if(dualScale < 1) 1261 dualScale = 1; 1262 else 1263 { 1264 SPX_MSG_INFO2(spxout, spxout << "Scaling dual by " << dualScale.str() << ".\n"); 1265 1266 // perform dual scaling 1267 ///@todo remove _modObj and use dualScale * sol._redCost directly 1268 _modObj *= dualScale; 1269 } 1270 } 1271 1272 1273 1274 /// applies scaled bounds 1275 ///@param primalScale in & out 1276 template <class R> 1277 template <typename T> 1278 void SoPlexBase<R>::_applyScaledBounds(SPxSolverBase<T>& solver, Rational& primalScale) 1279 { 1280 if(primalScale < 1) 1281 primalScale = 1; 1282 1283 if(primalScale > 1) 1284 SPX_MSG_INFO2(spxout, spxout << "Scaling primal by " << primalScale.str() << ".\n"); 1285 1286 for(int c = numColsRational() - 1; c >= 0; c--) 1287 { 1288 if(_lowerFinite(_colTypes[c])) 1289 { 1290 if(primalScale > 1) 1291 _modLower[c] *= primalScale; 1292 1293 if(_modLower[c] <= _rationalNegInfty) 1294 solver.changeLower(c, -realParam(SoPlexBase<R>::INFTY)); 1295 else if(primalScale > 1) 1296 solver.changeLower(c, T(_modLower[c])); 1297 else 1298 solver.changeLower(c, T(_modLower[c])); 1299 } 1300 1301 if(_upperFinite(_colTypes[c])) 1302 { 1303 if(primalScale > 1) 1304 _modUpper[c] *= primalScale; 1305 1306 if(_modUpper[c] >= _rationalPosInfty) 1307 solver.changeUpper(c, realParam(SoPlexBase<R>::INFTY)); 1308 else 1309 solver.changeUpper(c, T(_modUpper[c])); 1310 } 1311 } 1312 } 1313 1314 1315 1316 /// applies scaled sides, updates: _modLhs, _modRhs, Lhs (resp. Rhs) vector of _solver, 1317 ///@param primalScale in & out 1318 template <class R> 1319 template <typename T> 1320 void SoPlexBase<R>::_applyScaledSides(SPxSolverBase<T>& solver, Rational& primalScale) 1321 { 1322 assert(primalScale >= 1); 1323 1324 for(int r = numRowsRational() - 1; r >= 0; r--) 1325 { 1326 if(_lowerFinite(_rowTypes[r])) 1327 { 1328 if(primalScale != 1) 1329 _modLhs[r] *= primalScale; 1330 1331 if(_modLhs[r] <= _rationalNegInfty) 1332 solver.changeLhs(r, -realParam(SoPlexBase<R>::INFTY)); 1333 else 1334 solver.changeLhs(r, T(_modLhs[r])); 1335 } 1336 1337 if(_upperFinite(_rowTypes[r])) 1338 { 1339 if(primalScale != 1) 1340 _modRhs[r] *= primalScale; 1341 1342 if(_modRhs[r] >= _rationalPosInfty) 1343 solver.changeRhs(r, realParam(SoPlexBase<R>::INFTY)); 1344 else 1345 solver.changeRhs(r, T(_modRhs[r])); 1346 } 1347 } 1348 } 1349 1350 1351 1352 /// applies scaled objective function, updates: objective and row objective of _solver 1353 ///@param dualScale in 1354 ///@param sol in 1355 template <class R> 1356 template <typename T> 1357 void SoPlexBase<R>::_applyScaledObj(SPxSolverBase<T>& solver, Rational& dualScale, SolRational& sol) 1358 { 1359 for(int c = numColsRational() - 1; c >= 0; c--) 1360 { 1361 if(_modObj[c] >= _rationalPosInfty) 1362 solver.changeObj(c, realParam(SoPlexBase<R>::INFTY)); 1363 else if(_modObj[c] <= _rationalNegInfty) 1364 solver.changeObj(c, -realParam(SoPlexBase<R>::INFTY)); 1365 else 1366 solver.changeObj(c, T(_modObj[c])); 1367 } 1368 1369 for(int r = numRowsRational() - 1; r >= 0; r--) 1370 { 1371 Rational newRowObj; 1372 1373 if(_rowTypes[r] == RANGETYPE_FIXED) 1374 solver.changeRowObj(r, T(0.0)); 1375 else 1376 { 1377 newRowObj = sol._dual[r]; 1378 newRowObj *= dualScale; 1379 1380 if(newRowObj >= _rationalPosInfty) 1381 solver.changeRowObj(r, -realParam(SoPlexBase<R>::INFTY)); 1382 else if(newRowObj <= _rationalNegInfty) 1383 solver.changeRowObj(r, realParam(SoPlexBase<R>::INFTY)); 1384 else 1385 solver.changeRowObj(r, -T(newRowObj)); 1386 } 1387 } 1388 } 1389 1390 1391 1392 /// evaluates result of solve. Return true if the algorithm needs to stop, false otherwise. 1393 ///@param result in 1394 ///@param usingRefinedLP in 1395 ///@param dualReal in 1396 ///@param ... out 1397 template <class R> 1398 template <typename T> 1399 bool SoPlexBase<R>::_evaluateResult( 1400 SPxSolverBase<T>& solver, 1401 typename SPxSolverBase<T>::Status result, 1402 bool usingRefinedLP, 1403 SolRational& sol, 1404 VectorBase<T>& dualReal, 1405 bool& infeasible, 1406 bool& unbounded, 1407 bool& stoppedTime, 1408 bool& stoppedIter, 1409 bool& error) 1410 { 1411 // we always evaluate the result after a solve so the first time is a good place to 1412 // set the fp Time 1413 if(_statistics->fpTime == 0) 1414 _statistics->fpTime = _statistics->solvingTime->time(); 1415 1416 if(_statistics->iterationsFP == 0) 1417 _statistics->iterationsFP = _statistics->iterations; 1418 1419 switch(result) 1420 { 1421 case SPxSolverBase<T>::OPTIMAL: 1422 SPX_MSG_INFO1(spxout, spxout << "Floating-point optimal.\n"); 1423 return false; 1424 1425 case SPxSolverBase<T>::INFEASIBLE: 1426 SPX_MSG_INFO1(spxout, spxout << "Floating-point infeasible.\n"); 1427 1428 // when not using refined LP 1429 // the floating-point solve returns a Farkas ray if and only if the simplifier was not used, which is exactly 1430 // the case when a basis could be returned 1431 if(usingRefinedLP || _hasBasis) 1432 { 1433 sol._dualFarkas = dualReal; 1434 sol._hasDualFarkas = true; 1435 } 1436 else 1437 sol._hasDualFarkas = false; 1438 1439 if(usingRefinedLP) 1440 solver.clearRowObjs(); 1441 1442 infeasible = true; 1443 return true; 1444 1445 case SPxSolverBase<T>::UNBOUNDED: 1446 SPX_MSG_INFO1(spxout, spxout << "Floating-point unbounded.\n"); 1447 1448 if(usingRefinedLP) 1449 solver.clearRowObjs(); 1450 1451 unbounded = true; 1452 return true; 1453 1454 case SPxSolverBase<T>::ABORT_TIME: 1455 stoppedTime = true; 1456 return true; 1457 1458 case SPxSolverBase<T>::ABORT_ITER: 1459 if(usingRefinedLP) 1460 solver.clearRowObjs(); 1461 1462 stoppedIter = true; 1463 return true; 1464 1465 default: 1466 if(usingRefinedLP) 1467 solver.clearRowObjs(); 1468 1469 error = true; 1470 return true; 1471 } 1472 } 1473 1474 1475 1476 /// corrects primal solution and aligns with basis 1477 ///@param sol in & out 1478 ///@param primalSize out 1479 ///@param ... in 1480 template <class R> 1481 template <typename T> 1482 void SoPlexBase<R>::_correctPrimalSolution( 1483 SolRational& sol, 1484 Rational& primalScale, 1485 int& primalSize, 1486 const int& maxDimRational, 1487 VectorBase<T>& primalReal) 1488 { 1489 SPxOut::debug(this, "Correcting primal solution.\n"); 1490 1491 primalSize = 0; 1492 Rational primalScaleInverse = primalScale; 1493 invert(primalScaleInverse); 1494 _primalDualDiff.clear(); 1495 1496 for(int c = numColsRational() - 1; c >= 0; c--) 1497 { 1498 // force values of nonbasic variables to bounds 1499 typename SPxSolverBase<R>::VarStatus& basisStatusCol = _basisStatusCols[c]; 1500 1501 if(basisStatusCol == SPxSolverBase<R>::ON_LOWER) 1502 { 1503 if(sol._primal[c] != lowerRational(c)) 1504 { 1505 _forceNonbasicToBound(sol, c, maxDimRational, true); 1506 } 1507 } 1508 else if(basisStatusCol == SPxSolverBase<R>::ON_UPPER) 1509 { 1510 if(sol._primal[c] != upperRational(c)) 1511 { 1512 _forceNonbasicToBound(sol, c, maxDimRational, false); 1513 } 1514 } 1515 else if(basisStatusCol == SPxSolverBase<R>::FIXED) 1516 { 1517 // it may happen that lower and upper are only equal in the Real LP but different in the rational LP; we 1518 // do not check this to avoid rational comparisons, but simply switch the basis status to the lower 1519 // bound; this is necessary, because for fixed variables any reduced cost is feasible 1520 basisStatusCol = SPxSolverBase<R>::ON_LOWER; 1521 1522 if(sol._primal[c] != lowerRational(c)) 1523 { 1524 _forceNonbasicToBound(sol, c, maxDimRational, true); 1525 } 1526 } 1527 else if(basisStatusCol == SPxSolverBase<R>::ZERO) 1528 { 1529 if(sol._primal[c] != 0) 1530 { 1531 int i = _primalDualDiff.size(); 1532 _ensureDSVectorRationalMemory(_primalDualDiff, maxDimRational); 1533 _primalDualDiff.add(c); 1534 _primalDualDiff.value(i) = sol._primal[c]; 1535 _primalDualDiff.value(i) *= -1; 1536 sol._primal[c] = 0; 1537 } 1538 } 1539 else 1540 { 1541 if(primalReal[c] == 1.0) 1542 { 1543 int i = _primalDualDiff.size(); 1544 _ensureDSVectorRationalMemory(_primalDualDiff, maxDimRational); 1545 _primalDualDiff.add(c); 1546 _primalDualDiff.value(i) = primalScaleInverse; 1547 sol._primal[c] += _primalDualDiff.value(i); 1548 } 1549 else if(primalReal[c] == -1.0) 1550 { 1551 int i = _primalDualDiff.size(); 1552 _ensureDSVectorRationalMemory(_primalDualDiff, maxDimRational); 1553 _primalDualDiff.add(c); 1554 _primalDualDiff.value(i) = primalScaleInverse; 1555 _primalDualDiff.value(i) *= -1; 1556 sol._primal[c] += _primalDualDiff.value(i); 1557 } 1558 else if(primalReal[c] != 0.0) 1559 { 1560 int i = _primalDualDiff.size(); 1561 _ensureDSVectorRationalMemory(_primalDualDiff, maxDimRational); 1562 _primalDualDiff.add(c); 1563 _primalDualDiff.value(i).assign(primalReal[c]); 1564 _primalDualDiff.value(i) *= primalScaleInverse; 1565 sol._primal[c] += _primalDualDiff.value(i); 1566 } 1567 } 1568 1569 if(sol._primal[c] != 0) 1570 primalSize++; 1571 } 1572 } 1573 1574 1575 1576 /// updates or recomputes slacks depending on which looks faster 1577 ///@param sol in & out 1578 ///@param primalSize in 1579 template <class R> 1580 void SoPlexBase<R>::_updateSlacks(SolRational& sol, int& primalSize) 1581 { 1582 if(_primalDualDiff.size() < primalSize) 1583 { 1584 _rationalLP->addPrimalActivity(_primalDualDiff, sol._slacks); 1585 #ifndef NDEBUG 1586 { 1587 VectorRational activity(numRowsRational()); 1588 _rationalLP->computePrimalActivity(sol._primal, activity); 1589 assert(sol._slacks == activity); 1590 } 1591 #endif 1592 } 1593 else 1594 _rationalLP->computePrimalActivity(sol._primal, sol._slacks); 1595 } 1596 1597 1598 1599 /// corrects dual solution and aligns with basis 1600 ///@param sol in & out 1601 ///@param dualReal in & out 1602 ///@param dualSize out 1603 ///@param ... in 1604 template <class R> 1605 template <typename T> 1606 void SoPlexBase<R>::_correctDualSolution( 1607 SPxSolverBase<T>& solver, 1608 SolRational& sol, 1609 const bool& maximizing, 1610 VectorBase<T>& dualReal, 1611 Rational& dualScale, 1612 int& dualSize, 1613 const int& maxDimRational) 1614 { 1615 SPxOut::debug(this, "Correcting dual solution.\n"); 1616 1617 #ifndef NDEBUG 1618 { 1619 // compute reduced cost violation 1620 VectorRational debugRedCost(numColsRational()); 1621 debugRedCost = VectorRational(_realLP->maxObj()); 1622 debugRedCost *= -1; 1623 _rationalLP->subDualActivity(VectorRational(dualReal), debugRedCost); 1624 1625 Rational debugRedCostViolation = 0; 1626 1627 for(int c = numColsRational() - 1; c >= 0; c--) 1628 { 1629 if(_colTypes[c] == RANGETYPE_FIXED) 1630 continue; 1631 1632 const typename SPxSolverBase<R>::VarStatus& basisStatusCol = _basisStatusCols[c]; 1633 assert(basisStatusCol != SPxSolverBase<R>::FIXED); 1634 1635 if(((maximizing && basisStatusCol != SPxSolverBase<R>::ON_LOWER) || (!maximizing 1636 && basisStatusCol != SPxSolverBase<R>::ON_UPPER)) 1637 && debugRedCost[c] < -debugRedCostViolation) 1638 { 1639 SPxOut::debug(this, "basisStatusCol = {}, lower tight = {}, upper tight = {}, obj[c] = {}, " 1640 "debugRedCost[c] = {}\n", basisStatusCol, bool(sol._primal[c] <= lowerRational(c)), 1641 bool(sol._primal[c] >= upperRational(c)), _realLP->obj(c), debugRedCost[c].str()); 1642 debugRedCostViolation = -debugRedCost[c]; 1643 } 1644 1645 if(((maximizing && basisStatusCol != SPxSolverBase<R>::ON_UPPER) || (!maximizing 1646 && basisStatusCol != SPxSolverBase<R>::ON_LOWER)) 1647 && debugRedCost[c] > debugRedCostViolation) 1648 { 1649 SPxOut::debug(this, "basisStatusCol = {}, lower tight = {}, upper tight = {}, obj[c] = {}, " 1650 "debugRedCost[c] = {}\n", basisStatusCol, bool(sol._primal[c] <= lowerRational(c)), 1651 bool(sol._primal[c] >= upperRational(c)), _realLP->obj(c), debugRedCost[c].str()); 1652 debugRedCostViolation = debugRedCost[c]; 1653 } 1654 } 1655 1656 // compute dual violation 1657 Rational debugDualViolation = 0; 1658 Rational debugBasicDualViolation = 0; 1659 1660 for(int r = numRowsRational() - 1; r >= 0; r--) 1661 { 1662 if(_rowTypes[r] == RANGETYPE_FIXED) 1663 continue; 1664 1665 const typename SPxSolverBase<R>::VarStatus& basisStatusRow = _basisStatusRows[r]; 1666 assert(basisStatusRow != SPxSolverBase<R>::FIXED); 1667 1668 Rational val = (-dualScale * sol._dual[r]) - Rational(dualReal[r]); 1669 1670 if(((maximizing && basisStatusRow != SPxSolverBase<R>::ON_LOWER) || (!maximizing 1671 && basisStatusRow != SPxSolverBase<R>::ON_UPPER)) 1672 && val > debugDualViolation) 1673 { 1674 SPxOut::debug(this, "basisStatusRow = {} lower tight = {} upper tight = {} dualReal[r] = {} " 1675 "dualReal[r] = {}\n", basisStatusRow, bool(sol._slacks[r] <= lhsRational(r)), 1676 bool(sol._slacks[r] >= rhsRational(r)), val.str(), dualReal[r]); 1677 debugDualViolation = val; 1678 } 1679 1680 if(((maximizing && basisStatusRow != SPxSolverBase<R>::ON_UPPER) || (!maximizing 1681 && basisStatusRow != SPxSolverBase<R>::ON_LOWER)) 1682 && val < -debugDualViolation) 1683 { 1684 SPxOut::debug(this, "basisStatusRow = {} lower tight = {} upper tight = {} dualReal[r] = {} " 1685 "dualReal[r] = {}\n", basisStatusRow, bool(sol._slacks[r] <= lhsRational(r)), 1686 bool(sol._slacks[r] >= rhsRational(r)), val.str(), dualReal[r]); 1687 debugDualViolation = -val; 1688 } 1689 1690 if(basisStatusRow == SPxSolverBase<R>::BASIC && spxAbs(val) > debugBasicDualViolation) 1691 { 1692 SPxOut::debug(this, "basisStatusRow = {} lower tight = {} upper tight = {} dualReal[r] = {} " 1693 "dualReal[r] = {}\n", basisStatusRow, bool(sol._slacks[r] <= lhsRational(r)), 1694 bool(sol._slacks[r] >= rhsRational(r)), val.str(), dualReal[r]); 1695 debugBasicDualViolation = spxAbs(val); 1696 } 1697 } 1698 1699 if(R(debugRedCostViolation) > _solver.tolerances()->floatingPointOpttol() 1700 || R(debugDualViolation) > _solver.tolerances()->floatingPointOpttol() 1701 || debugBasicDualViolation > 1e-9) 1702 { 1703 SPX_MSG_WARNING(spxout, spxout << "Warning: floating-point dual solution with violation " 1704 << debugRedCostViolation.str() << " / " 1705 << debugDualViolation.str() << " / " 1706 << debugBasicDualViolation.str() 1707 << " (red. cost, dual, basic).\n"); 1708 } 1709 } 1710 #endif 1711 1712 Rational dualScaleInverseNeg = dualScale; 1713 invert(dualScaleInverseNeg); 1714 dualScaleInverseNeg *= -1; 1715 _primalDualDiff.clear(); 1716 dualSize = 0; 1717 1718 for(int r = numRowsRational() - 1; r >= 0; r--) 1719 { 1720 typename SPxSolverBase<R>::VarStatus& basisStatusRow = _basisStatusRows[r]; 1721 1722 // it may happen that left-hand and right-hand side are different in the rational, but equal in the Real LP, 1723 // leading to a fixed basis status; this is critical because rows with fixed basis status are ignored in the 1724 // computation of the dual violation; to avoid rational comparisons we do not check this but simply switch 1725 // to the left-hand side status 1726 if(basisStatusRow == SPxSolverBase<R>::FIXED) 1727 basisStatusRow = SPxSolverBase<R>::ON_LOWER; 1728 1729 { 1730 if(dualReal[r] != 0) 1731 { 1732 int i = _primalDualDiff.size(); 1733 _ensureDSVectorRationalMemory(_primalDualDiff, maxDimRational); 1734 _primalDualDiff.add(r); 1735 _primalDualDiff.value(i).assign(dualReal[r]); 1736 _primalDualDiff.value(i) *= dualScaleInverseNeg; 1737 sol._dual[r] -= _primalDualDiff.value(i); 1738 1739 dualSize++; 1740 } 1741 else 1742 { 1743 // we do not check whether the dual value is nonzero, because it probably is; this gives us an 1744 // overestimation of the number of nonzeros in the dual solution 1745 dualSize++; 1746 } 1747 } 1748 } 1749 } 1750 1751 1752 1753 /// updates or recomputes reduced cost values depending on which looks faster; adding one to the length of the 1754 /// dual vector accounts for the objective function vector 1755 ///@param sol in & out 1756 ///@param ... in 1757 template <class R> 1758 void SoPlexBase<R>::_updateReducedCosts(SolRational& sol, int& dualSize, 1759 const int& numCorrectedPrimals) 1760 { 1761 1762 if(_primalDualDiff.size() < dualSize + 1) 1763 { 1764 _rationalLP->addDualActivity(_primalDualDiff, sol._redCost); 1765 #ifndef NDEBUG 1766 { 1767 VectorRational activity(_rationalLP->maxObj()); 1768 activity *= -1; 1769 _rationalLP->subDualActivity(sol._dual, activity); 1770 } 1771 #endif 1772 } 1773 else 1774 { 1775 // we assume that the objective function vector has less nonzeros than the reduced cost vector, and so multiplying 1776 // with -1 first and subtracting the dual activity should be faster than adding the dual activity and negating 1777 // afterwards 1778 _rationalLP->getObj(sol._redCost); 1779 _rationalLP->subDualActivity(sol._dual, sol._redCost); 1780 } 1781 1782 const int numCorrectedDuals = _primalDualDiff.size(); 1783 1784 if(numCorrectedPrimals + numCorrectedDuals > 0) 1785 { 1786 SPX_MSG_INFO2(spxout, spxout << "Corrected " << numCorrectedPrimals << " primal variables and " << 1787 numCorrectedDuals << " dual values.\n"); 1788 } 1789 } 1790 1791 1792 1793 #ifdef SOPLEX_WITH_MPFR 1794 ///@todo precision-boosting move in a different file ? 1795 /// converts the given DataArray of VarStatus to boostedPrecision 1796 template <class R> 1797 void SoPlexBase<R>::_convertDataArrayVarStatusToBoosted( 1798 DataArray< typename SPxSolverBase<R>::VarStatus >& base, 1799 DataArray< typename SPxSolverBase<BP>::VarStatus >& copy 1800 ) 1801 { 1802 using boostedVStat = typename SPxSolverBase<BP>::VarStatus; 1803 1804 // resize copy vector 1805 copy.reSize(base.size()); 1806 1807 // copy while converting type 1808 for(int i = 0; i < base.size(); i++) 1809 { 1810 copy[i] = (boostedVStat)base[i]; 1811 } 1812 } 1813 1814 1815 1816 ///@todo precision-boosting move in a different file ? 1817 /// converts the given DataArray of VarStatus to R precision 1818 template <class R> 1819 void SoPlexBase<R>::_convertDataArrayVarStatusToRPrecision( 1820 DataArray< typename SPxSolverBase<BP>::VarStatus >& base, 1821 DataArray< typename SPxSolverBase<R>::VarStatus >& copy 1822 ) 1823 { 1824 using VStat = typename SPxSolverBase<R>::VarStatus; 1825 1826 // resize copy vector 1827 copy.reSize(base.size()); 1828 1829 // copy while converting type 1830 for(int i = 0; i < base.size(); i++) 1831 { 1832 copy[i] = (VStat)base[i]; 1833 } 1834 } 1835 #endif 1836 1837 1838 1839 /// solves current problem with floating-point solver 1840 template <class R> 1841 void SoPlexBase<R>::_solveRealForRationalStable( 1842 SolRational& sol, 1843 bool& primalFeasible, 1844 bool& dualFeasible, 1845 bool& infeasible, 1846 bool& unbounded, 1847 bool& stoppedTime, 1848 bool& stoppedIter, 1849 bool& error) 1850 { 1851 assert(!boolParam(SoPlexBase<R>::ITERATIVE_REFINEMENT)); 1852 1853 // start rational solving timing 1854 _statistics->rationalTime->start(); 1855 1856 primalFeasible = false; 1857 dualFeasible = false; 1858 infeasible = false; 1859 unbounded = false; 1860 stoppedTime = false; 1861 stoppedIter = false; 1862 error = false; 1863 1864 // declare vectors and variables 1865 typename SPxSolverBase<R>::Status result = SPxSolverBase<R>::UNKNOWN; 1866 1867 _modLower.reDim(numColsRational(), false); 1868 _modUpper.reDim(numColsRational(), false); 1869 _modLhs.reDim(numRowsRational(), false); 1870 _modRhs.reDim(numRowsRational(), false); 1871 _modObj.reDim(numColsRational(), false); 1872 1873 VectorBase<R> primalReal(numColsRational()); 1874 VectorBase<R> dualReal(numRowsRational()); 1875 1876 Rational boundsViolation; 1877 Rational sideViolation; 1878 Rational redCostViolation; 1879 Rational dualViolation; 1880 1881 // solve original LP 1882 SPX_MSG_INFO1(spxout, spxout << "Initial floating-point solve . . .\n"); 1883 1884 if(_hasBasis) 1885 { 1886 assert(_basisStatusRows.size() == numRowsRational()); 1887 assert(_basisStatusCols.size() == numColsRational()); 1888 _solver.setBasis(_basisStatusRows.get_const_ptr(), _basisStatusCols.get_const_ptr()); 1889 _hasBasis = (_solver.basis().status() > SPxBasisBase<R>::NO_PROBLEM); 1890 } 1891 1892 for(int r = numRowsRational() - 1; r >= 0; r--) 1893 { 1894 assert(_solver.maxRowObj(r) == 0.0); 1895 } 1896 1897 _statistics->rationalTime->stop(); 1898 1899 // if boosted solver is available, double precision solver is only used once. 1900 // otherwise use the expensive pipeline from _solveRealStable 1901 result = _solveRealForRational(false, primalReal, dualReal, _basisStatusRows, _basisStatusCols); 1902 1903 // evalute result 1904 if(_evaluateResult<R>(_solver, result, false, sol, dualReal, infeasible, unbounded, stoppedTime, 1905 stoppedIter, 1906 error)) 1907 return; 1908 1909 _statistics->rationalTime->start(); 1910 1911 int dualSize = 0; 1912 1913 // stores floating-point solution of original LP as current rational solution 1914 // ensure that solution vectors have right dimension; ensure that solution is aligned with basis 1915 _storeRealSolutionAsRational<R>(sol, primalReal, dualReal, dualSize); 1916 1917 // control progress 1918 Rational maxViolation; 1919 Rational bestViolation = _rationalPosInfty; 1920 const Rational violationImprovementFactor = 16; 1921 const Rational errorCorrectionFactor = 1.1; 1922 Rational errorCorrection = 2; 1923 1924 // refinement loop 1925 const bool maximizing = (intParam(SoPlexBase<R>::OBJSENSE) == SoPlexBase<R>::OBJSENSE_MAXIMIZE); 1926 1927 // used in _ratrecAndOrRatfac to order a break or a continue outside of the function. 1928 bool breakAfter; 1929 bool continueAfter; 1930 1931 int minIRRoundsRemaining = 0; 1932 1933 do 1934 { 1935 // decrement minIRRoundsRemaining counter 1936 minIRRoundsRemaining--; 1937 1938 SPxOut::debug(this, "Computing primal violations.\n"); 1939 1940 // computes violation of bounds 1941 _computeBoundsViolation(sol, boundsViolation); 1942 1943 // computes violation of sides 1944 _computeSidesViolation(sol, sideViolation); 1945 1946 SPxOut::debug(this, "Computing dual violations.\n"); 1947 1948 // compute reduced cost violation 1949 _computeReducedCostViolation(sol, redCostViolation, maximizing); 1950 1951 // compute dual violation 1952 _computeDualViolation(sol, dualViolation, maximizing); 1953 1954 _modObj = sol._redCost; 1955 1956 // output violations; the reduced cost violations for artificially introduced slack columns are actually violations of the dual multipliers 1957 SPX_MSG_INFO1(spxout, spxout 1958 << "Max. bound violation = " << boundsViolation.str() << ", ofm = 1e" << orderOfMagnitude( 1959 boundsViolation) << "\n" 1960 << "Max. row violation = " << sideViolation.str() << ", ofm = 1e" << orderOfMagnitude( 1961 sideViolation) << "\n" 1962 << "Max. reduced cost violation = " << redCostViolation.str() << ", ofm = 1e" << orderOfMagnitude( 1963 redCostViolation) << "\n" 1964 << "Max. dual violation = " << dualViolation.str() << ", ofm = 1e" << orderOfMagnitude( 1965 dualViolation) << "\n"); 1966 1967 SPxOut::debug(this, "Progress table: {} & {} & {} & {} & {} & {}\n", _statistics->refinements, 1968 _statistics->iterations, _statistics->solvingTime->time(), _statistics->rationalTime->time(), 1969 boundsViolation > sideViolation ? boundsViolation : sideViolation, 1970 redCostViolation > dualViolation ? redCostViolation : dualViolation); 1971 1972 // terminate if tolerances are satisfied 1973 primalFeasible = (boundsViolation <= _rationalFeastol && sideViolation <= _rationalFeastol); 1974 dualFeasible = (redCostViolation <= _rationalOpttol && dualViolation <= _rationalOpttol); 1975 1976 if(primalFeasible && dualFeasible) 1977 { 1978 SPX_MSG_INFO1(spxout, spxout << "Tolerances reached.\n"); 1979 return; 1980 } 1981 1982 1983 // terminate if some limit is reached 1984 if(_isSolveStopped(stoppedTime, stoppedIter)) 1985 return; 1986 1987 // compute maximal violation 1988 maxViolation = boundsViolation; 1989 1990 if(sideViolation > maxViolation) 1991 maxViolation = sideViolation; 1992 1993 if(redCostViolation > maxViolation) 1994 maxViolation = redCostViolation; 1995 1996 if(dualViolation > maxViolation) 1997 maxViolation = dualViolation; 1998 1999 2000 // perform rational reconstruction and/or factorization 2001 _ratrecAndOrRatfac(minIRRoundsRemaining, _lastStallPrecBoosts, _statistics->precBoosts, 2002 _factorSolNewBasisPrecBoost, _nextRatrecPrecBoost, 2003 errorCorrectionFactor, errorCorrection, maxViolation, 2004 sol, primalFeasible, dualFeasible, 2005 stoppedTime, stoppedIter, error, breakAfter, continueAfter); 2006 2007 if(!continueAfter) 2008 break; 2009 2010 } 2011 while(true); 2012 2013 // correct basis status for restricted inequalities 2014 if(_hasBasis) 2015 { 2016 for(int r = numRowsRational() - 1; r >= 0; r--) 2017 { 2018 assert((lhsRational(r) == rhsRational(r)) == (_rowTypes[r] == RANGETYPE_FIXED)); 2019 2020 if(_rowTypes[r] != RANGETYPE_FIXED && _basisStatusRows[r] == SPxSolverBase<R>::FIXED) 2021 _basisStatusRows[r] = (maximizing == (sol._dual[r] < 0)) 2022 ? SPxSolverBase<R>::ON_LOWER 2023 : SPxSolverBase<R>::ON_UPPER; 2024 } 2025 } 2026 2027 // compute objective function values 2028 assert(sol._isPrimalFeasible == sol._isDualFeasible); 2029 2030 if(sol._isPrimalFeasible) 2031 { 2032 sol._objVal = sol._primal * _rationalLP->maxObj(); 2033 2034 if(intParam(SoPlexBase<R>::OBJSENSE) == SoPlexBase<R>::OBJSENSE_MINIMIZE) 2035 sol._objVal *= -1; 2036 } 2037 2038 // set objective coefficients for all rows to zero 2039 _solver.clearRowObjs(); 2040 2041 // stop rational solving time 2042 _statistics->rationalTime->stop(); 2043 } 2044 2045 /// solves current problem with iterative refinement and recovery mechanism 2046 template <class R> 2047 void SoPlexBase<R>::_performOptIRStable( 2048 SolRational& sol, 2049 bool acceptUnbounded, 2050 bool acceptInfeasible, 2051 int minIRRoundsRemaining, 2052 bool& primalFeasible, 2053 bool& dualFeasible, 2054 bool& infeasible, 2055 bool& unbounded, 2056 bool& stoppedTime, 2057 bool& stoppedIter, 2058 bool& error) 2059 { 2060 assert(boolParam(SoPlexBase<R>::ITERATIVE_REFINEMENT)); 2061 2062 // start rational solving timing 2063 _statistics->rationalTime->start(); 2064 2065 primalFeasible = false; 2066 dualFeasible = false; 2067 infeasible = false; 2068 unbounded = false; 2069 stoppedTime = false; 2070 stoppedIter = false; 2071 error = false; 2072 2073 // declare vectors and variables 2074 typename SPxSolverBase<R>::Status result = SPxSolverBase<R>::UNKNOWN; 2075 2076 _modLower.reDim(numColsRational(), false); 2077 _modUpper.reDim(numColsRational(), false); 2078 _modLhs.reDim(numRowsRational(), false); 2079 _modRhs.reDim(numRowsRational(), false); 2080 _modObj.reDim(numColsRational(), false); 2081 2082 VectorBase<R> primalReal(numColsRational()); 2083 VectorBase<R> dualReal(numRowsRational()); 2084 2085 Rational boundsViolation; 2086 Rational sideViolation; 2087 Rational redCostViolation; 2088 Rational dualViolation; 2089 Rational primalScale; 2090 Rational dualScale; 2091 Rational maxScale; 2092 2093 // solve original LP 2094 SPX_MSG_INFO1(spxout, spxout << "Initial floating-point solve . . .\n"); 2095 2096 if(_hasBasis) 2097 { 2098 assert(_basisStatusRows.size() == numRowsRational()); 2099 assert(_basisStatusCols.size() == numColsRational()); 2100 _solver.setBasis(_basisStatusRows.get_const_ptr(), _basisStatusCols.get_const_ptr()); 2101 _hasBasis = (_solver.basis().status() > SPxBasisBase<R>::NO_PROBLEM); 2102 } 2103 2104 for(int r = numRowsRational() - 1; r >= 0; r--) 2105 { 2106 assert(_solver.maxRowObj(r) == 0.0); 2107 } 2108 2109 _statistics->rationalTime->stop(); 2110 2111 // determine if we should use the recovery mechanism pipeline in case of fails 2112 if(!boolParam(SoPlexBase<R>::RECOVERY_MECHANISM)) 2113 result = _solveRealForRational(false, primalReal, dualReal, _basisStatusRows, _basisStatusCols); 2114 else 2115 result = _solveRealStable(acceptUnbounded, acceptInfeasible, primalReal, dualReal, _basisStatusRows, 2116 _basisStatusCols); 2117 2118 // evalute result 2119 if(_evaluateResult<R>(_solver, result, false, sol, dualReal, infeasible, unbounded, stoppedTime, 2120 stoppedIter, 2121 error)) 2122 return; 2123 2124 _statistics->rationalTime->start(); 2125 2126 int dualSize = 0; 2127 2128 // stores floating-point solution of original LP as current rational solution 2129 // ensure that solution vectors have right dimension; ensure that solution is aligned with basis 2130 _storeRealSolutionAsRational<R>(sol, primalReal, dualReal, dualSize); 2131 2132 // initial scaling factors are one 2133 primalScale = _rationalPosone; 2134 dualScale = _rationalPosone; 2135 2136 // control progress 2137 Rational maxViolation; 2138 Rational bestViolation = _rationalPosInfty; 2139 const Rational violationImprovementFactor = 16; 2140 const Rational errorCorrectionFactor = 1.1; 2141 Rational errorCorrection = 2; 2142 int numFailedRefinements = 0; 2143 2144 // store basis status in case solving modified problem failed 2145 DataArray< typename SPxSolverBase<R>::VarStatus > basisStatusRowsFirst; 2146 DataArray< typename SPxSolverBase<R>::VarStatus > basisStatusColsFirst; 2147 2148 // refinement loop 2149 const bool maximizing = (intParam(SoPlexBase<R>::OBJSENSE) == SoPlexBase<R>::OBJSENSE_MAXIMIZE); 2150 const int maxDimRational = numColsRational() > numRowsRational() ? numColsRational() : 2151 numRowsRational(); 2152 bool factorSolNewBasis = true; 2153 int lastStallRefinements = 0; 2154 int nextRatrecRefinement = 0; 2155 2156 // used in _ratrecAndOrRatfac to order a break or a continue outside of the function. 2157 bool breakAfter; 2158 bool continueAfter; 2159 2160 do 2161 { 2162 // decrement minIRRoundsRemaining counter 2163 minIRRoundsRemaining--; 2164 2165 SPxOut::debug(this, "Computing primal violations.\n"); 2166 2167 // computes violation of bounds 2168 _computeBoundsViolation(sol, boundsViolation); 2169 2170 // computes violation of sides 2171 _computeSidesViolation(sol, sideViolation); 2172 2173 SPxOut::debug(this, "Computing dual violations.\n"); 2174 2175 // compute reduced cost violation 2176 _computeReducedCostViolation(sol, redCostViolation, maximizing); 2177 2178 // compute dual violation 2179 _computeDualViolation(sol, dualViolation, maximizing); 2180 2181 _modObj = sol._redCost; 2182 2183 // output violations; the reduced cost violations for artificially introduced slack columns are actually violations of the dual multipliers 2184 SPX_MSG_INFO1(spxout, spxout 2185 << "Max. bound violation = " << boundsViolation.str() << ", ofm = 1e" << orderOfMagnitude( 2186 boundsViolation) << "\n" 2187 << "Max. row violation = " << sideViolation.str() << ", ofm = 1e" << orderOfMagnitude( 2188 sideViolation) << "\n" 2189 << "Max. reduced cost violation = " << redCostViolation.str() << ", ofm = 1e" << orderOfMagnitude( 2190 redCostViolation) << "\n" 2191 << "Max. dual violation = " << dualViolation.str() << ", ofm = 1e" << orderOfMagnitude( 2192 dualViolation) << "\n"); 2193 2194 2195 SPxOut::debug(this, "Progress table: {} & {} & {} & {} & {} & {}\n", _statistics->refinements, 2196 _statistics->iterations, _statistics->solvingTime->time(), _statistics->rationalTime->time(), 2197 boundsViolation > sideViolation ? boundsViolation : sideViolation, 2198 redCostViolation > dualViolation ? redCostViolation : dualViolation); 2199 2200 // check termination criteria for refinement loop 2201 if(_isRefinementOver(primalFeasible, dualFeasible, boundsViolation, sideViolation, redCostViolation, 2202 dualViolation, 2203 minIRRoundsRemaining, stoppedTime, stoppedIter, numFailedRefinements)) 2204 break; 2205 2206 // check refinement progress 2207 _checkRefinementProgress(boundsViolation, sideViolation, redCostViolation, dualViolation, 2208 maxViolation, bestViolation, violationImprovementFactor, numFailedRefinements); 2209 2210 // perform rational reconstruction and/or factorization 2211 _ratrecAndOrRatfac(minIRRoundsRemaining, lastStallRefinements, _statistics->refinements, 2212 factorSolNewBasis, 2213 nextRatrecRefinement, 2214 errorCorrectionFactor, errorCorrection, maxViolation, 2215 sol, primalFeasible, dualFeasible, 2216 stoppedTime, stoppedIter, error, breakAfter, continueAfter); 2217 2218 if(breakAfter) 2219 break; 2220 else if(continueAfter) 2221 continue; 2222 2223 // start refinement 2224 2225 // compute primal scaling factor; limit increase in scaling by tolerance used in floating point solve 2226 _computePrimalScalingFactor(maxScale, primalScale, boundsViolation, sideViolation, 2227 redCostViolation); 2228 2229 // apply scaled bounds and scaled sides 2230 _applyScaledBounds<R>(_solver, primalScale); 2231 _applyScaledSides<R>(_solver, primalScale); 2232 2233 // compute dual scaling factor; limit increase in scaling by tolerance used in floating point solve 2234 _computeDualScalingFactor(maxScale, primalScale, dualScale, redCostViolation, dualViolation); 2235 2236 // apply scaled objective function 2237 _applyScaledObj<R>(_solver, dualScale, sol); 2238 2239 SPX_MSG_INFO1(spxout, spxout << "Refined floating-point solve . . .\n"); 2240 2241 // ensure that artificial slack columns are basic and inequality constraints are nonbasic; otherwise we may end 2242 // up with dual violation on inequality constraints after removing the slack columns; do not change this in the 2243 // floating-point solver, though, because the solver may require its original basis to detect optimality 2244 if(_slackCols.num() > 0 && _hasBasis) 2245 { 2246 int numOrigCols = numColsRational() - _slackCols.num(); 2247 assert(_slackCols.num() <= 0 || boolParam(SoPlexBase<R>::EQTRANS)); 2248 2249 for(int i = 0; i < _slackCols.num(); i++) 2250 { 2251 int row = _slackCols.colVector(i).index(0); 2252 int col = numOrigCols + i; 2253 2254 assert(row >= 0); 2255 assert(row < numRowsRational()); 2256 2257 if(_basisStatusRows[row] == SPxSolverBase<R>::BASIC 2258 && _basisStatusCols[col] != SPxSolverBase<R>::BASIC) 2259 { 2260 _basisStatusRows[row] = _basisStatusCols[col]; 2261 _basisStatusCols[col] = SPxSolverBase<R>::BASIC; 2262 _rationalLUSolver.clear(); 2263 } 2264 } 2265 } 2266 2267 // load basis 2268 if(_hasBasis && _solver.basis().status() < SPxBasisBase<R>::REGULAR) 2269 { 2270 SPxOut::debug(this, "basis (status = {}) desc before set:\n", _solver.basis().status()); 2271 #ifdef SOPLEX_DEBUG 2272 _solver.basis().desc().dump() 2273 #endif 2274 _solver.setBasis(_basisStatusRows.get_const_ptr(), _basisStatusCols.get_const_ptr()); 2275 SPxOut::debug(this, "basis (status = {}) desc after set:\n", _solver.basis().status()); 2276 #ifdef SOPLEX_DEBUG 2277 _solver.basis().desc().dump() 2278 #endif 2279 _hasBasis = _solver.basis().status() > SPxBasisBase<R>::NO_PROBLEM; 2280 SPxOut::debug(this, "setting basis in solver {} (3)\n", (_hasBasis ? "successful" : "failed")); 2281 } 2282 2283 // solve modified problem 2284 int prevIterations = _statistics->iterations; 2285 _statistics->rationalTime->stop(); 2286 2287 // determine if recovery pipeline should be used 2288 if(!boolParam(SoPlexBase<R>::RECOVERY_MECHANISM)) 2289 { 2290 // turn off simplifier if scaling factors are too high 2291 int simplifier = intParam(SoPlexBase<R>::SIMPLIFIER); 2292 2293 if(primalScale > 1e20 || dualScale > 1e20) 2294 setIntParam(SoPlexBase<R>::SIMPLIFIER, SoPlexBase<R>::SIMPLIFIER_OFF); 2295 2296 result = _solveRealForRational(false, primalReal, dualReal, _basisStatusRows, _basisStatusCols); 2297 // reset simplifier param to previous value 2298 setIntParam(SoPlexBase<R>::SIMPLIFIER, simplifier); 2299 } 2300 else 2301 result = _solveRealStable(acceptUnbounded, acceptInfeasible, primalReal, dualReal, _basisStatusRows, 2302 _basisStatusCols, primalScale > 1e20 || dualScale > 1e20); 2303 2304 // count refinements and remember whether we moved to a new basis 2305 _statistics->refinements++; 2306 2307 if(_statistics->iterations <= prevIterations) 2308 { 2309 lastStallRefinements++; 2310 _statistics->stallRefinements++; 2311 } 2312 else 2313 { 2314 factorSolNewBasis = true; 2315 lastStallRefinements = 0; 2316 _statistics->pivotRefinements = _statistics->refinements; 2317 } 2318 2319 // evaluate result; if modified problem was not solved to optimality, stop refinement; 2320 if(_evaluateResult<R>(_solver, result, true, sol, dualReal, infeasible, unbounded, stoppedTime, 2321 stoppedIter, 2322 error)) 2323 return; 2324 2325 _statistics->rationalTime->start(); 2326 2327 int primalSize; 2328 2329 // correct primal solution and align with basis 2330 _correctPrimalSolution<R>(sol, primalScale, primalSize, maxDimRational, primalReal); 2331 2332 // update or recompute slacks depending on which looks faster 2333 _updateSlacks(sol, primalSize); 2334 2335 const int numCorrectedPrimals = _primalDualDiff.size(); 2336 2337 // correct dual solution and align with basis 2338 _correctDualSolution<R>(_solver, sol, maximizing, dualReal, dualScale, dualSize, maxDimRational); 2339 2340 // update or recompute reduced cost values depending on which looks faster 2341 // adding one to the length of the dual vector accounts for the objective function vector 2342 _updateReducedCosts(sol, dualSize, numCorrectedPrimals); 2343 } 2344 while(true); 2345 2346 // correct basis status for restricted inequalities 2347 if(_hasBasis) 2348 { 2349 for(int r = numRowsRational() - 1; r >= 0; r--) 2350 { 2351 assert((lhsRational(r) == rhsRational(r)) == (_rowTypes[r] == RANGETYPE_FIXED)); 2352 2353 if(_rowTypes[r] != RANGETYPE_FIXED && _basisStatusRows[r] == SPxSolverBase<R>::FIXED) 2354 _basisStatusRows[r] = (maximizing == (sol._dual[r] < 0)) 2355 ? SPxSolverBase<R>::ON_LOWER 2356 : SPxSolverBase<R>::ON_UPPER; 2357 } 2358 } 2359 2360 // compute objective function values 2361 assert(sol._isPrimalFeasible == sol._isDualFeasible); 2362 2363 if(sol._isPrimalFeasible) 2364 { 2365 sol._objVal = sol._primal * _rationalLP->maxObj(); 2366 2367 if(intParam(SoPlexBase<R>::OBJSENSE) == SoPlexBase<R>::OBJSENSE_MINIMIZE) 2368 sol._objVal *= -1; 2369 } 2370 2371 // set objective coefficients for all rows to zero 2372 _solver.clearRowObjs(); 2373 2374 // stop rational solving time 2375 _statistics->rationalTime->stop(); 2376 } 2377 2378 2379 2380 /// disable initial precision solver and switch to boosted solver 2381 template <class R> 2382 void SoPlexBase<R>::_switchToBoosted() 2383 { 2384 assert(boolParam(SoPlexBase<R>::PRECISION_BOOSTING)); 2385 2386 if(!_switchedToBoosted) 2387 { 2388 SPX_MSG_INFO1(spxout, spxout << 2389 "Numerical troubles with initial precision solver. Disabling it and switching to multiprecision.\n"); 2390 2391 // switch timer 2392 _statistics->initialPrecisionTime->stop(); 2393 _statistics->extendedPrecisionTime->start(); 2394 2395 _switchedToBoosted = true; 2396 _hasBasis = (_boostedSolver.basis().status() > SPxBasisBase<BP>::NO_PROBLEM); 2397 2398 // invalidate rational basis factorization 2399 _rationalLUSolver.clear(); 2400 } 2401 else 2402 { 2403 SPX_MSG_INFO1(spxout, spxout << 2404 "Numerical troubles with multiprecision solver. Increase precision.\n"); 2405 } 2406 } 2407 2408 2409 2410 /// setup boosted solver before launching iteration 2411 template <class R> 2412 void SoPlexBase<R>::_setupBoostedSolver() 2413 { 2414 assert(boolParam(SoPlexBase<R>::PRECISION_BOOSTING)); 2415 2416 _statistics->boostingStepTime->start(); 2417 2418 // load the data from the rationalLP 2419 _boostedSolver.loadLP(*_rationalLP, false); 2420 2421 // if no warm-start, load available advanced basis 2422 if(!_isBoostedStartingFromSlack()) 2423 _loadBasisFromOldBasis(true); 2424 2425 _hasBasis = _boostedSolver.basis().status() > SPxBasisBase<BP>::NO_PROBLEM; 2426 2427 if(_hasBasis) 2428 { 2429 _tmpBasisStatusRows.reSize(_boostedSolver.nRows()); 2430 _tmpBasisStatusCols.reSize(_boostedSolver.nCols()); 2431 2432 _boostedSolver.getBasis(_tmpBasisStatusRows.get_ptr(), _tmpBasisStatusCols.get_ptr(), 2433 _tmpBasisStatusRows.size(), 2434 _tmpBasisStatusCols.size()); 2435 _convertDataArrayVarStatusToRPrecision(_tmpBasisStatusRows, _basisStatusRows); 2436 _convertDataArrayVarStatusToRPrecision(_tmpBasisStatusCols, _basisStatusCols); 2437 } 2438 2439 _statistics->boostingStepTime->stop(); 2440 } 2441 2442 2443 2444 /// increase the multiprecision, return false if maximum precision is reached, true otherwise 2445 template <class R> 2446 bool SoPlexBase<R>::_boostPrecision() 2447 { 2448 assert(boolParam(SoPlexBase<R>::PRECISION_BOOSTING)); 2449 assert(_switchedToBoosted); 2450 2451 #ifdef SOPLEX_WITH_MPFR 2452 2453 _statistics->boostingStepTime->start(); 2454 2455 _statistics->precBoosts++; 2456 // remember the number of iterations for the next comparison 2457 _prevIterations = _statistics->iterations; 2458 2459 // leave out quad precision and immediately go to 167 bits 2460 if(_statistics->precBoosts == 1) 2461 { 2462 BP nbDigitsSecondBoost = boost::multiprecision::floor(boost::multiprecision::log10( 2463 boost::multiprecision::pow(BP(2), 192))); 2464 2465 if(intParam(SoPlexBase<R>::MULTIPRECISION_LIMIT) >= nbDigitsSecondBoost) 2466 { 2467 BP::default_precision((int)nbDigitsSecondBoost); 2468 } 2469 else 2470 { 2471 SPX_MSG_INFO1(spxout, spxout << "Maximum number of digits for the multiprecision type reached.\n" 2472 << "To increase this limit, modify the parameter multiprecision_limit.\n" 2473 << "Giving up.\n"); 2474 _boostingLimitReached = true; 2475 } 2476 } 2477 else if(_statistics->precBoosts >= 2) 2478 { 2479 // general case: increase the number of decimal digits by 3/2, 2480 int newNbDigits = (int)floor(BP::default_precision() * realParam( 2481 SoPlexBase<R>::PRECISION_BOOSTING_FACTOR)); 2482 2483 if(intParam(SoPlexBase<R>::MULTIPRECISION_LIMIT) >= newNbDigits) 2484 { 2485 BP::default_precision(newNbDigits); 2486 } 2487 else 2488 { 2489 SPX_MSG_INFO1(spxout, spxout << "Maximum number of digits for the multiprecision type reached.\n" 2490 << "To increase this limit, modify the parameter multiprecision_limit.\n" 2491 << "Giving up.\n"); 2492 _boostingLimitReached = true; 2493 } 2494 } 2495 2496 _statistics->boostingStepTime->stop(); 2497 return !_boostingLimitReached; 2498 #else 2499 SPX_MSG_ERROR(std::cerr << 2500 "SoPlex was not compiled with MPFR, which is needed for precision changes \n"); 2501 return false; 2502 #endif 2503 } 2504 2505 /// increase the multiprecision, return false if maximum precision is reached, true otherwise 2506 template <class R> 2507 void SoPlexBase<R>::_resetBoostedPrecision() 2508 { 2509 _statistics->precBoosts = 0; 2510 #ifdef SOPLEX_WITH_MPFR 2511 BP::default_precision(50); 2512 #endif 2513 } 2514 2515 2516 2517 /// setup recovery mecanism using multiprecision, return false if maximum precision reached, true otherwise 2518 template <class R> 2519 bool SoPlexBase<R>::_setupBoostedSolverAfterRecovery() 2520 { 2521 assert(boolParam(SoPlexBase<R>::PRECISION_BOOSTING)); 2522 2523 _switchToBoosted(); 2524 2525 if(!_boostPrecision()) 2526 return false; 2527 2528 return true; 2529 } 2530 2531 2532 2533 /// return true if slack basis has to be loaded for boosted solver 2534 template <class R> 2535 bool SoPlexBase<R>::_isBoostedStartingFromSlack(bool initialSolve) 2536 { 2537 return (!boolParam(SoPlexBase<R>::BOOSTED_WARM_START) && initialSolve); 2538 } 2539 2540 2541 2542 // store last basis met in double precision solver for a potential iteration with increased precision 2543 template <class R> 2544 void SoPlexBase<R>::_storeBasisAsOldBasis(DataArray< typename SPxSolverBase<R>::VarStatus >& rows, 2545 DataArray< typename SPxSolverBase<R>::VarStatus >& cols) 2546 { 2547 if(_inStandardMode()) 2548 { 2549 SPX_MSG_INFO1(spxout, spxout << "Store basis as old basis (from solver)" << "\n"); 2550 _oldBasisStatusRows = rows; 2551 _oldBasisStatusCols = cols; 2552 _hasOldBasis = true; 2553 } 2554 else if(_inFeasMode()) 2555 { 2556 SPX_MSG_INFO1(spxout, spxout << "Store basis as old basis (from solver - testing feasibility)" << 2557 "\n"); 2558 _oldFeasBasisStatusRows = rows; 2559 _oldFeasBasisStatusCols = cols; 2560 _hasOldFeasBasis = true; 2561 } 2562 else if(_inUnbdMode()) 2563 { 2564 SPX_MSG_INFO1(spxout, spxout << "Store basis as old basis (from solver - testing unboundedness)" << 2565 "\n"); 2566 _oldUnbdBasisStatusRows = rows; 2567 _oldUnbdBasisStatusCols = cols; 2568 _hasOldUnbdBasis = true; 2569 } 2570 } 2571 2572 2573 2574 // store last basis met in boosted solver for a potential iteration with increased precision 2575 template <class R> 2576 void SoPlexBase<R>::_storeBasisAsOldBasisBoosted(DataArray< typename SPxSolverBase<BP>::VarStatus >& 2577 rows, DataArray< typename SPxSolverBase<BP>::VarStatus >& cols) 2578 { 2579 if(_inStandardMode()) 2580 { 2581 SPX_MSG_INFO1(spxout, spxout << "Store basis as old basis (from boosted solver)" << "\n"); 2582 _convertDataArrayVarStatusToRPrecision(rows, _oldBasisStatusRows); 2583 _convertDataArrayVarStatusToRPrecision(cols, _oldBasisStatusCols); 2584 _hasOldBasis = true; 2585 } 2586 else if(_inFeasMode()) 2587 { 2588 SPX_MSG_INFO1(spxout, spxout << 2589 "Store basis as old basis (from boosted solver - testing feasibility)" << "\n"); 2590 _convertDataArrayVarStatusToRPrecision(rows, _oldFeasBasisStatusRows); 2591 _convertDataArrayVarStatusToRPrecision(cols, _oldFeasBasisStatusCols); 2592 _hasOldFeasBasis = true; 2593 } 2594 else if(_inUnbdMode()) 2595 { 2596 SPX_MSG_INFO1(spxout, spxout << 2597 "Store basis as old basis (from boosted solver - testing unboundedness)" << "\n"); 2598 _convertDataArrayVarStatusToRPrecision(rows, _oldUnbdBasisStatusRows); 2599 _convertDataArrayVarStatusToRPrecision(cols, _oldUnbdBasisStatusCols); 2600 _hasOldUnbdBasis = true; 2601 } 2602 } 2603 2604 2605 2606 // load last basis if there is one available, else return false 2607 template <class R> 2608 bool SoPlexBase<R>::_loadBasisFromOldBasis(bool boosted) 2609 { 2610 if(!boosted) 2611 { 2612 if(_inStandardMode() && _hasOldBasis) 2613 { 2614 SPX_MSG_INFO1(spxout, spxout << "Load basis from old basis (in solver)" << "\n"); 2615 _solver.setBasis(_oldBasisStatusRows.get_const_ptr(), _oldBasisStatusCols.get_const_ptr()); 2616 } 2617 else if(_inFeasMode() && _hasOldFeasBasis) 2618 { 2619 SPX_MSG_INFO1(spxout, spxout << "Load basis from old basis (in solver - testing feasibility)" << 2620 "\n"); 2621 _solver.setBasis(_oldFeasBasisStatusRows.get_const_ptr(), _oldFeasBasisStatusCols.get_const_ptr()); 2622 } 2623 else if(_inUnbdMode() && _hasOldUnbdBasis) 2624 { 2625 SPX_MSG_INFO1(spxout, spxout << "Load basis from old basis (in solver - testing unboundedness)" << 2626 "\n"); 2627 _solver.setBasis(_oldUnbdBasisStatusRows.get_const_ptr(), _oldUnbdBasisStatusCols.get_const_ptr()); 2628 } 2629 else 2630 { 2631 SPX_MSG_INFO1(spxout, spxout << "No old basis available" << "\n"); 2632 return false; 2633 } 2634 } 2635 else 2636 { 2637 if(_inStandardMode() && _hasOldBasis) 2638 { 2639 SPX_MSG_INFO1(spxout, spxout << "Load basis from old basis (in boosted solver)" << "\n"); 2640 _convertDataArrayVarStatusToBoosted(_oldBasisStatusRows, _tmpBasisStatusRows); 2641 _convertDataArrayVarStatusToBoosted(_oldBasisStatusCols, _tmpBasisStatusCols); 2642 _boostedSolver.setBasis(_tmpBasisStatusRows.get_const_ptr(), _tmpBasisStatusCols.get_const_ptr()); 2643 } 2644 else if(_inFeasMode() && _hasOldFeasBasis) 2645 { 2646 SPX_MSG_INFO1(spxout, spxout << 2647 "Load basis from old basis (in boosted solver - testing feasibility)" << "\n"); 2648 _convertDataArrayVarStatusToBoosted(_oldFeasBasisStatusRows, _tmpBasisStatusRows); 2649 _convertDataArrayVarStatusToBoosted(_oldFeasBasisStatusCols, _tmpBasisStatusCols); 2650 _boostedSolver.setBasis(_tmpBasisStatusRows.get_const_ptr(), _tmpBasisStatusCols.get_const_ptr()); 2651 } 2652 else if(_inUnbdMode() && _hasOldUnbdBasis) 2653 { 2654 SPX_MSG_INFO1(spxout, spxout << 2655 "Load basis from old basis (in boosted solver - testing unboundedness)" << "\n"); 2656 _convertDataArrayVarStatusToBoosted(_oldUnbdBasisStatusRows, _tmpBasisStatusRows); 2657 _convertDataArrayVarStatusToBoosted(_oldUnbdBasisStatusCols, _tmpBasisStatusCols); 2658 _boostedSolver.setBasis(_tmpBasisStatusRows.get_const_ptr(), _tmpBasisStatusCols.get_const_ptr()); 2659 } 2660 else 2661 { 2662 SPX_MSG_INFO1(spxout, spxout << "No old basis available" << "\n"); 2663 return false; 2664 } 2665 } 2666 2667 return true; 2668 } 2669 2670 2671 2672 // update statistics for precision boosting 2673 template <class R> 2674 void SoPlexBase<R>::_updateBoostingStatistics() 2675 { 2676 // update statistics 2677 if(_statistics->iterations <= _prevIterations) 2678 { 2679 _lastStallPrecBoosts++; 2680 _statistics->stallPrecBoosts++; 2681 } 2682 else 2683 { 2684 _factorSolNewBasisPrecBoost = true; 2685 _lastStallPrecBoosts = 0; 2686 _statistics->pivotPrecBoosts = _statistics->precBoosts; 2687 } 2688 2689 if(_inFeasMode()) 2690 _statistics->feasPrecBoosts++; 2691 else if(_inUnbdMode()) 2692 _statistics->unbdPrecBoosts++; 2693 } 2694 2695 2696 2697 // indicate we are solving the original LP 2698 template <class R> 2699 void SoPlexBase<R>::_switchToStandardMode() 2700 { 2701 _certificateMode = 0; 2702 } 2703 2704 // indicate we are testing infeasibility 2705 template <class R> 2706 void SoPlexBase<R>::_switchToFeasMode() 2707 { 2708 _certificateMode = 1; 2709 } 2710 2711 // indicate we are testing unboundedness 2712 template <class R> 2713 void SoPlexBase<R>::_switchToUnbdMode() 2714 { 2715 _certificateMode = 2; 2716 } 2717 2718 // check we are solving the original LP 2719 template <class R> 2720 bool SoPlexBase<R>::_inStandardMode() 2721 { 2722 return _certificateMode == 0; 2723 } 2724 2725 // check we are testing infeasibility 2726 template <class R> 2727 bool SoPlexBase<R>::_inFeasMode() 2728 { 2729 return _certificateMode == 1; 2730 } 2731 2732 // check we are testing unboundedness 2733 template <class R> 2734 bool SoPlexBase<R>::_inUnbdMode() 2735 { 2736 return _certificateMode == 2; 2737 } 2738 2739 2740 2741 /// solves current problem using multiprecision floating-point solver 2742 /// return false if a new boosted iteration is necessary, true otherwise 2743 template <class R> 2744 void SoPlexBase<R>::_solveRealForRationalBoostedStable( 2745 SolRational& sol, 2746 bool& primalFeasible, 2747 bool& dualFeasible, 2748 bool& infeasible, 2749 bool& unbounded, 2750 bool& stoppedTime, 2751 bool& stoppedIter, 2752 bool& error, 2753 bool& needNewBoostedIt 2754 ) 2755 { 2756 assert(boolParam(SoPlexBase<R>::PRECISION_BOOSTING)); 2757 assert(!boolParam(SoPlexBase<R>::ITERATIVE_REFINEMENT)); 2758 2759 #ifdef SOPLEX_WITH_MPFR 2760 2761 // start rational solving timing 2762 _statistics->rationalTime->start(); 2763 2764 SPX_MSG_INFO1(spxout, spxout << "Current precision = 1e-" << BP::default_precision() << ", "); 2765 2766 primalFeasible = false; 2767 dualFeasible = false; 2768 infeasible = false; 2769 unbounded = false; 2770 stoppedTime = false; 2771 stoppedIter = false; 2772 error = false; 2773 needNewBoostedIt = false; 2774 2775 // initialize boosted solver 2776 _boostedSolver.init(); 2777 2778 _statistics->boostingStepTime->start(); 2779 2780 BP tolerance = boost::multiprecision::pow(BP(10), 2781 -(int)(BP::default_precision() * _tolPrecisionRatio)); 2782 2783 BP epsilonZero = boost::multiprecision::pow(BP(10), 2784 -(int)(BP::default_precision() * _epsZeroPrecisionRatio)); 2785 BP epsilonFactor = boost::multiprecision::pow(BP(10), 2786 -(int)(BP::default_precision() * _epsFactorPrecisionRatio)); 2787 BP epsilonUpdate = boost::multiprecision::pow(BP(10), 2788 -(int)(BP::default_precision() * _epsUpdatePrecisionRatio)); 2789 BP epsilonPivot = boost::multiprecision::pow(BP(10), 2790 -(int)(BP::default_precision() * _epsPivotPrecisionRatio)); 2791 2792 this->_tolerances->setEpsilon((Real) epsilonZero); 2793 this->_tolerances->setEpsilonFactorization((Real) epsilonFactor); 2794 this->_tolerances->setEpsilonUpdate((Real) epsilonUpdate); 2795 this->_tolerances->setEpsilonPivot((Real) epsilonPivot); 2796 2797 // set tolerances of the boosted solver 2798 if(boolParam(SoPlexBase<R>::ADAPT_TOLS_TO_MULTIPRECISION)) 2799 { 2800 this->_tolerances->setFloatingPointFeastol((Real) tolerance); 2801 this->_tolerances->setFloatingPointOpttol((Real) tolerance); 2802 } 2803 2804 _statistics->boostingStepTime->stop(); 2805 2806 typename SPxSolverBase<BP>::Status boostedResult = SPxSolverBase<BP>::UNKNOWN; 2807 2808 _modLower.reDim(numColsRational(), false); 2809 _modUpper.reDim(numColsRational(), false); 2810 _modLhs.reDim(numRowsRational(), false); 2811 _modRhs.reDim(numRowsRational(), false); 2812 _modObj.reDim(numColsRational(), false); 2813 2814 // declare real vectors after boosting precision 2815 VectorBase<BP> boostedPrimalReal(numColsRational()); 2816 VectorBase<BP> boostedDualReal(numRowsRational()); 2817 2818 Rational boundsViolation; 2819 Rational sideViolation; 2820 Rational redCostViolation; 2821 Rational dualViolation; 2822 2823 // solve original LP 2824 SPX_MSG_INFO1(spxout, spxout << "Initial floating-point boosted solve . . .\n"); 2825 2826 SPX_MSG_INFO1(spxout, spxout << "Boosted iteration " << _statistics->precBoosts << "\n"); 2827 2828 for(int r = numRowsRational() - 1; r >= 0; r--) 2829 assert(_boostedSolver.maxRowObj(r) == 0.0); 2830 2831 // solve original LP with boosted solver 2832 _statistics->rationalTime->stop(); 2833 _solveRealForRationalBoosted(boostedPrimalReal, boostedDualReal, _basisStatusRows, _basisStatusCols, 2834 boostedResult, true); 2835 2836 // special case on the first boosted iteration 2837 // otherwise attributs such as _factorSolNewBasisPrecBoost are not updated correctly 2838 if(_statistics->precBoosts == 1) 2839 _updateBoostingStatistics(); 2840 2841 // evalute result 2842 if(_evaluateResult<BP>(_boostedSolver, boostedResult, false, sol, boostedDualReal, infeasible, 2843 unbounded, stoppedTime, stoppedIter, 2844 error)) 2845 return; 2846 2847 // stop rational solving time 2848 _statistics->rationalTime->start(); 2849 2850 int dualSize = 0; 2851 2852 // stores floating-point solution of original LP as current rational solution 2853 // ensure that solution vectors have right dimension; ensure that solution is aligned with basis 2854 _storeRealSolutionAsRational<BP>(sol, boostedPrimalReal, boostedDualReal, dualSize); 2855 2856 // control progress 2857 Rational maxViolation; 2858 Rational bestViolation = _rationalPosInfty; 2859 const Rational violationImprovementFactor = 16; 2860 const Rational errorCorrectionFactor = 1.1; 2861 Rational errorCorrection = 2; 2862 2863 const bool maximizing = (intParam(SoPlexBase<R>::OBJSENSE) == SoPlexBase<R>::OBJSENSE_MAXIMIZE); 2864 2865 // used in _ratrecAndOrRatfac to order a break or a continue outside of the function. 2866 bool breakAfter; 2867 bool continueAfter; 2868 2869 int minIRRoundsRemaining = 0; 2870 2871 do 2872 { 2873 SPxOut::debug(this, "Computing primal violations.\n"); 2874 2875 // computes violation of bounds 2876 _computeBoundsViolation(sol, boundsViolation); 2877 2878 // computes violation of sides 2879 _computeSidesViolation(sol, sideViolation); 2880 2881 SPxOut::debug(this, "Computing dual violations.\n"); 2882 2883 // compute reduced cost violation 2884 _computeReducedCostViolation(sol, redCostViolation, maximizing); 2885 2886 // compute dual violation 2887 _computeDualViolation(sol, dualViolation, maximizing); 2888 2889 _modObj = sol._redCost; 2890 2891 // output violations; the reduced cost violations for artificially introduced slack columns are actually violations of the dual multipliers 2892 SPX_MSG_INFO1(spxout, spxout 2893 << "Max. bound violation = " << boundsViolation.str() << ", ofm = 1e" << orderOfMagnitude( 2894 boundsViolation) << "\n" 2895 << "Max. row violation = " << sideViolation.str() << ", ofm = 1e" << orderOfMagnitude( 2896 sideViolation) << "\n" 2897 << "Max. reduced cost violation = " << redCostViolation.str() << ", ofm = 1e" << orderOfMagnitude( 2898 redCostViolation) << "\n" 2899 << "Max. dual violation = " << dualViolation.str() << ", ofm = 1e" << orderOfMagnitude( 2900 dualViolation) << "\n"); 2901 2902 SPxOut::debug(this, "Progress table: {} & {} & {} & {} & {} & {}\n", _statistics->refinements, 2903 _statistics->iterations, _statistics->solvingTime->time(), _statistics->rationalTime->time(), 2904 boundsViolation > sideViolation ? boundsViolation : sideViolation, 2905 redCostViolation > dualViolation ? redCostViolation : dualViolation); 2906 2907 // terminate if tolerances are satisfied 2908 primalFeasible = (boundsViolation <= _rationalFeastol && sideViolation <= _rationalFeastol); 2909 dualFeasible = (redCostViolation <= _rationalOpttol && dualViolation <= _rationalOpttol); 2910 2911 if(primalFeasible && dualFeasible) 2912 { 2913 SPX_MSG_INFO1(spxout, spxout << "Tolerances reached.\n"); 2914 return; 2915 } 2916 2917 // terminate if some limit is reached 2918 if(_isSolveStopped(stoppedTime, stoppedIter)) 2919 return; 2920 2921 // compute maximal violation 2922 maxViolation = boundsViolation; 2923 2924 if(sideViolation > maxViolation) 2925 maxViolation = sideViolation; 2926 2927 if(redCostViolation > maxViolation) 2928 maxViolation = redCostViolation; 2929 2930 if(dualViolation > maxViolation) 2931 maxViolation = dualViolation; 2932 2933 2934 // perform rational reconstruction and/or factorization 2935 _ratrecAndOrRatfac(minIRRoundsRemaining, _lastStallPrecBoosts, _statistics->precBoosts, 2936 _factorSolNewBasisPrecBoost, _nextRatrecPrecBoost, 2937 errorCorrectionFactor, errorCorrection, maxViolation, 2938 sol, primalFeasible, dualFeasible, 2939 stoppedTime, stoppedIter, error, breakAfter, continueAfter); 2940 2941 if(!continueAfter) 2942 break; 2943 2944 } 2945 while(true); 2946 2947 if(primalFeasible && dualFeasible) 2948 { 2949 needNewBoostedIt = false; 2950 2951 // correct basis status for restricted inequalities 2952 if(_hasBasis) 2953 { 2954 for(int r = numRowsRational() - 1; r >= 0; r--) 2955 { 2956 assert((lhsRational(r) == rhsRational(r)) == (_rowTypes[r] == RANGETYPE_FIXED)); 2957 2958 if(_rowTypes[r] != RANGETYPE_FIXED && _basisStatusRows[r] == SPxSolverBase<R>::FIXED) 2959 _basisStatusRows[r] = (maximizing == (sol._dual[r] < 0)) 2960 ? SPxSolverBase<R>::ON_LOWER 2961 : SPxSolverBase<R>::ON_UPPER; 2962 } 2963 } 2964 2965 // compute objective function values 2966 assert(sol._isPrimalFeasible == sol._isDualFeasible); 2967 2968 if(sol._isPrimalFeasible) 2969 { 2970 sol._objVal = sol._primal * _rationalLP->maxObj(); 2971 2972 if(intParam(SoPlexBase<R>::OBJSENSE) == SoPlexBase<R>::OBJSENSE_MINIMIZE) 2973 sol._objVal *= -1; 2974 } 2975 2976 // set objective coefficients for all rows to zero 2977 _boostedSolver.clearRowObjs(); 2978 2979 // stop rational solving time 2980 _statistics->rationalTime->stop(); 2981 } 2982 else 2983 { 2984 SPX_MSG_INFO1(spxout, spxout << "\nNo success. Launch new boosted solve . . .\n"); 2985 needNewBoostedIt = true; 2986 2987 // stop rational solving time 2988 _statistics->rationalTime->stop(); 2989 } 2990 2991 #else 2992 SPX_MSG_ERROR(std::cerr << 2993 "SoPlex was not compiled with MPFR, which is needed for precision changes \n"); 2994 #endif 2995 } 2996 2997 2998 2999 /// solves current problem with iterative refinement and recovery mechanism using boosted solver 3000 /// return false if a new boosted iteration is necessary, true otherwise 3001 template <class R> 3002 void SoPlexBase<R>::_performOptIRStableBoosted( 3003 SolRational& sol, 3004 bool acceptUnbounded, 3005 bool acceptInfeasible, 3006 int minIRRoundsRemaining, 3007 bool& primalFeasible, 3008 bool& dualFeasible, 3009 bool& infeasible, 3010 bool& unbounded, 3011 bool& stoppedTime, 3012 bool& stoppedIter, 3013 bool& error, 3014 bool& needNewBoostedIt 3015 ) 3016 { 3017 assert(boolParam(SoPlexBase<R>::ITERATIVE_REFINEMENT)); 3018 assert(boolParam(SoPlexBase<R>::PRECISION_BOOSTING)); 3019 3020 #ifdef SOPLEX_WITH_MPFR 3021 3022 // start rational solving timing 3023 _statistics->rationalTime->start(); 3024 3025 SPX_MSG_INFO1(spxout, spxout << "Current precision = 1e-" << BP::default_precision() << ", "); 3026 3027 typename SPxSolverBase<BP>::Status boostedResult = SPxSolverBase<BP>::UNKNOWN; 3028 3029 _modLower.reDim(numColsRational(), false); 3030 _modUpper.reDim(numColsRational(), false); 3031 _modLhs.reDim(numRowsRational(), false); 3032 _modRhs.reDim(numRowsRational(), false); 3033 _modObj.reDim(numColsRational(), false); 3034 3035 primalFeasible = false; 3036 dualFeasible = false; 3037 infeasible = false; 3038 unbounded = false; 3039 stoppedTime = false; 3040 stoppedIter = false; 3041 error = false; 3042 needNewBoostedIt = false; 3043 3044 // initialize boosted solver 3045 _boostedSolver.init(); 3046 3047 _statistics->boostingStepTime->start(); 3048 3049 BP tolerance = boost::multiprecision::pow(BP(10), 3050 -(int)(BP::default_precision() * _tolPrecisionRatio)); 3051 3052 BP epsilonZero = boost::multiprecision::pow(BP(10), 3053 -(int)(BP::default_precision() * _epsZeroPrecisionRatio)); 3054 BP epsilonFactor = boost::multiprecision::pow(BP(10), 3055 -(int)(BP::default_precision() * _epsFactorPrecisionRatio)); 3056 BP epsilonUpdate = boost::multiprecision::pow(BP(10), 3057 -(int)(BP::default_precision() * _epsUpdatePrecisionRatio)); 3058 BP epsilonPivot = boost::multiprecision::pow(BP(10), 3059 -(int)(BP::default_precision() * _epsPivotPrecisionRatio)); 3060 3061 this->_tolerances->setEpsilon((Real) epsilonZero); 3062 this->_tolerances->setEpsilonFactorization((Real) epsilonFactor); 3063 this->_tolerances->setEpsilonUpdate((Real) epsilonUpdate); 3064 this->_tolerances->setEpsilonPivot((Real) epsilonPivot); 3065 3066 // set tolerances of the boosted solver 3067 if(boolParam(SoPlexBase<R>::ADAPT_TOLS_TO_MULTIPRECISION)) 3068 { 3069 this->_tolerances->setFloatingPointFeastol((Real) tolerance); 3070 this->_tolerances->setFloatingPointOpttol((Real) tolerance); 3071 } 3072 3073 _statistics->boostingStepTime->stop(); 3074 3075 // declare real vectors after boosting precision 3076 VectorBase<BP> boostedPrimalReal(numColsRational()); 3077 VectorBase<BP> boostedDualReal(numRowsRational()); 3078 3079 Rational boundsViolation; 3080 Rational sideViolation; 3081 Rational redCostViolation; 3082 Rational dualViolation; 3083 Rational primalScale; 3084 Rational dualScale; 3085 Rational maxScale; 3086 3087 // solve original LP 3088 SPX_MSG_INFO1(spxout, spxout << "Initial floating-point boosted solve . . .\n"); 3089 3090 SPX_MSG_INFO1(spxout, spxout << "Boosted iteration " << _statistics->precBoosts << "\n"); 3091 3092 for(int r = numRowsRational() - 1; r >= 0; r--) 3093 assert(_boostedSolver.maxRowObj(r) == 0.0); 3094 3095 // solve original LP with boosted solver 3096 _statistics->rationalTime->stop(); 3097 _solveRealForRationalBoosted(boostedPrimalReal, boostedDualReal, _basisStatusRows, _basisStatusCols, 3098 boostedResult, true); 3099 3100 // evalute result 3101 if(_evaluateResult<BP>(_boostedSolver, boostedResult, false, sol, boostedDualReal, infeasible, 3102 unbounded, stoppedTime, stoppedIter, 3103 error)) 3104 return; 3105 3106 _statistics->rationalTime->start(); 3107 3108 int dualSize = 0; 3109 3110 // stores floating-point solution of original LP as current rational solution 3111 // ensure that solution vectors have right dimension; ensure that solution is aligned with basis 3112 _storeRealSolutionAsRational<BP>(sol, boostedPrimalReal, boostedDualReal, dualSize); 3113 3114 // initial scaling factors are one 3115 primalScale = _rationalPosone; 3116 dualScale = _rationalPosone; 3117 3118 // control progress 3119 Rational maxViolation; 3120 Rational bestViolation = _rationalPosInfty; 3121 const Rational violationImprovementFactor = 16; 3122 const Rational errorCorrectionFactor = 1.1; 3123 Rational errorCorrection = 2; 3124 int numFailedRefinements = 0; 3125 3126 // refinement loop 3127 const bool maximizing = (intParam(SoPlexBase<R>::OBJSENSE) == SoPlexBase<R>::OBJSENSE_MAXIMIZE); 3128 const int maxDimRational = numColsRational() > numRowsRational() ? numColsRational() : 3129 numRowsRational(); 3130 bool factorSolNewBasis = true; 3131 int lastStallRefinements = 0; 3132 int nextRatrecRefinement = 0; 3133 3134 // used in _ratrecAndOrRatfac to order a break or a continue outside of the function. 3135 bool breakAfter; 3136 bool continueAfter; 3137 3138 do 3139 { 3140 // decrement minIRRoundsRemaining counter 3141 minIRRoundsRemaining--; 3142 3143 SPxOut::debug(this, "Computing primal violations.\n"); 3144 3145 // computes violation of bounds 3146 _computeBoundsViolation(sol, boundsViolation); 3147 3148 // computes violation of sides 3149 _computeSidesViolation(sol, sideViolation); 3150 3151 SPxOut::debug(this, "Computing dual violations.\n"); 3152 3153 // compute reduced cost violation 3154 _computeReducedCostViolation(sol, redCostViolation, maximizing); 3155 3156 // compute dual violation 3157 _computeDualViolation(sol, dualViolation, maximizing); 3158 3159 _modObj = sol._redCost; 3160 3161 // output violations; the reduced cost violations for artificially introduced slack columns are actually violations of the dual multipliers 3162 SPX_MSG_INFO1(spxout, spxout 3163 << "Max. bound violation = " << boundsViolation.str() << ", ofm = 1e" << orderOfMagnitude( 3164 boundsViolation) << "\n" 3165 << "Max. row violation = " << sideViolation.str() << ", ofm = 1e" << orderOfMagnitude( 3166 sideViolation) << "\n" 3167 << "Max. reduced cost violation = " << redCostViolation.str() << ", ofm = 1e" << orderOfMagnitude( 3168 redCostViolation) << "\n" 3169 << "Max. dual violation = " << dualViolation.str() << ", ofm = 1e" << orderOfMagnitude( 3170 dualViolation) << "\n"); 3171 3172 SPxOut::debug(this, "Progress table: {} & {} & {} & {} & {} & {}\n", _statistics->refinements, 3173 _statistics->iterations, _statistics->solvingTime->time(), _statistics->rationalTime->time(), 3174 boundsViolation > sideViolation ? boundsViolation : sideViolation, 3175 redCostViolation > dualViolation ? redCostViolation : dualViolation); 3176 3177 // check termination criteria for refinement loop 3178 if(_isRefinementOver(primalFeasible, dualFeasible, boundsViolation, sideViolation, redCostViolation, 3179 dualViolation, minIRRoundsRemaining, stoppedTime, stoppedIter, numFailedRefinements)) 3180 break; 3181 3182 // check refinement progress 3183 _checkRefinementProgress(boundsViolation, sideViolation, redCostViolation, dualViolation, 3184 maxViolation, bestViolation, violationImprovementFactor, numFailedRefinements); 3185 3186 // perform rational reconstruction and/or factorization 3187 _ratrecAndOrRatfac(minIRRoundsRemaining, lastStallRefinements, _statistics->refinements, 3188 factorSolNewBasis, 3189 nextRatrecRefinement, 3190 errorCorrectionFactor, errorCorrection, maxViolation, 3191 sol, primalFeasible, dualFeasible, 3192 stoppedTime, stoppedIter, error, breakAfter, continueAfter); 3193 3194 if(breakAfter) 3195 break; 3196 else if(continueAfter) 3197 continue; 3198 3199 // start refinement 3200 3201 // compute primal scaling factor; limit increase in scaling by tolerance used in floating point solve 3202 _computePrimalScalingFactor(maxScale, primalScale, boundsViolation, sideViolation, 3203 redCostViolation); 3204 3205 // apply scaled bounds and scaled sides 3206 _applyScaledBounds<BP>(_boostedSolver, primalScale); 3207 _applyScaledSides<BP>(_boostedSolver, primalScale); 3208 3209 // compute dual scaling factor; limit increase in scaling by tolerance used in floating point solve 3210 _computeDualScalingFactor(maxScale, primalScale, dualScale, redCostViolation, dualViolation); 3211 3212 // apply scaled objective function 3213 _applyScaledObj<BP>(_boostedSolver, dualScale, sol); 3214 3215 SPX_MSG_INFO1(spxout, spxout << "Refined floating-point solve . . .\n"); 3216 3217 // ensure that artificial slack columns are basic and inequality constraints are nonbasic; otherwise we may end 3218 // up with dual violation on inequality constraints after removing the slack columns; do not change this in the 3219 // floating-point solver, though, because the solver may require its original basis to detect optimality 3220 if(_slackCols.num() > 0 && _hasBasis) 3221 { 3222 int numOrigCols = numColsRational() - _slackCols.num(); 3223 assert(_slackCols.num() <= 0 || boolParam(SoPlexBase<R>::EQTRANS)); 3224 3225 for(int i = 0; i < _slackCols.num(); i++) 3226 { 3227 int row = _slackCols.colVector(i).index(0); 3228 int col = numOrigCols + i; 3229 3230 assert(row >= 0); 3231 assert(row < numRowsRational()); 3232 3233 if(_basisStatusRows[row] == SPxSolverBase<R>::BASIC 3234 && _basisStatusCols[col] != SPxSolverBase<R>::BASIC) 3235 { 3236 _basisStatusRows[row] = _basisStatusCols[col]; 3237 _basisStatusCols[col] = SPxSolverBase<R>::BASIC; 3238 _rationalLUSolver.clear(); 3239 } 3240 } 3241 } 3242 3243 // load basis 3244 if(_hasBasis && _boostedSolver.basis().status() < SPxBasisBase<BP>::REGULAR) 3245 { 3246 SPxOut::debug(this, "basis (status = {}) desc before set:\n", _boostedSolver.basis().status()); 3247 #ifdef SOPLEX_DEBUG 3248 _boostedSolver.basis().desc().dump() 3249 #endif 3250 _boostedSolver.setBasis(_tmpBasisStatusRows.get_const_ptr(), _tmpBasisStatusCols.get_const_ptr()); 3251 SPxOut::debug(this, "basis (status = {}) desc after set:\n", _boostedSolver.basis().status()); 3252 #ifdef SOPLEX_DEBUG 3253 _boostedSolver.basis().desc().dump() 3254 #endif 3255 _hasBasis = _boostedSolver.basis().status() > SPxBasisBase<BP>::NO_PROBLEM; 3256 SPxOut::debug(this, "setting basis in solver {} (3)\n", (_hasBasis ? "successful" : "failed")); 3257 } 3258 3259 // solve modified problem 3260 int prevIterations = _statistics->iterations; 3261 _statistics->rationalTime->stop(); 3262 // turn off simplifier if scaling factors are too high 3263 int simplifier = intParam(SoPlexBase<R>::SIMPLIFIER); 3264 3265 if(primalScale > 1e20 || dualScale > 1e20) 3266 setIntParam(SoPlexBase<R>::SIMPLIFIER, SoPlexBase<R>::SIMPLIFIER_OFF); 3267 3268 _solveRealForRationalBoosted(boostedPrimalReal, boostedDualReal, _basisStatusRows, _basisStatusCols, 3269 boostedResult, false); 3270 3271 // reset simplifier param to previous value 3272 setIntParam(SoPlexBase<R>::SIMPLIFIER, simplifier); 3273 3274 // count refinements and remember whether we moved to a new basis 3275 _statistics->refinements++; 3276 3277 if(_statistics->iterations <= prevIterations) 3278 { 3279 lastStallRefinements++; 3280 _statistics->stallRefinements++; 3281 } 3282 else 3283 { 3284 factorSolNewBasis = true; 3285 lastStallRefinements = 0; 3286 _statistics->pivotRefinements = _statistics->refinements; 3287 } 3288 3289 // evaluate result; if modified problem was not solved to optimality, stop refinement; 3290 if(_evaluateResult<BP>(_boostedSolver, boostedResult, true, sol, boostedDualReal, infeasible, 3291 unbounded, stoppedTime, stoppedIter, 3292 error)) 3293 return; 3294 3295 _statistics->rationalTime->start(); 3296 3297 int primalSize; 3298 3299 // correct primal solution and align with basis 3300 _correctPrimalSolution<BP>(sol, primalScale, primalSize, maxDimRational, boostedPrimalReal); 3301 3302 // update or recompute slacks depending on which looks faster 3303 _updateSlacks(sol, primalSize); 3304 3305 const int numCorrectedPrimals = _primalDualDiff.size(); 3306 3307 // correct dual solution and align with basis 3308 _correctDualSolution<BP>(_boostedSolver, sol, maximizing, boostedDualReal, dualScale, dualSize, 3309 maxDimRational); 3310 3311 // update or recompute reduced cost values depending on which looks faster 3312 // adding one to the length of the dual vector accounts for the objective function vector 3313 _updateReducedCosts(sol, dualSize, numCorrectedPrimals); 3314 } 3315 while(true); 3316 3317 if(primalFeasible && dualFeasible) 3318 { 3319 needNewBoostedIt = false; 3320 3321 // correct basis status for restricted inequalities 3322 if(_hasBasis) 3323 { 3324 for(int r = numRowsRational() - 1; r >= 0; r--) 3325 { 3326 assert((lhsRational(r) == rhsRational(r)) == (_rowTypes[r] == RANGETYPE_FIXED)); 3327 3328 if(_rowTypes[r] != RANGETYPE_FIXED && _basisStatusRows[r] == SPxSolverBase<R>::FIXED) 3329 _basisStatusRows[r] = (maximizing == (sol._dual[r] < 0)) 3330 ? SPxSolverBase<R>::ON_LOWER 3331 : SPxSolverBase<R>::ON_UPPER; 3332 } 3333 } 3334 3335 // compute objective function values 3336 assert(sol._isPrimalFeasible == sol._isDualFeasible); 3337 3338 if(sol._isPrimalFeasible) 3339 { 3340 sol._objVal = sol._primal * _rationalLP->maxObj(); 3341 3342 if(intParam(SoPlexBase<R>::OBJSENSE) == SoPlexBase<R>::OBJSENSE_MINIMIZE) 3343 sol._objVal *= -1; 3344 } 3345 3346 // set objective coefficients for all rows to zero 3347 _boostedSolver.clearRowObjs(); 3348 3349 // stop rational solving time 3350 _statistics->rationalTime->stop(); 3351 } 3352 else 3353 { 3354 SPX_MSG_INFO1(spxout, spxout << "\nNo success. Launch new boosted solve . . .\n"); 3355 needNewBoostedIt = true; 3356 3357 // stop rational solving time 3358 _statistics->rationalTime->stop(); 3359 } 3360 3361 #else 3362 SPX_MSG_ERROR(std::cerr << 3363 "SoPlex was not compiled with MPFR, which is needed for precision changes \n"); 3364 #endif 3365 } 3366 3367 3368 /// performs iterative refinement on the auxiliary problem for testing unboundedness 3369 template <class R> 3370 void SoPlexBase<R>::_performUnboundedIRStable( 3371 SolRational& sol, 3372 bool& hasUnboundedRay, 3373 bool& stoppedTime, 3374 bool& stoppedIter, 3375 bool& error) 3376 { 3377 bool primalFeasible; 3378 bool dualFeasible; 3379 bool infeasible; 3380 bool unbounded; 3381 3382 // move objective function to constraints and adjust sides and bounds 3383 _transformUnbounded(); 3384 3385 // invalidate solution 3386 sol.invalidate(); 3387 3388 // remember current number of refinements and precision boosts 3389 int oldRefinements = _statistics->refinements; 3390 3391 // perform iterative refinement 3392 _performOptIRWrapper(sol, false, false, 0, primalFeasible, dualFeasible, infeasible, unbounded, 3393 stoppedTime, stoppedIter, error); 3394 3395 // update unbounded refinement counter and unbounded precision boosts 3396 _statistics->unbdRefinements += _statistics->refinements - oldRefinements; 3397 3398 // stopped due to some limit 3399 if(stoppedTime || stoppedIter) 3400 { 3401 sol.invalidate(); 3402 hasUnboundedRay = false; 3403 error = false; 3404 } 3405 // the unbounded problem should always be solved to optimality 3406 else if(error || unbounded || infeasible || !primalFeasible || !dualFeasible) 3407 { 3408 sol.invalidate(); 3409 hasUnboundedRay = false; 3410 error = true; 3411 } 3412 else 3413 { 3414 const Rational& tau = sol._primal[numColsRational() - 1]; 3415 3416 SPxOut::debug(this, "tau = {} (roughly )\n", tau, tau.str()); 3417 3418 assert(tau <= 1.0 + 2.0 * realParam(SoPlexBase<R>::FEASTOL)); 3419 assert(tau >= -realParam(SoPlexBase<R>::FEASTOL)); 3420 3421 // because the right-hand side and all bounds (but tau's upper bound) are zero, tau should be approximately 3422 // zero if basic; otherwise at its upper bound 1 3423 error = !(tau >= _rationalPosone || tau <= _rationalFeastol); 3424 assert(!error); 3425 3426 hasUnboundedRay = (tau >= 1); 3427 } 3428 3429 // restore problem 3430 _untransformUnbounded(sol, hasUnboundedRay); 3431 } 3432 3433 3434 3435 /// performs iterative refinement on the auxiliary problem for testing feasibility 3436 template <class R> 3437 void SoPlexBase<R>::_performFeasIRStable( 3438 SolRational& sol, 3439 bool& withDualFarkas, 3440 bool& stoppedTime, 3441 bool& stoppedIter, 3442 bool& error) 3443 { 3444 bool primalFeasible; 3445 bool dualFeasible; 3446 bool infeasible; 3447 bool unbounded; 3448 bool success = false; 3449 error = false; 3450 3451 #if 0 3452 // if the problem has been found to be infeasible and an approximate Farkas proof is available, we compute a 3453 // scaled unit box around the origin that provably contains no feasible solution; this currently only works for 3454 // equality form 3455 ///@todo check whether approximate Farkas proof can be used 3456 _computeInfeasBox(_solRational, false); 3457 ///@todo if approx Farkas proof is good enough then exit without doing any transformation 3458 #endif 3459 3460 // remove objective function, shift, homogenize 3461 _transformFeasibility(); 3462 3463 // invalidate solution 3464 sol.invalidate(); 3465 3466 do 3467 { 3468 // remember current number of refinements and precision boosts 3469 int oldRefinements = _statistics->refinements; 3470 3471 // perform iterative refinement 3472 _performOptIRWrapper(sol, false, false, 0, primalFeasible, dualFeasible, infeasible, unbounded, 3473 stoppedTime, stoppedIter, error); 3474 3475 // update feasible refinement counter and precision boosts counter 3476 _statistics->feasRefinements += _statistics->refinements - oldRefinements; 3477 3478 // stopped due to some limit 3479 if(stoppedTime || stoppedIter) 3480 { 3481 sol.invalidate(); 3482 withDualFarkas = false; 3483 error = false; 3484 } 3485 // the feasibility problem should always be solved to optimality 3486 else if(error || unbounded || infeasible || !primalFeasible || !dualFeasible) 3487 { 3488 sol.invalidate(); 3489 withDualFarkas = false; 3490 error = true; 3491 } 3492 // else we should have either a refined Farkas proof or an approximate feasible solution to the original 3493 else 3494 { 3495 const Rational& tau = sol._primal[numColsRational() - 1]; 3496 3497 SPxOut::debug(this, "tau = {} (roughly )\n", tau, tau.str()); 3498 3499 assert(tau >= -realParam(SoPlexBase<R>::FEASTOL)); 3500 assert(tau <= 1.0 + realParam(SoPlexBase<R>::FEASTOL)); 3501 3502 error = (tau < -_rationalFeastol || tau > _rationalPosone + _rationalFeastol); 3503 withDualFarkas = (tau < _rationalPosone); 3504 3505 if(withDualFarkas) 3506 { 3507 _solRational._hasDualFarkas = true; 3508 _solRational._dualFarkas = _solRational._dual; 3509 3510 #if 0 3511 // check if we can compute sufficiently large Farkas box 3512 _computeInfeasBox(_solRational, true); 3513 #endif 3514 3515 if(true) //@todo check if computeInfeasBox found a sufficient box 3516 { 3517 3518 success = true; 3519 sol._isPrimalFeasible = false; 3520 } 3521 } 3522 else 3523 { 3524 sol._isDualFeasible = false; 3525 success = true; //successfully found approximate feasible solution 3526 } 3527 } 3528 } 3529 while(!error && !success && !(stoppedTime || stoppedIter)); 3530 3531 // restore problem 3532 _untransformFeasibility(sol, withDualFarkas); 3533 } 3534 3535 3536 3537 /// reduces matrix coefficient in absolute value by the lifting procedure of Thiele et al. 2013 3538 template <class R> 3539 void SoPlexBase<R>::_lift() 3540 { 3541 SPxOut::debug(this, "Reducing matrix coefficients by lifting.\n"); 3542 3543 // start timing 3544 _statistics->transformTime->start(); 3545 3546 SPX_DEBUG(_realLP->writeFileLPBase("beforeLift.lp", 0, 0, 0)); 3547 3548 // remember unlifted state 3549 _beforeLiftCols = numColsRational(); 3550 _beforeLiftRows = numRowsRational(); 3551 3552 // allocate vector memory 3553 DSVectorRational colVector; 3554 SVectorRational::Element liftingRowMem[2]; 3555 SVectorRational liftingRowVector(2, liftingRowMem); 3556 3557 // search each column for large nonzeros entries 3558 // @todo: rethink about the static_cast TODO 3559 const Rational maxValue = static_cast<Rational>(realParam(SoPlexBase<R>::LIFTMAXVAL)); 3560 3561 for(int i = 0; i < numColsRational(); i++) 3562 { 3563 SPxOut::debug(this, "in lifting: examining column {}\n", i); 3564 3565 // get column vector 3566 colVector = colVectorRational(i); 3567 3568 bool addedLiftingRow = false; 3569 int liftingColumnIndex = -1; 3570 3571 // go through nonzero entries of the column 3572 for(int k = colVector.size() - 1; k >= 0; k--) 3573 { 3574 const Rational& value = colVector.value(k); 3575 3576 if(spxAbs(value) > maxValue) 3577 { 3578 SPxOut::debug(this, " --> nonzero {} has value {} in row {}\n", k, value.str(), 3579 colVector.index(k)); 3580 3581 // add new column equal to maxValue times original column 3582 if(!addedLiftingRow) 3583 { 3584 SPxOut::debug(this, " --> adding lifting row\n"); 3585 3586 assert(liftingRowVector.size() == 0); 3587 3588 liftingColumnIndex = numColsRational(); 3589 liftingRowVector.add(i, maxValue); 3590 liftingRowVector.add(liftingColumnIndex, -1); 3591 3592 _rationalLP->addRow(LPRowRational(0, liftingRowVector, 0)); 3593 _realLP->addRow(LPRowBase<R>(0.0, DSVectorBase<R>(liftingRowVector), 0.0)); 3594 3595 assert(liftingColumnIndex == numColsRational() - 1); 3596 assert(liftingColumnIndex == numCols() - 1); 3597 3598 _rationalLP->changeBounds(liftingColumnIndex, _rationalNegInfty, _rationalPosInfty); 3599 _realLP->changeBounds(liftingColumnIndex, -realParam(SoPlexBase<R>::INFTY), 3600 realParam(SoPlexBase<R>::INFTY)); 3601 3602 liftingRowVector.clear(); 3603 addedLiftingRow = true; 3604 } 3605 3606 // get row index 3607 int rowIndex = colVector.index(k); 3608 assert(rowIndex >= 0); 3609 assert(rowIndex < _beforeLiftRows); 3610 assert(liftingColumnIndex == numColsRational() - 1); 3611 3612 SPxOut::debug(this, " --> changing matrix\n"); 3613 3614 // remove nonzero from original column 3615 _rationalLP->changeElement(rowIndex, i, 0); 3616 _realLP->changeElement(rowIndex, i, 0.0); 3617 3618 // add nonzero divided by maxValue to new column 3619 Rational newValue(value); 3620 newValue /= maxValue; 3621 _rationalLP->changeElement(rowIndex, liftingColumnIndex, newValue); 3622 _realLP->changeElement(rowIndex, liftingColumnIndex, R(newValue)); 3623 } 3624 } 3625 } 3626 3627 // search each column for small nonzeros entries 3628 const Rational minValue = Rational(realParam(SoPlexBase<R>::LIFTMINVAL)); 3629 3630 for(int i = 0; i < numColsRational(); i++) 3631 { 3632 SPxOut::debug(this, "in lifting: examining column {}\n", i); 3633 3634 // get column vector 3635 colVector = colVectorRational(i); 3636 3637 bool addedLiftingRow = false; 3638 int liftingColumnIndex = -1; 3639 3640 // go through nonzero entries of the column 3641 for(int k = colVector.size() - 1; k >= 0; k--) 3642 { 3643 const Rational& value = colVector.value(k); 3644 3645 if(spxAbs(value) < minValue) 3646 { 3647 SPxOut::debug(this, " --> nonzero {} has value {} in row {}\n", k, value.str(), 3648 colVector.index(k)); 3649 3650 // add new column equal to maxValue times original column 3651 if(!addedLiftingRow) 3652 { 3653 SPxOut::debug(this, " --> adding lifting row\n"); 3654 3655 assert(liftingRowVector.size() == 0); 3656 3657 liftingColumnIndex = numColsRational(); 3658 liftingRowVector.add(i, minValue); 3659 liftingRowVector.add(liftingColumnIndex, -1); 3660 3661 _rationalLP->addRow(LPRowRational(0, liftingRowVector, 0)); 3662 _realLP->addRow(LPRowBase<R>(0.0, DSVectorBase<R>(liftingRowVector), 0.0)); 3663 3664 assert(liftingColumnIndex == numColsRational() - 1); 3665 assert(liftingColumnIndex == numCols() - 1); 3666 3667 _rationalLP->changeBounds(liftingColumnIndex, _rationalNegInfty, _rationalPosInfty); 3668 _realLP->changeBounds(liftingColumnIndex, -realParam(SoPlexBase<R>::INFTY), 3669 realParam(SoPlexBase<R>::INFTY)); 3670 3671 liftingRowVector.clear(); 3672 addedLiftingRow = true; 3673 } 3674 3675 // get row index 3676 int rowIndex = colVector.index(k); 3677 assert(rowIndex >= 0); 3678 assert(rowIndex < _beforeLiftRows); 3679 assert(liftingColumnIndex == numColsRational() - 1); 3680 3681 SPxOut::debug(this, " --> changing matrix\n"); 3682 3683 // remove nonzero from original column 3684 _rationalLP->changeElement(rowIndex, i, 0); 3685 _realLP->changeElement(rowIndex, i, 0.0); 3686 3687 // add nonzero divided by maxValue to new column 3688 Rational newValue(value); 3689 newValue /= minValue; 3690 _rationalLP->changeElement(rowIndex, liftingColumnIndex, newValue); 3691 _realLP->changeElement(rowIndex, liftingColumnIndex, R(newValue)); 3692 } 3693 } 3694 } 3695 3696 // adjust basis 3697 if(_hasBasis) 3698 { 3699 assert(numColsRational() >= _beforeLiftCols); 3700 assert(numRowsRational() >= _beforeLiftRows); 3701 3702 _basisStatusCols.append(numColsRational() - _beforeLiftCols, SPxSolverBase<R>::BASIC); 3703 _basisStatusRows.append(numRowsRational() - _beforeLiftRows, SPxSolverBase<R>::FIXED); 3704 _rationalLUSolver.clear(); 3705 } 3706 3707 SPX_DEBUG(_realLP->writeFileLPBase("afterLift.lp", 0, 0, 0)); 3708 3709 // stop timing 3710 _statistics->transformTime->stop(); 3711 3712 if(numColsRational() > _beforeLiftCols || numRowsRational() > _beforeLiftRows) 3713 { 3714 SPX_MSG_INFO1(spxout, spxout << "Added " << numColsRational() - _beforeLiftCols << " columns and " 3715 << numRowsRational() - _beforeLiftRows << " rows to reduce large matrix coefficients\n."); 3716 } 3717 } 3718 3719 3720 3721 /// undoes lifting 3722 template <class R> 3723 void SoPlexBase<R>::_project(SolRational& sol) 3724 { 3725 // start timing 3726 _statistics->transformTime->start(); 3727 3728 // print LP if in debug mode 3729 SPX_DEBUG(_realLP->writeFileLPBase("beforeProject.lp", 0, 0, 0)); 3730 3731 assert(numColsRational() >= _beforeLiftCols); 3732 assert(numRowsRational() >= _beforeLiftRows); 3733 3734 // shrink rational LP to original size 3735 _rationalLP->removeColRange(_beforeLiftCols, numColsRational() - 1); 3736 _rationalLP->removeRowRange(_beforeLiftRows, numRowsRational() - 1); 3737 3738 // shrink real LP to original size 3739 _realLP->removeColRange(_beforeLiftCols, numColsReal() - 1); 3740 _realLP->removeRowRange(_beforeLiftRows, numRowsReal() - 1); 3741 3742 // adjust solution 3743 if(sol.isPrimalFeasible()) 3744 { 3745 sol._primal.reDim(_beforeLiftCols); 3746 sol._slacks.reDim(_beforeLiftRows); 3747 } 3748 3749 if(sol.hasPrimalRay()) 3750 { 3751 sol._primalRay.reDim(_beforeLiftCols); 3752 } 3753 3754 ///@todo if we know the mapping between original and lifting columns, we simply need to add the reduced cost of 3755 /// the lifting column to the reduced cost of the original column; this is not implemented now, because for 3756 /// optimal solutions the reduced costs of the lifting columns are zero 3757 const Rational maxValue = Rational(realParam(SoPlexBase<R>::LIFTMAXVAL)); 3758 3759 for(int i = _beforeLiftCols; i < numColsRational() && sol._isDualFeasible; i++) 3760 { 3761 if(spxAbs(Rational(maxValue * sol._redCost[i])) > _rationalOpttol) 3762 { 3763 SPX_MSG_INFO1(spxout, spxout << "Warning: lost dual solution during project phase.\n"); 3764 sol._isDualFeasible = false; 3765 } 3766 } 3767 3768 if(sol.isDualFeasible()) 3769 { 3770 sol._redCost.reDim(_beforeLiftCols); 3771 sol._dual.reDim(_beforeLiftRows); 3772 } 3773 3774 if(sol.hasDualFarkas()) 3775 { 3776 sol._dualFarkas.reDim(_beforeLiftRows); 3777 } 3778 3779 // adjust basis 3780 for(int i = _beforeLiftCols; i < numColsRational() && _hasBasis; i++) 3781 { 3782 if(_basisStatusCols[i] != SPxSolverBase<R>::BASIC) 3783 { 3784 SPX_MSG_INFO1(spxout, spxout << 3785 "Warning: lost basis during project phase because of nonbasic lifting column.\n"); 3786 _hasBasis = false; 3787 _rationalLUSolver.clear(); 3788 } 3789 } 3790 3791 for(int i = _beforeLiftRows; i < numRowsRational() && _hasBasis; i++) 3792 { 3793 if(_basisStatusRows[i] == SPxSolverBase<R>::BASIC) 3794 { 3795 SPX_MSG_INFO1(spxout, spxout << 3796 "Warning: lost basis during project phase because of basic lifting row.\n"); 3797 _hasBasis = false; 3798 _rationalLUSolver.clear(); 3799 } 3800 } 3801 3802 if(_hasBasis) 3803 { 3804 _basisStatusCols.reSize(_beforeLiftCols); 3805 _basisStatusRows.reSize(_beforeLiftRows); 3806 _rationalLUSolver.clear(); 3807 } 3808 3809 // print LP if in debug mode 3810 SPX_DEBUG(_realLP->writeFileLPBase("afterProject.lp", 0, 0, 0)); 3811 3812 // stop timing 3813 _statistics->transformTime->stop(); 3814 } 3815 3816 3817 3818 /// stores objective, bounds, and sides of real LP 3819 template <class R> 3820 void SoPlexBase<R>::_storeLPReal() 3821 { 3822 #ifndef SOPLEX_MANUAL_ALT 3823 3824 if(intParam(SoPlexBase<R>::SYNCMODE) == SYNCMODE_MANUAL) 3825 { 3826 _manualRealLP = *_realLP; 3827 return; 3828 } 3829 3830 #endif 3831 3832 _manualLower = _realLP->lower(); 3833 _manualUpper = _realLP->upper(); 3834 _manualLhs = _realLP->lhs(); 3835 _manualRhs = _realLP->rhs(); 3836 _manualObj.reDim(_realLP->nCols()); 3837 _realLP->getObj(_manualObj); 3838 } 3839 3840 3841 3842 /// restores objective, bounds, and sides of real LP 3843 template <class R> 3844 void SoPlexBase<R>::_restoreLPReal() 3845 { 3846 if(intParam(SoPlexBase<R>::SYNCMODE) == SYNCMODE_MANUAL) 3847 { 3848 #ifndef SOPLEX_MANUAL_ALT 3849 _solver.loadLP(_manualRealLP); 3850 #else 3851 _realLP->changeLower(_manualLower); 3852 _realLP->changeUpper(_manualUpper); 3853 _realLP->changeLhs(_manualLhs); 3854 _realLP->changeRhs(_manualRhs); 3855 _realLP->changeObj(_manualObj); 3856 #endif 3857 3858 if(_hasBasis) 3859 { 3860 // in manual sync mode, if the right-hand side of an equality constraint is not floating-point 3861 // representable, the user might have constructed the constraint in the real LP by rounding down the 3862 // left-hand side and rounding up the right-hand side; if the basis status is fixed, we need to adjust it 3863 for(int i = 0; i < _solver.nRows(); i++) 3864 { 3865 if(_basisStatusRows[i] == SPxSolverBase<R>::FIXED && _solver.lhs(i) != _solver.rhs(i)) 3866 { 3867 assert(_solver.rhs(i) == spxNextafter(_solver.lhs(i), R(infinity))); 3868 3869 if(_hasSolRational && _solRational.isDualFeasible() 3870 && ((intParam(SoPlexBase<R>::OBJSENSE) == SoPlexBase<R>::OBJSENSE_MAXIMIZE 3871 && _solRational._dual[i] > 0) 3872 || (intParam(SoPlexBase<R>::OBJSENSE) == SoPlexBase<R>::OBJSENSE_MINIMIZE 3873 && _solRational._dual[i] < 0))) 3874 { 3875 _basisStatusRows[i] = SPxSolverBase<R>::ON_UPPER; 3876 } 3877 else 3878 { 3879 _basisStatusRows[i] = SPxSolverBase<R>::ON_LOWER; 3880 } 3881 } 3882 } 3883 3884 _solver.setBasis(_basisStatusRows.get_const_ptr(), _basisStatusCols.get_const_ptr()); 3885 _hasBasis = (_solver.basis().status() > SPxBasisBase<R>::NO_PROBLEM); 3886 } 3887 } 3888 else 3889 { 3890 _realLP->changeLower(_manualLower); 3891 _realLP->changeUpper(_manualUpper); 3892 _realLP->changeLhs(_manualLhs); 3893 _realLP->changeRhs(_manualRhs); 3894 _realLP->changeObj(_manualObj); 3895 } 3896 } 3897 3898 3899 3900 /// introduces slack variables to transform inequality constraints into equations for both rational and real LP, 3901 /// which should be in sync 3902 template <class R> 3903 void SoPlexBase<R>::_transformEquality() 3904 { 3905 SPxOut::debug(this, "Transforming rows to equation form.\n"); 3906 3907 // start timing 3908 _statistics->transformTime->start(); 3909 3910 SPX_DEBUG(_realLP->writeFileLPBase("beforeTransEqu.lp", 0, 0, 0)); 3911 3912 // clear array of slack columns 3913 _slackCols.clear(); 3914 3915 // add artificial slack variables to convert inequality to equality constraints 3916 for(int i = 0; i < numRowsRational(); i++) 3917 { 3918 assert((lhsRational(i) == rhsRational(i)) == (_rowTypes[i] == RANGETYPE_FIXED)); 3919 3920 if(_rowTypes[i] != RANGETYPE_FIXED) 3921 { 3922 _slackCols.add(_rationalZero, -rhsRational(i), *_unitVectorRational(i), -lhsRational(i)); 3923 3924 if(_rationalLP->lhs(i) != 0) 3925 _rationalLP->changeLhs(i, _rationalZero); 3926 3927 if(_rationalLP->rhs(i) != 0) 3928 _rationalLP->changeRhs(i, _rationalZero); 3929 3930 assert(_rationalLP->lhs(i) == 0); 3931 assert(_rationalLP->rhs(i) == 0); 3932 _realLP->changeRange(i, R(0.0), R(0.0)); 3933 _colTypes.append(_switchRangeType(_rowTypes[i])); 3934 _rowTypes[i] = RANGETYPE_FIXED; 3935 } 3936 } 3937 3938 _rationalLP->addCols(_slackCols); 3939 _realLP->addCols(_slackCols); 3940 3941 // adjust basis 3942 if(_hasBasis) 3943 { 3944 for(int i = 0; i < _slackCols.num(); i++) 3945 { 3946 int row = _slackCols.colVector(i).index(0); 3947 3948 assert(row >= 0); 3949 assert(row < numRowsRational()); 3950 3951 switch(_basisStatusRows[row]) 3952 { 3953 case SPxSolverBase<R>::ON_LOWER: 3954 _basisStatusCols.append(SPxSolverBase<R>::ON_UPPER); 3955 break; 3956 3957 case SPxSolverBase<R>::ON_UPPER: 3958 _basisStatusCols.append(SPxSolverBase<R>::ON_LOWER); 3959 break; 3960 3961 case SPxSolverBase<R>::BASIC: 3962 case SPxSolverBase<R>::FIXED: 3963 default: 3964 _basisStatusCols.append(_basisStatusRows[row]); 3965 break; 3966 } 3967 3968 _basisStatusRows[row] = SPxSolverBase<R>::FIXED; 3969 } 3970 3971 _rationalLUSolver.clear(); 3972 } 3973 3974 SPX_DEBUG(_realLP->writeFileLPBase("afterTransEqu.lp", 0, 0, 0)); 3975 3976 // stop timing 3977 _statistics->transformTime->stop(); 3978 3979 if(_slackCols.num() > 0) 3980 { 3981 SPX_MSG_INFO1(spxout, spxout << "Added " << _slackCols.num() << 3982 " slack columns to transform rows to equality form.\n"); 3983 } 3984 } 3985 3986 3987 3988 /// restores original problem 3989 template <class R> 3990 void SoPlexBase<R>::_untransformEquality(SolRational& sol) 3991 { 3992 // start timing 3993 _statistics->transformTime->start(); 3994 3995 // print LP if in debug mode 3996 SPX_DEBUG(_realLP->writeFileLPBase("beforeUntransEqu.lp", 0, 0, 0)); 3997 3998 int numCols = numColsRational(); 3999 int numOrigCols = numColsRational() - _slackCols.num(); 4000 4001 // adjust solution 4002 if(sol.isPrimalFeasible()) 4003 { 4004 for(int i = 0; i < _slackCols.num(); i++) 4005 { 4006 int col = numOrigCols + i; 4007 int row = _slackCols.colVector(i).index(0); 4008 4009 assert(row >= 0); 4010 assert(row < numRowsRational()); 4011 4012 sol._slacks[row] -= sol._primal[col]; 4013 } 4014 4015 sol._primal.reDim(numOrigCols); 4016 } 4017 4018 if(sol.hasPrimalRay()) 4019 { 4020 sol._primalRay.reDim(numOrigCols); 4021 } 4022 4023 // adjust basis 4024 if(_hasBasis) 4025 { 4026 for(int i = 0; i < _slackCols.num(); i++) 4027 { 4028 int col = numOrigCols + i; 4029 int row = _slackCols.colVector(i).index(0); 4030 4031 assert(row >= 0); 4032 assert(row < numRowsRational()); 4033 assert(_basisStatusRows[row] != SPxSolverBase<R>::UNDEFINED); 4034 assert(_basisStatusRows[row] != SPxSolverBase<R>::ZERO || lhsRational(row) == 0); 4035 assert(_basisStatusRows[row] != SPxSolverBase<R>::ZERO || rhsRational(row) == 0); 4036 assert(_basisStatusRows[row] != SPxSolverBase<R>::BASIC 4037 || _basisStatusCols[col] != SPxSolverBase<R>::BASIC); 4038 4039 SPxOut::debug(this, 4040 "slack column {} for row {}: col status={}, row status={}, redcost={}, dual={}\n", 4041 col, row, _basisStatusCols[col], _basisStatusRows[row], 4042 sol._redCost[col].str(), sol._dual[row].str()); 4043 4044 if(_basisStatusRows[row] != SPxSolverBase<R>::BASIC) 4045 { 4046 switch(_basisStatusCols[col]) 4047 { 4048 case SPxSolverBase<R>::ON_LOWER: 4049 _basisStatusRows[row] = SPxSolverBase<R>::ON_UPPER; 4050 break; 4051 4052 case SPxSolverBase<R>::ON_UPPER: 4053 _basisStatusRows[row] = SPxSolverBase<R>::ON_LOWER; 4054 break; 4055 4056 case SPxSolverBase<R>::BASIC: 4057 case SPxSolverBase<R>::FIXED: 4058 default: 4059 _basisStatusRows[row] = _basisStatusCols[col]; 4060 break; 4061 } 4062 } 4063 } 4064 4065 _basisStatusCols.reSize(numOrigCols); 4066 4067 if(_slackCols.num() > 0) 4068 _rationalLUSolver.clear(); 4069 } 4070 4071 // not earlier because of debug message 4072 if(sol.isDualFeasible()) 4073 { 4074 sol._redCost.reDim(numOrigCols); 4075 } 4076 4077 // restore sides and remove slack columns 4078 for(int i = 0; i < _slackCols.num(); i++) 4079 { 4080 int col = numOrigCols + i; 4081 int row = _slackCols.colVector(i).index(0); 4082 4083 if(upperRational(col) != 0) 4084 _rationalLP->changeLhs(row, -upperRational(col)); 4085 4086 if(lowerRational(col) != 0) 4087 _rationalLP->changeRhs(row, -lowerRational(col)); 4088 4089 assert(_rationalLP->lhs(row) == -upperRational(col)); 4090 assert(_rationalLP->rhs(row) == -lowerRational(col)); 4091 _rowTypes[row] = _switchRangeType(_colTypes[col]); 4092 } 4093 4094 _rationalLP->removeColRange(numOrigCols, numCols - 1); 4095 _realLP->removeColRange(numOrigCols, numCols - 1); 4096 _colTypes.reSize(numOrigCols); 4097 4098 // objective, bounds, and sides of real LP are restored only after _solveRational() 4099 4100 // print LP if in debug mode 4101 SPX_DEBUG(_realLP->writeFileLPBase("afterUntransEqu.lp", 0, 0, 0)); 4102 4103 // stop timing 4104 _statistics->transformTime->stop(); 4105 } 4106 4107 4108 4109 /// transforms LP to unboundedness problem by moving the objective function to the constraints, changing right-hand 4110 /// side and bounds to zero, and adding an auxiliary variable for the decrease in the objective function 4111 template <class R> 4112 void SoPlexBase<R>::_transformUnbounded() 4113 { 4114 SPX_MSG_INFO1(spxout, spxout << "Setting up LP to compute primal unbounded ray.\n"); 4115 4116 // start timing 4117 _statistics->transformTime->start(); 4118 4119 // print LP if in debug mode 4120 SPX_DEBUG(_realLP->writeFileLPBase("beforeTransUnbounded.lp", 0, 0, 0)); 4121 4122 // store bounds 4123 _unboundedLower.reDim(numColsRational()); 4124 _unboundedUpper.reDim(numColsRational()); 4125 4126 for(int c = numColsRational() - 1; c >= 0; c--) 4127 { 4128 if(_lowerFinite(_colTypes[c])) 4129 _unboundedLower[c] = lowerRational(c); 4130 4131 if(_upperFinite(_colTypes[c])) 4132 _unboundedUpper[c] = upperRational(c); 4133 } 4134 4135 // store sides 4136 _unboundedLhs.reDim(numRowsRational()); 4137 _unboundedRhs.reDim(numRowsRational()); 4138 4139 for(int r = numRowsRational() - 1; r >= 0; r--) 4140 { 4141 if(_lowerFinite(_rowTypes[r])) 4142 _unboundedLhs[r] = lhsRational(r); 4143 4144 if(_upperFinite(_rowTypes[r])) 4145 _unboundedRhs[r] = rhsRational(r); 4146 } 4147 4148 // make right-hand side zero 4149 for(int r = numRowsRational() - 1; r >= 0; r--) 4150 { 4151 assert((lhsRational(r) > _rationalNegInfty) == _lowerFinite(_rowTypes[r])); 4152 4153 if(_lowerFinite(_rowTypes[r])) 4154 { 4155 _rationalLP->changeLhs(r, Rational(0)); 4156 _realLP->changeLhs(r, R(0.0)); 4157 } 4158 else if(_realLP->lhs(r) > -realParam(SoPlexBase<R>::INFTY)) 4159 _realLP->changeLhs(r, -realParam(SoPlexBase<R>::INFTY)); 4160 4161 assert((rhsRational(r) < _rationalPosInfty) == _upperFinite(_rowTypes[r])); 4162 4163 if(_upperFinite(_rowTypes[r])) 4164 { 4165 _rationalLP->changeRhs(r, Rational(0)); 4166 _realLP->changeRhs(r, R(0.0)); 4167 } 4168 else if(_realLP->rhs(r) < realParam(SoPlexBase<R>::INFTY)) 4169 _realLP->changeRhs(r, realParam(SoPlexBase<R>::INFTY)); 4170 } 4171 4172 // transform objective function to constraint and add auxiliary variable 4173 int numOrigCols = numColsRational(); 4174 DSVectorRational obj(numOrigCols + 1); 4175 ///@todo implement this without copying the objective function 4176 obj = _rationalLP->maxObj(); 4177 obj.add(numOrigCols, -1); 4178 _rationalLP->addRow(LPRowRational(0, obj, 0)); 4179 _realLP->addRow(LPRowBase<R>(0, DSVectorBase<R>(obj), 0)); 4180 _rowTypes.append(RANGETYPE_FIXED); 4181 4182 assert(numColsRational() == numOrigCols + 1); 4183 4184 // set objective coefficient and bounds for auxiliary variable 4185 _rationalLP->changeMaxObj(numOrigCols, Rational(1)); 4186 _realLP->changeMaxObj(numOrigCols, R(1.0)); 4187 4188 _rationalLP->changeBounds(numOrigCols, _rationalNegInfty, 1); 4189 _realLP->changeBounds(numOrigCols, -realParam(SoPlexBase<R>::INFTY), 1.0); 4190 _colTypes.append(RANGETYPE_UPPER); 4191 4192 // set objective coefficients to zero and adjust bounds for problem variables 4193 for(int c = numColsRational() - 2; c >= 0; c--) 4194 { 4195 _rationalLP->changeObj(c, Rational(0)); 4196 _realLP->changeObj(c, R(0.0)); 4197 4198 assert((lowerRational(c) > _rationalNegInfty) == _lowerFinite(_colTypes[c])); 4199 4200 if(_lowerFinite(_colTypes[c])) 4201 { 4202 _rationalLP->changeLower(c, Rational(0)); 4203 _realLP->changeLower(c, R(0.0)); 4204 } 4205 else if(_realLP->lower(c) > -realParam(SoPlexBase<R>::INFTY)) 4206 _realLP->changeLower(c, -realParam(SoPlexBase<R>::INFTY)); 4207 4208 assert((upperRational(c) < _rationalPosInfty) == _upperFinite(_colTypes[c])); 4209 4210 if(_upperFinite(_colTypes[c])) 4211 { 4212 _rationalLP->changeUpper(c, Rational(0)); 4213 _realLP->changeUpper(c, R(0.0)); 4214 } 4215 else if(_realLP->upper(c) < realParam(SoPlexBase<R>::INFTY)) 4216 _realLP->changeUpper(c, realParam(SoPlexBase<R>::INFTY)); 4217 } 4218 4219 // adjust basis 4220 if(_hasBasis) 4221 { 4222 _basisStatusCols.append(SPxSolverBase<R>::ON_UPPER); 4223 _basisStatusRows.append(SPxSolverBase<R>::BASIC); 4224 _rationalLUSolver.clear(); 4225 } 4226 4227 // print LP if in debug mode 4228 SPX_DEBUG(_realLP->writeFileLPBase("afterTransUnbounded.lp", 0, 0, 0)); 4229 4230 // stop timing 4231 _statistics->transformTime->stop(); 4232 } 4233 4234 4235 4236 /// undoes transformation to unboundedness problem 4237 template <class R> 4238 void SoPlexBase<R>::_untransformUnbounded(SolRational& sol, bool unbounded) 4239 { 4240 // start timing 4241 _statistics->transformTime->start(); 4242 4243 // print LP if in debug mode 4244 SPX_DEBUG(_realLP->writeFileLPBase("beforeUntransUnbounded.lp", 0, 0, 0)); 4245 4246 int numOrigCols = numColsRational() - 1; 4247 int numOrigRows = numRowsRational() - 1; 4248 const Rational& tau = sol._primal[numOrigCols]; 4249 4250 // adjust solution and basis 4251 if(unbounded) 4252 { 4253 assert(tau >= _rationalPosone); 4254 4255 sol._isPrimalFeasible = false; 4256 sol._hasPrimalRay = true; 4257 sol._isDualFeasible = false; 4258 sol._hasDualFarkas = false; 4259 4260 if(tau != 1) 4261 sol._primal /= tau; 4262 4263 sol._primalRay = sol._primal; 4264 sol._primalRay.reDim(numOrigCols); 4265 4266 _hasBasis = (_basisStatusCols[numOrigCols] != SPxSolverBase<R>::BASIC 4267 && _basisStatusRows[numOrigRows] == SPxSolverBase<R>::BASIC); 4268 _basisStatusCols.reSize(numOrigCols); 4269 _basisStatusRows.reSize(numOrigRows); 4270 } 4271 else if(boolParam(SoPlexBase<R>::TESTDUALINF) && tau < _rationalFeastol) 4272 { 4273 const Rational& alpha = sol._dual[numOrigRows]; 4274 4275 assert(sol._isDualFeasible); 4276 assert(alpha <= _rationalFeastol - _rationalPosone); 4277 4278 sol._isPrimalFeasible = false; 4279 sol._hasPrimalRay = false; 4280 sol._hasDualFarkas = false; 4281 4282 if(alpha != -1) 4283 { 4284 sol._dual /= -alpha; 4285 sol._redCost /= -alpha; 4286 } 4287 4288 sol._dual.reDim(numOrigRows); 4289 sol._redCost.reDim(numOrigCols); 4290 } 4291 else 4292 { 4293 sol.invalidate(); 4294 _hasBasis = false; 4295 _basisStatusCols.reSize(numOrigCols); 4296 _basisStatusCols.reSize(numOrigRows); 4297 } 4298 4299 // recover objective function 4300 const SVectorRational& objRowVector = _rationalLP->rowVector(numOrigRows); 4301 4302 for(int i = objRowVector.size() - 1; i >= 0; i--) 4303 { 4304 _rationalLP->changeMaxObj(objRowVector.index(i), objRowVector.value(i)); 4305 _realLP->changeMaxObj(objRowVector.index(i), R(objRowVector.value(i))); 4306 } 4307 4308 // remove objective function constraint and auxiliary variable 4309 _rationalLP->removeRow(numOrigRows); 4310 _realLP->removeRow(numOrigRows); 4311 _rowTypes.reSize(numOrigRows); 4312 4313 _rationalLP->removeCol(numOrigCols); 4314 _realLP->removeCol(numOrigCols); 4315 _colTypes.reSize(numOrigCols); 4316 4317 // restore objective, sides and bounds 4318 for(int r = numRowsRational() - 1; r >= 0; r--) 4319 { 4320 if(_lowerFinite(_rowTypes[r])) 4321 { 4322 _rationalLP->changeLhs(r, _unboundedLhs[r]); 4323 _realLP->changeLhs(r, R(_unboundedLhs[r])); 4324 } 4325 4326 if(_upperFinite(_rowTypes[r])) 4327 { 4328 _rationalLP->changeRhs(r, _unboundedRhs[r]); 4329 _realLP->changeRhs(r, R(_unboundedRhs[r])); 4330 } 4331 4332 assert((lhsRational(r) > _rationalNegInfty) == _lowerFinite(_rowTypes[r])); 4333 assert((rhsRational(r) < _rationalPosInfty) == _upperFinite(_rowTypes[r])); 4334 assert((lhsReal(r) > -realParam(SoPlexBase<R>::INFTY)) == _lowerFinite(_rowTypes[r])); 4335 assert((rhsReal(r) < realParam(SoPlexBase<R>::INFTY)) == _upperFinite(_rowTypes[r])); 4336 } 4337 4338 for(int c = numColsRational() - 1; c >= 0; c--) 4339 { 4340 if(_lowerFinite(_colTypes[c])) 4341 { 4342 _rationalLP->changeLower(c, _unboundedLower[c]); 4343 _realLP->changeLower(c, R(_unboundedLower[c])); 4344 } 4345 4346 if(_upperFinite(_colTypes[c])) 4347 { 4348 _rationalLP->changeUpper(c, _unboundedUpper[c]); 4349 _realLP->changeUpper(c, R(_unboundedUpper[c])); 4350 } 4351 4352 assert((lowerRational(c) > _rationalNegInfty) == _lowerFinite(_colTypes[c])); 4353 assert((upperRational(c) < _rationalPosInfty) == _upperFinite(_colTypes[c])); 4354 assert((lowerReal(c) > -realParam(SoPlexBase<R>::INFTY)) == _lowerFinite(_colTypes[c])); 4355 assert((upperReal(c) < realParam(SoPlexBase<R>::INFTY)) == _upperFinite(_colTypes[c])); 4356 } 4357 4358 // invalidate rational basis factorization 4359 _rationalLUSolver.clear(); 4360 4361 // print LP if in debug mode 4362 SPX_DEBUG(_realLP->writeFileLPBase("afterUntransUnbounded.lp", 0, 0, 0)); 4363 4364 // stop timing 4365 _statistics->transformTime->stop(); 4366 } 4367 4368 4369 4370 /// store basis 4371 template <class R> 4372 void SoPlexBase<R>::_storeBasis() 4373 { 4374 assert(!_storedBasis); 4375 4376 if(_hasBasis) 4377 { 4378 _storedBasis = true; 4379 _storedBasisStatusCols = _basisStatusCols; 4380 _storedBasisStatusRows = _basisStatusRows; 4381 } 4382 else 4383 _storedBasis = false; 4384 } 4385 4386 4387 4388 /// restore basis 4389 template <class R> 4390 void SoPlexBase<R>::_restoreBasis() 4391 { 4392 if(_storedBasis) 4393 { 4394 _hasBasis = true; 4395 _basisStatusCols = _storedBasisStatusCols; 4396 _basisStatusRows = _storedBasisStatusRows; 4397 _storedBasis = false; 4398 } 4399 } 4400 4401 4402 4403 /// transforms LP to feasibility problem by removing the objective function, shifting variables, and homogenizing the 4404 /// right-hand side 4405 template <class R> 4406 void SoPlexBase<R>::_transformFeasibility() 4407 { 4408 SPX_MSG_INFO1(spxout, spxout << "Setting up LP to test for feasibility.\n"); 4409 4410 // start timing 4411 _statistics->transformTime->start(); 4412 4413 // print LP if in debug mode 4414 SPX_DEBUG(_realLP->writeFileLPBase("beforeTransFeas.lp", 0, 0, 0)); 4415 4416 // store objective function 4417 _feasObj.reDim(numColsRational()); 4418 4419 for(int c = numColsRational() - 1; c >= 0; c--) 4420 _feasObj[c] = _rationalLP->maxObj(c); 4421 4422 // store bounds 4423 _feasLower.reDim(numColsRational()); 4424 _feasUpper.reDim(numColsRational()); 4425 4426 for(int c = numColsRational() - 1; c >= 0; c--) 4427 { 4428 if(_lowerFinite(_colTypes[c])) 4429 _feasLower[c] = lowerRational(c); 4430 4431 if(_upperFinite(_colTypes[c])) 4432 _feasUpper[c] = upperRational(c); 4433 } 4434 4435 // store sides 4436 _feasLhs.reDim(numRowsRational()); 4437 _feasRhs.reDim(numRowsRational()); 4438 4439 for(int r = numRowsRational() - 1; r >= 0; r--) 4440 { 4441 if(_lowerFinite(_rowTypes[r])) 4442 _feasLhs[r] = lhsRational(r); 4443 4444 if(_upperFinite(_rowTypes[r])) 4445 _feasRhs[r] = rhsRational(r); 4446 } 4447 4448 // set objective coefficients to zero; shift primal space such as to guarantee that the zero solution is within 4449 // the bounds 4450 Rational shiftValue; 4451 Rational shiftValue2; 4452 4453 for(int c = numColsRational() - 1; c >= 0; c--) 4454 { 4455 _rationalLP->changeMaxObj(c, Rational(0)); 4456 _realLP->changeMaxObj(c, R(0.0)); 4457 4458 if(lowerRational(c) > 0) 4459 { 4460 const SVectorRational& colVector = colVectorRational(c); 4461 4462 for(int i = 0; i < colVector.size(); i++) 4463 { 4464 shiftValue = colVector.value(i); 4465 shiftValue *= lowerRational(c); 4466 int r = colVector.index(i); 4467 4468 assert((lhsRational(r) > _rationalNegInfty) == _lowerFinite(_rowTypes[r])); 4469 assert((rhsRational(r) < _rationalPosInfty) == _upperFinite(_rowTypes[r])); 4470 4471 if(_lowerFinite(_rowTypes[r]) && _upperFinite(_rowTypes[r])) 4472 { 4473 shiftValue2 = lhsRational(r); 4474 shiftValue2 -= shiftValue; 4475 _rationalLP->changeLhs(r, shiftValue2); 4476 _realLP->changeLhs(r, R(shiftValue2)); 4477 4478 shiftValue -= rhsRational(r); 4479 shiftValue *= -1; 4480 _rationalLP->changeRhs(r, shiftValue); 4481 _realLP->changeRhs(r, R(shiftValue)); 4482 } 4483 else if(_lowerFinite(_rowTypes[r])) 4484 { 4485 shiftValue -= lhsRational(r); 4486 shiftValue *= -1; 4487 _rationalLP->changeLhs(r, shiftValue); 4488 _realLP->changeLhs(r, R(shiftValue)); 4489 } 4490 else if(_upperFinite(_rowTypes[r])) 4491 { 4492 shiftValue -= rhsRational(r); 4493 shiftValue *= -1; 4494 _rationalLP->changeRhs(r, shiftValue); 4495 _realLP->changeRhs(r, R(shiftValue)); 4496 } 4497 } 4498 4499 assert((upperRational(c) < _rationalPosInfty) == _upperFinite(_colTypes[c])); 4500 4501 if(_upperFinite(_colTypes[c])) 4502 { 4503 _rationalLP->changeBounds(c, 0, upperRational(c) - lowerRational(c)); 4504 _realLP->changeBounds(c, 0.0, R(upperRational(c))); 4505 } 4506 else if(_realLP->upper(c) < realParam(SoPlexBase<R>::INFTY)) 4507 { 4508 _rationalLP->changeLower(c, Rational(0)); 4509 _realLP->changeBounds(c, 0.0, realParam(SoPlexBase<R>::INFTY)); 4510 } 4511 else 4512 { 4513 _rationalLP->changeLower(c, Rational(0)); 4514 _realLP->changeLower(c, R(0.0)); 4515 } 4516 } 4517 else if(upperRational(c) < 0) 4518 { 4519 const SVectorRational& colVector = colVectorRational(c); 4520 4521 for(int i = 0; i < colVector.size(); i++) 4522 { 4523 shiftValue = colVector.value(i); 4524 shiftValue *= upperRational(c); 4525 int r = colVector.index(i); 4526 4527 assert((lhsRational(r) > _rationalNegInfty) == _lowerFinite(_rowTypes[r])); 4528 assert((rhsRational(r) < _rationalPosInfty) == _upperFinite(_rowTypes[r])); 4529 4530 if(_lowerFinite(_rowTypes[r]) && _upperFinite(_rowTypes[r])) 4531 { 4532 shiftValue2 = lhsRational(r); 4533 shiftValue2 -= shiftValue; 4534 _rationalLP->changeLhs(r, shiftValue2); 4535 _realLP->changeLhs(r, R(shiftValue2)); 4536 4537 shiftValue -= rhsRational(r); 4538 shiftValue *= -1; 4539 _rationalLP->changeRhs(r, shiftValue); 4540 _realLP->changeRhs(r, R(shiftValue)); 4541 } 4542 else if(_lowerFinite(_rowTypes[r])) 4543 { 4544 shiftValue -= lhsRational(r); 4545 shiftValue *= -1; 4546 _rationalLP->changeLhs(r, shiftValue); 4547 _realLP->changeLhs(r, R(shiftValue)); 4548 } 4549 else if(_upperFinite(_rowTypes[r])) 4550 { 4551 shiftValue -= rhsRational(r); 4552 shiftValue *= -1; 4553 _rationalLP->changeRhs(r, shiftValue); 4554 _realLP->changeRhs(r, R(shiftValue)); 4555 } 4556 } 4557 4558 assert((lowerRational(c) > _rationalNegInfty) == _lowerFinite(_colTypes[c])); 4559 4560 if(_lowerFinite(_colTypes[c])) 4561 { 4562 _rationalLP->changeBounds(c, lowerRational(c) - upperRational(c), 0); 4563 _realLP->changeBounds(c, R(lowerRational(c)), 0.0); 4564 } 4565 else if(_realLP->lower(c) > -realParam(SoPlexBase<R>::INFTY)) 4566 { 4567 _rationalLP->changeUpper(c, Rational(0)); 4568 _realLP->changeBounds(c, -realParam(SoPlexBase<R>::INFTY), 0.0); 4569 } 4570 else 4571 { 4572 _rationalLP->changeUpper(c, Rational(0)); 4573 _realLP->changeUpper(c, R(0.0)); 4574 } 4575 } 4576 else 4577 { 4578 if(_lowerFinite(_colTypes[c])) 4579 _realLP->changeLower(c, R(lowerRational(c))); 4580 else if(_realLP->lower(c) > -realParam(SoPlexBase<R>::INFTY)) 4581 _realLP->changeLower(c, -realParam(SoPlexBase<R>::INFTY)); 4582 4583 if(_upperFinite(_colTypes[c])) 4584 _realLP->changeUpper(c, R(upperRational(c))); 4585 else if(_realLP->upper(c) < realParam(SoPlexBase<R>::INFTY)) 4586 _realLP->changeUpper(c, realParam(SoPlexBase<R>::INFTY)); 4587 } 4588 4589 assert(lowerReal(c) <= upperReal(c)); 4590 } 4591 4592 // homogenize sides 4593 _tauColVector.clear(); 4594 4595 for(int r = numRowsRational() - 1; r >= 0; r--) 4596 { 4597 if(lhsRational(r) > 0) 4598 { 4599 _tauColVector.add(r, lhsRational(r)); 4600 assert((rhsRational(r) < _rationalPosInfty) == _upperFinite(_rowTypes[r])); 4601 4602 if(_upperFinite(_rowTypes[r])) 4603 { 4604 _rationalLP->changeRange(r, 0, rhsRational(r) - lhsRational(r)); 4605 _realLP->changeRange(r, 0.0, R(rhsRational(r))); 4606 } 4607 else 4608 { 4609 _rationalLP->changeLhs(r, Rational(0)); 4610 _realLP->changeLhs(r, R(0.0)); 4611 4612 if(_realLP->rhs(r) < realParam(SoPlexBase<R>::INFTY)) 4613 _realLP->changeRhs(r, realParam(SoPlexBase<R>::INFTY)); 4614 } 4615 } 4616 else if(rhsRational(r) < 0) 4617 { 4618 _tauColVector.add(r, rhsRational(r)); 4619 assert((lhsRational(r) > _rationalNegInfty) == _lowerFinite(_rowTypes[r])); 4620 4621 if(_lowerFinite(_rowTypes[r])) 4622 { 4623 _rationalLP->changeRange(r, lhsRational(r) - rhsRational(r), 0); 4624 _realLP->changeRange(r, R(lhsRational(r)), 0.0); 4625 } 4626 else 4627 { 4628 _rationalLP->changeRhs(r, Rational(0)); 4629 _realLP->changeRhs(r, R(0.0)); 4630 4631 if(_realLP->lhs(r) > -realParam(SoPlexBase<R>::INFTY)) 4632 _realLP->changeLhs(r, -realParam(SoPlexBase<R>::INFTY)); 4633 } 4634 } 4635 else 4636 { 4637 if(_lowerFinite(_rowTypes[r])) 4638 _realLP->changeLhs(r, R(lhsRational(r))); 4639 else if(_realLP->lhs(r) > -realParam(SoPlexBase<R>::INFTY)) 4640 _realLP->changeLhs(r, -realParam(SoPlexBase<R>::INFTY)); 4641 4642 if(_upperFinite(_rowTypes[r])) 4643 _realLP->changeRhs(r, R(rhsRational(r))); 4644 else if(_realLP->rhs(r) < realParam(SoPlexBase<R>::INFTY)) 4645 _realLP->changeRhs(r, realParam(SoPlexBase<R>::INFTY)); 4646 } 4647 4648 assert(rhsReal(r) <= rhsReal(r)); 4649 } 4650 4651 ///@todo exploit this case by returning without LP solving 4652 if(_tauColVector.size() == 0) 4653 { 4654 SPX_MSG_INFO3(spxout, spxout << "LP is trivially feasible.\n"); 4655 } 4656 4657 // add artificial column 4658 SPxColId id; 4659 _tauColVector *= -1; 4660 _rationalLP->addCol(id, 4661 LPColRational((intParam(SoPlexBase<R>::OBJSENSE) == SoPlexBase<R>::OBJSENSE_MAXIMIZE ? 4662 _rationalPosone : _rationalNegone), 4663 _tauColVector, 1, 0)); 4664 _realLP->addCol(id, 4665 LPColBase<R>((intParam(SoPlexBase<R>::OBJSENSE) == SoPlexBase<R>::OBJSENSE_MAXIMIZE ? 1.0 : -1.0), 4666 DSVectorBase<R>(_tauColVector), 1.0, 0.0)); 4667 _colTypes.append(RANGETYPE_BOXED); 4668 4669 // adjust basis 4670 if(_hasBasis) 4671 { 4672 _basisStatusCols.append(SPxSolverBase<R>::ON_UPPER); 4673 } 4674 4675 // invalidate rational basis factorization 4676 _rationalLUSolver.clear(); 4677 4678 // print LP if in debug mode 4679 SPX_DEBUG(_realLP->writeFileLPBase("afterTransFeas.lp", 0, 0, 0)); 4680 4681 // stop timing 4682 _statistics->transformTime->stop(); 4683 } 4684 4685 4686 4687 /// undoes transformation to feasibility problem 4688 template <class R> 4689 void SoPlexBase<R>::_untransformFeasibility(SolRational& sol, bool infeasible) 4690 { 4691 // start timing 4692 _statistics->transformTime->start(); 4693 4694 // print LP if in debug mode 4695 SPX_DEBUG(_realLP->writeFileLPBase("beforeUntransFeas.lp", 0, 0, 0)); 4696 4697 int numOrigCols = numColsRational() - 1; 4698 4699 // adjust solution and basis 4700 if(infeasible) 4701 { 4702 assert(sol._isDualFeasible); 4703 assert(sol._primal[numOrigCols] < 1); 4704 4705 sol._isPrimalFeasible = false; 4706 sol._hasPrimalRay = false; 4707 sol._isDualFeasible = false; 4708 sol._hasDualFarkas = true; 4709 4710 sol._dualFarkas = sol._dual; 4711 4712 _hasBasis = false; 4713 _basisStatusCols.reSize(numOrigCols); 4714 } 4715 else if(sol._isPrimalFeasible) 4716 { 4717 assert(sol._primal[numOrigCols] >= 1); 4718 4719 sol._hasPrimalRay = false; 4720 sol._isDualFeasible = false; 4721 sol._hasDualFarkas = false; 4722 4723 if(sol._primal[numOrigCols] != 1) 4724 { 4725 sol._slacks /= sol._primal[numOrigCols]; 4726 4727 for(int i = 0; i < numOrigCols; i++) 4728 sol._primal[i] /= sol._primal[numOrigCols]; 4729 4730 sol._primal[numOrigCols] = 1; 4731 } 4732 4733 sol._primal.reDim(numOrigCols); 4734 sol._slacks -= _rationalLP->colVector(numOrigCols); 4735 4736 _hasBasis = (_basisStatusCols[numOrigCols] != SPxSolverBase<R>::BASIC); 4737 _basisStatusCols.reSize(numOrigCols); 4738 } 4739 else 4740 { 4741 _hasBasis = false; 4742 _basisStatusCols.reSize(numOrigCols); 4743 } 4744 4745 // restore right-hand side 4746 for(int r = numRowsRational() - 1; r >= 0; r--) 4747 { 4748 assert(rhsRational(r) >= _rationalPosInfty || lhsRational(r) <= _rationalNegInfty 4749 || _feasLhs[r] - lhsRational(r) == _feasRhs[r] - rhsRational(r)); 4750 4751 if(_lowerFinite(_rowTypes[r])) 4752 { 4753 _rationalLP->changeLhs(r, _feasLhs[r]); 4754 _realLP->changeLhs(r, R(_feasLhs[r])); 4755 } 4756 else if(_realLP->lhs(r) > -realParam(SoPlexBase<R>::INFTY)) 4757 _realLP->changeLhs(r, -realParam(SoPlexBase<R>::INFTY)); 4758 4759 assert(_lowerFinite(_rowTypes[r]) == (lhsRational(r) > _rationalNegInfty)); 4760 assert(_lowerFinite(_rowTypes[r]) == (lhsReal(r) > -realParam(SoPlexBase<R>::INFTY))); 4761 4762 if(_upperFinite(_rowTypes[r])) 4763 { 4764 _rationalLP->changeRhs(r, _feasRhs[r]); 4765 _realLP->changeRhs(r, R(_feasRhs[r])); 4766 } 4767 else if(_realLP->rhs(r) < realParam(SoPlexBase<R>::INFTY)) 4768 _realLP->changeRhs(r, realParam(SoPlexBase<R>::INFTY)); 4769 4770 assert(_upperFinite(_rowTypes[r]) == (rhsRational(r) < _rationalPosInfty)); 4771 assert(_upperFinite(_rowTypes[r]) == (rhsReal(r) < realParam(SoPlexBase<R>::INFTY))); 4772 4773 assert(lhsReal(r) <= rhsReal(r)); 4774 } 4775 4776 // unshift primal space and restore objective coefficients 4777 Rational shiftValue; 4778 4779 for(int c = numOrigCols - 1; c >= 0; c--) 4780 { 4781 bool shifted = (_lowerFinite(_colTypes[c]) && _feasLower[c] > 0) || (_upperFinite(_colTypes[c]) 4782 && _feasUpper[c] < 0); 4783 assert(shifted || !_lowerFinite(_colTypes[c]) || _feasLower[c] == lowerRational(c)); 4784 assert(shifted || !_upperFinite(_colTypes[c]) || _feasUpper[c] == upperRational(c)); 4785 assert(upperRational(c) >= _rationalPosInfty || lowerRational(c) <= _rationalNegInfty 4786 || _feasLower[c] - lowerRational(c) == _feasUpper[c] - upperRational(c)); 4787 4788 if(shifted) 4789 { 4790 if(_lowerFinite(_colTypes[c])) 4791 { 4792 shiftValue = _feasLower[c]; 4793 shiftValue -= lowerRational(c); 4794 } 4795 else if(_upperFinite(_colTypes[c])) 4796 { 4797 shiftValue = _feasUpper[c]; 4798 shiftValue -= upperRational(c); 4799 } 4800 4801 if(sol._isPrimalFeasible) 4802 { 4803 sol._primal[c] += shiftValue; 4804 sol._slacks.multAdd(shiftValue, _rationalLP->colVector(c)); 4805 } 4806 } 4807 4808 if(_lowerFinite(_colTypes[c])) 4809 { 4810 if(shifted) 4811 _rationalLP->changeLower(c, _feasLower[c]); 4812 4813 _realLP->changeLower(c, R(_feasLower[c])); 4814 } 4815 else if(_realLP->lower(c) > -realParam(SoPlexBase<R>::INFTY)) 4816 _realLP->changeLower(c, -realParam(SoPlexBase<R>::INFTY)); 4817 4818 assert(_lowerFinite(_colTypes[c]) == (lowerRational(c) > -_rationalPosInfty)); 4819 assert(_lowerFinite(_colTypes[c]) == (lowerReal(c) > -realParam(SoPlexBase<R>::INFTY))); 4820 4821 if(_upperFinite(_colTypes[c])) 4822 { 4823 if(shifted) 4824 _rationalLP->changeUpper(c, _feasUpper[c]); 4825 4826 _realLP->changeUpper(c, R(upperRational(c))); 4827 } 4828 else if(_realLP->upper(c) < realParam(SoPlexBase<R>::INFTY)) 4829 _realLP->changeUpper(c, realParam(SoPlexBase<R>::INFTY)); 4830 4831 assert(_upperFinite(_colTypes[c]) == (upperRational(c) < _rationalPosInfty)); 4832 assert(_upperFinite(_colTypes[c]) == (upperReal(c) < realParam(SoPlexBase<R>::INFTY))); 4833 4834 _rationalLP->changeMaxObj(c, _feasObj[c]); 4835 _realLP->changeMaxObj(c, R(_feasObj[c])); 4836 4837 assert(lowerReal(c) <= upperReal(c)); 4838 } 4839 4840 // remove last column 4841 _rationalLP->removeCol(numOrigCols); 4842 _realLP->removeCol(numOrigCols); 4843 _colTypes.reSize(numOrigCols); 4844 4845 // invalidate rational basis factorization 4846 _rationalLUSolver.clear(); 4847 4848 // print LP if in debug mode 4849 SPX_DEBUG(_realLP->writeFileLPBase("afterUntransFeas.lp", 0, 0, 0)); 4850 4851 // stop timing 4852 _statistics->transformTime->stop(); 4853 4854 #ifndef NDEBUG 4855 4856 if(sol._isPrimalFeasible) 4857 { 4858 VectorRational activity(numRowsRational()); 4859 _rationalLP->computePrimalActivity(sol._primal, activity); 4860 assert(sol._slacks == activity); 4861 } 4862 4863 #endif 4864 } 4865 4866 /** computes radius of infeasibility box implied by an approximate Farkas' proof 4867 4868 Given constraints of the form \f$ lhs <= Ax <= rhs \f$, a farkas proof y should satisfy \f$ y^T A = 0 \f$ and 4869 \f$ y_+^T lhs - y_-^T rhs > 0 \f$, where \f$ y_+, y_- \f$ denote the positive and negative parts of \f$ y \f$. 4870 If \f$ y \f$ is approximate, it may not satisfy \f$ y^T A = 0 \f$ exactly, but the proof is still valid as long 4871 as the following holds for all potentially feasible \f$ x \f$: 4872 4873 \f[ 4874 y^T Ax < (y_+^T lhs - y_-^T rhs) (*) 4875 \f] 4876 4877 we may therefore calculate \f$ y^T A \f$ and \f$ y_+^T lhs - y_-^T rhs \f$ exactly and check if the upper and lower 4878 bounds on \f$ x \f$ imply that all feasible \f$ x \f$ satisfy (*), and if not then compute bounds on \f$ x \f$ to 4879 guarantee (*). The simplest way to do this is to compute 4880 4881 \f[ 4882 B = (y_+^T lhs - y_-^T rhs) / \sum_i(|(y^T A)_i|) 4883 \f] 4884 4885 noting that if every component of \f$ x \f$ has \f$ |x_i| < B \f$, then (*) holds. 4886 4887 \f$ B \f$ can be increased by iteratively including variable bounds smaller than \f$ B \f$. The speed of this 4888 method can be further improved by using interval arithmetic for all computations. For related information see 4889 Sec. 4 of Neumaier and Shcherbina, Mathematical Programming A, 2004. 4890 4891 Set transformed to true if this method is called after _transformFeasibility(). 4892 */ 4893 template <class R> 4894 void SoPlexBase<R>::_computeInfeasBox(SolRational& sol, bool transformed) 4895 { 4896 assert(sol.hasDualFarkas()); 4897 4898 const VectorRational& lower = transformed ? _feasLower : lowerRational(); 4899 const VectorRational& upper = transformed ? _feasUpper : upperRational(); 4900 const VectorRational& lhs = transformed ? _feasLhs : lhsRational(); 4901 const VectorRational& rhs = transformed ? _feasRhs : rhsRational(); 4902 const VectorRational& y = sol._dualFarkas; 4903 4904 const int numRows = numRowsRational(); 4905 const int numCols = transformed ? numColsRational() - 1 : numColsRational(); 4906 4907 SSVectorRational ytransA(numColsRational()); 4908 Rational ytransb; 4909 Rational temp; 4910 4911 // prepare ytransA and ytransb; since we want exact arithmetic, we set the zero threshold of the semi-sparse 4912 // vector to zero 4913 ytransA.setTolerances(0); 4914 ytransA.clear(); 4915 ytransb = 0; 4916 4917 ///@todo this currently works only if all constraints are equations aggregate rows and sides using the multipliers of the Farkas ray 4918 for(int r = 0; r < numRows; r++) 4919 { 4920 ytransA += y[r] * _rationalLP->rowVector(r); 4921 ytransb += y[r] * (y[r] > 0 ? lhs[r] : rhs[r]); 4922 } 4923 4924 // if we work on the feasibility problem, we ignore the last column 4925 if(transformed) 4926 ytransA.reDim(numCols); 4927 4928 SPxOut::debug(this, "ytransb = {}\n", ytransb.str()); 4929 4930 // if we choose minus ytransb as vector of multipliers for the bound constraints on the variables, we obtain an 4931 // exactly feasible dual solution for the LP with zero objective function; we aggregate the bounds of the 4932 // variables accordingly and store its negation in temp 4933 temp = 0; 4934 bool isTempFinite = true; 4935 4936 for(int c = 0; c < numCols && isTempFinite; c++) 4937 { 4938 const Rational& minusRedCost = ytransA[c]; 4939 4940 if(minusRedCost > 0) 4941 { 4942 assert((upper[c] < _rationalPosInfty) == _upperFinite(_colTypes[c])); 4943 4944 if(_upperFinite(_colTypes[c])) 4945 temp += minusRedCost * upper[c]; 4946 else 4947 isTempFinite = false; 4948 } 4949 else if(minusRedCost < 0) 4950 { 4951 assert((lower[c] > _rationalNegInfty) == _lowerFinite(_colTypes[c])); 4952 4953 if(_lowerFinite(_colTypes[c])) 4954 temp += minusRedCost * lower[c]; 4955 else 4956 isTempFinite = false; 4957 } 4958 } 4959 4960 SPxOut::debug(this, "max ytransA*[x_l,x_u] = {}\n", (isTempFinite ? temp.str() : "infinite")); 4961 4962 // ytransb - temp is the increase in the dual objective along the Farkas ray; if this is positive, the dual is 4963 // unbounded and certifies primal infeasibility 4964 if(isTempFinite && temp < ytransb) 4965 { 4966 SPX_MSG_INFO1(spxout, spxout << "Farkas infeasibility proof verified exactly. (1)\n"); 4967 return; 4968 } 4969 4970 // ensure that array of nonzero elements in ytransA is available 4971 assert(ytransA.isSetup()); 4972 ytransA.setup(); 4973 4974 // if ytransb is negative, try to make it zero by including a positive lower bound or a negative upper bound 4975 if(ytransb < 0) 4976 { 4977 for(int c = 0; c < numCols; c++) 4978 { 4979 if(lower[c] > 0) 4980 { 4981 ytransA.setValue(c, ytransA[c] - ytransb / lower[c]); 4982 ytransb = 0; 4983 break; 4984 } 4985 else if(upper[c] < 0) 4986 { 4987 ytransA.setValue(c, ytransA[c] - ytransb / upper[c]); 4988 ytransb = 0; 4989 break; 4990 } 4991 } 4992 } 4993 4994 // if ytransb is still zero then the zero solution is inside the bounds and cannot be cut off by the Farkas 4995 // constraint; in this case, we cannot compute a Farkas box 4996 if(ytransb < 0) 4997 { 4998 SPX_MSG_INFO1(spxout, spxout << 4999 "Approximate Farkas proof to weak. Could not compute Farkas box. (1)\n"); 5000 return; 5001 } 5002 5003 // compute the one norm of ytransA 5004 temp = 0; 5005 const int size = ytransA.size(); 5006 5007 for(int n = 0; n < size; n++) 5008 temp += spxAbs(ytransA.value(n)); 5009 5010 // if the one norm is zero then ytransA is zero the Farkas proof should have been verified above 5011 assert(temp != 0); 5012 5013 // initialize variables in loop: size of Farkas box B, flag whether B has been increased, and number of current 5014 // nonzero in ytransA 5015 Rational B = ytransb / temp; 5016 bool success = false; 5017 int n = 0; 5018 5019 // loop through nonzeros of ytransA 5020 SPxOut::debug(this, "B = {}\n", B.str()); 5021 assert(ytransb >= 0); 5022 5023 while(true) 5024 { 5025 // if all nonzeros have been inspected once without increasing B, we abort; otherwise, we start another round 5026 if(n >= ytransA.size()) 5027 { 5028 if(!success) 5029 break; 5030 5031 success = false; 5032 n = 0; 5033 } 5034 5035 // get Farkas multiplier of the bound constraint as minus the value in ytransA 5036 const Rational& minusRedCost = ytransA.value(n); 5037 int colIdx = ytransA.index(n); 5038 5039 // if the multiplier is positive we inspect the lower bound: if it is finite and within the Farkas box, we can 5040 // increase B by including it in the Farkas proof 5041 assert((upper[colIdx] < _rationalPosInfty) == _upperFinite(_colTypes[colIdx])); 5042 assert((lower[colIdx] > _rationalNegInfty) == _lowerFinite(_colTypes[colIdx])); 5043 5044 if(minusRedCost < 0 && lower[colIdx] > -B && _lowerFinite(_colTypes[colIdx])) 5045 { 5046 ytransA.clearNum(n); 5047 ytransb -= minusRedCost * lower[colIdx]; 5048 temp += minusRedCost; 5049 5050 assert(ytransb >= 0); 5051 assert(temp >= 0); 5052 assert(temp == 0 || ytransb / temp > B); 5053 5054 // if ytransA and ytransb are zero, we have 0^T x >= 0 and cannot compute a Farkas box 5055 if(temp == 0 && ytransb == 0) 5056 { 5057 SPX_MSG_INFO1(spxout, spxout << 5058 "Approximate Farkas proof to weak. Could not compute Farkas box. (2)\n"); 5059 return; 5060 } 5061 // if ytransb is positive and ytransA is zero, we have 0^T x > 0, proving infeasibility 5062 else if(temp == 0) 5063 { 5064 assert(ytransb > 0); 5065 SPX_MSG_INFO1(spxout, spxout << "Farkas infeasibility proof verified exactly. (2)\n"); 5066 return; 5067 } 5068 else 5069 { 5070 B = ytransb / temp; 5071 SPxOut::debug(this, "B = {}\n", B.str()); 5072 } 5073 5074 success = true; 5075 } 5076 // if the multiplier is negative we inspect the upper bound: if it is finite and within the Farkas box, we can 5077 // increase B by including it in the Farkas proof 5078 else if(minusRedCost > 0 && upper[colIdx] < B && _upperFinite(_colTypes[colIdx])) 5079 { 5080 ytransA.clearNum(n); 5081 ytransb -= minusRedCost * upper[colIdx]; 5082 temp -= minusRedCost; 5083 5084 assert(ytransb >= 0); 5085 assert(temp >= 0); 5086 assert(temp == 0 || ytransb / temp > B); 5087 5088 // if ytransA and ytransb are zero, we have 0^T x >= 0 and cannot compute a Farkas box 5089 if(temp == 0 && ytransb == 0) 5090 { 5091 SPX_MSG_INFO1(spxout, spxout << 5092 "Approximate Farkas proof to weak. Could not compute Farkas box. (2)\n"); 5093 return; 5094 } 5095 // if ytransb is positive and ytransA is zero, we have 0^T x > 0, proving infeasibility 5096 else if(temp == 0) 5097 { 5098 assert(ytransb > 0); 5099 SPX_MSG_INFO1(spxout, spxout << "Farkas infeasibility proof verified exactly. (2)\n"); 5100 return; 5101 } 5102 else 5103 { 5104 B = ytransb / temp; 5105 SPxOut::debug(this, "B = {}\n", B.str()); 5106 } 5107 5108 success = true; 5109 } 5110 // the multiplier is zero, we can ignore the bound constraints on this variable 5111 else if(minusRedCost == 0) 5112 ytransA.clearNum(n); 5113 // currently this bound cannot be used to increase B; we will check it again in the next round, because B might 5114 // have increased by then 5115 else 5116 n++; 5117 } 5118 5119 if(B > 0) 5120 { 5121 SPX_MSG_INFO1(spxout, spxout << 5122 "Computed Farkas box: provably no feasible solutions with components less than " 5123 << B.str() << " in absolute value.\n"); 5124 } 5125 } 5126 5127 5128 5129 // get the last stable basis from the initial solver and store it as old basis, unsimplify basis if simplifier activated 5130 template <class R> 5131 void SoPlexBase<R>::_storeLastStableBasisBoosted(bool vanished) 5132 { 5133 if(_boostedSimplifier != 0) 5134 { 5135 // get solution vectors for transformed problem 5136 VectorBase<BP> tmpPrimal(vanished ? 0 : _boostedSolver.nCols()); 5137 VectorBase<BP> tmpSlacks(vanished ? 0 : _boostedSolver.nRows()); 5138 VectorBase<BP> tmpDual(vanished ? 0 : _boostedSolver.nRows()); 5139 VectorBase<BP> tmpRedCost(vanished ? 0 : _boostedSolver.nCols()); 5140 5141 if(!vanished) 5142 { 5143 _boostedSolver.getPrimalSol(tmpPrimal); 5144 _boostedSolver.getSlacks(tmpSlacks); 5145 _boostedSolver.getDualSol(tmpDual); 5146 _boostedSolver.getRedCostSol(tmpRedCost); 5147 5148 try 5149 { 5150 _boostedSimplifier->unsimplify(tmpPrimal, tmpDual, tmpSlacks, tmpRedCost, 5151 _boostedSolver.getOldBasisStatusRows().get_ptr(), _boostedSolver.getOldBasisStatusCols().get_ptr()); 5152 } 5153 catch(const SPxInternalCodeException& E) 5154 { 5155 // error message has already been printed 5156 throw SPxInternalCodeException("Error storing the basis"); 5157 } 5158 5159 // store basis for original problem 5160 _boostedSolver.getOldBasisStatusRows().reSize(numRowsRational()); 5161 _boostedSolver.getOldBasisStatusCols().reSize(numColsRational()); 5162 _boostedSimplifier->getBasis(_boostedSolver.getOldBasisStatusRows().get_ptr(), 5163 _boostedSolver.getOldBasisStatusCols().get_ptr(), 5164 _boostedSolver.getOldBasisStatusRows().size(), _boostedSolver.getOldBasisStatusCols().size()); 5165 } 5166 } 5167 5168 // store last basis as old basis 5169 _storeBasisAsOldBasisBoosted(_boostedSolver.getOldBasisStatusRows(), 5170 _boostedSolver.getOldBasisStatusCols()); 5171 } 5172 5173 // get the last stable basis from the initial solver and store it as old basis, unsimplify basis if simplifier activated 5174 template <class R> 5175 void SoPlexBase<R>::_storeLastStableBasis(bool vanished) 5176 { 5177 if(_simplifier != 0) 5178 { 5179 // get solution vectors for transformed problem 5180 VectorBase<R> tmpPrimal(vanished ? 0 : _solver.nCols()); 5181 VectorBase<R> tmpSlacks(vanished ? 0 : _solver.nRows()); 5182 VectorBase<R> tmpDual(vanished ? 0 : _solver.nRows()); 5183 VectorBase<R> tmpRedCost(vanished ? 0 : _solver.nCols()); 5184 5185 if(!vanished) 5186 { 5187 _solver.getPrimalSol(tmpPrimal); 5188 _solver.getSlacks(tmpSlacks); 5189 _solver.getDualSol(tmpDual); 5190 _solver.getRedCostSol(tmpRedCost); 5191 5192 // unscale vectors 5193 if(_scaler && _solver.isScaled()) 5194 { 5195 _scaler->unscalePrimal(_solver, tmpPrimal); 5196 _scaler->unscaleSlacks(_solver, tmpSlacks); 5197 _scaler->unscaleDual(_solver, tmpDual); 5198 _scaler->unscaleRedCost(_solver, tmpRedCost); 5199 } 5200 5201 _solver.getBasis(_solver.getOldBasisStatusRows().get_ptr(), 5202 _solver.getOldBasisStatusCols().get_ptr(), _solver.getOldBasisStatusRows().size(), 5203 _solver.getOldBasisStatusCols().size()); 5204 5205 try 5206 { 5207 _simplifier->unsimplify(tmpPrimal, tmpDual, tmpSlacks, tmpRedCost, 5208 _solver.getOldBasisStatusRows().get_ptr(), _solver.getOldBasisStatusCols().get_ptr()); 5209 } 5210 catch(const SPxInternalCodeException& E) 5211 { 5212 // error message has already been printed 5213 throw SPxInternalCodeException("Error storing the basis"); 5214 } 5215 5216 // store basis for original problem 5217 _solver.getOldBasisStatusRows().reSize(numRowsRational()); 5218 _solver.getOldBasisStatusCols().reSize(numColsRational()); 5219 _simplifier->getBasis(_solver.getOldBasisStatusRows().get_ptr(), 5220 _solver.getOldBasisStatusCols().get_ptr(), 5221 _solver.getOldBasisStatusRows().size(), _solver.getOldBasisStatusCols().size()); 5222 } 5223 } 5224 5225 // store last basis as old basis 5226 _storeBasisAsOldBasis(_solver.getOldBasisStatusRows(), _solver.getOldBasisStatusCols()); 5227 } 5228 5229 5230 5231 // General specializations 5232 /// solves real LP during iterative refinement 5233 ///@param fromscratch in 5234 ///@param primal out 5235 ///@param dual out 5236 ///@param basisStatusRows out 5237 ///@param basisStatusRows out 5238 template <class R> 5239 typename SPxSolverBase<R>::Status SoPlexBase<R>::_solveRealForRational(bool fromscratch, 5240 VectorBase<R>& primal, VectorBase<R>& dual, 5241 DataArray< typename SPxSolverBase<R>::VarStatus >& basisStatusRows, 5242 DataArray< typename SPxSolverBase<R>::VarStatus >& basisStatusCols) 5243 { 5244 assert(_isConsistent()); 5245 5246 assert(_solver.nRows() == numRowsRational()); 5247 assert(_solver.nCols() == numColsRational()); 5248 assert(primal.dim() == numColsRational()); 5249 assert(dual.dim() == numRowsRational()); 5250 5251 typename SPxSolverBase<R>::Status result = SPxSolverBase<R>::UNKNOWN; 5252 5253 #ifndef SOPLEX_MANUAL_ALT 5254 5255 if(fromscratch || !_hasBasis) 5256 _enableSimplifierAndScaler(); 5257 else 5258 _disableSimplifierAndScaler(); 5259 5260 #else 5261 _disableSimplifierAndScaler(); 5262 #endif 5263 5264 // reset basis to slack basis when solving from scratch 5265 if(fromscratch) 5266 _solver.reLoad(); 5267 5268 // start timing 5269 _statistics->syncTime->start(); 5270 5271 // if preprocessing is applied, we need to restore the original LP at the end 5272 SPxLPRational* rationalLP = 0; 5273 5274 if(_simplifier != 0 || _scaler != nullptr) 5275 { 5276 spx_alloc(rationalLP); 5277 rationalLP = new(rationalLP) SPxLPRational(_solver); 5278 } 5279 5280 // with preprocessing or solving from scratch, the basis may change, hence invalidate the 5281 // rational basis factorization 5282 if(_simplifier != nullptr || _scaler != nullptr || fromscratch) 5283 _rationalLUSolver.clear(); 5284 5285 // stop timing 5286 _statistics->syncTime->stop(); 5287 5288 try 5289 { 5290 // apply problem simplification 5291 typename SPxSimplifier<R>::Result simplificationStatus = SPxSimplifier<R>::OKAY; 5292 5293 int oldIterations = _statistics->iterations; // save previous number of iterations 5294 5295 if(_simplifier != 0) 5296 { 5297 // do not remove bounds of boxed variables or sides of ranged rows if bound flipping is used 5298 bool keepbounds = intParam(SoPlexBase<R>::RATIOTESTER) == SoPlexBase<R>::RATIOTESTER_BOUNDFLIPPING; 5299 Real remainingTime = _solver.getMaxTime() - _solver.time(); 5300 simplificationStatus = _simplifier->simplify(_solver, remainingTime, keepbounds, 5301 _solver.random.getSeed()); 5302 } 5303 5304 // apply scaling after the simplification 5305 if(_scaler != nullptr && simplificationStatus == SPxSimplifier<R>::OKAY) 5306 _scaler->scale(_solver, false); 5307 5308 // run the simplex method if problem has not been solved by the simplifier 5309 if(simplificationStatus == SPxSimplifier<R>::OKAY) 5310 { 5311 SPX_MSG_INFO1(spxout, spxout << std::endl); 5312 5313 _solveRealLPAndRecordStatistics(); 5314 5315 SPX_MSG_INFO1(spxout, spxout << std::endl); 5316 } 5317 5318 ///@todo move to private helper methods 5319 // evaluate status flag 5320 if(simplificationStatus == SPxSimplifier<R>::INFEASIBLE) 5321 result = SPxSolverBase<R>::INFEASIBLE; 5322 else if(simplificationStatus == SPxSimplifier<R>::DUAL_INFEASIBLE) 5323 result = SPxSolverBase<R>::INForUNBD; 5324 else if(simplificationStatus == SPxSimplifier<R>::UNBOUNDED) 5325 result = SPxSolverBase<R>::UNBOUNDED; 5326 else if(simplificationStatus == SPxSimplifier<R>::VANISHED 5327 || simplificationStatus == SPxSimplifier<R>::OKAY) 5328 { 5329 result = simplificationStatus == SPxSimplifier<R>::VANISHED ? SPxSolverBase<R>::OPTIMAL : 5330 _solver.status(); 5331 5332 ///@todo move to private helper methods 5333 // process result 5334 switch(result) 5335 { 5336 case SPxSolverBase<R>::OPTIMAL: 5337 5338 // unsimplify if simplifier is active and LP is solved to optimality; this must be done here and not at solution 5339 // query, because we want to have the basis for the original problem 5340 if(_simplifier != 0) 5341 { 5342 assert(!_simplifier->isUnsimplified()); 5343 assert(simplificationStatus == SPxSimplifier<R>::VANISHED 5344 || simplificationStatus == SPxSimplifier<R>::OKAY); 5345 5346 bool vanished = simplificationStatus == SPxSimplifier<R>::VANISHED; 5347 5348 // get solution vectors for transformed problem 5349 VectorBase<R> tmpPrimal(vanished ? 0 : _solver.nCols()); 5350 VectorBase<R> tmpSlacks(vanished ? 0 : _solver.nRows()); 5351 VectorBase<R> tmpDual(vanished ? 0 : _solver.nRows()); 5352 VectorBase<R> tmpRedCost(vanished ? 0 : _solver.nCols()); 5353 5354 if(!vanished) 5355 { 5356 assert(_solver.status() == SPxSolverBase<R>::OPTIMAL); 5357 5358 _solver.getPrimalSol(tmpPrimal); 5359 _solver.getSlacks(tmpSlacks); 5360 _solver.getDualSol(tmpDual); 5361 _solver.getRedCostSol(tmpRedCost); 5362 5363 // unscale vectors 5364 if(_scaler != nullptr) 5365 { 5366 _scaler->unscalePrimal(_solver, tmpPrimal); 5367 _scaler->unscaleSlacks(_solver, tmpSlacks); 5368 _scaler->unscaleDual(_solver, tmpDual); 5369 _scaler->unscaleRedCost(_solver, tmpRedCost); 5370 } 5371 5372 // get basis of transformed problem 5373 basisStatusRows.reSize(_solver.nRows()); 5374 basisStatusCols.reSize(_solver.nCols()); 5375 _solver.getBasis(basisStatusRows.get_ptr(), basisStatusCols.get_ptr(), basisStatusRows.size(), 5376 basisStatusCols.size()); 5377 } 5378 5379 try 5380 { 5381 _simplifier->unsimplify(tmpPrimal, tmpDual, tmpSlacks, tmpRedCost, basisStatusRows.get_ptr(), 5382 basisStatusCols.get_ptr()); 5383 } 5384 catch(const SPxInternalCodeException& E) 5385 { 5386 // error message has already been printed 5387 throw SPxInternalCodeException("Error storing the basis"); 5388 } 5389 5390 // store basis for original problem 5391 basisStatusRows.reSize(numRowsRational()); 5392 basisStatusCols.reSize(numColsRational()); 5393 _simplifier->getBasis(basisStatusRows.get_ptr(), basisStatusCols.get_ptr(), basisStatusRows.size(), 5394 basisStatusCols.size()); 5395 _hasBasis = true; 5396 5397 primal = _simplifier->unsimplifiedPrimal(); 5398 dual = _simplifier->unsimplifiedDual(); 5399 } 5400 else 5401 { 5402 _solver.getPrimalSol(primal); 5403 _solver.getDualSol(dual); 5404 5405 // unscale vectors 5406 if(_scaler != nullptr) 5407 { 5408 _scaler->unscalePrimal(_solver, primal); 5409 _scaler->unscaleDual(_solver, dual); 5410 } 5411 5412 // get basis of transformed problem 5413 basisStatusRows.reSize(_solver.nRows()); 5414 basisStatusCols.reSize(_solver.nCols()); 5415 _solver.getBasis(basisStatusRows.get_ptr(), basisStatusCols.get_ptr(), basisStatusRows.size(), 5416 basisStatusCols.size()); 5417 _hasBasis = true; 5418 } 5419 5420 if(oldIterations < _statistics->iterations || _oldBasisStatusRows.size() != basisStatusRows.size() 5421 || _oldBasisStatusCols.size() != basisStatusCols.size()) 5422 _storeBasisAsOldBasis(basisStatusRows, basisStatusCols); 5423 5424 break; 5425 5426 case SPxSolverBase<R>::ABORT_CYCLING: 5427 if(_simplifier == 0 && boolParam(SoPlexBase<R>::ACCEPTCYCLING)) 5428 { 5429 _solver.getPrimalSol(primal); 5430 _solver.getDualSol(dual); 5431 5432 // unscale vectors 5433 if(_scaler != nullptr) 5434 { 5435 _scaler->unscalePrimal(_solver, primal); 5436 _scaler->unscaleDual(_solver, dual); 5437 } 5438 } 5439 5440 // get the last stable basis. The hope is that precision boosting will get rid of cycling. 5441 if(boolParam(SoPlexBase<R>::STORE_BASIS_BEFORE_SIMPLEX_PIVOT)) 5442 { 5443 try 5444 { 5445 if(oldIterations < _statistics->iterations || _oldBasisStatusRows.size() != basisStatusRows.size() 5446 || _oldBasisStatusCols.size() != basisStatusCols.size()) 5447 _storeLastStableBasis(simplificationStatus == SPxSimplifier<R>::VANISHED); 5448 } 5449 catch(const SPxInternalCodeException& E) 5450 { 5451 SPX_MSG_INFO1(spxout, spxout << "Caught exception <" << E.what() << 5452 "> while processing the result of the solve.\n"); 5453 SPX_MSG_INFO1(spxout, spxout << "Storage of last basis failed. Keep going.\n"); 5454 } 5455 } 5456 5457 break; 5458 5459 // intentional fallthrough 5460 case SPxSolverBase<R>::ABORT_TIME: 5461 case SPxSolverBase<R>::ABORT_ITER: 5462 case SPxSolverBase<R>::ABORT_VALUE: 5463 case SPxSolverBase<R>::REGULAR: 5464 case SPxSolverBase<R>::RUNNING: 5465 case SPxSolverBase<R>::UNBOUNDED: 5466 _hasBasis = (_solver.basis().status() > SPxBasisBase<R>::NO_PROBLEM); 5467 5468 if(_hasBasis && _simplifier == 0) 5469 { 5470 basisStatusRows.reSize(_solver.nRows()); 5471 basisStatusCols.reSize(_solver.nCols()); 5472 _solver.getBasis(basisStatusRows.get_ptr(), basisStatusCols.get_ptr(), basisStatusRows.size(), 5473 basisStatusCols.size()); 5474 5475 if(oldIterations < _statistics->iterations || _oldBasisStatusRows.size() != basisStatusRows.size() 5476 || _oldBasisStatusCols.size() != basisStatusCols.size()) 5477 _storeBasisAsOldBasis(basisStatusRows, basisStatusCols); 5478 } 5479 else 5480 { 5481 _hasBasis = false; 5482 _rationalLUSolver.clear(); 5483 } 5484 5485 // get the last stable basis. 5486 if(boolParam(SoPlexBase<R>::STORE_BASIS_BEFORE_SIMPLEX_PIVOT)) 5487 { 5488 try 5489 { 5490 if(oldIterations < _statistics->iterations || _oldBasisStatusRows.size() != basisStatusRows.size() 5491 || _oldBasisStatusCols.size() != basisStatusCols.size()) 5492 _storeLastStableBasis(simplificationStatus == SPxSimplifier<R>::VANISHED); 5493 } 5494 catch(const SPxInternalCodeException& E) 5495 { 5496 SPX_MSG_INFO1(spxout, spxout << "Caught exception <" << E.what() << 5497 "> while processing the result of the solve.\n"); 5498 SPX_MSG_INFO1(spxout, spxout << "Storage of last basis failed. Keep going.\n"); 5499 } 5500 } 5501 5502 break; 5503 5504 case SPxSolverBase<R>::INFEASIBLE: 5505 5506 // if simplifier is active we can currently not return a Farkas ray or basis 5507 if(_simplifier != 0) 5508 { 5509 _hasBasis = false; 5510 _rationalLUSolver.clear(); 5511 5512 // get the last stable basis. 5513 if(boolParam(SoPlexBase<R>::STORE_BASIS_BEFORE_SIMPLEX_PIVOT)) 5514 { 5515 try 5516 { 5517 if(oldIterations < _statistics->iterations || _oldBasisStatusRows.size() != basisStatusRows.size() 5518 || _oldBasisStatusCols.size() != basisStatusCols.size()) 5519 _storeLastStableBasis(simplificationStatus == SPxSimplifier<R>::VANISHED); 5520 } 5521 catch(const SPxInternalCodeException& E) 5522 { 5523 SPX_MSG_INFO1(spxout, spxout << "Caught exception <" << E.what() << 5524 "> while processing the result of the solve.\n"); 5525 SPX_MSG_INFO1(spxout, spxout << "Storage of last basis failed. Keep going.\n"); 5526 } 5527 } 5528 5529 break; 5530 } 5531 5532 // return Farkas ray as dual solution 5533 _solver.getDualfarkas(dual); 5534 5535 // unscale vectors 5536 if(_scaler != nullptr) 5537 _scaler->unscaleDual(_solver, dual); 5538 5539 // get basis of transformed problem 5540 basisStatusRows.reSize(_solver.nRows()); 5541 basisStatusCols.reSize(_solver.nCols()); 5542 _solver.getBasis(basisStatusRows.get_ptr(), basisStatusCols.get_ptr(), basisStatusRows.size(), 5543 basisStatusCols.size()); 5544 _hasBasis = true; 5545 5546 // if possible, get the last stable basis, otherwise store the infeasible basis. 5547 if(boolParam(SoPlexBase<R>::STORE_BASIS_BEFORE_SIMPLEX_PIVOT)) 5548 { 5549 try 5550 { 5551 if(oldIterations < _statistics->iterations || _oldBasisStatusRows.size() != basisStatusRows.size() 5552 || _oldBasisStatusCols.size() != basisStatusCols.size()) 5553 _storeLastStableBasis(simplificationStatus == SPxSimplifier<R>::VANISHED); 5554 } 5555 catch(const SPxInternalCodeException& E) 5556 { 5557 SPX_MSG_INFO1(spxout, spxout << "Caught exception <" << E.what() << 5558 "> while processing the result of the solve.\n"); 5559 SPX_MSG_INFO1(spxout, spxout << "Storage of last basis failed. Keep going.\n"); 5560 } 5561 } 5562 else 5563 { 5564 if(oldIterations < _statistics->iterations || _oldBasisStatusRows.size() != basisStatusRows.size() 5565 || _oldBasisStatusCols.size() != basisStatusCols.size()) 5566 _storeBasisAsOldBasis(basisStatusRows, basisStatusCols); 5567 } 5568 5569 break; 5570 5571 case SPxSolverBase<R>::INForUNBD: 5572 case SPxSolverBase<R>::SINGULAR: 5573 default: 5574 _hasBasis = false; 5575 _rationalLUSolver.clear(); 5576 5577 // get the last stable basis. 5578 if(boolParam(SoPlexBase<R>::STORE_BASIS_BEFORE_SIMPLEX_PIVOT)) 5579 { 5580 try 5581 { 5582 if(oldIterations < _statistics->iterations || _oldBasisStatusRows.size() != basisStatusRows.size() 5583 || _oldBasisStatusCols.size() != basisStatusCols.size()) 5584 _storeLastStableBasis(simplificationStatus == SPxSimplifier<R>::VANISHED); 5585 } 5586 catch(const SPxInternalCodeException& E) 5587 { 5588 SPX_MSG_INFO1(spxout, spxout << "Caught exception <" << E.what() << 5589 "> while processing the result of the solve.\n"); 5590 SPX_MSG_INFO1(spxout, spxout << "Storage of last basis failed. Keep going.\n"); 5591 } 5592 } 5593 5594 break; 5595 } 5596 } 5597 } 5598 catch(...) 5599 { 5600 SPX_MSG_INFO1(spxout, spxout << "Exception thrown during floating-point solve.\n"); 5601 result = SPxSolverBase<R>::ERROR; 5602 _hasBasis = false; 5603 _rationalLUSolver.clear(); 5604 5605 } 5606 5607 // restore original LP if necessary 5608 if(_simplifier != 0 || _scaler != nullptr) 5609 { 5610 assert(rationalLP != 0); 5611 _solver.loadLP((SPxLPBase<R>)(*rationalLP)); 5612 rationalLP->~SPxLPRational(); 5613 spx_free(rationalLP); 5614 5615 if(_hasBasis) 5616 _solver.setBasis(basisStatusRows.get_ptr(), basisStatusCols.get_ptr()); 5617 } 5618 5619 return result; 5620 } 5621 5622 5623 5624 /// solves real LP with recovery mechanism 5625 template <class R> 5626 typename SPxSolverBase<R>::Status SoPlexBase<R>::_solveRealStable(bool acceptUnbounded, 5627 bool acceptInfeasible, VectorBase<R>& primal, VectorBase<R>& dual, 5628 DataArray< typename SPxSolverBase<R>::VarStatus >& basisStatusRows, 5629 DataArray< typename SPxSolverBase<R>::VarStatus >& basisStatusCols, 5630 const bool forceNoSimplifier) 5631 { 5632 typename SPxSolverBase<R>::Status result = SPxSolverBase<R>::UNKNOWN; 5633 5634 bool fromScratch = false; 5635 bool solved = false; 5636 bool solvedFromScratch = false; 5637 bool initialSolve = true; 5638 bool increasedMarkowitz = false; 5639 bool relaxedTolerances = false; 5640 bool tightenedTolerances = false; 5641 bool switchedScaler = false; 5642 bool switchedSimplifier = false; 5643 bool switchedRatiotester = false; 5644 bool switchedPricer = false; 5645 bool turnedoffPre = false; 5646 5647 R markowitz = _slufactor.markowitz(); 5648 int ratiotester = intParam(SoPlexBase<R>::RATIOTESTER); 5649 int pricer = intParam(SoPlexBase<R>::PRICER); 5650 int simplifier = intParam(SoPlexBase<R>::SIMPLIFIER); 5651 int scaler = intParam(SoPlexBase<R>::SCALER); 5652 int type = intParam(SoPlexBase<R>::ALGORITHM); 5653 5654 if(forceNoSimplifier) 5655 setIntParam(SoPlexBase<R>::SIMPLIFIER, SoPlexBase<R>::SIMPLIFIER_OFF); 5656 5657 while(true) 5658 { 5659 assert(!increasedMarkowitz || GE(_slufactor.markowitz(), R(0.9), this->tolerances()->epsilon())); 5660 5661 result = _solveRealForRational(fromScratch, primal, dual, basisStatusRows, basisStatusCols); 5662 5663 solved = (result == SPxSolverBase<R>::OPTIMAL) 5664 || (result == SPxSolverBase<R>::INFEASIBLE && acceptInfeasible) 5665 || (result == SPxSolverBase<R>::UNBOUNDED && acceptUnbounded); 5666 5667 if(solved || result == SPxSolverBase<R>::ABORT_TIME || result == SPxSolverBase<R>::ABORT_ITER) 5668 break; 5669 5670 if(initialSolve) 5671 { 5672 SPX_MSG_INFO1(spxout, spxout << "Numerical troubles during floating-point solve." << std::endl); 5673 initialSolve = false; 5674 } 5675 5676 if(!turnedoffPre 5677 && (intParam(SoPlexBase<R>::SIMPLIFIER) != SoPlexBase<R>::SIMPLIFIER_OFF 5678 || intParam(SoPlexBase<R>::SCALER) != SoPlexBase<R>::SCALER_OFF)) 5679 { 5680 SPX_MSG_INFO1(spxout, spxout << "Turning off preprocessing." << std::endl); 5681 5682 turnedoffPre = true; 5683 5684 setIntParam(SoPlexBase<R>::SCALER, SoPlexBase<R>::SCALER_OFF); 5685 setIntParam(SoPlexBase<R>::SIMPLIFIER, SoPlexBase<R>::SIMPLIFIER_OFF); 5686 5687 fromScratch = true; 5688 solvedFromScratch = true; 5689 continue; 5690 } 5691 5692 setIntParam(SoPlexBase<R>::SCALER, scaler); 5693 setIntParam(SoPlexBase<R>::SIMPLIFIER, simplifier); 5694 5695 if(!increasedMarkowitz) 5696 { 5697 SPX_MSG_INFO1(spxout, spxout << "Increasing Markowitz threshold." << std::endl); 5698 5699 _slufactor.setMarkowitz(0.9); 5700 increasedMarkowitz = true; 5701 5702 try 5703 { 5704 _solver.factorize(); 5705 continue; 5706 } 5707 catch(...) 5708 { 5709 SPxOut::debug(this, "Factorization failed.\n"); 5710 } 5711 } 5712 5713 if(!solvedFromScratch) 5714 { 5715 SPX_MSG_INFO1(spxout, spxout << "Solving from scratch." << std::endl); 5716 5717 fromScratch = true; 5718 solvedFromScratch = true; 5719 continue; 5720 } 5721 5722 setIntParam(SoPlexBase<R>::RATIOTESTER, ratiotester); 5723 setIntParam(SoPlexBase<R>::PRICER, pricer); 5724 5725 if(!switchedScaler) 5726 { 5727 SPX_MSG_INFO1(spxout, spxout << "Switching scaling." << std::endl); 5728 5729 if(scaler == int(SoPlexBase<R>::SCALER_OFF)) 5730 setIntParam(SoPlexBase<R>::SCALER, SoPlexBase<R>::SCALER_BIEQUI); 5731 else 5732 setIntParam(SoPlexBase<R>::SCALER, SoPlexBase<R>::SCALER_OFF); 5733 5734 fromScratch = true; 5735 solvedFromScratch = true; 5736 switchedScaler = true; 5737 continue; 5738 } 5739 5740 if(!switchedSimplifier && !forceNoSimplifier) 5741 { 5742 SPX_MSG_INFO1(spxout, spxout << "Switching simplification." << std::endl); 5743 5744 if(simplifier == int(SoPlexBase<R>::SIMPLIFIER_OFF)) 5745 setIntParam(SoPlexBase<R>::SIMPLIFIER, SoPlexBase<R>::SIMPLIFIER_INTERNAL); 5746 else 5747 setIntParam(SoPlexBase<R>::SIMPLIFIER, SoPlexBase<R>::SIMPLIFIER_OFF); 5748 5749 fromScratch = true; 5750 solvedFromScratch = true; 5751 switchedSimplifier = true; 5752 continue; 5753 } 5754 5755 setIntParam(SoPlexBase<R>::SIMPLIFIER, SoPlexBase<R>::SIMPLIFIER_OFF); 5756 5757 if(!relaxedTolerances) 5758 { 5759 SPX_MSG_INFO1(spxout, spxout << "Relaxing tolerances." << std::endl); 5760 5761 setIntParam(SoPlexBase<R>::ALGORITHM, SoPlexBase<R>::ALGORITHM_PRIMAL); 5762 5763 // scale tols by up 1e3 but so that they do not exceed 1e-3 5764 if(_solver.entertol() > 1e-6) 5765 _solver.scaleEntertol(1e-3 / _solver.entertol()); 5766 else 5767 _solver.scaleEntertol(1e3); 5768 5769 if(_solver.leavetol() > 1e-6) 5770 _solver.scaleLeavetol(1e-3 / _solver.leavetol()); 5771 else 5772 _solver.scaleLeavetol(1e3); 5773 5774 relaxedTolerances = _solver.entertol() >= 1e-3; 5775 solvedFromScratch = false; 5776 continue; 5777 } 5778 5779 if(!tightenedTolerances && result != SPxSolverBase<R>::INFEASIBLE) 5780 { 5781 SPX_MSG_INFO1(spxout, spxout << "Tightening tolerances." << std::endl); 5782 5783 setIntParam(SoPlexBase<R>::ALGORITHM, SoPlexBase<R>::ALGORITHM_DUAL); 5784 5785 // scale tols by up 1e-3 but so that they are not smaller than 1e-9 5786 if(_solver.entertol() < 1e-6) 5787 _solver.scaleEntertol(1e-9 / _solver.entertol()); 5788 else 5789 _solver.scaleEntertol(1e-3); 5790 5791 tightenedTolerances = _solver.entertol() <= 1e-9; 5792 solvedFromScratch = false; 5793 continue; 5794 } 5795 5796 setIntParam(SoPlexBase<R>::ALGORITHM, type); 5797 5798 if(!switchedRatiotester) 5799 { 5800 SPX_MSG_INFO1(spxout, spxout << "Switching ratio test." << std::endl); 5801 5802 _solver.setType(_solver.type() == SPxSolverBase<R>::LEAVE ? SPxSolverBase<R>::ENTER : 5803 SPxSolverBase<R>::LEAVE); 5804 5805 if(_solver.ratiotester() != (SPxRatioTester<R>*)&_ratiotesterTextbook) 5806 setIntParam(SoPlexBase<R>::RATIOTESTER, RATIOTESTER_TEXTBOOK); 5807 else 5808 setIntParam(SoPlexBase<R>::RATIOTESTER, RATIOTESTER_FAST); 5809 5810 switchedRatiotester = true; 5811 solvedFromScratch = false; 5812 continue; 5813 } 5814 5815 if(!switchedPricer) 5816 { 5817 SPX_MSG_INFO1(spxout, spxout << "Switching pricer." << std::endl); 5818 5819 _solver.setType(_solver.type() == SPxSolverBase<R>::LEAVE ? SPxSolverBase<R>::ENTER : 5820 SPxSolverBase<R>::LEAVE); 5821 5822 if(_solver.pricer() != (SPxPricer<R>*)&_pricerDevex) 5823 setIntParam(SoPlexBase<R>::PRICER, PRICER_DEVEX); 5824 else 5825 setIntParam(SoPlexBase<R>::PRICER, PRICER_STEEP); 5826 5827 switchedPricer = true; 5828 solvedFromScratch = false; 5829 continue; 5830 } 5831 5832 SPX_MSG_INFO1(spxout, spxout << "Giving up." << std::endl); 5833 5834 break; 5835 } 5836 5837 _slufactor.setMarkowitz(markowitz); 5838 5839 setIntParam(SoPlexBase<R>::RATIOTESTER, ratiotester); 5840 setIntParam(SoPlexBase<R>::PRICER, pricer); 5841 setIntParam(SoPlexBase<R>::SIMPLIFIER, simplifier); 5842 setIntParam(SoPlexBase<R>::SCALER, scaler); 5843 setIntParam(SoPlexBase<R>::ALGORITHM, type); 5844 5845 return result; 5846 } 5847 5848 5849 5850 // General specializations 5851 /// solves real LP during iterative refinement 5852 ///@param initialSolve in 5853 ///@param primal out 5854 ///@param dual out 5855 ///@param basisStatusRows out 5856 ///@param basisStatusRows out 5857 ///@param boostedResult out 5858 template <class R> 5859 void SoPlexBase<R>::_solveRealForRationalBoosted( 5860 VectorBase<BP>& primal, VectorBase<BP>& dual, 5861 DataArray< typename SPxSolverBase<R>::VarStatus >& basisStatusRows, 5862 DataArray< typename SPxSolverBase<R>::VarStatus >& basisStatusCols, 5863 typename SPxSolverBase<BP>::Status& boostedResult, bool initialSolve) 5864 { 5865 assert(boolParam(SoPlexBase<R>::PRECISION_BOOSTING)); 5866 assert(_isBoostedConsistent()); 5867 5868 assert(_boostedSolver.nRows() == numRowsRational()); 5869 assert(_boostedSolver.nCols() == numColsRational()); 5870 assert(primal.dim() == numColsRational()); 5871 assert(dual.dim() == numRowsRational()); 5872 5873 #ifdef SOPLEX_WITH_MPFR 5874 5875 #ifndef SOPLEX_MANUAL_ALT 5876 5877 if(_isBoostedStartingFromSlack(initialSolve)) 5878 _enableSimplifierAndScaler(); 5879 else 5880 _disableSimplifierAndScaler(); 5881 5882 #else 5883 _disableSimplifierAndScaler(); 5884 #endif 5885 5886 // start timing 5887 _statistics->syncTime->start(); 5888 5889 // if preprocessing is applied, we need to restore the original LP at the end 5890 SPxLPRational* rationalLP = 0; 5891 5892 if(_boostedSimplifier != 0 || _boostedScaler != nullptr) 5893 { 5894 spx_alloc(rationalLP); 5895 rationalLP = new(rationalLP) SPxLPRational(_boostedSolver); 5896 } 5897 5898 // with preprocessing or solving from scratch, the basis may change, hence invalidate the 5899 // rational basis factorization 5900 if(_boostedSimplifier != nullptr || _boostedScaler != nullptr) 5901 _rationalLUSolver.clear(); 5902 5903 // stop timing 5904 _statistics->syncTime->stop(); 5905 5906 try 5907 { 5908 // apply problem simplification 5909 typename SPxSimplifier<BP>::Result simplificationStatus = SPxSimplifier<BP>::OKAY; 5910 5911 int oldIterations = _statistics->iterations; // save previous number of iterations 5912 5913 if(_boostedSimplifier != 0) 5914 { 5915 // do not remove bounds of boxed variables or sides of ranged rows if bound flipping is used 5916 bool keepbounds = intParam(SoPlexBase<R>::RATIOTESTER) == SoPlexBase<R>::RATIOTESTER_BOUNDFLIPPING; 5917 Real remainingTime = _boostedSolver.getMaxTime() - _boostedSolver.time(); 5918 BP tol = pow(10, -(int)(BP::default_precision() * _tolPrecisionRatio)); 5919 simplificationStatus = _boostedSimplifier->simplify(_boostedSolver, remainingTime, keepbounds, 5920 _boostedSolver.random.getSeed()); 5921 } 5922 5923 // apply scaling after the simplification 5924 if(_boostedScaler != nullptr && simplificationStatus == SPxSimplifier<BP>::OKAY) 5925 _boostedScaler->scale(_boostedSolver, false); 5926 5927 // run the simplex method if problem has not been solved by the simplifier 5928 if(simplificationStatus == SPxSimplifier<BP>::OKAY) 5929 { 5930 SPX_MSG_INFO1(spxout, spxout << std::endl); 5931 5932 _solveBoostedRealLPAndRecordStatistics(); 5933 5934 SPX_MSG_INFO1(spxout, spxout << std::endl); 5935 } 5936 5937 ///@todo move to private helper methods 5938 // evaluate status flag 5939 if(simplificationStatus == SPxSimplifier<BP>::INFEASIBLE) 5940 boostedResult = SPxSolverBase<BP>::INFEASIBLE; 5941 else if(simplificationStatus == SPxSimplifier<BP>::DUAL_INFEASIBLE) 5942 boostedResult = SPxSolverBase<BP>::INForUNBD; 5943 else if(simplificationStatus == SPxSimplifier<BP>::UNBOUNDED) 5944 boostedResult = SPxSolverBase<BP>::UNBOUNDED; 5945 else if(simplificationStatus == SPxSimplifier<BP>::VANISHED 5946 || simplificationStatus == SPxSimplifier<BP>::OKAY) 5947 { 5948 boostedResult = simplificationStatus == SPxSimplifier<BP>::VANISHED ? SPxSolverBase<BP>::OPTIMAL : 5949 _boostedSolver.status(); 5950 5951 ///@todo move to private helper methods 5952 // process boostedResult 5953 switch(boostedResult) 5954 { 5955 case SPxSolverBase<BP>::OPTIMAL: 5956 5957 // unsimplify if simplifier is active and LP is solved to optimality; this must be done here and not at solution 5958 // query, because we want to have the basis for the original problem 5959 if(_boostedSimplifier != 0) 5960 { 5961 assert(!_boostedSimplifier->isUnsimplified()); 5962 assert(simplificationStatus == SPxSimplifier<BP>::VANISHED 5963 || simplificationStatus == SPxSimplifier<BP>::OKAY); 5964 5965 bool vanished = simplificationStatus == SPxSimplifier<BP>::VANISHED; 5966 5967 // get solution vectors for transformed problem 5968 VectorBase<BP> tmpPrimal(vanished ? 0 : _boostedSolver.nCols()); 5969 VectorBase<BP> tmpSlacks(vanished ? 0 : _boostedSolver.nRows()); 5970 VectorBase<BP> tmpDual(vanished ? 0 : _boostedSolver.nRows()); 5971 VectorBase<BP> tmpRedCost(vanished ? 0 : _boostedSolver.nCols()); 5972 5973 if(!vanished) 5974 { 5975 assert(_boostedSolver.status() == SPxSolverBase<BP>::OPTIMAL); 5976 5977 _boostedSolver.getPrimalSol(tmpPrimal); 5978 _boostedSolver.getSlacks(tmpSlacks); 5979 _boostedSolver.getDualSol(tmpDual); 5980 _boostedSolver.getRedCostSol(tmpRedCost); 5981 5982 // unscale vectors 5983 if(_boostedScaler != nullptr) 5984 { 5985 _boostedScaler->unscalePrimal(_boostedSolver, tmpPrimal); 5986 _boostedScaler->unscaleSlacks(_boostedSolver, tmpSlacks); 5987 _boostedScaler->unscaleDual(_boostedSolver, tmpDual); 5988 _boostedScaler->unscaleRedCost(_boostedSolver, tmpRedCost); 5989 } 5990 5991 // get basis of transformed problem 5992 _tmpBasisStatusRows.reSize(_boostedSolver.nRows()); 5993 _tmpBasisStatusCols.reSize(_boostedSolver.nCols()); 5994 5995 _boostedSolver.getBasis(_tmpBasisStatusRows.get_ptr(), _tmpBasisStatusCols.get_ptr(), 5996 _tmpBasisStatusRows.size(), 5997 _tmpBasisStatusCols.size()); 5998 _convertDataArrayVarStatusToRPrecision(_tmpBasisStatusRows, basisStatusRows); 5999 _convertDataArrayVarStatusToRPrecision(_tmpBasisStatusCols, basisStatusCols); 6000 } 6001 6002 _convertDataArrayVarStatusToBoosted(basisStatusRows, _tmpBasisStatusRows); 6003 _convertDataArrayVarStatusToBoosted(basisStatusCols, _tmpBasisStatusCols); 6004 6005 try 6006 { 6007 _boostedSimplifier->unsimplify(tmpPrimal, tmpDual, tmpSlacks, tmpRedCost, 6008 _tmpBasisStatusRows.get_ptr(), _tmpBasisStatusCols.get_ptr()); 6009 } 6010 catch(const SPxInternalCodeException& E) 6011 { 6012 // error message has already been printed 6013 throw SPxInternalCodeException("Error storing the basis"); 6014 } 6015 6016 // store basis for original problem 6017 _tmpBasisStatusRows.reSize(numRowsRational()); 6018 _tmpBasisStatusCols.reSize(numColsRational()); 6019 6020 _boostedSimplifier->getBasis(_tmpBasisStatusRows.get_ptr(), _tmpBasisStatusCols.get_ptr(), 6021 _tmpBasisStatusRows.size(), 6022 _tmpBasisStatusCols.size()); 6023 _convertDataArrayVarStatusToRPrecision(_tmpBasisStatusRows, basisStatusRows); 6024 _convertDataArrayVarStatusToRPrecision(_tmpBasisStatusCols, basisStatusCols); 6025 _hasBasis = true; 6026 6027 primal = _boostedSimplifier->unsimplifiedPrimal(); 6028 dual = _boostedSimplifier->unsimplifiedDual(); 6029 } 6030 else 6031 { 6032 _boostedSolver.getPrimalSol(primal); 6033 _boostedSolver.getDualSol(dual); 6034 6035 // unscale vectors 6036 if(_boostedScaler != nullptr) 6037 { 6038 _boostedScaler->unscalePrimal(_boostedSolver, primal); 6039 _boostedScaler->unscaleDual(_boostedSolver, dual); 6040 } 6041 6042 // get basis of transformed problem 6043 _tmpBasisStatusRows.reSize(_boostedSolver.nRows()); 6044 _tmpBasisStatusCols.reSize(_boostedSolver.nCols()); 6045 6046 _boostedSolver.getBasis(_tmpBasisStatusRows.get_ptr(), _tmpBasisStatusCols.get_ptr(), 6047 _tmpBasisStatusRows.size(), 6048 _tmpBasisStatusCols.size()); 6049 _convertDataArrayVarStatusToRPrecision(_tmpBasisStatusRows, basisStatusRows); 6050 _convertDataArrayVarStatusToRPrecision(_tmpBasisStatusCols, basisStatusCols); 6051 _hasBasis = true; 6052 } 6053 6054 if(oldIterations < _statistics->iterations) // store optimal basis as old basis 6055 _storeBasisAsOldBasisBoosted(_tmpBasisStatusRows, _tmpBasisStatusCols); 6056 6057 break; 6058 6059 case SPxSolverBase<BP>::ABORT_CYCLING: 6060 if(_boostedSimplifier == 0 && boolParam(SoPlexBase<R>::ACCEPTCYCLING)) 6061 { 6062 _boostedSolver.getPrimalSol(primal); 6063 _boostedSolver.getDualSol(dual); 6064 6065 // unscale vectors 6066 if(_boostedScaler != nullptr) 6067 { 6068 _boostedScaler->unscalePrimal(_boostedSolver, primal); 6069 _boostedScaler->unscaleDual(_boostedSolver, dual); 6070 } 6071 } 6072 6073 // get the last stable basis. The hope is that precision boosting will get rid of cycling. 6074 if(boolParam(SoPlexBase<R>::STORE_BASIS_BEFORE_SIMPLEX_PIVOT)) 6075 { 6076 try 6077 { 6078 if(oldIterations < _statistics->iterations) 6079 _storeLastStableBasisBoosted(simplificationStatus == SPxSimplifier<BP>::VANISHED); 6080 } 6081 catch(const SPxInternalCodeException& E) 6082 { 6083 SPX_MSG_INFO1(spxout, spxout << "Caught exception <" << E.what() << 6084 "> while processing the result of the solve.\n"); 6085 SPX_MSG_INFO1(spxout, spxout << "Storage of last basis failed. Keep going.\n"); 6086 } 6087 } 6088 6089 break; 6090 6091 // intentional fallthrough 6092 case SPxSolverBase<BP>::ABORT_TIME: 6093 case SPxSolverBase<BP>::ABORT_ITER: 6094 case SPxSolverBase<BP>::ABORT_VALUE: 6095 case SPxSolverBase<BP>::REGULAR: 6096 case SPxSolverBase<BP>::RUNNING: 6097 case SPxSolverBase<BP>::UNBOUNDED: 6098 _hasBasis = (_boostedSolver.basis().status() > SPxBasisBase<BP>::NO_PROBLEM); 6099 6100 if(_hasBasis && _simplifier == 0) 6101 { 6102 _tmpBasisStatusRows.reSize(_boostedSolver.nRows()); 6103 _tmpBasisStatusCols.reSize(_boostedSolver.nCols()); 6104 _boostedSolver.getBasis(_tmpBasisStatusRows.get_ptr(), _tmpBasisStatusCols.get_ptr(), 6105 _tmpBasisStatusRows.size(), 6106 _tmpBasisStatusCols.size()); 6107 _convertDataArrayVarStatusToRPrecision(_tmpBasisStatusRows, basisStatusRows); 6108 _convertDataArrayVarStatusToRPrecision(_tmpBasisStatusCols, basisStatusCols); 6109 6110 if(oldIterations < _statistics->iterations) 6111 _storeBasisAsOldBasisBoosted(_tmpBasisStatusRows, _tmpBasisStatusCols); 6112 } 6113 else 6114 { 6115 _hasBasis = false; 6116 _rationalLUSolver.clear(); 6117 } 6118 6119 // get the last stable basis. 6120 if(boolParam(SoPlexBase<R>::STORE_BASIS_BEFORE_SIMPLEX_PIVOT)) 6121 { 6122 try 6123 { 6124 if(oldIterations < _statistics->iterations) 6125 _storeLastStableBasisBoosted(simplificationStatus == SPxSimplifier<BP>::VANISHED); 6126 } 6127 catch(const SPxInternalCodeException& E) 6128 { 6129 SPX_MSG_INFO1(spxout, spxout << "Caught exception <" << E.what() << 6130 "> while processing the result of the solve.\n"); 6131 SPX_MSG_INFO1(spxout, spxout << "Storage of last basis failed. Keep going.\n"); 6132 } 6133 } 6134 6135 break; 6136 6137 case SPxSolverBase<BP>::INFEASIBLE: 6138 6139 // if simplifier is active we can currently not return a Farkas ray or basis 6140 if(_boostedSimplifier != 0) 6141 { 6142 _hasBasis = false; 6143 _rationalLUSolver.clear(); 6144 6145 // get the last stable basis. 6146 if(boolParam(SoPlexBase<R>::STORE_BASIS_BEFORE_SIMPLEX_PIVOT)) 6147 { 6148 try 6149 { 6150 if(oldIterations < _statistics->iterations) 6151 _storeLastStableBasisBoosted(simplificationStatus == SPxSimplifier<BP>::VANISHED); 6152 } 6153 catch(const SPxInternalCodeException& E) 6154 { 6155 SPX_MSG_INFO1(spxout, spxout << "Caught exception <" << E.what() << 6156 "> while processing the result of the solve.\n"); 6157 SPX_MSG_INFO1(spxout, spxout << "Storage of last basis failed. Keep going.\n"); 6158 } 6159 } 6160 6161 break; 6162 } 6163 6164 // return Farkas ray as dual solution 6165 _boostedSolver.getDualfarkas(dual); 6166 6167 // unscale vectors 6168 if(_boostedScaler != nullptr) 6169 _boostedScaler->unscaleDual(_boostedSolver, dual); 6170 6171 // get basis of transformed problem 6172 _tmpBasisStatusRows.reSize(_boostedSolver.nRows()); 6173 _tmpBasisStatusCols.reSize(_boostedSolver.nCols()); 6174 _boostedSolver.getBasis(_tmpBasisStatusRows.get_ptr(), _tmpBasisStatusCols.get_ptr(), 6175 _tmpBasisStatusRows.size(), 6176 _tmpBasisStatusCols.size()); 6177 _convertDataArrayVarStatusToRPrecision(_tmpBasisStatusRows, basisStatusRows); 6178 _convertDataArrayVarStatusToRPrecision(_tmpBasisStatusCols, basisStatusCols); 6179 _hasBasis = true; 6180 6181 // if possible, get the last stable basis, otherwise store the infeasible basis. 6182 if(boolParam(SoPlexBase<R>::STORE_BASIS_BEFORE_SIMPLEX_PIVOT)) 6183 { 6184 try 6185 { 6186 if(oldIterations < _statistics->iterations) 6187 _storeLastStableBasisBoosted(simplificationStatus == SPxSimplifier<BP>::VANISHED); 6188 } 6189 catch(const SPxInternalCodeException& E) 6190 { 6191 SPX_MSG_INFO1(spxout, spxout << "Caught exception <" << E.what() << 6192 "> while processing the result of the solve.\n"); 6193 SPX_MSG_INFO1(spxout, spxout << "Storage of last basis failed. Keep going.\n"); 6194 } 6195 } 6196 else 6197 { 6198 if(oldIterations < _statistics->iterations) 6199 _storeBasisAsOldBasisBoosted(_tmpBasisStatusRows, _tmpBasisStatusCols); 6200 } 6201 6202 break; 6203 6204 case SPxSolverBase<BP>::INForUNBD: 6205 case SPxSolverBase<BP>::SINGULAR: 6206 default: 6207 _hasBasis = false; 6208 _rationalLUSolver.clear(); 6209 6210 // get the last stable basis. 6211 if(boolParam(SoPlexBase<R>::STORE_BASIS_BEFORE_SIMPLEX_PIVOT)) 6212 { 6213 try 6214 { 6215 if(oldIterations < _statistics->iterations) 6216 _storeLastStableBasisBoosted(simplificationStatus == SPxSimplifier<BP>::VANISHED); 6217 } 6218 catch(const SPxInternalCodeException& E) 6219 { 6220 SPX_MSG_INFO1(spxout, spxout << "Caught exception <" << E.what() << 6221 "> while processing the result of the solve.\n"); 6222 SPX_MSG_INFO1(spxout, spxout << "Storage of last basis failed. Keep going.\n"); 6223 } 6224 } 6225 6226 6227 break; 6228 } 6229 } 6230 } 6231 catch(...) 6232 { 6233 SPX_MSG_INFO1(spxout, spxout << "Exception thrown during floating-point solve.\n"); 6234 boostedResult = SPxSolverBase<BP>::ERROR; 6235 _hasBasis = false; 6236 _rationalLUSolver.clear(); 6237 6238 } 6239 6240 // restore original LP if necessary 6241 if(_boostedSimplifier != 0 || _boostedScaler != nullptr) 6242 { 6243 assert(rationalLP != 0); 6244 _boostedSolver.loadLP((SPxLPBase<BP>)(*rationalLP)); 6245 rationalLP->~SPxLPRational(); 6246 spx_free(rationalLP); 6247 6248 if(_hasBasis) 6249 { 6250 _convertDataArrayVarStatusToBoosted(basisStatusRows, _tmpBasisStatusRows); 6251 _convertDataArrayVarStatusToBoosted(basisStatusCols, _tmpBasisStatusCols); 6252 _boostedSolver.setBasis(_tmpBasisStatusRows.get_ptr(), _tmpBasisStatusCols.get_ptr()); 6253 } 6254 } 6255 6256 #else 6257 SPX_MSG_ERROR(std::cerr << "No support for precision boosting without GMP/MPFR.\n"); 6258 #endif 6259 } 6260 6261 6262 6263 /// computes rational inverse of basis matrix as defined by _rationalLUSolverBind 6264 template <class R> 6265 void SoPlexBase<R>::_computeBasisInverseRational() 6266 { 6267 assert(_rationalLUSolver.status() == SLinSolverRational::UNLOADED 6268 || _rationalLUSolver.status() == SLinSolverRational::TIME); 6269 6270 const int matrixdim = numRowsRational(); 6271 assert(_rationalLUSolverBind.size() == matrixdim); 6272 6273 Array< const SVectorRational* > matrix(matrixdim); 6274 _rationalLUSolverBind.reSize(matrixdim); 6275 6276 for(int i = 0; i < matrixdim; i++) 6277 { 6278 if(_rationalLUSolverBind[i] >= 0) 6279 { 6280 assert(_rationalLUSolverBind[i] < numColsRational()); 6281 matrix[i] = &colVectorRational(_rationalLUSolverBind[i]); 6282 } 6283 else 6284 { 6285 assert(-1 - _rationalLUSolverBind[i] >= 0); 6286 assert(-1 - _rationalLUSolverBind[i] < numRowsRational()); 6287 matrix[i] = _unitVectorRational(-1 - _rationalLUSolverBind[i]); 6288 } 6289 } 6290 6291 // load and factorize rational basis matrix 6292 if(realParam(SoPlexBase<R>::TIMELIMIT) < realParam(SoPlexBase<R>::INFTY)) 6293 _rationalLUSolver.setTimeLimit((double)realParam(SoPlexBase<R>::TIMELIMIT) - 6294 _statistics->solvingTime->time()); 6295 else 6296 _rationalLUSolver.setTimeLimit(-1.0); 6297 6298 _rationalLUSolver.load(matrix.get_ptr(), matrixdim); 6299 6300 // record statistics 6301 _statistics->luFactorizationTimeRational += _rationalLUSolver.getFactorTime(); 6302 _statistics->luFactorizationsRational += _rationalLUSolver.getFactorCount(); 6303 _rationalLUSolver.resetCounters(); 6304 6305 if(_rationalLUSolver.status() == SLinSolverRational::TIME) 6306 { 6307 SPX_MSG_INFO2(spxout, spxout << "Rational factorization hit time limit.\n"); 6308 } 6309 else if(_rationalLUSolver.status() != SLinSolverRational::OK) 6310 { 6311 SPX_MSG_INFO1(spxout, spxout << "Error performing rational LU factorization.\n"); 6312 } 6313 6314 return; 6315 } 6316 6317 6318 6319 /// factorizes rational basis matrix in column representation 6320 template <class R> 6321 void SoPlexBase<R>::_factorizeColumnRational(SolRational& sol, 6322 DataArray< typename SPxSolverBase<R>::VarStatus >& basisStatusRows, 6323 DataArray< typename SPxSolverBase<R>::VarStatus >& basisStatusCols, bool& stoppedTime, 6324 bool& stoppedIter, bool& error, bool& optimal) 6325 { 6326 // start rational solving time 6327 _statistics->rationalTime->start(); 6328 6329 stoppedTime = false; 6330 stoppedIter = false; 6331 error = false; 6332 optimal = false; 6333 6334 const bool maximizing = (intParam(SoPlexBase<R>::OBJSENSE) == SoPlexBase<R>::OBJSENSE_MAXIMIZE); 6335 const int matrixdim = numRowsRational(); 6336 bool loadMatrix = (_rationalLUSolver.status() == SLinSolverRational::UNLOADED 6337 || _rationalLUSolver.status() == SLinSolverRational::TIME); 6338 int numBasicRows; 6339 6340 assert(loadMatrix || matrixdim == _rationalLUSolver.dim()); 6341 assert(loadMatrix || matrixdim == _rationalLUSolverBind.size()); 6342 6343 if(!loadMatrix && (matrixdim != _rationalLUSolver.dim() 6344 || matrixdim != _rationalLUSolverBind.size())) 6345 { 6346 SPX_MSG_WARNING(spxout, spxout << 6347 "Warning: dimensioning error in rational matrix factorization (recovered).\n"); 6348 loadMatrix = true; 6349 } 6350 6351 _workSol._primal.reDim(matrixdim); 6352 _workSol._slacks.reDim(matrixdim); 6353 _workSol._dual.reDim(matrixdim); 6354 _workSol._redCost.reDim(numColsRational() > matrixdim ? numColsRational() : matrixdim); 6355 6356 if(loadMatrix) 6357 _rationalLUSolverBind.reSize(matrixdim); 6358 6359 VectorRational& basicPrimalRhs = _workSol._slacks; 6360 VectorRational& basicDualRhs = _workSol._redCost; 6361 VectorRational& basicPrimal = _workSol._primal; 6362 VectorRational& basicDual = _workSol._dual; 6363 6364 Rational violation; 6365 Rational primalViolation; 6366 Rational dualViolation; 6367 bool primalFeasible = false; 6368 bool dualFeasible = false; 6369 6370 assert(basisStatusCols.size() == numColsRational()); 6371 assert(basisStatusRows.size() == numRowsRational()); 6372 6373 int j = 0; 6374 6375 for(int i = 0; i < basisStatusRows.size(); i++) 6376 { 6377 if(basisStatusRows[i] == SPxSolverBase<R>::BASIC && j < matrixdim) 6378 { 6379 basicPrimalRhs[i] = 0; 6380 basicDualRhs[j] = 0; 6381 6382 if(loadMatrix) 6383 _rationalLUSolverBind[j] = -1 - i; 6384 6385 j++; 6386 } 6387 else if(basisStatusRows[i] == SPxSolverBase<R>::ON_LOWER) 6388 basicPrimalRhs[i] = lhsRational(i); 6389 else if(basisStatusRows[i] == SPxSolverBase<R>::ON_UPPER) 6390 basicPrimalRhs[i] = rhsRational(i); 6391 else if(basisStatusRows[i] == SPxSolverBase<R>::ZERO) 6392 basicPrimalRhs[i] = 0; 6393 else if(basisStatusRows[i] == SPxSolverBase<R>::FIXED) 6394 { 6395 assert(lhsRational(i) == rhsRational(i)); 6396 basicPrimalRhs[i] = lhsRational(i); 6397 } 6398 else if(basisStatusRows[i] == SPxSolverBase<R>::UNDEFINED) 6399 { 6400 SPX_MSG_INFO1(spxout, spxout << "Undefined basis status of row in rational factorization.\n"); 6401 error = true; 6402 goto TERMINATE; 6403 } 6404 else 6405 { 6406 assert(basisStatusRows[i] == SPxSolverBase<R>::BASIC); 6407 SPX_MSG_INFO1(spxout, spxout << "Too many basic rows in rational factorization.\n"); 6408 error = true; 6409 goto TERMINATE; 6410 } 6411 } 6412 6413 numBasicRows = j; 6414 6415 for(int i = 0; i < basisStatusCols.size(); i++) 6416 { 6417 if(basisStatusCols[i] == SPxSolverBase<R>::BASIC && j < matrixdim) 6418 { 6419 basicDualRhs[j] = objRational(i); 6420 6421 if(loadMatrix) 6422 _rationalLUSolverBind[j] = i; 6423 6424 j++; 6425 } 6426 else if(basisStatusCols[i] == SPxSolverBase<R>::ON_LOWER) 6427 basicPrimalRhs.multAdd(-lowerRational(i), colVectorRational(i)); 6428 else if(basisStatusCols[i] == SPxSolverBase<R>::ON_UPPER) 6429 basicPrimalRhs.multAdd(-upperRational(i), colVectorRational(i)); 6430 else if(basisStatusCols[i] == SPxSolverBase<R>::ZERO) 6431 {} 6432 else if(basisStatusCols[i] == SPxSolverBase<R>::FIXED) 6433 { 6434 assert(lowerRational(i) == upperRational(i)); 6435 basicPrimalRhs.multAdd(-lowerRational(i), colVectorRational(i)); 6436 } 6437 else if(basisStatusCols[i] == SPxSolverBase<R>::UNDEFINED) 6438 { 6439 SPX_MSG_INFO1(spxout, spxout << "Undefined basis status of column in rational factorization.\n"); 6440 error = true; 6441 goto TERMINATE; 6442 } 6443 else 6444 { 6445 assert(basisStatusCols[i] == SPxSolverBase<R>::BASIC); 6446 SPX_MSG_INFO1(spxout, spxout << "Too many basic columns in rational factorization.\n"); 6447 error = true; 6448 goto TERMINATE; 6449 } 6450 } 6451 6452 if(j != matrixdim) 6453 { 6454 SPX_MSG_INFO1(spxout, spxout << "Too few basic entries in rational factorization.\n"); 6455 error = true; 6456 goto TERMINATE; 6457 } 6458 6459 // load and factorize rational basis matrix 6460 if(loadMatrix) 6461 _computeBasisInverseRational(); 6462 6463 if(_rationalLUSolver.status() == SLinSolverRational::TIME) 6464 { 6465 stoppedTime = true; 6466 return; 6467 } 6468 else if(_rationalLUSolver.status() != SLinSolverRational::OK) 6469 { 6470 error = true; 6471 return; 6472 } 6473 6474 assert(_rationalLUSolver.status() == SLinSolverRational::OK); 6475 6476 // solve for primal solution 6477 if(realParam(SoPlexBase<R>::TIMELIMIT) < realParam(SoPlexBase<R>::INFTY)) 6478 _rationalLUSolver.setTimeLimit(Real(realParam(SoPlexBase<R>::TIMELIMIT)) - 6479 _statistics->solvingTime->time()); 6480 else 6481 _rationalLUSolver.setTimeLimit(-1.0); 6482 6483 _rationalLUSolver.solveRight(basicPrimal, basicPrimalRhs); 6484 6485 // record statistics 6486 _statistics->luSolveTimeRational += _rationalLUSolver.getSolveTime(); 6487 _rationalLUSolver.resetCounters(); 6488 6489 if(_isSolveStopped(stoppedTime, stoppedIter)) 6490 { 6491 SPX_MSG_INFO2(spxout, spxout << 6492 "Rational factorization hit time limit while solving for primal.\n"); 6493 return; 6494 } 6495 6496 // check bound violation on basic rows and columns 6497 j = 0; 6498 primalViolation = 0; 6499 primalFeasible = true; 6500 6501 for(int i = 0; i < basisStatusRows.size(); i++) 6502 { 6503 if(basisStatusRows[i] == SPxSolverBase<R>::BASIC) 6504 { 6505 assert(j < matrixdim); 6506 assert(_rationalLUSolverBind[j] == -1 - i); 6507 violation = lhsRational(i); 6508 violation += basicPrimal[j]; 6509 6510 if(violation > primalViolation) 6511 { 6512 primalFeasible = false; 6513 primalViolation = violation; 6514 } 6515 6516 violation = rhsRational(i); 6517 violation += basicPrimal[j]; 6518 violation *= -1; 6519 6520 if(violation > primalViolation) 6521 { 6522 primalFeasible = false; 6523 primalViolation = violation; 6524 } 6525 6526 j++; 6527 } 6528 } 6529 6530 for(int i = 0; i < basisStatusCols.size(); i++) 6531 { 6532 if(basisStatusCols[i] == SPxSolverBase<R>::BASIC) 6533 { 6534 assert(j < matrixdim); 6535 assert(_rationalLUSolverBind[j] == i); 6536 6537 if(basicPrimal[j] < lowerRational(i)) 6538 { 6539 violation = lowerRational(i); 6540 violation -= basicPrimal[j]; 6541 6542 if(violation > primalViolation) 6543 { 6544 primalFeasible = false; 6545 primalViolation = violation; 6546 } 6547 } 6548 6549 if(basicPrimal[j] > upperRational(i)) 6550 { 6551 violation = basicPrimal[j]; 6552 violation -= upperRational(i); 6553 6554 if(violation > primalViolation) 6555 { 6556 primalFeasible = false; 6557 primalViolation = violation; 6558 } 6559 } 6560 6561 j++; 6562 } 6563 } 6564 6565 if(!primalFeasible) 6566 { 6567 SPX_MSG_INFO1(spxout, spxout << "Rational solution primal infeasible.\n"); 6568 } 6569 6570 // solve for dual solution 6571 if(realParam(SoPlexBase<R>::TIMELIMIT) < realParam(SoPlexBase<R>::INFTY)) 6572 _rationalLUSolver.setTimeLimit(Real(realParam(SoPlexBase<R>::TIMELIMIT)) - 6573 _statistics->solvingTime->time()); 6574 else 6575 _rationalLUSolver.setTimeLimit(-1.0); 6576 6577 _rationalLUSolver.solveLeft(basicDual, basicDualRhs); 6578 6579 // record statistics 6580 _statistics->luSolveTimeRational += _rationalLUSolver.getSolveTime(); 6581 _rationalLUSolver.resetCounters(); 6582 6583 if(_isSolveStopped(stoppedTime, stoppedIter)) 6584 { 6585 SPX_MSG_INFO2(spxout, spxout << "Rational factorization hit time limit while solving for dual.\n"); 6586 return; 6587 } 6588 6589 // check dual violation on nonbasic rows 6590 dualViolation = 0; 6591 dualFeasible = true; 6592 6593 for(int i = 0; i < basisStatusRows.size(); i++) 6594 { 6595 if(_rowTypes[i] == RANGETYPE_FIXED 6596 && (basisStatusRows[i] == SPxSolverBase<R>::ON_LOWER 6597 || basisStatusRows[i] == SPxSolverBase<R>::ON_UPPER)) 6598 { 6599 assert(lhsRational(i) == rhsRational(i)); 6600 basisStatusRows[i] = SPxSolverBase<R>::FIXED; 6601 } 6602 6603 assert(basisStatusRows[i] != SPxSolverBase<R>::BASIC || basicDual[i] == 0); 6604 6605 if(basisStatusRows[i] == SPxSolverBase<R>::BASIC || basisStatusRows[i] == SPxSolverBase<R>::FIXED) 6606 continue; 6607 else if(basicDual[i] < 0) 6608 { 6609 if(((maximizing && basisStatusRows[i] != SPxSolverBase<R>::ON_LOWER) || (!maximizing 6610 && basisStatusRows[i] != SPxSolverBase<R>::ON_UPPER)) 6611 && (basisStatusRows[i] != SPxSolverBase<R>::ZERO || rhsRational(i) != 0)) 6612 { 6613 dualFeasible = false; 6614 violation = -basicDual[i]; 6615 6616 if(violation > dualViolation) 6617 dualViolation = violation; 6618 6619 SPxOut::debug(this, 6620 "negative dual multliplier for row {} with dual {} and status {} and [lhs,rhs] = [{},{}]\n", 6621 i, basicDual[i].str(), basisStatusRows[i], lhsRational(i).str(), rhsRational(i).str()); 6622 } 6623 } 6624 else if(basicDual[i] > 0) 6625 { 6626 if(((maximizing && basisStatusRows[i] != SPxSolverBase<R>::ON_UPPER) || (!maximizing 6627 && basisStatusRows[i] != SPxSolverBase<R>::ON_LOWER)) 6628 && (basisStatusRows[i] != SPxSolverBase<R>::ZERO || lhsRational(i) == 0)) 6629 { 6630 dualFeasible = false; 6631 6632 if(basicDual[i] > dualViolation) 6633 dualViolation = basicDual[i]; 6634 6635 SPxOut::debug(this, 6636 "positive dual multliplier for row {} with dual {} and status {} and [lhs,rhs] = [{},{}]\n", 6637 i, basicDual[i].str(), basisStatusRows[i], lhsRational(i).str(), rhsRational(i).str()); 6638 } 6639 } 6640 } 6641 6642 // check reduced cost violation on nonbasic columns 6643 for(int i = 0; i < basisStatusCols.size(); i++) 6644 { 6645 if(_colTypes[i] == RANGETYPE_FIXED 6646 && (basisStatusCols[i] == SPxSolverBase<R>::ON_LOWER 6647 || basisStatusCols[i] == SPxSolverBase<R>::ON_UPPER)) 6648 { 6649 assert(lowerRational(i) == upperRational(i)); 6650 basisStatusCols[i] = SPxSolverBase<R>::FIXED; 6651 } 6652 6653 assert(basisStatusCols[i] != SPxSolverBase<R>::BASIC 6654 || basicDual * colVectorRational(i) == objRational(i)); 6655 6656 if(basisStatusCols[i] == SPxSolverBase<R>::BASIC || basisStatusCols[i] == SPxSolverBase<R>::FIXED) 6657 continue; 6658 else 6659 { 6660 _workSol._redCost[i] = basicDual * colVectorRational(i); 6661 _workSol._redCost[i] -= objRational(i); 6662 6663 if(_workSol._redCost[i] > 0) 6664 { 6665 if(((maximizing && basisStatusCols[i] != SPxSolverBase<R>::ON_LOWER) || (!maximizing 6666 && basisStatusCols[i] != SPxSolverBase<R>::ON_UPPER)) 6667 && (basisStatusCols[i] != SPxSolverBase<R>::ZERO || upperRational(i) != 0)) 6668 { 6669 dualFeasible = false; 6670 6671 if(_workSol._redCost[i] > dualViolation) 6672 dualViolation = _workSol._redCost[i]; 6673 } 6674 6675 _workSol._redCost[i] *= -1; 6676 } 6677 else if(_workSol._redCost[i] < 0) 6678 { 6679 _workSol._redCost[i] *= -1; 6680 6681 if(((maximizing && basisStatusCols[i] != SPxSolverBase<R>::ON_UPPER) || (!maximizing 6682 && basisStatusCols[i] != SPxSolverBase<R>::ON_LOWER)) 6683 && (basisStatusCols[i] != SPxSolverBase<R>::ZERO || lowerRational(i) != 0)) 6684 { 6685 dualFeasible = false; 6686 6687 if(_workSol._redCost[i] > dualViolation) 6688 dualViolation = _workSol._redCost[i]; 6689 } 6690 } 6691 else 6692 _workSol._redCost[i] *= -1; 6693 } 6694 } 6695 6696 if(!dualFeasible) 6697 { 6698 SPX_MSG_INFO1(spxout, spxout << "Rational solution dual infeasible.\n"); 6699 } 6700 6701 // store solution 6702 optimal = primalFeasible && dualFeasible; 6703 6704 if(optimal || boolParam(SoPlexBase<R>::RATFACJUMP)) 6705 { 6706 _hasBasis = true; 6707 6708 if(&basisStatusRows != &_basisStatusRows) 6709 _basisStatusRows = basisStatusRows; 6710 6711 if(&basisStatusCols != &_basisStatusCols) 6712 _basisStatusCols = basisStatusCols; 6713 6714 sol._primal.reDim(numColsRational()); 6715 j = numBasicRows; 6716 6717 for(int i = 0; i < basisStatusCols.size(); i++) 6718 { 6719 if(basisStatusCols[i] == SPxSolverBase<R>::BASIC) 6720 { 6721 assert(j < matrixdim); 6722 assert(_rationalLUSolverBind[j] == i); 6723 sol._primal[i] = basicPrimal[j]; 6724 j++; 6725 } 6726 else if(basisStatusCols[i] == SPxSolverBase<R>::ON_LOWER) 6727 sol._primal[i] = lowerRational(i); 6728 else if(basisStatusCols[i] == SPxSolverBase<R>::ON_UPPER) 6729 sol._primal[i] = upperRational(i); 6730 else if(basisStatusCols[i] == SPxSolverBase<R>::ZERO) 6731 sol._primal[i] = 0; 6732 else if(basisStatusCols[i] == SPxSolverBase<R>::FIXED) 6733 { 6734 assert(lowerRational(i) == upperRational(i)); 6735 sol._primal[i] = lowerRational(i); 6736 } 6737 else 6738 { 6739 assert(basisStatusCols[i] == SPxSolverBase<R>::UNDEFINED); 6740 SPX_MSG_INFO1(spxout, spxout << "Undefined basis status of column in rational factorization.\n"); 6741 error = true; 6742 goto TERMINATE; 6743 } 6744 } 6745 6746 sol._slacks.reDim(numRowsRational()); 6747 _rationalLP->computePrimalActivity(sol._primal, sol._slacks); 6748 sol._isPrimalFeasible = true; 6749 6750 sol._dual = basicDual; 6751 6752 for(int i = 0; i < numColsRational(); i++) 6753 { 6754 if(basisStatusCols[i] == SPxSolverBase<R>::BASIC) 6755 sol._redCost[i] = 0; 6756 else if(basisStatusCols[i] == SPxSolverBase<R>::FIXED) 6757 { 6758 sol._redCost[i] = basicDual * colVectorRational(i); 6759 sol._redCost[i] -= objRational(i); 6760 sol._redCost[i] *= -1; 6761 } 6762 else 6763 sol._redCost[i] = _workSol._redCost[i]; 6764 } 6765 6766 sol._isDualFeasible = true; 6767 } 6768 else 6769 { 6770 _rationalLUSolver.clear(); 6771 } 6772 6773 6774 TERMINATE: 6775 // stop rational solving time 6776 _statistics->rationalTime->stop(); 6777 return; 6778 } 6779 6780 /// attempts rational reconstruction of primal-dual solution 6781 template <class R> 6782 bool SoPlexBase<R>::_reconstructSolutionRational(SolRational& sol, 6783 DataArray< typename SPxSolverBase<R>::VarStatus >& basisStatusRows, 6784 DataArray< typename SPxSolverBase<R>::VarStatus >& basisStatusCols, 6785 const Rational& denomBoundSquared) 6786 { 6787 bool success; 6788 bool isSolBasic; 6789 DIdxSet basicIndices(numColsRational()); 6790 6791 success = false; 6792 isSolBasic = true; 6793 6794 if(!sol.isPrimalFeasible() || !sol.isDualFeasible()) 6795 return success; 6796 6797 // start timing and increment statistics counter 6798 _statistics->reconstructionTime->start(); 6799 _statistics->rationalReconstructions++; 6800 6801 // reconstruct primal vector 6802 _workSol._primal = sol._primal; 6803 6804 for(int j = 0; j < numColsRational(); ++j) 6805 { 6806 if(basisStatusCols[j] == SPxSolverBase<R>::BASIC) 6807 basicIndices.addIdx(j); 6808 } 6809 6810 success = reconstructVector(_workSol._primal, denomBoundSquared, &basicIndices); 6811 6812 if(!success) 6813 { 6814 SPX_MSG_INFO1(spxout, spxout << "Rational reconstruction of primal solution failed.\n"); 6815 _statistics->reconstructionTime->stop(); 6816 return success; 6817 } 6818 6819 SPxOut::debug(this, "Rational reconstruction of primal solution successful.\n"); 6820 6821 // check violation of bounds 6822 for(int c = numColsRational() - 1; c >= 0; c--) 6823 { 6824 // we want to notify the user whether the reconstructed solution is basic; otherwise, this would be redundant 6825 typename SPxSolverBase<R>::VarStatus& basisStatusCol = _basisStatusCols[c]; 6826 6827 if((basisStatusCol == SPxSolverBase<R>::FIXED && _workSol._primal[c] != lowerRational(c)) 6828 || (basisStatusCol == SPxSolverBase<R>::ON_LOWER && _workSol._primal[c] != lowerRational(c)) 6829 || (basisStatusCol == SPxSolverBase<R>::ON_UPPER && _workSol._primal[c] != upperRational(c)) 6830 || (basisStatusCol == SPxSolverBase<R>::ZERO && _workSol._primal[c] != 0) 6831 || (basisStatusCol == SPxSolverBase<R>::UNDEFINED)) 6832 { 6833 isSolBasic = false; 6834 } 6835 6836 if(_lowerFinite(_colTypes[c]) && _workSol._primal[c] < lowerRational(c)) 6837 { 6838 SPxOut::debug(this, "Lower bound of variable {} violated by {}\n", 6839 c, (static_cast<Rational>(lowerRational(c) - _workSol._primal[c])).str()); 6840 SPX_MSG_INFO1(spxout, spxout << "Reconstructed solution primal infeasible (1).\n"); 6841 _statistics->reconstructionTime->stop(); 6842 return false; 6843 } 6844 6845 if(_upperFinite(_colTypes[c]) && _workSol._primal[c] > upperRational(c)) 6846 { 6847 SPxOut::debug(this, "Upper bound of variable {} violated by {}\n", 6848 c, (static_cast<Rational>(_workSol._primal[c] - upperRational(c))).str()); 6849 SPX_MSG_INFO1(spxout, spxout << "Reconstructed solution primal infeasible (2).\n"); 6850 _statistics->reconstructionTime->stop(); 6851 return false; 6852 } 6853 } 6854 6855 // compute slacks 6856 ///@todo we should compute them one by one so we can abort when encountering an infeasibility 6857 _workSol._slacks.reDim(numRowsRational()); 6858 _rationalLP->computePrimalActivity(_workSol._primal, _workSol._slacks); 6859 6860 // check violation of sides 6861 for(int r = numRowsRational() - 1; r >= 0; r--) 6862 { 6863 // we want to notify the user whether the reconstructed solution is basic; otherwise, this would be redundant 6864 typename SPxSolverBase<R>::VarStatus& basisStatusRow = _basisStatusRows[r]; 6865 6866 if((basisStatusRow == SPxSolverBase<R>::FIXED && _workSol._slacks[r] != lhsRational(r)) 6867 || (basisStatusRow == SPxSolverBase<R>::ON_LOWER && _workSol._slacks[r] != lhsRational(r)) 6868 || (basisStatusRow == SPxSolverBase<R>::ON_UPPER && _workSol._slacks[r] != rhsRational(r)) 6869 || (basisStatusRow == SPxSolverBase<R>::ZERO && _workSol._slacks[r] != 0) 6870 || (basisStatusRow == SPxSolverBase<R>::UNDEFINED)) 6871 { 6872 isSolBasic = false; 6873 } 6874 6875 if(_lowerFinite(_rowTypes[r]) && _workSol._slacks[r] < lhsRational(r)) 6876 { 6877 SPxOut::debug(this, "Lhs of row {} violated by {}\n", 6878 r, (static_cast<Rational>(lhsRational(r) - _workSol._slacks[r])).str()); 6879 SPX_MSG_INFO1(spxout, spxout << "Reconstructed solution primal infeasible (3).\n"); 6880 _statistics->reconstructionTime->stop(); 6881 return false; 6882 } 6883 6884 if(_upperFinite(_rowTypes[r]) && _workSol._slacks[r] > rhsRational(r)) 6885 { 6886 SPxOut::debug(this, "Rhs of row {} violated by {}\n", 6887 r, (static_cast<Rational>(_workSol._slacks[r] - rhsRational(r))).str()); 6888 SPX_MSG_INFO1(spxout, spxout << "Reconstructed solution primal infeasible (4).\n"); 6889 _statistics->reconstructionTime->stop(); 6890 return false; 6891 } 6892 } 6893 6894 // reconstruct dual vector 6895 _workSol._dual = sol._dual; 6896 6897 success = reconstructVector(_workSol._dual, denomBoundSquared); 6898 6899 if(!success) 6900 { 6901 SPX_MSG_INFO1(spxout, spxout << "Rational reconstruction of dual solution failed.\n"); 6902 _statistics->reconstructionTime->stop(); 6903 return success; 6904 } 6905 6906 SPxOut::debug(this, "Rational reconstruction of dual vector successful.\n"); 6907 6908 // check dual multipliers before reduced costs because this check is faster since it does not require the 6909 // computation of reduced costs 6910 const bool maximizing = (intParam(SoPlexBase<R>::OBJSENSE) == SoPlexBase<R>::OBJSENSE_MAXIMIZE); 6911 6912 for(int r = numRowsRational() - 1; r >= 0; r--) 6913 { 6914 int sig = sign(_workSol._dual[r]); 6915 6916 if((!maximizing && sig > 0) || (maximizing && sig < 0)) 6917 { 6918 if(!_lowerFinite(_rowTypes[r]) || _workSol._slacks[r] > lhsRational(r)) 6919 { 6920 SPxOut::debug(this, 6921 "complementary slackness violated by row {} with dual {} and slack {} not at lhs {}\n", 6922 r, _workSol._dual[r].str(), _workSol._slacks[r].str(), lhsRational(r).str()); 6923 SPX_MSG_INFO1(spxout, spxout << "Reconstructed solution dual infeasible (1).\n"); 6924 _statistics->reconstructionTime->stop(); 6925 return false; 6926 } 6927 6928 if(_basisStatusRows[r] != SPxSolverBase<R>::ON_LOWER 6929 && _basisStatusRows[r] != SPxSolverBase<R>::FIXED) 6930 { 6931 if(_basisStatusRows[r] == SPxSolverBase<R>::BASIC 6932 || _basisStatusRows[r] == SPxSolverBase<R>::UNDEFINED) 6933 isSolBasic = false; 6934 else 6935 _basisStatusRows[r] = SPxSolverBase<R>::ON_LOWER; 6936 } 6937 } 6938 else if((!maximizing && sig < 0) || (maximizing && sig > 0)) 6939 { 6940 if(!_upperFinite(_rowTypes[r]) || _workSol._slacks[r] < rhsRational(r)) 6941 { 6942 SPxOut::debug(this, 6943 "complementary slackness violated by row {} with dual {} and slack {} not at rhs {}\n", 6944 r, _workSol._dual[r].str(), _workSol._slacks[r].str(), rhsRational(r).str()); 6945 SPX_MSG_INFO1(spxout, spxout << "Reconstructed solution dual infeasible (2).\n"); 6946 _statistics->reconstructionTime->stop(); 6947 return false; 6948 } 6949 6950 if(_basisStatusRows[r] != SPxSolverBase<R>::ON_UPPER 6951 && _basisStatusRows[r] != SPxSolverBase<R>::FIXED) 6952 { 6953 if(_basisStatusRows[r] == SPxSolverBase<R>::BASIC 6954 || _basisStatusRows[r] == SPxSolverBase<R>::UNDEFINED) 6955 isSolBasic = false; 6956 else 6957 _basisStatusRows[r] = SPxSolverBase<R>::ON_UPPER; 6958 } 6959 } 6960 } 6961 6962 // compute reduced cost vector; we assume that the objective function vector has less nonzeros than the reduced 6963 // cost vector, and so multiplying with -1 first and subtracting the dual activity should be faster than adding 6964 // the dual activity and negating afterwards 6965 ///@todo we should compute them one by one so we can abort when encountering an infeasibility 6966 _workSol._redCost.reDim(numColsRational()); 6967 _rationalLP->getObj(_workSol._redCost); 6968 _rationalLP->subDualActivity(_workSol._dual, _workSol._redCost); 6969 6970 // check reduced cost violation 6971 for(int c = numColsRational() - 1; c >= 0; c--) 6972 { 6973 int sig = sign(_workSol._redCost[c]); 6974 6975 if((!maximizing && sig > 0) || (maximizing && sig < 0)) 6976 { 6977 if(!_lowerFinite(_colTypes[c]) || _workSol._primal[c] > lowerRational(c)) 6978 { 6979 SPxOut::debug(this, 6980 "complementary slackness violated by column {} with reduced cost {} and value {} not at lower bound {}\n", 6981 c, _workSol._redCost[c].str(), _workSol._primal[c].str(), lowerRational(c).str()); 6982 SPX_MSG_INFO1(spxout, spxout << "Reconstructed solution dual infeasible (3).\n"); 6983 _statistics->reconstructionTime->stop(); 6984 return false; 6985 } 6986 6987 if(_basisStatusCols[c] != SPxSolverBase<R>::ON_LOWER 6988 && _basisStatusCols[c] != SPxSolverBase<R>::FIXED) 6989 { 6990 if(_basisStatusCols[c] == SPxSolverBase<R>::BASIC 6991 || _basisStatusCols[c] == SPxSolverBase<R>::UNDEFINED) 6992 isSolBasic = false; 6993 else 6994 _basisStatusCols[c] = SPxSolverBase<R>::ON_LOWER; 6995 } 6996 } 6997 else if((!maximizing && sig < 0) || (maximizing && sig > 0)) 6998 { 6999 if(!_upperFinite(_colTypes[c]) || _workSol._primal[c] < upperRational(c)) 7000 { 7001 SPxOut::debug(this, 7002 "complementary slackness violated by column {} with reduced cost {} and value {} not at upper bound {}\n", 7003 c, _workSol._redCost[c].str(), _workSol._primal[c].str(), upperRational(c).str()); 7004 SPX_MSG_INFO1(spxout, spxout << "Reconstructed solution dual infeasible (4).\n"); 7005 _statistics->reconstructionTime->stop(); 7006 return false; 7007 } 7008 7009 if(_basisStatusCols[c] != SPxSolverBase<R>::ON_UPPER 7010 && _basisStatusCols[c] != SPxSolverBase<R>::FIXED) 7011 { 7012 if(_basisStatusCols[c] == SPxSolverBase<R>::BASIC 7013 || _basisStatusCols[c] == SPxSolverBase<R>::UNDEFINED) 7014 isSolBasic = false; 7015 else 7016 _basisStatusCols[c] = SPxSolverBase<R>::ON_UPPER; 7017 } 7018 } 7019 } 7020 7021 // update solution 7022 sol._primal = _workSol._primal; 7023 sol._slacks = _workSol._slacks; 7024 sol._dual = _workSol._dual; 7025 sol._redCost = _workSol._redCost; 7026 7027 if(!isSolBasic) 7028 { 7029 SPX_MSG_WARNING(spxout, spxout << "Warning: Reconstructed solution not basic.\n"); 7030 _hasBasis = false; 7031 } 7032 7033 // stop timing 7034 _statistics->reconstructionTime->stop(); 7035 7036 return success; 7037 } 7038 } // namespace soplex 7039