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