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