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 /* Updating the Basis for Leaving Variables 26 */ 27 #include <assert.h> 28 #include <stdio.h> 29 30 #include "soplex/spxdefines.h" 31 #include "soplex/spxpricer.h" 32 #include "soplex/spxsolver.h" 33 #include "soplex/spxratiotester.h" 34 #include "soplex/spxout.h" 35 #include "soplex/exceptions.h" 36 37 namespace soplex 38 { 39 static const Real default_reject_leave_tol = 1e-10; // = SOPLEX_LOWSTAB as defined in spxfastrt.hpp 40 41 /* 42 VectorBase<R> |fTest| gives the feasibility test of all basic variables. For its 43 computation |fVec|, |theUBbound| and |theLBbound| must be setup correctly. 44 Values of |fTest| $<0$ represent infeasible variables, which are eligible 45 for leaving the basis in the simplex loop. 46 */ 47 template <class R> 48 void SPxSolverBase<R>::computeFtest() 49 { 50 51 assert(type() == LEAVE); 52 53 R theeps = entertol(); 54 m_pricingViolUpToDate = true; 55 m_pricingViolCoUpToDate = true; 56 m_pricingViol = 0; 57 m_pricingViolCo = 0; 58 m_numViol = 0; 59 infeasibilities.clear(); 60 int sparsitythreshold = (int)(sparsePricingFactor * dim()); 61 62 for(int i = 0; i < dim(); ++i) 63 { 64 theCoTest[i] = ((*theFvec)[i] > theUBbound[i]) 65 ? theUBbound[i] - (*theFvec)[i] 66 : (*theFvec)[i] - theLBbound[i]; 67 68 if(remainingRoundsLeave == 0) 69 { 70 if(theCoTest[i] < -theeps) 71 { 72 m_pricingViol -= theCoTest[i]; 73 infeasibilities.addIdx(i); 74 isInfeasible[i] = SPxPricer<R>::VIOLATED; 75 ++m_numViol; 76 } 77 else 78 isInfeasible[i] = SPxPricer<R>::NOT_VIOLATED; 79 80 if(infeasibilities.size() > sparsitythreshold) 81 { 82 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << " --- using dense pricing" 83 << std::endl;) 84 remainingRoundsLeave = SOPLEX_DENSEROUNDS; 85 sparsePricingLeave = false; 86 infeasibilities.clear(); 87 } 88 } 89 else if(theCoTest[i] < -theeps) 90 { 91 m_pricingViol -= theCoTest[i]; 92 m_numViol++; 93 } 94 } 95 96 if(infeasibilities.size() == 0 && !sparsePricingLeave) 97 { 98 --remainingRoundsLeave; 99 } 100 else if(infeasibilities.size() <= sparsitythreshold && !sparsePricingLeave) 101 { 102 SPX_MSG_INFO2((*this->spxout), 103 std::streamsize prec = spxout->precision(); 104 105 if(hyperPricingLeave) 106 (*this->spxout) << " --- using hypersparse pricing, "; 107 else 108 (*this->spxout) << " --- using sparse pricing, "; 109 (*this->spxout) << "sparsity: " 110 << std::setw(6) << std::fixed << std::setprecision(4) 111 << (R) m_numViol / dim() 112 << std::scientific << std::setprecision(int(prec)) 113 << std::endl;) 114 sparsePricingLeave = true; 115 } 116 } 117 118 template <class R> 119 void SPxSolverBase<R>::updateFtest() 120 { 121 const IdxSet& idx = theFvec->idx(); 122 VectorBase<R>& ftest = theCoTest; // |== fTest()| 123 assert(&ftest == &fTest()); 124 125 assert(type() == LEAVE); 126 127 updateViols.clear(); 128 R theeps = entertol(); 129 130 for(int j = idx.size() - 1; j >= 0; --j) 131 { 132 int i = idx.index(j); 133 134 if(m_pricingViolUpToDate && ftest[i] < -theeps) 135 // violation was present before this iteration 136 m_pricingViol += ftest[i]; 137 138 ftest[i] = ((*theFvec)[i] > theUBbound[i]) 139 ? theUBbound[i] - (*theFvec)[i] 140 : (*theFvec)[i] - theLBbound[i]; 141 142 if(sparsePricingLeave && ftest[i] < -theeps) 143 { 144 assert(remainingRoundsLeave == 0); 145 146 if(m_pricingViolUpToDate) 147 m_pricingViol -= ftest[i]; 148 149 if(isInfeasible[i] == SPxPricer<R>::NOT_VIOLATED) 150 { 151 // this can cause problems - we cannot keep on adding indeces to infeasibilities, 152 // because they are not deleted in hyper mode... 153 // if( !hyperPricingLeave ) 154 infeasibilities.addIdx(i); 155 isInfeasible[i] = SPxPricer<R>::VIOLATED; 156 } 157 158 if(hyperPricingLeave) 159 updateViols.addIdx(i); 160 } 161 else if(m_pricingViolUpToDate && ftest[i] < -theeps) 162 m_pricingViol -= ftest[i]; 163 } 164 165 // if boundflips were performed, we need to update these indices as well 166 if(boundflips > 0) 167 { 168 R eps = epsilon(); 169 170 for(int j = 0; j < solveVector3->size(); ++j) 171 { 172 if(spxAbs(solveVector3->value(j)) > eps) 173 { 174 int i = solveVector3->index(j); 175 176 if(m_pricingViolUpToDate && ftest[i] < -theeps) 177 m_pricingViol += ftest[i]; 178 179 ftest[i] = ((*theFvec)[i] > theUBbound[i]) ? theUBbound[i] - (*theFvec)[i] : 180 (*theFvec)[i] - theLBbound[i]; 181 182 if(sparsePricingLeave && ftest[i] < -theeps) 183 { 184 assert(remainingRoundsLeave == 0); 185 186 if(m_pricingViolUpToDate) 187 m_pricingViol -= ftest[i]; 188 189 if(!isInfeasible[i]) 190 { 191 infeasibilities.addIdx(i); 192 isInfeasible[i] = true; 193 } 194 } 195 else if(m_pricingViolUpToDate && ftest[i] < -theeps) 196 m_pricingViol -= ftest[i]; 197 } 198 } 199 } 200 } 201 202 203 /* compute statistics on leaving variable 204 Compute a set of statistical values on the variable selected for leaving the 205 basis. 206 */ 207 template <class R> 208 void SPxSolverBase<R>::getLeaveVals( 209 int leaveIdx, 210 typename SPxBasisBase<R>::Desc::Status& leaveStat, 211 SPxId& leaveId, 212 R& leaveMax, 213 R& leavebound, 214 int& leaveNum, 215 StableSum<R>& objChange) 216 { 217 typename SPxBasisBase<R>::Desc& ds = this->desc(); 218 leaveId = this->baseId(leaveIdx); 219 220 if(leaveId.isSPxRowId()) 221 { 222 leaveNum = this->number(SPxRowId(leaveId)); 223 leaveStat = ds.rowStatus(leaveNum); 224 225 assert(isBasic(leaveStat)); 226 227 switch(leaveStat) 228 { 229 case SPxBasisBase<R>::Desc::P_ON_UPPER : 230 assert(rep() == ROW); 231 ds.rowStatus(leaveNum) = this->dualRowStatus(leaveNum); 232 leavebound = 0; 233 leaveMax = R(-infinity); 234 break; 235 236 case SPxBasisBase<R>::Desc::P_ON_LOWER : 237 assert(rep() == ROW); 238 ds.rowStatus(leaveNum) = this->dualRowStatus(leaveNum); 239 leavebound = 0; 240 leaveMax = R(infinity); 241 break; 242 243 case SPxBasisBase<R>::Desc::P_FREE : 244 assert(rep() == ROW); 245 throw SPxInternalCodeException("XLEAVE01 This should never happen."); 246 247 case SPxBasisBase<R>::Desc::D_FREE : 248 assert(rep() == COLUMN); 249 ds.rowStatus(leaveNum) = SPxBasisBase<R>::Desc::P_FIXED; 250 assert(this->lhs(leaveNum) == this->rhs(leaveNum)); 251 leavebound = -this->rhs(leaveNum); 252 253 if((*theFvec)[leaveIdx] < theLBbound[leaveIdx]) 254 leaveMax = R(infinity); 255 else 256 leaveMax = R(-infinity); 257 258 break; 259 260 case SPxBasisBase<R>::Desc::D_ON_LOWER : 261 assert(rep() == COLUMN); 262 ds.rowStatus(leaveNum) = SPxBasisBase<R>::Desc::P_ON_UPPER; 263 leavebound = -this->rhs(leaveNum); // slack !! 264 leaveMax = R(infinity); 265 objChange += theLRbound[leaveNum] * this->rhs(leaveNum); 266 break; 267 268 case SPxBasisBase<R>::Desc::D_ON_UPPER : 269 assert(rep() == COLUMN); 270 ds.rowStatus(leaveNum) = SPxBasisBase<R>::Desc::P_ON_LOWER; 271 leavebound = -this->lhs(leaveNum); // slack !! 272 leaveMax = R(-infinity); 273 objChange += theURbound[leaveNum] * this->lhs(leaveNum); 274 break; 275 276 case SPxBasisBase<R>::Desc::D_ON_BOTH : 277 assert(rep() == COLUMN); 278 279 if((*theFvec)[leaveIdx] > theLBbound[leaveIdx]) 280 { 281 ds.rowStatus(leaveNum) = SPxBasisBase<R>::Desc::P_ON_LOWER; 282 theLRbound[leaveNum] = R(-infinity); 283 leavebound = -this->lhs(leaveNum); // slack !! 284 leaveMax = R(-infinity); 285 objChange += theURbound[leaveNum] * this->lhs(leaveNum); 286 } 287 else 288 { 289 ds.rowStatus(leaveNum) = SPxBasisBase<R>::Desc::P_ON_UPPER; 290 theURbound[leaveNum] = R(infinity); 291 leavebound = -this->rhs(leaveNum); // slack !! 292 leaveMax = R(infinity); 293 objChange += theLRbound[leaveNum] * this->rhs(leaveNum); 294 } 295 296 break; 297 298 default: 299 throw SPxInternalCodeException("XLEAVE02 This should never happen."); 300 } 301 302 SPxOut::debug(this, "DLEAVE51 SPxSolverBase<R>::getLeaveVals() : row {} : {} -> {} objChange: {}\n", 303 leaveNum, leaveStat, ds.rowStatus(leaveNum), objChange); 304 } 305 306 else 307 { 308 assert(leaveId.isSPxColId()); 309 leaveNum = this->number(SPxColId(leaveId)); 310 leaveStat = ds.colStatus(leaveNum); 311 312 assert(isBasic(leaveStat)); 313 314 switch(leaveStat) 315 { 316 case SPxBasisBase<R>::Desc::P_ON_UPPER : 317 assert(rep() == ROW); 318 ds.colStatus(leaveNum) = this->dualColStatus(leaveNum); 319 leavebound = 0; 320 leaveMax = R(-infinity); 321 break; 322 323 case SPxBasisBase<R>::Desc::P_ON_LOWER : 324 assert(rep() == ROW); 325 ds.colStatus(leaveNum) = this->dualColStatus(leaveNum); 326 leavebound = 0; 327 leaveMax = R(infinity); 328 break; 329 330 case SPxBasisBase<R>::Desc::P_FREE : 331 assert(rep() == ROW); 332 ds.colStatus(leaveNum) = this->dualColStatus(leaveNum); 333 334 if((*theFvec)[leaveIdx] < theLBbound[leaveIdx]) 335 { 336 leavebound = theLBbound[leaveIdx]; 337 leaveMax = R(-infinity); 338 } 339 else 340 { 341 leavebound = theUBbound[leaveIdx]; 342 leaveMax = R(infinity); 343 } 344 345 break; 346 347 case SPxBasisBase<R>::Desc::D_FREE : 348 assert(rep() == COLUMN); 349 assert(SPxLPBase<R>::upper(leaveNum) == SPxLPBase<R>::lower(leaveNum)); 350 ds.colStatus(leaveNum) = SPxBasisBase<R>::Desc::P_FIXED; 351 leavebound = SPxLPBase<R>::upper(leaveNum); 352 objChange += this->maxObj(leaveNum) * leavebound; 353 354 if((*theFvec)[leaveIdx] < theLBbound[leaveIdx]) 355 leaveMax = R(infinity); 356 else 357 leaveMax = R(-infinity); 358 359 break; 360 361 case SPxBasisBase<R>::Desc::D_ON_LOWER : 362 assert(rep() == COLUMN); 363 ds.colStatus(leaveNum) = SPxBasisBase<R>::Desc::P_ON_UPPER; 364 leavebound = SPxLPBase<R>::upper(leaveNum); 365 objChange += theUCbound[leaveNum] * leavebound; 366 leaveMax = R(-infinity); 367 break; 368 369 case SPxBasisBase<R>::Desc::D_ON_UPPER : 370 assert(rep() == COLUMN); 371 ds.colStatus(leaveNum) = SPxBasisBase<R>::Desc::P_ON_LOWER; 372 leavebound = SPxLPBase<R>::lower(leaveNum); 373 objChange += theLCbound[leaveNum] * leavebound; 374 leaveMax = R(infinity); 375 break; 376 377 case SPxBasisBase<R>::Desc::D_ON_BOTH : 378 assert(rep() == COLUMN); 379 380 if((*theFvec)[leaveIdx] > theUBbound[leaveIdx]) 381 { 382 leaveMax = R(-infinity); 383 leavebound = SPxLPBase<R>::upper(leaveNum); 384 objChange += theUCbound[leaveNum] * leavebound; 385 theLCbound[leaveNum] = R(-infinity); 386 ds.colStatus(leaveNum) = SPxBasisBase<R>::Desc::P_ON_UPPER; 387 } 388 else 389 { 390 leaveMax = R(infinity); 391 leavebound = SPxLPBase<R>::lower(leaveNum); 392 objChange += theLCbound[leaveNum] * leavebound; 393 theUCbound[leaveNum] = R(infinity); 394 ds.colStatus(leaveNum) = SPxBasisBase<R>::Desc::P_ON_LOWER; 395 } 396 397 break; 398 399 default: 400 throw SPxInternalCodeException("XLEAVE03 This should never happen."); 401 } 402 403 SPxOut::debug(this, "DLEAVE52 SPxSolverBase<R>::getLeaveVals() : col {} : {} -> {} objChange: {}\n", 404 leaveNum, leaveStat, ds.colStatus(leaveNum), objChange); 405 } 406 } 407 408 template <class R> 409 void SPxSolverBase<R>::getLeaveVals2( 410 R leaveMax, 411 SPxId enterId, 412 R& enterBound, 413 R& newUBbound, 414 R& newLBbound, 415 R& newCoPrhs, 416 StableSum<R>& objChange 417 ) 418 { 419 typename SPxBasisBase<R>::Desc& ds = this->desc(); 420 421 enterBound = 0; 422 423 if(enterId.isSPxRowId()) 424 { 425 int idx = this->number(SPxRowId(enterId)); 426 typename SPxBasisBase<R>::Desc::Status enterStat = ds.rowStatus(idx); 427 428 // coverity[switch_selector_expr_is_constant] 429 switch(enterStat) 430 { 431 case SPxBasisBase<R>::Desc::D_FREE : 432 assert(rep() == ROW); 433 434 if(thePvec->delta()[idx] * leaveMax < 0) 435 newCoPrhs = theLRbound[idx]; 436 else 437 newCoPrhs = theURbound[idx]; 438 439 newUBbound = R(infinity); 440 newLBbound = R(-infinity); 441 ds.rowStatus(idx) = SPxBasisBase<R>::Desc::P_FIXED; 442 break; 443 444 case SPxBasisBase<R>::Desc::D_ON_UPPER : 445 assert(rep() == ROW); 446 newUBbound = 0; 447 newLBbound = R(-infinity); 448 ds.rowStatus(idx) = SPxBasisBase<R>::Desc::P_ON_LOWER; 449 newCoPrhs = theLRbound[idx]; 450 break; 451 452 case SPxBasisBase<R>::Desc::D_ON_LOWER : 453 assert(rep() == ROW); 454 newUBbound = R(infinity); 455 newLBbound = 0; 456 ds.rowStatus(idx) = SPxBasisBase<R>::Desc::P_ON_UPPER; 457 newCoPrhs = theURbound[idx]; 458 break; 459 460 case SPxBasisBase<R>::Desc::D_ON_BOTH : 461 assert(rep() == ROW); 462 463 if(leaveMax * thePvec->delta()[idx] < 0) 464 { 465 newUBbound = 0; 466 newLBbound = R(-infinity); 467 ds.rowStatus(idx) = SPxBasisBase<R>::Desc::P_ON_LOWER; 468 newCoPrhs = theLRbound[idx]; 469 } 470 else 471 { 472 newUBbound = R(infinity); 473 newLBbound = 0; 474 ds.rowStatus(idx) = SPxBasisBase<R>::Desc::P_ON_UPPER; 475 newCoPrhs = theURbound[idx]; 476 } 477 478 break; 479 480 case SPxBasisBase<R>::Desc::P_ON_UPPER : 481 assert(rep() == COLUMN); 482 ds.rowStatus(idx) = this->dualRowStatus(idx); 483 484 if(this->lhs(idx) > R(-infinity)) 485 theURbound[idx] = theLRbound[idx]; 486 487 newCoPrhs = theLRbound[idx]; // slack !! 488 newUBbound = -this->lhs(idx); 489 newLBbound = -this->rhs(idx); 490 enterBound = -this->rhs(idx); 491 objChange -= newCoPrhs * this->rhs(idx); 492 break; 493 494 case SPxBasisBase<R>::Desc::P_ON_LOWER : 495 assert(rep() == COLUMN); 496 ds.rowStatus(idx) = this->dualRowStatus(idx); 497 498 if(this->rhs(idx) < R(infinity)) 499 theLRbound[idx] = theURbound[idx]; 500 501 newCoPrhs = theURbound[idx]; // slack !! 502 newLBbound = -this->rhs(idx); 503 newUBbound = -this->lhs(idx); 504 enterBound = -this->lhs(idx); 505 objChange -= newCoPrhs * this->lhs(idx); 506 break; 507 508 case SPxBasisBase<R>::Desc::P_FREE : 509 assert(rep() == COLUMN); 510 #if 1 511 throw SPxInternalCodeException("XLEAVE04 This should never happen."); 512 #else 513 SPX_MSG_ERROR(std::cerr << "ELEAVE53 ERROR: not yet debugged!" << std::endl;) 514 ds.rowStatus(idx) = this->dualRowStatus(idx); 515 newCoPrhs = theURbound[idx]; // slack !! 516 newUBbound = R(infinity); 517 newLBbound = R(-infinity); 518 enterBound = 0; 519 #endif 520 break; 521 522 case SPxBasisBase<R>::Desc::P_FIXED : 523 assert(rep() == COLUMN); 524 SPX_MSG_ERROR(std::cerr << "ELEAVE54 " 525 << "ERROR! Tried to put a fixed row variable into the basis: " 526 << "idx=" << idx 527 << ", lhs=" << this->lhs(idx) 528 << ", rhs=" << this->rhs(idx) << std::endl;) 529 throw SPxInternalCodeException("XLEAVE05 This should never happen."); 530 531 default: 532 throw SPxInternalCodeException("XLEAVE06 This should never happen."); 533 } 534 535 SPxOut::debug(this, "DLEAVE55 SPxSolverBase<R>::getLeaveVals2(): row {} : {} -> {} objChange: {}\n", 536 idx, enterStat, ds.rowStatus(idx), objChange); 537 } 538 539 else 540 { 541 assert(enterId.isSPxColId()); 542 int idx = this->number(SPxColId(enterId)); 543 typename SPxBasisBase<R>::Desc::Status enterStat = ds.colStatus(idx); 544 545 // coverity[switch_selector_expr_is_constant] 546 switch(enterStat) 547 { 548 case SPxBasisBase<R>::Desc::D_ON_UPPER : 549 assert(rep() == ROW); 550 newUBbound = 0; 551 newLBbound = R(-infinity); 552 ds.colStatus(idx) = SPxBasisBase<R>::Desc::P_ON_LOWER; 553 newCoPrhs = theLCbound[idx]; 554 break; 555 556 case SPxBasisBase<R>::Desc::D_ON_LOWER : 557 assert(rep() == ROW); 558 newUBbound = R(infinity); 559 newLBbound = 0; 560 ds.colStatus(idx) = SPxBasisBase<R>::Desc::P_ON_UPPER; 561 newCoPrhs = theUCbound[idx]; 562 break; 563 564 case SPxBasisBase<R>::Desc::D_FREE : 565 assert(rep() == ROW); 566 newUBbound = R(infinity); 567 newLBbound = R(-infinity); 568 newCoPrhs = theLCbound[idx]; 569 ds.colStatus(idx) = SPxBasisBase<R>::Desc::P_FIXED; 570 break; 571 572 case SPxBasisBase<R>::Desc::D_ON_BOTH : 573 assert(rep() == ROW); 574 575 if(leaveMax * theCoPvec->delta()[idx] < 0) 576 { 577 newUBbound = 0; 578 newLBbound = R(-infinity); 579 ds.colStatus(idx) = SPxBasisBase<R>::Desc::P_ON_LOWER; 580 newCoPrhs = theLCbound[idx]; 581 } 582 else 583 { 584 newUBbound = R(infinity); 585 newLBbound = 0; 586 ds.colStatus(idx) = SPxBasisBase<R>::Desc::P_ON_UPPER; 587 newCoPrhs = theUCbound[idx]; 588 } 589 590 break; 591 592 case SPxBasisBase<R>::Desc::P_ON_UPPER : 593 assert(rep() == COLUMN); 594 ds.colStatus(idx) = this->dualColStatus(idx); 595 596 if(SPxLPBase<R>::lower(idx) > R(-infinity)) 597 theLCbound[idx] = theUCbound[idx]; 598 599 newCoPrhs = theUCbound[idx]; 600 newUBbound = SPxLPBase<R>::upper(idx); 601 newLBbound = SPxLPBase<R>::lower(idx); 602 enterBound = SPxLPBase<R>::upper(idx); 603 objChange -= newCoPrhs * enterBound; 604 break; 605 606 case SPxBasisBase<R>::Desc::P_ON_LOWER : 607 assert(rep() == COLUMN); 608 ds.colStatus(idx) = this->dualColStatus(idx); 609 610 if(SPxLPBase<R>::upper(idx) < R(infinity)) 611 theUCbound[idx] = theLCbound[idx]; 612 613 newCoPrhs = theLCbound[idx]; 614 newUBbound = SPxLPBase<R>::upper(idx); 615 newLBbound = SPxLPBase<R>::lower(idx); 616 enterBound = SPxLPBase<R>::lower(idx); 617 objChange -= newCoPrhs * enterBound; 618 break; 619 620 case SPxBasisBase<R>::Desc::P_FREE : 621 assert(rep() == COLUMN); 622 ds.colStatus(idx) = this->dualColStatus(idx); 623 624 if(thePvec->delta()[idx] * leaveMax > 0) 625 newCoPrhs = theUCbound[idx]; 626 else 627 newCoPrhs = theLCbound[idx]; 628 629 newUBbound = SPxLPBase<R>::upper(idx); 630 newLBbound = SPxLPBase<R>::lower(idx); 631 enterBound = 0; 632 break; 633 634 case SPxBasisBase<R>::Desc::P_FIXED : 635 assert(rep() == COLUMN); 636 SPX_MSG_ERROR(std::cerr << "ELEAVE56 " 637 << "ERROR! Tried to put a fixed column variable into the basis. " 638 << "idx=" << idx 639 << ", lower=" << this->lower(idx) 640 << ", upper=" << this->upper(idx) << std::endl;) 641 throw SPxInternalCodeException("XLEAVE07 This should never happen."); 642 643 default: 644 throw SPxInternalCodeException("XLEAVE08 This should never happen."); 645 } 646 647 SPxOut::debug(this, "DLEAVE57 SPxSolverBase<R>::getLeaveVals2(): col {} : {} -> {} objChange: {}\n", 648 idx, enterStat, ds.colStatus(idx), objChange); 649 } 650 651 } 652 653 template <class R> 654 void SPxSolverBase<R>::rejectLeave( 655 int leaveNum, 656 SPxId leaveId, 657 typename SPxBasisBase<R>::Desc::Status leaveStat, 658 const SVectorBase<R>* //newVec 659 ) 660 { 661 typename SPxBasisBase<R>::Desc& ds = this->desc(); 662 663 if(leaveId.isSPxRowId()) 664 { 665 SPxOut::debug(this, "DLEAVE58 rejectLeave() : row {}: {} -> {}\n", leaveNum, 666 ds.rowStatus(leaveNum), leaveStat); 667 668 if(leaveStat == SPxBasisBase<R>::Desc::D_ON_BOTH) 669 { 670 if(ds.rowStatus(leaveNum) == SPxBasisBase<R>::Desc::P_ON_LOWER) 671 theLRbound[leaveNum] = theURbound[leaveNum]; 672 else 673 theURbound[leaveNum] = theLRbound[leaveNum]; 674 } 675 676 ds.rowStatus(leaveNum) = leaveStat; 677 } 678 else 679 { 680 SPxOut::debug(this, "DLEAVE59 rejectLeave() : col {}: {} -> {}\n", leaveNum, 681 ds.colStatus(leaveNum), leaveStat); 682 683 if(leaveStat == SPxBasisBase<R>::Desc::D_ON_BOTH) 684 { 685 if(ds.colStatus(leaveNum) == SPxBasisBase<R>::Desc::P_ON_UPPER) 686 theLCbound[leaveNum] = theUCbound[leaveNum]; 687 else 688 theUCbound[leaveNum] = theLCbound[leaveNum]; 689 } 690 691 ds.colStatus(leaveNum) = leaveStat; 692 } 693 } 694 695 696 template <class R> 697 void SPxSolverBase<R>::computePrimalray4Row(R direction) 698 { 699 R sign = (direction > 0 ? 1.0 : -1.0); 700 701 primalRay.clear(); 702 primalRay.setMax(coPvec().delta().size()); 703 704 for(int i = 0; i < coPvec().delta().size(); ++i) 705 primalRay.add(coPvec().delta().index(i), sign * coPvec().delta().value(i)); 706 } 707 708 template <class R> 709 void SPxSolverBase<R>::computeDualfarkas4Col(R direction) 710 { 711 R sign = (direction > 0 ? -1.0 : 1.0); 712 713 dualFarkas.clear(); 714 dualFarkas.setMax(coPvec().delta().size()); 715 716 for(int i = 0; i < coPvec().delta().size(); ++i) 717 dualFarkas.add(coPvec().delta().index(i), sign * coPvec().delta().value(i)); 718 } 719 720 template <class R> 721 bool SPxSolverBase<R>::leave(int leaveIdx, bool polish) 722 { 723 assert(leaveIdx < dim() && leaveIdx >= 0); 724 assert(type() == LEAVE); 725 assert(initialized); 726 727 bool instable = instableLeave; 728 assert(!instable || instableLeaveNum >= 0); 729 730 /* 731 Before performing the actual basis update, we must determine, how this 732 is to be accomplished. 733 When using steepest edge pricing this solve is already performed by the pricer 734 */ 735 if(theCoPvec->delta().isSetup() && theCoPvec->delta().size() == 0) 736 { 737 this->coSolve(theCoPvec->delta(), unitVecs[leaveIdx]); 738 } 739 740 #ifdef ENABLE_ADDITIONAL_CHECKS 741 else 742 { 743 SSVectorBase<R> tmp(dim(), this->tolerances()); 744 tmp.clear(); 745 this->coSolve(tmp, unitVecs[leaveIdx]); 746 tmp -= theCoPvec->delta(); 747 748 if(tmp.length() > leavetol()) 749 { 750 // This happens very frequently and does usually not hurt, so print 751 // these warnings only with verbose level INFO2 and higher. 752 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "WLEAVE60 iteration=" << basis().iteration() 753 << ": coPvec.delta error = " << tmp.length() 754 << std::endl;) 755 } 756 } 757 758 #endif // ENABLE_ADDITIONAL_CHECKS 759 760 setupPupdate(); 761 762 assert(thePvec->isConsistent()); 763 assert(theCoPvec->isConsistent()); 764 765 typename SPxBasisBase<R>::Desc::Status leaveStat; // status of leaving var 766 SPxId leaveId; // id of leaving var 767 SPxId none; // invalid id used if leave fails 768 R leaveMax; // maximium lambda of leaving var 769 R leavebound; // current fVec value of leaving var 770 int leaveNum; // number of leaveId in bounds 771 StableSum<R> objChange; // amount of change in the objective function 772 773 getLeaveVals(leaveIdx, leaveStat, leaveId, leaveMax, leavebound, leaveNum, objChange); 774 775 if(!polish && m_numCycle > m_maxCycle) 776 { 777 if(leaveMax > 0) 778 perturbMaxLeave(); 779 else 780 perturbMinLeave(); 781 782 //@ m_numCycle /= 2; 783 // perturbation invalidates the currently stored nonbasic value 784 forceRecompNonbasicValue(); 785 } 786 787 //@ testBounds(); 788 789 R enterVal = leaveMax; 790 boundflips = 0; 791 R oldShift = theShift; 792 SPxId enterId = theratiotester->selectEnter(enterVal, leaveIdx, polish); 793 794 if(NE(theShift, oldShift, this->tolerances()->epsilon())) 795 { 796 SPxOut::debug(this, 797 "DLEAVE71 trigger recomputation of nonbasic value due to shifts in ratiotest\n"); 798 forceRecompNonbasicValue(); 799 } 800 801 assert(!enterId.isValid() || !isBasic(enterId)); 802 803 instableLeaveNum = -1; 804 instableLeave = false; 805 806 /* 807 No variable could be selected to enter the basis and even the leaving 808 variable is unbounded. 809 */ 810 if(!enterId.isValid()) 811 { 812 /* the following line originally was below in "rejecting leave" case; 813 we need it in the unbounded/infeasible case, too, to have the 814 correct basis size */ 815 rejectLeave(leaveNum, leaveId, leaveStat); 816 this->change(-1, none, 0); 817 objChange = R(0.0); // the nonbasicValue is not supposed to be updated in this case 818 819 if(polish) 820 return false; 821 822 if(NE(enterVal, leaveMax, this->tolerances()->epsilon())) 823 { 824 SPxOut::debug(this, "DLEAVE61 rejecting leave A (leaveIdx={}, theCoTest={})\n", leaveIdx, 825 theCoTest[leaveIdx]); 826 827 /* In the LEAVE algorithm, when for a selected leaving variable we find only 828 an instable entering variable, then the basis change is not conducted. 829 Instead, we save the leaving variable's index in instableLeaveNum and scale 830 theCoTest[leaveIdx] down by some factor, hoping to find a different leaving 831 variable with a stable entering variable. 832 If this fails, however, and no more leaving variable is found, we have to 833 perform the instable basis change using instableLeaveNum. In this (and only 834 in this) case, the flag instableLeave is set to true. 835 836 enterVal != leaveMax is the case that selectEnter has found only an instable entering 837 variable. We store this leaving variable for later -- if we are not already in the 838 instable case: then we continue and conclude unboundedness/infeasibility */ 839 if(!instable) 840 { 841 instableLeaveNum = leaveIdx; 842 843 // Note: These changes do not survive a refactorization 844 instableLeaveVal = theCoTest[leaveIdx]; 845 theCoTest[leaveIdx] = instableLeaveVal / 10.0; 846 847 return true; 848 } 849 } 850 851 if(this->lastUpdate() > 1) 852 { 853 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "ILEAVE01 factorization triggered in " 854 << "leave() for feasibility test" << std::endl;) 855 856 try 857 { 858 factorize(); 859 } 860 catch(const SPxStatusException& E) 861 { 862 // don't exit immediately but handle the singularity correctly 863 assert(SPxBasisBase<R>::status() == SPxBasisBase<R>::SINGULAR); 864 SPX_MSG_INFO3((*this->spxout), 865 (*this->spxout) << "Caught exception in factorization: " << E.what() << 866 std::endl;) 867 } 868 869 /* after a factorization, the leaving column/row might not be infeasible or suboptimal anymore, hence we do 870 * not try to call leave(leaveIdx), but rather return to the main solving loop and call the pricer again 871 */ 872 return true; 873 } 874 875 /* do not exit with status infeasible or unbounded if there is only a very small violation */ 876 if(!recomputedVectors && spxAbs(enterVal) < leavetol()) 877 { 878 SPX_MSG_INFO3((*this->spxout), 879 (*this->spxout) << "ILEAVE11 clean up step to reduce numerical errors" << 880 std::endl;) 881 882 computeFrhs(); 883 SPxBasisBase<R>::solve(*theFvec, *theFrhs); 884 computeFtest(); 885 886 /* only do this once per solve */ 887 recomputedVectors = true; 888 889 return true; 890 } 891 892 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "ILEAVE02 unboundedness/infeasibility found " 893 << "in leave()" << std::endl;) 894 895 if(rep() != COLUMN) 896 { 897 computePrimalray4Row(enterVal); 898 setBasisStatus(SPxBasisBase<R>::UNBOUNDED); 899 } 900 else 901 { 902 computeDualfarkas4Col(enterVal); 903 setBasisStatus(SPxBasisBase<R>::INFEASIBLE); 904 } 905 906 return false; 907 } 908 else 909 { 910 /* 911 If an entering variable has been found, a regular basis update is to 912 be performed. 913 */ 914 if(enterId != this->baseId((leaveIdx))) 915 { 916 const SVectorBase<R>& newVector = *enterVector(enterId); 917 918 // update feasibility vectors 919 if(solveVector2 != NULL && solveVector3 != NULL) 920 { 921 assert(solveVector2->isConsistent()); 922 assert(solveVector2rhs->isSetup()); 923 assert(solveVector3->isConsistent()); 924 assert(solveVector3rhs->isSetup()); 925 assert(boundflips > 0); 926 SPxBasisBase<R>::solve4update(theFvec->delta(), 927 *solveVector2, 928 *solveVector3, 929 newVector, 930 *solveVector2rhs, 931 *solveVector3rhs); 932 933 // perform update of basic solution 934 primVec -= (*solveVector3); 935 SPxOut::debug(this, "ILBFRT02 breakpoints passed / bounds flipped = {}\n", boundflips); 936 totalboundflips += boundflips; 937 } 938 else if(solveVector2 != NULL) 939 { 940 assert(solveVector2->isConsistent()); 941 assert(solveVector2rhs->isSetup()); 942 943 SPxBasisBase<R>::solve4update(theFvec->delta(), 944 *solveVector2, 945 newVector, 946 *solveVector2rhs); 947 } 948 else if(solveVector3 != NULL) 949 { 950 assert(solveVector3->isConsistent()); 951 assert(solveVector3rhs->isSetup()); 952 assert(boundflips > 0); 953 SPxBasisBase<R>::solve4update(theFvec->delta(), 954 *solveVector3, 955 newVector, 956 *solveVector3rhs); 957 958 // perform update of basic solution 959 primVec -= (*solveVector3); 960 SPxOut::debug(this, "ILBFRT02 breakpoints passed / bounds flipped = {}\n", boundflips); 961 totalboundflips += boundflips; 962 } 963 else 964 SPxBasisBase<R>::solve4update(theFvec->delta(), newVector); 965 966 #ifdef ENABLE_ADDITIONAL_CHECKS 967 { 968 SSVectorBase<R> tmp(dim(), this->tolerances()); 969 SPxBasisBase<R>::solve(tmp, newVector); 970 tmp -= fVec().delta(); 971 972 if(tmp.length() > entertol()) 973 { 974 // This happens very frequently and does usually not hurt, so print 975 // these warnings only with verbose level INFO2 and higher. 976 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "WLEAVE62\t(" << tmp.length() << ")\n";) 977 } 978 } 979 #endif // ENABLE_ADDITIONAL_CHECKS 980 981 R reject_treshold = this->tolerances()->scaleAccordingToEpsilon(default_reject_leave_tol); 982 983 if(spxAbs(theFvec->delta()[leaveIdx]) < reject_treshold) 984 { 985 if(instable) 986 { 987 /* We are in the case that for all leaving variables only instable entering 988 variables were found: Thus, above we already accepted such an instable 989 entering variable. Now even this seems to be impossible, thus we conclude 990 unboundedness/infeasibility. */ 991 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "ILEAVE03 unboundedness/infeasibility found " 992 << "in leave()" << std::endl;) 993 994 rejectLeave(leaveNum, leaveId, leaveStat); 995 this->change(-1, none, 0); 996 objChange = R(0.0); // the nonbasicValue is not supposed to be updated in this case 997 998 /**@todo if shift() is not zero we must not conclude unboundedness */ 999 if(rep() == ROW) 1000 { 1001 computePrimalray4Row(enterVal); 1002 setBasisStatus(SPxBasisBase<R>::UNBOUNDED); 1003 } 1004 else 1005 { 1006 computeDualfarkas4Col(enterVal); 1007 setBasisStatus(SPxBasisBase<R>::INFEASIBLE); 1008 } 1009 1010 return false; 1011 } 1012 else 1013 { 1014 theFvec->delta().clear(); 1015 rejectLeave(leaveNum, leaveId, leaveStat, &newVector); 1016 this->change(-1, none, 0); 1017 objChange = R(0.0); // the nonbasicValue is not supposed to be updated in this case 1018 1019 SPxOut::debug(this, "DLEAVE63 rejecting leave B (leaveIdx={}, theCoTest={})\n", leaveIdx, 1020 theCoTest[leaveIdx]); 1021 1022 // Note: These changes do not survive a refactorization 1023 theCoTest[leaveIdx] *= 0.01; 1024 1025 return true; 1026 } 1027 } 1028 1029 // process leaving variable 1030 if(leavebound > epsilon() || leavebound < -epsilon()) 1031 theFrhs->multAdd(-leavebound, this->baseVec(leaveIdx)); 1032 1033 // process entering variable 1034 R enterBound; 1035 R newUBbound; 1036 R newLBbound; 1037 R newCoPrhs; 1038 1039 try 1040 { 1041 getLeaveVals2(leaveMax, enterId, enterBound, newUBbound, newLBbound, newCoPrhs, objChange); 1042 } 1043 catch(const SPxException& F) 1044 { 1045 rejectLeave(leaveNum, leaveId, leaveStat); 1046 this->change(-1, none, 0); 1047 objChange = R(0.0); // the nonbasicValue is not supposed to be updated in this case 1048 throw F; 1049 } 1050 1051 theUBbound[leaveIdx] = newUBbound; 1052 theLBbound[leaveIdx] = newLBbound; 1053 (*theCoPrhs)[leaveIdx] = newCoPrhs; 1054 1055 if(enterBound > epsilon() || enterBound < -epsilon()) 1056 theFrhs->multAdd(enterBound, newVector); 1057 1058 // update pricing vectors 1059 theCoPvec->value() = enterVal; 1060 thePvec->value() = enterVal; 1061 1062 if(enterVal > epsilon() || enterVal < -epsilon()) 1063 doPupdate(); 1064 1065 // update feasibility vector 1066 theFvec->value() = -((*theFvec)[leaveIdx] - leavebound) 1067 / theFvec->delta()[leaveIdx]; 1068 theFvec->update(); 1069 (*theFvec)[leaveIdx] = enterBound - theFvec->value(); 1070 updateFtest(); 1071 1072 // update objective funtion value 1073 updateNonbasicValue(objChange); 1074 1075 // change basis matrix 1076 this->change(leaveIdx, enterId, &newVector, &(theFvec->delta())); 1077 } 1078 1079 /* 1080 No entering vector has been selected from the basis. However, if the 1081 shift amount for |coPvec| is bounded, we are in the case, that the 1082 entering variable is moved from one bound to its other, before any of 1083 the basis feasibility variables reaches their bound. This may only 1084 happen in primal/columnwise case with upper and lower bounds on 1085 variables. 1086 */ 1087 else 1088 { 1089 // @todo update obj function value here!!! 1090 assert(rep() == ROW); 1091 typename SPxBasisBase<R>::Desc& ds = this->desc(); 1092 1093 this->change(leaveIdx, none, 0); 1094 1095 if(leaveStat == SPxBasisBase<R>::Desc::P_ON_UPPER) 1096 { 1097 if(leaveId.isSPxRowId()) 1098 { 1099 ds.rowStatus(leaveNum) = SPxBasisBase<R>::Desc::P_ON_LOWER; 1100 (*theCoPrhs)[leaveIdx] = theLRbound[leaveNum]; 1101 } 1102 else 1103 { 1104 ds.colStatus(leaveNum) = SPxBasisBase<R>::Desc::P_ON_LOWER; 1105 (*theCoPrhs)[leaveIdx] = theLCbound[leaveNum]; 1106 } 1107 1108 theUBbound[leaveIdx] = 0; 1109 theLBbound[leaveIdx] = R(-infinity); 1110 } 1111 else 1112 { 1113 assert(leaveStat == SPxBasisBase<R>::Desc::P_ON_LOWER); 1114 1115 if(leaveId.isSPxRowId()) 1116 { 1117 ds.rowStatus(leaveNum) = SPxBasisBase<R>::Desc::P_ON_UPPER; 1118 (*theCoPrhs)[leaveIdx] = theURbound[leaveNum]; 1119 } 1120 else 1121 { 1122 ds.colStatus(leaveNum) = SPxBasisBase<R>::Desc::P_ON_UPPER; 1123 (*theCoPrhs)[leaveIdx] = theUCbound[leaveNum]; 1124 } 1125 1126 theUBbound[leaveIdx] = R(infinity); 1127 theLBbound[leaveIdx] = 0; 1128 } 1129 1130 // update copricing vector 1131 theCoPvec->value() = enterVal; 1132 thePvec->value() = enterVal; 1133 1134 if(enterVal > epsilon() || enterVal < -epsilon()) 1135 doPupdate(); 1136 1137 // update feasibility vectors 1138 theFvec->value() = 0; 1139 assert(theCoTest[leaveIdx] < 0.0); 1140 m_pricingViol += theCoTest[leaveIdx]; 1141 theCoTest[leaveIdx] *= -1; 1142 } 1143 1144 if((leaveMax > entertol() && enterVal <= entertol()) || (leaveMax < -entertol() 1145 && enterVal >= -entertol())) 1146 { 1147 if((theUBbound[leaveIdx] < R(infinity) || theLBbound[leaveIdx] > R(-infinity)) 1148 && leaveStat != SPxBasisBase<R>::Desc::P_FREE 1149 && leaveStat != SPxBasisBase<R>::Desc::D_FREE) 1150 { 1151 m_numCycle++; 1152 leaveCycles++; 1153 } 1154 } 1155 else 1156 m_numCycle /= 2; 1157 1158 #ifdef ENABLE_ADDITIONAL_CHECKS 1159 { 1160 VectorBase<R> tmp = fVec(); 1161 this->multBaseWith(tmp); 1162 tmp -= fRhs(); 1163 1164 if(tmp.length() > entertol()) 1165 { 1166 // This happens very frequently and does usually not hurt, so print 1167 // these warnings only with verbose level INFO2 and higher. 1168 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "WLEAVE64\t" << basis().iteration() 1169 << ": fVec error = " << tmp.length() << std::endl;) 1170 SPxBasisBase<R>::solve(tmp, fRhs()); 1171 tmp -= fVec(); 1172 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "WLEAVE65\t(" << tmp.length() << ")\n";) 1173 } 1174 } 1175 #endif // ENABLE_ADDITIONAL_CHECKS 1176 1177 return true; 1178 } 1179 } 1180 } // namespace soplex 1181