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 <iostream> 26 #include <assert.h> 27 28 #include "soplex/spxdefines.h" 29 #include "soplex.h" 30 #include "soplex/statistics.h" 31 #include "soplex/sorter.h" 32 33 #define SOPLEX_USE_FEASTOL 34 #define SOPLEX_MAX_DEGENCHECK 20 /**< the maximum number of degen checks that are performed before the DECOMP is abandoned */ 35 #define SOPLEX_DEGENCHECK_OFFSET 50 /**< the number of iteration before the degeneracy check is reperformed */ 36 #define SOPLEX_SLACKCOEFF 1.0 /**< the coefficient of the slack variable in the incompatible rows. */ 37 #define SOPLEX_TIMELIMIT_FRAC 0.5 /**< the fraction of the total time limit given to the setup of the reduced problem */ 38 39 /* This file contains the private functions for the Decomposition Based Dual Simplex (DBDS) 40 * 41 * An important note about the implementation of the DBDS is the reliance on the row representation of the basis matrix. 42 * The forming of the reduced and complementary problems is involves identifying rows in the row-form of the basis 43 * matrix that have zero dual multipliers. 44 * 45 * Ideally, this work will be extended such that the DBDS can be executed using the column-form of the basis matrix. */ 46 47 namespace soplex 48 { 49 50 /// solves LP using the decomposition dual simplex 51 template <class R> 52 void SoPlexBase<R>::_solveDecompositionDualSimplex() 53 { 54 assert(_solver.rep() == SPxSolverBase<R>::ROW); 55 assert(_solver.type() == SPxSolverBase<R>::LEAVE); 56 57 // flag to indicate that the algorithm must terminate 58 bool stop = false; 59 60 // setting the initial status of the reduced problem 61 _statistics->redProbStatus = SPxSolverBase<R>::NO_PROBLEM; 62 63 // start timing 64 _statistics->solvingTime->start(); 65 66 // remember that last solve was in floating-point 67 _lastSolveMode = SOLVEMODE_REAL; 68 69 // setting the current solving mode. 70 _currentProb = DECOMP_ORIG; 71 72 // setting the sense to maximise. This is to make all matrix operations more consistent. 73 // @todo if the objective sense is changed, the output of the current objective value is the negative of the 74 // actual value. This needs to be corrected in future versions. 75 if(intParam(SoPlexBase<R>::OBJSENSE) == SoPlexBase<R>::OBJSENSE_MINIMIZE) 76 { 77 assert(_solver.spxSense() == SPxLPBase<R>::MINIMIZE); 78 79 _solver.changeObj(_solver.maxObj()); 80 _solver.changeSense(SPxLPBase<R>::MAXIMIZE); 81 } 82 83 // it is necessary to solve the initial problem to find a starting basis 84 _solver.setDecompStatus(SPxSolverBase<R>::FINDSTARTBASIS); 85 86 // setting the decomposition iteration limit to the parameter setting 87 _solver.setDecompIterationLimit(intParam(SoPlexBase<R>::DECOMP_ITERLIMIT)); 88 89 // variables used in the initialisation phase of the decomposition solve. 90 int numDegenCheck = 0; 91 R degeneracyLevel = 0; 92 _decompFeasVector.reDim(_solver.nCols()); 93 bool initSolveFromScratch = true; 94 95 96 /************************/ 97 /* Initialisation phase */ 98 /************************/ 99 100 // arrays to store the basis status for the rows and columns at the interruption of the original problem solve. 101 DataArray< typename SPxSolverBase<R>::VarStatus > basisStatusRows; 102 DataArray< typename SPxSolverBase<R>::VarStatus > basisStatusCols; 103 104 // since the original LP may have been shifted, the dual multiplier will not be correct for the original LP. This 105 // loop will recheck the degeneracy level and compute the proper dual multipliers. 106 do 107 { 108 // solving the instance. 109 _decompSimplifyAndSolve(_solver, _slufactor, initSolveFromScratch, initSolveFromScratch); 110 initSolveFromScratch = false; 111 112 // checking whether the initialisation must terminate and the original problem is solved using the dual simplex. 113 if(_solver.type() == SPxSolverBase<R>::LEAVE || _solver.status() >= SPxSolverBase<R>::OPTIMAL 114 || _solver.status() == SPxSolverBase<R>::ABORT_EXDECOMP || numDegenCheck > SOPLEX_MAX_DEGENCHECK) 115 { 116 // decomposition is deemed not useful. Solving the original problem using regular SoPlexBase. 117 118 // returning the sense to minimise 119 if(intParam(SoPlexBase<R>::OBJSENSE) == SoPlexBase<R>::OBJSENSE_MINIMIZE) 120 { 121 assert(_solver.spxSense() == SPxLPBase<R>::MAXIMIZE); 122 123 _solver.changeObj(-(_solver.maxObj())); 124 _solver.changeSense(SPxLPBase<R>::MINIMIZE); 125 } 126 127 // switching off the starting basis check. This is only required to initialise the decomposition simplex. 128 _solver.setDecompStatus(SPxSolverBase<R>::DONTFINDSTARTBASIS); 129 130 // the basis is not computed correctly is the problem was unfeasible or unbounded. 131 if(_solver.status() == SPxSolverBase<R>::UNBOUNDED 132 || _solver.status() == SPxSolverBase<R>::INFEASIBLE) 133 _hasBasis = false; 134 135 136 // resolving the problem to update the real lp and solve with the correct objective. 137 // TODO: With some infeasible problem (e.g. refinery) the dual is violated. Setting fromScratch to true 138 // avoids this issue. Need to check what the problem is. 139 // @TODO: Should this be _preprocessAndSolveReal similar to the resolve at the end of the algorithm? 140 _decompSimplifyAndSolve(_solver, _slufactor, true, false); 141 142 // retreiving the original problem statistics prior to destroying it. 143 getOriginalProblemStatistics(); 144 145 // storing the solution from the reduced problem 146 _storeSolutionReal(); 147 148 // stop timing 149 _statistics->solvingTime->stop(); 150 return; 151 } 152 else if(_solver.status() == SPxSolverBase<R>::ABORT_TIME 153 || _solver.status() == SPxSolverBase<R>::ABORT_ITER 154 || _solver.status() == SPxSolverBase<R>::ABORT_VALUE) 155 { 156 // This cleans up the problem is an abort is reached. 157 158 // at this point, the _decompSimplifyAndSolve does not store the realLP. It stores the _decompLP. As such, it is 159 // important to reinstall the _realLP to the _solver. 160 if(!_isRealLPLoaded) 161 { 162 _solver.loadLP(*_decompLP); 163 spx_free(_decompLP); 164 _decompLP = &_solver; 165 _isRealLPLoaded = true; 166 } 167 168 // returning the sense to minimise 169 if(intParam(SoPlexBase<R>::OBJSENSE) == SoPlexBase<R>::OBJSENSE_MINIMIZE) 170 { 171 assert(_solver.spxSense() == SPxLPBase<R>::MAXIMIZE); 172 173 _solver.changeObj(-(_solver.maxObj())); 174 _solver.changeSense(SPxLPBase<R>::MINIMIZE); 175 } 176 177 // resolving the problem to set up all of the solution vectors. This is required because data from the 178 // initialisation solve may remain at termination causing an infeasible solution to be reported. 179 _preprocessAndSolveReal(false); 180 181 // storing the solution from the reduced problem 182 _storeSolutionReal(); 183 184 // stop timing 185 _statistics->solvingTime->stop(); 186 return; 187 } 188 189 // checking the degeneracy level 190 _solver.basis().solve(_decompFeasVector, _solver.maxObj()); 191 degeneracyLevel = _solver.getDegeneracyLevel(_decompFeasVector); 192 193 _solver.setDegenCompOffset(SOPLEX_DEGENCHECK_OFFSET); 194 195 numDegenCheck++; 196 } 197 while((degeneracyLevel > 0.9 || degeneracyLevel < 0.1) 198 || !checkBasisDualFeasibility(_decompFeasVector)); 199 200 // decomposition simplex will commence, so the start basis does not need to be found 201 _solver.setDecompStatus(SPxSolverBase<R>::DONTFINDSTARTBASIS); 202 203 // if the original problem has a basis, this will be stored 204 if(_hasBasis) 205 { 206 basisStatusRows.reSize(numRows()); 207 basisStatusCols.reSize(numCols()); 208 _solver.getBasis(basisStatusRows.get_ptr(), basisStatusCols.get_ptr(), basisStatusRows.size(), 209 basisStatusCols.size()); 210 } 211 212 // updating the algorithm iterations statistic 213 _statistics->callsReducedProb++; 214 215 // setting the statistic information 216 numDecompIter = 0; 217 numRedProbIter = _solver.iterations(); 218 numCompProbIter = 0; 219 220 // setting the display information 221 _decompDisplayLine = 0; 222 223 /************************/ 224 /* Decomposition phase */ 225 /************************/ 226 227 SPX_MSG_INFO1(spxout, 228 spxout << "======== Degeneracy Detected ========" << std::endl 229 << std::endl 230 << "======== Commencing decomposition solve ========" << std::endl 231 ); 232 233 //spxout.setVerbosity( SPxOut::DEBUG ); 234 SPX_MSG_INFO2(spxout, spxout << "Creating the Reduced and Complementary problems." << std::endl); 235 236 // setting the verbosity level 237 const SPxOut::Verbosity orig_verbosity = spxout.getVerbosity(); 238 const SPxOut::Verbosity decomp_verbosity = (SPxOut::Verbosity)intParam( 239 SoPlexBase<R>::DECOMP_VERBOSITY); 240 241 if(decomp_verbosity < orig_verbosity) 242 spxout.setVerbosity(decomp_verbosity); 243 244 // creating copies of the original problem that will be manipulated to form the reduced and complementary 245 // problems. 246 _createDecompReducedAndComplementaryProblems(); 247 248 // creating the initial reduced problem from the basis information 249 _formDecompReducedProblem(stop); 250 251 // setting flags for the decomposition solve 252 _hasBasis = false; 253 bool hasRedBasis = false; 254 bool redProbError = false; 255 bool noRedprobIter = false; 256 bool explicitviol = boolParam(SoPlexBase<R>::EXPLICITVIOL); 257 int algIterCount = 0; 258 259 260 // a stop will be triggered if the reduced problem exceeded a specified fraction of the time limit 261 if(stop) 262 { 263 SPX_MSG_INFO1(spxout, spxout << "==== Error constructing the reduced problem ====" << std::endl); 264 redProbError = true; 265 _statistics->redProbStatus = SPxSolverBase<R>::NOT_INIT; 266 } 267 268 // the main solving loop of the decomposition simplex. 269 // This loop solves the Reduced problem, and if the problem is feasible, the complementary problem is solved. 270 while(!stop) 271 { 272 int previter = _statistics->iterations; 273 274 // setting the current solving mode. 275 _currentProb = DECOMP_RED; 276 277 // solve the reduced problem 278 279 SPX_MSG_INFO2(spxout, 280 spxout << std::endl 281 << "=========================" << std::endl 282 << "Solving: Reduced Problem." << std::endl 283 << "=========================" << std::endl 284 << std::endl); 285 286 _hasBasis = hasRedBasis; 287 // solving the reduced problem 288 _decompSimplifyAndSolve(_solver, _slufactor, !algIterCount, !algIterCount); 289 290 stop = decompTerminate(realParam( 291 SoPlexBase<R>::TIMELIMIT)); // checking whether the algorithm should terminate 292 293 // updating the algorithm iterations statistics 294 _statistics->callsReducedProb++; 295 _statistics->redProbStatus = _solver.status(); 296 297 assert(_isRealLPLoaded); 298 hasRedBasis = _hasBasis; 299 _hasBasis = false; 300 301 // checking the status of the reduced problem 302 // It is expected that infeasibility and unboundedness will be found in the initialisation solve. If this is 303 // not the case and the reduced problem is infeasible or unbounded the decomposition simplex will terminate. 304 // If the decomposition simplex terminates, then the original problem will be solved using the stored basis. 305 if(_solver.status() != SPxSolverBase<R>::OPTIMAL) 306 { 307 if(_solver.status() == SPxSolverBase<R>::UNBOUNDED) 308 SPX_MSG_INFO2(spxout, spxout << "Unbounded reduced problem." << std::endl); 309 310 if(_solver.status() == SPxSolverBase<R>::INFEASIBLE) 311 SPX_MSG_INFO2(spxout, spxout << "Infeasible reduced problem." << std::endl); 312 313 SPX_MSG_INFO2(spxout, spxout << "Reduced problem status: " << static_cast<int> 314 (_solver.status()) << std::endl); 315 316 redProbError = true; 317 break; 318 } 319 320 // as a final check, if no iterations were performed with the reduced problem, then the algorithm terminates 321 if(_statistics->iterations == previter) 322 { 323 SPX_MSG_WARNING(spxout, 324 spxout << "WIMDSM02: reduced problem performed zero iterations. Terminating." << std::endl;); 325 326 noRedprobIter = true; 327 stop = true; 328 break; 329 } 330 331 // printing display line 332 printDecompDisplayLine(_solver, orig_verbosity, !algIterCount, !algIterCount); 333 334 // get the dual solutions from the reduced problem 335 VectorBase<R> reducedLPDualVector(_solver.nRows()); 336 VectorBase<R> reducedLPRedcostVector(_solver.nCols()); 337 _solver.getDualSol(reducedLPDualVector); 338 _solver.getRedCostSol(reducedLPRedcostVector); 339 340 341 // Using the solution to the reduced problem: 342 // In the first iteration of the algorithm create the complementary problem. 343 // In the subsequent iterations of the algorithm update the complementary problem. 344 if(algIterCount == 0) 345 _formDecompComplementaryProblem(); 346 else 347 { 348 if(boolParam(SoPlexBase<R>::USECOMPDUAL)) 349 _updateDecompComplementaryDualProblem(false); 350 else 351 _updateDecompComplementaryPrimalProblem(false); 352 } 353 354 // setting the current solving mode. 355 _currentProb = DECOMP_COMP; 356 357 // solve the complementary problem 358 SPX_MSG_INFO2(spxout, 359 spxout << std::endl 360 << "=========================" << std::endl 361 << "Solving: Complementary Problem." << std::endl 362 << "=========================" << std::endl 363 << std::endl); 364 365 if(!explicitviol) 366 { 367 //_compSolver.writeFileLPBase("comp.lp"); 368 369 _decompSimplifyAndSolve(_compSolver, _compSlufactor, true, true); 370 371 SPX_MSG_INFO2(spxout, spxout << "Iteration " << algIterCount 372 << "Objective Value: " << std::setprecision(10) << _compSolver.objValue() 373 << std::endl); 374 } 375 376 377 assert(_isRealLPLoaded); 378 379 380 // Check whether the complementary problem is solved with a non-negative objective function, is infeasible or 381 // unbounded. If true, then stop the algorithm. 382 if(!explicitviol && (GE(_compSolver.objValue(), R(0.0), R(1e-20)) 383 || _compSolver.status() == SPxSolverBase<R>::INFEASIBLE 384 || _compSolver.status() == SPxSolverBase<R>::UNBOUNDED)) 385 { 386 _statistics->compProbStatus = _compSolver.status(); 387 _statistics->finalCompObj = _compSolver.objValue(); 388 389 if(_compSolver.status() == SPxSolverBase<R>::UNBOUNDED) 390 SPX_MSG_INFO2(spxout, spxout << "Unbounded complementary problem." << std::endl); 391 392 if(_compSolver.status() == SPxSolverBase<R>::INFEASIBLE) 393 SPX_MSG_INFO2(spxout, spxout << "Infeasible complementary problem." << std::endl); 394 395 if(_compSolver.status() == SPxSolverBase<R>::INFEASIBLE 396 || _compSolver.status() == SPxSolverBase<R>::UNBOUNDED) 397 explicitviol = true; 398 399 stop = true; 400 } 401 402 if(!stop && !explicitviol) 403 { 404 // get the primal solutions from the complementary problem 405 VectorBase<R> compLPPrimalVector(_compSolver.nCols()); 406 _compSolver.getPrimalSol(compLPPrimalVector); 407 408 // get the dual solutions from the complementary problem 409 VectorBase<R> compLPDualVector(_compSolver.nRows()); 410 _compSolver.getDualSol(compLPDualVector); 411 412 // updating the reduced problem 413 _updateDecompReducedProblem(_compSolver.objValue(), reducedLPDualVector, reducedLPRedcostVector, 414 compLPPrimalVector, compLPDualVector); 415 } 416 // if the complementary problem is infeasible or unbounded, it is possible that the algorithm can continue. 417 // a check of the original problem is required to determine whether there are any violations. 418 else if(_compSolver.status() == SPxSolverBase<R>::INFEASIBLE 419 || _compSolver.status() == SPxSolverBase<R>::UNBOUNDED 420 || explicitviol) 421 { 422 // getting the primal vector from the reduced problem 423 VectorBase<R> reducedLPPrimalVector(_solver.nCols()); 424 _solver.getPrimalSol(reducedLPPrimalVector); 425 426 // checking the optimality of the reduced problem solution with the original problem 427 _checkOriginalProblemOptimality(reducedLPPrimalVector, true); 428 429 // if there are any violated rows or bounds then stop is reset and the algorithm continues. 430 if(_nDecompViolBounds > 0 || _nDecompViolRows > 0) 431 stop = false; 432 433 if(_nDecompViolBounds == 0 && _nDecompViolRows == 0) 434 stop = true; 435 436 // updating the reduced problem with the original problem violated rows 437 if(!stop) 438 _updateDecompReducedProblemViol(false); 439 } 440 441 442 // ============================================================================= 443 // Code check completed up to here 444 // ============================================================================= 445 446 numDecompIter++; 447 algIterCount++; 448 } 449 450 // if there is an error in solving the reduced problem, i.e. infeasible or unbounded, then we leave the 451 // decomposition solve and resolve the original. Infeasibility should be dealt with in the original problem. 452 if(!redProbError) 453 { 454 #ifndef NDEBUG 455 // computing the solution for the original variables 456 VectorBase<R> reducedLPPrimalVector(_solver.nCols()); 457 _solver.getPrimalSol(reducedLPPrimalVector); 458 459 // checking the optimality of the reduced problem solution with the original problem 460 _checkOriginalProblemOptimality(reducedLPPrimalVector, true); 461 462 // resetting the verbosity level 463 spxout.setVerbosity(orig_verbosity); 464 #endif 465 466 // This is an additional check to ensure that the complementary problem is solving correctly. 467 // This is only used for debugging 468 #ifdef PERFORM_COMPPROB_CHECK 469 470 // solving the complementary problem with the original objective function 471 if(boolParam(SoPlexBase<R>::USECOMPDUAL)) 472 _setComplementaryDualOriginalObjective(); 473 else 474 _setComplementaryPrimalOriginalObjective(); 475 476 // for clean up 477 // the real lp has to be set to not loaded. 478 //_isRealLPLoaded = false; 479 // the basis has to be set to false 480 _hasBasis = false; 481 482 if(boolParam(SoPlexBase<R>::USECOMPDUAL)) 483 { 484 SPxLPBase<R> compDualLP; 485 _compSolver.buildDualProblem(compDualLP, _decompPrimalRowIDs.get_ptr(), 486 _decompPrimalColIDs.get_ptr(), 487 _decompDualRowIDs.get_ptr(), _decompDualColIDs.get_ptr(), &_nPrimalRows, &_nPrimalCols, &_nDualRows, 488 &_nDualCols); 489 490 _compSolver.loadLP(compDualLP); 491 } 492 493 _decompSimplifyAndSolve(_compSolver, _compSlufactor, true, true); 494 495 _solReal._hasPrimal = true; 496 _hasSolReal = true; 497 // get the primal solutions from the reduced problem 498 VectorBase<R> testPrimalVector(_compSolver.nCols()); 499 _compSolver.getPrimalSol(testPrimalVector); 500 _solReal._primal.reDim(_compSolver.nCols()); 501 _solReal._primal = testPrimalVector; 502 503 R maxviol = 0; 504 R sumviol = 0; 505 506 // Since the original objective value has now been installed in the complementary problem, the solution will be 507 // a solution to the original problem. 508 // checking the bound violation of the solution from complementary problem 509 if(getDecompBoundViolation(maxviol, sumviol)) 510 SPX_MSG_INFO1(spxout, spxout << "Bound violation - " 511 << "Max: " << std::setprecision(20) << maxviol << " " 512 << "Sum: " << sumviol << std::endl); 513 514 // checking the row violation of the solution from the complementary problem 515 if(getDecompRowViolation(maxviol, sumviol)) 516 SPX_MSG_INFO1(spxout, spxout << "Row violation - " 517 << "Max: " << std::setprecision(21) << maxviol << " " 518 << "Sum: " << sumviol << std::endl); 519 520 SPX_MSG_INFO1(spxout, spxout << "Objective Value: " << _compSolver.objValue() << std::endl); 521 #endif 522 } 523 524 // resetting the verbosity level 525 spxout.setVerbosity(orig_verbosity); 526 527 SPX_MSG_INFO1(spxout, 528 spxout << "======== Decomposition solve completed ========" << std::endl 529 << std::endl 530 << "======== Resolving original problem ========" << std::endl 531 ); 532 533 // if there is a reduced problem error in the first iteration the complementary problme has not been 534 // set up. In this case no memory has been allocated for _decompCompProbColIDsIdx and _fixedOrigVars. 535 if((!redProbError && !noRedprobIter) || algIterCount > 0) 536 { 537 spx_free(_decompCompProbColIDsIdx); 538 spx_free(_fixedOrigVars); 539 } 540 541 spx_free(_decompViolatedRows); 542 spx_free(_decompViolatedBounds); 543 spx_free(_decompReducedProbCols); 544 spx_free(_decompReducedProbRows); 545 546 // retreiving the original problem statistics prior to destroying it. 547 getOriginalProblemStatistics(); 548 549 // setting the reduced problem statistics 550 _statistics->numRedProbRows = numIncludedRows; 551 _statistics->numRedProbCols = _solver.nCols(); 552 553 554 // printing display line for resolve of problem 555 _solver.printDisplayLine(false, true); 556 557 // if there is an error solving the reduced problem the LP must be solved from scratch 558 if(redProbError) 559 { 560 SPX_MSG_INFO1(spxout, spxout << "=== Reduced problem error - Solving original ===" << std::endl); 561 562 // the solver is loaded with the realLP to solve the problem from scratch 563 _solver.loadLP(*_realLP); 564 spx_free(_realLP); 565 _realLP = &_solver; 566 _isRealLPLoaded = true; 567 568 // returning the sense to minimise 569 if(intParam(SoPlexBase<R>::OBJSENSE) == SoPlexBase<R>::OBJSENSE_MINIMIZE) 570 { 571 assert(_solver.spxSense() == SPxLPBase<R>::MAXIMIZE); 572 573 _solver.changeObj(-(_solver.maxObj())); 574 _solver.changeSense(SPxLPBase<R>::MINIMIZE); 575 576 // Need to add commands to multiply the objective solution values by -1 577 } 578 579 //_solver.setBasis(basisStatusRows.get_const_ptr(), basisStatusCols.get_const_ptr()); 580 _preprocessAndSolveReal(true); 581 } 582 else 583 { 584 // resolving the problem to update the real lp and solve with the correct objective. 585 // the realLP is updated with the current solver. This is to resolve the problem to get the correct solution 586 // values 587 _realLP->~SPxLPBase<R>(); 588 spx_free(_realLP); 589 _realLP = &_solver; 590 _isRealLPLoaded = true; 591 592 // returning the sense to minimise 593 if(intParam(SoPlexBase<R>::OBJSENSE) == SoPlexBase<R>::OBJSENSE_MINIMIZE) 594 { 595 assert(_solver.spxSense() == SPxLPBase<R>::MAXIMIZE); 596 597 _solver.changeObj(-(_solver.maxObj())); 598 _solver.changeSense(SPxLPBase<R>::MINIMIZE); 599 600 // Need to add commands to multiply the objective solution values by -1 601 } 602 603 _preprocessAndSolveReal(false); 604 } 605 606 // stop timing 607 _statistics->solvingTime->stop(); 608 } 609 610 611 612 /// creating copies of the original problem that will be manipulated to form the reduced and complementary problems 613 template <class R> 614 void SoPlexBase<R>::_createDecompReducedAndComplementaryProblems() 615 { 616 // the reduced problem is formed from the current problem 617 // So, we copy the _solver to the _realLP and work on the _solver 618 // NOTE: there is no need to preprocess because we always have a starting basis. 619 _realLP = nullptr; 620 spx_alloc(_realLP); 621 _realLP = new(_realLP) SPxLPBase<R>(_solver); 622 623 // allocating memory for the reduced problem rows and cols flag array 624 _decompReducedProbRows = nullptr; 625 spx_alloc(_decompReducedProbRows, numRows()); 626 _decompReducedProbCols = nullptr; 627 spx_alloc(_decompReducedProbCols, numCols()); 628 629 // the complementary problem is formulated with all incompatible rows and those from the reduced problem that have 630 // a positive reduced cost. 631 _compSolver = _solver; 632 _compSolver.setOutstream(spxout); 633 _compSolver.setBasisSolver(&_compSlufactor); 634 635 // allocating memory for the violated bounds and rows arrays 636 _decompViolatedBounds = nullptr; 637 _decompViolatedRows = nullptr; 638 spx_alloc(_decompViolatedBounds, numCols()); 639 spx_alloc(_decompViolatedRows, numRows()); 640 _nDecompViolBounds = 0; 641 _nDecompViolRows = 0; 642 } 643 644 645 646 /// forms the reduced problem 647 template <class R> 648 void SoPlexBase<R>::_formDecompReducedProblem(bool& stop) 649 { 650 SPX_MSG_INFO2(spxout, spxout << "Forming the Reduced problem" << std::endl); 651 int* nonposind = 0; 652 int* compatind = 0; 653 int* rowsforremoval = 0; 654 int* colsforremoval = 0; 655 int nnonposind = 0; 656 int ncompatind = 0; 657 658 assert(_solver.nCols() == numCols()); 659 assert(_solver.nRows() == numRows()); 660 661 // capturing the basis used for the transformation 662 _decompTransBasis = _solver.basis(); 663 664 // setting row counter to zero 665 numIncludedRows = 0; 666 667 // the _decompLP is used as a helper LP object. This is used so that _realLP can still hold the original problem. 668 _decompLP = 0; 669 spx_alloc(_decompLP); 670 _decompLP = new(_decompLP) SPxLPBase<R>(_solver); 671 672 // retreiving the basis information 673 _basisStatusRows.reSize(numRows()); 674 _basisStatusCols.reSize(numCols()); 675 _solver.getBasis(_basisStatusRows.get_ptr(), _basisStatusCols.get_ptr()); 676 677 // get the indices of the rows with positive dual multipliers and columns with positive reduced costs. 678 spx_alloc(nonposind, numCols()); 679 spx_alloc(colsforremoval, numCols()); 680 681 if(!stop) 682 _getZeroDualMultiplierIndices(_decompFeasVector, nonposind, colsforremoval, &nnonposind, stop); 683 684 // get the compatible columns from the constraint matrix w.r.t the current basis matrix 685 SPX_MSG_INFO2(spxout, spxout << "Computing the compatible columns" << std::endl 686 << "Solving time: " << solveTime() << std::endl); 687 688 spx_alloc(compatind, _solver.nRows()); 689 spx_alloc(rowsforremoval, _solver.nRows()); 690 691 if(!stop) 692 _getCompatibleColumns(_decompFeasVector, nonposind, compatind, rowsforremoval, colsforremoval, 693 nnonposind, 694 &ncompatind, true, stop); 695 696 int* compatboundcons = 0; 697 int ncompatboundcons = 0; 698 spx_alloc(compatboundcons, numCols()); 699 700 LPRowSetBase<R> boundcons; 701 702 // identifying the compatible bound constraints 703 if(!stop) 704 _getCompatibleBoundCons(boundcons, compatboundcons, nonposind, &ncompatboundcons, nnonposind, stop); 705 706 // delete rows and columns from the LP to form the reduced problem 707 SPX_MSG_INFO2(spxout, spxout << "Deleting rows and columns to form the reduced problem" << std::endl 708 << "Solving time: " << solveTime() << std::endl); 709 710 // allocating memory to add bound constraints 711 SPxRowId* addedrowids = 0; 712 spx_alloc(addedrowids, ncompatboundcons); 713 714 // computing the reduced problem obj coefficient vector and updating the problem 715 if(!stop) 716 { 717 _computeReducedProbObjCoeff(stop); 718 719 _solver.loadLP(*_decompLP); 720 721 // the colsforremoval are the columns with a zero reduced cost. 722 // the rowsforremoval are the rows identified as incompatible. 723 _solver.removeRows(rowsforremoval); 724 // adding the rows for the compatible bound constraints 725 _solver.addRows(addedrowids, boundcons); 726 727 for(int i = 0; i < ncompatboundcons; i++) 728 _decompReducedProbColRowIDs[compatboundcons[i]] = addedrowids[i]; 729 } 730 731 732 // freeing allocated memory 733 spx_free(addedrowids); 734 spx_free(compatboundcons); 735 spx_free(rowsforremoval); 736 spx_free(compatind); 737 spx_free(colsforremoval); 738 spx_free(nonposind); 739 740 _decompLP->~SPxLPBase<R>(); 741 spx_free(_decompLP); 742 } 743 744 745 746 /// forms the complementary problem 747 template <class R> 748 void SoPlexBase<R>::_formDecompComplementaryProblem() 749 { 750 // delete rows and columns from the LP to form the reduced problem 751 _nElimPrimalRows = 0; 752 _decompElimPrimalRowIDs.reSize( 753 _realLP->nRows()); // the number of eliminated rows is less than the number of rows 754 // in the reduced problem 755 _decompCompProbColIDsIdx = 0; 756 spx_alloc(_decompCompProbColIDsIdx, _realLP->nCols()); 757 758 _fixedOrigVars = 0; 759 spx_alloc(_fixedOrigVars, _realLP->nCols()); 760 761 for(int i = 0; i < _realLP->nCols(); i++) 762 { 763 _decompCompProbColIDsIdx[i] = -1; 764 _fixedOrigVars[i] = -2; // this must be set to a value differet from -1, 0, 1 to ensure the 765 // _updateComplementaryFixedPrimalVars function updates all variables in the problem. 766 } 767 768 int naddedrows = 0; 769 DataArray < SPxRowId > rangedRowIds(numRows()); 770 _deleteAndUpdateRowsComplementaryProblem(rangedRowIds.get_ptr(), naddedrows); 771 772 if(boolParam( 773 SoPlexBase<R>::USECOMPDUAL)) // if we use the dual formulation of the complementary problem, we must 774 // perform the transformation 775 { 776 // initialising the arrays to store the row id's from the primal and the col id's from the dual 777 _decompPrimalRowIDs.reSize(3 * _compSolver.nRows()); 778 _decompPrimalColIDs.reSize(3 * _compSolver.nCols()); 779 _decompDualRowIDs.reSize(3 * _compSolver.nCols()); 780 _decompDualColIDs.reSize(3 * _compSolver.nRows()); 781 _nPrimalRows = 0; 782 _nPrimalCols = 0; 783 _nDualRows = 0; 784 _nDualCols = 0; 785 786 // convert complementary problem to dual problem 787 SPxLPBase<R> compDualLP; 788 _compSolver.buildDualProblem(compDualLP, _decompPrimalRowIDs.get_ptr(), 789 _decompPrimalColIDs.get_ptr(), 790 _decompDualRowIDs.get_ptr(), _decompDualColIDs.get_ptr(), &_nPrimalRows, &_nPrimalCols, &_nDualRows, 791 &_nDualCols); 792 compDualLP.setOutstream(spxout); 793 794 // setting the id index array for the complementary problem 795 for(int i = 0; i < _nPrimalRows; i++) 796 { 797 // a primal row may result in two dual columns. So the index array points to the first of the indicies for the 798 // primal row. 799 if(i + 1 < _nPrimalRows && 800 _realLP->number(SPxRowId(_decompPrimalRowIDs[i])) == _realLP->number(SPxRowId( 801 _decompPrimalRowIDs[i + 1]))) 802 i++; 803 } 804 805 for(int i = 0; i < _nPrimalCols; i++) 806 { 807 if(_decompPrimalColIDs[i].getIdx() != _compSlackColId.getIdx()) 808 _decompCompProbColIDsIdx[_realLP->number(_decompPrimalColIDs[i])] = i; 809 } 810 811 // retrieving the dual row id for the complementary slack column 812 // there should be a one to one relationship between the number of primal columns and the number of dual rows. 813 // hence, it should be possible to equate the dual row id to the related primal column. 814 assert(_nPrimalCols == _nDualRows); 815 assert(_compSolver.nCols() == compDualLP.nRows()); 816 _compSlackDualRowId = compDualLP.rId(_compSolver.number(_compSlackColId)); 817 818 _compSolver.loadLP(compDualLP); 819 820 _compSolver.init(); 821 822 _updateComplementaryDualSlackColCoeff(); 823 824 // initalising the array for the dual columns representing the fixed variables in the complementary problem. 825 _decompFixedVarDualIDs.reSize(_realLP->nCols()); 826 _decompVarBoundDualIDs.reSize(2 * _realLP->nCols()); 827 828 829 // updating the constraints for the complementary problem 830 _updateDecompComplementaryDualProblem(false); 831 } 832 else // when using the primal of the complementary problem, no transformation of the problem is required. 833 { 834 // initialising the arrays to store the row and column id's from the original problem the col id's from the dual 835 // if a row or column is not included in the complementary problem, the ID is set to INVALID in the 836 // _decompPrimalRowIDs or _decompPrimalColIDs respectively. 837 _decompPrimalRowIDs.reSize(_compSolver.nRows() * 2); 838 _decompPrimalColIDs.reSize(_compSolver.nCols()); 839 _decompCompPrimalRowIDs.reSize(_compSolver.nRows() * 2); 840 _decompCompPrimalColIDs.reSize(_compSolver.nCols()); 841 _nPrimalRows = 0; 842 _nPrimalCols = 0; 843 844 // at this point in time the complementary problem contains all rows and columns from the original problem 845 int addedrangedrows = 0; 846 _nCompPrimalRows = 0; 847 _nCompPrimalCols = 0; 848 849 for(int i = 0; i < _realLP->nRows(); i++) 850 { 851 assert(_nCompPrimalRows < _compSolver.nRows()); 852 _decompPrimalRowIDs[_nPrimalRows] = _realLP->rId(i); 853 _decompCompPrimalRowIDs[_nCompPrimalRows] = _compSolver.rId(i); 854 _nPrimalRows++; 855 _nCompPrimalRows++; 856 857 if(_realLP->rowType(i) == LPRowBase<R>::RANGE || _realLP->rowType(i) == LPRowBase<R>::EQUAL) 858 { 859 assert(EQ(_compSolver.lhs(rangedRowIds[addedrangedrows]), _realLP->lhs(i), _compSolver.epsilon())); 860 assert(LE(_compSolver.rhs(rangedRowIds[addedrangedrows]), R(infinity), _compSolver.epsilon())); 861 assert(LE(_compSolver.lhs(i), R(-infinity), _compSolver.epsilon())); 862 assert(LT(_compSolver.rhs(i), R(infinity), _compSolver.epsilon())); 863 864 _decompPrimalRowIDs[_nPrimalRows] = _realLP->rId(i); 865 _decompCompPrimalRowIDs[_nCompPrimalRows] = rangedRowIds[addedrangedrows]; 866 _nPrimalRows++; 867 _nCompPrimalRows++; 868 addedrangedrows++; 869 } 870 } 871 872 for(int i = 0; i < _compSolver.nCols(); i++) 873 { 874 875 if(_compSolver.cId(i).getIdx() != _compSlackColId.getIdx()) 876 { 877 _decompPrimalColIDs[i] = _realLP->cId(i); 878 _decompCompPrimalColIDs[i] = _compSolver.cId(i); 879 _nPrimalCols++; 880 _nCompPrimalCols++; 881 882 _decompCompProbColIDsIdx[_realLP->number(_decompPrimalColIDs[i])] = i; 883 } 884 } 885 886 // updating the constraints for the complementary problem 887 _updateDecompComplementaryPrimalProblem(false); 888 } 889 } 890 891 892 893 /// simplifies the problem and solves 894 template <class R> 895 void SoPlexBase<R>::_decompSimplifyAndSolve(SPxSolverBase<R>& solver, SLUFactor<R>& sluFactor, 896 bool fromScratch, bool applyPreprocessing) 897 { 898 if(realParam(SoPlexBase<R>::TIMELIMIT) < realParam(SoPlexBase<R>::INFTY)) 899 solver.setTerminationTime(Real(realParam(SoPlexBase<R>::TIMELIMIT)) - 900 _statistics->solvingTime->time()); 901 902 solver.changeObjOffset(realParam(SoPlexBase<R>::OBJ_OFFSET)); 903 _statistics->preprocessingTime->start(); 904 905 typename SPxSimplifier<R>::Result result = SPxSimplifier<R>::OKAY; 906 907 /* delete starting basis if solving from scratch */ 908 if(fromScratch) 909 { 910 try 911 { 912 solver.reLoad(); 913 } 914 catch(const SPxException& E) 915 { 916 SPX_MSG_ERROR(spxout << "Caught exception <" << E.what() << "> during simplify and solve.\n"); 917 } 918 } 919 920 //assert(!fromScratch || solver.status() == SPxSolverBase<R>::NO_PROBLEM); 921 922 if(applyPreprocessing) 923 { 924 _enableSimplifierAndScaler(); 925 solver.setTerminationValue(realParam(SoPlexBase<R>::INFTY)); 926 } 927 else 928 { 929 _disableSimplifierAndScaler(); 930 ///@todo implement for both objective senses 931 solver.setTerminationValue(solver.spxSense() == SPxLPBase<R>::MINIMIZE 932 ? realParam(SoPlexBase<R>::OBJLIMIT_UPPER) : realParam(SoPlexBase<R>::INFTY)); 933 } 934 935 /* store original lp */ 936 applyPreprocessing = (_scaler != nullptr || _simplifier != NULL); 937 938 // @TODO The use of _isRealLPLoaded is not correct. An additional parameter would be useful for this algorithm. 939 // Maybe a parameter _isDecompLPLoaded? 940 if(_isRealLPLoaded) 941 { 942 //assert(_decompLP == &solver); // there should be an assert here, but I don't know what to use. 943 944 // preprocessing is always applied to the LP in the solver; hence we have to create a copy of the original LP 945 // if preprocessing is turned on 946 if(applyPreprocessing) 947 { 948 _decompLP = 0; 949 spx_alloc(_decompLP); 950 _decompLP = new(_decompLP) SPxLPBase<R>(solver); 951 _isRealLPLoaded = false; 952 } 953 else 954 _decompLP = &solver; 955 } 956 else 957 { 958 assert(_decompLP != &solver); 959 960 // ensure that the solver has the original problem 961 solver.loadLP(*_decompLP); 962 963 // load basis if available 964 if(_hasBasis) 965 { 966 assert(_basisStatusRows.size() == solver.nRows()); 967 assert(_basisStatusCols.size() == solver.nCols()); 968 969 ///@todo this should not fail even if the basis is invalid (wrong dimension or wrong number of basic 970 /// entries); fix either in SPxSolverBase or in SPxBasisBase 971 solver.setBasis(_basisStatusRows.get_const_ptr(), _basisStatusCols.get_const_ptr()); 972 } 973 974 // if there is no preprocessing, then the original and the transformed problem are identical and it is more 975 // memory-efficient to keep only the problem in the solver 976 if(!applyPreprocessing) 977 { 978 _decompLP->~SPxLPBase<R>(); 979 spx_free(_decompLP); 980 _decompLP = &solver; 981 _isRealLPLoaded = true; 982 } 983 else 984 { 985 _decompLP = 0; 986 spx_alloc(_decompLP); 987 _decompLP = new(_decompLP) SPxLPBase<R>(solver); 988 } 989 } 990 991 // assert that we have two problems if and only if we apply preprocessing 992 assert(_decompLP == &solver || applyPreprocessing); 993 assert(_decompLP != &solver || !applyPreprocessing); 994 995 // apply problem simplification 996 if(_simplifier != 0) 997 { 998 Real remainingTime = _solver.getMaxTime() - _solver.time(); 999 result = _simplifier->simplify(solver, remainingTime); 1000 solver.changeObjOffset(_simplifier->getObjoffset() + realParam(SoPlexBase<R>::OBJ_OFFSET)); 1001 } 1002 1003 _statistics->preprocessingTime->stop(); 1004 1005 // run the simplex method if problem has not been solved by the simplifier 1006 if(result == SPxSimplifier<R>::OKAY) 1007 { 1008 if(_scaler != nullptr) 1009 _scaler->scale(solver); 1010 1011 bool _hadBasis = _hasBasis; 1012 1013 _statistics->simplexTime->start(); 1014 1015 try 1016 { 1017 solver.solve(); 1018 } 1019 catch(const SPxException& E) 1020 { 1021 SPX_MSG_ERROR(std::cerr << "Caught exception <" << E.what() << "> while solving real LP.\n"); 1022 _status = SPxSolverBase<R>::ERROR; 1023 } 1024 catch(...) 1025 { 1026 SPX_MSG_ERROR(std::cerr << "Caught unknown exception while solving real LP.\n"); 1027 _status = SPxSolverBase<R>::ERROR; 1028 } 1029 1030 _statistics->simplexTime->stop(); 1031 1032 // record statistics 1033 // only record the main statistics for the original problem and reduced problem. 1034 // the complementary problem is a pivoting rule, so these will be held in other statistics. 1035 // statistics on the number of iterations for each problem is stored in individual variables. 1036 if(_currentProb == DECOMP_ORIG || _currentProb == DECOMP_RED) 1037 { 1038 _statistics->iterations += solver.iterations(); 1039 _statistics->iterationsPrimal += solver.primalIterations(); 1040 _statistics->iterationsFromBasis += _hadBasis ? solver.iterations() : 0; 1041 _statistics->boundflips += solver.boundFlips(); 1042 _statistics->luFactorizationTimeReal += sluFactor.getFactorTime(); 1043 _statistics->luSolveTimeReal += sluFactor.getSolveTime(); 1044 _statistics->luFactorizationsReal += sluFactor.getFactorCount(); 1045 _statistics->luSolvesReal += sluFactor.getSolveCount(); 1046 sluFactor.resetCounters(); 1047 1048 _statistics->degenPivotsPrimal += solver.primalDegeneratePivots(); 1049 _statistics->degenPivotsDual += solver.dualDegeneratePivots(); 1050 1051 _statistics->sumDualDegen += _solver.sumDualDegeneracy(); 1052 _statistics->sumPrimalDegen += _solver.sumPrimalDegeneracy(); 1053 1054 if(_currentProb == DECOMP_ORIG) 1055 _statistics->iterationsInit += solver.iterations(); 1056 else 1057 _statistics->iterationsRedProb += solver.iterations(); 1058 } 1059 1060 if(_currentProb == DECOMP_COMP) 1061 _statistics->iterationsCompProb += solver.iterations(); 1062 1063 } 1064 1065 // check the result and run again without preprocessing if necessary 1066 _evaluateSolutionDecomp(solver, sluFactor, result); 1067 } 1068 1069 1070 1071 /// loads original problem into solver and solves again after it has been solved to optimality with preprocessing 1072 template <class R> 1073 void SoPlexBase<R>::_decompResolveWithoutPreprocessing(SPxSolverBase<R>& solver, 1074 SLUFactor<R>& sluFactor, typename SPxSimplifier<R>::Result result) 1075 { 1076 assert(_simplifier != 0 || _scaler != nullptr); 1077 assert(result == SPxSimplifier<R>::VANISHED 1078 || (result == SPxSimplifier<R>::OKAY 1079 && (solver.status() == SPxSolverBase<R>::OPTIMAL 1080 || solver.status() == SPxSolverBase<R>::ABORT_DECOMP 1081 || solver.status() == SPxSolverBase<R>::ABORT_EXDECOMP))); 1082 1083 // if simplifier is active and LP is solved in presolving or to optimality, then we unsimplify to get the basis 1084 if(_simplifier != 0) 1085 { 1086 assert(!_simplifier->isUnsimplified()); 1087 1088 bool vanished = result == SPxSimplifier<R>::VANISHED; 1089 1090 // get solution vectors for transformed problem 1091 VectorBase<R> primal(vanished ? 0 : solver.nCols()); 1092 VectorBase<R> slacks(vanished ? 0 : solver.nRows()); 1093 VectorBase<R> dual(vanished ? 0 : solver.nRows()); 1094 VectorBase<R> redCost(vanished ? 0 : solver.nCols()); 1095 1096 assert(!_isRealLPLoaded); 1097 _basisStatusRows.reSize(_decompLP->nRows()); 1098 _basisStatusCols.reSize(_decompLP->nCols()); 1099 assert(vanished || _basisStatusRows.size() >= solver.nRows()); 1100 assert(vanished || _basisStatusCols.size() >= solver.nCols()); 1101 1102 if(!vanished) 1103 { 1104 assert(solver.status() == SPxSolverBase<R>::OPTIMAL 1105 || solver.status() == SPxSolverBase<R>::ABORT_DECOMP 1106 || solver.status() == SPxSolverBase<R>::ABORT_EXDECOMP); 1107 1108 solver.getPrimalSol(primal); 1109 solver.getSlacks(slacks); 1110 solver.getDualSol(dual); 1111 solver.getRedCostSol(redCost); 1112 1113 // unscale vectors 1114 if(_scaler && solver.isScaled()) 1115 { 1116 _scaler->unscalePrimal(solver, primal); 1117 _scaler->unscaleSlacks(solver, slacks); 1118 _scaler->unscaleDual(solver, dual); 1119 _scaler->unscaleRedCost(solver, redCost); 1120 } 1121 1122 // get basis of transformed problem 1123 solver.getBasis(_basisStatusRows.get_ptr(), _basisStatusCols.get_ptr()); 1124 } 1125 1126 try 1127 { 1128 _simplifier->unsimplify(primal, dual, slacks, redCost, _basisStatusRows.get_ptr(), 1129 _basisStatusCols.get_ptr(), solver.status() == SPxSolverBase<R>::OPTIMAL); 1130 _simplifier->getBasis(_basisStatusRows.get_ptr(), _basisStatusCols.get_ptr()); 1131 _hasBasis = true; 1132 } 1133 catch(const SPxException& E) 1134 { 1135 SPX_MSG_ERROR(spxout << "Caught exception <" << E.what() << 1136 "> during unsimplification. Resolving without simplifier and scaler.\n"); 1137 } 1138 catch(...) 1139 { 1140 SPX_MSG_ERROR(spxout << 1141 "Caught unknown exception during unsimplification. Resolving without simplifier and scaler.\n"); 1142 _status = SPxSolverBase<R>::ERROR; 1143 } 1144 } 1145 // if the original problem is not in the solver because of scaling, we also need to store the basis 1146 else if(_scaler != nullptr) 1147 { 1148 _basisStatusRows.reSize(numRows()); 1149 _basisStatusCols.reSize(numCols()); 1150 assert(_basisStatusRows.size() == solver.nRows()); 1151 assert(_basisStatusCols.size() == solver.nCols()); 1152 1153 solver.getBasis(_basisStatusRows.get_ptr(), _basisStatusCols.get_ptr()); 1154 _hasBasis = true; 1155 } 1156 1157 // resolve the original problem 1158 _decompSimplifyAndSolve(solver, sluFactor, false, false); 1159 return; 1160 } 1161 1162 1163 /// updates the reduced problem with additional rows using the solution to the complementary problem 1164 template <class R> 1165 void SoPlexBase<R>::_updateDecompReducedProblem(R objValue, VectorBase<R> dualVector, 1166 VectorBase<R> redcostVector, 1167 VectorBase<R> compPrimalVector, VectorBase<R> compDualVector) 1168 { 1169 R feastol = realParam(SoPlexBase<R>::FEASTOL); 1170 1171 R maxDualRatio = R(infinity); 1172 1173 bool usecompdual = boolParam(SoPlexBase<R>::USECOMPDUAL); 1174 1175 for(int i = 0; i < _nPrimalRows; i++) 1176 { 1177 R reducedProbDual = 0; 1178 R compProbPrimal = 0; 1179 R dualRatio = 0; 1180 int rownumber = _realLP->number(SPxRowId(_decompPrimalRowIDs[i])); 1181 1182 if(_decompReducedProbRows[rownumber] && _solver.isBasic(_decompReducedProbRowIDs[rownumber])) 1183 { 1184 int solverRowNum = _solver.number(_decompReducedProbRowIDs[rownumber]); 1185 1186 // retreiving the reduced problem dual solutions and the complementary problem primal solutions 1187 reducedProbDual = dualVector[solverRowNum]; // this is y 1188 1189 if(usecompdual) 1190 compProbPrimal = compPrimalVector[_compSolver.number(SPxColId(_decompDualColIDs[i]))]; // this is u 1191 else 1192 compProbPrimal = compDualVector[_compSolver.number(SPxRowId(_decompCompPrimalRowIDs[i]))]; 1193 1194 // the variable in the basis is degenerate. 1195 if(EQ(reducedProbDual, R(0.0), feastol)) 1196 { 1197 SPX_MSG_WARNING(spxout, 1198 spxout << "WIMDSM01: reduced problem dual value is very close to zero." << std::endl;); 1199 continue; 1200 } 1201 1202 // the translation of the complementary primal problem to the dual some rows resulted in two columns. 1203 if(usecompdual && i < _nPrimalRows - 1 && 1204 _realLP->number(SPxRowId(_decompPrimalRowIDs[i])) == _realLP->number(SPxRowId( 1205 _decompPrimalRowIDs[i + 1]))) 1206 { 1207 i++; 1208 compProbPrimal += compPrimalVector[_compSolver.number(SPxColId(_decompDualColIDs[i]))]; // this is u 1209 } 1210 1211 1212 1213 // updating the ratio 1214 SoPlexBase<R>::DualSign varSign = getExpectedDualVariableSign(solverRowNum); 1215 1216 if(varSign == SoPlexBase<R>::IS_FREE || (varSign == SoPlexBase<R>::IS_POS 1217 && LE(compProbPrimal, (R)0, feastol)) || 1218 (varSign == SoPlexBase<R>::IS_NEG && GE(compProbPrimal, (R)0, feastol))) 1219 { 1220 dualRatio = R(infinity); 1221 } 1222 else 1223 { 1224 dualRatio = reducedProbDual / compProbPrimal; 1225 } 1226 1227 if(LT(dualRatio, maxDualRatio, feastol)) 1228 maxDualRatio = dualRatio; 1229 } 1230 else 1231 { 1232 if(usecompdual) 1233 compProbPrimal = compPrimalVector[_compSolver.number(SPxColId(_decompDualColIDs[i]))]; // this is v 1234 else 1235 compProbPrimal = compDualVector[_compSolver.number(SPxRowId(_decompCompPrimalRowIDs[i]))]; 1236 } 1237 1238 } 1239 1240 // setting the maxDualRatio to a maximum of 1 1241 if(maxDualRatio > 1.0) 1242 maxDualRatio = 1.0; 1243 1244 VectorBase<R> compProbRedcost( 1245 _compSolver.nCols()); // the reduced costs of the complementary problem 1246 1247 // Retrieving the reduced costs for each original problem row. 1248 _compSolver.getRedCostSol(compProbRedcost); 1249 1250 LPRowSetBase<R> updaterows; 1251 1252 1253 // Identifying the violated rows 1254 Array<RowViolation> violatedrows; 1255 int nviolatedrows = 0; 1256 int* newrowidx = 0; 1257 int nnewrowidx = 0; 1258 spx_alloc(newrowidx, _nPrimalRows); 1259 1260 violatedrows.reSize(_nPrimalRows); 1261 1262 bool ratioTest = true; 1263 1264 //ratioTest = false; 1265 for(int i = 0; i < _nPrimalRows; i++) 1266 { 1267 LPRowBase<R> origlprow; 1268 DSVectorBase<R> rowtoaddVec(_realLP->nCols()); 1269 R compProbPrimal = 0; 1270 R compRowRedcost = 0; 1271 int rownumber = _realLP->number(SPxRowId(_decompPrimalRowIDs[i])); 1272 1273 if(!_decompReducedProbRows[_realLP->number(SPxRowId(_decompPrimalRowIDs[i]))]) 1274 { 1275 int compRowNumber; 1276 1277 if(usecompdual) 1278 { 1279 compRowNumber = _compSolver.number(_decompDualColIDs[i]); 1280 // retreiving the complementary problem primal solutions 1281 compProbPrimal = compPrimalVector[compRowNumber]; // this is v 1282 compRowRedcost = compProbRedcost[compRowNumber]; 1283 } 1284 else 1285 { 1286 compRowNumber = _compSolver.number(_decompCompPrimalRowIDs[i]); 1287 // retreiving the complementary problem primal solutions 1288 compProbPrimal = compDualVector[compRowNumber]; // this is v 1289 } 1290 1291 // the translation of the complementary primal problem to the dual some rows resulted in two columns. 1292 if(usecompdual && i < _nPrimalRows - 1 && 1293 _realLP->number(SPxRowId(_decompPrimalRowIDs[i])) == _realLP->number(SPxRowId( 1294 _decompPrimalRowIDs[i + 1]))) 1295 { 1296 i++; 1297 compRowNumber = _compSolver.number(SPxColId(_decompDualColIDs[i])); 1298 compProbPrimal += compPrimalVector[compRowNumber]; // this is v 1299 compRowRedcost += compProbRedcost[compRowNumber]; 1300 } 1301 1302 SoPlexBase<R>::DualSign varSign = getOrigProbDualVariableSign(rownumber); 1303 1304 // add row to the reduced problem if the computed dual is of the correct sign for a feasible dual solution 1305 if(ratioTest && ((varSign == SoPlexBase<R>::IS_FREE && !isZero(compProbPrimal, feastol)) || 1306 (varSign == SoPlexBase<R>::IS_POS && GT(compProbPrimal * maxDualRatio, (R) 0, feastol)) || 1307 (varSign == SoPlexBase<R>::IS_NEG && LT(compProbPrimal * maxDualRatio, (R) 0, feastol)))) 1308 { 1309 //this set of statements are required to add rows to the reduced problem. This will add a row for every 1310 //violated constraint. I only want to add a row for the most violated constraint. That is why the row 1311 //adding functionality is put outside the for loop. 1312 if(!_decompReducedProbRows[rownumber]) 1313 { 1314 numIncludedRows++; 1315 assert(numIncludedRows <= _realLP->nRows()); 1316 } 1317 1318 violatedrows[nviolatedrows].idx = rownumber; 1319 violatedrows[nviolatedrows].violation = spxAbs(compProbPrimal * maxDualRatio); 1320 nviolatedrows++; 1321 } 1322 } 1323 } 1324 1325 1326 // sorting the violated rows by the absolute violation. 1327 // Only a predefined number of rows will be added to the reduced problem 1328 RowViolationCompare compare; 1329 compare.entry = violatedrows.get_const_ptr(); 1330 1331 // if no rows are identified by the pricing rule, we add rows based upon the constraint violations 1332 if(!ratioTest || nviolatedrows == 0) 1333 { 1334 _findViolatedRows(objValue, violatedrows, nviolatedrows); 1335 } 1336 1337 int sorted = 0; 1338 int sortsize = SOPLEX_MIN(intParam(SoPlexBase<R>::DECOMP_MAXADDEDROWS), nviolatedrows); 1339 1340 // only sorting if the sort size is less than the number of violated rows. 1341 if(sortsize > 0 && sortsize < nviolatedrows) 1342 sorted = SPxQuicksortPart(violatedrows.get_ptr(), compare, sorted + 1, nviolatedrows, sortsize); 1343 1344 // adding the violated rows. 1345 for(int i = 0; i < sortsize; i++) 1346 { 1347 updaterows.add(_transformedRows.lhs(violatedrows[i].idx), 1348 _transformedRows.rowVector(violatedrows[i].idx), 1349 _transformedRows.rhs(violatedrows[i].idx)); 1350 1351 _decompReducedProbRows[violatedrows[i].idx] = true; 1352 newrowidx[nnewrowidx] = violatedrows[i].idx; 1353 nnewrowidx++; 1354 } 1355 1356 SPxRowId* addedrowids = 0; 1357 spx_alloc(addedrowids, nnewrowidx); 1358 _solver.addRows(addedrowids, updaterows); 1359 1360 for(int i = 0; i < nnewrowidx; i++) 1361 _decompReducedProbRowIDs[newrowidx[i]] = addedrowids[i]; 1362 1363 // freeing allocated memory 1364 spx_free(addedrowids); 1365 spx_free(newrowidx); 1366 } 1367 1368 1369 1370 /// update the reduced problem with additional columns and rows based upon the violated original bounds and rows 1371 // TODO: Allow for the case that no rows are added. This should terminate the algorithm. 1372 // TODO: Check to make sure that only rows added to the problem do not currently exist in the reduced problem. 1373 template <class R> 1374 void SoPlexBase<R>::_updateDecompReducedProblemViol(bool allrows) 1375 { 1376 #ifdef NO_TOL 1377 R feastol = 0.0; 1378 #else 1379 #ifdef SOPLEX_USE_FEASTOL 1380 R feastol = realParam(SoPlexBase<R>::FEASTOL); 1381 #else 1382 R feastol = realParam(SoPlexBase<R>::EPSILON_ZERO); 1383 #endif 1384 #endif 1385 LPRowSetBase<R> updaterows; 1386 1387 int* newrowidx = 0; 1388 int nnewrowidx = 0; 1389 spx_alloc(newrowidx, _nPrimalRows); 1390 1391 int rowNumber; 1392 int bestrow = -1; 1393 R bestrownorm = R(infinity); 1394 R percenttoadd = 1; 1395 1396 int nrowstoadd = SOPLEX_MIN(intParam(SoPlexBase<R>::DECOMP_MAXADDEDROWS), _nDecompViolRows); 1397 1398 if(allrows) 1399 nrowstoadd = _nDecompViolRows; // adding all violated rows 1400 1401 SSVectorBase<R> y(_solver.nCols(), _solver.tolerances()); 1402 y.unSetup(); 1403 1404 // identifying the rows not included in the reduced problem that are violated by the current solution. 1405 for(int i = 0; i < nrowstoadd; i++) 1406 { 1407 rowNumber = _decompViolatedRows[i]; 1408 1409 if(!allrows) 1410 { 1411 R norm = 0; 1412 1413 // the rhs of this calculation are the rows of the constraint matrix 1414 // so we are solving y B = A_{i,.} 1415 try 1416 { 1417 _solver.basis().solve(y, _solver.vector(rowNumber)); 1418 } 1419 catch(const SPxException& E) 1420 { 1421 SPX_MSG_ERROR(spxout << "Caught exception <" << E.what() << "> while computing compatability.\n"); 1422 } 1423 1424 1425 // comparing the constraints based upon the row norm 1426 if(y.isSetup()) 1427 { 1428 for(int j = 0; j < y.size(); j++) 1429 { 1430 if(isZero(_solver.fVec()[i], feastol)) 1431 norm += spxAbs(y.value(j)) * spxAbs(y.value(j)); 1432 } 1433 } 1434 else 1435 { 1436 for(int j = 0; j < numCols(); j++) 1437 { 1438 if(isZero(_solver.fVec()[i], feastol)) 1439 norm += spxAbs(y[j]) * spxAbs(y[j]); 1440 } 1441 } 1442 1443 1444 // the best row is based upon the row norm 1445 // the best row is added if no violated row is found 1446 norm = soplex::spxSqrt(norm); 1447 1448 if(LT(norm, bestrownorm, _solver.epsilon())) 1449 { 1450 bestrow = rowNumber; 1451 bestrownorm = norm; 1452 } 1453 1454 1455 // adding the violated row 1456 if(isZero(norm, feastol) && LT(nnewrowidx / R(numRows()), percenttoadd, _solver.epsilon())) 1457 { 1458 updaterows.add(_transformedRows.lhs(rowNumber), _transformedRows.rowVector(rowNumber), 1459 _transformedRows.rhs(rowNumber)); 1460 1461 _decompReducedProbRows[rowNumber] = true; 1462 newrowidx[nnewrowidx] = rowNumber; 1463 nnewrowidx++; 1464 } 1465 1466 } 1467 else 1468 { 1469 // adding all violated rows 1470 updaterows.add(_transformedRows.lhs(rowNumber), _transformedRows.rowVector(rowNumber), 1471 _transformedRows.rhs(rowNumber)); 1472 1473 _decompReducedProbRows[rowNumber] = true; 1474 newrowidx[nnewrowidx] = rowNumber; 1475 nnewrowidx++; 1476 } 1477 } 1478 1479 // if no violated row is found during the first pass of the available rows, then we add all violated rows. 1480 if(nnewrowidx == 0) 1481 { 1482 for(int i = 0; i < nrowstoadd; i++) 1483 { 1484 rowNumber = _decompViolatedRows[i]; 1485 1486 updaterows.add(_transformedRows.lhs(rowNumber), _transformedRows.rowVector(rowNumber), 1487 _transformedRows.rhs(rowNumber)); 1488 1489 _decompReducedProbRows[rowNumber] = true; 1490 newrowidx[nnewrowidx] = rowNumber; 1491 nnewrowidx++; 1492 } 1493 } 1494 1495 // we will always add the row that is deemed best based upon the row norm. 1496 // TODO: check whether this should be skipped if a row has already been added. 1497 if(!allrows && bestrow >= 0) 1498 { 1499 updaterows.add(_transformedRows.lhs(bestrow), _transformedRows.rowVector(bestrow), 1500 _transformedRows.rhs(bestrow)); 1501 1502 _decompReducedProbRows[bestrow] = true; 1503 newrowidx[nnewrowidx] = bestrow; 1504 nnewrowidx++; 1505 } 1506 1507 1508 SPxRowId* addedrowids = 0; 1509 spx_alloc(addedrowids, nnewrowidx); 1510 _solver.addRows(addedrowids, updaterows); 1511 1512 for(int i = 0; i < nnewrowidx; i++) 1513 _decompReducedProbRowIDs[newrowidx[i]] = addedrowids[i]; 1514 1515 // freeing allocated memory 1516 spx_free(addedrowids); 1517 spx_free(newrowidx); 1518 } 1519 1520 1521 1522 /// builds the update rows with those violated in the complmentary problem 1523 // A row is violated in the constraint matrix Ax <= b, if b - A_{i}x < 0 1524 // To aid the computation, all violations are translated to <= constraints 1525 template <class R> 1526 void SoPlexBase<R>::_findViolatedRows(R compObjValue, Array<RowViolation>& violatedrows, 1527 int& nviolatedrows) 1528 { 1529 R feastol = realParam(SoPlexBase<R>::FEASTOL); 1530 VectorBase<R> compProbRedcost( 1531 _compSolver.nCols()); // the reduced costs of the complementary problem 1532 VectorBase<R> compProbPrimal( 1533 _compSolver.nCols()); // the primal vector of the complementary problem 1534 VectorBase<R> compProbActivity( 1535 _compSolver.nRows()); // the activity vector of the complementary problem 1536 R compProbSlackVal = 0; 1537 1538 if(boolParam(SoPlexBase<R>::USECOMPDUAL)) 1539 { 1540 // Retrieving the slacks for each row. 1541 _compSolver.getRedCostSol(compProbRedcost); 1542 } 1543 else 1544 { 1545 // Retrieving the primal vector for the complementary problem 1546 _compSolver.getPrimalSol(compProbPrimal); 1547 _compSolver.computePrimalActivity(compProbPrimal, compProbActivity); 1548 1549 compProbSlackVal = compProbPrimal[_compSolver.number(_compSlackColId)]; 1550 } 1551 1552 // scanning all rows of the complementary problem for violations. 1553 for(int i = 0; i < _nPrimalRows; i++) 1554 { 1555 LPRowBase<R> origlprow; 1556 DSVectorBase<R> rowtoaddVec(_realLP->nCols()); 1557 R compProbViol = 0; 1558 R compSlackCoeff = 0; 1559 int rownumber = _realLP->number(SPxRowId(_decompPrimalRowIDs[i])); 1560 int comprownum = _compSolver.number(SPxRowId(_decompPrimalRowIDs[i])); 1561 1562 if(!_decompReducedProbRows[rownumber]) 1563 { 1564 // retreiving the violation of the complementary problem primal constraints 1565 if(boolParam(SoPlexBase<R>::USECOMPDUAL)) 1566 { 1567 compSlackCoeff = getCompSlackVarCoeff(i); 1568 compProbViol = compProbRedcost[_compSolver.number(SPxColId( 1569 _decompDualColIDs[i]))]; // this is b - Ax 1570 // subtracting the slack variable value 1571 compProbViol += compObjValue * compSlackCoeff; // must add on the slack variable value. 1572 compProbViol *= compSlackCoeff; // translating the violation to a <= constraint 1573 } 1574 else 1575 { 1576 R viol = _compSolver.rhs(comprownum) - (compProbActivity[comprownum] + compProbSlackVal); 1577 1578 if(viol < 0.0) 1579 compProbViol = viol; 1580 1581 viol = (compProbActivity[comprownum] - compProbSlackVal) - _compSolver.lhs(comprownum); 1582 1583 if(viol < 0.0) 1584 compProbViol = viol; 1585 1586 } 1587 1588 // NOTE: if the row was originally a ranged constraint, we are only interest in one of the inequalities. 1589 // If one inequality of the range violates the bounds, then we will add the row. 1590 1591 // the translation of the complementary primal problem to the dual some rows resulted in two columns. 1592 if(boolParam(SoPlexBase<R>::USECOMPDUAL) && i < _nPrimalRows - 1 && 1593 _realLP->number(SPxRowId(_decompPrimalRowIDs[i])) == _realLP->number(SPxRowId( 1594 _decompPrimalRowIDs[i + 1]))) 1595 { 1596 i++; 1597 compSlackCoeff = getCompSlackVarCoeff(i); 1598 R tempViol = compProbRedcost[_compSolver.number(SPxColId(_decompDualColIDs[i]))]; // this is b - Ax 1599 tempViol += compObjValue * compSlackCoeff; 1600 tempViol *= compSlackCoeff; 1601 1602 // if the other side of the range constraint has a larger violation, then this is used for the 1603 // computation. 1604 if(tempViol < compProbViol) 1605 compProbViol = tempViol; 1606 } 1607 1608 1609 // checking the violation of the row. 1610 if(LT(compProbViol, (R) 0, feastol)) 1611 { 1612 if(!_decompReducedProbRows[rownumber]) 1613 { 1614 numIncludedRows++; 1615 assert(numIncludedRows <= _realLP->nRows()); 1616 } 1617 1618 violatedrows[nviolatedrows].idx = rownumber; 1619 violatedrows[nviolatedrows].violation = spxAbs(compProbViol); 1620 nviolatedrows++; 1621 } 1622 } 1623 } 1624 } 1625 1626 1627 1628 /// identifies the columns of the row-form basis that correspond to rows with zero dual multipliers. 1629 // This function assumes that the basis is in the row form. 1630 // @todo extend this to the case when the basis is in the column form. 1631 template <class R> 1632 void SoPlexBase<R>::_getZeroDualMultiplierIndices(VectorBase<R> feasVector, int* nonposind, 1633 int* colsforremoval, 1634 int* nnonposind, bool& stop) 1635 { 1636 assert(_solver.rep() == SPxSolverBase<R>::ROW); 1637 1638 #ifdef NO_TOL 1639 R feastol = 0.0; 1640 #else 1641 #ifdef SOPLEX_USE_FEASTOL 1642 R feastol = realParam(SoPlexBase<R>::FEASTOL); 1643 #else 1644 R feastol = realParam(SoPlexBase<R>::EPSILON_ZERO); 1645 #endif 1646 #endif 1647 1648 bool delCol; 1649 1650 _decompReducedProbColIDs.reSize(numCols()); 1651 *nnonposind = 0; 1652 1653 // iterating over all columns in the basis matrix 1654 // this identifies the basis indices and the indicies that are positive. 1655 for(int i = 0; i < _solver.nCols(); ++i) 1656 { 1657 _decompReducedProbCols[i] = true; 1658 _decompReducedProbColIDs[i].inValidate(); 1659 colsforremoval[i] = i; 1660 delCol = false; 1661 1662 if(_solver.basis().baseId(i).isSPxRowId()) // find the row id's for rows in the basis 1663 { 1664 // record the row if the dual multiple is zero. 1665 if(isZero(feasVector[i], feastol)) 1666 { 1667 nonposind[*nnonposind] = i; 1668 (*nnonposind)++; 1669 1670 // NOTE: commenting out the delCol flag at this point. The colsforremoval array should indicate the 1671 // columns that have a zero reduced cost. Hence, the delCol flag should only be set in the isSPxColId 1672 // branch of the if statement. 1673 //delCol = true; 1674 } 1675 } 1676 else if(_solver.basis().baseId( 1677 i).isSPxColId()) // get the column id's for the columns in the basis 1678 { 1679 if(isZero(feasVector[i], feastol)) 1680 { 1681 nonposind[*nnonposind] = i; 1682 (*nnonposind)++; 1683 1684 delCol = true; 1685 } 1686 } 1687 1688 // setting an array to identify the columns to be removed from the LP to form the reduced problem 1689 if(delCol) 1690 { 1691 colsforremoval[i] = -1; 1692 _decompReducedProbCols[i] = false; 1693 } 1694 else if(_solver.basis().baseId(i).isSPxColId()) 1695 { 1696 if(_solver.basis().baseId(i).isSPxColId()) 1697 _decompReducedProbColIDs[_solver.number(_solver.basis().baseId(i))] = SPxColId( 1698 _solver.basis().baseId(i)); 1699 else 1700 _decompReducedProbCols[_solver.number(_solver.basis().baseId(i))] = false; 1701 } 1702 } 1703 1704 stop = decompTerminate(realParam(SoPlexBase<R>::TIMELIMIT) * SOPLEX_TIMELIMIT_FRAC); 1705 } 1706 1707 1708 1709 /// retrieves the compatible columns from the constraint matrix 1710 // This function also updates the constraint matrix of the reduced problem. It is efficient to perform this in the 1711 // following function because the required linear algebra has been performed. 1712 template <class R> 1713 void SoPlexBase<R>::_getCompatibleColumns(VectorBase<R> feasVector, int* nonposind, int* compatind, 1714 int* rowsforremoval, 1715 int* colsforremoval, int nnonposind, int* ncompatind, bool formRedProb, bool& stop) 1716 { 1717 #ifdef NO_TOL 1718 R feastol = 0.0; 1719 #else 1720 #ifdef SOPLEX_USE_FEASTOL 1721 R feastol = realParam(SoPlexBase<R>::FEASTOL); 1722 #else 1723 R feastol = realParam(SoPlexBase<R>::EPSILON_ZERO); 1724 #endif 1725 #endif 1726 1727 bool compatible; 1728 SSVectorBase<R> y(_solver.nCols(), _solver.tolerances()); 1729 y.unSetup(); 1730 1731 *ncompatind = 0; 1732 1733 #ifndef NDEBUG 1734 int bind = 0; 1735 bool* activerows = 0; 1736 spx_alloc(activerows, numRows()); 1737 1738 for(int i = 0; i < numRows(); ++i) 1739 activerows[i] = false; 1740 1741 for(int i = 0; i < numCols(); ++i) 1742 { 1743 if(_solver.basis().baseId(i).isSPxRowId()) // find the row id's for rows in the basis 1744 { 1745 if(!isZero(feasVector[i], feastol)) 1746 { 1747 bind = _realLP->number(SPxRowId(_solver.basis().baseId(i))); // getting the corresponding row 1748 // for the original LP. 1749 assert(bind >= 0 && bind < numRows()); 1750 activerows[bind] = true; 1751 } 1752 } 1753 } 1754 1755 #endif 1756 1757 // this function is called from other places where to identify the columns for removal. In these places the 1758 // reduced problem is not formed. 1759 if(formRedProb) 1760 { 1761 _decompReducedProbRowIDs.reSize(_solver.nRows()); 1762 _decompReducedProbColRowIDs.reSize(_solver.nRows()); 1763 } 1764 1765 for(int i = 0; i < numRows(); ++i) 1766 { 1767 rowsforremoval[i] = i; 1768 1769 if(formRedProb) 1770 _decompReducedProbRows[i] = true; 1771 1772 numIncludedRows++; 1773 1774 // the rhs of this calculation are the rows of the constraint matrix 1775 // so we are solving y B = A_{i,.} 1776 // @NOTE: This may not actually be necessary. The solve process is very time consuming and is a point where the 1777 // approach breaks down. It could be simplier if we use a faster solve. Maybe something like: 1778 // Omer, J.; Towhidi, M. & Soumis, F., "The positive edge pricing rule for the dual simplex", 1779 // Computers & Operations Research , 2015, 61, 135-142 1780 try 1781 { 1782 _solver.basis().solve(y, _solver.vector(i)); 1783 } 1784 catch(const SPxException& E) 1785 { 1786 SPX_MSG_ERROR(spxout << "Caught exception <" << E.what() << "> while computing compatability.\n"); 1787 } 1788 1789 compatible = true; 1790 1791 // a compatible row is given by zeros in all columns related to the nonpositive indices 1792 for(int j = 0; j < nnonposind; ++j) 1793 { 1794 // @TODO: getting a tolerance issue with this check. Don't know how to fix it. 1795 if(!isZero(y[nonposind[j]], feastol)) 1796 { 1797 compatible = false; 1798 break; 1799 } 1800 } 1801 1802 // checking that the active rows are compatible 1803 assert(!activerows[i] || compatible); 1804 1805 // changing the matrix coefficients 1806 DSVectorBase<R> newRowVector; 1807 LPRowBase<R> rowtoupdate; 1808 1809 if(y.isSetup()) 1810 { 1811 for(int j = 0; j < y.size(); j++) 1812 // coverity[var_deref_model] 1813 newRowVector.add(y.index(j), y.value(j)); 1814 } 1815 else 1816 { 1817 for(int j = 0; j < numCols(); j++) 1818 { 1819 if(!isZero(y[j], feastol)) 1820 newRowVector.add(j, y[j]); 1821 } 1822 } 1823 1824 // transforming the original problem rows 1825 _solver.getRow(i, rowtoupdate); 1826 1827 #ifndef NO_TRANSFORM 1828 rowtoupdate.setRowVector(newRowVector); 1829 #endif 1830 1831 if(formRedProb) 1832 _transformedRows.add(rowtoupdate); 1833 1834 1835 // Making all equality constraints compatible, i.e. they are included in the reduced problem 1836 if(EQ(rowtoupdate.lhs(), rowtoupdate.rhs(), _solver.epsilon())) 1837 compatible = true; 1838 1839 if(compatible) 1840 { 1841 compatind[*ncompatind] = i; 1842 (*ncompatind)++; 1843 1844 if(formRedProb) 1845 { 1846 _decompReducedProbRowIDs[i] = _solver.rowId(i); 1847 1848 // updating the compatible row 1849 _decompLP->changeRow(i, rowtoupdate); 1850 } 1851 } 1852 else 1853 { 1854 // setting an array to identify the rows to be removed from the LP to form the reduced problem 1855 rowsforremoval[i] = -1; 1856 numIncludedRows--; 1857 1858 if(formRedProb) 1859 _decompReducedProbRows[i] = false; 1860 } 1861 1862 // determine whether the reduced problem setup should be terminated 1863 stop = decompTerminate(realParam(SoPlexBase<R>::TIMELIMIT) * SOPLEX_TIMELIMIT_FRAC); 1864 1865 if(stop) 1866 break; 1867 } 1868 1869 assert(numIncludedRows <= _solver.nRows()); 1870 1871 #ifndef NDEBUG 1872 spx_free(activerows); 1873 #endif 1874 } 1875 1876 1877 1878 /// computes the reduced problem objective coefficients 1879 template <class R> 1880 void SoPlexBase<R>::_computeReducedProbObjCoeff(bool& stop) 1881 { 1882 #ifdef NO_TOL 1883 R feastol = 0.0; 1884 #else 1885 #ifdef SOPLEX_USE_FEASTOL 1886 R feastol = realParam(SoPlexBase<R>::FEASTOL); 1887 #else 1888 R feastol = realParam(SoPlexBase<R>::EPSILON_ZERO); 1889 #endif 1890 #endif 1891 1892 SSVectorBase<R> y(numCols(), _solver.tolerances()); 1893 y.unSetup(); 1894 1895 // the rhs of this calculation is the original objective coefficient vector 1896 // so we are solving y B = c 1897 try 1898 { 1899 _solver.basis().solve(y, _solver.maxObj()); 1900 } 1901 catch(const SPxException& E) 1902 { 1903 SPX_MSG_ERROR(spxout << "Caught exception <" << E.what() << "> while computing compatability.\n"); 1904 } 1905 1906 _transformedObj.reDim(numCols()); 1907 1908 if(y.isSetup()) 1909 { 1910 int ycount = 0; 1911 1912 for(int i = 0; i < numCols(); i++) 1913 { 1914 if(ycount < y.size() && i == y.index(ycount)) 1915 { 1916 _transformedObj[i] = y.value(ycount); 1917 ycount++; 1918 } 1919 else 1920 _transformedObj[i] = 0.0; 1921 } 1922 } 1923 else 1924 { 1925 for(int i = 0; i < numCols(); i++) 1926 { 1927 if(isZero(y[i], feastol)) 1928 _transformedObj[i] = 0.0; 1929 else 1930 _transformedObj[i] = y[i]; 1931 } 1932 } 1933 1934 // setting the updated objective vector 1935 #ifndef NO_TRANSFORM 1936 _decompLP->changeObj(_transformedObj); 1937 #endif 1938 1939 // determine whether the reduced problem setup should be terminated 1940 stop = decompTerminate(realParam(SoPlexBase<R>::TIMELIMIT) * SOPLEX_TIMELIMIT_FRAC); 1941 } 1942 1943 1944 1945 /// computes the compatible bound constraints and adds them to the reduced problem 1946 // NOTE: No columns are getting removed from the reduced problem. Only the bound constraints are being removed. 1947 // So in the reduced problem, all variables are free unless the bound constraints are selected as compatible. 1948 template <class R> 1949 void SoPlexBase<R>::_getCompatibleBoundCons(LPRowSetBase<R>& boundcons, int* compatboundcons, 1950 int* nonposind, 1951 int* ncompatboundcons, int nnonposind, bool& stop) 1952 { 1953 #ifdef NO_TOL 1954 R feastol = 0.0; 1955 #else 1956 #ifdef SOPLEX_USE_FEASTOL 1957 R feastol = realParam(SoPlexBase<R>::FEASTOL); 1958 #else 1959 R feastol = realParam(SoPlexBase<R>::EPSILON_ZERO); 1960 #endif 1961 #endif 1962 1963 bool compatible; 1964 int ncols = numCols(); 1965 SSVectorBase<R> y(ncols, _solver.tolerances()); 1966 y.unSetup(); 1967 1968 _decompReducedProbColRowIDs.reSize(numCols()); 1969 1970 1971 // identifying the compatible bound constraints 1972 for(int i = 0; i < ncols; i++) 1973 { 1974 _decompReducedProbColRowIDs[i].inValidate(); 1975 1976 // setting all variables to free variables. 1977 // Bound constraints will only be added to the variables with compatible bound constraints. 1978 _decompLP->changeUpper(i, R(infinity)); 1979 _decompLP->changeLower(i, R(-infinity)); 1980 1981 // the rhs of this calculation are the unit vectors of the bound constraints 1982 // so we are solving y B = I_{i,.} 1983 // this solve should not be required. We only need the column of the basis inverse. 1984 try 1985 { 1986 _solver.basis().solve(y, _solver.unitVector(i)); 1987 } 1988 catch(const SPxException& E) 1989 { 1990 SPX_MSG_ERROR(spxout << "Caught exception <" << E.what() << "> while computing compatability.\n"); 1991 } 1992 1993 compatible = true; 1994 1995 // a compatible row is given by zeros in all columns related to the nonpositive indices 1996 for(int j = 0; j < nnonposind; ++j) 1997 { 1998 // @todo really need to check this part of the code. Run through this with Ambros or Matthias. 1999 if(!isZero(y[nonposind[j]], _solver.epsilon())) 2000 { 2001 compatible = false; 2002 break; 2003 } 2004 } 2005 2006 // changing the matrix coefficients 2007 DSVectorBase<R> newRowVector; 2008 LPRowBase<R> rowtoupdate; 2009 2010 #ifndef NO_TRANSFORM 2011 2012 if(y.isSetup()) 2013 { 2014 for(int j = 0; j < y.size(); j++) 2015 newRowVector.add(y.index(j), y.value(j)); 2016 } 2017 else 2018 { 2019 for(int j = 0; j < ncols; j++) 2020 { 2021 if(!isZero(y[j], feastol)) 2022 { 2023 newRowVector.add(j, y[j]); 2024 } 2025 } 2026 } 2027 2028 #else 2029 newRowVector.add(i, 1.0); 2030 #endif 2031 2032 // will probably need to map this set of rows. 2033 _transformedRows.add(_solver.lower(i), newRowVector, _solver.upper(i)); 2034 2035 // variable bounds are compatible in the current implementation. 2036 // this is still function is still required to transform the bound constraints with respect to the basis matrix. 2037 compatible = true; 2038 2039 // if the bound constraint is compatible, it remains in the reduced problem. 2040 // Otherwise the bound is removed and the variables are free. 2041 if(compatible) 2042 { 2043 R lhs = R(-infinity); 2044 R rhs = R(infinity); 2045 2046 if(_solver.lower(i) > R(-infinity)) 2047 lhs = _solver.lower(i); 2048 2049 if(_solver.upper(i) < R(infinity)) 2050 rhs = _solver.upper(i); 2051 2052 if(lhs > R(-infinity) || rhs < R(infinity)) 2053 { 2054 compatboundcons[(*ncompatboundcons)] = i; 2055 (*ncompatboundcons)++; 2056 2057 boundcons.add(lhs, newRowVector, rhs); 2058 } 2059 2060 2061 } 2062 2063 // determine whether the reduced problem setup should be terminated 2064 stop = decompTerminate(realParam(SoPlexBase<R>::TIMELIMIT) * SOPLEX_TIMELIMIT_FRAC); 2065 2066 if(stop) 2067 break; 2068 } 2069 2070 } 2071 2072 2073 2074 /// computes the rows to remove from the complementary problem 2075 template <class R> 2076 void SoPlexBase<R>::_getRowsForRemovalComplementaryProblem(int* nonposind, int* bind, 2077 int* rowsforremoval, 2078 int* nrowsforremoval, int nnonposind) 2079 { 2080 *nrowsforremoval = 0; 2081 2082 for(int i = 0; i < nnonposind; i++) 2083 { 2084 if(bind[nonposind[i]] < 0) 2085 { 2086 rowsforremoval[*nrowsforremoval] = -1 - bind[nonposind[i]]; 2087 (*nrowsforremoval)++; 2088 } 2089 } 2090 } 2091 2092 2093 2094 /// removing rows from the complementary problem. 2095 // the rows that are removed from decompCompLP are the rows from the reduced problem that have a non-positive dual 2096 // multiplier in the optimal solution. 2097 template <class R> 2098 void SoPlexBase<R>::_deleteAndUpdateRowsComplementaryProblem(SPxRowId rangedRowIds[], 2099 int& naddedrows) 2100 { 2101 DSVectorBase<R> slackColCoeff; 2102 2103 // setting the objective coefficients of the original variables to zero 2104 VectorBase<R> newObjCoeff(numCols()); 2105 2106 for(int i = 0; i < numCols(); i++) 2107 { 2108 _compSolver.changeBounds(_realLP->cId(i), R(-infinity), R(infinity)); 2109 newObjCoeff[i] = 0; 2110 } 2111 2112 _compSolver.changeObj(newObjCoeff); 2113 2114 // adding the slack column to the complementary problem 2115 SPxColId* addedcolid = 0; 2116 2117 if(boolParam(SoPlexBase<R>::USECOMPDUAL)) 2118 { 2119 spx_alloc(addedcolid, 1); 2120 LPColSetBase<R> compSlackCol; 2121 compSlackCol.add(R(1.0), R(-infinity), slackColCoeff, R(infinity)); 2122 _compSolver.addCols(addedcolid, compSlackCol); 2123 _compSlackColId = addedcolid[0]; 2124 } 2125 else 2126 { 2127 LPRowSetBase<R> 2128 addrangedrows; // the row set of ranged and equality rows that must be added to the complementary problem. 2129 naddedrows = 0; 2130 2131 // finding all of the ranged and equality rows and creating two <= constraints. 2132 for(int i = 0; i < numRows(); i++) 2133 { 2134 if(_realLP->rowType(i) == LPRowBase<R>::RANGE || _realLP->rowType(i) == LPRowBase<R>::EQUAL) 2135 { 2136 assert(_compSolver.lhs(i) > R(-infinity)); 2137 assert(_compSolver.rowType(i) == LPRowBase<R>::RANGE 2138 || _compSolver.rowType(i) == LPRowBase<R>::EQUAL); 2139 2140 _compSolver.changeLhs(i, R(-infinity)); 2141 addrangedrows.add(_realLP->lhs(i), _realLP->rowVector(i), R(infinity)); 2142 naddedrows++; 2143 } 2144 } 2145 2146 // adding the rows for the ranged rows to make <= conatraints 2147 SPxRowId* addedrowid = 0; 2148 spx_alloc(addedrowid, naddedrows); 2149 _compSolver.addRows(addedrowid, addrangedrows); 2150 2151 // collecting the row ids 2152 for(int i = 0; i < naddedrows; i++) 2153 rangedRowIds[i] = addedrowid[i]; 2154 2155 spx_free(addedrowid); 2156 2157 2158 // adding the slack column 2159 spx_alloc(addedcolid, 1); 2160 LPColSetBase<R> compSlackCol; 2161 compSlackCol.add(R(-1.0), R(0.0), slackColCoeff, R(infinity)); 2162 2163 _compSolver.addCols(addedcolid, compSlackCol); 2164 _compSlackColId = addedcolid[0]; 2165 } 2166 2167 // freeing allocated memory 2168 spx_free(addedcolid); 2169 } 2170 2171 2172 2173 /// update the dual complementary problem with additional columns and rows 2174 // Given the solution to the updated reduced problem, the complementary problem will be updated with modifications to 2175 // the constraints and the removal of variables 2176 template <class R> 2177 void SoPlexBase<R>::_updateDecompComplementaryDualProblem(bool origObj) 2178 { 2179 R feastol = realParam(SoPlexBase<R>::FEASTOL); 2180 2181 int prevNumCols = 2182 _compSolver.nCols(); // number of columns in the previous formulation of the complementary prob 2183 int prevNumRows = _compSolver.nRows(); 2184 int prevPrimalRowIds = _nPrimalRows; 2185 int prevDualColIds = _nDualCols; 2186 2187 LPColSetBase<R> addElimCols(_nElimPrimalRows); // columns previously eliminated from the 2188 // complementary problem that must be added 2189 int numElimColsAdded = 0; 2190 int currnumrows = prevNumRows; 2191 // looping over all rows from the original LP that were eliminated during the formation of the complementary 2192 // problem. The eliminated rows will be added if they are basic in the reduced problem. 2193 2194 for(int i = 0; i < _nElimPrimalRows; i++) 2195 { 2196 int rowNumber = _realLP->number(_decompElimPrimalRowIDs[i]); 2197 2198 int solverRowNum = _solver.number(_decompReducedProbRowIDs[rowNumber]); 2199 assert(solverRowNum >= 0 && solverRowNum < _solver.nRows()); 2200 2201 // checking for the basic rows in the reduced problem 2202 if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_UPPER || 2203 _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_LOWER || 2204 _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_FIXED || 2205 _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_FREE || 2206 (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_LOWER && 2207 LE(_solver.rhs(solverRowNum) - _solver.pVec()[solverRowNum], R(0.0), feastol)) || 2208 (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_UPPER && 2209 LE(_solver.pVec()[solverRowNum] - _solver.lhs(solverRowNum), R(0.0), feastol))) 2210 { 2211 LPRowBase<R> origlprow; 2212 DSVectorBase<R> coltoaddVec(_realLP->nCols()); 2213 2214 LPRowSetBase<R> additionalrows; 2215 int nnewrows = 0; 2216 2217 _realLP->getRow(rowNumber, origlprow); 2218 2219 for(int j = 0; j < origlprow.rowVector().size(); j++) 2220 { 2221 // the column of the new row may not exist in the current complementary problem. 2222 // if the column does not exist, then it is necessary to create the column. 2223 int colNumber = origlprow.rowVector().index(j); 2224 2225 if(_decompCompProbColIDsIdx[colNumber] == -1) 2226 { 2227 assert(!_decompReducedProbColIDs[colNumber].isValid()); 2228 _decompPrimalColIDs[_nPrimalCols] = _realLP->cId(colNumber); 2229 _decompCompProbColIDsIdx[colNumber] = _nPrimalCols; 2230 _fixedOrigVars[colNumber] = -2; 2231 _nPrimalCols++; 2232 2233 // all columns for the complementary problem are converted to unrestricted. 2234 additionalrows.create(1, _realLP->maxObj(colNumber), _realLP->maxObj(colNumber)); 2235 nnewrows++; 2236 2237 coltoaddVec.add(currnumrows, origlprow.rowVector().value(j)); 2238 currnumrows++; 2239 } 2240 else 2241 coltoaddVec.add(_compSolver.number(_decompDualRowIDs[_decompCompProbColIDsIdx[colNumber]]), 2242 origlprow.rowVector().value(j)); 2243 } 2244 2245 SPxRowId* addedrowids = 0; 2246 spx_alloc(addedrowids, nnewrows); 2247 _compSolver.addRows(addedrowids, additionalrows); 2248 2249 for(int j = 0; j < nnewrows; j++) 2250 { 2251 _decompDualRowIDs[_nDualRows] = addedrowids[j]; 2252 _nDualRows++; 2253 } 2254 2255 spx_free(addedrowids); 2256 2257 2258 2259 if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_UPPER || 2260 _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_FIXED || 2261 _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_FREE || 2262 (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_LOWER && 2263 LE(_solver.rhs(solverRowNum) - _solver.pVec()[solverRowNum], R(0.0), feastol))) 2264 { 2265 assert(_realLP->rhs(_decompElimPrimalRowIDs[i]) < R(infinity)); 2266 addElimCols.add(_realLP->rhs(_decompElimPrimalRowIDs[i]), R(-infinity), coltoaddVec, R(infinity)); 2267 2268 if(_nPrimalRows >= _decompPrimalRowIDs.size()) 2269 { 2270 _decompPrimalRowIDs.reSize(_nPrimalRows * 2); 2271 _decompDualColIDs.reSize(_nPrimalRows * 2); 2272 } 2273 2274 _decompPrimalRowIDs[_nPrimalRows] = _decompElimPrimalRowIDs[i]; 2275 _nPrimalRows++; 2276 2277 _decompElimPrimalRowIDs.remove(i); 2278 _nElimPrimalRows--; 2279 i--; 2280 2281 numElimColsAdded++; 2282 } 2283 else if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_LOWER || 2284 (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_UPPER && 2285 LE(_solver.pVec()[solverRowNum] - _solver.lhs(solverRowNum), R(0.0), feastol))) 2286 { 2287 // this assert should stay, but there is an issue with the status and the dual vector 2288 //assert(LT(dualVector[_solver.number(_decompReducedProbRowIDs[rowNumber])], 0.0)); 2289 assert(_realLP->lhs(_decompElimPrimalRowIDs[i]) > R(-infinity)); 2290 addElimCols.add(_realLP->lhs(_decompElimPrimalRowIDs[i]), R(-infinity), coltoaddVec, R(infinity)); 2291 2292 _decompPrimalRowIDs[_nPrimalRows] = _decompElimPrimalRowIDs[i]; 2293 _nPrimalRows++; 2294 2295 _decompElimPrimalRowIDs.remove(i); 2296 _nElimPrimalRows--; 2297 i--; 2298 2299 numElimColsAdded++; 2300 } 2301 } 2302 } 2303 2304 SPX_MSG_INFO2(spxout, spxout << "Number of eliminated columns added to complementary problem: " 2305 << numElimColsAdded << std::endl); 2306 2307 // updating the _decompDualColIDs with the additional columns from the eliminated rows. 2308 _compSolver.addCols(addElimCols); 2309 2310 for(int i = prevNumCols; i < _compSolver.nCols(); i++) 2311 { 2312 _decompDualColIDs[prevDualColIds + i - prevNumCols] = _compSolver.colId(i); 2313 _nDualCols++; 2314 } 2315 2316 assert(_nDualCols == _nPrimalRows); 2317 2318 // looping over all rows from the original problem that were originally contained in the complementary problem. 2319 // The basic rows will be set as free variables, the non-basic rows will be eliminated from the complementary 2320 // problem. 2321 DSVectorBase<R> slackRowCoeff(_compSolver.nCols()); 2322 2323 int* colsforremoval = 0; 2324 int ncolsforremoval = 0; 2325 spx_alloc(colsforremoval, prevPrimalRowIds); 2326 2327 for(int i = 0; i < prevPrimalRowIds; i++) 2328 { 2329 int rowNumber = _realLP->number(_decompPrimalRowIDs[i]); 2330 2331 // this loop runs over all rows previously in the complementary problem. If rows are added to the reduced 2332 // problem, they will be transfered from the incompatible set to the compatible set in the following if 2333 // statement. 2334 if(_decompReducedProbRows[rowNumber]) 2335 { 2336 // rows added to the reduced problem may have been equality constriants. The equality constraints from the 2337 // original problem are converted into <= and >= constraints. Upon adding these constraints to the reduced 2338 // problem, only a single dual column is needed in the complementary problem. Hence, one of the dual columns 2339 // is removed. 2340 // 2341 // 22.06.2015 Testing keeping all constraints in the complementary problem. 2342 // This requires a dual column to be fixed to zero for the range and equality rows. 2343 #ifdef KEEP_ALL_ROWS_IN_COMP_PROB 2344 bool incrementI = false; 2345 #endif 2346 2347 if(i + 1 < prevPrimalRowIds 2348 && _realLP->number(_decompPrimalRowIDs[i]) == _realLP->number(_decompPrimalRowIDs[i + 1])) 2349 { 2350 assert(_decompPrimalRowIDs[i].idx == _decompPrimalRowIDs[i + 1].idx); 2351 2352 #ifdef KEEP_ALL_ROWS_IN_COMP_PROB // 22.06.2015 2353 2354 if(_realLP->rowType(_decompPrimalRowIDs[i]) == LPRowBase<R>::RANGE) 2355 { 2356 _compSolver.changeObj(_decompDualColIDs[i + 1], 0.0); 2357 _compSolver.changeBounds(_decompDualColIDs[i + 1], 0.0, 0.0); 2358 } 2359 2360 incrementI = true; 2361 #else 2362 colsforremoval[ncolsforremoval] = _compSolver.number(SPxColId(_decompDualColIDs[i + 1])); 2363 ncolsforremoval++; 2364 2365 _decompPrimalRowIDs.remove(i + 1); 2366 _nPrimalRows--; 2367 _decompDualColIDs.remove(i + 1); 2368 _nDualCols--; 2369 2370 prevPrimalRowIds--; 2371 #endif 2372 } 2373 2374 //assert(i + 1 == prevPrimalRowIds || _decompPrimalRowIDs[i].idx != _decompPrimalRowIDs[i+1].idx); 2375 2376 int solverRowNum = _solver.number(_decompReducedProbRowIDs[rowNumber]); 2377 assert(solverRowNum >= 0 && solverRowNum < _solver.nRows()); 2378 2379 if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_UPPER || 2380 _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_FIXED || 2381 _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_FREE || 2382 (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_LOWER && 2383 LE(_solver.rhs(solverRowNum) - _solver.pVec()[solverRowNum], R(0.0), feastol))) 2384 { 2385 _compSolver.changeObj(_decompDualColIDs[i], _realLP->rhs(SPxRowId(_decompPrimalRowIDs[i]))); 2386 _compSolver.changeBounds(_decompDualColIDs[i], R(-infinity), R(infinity)); 2387 } 2388 else if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_LOWER || 2389 (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_UPPER && 2390 LE(_solver.pVec()[solverRowNum] - _solver.lhs(solverRowNum), R(0.0), feastol))) 2391 { 2392 //assert(LT(dualVector[solverRowNum], 0.0)); 2393 _compSolver.changeObj(_decompDualColIDs[i], _realLP->lhs(SPxRowId(_decompPrimalRowIDs[i]))); 2394 _compSolver.changeBounds(_decompDualColIDs[i], R(-infinity), R(infinity)); 2395 } 2396 else //if ( _solver.basis().desc().rowStatus(solverRowNum) != SPxBasisBase<R>::Desc::D_FREE ) 2397 { 2398 //assert(isZero(dualVector[solverRowNum], 0.0)); 2399 2400 // 22.06.2015 Testing keeping all rows in the complementary problem 2401 #ifdef KEEP_ALL_ROWS_IN_COMP_PROB 2402 switch(_realLP->rowType(_decompPrimalRowIDs[i])) 2403 { 2404 case LPRowBase<R>::RANGE: 2405 assert(_realLP->number(SPxColId(_decompPrimalRowIDs[i])) == 2406 _realLP->number(SPxColId(_decompPrimalRowIDs[i + 1]))); 2407 assert(_compSolver.obj(_compSolver.number(SPxColId(_decompDualColIDs[i]))) != 2408 _compSolver.obj(_compSolver.number(SPxColId(_decompDualColIDs[i + 1])))); 2409 2410 _compSolver.changeObj(_decompDualColIDs[i], _realLP->rhs(_decompPrimalRowIDs[i])); 2411 _compSolver.changeBounds(_decompDualColIDs[i], 0.0, R(infinity)); 2412 _compSolver.changeObj(_decompDualColIDs[i + 1], _realLP->lhs(_decompPrimalRowIDs[i])); 2413 _compSolver.changeBounds(_decompDualColIDs[i + 1], R(-infinity), 0.0); 2414 2415 i++; 2416 break; 2417 2418 case LPRowBase<R>::GREATER_EQUAL: 2419 _compSolver.changeObj(_decompDualColIDs[i], _realLP->lhs(_decompPrimalRowIDs[i])); 2420 _compSolver.changeBounds(_decompDualColIDs[i], R(-infinity), 0.0); 2421 break; 2422 2423 case LPRowBase<R>::LESS_EQUAL: 2424 _compSolver.changeObj(_decompDualColIDs[i], _realLP->rhs(_decompPrimalRowIDs[i])); 2425 _compSolver.changeBounds(_decompDualColIDs[i], 0.0, R(infinity)); 2426 break; 2427 2428 default: 2429 throw SPxInternalCodeException("XDECOMPSL01 This should never happen."); 2430 } 2431 2432 #else // 22.06.2015 testing keeping all rows in the complementary problem 2433 colsforremoval[ncolsforremoval] = _compSolver.number(SPxColId(_decompDualColIDs[i])); 2434 ncolsforremoval++; 2435 2436 if(_nElimPrimalRows >= _decompElimPrimalRowIDs.size()) 2437 _decompElimPrimalRowIDs.reSize(_realLP->nRows()); 2438 2439 _decompElimPrimalRowIDs[_nElimPrimalRows] = _decompPrimalRowIDs[i]; 2440 _nElimPrimalRows++; 2441 _decompPrimalRowIDs.remove(i); 2442 _nPrimalRows--; 2443 _decompDualColIDs.remove(i); 2444 _nDualCols--; 2445 2446 i--; 2447 prevPrimalRowIds--; 2448 #endif 2449 } 2450 2451 #ifdef KEEP_ALL_ROWS_IN_COMP_PROB 2452 2453 if(incrementI) 2454 i++; 2455 2456 #endif 2457 } 2458 else 2459 { 2460 switch(_realLP->rowType(_decompPrimalRowIDs[i])) 2461 { 2462 case LPRowBase<R>::RANGE: 2463 assert(_realLP->number(SPxColId(_decompPrimalRowIDs[i])) == 2464 _realLP->number(SPxColId(_decompPrimalRowIDs[i + 1]))); 2465 assert(_compSolver.obj(_compSolver.number(SPxColId(_decompDualColIDs[i]))) != 2466 _compSolver.obj(_compSolver.number(SPxColId(_decompDualColIDs[i + 1])))); 2467 2468 if(_compSolver.obj(_compSolver.number(SPxColId(_decompDualColIDs[i]))) < 2469 _compSolver.obj(_compSolver.number(SPxColId(_decompDualColIDs[i + 1])))) 2470 { 2471 slackRowCoeff.add(_compSolver.number(SPxColId(_decompDualColIDs[i])), -SOPLEX_SLACKCOEFF); 2472 slackRowCoeff.add(_compSolver.number(SPxColId(_decompDualColIDs[i + 1])), SOPLEX_SLACKCOEFF); 2473 } 2474 else 2475 { 2476 slackRowCoeff.add(_compSolver.number(SPxColId(_decompDualColIDs[i])), SOPLEX_SLACKCOEFF); 2477 slackRowCoeff.add(_compSolver.number(SPxColId(_decompDualColIDs[i + 1])), -SOPLEX_SLACKCOEFF); 2478 } 2479 2480 2481 i++; 2482 break; 2483 2484 case LPRowBase<R>::EQUAL: 2485 assert(_realLP->number(SPxColId(_decompPrimalRowIDs[i])) == 2486 _realLP->number(SPxColId(_decompPrimalRowIDs[i + 1]))); 2487 2488 slackRowCoeff.add(_compSolver.number(SPxColId(_decompDualColIDs[i])), SOPLEX_SLACKCOEFF); 2489 slackRowCoeff.add(_compSolver.number(SPxColId(_decompDualColIDs[i + 1])), SOPLEX_SLACKCOEFF); 2490 2491 i++; 2492 break; 2493 2494 case LPRowBase<R>::GREATER_EQUAL: 2495 slackRowCoeff.add(_compSolver.number(SPxColId(_decompDualColIDs[i])), -SOPLEX_SLACKCOEFF); 2496 break; 2497 2498 case LPRowBase<R>::LESS_EQUAL: 2499 slackRowCoeff.add(_compSolver.number(SPxColId(_decompDualColIDs[i])), SOPLEX_SLACKCOEFF); 2500 break; 2501 2502 default: 2503 throw SPxInternalCodeException("XDECOMPSL01 This should never happen."); 2504 } 2505 2506 if(origObj) 2507 { 2508 int numRemove = 1; 2509 int removeCount = 0; 2510 2511 if(_realLP->number(SPxColId(_decompPrimalRowIDs[i])) == 2512 _realLP->number(SPxColId(_decompPrimalRowIDs[i + 1]))) 2513 numRemove++; 2514 2515 do 2516 { 2517 colsforremoval[ncolsforremoval] = _compSolver.number(SPxColId(_decompDualColIDs[i])); 2518 ncolsforremoval++; 2519 2520 if(_nElimPrimalRows >= _decompElimPrimalRowIDs.size()) 2521 _decompElimPrimalRowIDs.reSize(_realLP->nRows()); 2522 2523 _decompElimPrimalRowIDs[_nElimPrimalRows] = _decompPrimalRowIDs[i]; 2524 _nElimPrimalRows++; 2525 _decompPrimalRowIDs.remove(i); 2526 _nPrimalRows--; 2527 _decompDualColIDs.remove(i); 2528 _nDualCols--; 2529 2530 i--; 2531 prevPrimalRowIds--; 2532 2533 removeCount++; 2534 } 2535 while(removeCount < numRemove); 2536 } 2537 } 2538 } 2539 2540 // updating the slack column in the complementary problem 2541 R lhs = 1.0; 2542 R rhs = 1.0; 2543 2544 // it is possible that all rows are included in the reduced problem. In this case, the slack row will be empty. To 2545 // avoid infeasibility, the lhs and rhs are set to 0. 2546 if(slackRowCoeff.size() == 0) 2547 { 2548 lhs = 0.0; 2549 rhs = 0.0; 2550 } 2551 2552 LPRowBase<R> compSlackRow(lhs, slackRowCoeff, rhs); 2553 _compSolver.changeRow(_compSlackDualRowId, compSlackRow); 2554 2555 2556 // if the original objective is used, then all dual columns related to primal rows not in the reduced problem are 2557 // removed from the complementary problem. 2558 // As a result, the slack row becomes an empty row. 2559 int* perm = 0; 2560 spx_alloc(perm, _compSolver.nCols() + numElimColsAdded); 2561 _compSolver.removeCols(colsforremoval, ncolsforremoval, perm); 2562 2563 2564 // updating the dual columns to represent the fixed primal variables. 2565 int* currFixedVars = 0; 2566 spx_alloc(currFixedVars, _realLP->nCols()); 2567 _identifyComplementaryDualFixedPrimalVars(currFixedVars); 2568 _removeComplementaryDualFixedPrimalVars(currFixedVars); 2569 _updateComplementaryDualFixedPrimalVars(currFixedVars); 2570 2571 2572 // freeing allocated memory 2573 spx_free(currFixedVars); 2574 spx_free(perm); 2575 spx_free(colsforremoval); 2576 } 2577 2578 2579 2580 /// update the primal complementary problem with additional columns and rows 2581 // Given the solution to the updated reduced problem, the complementary problem will be updated with modifications to 2582 // the constraints and the removal of variables 2583 template <class R> 2584 void SoPlexBase<R>::_updateDecompComplementaryPrimalProblem(bool origObj) 2585 { 2586 R feastol = realParam(SoPlexBase<R>::FEASTOL); 2587 2588 int prevNumRows = _compSolver.nRows(); 2589 int prevPrimalRowIds = _nPrimalRows; 2590 2591 assert(_nPrimalRows == _nCompPrimalRows); 2592 2593 LPRowSetBase<R> addElimRows(_nElimPrimalRows); // rows previously eliminated from the 2594 // complementary problem that must be added 2595 int numElimRowsAdded = 0; 2596 // looping over all rows from the original LP that were eliminated during the formation of the complementary 2597 // problem. The eliminated rows will be added if they are basic in the reduced problem. 2598 2599 for(int i = 0; i < _nElimPrimalRows; i++) 2600 { 2601 int rowNumber = _realLP->number(_decompElimPrimalRowIDs[i]); 2602 2603 int solverRowNum = _solver.number(_decompReducedProbRowIDs[rowNumber]); 2604 assert(solverRowNum >= 0 && solverRowNum < _solver.nRows()); 2605 2606 // checking the rows that are basic in the reduced problem that should be added to the complementary problem 2607 if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_UPPER 2608 || _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_LOWER 2609 || _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_FIXED 2610 || _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_FREE 2611 || (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_LOWER && 2612 EQ(_solver.rhs(solverRowNum) - _solver.pVec()[solverRowNum], R(0.0), feastol)) 2613 || (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_UPPER && 2614 EQ(_solver.pVec()[solverRowNum] - _solver.lhs(solverRowNum), R(0.0), feastol))) 2615 { 2616 LPRowBase<R> origlprow; 2617 _realLP->getRow(rowNumber, origlprow); 2618 2619 // NOTE: 11.02.2016 I am assuming that all columns from the original problem are contained in the 2620 // complementary problem. Will need to check this. Since nothing is happenning in the 2621 // _deleteAndUpdateRowsComplementaryProblem function, I am feeling confident that all columns remain. 2622 2623 2624 if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_UPPER 2625 || _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_FIXED 2626 || _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_FREE 2627 || (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_LOWER && 2628 EQ(_solver.rhs(solverRowNum) - _solver.pVec()[solverRowNum], R(0.0), feastol))) 2629 { 2630 assert(_realLP->rhs(_decompElimPrimalRowIDs[i]) < R(infinity)); 2631 2632 if(_nPrimalRows >= _decompPrimalRowIDs.size()) 2633 { 2634 _decompPrimalRowIDs.reSize(_nPrimalRows * 2); 2635 _decompCompPrimalRowIDs.reSize(_nPrimalRows * 2); 2636 } 2637 2638 addElimRows.add(_realLP->rhs(_decompElimPrimalRowIDs[i]), origlprow.rowVector(), 2639 _realLP->rhs(_decompElimPrimalRowIDs[i])); 2640 2641 _decompPrimalRowIDs[_nPrimalRows] = _decompElimPrimalRowIDs[i]; 2642 _nPrimalRows++; 2643 2644 _decompElimPrimalRowIDs.remove(i); 2645 _nElimPrimalRows--; 2646 i--; 2647 2648 numElimRowsAdded++; 2649 } 2650 else if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_LOWER 2651 || (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_UPPER && 2652 EQ(_solver.pVec()[solverRowNum] - _solver.lhs(solverRowNum), R(0.0), feastol))) 2653 { 2654 assert(_realLP->lhs(_decompElimPrimalRowIDs[i]) > R(-infinity)); 2655 2656 if(_nPrimalRows >= _decompPrimalRowIDs.size()) 2657 { 2658 _decompPrimalRowIDs.reSize(_nPrimalRows * 2); 2659 _decompCompPrimalRowIDs.reSize(_nPrimalRows * 2); 2660 } 2661 2662 addElimRows.add(_realLP->lhs(_decompElimPrimalRowIDs[i]), origlprow.rowVector(), 2663 _realLP->lhs(_decompElimPrimalRowIDs[i])); 2664 2665 _decompPrimalRowIDs[_nPrimalRows] = _decompElimPrimalRowIDs[i]; 2666 _nPrimalRows++; 2667 2668 _decompElimPrimalRowIDs.remove(i); 2669 _nElimPrimalRows--; 2670 i--; 2671 2672 numElimRowsAdded++; 2673 } 2674 } 2675 } 2676 2677 SPX_MSG_INFO2(spxout, spxout << "Number of eliminated rows added to the complementary problem: " 2678 << numElimRowsAdded << std::endl); 2679 2680 // adding the eliminated rows to the complementary problem. 2681 _compSolver.addRows(addElimRows); 2682 2683 for(int i = prevNumRows; i < _compSolver.nRows(); i++) 2684 { 2685 _decompCompPrimalRowIDs[prevPrimalRowIds + i - prevNumRows] = _compSolver.rowId(i); 2686 _nCompPrimalRows++; 2687 } 2688 2689 assert(_nPrimalRows == _nCompPrimalRows); 2690 2691 2692 // looping over all rows from the original problem that were originally contained in the complementary problem. 2693 // The basic rows will be set as equalities, the non-basic rows will be eliminated from the complementary 2694 // problem. 2695 DSVectorBase<R> slackColCoeff(_compSolver.nRows()); 2696 2697 int* rowsforremoval = 0; 2698 int nrowsforremoval = 0; 2699 spx_alloc(rowsforremoval, prevPrimalRowIds); 2700 2701 for(int i = 0; i < prevPrimalRowIds; i++) 2702 { 2703 int rowNumber = _realLP->number(_decompPrimalRowIDs[i]); 2704 2705 // this loop runs over all rows previously in the complementary problem. If rows are added to the reduced 2706 // problem, they will be transfered from the incompatible set to the compatible set in the following if 2707 // statement. 2708 if(_decompReducedProbRows[rowNumber]) 2709 { 2710 // rows added to the reduced problem may have been equality constriants. The equality constraints from the 2711 // original problem are converted into <= and >= constraints. Upon adding these constraints to the reduced 2712 // problem, only a single dual column is needed in the complementary problem. Hence, one of the dual columns 2713 // is removed. 2714 // 2715 2716 int solverRowNum = _solver.number(_decompReducedProbRowIDs[rowNumber]); 2717 assert(solverRowNum >= 0 && solverRowNum < _solver.nRows()); 2718 2719 if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_UPPER 2720 || _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_FIXED 2721 || _solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_FREE 2722 || (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_LOWER && 2723 EQ(_solver.rhs(solverRowNum) - _solver.pVec()[solverRowNum], R(0.0), feastol))) 2724 { 2725 _compSolver.changeLhs(_decompCompPrimalRowIDs[i], _realLP->rhs(SPxRowId(_decompPrimalRowIDs[i]))); 2726 // need to also update the RHS because a ranged row could have previously been fixed to LOWER 2727 _compSolver.changeRhs(_decompCompPrimalRowIDs[i], _realLP->rhs(SPxRowId(_decompPrimalRowIDs[i]))); 2728 } 2729 else if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_LOWER 2730 || (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_UPPER && 2731 EQ(_solver.pVec()[solverRowNum] - _solver.lhs(solverRowNum), R(0.0), feastol))) 2732 { 2733 _compSolver.changeRhs(_decompCompPrimalRowIDs[i], _realLP->lhs(SPxRowId(_decompPrimalRowIDs[i]))); 2734 // need to also update the LHS because a ranged row could have previously been fixed to UPPER 2735 _compSolver.changeLhs(_decompCompPrimalRowIDs[i], _realLP->lhs(SPxRowId(_decompPrimalRowIDs[i]))); 2736 } 2737 else //if ( _solver.basis().desc().rowStatus(solverRowNum) != SPxBasisBase<R>::Desc::D_FREE ) 2738 { 2739 rowsforremoval[nrowsforremoval] = _compSolver.number(SPxRowId(_decompCompPrimalRowIDs[i])); 2740 nrowsforremoval++; 2741 2742 if(_nElimPrimalRows >= _decompElimPrimalRowIDs.size()) 2743 _decompElimPrimalRowIDs.reSize(_realLP->nRows()); 2744 2745 _decompElimPrimalRowIDs[_nElimPrimalRows] = _decompPrimalRowIDs[i]; 2746 _nElimPrimalRows++; 2747 _decompPrimalRowIDs.remove(i); 2748 _nPrimalRows--; 2749 _decompCompPrimalRowIDs.remove(i); 2750 _nCompPrimalRows--; 2751 2752 i--; 2753 prevPrimalRowIds--; 2754 } 2755 } 2756 else 2757 { 2758 switch(_compSolver.rowType(_decompCompPrimalRowIDs[i])) 2759 { 2760 case LPRowBase<R>::RANGE: 2761 assert(false); 2762 break; 2763 2764 case LPRowBase<R>::EQUAL: 2765 assert(false); 2766 break; 2767 2768 case LPRowBase<R>::LESS_EQUAL: 2769 slackColCoeff.add(_compSolver.number(SPxRowId(_decompCompPrimalRowIDs[i])), -SOPLEX_SLACKCOEFF); 2770 break; 2771 2772 case LPRowBase<R>::GREATER_EQUAL: 2773 slackColCoeff.add(_compSolver.number(SPxRowId(_decompCompPrimalRowIDs[i])), SOPLEX_SLACKCOEFF); 2774 break; 2775 2776 default: 2777 throw SPxInternalCodeException("XDECOMPSL01 This should never happen."); 2778 } 2779 2780 // this is used as a check at the end of the algorithm. If the original objective function is used, then we 2781 // need to remove all unfixed variables. 2782 if(origObj) 2783 { 2784 rowsforremoval[nrowsforremoval] = _compSolver.number(SPxRowId(_decompCompPrimalRowIDs[i])); 2785 nrowsforremoval++; 2786 2787 if(_nElimPrimalRows >= _decompElimPrimalRowIDs.size()) 2788 _decompElimPrimalRowIDs.reSize(_realLP->nRows()); 2789 2790 _decompElimPrimalRowIDs[_nElimPrimalRows] = _decompPrimalRowIDs[i]; 2791 _nElimPrimalRows++; 2792 _decompPrimalRowIDs.remove(i); 2793 _nPrimalRows--; 2794 _decompCompPrimalRowIDs.remove(i); 2795 _nCompPrimalRows--; 2796 2797 i--; 2798 prevPrimalRowIds--; 2799 } 2800 } 2801 } 2802 2803 // updating the slack column in the complementary problem 2804 LPColBase<R> compSlackCol(-1, slackColCoeff, R(infinity), 0.0); 2805 _compSolver.changeCol(_compSlackColId, compSlackCol); 2806 2807 2808 // if the original objective is used, then all complementary rows related to primal rows not in the reduced problem are 2809 // removed from the complementary problem. 2810 // As a result, the slack column becomes an empty column. 2811 int* perm = 0; 2812 spx_alloc(perm, _compSolver.nRows() + numElimRowsAdded); 2813 _compSolver.removeRows(rowsforremoval, nrowsforremoval, perm); 2814 2815 // updating the dual columns to represent the fixed primal variables. 2816 int* currFixedVars = 0; 2817 spx_alloc(currFixedVars, _realLP->nCols()); 2818 _identifyComplementaryPrimalFixedPrimalVars(currFixedVars); 2819 _updateComplementaryPrimalFixedPrimalVars(currFixedVars); 2820 2821 2822 // freeing allocated memory 2823 spx_free(currFixedVars); 2824 spx_free(perm); 2825 spx_free(rowsforremoval); 2826 } 2827 2828 2829 2830 /// checking the optimality of the original problem. 2831 // this function is called if the complementary problem is solved with a non-negative objective value. This implies 2832 // that the rows currently included in the reduced problem are sufficient to identify the optimal solution to the 2833 // original problem. 2834 template <class R> 2835 void SoPlexBase<R>::_checkOriginalProblemOptimality(VectorBase<R> primalVector, bool printViol) 2836 { 2837 SSVectorBase<R> x(_solver.nCols(), this->tolerances()); 2838 x.unSetup(); 2839 2840 // multiplying the solution vector of the reduced problem with the transformed basis to identify the original 2841 // solution vector. 2842 _decompTransBasis.coSolve(x, primalVector); 2843 2844 if(printViol) 2845 { 2846 SPX_MSG_INFO1(spxout, spxout << std::endl 2847 << "Checking consistency between the reduced problem and the original problem." << std::endl); 2848 } 2849 2850 2851 // checking the objective function values of the reduced problem and the original problem. 2852 R redObjVal = 0; 2853 R objectiveVal = 0; 2854 2855 for(int i = 0; i < _solver.nCols(); i++) 2856 { 2857 redObjVal += _solver.maxObj(i) * primalVector[i]; 2858 objectiveVal += _realLP->maxObj(i) * x[i]; 2859 } 2860 2861 if(printViol) 2862 { 2863 SPX_MSG_INFO1(spxout, spxout << "Reduced Problem Objective Value: " << redObjVal << std::endl 2864 << "Original Problem Objective Value: " << objectiveVal << std::endl); 2865 } 2866 2867 _solReal._isPrimalFeasible = true; 2868 _hasSolReal = true; 2869 // get the primal solutions from the reduced problem 2870 _solReal._primal.reDim(_solver.nCols()); 2871 _solReal._primal = x; 2872 2873 R maxviol = 0; 2874 R sumviol = 0; 2875 2876 // checking the bound violations 2877 if(getDecompBoundViolation(maxviol, sumviol)) 2878 { 2879 if(printViol) 2880 SPX_MSG_INFO1(spxout, spxout << "Bound violation - " 2881 << "Max violation: " << maxviol << " Sum violation: " << sumviol << std::endl); 2882 } 2883 2884 _statistics->totalBoundViol = sumviol; 2885 _statistics->maxBoundViol = maxviol; 2886 2887 // checking the row violations 2888 if(getDecompRowViolation(maxviol, sumviol)) 2889 { 2890 if(printViol) 2891 SPX_MSG_INFO1(spxout, spxout << "Row violation - " 2892 << "Max violation: " << maxviol << " Sum violation: " << sumviol << std::endl); 2893 } 2894 2895 _statistics->totalRowViol = sumviol; 2896 _statistics->maxRowViol = maxviol; 2897 2898 if(printViol) 2899 SPX_MSG_INFO1(spxout, spxout << std::endl); 2900 } 2901 2902 2903 2904 /// updating the slack column coefficients to adjust for equality constraints 2905 template <class R> 2906 void SoPlexBase<R>::_updateComplementaryDualSlackColCoeff() 2907 { 2908 // the slack column for the equality constraints is not handled correctly in the dual conversion. Hence, it is 2909 // necessary to change the equality coefficients of the dual row related to the slack column. 2910 for(int i = 0; i < _nPrimalRows; i++) 2911 { 2912 int rowNumber = _realLP->number(SPxRowId(_decompPrimalRowIDs[i])); 2913 2914 if(!_decompReducedProbRows[rowNumber]) 2915 { 2916 if(_realLP->rowType(_decompPrimalRowIDs[i]) == LPRowBase<R>::EQUAL) 2917 { 2918 assert(_realLP->lhs(_decompPrimalRowIDs[i]) == _realLP->rhs(_decompPrimalRowIDs[i])); 2919 _compSolver.changeLower(_decompDualColIDs[i], 2920 0.0); // setting the lower bound of the dual column to zero. 2921 2922 LPColBase<R> addEqualityCol(-_realLP->rhs(_decompPrimalRowIDs[i]), 2923 R(-1.0)*_compSolver.colVector(_decompDualColIDs[i]), R(infinity), 2924 0.0); // adding a new column to the dual 2925 2926 SPxColId newDualCol; 2927 _compSolver.addCol(newDualCol, addEqualityCol); 2928 2929 // inserting the row and col ids for the added column. This is to be next to the original column that has 2930 // been duplicated. 2931 _decompPrimalRowIDs.insert(i + 1, 1, _decompPrimalRowIDs[i]); 2932 _decompDualColIDs.insert(i + 1, 1, newDualCol); 2933 assert(_realLP->number(_decompPrimalRowIDs[i]) == _realLP->number(_decompPrimalRowIDs[i + 1])); 2934 2935 i++; 2936 _nPrimalRows++; 2937 _nDualCols++; 2938 } 2939 } 2940 } 2941 } 2942 2943 2944 2945 /// identify the dual columns related to the fixed variables 2946 template <class R> 2947 void SoPlexBase<R>::_identifyComplementaryDualFixedPrimalVars(int* currFixedVars) 2948 { 2949 R feastol = realParam(SoPlexBase<R>::FEASTOL); 2950 2951 int numFixedVar = 0; 2952 2953 for(int i = 0; i < _realLP->nCols(); i++) 2954 { 2955 currFixedVars[i] = 0; 2956 2957 if(!_decompReducedProbColRowIDs[i].isValid()) 2958 continue; 2959 2960 int rowNumber = _solver.number(_decompReducedProbColRowIDs[i]); 2961 2962 if(_decompReducedProbColRowIDs[i].isValid()) 2963 { 2964 if(_solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::P_ON_UPPER || 2965 _solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::P_ON_LOWER || 2966 _solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::P_FIXED || 2967 _solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::D_FREE) 2968 { 2969 // setting the value of the _fixedOrigVars array to indicate which variables are at their bounds. 2970 currFixedVars[i] = getOrigVarFixedDirection(i); 2971 2972 numFixedVar++; 2973 } 2974 else 2975 { 2976 // the dual flags do not imply anything about the primal status of the rows. 2977 if(_solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::D_ON_LOWER && 2978 EQ(_solver.rhs(rowNumber) - _solver.pVec()[rowNumber], R(0.0), feastol)) 2979 currFixedVars[i] = 1; 2980 else if(_solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::D_ON_UPPER && 2981 EQ(_solver.pVec()[rowNumber] - _solver.lhs(rowNumber), R(0.0), feastol)) 2982 currFixedVars[i] = -1; 2983 } 2984 } 2985 } 2986 2987 SPX_MSG_INFO3(spxout, spxout << 2988 "Number of fixed primal variables in the complementary (dual) problem: " 2989 << numFixedVar << std::endl); 2990 } 2991 2992 2993 2994 /// removing the dual columns related to the fixed variables 2995 template <class R> 2996 void SoPlexBase<R>::_removeComplementaryDualFixedPrimalVars(int* currFixedVars) 2997 { 2998 SPxColId tempId; 2999 int ncolsforremoval = 0; 3000 int* colsforremoval = 0; 3001 spx_alloc(colsforremoval, _realLP->nCols() * 2); 3002 3003 tempId.inValidate(); 3004 3005 for(int i = 0; i < _realLP->nCols(); i++) 3006 { 3007 assert(_decompCompProbColIDsIdx[i] != -1); // this should be true in the current implementation 3008 3009 if(_decompCompProbColIDsIdx[i] != -1 3010 && _fixedOrigVars[i] != -2) //&& _fixedOrigVars[i] != currFixedVars[i] ) 3011 { 3012 if(_fixedOrigVars[i] != 0) 3013 { 3014 assert(_compSolver.number(SPxColId(_decompFixedVarDualIDs[i])) >= 0); 3015 assert(_fixedOrigVars[i] == -1 || _fixedOrigVars[i] == 1); 3016 assert(_decompFixedVarDualIDs[i].isValid()); 3017 3018 colsforremoval[ncolsforremoval] = _compSolver.number(SPxColId(_decompFixedVarDualIDs[i])); 3019 ncolsforremoval++; 3020 3021 _decompFixedVarDualIDs[i] = tempId; 3022 } 3023 else //if( false && !_decompReducedProbColRowIDs[i].isValid() ) // we want to remove all valid columns 3024 // in the current implementation, the only columns not included in the reduced problem are free columns. 3025 { 3026 assert((_realLP->lower(i) <= R(-infinity) && _realLP->upper(i) >= R(infinity)) || 3027 _compSolver.number(SPxColId(_decompVarBoundDualIDs[i * 2])) >= 0); 3028 int varcount = 0; 3029 3030 if(GT(_realLP->lower(i), R(-infinity), this->tolerances()->epsilon())) 3031 { 3032 colsforremoval[ncolsforremoval] = _compSolver.number(SPxColId(_decompVarBoundDualIDs[i * 2 + 3033 varcount])); 3034 ncolsforremoval++; 3035 3036 _decompVarBoundDualIDs[i * 2 + varcount] = tempId; 3037 varcount++; 3038 } 3039 3040 if(_realLP->upper(i) < R(infinity)) 3041 { 3042 colsforremoval[ncolsforremoval] = _compSolver.number(SPxColId(_decompVarBoundDualIDs[i * 2 + 3043 varcount])); 3044 ncolsforremoval++; 3045 3046 _decompVarBoundDualIDs[i * 2 + varcount] = tempId; 3047 } 3048 3049 } 3050 } 3051 } 3052 3053 int* perm = 0; 3054 spx_alloc(perm, _compSolver.nCols()); 3055 _compSolver.removeCols(colsforremoval, ncolsforremoval, perm); 3056 3057 // freeing allocated memory 3058 spx_free(perm); 3059 spx_free(colsforremoval); 3060 } 3061 3062 3063 3064 /// updating the dual columns related to the fixed primal variables. 3065 template <class R> 3066 void SoPlexBase<R>::_updateComplementaryDualFixedPrimalVars(int* currFixedVars) 3067 { 3068 DSVectorBase<R> col(1); 3069 LPColSetBase<R> boundConsCols; 3070 LPColSetBase<R> fixedVarsDualCols(_nPrimalCols); 3071 int numFixedVars = 0; 3072 // the solution to the reduced problem results in a number of variables at their bounds. If such variables exist 3073 // it is necessary to include a dual column to the complementary problem related to a variable fixing. This is 3074 // equivalent to the tight constraints being converted to equality constraints. 3075 int numBoundConsCols = 0; 3076 int* boundConsColsAdded = 0; 3077 spx_alloc(boundConsColsAdded, _realLP->nCols()); 3078 3079 // NOTE: this loop only goes over the primal columns that are included in the complementary problem, i.e. the 3080 // columns from the original problem. 3081 // 29.04.15 in the current implementation, all bound constraints are included in the reduced problem. So, all 3082 // variables (columns) are included in the reduced problem. 3083 for(int i = 0; i < _realLP->nCols(); i++) 3084 { 3085 boundConsColsAdded[i] = 0; 3086 assert(_decompCompProbColIDsIdx[i] != -1); 3087 3088 if(_decompCompProbColIDsIdx[i] != -1) 3089 { 3090 int idIndex = _decompCompProbColIDsIdx[i]; 3091 assert(_compSolver.number(SPxRowId(_decompDualRowIDs[idIndex])) >= 0); 3092 col.add(_compSolver.number(SPxRowId(_decompDualRowIDs[idIndex])), 1.0); 3093 3094 if(currFixedVars[i] != 0) 3095 { 3096 assert(currFixedVars[i] == -1 || currFixedVars[i] == 1); 3097 3098 assert(_realLP->lower(i) == _solver.lhs(_decompReducedProbColRowIDs[i])); 3099 assert(_realLP->upper(i) == _solver.rhs(_decompReducedProbColRowIDs[i])); 3100 R colObjCoeff = 0; 3101 3102 if(currFixedVars[i] == -1) 3103 colObjCoeff = _solver.lhs(_decompReducedProbColRowIDs[i]); 3104 else 3105 colObjCoeff = _solver.rhs(_decompReducedProbColRowIDs[i]); 3106 3107 fixedVarsDualCols.add(colObjCoeff, R(-infinity), col, R(infinity)); 3108 numFixedVars++; 3109 } 3110 // 09.02.15 I think that the else should only be entered if the column does not exist in the reduced 3111 // prob. I have tested by just leaving this as an else (without the if), but I think that this is wrong. 3112 //else if( !_decompReducedProbColRowIDs[i].isValid() ) 3113 3114 // NOTE: in the current implementation all columns, except the free columns, are included int the reduced 3115 // problem. There in no need to include the variable bounds in the complementary problem. 3116 else //if( false && _fixedOrigVars[i] == -2 ) 3117 { 3118 bool isRedProbCol = _decompReducedProbColRowIDs[i].isValid(); 3119 3120 // 29.04.15 in the current implementation only free variables are not included in the reduced problem 3121 if(_realLP->lower(i) > R(-infinity)) 3122 { 3123 if(!isRedProbCol) 3124 col.add(_compSolver.number(SPxRowId(_compSlackDualRowId)), -SOPLEX_SLACKCOEFF); 3125 3126 boundConsCols.add(_realLP->lower(i), R(-infinity), col, 0.0); 3127 3128 if(!isRedProbCol) 3129 col.remove(col.size() - 1); 3130 3131 boundConsColsAdded[i]++; 3132 numBoundConsCols++; 3133 } 3134 3135 if(_realLP->upper(i) < R(infinity)) 3136 { 3137 if(!isRedProbCol) 3138 col.add(_compSolver.number(SPxRowId(_compSlackDualRowId)), SOPLEX_SLACKCOEFF); 3139 3140 boundConsCols.add(_realLP->upper(i), 0.0, col, R(infinity)); 3141 3142 if(!isRedProbCol) 3143 col.remove(col.size() - 1); 3144 3145 boundConsColsAdded[i]++; 3146 numBoundConsCols++; 3147 } 3148 } 3149 3150 col.clear(); 3151 _fixedOrigVars[i] = currFixedVars[i]; 3152 } 3153 } 3154 3155 // adding the fixed var dual columns to the complementary problem 3156 SPxColId* addedcolids = 0; 3157 spx_alloc(addedcolids, numFixedVars); 3158 _compSolver.addCols(addedcolids, fixedVarsDualCols); 3159 3160 SPxColId tempId; 3161 int addedcolcount = 0; 3162 3163 tempId.inValidate(); 3164 3165 for(int i = 0; i < _realLP->nCols(); i++) 3166 { 3167 if(_fixedOrigVars[i] != 0) 3168 { 3169 assert(_fixedOrigVars[i] == -1 || _fixedOrigVars[i] == 1); 3170 _decompFixedVarDualIDs[i] = addedcolids[addedcolcount]; 3171 addedcolcount++; 3172 } 3173 else 3174 _decompFixedVarDualIDs[i] = tempId; 3175 } 3176 3177 // adding the bound cons dual columns to the complementary problem 3178 SPxColId* addedbndcolids = 0; 3179 spx_alloc(addedbndcolids, numBoundConsCols); 3180 _compSolver.addCols(addedbndcolids, boundConsCols); 3181 3182 addedcolcount = 0; 3183 3184 for(int i = 0; i < _realLP->nCols(); i++) 3185 { 3186 if(boundConsColsAdded[i] > 0) 3187 { 3188 for(int j = 0; j < boundConsColsAdded[i]; j++) 3189 { 3190 _decompVarBoundDualIDs[i * 2 + j] = addedbndcolids[addedcolcount]; 3191 addedcolcount++; 3192 } 3193 } 3194 3195 switch(boundConsColsAdded[i]) 3196 { 3197 case 0: 3198 _decompVarBoundDualIDs[i * 2] = tempId; 3199 3200 // FALLTHROUGH 3201 case 1: 3202 _decompVarBoundDualIDs[i * 2 + 1] = tempId; 3203 break; 3204 } 3205 } 3206 3207 // freeing allocated memory 3208 spx_free(addedbndcolids); 3209 spx_free(addedcolids); 3210 spx_free(boundConsColsAdded); 3211 } 3212 3213 3214 /// identify the dual columns related to the fixed variables 3215 template <class R> 3216 void SoPlexBase<R>::_identifyComplementaryPrimalFixedPrimalVars(int* currFixedVars) 3217 { 3218 int numFixedVar = 0; 3219 3220 for(int i = 0; i < _realLP->nCols(); i++) 3221 { 3222 currFixedVars[i] = 0; 3223 3224 if(!_decompReducedProbColRowIDs[i].isValid()) 3225 continue; 3226 3227 int rowNumber = _solver.number(_decompReducedProbColRowIDs[i]); 3228 3229 if(_decompReducedProbColRowIDs[i].isValid() && 3230 (_solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::P_ON_UPPER || 3231 _solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::P_ON_LOWER || 3232 _solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::P_FIXED)) 3233 { 3234 // setting the value of the _fixedOrigVars array to indicate which variables are at their bounds. 3235 currFixedVars[i] = getOrigVarFixedDirection(i); 3236 3237 numFixedVar++; 3238 } 3239 } 3240 3241 SPX_MSG_INFO3(spxout, spxout << 3242 "Number of fixed primal variables in the complementary (primal) problem: " 3243 << numFixedVar << std::endl); 3244 } 3245 3246 3247 3248 /// updating the dual columns related to the fixed primal variables. 3249 template <class R> 3250 void SoPlexBase<R>::_updateComplementaryPrimalFixedPrimalVars(int* currFixedVars) 3251 { 3252 // NOTE: this loop only goes over the primal columns that are included in the complementary problem, i.e. the 3253 // columns from the original problem. 3254 // 29.04.15 in the current implementation, all bound constraints are included in the reduced problem. So, all 3255 // variables (columns) are included in the reduced problem. 3256 for(int i = 0; i < _nCompPrimalCols; i++) 3257 { 3258 int colNumber = _compSolver.number(SPxColId(_decompCompPrimalColIDs[i])); 3259 3260 if(_fixedOrigVars[colNumber] != currFixedVars[colNumber]) 3261 { 3262 if(currFixedVars[colNumber] != 0) 3263 { 3264 assert(currFixedVars[colNumber] == -1 || currFixedVars[colNumber] == 1); 3265 3266 if(currFixedVars[colNumber] == -1) 3267 _compSolver.changeBounds(colNumber, _realLP->lower(SPxColId(_decompPrimalColIDs[i])), 3268 _realLP->lower(SPxColId(_decompPrimalColIDs[i]))); 3269 else 3270 _compSolver.changeBounds(colNumber, _realLP->upper(SPxColId(_decompPrimalColIDs[i])), 3271 _realLP->upper(SPxColId(_decompPrimalColIDs[i]))); 3272 } 3273 else 3274 { 3275 _compSolver.changeBounds(colNumber, R(-infinity), R(infinity)); 3276 } 3277 } 3278 3279 _fixedOrigVars[colNumber] = currFixedVars[colNumber]; 3280 } 3281 } 3282 3283 3284 3285 /// updating the complementary dual problem with the original objective function 3286 template <class R> 3287 void SoPlexBase<R>::_setComplementaryDualOriginalObjective() 3288 { 3289 for(int i = 0; i < _realLP->nCols(); i++) 3290 { 3291 assert(_decompCompProbColIDsIdx[i] != -1); // this should be true in the current implementation 3292 int idIndex = _decompCompProbColIDsIdx[i]; 3293 int compRowNumber = _compSolver.number(_decompDualRowIDs[idIndex]); 3294 3295 if(_decompCompProbColIDsIdx[i] != -1) 3296 { 3297 // In the dual conversion, when a variable has a non-standard bound it is converted to a free variable. 3298 if(LE(_realLP->lower(i), R(-infinity)) && GE(_realLP->upper(i), R(infinity))) 3299 { 3300 // unrestricted variable 3301 _compSolver.changeLhs(compRowNumber, _realLP->obj(i)); 3302 _compSolver.changeRhs(compRowNumber, _realLP->obj(i)); 3303 assert(LE(_compSolver.lhs(compRowNumber), _compSolver.rhs(compRowNumber))); 3304 } 3305 else if(LE(_realLP->lower(i), R(-infinity))) 3306 { 3307 // variable with a finite upper bound 3308 _compSolver.changeRhs(compRowNumber, _realLP->obj(i)); 3309 3310 if(isZero(_realLP->upper(i), this->tolerances()->epsilon())) 3311 _compSolver.changeLhs(compRowNumber, R(-infinity)); 3312 else 3313 _compSolver.changeLhs(compRowNumber, _realLP->obj(i)); 3314 } 3315 else if(GE(_realLP->upper(i), R(infinity))) 3316 { 3317 // variable with a finite lower bound 3318 _compSolver.changeLhs(compRowNumber, _realLP->obj(i)); 3319 3320 if(isZero(_realLP->upper(i), this->tolerances()->epsilon())) 3321 _compSolver.changeRhs(compRowNumber, R(infinity)); 3322 else 3323 _compSolver.changeRhs(compRowNumber, _realLP->obj(i)); 3324 } 3325 else if(NE(_realLP->lower(i), _realLP->upper(i))) 3326 { 3327 // variable with a finite upper and lower bound 3328 if(isZero(_realLP->upper(i), this->tolerances()->epsilon())) 3329 { 3330 _compSolver.changeLhs(compRowNumber, _realLP->obj(i)); 3331 _compSolver.changeRhs(compRowNumber, R(infinity)); 3332 } 3333 else if(isZero(_realLP->upper(i), this->tolerances->epsilon())) 3334 { 3335 _compSolver.changeLhs(compRowNumber, R(-infinity)); 3336 _compSolver.changeRhs(compRowNumber, _realLP->obj(i)); 3337 } 3338 else 3339 { 3340 _compSolver.changeLhs(compRowNumber, _realLP->obj(i)); 3341 _compSolver.changeRhs(compRowNumber, _realLP->obj(i)); 3342 } 3343 } 3344 else 3345 { 3346 // fixed variable 3347 _compSolver.changeLhs(compRowNumber, _realLP->obj(i)); 3348 _compSolver.changeRhs(compRowNumber, _realLP->obj(i)); 3349 } 3350 } 3351 } 3352 3353 // removing the complementary problem slack column dual row 3354 _compSolver.removeRow(_compSlackDualRowId); 3355 } 3356 3357 3358 3359 /// updating the complementary primal problem with the original objective function 3360 template <class R> 3361 void SoPlexBase<R>::_setComplementaryPrimalOriginalObjective() 3362 { 3363 // the comp solver has not removed any columns. Only the slack variables have been added. 3364 assert(_realLP->nCols() == _compSolver.nCols() - 1); 3365 3366 for(int i = 0; i < _realLP->nCols(); i++) 3367 { 3368 int colNumber = _realLP->number(_decompPrimalColIDs[i]); 3369 int compColNumber = _compSolver.number(_decompCompPrimalColIDs[i]); 3370 _compSolver.changeObj(compColNumber, _realLP->maxObj(colNumber)); 3371 } 3372 3373 // removing the complementary problem slack column dual row 3374 _compSolver.removeCol(_compSlackColId); 3375 } 3376 3377 3378 3379 /// determining which bound the primal variables will be fixed to. 3380 template <class R> 3381 int SoPlexBase<R>::getOrigVarFixedDirection(int colNum) 3382 { 3383 if(!_decompReducedProbColRowIDs[colNum].isValid()) 3384 return 0; 3385 3386 int rowNumber = _solver.number(_decompReducedProbColRowIDs[colNum]); 3387 3388 // setting the value of the _fixedOrigVars array to indicate which variables are at their bounds. 3389 if(_solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::P_ON_UPPER || 3390 _solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::P_FIXED || 3391 _solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::D_FREE) 3392 { 3393 assert(_solver.rhs(rowNumber) < R(infinity)); 3394 return 1; 3395 } 3396 else if(_solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::P_ON_LOWER) 3397 { 3398 assert(_solver.lhs(rowNumber) > R(-infinity)); 3399 return -1; 3400 } 3401 3402 return 0; 3403 } 3404 3405 3406 3407 // @todo update this function and related comments. It has only been hacked together. 3408 /// checks result of the solving process and solves again without preprocessing if necessary 3409 // @todo need to evaluate the solution to ensure that it is solved to optimality and then we are able to perform the 3410 // next steps in the algorithm. 3411 template <class R> 3412 void SoPlexBase<R>::_evaluateSolutionDecomp(SPxSolverBase<R>& solver, SLUFactor<R>& sluFactor, 3413 typename SPxSimplifier<R>::Result result) 3414 { 3415 typename SPxSolverBase<R>::Status solverStat = SPxSolverBase<R>::UNKNOWN; 3416 3417 if(result == SPxSimplifier<R>::INFEASIBLE) 3418 solverStat = SPxSolverBase<R>::INFEASIBLE; 3419 else if(result == SPxSimplifier<R>::DUAL_INFEASIBLE) 3420 solverStat = SPxSolverBase<R>::INForUNBD; 3421 else if(result == SPxSimplifier<R>::UNBOUNDED) 3422 solverStat = SPxSolverBase<R>::UNBOUNDED; 3423 else if(result == SPxSimplifier<R>::VANISHED) 3424 solverStat = SPxSolverBase<R>::OPTIMAL; 3425 else if(result == SPxSimplifier<R>::OKAY) 3426 solverStat = solver.status(); 3427 3428 // updating the status of SoPlexBase if the problem solved is the reduced problem. 3429 if(_currentProb == DECOMP_ORIG || _currentProb == DECOMP_RED) 3430 _status = solverStat; 3431 3432 // process result 3433 // coverity[switch_selector_expr_is_constant] 3434 switch(solverStat) 3435 { 3436 case SPxSolverBase<R>::OPTIMAL: 3437 if(!_isRealLPLoaded) 3438 { 3439 solver.changeObjOffset(realParam(SoPlexBase<R>::OBJ_OFFSET)); 3440 _decompResolveWithoutPreprocessing(solver, sluFactor, result); 3441 // Need to solve the complementary problem 3442 return; 3443 } 3444 else 3445 _hasBasis = true; 3446 3447 break; 3448 3449 case SPxSolverBase<R>::UNBOUNDED: 3450 case SPxSolverBase<R>::INFEASIBLE: 3451 case SPxSolverBase<R>::INForUNBD: 3452 3453 // in case of infeasibility or unboundedness, we currently can not unsimplify, but have to solve the original LP again 3454 if(!_isRealLPLoaded) 3455 { 3456 solver.changeObjOffset(realParam(SoPlexBase<R>::OBJ_OFFSET)); 3457 _decompSimplifyAndSolve(solver, sluFactor, false, false); 3458 return; 3459 } 3460 else 3461 _hasBasis = (solver.basis().status() > SPxBasisBase<R>::NO_PROBLEM); 3462 3463 break; 3464 3465 case SPxSolverBase<R>::ABORT_DECOMP: 3466 case SPxSolverBase<R>::ABORT_EXDECOMP: 3467 3468 // in the initialisation of the decomposition simplex, we want to keep the current basis. 3469 if(!_isRealLPLoaded) 3470 { 3471 solver.changeObjOffset(realParam(SoPlexBase<R>::OBJ_OFFSET)); 3472 _decompResolveWithoutPreprocessing(solver, sluFactor, result); 3473 //_decompSimplifyAndSolve(solver, sluFactor, false, false); 3474 return; 3475 } 3476 else 3477 _hasBasis = (solver.basis().status() > SPxBasisBase<R>::NO_PROBLEM); 3478 3479 break; 3480 3481 case SPxSolverBase<R>::SINGULAR: 3482 case SPxSolverBase<R>::ABORT_CYCLING: 3483 3484 // in case of infeasibility or unboundedness, we currently can not unsimplify, but have to solve the original LP again 3485 if(!_isRealLPLoaded) 3486 { 3487 solver.changeObjOffset(realParam(SoPlexBase<R>::OBJ_OFFSET)); 3488 _decompSimplifyAndSolve(solver, sluFactor, false, false); 3489 return; 3490 } 3491 3492 break; 3493 3494 3495 // else fallthrough 3496 case SPxSolverBase<R>::ABORT_TIME: 3497 case SPxSolverBase<R>::ABORT_ITER: 3498 case SPxSolverBase<R>::ABORT_VALUE: 3499 case SPxSolverBase<R>::REGULAR: 3500 case SPxSolverBase<R>::RUNNING: 3501 3502 // store regular basis if there is no simplifier and the original problem is not in the solver because of 3503 // scaling; non-optimal bases should currently not be unsimplified 3504 if(_simplifier == 0 && solver.basis().status() > SPxBasisBase<R>::NO_PROBLEM) 3505 { 3506 _basisStatusRows.reSize(_decompLP->nRows()); 3507 _basisStatusCols.reSize(_decompLP->nCols()); 3508 assert(_basisStatusRows.size() == solver.nRows()); 3509 assert(_basisStatusCols.size() == solver.nCols()); 3510 3511 solver.getBasis(_basisStatusRows.get_ptr(), _basisStatusCols.get_ptr()); 3512 _hasBasis = true; 3513 } 3514 else 3515 _hasBasis = false; 3516 3517 break; 3518 3519 default: 3520 _hasBasis = false; 3521 break; 3522 } 3523 } 3524 3525 3526 3527 3528 /// checks the dual feasibility of the current basis 3529 template <class R> 3530 bool SoPlexBase<R>::checkBasisDualFeasibility(VectorBase<R> feasVec) 3531 { 3532 assert(_solver.rep() == SPxSolverBase<R>::ROW); 3533 assert(_solver.spxSense() == SPxLPBase<R>::MAXIMIZE); 3534 3535 R feastol = realParam(SoPlexBase<R>::FEASTOL); 3536 3537 // iterating through all elements in the basis to determine whether the dual feasibility is satisfied. 3538 for(int i = 0; i < _solver.nCols(); i++) 3539 { 3540 if(_solver.basis().baseId(i).isSPxRowId()) // find the row id's for rows in the basis 3541 { 3542 int rownumber = _solver.number(_solver.basis().baseId(i)); 3543 3544 if(_solver.basis().desc().rowStatus(rownumber) != SPxBasisBase<R>::Desc::P_ON_UPPER && 3545 _solver.basis().desc().rowStatus(rownumber) != SPxBasisBase<R>::Desc::P_FIXED) 3546 { 3547 if(GT(feasVec[i], (R) 0, feastol)) 3548 return false; 3549 } 3550 3551 if(_solver.basis().desc().rowStatus(rownumber) != SPxBasisBase<R>::Desc::P_ON_LOWER && 3552 _solver.basis().desc().rowStatus(rownumber) != SPxBasisBase<R>::Desc::P_FIXED) 3553 { 3554 if(LT(feasVec[i], (R) 0, feastol)) 3555 return false; 3556 } 3557 3558 } 3559 else if(_solver.basis().baseId( 3560 i).isSPxColId()) // get the column id's for the columns in the basis 3561 { 3562 int colnumber = _solver.number(_solver.basis().baseId(i)); 3563 3564 if(_solver.basis().desc().colStatus(colnumber) != SPxBasisBase<R>::Desc::P_ON_UPPER && 3565 _solver.basis().desc().colStatus(colnumber) != SPxBasisBase<R>::Desc::P_FIXED) 3566 { 3567 if(GT(feasVec[i], (R) 0, feastol)) 3568 return false; 3569 } 3570 3571 if(_solver.basis().desc().colStatus(colnumber) != SPxBasisBase<R>::Desc::P_ON_LOWER && 3572 _solver.basis().desc().colStatus(colnumber) != SPxBasisBase<R>::Desc::P_FIXED) 3573 { 3574 if(LT(feasVec[i], (R) 0, feastol)) 3575 return false; 3576 } 3577 } 3578 } 3579 3580 return true; 3581 } 3582 3583 3584 3585 /// returns the expected sign of the dual variables for the reduced problem 3586 template <class R> 3587 typename SoPlexBase<R>::DualSign SoPlexBase<R>::getExpectedDualVariableSign(int rowNumber) 3588 { 3589 if(_solver.isRowBasic(rowNumber)) 3590 { 3591 if(_solver.basis().desc().rowStatus(rowNumber) != SPxBasisBase<R>::Desc::P_ON_UPPER && 3592 _solver.basis().desc().rowStatus(rowNumber) != SPxBasisBase<R>::Desc::P_FIXED) 3593 return SoPlexBase<R>::IS_NEG; 3594 3595 if(_solver.basis().desc().rowStatus(rowNumber) != SPxBasisBase<R>::Desc::P_ON_LOWER && 3596 _solver.basis().desc().rowStatus(rowNumber) != SPxBasisBase<R>::Desc::P_FIXED) 3597 return SoPlexBase<R>::IS_POS; 3598 } 3599 3600 return SoPlexBase<R>::IS_FREE; 3601 } 3602 3603 3604 3605 /// returns the expected sign of the dual variables for the original problem 3606 template <class R> 3607 typename SoPlexBase<R>::DualSign SoPlexBase<R>::getOrigProbDualVariableSign(int rowNumber) 3608 { 3609 if(_realLP->rowType(rowNumber) == LPRowBase<R>::LESS_EQUAL) 3610 return IS_POS; 3611 3612 if(_realLP->rowType(rowNumber) == LPRowBase<R>::GREATER_EQUAL) 3613 return IS_NEG; 3614 3615 if(_realLP->rowType(rowNumber) == LPRowBase<R>::EQUAL) 3616 return IS_FREE; 3617 3618 // this final statement needs to be checked. 3619 if(_realLP->rowType(rowNumber) == LPRowBase<R>::RANGE) 3620 return IS_FREE; 3621 3622 return IS_FREE; 3623 } 3624 3625 3626 /// print display line of flying table 3627 template <class R> 3628 void SoPlexBase<R>::printDecompDisplayLine(SPxSolverBase<R>& solver, 3629 const SPxOut::Verbosity origVerb, bool force, bool forceHead) 3630 { 3631 // setting the verbosity level 3632 const SPxOut::Verbosity currVerb = spxout.getVerbosity(); 3633 spxout.setVerbosity(origVerb); 3634 3635 int displayFreq = intParam(SoPlexBase<R>::DECOMP_DISPLAYFREQ); 3636 3637 SPX_MSG_INFO1(spxout, 3638 3639 if(forceHead || (_decompDisplayLine % (displayFreq * 30) == 0)) 3640 { 3641 spxout << "type | time | iters | red iter | alg iter | rows | cols | shift | value\n"; 3642 } 3643 if(force || (_decompDisplayLine % displayFreq == 0)) 3644 { 3645 Real currentTime = _statistics->solvingTime->time(); 3646 (solver.type() == SPxSolverBase<R>::LEAVE) ? spxout << " L |" : spxout << " E |"; 3647 spxout << std::fixed << std::setw(7) << std::setprecision(1) << currentTime << " |"; 3648 spxout << std::scientific << std::setprecision(2); 3649 spxout << std::setw(8) << _statistics->iterations << " | "; 3650 spxout << std::scientific << std::setprecision(2); 3651 spxout << std::setw(8) << _statistics->iterationsRedProb << " | "; 3652 spxout << std::scientific << std::setprecision(2); 3653 spxout << std::setw(8) << _statistics->callsReducedProb << " | "; 3654 spxout << std::scientific << std::setprecision(2); 3655 spxout << std::setw(8) << numIncludedRows << " | "; 3656 spxout << std::scientific << std::setprecision(2); 3657 spxout << std::setw(8) << solver.nCols() << " | " 3658 << solver.shift() << " | " 3659 << std::setprecision(8) << solver.value() + solver.objOffset() 3660 << std::endl; 3661 3662 } 3663 _decompDisplayLine++; 3664 ); 3665 3666 spxout.setVerbosity(currVerb); 3667 } 3668 3669 3670 3671 /// stores the problem statistics of the original problem 3672 template <class R> 3673 void SoPlexBase<R>::getOriginalProblemStatistics() 3674 { 3675 numProbRows = _realLP->nRows(); 3676 numProbCols = _realLP->nCols(); 3677 nNonzeros = _realLP->nNzos(); 3678 minAbsNonzero = _realLP->minAbsNzo(); 3679 maxAbsNonzero = _realLP->maxAbsNzo(); 3680 3681 origCountLower = 0; 3682 origCountUpper = 0; 3683 origCountBoxed = 0; 3684 origCountFreeCol = 0; 3685 3686 origCountLhs = 0; 3687 origCountRhs = 0; 3688 origCountRanged = 0; 3689 origCountFreeRow = 0; 3690 3691 for(int i = 0; i < _realLP->nCols(); i++) 3692 { 3693 bool hasLower = false; 3694 bool hasUpper = false; 3695 3696 if(_realLP->lower(i) > R(-infinity)) 3697 { 3698 origCountLower++; 3699 hasLower = true; 3700 } 3701 3702 if(_realLP->upper(i) < R(infinity)) 3703 { 3704 origCountUpper++; 3705 hasUpper = true; 3706 } 3707 3708 if(hasUpper && hasLower) 3709 { 3710 origCountBoxed++; 3711 origCountUpper--; 3712 origCountLower--; 3713 } 3714 3715 if(!hasUpper && !hasLower) 3716 origCountFreeCol++; 3717 } 3718 3719 for(int i = 0; i < _realLP->nRows(); i++) 3720 { 3721 bool hasRhs = false; 3722 bool hasLhs = false; 3723 3724 if(_realLP->lhs(i) > R(-infinity)) 3725 { 3726 origCountLhs++; 3727 hasLhs = true; 3728 } 3729 3730 if(_realLP->rhs(i) < R(infinity)) 3731 { 3732 origCountRhs++; 3733 hasRhs = true; 3734 } 3735 3736 if(hasRhs && hasLhs) 3737 { 3738 if(EQ(_realLP->rhs(i), _realLP->lhs(i), _solver.epsilon())) 3739 origCountEqual++; 3740 else 3741 origCountRanged++; 3742 3743 origCountLhs--; 3744 origCountRhs--; 3745 } 3746 3747 if(!hasRhs && !hasLhs) 3748 origCountFreeRow++; 3749 } 3750 } 3751 3752 template <class R> 3753 void SoPlexBase<R>::printOriginalProblemStatistics(std::ostream& os) 3754 { 3755 os << " Columns : " << numProbCols << "\n" 3756 << " boxed : " << origCountBoxed << "\n" 3757 << " lower bound : " << origCountLower << "\n" 3758 << " upper bound : " << origCountUpper << "\n" 3759 << " free : " << origCountFreeCol << "\n" 3760 << " Rows : " << numProbRows << "\n" 3761 << " equal : " << origCountEqual << "\n" 3762 << " ranged : " << origCountRanged << "\n" 3763 << " lhs : " << origCountLhs << "\n" 3764 << " rhs : " << origCountRhs << "\n" 3765 << " free : " << origCountFreeRow << "\n" 3766 << " Nonzeros : " << nNonzeros << "\n" 3767 << " per column : " << R(nNonzeros) / R(numProbCols) << "\n" 3768 << " per row : " << R(nNonzeros) / R(numProbRows) << "\n" 3769 << " sparsity : " << R(nNonzeros) / R(numProbCols) / R(numProbRows) << "\n" 3770 << " min. abs. value : " << R(minAbsNonzero) << "\n" 3771 << " max. abs. value : " << R(maxAbsNonzero) << "\n"; 3772 } 3773 3774 3775 3776 /// gets the coefficient of the slack variable in the primal complementary problem 3777 template <class R> 3778 R SoPlexBase<R>::getCompSlackVarCoeff(int primalRowNum) 3779 { 3780 int indDir = 1; 3781 3782 switch(_realLP->rowType(_decompPrimalRowIDs[primalRowNum])) 3783 { 3784 // NOTE: check the sign of the slackCoeff for the Range constraints. This will depend on the method of 3785 // dual conversion. 3786 case LPRowBase<R>::RANGE: 3787 assert((primalRowNum < _nPrimalRows - 1 3788 && _realLP->number(SPxColId(_decompPrimalRowIDs[primalRowNum])) == 3789 _realLP->number(SPxColId(_decompPrimalRowIDs[primalRowNum + 1]))) || 3790 (primalRowNum > 0 && _realLP->number(SPxColId(_decompPrimalRowIDs[primalRowNum - 1])) == 3791 _realLP->number(SPxColId(_decompPrimalRowIDs[primalRowNum])))); 3792 3793 // determine with primalRowNum and primalRowNum+1 or primalRowNum-1 and primalRowNum have the same row id. 3794 if(_realLP->number(SPxColId(_decompPrimalRowIDs[primalRowNum - 1])) == 3795 _realLP->number(SPxColId(_decompPrimalRowIDs[primalRowNum]))) 3796 indDir = -1; 3797 3798 if(_compSolver.obj(_compSolver.number(SPxColId(_decompDualColIDs[primalRowNum]))) < 3799 _compSolver.obj(_compSolver.number(SPxColId(_decompDualColIDs[primalRowNum + indDir])))) 3800 return -SOPLEX_SLACKCOEFF; 3801 else 3802 return SOPLEX_SLACKCOEFF; 3803 3804 break; 3805 3806 case LPRowBase<R>::GREATER_EQUAL: 3807 return -SOPLEX_SLACKCOEFF; 3808 break; 3809 3810 case LPRowBase<R>::EQUAL: 3811 assert((primalRowNum < _nPrimalRows - 1 3812 && _realLP->number(SPxColId(_decompPrimalRowIDs[primalRowNum])) == 3813 _realLP->number(SPxColId(_decompPrimalRowIDs[primalRowNum + 1]))) || 3814 (primalRowNum > 0 && _realLP->number(SPxColId(_decompPrimalRowIDs[primalRowNum - 1])) == 3815 _realLP->number(SPxColId(_decompPrimalRowIDs[primalRowNum])))); 3816 3817 // FALLTHROUGH 3818 case LPRowBase<R>::LESS_EQUAL: 3819 return SOPLEX_SLACKCOEFF; 3820 break; 3821 3822 default: 3823 throw SPxInternalCodeException("XDECOMPSL01 This should never happen."); 3824 } 3825 } 3826 3827 3828 3829 /// gets violation of bounds; returns true on success 3830 template <class R> 3831 bool SoPlexBase<R>::getDecompBoundViolation(R& maxviol, R& sumviol) 3832 { 3833 R feastol = realParam(SoPlexBase<R>::FEASTOL); 3834 3835 VectorBase<R>& primal = _solReal._primal; 3836 assert(primal.dim() == _realLP->nCols()); 3837 3838 _nDecompViolBounds = 0; 3839 3840 maxviol = 0.0; 3841 sumviol = 0.0; 3842 3843 bool isViol = false; 3844 bool isMaxViol = false; 3845 3846 for(int i = _realLP->nCols() - 1; i >= 0; i--) 3847 { 3848 R currviol = 0.0; 3849 3850 R viol = _realLP->lower(i) - primal[i]; 3851 3852 isViol = false; 3853 isMaxViol = false; 3854 3855 if(viol > 0.0) 3856 { 3857 sumviol += viol; 3858 3859 if(viol > maxviol) 3860 { 3861 maxviol = viol; 3862 isMaxViol = true; 3863 } 3864 3865 if(currviol < viol) 3866 currviol = viol; 3867 } 3868 3869 if(GT(viol, R(0.0), feastol)) 3870 isViol = true; 3871 3872 viol = primal[i] - _realLP->upper(i); 3873 3874 if(viol > 0.0) 3875 { 3876 sumviol += viol; 3877 3878 if(viol > maxviol) 3879 { 3880 maxviol = viol; 3881 isMaxViol = true; 3882 } 3883 3884 if(currviol < viol) 3885 currviol = viol; 3886 } 3887 3888 if(GT(viol, R(0.0), feastol)) 3889 isViol = true; 3890 3891 if(isViol) 3892 { 3893 // updating the violated bounds list 3894 if(isMaxViol) 3895 { 3896 _decompViolatedBounds[_nDecompViolBounds] = _decompViolatedBounds[0]; 3897 _decompViolatedBounds[0] = i; 3898 } 3899 else 3900 _decompViolatedBounds[_nDecompViolBounds] = i; 3901 3902 _nDecompViolBounds++; 3903 } 3904 } 3905 3906 return true; 3907 } 3908 3909 3910 3911 /// gets violation of constraints; returns true on success 3912 template <class R> 3913 bool SoPlexBase<R>::getDecompRowViolation(R& maxviol, R& sumviol) 3914 { 3915 R feastol = realParam(SoPlexBase<R>::FEASTOL); 3916 3917 VectorBase<R>& primal = _solReal._primal; 3918 assert(primal.dim() == _realLP->nCols()); 3919 3920 VectorBase<R> activity(_realLP->nRows()); 3921 _realLP->computePrimalActivity(primal, activity); 3922 3923 _nDecompViolRows = 0; 3924 3925 maxviol = 0.0; 3926 sumviol = 0.0; 3927 3928 bool isViol = false; 3929 bool isMaxViol = false; 3930 3931 for(int i = _realLP->nRows() - 1; i >= 0; i--) 3932 { 3933 R currviol = 0.0; 3934 3935 isViol = false; 3936 isMaxViol = false; 3937 3938 R viol = _realLP->lhs(i) - activity[i]; 3939 3940 if(viol > 0.0) 3941 { 3942 sumviol += viol; 3943 3944 if(viol > maxviol) 3945 { 3946 maxviol = viol; 3947 isMaxViol = true; 3948 } 3949 3950 if(viol > currviol) 3951 currviol = viol; 3952 } 3953 3954 if(GT(viol, R(0.0), feastol)) 3955 isViol = true; 3956 3957 viol = activity[i] - _realLP->rhs(i); 3958 3959 if(viol > 0.0) 3960 { 3961 sumviol += viol; 3962 3963 if(viol > maxviol) 3964 { 3965 maxviol = viol; 3966 isMaxViol = true; 3967 } 3968 3969 if(viol > currviol) 3970 currviol = viol; 3971 } 3972 3973 if(GT(viol, R(0.0), feastol)) 3974 isViol = true; 3975 3976 if(isViol) 3977 { 3978 // updating the violated rows list 3979 if(isMaxViol) 3980 { 3981 _decompViolatedRows[_nDecompViolRows] = _decompViolatedRows[0]; 3982 _decompViolatedRows[0] = i; 3983 } 3984 else 3985 _decompViolatedRows[_nDecompViolRows] = i; 3986 3987 _nDecompViolRows++; 3988 } 3989 } 3990 3991 return true; 3992 } 3993 3994 3995 3996 /// function call to terminate the decomposition simplex 3997 template <class R> 3998 bool SoPlexBase<R>::decompTerminate(R timeLimit) 3999 { 4000 R maxTime = timeLimit; 4001 4002 // check if a time limit is actually set 4003 if(maxTime < 0 || maxTime >= R(infinity)) 4004 return false; 4005 4006 Real currentTime = _statistics->solvingTime->time(); 4007 4008 if(currentTime >= maxTime) 4009 { 4010 SPX_MSG_INFO2(spxout, spxout << " --- timelimit (" << _solver.getMaxTime() 4011 << ") reached" << std::endl;) 4012 _solver.setSolverStatus(SPxSolverBase<R>::ABORT_TIME); 4013 return true; 4014 } 4015 4016 return false; 4017 } 4018 4019 4020 4021 /// function to retrieve the original problem row basis status from the reduced and complementary problems 4022 template <class R> 4023 void SoPlexBase<R>::getOriginalProblemBasisRowStatus(DataArray< int >& degenerateRowNums, 4024 DataArray< typename SPxSolverBase<R>::VarStatus >& degenerateRowStatus, int& nDegenerateRows, 4025 int& nNonBasicRows) 4026 { 4027 R feastol = realParam(SoPlexBase<R>::FEASTOL); 4028 int basicRow = 0; 4029 4030 // looping over all rows not contained in the complementary problem 4031 for(int i = 0; i < _nElimPrimalRows; i++) 4032 { 4033 int rowNumber = _realLP->number(_decompElimPrimalRowIDs[i]); 4034 4035 int solverRowNum = _solver.number(_decompReducedProbRowIDs[rowNumber]); 4036 assert(solverRowNum >= 0 && solverRowNum < _solver.nRows()); 4037 4038 if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_UPPER) /*|| 4039 (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_LOWER && 4040 EQ(_solver.rhs(solverRowNum) - _solver.pVec()[solverRowNum], 0.0, feastol)) )*/ 4041 { 4042 _basisStatusRows[rowNumber] = SPxSolverBase<R>::ON_UPPER; 4043 } 4044 else if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_LOWER) /*|| 4045 (_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_UPPER && 4046 EQ(_solver.pVec()[solverRowNum] - _solver.lhs(solverRowNum), 0.0, feastol)) )*/ 4047 { 4048 _basisStatusRows[rowNumber] = SPxSolverBase<R>::ON_LOWER; 4049 } 4050 else if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_FIXED) 4051 { 4052 _basisStatusRows[rowNumber] = SPxSolverBase<R>::FIXED; 4053 } 4054 else if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_FREE) 4055 { 4056 _basisStatusRows[rowNumber] = SPxSolverBase<R>::ZERO; 4057 } 4058 else 4059 { 4060 _basisStatusRows[rowNumber] = SPxSolverBase<R>::BASIC; 4061 basicRow++; 4062 } 4063 } 4064 4065 for(int i = 0; i < _nPrimalRows; i++) 4066 { 4067 int rowNumber = _realLP->number(_decompPrimalRowIDs[i]); 4068 4069 // this loop runs over all rows previously in the complementary problem. 4070 if(_decompReducedProbRows[rowNumber]) 4071 { 4072 int solverRowNum = _solver.number(_decompReducedProbRowIDs[rowNumber]); 4073 assert(solverRowNum >= 0 && solverRowNum < _solver.nRows()); 4074 4075 if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_UPPER) 4076 { 4077 _basisStatusRows[rowNumber] = SPxSolverBase<R>::ON_UPPER; 4078 } 4079 else if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_LOWER && 4080 EQ(_solver.rhs(solverRowNum) - _solver.pVec()[solverRowNum], 0.0, feastol)) 4081 { 4082 // if a row is non basic, but is at its bound then the row number and status is stored 4083 degenerateRowNums[nDegenerateRows] = rowNumber; 4084 degenerateRowStatus[nDegenerateRows] = SPxSolverBase<R>::ON_UPPER; 4085 nDegenerateRows++; 4086 } 4087 else if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_ON_LOWER) 4088 { 4089 _basisStatusRows[rowNumber] = SPxSolverBase<R>::ON_LOWER; 4090 } 4091 else if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::D_ON_UPPER && 4092 EQ(_solver.pVec()[solverRowNum] - _solver.lhs(solverRowNum), 0.0, feastol)) 4093 { 4094 // if a row is non basic, but is at its bound then the row number and status is stored 4095 degenerateRowNums[nDegenerateRows] = rowNumber; 4096 degenerateRowStatus[nDegenerateRows] = SPxSolverBase<R>::ON_LOWER; 4097 nDegenerateRows++; 4098 } 4099 else if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_FIXED) 4100 { 4101 _basisStatusRows[rowNumber] = SPxSolverBase<R>::FIXED; 4102 } 4103 else if(_solver.basis().desc().rowStatus(solverRowNum) == SPxBasisBase<R>::Desc::P_FREE) 4104 { 4105 _basisStatusRows[rowNumber] = SPxSolverBase<R>::ZERO; 4106 } 4107 else 4108 { 4109 _basisStatusRows[rowNumber] = SPxSolverBase<R>::BASIC; 4110 basicRow++; 4111 } 4112 } 4113 else 4114 { 4115 _basisStatusRows[rowNumber] = SPxSolverBase<R>::BASIC; 4116 basicRow++; 4117 4118 switch(_realLP->rowType(_decompPrimalRowIDs[i])) 4119 { 4120 case LPRowBase<R>::RANGE: 4121 //assert(_realLP->number(SPxColId(_decompPrimalRowIDs[i])) == 4122 //_realLP->number(SPxColId(_decompPrimalRowIDs[i+1]))); 4123 //assert(_compSolver.obj(_compSolver.number(SPxColId(_decompDualColIDs[i]))) != 4124 //_compSolver.obj(_compSolver.number(SPxColId(_decompDualColIDs[i + 1])))); 4125 4126 //// NOTE: This is not correct at present. Need to modify the inequalities. Look at the dual conversion 4127 //// function to get this correct. 4128 //if( _compSolver.obj(_compSolver.number(SPxColId(_decompDualColIDs[i]))) < 4129 //_compSolver.obj(_compSolver.number(SPxColId(_decompDualColIDs[i + 1])))) 4130 //{ 4131 //if( _compSolver.basis().desc().rowStatus(_compSolver.number(SPxColId(_decompDualColIDs[i]))) == 4132 //SPxBasisBase<R>::Desc::D_ON_LOWER ) 4133 //{ 4134 //_basisStatusRows[rowNumber] = SPxSolverBase<R>::ON_LOWER; 4135 //} 4136 //else if( _compSolver.basis().desc().rowStatus(_compSolver.number(SPxColId(_decompDualColIDs[i + 1]))) == 4137 //SPxBasisBase<R>::Desc::D_ON_UPPER ) 4138 //{ 4139 //_basisStatusRows[rowNumber] = SPxSolverBase<R>::ON_UPPER; 4140 //} 4141 //else 4142 //{ 4143 //_basisStatusRows[rowNumber] = SPxSolverBase<R>::BASIC; 4144 //basicRow++; 4145 //} 4146 //} 4147 //else 4148 //{ 4149 //if( _compSolver.basis().desc().rowStatus(_compSolver.number(SPxColId(_decompDualColIDs[i]))) == 4150 //SPxBasisBase<R>::Desc::D_ON_UPPER ) 4151 //{ 4152 //_basisStatusRows[rowNumber] = SPxSolverBase<R>::ON_UPPER; 4153 //} 4154 //else if( _compSolver.basis().desc().rowStatus(_compSolver.number(SPxColId(_decompDualColIDs[i + 1]))) == 4155 //SPxBasisBase<R>::Desc::D_ON_LOWER ) 4156 //{ 4157 //_basisStatusRows[rowNumber] = SPxSolverBase<R>::ON_LOWER; 4158 //} 4159 //else 4160 //{ 4161 //_basisStatusRows[rowNumber] = SPxSolverBase<R>::BASIC; 4162 //basicRow++; 4163 //} 4164 //} 4165 4166 i++; 4167 break; 4168 4169 case LPRowBase<R>::EQUAL: 4170 //assert(_realLP->number(SPxColId(_decompPrimalRowIDs[i])) == 4171 //_realLP->number(SPxColId(_decompPrimalRowIDs[i+1]))); 4172 4173 //_basisStatusRows[rowNumber] = SPxSolverBase<R>::FIXED; 4174 4175 i++; 4176 break; 4177 4178 case LPRowBase<R>::GREATER_EQUAL: 4179 //if( _compSolver.basis().desc().rowStatus(_compSolver.number(SPxColId(_decompDualColIDs[i]))) == 4180 //SPxBasisBase<R>::Desc::D_ON_LOWER ) 4181 //{ 4182 //_basisStatusRows[rowNumber] = SPxSolverBase<R>::ON_LOWER; 4183 //} 4184 //else 4185 //{ 4186 //_basisStatusRows[rowNumber] = SPxSolverBase<R>::BASIC; 4187 //basicRow++; 4188 //} 4189 4190 break; 4191 4192 case LPRowBase<R>::LESS_EQUAL: 4193 //if( _compSolver.basis().desc().rowStatus(_compSolver.number(SPxColId(_decompDualColIDs[i]))) == 4194 //SPxBasisBase<R>::Desc::D_ON_UPPER ) 4195 //{ 4196 //_basisStatusRows[rowNumber] = SPxSolverBase<R>::ON_UPPER; 4197 //} 4198 //else 4199 //{ 4200 //_basisStatusRows[rowNumber] = SPxSolverBase<R>::BASIC; 4201 //basicRow++; 4202 //} 4203 4204 break; 4205 4206 default: 4207 throw SPxInternalCodeException("XDECOMPSL01 This should never happen."); 4208 } 4209 } 4210 } 4211 4212 nNonBasicRows = _realLP->nRows() - basicRow - nDegenerateRows; 4213 SPX_MSG_INFO2(spxout, spxout << "Number of non-basic rows: " << nNonBasicRows << " (from " 4214 << _realLP->nRows() << ")" << std::endl); 4215 } 4216 4217 4218 /// function to retrieve the column status for the original problem basis from the reduced and complementary problems 4219 // all columns are currently in the reduced problem, so it is only necessary to retrieve the status from there. 4220 // However, a transformation of the columns has been made, so it is only possible to retrieve the status from the 4221 // variable bound constraints. 4222 template <class R> 4223 void SoPlexBase<R>::getOriginalProblemBasisColStatus(int& nNonBasicCols) 4224 { 4225 R feastol = realParam(SoPlexBase<R>::FEASTOL); 4226 int basicCol = 0; 4227 int numZeroDual = 0; 4228 4229 // looping over all variable bound constraints 4230 for(int i = 0; i < _realLP->nCols(); i++) 4231 { 4232 if(!_decompReducedProbColRowIDs[i].isValid()) 4233 continue; 4234 4235 int rowNumber = _solver.number(_decompReducedProbColRowIDs[i]); 4236 4237 if(_decompReducedProbColRowIDs[i].isValid()) 4238 { 4239 if(_solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::P_ON_UPPER) /*|| 4240 (_solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::D_ON_LOWER && 4241 EQ(_solver.rhs(rowNumber) - _solver.pVec()[rowNumber], 0.0, feastol)) )*/ 4242 { 4243 _basisStatusCols[i] = SPxSolverBase<R>::ON_UPPER; 4244 } 4245 else if(_solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::P_ON_LOWER) /*|| 4246 (_solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::D_ON_UPPER && 4247 EQ(_solver.pVec()[rowNumber] - _solver.lhs(rowNumber), 0.0, feastol)) )*/ 4248 { 4249 _basisStatusCols[i] = SPxSolverBase<R>::ON_LOWER; 4250 } 4251 else if(_solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::P_FIXED) 4252 { 4253 _basisStatusCols[i] = SPxSolverBase<R>::FIXED; 4254 } 4255 else if(_solver.basis().desc().rowStatus(rowNumber) == SPxBasisBase<R>::Desc::P_FREE) 4256 { 4257 _basisStatusCols[i] = SPxSolverBase<R>::ZERO; 4258 } 4259 else 4260 { 4261 _basisStatusCols[i] = SPxSolverBase<R>::BASIC; 4262 basicCol++; 4263 } 4264 } 4265 4266 if(EQ(_solver.fVec()[i], 0.0, feastol)) 4267 numZeroDual++; 4268 } 4269 4270 nNonBasicCols = _realLP->nCols() - basicCol; 4271 SPX_MSG_INFO2(spxout, spxout << "Number of non-basic columns: " 4272 << nNonBasicCols << " (from " << _realLP->nCols() << ")" << std::endl 4273 << "Number of zero dual columns: " << numZeroDual << " (from " << _realLP->nCols() << ")" << 4274 std::endl); 4275 } 4276 4277 4278 4279 /// function to build a basis for the original problem as given by the solution to the reduced problem 4280 template <class R> 4281 void SoPlexBase<R>::_writeOriginalProblemBasis(const char* filename, NameSet* rowNames, 4282 NameSet* colNames, bool cpxFormat) 4283 { 4284 int numrows = _realLP->nRows(); 4285 int numcols = _realLP->nCols(); 4286 int nNonBasicRows = 0; 4287 int nNonBasicCols = 0; 4288 int nDegenerateRows = 0; 4289 DataArray< int > degenerateRowNums; // array to store the orig prob row number of degenerate rows 4290 DataArray< typename SPxSolverBase<R>::VarStatus > 4291 degenerateRowStatus; // possible status for the degenerate rows in the orig prob 4292 4293 // setting the basis status rows and columns to size of the original problem 4294 _basisStatusRows.reSize(numrows); 4295 _basisStatusCols.reSize(numcols); 4296 4297 degenerateRowNums.reSize(numrows); 4298 degenerateRowStatus.reSize(numrows); 4299 4300 // setting the row and column status for the original problem 4301 getOriginalProblemBasisRowStatus(degenerateRowNums, degenerateRowStatus, nDegenerateRows, 4302 nNonBasicRows); 4303 getOriginalProblemBasisColStatus(nNonBasicCols); 4304 4305 // checking the consistency of the non-basic rows and columns for printing out the basis. 4306 // it is necessary to have numcols of non-basic rows and columns. 4307 // if there are not enought non-basic rows and columns, then the degenerate rows are set as non-basic. 4308 // all degenerate rows that are not changed to non-basic must be set to basic. 4309 assert(nDegenerateRows + nNonBasicRows + nNonBasicCols >= numcols); 4310 int degenRowsSetNonBasic = 0; 4311 4312 for(int i = 0; i < nDegenerateRows; i++) 4313 { 4314 if(nNonBasicRows + nNonBasicCols + degenRowsSetNonBasic < numcols) 4315 { 4316 _basisStatusRows[degenerateRowNums[i]] = degenerateRowStatus[i]; 4317 degenRowsSetNonBasic++; 4318 } 4319 else 4320 _basisStatusRows[degenerateRowNums[i]] = SPxSolverBase<R>::BASIC; 4321 } 4322 4323 // writing the basis file for the original problem 4324 bool wasRealLPLoaded = _isRealLPLoaded; 4325 _isRealLPLoaded = false; 4326 writeBasisFile(filename, rowNames, colNames, cpxFormat); 4327 _isRealLPLoaded = wasRealLPLoaded; 4328 } 4329 } // namespace soplex 4330