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 soplex.h 26 * @brief Preconfigured SoPlex LP solver 27 */ 28 29 #ifndef _SOPLEX_H_ 30 #define _SOPLEX_H_ 31 32 #include <string.h> 33 34 #include "soplex/spxgithash.h" 35 #include "soplex/spxdefines.h" 36 #include "soplex/basevectors.h" 37 #include "soplex/spxsolver.h" 38 #include "soplex/slufactor.h" 39 #include "soplex/slufactor_rational.h" 40 41 ///@todo try to move to cpp file by forward declaration 42 #include "soplex/spxsimplifier.h" 43 #include "soplex/spxmainsm.h" 44 45 #include "soplex/spxscaler.h" 46 #include "soplex/spxequilisc.h" 47 #include "soplex/spxleastsqsc.h" 48 #include "soplex/spxgeometsc.h" 49 50 #include "soplex/spxstarter.h" 51 #include "soplex/spxweightst.h" 52 #include "soplex/spxsumst.h" 53 #include "soplex/spxvectorst.h" 54 55 #include "soplex/spxpricer.h" 56 #include "soplex/spxautopr.h" 57 #include "soplex/spxdantzigpr.h" 58 #include "soplex/spxparmultpr.h" 59 #include "soplex/spxdevexpr.h" 60 #include "soplex/spxsteeppr.h" 61 #include "soplex/spxsteepexpr.h" 62 #include "soplex/spxhybridpr.h" 63 64 #include "soplex/spxratiotester.h" 65 #include "soplex/spxdefaultrt.h" 66 #include "soplex/spxharrisrt.h" 67 #include "soplex/spxfastrt.h" 68 #include "soplex/spxboundflippingrt.h" 69 70 #include "soplex/solbase.h" 71 #include "soplex/sol.h" 72 73 #include "soplex/spxlpbase.h" 74 75 #include "soplex/spxpapilo.h" 76 77 #ifdef SOPLEX_WITH_GMP 78 #include <gmp.h> 79 #endif 80 81 #ifdef SOPLEX_WITH_BOOST 82 #ifdef SOPLEX_WITH_MPFR 83 // For multiple precision 84 #include <boost/multiprecision/mpfr.hpp> 85 #endif 86 #ifdef SOPLEX_WITH_CPPMPF 87 #include <boost/multiprecision/cpp_dec_float.hpp> 88 #endif 89 90 // An alias for boost multiprecision 91 namespace mpf = boost::multiprecision; 92 #endif 93 94 #define SOPLEX_DEFAULT_RANDOM_SEED 0 // used to suppress output when the seed was not changed 95 96 ///@todo implement automatic rep switch, based on row/col dim 97 ///@todo introduce status codes for SoPlex, especially for rational solving 98 99 ///@todo record and return "best" solutions found during IR (Ambros) 100 ///@todo implement main IR loop for primal and dual feasible case with fail otherwise (Ambros) 101 ///@todo implement statistical info (time, factor time, iters, ...) since last call to solveReal() or solveRational() (Ambros?) 102 ///@todo implement performInfeasibilityIR (Ambros?) 103 ///@todo extend IR loop to infeasible case (Dan?) 104 ///@todo extend IR loop to unbounded case (Dan?) 105 106 ///@todo interface rational reconstruction code for rational vectors 107 ///@todo integrate rational reconstruction into IR loop 108 ///@todo integrate rational SPxSolver and distinguish between original and transformed rational LP 109 ///@todo rational scalers 110 ///@todo rational simplifier 111 112 namespace soplex 113 { 114 115 /**@class SoPlex 116 * @brief Preconfigured SoPlex LP-solver. 117 * @ingroup Algo 118 */ 119 template <class R> 120 class SoPlexBase 121 { 122 public: 123 124 ///@name Construction and destruction 125 ///@{ 126 127 /// default constructor 128 SoPlexBase(); 129 130 /// assignment operator 131 SoPlexBase<R>& operator=(const SoPlexBase<R>& rhs); 132 133 /// copy constructor 134 SoPlexBase(const SoPlexBase<R>& rhs); 135 136 /// destructor 137 virtual ~SoPlexBase(); 138 139 ///@} 140 141 142 ///@name Access to the real LP 143 ///@{ 144 145 /// returns number of rows 146 int numRows() const; 147 int numRowsReal() const; /* For SCIP compatibility */ 148 int numRowsRational() const; 149 150 /// Templated function that 151 /// returns number of columns 152 int numCols() const; 153 int numColsReal() const; /* For SCIP compatibility */ 154 int numColsRational() const; 155 156 /// returns number of nonzeros 157 int numNonzeros() const; 158 159 int numNonzerosRational() const; 160 161 /// returns smallest non-zero element in absolute value 162 R minAbsNonzeroReal() const; 163 164 /// returns biggest non-zero element in absolute value 165 R maxAbsNonzeroReal() const; 166 167 /// returns (unscaled) coefficient 168 R coefReal(int row, int col) const; 169 170 /// returns vector of row \p i, ignoring scaling 171 const SVectorBase<R>& rowVectorRealInternal(int i) const; 172 173 /// gets vector of row \p i 174 void getRowVectorReal(int i, DSVectorBase<R>& row) const; 175 176 /// returns right-hand side vector, ignoring scaling 177 const VectorBase<R>& rhsRealInternal() const; 178 179 /// gets right-hand side vector 180 void getRhsReal(VectorBase<R>& rhs) const; 181 182 /// returns right-hand side of row \p i 183 R rhsReal(int i) const; 184 185 /// returns left-hand side vector, ignoring scaling 186 const VectorBase<R>& lhsRealInternal() const; 187 188 /// gets left-hand side vector 189 void getLhsReal(VectorBase<R>& lhs) const; 190 191 /// returns left-hand side of row \p i 192 R lhsReal(int i) const; 193 194 /// returns inequality type of row \p i 195 typename LPRowBase<R>::Type rowTypeReal(int i) const; 196 197 /// returns vector of col \p i, ignoring scaling 198 const SVectorBase<R>& colVectorRealInternal(int i) const; 199 200 /// gets vector of col \p i 201 void getColVectorReal(int i, DSVectorBase<R>& col) const; 202 203 /// returns upper bound vector 204 const VectorBase<R>& upperRealInternal() const; 205 206 /// returns upper bound of column \p i 207 R upperReal(int i) const; 208 209 /// gets upper bound vector 210 void getUpperReal(VectorBase<R>& upper) const; 211 212 /// returns lower bound vector 213 const VectorBase<R>& lowerRealInternal() const; 214 215 /// returns lower bound of column \p i 216 R lowerReal(int i) const; 217 218 /// gets lower bound vector 219 void getLowerReal(VectorBase<R>& lower) const; 220 221 /// gets objective function vector 222 void getObjReal(VectorBase<R>& obj) const; 223 224 /// returns objective value of column \p i 225 R objReal(int i) const; 226 227 /// returns objective function vector after transformation to a maximization problem; since this is how it is stored 228 /// internally, this is generally faster 229 const VectorBase<R>& maxObjRealInternal() const; 230 231 /// returns objective value of column \p i after transformation to a maximization problem; since this is how it is 232 /// stored internally, this is generally faster 233 R maxObjReal(int i) const; 234 235 /// gets number of available dual norms 236 void getNdualNorms(int& nnormsRow, int& nnormsCol) const; 237 238 /// gets steepest edge norms and returns false if they are not available 239 bool getDualNorms(int& nnormsRow, int& nnormsCol, R* norms) const; 240 241 /// sets steepest edge norms and returns false if that's not possible 242 bool setDualNorms(int nnormsRow, int nnormsCol, R* norms); 243 244 /// pass integrality information about the variables to the solver 245 void setIntegralityInformation(int ncols, int* intInfo); 246 247 ///@} 248 249 250 ///@name Access to the rational LP 251 ///@{ 252 253 /// returns smallest non-zero element in absolute value 254 Rational minAbsNonzeroRational() const; 255 256 /// returns biggest non-zero element in absolute value 257 Rational maxAbsNonzeroRational() const; 258 259 /// gets row \p i 260 void getRowRational(int i, LPRowRational& lprow) const; 261 262 /// gets rows \p start, ..., \p end. 263 void getRowsRational(int start, int end, LPRowSetRational& lprowset) const; 264 265 /// returns vector of row \p i 266 const SVectorRational& rowVectorRational(int i) const; 267 268 /// returns right-hand side vector 269 const VectorRational& rhsRational() const; 270 271 /// returns right-hand side of row \p i 272 const Rational& rhsRational(int i) const; 273 274 /// returns left-hand side vector 275 const VectorRational& lhsRational() const; 276 277 /// returns left-hand side of row \p i 278 const Rational& lhsRational(int i) const; 279 280 /// returns inequality type of row \p i 281 LPRowRational::Type rowTypeRational(int i) const; 282 283 /// gets column \p i 284 void getColRational(int i, LPColRational& lpcol) const; 285 286 /// gets columns \p start, ..., \p end 287 void getColsRational(int start, int end, LPColSetRational& lpcolset) const; 288 289 /// returns vector of column \p i 290 const SVectorRational& colVectorRational(int i) const; 291 292 /// returns upper bound vector 293 const VectorRational& upperRational() const; 294 295 /// returns upper bound of column \p i 296 const Rational& upperRational(int i) const; 297 298 /// returns lower bound vector 299 const VectorRational& lowerRational() const; 300 301 /// returns lower bound of column \p i 302 const Rational& lowerRational(int i) const; 303 304 /// gets objective function vector 305 void getObjRational(VectorRational& obj) const; 306 307 /// gets objective value of column \p i 308 void getObjRational(int i, Rational& obj) const; 309 310 /// returns objective value of column \p i 311 Rational objRational(int i) const; 312 313 /// returns objective function vector after transformation to a maximization problem; since this is how it is stored 314 /// internally, this is generally faster 315 const VectorRational& maxObjRational() const; 316 317 /// returns objective value of column \p i after transformation to a maximization problem; since this is how it is 318 /// stored internally, this is generally faster 319 const Rational& maxObjRational(int i) const; 320 321 ///@} 322 323 324 ///@name Modification of the real LP 325 ///@{ 326 327 /// adds a single row 328 void addRowReal(const LPRowBase<R>& lprow); 329 330 /// adds multiple rows 331 void addRowsReal(const LPRowSetBase<R>& lprowset); 332 333 /// adds a single column 334 void addColReal(const LPColBase<R>& lpcol); 335 336 /// adds multiple columns 337 void addColsReal(const LPColSetBase<R>& lpcolset); 338 339 /// replaces row \p i with \p lprow 340 void changeRowReal(int i, const LPRowBase<R>& lprow); 341 342 /// changes left-hand side vector for constraints to \p lhs 343 void changeLhsReal(const VectorBase<R>& lhs); 344 345 /// changes left-hand side of row \p i to \p lhs 346 void changeLhsReal(int i, const R& lhs); 347 348 /// changes right-hand side vector to \p rhs 349 void changeRhsReal(const VectorBase<R>& rhs); 350 351 /// changes right-hand side of row \p i to \p rhs 352 void changeRhsReal(int i, const R& rhs); 353 354 /// changes left- and right-hand side vectors 355 void changeRangeReal(const VectorBase<R>& lhs, const VectorBase<R>& rhs); 356 357 /// changes left- and right-hand side of row \p i 358 void changeRangeReal(int i, const R& lhs, const R& rhs); 359 360 /// replaces column \p i with \p lpcol 361 void changeColReal(int i, const LPColReal& lpcol); 362 363 /// changes vector of lower bounds to \p lower 364 void changeLowerReal(const VectorBase<R>& lower); 365 366 /// changes lower bound of column i to \p lower 367 void changeLowerReal(int i, const R& lower); 368 369 /// changes vector of upper bounds to \p upper 370 void changeUpperReal(const VectorBase<R>& upper); 371 372 /// changes \p i 'th upper bound to \p upper 373 void changeUpperReal(int i, const R& upper); 374 375 /// changes vectors of column bounds to \p lower and \p upper 376 void changeBoundsReal(const VectorBase<R>& lower, const VectorBase<R>& upper); 377 378 /// changes bounds of column \p i to \p lower and \p upper 379 void changeBoundsReal(int i, const R& lower, const R& upper); 380 381 /// changes objective function vector to \p obj 382 void changeObjReal(const VectorBase<R>& obj); 383 384 /// changes objective coefficient of column i to \p obj 385 void changeObjReal(int i, const R& obj); 386 387 /// changes matrix entry in row \p i and column \p j to \p val 388 void changeElementReal(int i, int j, const R& val); 389 390 /// removes row \p i 391 void removeRowReal(int i); 392 393 /// removes all rows with an index \p i such that \p perm[i] < 0; upon completion, \p perm[i] >= 0 indicates the 394 /// new index where row \p i has been moved to; note that \p perm must point to an array of size at least 395 /// #numRows() 396 void removeRowsReal(int perm[]); 397 398 /// remove all rows with indices in array \p idx of size \p n; an array \p perm of size #numRows() may be passed 399 /// as buffer memory 400 void removeRowsReal(int idx[], int n, int perm[] = 0); 401 402 /// removes rows \p start to \p end including both; an array \p perm of size #numRows() may be passed as buffer 403 /// memory 404 void removeRowRangeReal(int start, int end, int perm[] = 0); 405 406 /// removes column i 407 void removeColReal(int i); 408 409 /// removes all columns with an index \p i such that \p perm[i] < 0; upon completion, \p perm[i] >= 0 indicates the 410 /// new index where column \p i has been moved to; note that \p perm must point to an array of size at least 411 /// #numColsReal() 412 void removeColsReal(int perm[]); 413 414 /// remove all columns with indices in array \p idx of size \p n; an array \p perm of size #numColsReal() may be 415 /// passed as buffer memory 416 void removeColsReal(int idx[], int n, int perm[] = 0); 417 418 /// removes columns \p start to \p end including both; an array \p perm of size #numColsReal() may be passed as 419 /// buffer memory 420 void removeColRangeReal(int start, int end, int perm[] = 0); 421 422 /// clears the LP 423 void clearLPReal(); 424 425 /// synchronizes real LP with rational LP, i.e., copies (rounded) rational LP into real LP, if sync mode is manual 426 void syncLPReal(); 427 428 ///@} 429 430 431 ///@name Modification of the rational LP 432 ///@{ 433 434 /// adds a single row 435 void addRowRational(const LPRowRational& lprow); 436 437 #ifdef SOPLEX_WITH_GMP 438 /// adds a single row (GMP only method) 439 void addRowRational(const mpq_t* lhs, const mpq_t* rowValues, const int* rowIndices, 440 const int rowSize, const mpq_t* rhs); 441 442 /// adds a set of rows (GMP only method) 443 void addRowsRational(const mpq_t* lhs, const mpq_t* rowValues, const int* rowIndices, 444 const int* rowStarts, const int* rowLengths, const int numRows, const int numValues, 445 const mpq_t* rhs); 446 #endif 447 448 /// adds multiple rows 449 void addRowsRational(const LPRowSetRational& lprowset); 450 451 /// adds a single column 452 void addColRational(const LPColRational& lpcol); 453 454 #ifdef SOPLEX_WITH_GMP 455 /// adds a single column (GMP only method) 456 void addColRational(const mpq_t* obj, const mpq_t* lower, const mpq_t* colValues, 457 const int* colIndices, const int colSize, const mpq_t* upper); 458 459 /// adds a set of columns (GMP only method) 460 void addColsRational(const mpq_t* obj, const mpq_t* lower, const mpq_t* colValues, 461 const int* colIndices, const int* colStarts, const int* colLengths, const int numCols, 462 const int numValues, const mpq_t* upper); 463 #endif 464 465 /// adds multiple columns 466 void addColsRational(const LPColSetRational& lpcolset); 467 468 /// replaces row \p i with \p lprow 469 void changeRowRational(int i, const LPRowRational& lprow); 470 471 /// changes left-hand side vector for constraints to \p lhs 472 void changeLhsRational(const VectorRational& lhs); 473 474 /// changes left-hand side of row \p i to \p lhs 475 void changeLhsRational(int i, const Rational& lhs); 476 477 #ifdef SOPLEX_WITH_GMP 478 /// changes left-hand side of row \p i to \p lhs (GMP only method) 479 void changeLhsRational(int i, const mpq_t* lhs); 480 #endif 481 482 /// changes right-hand side vector to \p rhs 483 void changeRhsRational(const VectorRational& rhs); 484 485 #ifdef SOPLEX_WITH_GMP 486 /// changes right-hand side vector to \p rhs (GMP only method) 487 void changeRhsRational(const mpq_t* rhs, int rhsSize); 488 #endif 489 490 /// changes right-hand side of row \p i to \p rhs 491 void changeRhsRational(int i, const Rational& rhs); 492 493 /// changes left- and right-hand side vectors 494 void changeRangeRational(const VectorRational& lhs, const VectorRational& rhs); 495 496 /// changes left- and right-hand side of row \p i 497 void changeRangeRational(int i, const Rational& lhs, const Rational& rhs); 498 499 #ifdef SOPLEX_WITH_GMP 500 /// changes left- and right-hand side of row \p i (GMP only method) 501 void changeRangeRational(int i, const mpq_t* lhs, const mpq_t* rhs); 502 #endif 503 504 /// replaces column \p i with \p lpcol 505 void changeColRational(int i, const LPColRational& lpcol); 506 507 /// changes vector of lower bounds to \p lower 508 void changeLowerRational(const VectorRational& lower); 509 510 /// changes lower bound of column i to \p lower 511 void changeLowerRational(int i, const Rational& lower); 512 513 #ifdef SOPLEX_WITH_GMP 514 /// changes lower bound of column i to \p lower (GMP only method) 515 void changeLowerRational(int i, const mpq_t* lower); 516 #endif 517 518 /// changes vector of upper bounds to \p upper 519 void changeUpperRational(const VectorRational& upper); 520 521 /// changes \p i 'th upper bound to \p upper 522 void changeUpperRational(int i, const Rational& upper); 523 524 #ifdef SOPLEX_WITH_GMP 525 /// changes upper bound of column i to \p upper (GMP only method) 526 void changeUpperRational(int i, const mpq_t* upper); 527 #endif 528 529 /// changes vectors of column bounds to \p lower and \p upper 530 void changeBoundsRational(const VectorRational& lower, const VectorRational& upper); 531 532 /// changes bounds of column \p i to \p lower and \p upper 533 void changeBoundsRational(int i, const Rational& lower, const Rational& upper); 534 535 #ifdef SOPLEX_WITH_GMP 536 /// changes bounds of column \p i to \p lower and \p upper (GMP only method) 537 void changeBoundsRational(int i, const mpq_t* lower, const mpq_t* upper); 538 #endif 539 540 /// changes objective function vector to \p obj 541 void changeObjRational(const VectorRational& obj); 542 543 /// changes objective coefficient of column i to \p obj 544 void changeObjRational(int i, const Rational& obj); 545 546 #ifdef SOPLEX_WITH_GMP 547 /// changes objective coefficient of column i to \p obj (GMP only method) 548 void changeObjRational(int i, const mpq_t* obj); 549 #endif 550 551 /// changes matrix entry in row \p i and column \p j to \p val 552 void changeElementRational(int i, int j, const Rational& val); 553 554 #ifdef SOPLEX_WITH_GMP 555 /// changes matrix entry in row \p i and column \p j to \p val (GMP only method) 556 void changeElementRational(int i, int j, const mpq_t* val); 557 #endif 558 559 /// removes row \p i 560 void removeRowRational(int i); 561 562 /// removes all rows with an index \p i such that \p perm[i] < 0; upon completion, \p perm[i] >= 0 indicates the new 563 /// index where row \p i has been moved to; note that \p perm must point to an array of size at least 564 /// #numRowsRational() 565 void removeRowsRational(int perm[]); 566 567 /// remove all rows with indices in array \p idx of size \p n; an array \p perm of size #numRowsRational() may be 568 /// passed as buffer memory 569 void removeRowsRational(int idx[], int n, int perm[] = 0); 570 571 /// removes rows \p start to \p end including both; an array \p perm of size #numRowsRational() may be passed as 572 /// buffer memory 573 void removeRowRangeRational(int start, int end, int perm[] = 0); 574 575 /// removes column i 576 void removeColRational(int i); 577 578 /// removes all columns with an index \p i such that \p perm[i] < 0; upon completion, \p perm[i] >= 0 indicates the 579 /// new index where column \p i has been moved to; note that \p perm must point to an array of size at least 580 /// #numColsRational() 581 void removeColsRational(int perm[]); 582 583 /// remove all columns with indices in array \p idx of size \p n; an array \p perm of size #numColsRational() may be 584 /// passed as buffer memory 585 void removeColsRational(int idx[], int n, int perm[] = 0); 586 587 /// removes columns \p start to \p end including both; an array \p perm of size #numColsRational() may be passed as 588 /// buffer memory 589 void removeColRangeRational(int start, int end, int perm[] = 0); 590 591 /// clears the LP 592 void clearLPRational(); 593 594 /// synchronizes rational LP with real LP, i.e., copies real LP to rational LP, if sync mode is manual 595 void syncLPRational(); 596 597 ///@} 598 599 ///@name Solving and general solution query 600 ///@{ 601 602 /// optimize the given LP 603 typename SPxSolverBase<R>::Status optimize(volatile bool* interrupt = NULL); 604 605 // old name for backwards compatibility 606 typename SPxSolverBase<R>::Status solve(volatile bool* interrupt = NULL) 607 { 608 return optimize(interrupt); 609 } 610 611 /// returns the current solver status 612 typename SPxSolverBase<R>::Status status() const; 613 614 /// is stored primal solution feasible? 615 bool isPrimalFeasible() const; 616 617 /// is a solution available (not neccessarily feasible)? 618 bool hasSol() const; 619 620 /// deprecated: use #hasSol() instead 621 bool hasPrimal() const 622 { 623 return hasSol(); 624 } 625 626 /// deprecated: use #hasSol() instead 627 bool hasDual() const 628 { 629 return hasSol(); 630 } 631 632 /// is a primal unbounded ray available? 633 bool hasPrimalRay() const; 634 635 /// is stored dual solution feasible? 636 bool isDualFeasible() const; 637 638 /// is Farkas proof of infeasibility available? 639 bool hasDualFarkas() const; 640 641 /// sets the status to OPTIMAL in case the LP has been solved with unscaled violations 642 bool ignoreUnscaledViolations() 643 { 644 if(_status == SPxSolverBase<R>::OPTIMAL_UNSCALED_VIOLATIONS) 645 { 646 _status = SPxSolverBase<R>::OPTIMAL; 647 return true; 648 } 649 else 650 return false; 651 } 652 ///@} 653 654 655 ///@name Query for the real solution data 656 ///@{ 657 658 /// returns the objective value if a primal solution is available 659 R objValueReal(); 660 661 /// gets the primal solution vector if available; returns true on success 662 bool getPrimal(VectorBase<R>& vector); 663 bool getPrimalReal(R* p_vector, int size); // For SCIP compatibility 664 bool getPrimalRational(VectorRational& vector); 665 666 /// gets the vector of slack values if available; returns true on success 667 bool getSlacksReal(VectorBase<R>& vector); 668 bool getSlacksReal(R* p_vector, int dim); 669 670 /// gets the primal ray if available; returns true on success 671 bool getPrimalRay(VectorBase<R>& vector); 672 bool getPrimalRayReal(R* vector, int dim); /* For SCIP compatibility */ 673 bool getPrimalRayRational(VectorRational& vector); 674 675 /// gets the dual solution vector if available; returns true on success 676 bool getDual(VectorBase<R>& vector); 677 bool getDualReal(R* p_vector, int dim); /* For SCIP compatibility */ 678 bool getDualRational(VectorRational& vector); 679 680 /// gets the vector of reduced cost values if available; returns true on success 681 bool getRedCost(VectorBase<R>& vector); 682 bool getRedCostReal(R* vector, int dim); /* For SCIP compatibility */ 683 bool getRedCostRational(VectorRational& vector); 684 685 /// gets the Farkas proof if available; returns true on success 686 bool getDualFarkas(VectorBase<R>& vector); 687 bool getDualFarkasReal(R* vector, int dim); 688 bool getDualFarkasRational(VectorRational& vector); 689 690 /// gets violation of bounds; returns true on success 691 bool getBoundViolation(R& maxviol, R& sumviol); 692 bool getBoundViolationRational(Rational& maxviol, Rational& sumviol); 693 694 /// gets violation of constraints; returns true on success 695 bool getRowViolation(R& maxviol, R& sumviol); 696 bool getRowViolationRational(Rational& maxviol, Rational& sumviol); 697 698 /// gets violation of reduced costs; returns true on success 699 bool getRedCostViolation(R& maxviol, R& sumviol); 700 bool getRedCostViolationRational(Rational& maxviol, Rational& sumviol); 701 702 /// gets violation of dual multipliers; returns true on success 703 bool getDualViolation(R& maxviol, R& sumviol); 704 bool getDualViolationRational(Rational& maxviol, Rational& sumviol); 705 706 ///@} 707 708 709 ///@name Query for the rational solution data 710 ///@{ 711 712 /// returns the objective value if a primal solution is available 713 Rational objValueRational(); 714 715 /// gets the vector of slack values if available; returns true on success 716 bool getSlacksRational(VectorRational& vector); 717 718 #ifdef SOPLEX_WITH_GMP 719 /// gets the primal solution vector if available; returns true on success (GMP only method) 720 bool getPrimalRational(mpq_t* vector, const int size); 721 722 /// gets the vector of slack values if available; returns true on success (GMP only method) 723 bool getSlacksRational(mpq_t* vector, const int size); 724 725 /// gets the primal ray if LP is unbounded; returns true on success (GMP only method) 726 bool getPrimalRayRational(mpq_t* vector, const int size); 727 728 /// gets the dual solution vector if available; returns true on success (GMP only method) 729 bool getDualRational(mpq_t* vector, const int size); 730 731 /// gets the vector of reduced cost values if available; returns true on success (GMP only method) 732 bool getRedCostRational(mpq_t* vector, const int size); 733 734 /// gets the Farkas proof if LP is infeasible; returns true on success (GMP only method) 735 bool getDualFarkasRational(mpq_t* vector, const int size); 736 #endif 737 738 /// get size of primal solution 739 int totalSizePrimalRational(const int base = 2); 740 741 /// get size of dual solution 742 int totalSizeDualRational(const int base = 2); 743 744 /// get size of least common multiple of denominators in primal solution 745 int dlcmSizePrimalRational(const int base = 2); 746 747 /// get size of least common multiple of denominators in dual solution 748 int dlcmSizeDualRational(const int base = 2); 749 750 /// get size of largest denominator in primal solution 751 int dmaxSizePrimalRational(const int base = 2); 752 753 /// get size of largest denominator in dual solution 754 int dmaxSizeDualRational(const int base = 2); 755 756 ///@} 757 758 759 ///@name Access and modification of basis information 760 ///@{ 761 762 /// is an advanced starting basis available? 763 bool hasBasis() const; 764 765 /// returns the current basis status 766 typename SPxBasisBase<R>::SPxStatus basisStatus() const; 767 768 /// returns basis status for a single row 769 typename SPxSolverBase<R>::VarStatus basisRowStatus(int row) const; 770 771 /// returns basis status for a single column 772 typename SPxSolverBase<R>::VarStatus basisColStatus(int col) const; 773 774 /// gets current basis via arrays of statuses 775 void getBasis(typename SPxSolverBase<R>::VarStatus rows[], 776 typename SPxSolverBase<R>::VarStatus cols[]) const; 777 778 /// gets the indices of the basic columns and rows; basic column n gives value n, basic row m gives value -1-m 779 void getBasisInd(int* bind) const; 780 781 /** compute one of several matrix metrics based on the diagonal of the LU factorization 782 * type = 0: max/min ratio 783 * type = 1: trace of U (sum of diagonal elements) 784 * type = 2: determinant (product of diagonal elements) 785 */ 786 bool getBasisMetric(R& metric, int type = 0); 787 788 /// computes an estimated condition number for the current basis matrix using the power method; returns true on success 789 bool getEstimatedCondition(R& condition); 790 791 /// computes the exact condition number for the current basis matrix using the power method; returns true on success 792 bool getExactCondition(R& condition); 793 794 /// computes row \p r of basis inverse; returns true on success 795 /// @param r which row of the basis inverse is computed 796 /// @param coef values of result vector (not packed but scattered) 797 /// @param inds indices of result vector (NULL if not to be used) 798 /// @param ninds number of nonzeros in result vector 799 /// @param unscale determines whether the result should be unscaled according to the original LP data 800 bool getBasisInverseRowReal(int r, R* coef, int* inds = NULL, int* ninds = NULL, 801 bool unscale = true); 802 803 /// computes column \p c of basis inverse; returns true on success 804 /// @param c which column of the basis inverse is computed 805 /// @param coef values of result vector (not packed but scattered) 806 /// @param inds indices of result vector (NULL if not to be used) 807 /// @param ninds number of nonzeros in result vector 808 /// @param unscale determines whether the result should be unscaled according to the original LP data 809 bool getBasisInverseColReal(int c, R* coef, int* inds = NULL, int* ninds = NULL, 810 bool unscale = true); 811 812 /// computes dense solution of basis matrix B * \p sol = \p rhs; returns true on success 813 bool getBasisInverseTimesVecReal(R* rhs, R* sol, bool unscale = true); 814 815 /// multiply with basis matrix; B * \p vec (inplace) 816 /// @param vec (dense) vector to be multiplied with 817 /// @param unscale determines whether the result should be unscaled according to the original LP data 818 bool multBasis(R* vec, bool unscale = true); 819 820 /// multiply with transpose of basis matrix; \p vec * B^T (inplace) 821 /// @param vec (dense) vector to be multiplied with 822 /// @param unscale determines whether the result should be unscaled according to the original LP data 823 bool multBasisTranspose(R* vec, bool unscale = true); 824 825 /// compute rational basis inverse; returns true on success 826 bool computeBasisInverseRational(); 827 828 /// gets an array of indices for the columns of the rational basis matrix; bind[i] >= 0 means that the i-th column of 829 /// the basis matrix contains variable bind[i]; bind[i] < 0 means that the i-th column of the basis matrix contains 830 /// the slack variable for row -bind[i]-1; performs rational factorization if not available; returns true on success 831 bool getBasisIndRational(DataArray<int>& bind); 832 833 /// computes row r of basis inverse; performs rational factorization if not available; returns true on success 834 bool getBasisInverseRowRational(const int r, SSVectorRational& vec); 835 836 /// computes column c of basis inverse; performs rational factorization if not available; returns true on success 837 bool getBasisInverseColRational(const int c, SSVectorRational& vec); 838 839 /// computes solution of basis matrix B * sol = rhs; performs rational factorization if not available; returns true 840 /// on success 841 bool getBasisInverseTimesVecRational(const SVectorRational& rhs, SSVectorRational& sol); 842 843 /// sets starting basis via arrays of statuses 844 void setBasis(const typename SPxSolverBase<R>::VarStatus rows[], 845 const typename SPxSolverBase<R>::VarStatus cols[]); 846 847 /// clears starting basis 848 void clearBasis(); 849 850 ///@} 851 852 853 ///@name Statistical information 854 ///@{ 855 856 /// number of iterations since last call to solve 857 int numIterations() const; 858 859 /// number of precision boosts since last call to solve 860 int numPrecisionBoosts() const; 861 862 /// number of iterations in higher precision since last call to solve 863 int numIterationsBoosted() const; 864 865 /// time spen in higher precision since last call to solve 866 Real precisionBoostTime() const; 867 868 /// time spent in last call to solve 869 Real solveTime() const; 870 871 /// statistical information in form of a string 872 std::string statisticString() const; 873 874 /// name of starter 875 const char* getStarterName(); 876 877 /// name of simplifier 878 const char* getSimplifierName(); 879 880 /// name of scaling method 881 const char* getScalerName(); 882 883 /// name of currently loaded pricer 884 const char* getPricerName(); 885 886 /// name of currently loaded ratiotester 887 const char* getRatiotesterName(); 888 889 ///@} 890 891 892 ///@name File I/O 893 ///@{ 894 895 /// reads LP file in LP or MPS format according to READMODE parameter; gets row names, column names, and 896 /// integer variables if desired; returns true on success 897 bool readFile(const char* filename, NameSet* rowNames = 0, NameSet* colNames = 0, 898 DIdxSet* intVars = 0); 899 900 /// Templated write function 901 /// Real 902 /// writes real LP to file; LP or MPS format is chosen from the extension in \p filename; if \p rowNames and \p 903 /// colNames are \c NULL, default names are used; if \p intVars is not \c NULL, the variables contained in it are 904 /// marked as integer; returns true on success 905 /// Rational 906 /// writes rational LP to file; LP or MPS format is chosen from the extension in \p filename; if \p rowNames and \p 907 /// colNames are \c NULL, default names are used; if \p intVars is not \c NULL, the variables contained in it are 908 /// marked as integer; returns true on success 909 bool writeFile(const char* filename, const NameSet* rowNames = 0, const NameSet* colNames = 0, 910 const DIdxSet* intvars = 0, const bool unscale = true) const; 911 912 bool writeFileRational(const char* filename, const NameSet* rowNames = 0, 913 const NameSet* colNames = 0, const DIdxSet* intvars = 0) const; 914 915 /* For SCIP compatibility */ 916 bool writeFileReal(const char* filename, const NameSet* rowNames = 0, const NameSet* colNames = 0, 917 const DIdxSet* intvars = 0, const bool unscale = true) const; 918 919 /// writes the dual of the real LP to file; LP or MPS format is chosen from the extension in \p filename; 920 /// if \p rowNames and \p colNames are \c NULL, default names are used; if \p intVars is not \c NULL, 921 /// the variables contained in it are marked as integer; returns true on success 922 bool writeDualFileReal(const char* filename, const NameSet* rowNames = 0, 923 const NameSet* colNames = 0, const DIdxSet* intvars = 0) const; 924 925 /// reads basis information from \p filename and returns true on success; if \p rowNames and \p colNames are \c NULL, 926 /// default names are assumed; returns true on success 927 bool readBasisFile(const char* filename, const NameSet* rowNames = 0, const NameSet* colNames = 0); 928 929 /// writes basis information to \p filename; if \p rowNames and \p colNames are \c NULL, default names are used; 930 /// returns true on success 931 bool writeBasisFile(const char* filename, const NameSet* rowNames = 0, const NameSet* colNames = 0, 932 const bool cpxFormat = false) const; 933 934 /// writes internal LP, basis information, and parameter settings; if \p rowNames and \p colNames are \c NULL, 935 /// default names are used 936 void writeStateReal(const char* filename, const NameSet* rowNames = 0, const NameSet* colNames = 0, 937 const bool cpxFormat = false) const; 938 939 /// writes internal LP, basis information, and parameter settings; if \p rowNames and \p colNames are \c NULL, 940 /// default names are used 941 void writeStateRational(const char* filename, const NameSet* rowNames = 0, 942 const NameSet* colNames = 0, const bool cpxFormat = false) const; 943 944 ///@} 945 946 947 ///@name Parameters 948 ///@{ 949 950 /// boolean parameters 951 typedef enum 952 { 953 /// should lifting be used to reduce range of nonzero matrix coefficients? 954 LIFTING = 0, 955 956 /// should LP be transformed to equality form before a rational solve? 957 EQTRANS = 1, 958 959 /// should dual infeasibility be tested in order to try to return a dual solution even if primal infeasible? 960 TESTDUALINF = 2, 961 962 /// should a rational factorization be performed after iterative refinement? 963 RATFAC = 3, 964 965 /// should the decomposition based dual simplex be used to solve the LP? Setting this to true forces the solve mode to 966 /// SOLVEMODE_REAL and the basis representation to REPRESENTATION_ROW 967 USEDECOMPDUALSIMPLEX = 4, 968 969 /// should the degeneracy be computed for each basis? 970 COMPUTEDEGEN = 5, 971 972 /// should the dual of the complementary problem be used in the decomposition simplex? 973 USECOMPDUAL = 6, 974 975 /// should row and bound violations be computed explicitly in the update of reduced problem in the decomposition simplex 976 EXPLICITVIOL = 7, 977 978 /// should cycling solutions be accepted during iterative refinement? 979 ACCEPTCYCLING = 8, 980 981 /// apply rational reconstruction after each iterative refinement? 982 RATREC = 9, 983 984 /// round scaling factors for iterative refinement to powers of two? 985 POWERSCALING = 10, 986 987 /// continue iterative refinement with exact basic solution if not optimal? 988 RATFACJUMP = 11, 989 990 /// use bound flipping also for row representation? 991 ROWBOUNDFLIPS = 12, 992 993 /// use persistent scaling? 994 PERSISTENTSCALING = 13, 995 996 /// perturb the entire problem or only the relevant bounds of s single pivot? 997 FULLPERTURBATION = 14, 998 999 /// re-optimize the original problem to get a proof (ray) of infeasibility/unboundedness? 1000 ENSURERAY = 15, 1001 1002 /// try to enforce that the optimal solution is a basic solution 1003 FORCEBASIC = 16, 1004 1005 // enable presolver SingletonCols in PaPILO? 1006 SIMPLIFIER_SINGLETONCOLS = 17, 1007 1008 // enable presolver ConstraintPropagation in PaPILO? 1009 SIMPLIFIER_CONSTRAINTPROPAGATION = 18, 1010 1011 // enable presolver ParallelRowDetection in PaPILO? 1012 SIMPLIFIER_PARALLELROWDETECTION = 19, 1013 1014 // enable presolver ParallelColDetection in PaPILO? 1015 SIMPLIFIER_PARALLELCOLDETECTION = 20, 1016 1017 // enable presolver SingletonStuffing in PaPILO? 1018 SIMPLIFIER_SINGLETONSTUFFING = 21, 1019 1020 // enable presolver DualFix in PaPILO? 1021 SIMPLIFIER_DUALFIX = 22, 1022 1023 // enable presolver FixContinuous in PaPILO? 1024 SIMPLIFIER_FIXCONTINUOUS = 23, 1025 1026 // enable presolver DominatedCols in PaPILO? 1027 SIMPLIFIER_DOMINATEDCOLS = 24, 1028 1029 // enable iterative refinement ? 1030 ITERATIVE_REFINEMENT = 25, 1031 1032 /// adapt tolerances to the multiprecision used 1033 ADAPT_TOLS_TO_MULTIPRECISION = 26, 1034 1035 /// enable precision boosting ? 1036 PRECISION_BOOSTING = 27, 1037 1038 /// boosted solver start from last basis 1039 BOOSTED_WARM_START = 28, 1040 1041 /// try different settings when solve fails 1042 RECOVERY_MECHANISM = 29, 1043 1044 /// store advanced and stable basis met before each simplex iteration, to better warm start 1045 STORE_BASIS_BEFORE_SIMPLEX_PIVOT = 30, 1046 1047 /// number of boolean parameters 1048 BOOLPARAM_COUNT = 31 1049 } BoolParam; 1050 1051 /// integer parameters 1052 typedef enum 1053 { 1054 /// objective sense 1055 OBJSENSE = 0, 1056 1057 /// type of computational form, i.e., column or row representation 1058 REPRESENTATION = 1, 1059 1060 /// type of algorithm, i.e., primal or dual 1061 ALGORITHM = 2, 1062 1063 /// type of LU update 1064 FACTOR_UPDATE_TYPE = 3, 1065 1066 /// maximum number of updates without fresh factorization 1067 FACTOR_UPDATE_MAX = 4, 1068 1069 /// iteration limit (-1 if unlimited) 1070 ITERLIMIT = 5, 1071 1072 /// refinement limit (-1 if unlimited) 1073 REFLIMIT = 6, 1074 1075 /// stalling refinement limit (-1 if unlimited) 1076 STALLREFLIMIT = 7, 1077 1078 /// display frequency 1079 DISPLAYFREQ = 8, 1080 1081 /// verbosity level 1082 VERBOSITY = 9, 1083 1084 /// type of simplifier 1085 SIMPLIFIER = 10, 1086 1087 /// type of scaler 1088 SCALER = 11, 1089 1090 /// type of starter used to create crash basis 1091 STARTER = 12, 1092 1093 /// type of pricer 1094 PRICER = 13, 1095 1096 /// type of ratio test 1097 RATIOTESTER = 14, 1098 1099 /// mode for synchronizing real and rational LP 1100 SYNCMODE = 15, 1101 1102 /// mode for reading LP files 1103 READMODE = 16, 1104 1105 /// mode for iterative refinement strategy 1106 SOLVEMODE = 17, 1107 1108 /// mode for a posteriori feasibility checks 1109 CHECKMODE = 18, 1110 1111 /// type of timer 1112 TIMER = 19, 1113 1114 /// mode for hyper sparse pricing 1115 HYPER_PRICING = 20, 1116 1117 /// minimum number of stalling refinements since last pivot to trigger rational factorization 1118 RATFAC_MINSTALLS = 21, 1119 1120 /// maximum number of conjugate gradient iterations in least square scaling 1121 LEASTSQ_MAXROUNDS = 22, 1122 1123 /// mode for solution polishing 1124 SOLUTION_POLISHING = 23, 1125 1126 /// the number of iterations before the decomposition simplex initialisation is terminated. 1127 DECOMP_ITERLIMIT = 24, 1128 1129 /// the maximum number of rows that are added in each iteration of the decomposition based simplex 1130 DECOMP_MAXADDEDROWS = 25, 1131 1132 /// the iteration frequency at which the decomposition solve output is displayed. 1133 DECOMP_DISPLAYFREQ = 26, 1134 1135 /// the verbosity of the decomposition based simplex 1136 DECOMP_VERBOSITY = 27, 1137 1138 /// print condition number during the solve 1139 PRINTBASISMETRIC = 28, 1140 1141 /// type of timer for statistics 1142 STATTIMER = 29, 1143 1144 // maximum number of digits for the multiprecision type 1145 MULTIPRECISION_LIMIT = 30, 1146 1147 ///@todo precision-boosting find better parameter name 1148 /// after how many simplex pivots do we store the advanced and stable basis, 1 = every iterations 1149 STORE_BASIS_SIMPLEX_FREQ = 31, 1150 1151 /// number of integer parameters 1152 INTPARAM_COUNT = 32 1153 } IntParam; 1154 1155 /// values for parameter OBJSENSE 1156 enum 1157 { 1158 /// minimization 1159 OBJSENSE_MINIMIZE = -1, 1160 1161 /// maximization 1162 OBJSENSE_MAXIMIZE = 1 1163 }; 1164 1165 /// values for parameter REPRESENTATION 1166 enum 1167 { 1168 /// automatic choice according to number of rows and columns 1169 REPRESENTATION_AUTO = 0, 1170 1171 /// column representation Ax - s = 0, lower <= x <= upper, lhs <= s <= rhs 1172 REPRESENTATION_COLUMN = 1, 1173 1174 /// row representation (lower,lhs) <= (x,Ax) <= (upper,rhs) 1175 REPRESENTATION_ROW = 2 1176 }; 1177 1178 /// values for parameter ALGORITHM 1179 enum 1180 { 1181 /// primal simplex algorithm, i.e., entering for column and leaving for row representation 1182 ALGORITHM_PRIMAL = 0, 1183 1184 /// dual simplex algorithm, i.e., leaving for column and entering for row representation 1185 ALGORITHM_DUAL = 1 1186 }; 1187 1188 /// values for parameter FACTOR_UPDATE_TYPE 1189 enum 1190 { 1191 /// product form update 1192 FACTOR_UPDATE_TYPE_ETA = 0, 1193 1194 /// Forrest-Tomlin type update 1195 FACTOR_UPDATE_TYPE_FT = 1 1196 }; 1197 1198 /// values for parameter VERBOSITY 1199 enum 1200 { 1201 /// only error output 1202 VERBOSITY_ERROR = 0, 1203 1204 /// only error and warning output 1205 VERBOSITY_WARNING = 1, 1206 1207 /// only error, warning, and debug output 1208 VERBOSITY_DEBUG = 2, 1209 1210 /// standard verbosity level 1211 VERBOSITY_NORMAL = 3, 1212 1213 /// high verbosity level 1214 VERBOSITY_HIGH = 4, 1215 1216 /// full verbosity level 1217 VERBOSITY_FULL = 5 1218 }; 1219 1220 /// values for parameter SIMPLIFIER 1221 enum 1222 { 1223 /// disabling presolving 1224 SIMPLIFIER_OFF = 0, 1225 1226 /// using internal presolving methods 1227 SIMPLIFIER_INTERNAL = 3, 1228 1229 /// using the presolve lib papilo 1230 SIMPLIFIER_PAPILO = 2, 1231 1232 /// SoPlex chooses automatically (currently always "internal") 1233 SIMPLIFIER_AUTO = 1 1234 }; 1235 1236 /// values for parameter SCALER 1237 enum 1238 { 1239 /// no scaler 1240 SCALER_OFF = 0, 1241 1242 /// equilibrium scaling on rows or columns 1243 SCALER_UNIEQUI = 1, 1244 1245 /// equilibrium scaling on rows and columns 1246 SCALER_BIEQUI = 2, 1247 1248 /// geometric mean scaling on rows and columns, max 1 round 1249 SCALER_GEO1 = 3, 1250 1251 /// geometric mean scaling on rows and columns, max 8 rounds 1252 SCALER_GEO8 = 4, 1253 1254 /// least square scaling 1255 SCALER_LEASTSQ = 5, 1256 1257 /// geometric mean scaling (max 8 rounds) followed by equilibrium scaling (rows and columns) 1258 SCALER_GEOEQUI = 6 1259 }; 1260 1261 /// values for parameter STARTER 1262 enum 1263 { 1264 /// slack basis 1265 STARTER_OFF = 0, 1266 1267 /// greedy crash basis weighted by objective, bounds, and sides 1268 STARTER_WEIGHT = 1, 1269 1270 /// crash basis from a greedy solution 1271 STARTER_SUM = 2, 1272 1273 /// generic solution-based crash basis 1274 STARTER_VECTOR = 3 1275 }; 1276 1277 /// values for parameter PRICER 1278 enum 1279 { 1280 /// automatic pricer 1281 PRICER_AUTO = 0, 1282 1283 /// Dantzig pricer 1284 PRICER_DANTZIG = 1, 1285 1286 /// partial multiple pricer based on Dantzig pricing 1287 PRICER_PARMULT = 2, 1288 1289 /// devex pricer 1290 PRICER_DEVEX = 3, 1291 1292 /// steepest edge pricer with initialization to unit norms 1293 PRICER_QUICKSTEEP = 4, 1294 1295 /// steepest edge pricer with exact initialization of norms 1296 PRICER_STEEP = 5 1297 }; 1298 1299 /// values for parameter RATIOTESTER 1300 enum 1301 { 1302 /// textbook ratio test without stabilization 1303 RATIOTESTER_TEXTBOOK = 0, 1304 1305 /// standard Harris ratio test 1306 RATIOTESTER_HARRIS = 1, 1307 1308 /// modified Harris ratio test 1309 RATIOTESTER_FAST = 2, 1310 1311 /// bound flipping ratio test for long steps in the dual simplex 1312 RATIOTESTER_BOUNDFLIPPING = 3 1313 }; 1314 1315 /// values for parameter SYNCMODE 1316 enum 1317 { 1318 /// store only real LP 1319 SYNCMODE_ONLYREAL = 0, 1320 1321 /// automatic sync of real and rational LP 1322 SYNCMODE_AUTO = 1, 1323 1324 /// user sync of real and rational LP 1325 SYNCMODE_MANUAL = 2 1326 }; 1327 1328 /// values for parameter READMODE 1329 enum 1330 { 1331 /// standard floating-point parsing 1332 READMODE_REAL = 0, 1333 1334 /// rational parsing 1335 READMODE_RATIONAL = 1 1336 }; 1337 1338 /// values for parameter SOLVEMODE 1339 enum 1340 { 1341 /// apply standard floating-point algorithm 1342 SOLVEMODE_REAL = 0, 1343 1344 /// decide depending on tolerances whether to apply iterative refinement 1345 SOLVEMODE_AUTO = 1, 1346 1347 /// force iterative refinement 1348 SOLVEMODE_RATIONAL = 2 1349 }; 1350 1351 /// values for parameter CHECKMODE 1352 enum 1353 { 1354 /// floating-point check 1355 CHECKMODE_REAL = 0, 1356 1357 /// decide according to READMODE 1358 CHECKMODE_AUTO = 1, 1359 1360 /// rational check 1361 CHECKMODE_RATIONAL = 2 1362 }; 1363 1364 /// values for parameter TIMER 1365 enum 1366 { 1367 /// disable timing 1368 TIMER_OFF = 0, 1369 1370 /// cpu or user time 1371 TIMER_CPU = 1, 1372 1373 /// wallclock time 1374 TIMER_WALLCLOCK = 2 1375 }; 1376 1377 /// values for parameter HYPER_PRICING 1378 enum 1379 { 1380 /// never 1381 HYPER_PRICING_OFF = 0, 1382 1383 /// decide according to problem size 1384 HYPER_PRICING_AUTO = 1, 1385 1386 /// always 1387 HYPER_PRICING_ON = 2 1388 }; 1389 1390 /// values for parameter SOLUTION_POLISHING 1391 enum 1392 { 1393 /// no solution polishing 1394 POLISHING_OFF = 0, 1395 1396 /// maximize number of basic slack variables, i.e. more variables on bounds 1397 POLISHING_INTEGRALITY = 1, 1398 1399 /// minimize number of basic slack variables, i.e. more variables between bounds 1400 POLISHING_FRACTIONALITY = 2 1401 }; 1402 1403 /// real parameters 1404 typedef enum 1405 { 1406 /// primal feasibility tolerance 1407 FEASTOL = 0, 1408 1409 /// dual feasibility tolerance 1410 OPTTOL = 1, 1411 1412 /// general zero tolerance 1413 EPSILON_ZERO = 2, 1414 1415 /// zero tolerance used in factorization 1416 EPSILON_FACTORIZATION = 3, 1417 1418 /// zero tolerance used in update of the factorization 1419 EPSILON_UPDATE = 4, 1420 1421 /// pivot zero tolerance used in factorization 1422 EPSILON_PIVOT = 5, 1423 1424 /// infinity threshold 1425 INFTY = 6, 1426 1427 /// time limit in seconds (INFTY if unlimited) 1428 TIMELIMIT = 7, 1429 1430 /// lower limit on objective value 1431 OBJLIMIT_LOWER = 8, 1432 1433 /// upper limit on objective value 1434 OBJLIMIT_UPPER = 9, 1435 1436 /// working tolerance for feasibility in floating-point solver during iterative refinement 1437 FPFEASTOL = 10, 1438 1439 /// working tolerance for optimality in floating-point solver during iterative refinement 1440 FPOPTTOL = 11, 1441 1442 /// maximum increase of scaling factors between refinements 1443 MAXSCALEINCR = 12, 1444 1445 /// lower threshold in lifting (nonzero matrix coefficients with smaller absolute value will be reformulated) 1446 LIFTMINVAL = 13, 1447 1448 /// upper threshold in lifting (nonzero matrix coefficients with larger absolute value will be reformulated) 1449 LIFTMAXVAL = 14, 1450 1451 /// sparse pricing threshold (\#violations < dimension * SPARSITY_THRESHOLD activates sparse pricing) 1452 SPARSITY_THRESHOLD = 15, 1453 1454 /// threshold on number of rows vs. number of columns for switching from column to row representations in auto mode 1455 REPRESENTATION_SWITCH = 16, 1456 1457 /// geometric frequency at which to apply rational reconstruction 1458 RATREC_FREQ = 17, 1459 1460 /// minimal reduction (sum of removed rows/cols) to continue simplification 1461 MINRED = 18, 1462 1463 /// refactor threshold for nonzeros in last factorized basis matrix compared to updated basis matrix 1464 REFAC_BASIS_NNZ = 19, 1465 1466 /// refactor threshold for fill-in in current factor update compared to fill-in in last factorization 1467 REFAC_UPDATE_FILL = 20, 1468 1469 /// refactor threshold for memory growth in factorization since last refactorization 1470 REFAC_MEM_FACTOR = 21, 1471 1472 /// accuracy of conjugate gradient method in least squares scaling (higher value leads to more iterations) 1473 LEASTSQ_ACRCY = 22, 1474 1475 /// objective offset 1476 OBJ_OFFSET = 23, 1477 1478 /// minimal Markowitz threshold to control sparsity/stability in LU factorization 1479 MIN_MARKOWITZ = 24, 1480 1481 /// minimal modification threshold to apply presolve reductions 1482 SIMPLIFIER_MODIFYROWFAC = 25, 1483 1484 /// factor by which the precision of the floating-point solver is multiplied 1485 PRECISION_BOOSTING_FACTOR = 26, 1486 1487 /// number of real parameters 1488 REALPARAM_COUNT = 27 1489 } RealParam; 1490 1491 #ifdef SOPLEX_WITH_RATIONALPARAM 1492 /// rational parameters 1493 typedef enum 1494 { 1495 /// number of rational parameters 1496 RATIONALPARAM_COUNT = 0 1497 } RationalParam; 1498 #endif 1499 1500 /// class of parameter settings 1501 class Settings 1502 { 1503 public: 1504 static struct BoolParam 1505 { 1506 /// constructor 1507 BoolParam(); 1508 /// array of names for boolean parameters 1509 std::string name[SoPlexBase<R>::BOOLPARAM_COUNT]; 1510 /// array of descriptions for boolean parameters 1511 std::string description[SoPlexBase<R>::BOOLPARAM_COUNT]; 1512 /// array of default values for boolean parameters 1513 bool defaultValue[SoPlexBase<R>::BOOLPARAM_COUNT]; 1514 } boolParam; 1515 1516 static struct IntParam 1517 { 1518 /// constructor 1519 IntParam(); 1520 /// array of names for integer parameters 1521 std::string name[SoPlexBase<R>::INTPARAM_COUNT]; 1522 /// array of descriptions for integer parameters 1523 std::string description[SoPlexBase<R>::INTPARAM_COUNT]; 1524 /// array of default values for integer parameters 1525 int defaultValue[SoPlexBase<R>::INTPARAM_COUNT]; 1526 /// array of lower bounds for int parameter values 1527 int lower[SoPlexBase<R>::INTPARAM_COUNT]; 1528 /// array of upper bounds for int parameter values 1529 int upper[SoPlexBase<R>::INTPARAM_COUNT]; 1530 } intParam; 1531 1532 static struct RealParam 1533 { 1534 /// constructor 1535 RealParam(); 1536 /// array of names for real parameters 1537 std::string name[SoPlexBase<R>::REALPARAM_COUNT]; 1538 /// array of descriptions for real parameters 1539 std::string description[SoPlexBase<R>::REALPARAM_COUNT]; 1540 /// array of default values for real parameters 1541 Real defaultValue[SoPlexBase<R>::REALPARAM_COUNT]; 1542 /// array of lower bounds for real parameter values 1543 Real lower[SoPlexBase<R>::REALPARAM_COUNT]; 1544 /// array of upper bounds for real parameter values 1545 Real upper[SoPlexBase<R>::REALPARAM_COUNT]; 1546 } realParam; 1547 1548 #ifdef SOPLEX_WITH_RATIONALPARAM 1549 static struct RationalParam 1550 { 1551 /// constructor 1552 RationalParam(); 1553 /// array of names for rational parameters 1554 std::string name[SoPlexBase<R>::RATIONALPARAM_COUNT]; 1555 /// array of descriptions for rational parameters 1556 std::string description[SoPlexBase<R>::RATIONALPARAM_COUNT]; 1557 /// array of default values for rational parameters 1558 Rational defaultValue[SoPlexBase<R>::RATIONALPARAM_COUNT]; 1559 /// array of lower bounds for rational parameter values 1560 Rational lower[SoPlexBase<R>::RATIONALPARAM_COUNT]; 1561 /// array of upper bounds for rational parameter values 1562 Rational upper[SoPlexBase<R>::RATIONALPARAM_COUNT]; 1563 } rationalParam; 1564 #endif 1565 1566 /// array of current boolean parameter values 1567 bool _boolParamValues[SoPlexBase<R>::BOOLPARAM_COUNT]; 1568 1569 /// array of current integer parameter values 1570 int _intParamValues[SoPlexBase<R>::INTPARAM_COUNT]; 1571 1572 /// array of current real parameter values 1573 Real _realParamValues[SoPlexBase<R>::REALPARAM_COUNT]; 1574 1575 #ifdef SOPLEX_WITH_RATIONALPARAM 1576 /// array of current rational parameter values 1577 Rational _rationalParamValues[SoPlexBase<R>::RATIONALPARAM_COUNT]; 1578 #endif 1579 1580 /// default constructor initializing default settings 1581 Settings(); 1582 1583 /// copy constructor 1584 Settings(const Settings& settings); 1585 1586 /// assignment operator 1587 Settings& operator=(const Settings& settings); 1588 }; 1589 1590 mutable SPxOut spxout; 1591 1592 /// returns boolean parameter value 1593 bool boolParam(const BoolParam param) const; 1594 1595 /// returns integer parameter value 1596 int intParam(const IntParam param) const; 1597 1598 /// returns real parameter value 1599 Real realParam(const RealParam param) const; 1600 1601 #ifdef SOPLEX_WITH_RATIONALPARAM 1602 /// returns rational parameter value 1603 Rational rationalParam(const RationalParam param) const; 1604 #endif 1605 1606 /// returns current parameter settings 1607 const Settings& settings() const; 1608 1609 /// returns current tolerances 1610 const std::shared_ptr<Tolerances> tolerances() const; 1611 1612 /// sets boolean parameter value; returns true on success 1613 bool setBoolParam(const BoolParam param, const bool value, const bool init = true); 1614 1615 /// sets integer parameter value; returns true on success 1616 bool setIntParam(const IntParam param, const int value, const bool init = true); 1617 1618 /// sets real parameter value; returns true on success 1619 bool setRealParam(const RealParam param, const Real value, const bool init = true); 1620 1621 #ifdef SOPLEX_WITH_RATIONALPARAM 1622 /// sets rational parameter value; returns true on success 1623 bool setRationalParam(const RationalParam param, const Rational value, const bool init = true); 1624 #endif 1625 1626 /// sets parameter settings; returns true on success 1627 bool setSettings(const Settings& newSettings, const bool init = true); 1628 1629 /// resets default parameter settings 1630 void resetSettings(const bool quiet = false, const bool init = true); 1631 1632 /// print non-default parameter values 1633 void printUserSettings(); 1634 1635 /// writes settings file; returns true on success 1636 bool saveSettingsFile(const char* filename, const bool onlyChanged = false, 1637 int solvemode = 1) const; 1638 1639 /// reads settings file; returns true on success 1640 bool loadSettingsFile(const char* filename); 1641 1642 /// parses one setting string and returns true on success; note that string is modified 1643 bool parseSettingsString(char* str); 1644 1645 ///@} 1646 1647 1648 ///@name Statistics 1649 ///@{ 1650 1651 /// set statistic timers to a certain type 1652 void setTimings(const Timer::TYPE ttype); 1653 1654 /// prints solution statistics 1655 void printSolutionStatistics(std::ostream& os); 1656 1657 /// prints statistics on solving process 1658 void printSolvingStatistics(std::ostream& os); 1659 1660 /// prints short statistics 1661 void printShortStatistics(std::ostream& os); 1662 1663 /// prints complete statistics 1664 void printStatistics(std::ostream& os); 1665 1666 /// prints status 1667 1668 void printStatus(std::ostream& os, typename SPxSolverBase<R>::Status status); 1669 1670 ///@} 1671 1672 1673 ///@name Miscellaneous 1674 ///@{ 1675 1676 /// prints version and compilation options 1677 void printVersion() const; 1678 1679 /// checks if real LP and rational LP are in sync; dimensions will always be compared, 1680 /// vector and matrix values only if the respective parameter is set to true. 1681 /// If quiet is set to true the function will only display which vectors are different. 1682 bool areLPsInSync(const bool checkVecVals = true, const bool checkMatVals = false, 1683 const bool quiet = false) const; 1684 1685 /// set the random seeds of the solver instance 1686 void setRandomSeed(unsigned int seed); 1687 1688 /// returns the current random seed of the solver instance 1689 unsigned int randomSeed() const; 1690 1691 ///@} 1692 1693 private: 1694 1695 ///@name Statistics on solving process 1696 ///@{ 1697 1698 /// class of statistics 1699 class Statistics; 1700 1701 /// statistics since last call to solveReal() or solveRational() 1702 Statistics* _statistics; 1703 1704 ///@} 1705 1706 1707 ///@name Parameter settings 1708 ///@{ 1709 1710 Settings* _currentSettings; 1711 1712 std::shared_ptr<Tolerances> _tolerances; 1713 1714 Rational _rationalPosInfty; 1715 Rational _rationalNegInfty; 1716 Rational _rationalFeastol; 1717 Rational _rationalOpttol; 1718 Rational _rationalMaxscaleincr; 1719 1720 ///@} 1721 1722 1723 ///@name Data for the real LP 1724 ///@{ 1725 1726 SPxSolverBase<R> _solver; 1727 SLUFactor<R> _slufactor; 1728 SPxMainSM<R> _simplifierMainSM; 1729 Presol<R> _simplifierPaPILO; 1730 SPxEquiliSC<R> _scalerUniequi; 1731 SPxEquiliSC<R> _scalerBiequi; 1732 SPxGeometSC<R> _scalerGeo1; 1733 SPxGeometSC<R> _scalerGeo8; 1734 SPxGeometSC<R> _scalerGeoequi; 1735 SPxLeastSqSC<R> _scalerLeastsq; 1736 SPxWeightST<R> _starterWeight; 1737 SPxSumST<R> _starterSum; 1738 SPxVectorST<R> _starterVector; 1739 SPxAutoPR<R> _pricerAuto; 1740 SPxDantzigPR<R> _pricerDantzig; 1741 SPxParMultPR<R> _pricerParMult; 1742 SPxDevexPR<R> _pricerDevex; 1743 SPxSteepPR<R> _pricerQuickSteep; 1744 SPxSteepExPR<R> _pricerSteep; 1745 SPxDefaultRT<R> _ratiotesterTextbook; 1746 SPxHarrisRT<R> _ratiotesterHarris; 1747 SPxFastRT<R> _ratiotesterFast; 1748 SPxBoundFlippingRT<R> _ratiotesterBoundFlipping; 1749 1750 SPxLPBase<R>* 1751 _realLP; // the real LP is also used as the original LP for the decomposition dual simplex 1752 SPxLPBase<R>* _decompLP; // used to store the original LP for the decomposition dual simplex 1753 SPxSimplifier<R>* _simplifier; 1754 SPxScaler<R>* _scaler; 1755 SPxStarter<R>* _starter; 1756 1757 #ifdef SOPLEX_WITH_BOOST 1758 #ifdef SOPLEX_WITH_MPFR 1759 //----------------------------- BOOSTED SOLVER ----------------------------- 1760 // multiprecision type used for the boosted solver 1761 using BP = number<mpfr_float_backend<0>, et_off>; 1762 #else 1763 #ifdef SOPLEX_WITH_GMP 1764 using BP = number<gmp_float<50>, et_off>; 1765 #else 1766 using BP = number<cpp_dec_float<50>, et_off>; 1767 #endif 1768 #endif 1769 #else 1770 using BP = double; 1771 #endif 1772 1773 // boosted solver object 1774 SPxSolverBase<BP> _boostedSolver; 1775 1776 // ------------- Main attributes for precision boosting 1777 1778 int _initialPrecision = 50; // initial number of digits for multiprecision 1779 bool _boostingLimitReached; // true if BP::default_precision() > max authorized number of digits 1780 bool _switchedToBoosted; // true if _boostedSolver is used instead of _solver to cope with the numerical failure of _solver 1781 // this attribute remembers wether we are testing feasibility (1), unboundedness (2) or neither (0) 1782 // it is used when storing/loading the right basis in precision boosting. 1783 // example: if _certificateMode == 1, it is the basis for the feasibility LP that should be stored/loaded. 1784 int _certificateMode; 1785 1786 // ------------- Buffers for statistics of precision boosting 1787 1788 // ideally these four attributes would be local variables, however the precision boosting loop 1789 // wraps the solve in a way that it is complicated to declare these variables locally. 1790 int _lastStallPrecBoosts; // number of previous stalling precision boosts 1791 bool _factorSolNewBasisPrecBoost; // false if the current basis has already been factorized (no new iterations have been done) 1792 int _nextRatrecPrecBoost; // the iteration during or after which rational reconstruction can be performed 1793 // buffer storing the number of iterations before a given precision boost 1794 // used to detect stalling (_prevIterations < _statistics->iterations) 1795 int _prevIterations; 1796 1797 // ------------- Tolerances Ratios 1798 1799 /// ratios for computing the tolerances for precision boosting 1800 /// ratio denotes the proportion of precision used by the tolerance 1801 /// e.g. ratio = 0.65, precision = 100 digits, new tol = 10^(0.65*100) 1802 Real _tolPrecisionRatio = 0.65; 1803 Real _epsZeroPrecisionRatio = 1.0; 1804 Real _epsFactorPrecisionRatio = 1.25; 1805 Real _epsUpdatePrecisionRatio = 1.0; 1806 Real _epsPivotPrecisionRatio = 0.625; 1807 1808 // ------------- [Boosted] SLUFactor, Pricers, RatioTesters, Scalers, Simplifiers 1809 1810 SLUFactor<BP> _boostedSlufactor; 1811 1812 SPxAutoPR<BP> _boostedPricerAuto; 1813 SPxDantzigPR<BP> _boostedPricerDantzig; 1814 SPxParMultPR<BP> _boostedPricerParMult; 1815 SPxDevexPR<BP> _boostedPricerDevex; 1816 SPxSteepPR<BP> _boostedPricerQuickSteep; 1817 SPxSteepExPR<BP> _boostedPricerSteep; 1818 1819 SPxDefaultRT<BP> _boostedRatiotesterTextbook; 1820 SPxHarrisRT<BP> _boostedRatiotesterHarris; 1821 SPxFastRT<BP> _boostedRatiotesterFast; 1822 SPxBoundFlippingRT<BP> _boostedRatiotesterBoundFlipping; 1823 1824 SPxScaler<BP>* _boostedScaler; 1825 SPxSimplifier<BP>* _boostedSimplifier; 1826 1827 SPxEquiliSC<BP> _boostedScalerUniequi; 1828 SPxEquiliSC<BP> _boostedScalerBiequi; 1829 SPxGeometSC<BP> _boostedScalerGeo1; 1830 SPxGeometSC<BP> _boostedScalerGeo8; 1831 SPxGeometSC<BP> _boostedScalerGeoequi; 1832 SPxLeastSqSC<BP> _boostedScalerLeastsq; 1833 1834 SPxMainSM<BP> _boostedSimplifierMainSM; 1835 Presol<BP> _boostedSimplifierPaPILO; 1836 1837 //-------------------------------------------------------------------------- 1838 1839 bool _isRealLPLoaded; // true indicates that the original LP is loaded in the _solver variable, hence all actions 1840 // are performed on the original LP. 1841 bool _isRealLPScaled; 1842 bool _applyPolishing; 1843 1844 VectorBase<R> _manualLower; 1845 VectorBase<R> _manualUpper; 1846 VectorBase<R> _manualLhs; 1847 VectorBase<R> _manualRhs; 1848 VectorBase<R> _manualObj; 1849 SPxLPBase<R> _manualRealLP; 1850 1851 ///@} 1852 1853 1854 ///@name Data for the rational LP 1855 ///@{ 1856 1857 SPxLPRational* _rationalLP; 1858 SLUFactorRational _rationalLUSolver; 1859 DataArray<int> _rationalLUSolverBind; 1860 1861 LPColSetRational _slackCols; 1862 VectorRational _unboundedLower; 1863 VectorRational _unboundedUpper; 1864 VectorRational _unboundedLhs; 1865 VectorRational _unboundedRhs; 1866 DSVectorRational _tauColVector; 1867 VectorRational _feasObj; 1868 VectorRational _feasLhs; 1869 VectorRational _feasRhs; 1870 VectorRational _feasLower; 1871 VectorRational _feasUpper; 1872 VectorRational _modLower; 1873 VectorRational _modUpper; 1874 VectorRational _modLhs; 1875 VectorRational _modRhs; 1876 VectorRational _modObj; 1877 DSVectorRational _primalDualDiff; 1878 DataArray< typename SPxSolverBase<R>::VarStatus > _storedBasisStatusRows; 1879 DataArray< typename SPxSolverBase<R>::VarStatus > _storedBasisStatusCols; 1880 Array< UnitVectorRational* > _unitMatrixRational; 1881 bool _storedBasis; 1882 int _beforeLiftRows; 1883 int _beforeLiftCols; 1884 1885 /// type of bounds and sides 1886 typedef enum 1887 { 1888 /// both bounds are infinite 1889 RANGETYPE_FREE = 0, 1890 1891 /// lower bound is finite, upper bound is infinite 1892 RANGETYPE_LOWER = 1, 1893 1894 /// upper bound is finite, lower bound is infinite 1895 RANGETYPE_UPPER = 2, 1896 1897 /// lower and upper bound finite, but different 1898 RANGETYPE_BOXED = 3, 1899 1900 /// lower bound equals upper bound 1901 RANGETYPE_FIXED = 4 1902 } RangeType; 1903 1904 DataArray< RangeType > _colTypes; 1905 DataArray< RangeType > _rowTypes; 1906 1907 ///@} 1908 1909 1910 ///@name Data for the Decomposition Based Dual Simplex 1911 ///@{ 1912 1913 /** row violation structure 1914 */ 1915 struct RowViolation 1916 { 1917 R violation; /**< the violation of the row */ 1918 int idx; /**< index of corresponding row */ 1919 }; 1920 1921 /** Compare class for row violations 1922 */ 1923 struct RowViolationCompare 1924 { 1925 public: 1926 /** constructor 1927 */ 1928 RowViolationCompare() 1929 : entry(0) 1930 { 1931 } 1932 1933 const RowViolation* entry; 1934 1935 R operator()( 1936 RowViolation i, 1937 RowViolation j 1938 ) const 1939 { 1940 return i.violation - j.violation; 1941 } 1942 }; 1943 1944 1945 typedef enum 1946 { 1947 // is the original problem currently being solved. 1948 DECOMP_ORIG = 0, 1949 1950 // is the reduced problem currently being solved. 1951 DECOMP_RED = 1, 1952 1953 // is the complementary problem currently being solved. 1954 DECOMP_COMP = 2 1955 } decompStatus; 1956 1957 // the expected sign of the dual variables. 1958 enum DualSign 1959 { 1960 IS_POS = 0, 1961 IS_NEG = 1, 1962 IS_FREE = 2 1963 }; 1964 1965 SPxSolverBase<R> 1966 _compSolver; // adding a solver to contain the complementary problem. It is too confusing to switch 1967 // the LP for the reduced and complementary problem in the one solver variable. The reduced 1968 // problem will be stored in _solver and the complementary problem will be stored in 1969 // _compSolver. 1970 SLUFactor<R> _compSlufactor; 1971 1972 SPxBasisBase<R> 1973 _decompTransBasis; // the basis required for the transformation to form the reduced problem 1974 1975 VectorBase<R> _transformedObj; // the objective coefficients of the transformed problem 1976 VectorBase<R> _decompFeasVector; // feasibility vector calculated using unshifted bounds. 1977 LPRowSetBase<R> 1978 _transformedRows; // a set of the original rows that have been transformed using the original basis. 1979 SPxColId _compSlackColId; // column id of the primal complementary problem slack column. 1980 SPxRowId _compSlackDualRowId; // row id in the dual of complementary problem related to the slack column. 1981 bool* _decompReducedProbRows; // flag to indicate the inclusion of a row in the reduced problem. 1982 bool* _decompReducedProbCols; // flag to indicate the inclusion of a col in the reduced problem. 1983 int* _decompRowStatus; 1984 int* _decompColStatus; 1985 int* _decompCompProbColIDsIdx; // the index to _decompPrimalColIDs for a given original col. 1986 DataArray < SPxRowId > 1987 _decompReducedProbRowIDs; // the row IDs for the related rows in the reduced problem 1988 DataArray < SPxRowId > 1989 _decompReducedProbColRowIDs;// the row IDs for the related cols in the reduced problem 1990 DataArray < SPxColId > 1991 _decompReducedProbColIDs; // the col IDs for the related cols in the reduced problem 1992 DataArray < SPxRowId > _decompPrimalRowIDs; // the primal row IDs from the original problem 1993 DataArray < SPxColId > _decompPrimalColIDs; // the primal col IDs from the original problem 1994 DataArray < SPxRowId > 1995 _decompElimPrimalRowIDs; // the primal row IDs eliminated in the complementary problem 1996 DataArray < SPxRowId > 1997 _decompDualRowIDs; // the dual row IDs from the complementary problem 1998 DataArray < SPxColId > 1999 _decompDualColIDs; // the dual col IDs from the complementary problem 2000 DataArray < SPxColId > _decompFixedVarDualIDs; // the column ids related to the fixed variables. 2001 DataArray < SPxColId > 2002 _decompVarBoundDualIDs; // the column ids related to the variable bound constraints. 2003 2004 DataArray < SPxColId > 2005 _decompCompPrimalFixedVarIDs; // the column ids related to the fixed variables in the complementary primal. 2006 DataArray < SPxColId > 2007 _decompCompPrimalVarBoundIDs; // the column ids related to the variable bound constraints in the complementary primal. 2008 2009 DataArray < SPxRowId > 2010 _decompCompPrimalRowIDs; // the primal row IDs from the complementary problem 2011 DataArray < SPxColId > 2012 _decompCompPrimalColIDs; // the primal col IDs from the complementary problem 2013 2014 int _nDecompViolBounds; // the number of violated bound constraints 2015 int _nDecompViolRows; // the number of violated rows 2016 int* _decompViolatedBounds; // the violated bounds given by the solution to the IDS reduced problem 2017 int* _decompViolatedRows; // the violated rows given by the solution to the IDS reduced problem 2018 2019 2020 int* _fixedOrigVars; // the original variables that are at their bounds in the reduced problem. 2021 // 1: fixed to upper, -1: fixed to lower, 0: unfixed. 2022 int _nPrimalRows; // the number of original problem rows included in the complementary problem 2023 int _nPrimalCols; // the number of original problem columns included in the complementary problem 2024 int _nElimPrimalRows; // the number of primal rows from the original problem eliminated from the complementary prob 2025 int _nDualRows; // the number of dual rows in the complementary problem. NOTE: _nPrimalRows = _nDualCols 2026 int _nDualCols; // the number of dual columns in the complementary problem. NOTE: _nPrimalRows = _nDualCols 2027 int _nCompPrimalRows; // the number of rows in the complementary primal problem. NOTE: _nPrimalRows = _nCompPrimalRows 2028 int _nCompPrimalCols; // the number of dual columns in the complementary problem. NOTE: _nPrimalCols = _nCompPrimalCols 2029 2030 int _decompDisplayLine; // the count for the display line 2031 2032 NameSet* _rowNames; // the row names from the input file 2033 NameSet* _colNames; // the col names from the input file 2034 2035 // Statistic information 2036 int numIncludedRows; // the number of rows currently included in the reduced problem. 2037 int numDecompIter; // the number of iterations of the decomposition dual simplex algorithm. 2038 int numRedProbIter; // the number of simplex iterations performed in the reduced problem. 2039 int numCompProbIter; // the number of iterations of the complementary problem. 2040 2041 // problem statistics 2042 int numProbRows; 2043 int numProbCols; 2044 int nNonzeros; 2045 R minAbsNonzero; 2046 R maxAbsNonzero; 2047 2048 int origCountLower; 2049 int origCountUpper; 2050 int origCountBoxed; 2051 int origCountFreeCol; 2052 2053 int origCountEqual; 2054 int origCountLhs; 2055 int origCountRhs; 2056 int origCountRanged; 2057 int origCountFreeRow; 2058 2059 decompStatus _currentProb; 2060 2061 ///@} 2062 2063 2064 ///@name Solution data 2065 ///@{ 2066 2067 typename SPxSolverBase<R>::Status _status; 2068 int _lastSolveMode; 2069 2070 DataArray<typename SPxSolverBase<R>::VarStatus > _basisStatusRows; 2071 DataArray<typename SPxSolverBase<R>::VarStatus > _basisStatusCols; 2072 2073 // indicates wether an old basis is currently stored for warm start 2074 bool _hasOldBasis; 2075 bool _hasOldFeasBasis; // basis for testing feasibility 2076 bool _hasOldUnbdBasis; // basis for testing unboundedness 2077 2078 // these vectors store the last basis met in precision boosting when not testing feasibility or unboundedness. 2079 DataArray<typename SPxSolverBase<R>::VarStatus > _oldBasisStatusRows; 2080 DataArray<typename SPxSolverBase<R>::VarStatus > _oldBasisStatusCols; 2081 2082 // these vectors store the last basis met when testing feasibility in precision boosting. 2083 DataArray<typename SPxSolverBase<R>::VarStatus > _oldFeasBasisStatusRows; 2084 DataArray<typename SPxSolverBase<R>::VarStatus > _oldFeasBasisStatusCols; 2085 2086 // these vectors store the last basis met when testing unboundedness in precision boosting. 2087 DataArray<typename SPxSolverBase<R>::VarStatus > _oldUnbdBasisStatusRows; 2088 DataArray<typename SPxSolverBase<R>::VarStatus > _oldUnbdBasisStatusCols; 2089 2090 // these vectors don't replace _basisStatusRows and _basisStatusCols 2091 // they aim to overcome the issue of having the enum VarStatus inside SPxSolverBase. 2092 // When calling setBasis or getBasis (from SPxSolverBase class), a specific conversion is needed. 2093 // Function: SPxSolverBase<BP>::setBasis(...) 2094 // Usage: copy _basisStatusRows(Cols) to _tmpBasisStatusRows(Cols) before calling 2095 // mysolver.setBasis(_tmpBasisStatusRows, _tmpBasisStatusCols) 2096 // Function: SPxSolverBase<BP>::getBasis(...) 2097 // Usage: copy _tmpBasisStatusRows(Cols) to _basisStatusRows(Cols) after calling 2098 // mysolver.getBasis(_tmpBasisStatusRows, _tmpBasisStatusCols, _basisStatusRows.size(), _basisStatusCols.size()) 2099 DataArray<typename SPxSolverBase<BP>::VarStatus > _tmpBasisStatusRows; 2100 DataArray<typename SPxSolverBase<BP>::VarStatus > _tmpBasisStatusCols; 2101 2102 SolBase<R> _solReal; 2103 SolRational _solRational; 2104 SolRational _workSol; 2105 2106 bool _hasBasis; 2107 bool _hasSolReal; 2108 bool _hasSolRational; 2109 2110 ///@} 2111 2112 ///@name Miscellaneous 2113 ///@{ 2114 2115 int _optimizeCalls; 2116 int _unscaleCalls; 2117 2118 Rational _rationalPosone; 2119 Rational _rationalNegone; 2120 Rational _rationalZero; 2121 2122 ///@} 2123 2124 ///@name Constant helper methods 2125 ///@{ 2126 2127 /// extends sparse vector to hold newmax entries if and only if it holds no more free entries 2128 void _ensureDSVectorRationalMemory(DSVectorRational& vec, const int newmax) const; 2129 2130 /// creates a permutation for removing rows/columns from an array of indices 2131 void _idxToPerm(int* idx, int idxSize, int* perm, int permSize) const; 2132 2133 /// creates a permutation for removing rows/columns from a range of indices 2134 void _rangeToPerm(int start, int end, int* perm, int permSize) const; 2135 2136 /// checks consistency for the boosted solver 2137 bool _isBoostedConsistent() const; 2138 2139 /// checks consistency 2140 bool _isConsistent() const; 2141 2142 /// should solving process be stopped? 2143 bool _isSolveStopped(bool& stoppedTime, bool& stoppedIter) const; 2144 2145 /// determines RangeType from real bounds 2146 RangeType _rangeTypeReal(const R& lower, const R& upper) const; 2147 2148 /// determines RangeType from rational bounds 2149 RangeType _rangeTypeRational(const Rational& lower, const Rational& upper) const; 2150 2151 /// switches RANGETYPE_LOWER to RANGETYPE_UPPER and vice versa 2152 RangeType _switchRangeType(const RangeType& rangeType) const; 2153 2154 /// checks whether RangeType corresponds to finite lower bound 2155 bool _lowerFinite(const RangeType& rangeType) const; 2156 2157 /// checks whether RangeType corresponds to finite upper bound 2158 bool _upperFinite(const RangeType& rangeType) const; 2159 2160 ///@} 2161 2162 2163 ///@name Non-constant helper methods 2164 ///@{ 2165 2166 /// adds a single row to the real LP and adjusts basis 2167 void _addRowReal(const LPRowBase<R>& lprow); 2168 2169 /// adds a single row to the real LP and adjusts basis 2170 void _addRowReal(R lhs, const SVectorBase<R>& lprow, R rhs); 2171 2172 /// adds multiple rows to the real LP and adjusts basis 2173 void _addRowsReal(const LPRowSetBase<R>& lprowset); 2174 2175 /// adds a single column to the real LP and adjusts basis 2176 void _addColReal(const LPColReal& lpcol); 2177 2178 /// adds a single column to the real LP and adjusts basis 2179 void _addColReal(R obj, R lower, const SVectorBase<R>& lpcol, R upper); 2180 2181 /// adds multiple columns to the real LP and adjusts basis 2182 void _addColsReal(const LPColSetReal& lpcolset); 2183 2184 /// replaces row \p i with \p lprow and adjusts basis 2185 void _changeRowReal(int i, const LPRowBase<R>& lprow); 2186 2187 /// changes left-hand side vector for constraints to \p lhs and adjusts basis 2188 void _changeLhsReal(const VectorBase<R>& lhs); 2189 2190 /// changes left-hand side of row \p i to \p lhs and adjusts basis 2191 void _changeLhsReal(int i, const R& lhs); 2192 2193 /// changes right-hand side vector to \p rhs and adjusts basis 2194 void _changeRhsReal(const VectorBase<R>& rhs); 2195 2196 /// changes right-hand side of row \p i to \p rhs and adjusts basis 2197 void _changeRhsReal(int i, const R& rhs); 2198 2199 /// changes left- and right-hand side vectors and adjusts basis 2200 void _changeRangeReal(const VectorBase<R>& lhs, const VectorBase<R>& rhs); 2201 2202 /// changes left- and right-hand side of row \p i and adjusts basis 2203 void _changeRangeReal(int i, const R& lhs, const R& rhs); 2204 2205 /// replaces column \p i with \p lpcol and adjusts basis 2206 void _changeColReal(int i, const LPColReal& lpcol); 2207 2208 /// changes vector of lower bounds to \p lower and adjusts basis 2209 void _changeLowerReal(const VectorBase<R>& lower); 2210 2211 /// changes lower bound of column i to \p lower and adjusts basis 2212 void _changeLowerReal(int i, const R& lower); 2213 2214 /// changes vector of upper bounds to \p upper and adjusts basis 2215 void _changeUpperReal(const VectorBase<R>& upper); 2216 2217 /// changes \p i 'th upper bound to \p upper and adjusts basis 2218 void _changeUpperReal(int i, const R& upper); 2219 2220 /// changes vectors of column bounds to \p lower and \p upper and adjusts basis 2221 void _changeBoundsReal(const VectorBase<R>& lower, const VectorBase<R>& upper); 2222 2223 /// changes bounds of column \p i to \p lower and \p upper and adjusts basis 2224 void _changeBoundsReal(int i, const R& lower, const R& upper); 2225 2226 /// changes matrix entry in row \p i and column \p j to \p val and adjusts basis 2227 void _changeElementReal(int i, int j, const R& val); 2228 2229 /// removes row \p i and adjusts basis 2230 void _removeRowReal(int i); 2231 2232 /// removes all rows with an index \p i such that \p perm[i] < 0; upon completion, \p perm[i] >= 0 indicates the 2233 /// new index where row \p i has been moved to; note that \p perm must point to an array of size at least 2234 /// #numRows() 2235 void _removeRowsReal(int perm[]); 2236 2237 /// remove all rows with indices in array \p idx of size \p n; an array \p perm of size #numRows() may be passed 2238 /// as buffer memory 2239 void _removeRowsReal(int idx[], int n, int perm[]); 2240 2241 /// removes rows \p start to \p end including both; an array \p perm of size #numRows() may be passed as buffer 2242 /// memory 2243 void _removeRowRangeReal(int start, int end, int perm[]); 2244 2245 /// removes column i 2246 void _removeColReal(int i); 2247 2248 /// removes all columns with an index \p i such that \p perm[i] < 0; upon completion, \p perm[i] >= 0 indicates the 2249 /// new index where column \p i has been moved to; note that \p perm must point to an array of size at least 2250 /// #numColsReal() 2251 void _removeColsReal(int perm[]); 2252 2253 /// remove all columns with indices in array \p idx of size \p n; an array \p perm of size #numColsReal() may be 2254 /// passed as buffer memory 2255 void _removeColsReal(int idx[], int n, int perm[]); 2256 2257 /// removes columns \p start to \p end including both; an array \p perm of size #numColsReal() may be passed as 2258 /// buffer memory 2259 void _removeColRangeReal(int start, int end, int perm[]); 2260 2261 /// invalidates solution 2262 void _invalidateSolution(); 2263 2264 /// enables simplifier and scaler according to current parameters 2265 void _enableSimplifierAndScaler(); 2266 2267 /// disables simplifier and scaler 2268 void _disableSimplifierAndScaler(); 2269 2270 /// ensures that the rational LP is available; performs no sync 2271 void _ensureRationalLP(); 2272 2273 /// ensures that the real LP and the basis are loaded in the solver; performs no sync 2274 void _ensureRealLPLoaded(); 2275 2276 /// call floating-point solver and update statistics on iterations etc. 2277 void _solveBoostedRealLPAndRecordStatistics(volatile bool* interrupt = NULL); 2278 2279 /// call floating-point solver and update statistics on iterations etc. 2280 void _solveRealLPAndRecordStatistics(volatile bool* interrupt = NULL); 2281 2282 /// reads real LP in LP or MPS format from file and returns true on success; gets row names, column names, and 2283 /// integer variables if desired 2284 bool _readFileReal(const char* filename, NameSet* rowNames = 0, NameSet* colNames = 0, 2285 DIdxSet* intVars = 0); 2286 2287 /// reads rational LP in LP or MPS format from file and returns true on success; gets row names, column names, and 2288 /// integer variables if desired 2289 bool _readFileRational(const char* filename, NameSet* rowNames = 0, NameSet* colNames = 0, 2290 DIdxSet* intVars = 0); 2291 2292 /// completes range type arrays after adding columns and/or rows 2293 void _completeRangeTypesRational(); 2294 2295 /// recomputes range types from scratch using real LP 2296 void _recomputeRangeTypesReal(); 2297 2298 /// recomputes range types from scratch using rational LP 2299 void _recomputeRangeTypesRational(); 2300 2301 /// synchronizes real LP with rational LP, i.e., copies (rounded) rational LP into real LP, without looking at the sync mode 2302 void _syncLPReal(bool time = true); 2303 2304 /// synchronizes rational LP with real LP, i.e., copies real LP to rational LP, without looking at the sync mode 2305 void _syncLPRational(bool time = true); 2306 2307 /// synchronizes rational solution with real solution, i.e., copies (rounded) rational solution to real solution 2308 void _syncRealSolution(); 2309 2310 /// synchronizes real solution with rational solution, i.e., copies real solution to rational solution 2311 void _syncRationalSolution(); 2312 2313 /// returns pointer to a constant unit vector available until destruction of the SoPlexBase class 2314 const UnitVectorRational* _unitVectorRational(const int i); 2315 2316 /// parses one line in a settings file and returns true on success; note that the string is modified 2317 bool _parseSettingsLine(char* line, const int lineNumber); 2318 2319 ///@} 2320 2321 2322 //**@name Private solving methods implemented in solverational.hpp */ 2323 ///@{ 2324 2325 /// stores floating-point solution of original LP as current rational solution and ensure that solution vectors have right dimension; ensure that solution is aligned with basis 2326 template <typename T> 2327 void _storeRealSolutionAsRational( 2328 SolRational& sol, 2329 VectorBase<T>& primalReal, 2330 VectorBase<T>& dualReal, 2331 int& dualSize); 2332 2333 /// computes violation of bounds during the refinement loop 2334 void _computeBoundsViolation(SolRational& sol, Rational& boundsViolation); 2335 2336 /// computes violation of sides during the refinement loop 2337 void _computeSidesViolation(SolRational& sol, Rational& sideViolation); 2338 2339 /// computes violation of reduced costs during the refinement loop 2340 void _computeReducedCostViolation( 2341 SolRational& sol, 2342 Rational& redCostViolation, 2343 const bool& maximizing); 2344 2345 /// computes dual violation during the refinement loop 2346 void _computeDualViolation( 2347 SolRational& sol, 2348 Rational& dualViolation, 2349 const bool& maximizing); 2350 2351 /// checks termination criteria for refinement loop 2352 bool _isRefinementOver( 2353 bool& primalFeasible, 2354 bool& dualFeasible, 2355 Rational& boundsViolation, 2356 Rational& sideViolation, 2357 Rational& redCostViolation, 2358 Rational& dualViolation, 2359 int minIRRoundsRemaining, 2360 bool& stoppedTime, 2361 bool& stoppedIter, 2362 int numFailedRefinements); 2363 2364 /// checks refinement loop progress 2365 void _checkRefinementProgress( 2366 Rational& boundsViolation, 2367 Rational& sideViolation, 2368 Rational& redCostViolation, 2369 Rational& dualViolation, 2370 Rational& maxViolation, 2371 Rational& bestViolation, 2372 const Rational& violationImprovementFactor, 2373 int& numFailedRefinements); 2374 2375 /// performs rational reconstruction and/or factorizationd 2376 void _ratrecAndOrRatfac( 2377 int& minIRRoundsRemaining, 2378 int& lastStallIterations, 2379 int& numberOfIterations, 2380 bool& factorSolNewBasis, 2381 int& nextRatrec, 2382 const Rational& errorCorrectionFactor, 2383 Rational& errorCorrection, 2384 Rational& maxViolation, 2385 SolRational& sol, 2386 bool& primalFeasible, 2387 bool& dualFeasible, 2388 bool& stoppedTime, 2389 bool& stoppedIter, 2390 bool& error, 2391 bool& breakAfter, 2392 bool& continueAfter); 2393 2394 /// forces value of given nonbasic variable to bound 2395 void _forceNonbasicToBound( 2396 SolRational& sol, 2397 int& c, 2398 const int& maxDimRational, 2399 bool toLower); 2400 2401 /// computes primal scaling factor; limit increase in scaling by tolerance used in floating point solve 2402 void _computePrimalScalingFactor( 2403 Rational& maxScale, 2404 Rational& primalScale, 2405 Rational& boundsViolation, 2406 Rational& sideViolation, 2407 Rational& redCostViolation); 2408 2409 /// computes dual scaling factor; limit increase in scaling by tolerance used in floating point solve 2410 void _computeDualScalingFactor( 2411 Rational& maxScale, 2412 Rational& primalScale, 2413 Rational& dualScale, 2414 Rational& redCostViolation, 2415 Rational& dualViolation); 2416 2417 /// applies scaled bounds 2418 template <typename T> 2419 void _applyScaledBounds(SPxSolverBase<T>& solver, Rational& primalScale); 2420 2421 /// applies scaled sides 2422 template <typename T> 2423 void _applyScaledSides(SPxSolverBase<T>& solver, Rational& primalScale); 2424 2425 /// applies scaled objective function 2426 template <typename T> 2427 void _applyScaledObj(SPxSolverBase<T>& solver, Rational& dualScale, SolRational& sol); 2428 2429 /// evaluates result of solve. Return true if the algorithm must to stopped, false otherwise. 2430 template <typename T> 2431 bool _evaluateResult( 2432 SPxSolverBase<T>& solver, 2433 typename SPxSolverBase<T>::Status result, 2434 bool usingRefinedLP, 2435 SolRational& sol, 2436 VectorBase<T>& dualReal, 2437 bool& infeasible, 2438 bool& unbounded, 2439 bool& stoppedTime, 2440 bool& stoppedIter, 2441 bool& error); 2442 2443 /// corrects primal solution and aligns with basis 2444 template <typename T> 2445 void _correctPrimalSolution( 2446 SolRational& sol, 2447 Rational& primalScale, 2448 int& primalSize, 2449 const int& maxDimRational, 2450 VectorBase<T>& primalReal); 2451 2452 /// updates or recomputes slacks depending on which looks faster 2453 void _updateSlacks(SolRational& sol, int& primalSize); 2454 2455 /// corrects dual solution and aligns with basis 2456 template <typename T> 2457 void _correctDualSolution( 2458 SPxSolverBase<T>& solver, 2459 SolRational& sol, 2460 const bool& maximizing, 2461 VectorBase<T>& dualReal, 2462 Rational& dualScale, 2463 int& dualSize, 2464 const int& maxDimRational); 2465 2466 /// updates or recomputes reduced cost values depending on which looks faster; adding one to the length of the 2467 /// dual vector accounts for the objective function vector 2468 void _updateReducedCosts(SolRational& sol, int& dualSize, const int& numCorrectedPrimals); 2469 2470 ///@todo precision-boosting move some place else 2471 /// converts the given DataArray of VarStatus to boostedPrecision 2472 void _convertDataArrayVarStatusToBoosted( 2473 DataArray< typename SPxSolverBase<R>::VarStatus >& base, 2474 DataArray< typename SPxSolverBase<BP>::VarStatus >& copy); 2475 2476 ///@todo precision-boosting move some place else 2477 /// converts the given DataArray of VarStatus to R precision 2478 void _convertDataArrayVarStatusToRPrecision( 2479 DataArray< typename SPxSolverBase<BP>::VarStatus >& base, 2480 DataArray< typename SPxSolverBase<R>::VarStatus >& copy); 2481 2482 /// disable initial precision solver and switch to boosted solver 2483 void _switchToBoosted(); 2484 2485 /// setup boosted solver before launching iteration 2486 void _setupBoostedSolver(); 2487 2488 /// increase the multiprecision, return false if maximum precision is reached, true otherwise 2489 bool _boostPrecision(); 2490 2491 /// reset the boosted precision to the default value 2492 void _resetBoostedPrecision(); 2493 2494 /// setup recovery mecanism using multiprecision, return false if maximum precision reached, true otherwise 2495 bool _setupBoostedSolverAfterRecovery(); 2496 2497 /// return true if slack basis has to be loaded for boosted solver 2498 bool _isBoostedStartingFromSlack(bool initialSolve = true); 2499 2500 /// indicate if we are testing feasibility, unboundedness or neither 2501 void _switchToStandardMode(); 2502 void _switchToFeasMode(); 2503 void _switchToUnbdMode(); 2504 2505 /// check if we are testing feasibility, unboundedness or neither 2506 bool _inStandardMode(); 2507 bool _inFeasMode(); 2508 bool _inUnbdMode(); 2509 2510 // stores given basis in old basis attributes: _oldBasisStatusRows, _oldFeasBasisStatusRows, _oldUnbdBasisStatusRows (and ...Cols) 2511 void _storeBasisAsOldBasis(DataArray< typename SPxSolverBase<R>::VarStatus >& rows, 2512 DataArray< typename SPxSolverBase<R>::VarStatus >& cols); 2513 2514 // stores given basis in old basis attributes: _oldBasisStatusRows, _oldFeasBasisStatusRows, _oldUnbdBasisStatusRows (and ...Cols) 2515 void _storeBasisAsOldBasisBoosted(DataArray< typename SPxSolverBase<BP>::VarStatus >& rows, 2516 DataArray< typename SPxSolverBase<BP>::VarStatus >& cols); 2517 2518 // get the last advanced and stable basis stored by the initial solver and store it as old basis, unsimplify basis if simplifier activated 2519 void _storeLastStableBasis(bool vanished); 2520 2521 // get the last advanced and stable basis stored by the boosted solver and store it as old basis, unsimplify basis if simplifier activated 2522 void _storeLastStableBasisBoosted(bool vanished); 2523 2524 // load old basis in solver. The old basis loaded depends on the certificate mode (feasibility, unboundedness, or neither) 2525 bool _loadBasisFromOldBasis(bool boosted); 2526 2527 // update statistics for precision boosting 2528 void _updateBoostingStatistics(); 2529 2530 /// solves current problem using multiprecision floating-point solver 2531 /// return false if a new boosted iteration is necessary, true otherwise 2532 void _solveRealForRationalBoostedStable( 2533 SolRational& sol, 2534 bool& primalFeasible, 2535 bool& dualFeasible, 2536 bool& infeasible, 2537 bool& unbounded, 2538 bool& stoppedTime, 2539 bool& stoppedIter, 2540 bool& error, 2541 bool& needNewBoostedIt); 2542 2543 /// solves current problem with iterative refinement and recovery mechanism using boosted solver 2544 void _performOptIRStableBoosted( 2545 SolRational& sol, 2546 bool acceptUnbounded, 2547 bool acceptInfeasible, 2548 int minIRRoundsRemaining, 2549 bool& primalFeasible, 2550 bool& dualFeasible, 2551 bool& infeasible, 2552 bool& unbounded, 2553 bool& stoppedTime, 2554 bool& stoppedIter, 2555 bool& error, 2556 bool& needNewBoostedIt); 2557 2558 /// perform iterative refinement using the right precision 2559 void _performOptIRWrapper( 2560 SolRational& sol, 2561 bool acceptUnbounded, 2562 bool acceptInfeasible, 2563 int minIRRoundsRemaining, 2564 bool& primalFeasible, 2565 bool& dualFeasible, 2566 bool& infeasible, 2567 bool& unbounded, 2568 bool& stoppedTime, 2569 bool& stoppedIter, 2570 bool& error 2571 ); 2572 2573 /// solves current problem using double floating-point solver 2574 void _solveRealForRationalStable( 2575 SolRational& sol, 2576 bool& primalFeasible, 2577 bool& dualFeasible, 2578 bool& infeasible, 2579 bool& unbounded, 2580 bool& stoppedTime, 2581 bool& stoppedIter, 2582 bool& error); 2583 2584 /// solves current problem with iterative refinement and recovery mechanism 2585 void _performOptIRStable(SolRational& sol, 2586 bool acceptUnbounded, 2587 bool acceptInfeasible, 2588 int minIRRoundsRemaining, 2589 bool& primalFeasible, 2590 bool& dualFeasible, 2591 bool& infeasible, 2592 bool& unbounded, 2593 bool& stoppedTime, 2594 bool& stoppedIter, 2595 bool& error); 2596 2597 /// performs iterative refinement on the auxiliary problem for testing unboundedness 2598 void _performUnboundedIRStable(SolRational& sol, bool& hasUnboundedRay, bool& stoppedTime, 2599 bool& stoppedIter, bool& error); 2600 2601 /// performs iterative refinement on the auxiliary problem for testing feasibility 2602 void _performFeasIRStable(SolRational& sol, bool& withDualFarkas, bool& stoppedTime, 2603 bool& stoppedIter, bool& error); 2604 2605 /// reduces matrix coefficient in absolute value by the lifting procedure of Thiele et al. 2013 2606 void _lift(); 2607 2608 /// undoes lifting 2609 void _project(SolRational& sol); 2610 2611 /// store basis 2612 void _storeBasis(); 2613 2614 /// restore basis 2615 void _restoreBasis(); 2616 2617 /// stores objective, bounds, and sides of real LP 2618 void _storeLPReal(); 2619 2620 /// restores objective, bounds, and sides of real LP 2621 void _restoreLPReal(); 2622 2623 /// introduces slack variables to transform inequality constraints into equations for both rational and real LP, 2624 /// which should be in sync 2625 void _transformEquality(); 2626 2627 /// undoes transformation to equality form 2628 void _untransformEquality(SolRational& sol); 2629 2630 /// transforms LP to unboundedness problem by moving the objective function to the constraints, changing right-hand 2631 /// side and bounds to zero, and adding an auxiliary variable for the decrease in the objective function 2632 void _transformUnbounded(); 2633 2634 /// undoes transformation to unboundedness problem 2635 void _untransformUnbounded(SolRational& sol, bool unbounded); 2636 2637 /// transforms LP to feasibility problem by removing the objective function, shifting variables, and homogenizing the 2638 /// right-hand side 2639 void _transformFeasibility(); 2640 2641 /// undoes transformation to feasibility problem 2642 void _untransformFeasibility(SolRational& sol, bool infeasible); 2643 2644 /** computes radius of infeasibility box implied by an approximate Farkas' proof 2645 2646 Given constraints of the form \f$ lhs <= Ax <= rhs \f$, a farkas proof y should satisfy \f$ y^T A = 0 \f$ and 2647 \f$ y_+^T lhs - y_-^T rhs > 0 \f$, where \f$ y_+, y_- \f$ denote the positive and negative parts of \f$ y \f$. 2648 If \f$ y \f$ is approximate, it may not satisfy \f$ y^T A = 0 \f$ exactly, but the proof is still valid as long 2649 as the following holds for all potentially feasible \f$ x \f$: 2650 2651 \f[ 2652 y^T Ax < (y_+^T lhs - y_-^T rhs) (*) 2653 \f] 2654 2655 we may therefore calculate \f$ y^T A \f$ and \f$ y_+^T lhs - y_-^T rhs \f$ exactly and check if the upper and lower 2656 bounds on \f$ x \f$ imply that all feasible \f$ x \f$ satisfy (*), and if not then compute bounds on \f$ x \f$ to 2657 guarantee (*). The simplest way to do this is to compute 2658 2659 \f[ 2660 B = (y_+^T lhs - y_-^T rhs) / \sum_i(|(y^T A)_i|) 2661 \f] 2662 2663 noting that if every component of \f$ x \f$ has \f$ |x_i| < B \f$, then (*) holds. 2664 2665 \f$ B \f$ can be increased by iteratively including variable bounds smaller than \f$ B \f$. The speed of this 2666 method can be further improved by using interval arithmetic for all computations. For related information see 2667 Sec. 4 of Neumaier and Shcherbina, Mathematical Programming A, 2004. 2668 2669 Set transformed to true if this method is called after _transformFeasibility(). 2670 */ 2671 void _computeInfeasBox(SolRational& sol, bool transformed); 2672 2673 /// solves real LP during iterative refinement 2674 typename SPxSolverBase<R>::Status _solveRealForRational(bool fromscratch, VectorBase<R>& primal, 2675 VectorBase<R>& dual, 2676 DataArray< typename SPxSolverBase<R>::VarStatus >& basisStatusRows, 2677 DataArray< typename SPxSolverBase<R>::VarStatus >& basisStatusCols); 2678 2679 /// solves real LP with recovery mechanism 2680 typename SPxSolverBase<R>::Status _solveRealStable(bool acceptUnbounded, bool acceptInfeasible, 2681 VectorBase<R>& primal, VectorBase<R>& dual, 2682 DataArray< typename SPxSolverBase<R>::VarStatus >& basisStatusRows, 2683 DataArray< typename SPxSolverBase<R>::VarStatus >& basisStatusCols, 2684 const bool forceNoSimplifier = false); 2685 2686 /// solves real LP during iterative refinement 2687 void _solveRealForRationalBoosted( 2688 VectorBase<BP>& primal, VectorBase<BP>& dual, 2689 DataArray< typename SPxSolverBase<R>::VarStatus >& basisStatusRows, 2690 DataArray< typename SPxSolverBase<R>::VarStatus >& basisStatusCols, 2691 typename SPxSolverBase<BP>::Status& boostedResult, bool initialSolve); 2692 2693 /// computes rational inverse of basis matrix as defined by _rationalLUSolverBind 2694 void _computeBasisInverseRational(); 2695 2696 /// factorizes rational basis matrix in column representation 2697 void _factorizeColumnRational(SolRational& sol, 2698 DataArray< typename SPxSolverBase<R>::VarStatus >& basisStatusRows, 2699 DataArray< typename SPxSolverBase<R>::VarStatus >& basisStatusCols, bool& stoppedTime, 2700 bool& stoppedIter, bool& error, bool& optimal); 2701 2702 /// attempts rational reconstruction of primal-dual solution 2703 bool _reconstructSolutionRational(SolRational& sol, 2704 DataArray< typename SPxSolverBase<R>::VarStatus >& basisStatusRows, 2705 DataArray< typename SPxSolverBase<R>::VarStatus >& basisStatusCols, 2706 const Rational& denomBoundSquared); 2707 ///@} 2708 2709 ///@name Private solving methods implemented in solvereal.cpp 2710 ///@{ 2711 2712 /// solves the templated LP 2713 void _optimize(volatile bool* interrupt = NULL); 2714 2715 /// temporary fix for Rational 2716 void _optimizeRational(volatile bool* interrupt = NULL); 2717 2718 /// checks result of the solving process and solves again without preprocessing if necessary 2719 void _evaluateSolutionReal(typename SPxSimplifier<R>::Result simplificationStatus); 2720 2721 /// solves real LP with/without preprocessing 2722 void _preprocessAndSolveReal(bool applyPreprocessing, volatile bool* interrupt = NULL); 2723 2724 /// loads original problem into solver and solves again after it has been solved to optimality with preprocessing 2725 void _resolveWithoutPreprocessing(typename SPxSimplifier<R>::Result simplificationStatus); 2726 2727 /// verify computed solution and resolve if necessary 2728 void _verifySolutionReal(); 2729 2730 /// verify computed obj stop and resolve if necessary 2731 void _verifyObjLimitReal(); 2732 2733 /// stores solution of the real LP; before calling this, the real LP must be loaded in the solver and solved (again) 2734 void _storeSolutionReal(bool verify = true); 2735 2736 /// stores solution from the simplifier because problem vanished in presolving step 2737 void _storeSolutionRealFromPresol(); 2738 2739 /// unscales stored solution to remove internal or external scaling of LP 2740 void _unscaleSolutionReal(SPxLPBase<R>& LP, bool persistent = true); 2741 2742 /// load original LP and possibly setup a slack basis 2743 void _loadRealLP(bool initBasis); 2744 2745 /// check scaling of LP 2746 void _checkScaling(SPxLPBase<R>* origLP) const; 2747 2748 /// check correctness of (un)scaled basis matrix operations 2749 void _checkBasisScaling(); 2750 2751 /// check whether persistent scaling is supposed to be reapplied again after unscaling 2752 bool _reapplyPersistentScaling() const; 2753 2754 /// solves LP using the decomposition based dual simplex 2755 void _solveDecompositionDualSimplex(); 2756 2757 /// creating copies of the original problem that will be manipulated to form the reduced and complementary problems 2758 void _createDecompReducedAndComplementaryProblems(); 2759 2760 /// forms the reduced problem 2761 void _formDecompReducedProblem(bool& stop); 2762 2763 /// solves the reduced problem 2764 void _solveDecompReducedProblem(); 2765 2766 /// forms the complementary problem 2767 void _formDecompComplementaryProblem(); 2768 2769 /// simplifies the problem and solves 2770 void _decompSimplifyAndSolve(SPxSolverBase<R>& solver, SLUFactor<R>& sluFactor, bool fromScratch, 2771 bool applyPreprocessing); 2772 2773 /// loads original problem into solver and solves again after it has been solved to optimality with preprocessing 2774 void _decompResolveWithoutPreprocessing(SPxSolverBase<R>& solver, SLUFactor<R>& sluFactor, 2775 typename SPxSimplifier<R>::Result result); 2776 2777 /// identifies the columns of the row-form basis that correspond to rows with zero dual multipliers. 2778 void _getZeroDualMultiplierIndices(VectorBase<R> feasVector, int* nonposind, int* colsforremoval, 2779 int* nnonposind, bool& stop); 2780 2781 /// retrieves the compatible columns from the constraint matrix 2782 void _getCompatibleColumns(VectorBase<R> feasVector, int* nonposind, int* compatind, 2783 int* rowsforremoval, int* colsforremoval, 2784 int nnonposind, int* ncompatind, bool formRedProb, bool& stop); 2785 2786 /// computes the reduced problem objective coefficients 2787 void _computeReducedProbObjCoeff(bool& stop); 2788 2789 /// computes the compatible bound constraints and adds them to the reduced problem 2790 void _getCompatibleBoundCons(LPRowSetBase<R>& boundcons, int* compatboundcons, int* nonposind, 2791 int* ncompatboundcons, 2792 int nnonposind, bool& stop); 2793 2794 /// computes the rows to remove from the complementary problem 2795 void _getRowsForRemovalComplementaryProblem(int* nonposind, int* bind, int* rowsforremoval, 2796 int* nrowsforremoval, 2797 int nnonposind); 2798 2799 /// removing rows from the complementary problem. 2800 void _deleteAndUpdateRowsComplementaryProblem(SPxRowId rangedRowIds[], int& naddedrows); 2801 2802 /// evaluates the solution of the reduced problem for the DBDS 2803 void _evaluateSolutionDecomp(SPxSolverBase<R>& solver, SLUFactor<R>& sluFactor, 2804 typename SPxSimplifier<R>::Result result); 2805 2806 /// update the reduced problem with additional columns and rows 2807 void _updateDecompReducedProblem(R objVal, VectorBase<R> dualVector, VectorBase<R> redcostVector, 2808 VectorBase<R> compPrimalVector, 2809 VectorBase<R> compDualVector); 2810 2811 /// update the reduced problem with additional columns and rows based upon the violated original bounds and rows 2812 void _updateDecompReducedProblemViol(bool allrows); 2813 2814 /// builds the update rows with those violated in the complmentary problem 2815 void _findViolatedRows(R compObjValue, Array<RowViolation>& violatedrows, int& nviolatedrows); 2816 2817 /// update the dual complementary problem with additional columns and rows 2818 void _updateDecompComplementaryDualProblem(bool origObj); 2819 2820 /// update the primal complementary problem with additional columns and rows 2821 void _updateDecompComplementaryPrimalProblem(bool origObj); 2822 2823 /// checking the optimality of the original problem. 2824 void _checkOriginalProblemOptimality(VectorBase<R> primalVector, bool printViol); 2825 2826 /// updating the slack column coefficients to adjust for equality constraints 2827 void _updateComplementaryDualSlackColCoeff(); 2828 2829 /// updating the slack column coefficients to adjust for equality constraints 2830 void _updateComplementaryPrimalSlackColCoeff(); 2831 2832 /// removing the dual columns related to the fixed variables 2833 void _removeComplementaryDualFixedPrimalVars(int* currFixedVars); 2834 2835 /// removing the dual columns related to the fixed variables 2836 void _identifyComplementaryDualFixedPrimalVars(int* currFixedVars); 2837 2838 /// updating the dual columns related to the fixed primal variables. 2839 void _updateComplementaryDualFixedPrimalVars(int* currFixedVars); 2840 2841 /// removing the dual columns related to the fixed variables 2842 void _identifyComplementaryPrimalFixedPrimalVars(int* currFixedVars); 2843 2844 /// updating the dual columns related to the fixed primal variables. 2845 void _updateComplementaryPrimalFixedPrimalVars(int* currFixedVars); 2846 2847 /// updating the complementary dual problem with the original objective function 2848 void _setComplementaryDualOriginalObjective(); 2849 2850 /// updating the complementary primal problem with the original objective function 2851 void _setComplementaryPrimalOriginalObjective(); 2852 2853 /// determining which bound the primal variables will be fixed to. 2854 int getOrigVarFixedDirection(int colNum); 2855 2856 /// checks the dual feasibility of the current basis 2857 bool checkBasisDualFeasibility(VectorBase<R> feasVec); 2858 2859 /// returns the expected sign of the dual variables for the reduced problem 2860 DualSign getExpectedDualVariableSign(int rowNumber); 2861 2862 /// returns the expected sign of the dual variables for the original problem 2863 DualSign getOrigProbDualVariableSign(int rowNumber); 2864 2865 /// prints a display line of the flying table for the DBDS 2866 void printDecompDisplayLine(SPxSolverBase<R>& solver, const SPxOut::Verbosity origVerb, bool force, 2867 bool forceHead); 2868 2869 /// stores the problem statistics of the original problem 2870 void getOriginalProblemStatistics(); 2871 2872 /// stores the problem statistics of the original problem 2873 void printOriginalProblemStatistics(std::ostream& os); 2874 2875 /// gets the coefficient of the slack variable in the primal complementary problem 2876 R getCompSlackVarCoeff(int primalRowNum); 2877 2878 /// gets violation of bounds; returns true on success 2879 bool getDecompBoundViolation(R& maxviol, R& sumviol); 2880 2881 /// gets violation of constraints; returns true on success 2882 bool getDecompRowViolation(R& maxviol, R& sumviol); 2883 2884 /// function call to terminate the decomposition simplex 2885 bool decompTerminate(R timeLimit); 2886 2887 /// function to build a basis for the original problem as given by the solution to the reduced problem 2888 void _writeOriginalProblemBasis(const char* filename, NameSet* rowNames, NameSet* colNames, 2889 bool cpxFormat); 2890 2891 /// function to retrieve the original problem row basis status from the reduced and complementary problems 2892 void getOriginalProblemBasisRowStatus(DataArray< int >& degenerateRowNums, 2893 DataArray< typename SPxSolverBase<R>::VarStatus >& degenerateRowStatus, int& nDegenerateRows, 2894 int& nNonBasicRows); 2895 2896 /// function to retrieve the column status for the original problem basis from the reduced and complementary problems 2897 void getOriginalProblemBasisColStatus(int& nNonBasicCols); 2898 2899 ///@} 2900 }; 2901 2902 /* Backwards compatibility */ 2903 typedef SoPlexBase<Real> SoPlex; 2904 // A header file containing all the general templated functions 2905 2906 } // namespace soplex 2907 2908 // General templated function 2909 #include "soplex/solvedbds.hpp" 2910 #include "soplex.hpp" 2911 #include "soplex/solverational.hpp" 2912 #include "soplex/testsoplex.hpp" 2913 #include "soplex/solvereal.hpp" 2914 2915 #endif // _SOPLEX_H_ 2916