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 "soplex/spxdefines.h" 27 #include "soplex/sorter.h" 28 #include "soplex/spxsolver.h" 29 #include "soplex/spxout.h" 30 #include "soplex/spxid.h" 31 32 namespace soplex 33 { 34 35 #define SOPLEX_LOWSTAB 1e-10 36 #define SOPLEX_MAX_RELAX_COUNT 2 37 #define SOPLEX_LONGSTEP_FREQ 100 38 39 40 /** perform necessary bound flips to restore dual feasibility */ 41 template <class R> 42 void SPxBoundFlippingRT<R>::flipAndUpdate( 43 int& nflips /**< number of bounds that should be flipped */ 44 ) 45 { 46 assert(nflips > 0); 47 48 // number of bound flips that are not performed 49 int skipped; 50 51 updPrimRhs.setup(); 52 updPrimRhs.reDim(this->thesolver->dim()); 53 updPrimVec.reDim(this->thesolver->dim()); 54 updPrimRhs.clear(); 55 updPrimVec.clear(); 56 57 skipped = 0; 58 59 for(int i = 0; i < nflips; ++i) 60 { 61 int idx; 62 idx = breakpoints[i].idx; 63 64 if(idx < 0) 65 { 66 ++skipped; 67 continue; 68 } 69 70 R range; 71 R upper; 72 R lower; 73 R objChange = 0.0; 74 typename SPxBasisBase<R>::Desc::Status stat; 75 typename SPxBasisBase<R>::Desc& ds = this->thesolver->basis().desc(); 76 77 range = 0; 78 79 if(breakpoints[i].src == PVEC) 80 { 81 assert(this->thesolver->rep() == SPxSolverBase<R>::COLUMN); 82 stat = ds.status(idx); 83 upper = this->thesolver->upper(idx); 84 lower = this->thesolver->lower(idx); 85 86 switch(stat) 87 { 88 case SPxBasisBase<R>::Desc::P_ON_UPPER : 89 ds.status(idx) = SPxBasisBase<R>::Desc::P_ON_LOWER; 90 range = lower - upper; 91 assert((*this->thesolver->theLbound)[idx] == R(-infinity)); 92 (*this->thesolver->theLbound)[idx] = (*this->thesolver->theUbound)[idx]; 93 (*this->thesolver->theUbound)[idx] = R(infinity); 94 objChange = range * (*this->thesolver->theLbound)[idx]; 95 break; 96 97 case SPxBasisBase<R>::Desc::P_ON_LOWER : 98 ds.status(idx) = SPxBasisBase<R>::Desc::P_ON_UPPER; 99 range = upper - lower; 100 assert((*this->thesolver->theUbound)[idx] == R(infinity)); 101 (*this->thesolver->theUbound)[idx] = (*this->thesolver->theLbound)[idx]; 102 (*this->thesolver->theLbound)[idx] = R(-infinity); 103 objChange = range * (*this->thesolver->theUbound)[idx]; 104 break; 105 106 default : 107 ++skipped; 108 SPX_MSG_WARNING((*this->thesolver->spxout), 109 (*this->thesolver->spxout) << "PVEC unexpected status: " << static_cast<int>(stat) 110 << " index: " << idx 111 << " val: " << this->thesolver->pVec()[idx] 112 << " upd: " << this->thesolver->pVec().delta()[idx] 113 << " lower: " << lower 114 << " upper: " << upper 115 << " bp.val: " << breakpoints[i].val 116 << std::endl;) 117 } 118 119 SPxOut::debug(this, 120 "PVEC flipped from: {} index: {} val: {} upd: {} lower: {} upper: {} bp.val: {} UCbound: {} LCbound: {}\n", 121 stat, idx, this->thesolver->pVec()[idx], this->thesolver->pVec().delta()[idx], lower, upper, 122 breakpoints[i].val, this->thesolver->theUCbound[idx], this->thesolver->theLCbound[idx]); 123 assert(spxAbs(range) < 1e20); 124 updPrimRhs.multAdd(range, this->thesolver->vector(idx)); 125 126 if(objChange != 0.0) 127 this->thesolver->updateNonbasicValue(objChange); 128 } 129 else if(breakpoints[i].src == COPVEC) 130 { 131 assert(this->thesolver->rep() == SPxSolverBase<R>::COLUMN); 132 stat = ds.coStatus(idx); 133 upper = this->thesolver->rhs(idx); 134 lower = this->thesolver->lhs(idx); 135 136 switch(stat) 137 { 138 case SPxBasisBase<R>::Desc::P_ON_UPPER : 139 ds.coStatus(idx) = SPxBasisBase<R>::Desc::P_ON_LOWER; 140 range = lower - upper; 141 assert((*this->thesolver->theCoUbound)[idx] == R(infinity)); 142 (*this->thesolver->theCoUbound)[idx] = -(*this->thesolver->theCoLbound)[idx]; 143 (*this->thesolver->theCoLbound)[idx] = R(-infinity); 144 objChange = range * (*this->thesolver->theCoUbound)[idx]; 145 break; 146 147 case SPxBasisBase<R>::Desc::P_ON_LOWER : 148 ds.coStatus(idx) = SPxBasisBase<R>::Desc::P_ON_UPPER; 149 range = upper - lower; 150 assert((*this->thesolver->theCoLbound)[idx] == R(-infinity)); 151 (*this->thesolver->theCoLbound)[idx] = -(*this->thesolver->theCoUbound)[idx]; 152 (*this->thesolver->theCoUbound)[idx] = R(infinity); 153 objChange = range * (*this->thesolver->theCoLbound)[idx]; 154 break; 155 156 default : 157 ++skipped; 158 SPX_MSG_WARNING((*this->thesolver->spxout), 159 (*this->thesolver->spxout) << "COPVEC unexpected status: " << static_cast<int>(stat) 160 << " index: " << idx 161 << " val: " << this->thesolver->coPvec()[idx] 162 << " upd: " << this->thesolver->coPvec().delta()[idx] 163 << " lower: " << lower 164 << " upper: " << upper 165 << " bp.val: " << breakpoints[i].val 166 << std::endl;) 167 } 168 169 SPxOut::debug(this, 170 "COPVEC flipped from: {} index: {} val: {} upd: {} lower: {} upper: {} bp.val: {} URbound: {} LRbound: {}\n", 171 stat, idx, this->thesolver->coPvec()[idx], this->thesolver->coPvec().delta()[idx], lower, upper, 172 breakpoints[i].val, this->thesolver->theURbound[idx], this->thesolver->theLRbound[idx]); 173 assert(spxAbs(range) < 1e20); 174 updPrimRhs.setValue(idx, updPrimRhs[idx] - range); 175 176 if(objChange != 0.0) 177 this->thesolver->updateNonbasicValue(objChange); 178 } 179 else if(breakpoints[i].src == FVEC) 180 { 181 assert(this->thesolver->rep() == SPxSolverBase<R>::ROW); 182 SPxId baseId = this->thesolver->basis().baseId(idx); 183 int IdNumber; 184 185 if(baseId.isSPxRowId()) 186 { 187 IdNumber = this->thesolver->number(SPxRowId(baseId)); 188 stat = ds.rowStatus(IdNumber); 189 upper = this->thesolver->rhs(IdNumber); 190 lower = this->thesolver->lhs(IdNumber); 191 192 switch(stat) 193 { 194 case SPxBasisBase<R>::Desc::P_ON_UPPER : 195 ds.rowStatus(IdNumber) = SPxBasisBase<R>::Desc::P_ON_LOWER; 196 range = upper - lower; 197 assert(this->thesolver->theUBbound[idx] == R(infinity)); 198 this->thesolver->theUBbound[idx] = -this->thesolver->theLBbound[idx]; 199 this->thesolver->theLBbound[idx] = R(-infinity); 200 break; 201 202 case SPxBasisBase<R>::Desc::P_ON_LOWER : 203 ds.rowStatus(IdNumber) = SPxBasisBase<R>::Desc::P_ON_UPPER; 204 range = lower - upper; 205 assert(this->thesolver->theLBbound[idx] == R(-infinity)); 206 this->thesolver->theLBbound[idx] = -this->thesolver->theUBbound[idx]; 207 this->thesolver->theUBbound[idx] = R(infinity); 208 break; 209 210 default : 211 ++skipped; 212 SPX_MSG_WARNING((*this->thesolver->spxout), 213 (*this->thesolver->spxout) << "unexpected basis status: " << static_cast<int>(stat) 214 << " index: " << idx 215 << " val: " << this->thesolver->fVec()[idx] 216 << " upd: " << this->thesolver->fVec().delta()[idx] 217 << " lower: " << lower 218 << " upper: " << upper 219 << " bp.val: " << breakpoints[i].val 220 << std::endl;) 221 } 222 } 223 else 224 { 225 assert(baseId.isSPxColId()); 226 IdNumber = this->thesolver->number(SPxColId(baseId)); 227 stat = ds.colStatus(IdNumber); 228 upper = this->thesolver->upper(IdNumber); 229 lower = this->thesolver->lower(IdNumber); 230 231 switch(stat) 232 { 233 case SPxBasisBase<R>::Desc::P_ON_UPPER : 234 ds.colStatus(IdNumber) = SPxBasisBase<R>::Desc::P_ON_LOWER; 235 range = upper - lower; 236 assert(this->thesolver->theUBbound[idx] == R(infinity)); 237 this->thesolver->theUBbound[idx] = -this->thesolver->theLBbound[idx]; 238 this->thesolver->theLBbound[idx] = R(-infinity); 239 break; 240 241 case SPxBasisBase<R>::Desc::P_ON_LOWER : 242 ds.colStatus(IdNumber) = SPxBasisBase<R>::Desc::P_ON_UPPER; 243 range = lower - upper; 244 assert(this->thesolver->theLBbound[idx] == R(-infinity)); 245 this->thesolver->theLBbound[idx] = -this->thesolver->theUBbound[idx]; 246 this->thesolver->theUBbound[idx] = R(infinity); 247 break; 248 249 default : 250 ++skipped; 251 SPX_MSG_WARNING((*this->thesolver->spxout), 252 (*this->thesolver->spxout) << "FVEC unexpected status: " << static_cast<int>(stat) 253 << " index: " << idx 254 << " val: " << this->thesolver->fVec()[idx] 255 << " upd: " << this->thesolver->fVec().delta()[idx] 256 << " lower: " << lower 257 << " upper: " << upper 258 << " bp.val: " << breakpoints[i].val 259 << std::endl;) 260 } 261 } 262 263 SPxOut::debug(this, 264 "basic row/col flipped from: {} index: {} val: {} upd: {} lower: {} upper: {} bp.val: {} \n", 265 stat, idx, this->thesolver->fVec()[idx], this->thesolver->fVec().delta()[idx], lower, upper, 266 breakpoints[i].val); 267 assert(spxAbs(range) < 1e20); 268 assert(updPrimRhs[idx] == 0); 269 updPrimRhs.add(idx, range); 270 } 271 } 272 273 nflips -= skipped; 274 275 if(nflips > 0) 276 { 277 if(this->thesolver->rep() == SPxSolverBase<R>::ROW) 278 { 279 assert(this->m_type == SPxSolverBase<R>::ENTER); 280 (*this->thesolver->theCoPrhs) -= updPrimRhs; 281 this->thesolver->setup4coSolve2(&updPrimVec, &updPrimRhs); 282 } 283 else 284 { 285 assert(this->thesolver->rep() == SPxSolverBase<R>::COLUMN); 286 assert(this->m_type == SPxSolverBase<R>::LEAVE); 287 (*this->thesolver->theFrhs) -= updPrimRhs; 288 this->thesolver->setup4solve2(&updPrimVec, &updPrimRhs); 289 } 290 } 291 292 return; 293 } 294 295 /** store all available pivots/breakpoints in an array (positive pivot search direction) */ 296 template <class R> 297 void SPxBoundFlippingRT<R>::collectBreakpointsMax( 298 int& nBp, /**< number of found breakpoints so far */ 299 int& minIdx, /**< index to current minimal breakpoint */ 300 const int* idx, /**< pointer to indices of current VectorBase<R> */ 301 int nnz, /**< number of nonzeros in current VectorBase<R> */ 302 const R* upd, /**< pointer to update values of current VectorBase<R> */ 303 const R* vec, /**< pointer to values of current VectorBase<R> */ 304 const R* upp, /**< pointer to upper bound/rhs of current VectorBase<R> */ 305 const R* low, /**< pointer to lower bound/lhs of current VectorBase<R> */ 306 BreakpointSource src /**< type of VectorBase<R> (pVec, coPvec or fVec)*/ 307 ) 308 { 309 R minVal; 310 R curVal; 311 const int* last; 312 313 minVal = (nBp == 0) ? R(infinity) : breakpoints[minIdx].val; 314 315 last = idx + nnz; 316 317 for(; idx < last; ++idx) 318 { 319 int i = *idx; 320 R x = upd[i]; 321 322 if(x > this->epsilonZero()) 323 { 324 if(upp[i] < R(infinity)) 325 { 326 R y = upp[i] - vec[i]; 327 curVal = (y <= 0) ? this->fastDelta / x : (y + this->fastDelta) / x; 328 assert(curVal > 0); 329 330 breakpoints[nBp].idx = i; 331 breakpoints[nBp].src = src; 332 breakpoints[nBp].val = curVal; 333 334 if(curVal < minVal) 335 { 336 minVal = curVal; 337 minIdx = nBp; 338 } 339 340 nBp++; 341 } 342 } 343 else if(x < -this->epsilonZero()) 344 { 345 if(low[i] > R(-infinity)) 346 { 347 R y = low[i] - vec[i]; 348 curVal = (y >= 0) ? -this->fastDelta / x : (y - this->fastDelta) / x; 349 assert(curVal > 0); 350 351 breakpoints[nBp].idx = i; 352 breakpoints[nBp].src = src; 353 breakpoints[nBp].val = curVal; 354 355 if(curVal < minVal) 356 { 357 minVal = curVal; 358 minIdx = nBp; 359 } 360 361 nBp++; 362 } 363 } 364 365 if(nBp >= breakpoints.size()) 366 breakpoints.reSize(nBp * 2); 367 } 368 369 return; 370 } 371 372 /** store all available pivots/breakpoints in an array (negative pivot search direction) */ 373 template <class R> 374 void SPxBoundFlippingRT<R>::collectBreakpointsMin( 375 int& nBp, /**< number of found breakpoints so far */ 376 int& minIdx, /**< index to current minimal breakpoint */ 377 const int* idx, /**< pointer to indices of current VectorBase<R> */ 378 int nnz, /**< number of nonzeros in current VectorBase<R> */ 379 const R* upd, /**< pointer to update values of current VectorBase<R> */ 380 const R* vec, /**< pointer to values of current VectorBase<R> */ 381 const R* upp, /**< pointer to upper bound/rhs of current VectorBase<R> */ 382 const R* low, /**< pointer to lower bound/lhs of current VectorBase<R> */ 383 BreakpointSource src /**< type of VectorBase<R> (pVec, coPvec or fVec)*/ 384 ) 385 { 386 R minVal; 387 R curVal; 388 const int* last; 389 390 minVal = (nBp == 0) ? R(infinity) : breakpoints[minIdx].val; 391 392 last = idx + nnz; 393 394 for(; idx < last; ++idx) 395 { 396 int i = *idx; 397 R x = upd[i]; 398 399 if(x > this->epsilonZero()) 400 { 401 if(low[i] > R(-infinity)) 402 { 403 R y = low[i] - vec[i]; 404 405 curVal = (y >= 0) ? this->fastDelta / x : (this->fastDelta - y) / x; 406 assert(curVal > 0); 407 408 breakpoints[nBp].idx = i; 409 breakpoints[nBp].src = src; 410 breakpoints[nBp].val = curVal; 411 412 if(curVal < minVal) 413 { 414 minVal = curVal; 415 minIdx = nBp; 416 } 417 418 nBp++; 419 } 420 } 421 else if(x < -this->epsilonZero()) 422 { 423 if(upp[i] < R(infinity)) 424 { 425 R y = upp[i] - vec[i]; 426 curVal = (y <= 0) ? -this->fastDelta / x : -(y + this->fastDelta) / x; 427 assert(curVal > 0); 428 429 breakpoints[nBp].idx = i; 430 breakpoints[nBp].src = src; 431 breakpoints[nBp].val = curVal; 432 433 if(curVal < minVal) 434 { 435 minVal = curVal; 436 minIdx = nBp; 437 } 438 439 nBp++; 440 } 441 } 442 443 if(nBp >= breakpoints.size()) 444 breakpoints.reSize(nBp * 2); 445 } 446 447 return; 448 } 449 450 /** get values for entering index and perform shifts if necessary */ 451 template <class R> 452 bool SPxBoundFlippingRT<R>::getData( 453 R& val, 454 SPxId& enterId, 455 int idx, 456 R stab, 457 R degeneps, 458 const R* upd, 459 const R* vec, 460 const R* low, 461 const R* upp, 462 BreakpointSource src, 463 R max 464 ) 465 { 466 if(src == PVEC) 467 { 468 this->thesolver->pVec()[idx] = this->thesolver->vector(idx) * this->thesolver->coPvec(); 469 R x = upd[idx]; 470 471 // skip breakpoint if it is too small 472 if(spxAbs(x) < stab) 473 { 474 return false; 475 } 476 477 enterId = this->thesolver->id(idx); 478 val = (max * x > 0) ? upp[idx] : low[idx]; 479 val = (val - vec[idx]) / x; 480 481 if(upp[idx] == low[idx]) 482 { 483 val = 0.0; 484 485 if(vec[idx] > upp[idx]) 486 this->thesolver->theShift += vec[idx] - upp[idx]; 487 else 488 this->thesolver->theShift += low[idx] - vec[idx]; 489 490 this->thesolver->upBound()[idx] = this->thesolver->lpBound()[idx] = vec[idx]; 491 } 492 else if((max > 0 && val < -degeneps) || (max < 0 && val > degeneps)) 493 { 494 val = 0.0; 495 496 if(max * x > 0) 497 this->thesolver->shiftUPbound(idx, vec[idx]); 498 else 499 this->thesolver->shiftLPbound(idx, vec[idx]); 500 } 501 } 502 else // src == COPVEC 503 { 504 R x = upd[idx]; 505 506 if(spxAbs(x) < stab) 507 { 508 return false; 509 } 510 511 enterId = this->thesolver->coId(idx); 512 val = (max * x > 0.0) ? upp[idx] : low[idx]; 513 val = (val - vec[idx]) / x; 514 515 if(upp[idx] == low[idx]) 516 { 517 val = 0.0; 518 519 if(vec[idx] > upp[idx]) 520 this->thesolver->theShift += vec[idx] - upp[idx]; 521 else 522 this->thesolver->theShift += low[idx] - vec[idx]; 523 524 this->thesolver->ucBound()[idx] = this->thesolver->lcBound()[idx] = vec[idx]; 525 } 526 else if((max > 0 && val < -degeneps) || (max < 0 && val > degeneps)) 527 { 528 val = 0.0; 529 530 if(max * x > 0) 531 this->thesolver->shiftUCbound(idx, vec[idx]); 532 else 533 this->thesolver->shiftLCbound(idx, vec[idx]); 534 } 535 } 536 537 return true; 538 } 539 540 /** get values for leaving index and perform shifts if necessary */ 541 template <class R> 542 bool SPxBoundFlippingRT<R>::getData( 543 R& val, 544 int& leaveIdx, 545 int idx, 546 R stab, 547 R degeneps, 548 const R* upd, 549 const R* vec, 550 const R* low, 551 const R* upp, 552 BreakpointSource src, 553 R max 554 ) 555 { 556 assert(src == FVEC); 557 558 R x = upd[idx]; 559 560 // skip breakpoint if it is too small 561 if(spxAbs(x) < stab) 562 { 563 return false; 564 } 565 566 leaveIdx = idx; 567 val = (max * x > 0) ? upp[idx] : low[idx]; 568 val = (val - vec[idx]) / x; 569 570 if(upp[idx] == low[idx]) 571 { 572 val = 0.0; 573 this->thesolver->shiftLBbound(idx, vec[idx]); 574 this->thesolver->shiftUBbound(idx, vec[idx]); 575 } 576 else if((max > 0 && val < -degeneps) || (max < 0 && val > degeneps)) 577 { 578 val = 0.0; 579 580 if(this->thesolver->dualStatus(this->thesolver->baseId(idx)) != SPxBasisBase<R>::Desc::D_ON_BOTH) 581 { 582 if(max * x > 0) 583 this->thesolver->shiftUBbound(idx, vec[idx]); 584 else 585 this->thesolver->shiftLBbound(idx, vec[idx]); 586 } 587 } 588 589 return true; 590 } 591 592 /** determine entering row/column */ 593 template <class R> 594 SPxId SPxBoundFlippingRT<R>::selectEnter( 595 R& val, 596 int leaveIdx, 597 bool polish 598 ) 599 { 600 assert(this->m_type == SPxSolverBase<R>::LEAVE); 601 assert(this->thesolver->boundflips == 0); 602 603 // reset the history and try again to do some long steps 604 if(this->thesolver->leaveCount % SOPLEX_LONGSTEP_FREQ == 0) 605 { 606 SPxOut::debug(this, "DLBFRT06 resetting long step history\n"); 607 flipPotential = 1; 608 } 609 610 if(!enableBoundFlips || polish || this->thesolver->rep() == SPxSolverBase<R>::ROW 611 || flipPotential <= 0) 612 { 613 SPxOut::debug(this, "DLBFRT07 switching to fast ratio test\n"); 614 return SPxFastRT<R>::selectEnter(val, leaveIdx, polish); 615 } 616 617 const R* pvec = this->thesolver->pVec().get_const_ptr(); 618 const R* pupd = this->thesolver->pVec().delta().values(); 619 const int* pidx = this->thesolver->pVec().delta().indexMem(); 620 int pupdnnz = this->thesolver->pVec().delta().size(); 621 const R* lpb = this->thesolver->lpBound().get_const_ptr(); 622 const R* upb = this->thesolver->upBound().get_const_ptr(); 623 624 const R* cvec = this->thesolver->coPvec().get_const_ptr(); 625 const R* cupd = this->thesolver->coPvec().delta().values(); 626 const int* cidx = this->thesolver->coPvec().delta().indexMem(); 627 int cupdnnz = this->thesolver->coPvec().delta().size(); 628 const R* lcb = this->thesolver->lcBound().get_const_ptr(); 629 const R* ucb = this->thesolver->ucBound().get_const_ptr(); 630 631 R max; 632 633 // index in breakpoint array of minimal value (i.e. choice of normal RT) 634 int minIdx; 635 636 // temporary breakpoint data structure to make swaps possible 637 Breakpoint tmp; 638 639 // most stable pivot value in candidate set 640 R moststable; 641 642 // initialize invalid enterId 643 SPxId enterId; 644 645 // slope of objective function improvement 646 R slope; 647 648 // number of found breakpoints 649 int nBp; 650 651 // number of passed breakpoints 652 int npassedBp; 653 654 R degeneps; 655 R stab; 656 bool instable; 657 658 max = val; 659 val = 0.0; 660 moststable = 0.0; 661 nBp = 0; 662 minIdx = -1; 663 664 // get breakpoints and and determine the index of the minimal value 665 if(max > 0) 666 { 667 collectBreakpointsMax(nBp, minIdx, pidx, pupdnnz, pupd, pvec, upb, lpb, PVEC); 668 // coverity[negative_returns] 669 collectBreakpointsMax(nBp, minIdx, cidx, cupdnnz, cupd, cvec, ucb, lcb, COPVEC); 670 } 671 else 672 { 673 collectBreakpointsMin(nBp, minIdx, pidx, pupdnnz, pupd, pvec, upb, lpb, PVEC); 674 // coverity[negative_returns] 675 collectBreakpointsMin(nBp, minIdx, cidx, cupdnnz, cupd, cvec, ucb, lcb, COPVEC); 676 } 677 678 if(nBp == 0) 679 { 680 val = max; 681 return enterId; 682 } 683 684 assert(minIdx >= 0); 685 686 // swap smallest breakpoint to the front to skip the sorting phase if no bound flip is possible 687 tmp = breakpoints[minIdx]; 688 breakpoints[minIdx] = breakpoints[0]; 689 breakpoints[0] = tmp; 690 691 // get initial slope 692 slope = spxAbs(this->thesolver->fTest()[leaveIdx]); 693 694 if(slope == 0) 695 { 696 // this may only happen if SoPlex decides to make an instable pivot 697 assert(this->thesolver->instableLeaveNum >= 0); 698 // restore original slope 699 slope = spxAbs(this->thesolver->instableLeaveVal); 700 } 701 702 // set up structures for the quicksort implementation 703 BreakpointCompare compare; 704 compare.entry = breakpoints.get_const_ptr(); 705 706 // pointer to end of sorted part of breakpoints 707 int sorted = 0; 708 // minimum number of entries that are supposed to be sorted by partial sort 709 int sortsize = 4; 710 711 // get all skipable breakpoints 712 for(npassedBp = 0; npassedBp < nBp && slope > 0; ++npassedBp) 713 { 714 // sort breakpoints only partially to save time 715 if(npassedBp > sorted) 716 { 717 sorted = SPxQuicksortPart(breakpoints.get_ptr(), compare, sorted + 1, nBp, sortsize); 718 } 719 720 int i = breakpoints[npassedBp].idx; 721 722 // compute new slope 723 if(breakpoints[npassedBp].src == PVEC) 724 { 725 if(this->thesolver->isBasic(i)) 726 { 727 // mark basic indices 728 breakpoints[npassedBp].idx = -1; 729 this->thesolver->pVec().delta().clearIdx(i); 730 } 731 else 732 { 733 R absupd = spxAbs(pupd[i]); 734 slope -= (this->thesolver->upper(i) * absupd) - (this->thesolver->lower(i) * absupd); 735 736 // get most stable pivot 737 if(absupd > moststable) 738 moststable = absupd; 739 } 740 } 741 else 742 { 743 assert(breakpoints[npassedBp].src == COPVEC); 744 745 if(this->thesolver->isCoBasic(i)) 746 { 747 // mark basic indices 748 breakpoints[npassedBp].idx = -1; 749 this->thesolver->coPvec().delta().clearIdx(i); 750 } 751 else 752 { 753 R absupd = spxAbs(cupd[i]); 754 slope -= (this->thesolver->rhs(i) * absupd) - (this->thesolver->lhs(i) * absupd); 755 756 if(absupd > moststable) 757 moststable = absupd; 758 } 759 } 760 } 761 762 --npassedBp; 763 assert(npassedBp >= 0); 764 765 // check for unboundedness/infeasibility 766 if(slope > this->delta && npassedBp >= nBp - 1) 767 { 768 SPxOut::debug(this, "DLBFRT02 {}: unboundedness in ratio test\n", 769 this->thesolver->basis().iteration()); 770 flipPotential -= 0.5; 771 val = max; 772 return SPxFastRT<R>::selectEnter(val, leaveIdx); 773 } 774 775 SPxOut::debug(this, "DLBFRT01 {}: number of flip candidates: {}\n", 776 this->thesolver->basis().iteration(), npassedBp); 777 778 // try to get a more stable pivot by looking at those with similar step length 779 int stableBp; // index to walk over additional breakpoints (after slope change) 780 int bestBp = -1; // breakpoints index with best possible stability 781 R bestDelta = breakpoints[npassedBp].val; // best step length (after bound flips) 782 783 for(stableBp = npassedBp + 1; stableBp < nBp; ++stableBp) 784 { 785 R stableDelta = 0; 786 787 // get next breakpoints in increasing order 788 if(stableBp > sorted) 789 { 790 sorted = SPxQuicksortPart(breakpoints.get_ptr(), compare, sorted + 1, nBp, sortsize); 791 } 792 793 int idx = breakpoints[stableBp].idx; 794 795 if(breakpoints[stableBp].src == PVEC) 796 { 797 if(this->thesolver->isBasic(idx)) 798 { 799 // mark basic indices 800 breakpoints[stableBp].idx = -1; 801 this->thesolver->pVec().delta().clearIdx(idx); 802 continue; 803 } 804 805 R x = pupd[idx]; 806 807 if(spxAbs(x) > moststable) 808 { 809 this->thesolver->pVec()[idx] = this->thesolver->vector(idx) * this->thesolver->coPvec(); 810 stableDelta = (x > 0.0) ? upb[idx] : lpb[idx]; 811 stableDelta = (stableDelta - pvec[idx]) / x; 812 813 if(stableDelta <= bestDelta) 814 { 815 moststable = spxAbs(x); 816 bestBp = stableBp; 817 } 818 } 819 } 820 else 821 { 822 if(this->thesolver->isCoBasic(idx)) 823 { 824 // mark basic indices 825 breakpoints[stableBp].idx = -1; 826 this->thesolver->coPvec().delta().clearIdx(idx); 827 continue; 828 } 829 830 R x = cupd[idx]; 831 832 if(spxAbs(x) > moststable) 833 { 834 stableDelta = (x > 0.0) ? ucb[idx] : lcb[idx]; 835 stableDelta = (stableDelta - cvec[idx]) / x; 836 837 if(stableDelta <= bestDelta) 838 { 839 moststable = spxAbs(x); 840 bestBp = stableBp; 841 } 842 } 843 } 844 845 // stop searching if the step length is too big 846 if(stableDelta > this->delta + bestDelta) 847 break; 848 } 849 850 degeneps = this->fastDelta / moststable; /* as in SPxFastRT */ 851 // get stability requirements 852 instable = this->thesolver->instableLeave; 853 assert(!instable || this->thesolver->instableLeaveNum >= 0); 854 R lowstab = this->tolerances()->scaleAccordingToEpsilon(SOPLEX_LOWSTAB); 855 stab = instable ? lowstab : SPxFastRT<R>::minStability(moststable); 856 857 bool foundStable = false; 858 859 if(bestBp >= 0) 860 { 861 // found a more stable pivot 862 if(moststable > stab) 863 { 864 // stability requirements are satisfied 865 int idx = breakpoints[bestBp].idx; 866 assert(idx >= 0); 867 868 if(breakpoints[bestBp].src == PVEC) 869 foundStable = getData(val, enterId, idx, stab, degeneps, pupd, pvec, lpb, upb, PVEC, max); 870 else 871 foundStable = getData(val, enterId, idx, stab, degeneps, cupd, cvec, lcb, ucb, COPVEC, max); 872 } 873 } 874 875 else 876 { 877 // scan passed breakpoints from back to front and stop as soon as a good one is found 878 while(!foundStable && npassedBp >= 0) 879 { 880 int idx = breakpoints[npassedBp].idx; 881 882 // only look for non-basic variables 883 if(idx >= 0) 884 { 885 if(breakpoints[npassedBp].src == PVEC) 886 foundStable = getData(val, enterId, idx, stab, degeneps, pupd, pvec, lpb, upb, PVEC, max); 887 else 888 foundStable = getData(val, enterId, idx, stab, degeneps, cupd, cvec, lcb, ucb, COPVEC, max); 889 } 890 891 --npassedBp; 892 } 893 894 ++npassedBp; 895 } 896 897 if(!foundStable) 898 { 899 assert(!enterId.isValid()); 900 901 if(relax_count < SOPLEX_MAX_RELAX_COUNT) 902 { 903 SPxOut::debug(this, "DLBFRT04 {}: no valid enterId found - relaxing...\n", 904 this->thesolver->basis().iteration()); 905 this->relax(); 906 ++relax_count; 907 // restore original value 908 val = max; 909 // try again with relaxed delta 910 return SPxBoundFlippingRT<R>::selectEnter(val, leaveIdx); 911 } 912 else 913 { 914 SPxOut::debug(this, "DLBFRT05 {}: no valid enterId found - breaking...\n", 915 this->thesolver->basis().iteration()); 916 return enterId; 917 } 918 } 919 else 920 { 921 relax_count = 0; 922 this->tighten(); 923 } 924 925 // flip bounds of skipped breakpoints only if a nondegenerate step is to be performed 926 if(npassedBp > 0 && spxAbs(breakpoints[npassedBp].val) > this->fastDelta) 927 { 928 flipAndUpdate(npassedBp); 929 this->thesolver->boundflips = npassedBp; 930 931 if(npassedBp >= 10) 932 flipPotential = 1; 933 else 934 flipPotential -= 0.05; 935 } 936 else 937 { 938 this->thesolver->boundflips = 0; 939 flipPotential -= 0.1; 940 } 941 942 SPxOut::debug(this, "DLBFRT06 {}: selected Id: {} number of candidates: {}\n", 943 this->thesolver->basis().iteration(), enterId, nBp); 944 return enterId; 945 } 946 947 /** determine leaving row/column */ 948 template <class R> 949 int SPxBoundFlippingRT<R>::selectLeave( 950 R& val, 951 R enterTest, 952 bool polish 953 ) 954 { 955 assert(this->m_type == SPxSolverBase<R>::ENTER); 956 assert(this->thesolver->boundflips == 0); 957 958 // reset the history and try again to do some long steps 959 if(this->thesolver->enterCount % SOPLEX_LONGSTEP_FREQ == 0) 960 { 961 SPxOut::debug(this, "DEBFRT06 resetting long step history\n"); 962 flipPotential = 1; 963 } 964 965 if(polish || !enableBoundFlips || !enableRowBoundFlips 966 || this->thesolver->rep() == SPxSolverBase<R>::COLUMN || flipPotential <= 0) 967 { 968 SPxOut::debug(this, "DEBFRT07 switching to fast ratio test\n"); 969 return SPxFastRT<R>::selectLeave(val, enterTest, polish); 970 } 971 972 const R* vec = 973 this->thesolver->fVec().get_const_ptr(); /**< pointer to values of current VectorBase<R> */ 974 const R* upd = 975 this->thesolver->fVec().delta().values(); /**< pointer to update values of current VectorBase<R> */ 976 const int* idx = 977 this->thesolver->fVec().delta().indexMem(); /**< pointer to indices of current VectorBase<R> */ 978 int updnnz = 979 this->thesolver->fVec().delta().size(); /**< number of nonzeros in update VectorBase<R> */ 980 const R* lb = 981 this->thesolver->lbBound().get_const_ptr(); /**< pointer to lower bound/lhs of current VectorBase<R> */ 982 const R* ub = 983 this->thesolver->ubBound().get_const_ptr(); /**< pointer to upper bound/rhs of current VectorBase<R> */ 984 985 R max; 986 987 // index in breakpoint array of minimal value (i.e. choice of normal RT) 988 int minIdx; 989 990 // temporary breakpoint data structure to make swaps possible 991 Breakpoint tmp; 992 993 // most stable pivot value in candidate set 994 R moststable; 995 996 // initialize invalid leaving index 997 int leaveIdx = -1; 998 999 // slope of objective function improvement 1000 R slope; 1001 1002 // number of found breakpoints 1003 int nBp; 1004 1005 // number of passed breakpoints 1006 int npassedBp; 1007 1008 R degeneps; 1009 R stab; 1010 bool instable; 1011 1012 max = val; 1013 val = 0.0; 1014 moststable = 0.0; 1015 nBp = 0; 1016 minIdx = -1; 1017 1018 assert(this->thesolver->fVec().delta().isSetup()); 1019 1020 // get breakpoints and and determine the index of the minimal value 1021 if(max > 0) 1022 { 1023 collectBreakpointsMax(nBp, minIdx, idx, updnnz, upd, vec, ub, lb, FVEC); 1024 } 1025 else 1026 { 1027 collectBreakpointsMin(nBp, minIdx, idx, updnnz, upd, vec, ub, lb, FVEC); 1028 } 1029 1030 // return -1 if no BP was found 1031 if(nBp == 0) 1032 { 1033 val = max; 1034 return leaveIdx; 1035 } 1036 1037 assert(minIdx >= 0); 1038 1039 // swap smallest breakpoint to the front to skip the sorting phase if no bound flip is possible 1040 tmp = breakpoints[minIdx]; 1041 breakpoints[minIdx] = breakpoints[0]; 1042 breakpoints[0] = tmp; 1043 1044 // get initial slope 1045 slope = spxAbs(enterTest); 1046 1047 if(slope == 0) 1048 { 1049 // this may only happen if SoPlex decides to make an instable pivot 1050 assert(this->thesolver->instableEnterId.isValid()); 1051 // restore original slope 1052 slope = this->thesolver->instableEnterVal; 1053 } 1054 1055 // set up structures for the quicksort implementation 1056 BreakpointCompare compare; 1057 compare.entry = breakpoints.get_const_ptr(); 1058 1059 // pointer to end of sorted part of breakpoints 1060 int sorted = 0; 1061 // minimum number of entries that are supposed to be sorted by partial sort 1062 int sortsize = 4; 1063 1064 // get all skipable breakpoints 1065 for(npassedBp = 0; npassedBp < nBp && slope > 0; ++npassedBp) 1066 { 1067 // sort breakpoints only partially to save time 1068 if(npassedBp > sorted) 1069 { 1070 sorted = SPxQuicksortPart(breakpoints.get_ptr(), compare, sorted + 1, nBp, sortsize); 1071 } 1072 1073 assert(breakpoints[npassedBp].src == FVEC); 1074 int breakpointidx = breakpoints[npassedBp].idx; 1075 // compute new slope 1076 R upper; 1077 R lower; 1078 R absupd = spxAbs(upd[breakpointidx]); 1079 SPxId baseId = this->thesolver->baseId(breakpointidx); 1080 int i = this->thesolver->number(baseId); 1081 1082 if(baseId.isSPxColId()) 1083 { 1084 upper = this->thesolver->upper(i); 1085 lower = this->thesolver->lower(i); 1086 } 1087 else 1088 { 1089 assert(baseId.isSPxRowId()); 1090 upper = this->thesolver->rhs(i); 1091 lower = this->thesolver->lhs(i); 1092 } 1093 1094 slope -= (upper * absupd) - (lower * absupd); 1095 1096 // get most stable pivot 1097 if(absupd > moststable) 1098 moststable = absupd; 1099 } 1100 1101 --npassedBp; 1102 assert(npassedBp >= 0); 1103 1104 // check for unboundedness/infeasibility 1105 if(slope > this->delta && npassedBp >= nBp - 1) 1106 { 1107 SPxOut::debug(this, "DEBFRT02 {}: unboundedness in ratio test\n", 1108 this->thesolver->basis().iteration()); 1109 flipPotential -= 0.5; 1110 val = max; 1111 return SPxFastRT<R>::selectLeave(val, enterTest); 1112 } 1113 1114 SPxOut::debug(this, "DEBFRT01 {}: number of flip candidates: {}\n", 1115 this->thesolver->basis().iteration(), npassedBp); 1116 1117 // try to get a more stable pivot by looking at those with similar step length 1118 int stableBp; // index to walk over additional breakpoints (after slope change) 1119 int bestBp = -1; // breakpoints index with best possible stability 1120 R bestDelta = breakpoints[npassedBp].val; // best step length (after bound flips) 1121 1122 for(stableBp = npassedBp + 1; stableBp < nBp; ++stableBp) 1123 { 1124 R stableDelta = 0; 1125 1126 // get next breakpoints in increasing order 1127 if(stableBp > sorted) 1128 { 1129 sorted = SPxQuicksortPart(breakpoints.get_ptr(), compare, sorted + 1, nBp, sortsize); 1130 } 1131 1132 int breakpointidx = breakpoints[stableBp].idx; 1133 assert(breakpoints[stableBp].src == FVEC); 1134 R x = upd[breakpointidx]; 1135 1136 if(spxAbs(x) > moststable) 1137 { 1138 stableDelta = (x > 0.0) ? ub[breakpointidx] : lb[breakpointidx]; 1139 stableDelta = (stableDelta - vec[breakpointidx]) / x; 1140 1141 if(stableDelta <= bestDelta) 1142 { 1143 moststable = spxAbs(x); 1144 bestBp = stableBp; 1145 } 1146 } 1147 // stop searching if the step length is too big 1148 else if(stableDelta > this->delta + bestDelta) 1149 break; 1150 } 1151 1152 degeneps = this->fastDelta / moststable; /* as in SPxFastRT */ 1153 // get stability requirements 1154 instable = this->thesolver->instableEnter; 1155 assert(!instable || this->thesolver->instableEnterId.isValid()); 1156 R lowstab = this->tolerances()->scaleAccordingToEpsilon(SOPLEX_LOWSTAB); 1157 stab = instable ? lowstab : SPxFastRT<R>::minStability(moststable); 1158 1159 bool foundStable = false; 1160 1161 if(bestBp >= 0) 1162 { 1163 // found a more stable pivot 1164 if(moststable > stab) 1165 { 1166 // stability requirements are satisfied 1167 int breakpointidx = breakpoints[bestBp].idx; 1168 assert(breakpointidx >= 0); 1169 foundStable = getData(val, leaveIdx, breakpointidx, moststable, degeneps, upd, vec, lb, ub, FVEC, 1170 max); 1171 } 1172 } 1173 1174 else 1175 { 1176 // scan passed breakpoints from back to front and stop as soon as a good one is found 1177 while(!foundStable && npassedBp >= 0) 1178 { 1179 int breakpointidx = breakpoints[npassedBp].idx; 1180 1181 // only look for non-basic variables 1182 if(breakpointidx >= 0) 1183 { 1184 foundStable = getData(val, leaveIdx, breakpointidx, moststable, degeneps, upd, vec, lb, ub, FVEC, 1185 max); 1186 } 1187 1188 --npassedBp; 1189 } 1190 1191 ++npassedBp; 1192 } 1193 1194 if(!foundStable) 1195 { 1196 assert(leaveIdx < 0); 1197 1198 if(relax_count < SOPLEX_MAX_RELAX_COUNT) 1199 { 1200 SPxOut::debug(this, "DEBFRT04 {}: no valid leaveIdx found - relaxing...\n", 1201 this->thesolver->basis().iteration()); 1202 this->relax(); 1203 ++relax_count; 1204 // restore original value 1205 val = max; 1206 // try again with relaxed delta 1207 return SPxBoundFlippingRT<R>::selectLeave(val, enterTest); 1208 } 1209 else 1210 { 1211 SPxOut::debug(this, "DEBFRT05 {}: no valid leaveIdx found - breaking...\n", 1212 this->thesolver->basis().iteration()); 1213 return leaveIdx; 1214 } 1215 } 1216 else 1217 { 1218 relax_count = 0; 1219 this->tighten(); 1220 } 1221 1222 // flip bounds of skipped breakpoints only if a nondegenerate step is to be performed 1223 if(npassedBp > 0 && spxAbs(breakpoints[npassedBp].val) > this->fastDelta) 1224 { 1225 flipAndUpdate(npassedBp); 1226 this->thesolver->boundflips = npassedBp; 1227 1228 if(npassedBp >= 10) 1229 flipPotential = 1; 1230 else 1231 flipPotential -= 0.05; 1232 } 1233 else 1234 { 1235 this->thesolver->boundflips = 0; 1236 flipPotential -= 0.1; 1237 } 1238 1239 SPxOut::debug(this, "DEBFRT06 {}: selected Index: {} number of candidates: {}\n", 1240 this->thesolver->basis().iteration(), leaveIdx, nBp); 1241 1242 return leaveIdx; 1243 } 1244 1245 1246 } // namespace soplex 1247