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