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 //TODO may be faster to have a greater zero tolerance for sparse pricing vectors 26 // to reduce the number of nonzero entries, e.g. for workVec 27 28 #include <assert.h> 29 #include <iostream> 30 31 #include "soplex/spxdefines.h" 32 #include "soplex/spxsteeppr.h" 33 #include "soplex/random.h" 34 35 #define SOPLEX_STEEP_REFINETOL 2.0 36 37 namespace soplex 38 { 39 40 template <class R> 41 void SPxSteepPR<R>::clear() 42 { 43 this->thesolver = 0; 44 } 45 46 template <class R> 47 void SPxSteepPR<R>::load(SPxSolverBase<R>* base) 48 { 49 this->thesolver = base; 50 51 if(base) 52 { 53 workVec.clear(); 54 workVec.reDim(base->dim()); 55 workRhs.clear(); 56 workRhs.reDim(base->dim()); 57 } 58 } 59 60 template <class R> 61 void SPxSteepPR<R>::setType(typename SPxSolverBase<R>::Type type) 62 { 63 workRhs.setTolerances(this->_tolerances); 64 65 setupWeights(type); 66 workVec.clear(); 67 workRhs.clear(); 68 refined = false; 69 70 bestPrices.clear(); 71 bestPrices.setMax(this->thesolver->dim()); 72 prices.reSize(this->thesolver->dim()); 73 74 if(type == SPxSolverBase<R>::ENTER) 75 { 76 bestPricesCo.clear(); 77 bestPricesCo.setMax(this->thesolver->coDim()); 78 pricesCo.reSize(this->thesolver->coDim()); 79 } 80 } 81 82 template <class R> 83 void SPxSteepPR<R>::setupWeights(typename SPxSolverBase<R>::Type type) 84 { 85 int i; 86 int endDim = 0; 87 int endCoDim = 0; 88 VectorBase<R>& weights = this->thesolver->weights; 89 VectorBase<R>& coWeights = this->thesolver->coWeights; 90 91 if(setup == DEFAULT) 92 { 93 if(type == SPxSolverBase<R>::ENTER) 94 { 95 if(this->thesolver->weightsAreSetup) 96 { 97 // check for added/removed rows and adapt norms accordingly 98 if(coWeights.dim() < this->thesolver->dim()) 99 endDim = coWeights.dim(); 100 else 101 endDim = this->thesolver->dim(); 102 103 if(weights.dim() < this->thesolver->coDim()) 104 endCoDim = weights.dim(); 105 else 106 endCoDim = this->thesolver->coDim(); 107 } 108 109 coWeights.reDim(this->thesolver->dim(), false); 110 111 for(i = this->thesolver->dim() - 1; i >= endDim; --i) 112 coWeights[i] = 2.0; 113 114 weights.reDim(this->thesolver->coDim(), false); 115 116 for(i = this->thesolver->coDim() - 1; i >= endCoDim; --i) 117 weights[i] = 1.0; 118 } 119 else 120 { 121 assert(type == SPxSolverBase<R>::LEAVE); 122 123 if(this->thesolver->weightsAreSetup) 124 { 125 // check for added/removed rows and adapt norms accordingly 126 if(coWeights.dim() < this->thesolver->dim()) 127 endDim = coWeights.dim(); 128 else 129 endDim = this->thesolver->dim(); 130 } 131 132 coWeights.reDim(this->thesolver->dim(), false); 133 134 for(i = this->thesolver->dim() - 1; i >= endDim; --i) 135 coWeights[i] = 1.0; 136 } 137 } 138 else 139 { 140 SPX_MSG_INFO1((*this->thesolver->spxout), 141 (*this->thesolver->spxout) << " --- initializing steepest edge multipliers" << std::endl;) 142 143 if(type == SPxSolverBase<R>::ENTER) 144 { 145 coWeights.reDim(this->thesolver->dim(), false); 146 147 for(i = this->thesolver->dim() - 1; i >= endDim; --i) 148 coWeights[i] = 1.0; 149 150 weights.reDim(this->thesolver->coDim(), false); 151 152 for(i = this->thesolver->coDim() - 1; i >= endCoDim; --i) 153 weights[i] = 1.0 + this->thesolver->vector(i).length2(); 154 } 155 else 156 { 157 assert(type == SPxSolverBase<R>::LEAVE); 158 coWeights.reDim(this->thesolver->dim(), false); 159 SSVectorBase<R> tmp(this->thesolver->dim(), this->thesolver->tolerances()); 160 161 for(i = this->thesolver->dim() - 1; i >= endDim && !this->thesolver->isTimeLimitReached(); --i) 162 { 163 this->thesolver->basis().coSolve(tmp, this->thesolver->unitVector(i)); 164 coWeights[i] = tmp.length2(); 165 } 166 } 167 } 168 169 this->thesolver->weightsAreSetup = true; 170 } 171 172 template <class R> 173 void SPxSteepPR<R>::setRep(typename SPxSolverBase<R>::Representation) 174 { 175 if(workVec.dim() != this->thesolver->dim()) 176 { 177 VectorBase<R> tmp = this->thesolver->weights; 178 this->thesolver->weights = this->thesolver->coWeights; 179 this->thesolver->coWeights = tmp; 180 181 workVec.clear(); 182 workVec.reDim(this->thesolver->dim()); 183 } 184 } 185 186 template <class R> 187 void SPxSteepPR<R>::left4(int n, SPxId id) 188 { 189 assert(this->thesolver->type() == SPxSolverBase<R>::LEAVE); 190 191 if(id.isValid()) 192 { 193 R delta = 0.1 + 1.0 / this->thesolver->basis().iteration(); 194 R* coWeights_ptr = this->thesolver->coWeights.get_ptr(); 195 const R* workVec_ptr = workVec.get_const_ptr(); 196 const R* rhoVec = this->thesolver->fVec().delta().values(); 197 R rhov_1 = 1.0 / rhoVec[n]; 198 R beta_q = this->thesolver->coPvec().delta().length2() * rhov_1 * rhov_1; 199 200 #ifndef NDEBUG 201 202 if(spxAbs(rhoVec[n]) < this->thetolerance * 0.5) 203 { 204 SPX_MSG_INFO3((*this->thesolver->spxout), (*this->thesolver->spxout) << "WSTEEP04: rhoVec = " 205 << rhoVec[n] << " with smaller absolute value than 0.5*thetolerance = " << 0.5 * this->thetolerance 206 << std::endl;) 207 } 208 209 #endif 210 211 const IdxSet& rhoIdx = this->thesolver->fVec().idx(); 212 int len = this->thesolver->fVec().idx().size(); 213 214 for(int i = 0; i < len; ++i) 215 { 216 int j = rhoIdx.index(i); 217 coWeights_ptr[j] += rhoVec[j] * (beta_q * rhoVec[j] - 2.0 * rhov_1 * workVec_ptr[j]); 218 219 if(coWeights_ptr[j] < delta) 220 coWeights_ptr[j] = delta; // coWeights_ptr[j] = delta / (1+delta-x); 221 else if(coWeights_ptr[j] >= R(infinity)) 222 coWeights_ptr[j] = 1.0 / this->thetolerance; 223 } 224 225 coWeights_ptr[n] = beta_q; 226 } 227 } 228 229 230 namespace steeppr 231 { 232 template <class R> 233 R computePrice(R viol, R weight, R tol) 234 { 235 if(weight < tol) 236 return viol * viol / tol; 237 else 238 return viol * viol / weight; 239 } 240 } 241 242 template <class R> 243 int SPxSteepPR<R>::buildBestPriceVectorLeave(R feastol) 244 { 245 int idx; 246 int nsorted; 247 R x; 248 const R* fTest = this->thesolver->fTest().get_const_ptr(); 249 const R* cpen = this->thesolver->coWeights.get_const_ptr(); 250 typename SPxPricer<R>::IdxElement price; 251 prices.clear(); 252 bestPrices.clear(); 253 254 // construct vector of all prices 255 for(int i = this->thesolver->infeasibilities.size() - 1; i >= 0; --i) 256 { 257 idx = this->thesolver->infeasibilities.index(i); 258 x = fTest[idx]; 259 260 if(x < -feastol) 261 { 262 // it might happen that we call the pricer with a tighter tolerance than what was used when computing the violations 263 this->thesolver->isInfeasible[idx] = this->VIOLATED; 264 price.val = steeppr::computePrice(x, cpen[idx], feastol); 265 price.idx = idx; 266 prices.append(price); 267 } 268 } 269 270 // set up structures for the quicksort implementation 271 this->compare.elements = prices.get_const_ptr(); 272 // do a partial sort to move the best ones to the front 273 // TODO this can be done more efficiently, since we only need the indices 274 nsorted = SPxQuicksortPart(prices.get_ptr(), this->compare, 0, prices.size(), 275 SOPLEX_HYPERPRICINGSIZE); 276 277 // copy indices of best values to bestPrices 278 for(int i = 0; i < nsorted; ++i) 279 { 280 bestPrices.addIdx(prices[i].idx); 281 this->thesolver->isInfeasible[prices[i].idx] = this->VIOLATED_AND_CHECKED; 282 } 283 284 if(nsorted > 0) 285 return prices[0].idx; 286 else 287 return -1; 288 } 289 290 291 template <class R> 292 int SPxSteepPR<R>::selectLeave() 293 { 294 assert(isConsistent()); 295 296 int retid; 297 298 if(this->thesolver->hyperPricingLeave && this->thesolver->sparsePricingLeave) 299 { 300 if(bestPrices.size() < 2 || this->thesolver->basis().lastUpdate() == 0) 301 { 302 // call init method to build up price-vector and return index of largest price 303 retid = buildBestPriceVectorLeave(this->thetolerance); 304 } 305 else 306 retid = selectLeaveHyper(this->thetolerance); 307 } 308 else if(this->thesolver->sparsePricingLeave) 309 retid = selectLeaveSparse(this->thetolerance); 310 else 311 retid = selectLeaveX(this->thetolerance); 312 313 if(retid < 0 && !refined) 314 { 315 refined = true; 316 SPX_MSG_INFO3((*this->thesolver->spxout), 317 (*this->thesolver->spxout) << "WSTEEP03 trying refinement step..\n";) 318 retid = selectLeaveX(this->thetolerance / SOPLEX_STEEP_REFINETOL); 319 } 320 321 if(retid >= 0) 322 { 323 assert(this->thesolver->coPvec().delta().isConsistent()); 324 // coPvec().delta() might be not setup after the solve when it contains too many nonzeros. 325 // This is intended and forcing to keep the sparsity information leads to a slowdown 326 // TODO implement a dedicated solve method for unitvectors 327 this->thesolver->basis().coSolve(this->thesolver->coPvec().delta(), 328 this->thesolver->unitVector(retid)); 329 assert(this->thesolver->coPvec().delta().isConsistent()); 330 workRhs.setup_and_assign(this->thesolver->coPvec().delta()); 331 this->thesolver->setup4solve(&workVec, &workRhs); 332 } 333 334 return retid; 335 } 336 337 template <class R> 338 int SPxSteepPR<R>::selectLeaveX(R tol) 339 { 340 const R* coWeights_ptr = this->thesolver->coWeights.get_const_ptr(); 341 const R* fTest = this->thesolver->fTest().get_const_ptr(); 342 R best = R(-infinity); 343 R x; 344 int lastIdx = -1; 345 346 for(int i = this->thesolver->dim() - 1; i >= 0; --i) 347 { 348 x = fTest[i]; 349 350 if(x < -tol) 351 { 352 x = steeppr::computePrice(x, coWeights_ptr[i], tol); 353 354 if(x > best) 355 { 356 best = x; 357 lastIdx = i; 358 } 359 } 360 } 361 362 return lastIdx; 363 } 364 365 template <class R> 366 int SPxSteepPR<R>::selectLeaveSparse(R tol) 367 { 368 const R* coWeights_ptr = this->thesolver->coWeights.get_const_ptr(); 369 const R* fTest = this->thesolver->fTest().get_const_ptr(); 370 R best = R(-infinity); 371 R x; 372 int lastIdx = -1; 373 int idx; 374 375 for(int i = this->thesolver->infeasibilities.size() - 1; i >= 0; --i) 376 { 377 idx = this->thesolver->infeasibilities.index(i); 378 x = fTest[idx]; 379 380 if(x < -tol) 381 { 382 x = steeppr::computePrice(x, coWeights_ptr[idx], tol); 383 384 if(x > best) 385 { 386 best = x; 387 lastIdx = idx; 388 } 389 } 390 else 391 { 392 this->thesolver->infeasibilities.remove(i); 393 assert(this->thesolver->isInfeasible[idx] == this->VIOLATED 394 || this->thesolver->isInfeasible[idx] == this->VIOLATED_AND_CHECKED); 395 this->thesolver->isInfeasible[idx] = this->NOT_VIOLATED; 396 } 397 } 398 399 return lastIdx; 400 } 401 402 template <class R> 403 int SPxSteepPR<R>::selectLeaveHyper(R tol) 404 { 405 const R* coPen = this->thesolver->coWeights.get_const_ptr(); 406 const R* fTest = this->thesolver->fTest().get_const_ptr(); 407 R leastBest = -1; 408 R best = R(-infinity); 409 R x; 410 int bestIdx = -1; 411 int idx = 0; 412 413 // find the best price from the short candidate list 414 for(int i = bestPrices.size() - 1; i >= 0; --i) 415 { 416 idx = bestPrices.index(i); 417 x = fTest[idx]; 418 419 if(x < -tol) 420 { 421 assert(this->thesolver->isInfeasible[idx] == this->VIOLATED 422 || this->thesolver->isInfeasible[idx] == this->VIOLATED_AND_CHECKED); 423 x = steeppr::computePrice(x, coPen[idx], tol); 424 425 assert(x >= 0); 426 427 // update the best price of candidate list 428 if(x > best) 429 { 430 best = x; 431 bestIdx = idx; 432 } 433 434 // update the smallest price of candidate list 435 if(x < leastBest || leastBest < 0) 436 leastBest = x; 437 } 438 else 439 { 440 bestPrices.remove(i); 441 this->thesolver->isInfeasible[idx] = this->NOT_VIOLATED; 442 } 443 } 444 445 // scan the updated indices for a better price 446 for(int i = this->thesolver->updateViols.size() - 1; i >= 0; --i) 447 { 448 idx = this->thesolver->updateViols.index(i); 449 450 // is this index a candidate for bestPrices? 451 if(this->thesolver->isInfeasible[idx] == this->VIOLATED) 452 { 453 x = fTest[idx]; 454 assert(x < -tol); 455 x = steeppr::computePrice(x, coPen[idx], tol); 456 457 if(x > leastBest) 458 { 459 if(x > best) 460 { 461 best = x; 462 bestIdx = idx; 463 } 464 465 this->thesolver->isInfeasible[idx] = this->VIOLATED_AND_CHECKED; 466 bestPrices.addIdx(idx); 467 } 468 } 469 } 470 471 return bestIdx; 472 } 473 474 /* Entering Simplex 475 */ 476 template <class R> 477 void SPxSteepPR<R>::entered4(SPxId /* id */, int n) 478 { 479 assert(this->thesolver->type() == SPxSolverBase<R>::ENTER); 480 481 if(n >= 0 && n < this->thesolver->dim()) 482 { 483 R delta = 2 + 1.0 / this->thesolver->basis().iteration(); 484 R* coWeights_ptr = this->thesolver->coWeights.get_ptr(); 485 R* weights_ptr = this->thesolver->weights.get_ptr(); 486 const R* workVec_ptr = workVec.get_const_ptr(); 487 const R* pVec = this->thesolver->pVec().delta().values(); 488 const IdxSet& pIdx = this->thesolver->pVec().idx(); 489 const R* coPvec = this->thesolver->coPvec().delta().values(); 490 const IdxSet& coPidx = this->thesolver->coPvec().idx(); 491 R xi_p = 1 / this->thesolver->fVec().delta()[n]; 492 int i, j; 493 R xi_ip; 494 495 assert(this->thesolver->fVec().delta()[n] > this->thesolver->epsilon() 496 || this->thesolver->fVec().delta()[n] < -this->thesolver->epsilon()); 497 498 for(j = coPidx.size() - 1; j >= 0; --j) 499 { 500 i = coPidx.index(j); 501 xi_ip = xi_p * coPvec[i]; 502 coWeights_ptr[i] += xi_ip * (xi_ip * pi_p - 2.0 * workVec_ptr[i]); 503 504 /* 505 if(coWeights_ptr[i] < 1) 506 coWeights_ptr[i] = 1 / (2-x); 507 */ 508 if(coWeights_ptr[i] < delta) 509 coWeights_ptr[i] = delta; 510 // coWeights_ptr[i] = 1; 511 else if(coWeights_ptr[i] > R(infinity)) 512 coWeights_ptr[i] = 1 / this->thesolver->epsilon(); 513 } 514 515 for(j = pIdx.size() - 1; j >= 0; --j) 516 { 517 i = pIdx.index(j); 518 xi_ip = xi_p * pVec[i]; 519 weights_ptr[i] += xi_ip * (xi_ip * pi_p - 2.0 * (this->thesolver->vector(i) * workVec)); 520 521 /* 522 if(weights_ptr[i] < 1) 523 weights_ptr[i] = 1 / (2-x); 524 */ 525 if(weights_ptr[i] < delta) 526 weights_ptr[i] = delta; 527 // weights_ptr[i] = 1; 528 else if(weights_ptr[i] > R(infinity)) 529 weights_ptr[i] = 1.0 / this->thesolver->epsilon(); 530 } 531 } 532 533 /*@ 534 if(this->thesolver->isId(id)) 535 weights[ this->thesolver->number(id) ] *= 1.0001; 536 else if(this->thesolver->isCoId(id)) 537 coWeights[ this->thesolver->number(id) ] *= 1.0001; 538 */ 539 540 } 541 542 543 template <class R> 544 SPxId SPxSteepPR<R>::buildBestPriceVectorEnterDim(R& best, R feastol) 545 { 546 const R* coTest = this->thesolver->coTest().get_const_ptr(); 547 const R* coWeights_ptr = this->thesolver->coWeights.get_const_ptr(); 548 int idx; 549 int nsorted; 550 R x; 551 typename SPxPricer<R>::IdxElement price; 552 553 prices.clear(); 554 bestPrices.clear(); 555 556 // construct vector of all prices 557 for(int i = this->thesolver->infeasibilities.size() - 1; i >= 0; --i) 558 { 559 idx = this->thesolver->infeasibilities.index(i); 560 x = coTest[idx]; 561 562 if(x < -feastol) 563 { 564 // it might happen that we call the pricer with a tighter tolerance than what was used when computing the violations 565 this->thesolver->isInfeasible[idx] = this->VIOLATED; 566 price.val = steeppr::computePrice(x, coWeights_ptr[idx], feastol); 567 price.idx = idx; 568 prices.append(price); 569 } 570 else 571 { 572 this->thesolver->infeasibilities.remove(i); 573 this->thesolver->isInfeasible[idx] = this->NOT_VIOLATED; 574 } 575 } 576 577 // set up structures for the quicksort implementation 578 this->compare.elements = prices.get_const_ptr(); 579 // do a partial sort to move the best ones to the front 580 // TODO this can be done more efficiently, since we only need the indices 581 nsorted = SPxQuicksortPart(prices.get_ptr(), this->compare, 0, prices.size(), 582 SOPLEX_HYPERPRICINGSIZE); 583 584 // copy indices of best values to bestPrices 585 for(int i = 0; i < nsorted; ++i) 586 { 587 bestPrices.addIdx(prices[i].idx); 588 this->thesolver->isInfeasible[prices[i].idx] = this->VIOLATED_AND_CHECKED; 589 } 590 591 if(nsorted > 0) 592 { 593 best = prices[0].val; 594 return this->thesolver->coId(prices[0].idx); 595 } 596 else 597 return SPxId(); 598 } 599 600 601 template <class R> 602 SPxId SPxSteepPR<R>::buildBestPriceVectorEnterCoDim(R& best, R feastol) 603 { 604 const R* test = this->thesolver->test().get_const_ptr(); 605 const R* weights_ptr = this->thesolver->weights.get_const_ptr(); 606 int idx; 607 int nsorted; 608 R x; 609 typename SPxPricer<R>::IdxElement price; 610 611 pricesCo.clear(); 612 bestPricesCo.clear(); 613 614 // construct vector of all prices 615 for(int i = this->thesolver->infeasibilitiesCo.size() - 1; i >= 0; --i) 616 { 617 idx = this->thesolver->infeasibilitiesCo.index(i); 618 x = test[idx]; 619 620 if(x < -feastol) 621 { 622 // it might happen that we call the pricer with a tighter tolerance than what was used when computing the violations 623 this->thesolver->isInfeasibleCo[idx] = this->VIOLATED; 624 price.val = steeppr::computePrice(x, weights_ptr[idx], feastol); 625 price.idx = idx; 626 pricesCo.append(price); 627 } 628 else 629 { 630 this->thesolver->infeasibilitiesCo.remove(i); 631 this->thesolver->isInfeasibleCo[idx] = this->NOT_VIOLATED; 632 } 633 } 634 635 // set up structures for the quicksort implementation 636 this->compare.elements = pricesCo.get_const_ptr(); 637 // do a partial sort to move the best ones to the front 638 // TODO this can be done more efficiently, since we only need the indices 639 nsorted = SPxQuicksortPart(pricesCo.get_ptr(), this->compare, 0, pricesCo.size(), 640 SOPLEX_HYPERPRICINGSIZE); 641 642 // copy indices of best values to bestPrices 643 for(int i = 0; i < nsorted; ++i) 644 { 645 bestPricesCo.addIdx(pricesCo[i].idx); 646 this->thesolver->isInfeasibleCo[pricesCo[i].idx] = this->VIOLATED_AND_CHECKED; 647 } 648 649 if(nsorted > 0) 650 { 651 best = pricesCo[0].val; 652 return this->thesolver->id(pricesCo[0].idx); 653 } 654 else 655 return SPxId(); 656 } 657 658 659 template <class R> 660 SPxId SPxSteepPR<R>::selectEnter() 661 { 662 assert(this->thesolver != 0); 663 SPxId enterId; 664 665 enterId = selectEnterX(this->thetolerance); 666 667 if(!enterId.isValid() && !refined) 668 { 669 refined = true; 670 SPX_MSG_INFO3((*this->thesolver->spxout), 671 (*this->thesolver->spxout) << "WSTEEP05 trying refinement step..\n";) 672 enterId = selectEnterX(this->thetolerance / SOPLEX_STEEP_REFINETOL); 673 } 674 675 assert(isConsistent()); 676 677 if(enterId.isValid()) 678 { 679 SSVectorBase<R>& delta = this->thesolver->fVec().delta(); 680 681 this->thesolver->basis().solve4update(delta, this->thesolver->vector(enterId)); 682 683 workRhs.setup_and_assign(delta); 684 pi_p = 1 + delta.length2(); 685 686 this->thesolver->setup4coSolve(&workVec, &workRhs); 687 } 688 689 return enterId; 690 } 691 692 template <class R> 693 SPxId SPxSteepPR<R>::selectEnterX(R tol) 694 { 695 SPxId enterId; 696 SPxId enterCoId; 697 R best; 698 R bestCo; 699 700 best = R(-infinity); 701 bestCo = R(-infinity); 702 703 if(this->thesolver->hyperPricingEnter && !refined) 704 { 705 if(bestPrices.size() < 2 || this->thesolver->basis().lastUpdate() == 0) 706 enterCoId = (this->thesolver->sparsePricingEnter) ? buildBestPriceVectorEnterDim(best, 707 tol) : selectEnterDenseDim(best, tol); 708 else 709 enterCoId = (this->thesolver->sparsePricingEnter) ? selectEnterHyperDim(best, 710 tol) : selectEnterDenseDim(best, tol); 711 712 if(bestPricesCo.size() < 2 || this->thesolver->basis().lastUpdate() == 0) 713 enterId = (this->thesolver->sparsePricingEnterCo) ? buildBestPriceVectorEnterCoDim(bestCo, 714 tol) : selectEnterDenseCoDim(bestCo, tol); 715 else 716 enterId = (this->thesolver->sparsePricingEnterCo) ? selectEnterHyperCoDim(bestCo, 717 tol) : selectEnterDenseCoDim(bestCo, tol); 718 } 719 else 720 { 721 enterCoId = (this->thesolver->sparsePricingEnter 722 && !refined) ? selectEnterSparseDim(best, tol) : selectEnterDenseDim(best, tol); 723 enterId = (this->thesolver->sparsePricingEnterCo 724 && !refined) ? selectEnterSparseCoDim(bestCo, tol) : selectEnterDenseCoDim(bestCo, tol); 725 } 726 727 // prefer slack indices to reduce nonzeros in basis matrix 728 if(enterCoId.isValid() && (best > SOPLEX_SPARSITY_TRADEOFF * bestCo || !enterId.isValid())) 729 return enterCoId; 730 else 731 return enterId; 732 } 733 734 735 template <class R> 736 SPxId SPxSteepPR<R>::selectEnterHyperDim(R& best, R tol) 737 { 738 const R* coTest = this->thesolver->coTest().get_const_ptr(); 739 const R* coWeights_ptr = this->thesolver->coWeights.get_const_ptr(); 740 741 R leastBest = -1; 742 R x; 743 int enterIdx = -1; 744 int idx; 745 746 // find the best price from short candidate list 747 for(int i = bestPrices.size() - 1; i >= 0; --i) 748 { 749 idx = bestPrices.index(i); 750 x = coTest[idx]; 751 752 if(x < -tol) 753 { 754 x = steeppr::computePrice(x, coWeights_ptr[idx], tol); 755 756 assert(x >= 0); 757 758 // update the best price of candidate list 759 if(x > best) 760 { 761 best = x; 762 enterIdx = idx; 763 } 764 765 // update the smallest price of candidate list 766 if(x < leastBest || leastBest < 0) 767 leastBest = x; 768 } 769 else 770 { 771 bestPrices.remove(i); 772 this->thesolver->isInfeasible[idx] = this->NOT_VIOLATED; 773 } 774 } 775 776 // scan the updated indices for a better price 777 for(int i = this->thesolver->updateViols.size() - 1; i >= 0; --i) 778 { 779 idx = this->thesolver->updateViols.index(i); 780 781 // only look at indices that were not checked already 782 if(this->thesolver->isInfeasible[idx] == this->VIOLATED) 783 { 784 x = coTest[idx]; 785 786 if(x < -tol) 787 { 788 x = steeppr::computePrice(x, coWeights_ptr[idx], tol); 789 790 if(x > leastBest) 791 { 792 if(x > best) 793 { 794 best = x; 795 enterIdx = idx; 796 } 797 798 // put index into candidate list 799 this->thesolver->isInfeasible[idx] = this->VIOLATED_AND_CHECKED; 800 bestPrices.addIdx(idx); 801 } 802 } 803 else 804 { 805 this->thesolver->isInfeasible[idx] = this->NOT_VIOLATED; 806 } 807 } 808 } 809 810 if(enterIdx >= 0) 811 return this->thesolver->coId(enterIdx); 812 else 813 return SPxId(); 814 } 815 816 817 template <class R> 818 SPxId SPxSteepPR<R>::selectEnterHyperCoDim(R& best, R tol) 819 { 820 const R* test = this->thesolver->test().get_const_ptr(); 821 const R* weights_ptr = this->thesolver->weights.get_const_ptr(); 822 823 R leastBest = -1; 824 R x; 825 int enterIdx = -1; 826 int idx; 827 828 // find the best price from short candidate list 829 for(int i = bestPricesCo.size() - 1; i >= 0; --i) 830 { 831 idx = bestPricesCo.index(i); 832 x = test[idx]; 833 834 if(x < -tol) 835 { 836 x = steeppr::computePrice(x, weights_ptr[idx], tol); 837 838 assert(x >= 0); 839 840 // update the best price of candidate list 841 if(x > best) 842 { 843 best = x; 844 enterIdx = idx; 845 } 846 847 // update the smallest price of candidate list 848 if(x < leastBest || leastBest < 0) 849 leastBest = x; 850 } 851 else 852 { 853 bestPricesCo.remove(i); 854 this->thesolver->isInfeasibleCo[idx] = this->NOT_VIOLATED; 855 } 856 } 857 858 // scan the updated indices for a better price 859 for(int i = this->thesolver->updateViolsCo.size() - 1; i >= 0; --i) 860 { 861 idx = this->thesolver->updateViolsCo.index(i); 862 863 // only look at indices that were not checked already 864 if(this->thesolver->isInfeasibleCo[idx] == this->VIOLATED) 865 { 866 x = test[idx]; 867 868 if(x < -tol) 869 { 870 x = steeppr::computePrice(x, weights_ptr[idx], tol); 871 872 if(x > leastBest) 873 { 874 if(x > best) 875 { 876 best = x; 877 enterIdx = idx; 878 } 879 880 // put index into candidate list 881 this->thesolver->isInfeasibleCo[idx] = this->VIOLATED_AND_CHECKED; 882 bestPricesCo.addIdx(idx); 883 } 884 } 885 else 886 { 887 this->thesolver->isInfeasibleCo[idx] = this->NOT_VIOLATED; 888 } 889 } 890 } 891 892 if(enterIdx >= 0) 893 return this->thesolver->id(enterIdx); 894 else 895 return SPxId(); 896 } 897 898 899 template <class R> 900 SPxId SPxSteepPR<R>::selectEnterSparseDim(R& best, R tol) 901 { 902 SPxId enterId; 903 const R* coTest = this->thesolver->coTest().get_const_ptr(); 904 const R* coWeights_ptr = this->thesolver->coWeights.get_const_ptr(); 905 906 int idx; 907 R x; 908 909 for(int i = this->thesolver->infeasibilities.size() - 1; i >= 0; --i) 910 { 911 idx = this->thesolver->infeasibilities.index(i); 912 x = coTest[idx]; 913 914 if(x < -tol) 915 { 916 x = steeppr::computePrice(x, coWeights_ptr[idx], tol); 917 918 if(x > best) 919 { 920 best = x; 921 enterId = this->thesolver->coId(idx); 922 } 923 } 924 else 925 { 926 this->thesolver->infeasibilities.remove(i); 927 this->thesolver->isInfeasible[idx] = this->NOT_VIOLATED; 928 } 929 } 930 931 return enterId; 932 } 933 934 template <class R> 935 SPxId SPxSteepPR<R>::selectEnterSparseCoDim(R& best, R tol) 936 { 937 SPxId enterId; 938 const R* test = this->thesolver->test().get_const_ptr(); 939 const R* weights_ptr = this->thesolver->weights.get_const_ptr(); 940 941 int idx; 942 R x; 943 944 for(int i = this->thesolver->infeasibilitiesCo.size() - 1; i >= 0; --i) 945 { 946 idx = this->thesolver->infeasibilitiesCo.index(i); 947 x = test[idx]; 948 949 if(x < -tol) 950 { 951 x = steeppr::computePrice(x, weights_ptr[idx], tol); 952 953 if(x > best) 954 { 955 best = x; 956 enterId = this->thesolver->id(idx); 957 } 958 } 959 else 960 { 961 this->thesolver->infeasibilitiesCo.remove(i); 962 this->thesolver->isInfeasibleCo[idx] = this->NOT_VIOLATED; 963 } 964 } 965 966 return enterId; 967 } 968 969 template <class R> 970 SPxId SPxSteepPR<R>::selectEnterDenseDim(R& best, R tol) 971 { 972 SPxId enterId; 973 const R* coTest = this->thesolver->coTest().get_const_ptr(); 974 const R* coWeights_ptr = this->thesolver->coWeights.get_const_ptr(); 975 976 R x; 977 978 for(int i = 0, end = this->thesolver->dim(); i < end; ++i) 979 { 980 x = coTest[i]; 981 982 if(x < -tol) 983 { 984 x = steeppr::computePrice(x, coWeights_ptr[i], tol); 985 986 if(x > best) 987 { 988 best = x; 989 enterId = this->thesolver->coId(i); 990 } 991 } 992 } 993 994 return enterId; 995 } 996 997 template <class R> 998 SPxId SPxSteepPR<R>::selectEnterDenseCoDim(R& best, R tol) 999 { 1000 SPxId enterId; 1001 const R* test = this->thesolver->test().get_const_ptr(); 1002 const R* weights_ptr = this->thesolver->weights.get_const_ptr(); 1003 1004 R x; 1005 1006 for(int i = 0, end = this->thesolver->coDim(); i < end; ++i) 1007 { 1008 x = test[i]; 1009 1010 if(x < -tol) 1011 { 1012 x = steeppr::computePrice(x, weights_ptr[i], tol); 1013 1014 if(x > best) 1015 { 1016 best = x; 1017 enterId = this->thesolver->id(i); 1018 } 1019 } 1020 } 1021 1022 return enterId; 1023 } 1024 1025 1026 template <class R> 1027 void SPxSteepPR<R>::addedVecs(int n) 1028 { 1029 VectorBase<R>& weights = this->thesolver->weights; 1030 n = weights.dim(); 1031 weights.reDim(this->thesolver->coDim()); 1032 1033 if(this->thesolver->type() == SPxSolverBase<R>::ENTER) 1034 { 1035 for(; n < weights.dim(); ++n) 1036 weights[n] = 2; 1037 } 1038 } 1039 1040 template <class R> 1041 void SPxSteepPR<R>::addedCoVecs(int n) 1042 { 1043 VectorBase<R>& coWeights = this->thesolver->coWeights; 1044 n = coWeights.dim(); 1045 workVec.reDim(this->thesolver->dim()); 1046 coWeights.reDim(this->thesolver->dim()); 1047 1048 for(; n < coWeights.dim(); ++n) 1049 coWeights[n] = 1; 1050 } 1051 1052 template <class R> 1053 void SPxSteepPR<R>::removedVec(int i) 1054 { 1055 assert(this->thesolver != 0); 1056 VectorBase<R>& weights = this->thesolver->weights; 1057 weights[i] = weights[weights.dim()]; 1058 weights.reDim(this->thesolver->coDim()); 1059 } 1060 1061 template <class R> 1062 void SPxSteepPR<R>::removedVecs(const int perm[]) 1063 { 1064 assert(this->thesolver != 0); 1065 VectorBase<R>& weights = this->thesolver->weights; 1066 1067 if(this->thesolver->type() == SPxSolverBase<R>::ENTER) 1068 { 1069 int i; 1070 int j = weights.dim(); 1071 1072 for(i = 0; i < j; ++i) 1073 { 1074 if(perm[i] >= 0) 1075 weights[perm[i]] = weights[i]; 1076 } 1077 } 1078 1079 weights.reDim(this->thesolver->coDim()); 1080 } 1081 1082 template <class R> 1083 void SPxSteepPR<R>::removedCoVec(int i) 1084 { 1085 assert(this->thesolver != 0); 1086 VectorBase<R>& coWeights = this->thesolver->coWeights; 1087 coWeights[i] = coWeights[coWeights.dim()]; 1088 coWeights.reDim(this->thesolver->dim()); 1089 } 1090 1091 template <class R> 1092 void SPxSteepPR<R>::removedCoVecs(const int perm[]) 1093 { 1094 assert(this->thesolver != 0); 1095 VectorBase<R>& coWeights = this->thesolver->coWeights; 1096 int i; 1097 int j = coWeights.dim(); 1098 1099 for(i = 0; i < j; ++i) 1100 { 1101 if(perm[i] >= 0) 1102 coWeights[perm[i]] = coWeights[i]; 1103 } 1104 1105 coWeights.reDim(this->thesolver->dim()); 1106 } 1107 1108 template <class R> 1109 bool SPxSteepPR<R>::isConsistent() const 1110 { 1111 #ifdef ENABLE_CONSISTENCY_CHECKS 1112 VectorBase<R>& w = this->thesolver->weights; 1113 VectorBase<R>& coW = this->thesolver->coWeights; 1114 1115 if(this->thesolver != 0 && this->thesolver->type() == SPxSolverBase<R>::LEAVE && setup == EXACT) 1116 { 1117 int i; 1118 SSVectorBase<R> tmp(this->thesolver->dim(), this->thesolver->tolerances()); 1119 R x; 1120 1121 for(i = this->thesolver->dim() - 1; i >= 0; --i) 1122 { 1123 this->thesolver->basis().coSolve(tmp, this->thesolver->unitVector(i)); 1124 x = coW[i] - tmp.length2(); 1125 1126 if(x > this->thesolver->leavetol() || -x > this->thesolver->leavetol()) 1127 { 1128 SPX_MSG_ERROR(std::cerr << "ESTEEP03 x[" << i << "] = " << x << std::endl;) 1129 } 1130 } 1131 } 1132 1133 if(this->thesolver != 0 && this->thesolver->type() == SPxSolverBase<R>::ENTER) 1134 { 1135 int i; 1136 1137 for(i = this->thesolver->dim() - 1; i >= 0; --i) 1138 { 1139 if(coW[i] < this->thesolver->epsilon()) 1140 return SPX_MSG_INCONSISTENT("SPxSteepPR"); 1141 } 1142 1143 for(i = this->thesolver->coDim() - 1; i >= 0; --i) 1144 { 1145 if(w[i] < this->thesolver->epsilon()) 1146 return SPX_MSG_INCONSISTENT("SPxSteepPR"); 1147 } 1148 } 1149 1150 #endif 1151 1152 return true; 1153 } 1154 } // namespace soplex 1155