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.h 26 * @brief Implementation of Sparse Linear Solver. 27 */ 28 #ifndef _SLUFACTOR_H_ 29 #define _SLUFACTOR_H_ 30 31 #include <assert.h> 32 33 #include "soplex/spxdefines.h" 34 #include "soplex/timerfactory.h" 35 #include "soplex/slinsolver.h" 36 #include "soplex/clufactor.h" 37 38 namespace soplex 39 { 40 /// maximum nr. of factorization updates allowed before refactorization. 41 #define SOPLEX_MAXUPDATES 1000 42 43 /**@brief Implementation of Sparse Linear Solver. 44 * @ingroup Algo 45 * 46 * This class implements a SLinSolver interface by using the sparse LU 47 * factorization implemented in CLUFactor. 48 */ 49 template <class R> 50 class SLUFactor : public SLinSolver<R>, protected CLUFactor<R> 51 { 52 public: 53 54 //-------------------------------- 55 /**@name Types */ 56 ///@{ 57 /// Specifies how to perform \ref soplex::SLUFactor<R>::change "change" method. 58 enum UpdateType 59 { 60 ETA = 0, ///< 61 FOREST_TOMLIN ///< 62 }; 63 /// for convenience 64 using Status = typename SLinSolver<R>::Status; 65 ///@} 66 67 private: 68 69 //-------------------------------- 70 /**@name Private data */ 71 ///@{ 72 VectorBase<R> vec; ///< Temporary VectorBase<R> 73 SSVectorBase<R> ssvec; ///< Temporary semi-sparse VectorBase<R> 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 SSVectorBase<R> eta; ///< 84 SSVectorBase<R> 85 forest; ///< ? Update VectorBase<R> set up by solveRight4update() and solve2right4update() 86 R lastThreshold; ///< pivoting threshold of last factorization 87 ///@} 88 89 //-------------------------------- 90 /**@name Control Parameters */ 91 ///@{ 92 /// minimum threshold to use. 93 R minThreshold; 94 /// minimum stability to achieve by setting threshold. 95 R 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, SSVectorBase<R>& eta); 112 ///@} 113 114 115 public: 116 117 //-------------------------------- 118 /**@name Update type */ 119 ///@{ 120 /// returns the current update type uptype. 121 UpdateType utype() const 122 { 123 return uptype; 124 } 125 126 /// sets update type. 127 /** The new UpdateType becomes valid only after the next call to 128 method load(). 129 */ 130 void setUtype(UpdateType tp) 131 { 132 uptype = tp; 133 } 134 135 /// sets minimum Markowitz threshold. 136 void setMarkowitz(R m) 137 { 138 if(m < 0.0001) 139 m = 0.0001; 140 141 if(m > 0.9999) 142 m = 0.9999; 143 144 minThreshold = m; 145 lastThreshold = m; 146 } 147 148 /// returns Markowitz threshold. 149 R markowitz() 150 { 151 return lastThreshold; 152 } 153 ///@} 154 155 //-------------------------------- 156 /**@name Derived from SLinSolver 157 See documentation of \ref soplex::SLinSolver "SLinSolver" for a 158 documentation of these methods. 159 */ 160 ///@{ 161 /// 162 void clear(); 163 /// 164 int dim() const 165 { 166 return this->thedim; 167 } 168 /// 169 int memory() const 170 { 171 return this->nzCnt + this->l.start[this->l.firstUnused]; 172 } 173 /// 174 const char* getName() const 175 { 176 return (uptype == SLUFactor<R>::ETA) ? "SLU-Eta" : "SLU-Forest-Tomlin"; 177 } 178 /// 179 Status status() const 180 { 181 return Status(this->stat); 182 } 183 /// 184 R stability() const; 185 /** return one of several matrix metrics based on the diagonal of U 186 * 0: condition number estimate by ratio of min/max 187 * 1: trace (sum of diagonal elements) 188 * 2: determinant (product of diagonal elements) 189 */ 190 R matrixMetric(int type = 0) const; 191 /// 192 std::string statistics() const; 193 /// 194 Status load(const SVectorBase<R>* vec[], int dim); 195 ///@} 196 197 public: 198 199 //-------------------------------- 200 /**@name Solve */ 201 ///@{ 202 /// Solves \f$Ax=b\f$. 203 void solveRight(VectorBase<R>& x, const VectorBase<R>& b); 204 void solveRight(SSVectorBase<R>& x, const SSVectorBase<R>& b) 205 { 206 x.unSetup(); 207 solveRight((VectorBase<R>&) x, (const VectorBase<R>&) b); 208 } 209 /// Solves \f$Ax=b\f$. 210 void solveRight(SSVectorBase<R>& x, const SVectorBase<R>& b); 211 /// Solves \f$Ax=b\f$. 212 void solveRight4update(SSVectorBase<R>& x, const SVectorBase<R>& b); 213 /// Solves \f$Ax=b\f$ and \f$Ay=d\f$. 214 void solve2right4update(SSVectorBase<R>& x, VectorBase<R>& y, const SVectorBase<R>& b, 215 SSVectorBase<R>& d); 216 /// Sparse version of solving two systems of equations 217 void solve2right4update(SSVectorBase<R>& x, SSVectorBase<R>& y, const SVectorBase<R>& b, 218 SSVectorBase<R>& d); 219 /// Solves \f$Ax=b\f$, \f$Ay=d\f$ and \f$Az=e\f$. 220 void solve3right4update(SSVectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& z, 221 const SVectorBase<R>& b, SSVectorBase<R>& d, SSVectorBase<R>& e); 222 /// sparse version of solving three systems of equations 223 void solve3right4update(SSVectorBase<R>& x, SSVectorBase<R>& y, SSVectorBase<R>& z, 224 const SVectorBase<R>& b, SSVectorBase<R>& d, SSVectorBase<R>& e); 225 /// sparse version of solving one system of equations with transposed basis matrix 226 void solveLeft(VectorBase<R>& x, const VectorBase<R>& b); 227 void solveLeft(SSVectorBase<R>& x, const SSVectorBase<R>& b) 228 { 229 x.unSetup(); 230 solveLeft((VectorBase<R>&) x, (const VectorBase<R>&) b); 231 } 232 /// Solves \f$Ax=b\f$. 233 void solveLeft(SSVectorBase<R>& x, const SVectorBase<R>& b); 234 /// Solves \f$Ax=b\f$ and \f$Ay=d\f$. 235 void solveLeft(SSVectorBase<R>& x, VectorBase<R>& y, const SVectorBase<R>& b, SSVectorBase<R>& d); 236 /// sparse version of solving two systems of equations with transposed basis matrix 237 void solveLeft(SSVectorBase<R>& x, SSVectorBase<R>& two, const SVectorBase<R>& b, 238 SSVectorBase<R>& rhs2); 239 /// Solves \f$Ax=b\f$, \f$Ay=d\f$ and \f$Az=e\f$. 240 void solveLeft(SSVectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& z, 241 const SVectorBase<R>& b, SSVectorBase<R>& d, SSVectorBase<R>& e); 242 /// sparse version of solving three systems of equations with transposed basis matrix 243 void solveLeft(SSVectorBase<R>& x, SSVectorBase<R>& y, SSVectorBase<R>& z, 244 const SVectorBase<R>& b, SSVectorBase<R>& d, SSVectorBase<R>& e); 245 /// 246 Status change(int idx, const SVectorBase<R>& subst, const SSVectorBase<R>* eta = 0); 247 ///@} 248 249 //-------------------------------- 250 /**@name Miscellaneous */ 251 ///@{ 252 /// time spent in factorizations 253 // @todo fix the return type from of the type form Real to a cpp time (Refactoring) TODO 254 Real getFactorTime() const 255 { 256 return this->factorTime->time(); 257 } 258 /// reset FactorTime 259 void resetFactorTime() 260 { 261 this->factorTime->reset(); 262 } 263 /// number of factorizations performed 264 int getFactorCount() const 265 { 266 return this->factorCount; 267 } 268 /// time spent in solves 269 // @todo fix the return type of time to a cpp time type TODO 270 Real getSolveTime() const 271 { 272 return solveTime->time(); 273 } 274 /// reset SolveTime 275 void resetSolveTime() 276 { 277 solveTime->reset(); 278 } 279 /// number of solves performed 280 int getSolveCount() const 281 { 282 return solveCount; 283 } 284 /// reset timers and counters 285 void resetCounters() 286 { 287 this->factorTime->reset(); 288 solveTime->reset(); 289 this->factorCount = 0; 290 this->hugeValues = 0; 291 solveCount = 0; 292 } 293 void changeTimer(const Timer::TYPE ttype) 294 { 295 solveTime = TimerFactory::switchTimer(solveTime, ttype); 296 this->factorTime = TimerFactory::switchTimer(this->factorTime, ttype); 297 timerType = ttype; 298 } 299 /// prints the LU factorization to stdout. 300 void dump() const; 301 302 /// consistency check. 303 bool isConsistent() const; 304 305 /// set tolerances 306 virtual void setTolerances(std::shared_ptr<Tolerances> tolerances) 307 { 308 this->_tolerances = tolerances; 309 this->eta.setTolerances(tolerances); 310 this->forest.setTolerances(tolerances); 311 this->ssvec.setTolerances(tolerances); 312 } 313 314 ///@} 315 316 //------------------------------------ 317 /**@name Constructors / Destructors */ 318 ///@{ 319 /// default constructor. 320 SLUFactor(); 321 /// assignment operator. 322 SLUFactor<R>& operator=(const SLUFactor<R>& old); 323 /// copy constructor. 324 SLUFactor(const SLUFactor<R>& old); 325 /// destructor. 326 virtual ~SLUFactor(); 327 /// clone function for polymorphism 328 inline virtual SLinSolver<R>* clone() const 329 { 330 return new SLUFactor<R>(*this); 331 } 332 ///@} 333 334 private: 335 336 //------------------------------------ 337 /**@name Private helpers */ 338 ///@{ 339 /// used to implement the assignment operator 340 void assign(const SLUFactor<R>& old); 341 ///@} 342 }; 343 344 } // namespace soplex 345 346 #include "slufactor.hpp" 347 #endif // _SLUFACTOR_H_ 348