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 /**@file spxsolver.h 26 * @brief main LP solver class 27 */ 28 #ifndef _SPXSOLVER_H_ 29 #define _SPXSOLVER_H_ 30 31 #include <assert.h> 32 #include <iostream> 33 #include <iomanip> 34 #include <sstream> 35 36 #include "soplex/spxdefines.h" 37 #include "soplex/timer.h" 38 #include "soplex/timerfactory.h" 39 #include "soplex/spxlp.h" 40 #include "soplex/spxbasis.h" 41 #include "soplex/array.h" 42 #include "soplex/random.h" 43 #include "soplex/unitvector.h" 44 #include "soplex/updatevector.h" 45 #include "soplex/stablesum.h" 46 47 #include "soplex/spxlpbase.h" 48 49 #define SOPLEX_HYPERPRICINGTHRESHOLD 5000 /**< do (auto) hyper pricing only if problem size (cols+rows) is larger than SOPLEX_HYPERPRICINGTHRESHOLD */ 50 #define SOPLEX_HYPERPRICINGSIZE 100 /**< size of initial candidate list for hyper pricing */ 51 #define SOPLEX_SPARSITYFACTOR 0.6 /**< percentage of infeasibilities that is considered sparse */ 52 #define SOPLEX_DENSEROUNDS 5 /**< number of refactorizations until sparsity is tested again */ 53 #define SOPLEX_SPARSITY_TRADEOFF 0.8 /**< threshold to decide whether Ids or coIds are preferred to enter the basis; 54 // * coIds are more likely to enter if SOPLEX_SPARSITY_TRADEOFF is close to 0 55 // */ 56 #define SOPLEX_MAXNCLCKSKIPS 32 /**< maximum number of clock skips (iterations without time measuring) */ 57 #define SOPLEX_SAFETYFACTOR 1e-2 /**< the probability to skip the clock when the time limit has been reached */ 58 #define SOPLEX_NINITCALLS 200 /**< the number of clock updates in isTimelimitReached() before clock skipping starts */ 59 namespace soplex 60 { 61 template <class R> 62 class SPxPricer; 63 template <class R> 64 class SPxRatioTester; 65 template <class R> 66 class SPxStarter; 67 template <class R> 68 class SPxFastRT; 69 template <class R> 70 class SPxBoundFlippingRT; 71 72 /**@brief Sequential object-oriented SimPlex. 73 @ingroup Algo 74 75 SPxSolverBase is an LP solver class using the revised Simplex algorithm. It 76 provides two basis representations, namely a column basis and a row basis 77 (see #Representation). For both representations, a primal and 78 dual algorithm is available (see \ref Type). 79 80 In addition, SPxSolverBase can be customized with various respects: 81 - pricing algorithms using SPxPricer 82 - ratio test using class SPxRatioTester 83 - computation of a start basis using class SPxStarter 84 - preprocessing of the LP using class SPxSimplifier 85 - termination criteria by overriding 86 87 SPxSolverBase is derived from SPxLPBase<R> that is used to store the LP to be solved. 88 Hence, the LPs solved with SPxSolverBase have the general format 89 90 \f[ 91 \begin{array}{rl} 92 \hbox{max} & \mbox{maxObj}^T x \\ 93 \hbox{s.t.} & \mbox{lhs} \le Ax \le \mbox{rhs} \\ 94 & \mbox{low} \le x \le \mbox{up} 95 \end{array} 96 \f] 97 98 Also, SPxLPBase<R> provide all manipulation methods for the LP. They allow 99 SPxSolverBase to be used within cutting plane algorithms. 100 */ 101 102 template <class R> 103 class SPxSolverBase : public SPxLPBase<R>, protected SPxBasisBase<R> 104 { 105 friend SPxFastRT<R>; 106 friend SPxBoundFlippingRT<R>; 107 108 public: 109 110 //----------------------------- 111 /**@name Data Types */ 112 ///@{ 113 /// LP basis representation. 114 /** Solving LPs with the Simplex algorithm requires the definition of a 115 * \em basis. A basis can be defined as a set of column vectors or a 116 * set of row vectors building a nonsingular matrix. We will refer to 117 * the first case as the \em columnwise representation and the latter 118 * case will be called the \em rowwise representation. 119 * 120 * Type Representation determines the representation of SPxSolverBase, i.e. 121 * a columnwise (#COLUMN == 1) or rowwise (#ROW == -1) one. 122 */ 123 enum Representation 124 { 125 ROW = -1, ///< rowwise representation. 126 COLUMN = 1 ///< columnwise representation. 127 }; 128 129 /// Algorithmic type. 130 /** SPxSolverBase uses the reviesed Simplex algorithm to solve LPs. 131 * Mathematically, one distinguishes the \em primal from the 132 * \em dual algorihm. Algorithmically, these relate to the two 133 * types #ENTER or #LEAVE. How they relate, depends on the chosen 134 * basis representation. This is desribed by the following table: 135 * 136 * <TABLE> 137 * <TR><TD> </TD><TD>ENTER </TD><TD>LEAVE </TD></TR> 138 * <TR><TD>ROW </TD><TD>DUAL </TD><TD>PRIMAL</TD></TR> 139 * <TR><TD>COLUMN</TD><TD>PRIMAL</TD><TD>DUAL </TD></TR> 140 * </TABLE> 141 */ 142 enum Type 143 { 144 /// Entering Simplex. 145 /** The Simplex loop for the entering Simplex can be sketched 146 * as follows: 147 * - \em Pricing : Select a variable to #ENTER the basis. 148 * - \em Ratio-Test : Select variable to #LEAVE the 149 * basis such that the basis remains feasible. 150 * - Perform the basis update. 151 */ 152 ENTER = -1, 153 /// Leaving Simplex. 154 /** The Simplex loop for the leaving Simplex can be sketched 155 * as follows: 156 * - \em Pricing: Select a variable to #LEAVE the basis. 157 * - \em Ratio-Test: Select variable to #ENTER the 158 * basis such that the basis remains priced. 159 * - Perform the basis update. 160 */ 161 LEAVE = 1 162 }; 163 164 /// Pricing type. 165 /** In case of the #ENTER%ing Simplex algorithm, for performance 166 * reasons it may be advisable not to compute and maintain up to 167 * date vectors #pVec() and #test() and instead compute only some 168 * of its elements explicitely. This is controled by the #Pricing type. 169 */ 170 enum Pricing 171 { 172 /// Full pricing. 173 /** If #FULL pricing in selected for the #ENTER%ing Simplex, 174 * vectors #pVec() and #test() are kept up to date by 175 * SPxSolverBase. An SPxPricer only needs to select an Id such 176 * that the #test() or #coTest() value is < 0. 177 */ 178 FULL, 179 /// Partial pricing. 180 /** When #PARTIAL pricing in selected for the #ENTER%ing 181 * Simplex, vectors #pVec() and #test() are not set up and 182 * updated by SPxSolverBase. However, vectors #coPvec() and 183 * #coTest() are still kept up to date by SPxSolverBase. 184 * An SPxPricer object needs to compute the values for 185 * #pVec() and #test() itself in order to select an 186 * appropriate pivot with #test() < 0. Methods \ref computePvec(int) 187 * "computePvec(i)" and \ref computeTest(int) "computeTest(i)" 188 * will assist the used to do so. Note 189 * that it may be feasible for a pricer to return an Id with 190 * #test() > 0; such will be rejected by SPxSolverBase. 191 */ 192 PARTIAL 193 }; 194 195 /// Improved dual simplex status 196 /** The improved dual simplex requires a starting basis to perform the problem partitioning. This flag sets the 197 * status of the improved dual simplex to indicate whether the starting basis must be found or not. 198 */ 199 enum DecompStatus 200 { 201 /// Starting basis has not been found yet 202 FINDSTARTBASIS = 0, 203 /// Starting basis has been found and the simplex can be executed as normal 204 DONTFINDSTARTBASIS = 1 205 }; 206 207 enum VarStatus 208 { 209 ON_UPPER, ///< variable set to its upper bound. 210 ON_LOWER, ///< variable set to its lower bound. 211 FIXED, ///< variable fixed to identical bounds. 212 ZERO, ///< free variable fixed to zero. 213 BASIC, ///< variable is basic. 214 UNDEFINED ///< nothing known about basis status (possibly due to a singular basis in transformed problem) 215 }; 216 217 /**@todo In spxchange, change the status to 218 if (m_status > 0) m_status = REGULAR; 219 */ 220 enum Status 221 { 222 ERROR = -15, ///< an error occured. 223 NO_RATIOTESTER = -14, ///< No ratiotester loaded 224 NO_PRICER = -13, ///< No pricer loaded 225 NO_SOLVER = -12, ///< No linear solver loaded 226 NOT_INIT = -11, ///< not initialised error 227 ABORT_EXDECOMP = -10, ///< solve() aborted to exit decomposition simplex 228 ABORT_DECOMP = -9, ///< solve() aborted due to commence decomposition simplex 229 ABORT_CYCLING = -8, ///< solve() aborted due to detection of cycling. 230 ABORT_TIME = -7, ///< solve() aborted due to time limit. 231 ABORT_ITER = -6, ///< solve() aborted due to iteration limit. 232 ABORT_VALUE = -5, ///< solve() aborted due to objective limit. 233 SINGULAR = -4, ///< Basis is singular, numerical troubles? 234 NO_PROBLEM = -3, ///< No Problem has been loaded. 235 REGULAR = -2, ///< LP has a usable Basis (maybe LP is changed). 236 RUNNING = -1, ///< algorithm is running 237 UNKNOWN = 0, ///< nothing known on loaded problem. 238 OPTIMAL = 1, ///< LP has been solved to optimality. 239 UNBOUNDED = 2, ///< LP has been proven to be primal unbounded. 240 INFEASIBLE = 3, ///< LP has been proven to be primal infeasible. 241 INForUNBD = 4, ///< LP is primal infeasible or unbounded. 242 OPTIMAL_UNSCALED_VIOLATIONS = 5 ///< LP has beed solved to optimality but unscaled solution contains violations. 243 }; 244 245 /// objective for solution polishing 246 enum SolutionPolish 247 { 248 POLISH_OFF, ///< don't perform modifications on optimal basis 249 POLISH_INTEGRALITY, ///< maximize number of basic slack variables, i.e. more variables on bounds 250 POLISH_FRACTIONALITY ///< minimize number of basic slack variables, i.e. more variables in between bounds 251 }; 252 253 254 ///@} 255 256 private: 257 258 //----------------------------- 259 /**@name Private data */ 260 ///@{ 261 Type theType; ///< entering or leaving algortihm. 262 Pricing thePricing; ///< full or partial pricing. 263 Representation theRep; ///< row or column representation. 264 SolutionPolish polishObj; ///< objective of solution polishing 265 Timer* theTime; ///< time spent in last call to method solve() 266 Timer::TYPE timerType; ///< type of timer (user or wallclock) 267 Real theCumulativeTime; ///< cumulative time spent in all calls to method solve() 268 int maxIters; ///< maximum allowed iterations. 269 Real maxTime; ///< maximum allowed time. 270 int nClckSkipsLeft; ///< remaining number of times the clock can be safely skipped 271 long nCallsToTimelim; /// < the number of calls to the method isTimeLimitReached() 272 R objLimit; ///< objective value limit. 273 Status m_status; ///< status of algorithm. 274 275 R m_nonbasicValue; ///< nonbasic part of current objective value 276 bool m_nonbasicValueUpToDate; ///< true, if the stored objValue is up to date 277 278 R m_pricingViol; ///< maximal feasibility violation of current solution 279 bool m_pricingViolUpToDate; ///< true, if the stored violation is up to date 280 281 R 282 m_pricingViolCo; ///< maximal feasibility violation of current solution in coDim 283 bool m_pricingViolCoUpToDate; ///< true, if the stored violation in coDim is up to date 284 int m_numViol; ///< number of violations of current solution 285 286 R entertolscale; ///< factor to temporarily decrease the entering tolerance 287 R leavetolscale; ///< factor to temporarily decrease the leaving tolerance 288 R theShift; ///< sum of all shifts applied to any bound. 289 R lastShift; ///< for forcing feasibility. 290 int m_maxCycle; ///< maximum steps before cycling is detected. 291 int m_numCycle; ///< actual number of degenerate steps so far. 292 bool initialized; ///< true, if all vectors are setup. 293 294 SSVectorBase<R>* 295 solveVector2; ///< when 2 systems are to be solved at a time; typically for speepest edge weights 296 SSVectorBase<R>* 297 solveVector2rhs; ///< when 2 systems are to be solved at a time; typically for speepest edge weights 298 SSVectorBase<R>* 299 solveVector3; ///< when 3 systems are to be solved at a time; typically reserved for bound flipping ratio test (basic solution will be modified!) 300 SSVectorBase<R>* 301 solveVector3rhs; ///< when 3 systems are to be solved at a time; typically reserved for bound flipping ratio test (basic solution will be modified!) 302 SSVectorBase<R>* 303 coSolveVector2; ///< when 2 systems are to be solved at a time; typically for speepest edge weights 304 SSVectorBase<R>* 305 coSolveVector2rhs; ///< when 2 systems are to be solved at a time; typically for speepest edge weights 306 SSVectorBase<R>* 307 coSolveVector3; ///< when 3 systems are to be solved at a time; typically reserved for bound flipping ratio test (basic solution will be modified!) 308 SSVectorBase<R>* 309 coSolveVector3rhs; ///< when 3 systems are to be solved at a time; typically reserved for bound flipping ratio test (basic solution will be modified!) 310 311 bool freePricer; ///< true iff thepricer should be freed inside of object 312 bool freeRatioTester; ///< true iff theratiotester should be freed inside of object 313 bool freeStarter; ///< true iff thestarter should be freed inside of object 314 315 /* Store the index of a leaving variable if only an instable entering variable has been found. 316 instableLeave == true iff this instable basis change should be performed. 317 (see spxsolve.hpp and leave.hpp) */ 318 int instableLeaveNum; 319 bool instableLeave; 320 R instableLeaveVal; 321 322 /* Store the id of an entering row or column if only an instable pivot has been found. 323 instableEnter == true iff this instable basis change should be performed. 324 (see spxsolve.hpp and enter.hpp) */ 325 SPxId instableEnterId; 326 bool instableEnter; 327 R instableEnterVal; 328 329 bool 330 recomputedVectors; ///< flag to perform clean up step to reduce numerical errors only once 331 332 int displayLine; 333 int displayFreq; 334 R sparsePricingFactor; ///< enable sparse pricing when viols < factor * dim() 335 336 bool 337 getStartingDecompBasis; ///< flag to indicate whether the simplex is solved to get the starting improved dual simplex basis 338 bool computeDegeneracy; 339 int 340 degenCompIterOffset; ///< the number of iterations performed before the degeneracy level is computed 341 int 342 decompIterationLimit; ///< the maximum number of iterations before the decomposition simplex is aborted. 343 344 DataArray<VarStatus> oldBasisStatusRows; 345 ///< stored stable basis met before a simplex pivot (used to warm start the solver) 346 DataArray<VarStatus> oldBasisStatusCols; 347 ///< They don't have setters because only the internal simplex method is meant to fill them 348 349 bool solvingForBoosted; ///< is this solver involved in a higher precision solving scheme? 350 int storeBasisSimplexFreq; ///< number of simplex pivots -1 to perform before storing stable basis 351 352 bool 353 fullPerturbation; ///< whether to perturb the entire problem or just the bounds relevant for the current pivot 354 int 355 printBasisMetric; ///< printing the current basis metric in the log (-1: off, 0: condition estimate, 1: trace, 2: determinant, 3: condition) 356 357 ///@} 358 359 protected: 360 361 //----------------------------- 362 /**@name Protected data */ 363 ///@{ 364 Array < UnitVectorBase<R> > unitVecs; ///< array of unit vectors 365 const SVSetBase<R>* thevectors; ///< the LP vectors according to representation 366 const SVSetBase<R>* thecovectors; ///< the LP coVectors according to representation 367 368 VectorBase<R> primRhs; ///< rhs VectorBase<R> for computing the primal vector 369 UpdateVector<R> primVec; ///< primal vector 370 VectorBase<R> dualRhs; ///< rhs VectorBase<R> for computing the dual vector 371 UpdateVector<R> dualVec; ///< dual vector 372 UpdateVector<R> addVec; ///< storage for thePvec = &addVec 373 374 VectorBase<R> theURbound; ///< Upper Row Feasibility bound 375 VectorBase<R> theLRbound; ///< Lower Row Feasibility bound 376 VectorBase<R> theUCbound; ///< Upper Column Feasibility bound 377 VectorBase<R> theLCbound; ///< Lower Column Feasibility bound 378 379 /** In entering Simplex algorithm, the ratio test must ensure that all 380 * \em basic variables remain within their feasibility bounds. To give fast 381 * acces to them, the bounds of basic variables are copied into the 382 * following two vectors. 383 */ 384 VectorBase<R> theUBbound; ///< Upper Basic Feasibility bound 385 VectorBase<R> theLBbound; ///< Lower Basic Feasibility bound 386 387 /** The values of the rhs corresponding to the current basis.*/ 388 VectorBase<R>* theFrhs; 389 /** The values of all basis variables. */ 390 UpdateVector<R>* theFvec; 391 392 /* The Copricing rhs and VectorBase<R> */ 393 VectorBase<R>* theCoPrhs; 394 UpdateVector<R>* theCoPvec; 395 /** The pricing VectorBase<R> */ 396 UpdateVector<R>* thePvec; 397 398 UpdateVector<R>* theRPvec; ///< row pricing vector 399 UpdateVector<R>* theCPvec; ///< column pricing vector 400 401 // The following vectors serve for the virtualization of shift bounds 402 //@todo In prinziple this schould be references. 403 VectorBase<R>* theUbound; ///< Upper bound for vars 404 VectorBase<R>* theLbound; ///< Lower bound for vars 405 VectorBase<R>* theCoUbound; ///< Upper bound for covars 406 VectorBase<R>* theCoLbound; ///< Lower bound for covars 407 408 // The following vectors serve for the virtualization of testing vectors 409 VectorBase<R> theCoTest; 410 VectorBase<R> theTest; 411 412 DSVectorBase<R> primalRay; ///< stores primal ray in case of unboundedness 413 DSVectorBase<R> dualFarkas; ///< stores dual farkas proof in case of infeasibility 414 415 int leaveCount; ///< number of LEAVE iterations 416 int enterCount; ///< number of ENTER iterations 417 int primalCount; ///< number of primal iterations 418 int polishCount; ///< number of solution polishing iterations 419 420 int boundflips; ///< number of performed bound flips 421 int totalboundflips; ///< total number of bound flips 422 423 int enterCycles; ///< the number of degenerate steps during the entering algorithm 424 int leaveCycles; ///< the number of degenerate steps during the leaving algorithm 425 int enterDegenCand; ///< the number of degenerate candidates in the entering algorithm 426 int leaveDegenCand; ///< the number of degenerate candidates in the leaving algorithm 427 R primalDegenSum; ///< the sum of the primal degeneracy percentage 428 R dualDegenSum; ///< the sum of the dual degeneracy percentage 429 430 SPxPricer<R>* thepricer; 431 SPxRatioTester<R>* theratiotester; 432 SPxStarter<R>* thestarter; 433 434 R boundrange; ///< absolute range of all bounds in the problem 435 R siderange; ///< absolute range of all side in the problem 436 R objrange; ///< absolute range of all objective coefficients in the problem 437 ///@} 438 439 //----------------------------- 440 /**@name Precision */ 441 ///@{ 442 /// is the solution precise enough, or should we increase delta() ? 443 virtual bool precisionReached(R& newpricertol) const; 444 445 /// determine ranges of problem values for bounds, sides and objective to assess numerical difficulties 446 void calculateProblemRanges(); 447 ///@} 448 449 public: 450 451 /// The random number generator used throughout the whole computation. Its seed can be modified. 452 Random random; 453 454 /** For the leaving Simplex algorithm this vector contains the indices of infeasible basic variables; 455 * for the entering Simplex algorithm this vector contains the indices of infeasible slack variables. 456 */ 457 DIdxSet infeasibilities; 458 /**For the entering Simplex algorithm these vectors contains the indices of infeasible basic variables. 459 */ 460 DIdxSet infeasibilitiesCo; 461 462 /// store indices that were changed in the previous iteration and must be checked in hyper pricing 463 DIdxSet updateViols; 464 DIdxSet updateViolsCo; 465 466 /** Binary vectors to store whether basic indices are infeasible 467 * the i-th entry equals false, if the i-th basic variable is not infeasible 468 * the i-th entry equals true, if the i-th basic variable is infeasible 469 */ 470 DataArray<int> 471 isInfeasible; ///< 0: index not violated, 1: index violated, 2: index violated and among candidate list 472 DataArray<int> 473 isInfeasibleCo; ///< 0: index not violated, 1: index violated, 2: index violated and among candidate list 474 475 /// These values enable or disable sparse pricing 476 bool sparsePricingLeave; ///< true if sparsePricing is turned on in the leaving Simplex 477 bool sparsePricingEnter; ///< true if sparsePricing is turned on in the entering Simplex for slack variables 478 bool sparsePricingEnterCo; ///< true if sparsePricing is turned on in the entering Simplex 479 bool hyperPricingLeave; ///< true if hyper sparse pricing is turned on in the leaving Simplex 480 bool hyperPricingEnter; ///< true if hyper sparse pricing is turned on in the entering Simplex 481 482 int remainingRoundsLeave; ///< number of dense rounds/refactorizations until sparsePricing is enabled again 483 int remainingRoundsEnter; 484 int remainingRoundsEnterCo; 485 486 /// dual pricing norms 487 VectorBase<R> weights; ///< store dual norms 488 VectorBase<R> coWeights; ///< store dual norms 489 bool weightsAreSetup; ///< are the dual norms already set up? 490 491 492 Timer* multTimeSparse; ///< time spent in setupPupdate() exploiting sparsity 493 Timer* multTimeFull; ///< time spent in setupPupdate() ignoring sparsity 494 Timer* multTimeColwise; ///< time spent in setupPupdate(), columnwise multiplication 495 Timer* multTimeUnsetup; ///< time spent in setupPupdate() w/o sparsity information 496 int multSparseCalls; ///< number of products exploiting sparsity 497 int multFullCalls; ///< number of products ignoring sparsity 498 int multColwiseCalls; ///< number of products, columnwise multiplication 499 int multUnsetupCalls; ///< number of products w/o sparsity information 500 501 SPxOut* spxout; ///< message handler 502 503 DataArray<int> 504 integerVariables; ///< supplementary variable information, 0: continous variable, 1: integer variable 505 506 //----------------------------- 507 void setOutstream(SPxOut& newOutstream) 508 { 509 spxout = &newOutstream; 510 SPxLPBase<R>::spxout = &newOutstream; 511 } 512 513 /// set the _tolerances member variable 514 virtual void setTolerances(std::shared_ptr<Tolerances> newTolerances) 515 { 516 this->_tolerances = newTolerances; 517 // set tolerances for all the UpdateVectors 518 this->primVec.setTolerances(newTolerances); 519 this->dualVec.setTolerances(newTolerances); 520 this->addVec.setTolerances(newTolerances); 521 this->theFvec->setTolerances(newTolerances); 522 this->theCoPvec->setTolerances(newTolerances); 523 this->thePvec->setTolerances(newTolerances); 524 this->theRPvec->setTolerances(newTolerances); 525 this->theCPvec->setTolerances(newTolerances); 526 } 527 528 /// returns current tolerances 529 const std::shared_ptr<Tolerances>& tolerances() const 530 { 531 return this->_tolerances; 532 } 533 534 /// set refactor threshold for nonzeros in last factorized basis matrix compared to updated basis matrix 535 void setNonzeroFactor(R f) 536 { 537 SPxBasisBase<R>::nonzeroFactor = f; 538 } 539 540 /// set refactor threshold for fill-in in current factor update compared to fill-in in last factorization 541 void setFillFactor(R f) 542 { 543 SPxBasisBase<R>::fillFactor = f; 544 } 545 546 /// set refactor threshold for memory growth in current factor update compared to the last factorization 547 void setMemFactor(R f) 548 { 549 SPxBasisBase<R>::memFactor = f; 550 } 551 552 /**@name Access */ 553 ///@{ 554 /// return the version of SPxSolverBase as number like 123 for 1.2.3 555 int version() const 556 { 557 return SOPLEX_VERSION; 558 } 559 /// return the internal subversion of SPxSolverBase as number 560 int subversion() const 561 { 562 return SOPLEX_SUBVERSION; 563 } 564 /// return the current basis representation. 565 Representation rep() const 566 { 567 return theRep; 568 } 569 570 /// return current Type. 571 Type type() const 572 { 573 return theType; 574 } 575 576 /// return current Pricing. 577 Pricing pricing() const 578 { 579 return thePricing; 580 } 581 582 /// return current starter. 583 SPxStarter<R>* starter() const 584 { 585 return thestarter; 586 } 587 ///@} 588 589 //----------------------------- 590 /**@name Setup 591 * Before solving an LP with an instance of SPxSolverBase, 592 * the following steps must be performed: 593 * 594 * -# Load the LP by copying an external LP or reading it from an 595 * input stream. 596 * -# Setup the pricer to use by loading an \ref soplex::SPxPricer 597 * "SPxPricer" object (if not already done in a previous call). 598 * -# Setup the ratio test method to use by loading an 599 * \ref soplex::SPxRatioTester "SPxRatioTester" object 600 * (if not already done in a previous call). 601 * -# Setup the linear system solver to use by loading an 602 * \ref soplex::SLinSolver "SLinSolver" object 603 * (if not already done in a previous call). 604 * -# Optionally setup an start basis generation method by loading an 605 * \ref soplex::SPxStarter "SPxStarter" object. 606 * -# Optionally setup a start basis by loading a 607 * \ref soplex::SPxBasisBase<R>::Desc "SPxBasisBase<R>::Desc" object. 608 * -# Optionally switch to another basis 609 * \ref soplex::SPxSolverBase<R>::Representation "Representation" 610 * by calling method \ref soplex::SPxSolverBase<R>::setRep() "setRep()". 611 * -# Optionally switch to another algorithm 612 * \ref soplex::SPxSolverBase<R>::Type "Type" 613 * by calling method \ref soplex::SPxSolverBase<R>::setType() "setType()". 614 * 615 * Now the solver is ready for execution. If the loaded LP is to be solved 616 * again from scratch, this can be done with method 617 * \ref soplex::SPxSolverBase<R>::reLoad() "reLoad()". Finally, 618 * \ref soplex::SPxSolverBase<R>::clear() "clear()" removes the LP from the solver. 619 */ 620 ///@{ 621 /// read LP from input stream. 622 virtual bool read(std::istream& in, NameSet* rowNames = 0, 623 NameSet* colNames = 0, DIdxSet* intVars = 0); 624 625 /// copy LP. 626 virtual void loadLP(const SPxLPBase<R>& LP, bool initSlackBasis = true); 627 /// setup linear solver to use. If \p destroy is true, \p slusolver will be freed in destructor. 628 virtual void setBasisSolver(SLinSolver<R>* slu, const bool destroy = false); 629 /// setup pricer to use. If \p destroy is true, \p pricer will be freed in destructor. 630 virtual void setPricer(SPxPricer<R>* pricer, const bool destroy = false); 631 /// setup ratio-tester to use. If \p destroy is true, \p tester will be freed in destructor. 632 virtual void setTester(SPxRatioTester<R>* tester, const bool destroy = false); 633 /// setup starting basis generator to use. If \p destroy is true, \p starter will be freed in destructor. 634 virtual void setStarter(SPxStarter<R>* starter, const bool destroy = false); 635 /// set a start basis. 636 virtual void loadBasis(const typename SPxBasisBase<R>::Desc&); 637 638 /// initialize #ROW or #COLUMN representation. 639 void initRep(Representation p_rep); 640 /// switch to #ROW or #COLUMN representation if not already used. 641 void setRep(Representation p_rep); 642 /// set \ref soplex::SPxSolverBase<R>::LEAVE "LEAVE" or \ref soplex::SPxSolverBase<R>::ENTER "ENTER" algorithm. 643 void setType(Type tp); 644 /// set \ref soplex::SPxSolverBase<R>::FULL "FULL" or \ref soplex::SPxSolverBase<R>::PARTIAL "PARTIAL" pricing. 645 void setPricing(Pricing pr); 646 /// turn on or off the improved dual simplex. 647 void setDecompStatus(DecompStatus decomp_stat); 648 649 /// reload LP. 650 virtual void reLoad(); 651 652 /// clear all data in solver. 653 virtual void clear(); 654 655 /// unscales the LP and reloads the basis 656 void unscaleLPandReloadBasis(); 657 658 /// invalidates the basis, triggers refactorization 659 void invalidateBasis(); 660 661 /** Load basis from \p filename in MPS format. If \p rowNames and \p 662 * colNames are \c NULL, default names are used for the constraints and 663 * variables. 664 */ 665 virtual bool readBasisFile(const char* filename, 666 const NameSet* rowNames, const NameSet* colNames); 667 668 /** Write basis to \p filename in MPS format. If \p rowNames and \p 669 * colNames are \c NULL, default names are used for the constraints and 670 * variables. 671 */ 672 virtual bool writeBasisFile(const char* filename, 673 const NameSet* rowNames, const NameSet* colNames, const bool cpxFormat = false) const; 674 675 /** Write current LP, basis, and parameter settings. 676 * LP is written in MPS format to "\p filename".mps, basis is written in "\p filename".bas, and parameters 677 * are written to "\p filename".set. If \p rowNames and \p colNames are \c NULL, default names are used for 678 * the constraints and variables. 679 */ 680 virtual bool writeState(const char* filename, 681 const NameSet* rowNames = NULL, const NameSet* colNames = NULL, const bool cpxFormat = false) const; 682 683 ///@} 684 685 /**@name Solving LPs */ 686 ///@{ 687 /// solve loaded LP. 688 /** Solves the loaded LP by processing the Simplex iteration until 689 * the termination criteria is fullfilled (see #terminate()). 690 * The SPxStatus of the solver will indicate the reason for termination. 691 * @param interrupt can be set externally to interrupt the solve 692 * @param polish should solution polishing be considered 693 * 694 * @throw SPxStatusException if either no problem, solver, pricer 695 * or ratiotester loaded or if solve is still running when it shouldn't be 696 */ 697 virtual Status solve(volatile bool* interrupt = NULL, bool polish = true); 698 699 /** Identify primal basic variables that have zero reduced costs and 700 * try to pivot them out of the basis to make them tight. 701 * This is supposed to decrease the number of fractional variables 702 * when solving LP relaxations of (mixed) integer programs. 703 * The objective must not be modified during this procedure. 704 * @return true, if objective was modified (due to numerics) and resolving is necessary 705 */ 706 bool performSolutionPolishing(); 707 708 /// set objective of solution polishing (0: off, 1: max_basic_slack, 2: min_basic_slack) 709 void setSolutionPolishing(SolutionPolish _polishObj) 710 { 711 polishObj = _polishObj; 712 } 713 714 /// return objective of solution polishing 715 SolutionPolish getSolutionPolishing() 716 { 717 return polishObj; 718 } 719 720 /// Status of solution process. 721 Status status() const; 722 723 /// current objective value. 724 /**@return Objective value of the current solution vector 725 * (see #getPrimalSol()). 726 */ 727 virtual R value(); 728 729 // update nonbasic part of the objective value by the given amount 730 /**@return whether nonbasic part of objective is reliable 731 */ 732 bool updateNonbasicValue(R objChange); 733 734 // trigger a recomputation of the nonbasic part of the objective value 735 void forceRecompNonbasicValue() 736 { 737 m_nonbasicValue = 0.0; 738 m_nonbasicValueUpToDate = false; 739 } 740 741 /// get solution vector for primal variables. 742 /** This method returns the Status of the basis. 743 * If it is #REGULAR or better, 744 * the primal solution vector of the current basis will be copied 745 * to the argument \p vector. Hence, \p vector must be of dimension 746 * #nCols(). 747 * 748 * @throw SPxStatusException if not initialized 749 */ 750 virtual Status getPrimalSol(VectorBase<R>& vector) const; 751 752 /// get VectorBase<R> of slack variables. 753 /** This method returns the Status of the basis. 754 * If it is #REGULAR or better, 755 * the slack variables of the current basis will be copied 756 * to the argument \p vector. Hence, \p VectorBase<R> must be of dimension 757 * #nRows(). 758 * 759 * @warning Because SPxSolverBase supports range constraints as its 760 * default, slack variables are defined in a nonstandard way: 761 * Let \em x be the current solution vector and \em A the constraint 762 * matrix. Then the vector of slack variables is defined as 763 * \f$s = Ax\f$. 764 * 765 * @throw SPxStatusException if no problem loaded 766 */ 767 virtual Status getSlacks(VectorBase<R>& vector) const; 768 769 /// get current solution VectorBase<R> for dual variables. 770 /** This method returns the Status of the basis. 771 * If it is #REGULAR or better, 772 * the VectorBase<R> of dual variables of the current basis will be copied 773 * to the argument \p vector. Hence, \p VectorBase<R> must be of dimension 774 * #nRows(). 775 * 776 * @warning Even though mathematically, each range constraint would 777 * account for two dual variables (one for each inequaility), only 778 * #nRows() dual variables are setup via the following 779 * construction: Given a range constraint, there are three possible 780 * situations: 781 * - None of its inequalities is tight: The dual variables 782 * for both are 0. However, when shifting (see below) 783 * occurs, it may be set to a value other than 0, which 784 * models a perturbed objective vector. 785 * - Both of its inequalities are tight: In this case the 786 * range constraint models an equality and we adopt the 787 * standard definition. 788 * - One of its inequalities is tight while the other is not: 789 * In this case only the dual variable for the tight 790 * constraint is given with the standard definition, while 791 * the other constraint is implicitely set to 0. 792 * 793 * @throw SPxStatusException if no problem loaded 794 */ 795 virtual Status getDualSol(VectorBase<R>& vector) const; 796 797 /// get vector of reduced costs. 798 /** This method returns the Status of the basis. 799 * If it is #REGULAR or better, 800 * the vector of reduced costs of the current basis will be copied 801 * to the argument \p vector. Hence, \p vector must be of dimension 802 * #nCols(). 803 * 804 * Let \em d denote the vector of dual variables, as defined above, 805 * and \em A the LPs constraint matrix. Then the reduced cost vector 806 * \em r is defined as \f$r^T = c^T - d^TA\f$. 807 * 808 * @throw SPxStatusException if no problem loaded 809 */ 810 virtual Status getRedCostSol(VectorBase<R>& vector) const; 811 812 /// get primal ray in case of unboundedness. 813 /// @throw SPxStatusException if no problem loaded 814 virtual Status getPrimalray(VectorBase<R>& vector) const; 815 816 /// get dual farkas proof of infeasibility. 817 /// @throw SPxStatusException if no problem loaded 818 virtual Status getDualfarkas(VectorBase<R>& vector) const; 819 820 /// print display line of flying table 821 virtual void printDisplayLine(const bool force = false, const bool forceHead = false); 822 823 /// Termination criterion. 824 /** This method is called in each Simplex iteration to determine, if 825 * the algorithm is to terminate. In this case a nonzero value is 826 * returned. 827 * 828 * This method is declared virtual to allow for implementation of 829 * other stopping criteria or using it as callback method within the 830 * Simplex loop, by overriding the method in a derived class. 831 * However, all implementations must terminate with the 832 * statement \c return SPxSolverBase<R>::#terminate(), if no own termination 833 * criteria is encountered. 834 * 835 * Note, that the Simplex loop stopped even when #terminate() 836 * returns 0, if the LP has been solved to optimality (i.e. no 837 * further pricing succeeds and no shift is present). 838 */ 839 virtual bool terminate(); 840 ///@} 841 842 //----------------------------- 843 /**@name Control Parameters */ 844 ///@{ 845 /// values \f$|x| < \epsilon\f$ are considered to be 0. 846 /** if you want another value for epsilon, use 847 * \ref soplex::Tolerances::setEpsilon() "Tolerances::setEpsilon()". 848 */ 849 R epsilon() const 850 { 851 return this->tolerances()->epsilon(); 852 } 853 /// feasibility tolerance maintained by ratio test during ENTER algorithm. 854 R entertol() const 855 { 856 if(theRep == COLUMN) 857 return this->tolerances()->floatingPointFeastol() * this->entertolscale; 858 else 859 return this->tolerances()->floatingPointOpttol() * this->entertolscale; 860 } 861 /// feasibility tolerance maintained by ratio test during LEAVE algorithm. 862 R leavetol() const 863 { 864 if(theRep == COLUMN) 865 return this->tolerances()->floatingPointOpttol() * this->leavetolscale; 866 else 867 return this->tolerances()->floatingPointFeastol() * this->leavetolscale; 868 } 869 /// scale the entering tolerance 870 void scaleEntertol(R d) 871 { 872 this->entertolscale = d; 873 } 874 /// scale the leaving tolerance 875 void scaleLeavetol(R d) 876 { 877 this->leavetolscale = d; 878 } 879 void scaleTolerances(R d) 880 { 881 this->scaleEntertol(d); 882 this->scaleLeavetol(d); 883 } 884 /// guaranteed primal and dual bound violation for optimal solution, returning the maximum of floatingPointFeastol() and floatingPointOpttol(). 885 R delta() const 886 { 887 return SOPLEX_MAX(this->tolerances()->floatingPointFeastol(), 888 this->tolerances()->floatingPointOpttol()); 889 } 890 891 /// set timing type 892 void setTiming(Timer::TYPE ttype) 893 { 894 theTime = TimerFactory::switchTimer(theTime, ttype); 895 multTimeSparse = TimerFactory::switchTimer(multTimeSparse, ttype); 896 multTimeFull = TimerFactory::switchTimer(multTimeFull, ttype); 897 multTimeColwise = TimerFactory::switchTimer(multTimeColwise, ttype); 898 multTimeUnsetup = TimerFactory::switchTimer(multTimeUnsetup, ttype); 899 timerType = ttype; 900 } 901 902 /// set timing type 903 Timer::TYPE getTiming() 904 { 905 assert(timerType == theTime->type()); 906 assert(timerType == multTimeSparse->type()); 907 assert(timerType == multTimeFull->type()); 908 assert(timerType == multTimeColwise->type()); 909 assert(timerType == multTimeUnsetup->type()); 910 return timerType; 911 } 912 913 /// set display frequency 914 void setDisplayFreq(int freq) 915 { 916 displayFreq = freq; 917 } 918 919 /// get display frequency 920 int getDisplayFreq() 921 { 922 return displayFreq; 923 } 924 925 /// print basis metric within the usual output 926 void setMetricInformation(int type) 927 { 928 printBasisMetric = type; 929 } 930 931 // enable sparse pricing when viols < fac * dim() 932 void setSparsePricingFactor(R fac) 933 { 934 sparsePricingFactor = fac; 935 } 936 /// enable or disable hyper sparse pricing 937 void hyperPricing(bool h); 938 939 // get old basis status rows 940 DataArray<VarStatus>& getOldBasisStatusRows() 941 { 942 return oldBasisStatusRows; 943 } 944 945 // get old basis status cols 946 DataArray<VarStatus>& getOldBasisStatusCols() 947 { 948 return oldBasisStatusCols; 949 } 950 951 // should the basis be stored for use in precision boosting? 952 void setSolvingForBoosted(bool value) 953 { 954 solvingForBoosted = value; 955 } 956 957 // set frequency of storing the basis for use in precision boosting 958 void setStoreBasisFreqForBoosting(int freq) 959 { 960 storeBasisSimplexFreq = freq; 961 } 962 963 /** SPxSolverBase considers a Simplex step as degenerate if the 964 * steplength does not exceed #epsilon(). Cycling occurs if only 965 * degenerate steps are taken. To prevent this situation, SPxSolverBase 966 * perturbs the problem such that nondegenerate steps are ensured. 967 * 968 * maxCycle() controls how agressive such perturbation is 969 * performed, since no more than maxCycle() degenerate steps are 970 * accepted before perturbing the LP. The current number of consecutive 971 * degenerate steps is counted by numCycle(). 972 */ 973 /// maximum number of degenerate simplex steps before we detect cycling. 974 int maxCycle() const 975 { 976 return m_maxCycle; 977 } 978 /// actual number of degenerate simplex steps encountered so far. 979 int numCycle() const 980 { 981 return m_numCycle; 982 } 983 984 /// perturb entire problem or only the bounds relevant to the current pivot 985 void useFullPerturbation(bool full) 986 { 987 fullPerturbation = full; 988 } 989 990 virtual R getBasisMetric(int type) 991 { 992 return basis().getMatrixMetric(type); 993 } 994 995 ///@} 996 997 private: 998 999 //----------------------------- 1000 /**@name Private helpers */ 1001 ///@{ 1002 /// 1003 void localAddRows(int start); 1004 /// 1005 void localAddCols(int start); 1006 /// 1007 void setPrimal(VectorBase<R>& p_vector); 1008 /// 1009 void setSlacks(VectorBase<R>& p_vector); 1010 /// 1011 void setDual(VectorBase<R>& p_vector); 1012 /// 1013 void setRedCost(VectorBase<R>& p_vector); 1014 ///@} 1015 1016 protected: 1017 1018 //----------------------------- 1019 /**@name Protected helpers */ 1020 ///@{ 1021 /// 1022 virtual void addedRows(int n); 1023 /// 1024 virtual void addedCols(int n); 1025 /// 1026 virtual void doRemoveRow(int i); 1027 /// 1028 virtual void doRemoveRows(int perm[]); 1029 /// 1030 virtual void doRemoveCol(int i); 1031 /// 1032 virtual void doRemoveCols(int perm[]); 1033 ///@} 1034 1035 public: 1036 1037 //----------------------------- 1038 /**@name Modification */ 1039 /// \p scale determines whether the new data needs to be scaled according to the existing LP (persistent scaling) 1040 ///@{ 1041 /// 1042 virtual void changeObj(const VectorBase<R>& newObj, bool scale = false); 1043 /// 1044 virtual void changeObj(int i, const R& newVal, bool scale = false); 1045 /// 1046 using SPxLPBase<R>::changeObj; /// overloading a virtual function 1047 virtual void changeObj(SPxColId p_id, const R& p_newVal, bool scale = false) 1048 { 1049 changeObj(this->number(p_id), p_newVal, scale); 1050 } 1051 /// 1052 virtual void changeMaxObj(const VectorBase<R>& newObj, bool scale = false); 1053 /// 1054 virtual void changeMaxObj(int i, const R& newVal, bool scale = false); 1055 /// 1056 using SPxLPBase<R>::changeMaxObj; /// overloading a virtual function 1057 virtual void changeMaxObj(SPxColId p_id, const R& p_newVal, bool scale = false) 1058 { 1059 changeMaxObj(this->number(p_id), p_newVal, scale); 1060 } 1061 /// 1062 virtual void changeRowObj(const VectorBase<R>& newObj, bool scale = false); 1063 /// 1064 virtual void changeRowObj(int i, const R& newVal, bool scale = false); 1065 /// 1066 using SPxLPBase<R>::changeRowObj; 1067 virtual void changeRowObj(SPxRowId p_id, const R& p_newVal, bool scale = false) 1068 { 1069 changeRowObj(this->number(p_id), p_newVal); 1070 } 1071 /// 1072 virtual void clearRowObjs() 1073 { 1074 SPxLPBase<R>::clearRowObjs(); 1075 unInit(); 1076 } 1077 /// 1078 virtual void changeLowerStatus(int i, R newLower, R oldLower = 0.0); 1079 /// 1080 virtual void changeLower(const VectorBase<R>& newLower, bool scale = false); 1081 /// 1082 virtual void changeLower(int i, const R& newLower, bool scale = false); 1083 /// 1084 using SPxLPBase<R>::changeLower; 1085 virtual void changeLower(SPxColId p_id, const R& p_newLower, bool scale = false) 1086 { 1087 changeLower(this->number(p_id), p_newLower, scale); 1088 } 1089 /// 1090 virtual void changeUpperStatus(int i, R newUpper, R oldLower = 0.0); 1091 /// 1092 virtual void changeUpper(const VectorBase<R>& newUpper, bool scale = false); 1093 /// 1094 virtual void changeUpper(int i, const R& newUpper, bool scale = false); 1095 /// 1096 using SPxLPBase<R>::changeUpper; /// overloading virtual function 1097 virtual void changeUpper(SPxColId p_id, const R& p_newUpper, bool scale = false) 1098 { 1099 changeUpper(this->number(p_id), p_newUpper, scale); 1100 } 1101 /// 1102 virtual void changeBounds(const VectorBase<R>& newLower, const VectorBase<R>& newUpper, 1103 bool scale = false); 1104 /// 1105 virtual void changeBounds(int i, const R& newLower, const R& newUpper, bool scale = false); 1106 /// 1107 using SPxLPBase<R>::changeBounds; 1108 virtual void changeBounds(SPxColId p_id, const R& p_newLower, const R& p_newUpper, 1109 bool scale = false) 1110 { 1111 changeBounds(this->number(p_id), p_newLower, p_newUpper, scale); 1112 } 1113 /// 1114 virtual void changeLhsStatus(int i, R newLhs, R oldLhs = 0.0); 1115 /// 1116 virtual void changeLhs(const VectorBase<R>& newLhs, bool scale = false); 1117 /// 1118 virtual void changeLhs(int i, const R& newLhs, bool scale = false); 1119 /// 1120 using SPxLPBase<R>::changeLhs; 1121 virtual void changeLhs(SPxRowId p_id, const R& p_newLhs, bool scale = false) 1122 { 1123 changeLhs(this->number(p_id), p_newLhs, scale); 1124 } 1125 /// 1126 virtual void changeRhsStatus(int i, R newRhs, R oldRhs = 0.0); 1127 /// 1128 virtual void changeRhs(const VectorBase<R>& newRhs, bool scale = false); 1129 /// 1130 virtual void changeRhs(int i, const R& newRhs, bool scale = false); 1131 /// 1132 using SPxLPBase<R>::changeRhs; 1133 virtual void changeRhs(SPxRowId p_id, const R& p_newRhs, bool scale = false) 1134 { 1135 changeRhs(this->number(p_id), p_newRhs, scale); 1136 } 1137 /// 1138 virtual void changeRange(const VectorBase<R>& newLhs, const VectorBase<R>& newRhs, 1139 bool scale = false); 1140 /// 1141 virtual void changeRange(int i, const R& newLhs, const R& newRhs, bool scale = false); 1142 /// 1143 using SPxLPBase<R>::changeRange; 1144 virtual void changeRange(SPxRowId p_id, const R& p_newLhs, const R& p_newRhs, bool scale = false) 1145 { 1146 changeRange(this->number(p_id), p_newLhs, p_newRhs, scale); 1147 } 1148 /// 1149 virtual void changeRow(int i, const LPRowBase<R>& newRow, bool scale = false); 1150 /// 1151 using SPxLPBase<R>::changeRow; 1152 virtual void changeRow(SPxRowId p_id, const LPRowBase<R>& p_newRow, bool scale = false) 1153 { 1154 changeRow(this->number(p_id), p_newRow, scale); 1155 } 1156 /// 1157 virtual void changeCol(int i, const LPColBase<R>& newCol, bool scale = false); 1158 /// 1159 using SPxLPBase<R>::changeCol; 1160 virtual void changeCol(SPxColId p_id, const LPColBase<R>& p_newCol, bool scale = false) 1161 { 1162 changeCol(this->number(p_id), p_newCol, scale); 1163 } 1164 /// 1165 virtual void changeElement(int i, int j, const R& val, bool scale = false); 1166 /// 1167 using SPxLPBase<R>::changeElement; 1168 virtual void changeElement(SPxRowId rid, SPxColId cid, const R& val, bool scale = false) 1169 { 1170 changeElement(this->number(rid), this->number(cid), val, scale); 1171 } 1172 /// 1173 virtual void changeSense(typename SPxLPBase<R>::SPxSense sns); 1174 ///@} 1175 1176 //------------------------------------ 1177 /**@name Dimension and codimension */ 1178 ///@{ 1179 /// dimension of basis matrix. 1180 int dim() const 1181 { 1182 return thecovectors->num(); 1183 } 1184 /// codimension. 1185 int coDim() const 1186 { 1187 return thevectors->num(); 1188 } 1189 ///@} 1190 1191 //------------------------------------ 1192 /**@name Variables and Covariables 1193 * Class SPxLPBase<R> introduces \ref soplex::SPxId "SPxIds" to identify 1194 * row or column data of an LP. SPxSolverBase uses this concept to 1195 * access data with respect to the chosen representation. 1196 */ 1197 ///@{ 1198 /// id of \p i 'th vector. 1199 /** The \p i 'th Id is the \p i 'th SPxRowId for a rowwise and the 1200 * \p i 'th SPxColId for a columnwise basis represenation. Hence, 1201 * 0 <= i < #coDim(). 1202 */ 1203 SPxId id(int i) const 1204 { 1205 if(rep() == ROW) 1206 { 1207 SPxRowId rid = SPxLPBase<R>::rId(i); 1208 return SPxId(rid); 1209 } 1210 else 1211 { 1212 SPxColId cid = SPxLPBase<R>::cId(i); 1213 return SPxId(cid); 1214 } 1215 } 1216 1217 /// id of \p i 'th covector. 1218 /** The \p i 'th #coId() is the \p i 'th SPxColId for a rowwise and the 1219 * \p i 'th SPxRowId for a columnwise basis represenation. Hence, 1220 * 0 <= i < #dim(). 1221 */ 1222 SPxId coId(int i) const 1223 { 1224 if(rep() == ROW) 1225 { 1226 SPxColId cid = SPxLPBase<R>::cId(i); 1227 return SPxId(cid); 1228 } 1229 else 1230 { 1231 SPxRowId rid = SPxLPBase<R>::rId(i); 1232 return SPxId(rid); 1233 } 1234 } 1235 1236 /// Is \p p_id an SPxId ? 1237 /** This method returns wheather or not \p p_id identifies a vector 1238 * with respect to the chosen representation. 1239 */ 1240 bool isId(const SPxId& p_id) const 1241 { 1242 return p_id.info * theRep > 0; 1243 } 1244 1245 /// Is \p p_id a CoId. 1246 /** This method returns wheather or not \p p_id identifies a coVector 1247 * with respect to the chosen representation. 1248 */ 1249 bool isCoId(const SPxId& p_id) const 1250 { 1251 return p_id.info * theRep < 0; 1252 } 1253 ///@} 1254 1255 //------------------------------------ 1256 /**@name Vectors and Covectors */ 1257 ///@{ 1258 /// \p i 'th vector. 1259 /**@return a reference to the \p i 'th, 0 <= i < #coDim(), vector of 1260 * the loaded LP (with respect to the chosen representation). 1261 */ 1262 const SVectorBase<R>& vector(int i) const 1263 { 1264 return (*thevectors)[i]; 1265 } 1266 1267 /// 1268 const SVectorBase<R>& vector(const SPxRowId& rid) const 1269 { 1270 assert(rid.isValid()); 1271 return (rep() == ROW) 1272 ? (*thevectors)[this->number(rid)] 1273 : static_cast<const SVectorBase<R>&>(unitVecs[this->number(rid)]); 1274 } 1275 /// 1276 const SVectorBase<R>& vector(const SPxColId& cid) const 1277 { 1278 assert(cid.isValid()); 1279 return (rep() == COLUMN) 1280 ? (*thevectors)[this->number(cid)] 1281 : static_cast<const SVectorBase<R>&>(unitVecs[this->number(cid)]); 1282 } 1283 1284 /// VectorBase<R> associated to \p p_id. 1285 /**@return Returns a reference to the VectorBase<R> of the loaded LP corresponding 1286 * to \p id (with respect to the chosen representation). If \p p_id is 1287 * an id, a vector of the constraint matrix is returned, otherwise 1288 * the corresponding unit vector (of the slack variable or bound 1289 * inequality) is returned. 1290 * @todo The implementation does not exactly look like it will do 1291 * what is promised in the describtion. 1292 */ 1293 const SVectorBase<R>& vector(const SPxId& p_id) const 1294 { 1295 assert(p_id.isValid()); 1296 1297 return p_id.isSPxRowId() 1298 ? vector(SPxRowId(p_id)) 1299 : vector(SPxColId(p_id)); 1300 } 1301 1302 /// \p i 'th covector of LP. 1303 /**@return a reference to the \p i 'th, 0 <= i < #dim(), covector of 1304 * the loaded LP (with respect to the chosen representation). 1305 */ 1306 const SVectorBase<R>& coVector(int i) const 1307 { 1308 return (*thecovectors)[i]; 1309 } 1310 /// 1311 const SVectorBase<R>& coVector(const SPxRowId& rid) const 1312 { 1313 assert(rid.isValid()); 1314 return (rep() == COLUMN) 1315 ? (*thecovectors)[this->number(rid)] 1316 : static_cast<const SVector&>(unitVecs[this->number(rid)]); 1317 } 1318 /// 1319 const SVectorBase<R>& coVector(const SPxColId& cid) const 1320 { 1321 assert(cid.isValid()); 1322 return (rep() == ROW) 1323 ? (*thecovectors)[this->number(cid)] 1324 : static_cast<const SVectorBase<R>&>(unitVecs[this->number(cid)]); 1325 } 1326 /// coVector associated to \p p_id. 1327 /**@return a reference to the covector of the loaded LP 1328 * corresponding to \p p_id (with respect to the chosen 1329 * representation). If \p p_id is a coid, a covector of the constraint 1330 * matrix is returned, otherwise the corresponding unit vector is 1331 * returned. 1332 */ 1333 const SVectorBase<R>& coVector(const SPxId& p_id) const 1334 { 1335 assert(p_id.isValid()); 1336 return p_id.isSPxRowId() 1337 ? coVector(SPxRowId(p_id)) 1338 : coVector(SPxColId(p_id)); 1339 } 1340 /// return \p i 'th unit vector. 1341 const SVectorBase<R>& unitVector(int i) const 1342 { 1343 return unitVecs[i]; 1344 } 1345 ///@} 1346 1347 //------------------------------------ 1348 /**@name Variable status 1349 * The Simplex basis assigns a \ref soplex::SPxBasisBase<R>::Desc::Status 1350 * "Status" to each variable and covariable. Depending on the 1351 * representation, the status indicates that the corresponding 1352 * vector is in the basis matrix or not. 1353 */ 1354 ///@{ 1355 /// Status of \p i 'th variable. 1356 typename SPxBasisBase<R>::Desc::Status varStatus(int i) const 1357 { 1358 return this->desc().status(i); 1359 } 1360 1361 /// Status of \p i 'th covariable. 1362 typename SPxBasisBase<R>::Desc::Status covarStatus(int i) const 1363 { 1364 return this->desc().coStatus(i); 1365 } 1366 1367 /// does \p stat describe a basic index ? 1368 bool isBasic(typename SPxBasisBase<R>::Desc::Status stat) const 1369 { 1370 return (stat * rep() > 0); 1371 } 1372 1373 /// is the \p p_id 'th vector basic ? 1374 bool isBasic(const SPxId& p_id) const 1375 { 1376 assert(p_id.isValid()); 1377 return p_id.isSPxRowId() 1378 ? isBasic(SPxRowId(p_id)) 1379 : isBasic(SPxColId(p_id)); 1380 } 1381 1382 /// is the \p rid 'th vector basic ? 1383 bool isBasic(const SPxRowId& rid) const 1384 { 1385 return isBasic(this->desc().rowStatus(this->number(rid))); 1386 } 1387 1388 /// is the \p cid 'th vector basic ? 1389 bool isBasic(const SPxColId& cid) const 1390 { 1391 return isBasic(this->desc().colStatus(this->number(cid))); 1392 } 1393 1394 /// is the \p i 'th row vector basic ? 1395 bool isRowBasic(int i) const 1396 { 1397 return isBasic(this->desc().rowStatus(i)); 1398 } 1399 1400 /// is the \p i 'th column vector basic ? 1401 bool isColBasic(int i) const 1402 { 1403 return isBasic(this->desc().colStatus(i)); 1404 } 1405 1406 /// is the \p i 'th vector basic ? 1407 bool isBasic(int i) const 1408 { 1409 return isBasic(this->desc().status(i)); 1410 } 1411 1412 /// is the \p i 'th covector basic ? 1413 bool isCoBasic(int i) const 1414 { 1415 return isBasic(this->desc().coStatus(i)); 1416 } 1417 ///@} 1418 1419 /// feasibility vector. 1420 /** This method return the \em feasibility vector. If it satisfies its 1421 * bound, the basis is called feasible (independently of the chosen 1422 * representation). The feasibility vector has dimension #dim(). 1423 * 1424 * For the entering Simplex, #fVec is kept within its bounds. In 1425 * contrast to this, the pricing of the leaving Simplex selects an 1426 * element of #fVec, that violates its bounds. 1427 */ 1428 UpdateVector<R>& fVec() const 1429 { 1430 return *theFvec; 1431 } 1432 /// right-hand side vector for \ref soplex::SPxSolverBase<R>::fVec "fVec" 1433 /** The feasibility vector is computed by solving a linear system with the 1434 * basis matrix. The right-hand side vector of this system is referred 1435 * to as \em feasibility, \em right-hand \em side \em vector #fRhs(). 1436 * 1437 * For a row basis, #fRhs() is the objective vector (ignoring shifts). 1438 * For a column basis, it is the sum of all nonbasic vectors scaled by 1439 * the factor of their bound. 1440 */ 1441 const VectorBase<R>& fRhs() const 1442 { 1443 return *theFrhs; 1444 } 1445 /// upper bound for \ref soplex::SPxSolverBase<R>::fVec "fVec". 1446 const VectorBase<R>& ubBound() const 1447 { 1448 return theUBbound; 1449 } 1450 /// upper bound for #fVec, writable. 1451 /** This method returns the upper bound for the feasibility vector. 1452 * It may only be called for the #ENTER%ing Simplex. 1453 * 1454 * For the #ENTER%ing Simplex algorithms, the feasibility vector is 1455 * maintained to fullfill its bounds. As #fVec itself, also its 1456 * bounds depend on the chosen representation. Further, they may 1457 * need to be shifted (see below). 1458 */ 1459 VectorBase<R>& ubBound() 1460 { 1461 return theUBbound; 1462 } 1463 /// lower bound for \ref soplex::SPxSolverBase<R>::fVec "fVec". 1464 const VectorBase<R>& lbBound() const 1465 { 1466 return theLBbound; 1467 } 1468 /// lower bound for #fVec, writable. 1469 /** This method returns the lower bound for the feasibility vector. 1470 * It may only be called for the #ENTER%ing Simplex. 1471 * 1472 * For the #ENTER%ing Simplex algorithms, the feasibility vector is 1473 * maintained to fullfill its bounds. As #fVec itself, also its 1474 * bound depend on the chosen representation. Further, they may 1475 * need to be shifted (see below). 1476 */ 1477 VectorBase<R>& lbBound() 1478 { 1479 return theLBbound; 1480 } 1481 1482 /// Violations of \ref soplex::SPxSolverBase<R>::fVec "fVec" 1483 /** For the leaving Simplex algorithm, pricing involves selecting a 1484 * variable from #fVec that violates its bounds that is to leave 1485 * the basis. When a SPxPricer is called to select such a 1486 * leaving variable, #fTest() contains the vector of violations: 1487 * For #fTest()[i] < 0, the \c i 'th basic variable violates one of 1488 * its bounds by the given value. Otherwise no bound is violated. 1489 */ 1490 const VectorBase<R>& fTest() const 1491 { 1492 assert(type() == LEAVE); 1493 return theCoTest; 1494 } 1495 1496 /// copricing vector. 1497 /** The copricing vector #coPvec along with the pricing vector 1498 * #pVec are used for pricing in the #ENTER%ing Simplex algorithm, 1499 * i.e. one variable is selected, that violates its bounds. In 1500 * contrast to this, the #LEAVE%ing Simplex algorithm keeps both 1501 * vectors within their bounds. 1502 */ 1503 UpdateVector<R>& coPvec() const 1504 { 1505 return *theCoPvec; 1506 } 1507 1508 /// Right-hand side vector for \ref soplex::SPxSolverBase<R>::coPvec "coPvec". 1509 /** The vector #coPvec is computed by solving a linear system with the 1510 * basis matrix and #coPrhs as the right-hand side vector. For 1511 * column basis representation, #coPrhs is build up of the 1512 * objective vector elements of all basic variables. For a row 1513 * basis, it consists of the tight bounds of all basic 1514 * constraints. 1515 */ 1516 const VectorBase<R>& coPrhs() const 1517 { 1518 return *theCoPrhs; 1519 } 1520 1521 /// 1522 const VectorBase<R>& ucBound() const 1523 { 1524 assert(theType == LEAVE); 1525 return *theCoUbound; 1526 } 1527 /// upper bound for #coPvec. 1528 /** This method returns the upper bound for #coPvec. It may only be 1529 * called for the leaving Simplex algorithm. 1530 * 1531 * For the leaving Simplex algorithms #coPvec is maintained to 1532 * fullfill its bounds. As #coPvec itself, also its bound depend 1533 * on the chosen representation. Further, they may need to be 1534 * shifted (see below). 1535 */ 1536 VectorBase<R>& ucBound() 1537 { 1538 assert(theType == LEAVE); 1539 return *theCoUbound; 1540 } 1541 1542 /// 1543 const VectorBase<R>& lcBound() const 1544 { 1545 assert(theType == LEAVE); 1546 return *theCoLbound; 1547 } 1548 /// lower bound for #coPvec. 1549 /** This method returns the lower bound for #coPvec. It may only be 1550 * called for the leaving Simplex algorithm. 1551 * 1552 * For the leaving Simplex algorithms #coPvec is maintained to 1553 * fullfill its bounds. As #coPvec itself, also its bound depend 1554 * on the chosen representation. Further, they may need to be 1555 * shifted (see below). 1556 */ 1557 VectorBase<R>& lcBound() 1558 { 1559 assert(theType == LEAVE); 1560 return *theCoLbound; 1561 } 1562 1563 /// violations of \ref soplex::SPxSolverBase<R>::coPvec "coPvec". 1564 /** In entering Simplex pricing selects checks vectors #coPvec() 1565 * and #pVec() for violation of its bounds. #coTest() contains 1566 * the violations for #coPvec() which are indicated by a negative 1567 * value. That is, if #coTest()[i] < 0, the \p i 'th element of #coPvec() 1568 * is violated by -#coTest()[i]. 1569 */ 1570 const VectorBase<R>& coTest() const 1571 { 1572 assert(type() == ENTER); 1573 return theCoTest; 1574 } 1575 /// pricing vector. 1576 /** The pricing vector #pVec is the product of #coPvec with the 1577 * constraint matrix. As #coPvec, also #pVec is maintained within 1578 * its bound for the leaving Simplex algorithm, while the bounds 1579 * are tested for the entering Simplex. #pVec is of dimension 1580 * #coDim(). Vector #pVec() is only up to date for #LEAVE%ing 1581 * Simplex or #FULL pricing in #ENTER%ing Simplex. 1582 */ 1583 UpdateVector<R>& pVec() const 1584 { 1585 return *thePvec; 1586 } 1587 /// 1588 const VectorBase<R>& upBound() const 1589 { 1590 assert(theType == LEAVE); 1591 return *theUbound; 1592 } 1593 /// upper bound for #pVec. 1594 /** This method returns the upper bound for #pVec. It may only be 1595 * called for the leaving Simplex algorithm. 1596 * 1597 * For the leaving Simplex algorithms #pVec is maintained to 1598 * fullfill its bounds. As #pVec itself, also its bound depend 1599 * on the chosen representation. Further, they may need to be 1600 * shifted (see below). 1601 */ 1602 VectorBase<R>& upBound() 1603 { 1604 assert(theType == LEAVE); 1605 return *theUbound; 1606 } 1607 1608 /// 1609 const VectorBase<R>& lpBound() const 1610 { 1611 assert(theType == LEAVE); 1612 return *theLbound; 1613 } 1614 /// lower bound for #pVec. 1615 /** This method returns the lower bound for #pVec. It may only be 1616 * called for the leaving Simplex algorithm. 1617 * 1618 * For the leaving Simplex algorithms #pVec is maintained to 1619 * fullfill its bounds. As #pVec itself, also its bound depend 1620 * on the chosen representation. Further, they may need to be 1621 * shifted (see below). 1622 */ 1623 VectorBase<R>& lpBound() 1624 { 1625 assert(theType == LEAVE); 1626 return *theLbound; 1627 } 1628 1629 /// Violations of \ref soplex::SPxSolverBase<R>::pVec "pVec". 1630 /** In entering Simplex pricing selects checks vectors #coPvec() 1631 * and #pVec() for violation of its bounds. Vector #test() 1632 * contains the violations for #pVec(), i.e., if #test()[i] < 0, 1633 * the i'th element of #pVec() is violated by #test()[i]. 1634 * Vector #test() is only up to date for #FULL pricing. 1635 */ 1636 const VectorBase<R>& test() const 1637 { 1638 assert(type() == ENTER); 1639 return theTest; 1640 } 1641 1642 /// compute and return \ref soplex::SPxSolverBase<R>::pVec() "pVec()"[i]. 1643 R computePvec(int i); 1644 /// compute entire \ref soplex::SPxSolverBase<R>::pVec() "pVec()". 1645 void computePvec(); 1646 /// compute and return \ref soplex::SPxSolverBase<R>::test() "test()"[i] in \ref soplex::SPxSolverBase<R>::ENTER "ENTER"ing Simplex. 1647 R computeTest(int i); 1648 /// compute test VectorBase<R> in \ref soplex::SPxSolverBase<R>::ENTER "ENTER"ing Simplex. 1649 void computeTest(); 1650 1651 //------------------------------------ 1652 /**@name Shifting 1653 * The task of the ratio test (implemented in SPxRatioTester classes) 1654 * is to select a variable for the basis update, such that the basis 1655 * remains priced (i.e. both, the pricing and copricing vectors satisfy 1656 * their bounds) or feasible (i.e. the feasibility vector satisfies its 1657 * bounds). However, this can lead to numerically instable basis matrices 1658 * or -- after accumulation of various errors -- even to a singular basis 1659 * matrix. 1660 * 1661 * The key to overcome this problem is to allow the basis to become "a 1662 * bit" infeasible or unpriced, in order provide a better choice for the 1663 * ratio test to select a stable variable. This is equivalent to enlarging 1664 * the bounds by a small amount. This is referred to as \em shifting. 1665 * 1666 * These methods serve for shifting feasibility bounds, either in order 1667 * to maintain numerical stability or initially for computation of 1668 * phase 1. The sum of all shifts applied to any bound is stored in 1669 * \ref soplex::SPxSolverBase<R>::theShift "theShift". 1670 * 1671 * The following methods are used to shift individual bounds. They are 1672 * mainly intended for stable implenentations of SPxRatioTester. 1673 */ 1674 ///@{ 1675 /// Perform initial shifting to optain an feasible or pricable basis. 1676 void shiftFvec(); 1677 /// Perform initial shifting to optain an feasible or pricable basis. 1678 void shiftPvec(); 1679 1680 /// shift \p i 'th \ref soplex::SPxSolver::ubBound "ubBound" to \p to. 1681 void shiftUBbound(int i, R to) 1682 { 1683 assert(theType == ENTER); 1684 // use maximum to not count tightened bounds in case of equality shifts 1685 theShift += SOPLEX_MAX(to - theUBbound[i], 0.0); 1686 theUBbound[i] = to; 1687 } 1688 /// shift \p i 'th \ref soplex::SPxSolver::lbBound "lbBound" to \p to. 1689 void shiftLBbound(int i, R to) 1690 { 1691 assert(theType == ENTER); 1692 // use maximum to not count tightened bounds in case of equality shifts 1693 theShift += SOPLEX_MAX(theLBbound[i] - to, 0.0); 1694 theLBbound[i] = to; 1695 } 1696 /// shift \p i 'th \ref soplex::SPxSolver::upBound "upBound" to \p to. 1697 void shiftUPbound(int i, R to) 1698 { 1699 assert(theType == LEAVE); 1700 // use maximum to not count tightened bounds in case of equality shifts 1701 theShift += SOPLEX_MAX(to - (*theUbound)[i], 0.0); 1702 (*theUbound)[i] = to; 1703 } 1704 /// shift \p i 'th \ref soplex::SPxSolver::lpBound "lpBound" to \p to. 1705 void shiftLPbound(int i, R to) 1706 { 1707 assert(theType == LEAVE); 1708 // use maximum to not count tightened bounds in case of equality shifts 1709 theShift += SOPLEX_MAX((*theLbound)[i] - to, 0.0); 1710 (*theLbound)[i] = to; 1711 } 1712 /// shift \p i 'th \ref soplex::SPxSolver::ucBound "ucBound" to \p to. 1713 void shiftUCbound(int i, R to) 1714 { 1715 assert(theType == LEAVE); 1716 // use maximum to not count tightened bounds in case of equality shifts 1717 theShift += SOPLEX_MAX(to - (*theCoUbound)[i], 0.0); 1718 (*theCoUbound)[i] = to; 1719 } 1720 /// shift \p i 'th \ref soplex::SPxSolver::lcBound "lcBound" to \p to. 1721 void shiftLCbound(int i, R to) 1722 { 1723 assert(theType == LEAVE); 1724 // use maximum to not count tightened bounds in case of equality shifts 1725 theShift += SOPLEX_MAX((*theCoLbound)[i] - to, 0.0); 1726 (*theCoLbound)[i] = to; 1727 } 1728 /// 1729 void testBounds() const; 1730 1731 /// total current shift amount. 1732 virtual R shift() const 1733 { 1734 return theShift; 1735 } 1736 /// remove shift as much as possible. 1737 virtual void unShift(void); 1738 1739 /// get violation of constraints. 1740 virtual void qualConstraintViolation(R& maxviol, R& sumviol) const; 1741 /// get violations of bounds. 1742 virtual void qualBoundViolation(R& maxviol, R& sumviol) const; 1743 /// get the residuum |Ax-b|. 1744 virtual void qualSlackViolation(R& maxviol, R& sumviol) const; 1745 /// get violation of optimality criterion. 1746 virtual void qualRedCostViolation(R& maxviol, R& sumviol) const; 1747 ///@} 1748 1749 private: 1750 1751 //------------------------------------ 1752 /**@name Perturbation */ 1753 ///@{ 1754 /// 1755 void perturbMin( 1756 const UpdateVector<R>& vec, VectorBase<R>& low, VectorBase<R>& up, R eps, R delta, 1757 int start = 0, int incr = 1); 1758 /// 1759 void perturbMax( 1760 const UpdateVector<R>& vec, VectorBase<R>& low, VectorBase<R>& up, R eps, R delta, 1761 int start = 0, int incr = 1); 1762 /// 1763 R perturbMin(const UpdateVector<R>& uvec, 1764 VectorBase<R>& low, VectorBase<R>& up, R eps, R delta, 1765 const typename SPxBasisBase<R>::Desc::Status* stat, int start, int incr); 1766 /// 1767 R perturbMax(const UpdateVector<R>& uvec, 1768 VectorBase<R>& low, VectorBase<R>& up, R eps, R delta, 1769 const typename SPxBasisBase<R>::Desc::Status* stat, int start, int incr); 1770 ///@} 1771 1772 //------------------------------------ 1773 /**@name The Simplex Loop 1774 * We now present a set of methods that may be usefull when implementing 1775 * own SPxPricer or SPxRatioTester classes. Here is, how 1776 * SPxSolverBase will call methods from its loaded SPxPricer and 1777 * SPxRatioTester. 1778 * 1779 * For the entering Simplex: 1780 * -# \ref soplex::SPxPricer::selectEnter() "SPxPricer::selectEnter()" 1781 * -# \ref soplex::SPxRatioTester::selectLeave() "SPxRatioTester::selectLeave()" 1782 * -# \ref soplex::SPxPricer::entered4() "SPxPricer::entered4()" 1783 * 1784 * For the leaving Simplex: 1785 * -# \ref soplex::SPxPricer::selectLeave() "SPxPricer::selectLeave()" 1786 * -# \ref soplex::SPxRatioTester::selectEnter() "SPxRatioTester::selectEnter()" 1787 * -# \ref soplex::SPxPricer::left4() "SPxPricer::left4()" 1788 */ 1789 ///@{ 1790 public: 1791 /// Setup vectors to be solved within Simplex loop. 1792 /** Load vector \p y to be #solve%d with the basis matrix during the 1793 * #LEAVE Simplex. The system will be solved after #SPxSolverBase%'s call 1794 * to SPxRatioTester. The system will be solved along with 1795 * another system. Solving two linear system at a time has 1796 * performance advantages over solving the two linear systems 1797 * seperately. 1798 */ 1799 void setup4solve(SSVectorBase<R>* p_y, SSVectorBase<R>* p_rhs) 1800 { 1801 assert(type() == LEAVE); 1802 solveVector2 = p_y; 1803 solveVector2rhs = p_rhs; 1804 } 1805 /// Setup vectors to be solved within Simplex loop. 1806 /** Load a second additional vector \p y2 to be #solve%d with the 1807 * basis matrix during the #LEAVE Simplex. The system will be 1808 * solved after #SPxSolverBase%'s call to SPxRatioTester. 1809 * The system will be solved along with at least one 1810 * other system. Solving several linear system at a time has 1811 * performance advantages over solving them seperately. 1812 */ 1813 void setup4solve2(SSVectorBase<R>* p_y2, SSVectorBase<R>* p_rhs2) 1814 { 1815 assert(type() == LEAVE); 1816 solveVector3 = p_y2; 1817 solveVector3rhs = p_rhs2; 1818 } 1819 /// Setup vectors to be cosolved within Simplex loop. 1820 /** Load vector \p y to be #coSolve%d with the basis matrix during 1821 * the #ENTER Simplex. The system will be solved after #SPxSolverBase%'s 1822 * call to SPxRatioTester. The system will be solved along 1823 * with another system. Solving two linear system at a time has 1824 * performance advantages over solving the two linear systems 1825 * seperately. 1826 */ 1827 void setup4coSolve(SSVectorBase<R>* p_y, SSVectorBase<R>* p_rhs) 1828 { 1829 assert(type() == ENTER); 1830 coSolveVector2 = p_y; 1831 coSolveVector2rhs = p_rhs; 1832 } 1833 /// Setup vectors to be cosolved within Simplex loop. 1834 /** Load a second vector \p z to be #coSolve%d with the basis matrix during 1835 * the #ENTER Simplex. The system will be solved after #SPxSolverBase%'s 1836 * call to SPxRatioTester. The system will be solved along 1837 * with two other systems. 1838 */ 1839 void setup4coSolve2(SSVectorBase<R>* p_z, SSVectorBase<R>* p_rhs) 1840 { 1841 assert(type() == ENTER); 1842 coSolveVector3 = p_z; 1843 coSolveVector3rhs = p_rhs; 1844 } 1845 1846 /// maximal infeasibility of basis 1847 /** This method is called before concluding optimality. Since it is 1848 * possible that some stable implementation of class 1849 * SPxRatioTester yielded a slightly infeasible (or unpriced) 1850 * basis, this must be checked before terminating with an optimal 1851 * solution. 1852 */ 1853 virtual R maxInfeas() const; 1854 1855 /// check for violations above tol and immediately return false w/o checking the remaining values 1856 /** This method is useful for verifying whether an objective limit can be used as termination criterion 1857 */ 1858 virtual bool noViols(R tol) const; 1859 1860 /// Return current basis. 1861 /**@note The basis can be used to solve linear systems or use 1862 * any other of its (const) methods. It is, however, encuraged 1863 * to use methods #setup4solve() and #setup4coSolve() for solving 1864 * systems, since this is likely to have perfomance advantages. 1865 */ 1866 const SPxBasisBase<R>& basis() const 1867 { 1868 return *this; 1869 } 1870 /// 1871 SPxBasisBase<R>& basis() 1872 { 1873 return *this; 1874 } 1875 /// return loaded SPxPricer. 1876 const SPxPricer<R>* pricer() const 1877 { 1878 return thepricer; 1879 } 1880 /// return loaded SLinSolver. 1881 const SLinSolver<R>* slinSolver() const 1882 { 1883 return SPxBasisBase<R>::factor; 1884 } 1885 /// return loaded SPxRatioTester. 1886 const SPxRatioTester<R>* ratiotester() const 1887 { 1888 return theratiotester; 1889 } 1890 1891 /// Factorize basis matrix. 1892 /// @throw SPxStatusException if loaded matrix is singular 1893 virtual void factorize(); 1894 1895 private: 1896 1897 /** let index \p i leave the basis and manage entering of another one. 1898 @returns \c false if LP is unbounded/infeasible. */ 1899 bool leave(int i, bool polish = false); 1900 /** let id enter the basis and manage leaving of another one. 1901 @returns \c false if LP is unbounded/infeasible. */ 1902 bool enter(SPxId& id, bool polish = false); 1903 1904 /// test coVector \p i with status \p stat. 1905 R coTest(int i, typename SPxBasisBase<R>::Desc::Status stat) const; 1906 /// compute coTest vector. 1907 void computeCoTest(); 1908 /// recompute coTest vector. 1909 void updateCoTest(); 1910 1911 /// test VectorBase<R> \p i with status \p stat. 1912 R test(int i, typename SPxBasisBase<R>::Desc::Status stat) const; 1913 /// recompute test vector. 1914 void updateTest(); 1915 1916 /// compute basis feasibility test vector. 1917 void computeFtest(); 1918 /// update basis feasibility test vector. 1919 void updateFtest(); 1920 1921 ///@} 1922 1923 //------------------------------------ 1924 /**@name Parallelization 1925 * In this section we present the methods, that are provided in order to 1926 * allow a parallel version to be implemented as a derived class, thereby 1927 * inheriting most of the code of SPxSolverBase. 1928 * 1929 * @par Initialization 1930 * These methods are used to setup all the vectors used in the Simplex 1931 * loop, that where described in the previous sectios. 1932 */ 1933 ///@{ 1934 public: 1935 /// intialize data structures. 1936 /** If SPxSolverBase is not \ref isInitialized() "initialized", the method 1937 * #solve() calls #init() to setup all vectors and internal data structures. 1938 * Most of the other methods within this section are called by #init(). 1939 * 1940 * Derived classes should add the initialization of additional 1941 * data structures by overriding this method. Don't forget, 1942 * however, to call SPxSolverBase<R>::init(). 1943 */ 1944 virtual void init(); 1945 1946 protected: 1947 1948 /// has the internal data been initialized? 1949 /** As long as an instance of SPxSolverBase is not initialized, no member 1950 * contains setup data. Initialization is performed via method 1951 * #init(). Afterwards all data structures are kept up to date (even 1952 * for all manipulation methods), until #unInit() is called. However, 1953 * some manipulation methods call #unInit() themselfs. 1954 */ 1955 bool isInitialized() const 1956 { 1957 return initialized; 1958 } 1959 1960 /// resets clock average statistics 1961 void resetClockStats(); 1962 1963 /// uninitialize data structures. 1964 virtual void unInit() 1965 { 1966 initialized = false; 1967 } 1968 /// setup all vecs fresh 1969 virtual void reinitializeVecs(); 1970 /// reset dimensions of vectors according to loaded LP. 1971 virtual void reDim(); 1972 /// compute feasibility vector from scratch. 1973 void computeFrhs(); 1974 /// 1975 virtual void computeFrhsXtra(); 1976 /// 1977 virtual void computeFrhs1(const VectorBase<R>&, const VectorBase<R>&); 1978 /// 1979 void computeFrhs2(VectorBase<R>&, VectorBase<R>&); 1980 /// compute \ref soplex::SPxSolverBase<R>::theCoPrhs "theCoPrhs" for entering Simplex. 1981 virtual void computeEnterCoPrhs(); 1982 /// 1983 void computeEnterCoPrhs4Row(int i, int n); 1984 /// 1985 void computeEnterCoPrhs4Col(int i, int n); 1986 /// compute \ref soplex::SPxSolverBase<R>::theCoPrhs "theCoPrhs" for leaving Simplex. 1987 virtual void computeLeaveCoPrhs(); 1988 /// 1989 void computeLeaveCoPrhs4Row(int i, int n); 1990 /// 1991 void computeLeaveCoPrhs4Col(int i, int n); 1992 1993 /// Compute part of objective value. 1994 /** This method is called from #value() in order to compute the part of 1995 * the objective value resulting form nonbasic variables for #COLUMN 1996 * Representation. 1997 */ 1998 R nonbasicValue(); 1999 2000 /// Get pointer to the \p id 'th vector 2001 virtual const SVectorBase<R>* enterVector(const SPxId& p_id) 2002 { 2003 assert(p_id.isValid()); 2004 return p_id.isSPxRowId() 2005 ? &vector(SPxRowId(p_id)) : &vector(SPxColId(p_id)); 2006 } 2007 /// 2008 virtual void getLeaveVals(int i, 2009 typename SPxBasisBase<R>::Desc::Status& leaveStat, SPxId& leaveId, 2010 R& leaveMax, R& leavebound, int& leaveNum, StableSum<R>& objChange); 2011 /// 2012 virtual void getLeaveVals2(R leaveMax, SPxId enterId, 2013 R& enterBound, R& newUBbound, 2014 R& newLBbound, R& newCoPrhs, StableSum<R>& objChange); 2015 /// 2016 virtual void getEnterVals(SPxId id, R& enterTest, 2017 R& enterUB, R& enterLB, R& enterVal, R& enterMax, 2018 R& enterPric, typename SPxBasisBase<R>::Desc::Status& enterStat, R& enterRO, 2019 StableSum<R>& objChange); 2020 /// 2021 virtual void getEnterVals2(int leaveIdx, 2022 R enterMax, R& leaveBound, StableSum<R>& objChange); 2023 /// 2024 virtual void ungetEnterVal(SPxId enterId, typename SPxBasisBase<R>::Desc::Status enterStat, 2025 R leaveVal, const SVectorBase<R>& vec, StableSum<R>& objChange); 2026 /// 2027 virtual void rejectEnter(SPxId enterId, 2028 R enterTest, typename SPxBasisBase<R>::Desc::Status enterStat); 2029 /// 2030 virtual void rejectLeave(int leaveNum, SPxId leaveId, 2031 typename SPxBasisBase<R>::Desc::Status leaveStat, const SVectorBase<R>* newVec = 0); 2032 /// 2033 virtual void setupPupdate(void); 2034 /// 2035 virtual void doPupdate(void); 2036 /// 2037 virtual void clearUpdateVecs(void); 2038 /// 2039 virtual void perturbMinEnter(void); 2040 /// perturb basis bounds. 2041 virtual void perturbMaxEnter(void); 2042 /// 2043 virtual void perturbMinLeave(void); 2044 /// perturb nonbasic bounds. 2045 virtual void perturbMaxLeave(void); 2046 ///@} 2047 2048 //------------------------------------ 2049 /** The following methods serve for initializing the bounds for dual or 2050 * primal Simplex algorithm of entering or leaving type. 2051 */ 2052 ///@{ 2053 /// 2054 void clearDualBounds(typename SPxBasisBase<R>::Desc::Status, R&, R&) const; 2055 /// 2056 void setDualColBounds(); 2057 /// 2058 void setDualRowBounds(); 2059 /// setup feasibility bounds for entering algorithm 2060 void setPrimalBounds(); 2061 /// 2062 void setEnterBound4Col(int, int); 2063 /// 2064 void setEnterBound4Row(int, int); 2065 /// 2066 virtual void setEnterBounds(); 2067 /// 2068 void setLeaveBound4Row(int i, int n); 2069 /// 2070 void setLeaveBound4Col(int i, int n); 2071 /// 2072 virtual void setLeaveBounds(); 2073 ///@} 2074 2075 //------------------------------------ 2076 /** Compute the primal ray or the farkas proof in case of unboundedness 2077 * or infeasibility. 2078 */ 2079 ///@{ 2080 /// 2081 void computePrimalray4Col(R direction, SPxId enterId); 2082 /// 2083 void computePrimalray4Row(R direction); 2084 /// 2085 void computeDualfarkas4Col(R direction); 2086 /// 2087 void computeDualfarkas4Row(R direction, SPxId enterId); 2088 ///@} 2089 2090 public: 2091 2092 //------------------------------------ 2093 /** Limits and status inquiry */ 2094 ///@{ 2095 /// set time limit. 2096 virtual void setTerminationTime(Real time = infinity); 2097 /// return time limit. 2098 virtual Real terminationTime() const; 2099 /// set iteration limit. 2100 virtual void setTerminationIter(int iteration = -1); 2101 /// return iteration limit. 2102 virtual int terminationIter() const; 2103 /// set objective limit. 2104 virtual void setTerminationValue(R value = R(infinity)); 2105 /// return objective limit. 2106 virtual R terminationValue() const; 2107 /// get objective value of current solution. 2108 virtual R objValue() 2109 { 2110 return value(); 2111 } 2112 /// get all results of last solve. 2113 Status 2114 getResult(R* value = 0, VectorBase<R>* primal = 0, 2115 VectorBase<R>* slacks = 0, VectorBase<R>* dual = 0, 2116 VectorBase<R>* reduCost = 0); 2117 2118 protected: 2119 2120 /**@todo put the following basis methods near the variable status methods!*/ 2121 /// converts basis status to VarStatus 2122 VarStatus basisStatusToVarStatus(typename SPxBasisBase<R>::Desc::Status stat) const; 2123 2124 /// converts VarStatus to basis status for rows 2125 typename SPxBasisBase<R>::Desc::Status varStatusToBasisStatusRow(int row, VarStatus stat) 2126 const; 2127 2128 /// converts VarStatus to basis status for columns 2129 typename SPxBasisBase<R>::Desc::Status varStatusToBasisStatusCol(int col, VarStatus stat) 2130 const; 2131 2132 public: 2133 2134 /// gets basis status for a single row 2135 VarStatus getBasisRowStatus(int row) const; 2136 2137 /// gets basis status for a single column 2138 VarStatus getBasisColStatus(int col) const; 2139 2140 /// get current basis, and return solver status. 2141 Status getBasis(VarStatus rows[], VarStatus cols[], const int rowsSize = -1, 2142 const int colsSize = -1) const; 2143 2144 /// gets basis status 2145 typename SPxBasisBase<R>::SPxStatus getBasisStatus() const 2146 { 2147 return SPxBasisBase<R>::status(); 2148 } 2149 2150 /// check a given basis for validity. 2151 bool isBasisValid(DataArray<VarStatus> rows, DataArray<VarStatus> cols); 2152 2153 /// set the lp solver's basis. 2154 void setBasis(const VarStatus rows[], const VarStatus cols[]); 2155 2156 /// set the lp solver's basis status. 2157 void setBasisStatus(typename SPxBasisBase<R>::SPxStatus stat) 2158 { 2159 if(m_status == OPTIMAL) 2160 m_status = UNKNOWN; 2161 2162 SPxBasisBase<R>::setStatus(stat); 2163 } 2164 2165 /// setting the solver status external from the solve loop. 2166 void setSolverStatus(typename SPxSolverBase<R>::Status stat) 2167 { 2168 m_status = stat; 2169 } 2170 2171 /// get level of dual degeneracy 2172 // this function is used for the improved dual simplex 2173 R getDegeneracyLevel(VectorBase<R> degenvec); 2174 2175 /// get number of dual norms 2176 void getNdualNorms(int& nnormsRow, int& nnormsCol) const; 2177 2178 /// get dual norms 2179 bool getDualNorms(int& nnormsRow, int& nnormsCol, R* norms) const; 2180 2181 /// set dual norms 2182 bool setDualNorms(int nnormsRow, int nnormsCol, R* norms); 2183 2184 /// pass integrality information about the variables to the solver 2185 void setIntegralityInformation(int ncols, int* intInfo); 2186 2187 /// reset cumulative time counter to zero. 2188 void resetCumulativeTime() 2189 { 2190 theCumulativeTime = 0.0; 2191 } 2192 2193 /// get number of bound flips. 2194 int boundFlips() const 2195 { 2196 return totalboundflips; 2197 } 2198 2199 /// get number of dual degenerate pivots 2200 int dualDegeneratePivots() 2201 { 2202 return (rep() == ROW) ? enterCycles : leaveCycles; 2203 } 2204 2205 /// get number of primal degenerate pivots 2206 int primalDegeneratePivots() 2207 { 2208 return (rep() == ROW) ? leaveCycles : enterCycles; 2209 } 2210 2211 /// get the sum of dual degeneracy 2212 R sumDualDegeneracy() 2213 { 2214 return dualDegenSum; 2215 } 2216 2217 /// get the sum of primal degeneracy 2218 R sumPrimalDegeneracy() 2219 { 2220 return primalDegenSum; 2221 } 2222 2223 /// get number of iterations of current solution. 2224 int iterations() const 2225 { 2226 return basis().iteration(); 2227 } 2228 2229 /// return number of iterations done with primal algorithm 2230 int primalIterations() 2231 { 2232 assert(iterations() == 0 || primalCount <= iterations()); 2233 return (iterations() == 0) ? 0 : primalCount; 2234 } 2235 2236 /// return number of iterations done with primal algorithm 2237 int dualIterations() 2238 { 2239 return iterations() - primalIterations(); 2240 } 2241 2242 /// return number of iterations done with primal algorithm 2243 int polishIterations() 2244 { 2245 return polishCount; 2246 } 2247 2248 /// time spent in last call to method solve(). 2249 Real time() const 2250 { 2251 return theTime->time(); 2252 } 2253 2254 /// returns whether current time limit is reached; call to time() may be skipped unless \p forceCheck is true 2255 /// 2256 bool isTimeLimitReached(const bool forceCheck = false); 2257 2258 /// the maximum runtime 2259 Real getMaxTime() 2260 { 2261 return maxTime; 2262 } 2263 2264 /// cumulative time spent in all calls to method solve(). 2265 Real cumulativeTime() const 2266 { 2267 return theCumulativeTime; 2268 } 2269 2270 /// the maximum number of iterations 2271 int getMaxIters() 2272 { 2273 return maxIters; 2274 } 2275 2276 /// return const lp's rows if available. 2277 const LPRowSetBase<R>& rows() const 2278 { 2279 return *this->lprowset(); 2280 } 2281 2282 /// return const lp's cols if available. 2283 const LPColSet& cols() const 2284 { 2285 return *this->lpcolset(); 2286 } 2287 2288 /// copy lower bound VectorBase<R> to \p p_low. 2289 void getLower(VectorBase<R>& p_low) const 2290 { 2291 p_low = SPxLPBase<R>::lower(); 2292 } 2293 /// copy upper bound VectorBase<R> to \p p_up. 2294 void getUpper(VectorBase<R>& p_up) const 2295 { 2296 p_up = SPxLPBase<R>::upper(); 2297 } 2298 2299 /// copy lhs value VectorBase<R> to \p p_lhs. 2300 void getLhs(VectorBase<R>& p_lhs) const 2301 { 2302 p_lhs = SPxLPBase<R>::lhs(); 2303 } 2304 2305 /// copy rhs value VectorBase<R> to \p p_rhs. 2306 void getRhs(VectorBase<R>& p_rhs) const 2307 { 2308 p_rhs = SPxLPBase<R>::rhs(); 2309 } 2310 2311 /// optimization sense. 2312 typename SPxLPBase<R>::SPxSense sense() const 2313 { 2314 return this->spxSense(); 2315 } 2316 2317 /// returns statistical information in form of a string. 2318 std::string statistics() const 2319 { 2320 std::stringstream s; 2321 s << basis().statistics() 2322 << "Solution time : " << std::setw(10) << std::fixed << std::setprecision( 2323 2) << time() << std::endl 2324 << "Iterations : " << std::setw(10) << iterations() << std::endl; 2325 2326 return s.str(); 2327 } 2328 2329 /// returns whether a basis needs to be found for the improved dual simplex 2330 DecompStatus getDecompStatus() const 2331 { 2332 if(getStartingDecompBasis) 2333 return FINDSTARTBASIS; 2334 else 2335 return DONTFINDSTARTBASIS; 2336 } 2337 2338 /// sets whether the degeneracy is computed at each iteration 2339 void setComputeDegenFlag(bool computeDegen) 2340 { 2341 computeDegeneracy = computeDegen; 2342 } 2343 2344 2345 /// returns whether the degeneracy is computed in each iteration 2346 bool getComputeDegeneracy() const 2347 { 2348 return computeDegeneracy; 2349 } 2350 2351 2352 /// sets the offset for the number of iterations before the degeneracy is computed 2353 void setDegenCompOffset(int iterOffset) 2354 { 2355 degenCompIterOffset = iterOffset; 2356 } 2357 2358 2359 /// gets the offset for the number of iterations before the degeneracy is computed 2360 int getDegenCompOffset() const 2361 { 2362 return degenCompIterOffset; 2363 } 2364 2365 /// sets the iteration limit for the decomposition simplex initialisation 2366 void setDecompIterationLimit(int iterationLimit) 2367 { 2368 decompIterationLimit = iterationLimit; 2369 } 2370 2371 /// returns the iteration limit for the decomposition simplex initialisation 2372 int getDecompIterationLimit() const 2373 { 2374 return decompIterationLimit; 2375 } 2376 ///@} 2377 2378 //------------------------------------ 2379 /** Mapping between numbers and Ids */ 2380 ///@{ 2381 /// RowId of \p i 'th inequality. 2382 SPxRowId rowId(int i) const 2383 { 2384 return this->rId(i); 2385 } 2386 /// ColId of \p i 'th column. 2387 SPxColId colId(int i) const 2388 { 2389 return this->cId(i); 2390 } 2391 ///@} 2392 2393 //------------------------------------ 2394 /** Constructors / destructors */ 2395 ///@{ 2396 /// default constructor. 2397 explicit 2398 SPxSolverBase(Type type = LEAVE, 2399 Representation rep = ROW, 2400 Timer::TYPE ttype = Timer::USER_TIME); 2401 // virtual destructor 2402 virtual ~SPxSolverBase(); 2403 ///@} 2404 2405 //------------------------------------ 2406 /** Miscellaneous */ 2407 ///@{ 2408 /// check consistency. 2409 bool isConsistent() const; 2410 ///@} 2411 2412 //------------------------------------ 2413 /** assignment operator and copy constructor */ 2414 ///@{ 2415 /// assignment operator 2416 SPxSolverBase<R>& operator=(const SPxSolverBase<R>& base); 2417 /// copy constructor 2418 SPxSolverBase(const SPxSolverBase<R>& base); 2419 ///@} 2420 2421 void testVecs(); 2422 }; 2423 2424 // 2425 // Auxiliary functions. 2426 // 2427 2428 /// Pretty-printing of variable status. 2429 template <class R> 2430 std::ostream& operator<<(std::ostream& os, 2431 const typename SPxSolverBase<R>::VarStatus& status); 2432 2433 /// Pretty-printing of solver status. 2434 template <class R> 2435 std::ostream& operator<<(std::ostream& os, 2436 const typename SPxSolverBase<R>::Status& status); 2437 2438 /// Pretty-printing of algorithm. 2439 template <class R> 2440 std::ostream& operator<<(std::ostream& os, 2441 const typename SPxSolverBase<R>::Type& status); 2442 2443 /// Pretty-printing of representation. 2444 template <class R> 2445 std::ostream& operator<<(std::ostream& os, 2446 const typename SPxSolverBase<R>::Representation& status); 2447 2448 /* For Backwards compatibility */ 2449 typedef SPxSolverBase<Real> SPxSolver; 2450 2451 } // namespace soplex 2452 2453 // For general templated functions 2454 #include "spxsolver.hpp" 2455 #include "spxsolve.hpp" 2456 #include "changesoplex.hpp" 2457 #include "leave.hpp" 2458 #include "enter.hpp" 2459 #include "spxshift.hpp" 2460 #include "spxbounds.hpp" 2461 #include "spxchangebasis.hpp" 2462 #include "spxvecs.hpp" 2463 #include "spxwritestate.hpp" 2464 #include "spxfileio.hpp" 2465 #include "spxquality.hpp" 2466 2467 #endif // _SPXSOLVER_H_ 2468