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 31 #define SOPLEX_ALLOWED_UNSCALE_PERCENTAGE 0.1 32 #define SOPLEX_MIN_OPT_CALLS_WITH_SCALING 10 33 34 namespace soplex 35 { 36 /// solves real LP 37 template <class R> 38 void SoPlexBase<R>::_optimize(volatile bool* interrupt) 39 { 40 assert(_realLP != 0); 41 assert(_realLP == &_solver); 42 43 _solReal.invalidate(); 44 ++_optimizeCalls; 45 46 // start timing 47 _statistics->solvingTime->start(); 48 49 if(boolParam(SoPlexBase<R>::PERSISTENTSCALING)) 50 { 51 // scale original problem; overwriting _realLP 52 if(_scaler && !_realLP->isScaled() && _reapplyPersistentScaling()) 53 { 54 #ifdef SOPLEX_DEBUG 55 SPxLPBase<R>* origLP = 0; 56 spx_alloc(origLP); 57 origLP = new(origLP) SPxLPBase<R>(*_realLP); 58 #endif 59 _scaler->scale(*_realLP, true); 60 _isRealLPScaled = _realLP->isScaled(); // a scaler might decide not to apply scaling 61 _solver.invalidateBasis(); 62 #ifdef SOPLEX_DEBUG 63 _checkScaling(origLP); 64 #endif 65 } 66 // unscale previously scaled problem, overwriting _realLP 67 else if(!_scaler && _realLP->isScaled()) 68 { 69 _solver.unscaleLPandReloadBasis(); 70 _isRealLPScaled = false; 71 ++_unscaleCalls; 72 } 73 } 74 75 // remember that last solve was in floating-point 76 _lastSolveMode = SOLVEMODE_REAL; 77 78 // solve and store solution; if we have a starting basis, do not apply preprocessing; if we are solving from 79 // scratch, apply preprocessing according to parameter settings 80 if(!_hasBasis && realParam(SoPlexBase<R>::OBJLIMIT_LOWER) == -realParam(SoPlexBase<R>::INFTY) 81 && realParam(SoPlexBase<R>::OBJLIMIT_UPPER) == realParam(SoPlexBase<R>::INFTY)) 82 _preprocessAndSolveReal(true, interrupt); 83 else 84 _preprocessAndSolveReal(false, interrupt); 85 86 _statistics->finalBasisCondition = _solver.getBasisMetric(0); 87 88 // stop timing 89 _statistics->solvingTime->stop(); 90 } 91 92 93 94 /// check whether persistent scaling is supposed to be reapplied again after unscaling 95 template <class R> 96 bool SoPlexBase<R>::_reapplyPersistentScaling() const 97 { 98 if((_unscaleCalls > _optimizeCalls * SOPLEX_ALLOWED_UNSCALE_PERCENTAGE) 99 && _optimizeCalls > SOPLEX_MIN_OPT_CALLS_WITH_SCALING) 100 return false; 101 else 102 return true; 103 } 104 105 106 107 /// checks result of the solving process and solves again without preprocessing if necessary 108 template <class R> 109 void SoPlexBase<R>::_evaluateSolutionReal(typename SPxSimplifier<R>::Result simplificationStatus) 110 { 111 // if the simplifier detected infeasibility or unboundedness we optimize again 112 // just to get the proof (primal or dual ray) 113 // todo get infeasibility proof from simplifier 114 switch(simplificationStatus) 115 { 116 case SPxSimplifier<R>::INFEASIBLE: 117 case SPxSimplifier<R>::DUAL_INFEASIBLE: 118 case SPxSimplifier<R>::UNBOUNDED: 119 _hasBasis = false; 120 121 if(boolParam(SoPlexBase<R>::ENSURERAY)) 122 { 123 SPX_MSG_INFO1(spxout, spxout << 124 "simplifier detected infeasibility or unboundedness - solve again without simplifying" << std::endl; 125 ) 126 _preprocessAndSolveReal(false); 127 } 128 else 129 { 130 if(simplificationStatus == SPxSimplifier<R>::INFEASIBLE) 131 _status = SPxSolverBase<R>::INFEASIBLE; 132 else if(simplificationStatus == SPxSimplifier<R>::UNBOUNDED) 133 _status = SPxSolverBase<R>::UNBOUNDED; 134 else 135 _status = SPxSolverBase<R>::INForUNBD; 136 137 // load original LP to restore clean solver state 138 _loadRealLP(false); 139 } 140 141 return; 142 143 case SPxSimplifier<R>::VANISHED: 144 _status = SPxSolverBase<R>::OPTIMAL; 145 _storeSolutionRealFromPresol(); 146 return; 147 148 case SPxSimplifier<R>::OKAY: 149 _status = _solver.status(); 150 } 151 152 // process result 153 switch(_status) 154 { 155 case SPxSolverBase<R>::OPTIMAL: 156 _storeSolutionReal(!_isRealLPLoaded || _isRealLPScaled); 157 158 // apply polishing on original problem 159 if(_applyPolishing) 160 { 161 int polishing = intParam(SoPlexBase<R>::SOLUTION_POLISHING); 162 setIntParam(SoPlexBase<R>::SOLUTION_POLISHING, polishing); 163 _preprocessAndSolveReal(false); 164 } 165 166 break; 167 168 case SPxSolverBase<R>::UNBOUNDED: 169 case SPxSolverBase<R>::INFEASIBLE: 170 case SPxSolverBase<R>::INForUNBD: 171 172 // in case of infeasibility or unboundedness, we currently can not unsimplify, but have to solve the original LP again 173 if(!_isRealLPLoaded && boolParam(SoPlexBase<R>::ENSURERAY)) 174 { 175 SPX_MSG_INFO1(spxout, spxout << " --- loading original problem" << std::endl;) 176 _solver.changeObjOffset(realParam(SoPlexBase<R>::OBJ_OFFSET)); 177 // we cannot do more to remove violations 178 _resolveWithoutPreprocessing(simplificationStatus); 179 } 180 else 181 { 182 _storeSolutionReal(false); 183 } 184 185 break; 186 187 case SPxSolverBase<R>::SINGULAR: 188 189 // if preprocessing was applied, try to run again without to avoid singularity 190 if(!_isRealLPLoaded) 191 { 192 SPX_MSG_INFO1(spxout, spxout << 193 "encountered singularity - trying to solve again without simplifying" << 194 std::endl;) 195 _preprocessAndSolveReal(false); 196 return; 197 } 198 199 _hasBasis = false; 200 break; 201 202 case SPxSolverBase<R>::ABORT_VALUE: 203 204 // If we aborted the solve for some reason and there is still a shift, ensure that the basis status is correct 205 if(_solver.shift() > _solver.epsilon()) 206 _solver.setBasisStatus(SPxBasisBase<R>::REGULAR); 207 208 _storeSolutionReal(true); 209 break; 210 211 case SPxSolverBase<R>::ABORT_CYCLING: 212 213 // if preprocessing or scaling was applied, try to run again without to 214 // avoid cycling 215 if(!_isRealLPLoaded || _isRealLPScaled) 216 { 217 SPX_MSG_INFO1(spxout, spxout << "encountered cycling - trying to solve again without simplifying" << 218 std::endl;) 219 // store and unsimplify sub-optimal solution and basis, may trigger re-solve 220 _storeSolutionReal(true); 221 return; 222 } 223 224 if(_solReal.isPrimalFeasible() || _solReal.isDualFeasible()) 225 _status = SPxSolverBase<R>::OPTIMAL_UNSCALED_VIOLATIONS; 226 227 // FALLTHROUGH 228 case SPxSolverBase<R>::ABORT_TIME: 229 case SPxSolverBase<R>::ABORT_ITER: 230 case SPxSolverBase<R>::REGULAR: 231 case SPxSolverBase<R>::RUNNING: 232 233 // If we aborted the solve for some reason and there is still a shift, ensure that the basis status is correct 234 if(_solver.shift() > _solver.epsilon()) 235 _solver.setBasisStatus(SPxBasisBase<R>::REGULAR); 236 237 _storeSolutionReal(false); 238 break; 239 240 default: 241 _hasBasis = false; 242 break; 243 } 244 } 245 246 247 248 /// solves real LP with/without preprocessing 249 template <class R> 250 void SoPlexBase<R>::_preprocessAndSolveReal(bool applySimplifier, volatile bool* interrupt) 251 { 252 _solver.changeObjOffset(realParam(SoPlexBase<R>::OBJ_OFFSET)); 253 _statistics->preprocessingTime->start(); 254 255 _applyPolishing = false; 256 257 if(applySimplifier) 258 _enableSimplifierAndScaler(); 259 else 260 _disableSimplifierAndScaler(); 261 262 // create a copy of the LP when simplifying or when using internal scaling, i.e. w/o persistent scaling 263 bool copyLP = (_simplifier != 0 || (_scaler && !_isRealLPScaled)); 264 265 _solver.setTerminationValue(intParam(SoPlexBase<R>::OBJSENSE) == SoPlexBase<R>::OBJSENSE_MINIMIZE 266 ? realParam(SoPlexBase<R>::OBJLIMIT_UPPER) : realParam(SoPlexBase<R>::OBJLIMIT_LOWER)); 267 268 if(_isRealLPLoaded) 269 { 270 assert(_realLP == &_solver); 271 272 // preprocessing is always applied to the LP in the solver; hence we have to create a copy of the original LP 273 // if simplifier is turned on 274 if(copyLP) 275 { 276 _realLP = nullptr; 277 spx_alloc(_realLP); 278 _realLP = new(_realLP) SPxLPBase<R>(_solver); 279 _isRealLPLoaded = false; 280 } 281 } 282 else 283 { 284 assert(_realLP != &_solver); 285 286 // load real LP and basis if available 287 if(_hasBasis) 288 { 289 assert(_basisStatusRows.size() == numRows()); 290 assert(_basisStatusCols.size() == this->numCols()); 291 292 _solver.loadLP(*_realLP, false); 293 _solver.setBasis(_basisStatusRows.get_const_ptr(), _basisStatusCols.get_const_ptr()); 294 } 295 // load real LP and set up slack basis 296 else 297 _solver.loadLP(*_realLP, true); 298 299 // if there is no simplifier, then the original and the transformed problem are identical and it is more 300 // memory-efficient to keep only the problem in the solver 301 if(!copyLP) 302 { 303 _realLP->~SPxLPBase<R>(); 304 spx_free(_realLP); 305 _realLP = &_solver; 306 _isRealLPLoaded = true; 307 } 308 } 309 310 // assert that we have two problems if and only if we apply the simplifier 311 assert(_realLP == &_solver || copyLP); 312 assert(_realLP != &_solver || !copyLP); 313 314 // apply problem simplification 315 typename SPxSimplifier<R>::Result simplificationStatus = SPxSimplifier<R>::OKAY; 316 317 if(_simplifier) 318 { 319 assert(!_isRealLPLoaded); 320 // do not remove bounds of boxed variables or sides of ranged rows if bound flipping is used; also respect row-boundflip parameter 321 bool keepbounds = intParam(SoPlexBase<R>::RATIOTESTER) == SoPlexBase<R>::RATIOTESTER_BOUNDFLIPPING; 322 323 if(intParam(SoPlexBase<R>::REPRESENTATION) == SoPlexBase<R>::REPRESENTATION_ROW 324 || (intParam(SoPlexBase<R>::REPRESENTATION) == SoPlexBase<R>::REPRESENTATION_AUTO 325 && (_solver.nCols() + 1) * realParam(SoPlexBase<R>::REPRESENTATION_SWITCH) < (_solver.nRows() + 1))) 326 keepbounds &= boolParam(SoPlexBase<R>::ROWBOUNDFLIPS); 327 328 Real remainingTime = _solver.getMaxTime() - _solver.time(); 329 simplificationStatus = _simplifier->simplify(_solver, remainingTime, keepbounds, 330 _solver.random.getSeed()); 331 _solver.changeObjOffset(_simplifier->getObjoffset() + realParam(SoPlexBase<R>::OBJ_OFFSET)); 332 _solver.setScalingInfo(false); 333 _applyPolishing = true; 334 335 _solver.setSolutionPolishing(SPxSolverBase<R>::POLISH_OFF); 336 } 337 338 _statistics->preprocessingTime->stop(); 339 340 // run the simplex method if problem has not been solved by the simplifier 341 if(simplificationStatus == SPxSimplifier<R>::OKAY) 342 { 343 if(_scaler && !_solver.isScaled()) 344 { 345 _scaler->scale(_solver, false); 346 _solver.invalidateBasis(); 347 } 348 349 _solveRealLPAndRecordStatistics(interrupt); 350 } 351 352 _evaluateSolutionReal(simplificationStatus); 353 } 354 355 356 357 /// loads original problem into solver and solves again after it has been solved to infeasibility or unboundedness with preprocessing 358 template <class R> 359 void SoPlexBase<R>::_resolveWithoutPreprocessing(typename SPxSimplifier<R>::Result 360 simplificationStatus) 361 { 362 assert(!_isRealLPLoaded || _scaler != nullptr); 363 assert(_simplifier != 0 || _scaler != nullptr); 364 assert(_status == SPxSolverBase<R>::INFEASIBLE || _status == SPxSolverBase<R>::INForUNBD 365 || _status == SPxSolverBase<R>::UNBOUNDED); 366 367 // if simplifier was active, then we unsimplify to get the basis 368 if(_simplifier) 369 { 370 assert(!_simplifier->isUnsimplified()); 371 assert(simplificationStatus == SPxSimplifier<R>::OKAY); 372 373 // get temporary solution vectors for transformed problem 374 VectorBase<R> primal(_solver.nCols()); 375 VectorBase<R> slacks(_solver.nRows()); 376 VectorBase<R> dual(_solver.nRows()); 377 VectorBase<R> redCost(_solver.nCols()); 378 379 _basisStatusRows.reSize(numRows()); 380 _basisStatusCols.reSize(numCols()); 381 assert(_basisStatusRows.size() >= _solver.nRows()); 382 assert(_basisStatusCols.size() >= _solver.nCols()); 383 384 // get solution data from transformed problem 385 _solver.getPrimalSol(primal); 386 _solver.getSlacks(slacks); 387 _solver.getDualSol(dual); 388 _solver.getRedCostSol(redCost); 389 390 // unscale vectors 391 if(_scaler && _solver.isScaled()) 392 { 393 _scaler->unscalePrimal(_solver, primal); 394 _scaler->unscaleSlacks(_solver, slacks); 395 _scaler->unscaleDual(_solver, dual); 396 _scaler->unscaleRedCost(_solver, redCost); 397 } 398 399 // get basis of transformed problem 400 _solver.getBasis(_basisStatusRows.get_ptr(), _basisStatusCols.get_ptr(), _basisStatusRows.size(), 401 _basisStatusCols.size()); 402 403 try 404 { 405 _simplifier->unsimplify(primal, dual, slacks, redCost, _basisStatusRows.get_ptr(), 406 _basisStatusCols.get_ptr(), false); 407 _simplifier->getBasis(_basisStatusRows.get_ptr(), _basisStatusCols.get_ptr(), 408 _basisStatusRows.size(), _basisStatusCols.size()); 409 _hasBasis = true; 410 } 411 catch(const SPxException& E) 412 { 413 SPX_MSG_INFO1(spxout, spxout << "Caught exception <" << E.what() << 414 "> during unsimplification. Resolving without simplifier and scaler.\n"); 415 _hasBasis = false; 416 } 417 } 418 // if the original problem is not in the solver because of scaling, we also need to store the basis 419 else if(_scaler != nullptr) 420 { 421 _basisStatusRows.reSize(numRows()); 422 _basisStatusCols.reSize(numCols()); 423 assert(_basisStatusRows.size() == _solver.nRows()); 424 assert(_basisStatusCols.size() == _solver.nCols()); 425 426 _solver.getBasis(_basisStatusRows.get_ptr(), _basisStatusCols.get_ptr(), _basisStatusRows.size(), 427 _basisStatusCols.size()); 428 _hasBasis = true; 429 } 430 431 // resolve the original problem 432 _preprocessAndSolveReal(false); 433 return; 434 } 435 436 437 438 /// verify computed solution based on status and resolve if claimed primal or dual feasibility is not fulfilled 439 template <class R> 440 void SoPlexBase<R>::_verifySolutionReal() 441 { 442 assert(_hasSolReal); 443 444 SPX_MSG_INFO1(spxout, spxout << " --- verifying computed solution" << std::endl;) 445 446 R sumviol = 0; 447 R boundviol = 0; 448 R rowviol = 0; 449 R dualviol = 0; 450 R redcostviol = 0; 451 452 (void) getBoundViolation(boundviol, sumviol); 453 (void) getRowViolation(rowviol, sumviol); 454 (void) getDualViolation(dualviol, sumviol); 455 (void) getRedCostViolation(redcostviol, sumviol); 456 457 if(boundviol >= _solver.tolerances()->floatingPointFeastol() 458 || rowviol >= _solver.tolerances()->floatingPointFeastol() 459 || dualviol >= _solver.tolerances()->floatingPointOpttol() 460 || redcostviol >= _solver.tolerances()->floatingPointOpttol()) 461 { 462 assert(&_solver == _realLP); 463 assert(_isRealLPLoaded); 464 SPX_MSG_INFO3(spxout, spxout << "bound violation: " << boundviol 465 << ", row violation: " << rowviol 466 << ", dual violation: " << dualviol 467 << ", redcost violation: " << redcostviol << std::endl;) 468 SPX_MSG_INFO1(spxout, spxout << 469 " --- detected violations in original problem space -- solve again without presolving/scaling" << 470 std::endl;) 471 472 if(_isRealLPScaled) 473 { 474 _solver.unscaleLPandReloadBasis(); 475 _isRealLPScaled = false; 476 ++_unscaleCalls; 477 } 478 479 _preprocessAndSolveReal(false); 480 } 481 } 482 483 /// verify computed solution based on status and resolve if claimed primal or dual feasibility is not fulfilled 484 template <class R> 485 void SoPlexBase<R>::_verifyObjLimitReal() 486 { 487 SPX_MSG_INFO1(spxout, spxout << " --- verifying objective limit" << std::endl;) 488 489 R sumviol = 0; 490 R dualviol = 0; 491 R redcostviol = 0; 492 493 bool dualfeasible = true; 494 495 if(!getDualViolation(dualviol, sumviol)) 496 dualfeasible = false; 497 498 if(!getRedCostViolation(redcostviol, sumviol)) 499 dualfeasible = false; 500 501 if(!dualfeasible || dualviol >= _solver.tolerances()->floatingPointOpttol() 502 || redcostviol >= _solver.tolerances()->floatingPointOpttol()) 503 { 504 assert(&_solver == _realLP); 505 assert(_isRealLPLoaded); 506 SPX_MSG_INFO3(spxout, spxout << ", dual violation: " << dualviol 507 << ", redcost violation: " << redcostviol << std::endl;) 508 SPX_MSG_INFO1(spxout, spxout << 509 " --- detected violations in original problem space -- solve again without presolving/scaling" << 510 std::endl;) 511 512 if(_isRealLPScaled) 513 { 514 _solver.unscaleLPandReloadBasis(); 515 _isRealLPScaled = false; 516 ++_unscaleCalls; 517 } 518 519 _preprocessAndSolveReal(false); 520 } 521 } 522 523 524 /// stores solution data from the solver, possibly after applying unscaling and unsimplifying 525 template <class R> 526 void SoPlexBase<R>::_storeSolutionReal(bool verify) 527 { 528 // prepare storage for basis (enough to fit the original basis) 529 _basisStatusRows.reSize(numRows()); 530 _basisStatusCols.reSize(numCols()); 531 532 // prepare storage for the solution data (only in transformed space due to unscaling), w/o setting it to zero 533 _solReal._primal.reDim(_solver.nCols(), false); 534 _solReal._slacks.reDim(_solver.nRows(), false); 535 _solReal._dual.reDim(_solver.nRows(), false); 536 _solReal._redCost.reDim(_solver.nCols(), false); 537 538 // check primal status consistency and query solution status 539 assert(_solver.basis().status() != SPxBasisBase<R>::PRIMAL || status() != SPxSolverBase<R>::ERROR); 540 assert(_solver.basis().status() != SPxBasisBase<R>::PRIMAL 541 || status() != SPxSolverBase<R>::NO_RATIOTESTER); 542 assert(_solver.basis().status() != SPxBasisBase<R>::PRIMAL 543 || status() != SPxSolverBase<R>::NO_PRICER); 544 assert(_solver.basis().status() != SPxBasisBase<R>::PRIMAL 545 || status() != SPxSolverBase<R>::NO_SOLVER); 546 assert(_solver.basis().status() != SPxBasisBase<R>::PRIMAL 547 || status() != SPxSolverBase<R>::NOT_INIT); 548 assert(_solver.basis().status() != SPxBasisBase<R>::PRIMAL 549 || status() != SPxSolverBase<R>::SINGULAR); 550 assert(_solver.basis().status() != SPxBasisBase<R>::PRIMAL 551 || status() != SPxSolverBase<R>::NO_PROBLEM); 552 assert(_solver.basis().status() != SPxBasisBase<R>::PRIMAL 553 || status() != SPxSolverBase<R>::UNBOUNDED); 554 assert(_solver.basis().status() != SPxBasisBase<R>::PRIMAL 555 || status() != SPxSolverBase<R>::INFEASIBLE); 556 assert(_solver.basis().status() != SPxBasisBase<R>::UNBOUNDED 557 || status() == SPxSolverBase<R>::UNBOUNDED); 558 assert(_solver.basis().status() == SPxBasisBase<R>::UNBOUNDED 559 || _solver.basis().status() == SPxBasisBase<R>::NO_PROBLEM 560 || status() != SPxSolverBase<R>::UNBOUNDED); 561 562 _solReal._isPrimalFeasible = (status() == SPxSolverBase<R>::OPTIMAL 563 || ((_solver.basis().status() == SPxBasisBase<R>::PRIMAL 564 || _solver.basis().status() == SPxBasisBase<R>::UNBOUNDED) 565 && _solver.shift() < 10.0 * realParam(SoPlexBase<R>::EPSILON_ZERO))); 566 567 _solReal._hasPrimalRay = (status() == SPxSolverBase<R>::UNBOUNDED && _isRealLPLoaded); 568 569 // check dual status consistency and query solution status 570 assert(_solver.basis().status() != SPxBasisBase<R>::DUAL || status() != SPxSolverBase<R>::ERROR); 571 assert(_solver.basis().status() != SPxBasisBase<R>::DUAL 572 || status() != SPxSolverBase<R>::NO_RATIOTESTER); 573 assert(_solver.basis().status() != SPxBasisBase<R>::DUAL 574 || status() != SPxSolverBase<R>::NO_PRICER); 575 assert(_solver.basis().status() != SPxBasisBase<R>::DUAL 576 || status() != SPxSolverBase<R>::NO_SOLVER); 577 assert(_solver.basis().status() != SPxBasisBase<R>::DUAL || status() != SPxSolverBase<R>::NOT_INIT); 578 assert(_solver.basis().status() != SPxBasisBase<R>::DUAL || status() != SPxSolverBase<R>::SINGULAR); 579 assert(_solver.basis().status() != SPxBasisBase<R>::DUAL 580 || status() != SPxSolverBase<R>::NO_PROBLEM); 581 assert(_solver.basis().status() != SPxBasisBase<R>::DUAL 582 || status() != SPxSolverBase<R>::UNBOUNDED); 583 assert(_solver.basis().status() != SPxBasisBase<R>::DUAL 584 || status() != SPxSolverBase<R>::INFEASIBLE); 585 assert(_solver.basis().status() != SPxBasisBase<R>::INFEASIBLE 586 || status() == SPxSolverBase<R>::INFEASIBLE); 587 assert(_solver.basis().status() == SPxBasisBase<R>::INFEASIBLE 588 || _solver.basis().status() == SPxBasisBase<R>::NO_PROBLEM 589 || status() != SPxSolverBase<R>::INFEASIBLE); 590 591 _solReal._isDualFeasible = (status() == SPxSolverBase<R>::OPTIMAL 592 || ((_solver.basis().status() == SPxBasisBase<R>::DUAL 593 || _solver.basis().status() == SPxBasisBase<R>::INFEASIBLE) 594 && _solver.shift() < 10.0 * realParam(SoPlexBase<R>::EPSILON_ZERO))); 595 596 _solReal._hasDualFarkas = (status() == SPxSolverBase<R>::INFEASIBLE && _isRealLPLoaded); 597 598 // get infeasibility or unboundedness proof if available 599 if(_solReal._hasPrimalRay) 600 { 601 _solReal._primalRay.reDim(_solver.nCols(), false); 602 _solver.getPrimalray(_solReal._primalRay); 603 } 604 605 if(_solReal._hasDualFarkas) 606 { 607 _solReal._dualFarkas.reDim(_solver.nRows(), false); 608 _solver.getDualfarkas(_solReal._dualFarkas); 609 } 610 611 // get solution data from the solver; independent of solution status 612 _solver.getBasis(_basisStatusRows.get_ptr(), _basisStatusCols.get_ptr(), 613 _basisStatusRows.size(), _basisStatusCols.size()); 614 615 _solver.getPrimalSol(_solReal._primal); 616 _solver.getSlacks(_solReal._slacks); 617 _solver.getDualSol(_solReal._dual); 618 _solver.getRedCostSol(_solReal._redCost); 619 620 _hasBasis = true; 621 622 // get primal and/or dual objective function value depending on status 623 _solver.forceRecompNonbasicValue(); 624 _solReal._objVal = _solver.objValue(); 625 626 // infeasible solutions shall also be stored and be accessible 627 _hasSolReal = true; 628 629 // unscale vectors 630 if(_solver.isScaled() && !_isRealLPLoaded) 631 _unscaleSolutionReal(_solver, false); 632 633 // get unsimplified solution data from simplifier 634 if(_simplifier) 635 { 636 assert(!_simplifier->isUnsimplified()); 637 assert(_simplifier->result() == SPxSimplifier<R>::OKAY); 638 assert(_realLP != &_solver); 639 640 typename SPxBasisBase<R>::SPxStatus simplifiedBasisStatus = _solver.getBasisStatus(); 641 642 try 643 { 644 // pass solution data of transformed problem to simplifier 645 _simplifier->unsimplify(_solReal._primal, _solReal._dual, _solReal._slacks, _solReal._redCost, 646 _basisStatusRows.get_ptr(), _basisStatusCols.get_ptr(), status() == SPxSolverBase<R>::OPTIMAL); 647 } 648 catch(const SPxException& E) 649 { 650 SPX_MSG_INFO1(spxout, spxout << "Caught exception <" << E.what() << 651 "> during unsimplification. Resolving without simplifier and scaler.\n"); 652 _hasBasis = false; 653 _preprocessAndSolveReal(false); 654 return; 655 } 656 657 // copy unsimplified solution data from simplifier (size and dimension is adapted during copy) 658 _solReal._primal = _simplifier->unsimplifiedPrimal(); 659 _solReal._slacks = _simplifier->unsimplifiedSlacks(); 660 _solReal._dual = _simplifier->unsimplifiedDual(); 661 _solReal._redCost = _simplifier->unsimplifiedRedCost(); 662 663 // overwrite the transformed basis with the unsimplified one 664 _simplifier->getBasis(_basisStatusRows.get_ptr(), _basisStatusCols.get_ptr(), 665 _basisStatusRows.size(), _basisStatusCols.size()); 666 667 // load original problem but don't setup a slack basis 668 _loadRealLP(false); 669 670 // since we presolved the problem, we should not allow access to the dual norms 671 _solver.weightsAreSetup = false; 672 673 assert(_realLP == &_solver); 674 675 // reset basis status 676 _solver.setBasisStatus(simplifiedBasisStatus); 677 678 _solver.setBasis(_basisStatusRows.get_const_ptr(), _basisStatusCols.get_const_ptr()); 679 // load unsimplified basis into solver 680 assert(_basisStatusRows.size() == numRows()); 681 assert(_basisStatusCols.size() == this->numCols()); 682 _hasBasis = true; 683 } 684 685 // load realLP into the solver again (internal scaling was applied) 686 else if(_realLP != &_solver) 687 { 688 assert(_solver.isScaled()); 689 _loadRealLP(false); 690 } 691 692 // unscale stored solution (removes persistent scaling) 693 if(_isRealLPScaled) 694 _unscaleSolutionReal(*_realLP, true); 695 696 // check solution for violations and solve again if necessary 697 if(verify) 698 { 699 if(_status == SPxSolverBase<R>::ABORT_VALUE) 700 _verifyObjLimitReal(); 701 else 702 _verifySolutionReal(); 703 } 704 705 assert(_solver.nCols() == this->numCols()); 706 assert(_solver.nRows() == numRows()); 707 } 708 709 710 711 template <class R> 712 void SoPlexBase<R>::_storeSolutionRealFromPresol() 713 { 714 assert(_simplifier); 715 assert(_simplifier->result() == SPxSimplifier<R>::VANISHED); 716 717 // prepare storage for basis (enough to fit the original basis) 718 _basisStatusRows.reSize(numRows()); 719 _basisStatusCols.reSize(numCols()); 720 721 // prepare storage for the solution data and initialize it to zero 722 _solReal._primal.reDim(numCols(), true); 723 _solReal._slacks.reDim(numRows(), true); 724 _solReal._dual.reDim(numRows(), true); 725 _solReal._redCost.reDim(numCols(), true); 726 727 // load original LP and setup slack basis for unsimplifying 728 _loadRealLP(true); 729 730 // store slack basis 731 _solver.getBasis(_basisStatusRows.get_ptr(), _basisStatusCols.get_ptr(), 732 _basisStatusRows.size(), _basisStatusCols.size()); 733 734 assert(!_simplifier->isUnsimplified()); 735 736 try 737 { 738 // unsimplify basis and solution data 739 _simplifier->unsimplify(_solReal._primal, _solReal._dual, _solReal._slacks, _solReal._redCost, 740 _basisStatusRows.get_ptr(), _basisStatusCols.get_ptr()); 741 742 } 743 catch(const SPxException& E) 744 { 745 SPX_MSG_INFO1(spxout, spxout << "Caught exception <" << E.what() << 746 "> during unsimplification. Resolving without simplifier and scaler.\n"); 747 _preprocessAndSolveReal(false); 748 return; 749 } 750 751 // copy unsimplified solution data from simplifier 752 _solReal._primal = _simplifier->unsimplifiedPrimal(); 753 _solReal._slacks = _simplifier->unsimplifiedSlacks(); 754 _solReal._dual = _simplifier->unsimplifiedDual(); 755 _solReal._redCost = _simplifier->unsimplifiedRedCost(); 756 757 // unscale stored solution (removes persistent scaling) 758 if(_isRealLPScaled) 759 _unscaleSolutionReal(*_realLP, true); 760 761 // compute the original objective function value 762 StableSum<R> objVal(realParam(SoPlexBase<R>::OBJ_OFFSET)); 763 764 for(int i = 0; i < numCols(); ++i) 765 objVal += _solReal._primal[i] * objReal(i); 766 767 _solReal._objVal = R(objVal); 768 769 // store the unsimplified basis 770 _simplifier->getBasis(_basisStatusRows.get_ptr(), _basisStatusCols.get_ptr(), 771 _basisStatusRows.size(), _basisStatusCols.size()); 772 _hasBasis = true; 773 _hasSolReal = true; 774 _solReal._isPrimalFeasible = true; 775 _solReal._isDualFeasible = true; 776 _solver.setBasisStatus(SPxBasisBase<R>::OPTIMAL); 777 778 // check solution for violations and solve again if necessary 779 _verifySolutionReal(); 780 } 781 782 783 784 /// load original LP and possibly setup a slack basis 785 template <class R> 786 void SoPlexBase<R>::_loadRealLP(bool initBasis) 787 { 788 _solver.loadLP(*_realLP, initBasis); 789 _isRealLPLoaded = true; 790 _realLP->~SPxLPBase<R>(); 791 spx_free(_realLP); 792 _realLP = &_solver; 793 794 if(initBasis) 795 _solver.init(); 796 } 797 798 799 800 /// unscales stored solution to remove internal or external scaling of LP 801 template <class R> 802 void SoPlexBase<R>::_unscaleSolutionReal(SPxLPBase<R>& LP, bool persistent) 803 { 804 SPX_MSG_INFO1(spxout, spxout << " --- unscaling " << (persistent ? "external" : "internal") << 805 " solution" << std::endl;) 806 assert(_scaler); 807 assert(!persistent || (boolParam(SoPlexBase<R>::PERSISTENTSCALING) && _isRealLPScaled)); 808 _scaler->unscalePrimal(LP, _solReal._primal); 809 _scaler->unscaleSlacks(LP, _solReal._slacks); 810 _scaler->unscaleDual(LP, _solReal._dual); 811 _scaler->unscaleRedCost(LP, _solReal._redCost); 812 813 if(_solReal.hasPrimalRay()) 814 _scaler->unscalePrimalray(LP, _solReal._primalRay); 815 816 if(_solReal.hasDualFarkas()) 817 _scaler->unscaleDualray(LP, _solReal._dualFarkas); 818 } 819 } // namespace soplex 820