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 <assert.h> 26 #include <iostream> 27 28 #include "soplex/spxdefines.h" 29 #include "soplex/rational.h" 30 #include "soplex/spxsolver.h" 31 #include "soplex/spxpricer.h" 32 #include "soplex/spxratiotester.h" 33 #include "soplex/spxdefaultrt.h" 34 #include "soplex/spxstarter.h" 35 #include "soplex/spxout.h" 36 37 #define SOPLEX_MAXCYCLES 400 38 #define SOPLEX_MAXSTALLS 10000 39 #define SOPLEX_MAXSTALLRECOVERS 10 40 #define SOPLEX_MAXREFACPIVOTS 10 41 42 namespace soplex 43 { 44 45 /**@todo check separately for ENTER and LEAVE algorithm */ 46 template <class R> 47 bool SPxSolverBase<R>::precisionReached(R& newpricertol) const 48 { 49 R maxViolRedCost; 50 R sumViolRedCost; 51 R maxViolBounds; 52 R sumViolBounds; 53 R maxViolConst; 54 R sumViolConst; 55 56 qualRedCostViolation(maxViolRedCost, sumViolRedCost); 57 qualBoundViolation(maxViolBounds, sumViolBounds); 58 qualConstraintViolation(maxViolConst, sumViolConst); 59 60 // is the solution good enough ? 61 bool reached = maxViolRedCost < tolerances()->floatingPointOpttol() 62 && maxViolBounds < tolerances()->floatingPointFeastol() 63 && maxViolConst < tolerances()->floatingPointFeastol(); 64 65 if(!reached) 66 { 67 newpricertol = thepricer->pricingTolerance() / 10.0; 68 69 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "Precision not reached: Pricer tolerance = " 70 << thepricer->pricingTolerance() 71 << " new tolerance = " << newpricertol 72 << std::endl 73 << " maxViolRedCost= " << maxViolRedCost 74 << " maxViolBounds= " << maxViolBounds 75 << " maxViolConst= " << maxViolConst 76 << std::endl 77 << " sumViolRedCost= " << sumViolRedCost 78 << " sumViolBounds= " << sumViolBounds 79 << " sumViolConst= " << sumViolConst 80 << std::endl;); 81 } 82 83 return reached; 84 } 85 86 template <class R> 87 void SPxSolverBase<R>::calculateProblemRanges() 88 { 89 // only collect absolute values 90 R minobj = R(infinity); 91 R maxobj = 0.0; 92 R minbound = R(infinity); 93 R maxbound = 0.0; 94 R minside = R(infinity); 95 R maxside = 0.0; 96 97 // get min and max absolute values of bounds and objective 98 for(int j = 0; j < this->nCols(); ++j) 99 { 100 R abslow = spxAbs(this->lower(j)); 101 R absupp = spxAbs(this->lower(j)); 102 R absobj = spxAbs(this->obj(j)); 103 104 if(abslow < R(infinity)) 105 { 106 minbound = SOPLEX_MIN(minbound, abslow); 107 maxbound = SOPLEX_MAX(maxbound, abslow); 108 } 109 110 if(absupp < R(infinity)) 111 { 112 minbound = SOPLEX_MIN(minbound, absupp); 113 maxbound = SOPLEX_MAX(maxbound, absupp); 114 } 115 116 minobj = SOPLEX_MIN(minobj, absobj); 117 maxobj = SOPLEX_MAX(maxobj, absobj); 118 } 119 120 // get min and max absoute values of sides 121 for(int i = 0; i < this->nRows(); ++i) 122 { 123 R abslhs = spxAbs(this->lhs(i)); 124 R absrhs = spxAbs(this->rhs(i)); 125 126 if(abslhs > R(infinity)) 127 { 128 minside = SOPLEX_MIN(minside, abslhs); 129 maxside = SOPLEX_MAX(maxside, abslhs); 130 } 131 132 if(absrhs < R(infinity)) 133 { 134 minside = SOPLEX_MIN(minside, absrhs); 135 maxside = SOPLEX_MAX(maxside, absrhs); 136 } 137 } 138 139 boundrange = maxbound - minbound; 140 siderange = maxside - minside; 141 objrange = maxobj - minobj; 142 } 143 144 template <class R> 145 typename SPxSolverBase<R>::Status SPxSolverBase<R>::solve(volatile bool* interrupt, bool polish) 146 { 147 SPxId enterId; 148 int leaveNum; 149 int loopCount = 0; 150 R minShift = R(infinity); 151 int cycleCount = 0; 152 bool priced = false; 153 R lastDelta = 1; 154 155 /* allow clean up step only once */ 156 recomputedVectors = false; 157 158 /* store the last (primal or dual) feasible objective value to recover/abort in case of stalling */ 159 R stallRefValue; 160 R stallRefShift; 161 int stallRefIter; 162 int stallNumRecovers; 163 164 int timesBasisWasStored = 0; 165 /// number of times the current basis was stored in oldBasisStatusRows oldBasisStatusCols 166 /* storeBasisFreqLog : if true, store basis if iterations() == 2^timesBasisWasStored 167 else, store basis if iterations() % storeBasisSimplexFreq == 0 */ 168 bool storeBasisFreqLog = true; 169 170 if(dim() <= 0 && coDim() <= 0) // no problem loaded 171 { 172 m_status = NO_PROBLEM; 173 throw SPxStatusException("XSOLVE01 No Problem loaded"); 174 } 175 176 if(slinSolver() == 0) // linear system solver is required. 177 { 178 m_status = NO_SOLVER; 179 throw SPxStatusException("XSOLVE02 No Solver loaded"); 180 } 181 182 if(thepricer == 0) // pricer is required. 183 { 184 m_status = NO_PRICER; 185 throw SPxStatusException("XSOLVE03 No Pricer loaded"); 186 } 187 188 if(theratiotester == 0) // ratiotester is required. 189 { 190 m_status = NO_RATIOTESTER; 191 throw SPxStatusException("XSOLVE04 No RatioTester loaded"); 192 } 193 194 theTime->reset(); 195 theTime->start(); 196 197 m_numCycle = 0; 198 this->iterCount = 0; 199 this->lastIterCount = 0; 200 this->iterDegenCheck = 0; 201 202 if(!isInitialized()) 203 { 204 /* 205 if(SPxBasisBase<R>::status() <= NO_PROBLEM) 206 SPxBasisBase<R>::load(this); 207 */ 208 /**@todo != REGULAR is not enough. Also OPTIMAL/DUAL/PRIMAL should 209 * be tested and acted accordingly. 210 */ 211 if(thestarter != 0 && status() != REGULAR 212 && this->theLP->status() == NO_PROBLEM) // no basis and no starter. 213 thestarter->generate(*this); // generate start basis. 214 215 init(); 216 217 // Inna/Tobi: init might fail, if the basis is singular 218 if(!isInitialized()) 219 { 220 assert(SPxBasisBase<R>::status() == SPxBasisBase<R>::SINGULAR); 221 m_status = UNKNOWN; 222 return status(); 223 } 224 } 225 226 //setType(type()); 227 228 if(!this->matrixIsSetup) 229 SPxBasisBase<R>::load(this); 230 231 //factorized = false; 232 233 assert(thepricer->solver() == this); 234 assert(theratiotester->solver() == this); 235 236 // maybe this should be done in init() ? 237 thepricer->setType(type()); 238 theratiotester->setType(type()); 239 240 SPX_MSG_INFO3((*this->spxout), 241 (*this->spxout) << "starting value = " << value() << std::endl 242 << "starting shift = " << shift() << std::endl; 243 ) 244 245 if(SPxBasisBase<R>::status() == SPxBasisBase<R>::OPTIMAL) 246 setBasisStatus(SPxBasisBase<R>::REGULAR); 247 248 m_status = RUNNING; 249 bool stop = terminate(); 250 leaveCount = 0; 251 enterCount = 0; 252 primalCount = 0; 253 polishCount = 0; 254 boundflips = 0; 255 totalboundflips = 0; 256 enterCycles = 0; 257 leaveCycles = 0; 258 primalDegenSum = 0; 259 dualDegenSum = 0; 260 261 multSparseCalls = 0; 262 multFullCalls = 0; 263 multColwiseCalls = 0; 264 multUnsetupCalls = 0; 265 266 stallNumRecovers = 0; 267 268 /* if we run into a singular basis, we will retry from regulardesc with tighter tolerance in the ratio test */ 269 typename SPxSolverBase<R>::Type tightenedtype = type(); 270 bool tightened = false; 271 272 oldBasisStatusRows.reSize(this->nRows()); 273 oldBasisStatusCols.reSize(this->nCols()); 274 275 while(!stop) 276 { 277 const typename SPxBasisBase<R>::Desc regulardesc = this->desc(); 278 279 // we need to reset these pointers to avoid unnecessary/wrong solves in leave() or enter() 280 solveVector2 = 0; 281 solveVector3 = 0; 282 coSolveVector2 = 0; 283 coSolveVector3 = 0; 284 285 updateViols.clear(); 286 updateViolsCo.clear(); 287 288 try 289 { 290 291 if(type() == ENTER) 292 { 293 forceRecompNonbasicValue(); 294 295 int enterCycleCount = 0; 296 int enterFacPivotCount = 0; 297 298 instableEnterVal = 0; 299 instableEnterId = SPxId(); 300 instableEnter = false; 301 302 stallRefIter = this->iteration() - 1; 303 stallRefShift = shift(); 304 stallRefValue = value(); 305 306 /* in the entering algorithm, entertol() should be maintained by the ratio test and leavetol() should be 307 * reached by the pricer 308 */ 309 R maxpricertol = leavetol(); 310 R minpricertol = 0.01 * maxpricertol; 311 312 thepricer->setPricingTolerance(maxpricertol); 313 priced = false; 314 315 // to avoid shifts we restrict tolerances in the ratio test 316 if(loopCount > 0) 317 { 318 lastDelta = (lastDelta < entertol()) ? lastDelta : entertol(); 319 lastDelta *= 0.01; 320 theratiotester->setDelta(lastDelta); 321 assert(theratiotester->getDelta() > 0); 322 SPxOut::debug(this, "decreased delta for ratiotest to: {}\n", theratiotester->getDelta()); 323 } 324 else 325 { 326 lastDelta = 1; 327 theratiotester->setDelta(entertol()); 328 } 329 330 printDisplayLine(true); 331 332 do 333 { 334 printDisplayLine(); 335 336 // if it is time to store the basis, store it (only used in rational solve) 337 if(solvingForBoosted) 338 { 339 if((storeBasisFreqLog && iterations() == pow(2, timesBasisWasStored) - 1) 340 || (!storeBasisFreqLog && iterations() % storeBasisSimplexFreq == 0)) 341 { 342 // switch off storeBasisFreqLog if 2^timesBasisWasStored becomes too big 343 // in order to avoid computing enormous powers of 2 344 if(storeBasisFreqLog && pow(2, timesBasisWasStored) > storeBasisSimplexFreq) 345 storeBasisFreqLog = false; 346 347 // store basis 348 getBasis(oldBasisStatusRows.get_ptr(), oldBasisStatusCols.get_ptr(), oldBasisStatusRows.size(), 349 oldBasisStatusCols.size()); 350 timesBasisWasStored++; 351 } 352 } 353 354 enterId = thepricer->selectEnter(); 355 356 if(!enterId.isValid() && instableEnterId.isValid() && this->lastUpdate() == 0) 357 { 358 /* no entering variable was found, but because of valid instableEnterId we know 359 that this is due to the scaling of the test values. Thus, we use 360 instableEnterId and SPxFastRT<R>::selectEnter shall accept even an instable 361 leaving variable. */ 362 SPX_MSG_INFO3((*this->spxout), 363 (*this->spxout) << " --- trying instable enter iteration" << std::endl;) 364 365 enterId = instableEnterId; 366 instableEnter = true; 367 // we also need to reset the test() or coTest() value for getEnterVals() 368 assert(instableEnterVal < 0); 369 370 if(enterId.isSPxColId()) 371 { 372 int idx = this->number(SPxColId(enterId)); 373 374 if(rep() == COLUMN) 375 { 376 theTest[idx] = instableEnterVal; 377 378 if(sparsePricingEnterCo && isInfeasibleCo[idx] == SPxPricer<R>::NOT_VIOLATED) 379 { 380 infeasibilitiesCo.addIdx(idx); 381 isInfeasibleCo[idx] = SPxPricer<R>::VIOLATED; 382 } 383 384 if(hyperPricingEnter) 385 updateViolsCo.addIdx(idx); 386 } 387 else 388 { 389 theCoTest[idx] = instableEnterVal; 390 391 if(sparsePricingEnter && isInfeasible[idx] == SPxPricer<R>::NOT_VIOLATED) 392 { 393 infeasibilities.addIdx(idx); 394 isInfeasible[idx] = SPxPricer<R>::VIOLATED; 395 } 396 397 if(hyperPricingEnter) 398 updateViols.addIdx(idx); 399 } 400 } 401 else 402 { 403 int idx = this->number(SPxRowId(enterId)); 404 405 if(rep() == COLUMN) 406 { 407 theCoTest[idx] = instableEnterVal; 408 409 if(sparsePricingEnter && isInfeasible[idx] == SPxPricer<R>::NOT_VIOLATED) 410 { 411 infeasibilities.addIdx(idx); 412 isInfeasible[idx] = SPxPricer<R>::VIOLATED; 413 } 414 415 if(hyperPricingEnter) 416 updateViols.addIdx(idx); 417 } 418 else 419 { 420 theTest[idx] = instableEnterVal; 421 422 if(sparsePricingEnterCo && isInfeasibleCo[idx] == SPxPricer<R>::NOT_VIOLATED) 423 { 424 infeasibilitiesCo.addIdx(idx); 425 isInfeasibleCo[idx] = SPxPricer<R>::VIOLATED; 426 } 427 428 if(hyperPricingEnter) 429 updateViolsCo.addIdx(idx); 430 } 431 } 432 } 433 else 434 { 435 instableEnter = false; 436 } 437 438 if(!enterId.isValid()) 439 { 440 // we are not infeasible and have no shift 441 if(shift() <= epsilon() 442 && (SPxBasisBase<R>::status() == SPxBasisBase<R>::REGULAR 443 || SPxBasisBase<R>::status() == SPxBasisBase<R>::DUAL 444 || SPxBasisBase<R>::status() == SPxBasisBase<R>::PRIMAL)) 445 { 446 R newpricertol = minpricertol; 447 448 // refactorize to eliminate accumulated errors from LU updates 449 if(this->lastUpdate() > 0) 450 factorize(); 451 452 // recompute Fvec, Pvec and CoPvec to get a more precise solution and obj value 453 computeFrhs(); 454 SPxBasisBase<R>::solve(*theFvec, *theFrhs); 455 456 computeEnterCoPrhs(); 457 SPxBasisBase<R>::coSolve(*theCoPvec, *theCoPrhs); 458 computePvec(); 459 460 forceRecompNonbasicValue(); 461 462 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << " --- checking feasibility and optimality\n") 463 computeCoTest(); 464 computeTest(); 465 466 // is the solution good enough ? 467 // max three times reduced 468 if((thepricer->pricingTolerance() > minpricertol) && !precisionReached(newpricertol)) 469 { 470 // no! 471 // we reduce the pricer tolerance. Note that if the pricer does not find a candiate 472 // with the reduced tolerance, we quit, regardless of the violations. 473 if(newpricertol < minpricertol) 474 newpricertol = minpricertol; 475 476 thepricer->setPricingTolerance(newpricertol); 477 478 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << " --- setting pricer tolerance to " 479 << thepricer->pricingTolerance() 480 << std::endl;) 481 } 482 } 483 484 // if the factorization is not fresh, we better refactorize and call the pricer again; however, this can 485 // create cycling, so it is performed only a limited number of times per ENTER round 486 if(this->lastUpdate() > 0 && enterFacPivotCount < SOPLEX_MAXREFACPIVOTS) 487 { 488 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << " --- solve(enter) triggers refactorization" << 489 std::endl;) 490 491 factorize(); 492 493 // if the factorization was found out to be singular, we have to quit 494 if(SPxBasisBase<R>::status() < SPxBasisBase<R>::REGULAR) 495 { 496 SPX_MSG_INFO1((*this->spxout), 497 (*this->spxout) << "Something wrong with factorization, Basis status: " 498 << static_cast<int>(SPxBasisBase<R>::status()) << std::endl;) 499 stop = true; 500 break; 501 } 502 503 // call pricer again 504 enterId = thepricer->selectEnter(); 505 506 // count how often the pricer has found something only after refactorizing 507 if(enterId.isValid()) 508 enterFacPivotCount++; 509 } 510 511 if(!enterId.isValid()) 512 { 513 priced = true; 514 break; 515 } 516 } 517 518 /* check if we have iterations left */ 519 if(maxIters >= 0 && iterations() >= maxIters) 520 { 521 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << " --- maximum number of iterations (" << maxIters 522 << ") reached" << std::endl;) 523 m_status = ABORT_ITER; 524 stop = true; 525 break; 526 } 527 528 if(interrupt != NULL && *interrupt) 529 { 530 SPX_MSG_INFO2((*this->spxout), 531 (*this->spxout) << " --- aborted due to interrupt signal" << std::endl;) 532 m_status = ABORT_TIME; 533 stop = true; 534 break; 535 } 536 537 enter(enterId); 538 assert((testBounds(), 1)); 539 thepricer->entered4(this->lastEntered(), this->lastIndex()); 540 stop = terminate(); 541 clearUpdateVecs(); 542 543 /* if a successful pivot was performed or a nonbasic variable was flipped to its other bound, we reset the 544 * cycle counter 545 */ 546 if(this->lastEntered().isValid()) 547 enterCycleCount = 0; 548 else if(basis().status() != SPxBasisBase<R>::INFEASIBLE 549 && basis().status() != SPxBasisBase<R>::UNBOUNDED) 550 { 551 enterCycleCount++; 552 553 if(enterCycleCount > SOPLEX_MAXCYCLES) 554 { 555 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << " --- abort solving due to cycling in " 556 << "entering algorithm" << std::endl;); 557 m_status = ABORT_CYCLING; 558 stop = true; 559 } 560 } 561 562 /* only if the basis has really changed, we increase the iterations counter; this is not the case when only 563 * a nonbasic variable was flipped to its other bound 564 */ 565 if(this->lastIndex() >= 0) 566 { 567 enterCount++; 568 assert(this->lastEntered().isValid()); 569 } 570 571 /* check every SOPLEX_MAXSTALLS iterations whether shift and objective value have not changed */ 572 if((this->iteration() - stallRefIter) % SOPLEX_MAXSTALLS == 0 573 && basis().status() != SPxBasisBase<R>::INFEASIBLE) 574 { 575 if(spxAbs(value() - stallRefValue) <= epsilon() && spxAbs(shift() - stallRefShift) <= epsilon()) 576 { 577 if(stallNumRecovers < SOPLEX_MAXSTALLRECOVERS) 578 { 579 /* try to recover by unshifting/switching algorithm up to SOPLEX_MAXSTALLRECOVERS times (just a number picked) */ 580 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << 581 " --- stalling detected - trying to recover by switching to LEAVING algorithm." << std::endl;) 582 583 ++stallNumRecovers; 584 break; 585 } 586 else 587 { 588 /* giving up */ 589 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << 590 " --- abort solving due to stalling in entering algorithm." << std::endl;); 591 592 m_status = ABORT_CYCLING; 593 stop = true; 594 } 595 } 596 else 597 { 598 /* merely update reference values */ 599 stallRefIter = this->iteration() - 1; 600 stallRefShift = shift(); 601 stallRefValue = value(); 602 } 603 } 604 605 //@ assert(isConsistent()); 606 } 607 while(!stop); 608 609 SPX_MSG_INFO3((*this->spxout), 610 (*this->spxout) << " --- enter finished. iteration: " << this->iteration() 611 << ", value: " << value() 612 << ", shift: " << shift() 613 << ", epsilon: " << epsilon() 614 << ", feastol: " << tolerances()->floatingPointFeastol() 615 << ", opttol: " << tolerances()->floatingPointOpttol() 616 << std::endl 617 << "ISOLVE56 stop: " << stop 618 << ", basis status: " << static_cast<int>(SPxBasisBase<R>::status()) << " (" << static_cast<int> 619 (SPxBasisBase<R>::status()) << ")" 620 << ", solver status: " << static_cast<int>(m_status) << " (" << static_cast<int> 621 (m_status) << ")" << std::endl; 622 ) 623 624 if(!stop) 625 { 626 /**@todo technically it would be ok to finish already when (priced && maxinfeas + shift() <= entertol()) is 627 * satisfied; maybe at least in the case when SoPlex keeps jumping back between ENTER and LEAVE always 628 * shifting (looping), we may relax this condition here; 629 * note also that unShift may increase shift() slightly due to roundoff errors 630 */ 631 if(shift() <= epsilon()) 632 { 633 // factorize(); 634 unShift(); 635 636 R maxinfeas = maxInfeas(); 637 638 SPX_MSG_INFO3((*this->spxout), 639 (*this->spxout) << " --- maxInfeas: " << maxinfeas 640 << ", shift: " << shift() 641 << ", entertol: " << entertol() << std::endl; 642 ) 643 644 if(priced && maxinfeas + shift() <= entertol()) 645 { 646 setBasisStatus(SPxBasisBase<R>::OPTIMAL); 647 m_status = OPTIMAL; 648 break; 649 } 650 else if(loopCount > 2) 651 { 652 // calculate problem ranges if not done already 653 if(boundrange == 0.0 || siderange == 0.0 || objrange == 0.0) 654 calculateProblemRanges(); 655 656 if(SOPLEX_MAX(SOPLEX_MAX(boundrange, siderange), objrange) >= 1e9 && !solvingForBoosted) 657 { 658 SPxOut::setScientific(spxout->getCurrentStream(), 0); 659 SPX_MSG_INFO1((*this->spxout), (*this->spxout) << 660 " --- termination despite violations (numerical difficulties," 661 << " bound range = " << boundrange 662 << ", side range = " << siderange 663 << ", obj range = " << objrange 664 << ")" << std::endl;) 665 setBasisStatus(SPxBasisBase<R>::OPTIMAL); 666 m_status = OPTIMAL; 667 break; 668 } 669 else 670 { 671 m_status = ABORT_CYCLING; 672 throw SPxStatusException("XSOLVE14 Abort solving due to looping"); 673 } 674 } 675 676 loopCount++; 677 } 678 679 setType(LEAVE); 680 init(); 681 thepricer->setType(type()); 682 theratiotester->setType(type()); 683 } 684 } 685 else 686 { 687 assert(type() == LEAVE); 688 689 forceRecompNonbasicValue(); 690 691 int leaveCycleCount = 0; 692 int leaveFacPivotCount = 0; 693 694 instableLeaveNum = -1; 695 instableLeave = false; 696 instableLeaveVal = 0; 697 698 stallRefIter = this->iteration() - 1; 699 stallRefShift = shift(); 700 stallRefValue = value(); 701 702 /* in the leaving algorithm, leavetol() should be maintained by the ratio test and entertol() should be reached 703 * by the pricer 704 */ 705 R maxpricertol = entertol(); 706 R minpricertol = 0.01 * maxpricertol; 707 708 thepricer->setPricingTolerance(maxpricertol); 709 priced = false; 710 711 // to avoid shifts we restrict tolerances in the ratio test 712 if(loopCount > 0) 713 { 714 lastDelta = (lastDelta < leavetol()) ? lastDelta : leavetol(); 715 lastDelta *= 0.01; 716 theratiotester->setDelta(lastDelta); 717 assert(theratiotester->getDelta() > 0); 718 SPxOut::debug(this, "decreased delta for ratiotest to: {}\n", theratiotester->getDelta()); 719 } 720 else 721 { 722 lastDelta = 1; 723 theratiotester->setDelta(leavetol()); 724 } 725 726 printDisplayLine(true); 727 728 do 729 { 730 printDisplayLine(); 731 732 // if it is time to store the basis, store it (only used in rational solve) 733 if(solvingForBoosted) 734 { 735 if((storeBasisFreqLog && iterations() == pow(2, timesBasisWasStored) - 1) 736 || (!storeBasisFreqLog && iterations() % storeBasisSimplexFreq == 0)) 737 { 738 // switch off storeBasisFreqLog if 2^timesBasisWasStored becomes too big 739 // in order to avoid computing enormous powers of 2 740 if(storeBasisFreqLog && pow(2, timesBasisWasStored) > storeBasisSimplexFreq) 741 storeBasisFreqLog = false; 742 743 // store basis 744 getBasis(oldBasisStatusRows.get_ptr(), oldBasisStatusCols.get_ptr(), oldBasisStatusRows.size(), 745 oldBasisStatusCols.size()); 746 timesBasisWasStored++; 747 } 748 } 749 750 leaveNum = thepricer->selectLeave(); 751 752 if(leaveNum < 0 && instableLeaveNum >= 0 && this->lastUpdate() == 0) 753 { 754 /* no leaving variable was found, but because of instableLeaveNum >= 0 we know 755 that this is due to the scaling of theCoTest[...]. Thus, we use 756 instableLeaveNum and SPxFastRT<R>::selectEnter shall accept even an instable 757 entering variable. */ 758 SPX_MSG_INFO3((*this->spxout), 759 (*this->spxout) << " --- trying instable leave iteration" << std::endl; 760 ) 761 762 leaveNum = instableLeaveNum; 763 instableLeave = true; 764 // we also need to reset the fTest() value for getLeaveVals() 765 assert(instableLeaveVal < 0); 766 theCoTest[instableLeaveNum] = instableLeaveVal; 767 768 if(sparsePricingLeave) 769 { 770 if(isInfeasible[instableLeaveNum] == SPxPricer<R>::NOT_VIOLATED) 771 { 772 infeasibilities.addIdx(instableLeaveNum); 773 isInfeasible[instableLeaveNum] = SPxPricer<R>::VIOLATED; 774 } 775 776 if(hyperPricingLeave) 777 updateViols.addIdx(instableLeaveNum); 778 } 779 } 780 else 781 { 782 instableLeave = false; 783 } 784 785 if(leaveNum < 0) 786 { 787 // we are not infeasible and have no shift 788 if(shift() <= epsilon() 789 && (SPxBasisBase<R>::status() == SPxBasisBase<R>::REGULAR 790 || SPxBasisBase<R>::status() == SPxBasisBase<R>::DUAL 791 || SPxBasisBase<R>::status() == SPxBasisBase<R>::PRIMAL)) 792 { 793 R newpricertol = minpricertol; 794 795 // refactorize to eliminate accumulated errors from LU updates 796 if(this->lastUpdate() > 0) 797 factorize(); 798 799 // recompute Fvec, Pvec and CoPvec to get a more precise solution and obj value 800 computeFrhs(); 801 SPxBasisBase<R>::solve(*theFvec, *theFrhs); 802 803 computeLeaveCoPrhs(); 804 SPxBasisBase<R>::coSolve(*theCoPvec, *theCoPrhs); 805 computePvec(); 806 807 forceRecompNonbasicValue(); 808 809 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << " --- checking feasibility and optimality\n") 810 computeFtest(); 811 812 // is the solution good enough ? 813 // max three times reduced 814 if((thepricer->pricingTolerance() > minpricertol) && !precisionReached(newpricertol)) 815 { 816 // no 817 // we reduce the pricer tolerance. Note that if the pricer does not find a candiate 818 // with the reduced pricer tolerance, we quit, regardless of the violations. 819 if(newpricertol < minpricertol) 820 newpricertol = minpricertol; 821 822 thepricer->setPricingTolerance(newpricertol); 823 824 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << " --- setting pricer tolerance to " 825 << thepricer->pricingTolerance() << std::endl;); 826 } 827 } 828 829 // if the factorization is not fresh, we better refactorize and call the pricer again; however, this can 830 // create cycling, so it is performed only a limited number of times per LEAVE round 831 if(this->lastUpdate() > 0 && leaveFacPivotCount < SOPLEX_MAXREFACPIVOTS) 832 { 833 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << " --- solve(leave) triggers refactorization" << 834 std::endl;) 835 836 factorize(); 837 838 // Inna/Tobi: if the factorization was found out to be singular, we have to quit 839 if(SPxBasisBase<R>::status() < SPxBasisBase<R>::REGULAR) 840 { 841 SPX_MSG_INFO1((*this->spxout), 842 (*this->spxout) << "Something wrong with factorization, Basis status: " 843 << static_cast<int>(SPxBasisBase<R>::status()) << std::endl;) 844 stop = true; 845 break; 846 } 847 848 // call pricer again 849 leaveNum = thepricer->selectLeave(); 850 851 // count how often the pricer has found something only after refactorizing 852 if(leaveNum >= 0) 853 leaveFacPivotCount++; 854 } 855 856 if(leaveNum < 0) 857 { 858 priced = true; 859 break; 860 } 861 } 862 863 /* check if we have iterations left */ 864 if(maxIters >= 0 && iterations() >= maxIters) 865 { 866 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << " --- maximum number of iterations (" << maxIters 867 << ") reached" << std::endl;) 868 m_status = ABORT_ITER; 869 stop = true; 870 break; 871 } 872 873 if(interrupt != NULL && *interrupt) 874 { 875 SPX_MSG_INFO2((*this->spxout), 876 (*this->spxout) << " --- aborted due to interrupt signal" << std::endl;) 877 m_status = ABORT_TIME; 878 stop = true; 879 break; 880 } 881 882 leave(leaveNum); 883 assert((testBounds(), 1)); 884 thepricer->left4(this->lastIndex(), this->lastLeft()); 885 stop = terminate(); 886 clearUpdateVecs(); 887 888 /* if a successful pivot was performed or a nonbasic variable was flipped to its other bound, we reset the 889 * cycle counter 890 */ 891 if(this->lastIndex() >= 0) 892 leaveCycleCount = 0; 893 else if(basis().status() != SPxBasisBase<R>::INFEASIBLE 894 && basis().status() != SPxBasisBase<R>::UNBOUNDED) 895 { 896 leaveCycleCount++; 897 898 if(leaveCycleCount > SOPLEX_MAXCYCLES) 899 { 900 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << 901 " --- abort solving due to cycling in leaving algorithm" << std::endl;); 902 m_status = ABORT_CYCLING; 903 stop = true; 904 } 905 } 906 907 /* only if the basis has really changed, we increase the iterations counter; this is not the case when only 908 * a nonbasic variable was flipped to its other bound 909 */ 910 if(this->lastEntered().isValid()) 911 { 912 leaveCount++; 913 assert(this->lastIndex() >= 0); 914 } 915 916 /* check every SOPLEX_MAXSTALLS iterations whether shift and objective value have not changed */ 917 if((this->iteration() - stallRefIter) % SOPLEX_MAXSTALLS == 0 918 && basis().status() != SPxBasisBase<R>::INFEASIBLE) 919 { 920 if(spxAbs(value() - stallRefValue) <= epsilon() && spxAbs(shift() - stallRefShift) <= epsilon()) 921 { 922 if(stallNumRecovers < SOPLEX_MAXSTALLRECOVERS) 923 { 924 /* try to recover by switching algorithm up to SOPLEX_MAXSTALLRECOVERS times */ 925 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << 926 " --- stalling detected - trying to recover by switching to ENTERING algorithm." << std::endl;) 927 928 ++stallNumRecovers; 929 break; 930 } 931 else 932 { 933 /* giving up */ 934 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << 935 " --- abort solving due to stalling in leaving algorithm" << std::endl;); 936 937 m_status = ABORT_CYCLING; 938 stop = true; 939 } 940 } 941 else 942 { 943 /* merely update reference values */ 944 stallRefIter = this->iteration() - 1; 945 stallRefShift = shift(); 946 stallRefValue = value(); 947 } 948 } 949 950 //@ assert(isConsistent()); 951 } 952 while(!stop); 953 954 SPX_MSG_INFO3((*this->spxout), 955 (*this->spxout) << " --- leave finished. iteration: " << this->iteration() 956 << ", value: " << value() 957 << ", shift: " << shift() 958 << ", epsilon: " << epsilon() 959 << ", feastol: " << tolerances()->floatingPointFeastol() 960 << ", opttol: " << tolerances()->floatingPointOpttol() 961 << std::endl 962 << "ISOLVE57 stop: " << stop 963 << ", basis status: " << static_cast<int>(SPxBasisBase<R>::status()) << " (" << static_cast<int> 964 (SPxBasisBase<R>::status()) << ")" 965 << ", solver status: " << static_cast<int>(m_status) << " (" << static_cast<int> 966 (m_status) << ")" << std::endl; 967 ) 968 969 if(!stop) 970 { 971 if(shift() < minShift) 972 { 973 minShift = shift(); 974 cycleCount = 0; 975 } 976 else 977 { 978 cycleCount++; 979 980 if(cycleCount > SOPLEX_MAXCYCLES) 981 { 982 m_status = ABORT_CYCLING; 983 throw SPxStatusException("XSOLVE13 Abort solving due to cycling"); 984 } 985 986 SPX_MSG_INFO3((*this->spxout), 987 (*this->spxout) << " --- maxInfeas: " << maxInfeas() 988 << ", shift: " << shift() 989 << ", leavetol: " << leavetol() 990 << ", cycle count: " << cycleCount << std::endl; 991 ) 992 } 993 994 /**@todo technically it would be ok to finish already when (priced && maxinfeas + shift() <= entertol()) is 995 * satisfied; maybe at least in the case when SoPlex keeps jumping back between ENTER and LEAVE always 996 * shifting (looping), we may relax this condition here; 997 * note also that unShift may increase shift() slightly due to roundoff errors 998 */ 999 if(shift() <= epsilon()) 1000 { 1001 cycleCount = 0; 1002 // factorize(); 1003 unShift(); 1004 1005 R maxinfeas = maxInfeas(); 1006 1007 SPX_MSG_INFO3((*this->spxout), 1008 (*this->spxout) << " --- maxInfeas: " << maxinfeas 1009 << ", shift: " << shift() 1010 << ", leavetol: " << leavetol() << std::endl; 1011 ) 1012 1013 // We stop if we are indeed optimal, or if we have already been 1014 // two times at this place. In this case it seems futile to 1015 // continue. 1016 if(priced && maxinfeas + shift() <= leavetol()) 1017 { 1018 setBasisStatus(SPxBasisBase<R>::OPTIMAL); 1019 m_status = OPTIMAL; 1020 break; 1021 } 1022 else if(loopCount > 2) 1023 { 1024 // calculate problem ranges if not done already 1025 if(boundrange == 0.0 || siderange == 0.0 || objrange == 0.0) 1026 calculateProblemRanges(); 1027 1028 if(SOPLEX_MAX(SOPLEX_MAX(boundrange, siderange), objrange) >= 1e9 && !solvingForBoosted) 1029 { 1030 SPxOut::setScientific(spxout->getCurrentStream(), 0); 1031 SPX_MSG_INFO1((*this->spxout), (*this->spxout) << 1032 " --- termination despite violations (numerical difficulties," 1033 << " bound range = " << boundrange 1034 << ", side range = " << siderange 1035 << ", obj range = " << objrange 1036 << ")" << std::endl;) 1037 setBasisStatus(SPxBasisBase<R>::OPTIMAL); 1038 m_status = OPTIMAL; 1039 break; 1040 } 1041 else 1042 { 1043 m_status = ABORT_CYCLING; 1044 throw SPxStatusException("XSOLVE14 Abort solving due to looping"); 1045 } 1046 } 1047 1048 loopCount++; 1049 } 1050 1051 setType(ENTER); 1052 init(); 1053 thepricer->setType(type()); 1054 theratiotester->setType(type()); 1055 } 1056 } 1057 1058 assert(m_status != SINGULAR); 1059 1060 } 1061 catch(const SPxException& E) 1062 { 1063 // if we stopped due to a singular basis, we reload the original basis and try again with tighter 1064 // tolerance (only once) 1065 if(m_status == SINGULAR && !tightened) 1066 { 1067 tightenedtype = type(); 1068 1069 if(tightenedtype == ENTER) 1070 { 1071 this->scaleEntertol(0.01); 1072 1073 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << 1074 " --- basis singular: reloading basis and solving with tighter ratio test tolerance " << 1075 this->entertol() << std::endl;) 1076 } 1077 else 1078 { 1079 this->scaleLeavetol(0.01); 1080 1081 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << 1082 " --- basis singular: reloading basis and solving with tighter ratio test tolerance " << 1083 this->leavetol() << std::endl;) 1084 } 1085 1086 // load original basis 1087 int niters = iterations(); 1088 loadBasis(regulardesc); 1089 1090 // remember iteration count 1091 this->iterCount = niters; 1092 1093 // try initializing basis (might fail if starting basis was already singular) 1094 try 1095 { 1096 init(); 1097 theratiotester->setType(type()); 1098 } 1099 catch(const SPxException& Ex) 1100 { 1101 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << 1102 " --- reloaded basis singular, resetting original tolerances" << std::endl;) 1103 1104 if(tightenedtype == ENTER) 1105 this->scaleEntertol(1); 1106 else 1107 this->scaleLeavetol(1); 1108 1109 theratiotester->setType(type()); 1110 1111 throw Ex; 1112 } 1113 1114 // reset status and counters 1115 m_status = RUNNING; 1116 m_numCycle = 0; 1117 leaveCount = 0; 1118 enterCount = 0; 1119 stallNumRecovers = 0; 1120 1121 // continue 1122 stop = false; 1123 tightened = true; 1124 } 1125 // reset tolerance to its original value and pass on the exception 1126 else if(tightened) 1127 { 1128 if(tightenedtype == ENTER) 1129 this->scaleEntertol(1); 1130 else 1131 this->scaleLeavetol(1); 1132 1133 theratiotester->setType(type()); 1134 1135 throw E; 1136 } 1137 // pass on the exception 1138 else 1139 throw E; 1140 } 1141 } 1142 1143 // reset tolerance to its original value 1144 if(tightened) 1145 { 1146 if(tightenedtype == ENTER) 1147 this->scaleEntertol(1); 1148 else 1149 this->scaleLeavetol(1); 1150 1151 theratiotester->setType(type()); 1152 } 1153 1154 theTime->stop(); 1155 theCumulativeTime += time(); 1156 1157 if(m_status == RUNNING) 1158 { 1159 m_status = ERROR; 1160 throw SPxStatusException("XSOLVE05 Status is still RUNNING when it shouldn't be"); 1161 } 1162 1163 SPX_MSG_INFO3((*this->spxout), 1164 (*this->spxout) << "Finished solving (status=" << static_cast<int>(status()) 1165 << ", iters=" << this->iterCount 1166 << ", leave=" << leaveCount 1167 << ", enter=" << enterCount 1168 << ", flips=" << totalboundflips; 1169 1170 if(status() == OPTIMAL) 1171 (*this->spxout) << ", objValue=" << value(); 1172 (*this->spxout) << ")" << std::endl; 1173 ) 1174 1175 #ifdef ENABLE_ADDITIONAL_CHECKS 1176 1177 /* check if solution is really feasible */ 1178 if(status() == OPTIMAL) 1179 { 1180 int c; 1181 R val; 1182 VectorBase<R> sol(this->nCols()); 1183 1184 getPrimalSol(sol); 1185 1186 for(int row = 0; row < this->nRows(); ++row) 1187 { 1188 const SVectorBase<R>& rowvec = this->rowVector(row); 1189 val = 0.0; 1190 1191 for(c = 0; c < rowvec.size(); ++c) 1192 val += rowvec.value(c) * sol[rowvec.index(c)]; 1193 1194 if(LT(val, this->lhs(row), tolerances()->floatingPointFeastol()) || 1195 GT(val, this->rhs(row), tolerances()->floatingPointFeastol())) 1196 { 1197 // Minor rhs violations happen frequently, so print these 1198 // warnings only with verbose level INFO2 and higher. 1199 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "WSOLVE88 Warning! Constraint " << row 1200 << " is violated by solution" << std::endl 1201 << " lhs:" << this->lhs(row) 1202 << " <= val:" << val 1203 << " <= rhs:" << this->rhs(row) << std::endl;) 1204 1205 if(type() == LEAVE && isRowBasic(row)) 1206 { 1207 // find basis variable 1208 for(c = 0; c < this->nRows(); ++c) 1209 if(basis().baseId(c).isSPxRowId() 1210 && (this->number(basis().baseId(c)) == row)) 1211 break; 1212 1213 assert(c < this->nRows()); 1214 1215 SPX_MSG_WARNING((*this->spxout), (*this->spxout) << "WSOLVE90 basis idx:" << c 1216 << " fVec:" << fVec()[c] 1217 << " fRhs:" << fRhs()[c] 1218 << " fTest:" << fTest()[c] << std::endl;) 1219 } 1220 } 1221 } 1222 1223 for(int col = 0; col < this->nCols(); ++col) 1224 { 1225 if(LT(sol[col], this->lower(col), tolerances()->floatingPointFeastol()) || 1226 GT(sol[col], this->upper(col), tolerances()->floatingPointFeastol())) 1227 { 1228 // Minor bound violations happen frequently, so print these 1229 // warnings only with verbose level INFO2 and higher. 1230 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "WSOLVE91 Warning! Bound for column " << col 1231 << " is violated by solution" << std::endl 1232 << " lower:" << this->lower(col) 1233 << " <= val:" << sol[col] 1234 << " <= upper:" << this->upper(col) << std::endl;) 1235 1236 if(type() == LEAVE && isColBasic(col)) 1237 { 1238 for(c = 0; c < this->nRows() ; ++c) 1239 if(basis().baseId(c).isSPxColId() 1240 && (this->number(basis().baseId(c)) == col)) 1241 break; 1242 1243 assert(c < this->nRows()); 1244 SPX_MSG_WARNING((*this->spxout), (*this->spxout) << "WSOLVE92 basis idx:" << c 1245 << " fVec:" << fVec()[c] 1246 << " fRhs:" << fRhs()[c] 1247 << " fTest:" << fTest()[c] << std::endl;) 1248 } 1249 } 1250 } 1251 } 1252 1253 #endif // ENABLE_ADDITIONAL_CHECKS 1254 1255 1256 primalCount = (rep() == SPxSolverBase<R>::COLUMN) 1257 ? enterCount 1258 : leaveCount; 1259 1260 printDisplayLine(true); 1261 1262 if(polish) 1263 { 1264 bool resolve; 1265 resolve = performSolutionPolishing(); 1266 1267 if(resolve) 1268 solve(interrupt, false); 1269 } 1270 1271 return status(); 1272 } 1273 1274 template <class R> 1275 bool SPxSolverBase<R>::performSolutionPolishing() 1276 { 1277 // catch rare case that the iteration limit is exactly reached at optimality 1278 bool stop = (maxIters >= 0 && iterations() >= maxIters && !isTimeLimitReached()); 1279 1280 // only polish an already optimal basis 1281 if(stop || polishObj == POLISH_OFF || status() != OPTIMAL) 1282 return false; 1283 1284 int nSuccessfulPivots; 1285 const typename SPxBasisBase<R>::Desc& ds = this->desc(); 1286 const typename SPxBasisBase<R>::Desc::Status* rowstatus = ds.rowStatus(); 1287 const typename SPxBasisBase<R>::Desc::Status* colstatus = ds.colStatus(); 1288 typename SPxBasisBase<R>::Desc::Status stat; 1289 SPxId polishId; 1290 bool success = false; 1291 1292 R alloweddeviation; 1293 R origval = value(); 1294 1295 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << " --- perform solution polishing" << std::endl;) 1296 1297 if(rep() == COLUMN) 1298 { 1299 setType(ENTER); // use primal simplex to preserve feasibility 1300 init(); 1301 alloweddeviation = entertol(); 1302 1303 instableEnter = false; 1304 theratiotester->setType(type()); 1305 1306 int nrows = this->nRows(); 1307 int ncols = this->nCols(); 1308 1309 if(polishObj == POLISH_INTEGRALITY) 1310 { 1311 DIdxSet slackcandidates(nrows); 1312 DIdxSet continuousvars(ncols); 1313 1314 // collect nonbasic slack variables that could be made basic 1315 for(int i = 0; i < nrows; ++i) 1316 { 1317 // only check nonbasic rows, skip equations 1318 if(rowstatus[i] == SPxBasisBase<R>::Desc::P_ON_LOWER 1319 || rowstatus[i] == SPxBasisBase<R>::Desc::P_ON_UPPER) 1320 { 1321 // only consider rows with zero dual multiplier to preserve optimality 1322 if(EQrel((*theCoPvec)[i], (R) 0, this->epsilon())) 1323 slackcandidates.addIdx(i); 1324 } 1325 } 1326 1327 // collect continuous variables that could be made basic 1328 if(integerVariables.size() == ncols) 1329 { 1330 for(int i = 0; i < ncols; ++i) 1331 { 1332 // skip fixed variables 1333 if(colstatus[i] == SPxBasisBase<R>::Desc::P_ON_LOWER 1334 || colstatus[i] == SPxBasisBase<R>::Desc::P_ON_UPPER) 1335 { 1336 // only consider continuous variables with zero dual multiplier to preserve optimality 1337 if(EQrel(this->maxObj(i) - (*thePvec)[i], (R) 0, this->epsilon()) && integerVariables[i] == 0) 1338 continuousvars.addIdx(i); 1339 } 1340 } 1341 } 1342 1343 while(!stop) 1344 { 1345 nSuccessfulPivots = 0; 1346 1347 // identify nonbasic slack variables, i.e. rows, that may be moved into the basis 1348 for(int i = slackcandidates.size() - 1; i >= 0 && !stop; --i) 1349 { 1350 polishId = coId(slackcandidates.index(i)); 1351 SPxOut::debug(this, "try pivoting: {} stat: {}\n", polishId, rowstatus[slackcandidates.index(i)]); 1352 success = enter(polishId, true); 1353 clearUpdateVecs(); 1354 1355 if(success) 1356 { 1357 SPxOut::debug(this, " -> success!"); 1358 ++nSuccessfulPivots; 1359 slackcandidates.remove(i); 1360 1361 if(maxIters >= 0 && iterations() >= maxIters) 1362 stop = true; 1363 } 1364 1365 SPxOut::debug(this, "\n"); 1366 1367 if(isTimeLimitReached()) 1368 stop = true; 1369 } 1370 1371 // identify nonbasic variables that may be moved into the basis 1372 for(int i = continuousvars.size() - 1; i >= 0 && !stop; --i) 1373 { 1374 polishId = id(continuousvars.index(i)); 1375 SPxOut::debug(this, "try pivoting: {} stat: {}\n", polishId, colstatus[continuousvars.index(i)]); 1376 success = enter(polishId, true); 1377 clearUpdateVecs(); 1378 1379 if(success) 1380 { 1381 SPxOut::debug(this, " -> success!"); 1382 ++nSuccessfulPivots; 1383 continuousvars.remove(i); 1384 1385 if(maxIters >= 0 && iterations() >= maxIters) 1386 stop = true; 1387 } 1388 1389 SPxOut::debug(this, "\n"); 1390 1391 if(isTimeLimitReached()) 1392 stop = true; 1393 } 1394 1395 // terminate if in the last round no more polishing steps were performed 1396 if(nSuccessfulPivots == 0) 1397 stop = true; 1398 1399 polishCount += nSuccessfulPivots; 1400 } 1401 } 1402 else 1403 { 1404 assert(polishObj == POLISH_FRACTIONALITY); 1405 DIdxSet candidates(dim()); 1406 1407 // identify nonbasic variables, i.e. columns, that may be moved into the basis 1408 for(int i = 0; i < this->nCols() && !stop; ++i) 1409 { 1410 if(colstatus[i] == SPxBasisBase<R>::Desc::P_ON_LOWER 1411 || colstatus[i] == SPxBasisBase<R>::Desc::P_ON_UPPER) 1412 { 1413 // only consider variables with zero reduced costs to preserve optimality 1414 if(EQrel(this->maxObj(i) - (*thePvec)[i], (R) 0, this->epsilon())) 1415 candidates.addIdx(i); 1416 } 1417 } 1418 1419 while(!stop) 1420 { 1421 nSuccessfulPivots = 0; 1422 1423 for(int i = candidates.size() - 1; i >= 0 && !stop; --i) 1424 { 1425 polishId = id(candidates.index(i)); 1426 SPxOut::debug(this, "try pivoting: {} stat: {}\n", polishId, colstatus[candidates.index(i)]); 1427 success = enter(polishId, true); 1428 clearUpdateVecs(); 1429 1430 if(success) 1431 { 1432 SPxOut::debug(this, " -> success!"); 1433 ++nSuccessfulPivots; 1434 candidates.remove(i); 1435 1436 if(maxIters >= 0 && iterations() >= maxIters) 1437 stop = true; 1438 } 1439 1440 SPxOut::debug(this, "\n"); 1441 1442 if(isTimeLimitReached()) 1443 stop = true; 1444 } 1445 1446 // terminate if in the last round no more polishing steps were performed 1447 if(nSuccessfulPivots == 0) 1448 stop = true; 1449 1450 polishCount += nSuccessfulPivots; 1451 } 1452 } 1453 } 1454 else 1455 { 1456 setType(LEAVE); // use primal simplex to preserve feasibility 1457 init(); 1458 alloweddeviation = leavetol(); 1459 instableLeave = false; 1460 theratiotester->setType(type()); 1461 bool useIntegrality = false; 1462 int ncols = this->nCols(); 1463 1464 if(integerVariables.size() == ncols) 1465 useIntegrality = true; 1466 1467 // in ROW rep: pivot slack out of the basis 1468 if(polishObj == POLISH_INTEGRALITY) 1469 { 1470 DIdxSet basiccandidates(dim()); 1471 1472 // collect basic candidates that may be moved out of the basis 1473 for(int i = 0; i < dim(); ++i) 1474 { 1475 polishId = this->baseId(i); 1476 1477 if(polishId.isSPxRowId()) 1478 stat = ds.rowStatus(this->number(polishId)); 1479 else 1480 { 1481 // skip (integer) variables 1482 if(!useIntegrality || integerVariables[this->number(SPxColId(polishId))] == 1) 1483 continue; 1484 1485 stat = ds.colStatus(this->number(polishId)); 1486 } 1487 1488 if(stat == SPxBasisBase<R>::Desc::P_ON_LOWER || stat == SPxBasisBase<R>::Desc::P_ON_UPPER) 1489 { 1490 if(EQrel((*theFvec)[i], (R) 0, this->epsilon())) 1491 basiccandidates.addIdx(i); 1492 } 1493 } 1494 1495 while(!stop) 1496 { 1497 nSuccessfulPivots = 0; 1498 1499 for(int i = basiccandidates.size() - 1; i >= 0 && !stop; --i) 1500 { 1501 1502 SPxOut::debug(this, "try pivoting: {}", this->baseId(basiccandidates.index(i))); 1503 success = leave(basiccandidates.index(i), true); 1504 clearUpdateVecs(); 1505 1506 if(success) 1507 { 1508 SPxOut::debug(this, " -> success!"); 1509 ++nSuccessfulPivots; 1510 basiccandidates.remove(i); 1511 1512 if(maxIters >= 0 && iterations() >= maxIters) 1513 stop = true; 1514 } 1515 1516 SPxOut::debug(this, "\n"); 1517 1518 if(isTimeLimitReached()) 1519 stop = true; 1520 } 1521 1522 // terminate if in the last round no more polishing steps were performed 1523 if(nSuccessfulPivots == 0) 1524 stop = true; 1525 1526 polishCount += nSuccessfulPivots; 1527 } 1528 } 1529 else 1530 { 1531 assert(polishObj == POLISH_FRACTIONALITY); 1532 DIdxSet basiccandidates(dim()); 1533 1534 // collect basic (integer) variables, that may be moved out of the basis 1535 for(int i = 0; i < dim(); ++i) 1536 { 1537 polishId = this->baseId(i); 1538 1539 if(polishId.isSPxRowId()) 1540 continue; 1541 else 1542 { 1543 if(useIntegrality && integerVariables[this->number(SPxColId(polishId))] == 0) 1544 continue; 1545 1546 stat = ds.colStatus(i); 1547 } 1548 1549 if(stat == SPxBasisBase<R>::Desc::P_ON_LOWER || stat == SPxBasisBase<R>::Desc::P_ON_UPPER) 1550 { 1551 if(EQrel((*theFvec)[i], (R) 0, this->epsilon())) 1552 basiccandidates.addIdx(i); 1553 } 1554 } 1555 1556 while(!stop) 1557 { 1558 nSuccessfulPivots = 0; 1559 1560 for(int i = basiccandidates.size() - 1; i >= 0 && !stop; --i) 1561 { 1562 SPxOut::debug(this, "try pivoting: {}", this->baseId(basiccandidates.index(i))); 1563 success = leave(basiccandidates.index(i), true); 1564 clearUpdateVecs(); 1565 1566 if(success) 1567 { 1568 SPxOut::debug(this, " -> success!"); 1569 ++nSuccessfulPivots; 1570 basiccandidates.remove(i); 1571 1572 if(maxIters >= 0 && iterations() >= maxIters) 1573 stop = true; 1574 } 1575 1576 SPxOut::debug(this, "\n"); 1577 1578 if(isTimeLimitReached()) 1579 stop = true; 1580 } 1581 1582 // terminate if in the last round no more polishing steps were performed 1583 if(nSuccessfulPivots == 0) 1584 stop = true; 1585 1586 polishCount += nSuccessfulPivots; 1587 } 1588 } 1589 } 1590 1591 SPX_MSG_INFO1((*this->spxout), 1592 (*this->spxout) << " --- finished solution polishing (" << polishCount << " pivots)" << std::endl;) 1593 1594 this->setStatus(SPxBasisBase<R>::OPTIMAL); 1595 1596 // if the value() changed significantly (due to numerics) reoptimize after polishing 1597 if(!EQrel(value(), origval, alloweddeviation)) 1598 return true; 1599 else 1600 return false; 1601 } 1602 1603 1604 template <class R> 1605 void SPxSolverBase<R>::testVecs() 1606 { 1607 1608 assert(SPxBasisBase<R>::status() > SPxBasisBase<R>::SINGULAR); 1609 1610 VectorBase<R> tmp(dim()); 1611 1612 tmp = *theCoPvec; 1613 this->multWithBase(tmp); 1614 tmp -= *theCoPrhs; 1615 1616 if(tmp.length() > leavetol()) 1617 { 1618 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "ISOLVE93 " << this->iteration() << 1619 ":\tcoP error = \t" 1620 << tmp.length() << std::endl;) 1621 1622 tmp.clear(); 1623 SPxBasisBase<R>::coSolve(tmp, *theCoPrhs); 1624 this->multWithBase(tmp); 1625 tmp -= *theCoPrhs; 1626 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "ISOLVE94\t\t" << tmp.length() << std::endl;) 1627 1628 tmp.clear(); 1629 SPxBasisBase<R>::coSolve(tmp, *theCoPrhs); 1630 tmp -= *theCoPvec; 1631 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "ISOLVE95\t\t" << tmp.length() << std::endl;) 1632 } 1633 1634 tmp = *theFvec; 1635 this->multBaseWith(tmp); 1636 tmp -= *theFrhs; 1637 1638 if(tmp.length() > entertol()) 1639 { 1640 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "ISOLVE96 " << this->iteration() << 1641 ":\t F error = \t" 1642 << tmp.length() << std::endl;) 1643 1644 tmp.clear(); 1645 SPxBasisBase<R>::solve(tmp, *theFrhs); 1646 tmp -= *theFvec; 1647 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "ISOLVE97\t\t" << tmp.length() << std::endl;) 1648 } 1649 1650 if(type() == ENTER) 1651 { 1652 for(int i = 0; i < dim(); ++i) 1653 { 1654 if(theCoTest[i] < -leavetol() && isCoBasic(i)) 1655 { 1656 /// @todo Error message "this shalt not be": shalt this be an assert (also below)? 1657 SPX_MSG_INFO1((*this->spxout), (*this->spxout) << "ESOLVE98 testVecs: theCoTest: this shalt not be!" 1658 << std::endl 1659 << " i=" << i 1660 << ", theCoTest[i]=" << theCoTest[i] 1661 << ", leavetol()=" << leavetol() << std::endl;) 1662 } 1663 } 1664 1665 for(int i = 0; i < coDim(); ++i) 1666 { 1667 if(theTest[i] < -leavetol() && isBasic(i)) 1668 { 1669 SPX_MSG_INFO1((*this->spxout), (*this->spxout) << "ESOLVE99 testVecs: theTest: this shalt not be!" 1670 << std::endl 1671 << " i=" << i 1672 << ", theTest[i]=" << theTest[i] 1673 << ", leavetol()=" << leavetol() << std::endl;) 1674 } 1675 } 1676 } 1677 } 1678 1679 1680 /// print display line of flying table 1681 template <class R> 1682 void SPxSolverBase<R>::printDisplayLine(const bool force, const bool forceHead) 1683 { 1684 SPX_MSG_INFO1((*this->spxout), 1685 1686 if(forceHead || displayLine % (displayFreq * 30) == 0) 1687 { 1688 (*this->spxout) 1689 << "type | time | iters | facts | shift | viol sum | viol num | obj value "; 1690 1691 if(printBasisMetric >= 0) 1692 (*this->spxout) << " | basis metric"; 1693 1694 (*this->spxout) << std::endl; 1695 } 1696 if((force || (displayLine % displayFreq == 0)) && !forceHead) 1697 { 1698 (type() == LEAVE) 1699 ? (*this->spxout) << " L |" : (*this->spxout) << " E |"; 1700 (*this->spxout) << std::fixed << std::setw(7) << std::setprecision(1) << time() << " |"; 1701 (*this->spxout) << std::scientific << std::setprecision(2); 1702 (*this->spxout) << std::setw(8) << this->iteration() << " | " 1703 << std::setw(5) << slinSolver()->getFactorCount() << " | " 1704 << shift() << " | " 1705 << SOPLEX_MAX(0.0, m_pricingViol + m_pricingViolCo) << " | " 1706 << std::setw(8) << SOPLEX_MAX(0, m_numViol) << " | " 1707 << std::setprecision(8) << value(); 1708 1709 if(getStartingDecompBasis && rep() == SPxSolverBase<R>::ROW) 1710 (*this->spxout) << " (" << std::fixed << std::setprecision(2) << getDegeneracyLevel(fVec()) << ")"; 1711 1712 if(printBasisMetric == 0) 1713 (*this->spxout) << " | " << std::scientific << std::setprecision(2) << getBasisMetric(0); 1714 1715 if(printBasisMetric == 1) 1716 (*this->spxout) << " | " << std::scientific << std::setprecision(2) << getBasisMetric(1); 1717 1718 if(printBasisMetric == 2) 1719 (*this->spxout) << " | " << std::scientific << std::setprecision(2) << getBasisMetric(2); 1720 1721 if(printBasisMetric == 3) 1722 (*this->spxout) << " | " << std::scientific << std::setprecision(2) << 1723 basis().getEstimatedCondition(); 1724 1725 (*this->spxout) << std::endl; 1726 } 1727 displayLine++; 1728 ); 1729 } 1730 1731 1732 template <class R> 1733 bool SPxSolverBase<R>::terminate() 1734 { 1735 #ifdef ENABLE_ADDITIONAL_CHECKS 1736 1737 if(SPxBasisBase<R>::status() > SPxBasisBase<R>::SINGULAR) 1738 testVecs(); 1739 1740 #endif 1741 1742 int redo = dim(); 1743 1744 if(redo < 1000) 1745 redo = 1000; 1746 1747 if(this->iteration() > 10 && this->iteration() % redo == 0) 1748 { 1749 #ifdef ENABLE_ADDITIONAL_CHECKS 1750 VectorBase<R> cr(*theCoPrhs); 1751 VectorBase<R> fr(*theFrhs); 1752 #endif 1753 1754 if(type() == ENTER) 1755 computeEnterCoPrhs(); 1756 else 1757 computeLeaveCoPrhs(); 1758 1759 computeFrhs(); 1760 1761 #ifdef ENABLE_ADDITIONAL_CHECKS 1762 cr -= *theCoPrhs; 1763 fr -= *theFrhs; 1764 1765 if(cr.length() > leavetol()) 1766 SPX_MSG_WARNING((*this->spxout), (*this->spxout) << "WSOLVE50 unexpected change of coPrhs " 1767 << cr.length() << std::endl;) 1768 if(fr.length() > entertol()) 1769 SPX_MSG_WARNING((*this->spxout), (*this->spxout) << "WSOLVE51 unexpected change of Frhs " 1770 << fr.length() << std::endl;) 1771 #endif 1772 1773 if(this->updateCount > 1) 1774 { 1775 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << " --- terminate triggers refactorization" 1776 << std::endl;) 1777 factorize(); 1778 } 1779 1780 SPxBasisBase<R>::coSolve(*theCoPvec, *theCoPrhs); 1781 SPxBasisBase<R>::solve(*theFvec, *theFrhs); 1782 1783 if(pricing() == FULL) 1784 { 1785 computePvec(); 1786 1787 if(type() == ENTER) 1788 { 1789 computeCoTest(); 1790 computeTest(); 1791 } 1792 } 1793 1794 if(shift() > 0.0) 1795 unShift(); 1796 } 1797 1798 // check time limit and objective limit only for non-terminal bases 1799 if(SPxBasisBase<R>::status() >= SPxBasisBase<R>::OPTIMAL || 1800 SPxBasisBase<R>::status() <= SPxBasisBase<R>::SINGULAR) 1801 { 1802 m_status = UNKNOWN; 1803 return true; 1804 } 1805 1806 if(isTimeLimitReached()) 1807 { 1808 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << " --- timelimit (" << maxTime 1809 << ") reached" << std::endl;) 1810 m_status = ABORT_TIME; 1811 return true; 1812 } 1813 1814 // objLimit is set and we are running DUAL: 1815 // - objLimit is set if objLimit < R(infinity) 1816 // - DUAL is running if rep() * type() > 0 == DUAL (-1 == PRIMAL) 1817 // 1818 // In this case we have given a objective value limit, e.g, through a 1819 // MIP solver, and we want stop solving the LP if we figure out that the 1820 // optimal value of the current LP can not be better then this objective 1821 // limit. More precisely: 1822 // - MINIMIZATION Problem 1823 // We want stop the solving process if 1824 // objLimit <= current objective value of the DUAL LP 1825 // - MAXIMIZATION Problem 1826 // We want stop the solving process if 1827 // objLimit >= current objective value of the DUAL LP 1828 if(objLimit < R(infinity) && type() * rep() > 0) 1829 { 1830 // We have no bound shifts; therefore, we can trust the current 1831 // objective value. 1832 // It might be even possible to use this termination value in case of 1833 // bound violations (shifting) but in this case it is quite difficult 1834 // to determine if we already reached the limit. 1835 if(shift() < epsilon() && noViols(tolerances()->floatingPointOpttol() - shift())) 1836 { 1837 // SPxSense::MINIMIZE == -1, so we have sign = 1 on minimizing 1838 if(int(this->spxSense()) * value() <= int(this->spxSense()) * objLimit) 1839 { 1840 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << " --- objective value limit (" << objLimit 1841 << ") reached" << std::endl;) 1842 SPxOut::debug(this, " --- objective value limit reached\n (value:{} limit:{})\n", value(), 1843 objLimit); 1844 SPxOut::debug(this, " (spxSense:{} rep:{} type:{})\n", int(this->spxSense()), int(rep()), 1845 int(type())); 1846 1847 m_status = ABORT_VALUE; 1848 return true; 1849 } 1850 } 1851 } 1852 1853 1854 1855 if(getComputeDegeneracy() && this->iteration() > this->prevIteration()) 1856 { 1857 VectorBase<R> degenvec(this->nCols()); 1858 1859 if(rep() == ROW) 1860 { 1861 if(type() == ENTER) // dual simplex 1862 dualDegenSum += getDegeneracyLevel(fVec()); 1863 else // primal simplex 1864 { 1865 getPrimalSol(degenvec); 1866 primalDegenSum += getDegeneracyLevel(degenvec); 1867 } 1868 } 1869 else 1870 { 1871 assert(rep() == COLUMN); 1872 1873 if(type() == LEAVE) // dual simplex 1874 dualDegenSum += getDegeneracyLevel(pVec()); 1875 else 1876 { 1877 getPrimalSol(degenvec); 1878 primalDegenSum += getDegeneracyLevel(degenvec); 1879 } 1880 } 1881 } 1882 1883 1884 // the improved dual simplex requires a starting basis 1885 // if the flag getStartingDecompBasis is set to true the simplex will terminate when a dual basis is found 1886 if(getStartingDecompBasis) 1887 { 1888 R iterationFrac = 0.6; 1889 1890 if(type() == ENTER && SPxBasisBase<R>::status() == SPxBasisBase<R>::DUAL && 1891 this->iteration() - this->lastDegenCheck() > getDegenCompOffset()/*iteration() % 10 == 0*/) 1892 { 1893 this->iterDegenCheck = this->iterCount; 1894 1895 if(SPxBasisBase<R>::status() >= SPxBasisBase<R>::OPTIMAL) 1896 { 1897 m_status = RUNNING; 1898 return true; 1899 } 1900 1901 R degeneracyLevel = 0; 1902 R degeneracyLB = 0.1; 1903 R degeneracyUB = 0.9; 1904 degeneracyLevel = getDegeneracyLevel(fVec()); 1905 1906 if((degeneracyLevel < degeneracyUB && degeneracyLevel > degeneracyLB) 1907 && this->iteration() > this->nRows() * 0.2) 1908 { 1909 m_status = ABORT_DECOMP; 1910 return true; 1911 } 1912 1913 if(degeneracyLevel < degeneracyLB 1914 && this->iteration() > SOPLEX_MIN(getDecompIterationLimit(), int(this->nCols()*iterationFrac))) 1915 { 1916 setDecompIterationLimit(0); 1917 setDegenCompOffset(0); 1918 m_status = ABORT_EXDECOMP; 1919 return true; 1920 } 1921 } 1922 else if(type() == LEAVE 1923 && this->iteration() > SOPLEX_MIN(getDecompIterationLimit(), int(this->nCols()*iterationFrac))) 1924 { 1925 setDecompIterationLimit(0); 1926 setDegenCompOffset(0); 1927 m_status = ABORT_EXDECOMP; 1928 return true; 1929 } 1930 } 1931 1932 this->lastIterCount = this->iterCount; 1933 1934 return false; 1935 } 1936 1937 template <class R> 1938 typename SPxSolverBase<R>::Status SPxSolverBase<R>::getPrimalSol(VectorBase<R>& p_vector) const 1939 { 1940 1941 if(!isInitialized()) 1942 { 1943 /* exit if presolving/simplifier cleared the problem */ 1944 if(status() == NO_PROBLEM) 1945 return status(); 1946 1947 throw SPxStatusException("XSOLVE06 Not Initialized"); 1948 } 1949 1950 if(rep() == ROW) 1951 p_vector = coPvec(); 1952 else 1953 { 1954 const typename SPxBasisBase<R>::Desc& ds = this->desc(); 1955 1956 for(int i = 0; i < this->nCols(); ++i) 1957 { 1958 switch(ds.colStatus(i)) 1959 { 1960 case SPxBasisBase<R>::Desc::P_ON_LOWER : 1961 p_vector[i] = SPxLPBase<R>::lower(i); 1962 break; 1963 1964 case SPxBasisBase<R>::Desc::P_ON_UPPER : 1965 case SPxBasisBase<R>::Desc::P_FIXED : 1966 p_vector[i] = SPxLPBase<R>::upper(i); 1967 break; 1968 1969 case SPxBasisBase<R>::Desc::P_FREE : 1970 p_vector[i] = 0; 1971 break; 1972 1973 case SPxBasisBase<R>::Desc::D_FREE : 1974 case SPxBasisBase<R>::Desc::D_ON_UPPER : 1975 case SPxBasisBase<R>::Desc::D_ON_LOWER : 1976 case SPxBasisBase<R>::Desc::D_ON_BOTH : 1977 case SPxBasisBase<R>::Desc::D_UNDEFINED : 1978 break; 1979 1980 default: 1981 throw SPxInternalCodeException("XSOLVE07 This should never happen."); 1982 } 1983 } 1984 1985 for(int j = 0; j < dim(); ++j) 1986 { 1987 if(this->baseId(j).isSPxColId()) 1988 p_vector[ this->number(SPxColId(this->baseId(j))) ] = fVec()[j]; 1989 } 1990 } 1991 1992 return status(); 1993 } 1994 1995 template <class R> 1996 typename SPxSolverBase<R>::Status SPxSolverBase<R>::getDualSol(VectorBase<R>& p_vector) const 1997 { 1998 1999 assert(isInitialized()); 2000 2001 if(!isInitialized()) 2002 { 2003 /* exit if presolving/simplifier cleared the problem */ 2004 if(status() == NO_PROBLEM) 2005 return status(); 2006 2007 throw SPxStatusException("XSOLVE08 No Problem loaded"); 2008 } 2009 2010 if(rep() == ROW) 2011 { 2012 int i; 2013 p_vector = this->maxRowObj(); 2014 2015 for(i = this->nCols() - 1; i >= 0; --i) 2016 { 2017 if(this->baseId(i).isSPxRowId()) 2018 p_vector[ this->number(SPxRowId(this->baseId(i))) ] = fVec()[i]; 2019 } 2020 } 2021 else 2022 { 2023 const typename SPxBasisBase<R>::Desc& ds = this->desc(); 2024 2025 for(int i = 0; i < this->nRows(); ++i) 2026 { 2027 switch(ds.rowStatus(i)) 2028 { 2029 case SPxBasisBase<R>::Desc::D_FREE: 2030 case SPxBasisBase<R>::Desc::D_ON_UPPER: 2031 case SPxBasisBase<R>::Desc::D_ON_LOWER: 2032 case SPxBasisBase<R>::Desc::D_ON_BOTH: 2033 case SPxBasisBase<R>::Desc::D_UNDEFINED: 2034 p_vector[i] = 0; 2035 break; 2036 2037 default: 2038 p_vector[i] = (*theCoPvec)[i]; 2039 } 2040 } 2041 } 2042 2043 p_vector *= Real(this->spxSense()); 2044 2045 return status(); 2046 } 2047 2048 template <class R> 2049 typename SPxSolverBase<R>::Status SPxSolverBase<R>::getRedCostSol(VectorBase<R>& p_vector) const 2050 { 2051 2052 assert(isInitialized()); 2053 2054 if(!isInitialized()) 2055 { 2056 throw SPxStatusException("XSOLVE09 No Problem loaded"); 2057 // return NOT_INIT; 2058 } 2059 2060 if(rep() == ROW) 2061 { 2062 int i; 2063 p_vector.clear(); 2064 2065 if(this->spxSense() == SPxLPBase<R>::MINIMIZE) 2066 { 2067 for(i = dim() - 1; i >= 0; --i) 2068 { 2069 if(this->baseId(i).isSPxColId()) 2070 p_vector[ this->number(SPxColId(this->baseId(i))) ] = -fVec()[i]; 2071 } 2072 } 2073 else 2074 { 2075 for(i = dim() - 1; i >= 0; --i) 2076 { 2077 if(this->baseId(i).isSPxColId()) 2078 p_vector[ this->number(SPxColId(this->baseId(i))) ] = fVec()[i]; 2079 } 2080 } 2081 } 2082 else 2083 { 2084 const typename SPxBasisBase<R>::Desc& ds = this->desc(); 2085 2086 for(int i = 0; i < this->nCols(); ++i) 2087 { 2088 switch(ds.colStatus(i)) 2089 { 2090 case SPxBasisBase<R>::Desc::D_FREE: 2091 case SPxBasisBase<R>::Desc::D_ON_UPPER: 2092 case SPxBasisBase<R>::Desc::D_ON_LOWER: 2093 case SPxBasisBase<R>::Desc::D_ON_BOTH: 2094 case SPxBasisBase<R>::Desc::D_UNDEFINED: 2095 p_vector[i] = 0; 2096 break; 2097 2098 default: 2099 p_vector[i] = this->maxObj()[i] - (*thePvec)[i]; 2100 } 2101 } 2102 2103 if(this->spxSense() == SPxLPBase<R>::MINIMIZE) 2104 p_vector *= -1.0; 2105 } 2106 2107 return status(); 2108 } 2109 2110 template <class R> 2111 typename SPxSolverBase<R>::Status SPxSolverBase<R>::getPrimalray(VectorBase<R>& p_vector) const 2112 { 2113 2114 assert(isInitialized()); 2115 2116 if(!isInitialized()) 2117 { 2118 throw SPxStatusException("XSOLVE10 No Problem loaded"); 2119 // return NOT_INIT; 2120 } 2121 2122 assert(SPxBasisBase<R>::status() == SPxBasisBase<R>::UNBOUNDED); 2123 p_vector.clear(); 2124 p_vector = primalRay; 2125 2126 return status(); 2127 } 2128 2129 template <class R> 2130 typename SPxSolverBase<R>::Status SPxSolverBase<R>::getDualfarkas(VectorBase<R>& p_vector) const 2131 { 2132 2133 assert(isInitialized()); 2134 2135 if(!isInitialized()) 2136 { 2137 throw SPxStatusException("XSOLVE10 No Problem loaded"); 2138 // return NOT_INIT; 2139 } 2140 2141 assert(SPxBasisBase<R>::status() == SPxBasisBase<R>::INFEASIBLE); 2142 p_vector.clear(); 2143 p_vector = dualFarkas; 2144 2145 return status(); 2146 } 2147 2148 template <class R> 2149 typename SPxSolverBase<R>::Status SPxSolverBase<R>::getSlacks(VectorBase<R>& p_vector) const 2150 { 2151 2152 assert(isInitialized()); 2153 2154 if(!isInitialized()) 2155 { 2156 throw SPxStatusException("XSOLVE11 No Problem loaded"); 2157 // return NOT_INIT; 2158 } 2159 2160 if(rep() == COLUMN) 2161 { 2162 int i; 2163 const typename SPxBasisBase<R>::Desc& ds = this->desc(); 2164 2165 for(i = this->nRows() - 1; i >= 0; --i) 2166 { 2167 switch(ds.rowStatus(i)) 2168 { 2169 case SPxBasisBase<R>::Desc::P_ON_LOWER : 2170 p_vector[i] = this->lhs(i); 2171 break; 2172 2173 case SPxBasisBase<R>::Desc::P_ON_UPPER : 2174 case SPxBasisBase<R>::Desc::P_FIXED : 2175 p_vector[i] = this->rhs(i); 2176 break; 2177 2178 case SPxBasisBase<R>::Desc::P_FREE : 2179 p_vector[i] = 0; 2180 break; 2181 2182 case SPxBasisBase<R>::Desc::D_FREE : 2183 case SPxBasisBase<R>::Desc::D_ON_UPPER : 2184 case SPxBasisBase<R>::Desc::D_ON_LOWER : 2185 case SPxBasisBase<R>::Desc::D_ON_BOTH : 2186 case SPxBasisBase<R>::Desc::D_UNDEFINED : 2187 break; 2188 2189 default: 2190 throw SPxInternalCodeException("XSOLVE12 This should never happen."); 2191 } 2192 } 2193 2194 for(i = dim() - 1; i >= 0; --i) 2195 { 2196 if(this->baseId(i).isSPxRowId()) 2197 p_vector[ this->number(SPxRowId(this->baseId(i))) ] = -(*theFvec)[i]; 2198 } 2199 } 2200 else 2201 p_vector = pVec(); 2202 2203 return status(); 2204 } 2205 2206 template <class R> 2207 void SPxSolverBase<R>::setPrimal(VectorBase<R>& p_vector) 2208 { 2209 2210 if(!isInitialized()) 2211 { 2212 throw SPxStatusException("XSOLVE20 Not Initialized"); 2213 } 2214 2215 if(rep() == ROW) 2216 coPvec() = p_vector; 2217 else 2218 { 2219 for(int j = 0; j < dim(); ++j) 2220 { 2221 if(this->baseId(j).isSPxColId()) 2222 fVec()[j] = p_vector[ this->number(SPxColId(this->baseId(j))) ]; 2223 } 2224 } 2225 } 2226 2227 template <class R> 2228 void SPxSolverBase<R>::setDual(VectorBase<R>& p_vector) 2229 { 2230 2231 assert(isInitialized()); 2232 2233 if(!isInitialized()) 2234 { 2235 throw SPxStatusException("XSOLVE21 Not Initialized"); 2236 } 2237 2238 if(rep() == ROW) 2239 { 2240 for(int i = this->nCols() - 1; i >= 0; --i) 2241 { 2242 if(this->baseId(i).isSPxRowId()) 2243 { 2244 if(this->spxSense() == SPxLPBase<R>::MAXIMIZE) 2245 fVec()[i] = p_vector[ this->number(SPxRowId(this->baseId(i))) ]; 2246 else 2247 fVec()[i] = -p_vector[ this->number(SPxRowId(this->baseId(i))) ]; 2248 } 2249 } 2250 } 2251 else 2252 { 2253 coPvec() = p_vector; 2254 2255 if(this->spxSense() == SPxLPBase<R>::MINIMIZE) 2256 coPvec() *= -1.0; 2257 } 2258 } 2259 2260 template <class R> 2261 void SPxSolverBase<R>::setRedCost(VectorBase<R>& p_vector) 2262 { 2263 2264 assert(isInitialized()); 2265 2266 if(!isInitialized()) 2267 { 2268 throw SPxStatusException("XSOLVE22 Not Initialized"); 2269 } 2270 2271 if(rep() == ROW) 2272 { 2273 for(int i = dim() - 1; i >= 0; --i) 2274 { 2275 if(this->baseId(i).isSPxColId()) 2276 { 2277 if(this->spxSense() == SPxLPBase<R>::MINIMIZE) 2278 fVec()[i] = -p_vector[ this->number(SPxColId(this->baseId(i))) ]; 2279 else 2280 fVec()[i] = p_vector[ this->number(SPxColId(this->baseId(i))) ]; 2281 } 2282 } 2283 } 2284 else 2285 { 2286 pVec() = this->maxObj(); 2287 2288 if(this->spxSense() == SPxLPBase<R>::MINIMIZE) 2289 pVec() += p_vector; 2290 else 2291 pVec() -= p_vector; 2292 } 2293 } 2294 2295 template <class R> 2296 void SPxSolverBase<R>::setSlacks(VectorBase<R>& p_vector) 2297 { 2298 2299 assert(isInitialized()); 2300 2301 if(!isInitialized()) 2302 { 2303 throw SPxStatusException("XSOLVE23 Not Initialized"); 2304 } 2305 2306 if(rep() == COLUMN) 2307 { 2308 for(int i = dim() - 1; i >= 0; --i) 2309 { 2310 if(this->baseId(i).isSPxRowId()) 2311 (*theFvec)[i] = -p_vector[ this->number(SPxRowId(this->baseId(i))) ]; 2312 } 2313 } 2314 else 2315 pVec() = p_vector; 2316 } 2317 2318 template <class R> 2319 typename SPxSolverBase<R>::Status SPxSolverBase<R>::status() const 2320 { 2321 switch(m_status) 2322 { 2323 case UNKNOWN : 2324 switch(SPxBasisBase<R>::status()) 2325 { 2326 case SPxBasisBase<R>::NO_PROBLEM : 2327 return NO_PROBLEM; 2328 2329 case SPxBasisBase<R>::SINGULAR : 2330 return SINGULAR; 2331 2332 case SPxBasisBase<R>::REGULAR : 2333 case SPxBasisBase<R>::DUAL : 2334 case SPxBasisBase<R>::PRIMAL : 2335 return UNKNOWN; 2336 2337 case SPxBasisBase<R>::OPTIMAL : 2338 return OPTIMAL; 2339 2340 case SPxBasisBase<R>::UNBOUNDED : 2341 return UNBOUNDED; 2342 2343 case SPxBasisBase<R>::INFEASIBLE : 2344 return INFEASIBLE; 2345 2346 default: 2347 return ERROR; 2348 } 2349 2350 case SINGULAR : 2351 return m_status; 2352 2353 case OPTIMAL : 2354 assert(SPxBasisBase<R>::status() == SPxBasisBase<R>::OPTIMAL); 2355 2356 /*lint -fallthrough*/ 2357 case ABORT_EXDECOMP : 2358 case ABORT_DECOMP : 2359 case ABORT_CYCLING : 2360 case ABORT_TIME : 2361 case ABORT_ITER : 2362 case ABORT_VALUE : 2363 case RUNNING : 2364 case REGULAR : 2365 case NOT_INIT : 2366 case NO_SOLVER : 2367 case NO_PRICER : 2368 case NO_RATIOTESTER : 2369 case ERROR: 2370 return m_status; 2371 2372 default: 2373 return ERROR; 2374 } 2375 } 2376 2377 template <class R> 2378 typename SPxSolverBase<R>::Status SPxSolverBase<R>::getResult( 2379 R* p_value, 2380 VectorBase<R>* p_primal, 2381 VectorBase<R>* p_slacks, 2382 VectorBase<R>* p_dual, 2383 VectorBase<R>* reduCosts) 2384 { 2385 if(p_value) 2386 *p_value = this->value(); 2387 2388 if(p_primal) 2389 getPrimalSol(*p_primal); 2390 2391 if(p_slacks) 2392 getSlacks(*p_slacks); 2393 2394 if(p_dual) 2395 getDualSol(*p_dual); 2396 2397 if(reduCosts) 2398 getRedCostSol(*reduCosts); 2399 2400 return status(); 2401 } 2402 } // namespace soplex 2403