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   	#include <assert.h>
26   	#include <cstdio>
27   	#include <iostream>
28   	#include <iomanip>
29   	#include <sstream>
30   	
31   	#include "soplex/spxdefines.h"
32   	#include "soplex/didxset.h"
33   	#include "soplex/mpsinput.h"
34   	#include "soplex/spxout.h"
35   	#include "soplex/exceptions.h"
36   	
37   	namespace soplex
38   	{
39   	
40   	template <class R> typename SPxBasisBase<R>::Desc::Status
41   	SPxBasisBase<R>::dualStatus(const SPxColId& id) const
42   	{
43   	   return dualColStatus(static_cast<SPxLPBase<R>*>(theLP)->number(id));
44   	}
45   	
46   	template <class R>
47   	typename SPxBasisBase<R>::Desc::Status
48   	SPxBasisBase<R>::dualStatus(const SPxRowId& id) const
49   	{
50   	   return dualRowStatus((static_cast<SPxLPBase<R>*>(theLP))->number(id));
51   	}
52   	
53   	template <class R>
54   	typename SPxBasisBase<R>::Desc::Status
55   	SPxBasisBase<R>::dualRowStatus(int i) const
56   	{
57   	   assert(theLP != 0);
58   	
59   	   if(theLP->rhs(i) < R(infinity))
60   	   {
61   	      if(theLP->lhs(i) > R(-infinity))
62   	      {
63   	         if(theLP->lhs(i) == theLP->rhs(i))
64   	            return Desc::D_FREE;
65   	         else
66   	            return Desc::D_ON_BOTH;
67   	      }
68   	      else
69   	         return Desc::D_ON_LOWER;
70   	   }
71   	   else if(theLP->lhs(i) > R(-infinity))
72   	      return Desc::D_ON_UPPER;
73   	   else
74   	      return Desc::D_UNDEFINED;
75   	}
76   	
77   	template <class R>
78   	typename SPxBasisBase<R>::Desc::Status
79   	SPxBasisBase<R>::dualColStatus(int i) const
80   	{
81   	   assert(theLP != 0);
82   	
83   	   if(theLP->SPxLPBase<R>::upper(i) < R(infinity))
84   	   {
85   	      if(theLP->SPxLPBase<R>::lower(i) > R(-infinity))
86   	      {
87   	         if(theLP->SPxLPBase<R>::lower(i) == theLP->SPxLPBase<R>::upper(i))
88   	            return Desc::D_FREE;
89   	         else
90   	            return Desc::D_ON_BOTH;
91   	      }
92   	      else
93   	         return Desc::D_ON_LOWER;
94   	   }
95   	   else if(theLP->SPxLPBase<R>::lower(i) > R(-infinity))
96   	      return Desc::D_ON_UPPER;
97   	   else
98   	      return Desc::D_UNDEFINED;
99   	}
100  	
101  	template <class R>
102  	void SPxBasisBase<R>::loadMatrixVecs()
103  	{
104  	   assert(theLP != 0);
105  	   assert(theLP->dim() == matrix.size());
106  	
107  	   SPX_MSG_INFO3((*this->spxout),
108  	                 (*this->spxout) << "IBASIS01 loadMatrixVecs() invalidates factorization"
109  	                 << std::endl;)
110  	
111  	   int i;
112  	   nzCount = 0;
113  	
114  	   for(i = theLP->dim() - 1; i >= 0; --i)
115  	   {
116  	      matrix[i] = &theLP->vector(baseId(i));
117  	      nzCount += matrix[i]->size();
118  	   }
119  	
120  	   matrixIsSetup = true;
121  	   factorized = false;
122  	
123  	   if(factor != 0)
124  	      factor->clear();
125  	}
126  	
127  	template <class R>
128  	bool SPxBasisBase<R>::isDescValid(const Desc& ds)
129  	{
130  	
131  	   assert(status() > NO_PROBLEM);
132  	   assert(theLP != 0);
133  	
134  	   int basisdim;
135  	
136  	   if(ds.nRows() != theLP->nRows() || ds.nCols() != theLP->nCols())
137  	   {
138  	      SPxOut::debug(this, "IBASIS20 Dimension mismatch\n");
139  	      return false;
140  	   }
141  	
142  	   basisdim = 0;
143  	
144  	   for(int row = ds.nRows() - 1; row >= 0; --row)
145  	   {
146  	      if(ds.rowstat[row] >= 0)
147  	      {
148  	         if(ds.rowstat[row] != dualRowStatus(row))
149  	         {
150  	            SPxOut::debug(this, "IBASIS21 Basic row {} with incorrect dual status {}\n", row,
151  	                          ds.rowstat[row]);
152  	            return false;
153  	         }
154  	      }
155  	      else
156  	      {
157  	         basisdim++;
158  	
159  	         if((ds.rowstat[row] == Desc::P_FIXED
160  	               && theLP->SPxLPBase<R>::lhs(row) != theLP->SPxLPBase<R>::rhs(row))
161  	               || (ds.rowstat[row] == Desc::P_ON_UPPER && theLP->SPxLPBase<R>::rhs(row) >= R(infinity))
162  	               || (ds.rowstat[row] == Desc::P_ON_LOWER && theLP->SPxLPBase<R>::lhs(row) <= R(-infinity)))
163  	         {
164  	            SPxOut::debug(this, "IBASIS22 Nonbasic row with incorrect status: lhs={}, rhs={}, stat={}\n",
165  	                          theLP->SPxLPBase<R>::lhs(row), theLP->SPxLPBase<R>::rhs(row),
166  	                          ds.rowstat[row]);
167  	            return false;
168  	         }
169  	      }
170  	   }
171  	
172  	   for(int col = ds.nCols() - 1; col >= 0; --col)
173  	   {
174  	      if(ds.colstat[col] >= 0)
175  	      {
176  	         if(ds.colstat[col] !=  dualColStatus(col))
177  	         {
178  	            SPxOut::debug(this, "IBASIS23 Basic column {} with incorrect dual status {}\n", col,
179  	                          ds.colstat[col]);
180  	            return false;
181  	         }
182  	      }
183  	      else
184  	      {
185  	         basisdim++;
186  	
187  	         if((ds.colstat[col] == Desc::P_FIXED
188  	               && theLP->SPxLPBase<R>::lower(col) != theLP->SPxLPBase<R>::upper(col))
189  	               || (ds.colstat[col] == Desc::P_ON_UPPER && theLP->SPxLPBase<R>::upper(col) >= R(infinity))
190  	               || (ds.colstat[col] == Desc::P_ON_LOWER && theLP->SPxLPBase<R>::lower(col) <= R(-infinity)))
191  	         {
192  	            SPxOut::debug(this,
193  	                          "IBASIS24 Nonbasic column {} with incorrect status: lower={}, upper={}, stat={}\n",
194  	                          col, theLP->SPxLPBase<R>::lower(col), theLP->SPxLPBase<R>::upper(col),
195  	                          ds.colstat[col]);
196  	            return false;
197  	         }
198  	      }
199  	   }
200  	
201  	   if(basisdim != theLP->nCols())
202  	   {
203  	      SPxOut::debug(this, "IBASIS25 Incorrect basis dimension {} != {}\n", basisdim, theLP->nCols());
204  	      return false;
205  	   }
206  	
207  	   // basis descriptor valid
208  	   return true;
209  	}
210  	
211  	
212  	/*
213  	  Loading a #Desc# into the basis can be done more efficiently, by
214  	  explicitely programming both cases, for the rowwise and for the columnwise
215  	  representation. This implementation hides this distinction in the use of
216  	  methods #isBasic()# and #vector()#.
217  	*/
218  	template <class R>
219  	void SPxBasisBase<R>::loadDesc(const Desc& ds)
220  	{
221  	   assert(status() > NO_PROBLEM);
222  	   assert(theLP != 0);
223  	   assert(ds.nRows() == theLP->nRows());
224  	   assert(ds.nCols() == theLP->nCols());
225  	
226  	   SPxId none;
227  	   int   i;
228  	   int   j;
229  	   bool consistent = true;
230  	
231  	   SPX_MSG_INFO3((*this->spxout),
232  	                 (*this->spxout) << "IBASIS02 loading of Basis invalidates factorization"
233  	                 << std::endl;)
234  	
235  	   lastin      = none;
236  	   lastout     = none;
237  	   lastidx     = -1;
238  	   iterCount   = 0;
239  	   updateCount = 0;
240  	
241  	   if(&ds != &thedesc)
242  	   {
243  	      thedesc = ds;
244  	      setRep();
245  	   }
246  	
247  	   assert(theLP->dim() == matrix.size());
248  	
249  	   nzCount = 0;
250  	
251  	   for(j = i = 0; i < theLP->nRows(); ++i)
252  	   {
253  	      /* for columns and rows with D_... status, the correct D_... status depends on bounds and sides; if a basis
254  	       * descriptor is loaded after changing bounds or sides, e.g. in the refine() method, we have to correct them
255  	       */
256  	      if(thedesc.rowStatus(i) >= 0)
257  	         thedesc.rowStatus(i) = dualRowStatus(i);
258  	      else if(thedesc.rowStatus(i) == SPxBasisBase<R>::Desc::P_FIXED
259  	              && theLP->SPxLPBase<R>::lhs(i) != theLP->SPxLPBase<R>::rhs(i))
260  	      {
261  	         if(theLP->SPxLPBase<R>::lhs(i) > R(-infinity) && theLP->SPxLPBase<R>::maxRowObj(i) < 0.0)
262  	            thedesc.rowStatus(i) = SPxBasisBase<R>::Desc::P_ON_LOWER;
263  	         else if(theLP->SPxLPBase<R>::rhs(i) < R(infinity))
264  	            thedesc.rowStatus(i) = SPxBasisBase<R>::Desc::P_ON_UPPER;
265  	         else
266  	            thedesc.rowStatus(i) = SPxBasisBase<R>::Desc::P_FREE;
267  	      }
268  	
269  	      if(theLP->isBasic(thedesc.rowStatus(i)))
270  	      {
271  	         assert(theLP->dim() == matrix.size());
272  	         assert(j <= matrix.size());
273  	
274  	         if(j == matrix.size())
275  	         {
276  	            // too many basic variables
277  	            consistent = false;
278  	            break;
279  	         }
280  	
281  	         SPxRowId id = theLP->SPxLPBase<R>::rId(i);
282  	         theBaseId[j] = id;
283  	         matrix[j] = &theLP->vector(id);
284  	         nzCount += matrix[j++]->size();
285  	      }
286  	   }
287  	
288  	   for(i = 0; i < theLP->nCols(); ++i)
289  	   {
290  	      /* for columns and rows with D_... status, the correct D_... status depends on bounds and sides; if a basis
291  	       * descriptor is loaded after changing bounds or sides, e.g. in the refine() method, we have to correct them
292  	       */
293  	      if(thedesc.colStatus(i) >= 0)
294  	         thedesc.colStatus(i) = dualColStatus(i);
295  	      else if(thedesc.colStatus(i) == SPxBasisBase<R>::Desc::P_FIXED
296  	              && theLP->SPxLPBase<R>::lower(i) != theLP->SPxLPBase<R>::upper(i))
297  	      {
298  	         if(theLP->SPxLPBase<R>::lower(i) <= R(-infinity) && theLP->SPxLPBase<R>::upper(i) >= R(infinity))
299  	            thedesc.colStatus(i) = SPxBasisBase<R>::Desc::P_FREE;
300  	         else if(theLP->SPxLPBase<R>::upper(i) >= R(infinity)
301  	                 || (theLP->SPxLPBase<R>::lower(i) > R(-infinity) && theLP->SPxLPBase<R>::maxObj(i) < 0.0))
302  	            thedesc.colStatus(i) = SPxBasisBase<R>::Desc::P_ON_LOWER;
303  	         else
304  	            thedesc.colStatus(i) = SPxBasisBase<R>::Desc::P_ON_UPPER;
305  	      }
306  	
307  	      if(theLP->isBasic(thedesc.colStatus(i)))
308  	      {
309  	         assert(theLP->dim() == matrix.size());
310  	         assert(j <= matrix.size());
311  	
312  	         if(j == matrix.size())
313  	         {
314  	            // too many basic variables
315  	            consistent = false;
316  	            break;
317  	         }
318  	
319  	         SPxColId id = theLP->SPxLPBase<R>::cId(i);
320  	         theBaseId[j] = id;
321  	         matrix[j] = &theLP->vector(id);
322  	         nzCount += matrix[j++]->size();
323  	      }
324  	   }
325  	
326  	   if(j < matrix.size())
327  	      consistent = false;  // not enough basic variables
328  	
329  	   /* if dimensions are inconsistent, restore slack basis
330  	    * if dimensions are consistent, then we have setup a correct matrix
331  	    */
332  	   if(!consistent)
333  	      restoreInitialBasis();
334  	   else
335  	      matrixIsSetup = true;
336  	
337  	   assert(isDescValid(thedesc));
338  	
339  	   factorized = false;
340  	
341  	   if(factor != 0)
342  	      factor->clear();
343  	}
344  	
345  	template <class R>
346  	void SPxBasisBase<R>::setRep()
347  	{
348  	   assert(theLP != 0);
349  	
350  	   reDim();
351  	   minStab = 0.0;
352  	
353  	   if(theLP->rep() == SPxSolverBase<R>::ROW)
354  	   {
355  	      thedesc.stat   = &thedesc.rowstat;
356  	      thedesc.costat = &thedesc.colstat;
357  	   }
358  	   else
359  	   {
360  	      thedesc.stat   = &thedesc.colstat;
361  	      thedesc.costat = &thedesc.rowstat;
362  	   }
363  	}
364  	
365  	template <class R>
366  	void SPxBasisBase<R>::load(SPxSolverBase<R>* lp, bool initSlackBasis)
367  	{
368  	   assert(lp != 0);
369  	   theLP = lp;
370  	
371  	   setOutstream(*theLP->spxout);
372  	
373  	   setRep();
374  	
375  	   if(initSlackBasis)
376  	   {
377  	      restoreInitialBasis();
378  	      loadDesc(thedesc);
379  	   }
380  	}
381  	
382  	template <class R>
383  	void SPxBasisBase<R>::loadBasisSolver(SLinSolver<R>* p_solver, const bool destroy)
384  	{
385  	   assert(!freeSlinSolver || factor != 0);
386  	
387  	   setOutstream(*p_solver->spxout);
388  	
389  	   SPX_MSG_INFO3((*this->spxout),
390  	                 (*this->spxout) << "IBASIS03 loading of Solver invalidates factorization"
391  	                 << std::endl;)
392  	
393  	   if(freeSlinSolver)
394  	   {
395  	      delete factor;
396  	      factor = 0;
397  	   }
398  	
399  	   factor = p_solver;
400  	   factorized = false;
401  	   factor->clear();
402  	   freeSlinSolver = destroy;
403  	}
404  	
405  	/**
406  	 *  The specification is taken from the
407  	 *
408  	 *  ILOG CPLEX 7.0 Reference Manual, Appendix E, Page 543.
409  	 *
410  	 *  This routine should read valid BAS format files.
411  	 *
412  	 *  @return true if the file was read correctly.
413  	 *
414  	 *  Here is a very brief outline of the format:
415  	 *
416  	 *  The format is in a form similar to an MPS file. The basic assumption is that all (column)
417  	 *  variables are nonbasic at their lower bound and all row (variables) are basic; only the
418  	 *  differences to this rule are given. Each data line contains an indicator, a variable name and
419  	 *  possibly a row/constraint name. The following meaning applies with respect to the indicators:
420  	 *
421  	 *  - XU: the variable is basic, the row is nonbasic at its upper bound
422  	 *  - XL: the variable is basic, the row is nonbasic at its lower bound
423  	 *  - UL: the variable is nonbasic and at its upper bound
424  	 *  - LL: the variable is nonbasic and at its lower bound
425  	 *
426  	 *  The CPLEX format contains an additional indicator 'BS', but this is unsupported here.
427  	 *
428  	 *  Nonbasic variables without lower bound have the following default status for SoPlex:
429  	 *  - at their upper bound if finite,
430  	 *  - at zero if free.
431  	 */
432  	template <class R>
433  	bool SPxBasisBase<R>::readBasis(
434  	   std::istream&  is,
435  	   const NameSet* rowNames,
436  	   const NameSet* colNames)
437  	{
438  	   assert(theLP != 0);
439  	
440  	   /* prepare names */
441  	   const NameSet* rNames = rowNames;
442  	   const NameSet* cNames = colNames;
443  	
444  	   NameSet* p_colNames = 0;
445  	   NameSet* p_rowNames = 0;
446  	
447  	   if(colNames == 0)
448  	   {
449  	      int nCols = theLP->nCols();
450  	      std::stringstream name;
451  	
452  	      spx_alloc(p_colNames);
453  	      p_colNames = new(p_colNames) NameSet();
454  	      p_colNames->reMax(nCols);
455  	
456  	      for(int j = 0; j < nCols; ++j)
457  	      {
458  	         name << "x" << j;
459  	         DataKey key = theLP->colId(j);
460  	         p_colNames->add(key, name.str().c_str());
461  	      }
462  	
463  	      cNames = p_colNames;
464  	   }
465  	
466  	   if(rNames == 0)
467  	   {
468  	      int nRows = theLP->nRows();
469  	      std::stringstream name;
470  	
471  	      spx_alloc(p_rowNames);
472  	      p_rowNames = new(p_rowNames) NameSet();
473  	      p_rowNames->reMax(nRows);
474  	
475  	      for(int i = 0; i < nRows; ++i)
476  	      {
477  	         name << "C" << i;
478  	         DataKey key = theLP->rowId(i);
479  	         p_rowNames->add(key, name.str().c_str());
480  	      }
481  	
482  	      rNames = p_rowNames;
483  	   }
484  	
485  	   /* load default basis if necessary */
486  	   if(status() == NO_PROBLEM)
487  	      load(theLP, false);
488  	
489  	   /* initialize with standard settings */
490  	   Desc l_desc(thedesc);
491  	
492  	   for(int i = 0; i < theLP->nRows(); i++)
493  	      l_desc.rowstat[i] = dualRowStatus(i);
494  	
495  	   for(int i = 0; i < theLP->nCols(); i++)
496  	   {
497  	      if(theLP->SPxLPBase<R>::lower(i) == theLP->SPxLPBase<R>::upper(i))
498  	         l_desc.colstat[i] = Desc::P_FIXED;
499  	      else if(theLP->SPxLPBase<R>::lower(i) <= R(-infinity)
500  	              && theLP->SPxLPBase<R>::upper(i) >= R(infinity))
501  	         l_desc.colstat[i] = Desc::P_FREE;
502  	      else if(theLP->SPxLPBase<R>::lower(i) <= R(-infinity))
503  	         l_desc.colstat[i] = Desc::P_ON_UPPER;
504  	      else
505  	         l_desc.colstat[i] = Desc::P_ON_LOWER;
506  	   }
507  	
508  	   MPSInput mps(is);
509  	
510  	   if(mps.readLine() && (mps.field0() != 0) && !strcmp(mps.field0(), "NAME"))
511  	   {
512  	      while(mps.readLine())
513  	      {
514  	         int c = -1;
515  	         int r = -1;
516  	
517  	         if((mps.field0() != 0) && !strcmp(mps.field0(), "ENDATA"))
518  	         {
519  	            mps.setSection(MPSInput::ENDATA);
520  	            break;
521  	         }
522  	
523  	         if((mps.field1() == 0) || (mps.field2() == 0))
524  	            break;
525  	
526  	         if((c = cNames->number(mps.field2())) < 0)
527  	            break;
528  	
529  	         if(*mps.field1() == 'X')
530  	            if(mps.field3() == 0 || (r = rNames->number(mps.field3())) < 0)
531  	               break;
532  	
533  	         if(!strcmp(mps.field1(), "XU"))
534  	         {
535  	            l_desc.colstat[c] = dualColStatus(c);
536  	
537  	            if(theLP->LPRowSetBase<R>::type(r) == LPRowBase<R>::GREATER_EQUAL)
538  	               l_desc.rowstat[r] = Desc::P_ON_LOWER;
539  	            else if(theLP->LPRowSetBase<R>::type(r) == LPRowBase<R>::EQUAL)
540  	               l_desc.rowstat[r] = Desc::P_FIXED;
541  	            else
542  	               l_desc.rowstat[r] = Desc::P_ON_UPPER;
543  	         }
544  	         else if(!strcmp(mps.field1(), "XL"))
545  	         {
546  	            l_desc.colstat[c] = dualColStatus(c);
547  	
548  	            if(theLP->LPRowSetBase<R>::type(r) == LPRowBase<R>::LESS_EQUAL)
549  	               l_desc.rowstat[r] = Desc::P_ON_UPPER;
550  	            else if(theLP->LPRowSetBase<R>::type(r) == LPRowBase<R>::EQUAL)
551  	               l_desc.rowstat[r] = Desc::P_FIXED;
552  	            else
553  	               l_desc.rowstat[r] = Desc::P_ON_LOWER;
554  	         }
555  	         else if(!strcmp(mps.field1(), "UL"))
556  	         {
557  	            l_desc.colstat[c] = Desc::P_ON_UPPER;
558  	         }
559  	         else if(!strcmp(mps.field1(), "LL"))
560  	         {
561  	            l_desc.colstat[c] = Desc::P_ON_LOWER;
562  	         }
563  	         else
564  	         {
565  	            mps.syntaxError();
566  	            break;
567  	         }
568  	      }
569  	   }
570  	
571  	   if(!mps.hasError())
572  	   {
573  	      if(mps.section() == MPSInput::ENDATA)
574  	      {
575  	         // force basis to be different from NO_PROBLEM
576  	         // otherwise the basis will be overwritten at later stages.
577  	         setStatus(REGULAR);
578  	         loadDesc(l_desc);
579  	      }
580  	      else
581  	         mps.syntaxError();
582  	   }
583  	
584  	   if(rowNames == 0)
585  	   {
586  	      p_rowNames->~NameSet();
587  	      spx_free(p_rowNames);
588  	   }
589  	
590  	   if(colNames == 0)
591  	   {
592  	      p_colNames->~NameSet();
593  	      spx_free(p_colNames);
594  	   }
595  	
596  	#ifndef NDEBUG
597  	   SPX_DEBUG(thedesc.dump());
598  	#endif
599  	
600  	   return !mps.hasError();
601  	}
602  	
603  	
604  	/* Get row name - copied from spxmpswrite.cpp
605  	 *
606  	 * @todo put this in a common file and unify accross different formats (mps, lp, basis).
607  	 */
608  	template <class R>
609  	static const char* getRowName(
610  	   const SPxLPBase<R>*   lp,
611  	   int            idx,
612  	   const NameSet* rnames,
613  	   char*          buf)
614  	{
615  	   assert(buf != 0);
616  	   assert(idx >= 0);
617  	   assert(idx < lp->nRows());
618  	
619  	   if(rnames != 0)
620  	   {
621  	      DataKey key = lp->rId(idx);
622  	
623  	      if(rnames->has(key))
624  	         return (*rnames)[key];
625  	   }
626  	
627  	   spxSnprintf(buf, 16, "C%d", idx);
628  	
629  	   return buf;
630  	}
631  	
632  	/* Get column name - copied from spxmpswrite.cpp
633  	 *
634  	 * @todo put this in a common file and unify accross different formats (mps, lp, basis).
635  	 */
636  	template <class R>
637  	static const char* getColName(
638  	   const SPxLPBase<R>*   lp,
639  	   int            idx,
640  	   const NameSet* cnames,
641  	   char*          buf)
642  	{
643  	   assert(buf != 0);
644  	   assert(idx >= 0);
645  	   assert(idx < lp->nCols());
646  	
647  	   if(cnames != 0)
648  	   {
649  	      DataKey key = lp->cId(idx);
650  	
651  	      if(cnames->has(key))
652  	         return (*cnames)[key];
653  	   }
654  	
655  	   spxSnprintf(buf, 16, "x%d", idx);
656  	
657  	   return buf;
658  	}
659  	
660  	/* writes a file in MPS basis format to \p os.
661  	 *
662  	 * See SPxBasisBase<R>::readBasis() for a short description of the format.
663  	 */
664  	template <class R>
665  	void SPxBasisBase<R>::writeBasis(
666  	   std::ostream&  os,
667  	   const NameSet* rowNames,
668  	   const NameSet* colNames,
669  	   const bool cpxFormat
670  	) const
671  	{
672  	   assert(theLP != 0);
673  	
674  	   os.setf(std::ios::left);
675  	   os << "NAME  soplex.bas\n";
676  	
677  	   /* do not write basis if there is none */
678  	   if(status() == NO_PROBLEM)
679  	   {
680  	      os << "ENDATA" << std::endl;
681  	      return;
682  	   }
683  	
684  	   /* start writing */
685  	   char buf[255];
686  	   int row = 0;
687  	
688  	   for(int col = 0; col < theLP->nCols(); col++)
689  	   {
690  	      if(thedesc.colStatus(col) > 0)
691  	      {
692  	         /* Find non basic row */
693  	         for(; row < theLP->nRows(); row++)
694  	         {
695  	            if(thedesc.rowStatus(row) < 0)
696  	               break;
697  	         }
698  	
699  	         assert(row != theLP->nRows());
700  	
701  	         if(thedesc.rowStatus(row) == Desc::P_ON_UPPER && (!cpxFormat
702  	               || theLP->LPRowSetBase<R>::type(row) == LPRowBase<R>::RANGE))
703  	            os << " XU ";
704  	         else
705  	            os << " XL ";
706  	
707  	         os << std::setw(8) << getColName(theLP, col, colNames, buf);
708  	
709  	         /* break in two parts since buf is reused */
710  	         os << "       "
711  	            << getRowName(theLP, row, rowNames, buf)
712  	            << std::endl;
713  	
714  	         row++;
715  	      }
716  	      else
717  	      {
718  	         if(thedesc.colStatus(col) == Desc::P_ON_UPPER)
719  	         {
720  	            os << " UL "
721  	               << getColName(theLP, col, colNames, buf)
722  	               << std::endl;
723  	         }
724  	         else
725  	         {
726  	            /* Default is all non-basic variables on lower bound (if finite) or at zero (if free).
727  	             * nothing to do in this case.
728  	             */
729  	            assert(theLP->lower(col) <= R(-infinity) || thedesc.colStatus(col) == Desc::P_ON_LOWER
730  	                   || thedesc.colStatus(col) == Desc::P_FIXED);
731  	            assert(theLP->lower(col) > R(-infinity) || theLP->upper(col) < R(infinity)
732  	                   || thedesc.colStatus(col) == Desc::P_FREE);
733  	         }
734  	      }
735  	   }
736  	
737  	#ifndef NDEBUG
738  	   SPX_DEBUG(thedesc.dump());
739  	
740  	   // Check that we covered all nonbasic rows - the remaining should be basic.
741  	   for(; row < theLP->nRows(); row++)
742  	   {
743  	      if(thedesc.rowStatus(row) < 0)
744  	         break;
745  	   }
746  	
747  	   assert(row == theLP->nRows());
748  	
749  	#endif // NDEBUG
750  	
751  	   os << "ENDATA" << std::endl;
752  	}
753  	
754  	template <class R>
755  	void SPxBasisBase<R>::printMatrix() const
756  	{
757  	
758  	   assert(matrixIsSetup);
759  	
760  	   for(int i = 0; i < matrix.size(); i++)
761  	   {
762  	      std::cout << "C" << i << "=" << *matrix[i] << std::endl;
763  	   }
764  	}
765  	
766  	template <class R>
767  	void SPxBasisBase<R>::printMatrixMTX(int number)
768  	{
769  	   int dim;
770  	   int nnz;
771  	   char filename[SPX_MAXSTRLEN];
772  	
773  	   dim = matrix.size();
774  	   nnz = nzCount;
775  	   spxSnprintf(filename, SPX_MAXSTRLEN, "basis/basis%d.mtx", number);
776  	   std::cout << "printing basis matrix to file " << filename << "\n";
777  	   FILE* basisfile;
778  	   basisfile = fopen(filename, "w");
779  	   // print marker necessary for reading the file in Matlab
780  	   fprintf(basisfile, "%%%%MatrixMarket matrix coordinate Real general\n");
781  	   // print matrix information
782  	   fprintf(basisfile, "%d %d %d\n", dim, dim, nnz);
783  	
784  	   // print matrix data
785  	   for(int i = 0; i < matrix.size(); ++i)
786  	   {
787  	      for(int j = 0; j < baseVec(i).size(); ++j)
788  	      {
789  	         int idx = baseVec(i).index(j);
790  	         R val = baseVec(i).value(j);
791  	         fprintf(basisfile, "%d %d %.13" SOPLEX_REAL_FORMAT "\n", i + 1, idx + 1, val);
792  	      }
793  	   }
794  	
795  	   fclose(basisfile);
796  	
797  	   return;
798  	}
799  	
800  	template <class R>
801  	void SPxBasisBase<R>::change(
802  	   int i,
803  	   SPxId& id,
804  	   const SVectorBase<R>* enterVec,
805  	   const SSVectorBase<R>* eta)
806  	{
807  	
808  	   assert(matrixIsSetup);
809  	   assert(!id.isValid() || (enterVec != 0));
810  	   assert(factor != 0);
811  	
812  	   lastidx = i;
813  	   lastin  = id;
814  	
815  	   if(id.isValid() && i >= 0)
816  	   {
817  	      assert(enterVec != 0);
818  	
819  	      // update the counter for nonzeros in the basis matrix
820  	      nzCount      = nzCount - matrix[i]->size() + enterVec->size();
821  	      // let the new id enter the basis
822  	      matrix[i]    = enterVec;
823  	      lastout      = theBaseId[i];
824  	      theBaseId[i] = id;
825  	
826  	      ++iterCount;
827  	      ++updateCount;
828  	
829  	      SPxOut::debug(this,
830  	                    "factor_stats: iteration= {}, update= {}, total_update= {}, nonzero_B= {}, nonzero_LU= {}, factor_fill= {}, time= {}\n",
831  	                    this->iteration(), updateCount, totalUpdateCount, nzCount, factor->memory(), lastFill,
832  	                    theLP->time());
833  	
834  	      // never factorize? Just do it !
835  	      if(!factorized)
836  	         factorize();
837  	
838  	      // too much memory growth ?
839  	      else if(R(factor->memory()) > 1000 + factor->dim() + lastMem * memFactor)
840  	      {
841  	         SPX_MSG_INFO3((*this->spxout), (*this->spxout) <<
842  	                       "IBASIS04 memory growth factor triggers refactorization"
843  	                       << " memory= " << factor->memory()
844  	                       << " lastMem= " << lastMem
845  	                       << " memFactor= " << memFactor
846  	                       << std::endl;)
847  	         factorize();
848  	      }
849  	
850  	      // relative fill too high ?
851  	      else if(R(factor->memory()) > lastFill * R(nzCount))
852  	      {
853  	         SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "IBASIS04 fill factor triggers refactorization"
854  	                       << " memory= " << factor->memory()
855  	                       << " nzCount= " << nzCount
856  	                       << " lastFill= " << lastFill
857  	                       << std::endl;)
858  	
859  	         factorize();
860  	      }
861  	      // absolute fill in basis matrix too high ?
862  	      else if(nzCount > lastNzCount)
863  	      {
864  	         SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "IBASIS05 nonzero factor triggers refactorization"
865  	                       << " nzCount= " << nzCount
866  	                       << " lastNzCount= " << lastNzCount
867  	                       << " nonzeroFactor= " << nonzeroFactor
868  	                       << std::endl;)
869  	         factorize();
870  	      }
871  	      // too many updates ?
872  	      else if(updateCount >= maxUpdates)
873  	      {
874  	         SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "IBASIS06 update count triggers refactorization"
875  	                       << " updateCount= " << updateCount
876  	                       << " maxUpdates= " << maxUpdates
877  	                       << std::endl;)
878  	         factorize();
879  	      }
880  	      else
881  	      {
882  	         try
883  	         {
884  	#ifdef MEASUREUPDATETIME
885  	            theTime.start();
886  	#endif
887  	            factor->change(i, *enterVec, eta);
888  	            totalUpdateCount++;
889  	#ifdef MEASUREUPDATETIME
890  	            theTime.stop();
891  	#endif
892  	         }
893  	         catch(...)
894  	         {
895  	            SPX_MSG_INFO3((*this->spxout), (*this->spxout) <<
896  	                          "IBASIS13 problems updating factorization; refactorizing basis"
897  	                          << std::endl;)
898  	
899  	#ifdef MEASUREUPDATETIME
900  	            theTime.stop();
901  	#endif
902  	
903  	            // singularity was detected in update; we refactorize
904  	            factorize();
905  	
906  	            // if factorize() detects singularity, an exception is thrown, hence at this point we have a regular basis
907  	            // and can try the update again
908  	            assert(status() >= SPxBasisBase<R>::REGULAR);
909  	
910  	            try
911  	            {
912  	#ifdef MEASUREUPDATETIME
913  	               theTime.start();
914  	#endif
915  	               factor->change(i, *enterVec, eta);
916  	               totalUpdateCount++;
917  	#ifdef MEASUREUPDATETIME
918  	               theTime.stop();
919  	#endif
920  	            }
921  	            // with a freshly factorized, regular basis, the update is unlikely to fail; if this happens nevertheless,
922  	            // we have to invalidate the basis to have the statuses correct
923  	            catch(const SPxException& F)
924  	            {
925  	               SPX_MSG_INFO3((*this->spxout), (*this->spxout) <<
926  	                             "IBASIS14 problems updating factorization; invalidating factorization"
927  	                             << std::endl;)
928  	
929  	#ifdef MEASUREUPDATETIME
930  	               theTime.stop();
931  	#endif
932  	
933  	               factorized = false;
934  	               throw F;
935  	            }
936  	         }
937  	
938  	         assert(minStab > 0.0);
939  	
940  	         if(factor->status() != SLinSolver<R>::OK || factor->stability() < minStab)
941  	         {
942  	            SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "IBASIS07 stability triggers refactorization"
943  	                          << " stability= " << factor->stability()
944  	                          << " minStab= " << minStab
945  	                          << std::endl;)
946  	            factorize();
947  	         }
948  	      }
949  	   }
950  	   else
951  	      lastout = id;
952  	}
953  	
954  	template <class R>
955  	void SPxBasisBase<R>::factorize()
956  	{
957  	
958  	   assert(factor != 0);
959  	
960  	   if(!matrixIsSetup)
961  	      loadDesc(thedesc);
962  	
963  	   assert(matrixIsSetup);
964  	
965  	   updateCount = 0;
966  	
967  	   switch(factor->load(matrix.get_ptr(), matrix.size()))
968  	   {
969  	   case SLinSolver<R>::OK :
970  	      if(status() == SINGULAR)
971  	         setStatus(REGULAR);
972  	
973  	      factorized = true;
974  	      minStab = factor->stability();
975  	
976  	      // This seems always to be about 1e-7
977  	      if(minStab > 1e-4)
978  	         minStab *= 0.001;
979  	
980  	      if(minStab > 1e-5)
981  	         minStab *= 0.01;
982  	
983  	      if(minStab > 1e-6)
984  	         minStab *= 0.1;
985  	
986  	      break;
987  	
988  	   case SLinSolver<R>::SINGULAR :
989  	      setStatus(SINGULAR);
990  	      factorized = false;
991  	      break;
992  	
993  	   default :
994  	      SPX_MSG_ERROR(std::cerr << "EBASIS08 error: unknown status of factorization.\n";)
995  	      factorized = false;
996  	      throw SPxInternalCodeException("XBASIS01 This should never happen.");
997  	   }
998  	
999  	   // get nonzero count of factorization
1000 	   lastMem    = factor->memory();
1001 	   // compute fill ratio between factorization and basis matrix (multiplied with tolerance)
1002 	   lastFill   = fillFactor * R(lastMem) / R(nzCount > 0 ? nzCount : 1);
1003 	   lastNzCount = int(nonzeroFactor * R(nzCount > 0 ? nzCount : 1));
1004 	
1005 	   if(status() == SINGULAR)
1006 	   {
1007 	      throw SPxStatusException("Cannot factorize singular matrix");
1008 	   }
1009 	}
1010 	
1011 	template <class R>
1012 	VectorBase<R>& SPxBasisBase<R>::multWithBase(VectorBase<R>& x) const
1013 	{
1014 	   assert(status() > SINGULAR);
1015 	   assert(theLP->dim() == x.dim());
1016 	
1017 	   int i;
1018 	   VectorBase<R> tmp(x);
1019 	
1020 	   if(!matrixIsSetup)
1021 	      (const_cast<SPxBasisBase<R>*>(this))->loadDesc(thedesc);
1022 	
1023 	   assert(matrixIsSetup);
1024 	
1025 	   for(i = x.dim() - 1; i >= 0; --i)
1026 	      x[i] = *(matrix[i]) * tmp;
1027 	
1028 	   return x;
1029 	}
1030 	
1031 	template <class R>
1032 	void SPxBasisBase<R>::multWithBase(SSVectorBase<R>& x, SSVectorBase<R>& result) const
1033 	{
1034 	   assert(status() > SINGULAR);
1035 	   assert(theLP->dim() == x.dim());
1036 	   assert(x.dim() == result.dim());
1037 	
1038 	   if(!matrixIsSetup)
1039 	      (const_cast<SPxBasisBase<R>*>(this))->loadDesc(thedesc);
1040 	
1041 	   result.clear();
1042 	
1043 	   assert(matrixIsSetup);
1044 	
1045 	   for(int i = 0; i < x.dim(); ++i)
1046 	      result.add(i, (*matrix[i]) * x);
1047 	
1048 	   return;
1049 	}
1050 	
1051 	template <class R>
1052 	VectorBase<R>& SPxBasisBase<R>::multBaseWith(VectorBase<R>& x) const
1053 	{
1054 	   assert(status() > SINGULAR);
1055 	   assert(theLP->dim() == x.dim());
1056 	
1057 	   int i;
1058 	   VectorBase<R> tmp(x);
1059 	
1060 	   if(!matrixIsSetup)
1061 	      (const_cast<SPxBasisBase<R>*>(this))->loadDesc(thedesc);
1062 	
1063 	   assert(matrixIsSetup);
1064 	
1065 	   x.clear();
1066 	
1067 	   for(i = x.dim() - 1; i >= 0; --i)
1068 	   {
1069 	      if(tmp[i] != 0.0)
1070 	         x.multAdd(tmp[i], *(matrix[i]));
1071 	   }
1072 	
1073 	   return x;
1074 	}
1075 	
1076 	template <class R>
1077 	void SPxBasisBase<R>::multBaseWith(SSVectorBase<R>& x, SSVectorBase<R>& result) const
1078 	{
1079 	   assert(status() > SINGULAR);
1080 	   assert(theLP->dim() == x.dim());
1081 	   assert(x.dim() == result.dim());
1082 	
1083 	   if(!matrixIsSetup)
1084 	      (const_cast<SPxBasisBase<R>*>(this))->loadDesc(thedesc);
1085 	
1086 	   assert(matrixIsSetup);
1087 	
1088 	   result.clear();
1089 	
1090 	   if(x.isSetup())
1091 	   {
1092 	      for(int i = 0; i < x.size(); ++i)
1093 	      {
1094 	         int idx = x.index(i);
1095 	         result.multAdd(x[idx], (*matrix[idx]));
1096 	      }
1097 	   }
1098 	   else
1099 	   {
1100 	      for(int i = 0; i < x.dim(); ++i)
1101 	         result.multAdd(x[i], (*matrix[i]));
1102 	   }
1103 	
1104 	   return;
1105 	}
1106 	
1107 	template <class R>
1108 	/* compute an estimated condition number for the current basis matrix
1109 	 * by computing estimates of the norms of B and B^-1 using the power method.
1110 	 * maxiters and tolerance control the accuracy of the estimate.
1111 	 */
1112 	R SPxBasisBase<R>::condition(int maxiters, R tolerance)
1113 	{
1114 	   int dimension = matrix.size();
1115 	   int miniters = 3;    // minimum number of power method iterations
1116 	   int i;
1117 	   int c;
1118 	   R norm;
1119 	   R norminv;
1120 	   R norm1;
1121 	   R norm2;
1122 	
1123 	   // catch corner case of empty matrix
(1) Event cond_false: Condition "dimension <= 0", taking false branch.
1124 	   if(dimension <= 0)
(2) Event if_end: End of if statement.
1125 	      return 1.0;
1126 	
(3) Event write_zero_model: "SSVectorBase" sets "x.idx" to "nullptr". [details]
Also see events: [var_deref_model]
1127 	   SSVectorBase<R> x(dimension, theLP->tolerances());
1128 	   SSVectorBase<R> y(dimension, theLP->tolerances());
1129 	
1130 	   // check whether a regular basis matrix is available
(4) Event cond_false: Condition "this->status() < soplex::SPxBasisBase<double>::REGULAR", taking false branch.
1131 	   if(status() < REGULAR)
(5) Event if_end: End of if statement.
1132 	      return 0;
1133 	
(6) Event cond_true: Condition "!this->matrixIsSetup", taking true branch.
1134 	   if(!matrixIsSetup)
1135 	      (const_cast<SPxBasisBase<R>*>(this))->loadDesc(thedesc);
1136 	
(7) Event cond_true: Condition "!this->factorized", taking true branch.
1137 	   if(!factorized)
1138 	      factorize();
1139 	
1140 	   // initialize vectors
1141 	   norm1 = 1.0 / (R) dimension;
1142 	
(8) Event cond_true: Condition "i < dimension", taking true branch.
1143 	   for(i = 0; i < dimension; i++)
(9) Event var_deref_model: "add" dereferences null "x->idx". [details]
Also see events: [write_zero_model]
1144 	      x.add(i, norm1);
1145 	
1146 	   y = x;
1147 	
1148 	   // compute norm of B
1149 	   for(c = 0; c < maxiters; ++c)
1150 	   {
1151 	      norm2 = norm1;
1152 	
1153 	      // y = B*x
1154 	      multBaseWith(x, y);
1155 	      norm1 = y.length();
1156 	
1157 	      // stop if converged
1158 	      if(c >= miniters && spxAbs(norm1 - norm2) < tolerance * norm1)
1159 	         break;
1160 	
1161 	      // x = B^T*y and normalize
1162 	      multWithBase(y, x);
1163 	      norm2 = 1.0 / x.length();
1164 	      x *= norm2;
1165 	   }
1166 	
1167 	   norm = norm1;
1168 	
1169 	   // reinitialize vectors
1170 	   x.clear();
1171 	   y.clear();
1172 	   norm1 = 1.0 / (R) dimension;
1173 	
1174 	   for(i = 0; i < dimension; i++)
1175 	      x.add(i, norm1);
1176 	
1177 	   y = x;
1178 	
1179 	   // compute norm of B^-1
1180 	   for(c = 0; c < maxiters; ++c)
1181 	   {
1182 	      norm2 = norm1;
1183 	
1184 	      // x = B^-1*y
1185 	      factor->solveRight(x, y);
1186 	      x.setup();
1187 	      norm1 = x.length();
1188 	
1189 	      // stop if converged
1190 	      if(c >= miniters && spxAbs(norm1 - norm2) < tolerance * norm1)
1191 	         break;
1192 	
1193 	      // y = B^-T*x and normalize
1194 	      factor->solveLeft(y, x);
1195 	      y.setup();
1196 	      norm2 = 1.0 / y.length();
1197 	      y *= norm2;
1198 	   }
1199 	
1200 	   norminv = norm1;
1201 	
1202 	   return norm * norminv;
1203 	}
1204 	
1205 	/* compute one of several matrix metrics based on the diagonal of the LU factorization */
1206 	template <class R>
1207 	R SPxBasisBase<R>::getMatrixMetric(int type)
1208 	{
1209 	   R metric = R(infinity);
1210 	
1211 	   if(factorized)
1212 	      metric = factor->matrixMetric(type);
1213 	
1214 	   return metric;
1215 	}
1216 	
1217 	template <class R>
1218 	void SPxBasisBase<R>::dump()
1219 	{
1220 	   assert(status() > NO_PROBLEM);
1221 	   assert(theLP != 0);
1222 	   assert(thedesc.nRows() == theLP->nRows());
1223 	   assert(thedesc.nCols() == theLP->nCols());
1224 	   assert(theLP->dim() == matrix.size());
1225 	
1226 	   int i, basesize;
1227 	
1228 	   // Dump regardless of the verbosity level if this method is called.
1229 	
1230 	   std::cout << "DBASIS09 Basis entries:";
1231 	   basesize = 0;
1232 	
1233 	   for(i = 0; i < theLP->nRows(); ++i)
1234 	   {
1235 	      if(theLP->isBasic(thedesc.rowStatus(i)))
1236 	      {
1237 	         if(basesize % 10 == 0)
1238 	            std::cout << std::endl << "DBASIS10 ";
1239 	
1240 	         SPxRowId id = theLP->SPxLPBase<R>::rId(i);
1241 	         std::cout << "\tR" << theLP->number(id);
1242 	         basesize++;
1243 	      }
1244 	   }
1245 	
1246 	   for(i = 0; i < theLP->nCols(); ++i)
1247 	   {
1248 	      if(theLP->isBasic(thedesc.colStatus(i)))
1249 	      {
1250 	         if(basesize % 10 == 0)
1251 	            std::cout << std::endl << "DBASIS11 ";
1252 	
1253 	         SPxColId id = theLP->SPxLPBase<R>::cId(i);
1254 	         std::cout << "\tC" << theLP->number(id);
1255 	         basesize++;
1256 	      }
1257 	   }
1258 	
1259 	   std::cout << std::endl;
1260 	
1261 	   assert(basesize == matrix.size());
1262 	}
1263 	
1264 	template <class R>
1265 	
1266 	bool SPxBasisBase<R>::isConsistent() const
1267 	{
1268 	#ifdef ENABLE_CONSISTENCY_CHECKS
1269 	   int primals = 0;
1270 	   int i;
1271 	
1272 	   if(status() > NO_PROBLEM)
1273 	   {
1274 	      if(theLP == 0)
1275 	         return SPX_MSG_INCONSISTENT("SPxBasisBase<R>");
1276 	
1277 	      if(theBaseId.size() != theLP->dim() || matrix.size() != theLP->dim())
1278 	         return SPX_MSG_INCONSISTENT("SPxBasisBase<R>");
1279 	
1280 	      if(thedesc.nCols() != theLP->nCols() || thedesc.nRows() != theLP->nRows())
1281 	         return SPX_MSG_INCONSISTENT("SPxBasisBase<R>");
1282 	
1283 	      for(i = 0; i < thedesc.nRows(); ++i)
1284 	      {
1285 	         if(thedesc.rowStatus(i) >= 0)
1286 	         {
1287 	            if(thedesc.rowStatus(i) != dualRowStatus(i))
1288 	               return SPX_MSG_INCONSISTENT("SPxBasisBase<R>");
1289 	         }
1290 	         else
1291 	            ++primals;
1292 	      }
1293 	
1294 	      for(i = 0; i < thedesc.nCols(); ++i)
1295 	      {
1296 	         if(thedesc.colStatus(i) >= 0)
1297 	         {
1298 	            if(thedesc.colStatus(i) != dualColStatus(i))
1299 	               return SPX_MSG_INCONSISTENT("SPxBasisBase<R>");
1300 	         }
1301 	         else
1302 	            ++primals;
1303 	      }
1304 	
1305 	      if(primals != thedesc.nCols())
1306 	         return SPX_MSG_INCONSISTENT("SPxBasisBase<R>");
1307 	   }
1308 	
1309 	   return thedesc.isConsistent() && theBaseId.isConsistent()
1310 	          && matrix.isConsistent() && factor->isConsistent();
1311 	#else
1312 	   return true;
1313 	#endif // CONSISTENCY_CHECKS
1314 	}
1315 	
1316 	template <class R>
1317 	SPxBasisBase<R>::SPxBasisBase(Timer::TYPE ttype)
1318 	   : theLP(0)
1319 	   , matrixIsSetup(false)
1320 	   , factor(0)
1321 	   , factorized(false)
1322 	   , maxUpdates(200)
1323 	   , nonzeroFactor(10.0)
1324 	   , fillFactor(5.0)
1325 	   , memFactor(1.5)
1326 	   , iterCount(0)
1327 	   , lastIterCount(0)
1328 	   , iterDegenCheck(0)
1329 	   , updateCount(0)
1330 	   , totalUpdateCount(0)
1331 	   , nzCount(1)
1332 	   , lastMem(0)
1333 	   , lastFill(0)
1334 	   , lastNzCount(0)
1335 	   , theTime(0)
1336 	   , timerType(ttype)
1337 	   , lastidx(0)
1338 	   , minStab(0.0)
1339 	   , thestatus(NO_PROBLEM)
1340 	   , freeSlinSolver(false)
1341 	   , spxout(0)
1342 	{
1343 	   // info: is not consistent at this moment, e.g. because theLP == 0
1344 	
1345 	   theTime = TimerFactory::createTimer(timerType);
1346 	}
1347 	
1348 	
1349 	/**@warning Do not change the LP object.
1350 	 *  Only pointer to that object is copied.
1351 	 *  Hint: no problem, we use this function for copy
1352 	 *   constructor of SPxSolverBase<R>
1353 	 */
1354 	template <class R>
1355 	SPxBasisBase<R>::SPxBasisBase(const SPxBasisBase<R>& old)
1356 	   : theLP(old.theLP)
1357 	   , theBaseId(old.theBaseId)
1358 	   , matrix(old.matrix)
1359 	   , matrixIsSetup(old.matrixIsSetup)
1360 	   , factor(old.factor)
1361 	   , factorized(old.factorized)
1362 	   , maxUpdates(old.maxUpdates)
1363 	   , nonzeroFactor(old.nonzeroFactor)
1364 	   , fillFactor(old.fillFactor)
1365 	   , memFactor(old.memFactor)
1366 	   , iterCount(old.iterCount)
1367 	   , lastIterCount(old.lastIterCount)
1368 	   , iterDegenCheck(old.iterDegenCheck)
1369 	   , updateCount(old.updateCount)
1370 	   , totalUpdateCount(old.totalUpdateCount)
1371 	   , nzCount(old.nzCount)
1372 	   , lastMem(old.lastMem)
1373 	   , lastFill(old.lastFill)
1374 	   , lastNzCount(old.lastNzCount)
1375 	   , theTime(old.theTime)
1376 	   , lastin(old.lastin)
1377 	   , lastout(old.lastout)
1378 	   , lastidx(old.lastidx)
1379 	   , minStab(old.minStab)
1380 	   , thestatus(old.thestatus)
1381 	   , thedesc(old.thedesc)
1382 	   , spxout(old.spxout)
1383 	{
1384 	   theTime = TimerFactory::createTimer(old.theTime->type());
1385 	
1386 	   this->factor = old.factor->clone();
1387 	   freeSlinSolver = true;
1388 	
1389 	   assert(SPxBasisBase<R>::isConsistent());
1390 	}
1391 	
1392 	template <class R>
1393 	SPxBasisBase<R>::~SPxBasisBase()
1394 	{
1395 	
1396 	   assert(!freeSlinSolver || factor != 0);
1397 	
1398 	   if(freeSlinSolver)
1399 	   {
1400 	      delete factor;
1401 	      factor = 0;
1402 	   }
1403 	
1404 	   theTime->~Timer();
1405 	   spx_free(theTime);
1406 	}
1407 	
1408 	template <class R>
1409 	
1410 	/**@warning  Note that we do not create a deep copy of the corresponding SPxSolverBase<R> object.
1411 	 *  Only the reference to that object is copied.
1412 	 */
1413 	SPxBasisBase<R>& SPxBasisBase<R>::operator=(const SPxBasisBase<R>& rhs)
1414 	{
1415 	
1416 	   assert(!freeSlinSolver || factor != 0);
1417 	
1418 	   if(this != &rhs)
1419 	   {
1420 	      theLP         = rhs.theLP;
1421 	      theBaseId     = rhs.theBaseId;
1422 	      matrix        = rhs.matrix;
1423 	      matrixIsSetup = rhs.matrixIsSetup;
1424 	
1425 	      if(freeSlinSolver)
1426 	      {
1427 	         delete factor;
1428 	         factor = 0;
1429 	      }
1430 	
1431 	      factor = rhs.factor->clone();
1432 	      freeSlinSolver = true;
1433 	
1434 	      factorized    = rhs.factorized;
1435 	      maxUpdates    = rhs.maxUpdates;
1436 	      nonzeroFactor = rhs.nonzeroFactor;
1437 	      fillFactor    = rhs.fillFactor;
1438 	      memFactor     = rhs.memFactor;
1439 	      iterCount     = rhs.iterCount;
1440 	      nzCount       = rhs.nzCount;
1441 	      lastFill      = rhs.lastFill;
1442 	      lastNzCount   = rhs.lastNzCount;
1443 	      lastin        = rhs.lastin;
1444 	      lastout       = rhs.lastout;
1445 	      lastidx       = rhs.lastidx;
1446 	      minStab       = rhs.minStab;
1447 	      thestatus     = rhs.thestatus;
1448 	      thedesc       = rhs.thedesc;
1449 	
1450 	      assert(SPxBasisBase<R>::isConsistent());
1451 	   }
1452 	
1453 	   return *this;
1454 	}
1455 	
1456 	
1457 	
1458 	//
1459 	// Auxiliary functions.
1460 	//
1461 	
1462 	// Pretty-printing of basis status.
1463 	template <class R> // Why can't we remove the class R and make it empty?
1464 	std::ostream& operator<<(std::ostream& os,
1465 	                         const typename SPxBasisBase<R>::SPxStatus& status)
1466 	{
1467 	   switch(status)
1468 	   {
1469 	   case SPxBasisBase<R>::NO_PROBLEM:
1470 	      os << "NO_PROBLEM";
1471 	      break;
1472 	
1473 	   case SPxBasisBase<R>::SINGULAR:
1474 	      os << "SINGULAR";
1475 	      break;
1476 	
1477 	   case SPxBasisBase<R>::REGULAR:
1478 	      os << "REGULAR";
1479 	      break;
1480 	
1481 	   case SPxBasisBase<R>::DUAL:
1482 	      os << "DUAL";
1483 	      break;
1484 	
1485 	   case SPxBasisBase<R>::PRIMAL:
1486 	      os << "PRIMAL";
1487 	      break;
1488 	
1489 	   case SPxBasisBase<R>::OPTIMAL:
1490 	      os << "OPTIMAL";
1491 	      break;
1492 	
1493 	   case SPxBasisBase<R>::UNBOUNDED:
1494 	      os << "UNBOUNDED";
1495 	      break;
1496 	
1497 	   case SPxBasisBase<R>::INFEASIBLE:
1498 	      os << "INFEASIBLE";
1499 	      break;
1500 	
1501 	   default:
1502 	      os << "?other?";
1503 	      break;
1504 	   }
1505 	
1506 	   return os;
1507 	}
1508 	
1509 	
1510 	} // namespace soplex
1511