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 spxbasis.h 26 * @brief Simplex basis. 27 */ 28 #ifndef _SPXBASIS_H_ 29 #define _SPXBASIS_H_ 30 31 /* undefine SOPLEX_DEBUG flag from including files; if SOPLEX_DEBUG should be defined in this file, do so below */ 32 #ifdef SOPLEX_DEBUG 33 #define SOPLEX_DEBUG_SPXBASIS 34 #undef SOPLEX_DEBUG 35 #endif 36 37 #include <assert.h> 38 #include <iostream> 39 #include <iomanip> 40 #include <string.h> 41 #include <sstream> 42 43 #include "soplex/spxdefines.h" 44 #include "soplex/spxlp.h" 45 #include "soplex/svector.h" 46 #include "soplex/ssvector.h" 47 #include "soplex/dataarray.h" 48 #include "soplex/slinsolver.h" 49 #include "soplex/nameset.h" 50 #include "soplex/spxout.h" 51 #include "soplex/timerfactory.h" 52 53 //#define SOPLEX_MEASUREUPDATETIME 54 55 namespace soplex 56 { 57 template <class R> 58 class SPxSolverBase; 59 60 /**@class SPxBasisBase 61 @brief Simplex basis. 62 @ingroup Algo 63 64 Consider the linear program as provided from class SPxLP: 65 \f[ 66 \begin{array}{rl} 67 \hbox{max} & c^T x \\ 68 \hbox{s.t.} & l_r \le Ax \le u_r \\ 69 & l_c \le x \le u_c 70 \end{array} 71 \f] 72 where \f$c, l_c, u_c, x \in {\bf R}^n\f$, \f$l_r, u_r \in {\bf R}^m\f$ and 73 \f$A \in {\bf R}^{m \times n}\f$. Solving this LP with the simplex algorithm 74 requires the definition of a \em basis. Such can be defined as a set of 75 column vectors or a set of row vectors building a non-singular matrix. We 76 will refer to the first case as the \em columnwise \em representation and 77 the latter case will be called the \em rowwise \em representation. In both 78 cases, a \em basis is a set of vectors forming a non-singular matrix. The 79 dimension of the vectors is referred to as the basis' \em dimension, 80 whereas the number of vectors belonging to the LP is called the basis' 81 \em codimension. 82 83 Class SPxBasisBase is designed to represent a generic simplex basis, suitable 84 for both representations. At any time the representation can be changed by 85 calling method setRep(). 86 87 Class SPxBasisBase provides methods for solving linear systems with the basis 88 matrix. However, SPxBasisBase does not provide a linear solver by its own. 89 Instead, a SLinSolver object must be #load%ed to a SPxBasisBase which will 90 be called for solving linear systems. 91 */ 92 template <class R> // theLP gets templated 93 class SPxBasisBase 94 { 95 public: 96 97 /// basis status. 98 /** Each SPxBasisBase is assigned a status flag, which can take on of the 99 above values. 100 */ 101 enum SPxStatus 102 { 103 NO_PROBLEM = -2, ///< No Problem has been loaded to the basis. 104 SINGULAR = -1, ///< Basis is singular. 105 REGULAR = 0, ///< Basis is not known to be dual nor primal feasible. 106 DUAL = 1, ///< Basis is dual feasible. 107 PRIMAL = 2, ///< Basis is primal feasible. 108 OPTIMAL = 3, ///< Basis is optimal, i.e. dual and primal feasible. 109 UNBOUNDED = 4, ///< LP has been proven to be primal unbounded. 110 INFEASIBLE = 5 ///< LP has been proven to be primal infeasible. 111 }; 112 113 114 /// Basis descriptor. 115 class Desc 116 { 117 public: 118 119 //------------------------------------ 120 ///@name Status 121 ///@{ 122 /// Status of a variable. 123 /** A basis is described by assigning a Status to all of the LP 124 variables and covariables. This assignment is maintained by the 125 basis #Desc%riptor. 126 127 Variables and covariables (slackvariables) may have a primal or dual Status. The 128 first type specifies that a variable is set on a primal bound, while 129 the latter type indicates a dual variable to be set on a bound. 130 If a row variable has a primal status, say #P_ON_UPPER, this means 131 that the upper bound of the inequality is set to be tight. Hence, 132 in this case the upper bound must not be infinity. 133 134 Equivalently, if the status of a variable is dual, say #D_ON_UPPER, 135 it means that the dual variable corresponding to the upper bound 136 inequality of this variable is set to 0. 137 138 For a column basis, primal #Status%es correspond to nonbasic 139 variables, while dual ones are basic. This is reversed for a row 140 basis. We will now reveal in more detail the significance of 141 variable #Status%es. 142 143 <b>Primal Variables</b> 144 145 Consider a range inequality \f$l_r \le a^T x \le u_r\f$ or bounds on 146 a variable \f$l_c \le x_c \le u_c\f$. The following table reveals 147 what is implied if the corresponding variable or covariable is 148 assigned to a primal #Status: 149 150 \f[ 151 \begin{array}{lcl} 152 l_c \le x_c \le u_c & \mbox{Status}(x_i) & l_r \le a^T x \le u_r \\ 153 \hline 154 x_c = u_c < \infty & \mbox{P\_ON\_UPPER} & a^T x = u_r < \infty \\ 155 x_c = l_c > -\infty & \mbox{P\_ON\_LOWER} & a^T x = l_r > -\infty \\ 156 -\infty < l_c = x_c = u_c < \infty 157 & \mbox{P\_FIXED} & 158 -\infty < l_r = a^T x = u_r < \infty \\ 159 -\infty = l_i < x_i=0 < u_i = \infty 160 & \mbox{P\_FREE} & 161 -\infty = l_r < a^T x = 0 < u_r = \infty \\ 162 \end{array} 163 \f] 164 165 Note that to determine whether a variable with #Status stat is set to 166 its upper bound, one can compute the test (-stat | -#P_ON_UPPER). 167 This will yield true even if the variable is fixed, i.e., sitting on 168 both bounds at the same time. 169 170 <b>Dual Variables</b> 171 172 In principle for implementing the Simplex algorithm it would suffice 173 to use only one dual #Status. However, for performance reasons it 174 is advisable to introduce various dual status types, reflecting 175 the structure of the bounds. Given an upper bound \f$u\f$ and a lower 176 bound \f$l\f$ of a constraint or variable, the following table 177 indicates the setting of the dual Status of this variable. 178 179 \f[ 180 \begin{array}{cl} 181 l \le ... \le u & \mbox{Status} \\ 182 \hline 183 -\infty < l \ne u < \infty & \mbox{D\_ON\_BOTH} \\ 184 -\infty < l \ne u = \infty & \mbox{D\_ON\_UPPER} \\ 185 -\infty = l \ne u < \infty & \mbox{D\_ON\_LOWER} \\ 186 -\infty < l = u < \infty & \mbox{D\_FREE} \\ 187 -\infty = l \ne u = \infty & \mbox{D\_UNDEFINED} \\ 188 \end{array} 189 \f] 190 191 Note that unbounded primal variables are reflected by an #D_UNDEFINED 192 dual variable, since no reduced costs exist for them. To facilitate 193 the assignment of dual #Status%es, class SPxBasisBase provides methods 194 #dualStatus(), #dualColStatus() and #dualRowStatus)(). 195 */ 196 enum Status 197 { 198 P_ON_LOWER = -4, ///< primal variable is set to its lower bound 199 P_ON_UPPER = -2, ///< primal variable is set to its upper bound 200 P_FREE = -1, ///< primal variable is left free, but unset 201 P_FIXED = P_ON_UPPER + P_ON_LOWER, ///< primal variable is fixed to both bounds 202 D_FREE = 1, ///< dual variable is left free, but unset 203 D_ON_UPPER = 2, ///< dual variable is set to its upper bound 204 D_ON_LOWER = 4, ///< dual variable is set to its lower bound 205 D_ON_BOTH = D_ON_LOWER + D_ON_UPPER, ///< dual variable has two bounds 206 D_UNDEFINED = 8 ///< primal or dual variable is undefined 207 }; 208 ///@} 209 210 friend SPxBasisBase<R>; 211 template <class T> friend std::ostream& operator<< (std::ostream& os, 212 const Status& stat); //@todo is the <> required here? 213 214 private: 215 216 //------------------------------------ 217 ///@name Data 218 ///@{ 219 DataArray < Status > rowstat; ///< status of rows. 220 DataArray < Status > colstat; ///< status of columns. 221 DataArray < Status >* stat; ///< basis' status. 222 DataArray < Status >* costat; ///< cobasis' status. 223 ///@} 224 225 public: 226 227 //------------------------------------ 228 ///@name Access / modification 229 ///@{ 230 /// returns number of columns. 231 int nCols() const 232 { 233 return colstat.size(); 234 } 235 /// returns number of rows. 236 int nRows() const 237 { 238 return rowstat.size(); 239 } 240 /// returns dimension. 241 int dim() const 242 { 243 return stat->size(); 244 } 245 /// returns codimension. 246 int coDim() const 247 { 248 return costat->size(); 249 } 250 /// 251 Status& rowStatus(int i) 252 { 253 return rowstat[i]; 254 } 255 /// returns status of row \p i. 256 Status rowStatus(int i) const 257 { 258 return rowstat[i]; 259 } 260 /// returns the array of row \ref soplex::SPxBasisBase<R>::Desc::Status "Status"es. 261 const Status* rowStatus(void) const 262 { 263 return rowstat.get_const_ptr(); 264 } 265 /// 266 Status& colStatus(int i) 267 { 268 return colstat[i]; 269 } 270 /// returns status of column \p i. 271 Status colStatus(int i) const 272 { 273 return colstat[i]; 274 } 275 /// returns the array of column \ref soplex::SPxBasisBase<R>::Desc::Status "Status"es. 276 const Status* colStatus(void) const 277 { 278 return colstat.get_const_ptr(); 279 } 280 /// 281 Status& status(int i) 282 { 283 return (*stat)[i]; 284 } 285 /// returns status of variable \p i. 286 Status status(int i) const 287 { 288 return (*stat)[i]; 289 } 290 /// returns the array of variable \ref soplex::SPxBasisBase<R>::Desc::Status "Status"es. 291 const Status* status(void) const 292 { 293 return stat->get_const_ptr(); 294 } 295 /// 296 Status& coStatus(int i) 297 { 298 return (*costat)[i]; 299 } 300 /// returns status of covariable \p i. 301 Status coStatus(int i) const 302 { 303 return (*costat)[i]; 304 } 305 /// returns the array of covariable \ref soplex::SPxBasisBase<R>::Desc::Status "Status"es. 306 const Status* coStatus(void) const 307 { 308 return costat->get_const_ptr(); 309 } 310 /// resets dimensions. 311 void reSize(int rowDim, int colDim); 312 ///@} 313 314 //------------------------------------ 315 ///@name Debugging 316 ///@{ 317 /// Prints out status. 318 void dump() const; 319 320 /// consistency check. 321 bool isConsistent() const; 322 ///@} 323 324 //------------------------------------ 325 ///@name Construction / destruction 326 ///@{ 327 /// default constructor 328 Desc() 329 : stat(0) 330 , costat(0) 331 {} 332 explicit Desc(const SPxSolverBase<R>& base); 333 334 /// copy constructor 335 Desc(const Desc& old); 336 /// assignment operator 337 Desc& operator=(const Desc& rhs); 338 ///@} 339 }; 340 341 protected: 342 343 //------------------------------------ 344 //**@name Protected data 345 /** 346 For storing the basis matrix we keep two arrays: Array #theBaseId 347 contains the SPxId%s of the basis vectors, and #matrix the pointers to 348 the vectors themselfes. Method #loadMatrixVecs() serves for loading 349 #matrix according to the SPxId%s stored in #theBaseId. This method must 350 be called whenever the VectorBase<R> pointers may have 351 changed due to manipulations of the LP. 352 */ 353 ///@{ 354 /// the LP 355 SPxSolverBase<R>* theLP; 356 /// SPxId%s of basic vectors. 357 DataArray < SPxId > theBaseId; 358 /// pointers to the vectors of the basis matrix. 359 DataArray < const SVectorBase<R>* > matrix; 360 /// \c true iff the pointers in \ref soplex::SPxBasisBase<R>::matrix "matrix" are set up correctly. 361 bool matrixIsSetup; 362 363 /* @brief LU factorization of basis matrix 364 The factorization of the matrix is stored in #factor if #factorized != 0. 365 Otherwise #factor is undefined. 366 */ 367 SLinSolver<R>* factor; 368 /// \c true iff \ref soplex::SPxBasisBase<R>::factor "factor" = \ref soplex::SPxBasisBase<R>::matrix "matrix" \f$^{-1}\f$. 369 bool factorized; 370 371 /// number of updates before refactorization. 372 /** When a vector of the basis matrix is exchanged by a call to method 373 #change(), the LU factorization of the matrix is updated 374 accordingly. However, after atmost #maxUpdates updates of the 375 factorization, it is recomputed in order to regain numerical 376 stability and reduce fill in. 377 */ 378 int maxUpdates; 379 380 /// allowed increase of nonzeros before refactorization. 381 /** When the number of nonzeros in LU factorization exceeds 382 #nonzeroFactor times the number of nonzeros in B, the 383 basis matrix is refactorized. 384 */ 385 R nonzeroFactor; 386 387 /// allowed increase in relative fill before refactorization 388 /** When the real relative fill is bigger than fillFactor times lastFill 389 * the Basis will be refactorized. 390 */ 391 R fillFactor; 392 393 /// allowed total increase in memory consumption before refactorization 394 R memFactor; 395 396 /* Rank-1-updates to the basis may be performed via method #change(). In 397 this case, the factorization is updated, and the following members are 398 reset. 399 */ 400 int iterCount; ///< number of calls to change() since last manipulation 401 int lastIterCount; ///< number of calls to change() before halting the simplex 402 int iterDegenCheck;///< number of calls to change() since last degeneracy check 403 int updateCount; ///< number of calls to change() since last factorize() 404 int totalUpdateCount; ///< number of updates 405 int nzCount; ///< number of nonzeros in basis matrix 406 int lastMem; ///< memory needed after last fresh factorization 407 R lastFill; ///< fill ratio that occured during last factorization 408 int lastNzCount; ///< number of nonzeros in basis matrix after last fresh factorization 409 410 Timer* theTime; ///< time spent in updates 411 Timer::TYPE timerType; ///< type of timer (user or wallclock) 412 413 SPxId lastin; ///< lastEntered(): variable entered the base last 414 SPxId lastout; ///< lastLeft(): variable left the base last 415 int lastidx; ///< lastIndex(): basis index where last update was done 416 R minStab; ///< minimum stability 417 ///@} 418 419 private: 420 421 //------------------------------------ 422 //**@name Private data */ 423 ///@{ 424 SPxStatus thestatus; ///< current status of the basis. 425 Desc thedesc; ///< the basis' Descriptor 426 bool freeSlinSolver; ///< true iff factor should be freed inside of this object 427 SPxOut* spxout; ///< message handler 428 429 ///@} 430 431 public: 432 433 //------------------------------------------------ 434 /**@name Status and Descriptor related Methods */ 435 ///@{ 436 /// returns current SPxStatus. 437 SPxStatus status() const 438 { 439 return thestatus; 440 } 441 442 /// sets basis SPxStatus to \p stat. 443 void setStatus(SPxStatus stat) 444 { 445 446 if(thestatus != stat) 447 { 448 SPxOut::debug(this, "DBSTAT01 SPxBasisBase<R>::setStatus(): status: {} ({}) -> {} ({})\n", 449 int(thestatus), thestatus, int(stat), stat); 450 451 thestatus = stat; 452 453 if(stat == NO_PROBLEM) 454 invalidate(); 455 } 456 } 457 458 // TODO control factorization frequency dynamically 459 /// change maximum number of iterations until a refactorization is performed 460 void setMaxUpdates(int maxUp) 461 { 462 assert(maxUp >= 0); 463 maxUpdates = maxUp; 464 } 465 466 /// returns maximum number of updates before a refactorization is performed 467 int getMaxUpdates() const 468 { 469 return maxUpdates; 470 } 471 472 /// 473 const Desc& desc() const 474 { 475 return thedesc; 476 } 477 /// returns current basis Descriptor. 478 Desc& desc() 479 { 480 return thedesc; 481 } 482 483 /// dual Status for the \p i'th column variable of the loaded LP. 484 typename Desc::Status dualColStatus(int i) const; 485 486 /// dual Status for the column variable with ID \p id of the loaded LP. 487 typename Desc::Status dualStatus(const SPxColId& id) const; 488 489 /// dual Status for the \p i'th row variable of the loaded LP. 490 typename Desc::Status dualRowStatus(int i) const; 491 492 /// dual Status for the row variable with ID \p id of the loaded LP. 493 typename Desc::Status dualStatus(const SPxRowId& id) const; 494 495 /// dual Status for the variable with ID \p id of the loaded LP. 496 /** It is automatically detected, whether the \p id is one of a 497 row or a column variable, and the correct row or column status 498 is returned. 499 */ 500 typename Desc::Status dualStatus(const SPxId& id) const 501 { 502 return id.isSPxRowId() 503 ? dualStatus(SPxRowId(id)) 504 : dualStatus(SPxColId(id)); 505 } 506 ///@} 507 508 509 //----------------------------------- 510 /**@name Inquiry Methods */ 511 ///@{ 512 /// 513 inline SPxId& baseId(int i) 514 { 515 return theBaseId[i]; 516 } 517 /// returns the Id of the \p i'th basis vector. 518 inline SPxId baseId(int i) const 519 { 520 return theBaseId[i]; 521 } 522 523 /// returns the \p i'th basic vector. 524 const SVectorBase<R>& baseVec(int i) const 525 { 526 assert(matrixIsSetup); 527 return *matrix[i]; 528 } 529 530 /// returns SPxId of last VectorBase<R> included to the basis. 531 inline SPxId lastEntered() const 532 { 533 return lastin; 534 } 535 536 /// returns SPxId of last vector that left the basis. 537 inline SPxId lastLeft() const 538 { 539 return lastout; 540 } 541 542 /// returns index in basis where last update was done. 543 inline int lastIndex() const 544 { 545 return lastidx; 546 } 547 548 /// returns number of basis changes since last refactorization. 549 inline int lastUpdate() const 550 { 551 return updateCount; 552 } 553 554 /// returns number of basis changes since last \ref soplex::SPxBasisBase<R>::load() "load()". 555 inline int iteration() const 556 { 557 return iterCount; 558 } 559 560 /// returns the number of iterations prior to the last break in execution 561 inline int prevIteration() const 562 { 563 return lastIterCount; 564 } 565 566 /// returns the number of iterations since the last degeneracy check 567 inline int lastDegenCheck() const 568 { 569 return iterDegenCheck; 570 } 571 572 /// returns loaded solver. 573 inline SPxSolverBase<R>* solver() const 574 { 575 return theLP; 576 } 577 ///@} 578 579 //----------------------------------- 580 /**@name Linear Algebra */ 581 ///@{ 582 /// Basis-vector product. 583 /** Depending on the representation, for an SPxBasisBase B, 584 B.multBaseWith(x) computes 585 - \f$x \leftarrow Bx\f$ in the columnwise case, and 586 - \f$x \leftarrow x^TB\f$ in the rowwise case. 587 588 Both can be seen uniformly as multiplying the basis matrix \p B with 589 a vector \p x aligned the same way as the \em vectors of \p B. 590 */ 591 VectorBase<R>& multBaseWith(VectorBase<R>& x) const; 592 593 /// Basis-vector product 594 void multBaseWith(SSVectorBase<R>& x, SSVectorBase<R>& result) const; 595 596 /// Vector-basis product. 597 /** Depending on the representation, for a #SPxBasisBase B, 598 B.multWithBase(x) computes 599 - \f$x \leftarrow x^TB\f$ in the columnwise case and 600 - \f$x \leftarrow Bx\f$ in the rowwise case. 601 602 Both can be seen uniformly as multiplying the basis matrix \p B with 603 a vector \p x aligned the same way as the \em covectors of \p B. 604 */ 605 VectorBase<R>& multWithBase(VectorBase<R>& x) const; 606 607 /// VectorBase<R>-basis product 608 void multWithBase(SSVectorBase<R>& x, SSVectorBase<R>& result) const; 609 610 /* compute an estimated condition number for the current basis matrix 611 * by computing estimates of the norms of B and B^-1 using the power method. 612 * maxiters and tolerance control the accuracy of the estimate. 613 */ 614 R condition(int maxiters = 10, R tolerance = 1e-6); 615 616 /* wrapper to compute an estimate of the condition number of the current basis matrix */ 617 R getEstimatedCondition() 618 { 619 return condition(20, 1e-6); 620 } 621 622 /* wrapper to compute the exact condition number of the current basis matrix */ 623 R getExactCondition() 624 { 625 return condition(1000, 1e-9); 626 } 627 628 /** compute one of several matrix metrics based on the diagonal of the LU factorization 629 * type = 0: max/min ratio 630 * type = 1: trace of U (sum of diagonal elements) 631 * type = 2: determinant (product of diagonal elements) 632 */ 633 R getMatrixMetric(int type = 0); 634 635 /// returns the stability of the basis matrix. 636 R stability() const 637 { 638 return factor->stability(); 639 } 640 /// 641 void solve(VectorBase<R>& x, const VectorBase<R>& rhs) 642 { 643 if(rhs.dim() == 0) 644 { 645 x.clear(); 646 return; 647 } 648 649 if(!factorized) 650 SPxBasisBase<R>::factorize(); 651 652 factor->solveRight(x, rhs); 653 } 654 /// 655 void solve(SSVectorBase<R>& x, const SVectorBase<R>& rhs) 656 { 657 if(rhs.size() == 0) 658 { 659 x.clear(); 660 return; 661 } 662 663 if(!factorized) 664 SPxBasisBase<R>::factorize(); 665 666 factor->solveRight(x, rhs); 667 } 668 /// solves linear system with basis matrix. 669 /** Depending on the representation, for a SPxBasisBase B, 670 B.solve(x) computes 671 - \f$x \leftarrow B^{-1}rhs\f$ in the columnwise case and 672 - \f$x \leftarrow rhs^TB^{-1}\f$ in the rowwise case. 673 674 Both can be seen uniformly as solving a linear system with the basis 675 matrix \p B and a right handside vector \p x aligned the same way as 676 the \em vectors of \p B. 677 */ 678 void solve4update(SSVectorBase<R>& x, const SVectorBase<R>& rhs) 679 { 680 if(rhs.size() == 0) 681 { 682 x.clear(); 683 return; 684 } 685 686 if(!factorized) 687 SPxBasisBase<R>::factorize(); 688 689 factor->solveRight4update(x, rhs); 690 } 691 /// solves two systems in one call. 692 void solve4update(SSVectorBase<R>& x, VectorBase<R>& y, const SVectorBase<R>& rhsx, 693 SSVectorBase<R>& rhsy) 694 { 695 if(!factorized) 696 SPxBasisBase<R>::factorize(); 697 698 factor->solve2right4update(x, y, rhsx, rhsy); 699 } 700 /// solves two systems in one call using only sparse data structures 701 void solve4update(SSVectorBase<R>& x, SSVectorBase<R>& y, const SVectorBase<R>& rhsx, 702 SSVectorBase<R>& rhsy) 703 { 704 if(!factorized) 705 SPxBasisBase<R>::factorize(); 706 707 factor->solve2right4update(x, y, rhsx, rhsy); 708 } 709 /// solves three systems in one call. 710 void solve4update(SSVectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& y2, 711 const SVectorBase<R>& rhsx, SSVectorBase<R>& rhsy, SSVectorBase<R>& rhsy2) 712 { 713 if(!factorized) 714 SPxBasisBase<R>::factorize(); 715 716 assert(rhsy.isSetup()); 717 assert(rhsy2.isSetup()); 718 factor->solve3right4update(x, y, y2, rhsx, rhsy, rhsy2); 719 } 720 /// solves three systems in one call using only sparse data structures 721 void solve4update(SSVectorBase<R>& x, SSVectorBase<R>& y, SSVectorBase<R>& y2, 722 const SVectorBase<R>& rhsx, SSVectorBase<R>& rhsy, SSVectorBase<R>& rhsy2) 723 { 724 if(!factorized) 725 SPxBasisBase<R>::factorize(); 726 727 assert(rhsy.isSetup()); 728 assert(rhsy2.isSetup()); 729 factor->solve3right4update(x, y, y2, rhsx, rhsy, rhsy2); 730 } 731 /// Cosolves linear system with basis matrix. 732 /** Depending on the representation, for a SPxBasisBase B, 733 B.coSolve(x) computes 734 - \f$x \leftarrow rhs^TB^{-1}\f$ in the columnwise case and 735 - \f$x \leftarrow B^{-1}rhs\f$ in the rowwise case. 736 737 Both can be seen uniformly as solving a linear system with the basis 738 matrix \p B and a right handside vector \p x aligned the same way as 739 the \em covectors of \p B. 740 */ 741 void coSolve(VectorBase<R>& x, const VectorBase<R>& rhs) 742 { 743 if(rhs.dim() == 0) 744 { 745 x.clear(); 746 return; 747 } 748 749 if(!factorized) 750 SPxBasisBase<R>::factorize(); 751 752 factor->solveLeft(x, rhs); 753 } 754 /// Sparse version of coSolve 755 void coSolve(SSVectorBase<R>& x, const SVectorBase<R>& rhs) 756 { 757 if(rhs.size() == 0) 758 { 759 x.clear(); 760 return; 761 } 762 763 if(!factorized) 764 SPxBasisBase<R>::factorize(); 765 766 factor->solveLeft(x, rhs); 767 } 768 /// solves two systems in one call. 769 void coSolve(SSVectorBase<R>& x, VectorBase<R>& y, const SVectorBase<R>& rhsx, 770 SSVectorBase<R>& rhsy) 771 { 772 if(!factorized) 773 SPxBasisBase<R>::factorize(); 774 775 factor->solveLeft(x, y, rhsx, rhsy); 776 } 777 /// Sparse version of solving two systems in one call. 778 void coSolve(SSVectorBase<R>& x, SSVectorBase<R>& y, const SVectorBase<R>& rhsx, 779 SSVectorBase<R>& rhsy) 780 { 781 if(!factorized) 782 SPxBasisBase<R>::factorize(); 783 784 factor->solveLeft(x, y, rhsx, rhsy); 785 } 786 /// solves three systems in one call. May be improved by using just one pass through the basis. 787 void coSolve(SSVectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& z, const SVectorBase<R>& rhsx, 788 SSVectorBase<R>& rhsy, SSVectorBase<R>& rhsz) 789 { 790 if(!factorized) 791 SPxBasisBase<R>::factorize(); 792 793 factor->solveLeft(x, y, z, rhsx, rhsy, rhsz); 794 } 795 /// Sparse version of solving three systems in one call. 796 void coSolve(SSVectorBase<R>& x, SSVectorBase<R>& y, SSVectorBase<R>& z, const SVectorBase<R>& rhsx, 797 SSVectorBase<R>& rhsy, SSVectorBase<R>& rhsz) 798 { 799 if(!factorized) 800 SPxBasisBase<R>::factorize(); 801 802 factor->solveLeft(x, y, z, rhsx, rhsy, rhsz); 803 } 804 ///@} 805 806 807 //------------------------------------ 808 /**@name Modification notification. 809 These methods must be called after the loaded LP has been modified. 810 */ 811 ///@{ 812 /// inform SPxBasisBase, that \p n new rows had been added. 813 void addedRows(int n); 814 /// inform SPxBasisBase that row \p i had been removed. 815 void removedRow(int i); 816 /// inform SPxBasisBase that rows in \p perm with negative entry were removed. 817 void removedRows(const int perm[]); 818 /// inform SPxBasisBase that \p n new columns had been added. 819 void addedCols(int n); 820 /// inform SPxBasisBase that column \p i had been removed. 821 void removedCol(int i); 822 /// inform SPxBasisBase that columns in \p perm with negative entry were removed. 823 void removedCols(const int perm[]); 824 /// inform SPxBasisBase that a row had been changed. 825 void changedRow(int); 826 /// inform SPxBasisBase that a column had been changed. 827 void changedCol(int); 828 /// inform SPxBasisBase that a matrix entry had been changed. 829 void changedElement(int, int); 830 ///@} 831 832 833 //-------------------------------- 834 /**@name Miscellaneous */ 835 ///@{ 836 /// performs basis update. 837 /** Changes the \p i 'th vector of the basis with the vector associated to 838 \p id. This includes: 839 - updating the factorization, or recomputing it from scratch by 840 calling \ref soplex::SPxSolverBase<R>::factorize() "factorize()", 841 - resetting \ref soplex::SPxSolverBase<R>::lastEntered() "lastEntered()", 842 - resetting \ref soplex::SPxSolverBase<R>::lastIndex() "lastIndex()", 843 - resetting \ref soplex::SPxSolverBase<R>::lastLeft() "lastLeft()", 844 - resetting \ref soplex::SPxSolverBase<R>::lastUpdate() "lastUpdate()", 845 - resetting \ref soplex::SPxSolverBase<R>::iterations() "iterations()". 846 847 The basis descriptor is \em not \em modified, since #factor() 848 cannot know about how to set up the status of the involved variables 849 correctly. 850 851 A vector \p enterVec may be passed for a fast ETA update of the LU 852 factorization associated to the basis. It must be initialized with 853 the solution vector \f$x\f$ of the right linear system \f$Bx = b\f$ 854 with the entering vector as right-hand side vector \f$b\f$, where \f$B\f$ 855 denotes the basis matrix. This can be computed using method #solve(). 856 When using FAST updates, a vector \p eta may be passed for 857 improved performance. It must be initialized by a call to 858 factor->solveRightUpdate() as described in SLinSolver. The 859 implementation hidden behind FAST updates depends on the 860 SLinSolver implementation class. 861 */ 862 virtual void change(int i, SPxId& id, 863 const SVectorBase<R>* enterVec, const SSVectorBase<R>* eta = 0); 864 865 /** Load basis from \p in in MPS format. If \p rowNames and \p colNames 866 * are \c NULL, default names are used for the constraints and variables. 867 */ 868 virtual bool readBasis(std::istream& in, 869 const NameSet* rowNames, const NameSet* colNames); 870 871 /** Write basis to \p os in MPS format. If \p rowNames and \p colNames are 872 * \c NULL, default names are used for the constraints and variables. 873 */ 874 virtual void writeBasis(std::ostream& os, 875 const NameSet* rownames, const NameSet* colnames, const bool cpxFormat = false) const; 876 877 virtual void printMatrix() const; 878 879 /** Prints current basis matrix to a file using the MatrixMarket format: 880 * row col value 881 * The filename is basis/basis[number].mtx where number is a parameter. 882 */ 883 void printMatrixMTX(int number); 884 885 /// checks if a Descriptor is valid for the current LP w.r.t. its bounds 886 virtual bool isDescValid(const Desc& ds); 887 888 /// sets up basis. 889 /** Loads a Descriptor to the basis and sets up the basis matrix and 890 all vectors accordingly. The Descriptor must have the same number of 891 rows and columns as the currently loaded LP. 892 */ 893 virtual void loadDesc(const Desc&); 894 895 /// sets up linear solver to use. 896 /** If destroy is true, solver will be freed inside this object, e.g. in the destructor. 897 */ 898 virtual void loadBasisSolver(SLinSolver<R>* solver, const bool destroy = false); 899 900 /// loads the LP \p lp to the basis. 901 /** This involves resetting all counters to 0 and setting up a regular 902 default basis consisting of slacks, artificial variables or bounds. 903 */ 904 virtual void load(SPxSolverBase<R>* lp, bool initSlackBasis = true); 905 906 /// unloads the LP from the basis. 907 virtual void unLoad() 908 { 909 theLP = 0; 910 setStatus(NO_PROBLEM); 911 } 912 913 /// invalidates actual basis. 914 /** This method makes the basis matrix and vectors invalid. The basis will 915 be reinitialized if needed. 916 */ 917 void invalidate(); 918 919 /// Restores initial basis. 920 /** This method changes the basis to that present just after loading the LP 921 (see addedRows() and addedCols()). This may be necessary if a row or a 922 column is changed, since then the current basis may become singular. 923 */ 924 void restoreInitialBasis(); 925 926 /// output basis entries. 927 void dump(); 928 929 /// consistency check. 930 bool isConsistent() const; 931 932 /// time spent in updates 933 Real getTotalUpdateTime() const 934 { 935 return theTime->time(); 936 } 937 /// number of updates performed 938 int getTotalUpdateCount() const 939 { 940 return totalUpdateCount; 941 } 942 943 /// returns statistical information in form of a string. 944 std::string statistics() const 945 { 946 std::stringstream s; 947 s << factor->statistics() 948 #ifdef SOPLEX_MEASUREUPDATETIME 949 << "Updates : " << std::setw(10) << getTotalUpdateCount() << std::endl 950 << " Time spent : " << std::setw(10) << getTotalUpdateTime() << std::endl 951 #endif 952 ; 953 954 return s.str(); 955 } 956 957 void setOutstream(SPxOut& newOutstream) 958 { 959 spxout = &newOutstream; 960 } 961 ///@} 962 963 //-------------------------------------- 964 /**@name Constructors / Destructors */ 965 ///@{ 966 /// default constructor. 967 SPxBasisBase(Timer::TYPE ttype = Timer::USER_TIME); 968 /// copy constructor 969 SPxBasisBase(const SPxBasisBase<R>& old); 970 /// assignment operator 971 SPxBasisBase<R>& operator=(const SPxBasisBase<R>& rhs); 972 /// destructor. 973 virtual ~SPxBasisBase(); 974 ///@} 975 976 977 protected: 978 979 //-------------------------------------- 980 /**@name Protected helpers */ 981 ///@{ 982 /// loads \ref soplex::SPxBasisBase<R>::matrix "matrix" according to the SPxId%s stored in \ref soplex::SPxBasisBase<R>::theBaseId "theBaseId". 983 /** This method must be called whenever there is a chance, that the vector 984 pointers may have changed due to manipulations of the LP. 985 */ 986 void loadMatrixVecs(); 987 988 /// resizes internal arrays. 989 /** When a new LP is loaded, the basis matrix and vectors become invalid 990 and possibly also of the wrong dimension. Hence, after loading an 991 LP, #reDim() is called to reset all arrays etc. accoriding to the 992 dimensions of the loaded LP. 993 */ 994 void reDim(); 995 996 /// factorizes the basis matrix. 997 virtual void factorize(); 998 999 /// sets descriptor representation according to loaded LP. 1000 void setRep(); 1001 ///@} 1002 1003 }; 1004 1005 1006 // 1007 // Auxiliary functions. 1008 // 1009 1010 /// Pretty-printing of basis status. 1011 template <class R> 1012 std::ostream& operator<<(std::ostream& os, 1013 const typename SPxBasisBase<R>::SPxStatus& status); 1014 1015 1016 /* For backwards compatibility */ 1017 typedef SPxBasisBase<Real> SPxBasis; 1018 1019 } // namespace soplex 1020 1021 // General templated definitions 1022 #include "spxbasis.hpp" 1023 #include "spxdesc.hpp" 1024 1025 /* reset the SOPLEX_DEBUG flag to its original value */ 1026 #undef SOPLEX_DEBUG 1027 #ifdef SOPLEX_DEBUG_SPXBASIS 1028 #define SOPLEX_DEBUG 1029 #undef SOPLEX_DEBUG_SPXBASIS 1030 #endif 1031 1032 #endif // _SPXBASIS_H_ 1033