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 /* \SubSection{Updating the Basis for Entering Variables} 26 */ 27 #include <assert.h> 28 29 #include "soplex/spxdefines.h" 30 #include "soplex/spxratiotester.h" 31 #include "soplex/spxpricer.h" 32 #include "soplex/spxout.h" 33 #include "soplex/exceptions.h" 34 35 namespace soplex 36 { 37 38 /* 39 In the entering simplex algorithms (i.e. iteratively a vector is selected to 40 \em enter the simplex basis as in the dual rowwise and primal columnwise case) 41 let $A$ denote the current basis, $x$ and entering vector and $f$ the 42 feasibility vector. For a feasible basis $l \le f \le u$ holds. 43 For the rowwise case $f$ is obtained by solving $f^T = c^T A^{-1}$, 44 wherease in columnwisecase $f = A^{-1} b$. 45 46 Let us further consider the rowwise case. Exchanging $x$ with the $i$-th 47 vector of $A$ yields 48 49 \begin{equation}\label{update.eq} 50 A^{(i)} = E_i A \hbox{, with } E_i = I + e_i (x^T A^{-1} - e_i^T). 51 \end{equation} 52 53 With $E_i^{-1} = I + e_i \frac{e_i^T - \delta^T}{\delta}$, 54 $\delta^T = x^T A^{-1}$ one gets the new feasibility vector 55 56 \begin{eqnarray*} 57 (f^{(i)})^T 58 &=& c^T (A^{(i)})^{-1} \\ 59 &=& c^T A^{-1} + c^T A^{-1} e_i \frac{e_i^T - \delta^T}{\delta_i} \\ 60 &=& f^T + \frac{f_i}{\delta_i} e_i^T - \frac{f_i}{\delta_i} \delta^T. \\ 61 \end{eqnarray*} 62 63 The selection of the leaving vector $i^*$ for the basis must ensure, that for 64 all $j \ne i^*$ $f^{(i^*)}_j$ remains within its bounds $l_j$ and $u_j$. 65 */ 66 67 68 /* 69 Testing all values of |pVec| against its bounds. If $i$, say, is violated 70 the violation is saved as negative value in |theTest[i]|. 71 */ 72 73 template <class R> 74 R SPxSolverBase<R>::test(int i, typename SPxBasisBase<R>::Desc::Status stat) const 75 { 76 assert(type() == ENTER); 77 assert(!isBasic(stat)); 78 79 R x; 80 81 switch(stat) 82 { 83 case SPxBasisBase<R>::Desc::D_FREE: 84 case SPxBasisBase<R>::Desc::D_ON_BOTH: 85 assert(rep() == ROW); 86 x = (*thePvec)[i] - this->lhs(i); 87 88 if(x < 0) 89 return x; 90 91 // no break: next is else case 92 //lint -fallthrough 93 case SPxBasisBase<R>::Desc::D_ON_LOWER: 94 assert(rep() == ROW); 95 return this->rhs(i) - (*thePvec)[i]; 96 97 case SPxBasisBase<R>::Desc::D_ON_UPPER: 98 assert(rep() == ROW); 99 return (*thePvec)[i] - this->lhs(i); 100 101 case SPxBasisBase<R>::Desc::P_ON_UPPER: 102 assert(rep() == COLUMN); 103 return this->maxObj(i) - (*thePvec)[i]; 104 105 case SPxBasisBase<R>::Desc::P_ON_LOWER: 106 assert(rep() == COLUMN); 107 return (*thePvec)[i] - this->maxObj(i); 108 109 case SPxBasisBase<R>::Desc::P_FREE : 110 x = this->maxObj(i) - (*thePvec)[i]; 111 return (x < 0) ? x : -x; 112 113 default: 114 return 0; 115 } 116 } 117 118 template <class R> 119 void SPxSolverBase<R>::computeTest() 120 { 121 122 const typename SPxBasisBase<R>::Desc& ds = this->desc(); 123 R pricingTol = leavetol(); 124 m_pricingViolCoUpToDate = true; 125 m_pricingViolCo = 0; 126 127 infeasibilitiesCo.clear(); 128 int sparsitythreshold = (int)(sparsePricingFactor * coDim()); 129 130 for(int i = 0; i < coDim(); ++i) 131 { 132 typename SPxBasisBase<R>::Desc::Status stat = ds.status(i); 133 134 if(isBasic(stat)) 135 { 136 theTest[i] = 0.0; 137 138 if(remainingRoundsEnterCo == 0) 139 isInfeasibleCo[i] = SPxPricer<R>::NOT_VIOLATED; 140 } 141 else 142 { 143 assert(!isBasic(stat)); 144 theTest[i] = test(i, stat); 145 146 if(remainingRoundsEnterCo == 0) 147 { 148 if(theTest[i] < -pricingTol) 149 { 150 assert(infeasibilitiesCo.size() < infeasibilitiesCo.max()); 151 m_pricingViolCo -= theTest[i]; 152 infeasibilitiesCo.addIdx(i); 153 isInfeasibleCo[i] = SPxPricer<R>::VIOLATED; 154 ++m_numViol; 155 } 156 else 157 isInfeasibleCo[i] = SPxPricer<R>::NOT_VIOLATED; 158 159 if(infeasibilitiesCo.size() > sparsitythreshold) 160 { 161 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << " --- using dense pricing" 162 << std::endl;) 163 remainingRoundsEnterCo = SOPLEX_DENSEROUNDS; 164 sparsePricingEnterCo = false; 165 infeasibilitiesCo.clear(); 166 } 167 } 168 else if(theTest[i] < -pricingTol) 169 { 170 m_pricingViolCo -= theTest[i]; 171 ++m_numViol; 172 } 173 } 174 } 175 176 if(infeasibilitiesCo.size() == 0 && !sparsePricingEnterCo) 177 --remainingRoundsEnterCo; 178 else if(infeasibilitiesCo.size() <= sparsitythreshold && !sparsePricingEnterCo) 179 { 180 SPX_MSG_INFO2((*this->spxout), 181 std::streamsize prec = spxout->precision(); 182 183 if(hyperPricingEnter) 184 (*this->spxout) << " --- using hypersparse pricing, "; 185 else 186 (*this->spxout) << " --- using sparse pricing, "; 187 (*this->spxout) << "sparsity: " 188 << std::setw(6) << std::fixed << std::setprecision(4) 189 << (R) infeasibilitiesCo.size() / coDim() 190 << std::scientific << std::setprecision(int(prec)) 191 << std::endl;) 192 sparsePricingEnterCo = true; 193 } 194 } 195 196 template <class R> 197 R SPxSolverBase<R>::computePvec(int i) 198 { 199 200 return (*thePvec)[i] = vector(i) * (*theCoPvec); 201 } 202 203 template <class R> 204 R SPxSolverBase<R>::computeTest(int i) 205 { 206 typename SPxBasisBase<R>::Desc::Status stat = this->desc().status(i); 207 208 if(isBasic(stat)) 209 return theTest[i] = 0; 210 else 211 return theTest[i] = test(i, stat); 212 } 213 214 /* 215 Testing all values of #coPvec# against its bounds. If $i$, say, is violated 216 the violation is saved as negative value in |theCoTest[i]|. 217 */ 218 template <class R> 219 R SPxSolverBase<R>::coTest(int i, typename SPxBasisBase<R>::Desc::Status stat) const 220 { 221 assert(type() == ENTER); 222 assert(!isBasic(stat)); 223 224 R x; 225 226 switch(stat) 227 { 228 case SPxBasisBase<R>::Desc::D_FREE: 229 case SPxBasisBase<R>::Desc::D_ON_BOTH : 230 assert(rep() == ROW); 231 x = (*theCoPvec)[i] - SPxLPBase<R>::lower(i); 232 233 if(x < 0) 234 return x; 235 236 // no break: next is else case 237 //lint -fallthrough 238 case SPxBasisBase<R>::Desc::D_ON_LOWER: 239 assert(rep() == ROW); 240 return SPxLPBase<R>::upper(i) - (*theCoPvec)[i]; 241 242 case SPxBasisBase<R>::Desc::D_ON_UPPER: 243 assert(rep() == ROW); 244 return (*theCoPvec)[i] - SPxLPBase<R>::lower(i); 245 246 case SPxBasisBase<R>::Desc::P_ON_UPPER: 247 assert(rep() == COLUMN); 248 return (*theCoPvec)[i] - this->maxRowObj(i); // slacks ! 249 250 case SPxBasisBase<R>::Desc::P_ON_LOWER: 251 assert(rep() == COLUMN); 252 return this->maxRowObj(i) - (*theCoPvec)[i]; // slacks ! 253 254 default: 255 return 0; 256 } 257 } 258 259 template <class R> 260 void SPxSolverBase<R>::computeCoTest() 261 { 262 int i; 263 R pricingTol = leavetol(); 264 m_pricingViolUpToDate = true; 265 m_pricingViol = 0; 266 m_numViol = 0; 267 infeasibilities.clear(); 268 int sparsitythreshold = (int)(sparsePricingFactor * dim()); 269 const typename SPxBasisBase<R>::Desc& ds = this->desc(); 270 271 for(i = dim() - 1; i >= 0; --i) 272 { 273 typename SPxBasisBase<R>::Desc::Status stat = ds.coStatus(i); 274 275 if(isBasic(stat)) 276 { 277 theCoTest[i] = 0; 278 279 if(remainingRoundsEnter == 0) 280 isInfeasible[i] = SPxPricer<R>::NOT_VIOLATED; 281 } 282 else 283 { 284 theCoTest[i] = coTest(i, stat); 285 286 if(remainingRoundsEnter == 0) 287 { 288 if(theCoTest[i] < -pricingTol) 289 { 290 assert(infeasibilities.size() < infeasibilities.max()); 291 m_pricingViol -= theCoTest[i]; 292 infeasibilities.addIdx(i); 293 isInfeasible[i] = SPxPricer<R>::VIOLATED; 294 ++m_numViol; 295 } 296 else 297 isInfeasible[i] = SPxPricer<R>::NOT_VIOLATED; 298 299 if(infeasibilities.size() > sparsitythreshold) 300 { 301 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << " --- using dense pricing" 302 << std::endl;) 303 remainingRoundsEnter = SOPLEX_DENSEROUNDS; 304 sparsePricingEnter = false; 305 infeasibilities.clear(); 306 } 307 } 308 else if(theCoTest[i] < -pricingTol) 309 { 310 m_pricingViol -= theCoTest[i]; 311 ++m_numViol; 312 } 313 } 314 } 315 316 if(infeasibilities.size() == 0 && !sparsePricingEnter) 317 --remainingRoundsEnter; 318 else if(infeasibilities.size() <= sparsitythreshold && !sparsePricingEnter) 319 { 320 SPX_MSG_INFO2((*this->spxout), 321 std::streamsize prec = spxout->precision(); 322 323 if(hyperPricingEnter) 324 (*this->spxout) << " --- using hypersparse pricing, "; 325 else 326 (*this->spxout) << " --- using sparse pricing, "; 327 (*this->spxout) << "sparsity: " 328 << std::setw(6) << std::fixed << std::setprecision(4) 329 << (R) infeasibilities.size() / dim() 330 << std::scientific << std::setprecision(int(prec)) 331 << std::endl;) 332 sparsePricingEnter = true; 333 } 334 } 335 336 /* 337 The following methods require propersy initialized vectors |fVec| and 338 #coPvec#. 339 */ 340 template <class R> 341 void SPxSolverBase<R>::updateTest() 342 { 343 thePvec->delta().setup(); 344 345 const IdxSet& idx = thePvec->idx(); 346 const typename SPxBasisBase<R>::Desc& ds = this->desc(); 347 R pricingTol = leavetol(); 348 349 int i; 350 updateViolsCo.clear(); 351 352 for(i = idx.size() - 1; i >= 0; --i) 353 { 354 int j = idx.index(i); 355 typename SPxBasisBase<R>::Desc::Status stat = ds.status(j); 356 357 if(!isBasic(stat)) 358 { 359 if(m_pricingViolCoUpToDate && theTest[j] < -pricingTol) 360 m_pricingViolCo += theTest[j]; 361 362 theTest[j] = test(j, stat); 363 364 if(sparsePricingEnterCo) 365 { 366 if(theTest[j] < -pricingTol) 367 { 368 assert(remainingRoundsEnterCo == 0); 369 m_pricingViolCo -= theTest[j]; 370 371 if(isInfeasibleCo[j] == SPxPricer<R>::NOT_VIOLATED) 372 { 373 infeasibilitiesCo.addIdx(j); 374 isInfeasibleCo[j] = SPxPricer<R>::VIOLATED; 375 } 376 377 if(hyperPricingEnter) 378 updateViolsCo.addIdx(j); 379 } 380 else 381 { 382 isInfeasibleCo[j] = SPxPricer<R>::NOT_VIOLATED; 383 } 384 } 385 else if(theTest[j] < -pricingTol) 386 m_pricingViolCo -= theTest[j]; 387 } 388 else 389 { 390 isInfeasibleCo[j] = SPxPricer<R>::NOT_VIOLATED; 391 theTest[j] = 0; 392 } 393 } 394 } 395 396 template <class R> 397 void SPxSolverBase<R>::updateCoTest() 398 { 399 theCoPvec->delta().setup(); 400 401 const IdxSet& idx = theCoPvec->idx(); 402 const typename SPxBasisBase<R>::Desc& ds = this->desc(); 403 R pricingTol = leavetol(); 404 405 int i; 406 updateViols.clear(); 407 408 for(i = idx.size() - 1; i >= 0; --i) 409 { 410 int j = idx.index(i); 411 typename SPxBasisBase<R>::Desc::Status stat = ds.coStatus(j); 412 413 if(!isBasic(stat)) 414 { 415 if(m_pricingViolUpToDate && theCoTest[j] < -pricingTol) 416 m_pricingViol += theCoTest[j]; 417 418 theCoTest[j] = coTest(j, stat); 419 420 if(sparsePricingEnter) 421 { 422 if(theCoTest[j] < -pricingTol) 423 { 424 assert(remainingRoundsEnter == 0); 425 m_pricingViol -= theCoTest[j]; 426 427 if(isInfeasible[j] == SPxPricer<R>::NOT_VIOLATED) 428 { 429 // if( !hyperPricingEnter ) 430 infeasibilities.addIdx(j); 431 isInfeasible[j] = SPxPricer<R>::VIOLATED; 432 } 433 434 if(hyperPricingEnter) 435 updateViols.addIdx(j); 436 } 437 else 438 { 439 // @todo do we need to remove index j from infeasibilitiesCo? 440 isInfeasible[j] = SPxPricer<R>::NOT_VIOLATED; 441 } 442 } 443 else if(theCoTest[j] < -pricingTol) 444 m_pricingViol -= theCoTest[j]; 445 } 446 else 447 { 448 isInfeasible[j] = SPxPricer<R>::NOT_VIOLATED; 449 theCoTest[j] = 0; 450 } 451 } 452 } 453 454 455 456 /* \Section{Compute statistics on entering variable} 457 Here is a list of variables relevant when including |Id| to the basis. 458 They are computed by |computeEnterStats()|. 459 */ 460 template <class R> 461 void SPxSolverBase<R>::getEnterVals 462 ( 463 SPxId enterId, 464 R& enterTest, 465 R& enterUB, 466 R& enterLB, 467 R& enterVal, 468 R& enterMax, 469 R& enterPric, 470 typename SPxBasisBase<R>::Desc::Status& enterStat, 471 R& enterRO, 472 StableSum<R>& objChange 473 ) 474 { 475 int enterIdx; 476 typename SPxBasisBase<R>::Desc& ds = this->desc(); 477 478 if(enterId.isSPxColId()) 479 { 480 enterIdx = this->number(SPxColId(enterId)); 481 enterStat = ds.colStatus(enterIdx); 482 assert(!isBasic(enterStat)); 483 484 /* For an #Id# to enter the basis we better recompute the Test value. 485 */ 486 if(rep() == COLUMN) 487 { 488 computePvec(enterIdx); 489 enterTest = computeTest(enterIdx); 490 theTest[enterIdx] = 0; 491 } 492 else 493 { 494 enterTest = coTest()[enterIdx]; 495 theCoTest[enterIdx] = 0; 496 } 497 498 switch(enterStat) 499 { 500 // primal/columnwise cases: 501 case SPxBasisBase<R>::Desc::P_ON_UPPER : 502 assert(rep() == COLUMN); 503 enterUB = theUCbound[enterIdx]; 504 enterLB = theLCbound[enterIdx]; 505 enterVal = enterUB; 506 enterMax = enterLB - enterUB; 507 enterPric = (*thePvec)[enterIdx]; 508 enterRO = this->maxObj(enterIdx); 509 objChange -= enterVal * enterRO; 510 511 if(enterLB <= R(-infinity)) 512 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::D_ON_LOWER; 513 else if(EQ(enterLB, enterUB, this->epsilon())) 514 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::D_FREE; 515 else 516 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::D_ON_BOTH; 517 518 break; 519 520 case SPxBasisBase<R>::Desc::P_ON_LOWER : 521 assert(rep() == COLUMN); 522 enterUB = theUCbound[enterIdx]; 523 enterLB = theLCbound[enterIdx]; 524 enterVal = enterLB; 525 enterMax = enterUB - enterLB; 526 enterPric = (*thePvec)[enterIdx]; 527 enterRO = this->maxObj(enterIdx); 528 objChange -= enterVal * enterRO; 529 530 if(enterUB >= R(infinity)) 531 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::D_ON_UPPER; 532 else if(EQ(enterLB, enterUB, this->epsilon())) 533 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::D_FREE; 534 else 535 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::D_ON_BOTH; 536 537 break; 538 539 case SPxBasisBase<R>::Desc::P_FREE : 540 assert(rep() == COLUMN); 541 enterUB = theUCbound[enterIdx]; 542 enterLB = theLCbound[enterIdx]; 543 enterVal = 0; 544 enterPric = (*thePvec)[enterIdx]; 545 enterRO = this->maxObj(enterIdx); 546 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::D_UNDEFINED; 547 enterMax = (enterRO - enterPric > 0) ? R(infinity) : R(-infinity); 548 break; 549 550 // dual/rowwise cases: 551 case SPxBasisBase<R>::Desc::D_ON_UPPER : 552 assert(rep() == ROW); 553 assert(theUCbound[enterIdx] < R(infinity)); 554 enterUB = theUCbound[enterIdx]; 555 enterLB = R(-infinity); 556 enterMax = R(-infinity); 557 enterVal = enterUB; 558 enterPric = (*theCoPvec)[enterIdx]; 559 enterRO = SPxLPBase<R>::lower(enterIdx); 560 objChange -= enterRO * enterVal; 561 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_LOWER; 562 break; 563 564 case SPxBasisBase<R>::Desc::D_ON_LOWER : 565 assert(rep() == ROW); 566 assert(theLCbound[enterIdx] > R(-infinity)); 567 enterLB = theLCbound[enterIdx]; 568 enterUB = R(infinity); 569 enterMax = R(infinity); 570 enterVal = enterLB; 571 enterPric = (*theCoPvec)[enterIdx]; 572 enterRO = SPxLPBase<R>::upper(enterIdx); 573 objChange -= enterRO * enterVal; 574 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_UPPER; 575 break; 576 577 case SPxBasisBase<R>::Desc::D_FREE: 578 assert(rep() == ROW); 579 assert(SPxLPBase<R>::lower(enterIdx) == SPxLPBase<R>::upper(enterIdx)); 580 enterUB = R(infinity); 581 enterLB = R(-infinity); 582 enterVal = 0; 583 enterRO = SPxLPBase<R>::upper(enterIdx); 584 enterPric = (*theCoPvec)[enterIdx]; 585 586 if(enterPric > enterRO) 587 enterMax = R(infinity); 588 else 589 enterMax = R(-infinity); 590 591 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::P_FIXED; 592 break; 593 594 case SPxBasisBase<R>::Desc::D_ON_BOTH : 595 assert(rep() == ROW); 596 enterPric = (*theCoPvec)[enterIdx]; 597 598 if(enterPric > SPxLPBase<R>::upper(enterIdx)) 599 { 600 enterLB = theLCbound[enterIdx]; 601 enterUB = R(infinity); 602 enterMax = R(infinity); 603 enterVal = enterLB; 604 enterRO = SPxLPBase<R>::upper(enterIdx); 605 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_UPPER; 606 } 607 else 608 { 609 enterUB = theUCbound[enterIdx]; 610 enterVal = enterUB; 611 enterRO = SPxLPBase<R>::lower(enterIdx); 612 enterLB = R(-infinity); 613 enterMax = R(-infinity); 614 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_LOWER; 615 } 616 617 objChange -= theLCbound[enterIdx] * SPxLPBase<R>::upper(enterIdx); 618 objChange -= theUCbound[enterIdx] * SPxLPBase<R>::lower(enterIdx); 619 break; 620 621 default: 622 throw SPxInternalCodeException("XENTER01 This should never happen."); 623 } 624 625 SPxOut::debug(this, "DENTER03 SPxSolverBase::getEnterVals() : col {}: {} -> {} objChange: {}\n", 626 enterIdx, enterStat, ds.colStatus(enterIdx), objChange); 627 } 628 629 else 630 { 631 assert(enterId.isSPxRowId()); 632 enterIdx = this->number(SPxRowId(enterId)); 633 enterStat = ds.rowStatus(enterIdx); 634 assert(!isBasic(enterStat)); 635 636 /* For an #Id# to enter the basis we better recompute the Test value. 637 */ 638 if(rep() == ROW) 639 { 640 computePvec(enterIdx); 641 enterTest = computeTest(enterIdx); 642 theTest[enterIdx] = 0; 643 } 644 else 645 { 646 enterTest = coTest()[enterIdx]; 647 theCoTest[enterIdx] = 0; 648 } 649 650 switch(enterStat) 651 { 652 // primal/columnwise cases: 653 case SPxBasisBase<R>::Desc::P_ON_UPPER : 654 assert(rep() == COLUMN); 655 enterUB = theURbound[enterIdx]; 656 enterLB = theLRbound[enterIdx]; 657 enterVal = enterLB; 658 enterMax = enterUB - enterLB; 659 enterPric = (*theCoPvec)[enterIdx]; 660 enterRO = this->maxRowObj(enterIdx); 661 objChange -= enterRO * enterVal; 662 663 if(enterUB >= R(infinity)) 664 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::D_ON_LOWER; 665 else if(EQ(enterLB, enterUB, this->epsilon())) 666 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::D_FREE; 667 else 668 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::D_ON_BOTH; 669 670 break; 671 672 case SPxBasisBase<R>::Desc::P_ON_LOWER : 673 assert(rep() == COLUMN); 674 enterUB = theURbound[enterIdx]; 675 enterLB = theLRbound[enterIdx]; 676 enterVal = enterUB; 677 enterMax = enterLB - enterUB; 678 enterPric = (*theCoPvec)[enterIdx]; 679 enterRO = this->maxRowObj(enterIdx); 680 objChange -= enterRO * enterVal; 681 682 if(enterLB <= R(-infinity)) 683 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::D_ON_UPPER; 684 else if(EQ(enterLB, enterUB, this->epsilon())) 685 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::D_FREE; 686 else 687 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::D_ON_BOTH; 688 689 break; 690 691 case SPxBasisBase<R>::Desc::P_FREE : 692 assert(rep() == COLUMN); 693 #if 1 694 throw SPxInternalCodeException("XENTER02 This should never happen."); 695 #else 696 SPX_MSG_ERROR(std::cerr << "EENTER99 ERROR: not yet debugged!" << std::endl;) 697 enterPric = (*theCoPvec)[enterIdx]; 698 enterRO = this->maxRowObj(enterIdx); 699 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::D_UNDEFINED; 700 #endif 701 break; 702 703 // dual/rowwise cases: 704 case SPxBasisBase<R>::Desc::D_ON_UPPER : 705 assert(rep() == ROW); 706 assert(theURbound[enterIdx] < R(infinity)); 707 enterUB = theURbound[enterIdx]; 708 enterLB = R(-infinity); 709 enterVal = enterUB; 710 enterMax = R(-infinity); 711 enterPric = (*thePvec)[enterIdx]; 712 enterRO = this->lhs(enterIdx); 713 objChange -= enterRO * enterVal; 714 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_LOWER; 715 break; 716 717 case SPxBasisBase<R>::Desc::D_ON_LOWER : 718 assert(rep() == ROW); 719 assert(theLRbound[enterIdx] > R(-infinity)); 720 enterLB = theLRbound[enterIdx]; 721 enterUB = R(infinity); 722 enterVal = enterLB; 723 enterMax = R(infinity); 724 enterPric = (*thePvec)[enterIdx]; 725 enterRO = this->rhs(enterIdx); 726 objChange -= enterRO * enterVal; 727 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_UPPER; 728 break; 729 730 case SPxBasisBase<R>::Desc::D_FREE: 731 assert(rep() == ROW); 732 assert(this->rhs(enterIdx) == this->lhs(enterIdx)); 733 enterUB = R(infinity); 734 enterLB = R(-infinity); 735 enterVal = 0; 736 enterPric = (*thePvec)[enterIdx]; 737 enterRO = this->rhs(enterIdx); 738 enterMax = (enterPric > enterRO) ? R(infinity) : R(-infinity); 739 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::P_FIXED; 740 break; 741 742 case SPxBasisBase<R>::Desc::D_ON_BOTH : 743 assert(rep() == ROW); 744 enterPric = (*thePvec)[enterIdx]; 745 746 if(enterPric > this->rhs(enterIdx)) 747 { 748 enterLB = theLRbound[enterIdx]; 749 enterVal = enterLB; 750 enterUB = R(infinity); 751 enterMax = R(infinity); 752 enterRO = this->rhs(enterIdx); 753 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_UPPER; 754 } 755 else 756 { 757 enterUB = theURbound[enterIdx]; 758 enterVal = enterUB; 759 enterLB = R(-infinity); 760 enterMax = R(-infinity); 761 enterRO = this->lhs(enterIdx); 762 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_LOWER; 763 } 764 765 objChange -= theLRbound[enterIdx] * this->rhs(enterIdx); 766 objChange -= theURbound[enterIdx] * this->lhs(enterIdx); 767 break; 768 769 default: 770 throw SPxInternalCodeException("XENTER03 This should never happen."); 771 } 772 773 SPxOut::debug(this, "DENTER05 SPxSolverBase::getEnterVals() : row {}: {} -> {} objChange: {}\n", 774 enterIdx, enterStat, ds.rowStatus(enterIdx), objChange); 775 } 776 } 777 778 /* process leaving variable 779 */ 780 template <class R> 781 void SPxSolverBase<R>::getEnterVals2 782 ( 783 int leaveIdx, 784 R enterMax, 785 R& leavebound, 786 StableSum<R>& objChange 787 ) 788 { 789 int idx; 790 typename SPxBasisBase<R>::Desc& ds = this->desc(); 791 SPxId leftId = this->baseId(leaveIdx); 792 793 if(leftId.isSPxRowId()) 794 { 795 idx = this->number(SPxRowId(leftId)); 796 typename SPxBasisBase<R>::Desc::Status leaveStat = ds.rowStatus(idx); 797 798 // coverity[switch_selector_expr_is_constant] 799 switch(leaveStat) 800 { 801 case SPxBasisBase<R>::Desc::P_FIXED : 802 assert(rep() == ROW); 803 throw SPxInternalCodeException("XENTER04 This should never happen."); 804 break; 805 806 case SPxBasisBase<R>::Desc::P_ON_UPPER : 807 assert(rep() == ROW); 808 leavebound = theLBbound[leaveIdx]; 809 theLRbound[idx] = leavebound; 810 ds.rowStatus(idx) = this->dualRowStatus(idx); 811 812 switch(ds.rowStatus(idx)) 813 { 814 case SPxBasisBase<R>::Desc::D_ON_UPPER : 815 objChange += theURbound[idx] * this->lhs(idx); 816 break; 817 818 case SPxBasisBase<R>::Desc::D_ON_LOWER : 819 objChange += theLRbound[idx] * this->rhs(idx); 820 break; 821 822 case SPxBasisBase<R>::Desc::D_ON_BOTH : 823 objChange += theURbound[idx] * this->lhs(idx); 824 objChange += theLRbound[idx] * this->rhs(idx); 825 break; 826 827 default: 828 break; 829 } 830 831 break; 832 833 case SPxBasisBase<R>::Desc::P_ON_LOWER : 834 assert(rep() == ROW); 835 leavebound = theUBbound[leaveIdx]; 836 theURbound[idx] = leavebound; 837 ds.rowStatus(idx) = this->dualRowStatus(idx); 838 839 switch(ds.rowStatus(idx)) 840 { 841 case SPxBasisBase<R>::Desc::D_ON_UPPER : 842 objChange += theURbound[idx] * this->lhs(idx); 843 break; 844 845 case SPxBasisBase<R>::Desc::D_ON_LOWER : 846 objChange += theLRbound[idx] * this->rhs(idx); 847 break; 848 849 case SPxBasisBase<R>::Desc::D_ON_BOTH : 850 objChange += theURbound[idx] * this->lhs(idx); 851 objChange += theLRbound[idx] * this->rhs(idx); 852 break; 853 854 default: 855 break; 856 } 857 858 break; 859 860 case SPxBasisBase<R>::Desc::P_FREE : 861 assert(rep() == ROW); 862 #if 1 863 throw SPxInternalCodeException("XENTER05 This should never happen."); 864 #else 865 SPX_MSG_ERROR(std::cerr << "EENTER98 ERROR: not yet debugged!" << std::endl;) 866 867 if((*theCoPvec)[leaveIdx] - theLBbound[leaveIdx] < 868 theUBbound[leaveIdx] - (*theCoPvec)[leaveIdx]) 869 { 870 leavebound = theLBbound[leaveIdx]; 871 theLRbound[idx] = leavebound; 872 } 873 else 874 { 875 leavebound = theUBbound[leaveIdx]; 876 theURbound[idx] = leavebound; 877 } 878 879 ds.rowStatus(idx) = SPxBasisBase<R>::Desc::D_UNDEFINED; 880 #endif 881 break; 882 883 // primal/columnwise cases: 884 case SPxBasisBase<R>::Desc::D_UNDEFINED : 885 assert(rep() == COLUMN); 886 throw SPxInternalCodeException("XENTER06 This should never happen."); 887 break; 888 889 case SPxBasisBase<R>::Desc::D_FREE : 890 assert(rep() == COLUMN); 891 892 if(theFvec->delta()[leaveIdx] * enterMax < 0) 893 leavebound = theUBbound[leaveIdx]; 894 else 895 leavebound = theLBbound[leaveIdx]; 896 897 theLRbound[idx] = leavebound; 898 theURbound[idx] = leavebound; 899 objChange += leavebound * this->maxRowObj(leaveIdx); 900 ds.rowStatus(idx) = SPxBasisBase<R>::Desc::P_FIXED; 901 break; 902 903 case SPxBasisBase<R>::Desc::D_ON_UPPER : 904 assert(rep() == COLUMN); 905 leavebound = theUBbound[leaveIdx]; 906 theURbound[idx] = leavebound; 907 objChange += leavebound * this->maxRowObj(leaveIdx); 908 ds.rowStatus(idx) = SPxBasisBase<R>::Desc::P_ON_LOWER; 909 break; 910 911 case SPxBasisBase<R>::Desc::D_ON_LOWER : 912 assert(rep() == COLUMN); 913 leavebound = theLBbound[leaveIdx]; 914 theLRbound[idx] = leavebound; 915 objChange += leavebound * this->maxRowObj(leaveIdx); 916 ds.rowStatus(idx) = SPxBasisBase<R>::Desc::P_ON_UPPER; 917 break; 918 919 case SPxBasisBase<R>::Desc::D_ON_BOTH : 920 assert(rep() == COLUMN); 921 922 if(enterMax * theFvec->delta()[leaveIdx] < 0) 923 { 924 leavebound = theUBbound[leaveIdx]; 925 theURbound[idx] = leavebound; 926 objChange += leavebound * this->maxRowObj(leaveIdx); 927 ds.rowStatus(idx) = SPxBasisBase<R>::Desc::P_ON_LOWER; 928 } 929 else 930 { 931 leavebound = theLBbound[leaveIdx]; 932 theLRbound[idx] = leavebound; 933 objChange += leavebound * this->maxRowObj(leaveIdx); 934 ds.rowStatus(idx) = SPxBasisBase<R>::Desc::P_ON_UPPER; 935 } 936 937 break; 938 939 default: 940 throw SPxInternalCodeException("XENTER07 This should never happen."); 941 } 942 943 SPxOut::debug(this, "DENTER06 SPxSolverBase::getEnterVals2(): row {}: {} -> {} objChange: {}\n", 944 idx, leaveStat, ds.rowStatus(idx), objChange); 945 } 946 947 else 948 { 949 assert(leftId.isSPxColId()); 950 idx = this->number(SPxColId(leftId)); 951 typename SPxBasisBase<R>::Desc::Status leaveStat = ds.colStatus(idx); 952 953 // coverity[switch_selector_expr_is_constant] 954 switch(leaveStat) 955 { 956 case SPxBasisBase<R>::Desc::P_ON_UPPER : 957 assert(rep() == ROW); 958 leavebound = theLBbound[leaveIdx]; 959 theLCbound[idx] = leavebound; 960 ds.colStatus(idx) = this->dualColStatus(idx); 961 962 switch(ds.colStatus(idx)) 963 { 964 case SPxBasisBase<R>::Desc::D_ON_UPPER : 965 objChange += theUCbound[idx] * this->lower(idx); 966 break; 967 968 case SPxBasisBase<R>::Desc::D_ON_LOWER : 969 objChange += theLCbound[idx] * this->upper(idx); 970 break; 971 972 case SPxBasisBase<R>::Desc::D_ON_BOTH : 973 objChange += theLCbound[idx] * this->upper(idx); 974 objChange += theUCbound[idx] * this->lower(idx); 975 break; 976 977 default: 978 break; 979 } 980 981 break; 982 983 case SPxBasisBase<R>::Desc::P_ON_LOWER : 984 assert(rep() == ROW); 985 leavebound = theUBbound[leaveIdx]; 986 theUCbound[idx] = leavebound; 987 ds.colStatus(idx) = this->dualColStatus(idx); 988 989 switch(ds.colStatus(idx)) 990 { 991 case SPxBasisBase<R>::Desc::D_ON_UPPER : 992 objChange += theUCbound[idx] * this->lower(idx); 993 break; 994 995 case SPxBasisBase<R>::Desc::D_ON_LOWER : 996 objChange += theLCbound[idx] * this->upper(idx); 997 break; 998 999 case SPxBasisBase<R>::Desc::D_ON_BOTH : 1000 objChange += theLCbound[idx] * this->upper(idx); 1001 objChange += theUCbound[idx] * this->lower(idx); 1002 break; 1003 1004 default: 1005 break; 1006 } 1007 1008 break; 1009 1010 case SPxBasisBase<R>::Desc::P_FREE : 1011 assert(rep() == ROW); 1012 1013 if(theFvec->delta()[leaveIdx] * enterMax > 0) 1014 { 1015 leavebound = theLBbound[leaveIdx]; 1016 theLCbound[idx] = leavebound; 1017 } 1018 else 1019 { 1020 leavebound = theUBbound[leaveIdx]; 1021 theUCbound[idx] = leavebound; 1022 } 1023 1024 ds.colStatus(idx) = SPxBasisBase<R>::Desc::D_UNDEFINED; 1025 break; 1026 1027 case SPxBasisBase<R>::Desc::P_FIXED: 1028 assert(rep() == ROW); 1029 throw SPxInternalCodeException("XENTER08 This should never happen."); 1030 break; 1031 1032 // primal/columnwise cases: 1033 case SPxBasisBase<R>::Desc::D_FREE : 1034 assert(rep() == COLUMN); 1035 1036 if(theFvec->delta()[leaveIdx] * enterMax > 0) 1037 leavebound = theLBbound[leaveIdx]; 1038 else 1039 leavebound = theUBbound[leaveIdx]; 1040 1041 theUCbound[idx] = 1042 theLCbound[idx] = leavebound; 1043 objChange += this->maxObj(idx) * leavebound; 1044 ds.colStatus(idx) = SPxBasisBase<R>::Desc::P_FIXED; 1045 break; 1046 1047 case SPxBasisBase<R>::Desc::D_ON_UPPER : 1048 assert(rep() == COLUMN); 1049 leavebound = theLBbound[leaveIdx]; 1050 theLCbound[idx] = leavebound; 1051 objChange += this->maxObj(idx) * leavebound; 1052 ds.colStatus(idx) = SPxBasisBase<R>::Desc::P_ON_LOWER; 1053 break; 1054 1055 case SPxBasisBase<R>::Desc::D_ON_LOWER : 1056 assert(rep() == COLUMN); 1057 leavebound = theUBbound[leaveIdx]; 1058 theUCbound[idx] = leavebound; 1059 objChange += this->maxObj(idx) * leavebound; 1060 ds.colStatus(idx) = SPxBasisBase<R>::Desc::P_ON_UPPER; 1061 break; 1062 1063 case SPxBasisBase<R>::Desc::D_ON_BOTH : 1064 case SPxBasisBase<R>::Desc::D_UNDEFINED : 1065 assert(rep() == COLUMN); 1066 1067 if(enterMax * theFvec->delta()[leaveIdx] < 0) 1068 { 1069 leavebound = theUBbound[leaveIdx]; 1070 theUCbound[idx] = leavebound; 1071 objChange += this->maxObj(idx) * leavebound; 1072 ds.colStatus(idx) = SPxBasisBase<R>::Desc::P_ON_UPPER; 1073 } 1074 else 1075 { 1076 leavebound = theLBbound[leaveIdx]; 1077 theLCbound[idx] = leavebound; 1078 objChange += this->maxObj(idx) * leavebound; 1079 ds.colStatus(idx) = SPxBasisBase<R>::Desc::P_ON_LOWER; 1080 } 1081 1082 break; 1083 1084 default: 1085 throw SPxInternalCodeException("XENTER09 This should never happen."); 1086 } 1087 1088 SPxOut::debug(this, "DENTER07 SPxSolverBase::getEnterVals2(): col {}: {} -> {} objChange: {}\n", 1089 idx, leaveStat, ds.colStatus(idx), objChange); 1090 } 1091 } 1092 1093 template <class R> 1094 void 1095 SPxSolverBase<R>::ungetEnterVal( 1096 SPxId enterId, 1097 typename SPxBasisBase<R>::Desc::Status enterStat, 1098 R leaveVal, 1099 const SVectorBase<R>& vec, 1100 StableSum<R>& objChange 1101 ) 1102 { 1103 assert(rep() == COLUMN); 1104 int enterIdx; 1105 typename SPxBasisBase<R>::Desc& ds = this->desc(); 1106 1107 if(enterId.isSPxColId()) 1108 { 1109 enterIdx = this->number(SPxColId(enterId)); 1110 1111 if(enterStat == SPxBasisBase<R>::Desc::P_ON_UPPER) 1112 { 1113 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_LOWER; 1114 objChange += theLCbound[enterIdx] * this->maxObj(enterIdx); 1115 } 1116 else 1117 { 1118 ds.colStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_UPPER; 1119 objChange += theUCbound[enterIdx] * this->maxObj(enterIdx); 1120 } 1121 1122 theFrhs->multAdd(leaveVal, vec); 1123 } 1124 else 1125 { 1126 enterIdx = this->number(SPxRowId(enterId)); 1127 assert(enterId.isSPxRowId()); 1128 1129 if(enterStat == SPxBasisBase<R>::Desc::P_ON_UPPER) 1130 { 1131 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_LOWER; 1132 objChange += (theURbound[enterIdx]) * this->maxRowObj(enterIdx); 1133 } 1134 else 1135 { 1136 ds.rowStatus(enterIdx) = SPxBasisBase<R>::Desc::P_ON_UPPER; 1137 objChange += (theLRbound[enterIdx]) * this->maxRowObj(enterIdx); 1138 } 1139 1140 (*theFrhs)[enterIdx] += leaveVal; 1141 } 1142 1143 if(isId(enterId)) 1144 { 1145 theTest[enterIdx] = 0; 1146 isInfeasibleCo[enterIdx] = SPxPricer<R>::NOT_VIOLATED; 1147 } 1148 else 1149 { 1150 theCoTest[enterIdx] = 0; 1151 isInfeasible[enterIdx] = SPxPricer<R>::NOT_VIOLATED; 1152 } 1153 } 1154 1155 template <class R> 1156 void SPxSolverBase<R>::rejectEnter( 1157 SPxId enterId, 1158 R enterTest, 1159 typename SPxBasisBase<R>::Desc::Status enterStat 1160 ) 1161 { 1162 int enterIdx = this->number(enterId); 1163 1164 if(isId(enterId)) 1165 { 1166 theTest[enterIdx] = enterTest; 1167 this->desc().status(enterIdx) = enterStat; 1168 } 1169 else 1170 { 1171 theCoTest[enterIdx] = enterTest; 1172 this->desc().coStatus(enterIdx) = enterStat; 1173 } 1174 } 1175 1176 template <class R> 1177 void SPxSolverBase<R>::computePrimalray4Col(R direction, SPxId enterId) 1178 { 1179 R sign = (direction > 0 ? 1.0 : -1.0); 1180 1181 primalRay.clear(); 1182 primalRay.setMax(fVec().delta().size() + 1); 1183 1184 for(int j = 0; j < fVec().delta().size(); ++j) 1185 { 1186 SPxId i = this->baseId(fVec().idx().index(j)); 1187 1188 if(i.isSPxColId()) 1189 primalRay.add(this->number(SPxColId(i)), sign * fVec().delta().value(j)); 1190 } 1191 1192 if(enterId.isSPxColId()) 1193 primalRay.add(this->number(SPxColId(enterId)), -sign); 1194 } 1195 1196 template <class R> 1197 void SPxSolverBase<R>::computeDualfarkas4Row(R direction, SPxId enterId) 1198 { 1199 R sign = (direction > 0 ? -1.0 : 1.0); 1200 1201 dualFarkas.clear(); 1202 dualFarkas.setMax(fVec().delta().size() + 1); 1203 1204 for(int j = 0; j < fVec().delta().size(); ++j) 1205 { 1206 SPxId spxid = this->baseId(fVec().idx().index(j)); 1207 1208 if(spxid.isSPxRowId()) 1209 dualFarkas.add(this->number(SPxRowId(spxid)), sign * fVec().delta().value(j)); 1210 } 1211 1212 if(enterId.isSPxRowId()) 1213 dualFarkas.add(this->number(SPxRowId(enterId)), -sign); 1214 } 1215 1216 template <class R> 1217 bool SPxSolverBase<R>::enter(SPxId& enterId, bool polish) 1218 { 1219 assert(enterId.isValid()); 1220 assert(type() == ENTER); 1221 assert(initialized); 1222 1223 SPxId none; // invalid id used when enter fails 1224 R enterTest; // correct test value of entering var 1225 R enterUB; // upper bound of entering variable 1226 R enterLB; // lower bound of entering variable 1227 R enterVal; // current value of entering variable 1228 R enterMax; // maximum value for entering shift 1229 R enterPric; // priced value of entering variable 1230 typename SPxBasisBase<R>::Desc::Status enterStat; // status of entering variable 1231 R enterRO; // rhs/obj of entering variable 1232 StableSum<R> objChange; 1233 const SVectorBase<R>* enterVec = enterVector(enterId); 1234 1235 bool instable = instableEnter; 1236 assert(!instable || instableEnterId.isValid()); 1237 1238 getEnterVals(enterId, enterTest, enterUB, enterLB, 1239 enterVal, enterMax, enterPric, enterStat, enterRO, objChange); 1240 1241 if(!polish && enterTest > -epsilon()) 1242 { 1243 rejectEnter(enterId, enterTest, enterStat); 1244 this->change(-1, none, 0); 1245 1246 SPxOut::debug(this, "DENTER08 rejecting false enter pivot\n"); 1247 1248 return false; 1249 } 1250 1251 /* Before performing the actual basis update, we must determine, how this 1252 is to be accomplished. 1253 */ 1254 // BH 2005-11-15: Obviously solve4update() is only called if theFvec.delta() 1255 // is setup (i.e. the indices of the NZEs are stored within it) and there are 1256 // 0 NZEs (???). 1257 // In that case theFvec->delta() is set such that 1258 // Base * theFvec->delta() = enterVec 1259 if(theFvec->delta().isSetup() && theFvec->delta().size() == 0) 1260 SPxBasisBase<R>::solve4update(theFvec->delta(), *enterVec); 1261 1262 #ifdef ENABLE_ADDITIONAL_CHECKS 1263 else 1264 { 1265 // BH 2005-11-29: This code block seems to check the assertion 1266 // || Base * theFvec->delta() - enterVec ||_2 <= entertol() 1267 VectorBase<R> tmp(dim()); 1268 // BH 2005-11-15: This cast is necessary since SSVectorBase<R> inherits protected from VectorBase<R>. 1269 tmp = reinterpret_cast<VectorBase<R>&>(theFvec->delta()); 1270 this->multBaseWith(tmp); 1271 tmp -= *enterVec; 1272 1273 if(tmp.length() > entertol()) 1274 { 1275 // This happens frequently and does usually not hurt, so print these 1276 // warnings only with verbose level INFO2 and higher. 1277 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "WENTER09 fVec updated error = " 1278 << tmp.length() << std::endl;) 1279 } 1280 } 1281 1282 #endif // ENABLE_ADDITIONAL_CHECKS 1283 1284 if(!polish && m_numCycle > m_maxCycle) 1285 { 1286 if(-enterMax > 0) 1287 perturbMaxEnter(); 1288 else 1289 perturbMinEnter(); 1290 } 1291 1292 R leaveVal = -enterMax; 1293 1294 boundflips = 0; 1295 int leaveIdx = theratiotester->selectLeave(leaveVal, enterTest, polish); 1296 1297 /* in row representation, fixed columns and rows should not leave the basis */ 1298 assert(leaveIdx < 0 || !this->baseId(leaveIdx).isSPxColId() 1299 || this->desc().colStatus(this->number(SPxColId(this->baseId(leaveIdx)))) != 1300 SPxBasisBase<R>::Desc::P_FIXED); 1301 assert(leaveIdx < 0 || !this->baseId(leaveIdx).isSPxRowId() 1302 || this->desc().rowStatus(this->number(SPxRowId(this->baseId(leaveIdx)))) != 1303 SPxBasisBase<R>::Desc::P_FIXED); 1304 1305 instableEnterVal = 0; 1306 instableEnterId = SPxId(); 1307 instableEnter = false; 1308 1309 /* 1310 We now tried to find a variable to leave the basis. If one has been 1311 found, a regular basis update is to be performed. 1312 */ 1313 if(leaveIdx >= 0) 1314 { 1315 if(spxAbs(leaveVal) < entertol()) 1316 { 1317 if(NE(theUBbound[leaveIdx], theLBbound[leaveIdx], this->tolerances()->epsilon()) 1318 && enterStat != SPxBasisBase<R>::Desc::P_FREE && enterStat != SPxBasisBase<R>::Desc::D_FREE) 1319 { 1320 m_numCycle++; 1321 enterCycles++; 1322 } 1323 } 1324 else 1325 m_numCycle /= 2; 1326 1327 // setup for updating the copricing vector 1328 if(coSolveVector3 && coSolveVector2) 1329 { 1330 assert(coSolveVector2->isConsistent()); 1331 assert(coSolveVector2rhs->isSetup()); 1332 assert(coSolveVector3->isConsistent()); 1333 assert(coSolveVector3rhs->isSetup()); 1334 assert(boundflips > 0); 1335 SPxBasisBase<R>::coSolve(theCoPvec->delta(), *coSolveVector2, *coSolveVector3 1336 , unitVecs[leaveIdx], *coSolveVector2rhs, *coSolveVector3rhs); 1337 (*theCoPvec) -= (*coSolveVector3); 1338 } 1339 else if(coSolveVector3) 1340 { 1341 assert(coSolveVector3->isConsistent()); 1342 assert(coSolveVector3rhs->isSetup()); 1343 assert(boundflips > 0); 1344 SPxBasisBase<R>::coSolve(theCoPvec->delta(), *coSolveVector3, unitVecs[leaveIdx], 1345 *coSolveVector3rhs); 1346 (*theCoPvec) -= (*coSolveVector3); 1347 } 1348 else if(coSolveVector2) 1349 SPxBasisBase<R>::coSolve(theCoPvec->delta(), *coSolveVector2, unitVecs[leaveIdx], 1350 *coSolveVector2rhs); 1351 else 1352 SPxBasisBase<R>::coSolve(theCoPvec->delta(), unitVecs[leaveIdx]); 1353 1354 if(boundflips > 0) 1355 { 1356 for(int i = coSolveVector3->dim() - 1; i >= 0; --i) 1357 { 1358 if(spxAbs((*coSolveVector3)[i]) > epsilon()) 1359 (*thePvec).multAdd(-(*coSolveVector3)[i], (*thecovectors)[i]); 1360 } 1361 1362 // we need to update enterPric in case it was changed by bound flips 1363 if(enterId.isSPxColId()) 1364 enterPric = (*theCoPvec)[this->number(SPxColId(enterId))]; 1365 else 1366 enterPric = (*thePvec)[this->number(SPxRowId(enterId))]; 1367 1368 SPxOut::debug(this, "IEBFRT02 breakpoints passed / bounds flipped = {}\n", boundflips); 1369 totalboundflips += boundflips; 1370 } 1371 1372 (*theCoPrhs)[leaveIdx] = enterRO; 1373 theCoPvec->value() = (enterRO - enterPric) / theFvec->delta()[leaveIdx]; 1374 1375 if(theCoPvec->value() > epsilon() || theCoPvec->value() < -epsilon()) 1376 { 1377 if(pricing() == FULL) 1378 { 1379 thePvec->value() = theCoPvec->value(); 1380 setupPupdate(); 1381 } 1382 1383 doPupdate(); 1384 } 1385 1386 assert(thePvec->isConsistent()); 1387 assert(theCoPvec->isConsistent()); 1388 1389 assert(!this->baseId(leaveIdx).isSPxRowId() 1390 || this->desc().rowStatus(this->number(SPxRowId(this->baseId(leaveIdx)))) != 1391 SPxBasisBase<R>::Desc::P_FIXED); 1392 assert(!this->baseId(leaveIdx).isSPxColId() 1393 || this->desc().colStatus(this->number(SPxColId(this->baseId(leaveIdx)))) != 1394 SPxBasisBase<R>::Desc::P_FIXED); 1395 1396 R leavebound; // bound on which leaving variable moves 1397 1398 try 1399 { 1400 getEnterVals2(leaveIdx, enterMax, leavebound, objChange); 1401 } 1402 catch(const SPxException& F) 1403 { 1404 rejectEnter(enterId, enterTest, enterStat); 1405 this->change(-1, none, 0); 1406 throw F; 1407 } 1408 1409 // process entering variable 1410 theUBbound[leaveIdx] = enterUB; 1411 theLBbound[leaveIdx] = enterLB; 1412 1413 // compute tests: 1414 updateCoTest(); 1415 1416 if(pricing() == FULL) 1417 updateTest(); 1418 1419 // update feasibility vectors 1420 theFvec->value() = leaveVal; 1421 theFvec->update(); 1422 (*theFvec)[leaveIdx] = enterVal - leaveVal; 1423 1424 if(leavebound > epsilon() || leavebound < -epsilon()) 1425 theFrhs->multAdd(-leavebound, this->baseVec(leaveIdx)); 1426 1427 if(enterVal > epsilon() || enterVal < -epsilon()) 1428 theFrhs->multAdd(enterVal, *enterVec); 1429 1430 // update objective funtion value 1431 updateNonbasicValue(objChange); 1432 1433 // change basis matrix 1434 this->change(leaveIdx, enterId, enterVec, &(theFvec->delta())); 1435 1436 return true; 1437 } 1438 /* No leaving vector could be found that would yield a stable pivot step. 1439 */ 1440 else if(NE(leaveVal, -enterMax, this->tolerances()->epsilon())) 1441 { 1442 /* In the ENTER algorithm, when for a selected entering variable we find only 1443 an instable leaving variable, then the basis change is not conducted. 1444 Instead, we save the entering variable's id in instableEnterId and set 1445 the test value to zero, hoping to find a different leaving 1446 variable with a stable leavingvariable. 1447 If this fails, however, and no more entering variable is found, we have to 1448 perform the instable basis change using instableEnterId. In this (and only 1449 in this) case, the flag instableEnter is set to true. 1450 1451 leaveVal != enterMax is the case that selectLeave has found only an instable leaving 1452 variable. We store this leaving variable for later if we are not already in the 1453 instable case */ 1454 1455 if(!instable) 1456 { 1457 instableEnterId = enterId; 1458 instableEnterVal = enterTest; 1459 1460 SPxOut::debug(this, "DENTER09 rejecting enter pivot and looking for others\n"); 1461 1462 rejectEnter(enterId, enterTest / 10.0, enterStat); 1463 this->change(-1, none, 0); 1464 } 1465 else 1466 { 1467 SPxOut::debug(this, "DENTER10 rejecting enter pivot in instable state, resetting values\n"); 1468 rejectEnter(enterId, enterTest, enterStat); 1469 this->change(-1, none, 0); 1470 } 1471 1472 return false; 1473 } 1474 /* No leaving vector has been selected from the basis. However, if the 1475 shift amount for |fVec| is bounded, we are in the case, that the 1476 entering variable is moved from one bound to its other, before any of 1477 the basis feasibility variables reaches their bound. This may only 1478 happen in primal/columnwise case with upper and lower bounds on 1479 variables. 1480 */ 1481 else if(!polish && leaveVal < R(infinity) && leaveVal > R(-infinity)) 1482 { 1483 assert(rep() == COLUMN); 1484 assert(EQ(leaveVal, -enterMax, this->epsilon())); 1485 1486 this->change(-1, enterId, enterVec); 1487 1488 theFvec->value() = leaveVal; 1489 theFvec->update(); 1490 1491 ungetEnterVal(enterId, enterStat, leaveVal, *enterVec, objChange); 1492 1493 // update objective funtion value 1494 updateNonbasicValue(objChange); 1495 1496 SPxOut::debug(this, "DENTER11 moving entering variable from one bound to the other\n"); 1497 1498 return false; 1499 } 1500 /* No variable could be selected to leave the basis and even the entering 1501 variable is unbounded --- this is a failure. 1502 */ 1503 else 1504 { 1505 /* The following line originally was in the "lastUpdate() > 1" case; 1506 we need it in the INFEASIBLE/UNBOUNDED case, too, to have the 1507 basis descriptor at the correct size. 1508 */ 1509 rejectEnter(enterId, enterTest, enterStat); 1510 this->change(-1, none, 0); 1511 1512 if(polish) 1513 return false; 1514 1515 else if(this->lastUpdate() > 1) 1516 { 1517 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "IENTER01 factorization triggered in " 1518 << "enter() for feasibility test" << std::endl;) 1519 1520 try 1521 { 1522 factorize(); 1523 } 1524 catch(const SPxStatusException& E) 1525 { 1526 // don't exit immediately but handle the singularity correctly 1527 assert(SPxBasisBase<R>::status() == SPxBasisBase<R>::SINGULAR); 1528 SPX_MSG_INFO3((*this->spxout), 1529 (*this->spxout) << "Caught exception in factorization: " << E.what() << 1530 std::endl;) 1531 } 1532 1533 /* after a factorization, the entering column/row might not be infeasible or suboptimal anymore, hence we do 1534 * not try to call leave(leaveIdx), but rather return to the main solving loop and call the pricer again 1535 */ 1536 return false; 1537 } 1538 1539 /* do not exit with status infeasible or unbounded if there is only a very small violation 1540 * ROW: recompute the primal variables and activities for another, more precise, round of pricing 1541 */ 1542 else if(spxAbs(enterTest) < entertol()) 1543 { 1544 SPX_MSG_INFO3((*this->spxout), 1545 (*this->spxout) << "IENTER11 clean up step to reduce numerical errors" << 1546 std::endl;) 1547 1548 SPxBasisBase<R>::coSolve(*theCoPvec, *theCoPrhs); 1549 computePvec(); 1550 computeCoTest(); 1551 computeTest(); 1552 1553 return false; 1554 } 1555 1556 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "IENTER02 unboundedness/infeasibility found in " 1557 << "enter()" << std::endl;) 1558 1559 if(rep() == ROW) 1560 { 1561 computeDualfarkas4Row(leaveVal, enterId); 1562 setBasisStatus(SPxBasisBase<R>::INFEASIBLE); 1563 } 1564 /**@todo if shift() is not zero, we must not conclude primal unboundedness */ 1565 else 1566 { 1567 computePrimalray4Col(leaveVal, enterId); 1568 setBasisStatus(SPxBasisBase<R>::UNBOUNDED); 1569 } 1570 1571 return false; 1572 } 1573 } 1574 } // namespace soplex 1575