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 slufactor_rational.h 26 * @brief Implementation of Sparse Linear Solver with Rational precision. 27 */ 28 #ifndef _SLUFACTOR_RATIONAL_H_ 29 #define _SLUFACTOR_RATIONAL_H_ 30 31 #include <assert.h> 32 33 #include "soplex/spxdefines.h" 34 #include "soplex/timerfactory.h" 35 #include "soplex/slinsolver_rational.h" 36 #include "soplex/clufactor_rational.h" 37 #include "soplex/rational.h" 38 39 namespace soplex 40 { 41 /// maximum nr. of factorization updates allowed before refactorization. 42 #define SOPLEX_MAXUPDATES 1000 43 44 /**@brief Implementation of Sparse Linear Solver with Rational precision. 45 * @ingroup Algo 46 * 47 * This class implements a SLinSolverRational interface by using the sparse LU 48 * factorization implemented in CLUFactorRational. 49 */ 50 class SLUFactorRational : public SLinSolverRational, protected CLUFactorRational 51 { 52 public: 53 54 //-------------------------------- 55 /**@name Types */ 56 ///@{ 57 /// Specifies how to perform \ref soplex::SLUFactorRational::change "change" method. 58 enum UpdateType 59 { 60 ETA = 0, ///< 61 FOREST_TOMLIN ///< 62 }; 63 /// for convenience 64 typedef SLinSolverRational::Status Status; 65 ///@} 66 67 private: 68 69 //-------------------------------- 70 /**@name Private data */ 71 ///@{ 72 VectorRational vec; ///< Temporary vector 73 SSVectorRational ssvec; ///< Temporary semi-sparse vector 74 ///@} 75 76 protected: 77 78 //-------------------------------- 79 /**@name Protected data */ 80 ///@{ 81 bool usetup; ///< TRUE iff update vector has been setup 82 UpdateType uptype; ///< the current \ref soplex::SLUFactor<R>::UpdateType "UpdateType". 83 SSVectorRational eta; ///< 84 SSVectorRational 85 forest; ///< ? Update vector set up by solveRight4update() and solve2right4update() 86 Rational lastThreshold; ///< pivoting threshold of last factorization 87 ///@} 88 89 //-------------------------------- 90 /**@name Control Parameters */ 91 ///@{ 92 /// minimum threshold to use. 93 Rational minThreshold; 94 /// minimum stability to achieve by setting threshold. 95 Rational minStability; 96 /// Time spent in solves 97 Timer* solveTime; 98 Timer::TYPE timerType; 99 /// Number of solves 100 int solveCount; 101 ///@} 102 103 protected: 104 105 //-------------------------------- 106 /**@name Protected helpers */ 107 ///@{ 108 /// 109 void freeAll(); 110 /// 111 void changeEta(int idx, SSVectorRational& eta); 112 /// 113 void init(); 114 ///@} 115 116 117 public: 118 119 //-------------------------------- 120 /**@name Update type */ 121 ///@{ 122 /// returns the current update type uptype. 123 UpdateType utype() const 124 { 125 return uptype; 126 } 127 128 /// sets update type. 129 /** The new UpdateType becomes valid only after the next call to 130 method load(). 131 */ 132 void setUtype(UpdateType tp) 133 { 134 uptype = tp; 135 } 136 137 /// sets minimum Markowitz threshold. 138 void setMarkowitz(const Rational& m) 139 { 140 if(m < 0.01) 141 { 142 minThreshold = 0.01; 143 lastThreshold = 0.01; 144 } 145 else if(m > 0.99) 146 { 147 minThreshold = 0.99; 148 lastThreshold = 0.99; 149 } 150 else 151 { 152 minThreshold = m; 153 lastThreshold = m; 154 } 155 } 156 157 /// returns Markowitz threshold. 158 Rational markowitz() 159 { 160 return lastThreshold; 161 } 162 ///@} 163 164 //-------------------------------- 165 /**@name Derived from SLinSolverRational 166 See documentation of \ref soplex::SLinSolverRational "SLinSolverRational" for a 167 documentation of these methods. 168 */ 169 ///@{ 170 /// 171 void clear(); 172 /// 173 int dim() const 174 { 175 return thedim; 176 } 177 /// 178 int memory() const 179 { 180 return nzCnt + l.start[l.firstUnused]; 181 } 182 /// 183 const char* getName() const 184 { 185 return (uptype == SLUFactorRational::ETA) ? "SLU-Eta" : "SLU-Forest-Tomlin"; 186 } 187 /// 188 Status status() const 189 { 190 return Status(stat); 191 } 192 /// 193 Rational stability() const; 194 /// 195 std::string statistics() const; 196 /// 197 Status load(const SVectorRational* vec[], int dim); 198 ///@} 199 200 public: 201 202 //-------------------------------- 203 /**@name Solve */ 204 ///@{ 205 /// Solves \f$Ax=b\f$. 206 void solveRight(VectorRational& x, const VectorRational& b); 207 /// Solves \f$Ax=b\f$. 208 void solveRight(SSVectorRational& x, const SVectorRational& b); 209 /// Solves \f$Ax=b\f$. 210 void solveRight4update(SSVectorRational& x, const SVectorRational& b); 211 /// Solves \f$Ax=b\f$ and \f$Ay=d\f$. 212 void solve2right4update(SSVectorRational& x, VectorRational& y, const SVectorRational& b, 213 SSVectorRational& d); 214 /// Solves \f$Ax=b\f$, \f$Ay=d\f$ and \f$Az=e\f$. 215 void solve3right4update(SSVectorRational& x, VectorRational& y, VectorRational& z, 216 const SVectorRational& b, SSVectorRational& d, SSVectorRational& e); 217 /// Solves \f$Ax=b\f$. 218 void solveLeft(VectorRational& x, const VectorRational& b); 219 /// Solves \f$Ax=b\f$. 220 void solveLeft(SSVectorRational& x, const SVectorRational& b); 221 /// Solves \f$Ax=b\f$ and \f$Ay=d\f$. 222 void solveLeft(SSVectorRational& x, VectorRational& y, const SVectorRational& b, 223 SSVectorRational& d); 224 /// Solves \f$Ax=b\f$, \f$Ay=d\f$ and \f$Az=e\f$. 225 void solveLeft(SSVectorRational& x, VectorRational& y, VectorRational& z, 226 const SVectorRational& b, SSVectorRational& d, SSVectorRational& e); 227 /// 228 Status change(int idx, const SVectorRational& subst, const SSVectorRational* eta = 0); 229 ///@} 230 231 //-------------------------------- 232 /**@name Miscellaneous */ 233 ///@{ 234 /// time spent in factorizations 235 Real getFactorTime() const 236 { 237 return factorTime->time(); 238 } 239 /// set time limit on factorization 240 void setTimeLimit(const Real limit) 241 { 242 timeLimit = limit; 243 } 244 /// reset FactorTime 245 void resetFactorTime() 246 { 247 factorTime->reset(); 248 } 249 /// number of factorizations performed 250 int getFactorCount() const 251 { 252 return factorCount; 253 } 254 /// time spent in solves 255 Real getSolveTime() const 256 { 257 return solveTime->time(); 258 } 259 /// reset SolveTime 260 void resetSolveTime() 261 { 262 solveTime->reset(); 263 } 264 /// number of solves performed 265 int getSolveCount() const 266 { 267 return solveCount; 268 } 269 /// reset timers and counters 270 void resetCounters() 271 { 272 factorTime->reset(); 273 solveTime->reset(); 274 factorCount = 0; 275 solveCount = 0; 276 } 277 /// prints the LU factorization to stdout. 278 void dump() const; 279 280 /// consistency check. 281 bool isConsistent() const; 282 ///@} 283 284 //------------------------------------ 285 /**@name Constructors / Destructors */ 286 ///@{ 287 /// default constructor. 288 SLUFactorRational() 289 : CLUFactorRational() 290 , vec(1) 291 , ssvec(1) 292 , usetup(false) 293 , uptype(FOREST_TOMLIN) 294 , eta(1) 295 , forest(1) 296 , minThreshold(0.01) 297 , timerType(Timer::USER_TIME) 298 { 299 row.perm = 0; 300 row.orig = 0; 301 col.perm = 0; 302 col.orig = 0; 303 u.row.elem = 0; 304 u.row.idx = 0; 305 u.row.start = 0; 306 u.row.len = 0; 307 u.row.max = 0; 308 u.col.elem = 0; 309 u.col.idx = 0; 310 u.col.start = 0; 311 u.col.len = 0; 312 u.col.max = 0; 313 l.idx = 0; 314 l.start = 0; 315 l.row = 0; 316 l.ridx = 0; 317 l.rbeg = 0; 318 l.rorig = 0; 319 l.rperm = 0; 320 321 nzCnt = 0; 322 thedim = 0; 323 324 try 325 { 326 solveTime = TimerFactory::createTimer(timerType); 327 factorTime = TimerFactory::createTimer(timerType); 328 spx_alloc(row.perm, thedim); 329 spx_alloc(row.orig, thedim); 330 spx_alloc(col.perm, thedim); 331 spx_alloc(col.orig, thedim); 332 diag.reDim(thedim); 333 334 work = vec.get_ptr(); 335 336 u.row.used = 0; 337 spx_alloc(u.row.elem, thedim); 338 u.row.val.reDim(1); 339 spx_alloc(u.row.idx, u.row.val.dim()); 340 spx_alloc(u.row.start, thedim + 1); 341 spx_alloc(u.row.len, thedim + 1); 342 spx_alloc(u.row.max, thedim + 1); 343 344 u.row.list.idx = thedim; 345 u.row.start[thedim] = 0; 346 u.row.max [thedim] = 0; 347 u.row.len [thedim] = 0; 348 349 u.col.size = 1; 350 u.col.used = 0; 351 spx_alloc(u.col.elem, thedim); 352 spx_alloc(u.col.idx, u.col.size); 353 spx_alloc(u.col.start, thedim + 1); 354 spx_alloc(u.col.len, thedim + 1); 355 spx_alloc(u.col.max, thedim + 1); 356 u.col.val.reDim(0); 357 358 u.col.list.idx = thedim; 359 u.col.start[thedim] = 0; 360 u.col.max[thedim] = 0; 361 u.col.len[thedim] = 0; 362 363 l.val.reDim(1); 364 spx_alloc(l.idx, l.val.dim()); 365 366 l.startSize = 1; 367 l.firstUpdate = 0; 368 l.firstUnused = 0; 369 370 spx_alloc(l.start, l.startSize); 371 spx_alloc(l.row, l.startSize); 372 } 373 catch(const SPxMemoryException& x) 374 { 375 freeAll(); 376 throw x; 377 } 378 379 l.rval.reDim(0); 380 l.ridx = 0; 381 l.rbeg = 0; 382 l.rorig = 0; 383 l.rperm = 0; 384 385 SLUFactorRational::init(); 386 387 factorCount = 0; 388 timeLimit = -1.0; 389 solveCount = 0; 390 391 assert(row.perm != 0); 392 assert(row.orig != 0); 393 assert(col.perm != 0); 394 assert(col.orig != 0); 395 396 assert(u.row.elem != 0); 397 assert(u.row.idx != 0); 398 assert(u.row.start != 0); 399 assert(u.row.len != 0); 400 assert(u.row.max != 0); 401 402 assert(u.col.elem != 0); 403 assert(u.col.idx != 0); 404 assert(u.col.start != 0); 405 assert(u.col.len != 0); 406 assert(u.col.max != 0); 407 408 assert(l.idx != 0); 409 assert(l.start != 0); 410 assert(l.row != 0); 411 412 assert(SLUFactorRational::isConsistent()); 413 } 414 /// assignment operator. 415 SLUFactorRational& operator=(const SLUFactorRational& old) 416 { 417 418 if(this != &old) 419 { 420 // we don't need to copy them, because they are temporary vectors 421 vec.clear(); 422 ssvec.clear(); 423 424 eta = old.eta; 425 forest = old.forest; 426 427 freeAll(); 428 429 try 430 { 431 assign(old); 432 } 433 catch(const SPxMemoryException& x) 434 { 435 freeAll(); 436 throw x; 437 } 438 439 assert(isConsistent()); 440 } 441 442 return *this; 443 } 444 445 /// copy constructor. 446 SLUFactorRational(const SLUFactorRational& old) 447 : SLinSolverRational(old) 448 , CLUFactorRational() 449 , vec(1) // we don't need to copy it, because they are temporary vectors 450 , ssvec(1) // we don't need to copy it, because they are temporary vectors 451 , usetup(old.usetup) 452 , eta(old.eta) 453 , forest(old.forest) 454 , timerType(old.timerType) 455 { 456 row.perm = 0; 457 row.orig = 0; 458 col.perm = 0; 459 col.orig = 0; 460 u.row.elem = 0; 461 u.row.idx = 0; 462 u.row.start = 0; 463 u.row.len = 0; 464 u.row.max = 0; 465 u.col.elem = 0; 466 u.col.idx = 0; 467 u.col.start = 0; 468 u.col.len = 0; 469 u.col.max = 0; 470 l.idx = 0; 471 l.start = 0; 472 l.row = 0; 473 l.ridx = 0; 474 l.rbeg = 0; 475 l.rorig = 0; 476 l.rperm = 0; 477 478 solveCount = 0; 479 solveTime = TimerFactory::createTimer(timerType); 480 factorTime = TimerFactory::createTimer(timerType); 481 482 try 483 { 484 assign(old); 485 } 486 catch(const SPxMemoryException& x) 487 { 488 freeAll(); 489 throw x; 490 } 491 492 assert(SLUFactorRational::isConsistent()); 493 } 494 495 /// destructor. 496 virtual ~SLUFactorRational(); 497 /// clone function for polymorphism 498 inline virtual SLinSolverRational* clone() const 499 { 500 return new SLUFactorRational(*this); 501 } 502 ///@} 503 504 private: 505 506 //------------------------------------ 507 /**@name Private helpers */ 508 ///@{ 509 /// used to implement the assignment operator 510 void assign(const SLUFactorRational& old); 511 ///@} 512 }; 513 514 } // namespace soplex 515 #include "slufactor_rational.hpp" 516 #endif // _SLUFACTOR_RATIONAL_H_ 517