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 "soplex/spxdefines.h" 26 27 #define SOPLEX_DEVEX_REFINETOL 2.0 28 29 namespace soplex 30 { 31 // Definition of signature to avoid the specialization after instantiation error 32 33 template <class R> 34 void SPxDevexPR<R>::load(SPxSolverBase<R>* base) 35 { 36 this->thesolver = base; 37 setRep(base->rep()); 38 assert(isConsistent()); 39 } 40 41 template <class R> 42 bool SPxDevexPR<R>::isConsistent() const 43 { 44 #ifdef ENABLE_CONSISTENCY_CHECKS 45 46 if(this->thesolver != 0) 47 if(this->thesolver->weights.dim() != this->thesolver->coDim() 48 || this->thesolver->coWeights.dim() != this->thesolver->dim()) 49 return SPX_MSG_INCONSISTENT("SPxDevexPR"); 50 51 #endif 52 53 return true; 54 } 55 56 template <class R> 57 void SPxDevexPR<R>::setupWeights(typename SPxSolverBase<R>::Type tp) 58 { 59 int i; 60 int coWeightSize = 0; 61 int weightSize = 0; 62 63 VectorBase<R>& weights = this->thesolver->weights; 64 VectorBase<R>& coWeights = this->thesolver->coWeights; 65 66 if(tp == SPxSolverBase<R>::ENTER) 67 { 68 coWeights.reDim(this->thesolver->dim(), false); 69 70 for(i = this->thesolver->dim() - 1; i >= coWeightSize; --i) 71 coWeights[i] = 2.0; 72 73 weights.reDim(this->thesolver->coDim(), false); 74 75 for(i = this->thesolver->coDim() - 1; i >= weightSize; --i) 76 weights[i] = 2.0; 77 } 78 else 79 { 80 coWeights.reDim(this->thesolver->dim(), false); 81 82 for(i = this->thesolver->dim() - 1; i >= coWeightSize; --i) 83 coWeights[i] = 1.0; 84 } 85 86 this->thesolver->weightsAreSetup = true; 87 } 88 89 template <class R> 90 void SPxDevexPR<R>::setType(typename SPxSolverBase<R>::Type tp) 91 { 92 setupWeights(tp); 93 refined = false; 94 95 bestPrices.clear(); 96 bestPrices.setMax(this->thesolver->dim()); 97 prices.reSize(this->thesolver->dim()); 98 99 if(tp == SPxSolverBase<R>::ENTER) 100 { 101 bestPricesCo.clear(); 102 bestPricesCo.setMax(this->thesolver->coDim()); 103 pricesCo.reSize(this->thesolver->coDim()); 104 } 105 106 assert(isConsistent()); 107 } 108 109 /**@todo suspicious: Shouldn't the relation between dim, coDim, Vecs, 110 * and CoVecs be influenced by the representation ? 111 */ 112 template <class R> 113 void SPxDevexPR<R>::setRep(typename SPxSolverBase<R>::Representation) 114 { 115 if(this->thesolver != 0) 116 { 117 // resize weights and initialize new entries 118 addedVecs(this->thesolver->coDim()); 119 addedCoVecs(this->thesolver->dim()); 120 assert(isConsistent()); 121 } 122 } 123 124 // A namespace to avoid collision with other pricers 125 namespace devexpr 126 { 127 template <class R> 128 R computePrice(R viol, R weight, R tol) 129 { 130 if(weight < tol) 131 return viol * viol / tol; 132 else 133 return viol * viol / weight; 134 } 135 } 136 137 template <class R> 138 int SPxDevexPR<R>::buildBestPriceVectorLeave(R feastol) 139 { 140 int idx; 141 int nsorted; 142 R fTesti; 143 const R* fTest = this->thesolver->fTest().get_const_ptr(); 144 const R* cpen = this->thesolver->coWeights.get_const_ptr(); 145 typename SPxPricer<R>::IdxElement price; 146 prices.clear(); 147 bestPrices.clear(); 148 149 // TODO we should check infeasiblities for duplicates or loop over dimension 150 // bestPrices may then also contain duplicates! 151 // construct vector of all prices 152 for(int i = this->thesolver->infeasibilities.size() - 1; i >= 0; --i) 153 { 154 idx = this->thesolver->infeasibilities.index(i); 155 fTesti = fTest[idx]; 156 157 if(fTesti < -feastol) 158 { 159 this->thesolver->isInfeasible[idx] = this->VIOLATED; 160 price.idx = idx; 161 price.val = devexpr::computePrice(fTesti, cpen[idx], feastol); 162 prices.append(price); 163 } 164 } 165 166 // set up structures for the quicksort implementation 167 this->compare.elements = prices.get_const_ptr(); 168 // do a partial sort to move the best ones to the front 169 // TODO this can be done more efficiently, since we only need the indices 170 nsorted = SPxQuicksortPart(prices.get_ptr(), this->compare, 0, prices.size(), 171 SOPLEX_HYPERPRICINGSIZE); 172 173 // copy indices of best values to bestPrices 174 for(int i = 0; i < nsorted; ++i) 175 { 176 bestPrices.addIdx(prices[i].idx); 177 this->thesolver->isInfeasible[prices[i].idx] = this->VIOLATED_AND_CHECKED; 178 } 179 180 if(nsorted > 0) 181 return prices[0].idx; 182 else 183 return -1; 184 } 185 186 template <class R> 187 int SPxDevexPR<R>::selectLeave() 188 { 189 int retid; 190 191 if(this->thesolver->hyperPricingLeave && this->thesolver->sparsePricingLeave) 192 { 193 if(bestPrices.size() < 2 || this->thesolver->basis().lastUpdate() == 0) 194 { 195 // call init method to build up price-vector and return index of largest price 196 retid = buildBestPriceVectorLeave(this->thetolerance); 197 } 198 else 199 retid = selectLeaveHyper(this->thetolerance); 200 } 201 else if(this->thesolver->sparsePricingLeave) 202 retid = selectLeaveSparse(this->thetolerance); 203 else 204 retid = selectLeaveX(this->thetolerance); 205 206 if(retid < 0 && !refined) 207 { 208 refined = true; 209 SPX_MSG_INFO3((*this->thesolver->spxout), 210 (*this->thesolver->spxout) << "WDEVEX02 trying refinement step..\n";) 211 retid = selectLeaveX(this->thetolerance / SOPLEX_DEVEX_REFINETOL); 212 } 213 214 assert(retid < this->thesolver->dim()); 215 216 return retid; 217 } 218 219 template <class R> 220 int SPxDevexPR<R>::selectLeaveX(R feastol, int start, int incr) 221 { 222 R x; 223 224 const R* fTest = this->thesolver->fTest().get_const_ptr(); 225 const R* cpen = this->thesolver->coWeights.get_const_ptr(); 226 R best = 0; 227 int bstI = -1; 228 int end = this->thesolver->coWeights.dim(); 229 230 for(; start < end; start += incr) 231 { 232 if(fTest[start] < -feastol) 233 { 234 x = devexpr::computePrice(fTest[start], cpen[start], feastol); 235 236 if(x > best) 237 { 238 best = x; 239 bstI = start; 240 last = cpen[start]; 241 } 242 } 243 } 244 245 return bstI; 246 } 247 248 template <class R> 249 int SPxDevexPR<R>::selectLeaveSparse(R feastol) 250 { 251 R x; 252 253 const R* fTest = this->thesolver->fTest().get_const_ptr(); 254 const R* cpen = this->thesolver->coWeights.get_const_ptr(); 255 R best = 0; 256 int bstI = -1; 257 int idx = -1; 258 259 for(int i = this->thesolver->infeasibilities.size() - 1; i >= 0; --i) 260 { 261 idx = this->thesolver->infeasibilities.index(i); 262 x = fTest[idx]; 263 264 if(x < -feastol) 265 { 266 x = devexpr::computePrice(x, cpen[idx], feastol); 267 268 if(x > best) 269 { 270 best = x; 271 bstI = idx; 272 last = cpen[idx]; 273 } 274 } 275 else 276 { 277 this->thesolver->infeasibilities.remove(i); 278 assert(this->thesolver->isInfeasible[idx] == this->VIOLATED 279 || this->thesolver->isInfeasible[idx] == this->VIOLATED_AND_CHECKED); 280 this->thesolver->isInfeasible[idx] = this->NOT_VIOLATED; 281 } 282 } 283 284 return bstI; 285 } 286 287 template <class R> 288 int SPxDevexPR<R>::selectLeaveHyper(R feastol) 289 { 290 R x; 291 292 const R* fTest = this->thesolver->fTest().get_const_ptr(); 293 const R* cpen = this->thesolver->coWeights.get_const_ptr(); 294 R best = 0; 295 R leastBest = -1; 296 int bstI = -1; 297 int idx = -1; 298 299 // find the best price from the short candidate list 300 for(int i = bestPrices.size() - 1; i >= 0; --i) 301 { 302 idx = bestPrices.index(i); 303 x = fTest[idx]; 304 305 if(x < -feastol) 306 { 307 x = devexpr::computePrice(x, cpen[idx], feastol); 308 309 assert(x >= 0); 310 311 // update the best price of candidate list 312 if(x > best) 313 { 314 best = x; 315 bstI = idx; 316 last = cpen[idx]; 317 } 318 319 // update the smallest price of candidate list 320 if(x < leastBest || leastBest < 0) 321 leastBest = x; 322 } 323 else 324 { 325 bestPrices.remove(i); 326 this->thesolver->isInfeasible[idx] = this->NOT_VIOLATED; 327 } 328 } 329 330 // scan the updated indices for a better price 331 for(int i = this->thesolver->updateViols.size() - 1; i >= 0; --i) 332 { 333 idx = this->thesolver->updateViols.index(i); 334 335 // only look at indeces that were not checked already 336 if(this->thesolver->isInfeasible[idx] == this->VIOLATED) 337 { 338 x = fTest[idx]; 339 assert(x < -feastol); 340 x = devexpr::computePrice(x, cpen[idx], feastol); 341 342 if(x > leastBest) 343 { 344 if(x > best) 345 { 346 best = x; 347 bstI = idx; 348 last = cpen[idx]; 349 } 350 351 // put index into candidate list 352 this->thesolver->isInfeasible[idx] = this->VIOLATED_AND_CHECKED; 353 bestPrices.addIdx(idx); 354 } 355 } 356 } 357 358 return bstI; 359 } 360 361 template <class R> 362 void SPxDevexPR<R>::left4(int n, SPxId id) 363 { 364 VectorBase<R>& coWeights = this->thesolver->coWeights; 365 366 if(id.isValid()) 367 { 368 int i, j; 369 R x; 370 const R* rhoVec = this->thesolver->fVec().delta().values(); 371 R rhov_1 = 1 / rhoVec[n]; 372 R beta_q = this->thesolver->coPvec().delta().length2() * rhov_1 * rhov_1; 373 374 #ifndef NDEBUG 375 376 if(spxAbs(rhoVec[n]) < this->thetolerance) 377 { 378 SPX_MSG_INFO3((*this->thesolver->spxout), (*this->thesolver->spxout) << "WDEVEX01: rhoVec = " 379 << rhoVec[n] << " with smaller absolute value than this->thetolerance = " << this->thetolerance << 380 std::endl;) 381 } 382 383 #endif // NDEBUG 384 385 // Update #coPenalty# vector 386 const IdxSet& rhoIdx = this->thesolver->fVec().idx(); 387 int len = this->thesolver->fVec().idx().size(); 388 389 for(i = len - 1; i >= 0; --i) 390 { 391 j = rhoIdx.index(i); 392 x = rhoVec[j] * rhoVec[j] * beta_q; 393 // if(x > coPenalty[j]) 394 coWeights[j] += x; 395 } 396 397 coWeights[n] = beta_q; 398 } 399 } 400 401 template <class R> 402 SPxId SPxDevexPR<R>::buildBestPriceVectorEnterDim(R& best, R feastol) 403 { 404 int idx; 405 int nsorted; 406 R x; 407 const R* coTest = this->thesolver->coTest().get_const_ptr(); 408 const R* cpen = this->thesolver->coWeights.get_const_ptr(); 409 typename SPxPricer<R>::IdxElement price; 410 prices.clear(); 411 bestPrices.clear(); 412 413 // construct vector of all prices 414 for(int i = this->thesolver->infeasibilities.size() - 1; i >= 0; --i) 415 { 416 idx = this->thesolver->infeasibilities.index(i); 417 x = coTest[idx]; 418 419 if(x < -feastol) 420 { 421 this->thesolver->isInfeasible[idx] = this->VIOLATED; 422 price.idx = idx; 423 price.val = devexpr::computePrice(x, cpen[idx], feastol); 424 prices.append(price); 425 } 426 else 427 { 428 this->thesolver->infeasibilities.remove(i); 429 this->thesolver->isInfeasible[idx] = this->NOT_VIOLATED; 430 } 431 } 432 433 // set up structures for the quicksort implementation 434 this->compare.elements = prices.get_const_ptr(); 435 // do a partial sort to move the best ones to the front 436 // TODO this can be done more efficiently, since we only need the indices 437 nsorted = SPxQuicksortPart(prices.get_ptr(), this->compare, 0, prices.size(), 438 SOPLEX_HYPERPRICINGSIZE); 439 440 // copy indices of best values to bestPrices 441 for(int i = 0; i < nsorted; ++i) 442 { 443 bestPrices.addIdx(prices[i].idx); 444 this->thesolver->isInfeasible[prices[i].idx] = this->VIOLATED_AND_CHECKED; 445 } 446 447 if(nsorted > 0) 448 { 449 best = prices[0].val; 450 return this->thesolver->coId(prices[0].idx); 451 } 452 else 453 return SPxId(); 454 } 455 456 template <class R> 457 SPxId SPxDevexPR<R>::buildBestPriceVectorEnterCoDim(R& best, R feastol) 458 { 459 int idx; 460 int nsorted; 461 R x; 462 const R* test = this->thesolver->test().get_const_ptr(); 463 const R* pen = this->thesolver->weights.get_const_ptr(); 464 typename SPxPricer<R>::IdxElement price; 465 pricesCo.clear(); 466 bestPricesCo.clear(); 467 468 // construct vector of all prices 469 for(int i = this->thesolver->infeasibilitiesCo.size() - 1; i >= 0; --i) 470 { 471 idx = this->thesolver->infeasibilitiesCo.index(i); 472 x = test[idx]; 473 474 if(x < -feastol) 475 { 476 this->thesolver->isInfeasibleCo[idx] = this->VIOLATED; 477 price.idx = idx; 478 price.val = devexpr::computePrice(x, pen[idx], feastol); 479 pricesCo.append(price); 480 } 481 else 482 { 483 this->thesolver->infeasibilitiesCo.remove(i); 484 this->thesolver->isInfeasibleCo[idx] = this->NOT_VIOLATED; 485 } 486 } 487 488 // set up structures for the quicksort implementation 489 this->compare.elements = pricesCo.get_const_ptr(); 490 // do a partial sort to move the best ones to the front 491 // TODO this can be done more efficiently, since we only need the indices 492 nsorted = SPxQuicksortPart(pricesCo.get_ptr(), this->compare, 0, pricesCo.size(), 493 SOPLEX_HYPERPRICINGSIZE); 494 495 // copy indices of best values to bestPrices 496 for(int i = 0; i < nsorted; ++i) 497 { 498 bestPricesCo.addIdx(pricesCo[i].idx); 499 this->thesolver->isInfeasibleCo[pricesCo[i].idx] = this->VIOLATED_AND_CHECKED; 500 } 501 502 if(nsorted > 0) 503 { 504 best = pricesCo[0].val; 505 return this->thesolver->id(pricesCo[0].idx); 506 } 507 else 508 return SPxId(); 509 } 510 511 template <class R> 512 SPxId SPxDevexPR<R>::selectEnter() 513 { 514 assert(this->thesolver != 0); 515 516 SPxId enterId; 517 518 enterId = selectEnterX(this->thetolerance); 519 520 if(enterId.isSPxColId() && this->thesolver->isBasic(SPxColId(enterId))) 521 enterId.info = 0; 522 523 if(enterId.isSPxRowId() && this->thesolver->isBasic(SPxRowId(enterId))) 524 enterId.info = 0; 525 526 if(!enterId.isValid() && !refined) 527 { 528 refined = true; 529 SPX_MSG_INFO3((*this->thesolver->spxout), 530 (*this->thesolver->spxout) << "WDEVEX02 trying refinement step..\n";) 531 enterId = selectEnterX(this->thetolerance / SOPLEX_DEVEX_REFINETOL); 532 533 if(enterId.isSPxColId() && this->thesolver->isBasic(SPxColId(enterId))) 534 enterId.info = 0; 535 536 if(enterId.isSPxRowId() && this->thesolver->isBasic(SPxRowId(enterId))) 537 enterId.info = 0; 538 } 539 540 return enterId; 541 } 542 543 // choose the best entering index among columns and rows but prefer sparsity 544 template <class R> 545 SPxId SPxDevexPR<R>::selectEnterX(R tol) 546 { 547 SPxId enterId; 548 SPxId enterCoId; 549 R best; 550 R bestCo; 551 552 best = 0; 553 bestCo = 0; 554 last = 1.0; 555 556 // avoid uninitialized value later on in entered4X() 557 last = 1.0; 558 559 if(this->thesolver->hyperPricingEnter && !refined) 560 { 561 if(bestPrices.size() < 2 || this->thesolver->basis().lastUpdate() == 0) 562 enterCoId = (this->thesolver->sparsePricingEnter) ? buildBestPriceVectorEnterDim(best, 563 tol) : selectEnterDenseDim(best, tol); 564 else 565 enterCoId = (this->thesolver->sparsePricingEnter) ? selectEnterHyperDim(best, 566 tol) : selectEnterDenseDim(best, tol); 567 568 if(bestPricesCo.size() < 2 || this->thesolver->basis().lastUpdate() == 0) 569 enterId = (this->thesolver->sparsePricingEnterCo) ? buildBestPriceVectorEnterCoDim(bestCo, 570 tol) : selectEnterDenseCoDim(bestCo, tol); 571 else 572 enterId = (this->thesolver->sparsePricingEnterCo) ? selectEnterHyperCoDim(bestCo, 573 tol) : selectEnterDenseCoDim(bestCo, tol); 574 } 575 else 576 { 577 enterCoId = (this->thesolver->sparsePricingEnter 578 && !refined) ? selectEnterSparseDim(best, tol) : selectEnterDenseDim(best, tol); 579 enterId = (this->thesolver->sparsePricingEnterCo 580 && !refined) ? selectEnterSparseCoDim(bestCo, tol) : selectEnterDenseCoDim(bestCo, tol); 581 } 582 583 // prefer coIds to increase the number of unit vectors in the basis matrix, i.e., rows in colrep and cols in rowrep 584 if(enterCoId.isValid() && (best > SOPLEX_SPARSITY_TRADEOFF * bestCo || !enterId.isValid())) 585 return enterCoId; 586 else 587 return enterId; 588 } 589 590 template <class R> 591 SPxId SPxDevexPR<R>::selectEnterHyperDim(R& best, R feastol) 592 { 593 const R* cTest = this->thesolver->coTest().get_const_ptr(); 594 const R* cpen = this->thesolver->coWeights.get_const_ptr(); 595 R leastBest = -1; 596 R x; 597 int enterIdx = -1; 598 int idx; 599 600 // find the best price from short candidate list 601 for(int i = bestPrices.size() - 1; i >= 0; --i) 602 { 603 idx = bestPrices.index(i); 604 x = cTest[idx]; 605 606 if(x < -feastol) 607 { 608 x = devexpr::computePrice(x, cpen[idx], feastol); 609 610 assert(x >= 0); 611 612 // update the best price of candidate list 613 if(x > best) 614 { 615 best = x; 616 enterIdx = idx; 617 last = cpen[idx]; 618 } 619 620 // update the smallest price of candidate list 621 if(x < leastBest || leastBest < 0) 622 leastBest = x; 623 } 624 else 625 { 626 bestPrices.remove(i); 627 this->thesolver->isInfeasible[idx] = this->NOT_VIOLATED; 628 } 629 } 630 631 // scan the updated indeces for a better price 632 for(int i = this->thesolver->updateViols.size() - 1; i >= 0; --i) 633 { 634 idx = this->thesolver->updateViols.index(i); 635 636 // only look at indeces that were not checked already 637 if(this->thesolver->isInfeasible[idx] == this->VIOLATED) 638 { 639 x = cTest[idx]; 640 641 if(x < -feastol) 642 { 643 x = devexpr::computePrice(x, cpen[idx], feastol); 644 645 if(x > leastBest) 646 { 647 if(x > best) 648 { 649 best = x; 650 enterIdx = idx; 651 last = cpen[idx]; 652 } 653 654 // put index into candidate list 655 this->thesolver->isInfeasible[idx] = this->VIOLATED_AND_CHECKED; 656 bestPrices.addIdx(idx); 657 } 658 } 659 else 660 { 661 this->thesolver->isInfeasible[idx] = this->NOT_VIOLATED; 662 } 663 } 664 } 665 666 if(enterIdx >= 0) 667 return this->thesolver->coId(enterIdx); 668 else 669 return SPxId(); 670 } 671 672 template <class R> 673 SPxId SPxDevexPR<R>::selectEnterHyperCoDim(R& best, R feastol) 674 { 675 const R* test = this->thesolver->test().get_const_ptr(); 676 const R* pen = this->thesolver->weights.get_const_ptr(); 677 R leastBest = -1; 678 R x; 679 int enterIdx = -1; 680 int idx; 681 682 // find the best price from short candidate list 683 for(int i = bestPricesCo.size() - 1; i >= 0; --i) 684 { 685 idx = bestPricesCo.index(i); 686 x = test[idx]; 687 688 if(x < -feastol) 689 { 690 x = devexpr::computePrice(x, pen[idx], feastol); 691 692 assert(x >= 0); 693 694 // update the best price of candidate list 695 if(x > best) 696 { 697 best = x; 698 enterIdx = idx; 699 last = pen[idx]; 700 } 701 702 // update the smallest price of candidate list 703 if(x < leastBest || leastBest < 0) 704 leastBest = x; 705 } 706 else 707 { 708 bestPricesCo.remove(i); 709 this->thesolver->isInfeasibleCo[idx] = this->NOT_VIOLATED; 710 } 711 } 712 713 //scan the updated indeces for a better price 714 for(int i = this->thesolver->updateViolsCo.size() - 1; i >= 0; --i) 715 { 716 idx = this->thesolver->updateViolsCo.index(i); 717 718 // only look at indeces that were not checked already 719 if(this->thesolver->isInfeasibleCo[idx] == this->VIOLATED) 720 { 721 x = test[idx]; 722 723 if(x < -feastol) 724 { 725 x = devexpr::computePrice(x, pen[idx], feastol); 726 727 if(x > leastBest) 728 { 729 if(x > best) 730 { 731 best = x; 732 enterIdx = idx; 733 last = pen[idx]; 734 } 735 736 // put index into candidate list 737 this->thesolver->isInfeasibleCo[idx] = this->VIOLATED_AND_CHECKED; 738 bestPricesCo.addIdx(idx); 739 } 740 } 741 else 742 { 743 this->thesolver->isInfeasibleCo[idx] = this->NOT_VIOLATED; 744 } 745 } 746 } 747 748 if(enterIdx >= 0) 749 return this->thesolver->id(enterIdx); 750 else 751 return SPxId(); 752 } 753 754 755 template <class R> 756 SPxId SPxDevexPR<R>::selectEnterSparseDim(R& best, R feastol) 757 { 758 const R* cTest = this->thesolver->coTest().get_const_ptr(); 759 const R* cpen = this->thesolver->coWeights.get_const_ptr(); 760 int enterIdx = -1; 761 int idx; 762 R x; 763 764 assert(this->thesolver->coWeights.dim() == this->thesolver->coTest().dim()); 765 766 for(int i = this->thesolver->infeasibilities.size() - 1; i >= 0; --i) 767 { 768 idx = this->thesolver->infeasibilities.index(i); 769 x = cTest[idx]; 770 771 if(x < -feastol) 772 { 773 x = devexpr::computePrice(x, cpen[idx], feastol); 774 775 if(x > best) 776 { 777 best = x; 778 enterIdx = idx; 779 last = cpen[idx]; 780 } 781 } 782 else 783 { 784 this->thesolver->infeasibilities.remove(i); 785 this->thesolver->isInfeasible[idx] = this->NOT_VIOLATED; 786 } 787 } 788 789 if(enterIdx >= 0) 790 return this->thesolver->coId(enterIdx); 791 792 return SPxId(); 793 } 794 795 796 template <class R> 797 SPxId SPxDevexPR<R>::selectEnterSparseCoDim(R& best, R feastol) 798 { 799 const R* test = this->thesolver->test().get_const_ptr(); 800 const R* pen = this->thesolver->weights.get_const_ptr(); 801 int enterIdx = -1; 802 int idx; 803 R x; 804 805 assert(this->thesolver->weights.dim() == this->thesolver->test().dim()); 806 807 for(int i = this->thesolver->infeasibilitiesCo.size() - 1; i >= 0; --i) 808 { 809 idx = this->thesolver->infeasibilitiesCo.index(i); 810 x = test[idx]; 811 812 if(x < -feastol) 813 { 814 x = devexpr::computePrice(x, pen[idx], feastol); 815 816 if(x > best) 817 { 818 best = x; 819 enterIdx = idx; 820 last = pen[idx]; 821 } 822 } 823 else 824 { 825 this->thesolver->infeasibilitiesCo.remove(i); 826 this->thesolver->isInfeasibleCo[idx] = this->NOT_VIOLATED; 827 } 828 } 829 830 if(enterIdx >= 0) 831 return this->thesolver->id(enterIdx); 832 833 return SPxId(); 834 } 835 836 837 template <class R> 838 SPxId SPxDevexPR<R>::selectEnterDenseDim(R& best, R feastol, int start, int incr) 839 { 840 const R* cTest = this->thesolver->coTest().get_const_ptr(); 841 const R* cpen = this->thesolver->coWeights.get_const_ptr(); 842 int end = this->thesolver->coWeights.dim(); 843 int enterIdx = -1; 844 R x; 845 846 assert(end == this->thesolver->coTest().dim()); 847 848 for(; start < end; start += incr) 849 { 850 x = cTest[start]; 851 852 if(x < -feastol) 853 { 854 x = devexpr::computePrice(x, cpen[start], feastol); 855 856 if(x > best) 857 { 858 best = x; 859 enterIdx = start; 860 last = cpen[start]; 861 } 862 } 863 } 864 865 if(enterIdx >= 0) 866 return this->thesolver->coId(enterIdx); 867 868 return SPxId(); 869 } 870 871 template <class R> 872 SPxId SPxDevexPR<R>::selectEnterDenseCoDim(R& best, R feastol, int start, int incr) 873 { 874 const R* test = this->thesolver->test().get_const_ptr(); 875 const R* pen = this->thesolver->weights.get_const_ptr(); 876 int end = this->thesolver->weights.dim(); 877 int enterIdx = -1; 878 R x; 879 880 assert(end == this->thesolver->test().dim()); 881 882 for(; start < end; start += incr) 883 { 884 x = test[start]; 885 886 if(test[start] < -feastol) 887 { 888 x = devexpr::computePrice(x, pen[start], feastol); 889 890 if(x > best) 891 { 892 best = x; 893 enterIdx = start; 894 last = pen[start]; 895 } 896 } 897 } 898 899 if(enterIdx >= 0) 900 return this->thesolver->id(enterIdx); 901 902 return SPxId(); 903 } 904 905 906 /**@todo suspicious: the pricer should be informed, that variable id 907 has entered the basis at position n, but the id is not used here 908 (this is true for all pricers) 909 */ 910 template <class R> 911 void SPxDevexPR<R>::entered4(SPxId /*id*/, int n) 912 { 913 VectorBase<R>& weights = this->thesolver->weights; 914 VectorBase<R>& coWeights = this->thesolver->coWeights; 915 916 if(n >= 0 && n < this->thesolver->dim()) 917 { 918 const R* pVec = this->thesolver->pVec().delta().values(); 919 const IdxSet& pIdx = this->thesolver->pVec().idx(); 920 const R* coPvec = this->thesolver->coPvec().delta().values(); 921 const IdxSet& coPidx = this->thesolver->coPvec().idx(); 922 R xi_p = 1 / this->thesolver->fVec().delta()[n]; 923 int i, j; 924 925 assert(this->thesolver->fVec().delta()[n] > this->thesolver->epsilon() 926 || this->thesolver->fVec().delta()[n] < -this->thesolver->epsilon()); 927 928 xi_p = xi_p * xi_p * last; 929 930 for(j = coPidx.size() - 1; j >= 0; --j) 931 { 932 i = coPidx.index(j); 933 coWeights[i] += xi_p * coPvec[i] * coPvec[i]; 934 935 if(coWeights[i] <= 1 || coWeights[i] > 1e+6) 936 { 937 setupWeights(SPxSolverBase<R>::ENTER); 938 return; 939 } 940 } 941 942 for(j = pIdx.size() - 1; j >= 0; --j) 943 { 944 i = pIdx.index(j); 945 weights[i] += xi_p * pVec[i] * pVec[i]; 946 947 if(weights[i] <= 1 || weights[i] > 1e+6) 948 { 949 setupWeights(SPxSolverBase<R>::ENTER); 950 return; 951 } 952 } 953 } 954 } 955 956 template <class R> 957 void SPxDevexPR<R>::addedVecs(int n) 958 { 959 int initval = (this->thesolver->type() == SPxSolverBase<R>::ENTER) ? 2 : 1; 960 VectorBase<R>& weights = this->thesolver->weights; 961 n = weights.dim(); 962 weights.reDim(this->thesolver->coDim()); 963 964 for(int i = weights.dim() - 1; i >= n; --i) 965 weights[i] = initval; 966 } 967 968 template <class R> 969 void SPxDevexPR<R>::addedCoVecs(int n) 970 { 971 int initval = (this->thesolver->type() == SPxSolverBase<R>::ENTER) ? 2 : 1; 972 VectorBase<R>& coWeights = this->thesolver->coWeights; 973 n = coWeights.dim(); 974 coWeights.reDim(this->thesolver->dim()); 975 976 for(int i = coWeights.dim() - 1; i >= n; --i) 977 coWeights[i] = initval; 978 } 979 980 } // namespace soplex 981