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 <iostream> 27 28 #include "soplex/spxdefines.h" 29 #include "soplex/spxsolver.h" 30 #include "soplex/exceptions.h" 31 32 namespace soplex 33 { 34 /** Initialize Vectors 35 36 Computing the right hand side vector for the feasibility vector depends on 37 the chosen representation and type of the basis. 38 39 In columnwise case, |theFvec| = \f$x_B = A_B^{-1} (- A_N x_N)\f$, where \f$x_N\f$ 40 are either upper or lower bounds for the nonbasic variables (depending on 41 the variables |Status|). If these values remain unchanged throughout the 42 simplex algorithm, they may be taken directly from LP. However, in the 43 entering type algorith they are changed and, hence, retreived from the 44 column or row upper or lower bound vectors. 45 46 In rowwise case, |theFvec| = \f$\pi^T_B = (c^T - 0^T A_N) A_B^{-1}\f$. However, 47 this applies only to leaving type algorithm, where no bounds on dual 48 variables are altered. In entering type algorithm they are changed and, 49 hence, retreived from the column or row upper or lower bound vectors. 50 */ 51 template <class R> 52 void SPxSolverBase<R>::computeFrhs() 53 { 54 55 if(rep() == COLUMN) 56 { 57 theFrhs->clear(); 58 59 if(type() == LEAVE) 60 { 61 computeFrhsXtra(); 62 63 for(int i = 0; i < this->nRows(); i++) 64 { 65 R x; 66 67 typename SPxBasisBase<R>::Desc::Status stat = this->desc().rowStatus(i); 68 69 if(!isBasic(stat)) 70 { 71 // coverity[switch_selector_expr_is_constant] 72 switch(stat) 73 { 74 // columnwise cases: 75 case SPxBasisBase<R>::Desc::P_FREE : 76 continue; 77 78 case(SPxBasisBase<R>::Desc::P_FIXED) : 79 assert(EQ(this->lhs(i), this->rhs(i), this->epsilon())); 80 81 //lint -fallthrough 82 case SPxBasisBase<R>::Desc::P_ON_UPPER : 83 x = this->rhs(i); 84 break; 85 86 case SPxBasisBase<R>::Desc::P_ON_LOWER : 87 x = this->lhs(i); 88 break; 89 90 default: 91 SPX_MSG_ERROR(std::cerr << "ESVECS01 ERROR: " 92 << "inconsistent basis must not happen!" 93 << std::endl;) 94 throw SPxInternalCodeException("XSVECS01 This should never happen."); 95 } 96 97 assert(x < R(infinity)); 98 assert(x > R(-infinity)); 99 (*theFrhs)[i] += x; // slack ! 100 } 101 } 102 } 103 else 104 { 105 computeFrhs1(*theUbound, *theLbound); 106 computeFrhs2(*theCoUbound, *theCoLbound); 107 } 108 } 109 else 110 { 111 assert(rep() == ROW); 112 113 if(type() == ENTER) 114 { 115 theFrhs->clear(); 116 computeFrhs1(*theUbound, *theLbound); 117 computeFrhs2(*theCoUbound, *theCoLbound); 118 *theFrhs += this->maxObj(); 119 } 120 else 121 { 122 ///@todo put this into a separate method 123 *theFrhs = this->maxObj(); 124 const typename SPxBasisBase<R>::Desc& ds = this->desc(); 125 126 for(int i = 0; i < this->nRows(); ++i) 127 { 128 typename SPxBasisBase<R>::Desc::Status stat = ds.rowStatus(i); 129 130 if(!isBasic(stat)) 131 { 132 R x; 133 134 switch(stat) 135 { 136 case SPxBasisBase<R>::Desc::D_FREE : 137 continue; 138 139 case SPxBasisBase<R>::Desc::D_ON_UPPER : 140 case SPxBasisBase<R>::Desc::D_ON_LOWER : 141 case(SPxBasisBase<R>::Desc::D_ON_BOTH) : 142 x = this->maxRowObj(i); 143 break; 144 145 default: 146 assert(this->lhs(i) <= R(-infinity) && this->rhs(i) >= R(infinity)); 147 x = 0.0; 148 break; 149 } 150 151 assert(x < R(infinity)); 152 assert(x > R(-infinity)); 153 // assert(x == 0.0); 154 155 if(x != 0.0) 156 theFrhs->multAdd(x, vector(i)); 157 } 158 } 159 } 160 } 161 } 162 163 template <class R> 164 void SPxSolverBase<R>::computeFrhsXtra() 165 { 166 167 assert(rep() == COLUMN); 168 assert(type() == LEAVE); 169 170 for(int i = 0; i < this->nCols(); ++i) 171 { 172 typename SPxBasisBase<R>::Desc::Status stat = this->desc().colStatus(i); 173 174 if(!isBasic(stat)) 175 { 176 R x; 177 178 // coverity[switch_selector_expr_is_constant] 179 switch(stat) 180 { 181 // columnwise cases: 182 case SPxBasisBase<R>::Desc::P_FREE : 183 continue; 184 185 case(SPxBasisBase<R>::Desc::P_FIXED) : 186 assert(EQ(SPxLPBase<R>::lower(i), SPxLPBase<R>::upper(i), this->epsilon())); 187 188 //lint -fallthrough 189 case SPxBasisBase<R>::Desc::P_ON_UPPER : 190 x = SPxLPBase<R>::upper(i); 191 break; 192 193 case SPxBasisBase<R>::Desc::P_ON_LOWER : 194 x = SPxLPBase<R>::lower(i); 195 break; 196 197 default: 198 SPX_MSG_ERROR(std::cerr << "ESVECS02 ERROR: " 199 << "inconsistent basis must not happen!" 200 << std::endl;) 201 throw SPxInternalCodeException("XSVECS02 This should never happen."); 202 } 203 204 assert(x < R(infinity)); 205 assert(x > R(-infinity)); 206 207 if(x != 0.0) 208 theFrhs->multAdd(-x, vector(i)); 209 } 210 } 211 } 212 213 214 /** This methods subtracts \f$A_N x_N\f$ or \f$\pi_N^T A_N\f$ from |theFrhs| as 215 specified by the |Status| of all nonbasic variables. The values of \f$x_N\f$ or 216 \f$\pi_N\f$ are taken from the passed arrays. 217 */ 218 template <class R> 219 void SPxSolverBase<R>::computeFrhs1( 220 const VectorBase<R>& ufb, ///< upper feasibility bound for variables 221 const VectorBase<R>& lfb) ///< lower feasibility bound for variables 222 { 223 224 const typename SPxBasisBase<R>::Desc& ds = this->desc(); 225 226 for(int i = 0; i < coDim(); ++i) 227 { 228 typename SPxBasisBase<R>::Desc::Status stat = ds.status(i); 229 230 if(!isBasic(stat)) 231 { 232 R x; 233 234 // coverity[switch_selector_expr_is_constant] 235 switch(stat) 236 { 237 case SPxBasisBase<R>::Desc::D_FREE : 238 case SPxBasisBase<R>::Desc::D_UNDEFINED : 239 case SPxBasisBase<R>::Desc::P_FREE : 240 continue; 241 242 case SPxBasisBase<R>::Desc::P_ON_UPPER : 243 case SPxBasisBase<R>::Desc::D_ON_UPPER : 244 x = ufb[i]; 245 break; 246 247 case SPxBasisBase<R>::Desc::P_ON_LOWER : 248 case SPxBasisBase<R>::Desc::D_ON_LOWER : 249 x = lfb[i]; 250 break; 251 252 case(SPxBasisBase<R>::Desc::P_FIXED) : 253 assert(EQ(lfb[i], ufb[i], this->epsilon())); 254 255 //lint -fallthrough 256 case(SPxBasisBase<R>::Desc::D_ON_BOTH) : 257 x = lfb[i]; 258 break; 259 260 default: 261 SPX_MSG_ERROR(std::cerr << "ESVECS03 ERROR: " 262 << "inconsistent basis must not happen!" 263 << std::endl;) 264 throw SPxInternalCodeException("XSVECS04 This should never happen."); 265 } 266 267 assert(x < R(infinity)); 268 assert(x > R(-infinity)); 269 270 if(x != 0.0) 271 theFrhs->multAdd(-x, vector(i)); 272 } 273 } 274 } 275 276 /** This methods subtracts \f$A_N x_N\f$ or \f$\pi_N^T A_N\f$ from |theFrhs| as 277 specified by the |Status| of all nonbasic variables. The values of 278 \f$x_N\f$ or \f$\pi_N\f$ are taken from the passed arrays. 279 */ 280 template <class R> 281 void SPxSolverBase<R>::computeFrhs2( 282 VectorBase<R>& coufb, ///< upper feasibility bound for covariables 283 VectorBase<R>& colfb) ///< lower feasibility bound for covariables 284 { 285 const typename SPxBasisBase<R>::Desc& ds = this->desc(); 286 287 for(int i = 0; i < dim(); ++i) 288 { 289 typename SPxBasisBase<R>::Desc::Status stat = ds.coStatus(i); 290 291 if(!isBasic(stat)) 292 { 293 R x; 294 295 // coverity[switch_selector_expr_is_constant] 296 switch(stat) 297 { 298 case SPxBasisBase<R>::Desc::D_FREE : 299 case SPxBasisBase<R>::Desc::D_UNDEFINED : 300 case SPxBasisBase<R>::Desc::P_FREE : 301 continue; 302 303 case SPxBasisBase<R>::Desc::P_ON_LOWER : // negative slack bounds! 304 case SPxBasisBase<R>::Desc::D_ON_UPPER : 305 x = coufb[i]; 306 break; 307 308 case SPxBasisBase<R>::Desc::P_ON_UPPER : // negative slack bounds! 309 case SPxBasisBase<R>::Desc::D_ON_LOWER : 310 x = colfb[i]; 311 break; 312 313 case SPxBasisBase<R>::Desc::P_FIXED : 314 case SPxBasisBase<R>::Desc::D_ON_BOTH : 315 316 if(colfb[i] != coufb[i]) 317 { 318 SPX_MSG_WARNING((*this->spxout), 319 (*this->spxout) << "WSVECS04 Frhs2[" << i << "]: " << static_cast<int> 320 (stat) << " " 321 << colfb[i] << " " << coufb[i] 322 << " shouldn't be" << std::endl;) 323 324 if(isZero(colfb[i], this->epsilon()) || isZero(coufb[i], this->epsilon())) 325 colfb[i] = coufb[i] = 0.0; 326 else 327 { 328 R mid = (colfb[i] + coufb[i]) / 2.0; 329 colfb[i] = coufb[i] = mid; 330 } 331 } 332 333 assert(EQ(colfb[i], coufb[i], this->epsilon())); 334 x = colfb[i]; 335 break; 336 337 default: 338 SPX_MSG_ERROR(std::cerr << "ESVECS05 ERROR: " 339 << "inconsistent basis must not happen!" 340 << std::endl;) 341 throw SPxInternalCodeException("XSVECS05 This should never happen."); 342 } 343 344 assert(x < R(infinity)); 345 assert(x > R(-infinity)); 346 347 (*theFrhs)[i] -= x; // This is a slack, so no need to multiply a vector. 348 } 349 } 350 } 351 352 /** Computing the right hand side vector for |theCoPvec| depends on 353 the type of the simplex algorithm. In entering algorithms, the 354 values are taken from the inequality's right handside or the 355 column's objective value. 356 357 In contrast to this leaving algorithms take the values from vectors 358 |theURbound| and so on. 359 360 We reflect this difference by providing two pairs of methods 361 |computeEnterCoPrhs(n, stat)| and |computeLeaveCoPrhs(n, stat)|. The first 362 pair operates for entering algorithms, while the second one is intended for 363 leaving algorithms. The return value of these methods is the right hand 364 side value for the \f$n\f$-th row or column id, respectively, if it had the 365 passed |Status| for both. 366 367 Both methods are again split up into two methods named |...4Row(i,n)| and 368 |...4Col(i,n)|, respectively. They do their job for the |i|-th basis 369 variable, being the |n|-th row or column. 370 */ 371 template <class R> 372 void SPxSolverBase<R>::computeEnterCoPrhs4Row(int i, int n) 373 { 374 assert(this->baseId(i).isSPxRowId()); 375 assert(this->number(SPxRowId(this->baseId(i))) == n); 376 377 switch(this->desc().rowStatus(n)) 378 { 379 // rowwise representation: 380 case SPxBasisBase<R>::Desc::P_FIXED : 381 assert(this->lhs(n) > R(-infinity)); 382 assert(EQ(this->rhs(n), this->lhs(n), this->epsilon())); 383 384 //lint -fallthrough 385 case SPxBasisBase<R>::Desc::P_ON_UPPER : 386 assert(rep() == ROW); 387 assert(this->rhs(n) < R(infinity)); 388 (*theCoPrhs)[i] = this->rhs(n); 389 break; 390 391 case SPxBasisBase<R>::Desc::P_ON_LOWER : 392 assert(rep() == ROW); 393 assert(this->lhs(n) > R(-infinity)); 394 (*theCoPrhs)[i] = this->lhs(n); 395 break; 396 397 // columnwise representation: 398 // slacks must be left 0! 399 default: 400 (*theCoPrhs)[i] = this->maxRowObj(n); 401 break; 402 } 403 } 404 405 template <class R> 406 void SPxSolverBase<R>::computeEnterCoPrhs4Col(int i, int n) 407 { 408 assert(this->baseId(i).isSPxColId()); 409 assert(this->number(SPxColId(this->baseId(i))) == n); 410 411 switch(this->desc().colStatus(n)) 412 { 413 // rowwise representation: 414 case SPxBasisBase<R>::Desc::P_FIXED : 415 assert(EQ(SPxLPBase<R>::upper(n), SPxLPBase<R>::lower(n), this->epsilon())); 416 assert(SPxLPBase<R>::lower(n) > R(-infinity)); 417 418 //lint -fallthrough 419 case SPxBasisBase<R>::Desc::P_ON_UPPER : 420 assert(rep() == ROW); 421 assert(SPxLPBase<R>::upper(n) < R(infinity)); 422 (*theCoPrhs)[i] = SPxLPBase<R>::upper(n); 423 break; 424 425 case SPxBasisBase<R>::Desc::P_ON_LOWER : 426 assert(rep() == ROW); 427 assert(SPxLPBase<R>::lower(n) > R(-infinity)); 428 (*theCoPrhs)[i] = SPxLPBase<R>::lower(n); 429 break; 430 431 // columnwise representation: 432 case SPxBasisBase<R>::Desc::D_UNDEFINED : 433 case SPxBasisBase<R>::Desc::D_ON_BOTH : 434 case SPxBasisBase<R>::Desc::D_ON_UPPER : 435 case SPxBasisBase<R>::Desc::D_ON_LOWER : 436 case SPxBasisBase<R>::Desc::D_FREE : 437 assert(rep() == COLUMN); 438 (*theCoPrhs)[i] = this->maxObj(n); 439 break; 440 441 default: // variable left 0 442 (*theCoPrhs)[i] = 0; 443 break; 444 } 445 } 446 447 template <class R> 448 void SPxSolverBase<R>::computeEnterCoPrhs() 449 { 450 assert(type() == ENTER); 451 452 for(int i = dim() - 1; i >= 0; --i) 453 { 454 SPxId l_id = this->baseId(i); 455 456 if(l_id.isSPxRowId()) 457 computeEnterCoPrhs4Row(i, this->number(SPxRowId(l_id))); 458 else 459 computeEnterCoPrhs4Col(i, this->number(SPxColId(l_id))); 460 } 461 } 462 463 template <class R> 464 void SPxSolverBase<R>::computeLeaveCoPrhs4Row(int i, int n) 465 { 466 assert(this->baseId(i).isSPxRowId()); 467 assert(this->number(SPxRowId(this->baseId(i))) == n); 468 469 switch(this->desc().rowStatus(n)) 470 { 471 case SPxBasisBase<R>::Desc::D_ON_BOTH : 472 case SPxBasisBase<R>::Desc::P_FIXED : 473 assert(theLRbound[n] > R(-infinity)); 474 assert(EQ(theURbound[n], theLRbound[n], this->epsilon())); 475 476 //lint -fallthrough 477 case SPxBasisBase<R>::Desc::D_ON_UPPER : 478 case SPxBasisBase<R>::Desc::P_ON_UPPER : 479 assert(theURbound[n] < R(infinity)); 480 (*theCoPrhs)[i] = theURbound[n]; 481 break; 482 483 case SPxBasisBase<R>::Desc::D_ON_LOWER : 484 case SPxBasisBase<R>::Desc::P_ON_LOWER : 485 assert(theLRbound[n] > R(-infinity)); 486 (*theCoPrhs)[i] = theLRbound[n]; 487 break; 488 489 default: 490 (*theCoPrhs)[i] = this->maxRowObj(n); 491 break; 492 } 493 } 494 495 template <class R> 496 void SPxSolverBase<R>::computeLeaveCoPrhs4Col(int i, int n) 497 { 498 assert(this->baseId(i).isSPxColId()); 499 assert(this->number(SPxColId(this->baseId(i))) == n); 500 501 switch(this->desc().colStatus(n)) 502 { 503 case SPxBasisBase<R>::Desc::D_UNDEFINED : 504 case SPxBasisBase<R>::Desc::D_ON_BOTH : 505 case SPxBasisBase<R>::Desc::P_FIXED : 506 assert(theLCbound[n] > R(-infinity)); 507 assert(EQ(theUCbound[n], theLCbound[n], this->epsilon())); 508 509 //lint -fallthrough 510 case SPxBasisBase<R>::Desc::D_ON_LOWER : 511 case SPxBasisBase<R>::Desc::P_ON_UPPER : 512 assert(theUCbound[n] < R(infinity)); 513 (*theCoPrhs)[i] = theUCbound[n]; 514 break; 515 516 case SPxBasisBase<R>::Desc::D_ON_UPPER : 517 case SPxBasisBase<R>::Desc::P_ON_LOWER : 518 assert(theLCbound[n] > R(-infinity)); 519 (*theCoPrhs)[i] = theLCbound[n]; 520 break; 521 522 default: 523 (*theCoPrhs)[i] = this->maxObj(n); 524 // (*theCoPrhs)[i] = 0; 525 break; 526 } 527 } 528 529 template <class R> 530 void SPxSolverBase<R>::computeLeaveCoPrhs() 531 { 532 assert(type() == LEAVE); 533 534 for(int i = dim() - 1; i >= 0; --i) 535 { 536 SPxId l_id = this->baseId(i); 537 538 if(l_id.isSPxRowId()) 539 computeLeaveCoPrhs4Row(i, this->number(SPxRowId(l_id))); 540 else 541 computeLeaveCoPrhs4Col(i, this->number(SPxColId(l_id))); 542 } 543 } 544 545 546 /** When computing the copricing vector, we expect the pricing vector to be 547 setup correctly. Then computing the copricing vector is nothing but 548 computing all scalar products of the pricing vector with the vectors of the 549 LPs constraint matrix. 550 */ 551 template <class R> 552 void SPxSolverBase<R>::computePvec() 553 { 554 int i; 555 556 for(i = coDim() - 1; i >= 0; --i) 557 (*thePvec)[i] = vector(i) * (*theCoPvec); 558 } 559 560 template <class R> 561 void SPxSolverBase<R>::setupPupdate(void) 562 { 563 SSVectorBase<R>& p = thePvec->delta(); 564 SSVectorBase<R>& c = theCoPvec->delta(); 565 566 if(c.isSetup()) 567 { 568 if(c.size() < 0.95 * theCoPvec->dim()) 569 p.assign2product4setup(*thecovectors, c, 570 multTimeSparse, multTimeFull, 571 multSparseCalls, multFullCalls); 572 else 573 { 574 multTimeColwise->start(); 575 p.assign2product(c, *thevectors); 576 multTimeColwise->stop(); 577 ++multColwiseCalls; 578 } 579 } 580 else 581 { 582 multTimeUnsetup->start(); 583 p.assign2productAndSetup(*thecovectors, c); 584 multTimeUnsetup->stop(); 585 ++multUnsetupCalls; 586 } 587 588 p.setup(); 589 } 590 591 template <class R> 592 void SPxSolverBase<R>::doPupdate(void) 593 { 594 theCoPvec->update(); 595 596 if(pricing() == FULL) 597 { 598 // thePvec->delta().setup(); 599 thePvec->update(); 600 } 601 } 602 } // namespace soplex 603