/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
	/*                                                                           */
	/*                  This file is part of the class library                   */
	/*       SoPlex --- the Sequential object-oriented simPlex.                  */
	/*                                                                           */
	/*    Copyright (C) 1996-2022 Konrad-Zuse-Zentrum                            */
	/*                            fuer Informationstechnik Berlin                */
	/*                                                                           */
	/*  SoPlex is distributed under the terms of the ZIB Academic Licence.       */
	/*                                                                           */
	/*  You should have received a copy of the ZIB Academic License              */
	/*  along with SoPlex; see the file COPYING. If not email to soplex@zib.de.  */
	/*                                                                           */
	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
	
	/**@file  spxlpbase.h
	 * @brief Saving LPs in a form suitable for SoPlex.
	 */
	#ifndef _SPXLPBASE_H_
	#define _SPXLPBASE_H_
	
	/* undefine SOPLEX_DEBUG flag from including files; if SOPLEX_DEBUG should be defined in this file, do so below */
	#ifdef SOPLEX_DEBUG
	#define SOPLEX_DEBUG_SPXLPBASE
	#undef SOPLEX_DEBUG
	#endif
	
	#include <assert.h>
	#include <iostream>
	#include <iomanip>
	#include <typeinfo>
	
	#include "soplex/spxdefines.h"
	#include "soplex/basevectors.h"
	#include "soplex/dataarray.h"
	#include "soplex/datakey.h"
	#include "soplex/spxid.h"
	#include "soplex/lprowbase.h"
	#include "soplex/lpcolbase.h"
	#include "soplex/lprowsetbase.h"
	#include "soplex/lpcolsetbase.h"
	#include "soplex/nameset.h"
	#include "soplex/didxset.h"
	#include "soplex/spxfileio.h"
	#include "soplex/spxscaler.h"
	#include "soplex/rational.h"
	
	namespace soplex
	{
	// Declarations to fix errors of the form "SPxMainSM is not a type"
	template <class R>
	class SPxSolverBase;
	template <class R>
	class SPxMainSM;
	template <class R>
	class SPxLPBase;
	template <class R>
	class SPxBasisBase;
	
	template <class R>
	class SPxEquiliSC;
	template <class R>
	class SPxLeastSqSC;
	template <class R>
	class SPxGeometSC;
	template <class R>
	class SPxMainSM;
	
	
	/**@brief   Saving LPs in a form suitable for SoPlex.
	 * @ingroup Algo
	 *
	 *  Class SPxLPBase provides the data structures required for saving a linear program in the form
	 *  \f[
	 *  \begin{array}{rl}
	 *      \hbox{max}  & c^T x              \\
	 *      \hbox{s.t.} & l_r \le Ax \le u_r \\
	 *                  & l_c \le x \le u_c
	 *  \end{array}
	 *  \f]
	 *  suitable for solving with SoPlex. This includes:
	 *  - SVSetBase%s for both columns and rows
	 *  - objective Vector
	 *  - upper and lower bound Vectors for variables (\f$l_c\f$ and \f$u_c\f$)
	 *  - upper and lower bound Vectors for inequalities (\f$l_r\f$ and \f$u_r\f$)
	 *
	 *  Note, that the optimization sense is not saved directly. Instead, the objective function are multiplied by -1 to
	 *  transform the LP to our standard form maximizing the objective function. However, the sense of the loaded LP can be
	 *  retrieved with method #spxSense().
	 *
	 *  Further, equality constraints are modeled by \f$l_r = u_r\f$.  Analogously, fixed variables have \f$l_c = u_c\f$.
	 *
	 *  #SPxLPBase%s are saved as an SVSet, both for columns and rows. Note that this is redundant but eases the access.
	 */
	
	
	template <class R>
	class SPxLPBase : protected LPRowSetBase<R>, protected LPColSetBase<R>
	{
	   template <class S> friend class SPxLPBase;
	   friend SPxBasisBase<R>;
	   friend SPxScaler<R>;
	   friend SPxEquiliSC<R>;
	   friend SPxLeastSqSC<R>;
	   friend SPxGeometSC<R>;
	   friend SPxMainSM<R>;
	
	public:
	
	   // ------------------------------------------------------------------------------------------------------------------
	   /**@name Types */
	   ///@{
	
	   /// Optimization sense.
	   enum SPxSense
	   {
	      MAXIMIZE = 1,
	      MINIMIZE = -1
	   };
	
	   ///@}
	
	private:
	
	   // ------------------------------------------------------------------------------------------------------------------
	   /**@name Data */
	   ///@{
	
	   SPxSense thesense;                ///< optimization sense.
	   R offset;                         ///< offset computed, e.g., in simplification step
	   bool _isScaled;                   ///< true, if scaling has been performed
	   SPxScaler<R>*
	   lp_scaler;             ///< points to the scaler if the lp has been scaled, to nullptr otherwise
	
	   ///@}
	
	public:
	
	   // message handler
	   SPxOut* spxout;
	
	public:
	
	   void setOutstream(SPxOut& newOutstream)
	   {
	      spxout = &newOutstream;
	   }
	
	   // ------------------------------------------------------------------------------------------------------------------
	
	   /// unscales the lp and clears basis
	   void unscaleLP();
	
	   /**@name Inquiry */
	   ///@{
	
	   /// Returns true if and only if the LP is scaled
	   bool isScaled() const
	   {
	      return _isScaled;
	   }
	
	   /// set whether the LP is scaled or not
	   void setScalingInfo(bool scaled)
	   {
	      _isScaled = scaled;
	   }
	
	   /// Returns number of rows in LP.
	   int nRows() const
	   {
	      return LPRowSetBase<R>::num();
	   }
	
	   /// Returns number of columns in LP.
	   int nCols() const
	   {
	      return LPColSetBase<R>::num();
	   }
	
	   /// Returns number of nonzeros in LP.
	   int nNzos() const
	   {
	
	      int n = 0;
	
	      for(int i = 0; i < nCols(); ++i)
	         n += colVector(i).size();
	
	      return n;
	   }
	
	   /// Absolute smallest non-zero element in (possibly scaled) LP.
	   virtual R minAbsNzo(bool unscaled = true) const;
	
	   /// Absolute biggest non-zero element in (in rational case possibly scaled) LP.
	   virtual R maxAbsNzo(bool unscaled = true) const;
	
	   /// Gets \p i 'th row.
	   void getRow(int i, LPRowBase<R>& row) const
	   {
	      row.setLhs(lhs(i));
	      row.setRhs(rhs(i));
	      row.setObj(rowObj(i));
	      row.setRowVector(DSVectorBase<R>(rowVector(i)));
	   }
	
	   /// Gets row with identifier \p id.
	   void getRow(const SPxRowId& id, LPRowBase<R>& row) const
	   {
	      getRow(number(id), row);
	   }
	
	   /// Gets rows \p start, ... \p end.
	   void getRows(int start, int end, LPRowSetBase<R>& set) const
	   {
	      set.clear();
	
	      for(int i = start; i <= end; i++)
	         set.add(lhs(i), rowVector(i), rhs(i), rowObj(i));
	   }
	
	   /// Gets row vector of row \p i.
	   const SVectorBase<R>& rowVector(int i) const
	   {
	      return LPRowSetBase<R>::rowVector(i);
	   }
	
	   /// Gets row vector of row with identifier \p id.
	   const SVectorBase<R>& rowVector(const SPxRowId& id) const
	   {
	      return LPRowSetBase<R>::rowVector(id);
	   }
	
	   /// Gets unscaled row vector of row \p i.
	   void getRowVectorUnscaled(int i, DSVectorBase<R>& vec) const;
	
	   /// Returns right hand side vector.
	   const VectorBase<R>& rhs() const
	   {
	      return LPRowSetBase<R>::rhs();
	   }
	
	   /// Returns right hand side of row number \p i.
	   const R& rhs(int i) const
	   {
	      return LPRowSetBase<R>::rhs(i);
	   }
	
	
	   /// Returns right hand side of row with identifier \p id.
	   const R& rhs(const SPxRowId& id) const
	   {
	      return LPRowSetBase<R>::rhs(id);
	   }
	
	   /// Gets (internal and possibly scaled) right hand side vector.
	   void getRhs(VectorBase<R>& vec) const
	   {
	      vec = LPRowSetBase<R>::rhs();
	   }
	
	   /// Gets unscaled right hand side vector.
	   void getRhsUnscaled(VectorBase<R>& vec) const;
	
	   /// Returns unscaled right hand side of row number \p i.
	   R rhsUnscaled(int i) const;
	
	   /// Returns unscaled right hand side of row with identifier \p id.
	   R rhsUnscaled(const SPxRowId& id) const;
	
	   /// Returns left hand side vector.
	   const VectorBase<R>& lhs() const
	   {
	      return LPRowSetBase<R>::lhs();
	   }
	
	   /// Returns left hand side of row number \p i.
	   const R& lhs(int i) const
	   {
	      return LPRowSetBase<R>::lhs(i);
	   }
	
	   /// Returns left hand side of row with identifier \p id.
	   const R& lhs(const SPxRowId& id) const
	   {
	      return LPRowSetBase<R>::lhs(id);
	   }
	
	   /// Gets row objective function vector.
	   void getRowObj(VectorBase<R>& prowobj) const
	   {
	      prowobj = LPRowSetBase<R>::obj();
	
	      if(spxSense() == MINIMIZE)
	         prowobj *= -1.0;
	   }
	
	   ///
	   R rowObj(int i) const
	   {
	      if(spxSense() == MINIMIZE)
	         return -maxRowObj(i);
	      else
	         return maxRowObj(i);
	   }
	
	   /// Returns row objective function value of row with identifier \p id.
	   R rowObj(const SPxRowId& id) const
	   {
	      if(spxSense() == MINIMIZE)
	         return -maxRowObj(id);
	      else
	         return maxRowObj(id);
	   }
	
	   ///
	   const VectorBase<R>& maxRowObj() const
	   {
	      return LPRowSetBase<R>::obj();
	   }
	
	   ///
	   const R& maxRowObj(int i) const
	   {
	      return LPRowSetBase<R>::obj(i);
	   }
	
	   /// Returns row objective function value of row with identifier \p id.
	   const R& maxRowObj(const SPxRowId& id) const
	   {
	      return LPRowSetBase<R>::obj(id);
	   }
	
	   /// Returns unscaled left hand side vector.
	   void getLhsUnscaled(VectorBase<R>& vec) const;
	
	   /// Returns unscaled left hand side of row number \p i.
	   R lhsUnscaled(int i) const;
	
	   /// Returns left hand side of row with identifier \p id.
	   R lhsUnscaled(const SPxRowId& id) const;
	
	   /// Returns the inequality type of the \p i'th LPRow.
	   typename LPRowBase<R>::Type rowType(int i) const
	   {
	      return LPRowSetBase<R>::type(i);
	   }
	
	   /// Returns the inequality type of the row with identifier \p key.
	   typename LPRowBase<R>::Type rowType(const SPxRowId& id) const
	   {
	      return LPRowSetBase<R>::type(id);
	   }
	
	   /// Gets \p i 'th column.
	   void getCol(int i, LPColBase<R>& col) const
	   {
	      col.setUpper(upper(i));
	      col.setLower(lower(i));
	      col.setObj(obj(i));
	      col.setColVector(colVector(i));
	   }
	
	   /// Gets column with identifier \p id.
	   void getCol(const SPxColId& id, LPColBase<R>& col) const
	   {
	      getCol(number(id), col);
	   }
	
	   /// Gets columns \p start, ..., \p end.
	   void getCols(int start, int end, LPColSetBase<R>& set) const
	   {
	      if(_isScaled)
	      {
	         LPColBase<R> lpcol;
	
	         for(int i = start; i <= end; i++)
	         {
	            getCol(i, lpcol);
	            set.add(lpcol);
	         }
	
	      }
	      else
	      {
	         set.clear();
	
	         for(int i = start; i <= end; i++)
	            set.add(obj(i), lower(i), colVector(i), upper(i));
	      }
	   }
	
	   /// Returns column vector of column \p i.
	   const SVectorBase<R>& colVector(int i) const
	   {
	      return LPColSetBase<R>::colVector(i);
	   }
	
	   /// Returns column vector of column with identifier \p id.
	   const SVectorBase<R>& colVector(const SPxColId& id) const
	   {
	      return LPColSetBase<R>::colVector(id);
	   }
	
	   /// Gets column vector of column \p i.
	   void getColVectorUnscaled(int i, DSVectorBase<R>& vec) const;
	
	   /// Gets column vector of column with identifier \p id.
	   void getColVectorUnscaled(const SPxColId& id, DSVectorBase<R>& vec) const;
	
	   /// Gets unscaled objective vector.
	   void getObjUnscaled(VectorBase<R>& pobj) const;
	
	   /// Gets objective vector.
	   void getObj(VectorBase<R>& pobj) const
	   {
	      pobj = LPColSetBase<R>::maxObj();
	
	      if(spxSense() == MINIMIZE)
	         pobj *= -1.0;
	   }
	
	   /// Returns objective value of column \p i.
	   R obj(int i) const
	   {
	      R res = maxObj(i);
	
	      if(spxSense() == MINIMIZE)
	         res *= -1;
	
	      return res;
	   }
	
	   /// Returns objective value of column with identifier \p id.
	   R obj(const SPxColId& id) const
	   {
	      return obj(number(id));
	   }
	
	   /// Returns unscaled objective value of column \p i.
	   R objUnscaled(int i) const;
	
	   /// Returns unscaled objective value of column with identifier \p id.
	   R objUnscaled(const SPxColId& id) const;
	
	   /// Returns objective vector for maximization problem.
	   /** Methods #maxObj() return the objective vector or its elements, after transformation to a maximization
	    *  problem. Since this is how SPxLPBase internally stores any LP these methods are generally faster. The following
	    *  condition holds: #obj() = #spxSense() * maxObj().
	    */
	   const VectorBase<R>& maxObj() const
	   {
	      return LPColSetBase<R>::maxObj();
	   }
	
	   /// Returns objective value of column \p i for maximization problem.
	   const R& maxObj(int i) const
	   {
	      return LPColSetBase<R>::maxObj(i);
	   }
	
	   /// Returns objective value of column with identifier \p id for maximization problem.
	   const R& maxObj(const SPxColId& id) const
	   {
	      return maxObj(number(id));
	   }
	
	   /// Returns unscaled objective vector for maximization problem.
	   void maxObjUnscaled(VectorBase<R>& vec) const;
	
	   /// Returns unscaled objective value of column \p i for maximization problem.
	   R maxObjUnscaled(int i) const;
	
	   /// Returns unscaled objective value of column with identifier \p id for maximization problem.
	   R maxObjUnscaled(const SPxColId& id) const;
	
	   /// Returns upper bound vector.
	   const VectorBase<R>& upper() const
	   {
	      return LPColSetBase<R>::upper();
	   }
	
	   /// Returns upper bound of column \p i.
	   const R& upper(int i) const
	   {
	      return LPColSetBase<R>::upper(i);
	   }
	
	   /// Returns upper bound of column with identifier \p id.
	   const R& upper(const SPxColId& id) const
	   {
	      return LPColSetBase<R>::upper(id);
	   }
	
	   /// Gets unscaled upper bound vector
	   void getUpperUnscaled(VectorBase<R>& vec) const;
	
	   /// Returns unscaled upper bound of column \p i.
	   R upperUnscaled(int i) const;
	
	   /// Returns unscaled upper bound of column with identifier \p id.
	   R upperUnscaled(const SPxColId& id) const;
	
	   /// Returns (internal and possibly scaled) lower bound vector.
	   const VectorBase<R>& lower() const
	   {
	      return LPColSetBase<R>::lower();
	   }
	
	   /// Returns (internal and possibly scaled) lower bound of column \p i.
	   const R& lower(int i) const
	   {
	      return LPColSetBase<R>::lower(i);
	   }
	
	   /// Returns (internal and possibly scaled) lower bound of column with identifier \p id.
	   const R& lower(const SPxColId& id) const
	   {
	      return LPColSetBase<R>::lower(id);
	   }
	
	   /// Gets unscaled lower bound vector.
	   void getLowerUnscaled(VectorBase<R>& vec) const;
	
	   /// Returns unscaled lower bound of column \p i.
	   R lowerUnscaled(int i) const;
	
	   /// Returns unscaled lower bound of column with identifier \p id.
	   R lowerUnscaled(const SPxColId& id) const;
	
	   /// Returns the optimization sense.
	   SPxSense spxSense() const
	   {
	      return thesense;
	   }
	
	   /// Returns the objective function value offset
	   const R& objOffset() const
	   {
	      return offset;
	   }
	
	   /// Returns the row number of the row with identifier \p id.
	   int number(const SPxRowId& id) const
	   {
	      return LPRowSetBase<R>::number(id);
	   }
	
	   /// Returns the column number of the column with identifier \p id.
	   int number(const SPxColId& id) const
	   {
	      return LPColSetBase<R>::number(id);
	   }
	
	   /// Returns the row or column number for identifier \p id.
	   int number(const SPxId& id) const
	   {
	      return (id.type() == SPxId::COL_ID)
	             ? LPColSetBase<R>::number(id)
	             : LPRowSetBase<R>::number(id);
	   }
	
	   /// Returns the row number of the row with identifier \p id.
	   bool has(const SPxRowId& id) const
	   {
	      return LPRowSetBase<R>::has(id);
	   }
	
	   /// Returns the column number of the column with identifier \p id.
	   bool has(const SPxColId& id) const
	   {
	      return LPColSetBase<R>::has(id);
	   }
	
	   /// Returns the row or column number for identifier \p id.
	   bool has(const SPxId& id) const
	   {
	      return (id.type() == SPxId::COL_ID)
	             ? LPColSetBase<R>::has(id)
	             : LPRowSetBase<R>::has(id);
	   }
	
	   /// Returns the row identifier for row \p n.
	   SPxRowId rId(int n) const
	   {
	      return SPxRowId(LPRowSetBase<R>::key(n));
	   }
	
	   /// Returns the column identifier for column \p n.
	   SPxColId cId(int n) const
	   {
	      return SPxColId(LPColSetBase<R>::key(n));
	   }
	
	   ///@}
	
	   // ------------------------------------------------------------------------------------------------------------------
	   /**@name Extension */
	   ///@{
	
	   ///
	   virtual void addRow(const LPRowBase<R>& row, bool scale = false)
	   {
	      doAddRow(row, scale);
	   }
	
	   ///
	   virtual void addRow(const R& lhsValue, const SVectorBase<R>& rowVec, const R& rhsValue,
	                       bool scale = false)
	   {
	      doAddRow(lhsValue, rowVec, rhsValue, scale);
	   }
	
	   ///
	   template < class S >
	   void addRow(const S* lhsValue, const S* rowValues, const int* rowIndices, int rowSize,
	               const S* rhsValue)
	   {
	      assert(lhsValue != 0);
	      assert(rowSize <= 0 || rowValues != 0);
	      assert(rowSize <= 0 || rowIndices != 0);
	      assert(rhsValue != 0);
	
	      int idx = nRows();
	      int oldColNumber = nCols();
	
	      LPRowSetBase<R>::add(lhsValue, rowValues, rowIndices, rowSize, rhsValue);
	
	      // now insert nonzeros to column file also
	      for(int j = rowSize - 1; j >= 0; --j)
	      {
	         const S& val = rowValues[j];
	         int i = rowIndices[j];
	
	         // create new columns if required
	         if(i >= nCols())
	         {
	            LPColBase<R> empty;
	
	            for(int k = nCols(); k <= i; ++k)
	               LPColSetBase<R>::add(empty);
	         }
	
	         assert(i < nCols());
	         LPColSetBase<R>::add2(i, 1, &idx, &val);
	      }
	
	      addedRows(1);
	      addedCols(nCols() - oldColNumber);
	   }
	
	   /// Adds \p row to LPRowSetBase.
	   virtual void addRow(SPxRowId& id, const LPRowBase<R>& row, bool scale = false)
	   {
	      addRow(row, scale);
	      id = rId(nRows() - 1);
	   }
	
	   ///
	   virtual void addRows(const LPRowSetBase<R>& pset, bool scale = false)
	   {
	      doAddRows(pset, scale);
	   }
	
	   ///
	   template < class S >
	   void addRows(const S* lhsValues, const S* rowValues, const int* rowIndices, const int* rowStarts,
	                const int* rowLengths, const int numRows, const int numValues, const S* rhsValues)
	   {
	      assert(lhsValues != 0);
	      assert(numValues <= 0 || rowValues != 0);
	      assert(numValues <= 0 || rowIndices != 0);
	      assert(numValues <= 0 || rowStarts != 0);
	      assert(numValues <= 0 || rowLengths != 0);
	      assert(rhsValues != 0);
	
	      int i, j, k, idx;
	      SVectorBase<R>* col;
	      DataArray < int > newCols(nCols());
	      int oldRowNumber = nRows();
	      int oldColNumber = nCols();
	
	      LPRowSetBase<R>::memRemax(oldRowNumber + numRows);
	
	      for(i = 0; i < numRows; i++)
	      {
	         assert(numValues <= 0 || rowStarts[i] + rowLengths[i] <= numValues);
	
	         if(numValues <= 0)
	            LPRowSetBase<R>::add(&(lhsValues[i]), (S*)0, (int*)0, 0, &(rhsValues[i]));
	         else
	            LPRowSetBase<R>::add(&(lhsValues[i]), &(rowValues[rowStarts[i]]), &(rowIndices[rowStarts[i]]),
	                                 rowLengths[i], &(rhsValues[i]));
	      }
	
	      assert(LPRowSetBase<R>::isConsistent());
	      assert(LPColSetBase<R>::isConsistent());
	
	      // count additional nonzeros per column
	      for(i = nCols() - 1; i >= 0; --i)
	         newCols[i] = 0;
	
	      if(numValues > 0)
	      {
	         for(i = 0; i < numRows; i++)
	         {
	            for(j = rowStarts[i]; j < rowStarts[i] + rowLengths[i]; j++)
	            {
	               ///@todo implement the addition of new columns as in doAddRows()
	               assert(rowIndices[j] >= 0);
	               assert(rowIndices[j] < oldColNumber);
	               newCols[rowIndices[j]]++;
	            }
	         }
	      }
	
	      // extend columns as required (backward because of memory efficiency reasons)
	      for(i = nCols() - 1; i >= 0; --i)
	      {
	         if(newCols[i] > 0)
	         {
	            int len = newCols[i] + colVector(i).size();
	            LPColSetBase<R>::xtend(i, len);
	
	            /* preset the sizes: beware that this can irritate a consistency check call from xtend(). We need to set the
	             * sizes here, because a possible garbage collection called from xtend might destroy the sizes again. */
	            colVector_w(i).set_size(len);
	         }
	      }
	
	      // insert new elements to column file
	      for(i = nRows() - 1; i >= oldRowNumber; --i)
	      {
	         const SVectorBase<R>& vec = rowVector(i);
	
	         for(j = vec.size() - 1; j >= 0; --j)
	         {
	            k = vec.index(j);
	            col = &colVector_w(k);
	            idx = col->size() - newCols[k];
	            assert(newCols[k] > 0);
	            assert(idx >= 0);
	            newCols[k]--;
	            col->index(idx) = i;
	            col->value(idx) = vec.value(j);
	         }
	      }
	
	#ifndef NDEBUG
	
	      for(i = 0; i < nCols(); ++i)
	         assert(newCols[i] == 0);
	
	#endif
	
	      assert(SPxLPBase<R>::isConsistent());
	
	      assert(numRows == nRows() - oldRowNumber);
	      addedRows(nRows() - oldRowNumber);
	      addedCols(nCols() - oldColNumber);
	   }
	
	   /// adds all LPRowBase%s of \p pset to LPRowSetBase.
	   virtual void addRows(SPxRowId id[], const LPRowSetBase<R>& set, bool scale = false)
	   {
	      int i = nRows();
	      addRows(set, scale);
	
	      for(int j = 0; i < nRows(); ++i, ++j)
	         id[j] = rId(i);
	   }
	
	   ///
	   virtual void addCol(const LPColBase<R>& col, bool scale = false)
	   {
	      doAddCol(col, scale);
	   }
	
	   ///
	   virtual void addCol(const R& objValue, const R& lowerValue, const SVectorBase<R>& colVec,
	                       const R& upperValue, bool scale = false)
	   {
	      doAddCol(objValue, lowerValue, colVec, upperValue, scale);
	   }
	
	   ///
	   template < class S >
	   void addCol(const S* objValue, const S* lowerValue, const S* colValues, const int* colIndices,
	               int colSize, const S* upperValue)
	   {
	      int idx = nCols();
	      int oldRowNumber = nRows();
	
	      LPColSetBase<R>::add(objValue, lowerValue, colValues, colIndices, colSize, upperValue);
	
	      if(thesense != MAXIMIZE)
	         LPColSetBase<R>::maxObj_w(idx) *= -1;
	
	      // now insert nonzeros to column file also
	      for(int j = colSize - 1; j >= 0; --j)
	      {
	         const S& val = colValues[j];
	         int i = colIndices[j];
	
	         // create new rows if required
	         if(i >= nRows())
	         {
	            LPRowBase<R> empty;
	
	            for(int k = nRows(); k <= i; ++k)
	               LPRowSetBase<R>::add(empty);
	         }
	
	         assert(i < nRows());
	         LPRowSetBase<R>::add2(i, 1, &idx, &val);
	      }
	
	      addedCols(1);
	      addedRows(nRows() - oldRowNumber);
	   }
	
	   /// Adds \p col to LPColSetVBase.
	   virtual void addCol(SPxColId& id, const LPColBase<R>& col, bool scale = false)
	   {
	      addCol(col, scale);
	      id = cId(nCols() - 1);
	   }
	
	   ///
	   virtual void addCols(const LPColSetBase<R>& pset, bool scale = false)
	   {
	      doAddCols(pset, scale);
	   }
	
	   ///
	   template < class S >
	   void addCols(const S* objValue, const S* lowerValues, const S* colValues, const int* colIndices,
	                const int* colStarts, const int* colLengths, const int numCols, const int numValues,
	                const S* upperValues)
	   {
	      assert(lowerValues != 0);
	      assert(numValues <= 0 || colValues != 0);
	      assert(numValues <= 0 || colIndices != 0);
	      assert(numValues <= 0 || colStarts != 0);
	      assert(numValues <= 0 || colLengths != 0);
	      assert(upperValues != 0);
	
	      int i, j, k, idx;
	      SVectorBase<R>* row;
	      DataArray < int > newRows(nRows());
	      int oldColNumber = nCols();
	      int oldRowNumber = nRows();
	      idx = nCols();
	
	      LPColSetBase<R>::memRemax(oldColNumber + numCols);
	
	      for(i = 0; i < numCols; i++)
	      {
	         assert(numValues <= 0 || colStarts[i] + colLengths[i] <= numValues);
	
	         if(numValues <= 0)
	            LPColSetBase<R>::add(&(objValue[i]), &(lowerValues[i]), (S*)0, (int*)0, 0, &(upperValues[i]));
	         else
	            LPColSetBase<R>::add(&(objValue[i]), &(lowerValues[i]), &(colValues[colStarts[i]]),
	                                 &(colIndices[colStarts[i]]), colLengths[i], &(upperValues[i]));
	
	         if(thesense != MAXIMIZE)
	            LPColSetBase<R>::maxObj_w(idx + i) *= -1;
	      }
	
	      assert(LPColSetBase<R>::isConsistent());
	      assert(LPRowSetBase<R>::isConsistent());
	
	      // count additional nonzeros per rows
	      for(i = nRows() - 1; i >= 0; --i)
	         newRows[i] = 0;
	
	      for(i = numValues - 1; i >= 0; --i)
	      {
	         ///@todo implement the addition of new rows as in doAddCols()
	         assert(colIndices[i] >= 0);
	         assert(colIndices[i] < oldRowNumber);
	         newRows[colIndices[i]]++;
	      }
	
	      // extend rows as required (backward because of memory efficiency reasons)
	      for(i = nRows() - 1; i >= 0; --i)
	      {
	         if(newRows[i] > 0)
	         {
	            int len = newRows[i] + rowVector(i).size();
	            LPRowSetBase<R>::xtend(i, len);
	
	            /* preset the sizes: beware that this can irritate a consistency check call from xtend(). We need to set the
	             * sizes here, because a possible garbage collection called from xtend might destroy the sizes again. */
	            rowVector_w(i).set_size(len);
	         }
	      }
	
	      // insert new elements to row file
	      for(i = nCols() - 1; i >= oldColNumber; --i)
	      {
	         const SVectorBase<R>& vec = colVector(i);
	
	         for(j = vec.size() - 1; j >= 0; --j)
	         {
	            k = vec.index(j);
	            row = &rowVector_w(k);
	            idx = row->size() - newRows[k];
	            assert(newRows[k] > 0);
	            assert(idx >= 0);
	            newRows[k]--;
	            row->index(idx) = i;
	            row->value(idx) = vec.value(j);
	         }
	      }
	
	#ifndef NDEBUG
	
	      for(i = 0; i < nRows(); ++i)
	         assert(newRows[i] == 0);
	
	#endif
	
	      assert(SPxLPBase<R>::isConsistent());
	
	      assert(numCols == nCols() - oldColNumber);
	      addedCols(nCols() - oldColNumber);
	      addedRows(nRows() - oldRowNumber);
	   }
	
	   /// Adds all LPColBase%s of \p set to LPColSetBase.
	   virtual void addCols(SPxColId id[], const LPColSetBase<R>& set, bool scale = false)
	   {
	
	      int i = nCols();
	      addCols(set, scale);
	
	      for(int j = 0; i < nCols(); ++i, ++j)
	         id[j] = cId(i);
	   }
	
	   ///@}
	
	   // ------------------------------------------------------------------------------------------------------------------
	   /**@name Shrinking */
	   ///@{
	
	   /// Removes \p i 'th row.
	   virtual void removeRow(int i)
	   {
	      if(i < 0)
	         return;
	
	      doRemoveRow(i);
	   }
	
	   /// Removes row with identifier \p id.
	   virtual void removeRow(SPxRowId id)
	   {
	      removeRow(number(id));
	   }
	
	   /// Removes multiple rows.
	   /** This method removes all LPRowBase%s from the SPxLPBase with an index \p i such that \p perm[i] < 0. Upon
	    *  completion, \p perm[i] >= 0 indicates the new index where the \p i'th LPRowBase<R> has been moved to due to this
	    *  removal. Note that \p perm must point to an array of at least #nRows() ints.
	    */
	   virtual void removeRows(int perm[])
	   {
	      doRemoveRows(perm);
	   }
	
	   ///
	   virtual void removeRows(SPxRowId id[], int n, int perm[] = 0)
	   {
	
	      if(perm == 0)
	      {
	         DataArray < int > p(nRows());
	         removeRows(id, n, p.get_ptr());
	         return;
	      }
	
	      for(int i = nRows() - 1; i >= 0; --i)
	         perm[i] = i;
	
	      while(n--)
	         perm[number(id[n])] = -1;
	
	      removeRows(perm);
	   }
	
	   /// Removes \p n LPRowBase%s.
	   /** Removing multiple rows with one method invocation is available in two flavours. An array \p perm can be passed as
	    *  third argument or not. If given, \p perm must be an array at least of size #nRows(). It is used to return the
	    *  permutations resulting from this removal: \p perm[i] < 0 indicates, that the element to index \p i has been
	    *  removed.  Otherwise, \p perm[i] is the new index of the element with index \p i before the removal.
	    */
	   virtual void removeRows(int nums[], int n, int perm[] = 0)
	   {
	
	      if(perm == 0)
	      {
	         DataArray < int > p(nRows());
	         removeRows(nums, n, p.get_ptr());
	         return;
	      }
	
	      for(int i = nRows() - 1; i >= 0; --i)
	         perm[i] = i;
	
	      while(n--)
	         perm[nums[n]] = -1;
	
	      removeRows(perm);
	   }
	
	   /// Removes rows from \p start to \p end (including both).
	   virtual void removeRowRange(int start, int end, int perm[] = 0)
	   {
	
	      if(perm == 0)
	      {
	         int i = end - start + 1;
	         DataArray < int > p(i);
	
	         while(--i >= 0)
	            p[i] = start + i;
	
	         removeRows(p.get_ptr(), end - start + 1);
	         return;
	      }
	
	      int i;
	
	      for(i = 0; i < start; ++i)
	         perm[i] = i;
	
	      for(; i <= end; ++i)
	         perm[i] = -1;
	
	      for(; i < nRows(); ++i)
	         perm[i] = i;
	
	      removeRows(perm);
	   }
	
	   /// Removes \p i 'th column.
	   virtual void removeCol(int i)
	   {
	      if(i < 0)
	         return;
	
	      doRemoveCol(i);
	   }
	
	   /// Removes column with identifier \p id.
	   virtual void removeCol(SPxColId id)
	   {
	      removeCol(number(id));
	   }
	
	   /// Removes multiple columns.
	   /** This method removes all LPColBase%s from the SPxLPBase with an index \p i such that \p perm[i] < 0. Upon
	    *  completion, \p perm[i] >= 0 indicates the new index where the \p i 'th LPColBase has been moved to due to this
	    *  removal. Note, that \p perm must point to an array of at least #nCols() ints.
	    */
	   virtual void removeCols(int perm[])
	   {
	      doRemoveCols(perm);
	   }
	
	   ///
	   virtual void removeCols(SPxColId id[], int n, int perm[] = 0)
	   {
	
	      if(perm == 0)
	      {
	         DataArray < int > p(nCols());
	         removeCols(id, n, p.get_ptr());
	         return;
	      }
	
	      for(int i = nCols() - 1; i >= 0; --i)
	         perm[i] = i;
	
	      while(n--)
	         perm[number(id[n])] = -1;
	
	      removeCols(perm);
	   }
	
	   /// Removes \p n LPCols.
	   /** Removing multiple columns with one method invocation is available in two flavours. An array \p perm can be passed
	    *  as third argument or not. If given, \p perm must be an array at least of size #nCols(). It is used to return the
	    *  permutations resulting from this removal: \p perm[i] < 0 indicates, that the element to index \p i has been
	    *  removed.  Otherwise, \p perm[i] is the new index of the element with index \p i before the removal.
	    */
	   virtual void removeCols(int nums[], int n, int perm[] = 0)
	   {
	
	      if(perm == 0)
	      {
	         DataArray < int > p(nCols());
	         removeCols(nums, n, p.get_ptr());
	         return;
	      }
	
	      for(int i = nCols() - 1; i >= 0; --i)
	         perm[i] = i;
	
	      while(n--)
	         perm[nums[n]] = -1;
	
	      removeCols(perm);
	   }
	
	   /// Removes columns from \p start to \p end (including both).
	   virtual void removeColRange(int start, int end, int perm[] = 0)
	   {
	
	      if(perm == 0)
	      {
	         int i = end - start + 1;
	         DataArray < int > p(i);
	
	         while(--i >= 0)
	            p[i] = start + i;
	
	         removeCols(p.get_ptr(), end - start + 1);
	         return;
	      }
	
	      int i;
	
	      for(i = 0; i < start; ++i)
	         perm[i] = i;
	
	      for(; i <= end; ++i)
	         perm[i] = -1;
	
	      for(; i < nCols(); ++i)
	         perm[i] = i;
	
	      removeCols(perm);
	   }
	
	   /// clears the LP.
	   virtual void clear()
	   {
	
	      LPRowSetBase<R>::clear();
	      LPColSetBase<R>::clear();
	      thesense = MAXIMIZE;
	      offset = 0;
	      _isScaled = false;
	      lp_scaler = nullptr;
	      LPColSetBase<R>::scaleExp.clear();
	      LPRowSetBase<R>::scaleExp.clear();
	   }
	
	   ///@}
	
	   // ------------------------------------------------------------------------------------------------------------------
	   /**@name IO */
	   ///@{
	
	   /// Reads LP in LP format from input stream \p in.
	   virtual bool readLPF(std::istream& in, NameSet* rowNames = 0, NameSet* colNames = 0,
	                        DIdxSet* intVars = 0);
	
	   /// Reads an LP in MPS format from input stream \p in.
	   virtual bool readMPS(std::istream& in, NameSet* rowNames = 0, NameSet* colNames = 0,
	                        DIdxSet* intVars = 0);
	
	   /// Reads LP in LP or MPS format from input stream \p in.
	   /**@param in       input stream.
	    * @param rowNames contains after the call the names of the constraints (rows) in the same order as the rows in the
	    *                 LP.  Constraints without a name (only possible with LPF files) are automatically assigned a name.
	    *                 Maybe 0 if the names are not needed.
	    * @param colNames contains after the call the names of the variables (columns) in the same order as the columns in
	    *                 the LP.  Maybe 0 if the names are not needed.
	    * @param intVars contains after the call the indices of those variables that where marked as beeing integer in the
	    *                 file.  Maybe 0 if the information is not needed.
	    * @todo Make sure the Id's in the NameSet%s are the same as in the LP.
	    */
	   virtual bool read(std::istream& in, NameSet* rowNames = 0, NameSet* colNames = 0,
	                     DIdxSet* intVars  = 0)
	   {
	      bool ok;
	      char c;
	
	      in.get(c);
	      in.putback(c);
	
	      /* MPS starts either with a comment mark '*' or with the keyword 'NAME' at the first column.  LPF starts either
	       * with blanks, a comment mark '\' or with the keyword "MAX" or "MIN" in upper or lower case.  There is no
	       * possible valid LPF file starting with a '*' or 'N'.
	       */
	      ok = ((c == '*') || (c == 'N'))
	           ? readMPS(in, rowNames, colNames, intVars)
	           : readLPF(in, rowNames, colNames, intVars);
	
	      return ok;
	   }
	
	   /// Reads LP from a file.
	   virtual bool readFile(const char* filename, NameSet* rowNames = 0, NameSet* colNames = 0,
	                         DIdxSet* intVars = 0)
	   {
	
	      spxifstream file(filename);
	
	      if(!file)
	         return false;
	
	      return read(file, rowNames, colNames, intVars);
	   }
	
	   /** Writes a file in LP format to \p out. If \p rowNames and \p colNames are \c NULL, default names are used for the
	    *  constraints and variables. If \p intVars is not \c NULL, the variables contained in it are marked as integer in
	    *  the output.
	    */
	   virtual void writeLPF(std::ostream&  out, const NameSet* rowNames, const NameSet* colNames,
	                         const DIdxSet* p_intvars = 0) const;
	
	   /// Writes a file in MPS format to \p out.
	   virtual void writeMPS(std::ostream&  out, const NameSet* rowNames, const NameSet* colNames,
	                         const DIdxSet* p_intvars = 0) const;
	
	   /// Write loaded LP to \p filename.
	   virtual void writeFileLPBase(const char* filename, const NameSet* rowNames = 0,
	                                const NameSet* colNames = 0, const DIdxSet* p_intvars = 0) const
	   {
	
	      std::ofstream tmp(filename);
	      size_t len_f = strlen(filename);
	
	      if(len_f > 4 && filename[len_f - 1] == 's' && filename[len_f - 2] == 'p'
	            && filename[len_f - 3] == 'm' && filename[len_f - 4] == '.')
	      {
	         writeMPS(tmp, rowNames, colNames, p_intvars);
	      }
	      else
	      {
	         writeLPF(tmp, rowNames, colNames, p_intvars);
	      }
	   }
	
	   /** prints problem statistics */
	   void printProblemStatistics(std::ostream& os)
	   {
	      int countLower = 0;
	      int countUpper = 0;
	      int countBoxed = 0;
	      int countFreeCol = 0;
	
	      int countEqual = 0;
	      int countLhs = 0;
	      int countRhs = 0;
	      int countRanged = 0;
	      int countFreeRow = 0;
	
	      for(int i = 0; i < nCols(); i++)
	      {
	         bool hasLower = false;
	         bool hasUpper = false;
	
	         if(lower(i) > R(-infinity))
	         {
	            countLower++;
	            hasLower = true;
	         }
	
	         if(upper(i) < R(infinity))
	         {
	            countUpper++;
	            hasUpper = true;
	         }
	
	         if(hasUpper && hasLower)
	         {
	            countBoxed++;
	            countLower--;
	            countUpper--;
	         }
	
	         if(!hasUpper && !hasLower)
	            countFreeCol++;
	      }
	
	      for(int i = 0; i < nRows(); i++)
	      {
	         bool hasRhs = false;
	         bool hasLhs = false;
	
	         if(lhs(i) > R(-infinity))
	         {
	            countLhs++;
	            hasLhs = true;
	         }
	
	         if(rhs(i) < R(infinity))
	         {
	            countRhs++;
	            hasRhs = true;
	         }
	
	         if(hasRhs && hasLhs)
	         {
	            if(EQ(lhs(i), rhs(i)))
	               countEqual++;
	            else
	               countRanged++;
	
	            countLhs--;
	            countRhs--;
	         }
	
	         if(!hasRhs && !hasLhs)
	            countFreeRow++;
	      }
	
	      SPxOut::setFixed(os);
	      os << "  Columns           : " << nCols() << "\n"
	         << "              boxed : " << countBoxed << "\n"
	         << "        lower bound : " << countLower << "\n"
	         << "        upper bound : " << countUpper << "\n"
	         << "               free : " << countFreeCol << "\n"
	         << "  Rows              : " << nRows() << "\n"
	         << "              equal : " << countEqual << "\n"
	         << "             ranged : " << countRanged << "\n"
	         << "                lhs : " << countLhs << "\n"
	         << "                rhs : " << countRhs << "\n"
	         << "               free : " << countFreeRow << "\n"
	         << "  Nonzeros          : " << nNzos() << "\n"
	         << "         per column : " << R(nNzos()) / R(nCols()) << "\n"
	         << "            per row : " << R(nNzos()) / R(nRows()) << "\n"
	         << "           sparsity : " << R(nNzos()) / R(nCols()) / R(nRows()) << "\n"
	         << "    min. abs. value : " << R(minAbsNzo()) << "\n"
	         << "    max. abs. value : " << R(maxAbsNzo()) << "\n";
	   }
	
	   ///@}
	
	   // ------------------------------------------------------------------------------------------------------------------
	   /**@name Manipulation */
	   ///@{
	
	   /// Changes objective vector to \p newObj. \p scale determines whether the new data should be scaled
	   virtual void changeObj(const VectorBase<R>& newObj, bool scale = false)
	   {
	      changeMaxObj(newObj, scale);
	
	      if(spxSense() == MINIMIZE)
	         LPColSetBase<R>::maxObj_w() *= -1;
	   }
	
	   /// changes \p i 'th objective vector element to \p newVal. \p scale determines whether the new data should be scaled
	   virtual void changeObj(int i, const R& newVal, bool scale = false)
	   {
	      changeMaxObj(i, newVal, scale);
	
	      if(spxSense() == MINIMIZE)
	         LPColSetBase<R>::maxObj_w(i) *= -1;
	   }
	
	   /// changes \p i 'th objective vector element to \p newVal.
	   template < class S >
	   void changeObj(int i, const S* newVal)
	   {
	      LPColSetBase<R>::maxObj_w(i) = *newVal;
	
	      if(spxSense() == MINIMIZE)
	         LPColSetBase<R>::maxObj_w(i) *= -1;
	
	      assert(isConsistent());
	   }
	
	   /// Changes objective value of column with identifier \p id to \p newVal. \p scale determines whether the new data should be scaled
	   virtual void changeObj(SPxColId id, const R& newVal, bool scale = false)
	   {
	      this->changeObj(number(id), newVal, scale);
	   }
	
	   /// Changes objective vector to \p newObj. \p scale determines whether the new data should be scaled
	   virtual void changeMaxObj(const VectorBase<R>& newObj, bool scale = false)
	   {
	      assert(maxObj().dim() == newObj.dim());
	
	      if(scale)
	      {
	         assert(_isScaled);
	         assert(lp_scaler);
	
	         for(int i = 0; i < maxObj().dim(); i++)

	            LPColSetBase<R>::maxObj_w(i) = lp_scaler->scaleObj(*this, i, newObj[i]);
	      }
	      else
	         LPColSetBase<R>::maxObj_w() = newObj;
	
	      assert(isConsistent());
	   }
	
	   /// changes \p i 'th objective vector element to \p newVal. \p scale determines whether the new data should be scaled
	   virtual void changeMaxObj(int i, const R& newVal, bool scale = false)
	   {
	      if(scale)
	      {
	         assert(_isScaled);
	         assert(lp_scaler);
	         LPColSetBase<R>::maxObj_w(i) = lp_scaler->scaleObj(*this, i, newVal);
	      }
	      else
	         LPColSetBase<R>::maxObj_w(i) = newVal;
	
	      assert(isConsistent());
	   }
	
	   /// changes \p i 'th objective vector element to \p newVal.
	   template < class S >
	   void changeMaxObj(int i, const S* newVal)
	   {
	      LPColSetBase<R>::maxObj_w(i) = *newVal;
	      assert(isConsistent());
	   }
	
	   /// Changes objective value of column with identifier \p id to \p newVal. \p scale determines whether the new data should be scaled
	   virtual void changeMaxObj(SPxColId id, const R& newVal, bool scale = false)
	   {
	      changeMaxObj(number(id), newVal, scale);
	   }
	
	   /// Changes vector of lower bounds to \p newLower. \p scale determines whether the new data should be scaled
	   virtual void changeLower(const VectorBase<R>& newLower, bool scale = false)
	   {
	      assert(lower().dim() == newLower.dim());
	
	      if(scale)
	      {
	         assert(_isScaled);
	         assert(lp_scaler);
	
	         for(int i = 0; i < lower().dim(); i++)

	            LPColSetBase<R>::lower_w(i) = lp_scaler->scaleLower(*this, i, newLower[i]);
	      }
	      else
	         LPColSetBase<R>::lower_w() = newLower;
	
	      assert(isConsistent());
	   }
	
	   /// changes \p i 'th lower bound to \p newLower. \p scale determines whether the new data should be scaled
	   virtual void changeLower(int i, const R& newLower, bool scale = false)
	   {
	      if(scale && newLower > R(-infinity))
	      {
	         assert(_isScaled);
	         assert(lp_scaler);
	         LPColSetBase<R>::lower_w(i) = lp_scaler->scaleLower(*this, i, newLower);
	      }
	      else
	         LPColSetBase<R>::lower_w(i) = newLower;
	
	      assert(isConsistent());
	   }
	
	   /// changes \p i 'th lower bound to \p newLower.
	   template < class S >
	   void changeLower(int i, const S* newLower)
	   {
	      LPColSetBase<R>::lower_w(i) = *newLower;
	      assert(isConsistent());
	   }
	
	   /// changes lower bound of column with identifier \p id to \p newLower. \p scale determines whether the new data should be scaled
	   virtual void changeLower(SPxColId id, const R& newLower, bool scale = false)
	   {
	      changeLower(number(id), newLower, scale);
	   }
	
	   /// Changes vector of upper bounds to \p newUpper. \p scale determines whether the new data should be scaled
	   virtual void changeUpper(const VectorBase<R>& newUpper, bool scale = false)
	   {
	      assert(upper().dim() == newUpper.dim());
	
	      if(scale)
	      {
	         assert(_isScaled);
	         assert(lp_scaler);
	
	         for(int i = 0; i < upper().dim(); i++)

	            LPColSetBase<R>::upper_w(i) = lp_scaler->scaleUpper(*this, i, newUpper[i]);
	      }
	      else
	         LPColSetBase<R>::upper_w() = newUpper;
	
	      assert(isConsistent());
	   }
	
	   /// Changes \p i 'th upper bound to \p newUpper. \p scale determines whether the new data should be scaled
	   virtual void changeUpper(int i, const R& newUpper, bool scale = false)
	   {
	      if(scale && newUpper < R(infinity))
	      {
	         assert(_isScaled);
	         assert(lp_scaler);
	         LPColSetBase<R>::upper_w(i) = lp_scaler->scaleUpper(*this, i, newUpper);
	      }
	      else
	         LPColSetBase<R>::upper_w(i) = newUpper;
	
	      assert(isConsistent());
	   }
	
	   /// Changes \p i 'th upper bound to \p newUpper.
	   template < class S >
	   void changeUpper(int i, const S* newUpper)
	   {
	      LPColSetBase<R>::upper_w(i) = *newUpper;
	      assert(isConsistent());
	   }
	
	   /// Changes upper bound of column with identifier \p id to \p newLower. \p scale determines whether the new data should be scaled
	   virtual void changeUpper(SPxColId id, const R& newUpper, bool scale = false)
	   {
	      changeUpper(number(id), newUpper, scale);
	   }
	
	   /// Changes variable bounds to \p newLower and \p newUpper. \p scale determines whether the new data should be scaled
	   virtual void changeBounds(const VectorBase<R>& newLower, const VectorBase<R>& newUpper,
	                             bool scale = false)
	   {
	      changeLower(newLower, scale);
	      changeUpper(newUpper, scale);
	      assert(isConsistent());
	   }
	
	   /// Changes bounds of column \p i to \p newLower and \p newUpper. \p scale determines whether the new data should be scaled
	   virtual void changeBounds(int i, const R& newLower, const R& newUpper, bool scale = false)
	   {
	      changeLower(i, newLower, scale);
	      changeUpper(i, newUpper, scale);
	      assert(isConsistent());
	   }
	
	   /// Changes bounds of column \p i to \p newLower and \p newUpper.
	   template < class S >
	   void changeBounds(int i, const S* newLower, const S* newUpper)
	   {
	      LPColSetBase<R>::lower_w(i) = *newLower;
	      LPColSetBase<R>::upper_w(i) = *newUpper;
	      assert(isConsistent());
	   }
	
	   /// Changes bounds of column with identifier \p id. \p scale determines whether the new data should be scaled
	   virtual void changeBounds(SPxColId id, const R& newLower, const R& newUpper, bool scale = false)
	   {
	      changeBounds(number(id), newLower, newUpper, scale);
	   }
	
	   /// Changes left hand side vector for constraints to \p newLhs. \p scale determines whether the new data should be scaled
	   virtual void changeLhs(const VectorBase<R>& newLhs, bool scale = false)
	   {
	      assert(lhs().dim() == newLhs.dim());
	
	      if(scale)
	      {
	         assert(_isScaled);
	         assert(lp_scaler);
	
	         for(int i = 0; i < lhs().dim(); i++)

	            LPRowSetBase<R>::lhs_w(i) = lp_scaler->scaleLhs(*this, i, newLhs[i]);
	      }
	      else
	         LPRowSetBase<R>::lhs_w() = newLhs;
	
	      assert(isConsistent());
	   }
	
	   /// Changes \p i 'th left hand side value to \p newLhs. \p scale determines whether the new data should be scaled
	   virtual void changeLhs(int i, const R& newLhs, bool scale = false)
	   {
	      if(scale && newLhs > R(-infinity))
	      {
	         assert(_isScaled);
	         assert(lp_scaler);
	         LPRowSetBase<R>::lhs_w(i) = lp_scaler->scaleLhs(*this, i, newLhs);
	      }
	      else
	         LPRowSetBase<R>::lhs_w(i) = newLhs;
	
	      assert(isConsistent());
	   }
	
	   /// Changes \p i 'th left hand side value to \p newLhs.
	   template < class S >
	   void changeLhs(int i, const S* newLhs)
	   {
	      LPRowSetBase<R>::lhs_w(i) = *newLhs;
	      assert(isConsistent());
	   }
	
	   /// Changes left hand side value for row with identifier \p id. \p scale determines whether the new data should be scaled
	   virtual void changeLhs(SPxRowId id, const R& newLhs, bool scale = false)
	   {
	      changeLhs(number(id), newLhs, scale);
	   }
	
	   /// Changes right hand side vector for constraints to \p newRhs. \p scale determines whether the new data should be scaled
	   virtual void changeRhs(const VectorBase<R>& newRhs, bool scale = false)
	   {
	      assert(rhs().dim() == newRhs.dim());
	
	      if(scale)
	      {
	         assert(_isScaled);
	         assert(lp_scaler);
	
	         for(int i = 0; i < rhs().dim(); i++)
	            LPRowSetBase<R>::rhs_w(i) = lp_scaler->scaleRhs(*this, i, newRhs[i]);
	      }
	      else
	         LPRowSetBase<R>::rhs_w() = newRhs;
	
	      assert(isConsistent());
	   }
	
	   /// Changes \p i 'th right hand side value to \p newRhs. \p scale determines whether the new data should be scaled
	   virtual void changeRhs(int i, const R& newRhs, bool scale = false)
	   {
	      if(scale && newRhs < R(infinity))
	      {
	         assert(_isScaled);
	         assert(lp_scaler);
	         LPRowSetBase<R>::rhs_w(i) = lp_scaler->scaleRhs(*this, i, newRhs);
	      }
	      else
	         LPRowSetBase<R>::rhs_w(i) = newRhs;
	
	      assert(isConsistent());
	   }
	
	   /// Changes right hand side value for row with identifier \p id. \p scale determines whether the new data should be scaled
	   virtual void changeRhs(SPxRowId id, const R& newRhs, bool scale = false)
	   {
	      changeRhs(number(id), newRhs, scale);
	   }
	
	   /// Changes left and right hand side vectors. \p scale determines whether the new data should be scaled
	   virtual void changeRange(const VectorBase<R>& newLhs, const VectorBase<R>& newRhs,
	                            bool scale = false)
	   {
	      changeLhs(newLhs, scale);
	      changeRhs(newRhs, scale);
	      assert(isConsistent());
	   }
	
	   /// Changes left and right hand side of row \p i. \p scale determines whether the new data should be scaled
	   virtual void changeRange(int i, const R& newLhs, const R& newRhs, bool scale = false)
	   {
	      changeLhs(i, newLhs, scale);
	      changeRhs(i, newRhs, scale);
	      assert(isConsistent());
	   }
	
	   /// Changes left and right hand side of row \p i.
	   template < class S >
	   void changeRange(int i, const S* newLhs, const S* newRhs)
	   {
	      LPRowSetBase<R>::lhs_w(i) = *newLhs;
	      LPRowSetBase<R>::rhs_w(i) = *newRhs;
	      assert(isConsistent());
	   }
	
	   /// Changes left and right hand side of row with identifier \p id. \p scale determines whether the new data should be scaled
	   virtual void changeRange(SPxRowId id, const R& newLhs, const R& newRhs, bool scale = false)
	   {
	      changeRange(number(id), newLhs, newRhs, scale);
	   }
	
	   /// Changes row objective function vector to \p newRowObj. \p scale determines whether the new data should be scaled
	   virtual void changeRowObj(const VectorBase<R>& newRowObj, bool scale = false)
	   {
	      assert(maxRowObj().dim() == newRowObj.dim());
	      LPRowSetBase<R>::obj_w() = newRowObj;
	
	      if(spxSense() == MINIMIZE)
	         LPRowSetBase<R>::obj_w() *= -1;
	
	      assert(isConsistent());
	   }
	
	   /// Changes \p i 'th row objective function value to \p newRowObj. \p scale determines whether the new data should be scaled
	   virtual void changeRowObj(int i, const R& newRowObj, bool scale = false)
	   {
	      LPRowSetBase<R>::obj_w(i) = newRowObj;
	
	      if(spxSense() == MINIMIZE)
	         LPRowSetBase<R>::obj_w(i) *= -1;
	
	      assert(isConsistent());
	   }
	
	   /// Changes row objective function value for row with identifier \p id. \p scale determines whether the new data should be scaled
	   virtual void changeRowObj(SPxRowId id, const R& newRowObj, bool scale = false)
	   {
	      changeRowObj(number(id), newRowObj, scale);
	   }
	
	   /// Clears row objective function values for all rows
	   virtual void clearRowObjs()
	   {
	      LPRowSetBase<R>::obj_w().clear();
	   }
	
	   /// Replaces \p i 'th row of LP with \p newRow. \p scale determines whether the new data should be scaled
	   virtual void changeRow(int n, const LPRowBase<R>& newRow, bool scale = false)
	   {
	      if(n < 0)
	         return;
	
	      int j;
	      SVectorBase<R>& row = rowVector_w(n);
	
	      for(j = row.size() - 1; j >= 0; --j)
	      {
	         SVectorBase<R>& col = colVector_w(row.index(j));
	         int position = col.pos(n);
	
	         assert(position != -1);
	
	         if(position >= 0)
	            col.remove(position);
	      }
	
	      row.clear();
	
	      changeLhs(n, newRow.lhs(), scale);
	      changeRhs(n, newRow.rhs(), scale);
	      changeRowObj(n, newRow.obj(), scale);
	
	      const SVectorBase<R>& newrow = newRow.rowVector();
	
	      for(j = newrow.size() - 1; j >= 0; --j)
	      {
	         int idx = newrow.index(j);
	         R val = newrow.value(j);
	
	         if(scale)
	            val = spxLdexp(val, LPRowSetBase<R>::scaleExp[n] + LPColSetBase<R>::scaleExp[idx]);
	
	         LPRowSetBase<R>::add2(n, 1, &idx, &val);
	         LPColSetBase<R>::add2(idx, 1, &n, &val);
	      }
	
	      assert(isConsistent());
	   }
	
	   /// Replaces row with identifier \p id with \p newRow. \p scale determines whether the new data should be scaled
	   virtual void changeRow(SPxRowId id, const LPRowBase<R>& newRow, bool scale = false)
	   {
	      changeRow(number(id), newRow, scale);
	   }
	
	   /// Replaces \p i 'th column of LP with \p newCol. \p scale determines whether the new data should be scaled
	   virtual void changeCol(int n, const LPColBase<R>& newCol, bool scale = false)
	   {
	      if(n < 0)
	         return;
	
	      int j;
	      SVectorBase<R>& col = colVector_w(n);
	
	      for(j = col.size() - 1; j >= 0; --j)
	      {
	         SVectorBase<R>& row = rowVector_w(col.index(j));
	         int position = row.pos(n);
	
	         assert(position != -1);
	
	         if(position >= 0)
	            row.remove(position);
	      }
	
	      col.clear();
	
	      changeUpper(n, newCol.upper(), scale);
	      changeLower(n, newCol.lower(), scale);
	      changeObj(n, newCol.obj(), scale);
	
	      const SVectorBase<R>& newcol = newCol.colVector();
	
	      for(j = newcol.size() - 1; j >= 0; --j)
	      {
	         int idx = newcol.index(j);
	         R val = newcol.value(j);
	
	         if(scale)
	            val = spxLdexp(val, LPColSetBase<R>::scaleExp[n] + LPRowSetBase<R>::scaleExp[idx]);
	
	         LPColSetBase<R>::add2(n, 1, &idx, &val);
	         LPRowSetBase<R>::add2(idx, 1, &n, &val);
	      }
	
	      assert(isConsistent());
	   }
	
	   /// Replaces column with identifier \p id with \p newCol. \p scale determines whether the new data should be scaled
	   virtual void changeCol(SPxColId id, const LPColBase<R>& newCol, bool scale = false)
	   {
	      changeCol(number(id), newCol, scale);
	   }
	
	   /// Changes LP element (\p i, \p j) to \p val. \p scale determines whether the new data should be scaled
	   virtual void changeElement(int i, int j, const R& val, bool scale = false)
	   {
	      if(i < 0 || j < 0)
	         return;
	
	      SVectorBase<R>& row = rowVector_w(i);
	      SVectorBase<R>& col = colVector_w(j);
	
	      if(isNotZero(val))
	      {
	         R newVal;
	
	         if(scale)
	         {
	            assert(_isScaled);
	            assert(lp_scaler);
	            newVal = lp_scaler->scaleElement(*this, i, j, val);
	         }
	         else
	            newVal = val;
	
	         if(row.pos(j) >= 0 && col.pos(i) >= 0)
	         {
	            row.value(row.pos(j)) = newVal;
	            col.value(col.pos(i)) = newVal;
	         }
	         else
	         {
	            LPRowSetBase<R>::add2(i, 1, &j, &newVal);
	            LPColSetBase<R>::add2(j, 1, &i, &newVal);
	         }
	      }
	      else if(row.pos(j) >= 0 && col.pos(i) >= 0)
	      {
	         row.remove(row.pos(j));
	         col.remove(col.pos(i));
	      }
	
	      assert(isConsistent());
	   }
	
	   /// Changes LP element (\p i, \p j) to \p val.
	   template < class S >
	   void changeElement(int i, int j, const S* val)
	   {
	      if(i < 0 || j < 0)
	         return;
	
	      SVectorBase<R>& row = rowVector_w(i);
	      SVectorBase<R>& col = colVector_w(j);
	
	      if(mpq_get_d(*val) != R(0))
	      {
	         if(row.pos(j) >= 0 && col.pos(i) >= 0)
	         {
	            row.value(row.pos(j)) = *val;
	            col.value(col.pos(i)) = *val;
	         }
	         else
	         {
	            LPRowSetBase<R>::add2(i, 1, &j, val);
	            LPColSetBase<R>::add2(j, 1, &i, val);
	         }
	      }
	      else if(row.pos(j) >= 0 && col.pos(i) >= 0)
	      {
	         row.remove(row.pos(j));
	         col.remove(col.pos(i));
	      }
	
	      assert(isConsistent());
	   }
	
	   /// Changes LP element identified by (\p rid, \p cid) to \p val. \p scale determines whether the new data should be scaled
	   virtual void changeElement(SPxRowId rid, SPxColId cid, const R& val, bool scale = false)
	   {
	      changeElement(number(rid), number(cid), val, scale);
	   }
	
	   /// Changes optimization sense to \p sns.
	   virtual void changeSense(SPxSense sns)
	   {
	      if(sns != thesense)
	      {
	         LPColSetBase<R>::maxObj_w() *= -1;
	         LPRowSetBase<R>::obj_w() *= -1;
	      }
	
	      thesense = sns;
	   }
	
	   template <typename T>
	   void changeObjOffset(const T& o)
	   {
	      offset = o;               // Converts o into type R. Example Rational into
	      // R
	   }
	
	   /// Computes activity of the rows for a given primal vector; activity does not need to be zero
	   /// @throw SPxInternalCodeException if the dimension of primal vector does not match number of columns or if the
	   ///        dimension of the activity vector does not match the number of rows
	   /// \p unscaled determines whether the returned data should be unscaled (if scaling was applied prior)
	   virtual void computePrimalActivity(const VectorBase<R>& primal, VectorBase<R>& activity,
	                                      const bool unscaled = true) const;
	
	   /// Updates activity of the rows for a given primal vector; activity does not need to be zero
	   /// @throw SPxInternalCodeException if the dimension of primal vector does not match number of columns or if the
	   ///        dimension of the activity vector does not match the number of rows
	   virtual void addPrimalActivity(const SVectorBase<R>& primal, VectorBase<R>& activity) const
	   {
	      if(activity.dim() != nRows())
	      {
	         throw SPxInternalCodeException("XSPXLP03 Activity vector computing row activity has wrong dimension");
	      }
	
	      for(int i = primal.size() - 1; i >= 0; i--)
	      {
	         assert(primal.index(i) >= 0);
	         assert(primal.index(i) < nCols());
	         activity.multAdd(primal.value(i), colVector(primal.index(i)));
	      }
	   }
	
	   /// Computes "dual" activity of the columns for a given dual vector, i.e., y^T A; activity does not need to be zero
	   /// @throw SPxInternalCodeException if dimension of dual vector does not match number of rows or if the dimension of
	   ///        the activity vector does not match the number of columns
	   virtual void computeDualActivity(const VectorBase<R>& dual, VectorBase<R>& activity,
	                                    const bool unscaled = true) const;
	
	   /// Updates "dual" activity of the columns for a given dual vector, i.e., y^T A; activity does not need to be zero
	   /// @throw SPxInternalCodeException if dimension of dual vector does not match number of rows or if the dimension of
	   ///        the activity vector does not match the number of columns
	   virtual void addDualActivity(const SVectorBase<R>& dual, VectorBase<R>& activity) const
	   {
	      if(activity.dim() != nCols())
	      {
	         throw SPxInternalCodeException("XSPXLP04 Activity vector computing dual activity has wrong dimension");
	      }
	
	      for(int i = dual.size() - 1; i >= 0; i--)
	      {
	         assert(dual.index(i) >= 0);
	         assert(dual.index(i) < nRows());
	         activity.multAdd(dual.value(i), rowVector(dual.index(i)));
	      }
	   }
	
	   /// Updates "dual" activity of the columns for a given dual vector, i.e., y^T A; activity does not need to be zero
	   /// @throw SPxInternalCodeException if dimension of dual vector does not match number of rows or if the dimension of
	   ///        the activity vector does not match the number of columns
	   virtual void subDualActivity(const VectorBase<R>& dual, VectorBase<R>& activity) const
	   {
	      if(dual.dim() != nRows())
	      {
	         throw SPxInternalCodeException("XSPXLP02 Dual vector for computing dual activity has wrong dimension");
	      }
	
	      if(activity.dim() != nCols())
	      {
	         throw SPxInternalCodeException("XSPXLP04 Activity vector computing dual activity has wrong dimension");
	      }
	
	      for(int r = 0; r < nRows(); r++)
	      {
	         if(dual[r] != 0)
	            activity.multSub(dual[r], rowVector(r));
	      }
	   }
	
	   ///@}
	
	   // ------------------------------------------------------------------------------------------------------------------
	   /**@name Construction of dual problem */
	   ///@{
	
	   /// Building the dual problem from a given LP
	   /// @note primalRows must be as large as the number of unranged primal rows + 2 * the number of ranged primal rows.
	   ///       dualCols must have the identical size to the primal rows.
	   virtual void buildDualProblem(SPxLPBase<R>& dualLP, SPxRowId primalRowIds[] = 0,
	                                 SPxColId primalColIds[] = 0,
	                                 SPxRowId dualRowIds[] = 0, SPxColId dualColIds[] = 0, int* nprimalrows = 0, int* nprimalcols = 0,
	                                 int* ndualrows = 0, int* ndualcols = 0);
	
	   ///@}
	
	   // ------------------------------------------------------------------------------------------------------------------
	   /**@name Miscellaneous */
	   ///@{
	
	   /// Consistency check.
	   bool isConsistent() const
	   {
	#ifdef ENABLE_CONSISTENCY_CHECKS
	
	      for(int i = nCols() - 1; i >= 0; --i)
	      {
	         const SVectorBase<R>& v = colVector(i);
	
	         for(int j = v.size() - 1; j >= 0; --j)
	         {
	            const SVectorBase<R>& w = rowVector(v.index(j));
	            int n = w.pos(i);
	
	            if(n < 0)
	               return MSGinconsistent("SPxLPBase");
	
	            if(v.value(j) != w.value(n))
	               return MSGinconsistent("SPxLPBase");
	         }
	      }
	
	      for(int i = nRows() - 1; i >= 0; --i)
	      {
	         const SVectorBase<R>& v = rowVector(i);
	
	         for(int j = v.size() - 1; j >= 0; --j)
	         {
	            const SVectorBase<R>& w = colVector(v.index(j));
	            int n = w.pos(i);
	
	            if(n < 0)
	               return MSGinconsistent("SPxLPBase");
	
	            if(v.value(j) != w.value(n))
	               return MSGinconsistent("SPxLPBase");
	         }
	      }
	
	      return LPRowSetBase<R>::isConsistent() && LPColSetBase<R>::isConsistent();
	#else
	      return true;
	#endif
	   }
	
	   ///@}
	
	protected:
	
	   // ------------------------------------------------------------------------------------------------------------------
	   /**@name Protected write access */
	   ///@{
	
	   /// Returns right hand side of row \p i.
	   R& rhs_w(int i)
	   {
	      return LPRowSetBase<R>::rhs_w(i);
	   }
	
	   /// Returns left hand side of row \p i.
	   R& lhs_w(int i)
	   {
	      return LPRowSetBase<R>::lhs_w(i);
	   }
	
	   /// Returns objective function value of row \p i.
	   R& maxRowObj_w(int i)
	   {
	      return LPRowSetBase<R>::obj_w(i);
	   }
	
	   /// Returns objective value of column \p i for maximization problem.
	   R& maxObj_w(int i)
	   {
	      return LPColSetBase<R>::maxObj_w(i);
	   }
	
	   /// Returns upper bound of column \p i.
	   R& upper_w(int i)
	   {
	      return LPColSetBase<R>::upper_w(i);
	   }
	
	   /// Returns lower bound of column \p i.
	   R& lower_w(int i)
	   {
	      return LPColSetBase<R>::lower_w(i);
	   }
	
	   ///@}
	
	   // ------------------------------------------------------------------------------------------------------------------
	   /**@name Protected helpers */
	   ///@{
	
	   /// Returns the LP as an LPRowSetBase.
	   const LPRowSetBase<R>* lprowset() const
	   {
	      return static_cast<const LPRowSetBase<R>*>(this);
	   }
	
	   /// Returns the LP as an LPColSetBase.
	   const LPColSetBase<R>* lpcolset() const
	   {
	      return static_cast<const LPColSetBase<R>*>(this);
	   }
	
	   /// Internal helper method.
	   virtual void doRemoveRow(int j)
	   {
	
	      const SVectorBase<R>& vec = rowVector(j);
	
	      // remove row vector from column file
	      for(int i = vec.size() - 1; i >= 0; --i)
	      {
	         SVectorBase<R>& remvec = colVector_w(vec.index(i));
	         int position = remvec.pos(j);
	
	         if(position >= 0)
	            remvec.remove(position);
	      }
	
	      // move last row to removed position
	      int idx = nRows() - 1;
	
	      if(j != idx)
	      {
	         const SVectorBase<R>& l_vec = rowVector(idx);
	
	         for(int i = l_vec.size() - 1; i >= 0; --i)
	         {
	            SVectorBase<R>& movevec = colVector_w(l_vec.index(i));
	            int position = movevec.pos(idx);
	
	            assert(position != -1);
	
	            if(position >= 0)
	               movevec.index(position) = j;
	         }
	      }
	
	      LPRowSetBase<R>::remove(j);
	   }
	
	   /// Internal helper method.
	   virtual void doRemoveRows(int perm[])
	   {
	      int j = nCols();
	
	      LPRowSetBase<R>::remove(perm);
	
	      for(int i = 0; i < j; ++i)
	      {
	         SVectorBase<R>& vec = colVector_w(i);
	
	         for(int k = vec.size() - 1; k >= 0; --k)
	         {
	            int idx = vec.index(k);
	
	            if(perm[idx] < 0)
	               vec.remove(k);
	            else
	               vec.index(k) = perm[idx];
	         }
	      }
	   }
	
	   /// Internal helper method.
	   virtual void doRemoveCol(int j)
	   {
	
	      const SVectorBase<R>& vec = colVector(j);
	      int i;
	
	      // remove column vector from row file
	      for(i = vec.size() - 1; i >= 0; --i)
	      {
	         SVectorBase<R>& remvec = rowVector_w(vec.index(i));
	         int position = remvec.pos(j);
	
	         assert(position != -1);
	
	         if(position >= 0)
	            remvec.remove(position);
	      }
	
	      // move last column to removed position
	      int idx = nCols() - 1;
	
	      if(j != idx)
	      {
	         const SVectorBase<R>& l_vec = colVector(idx);
	
	         for(i = l_vec.size() - 1; i >= 0; --i)
	         {
	            SVectorBase<R>& movevec = rowVector_w(l_vec.index(i));
	            int position = movevec.pos(idx);
	
	            assert(position != -1);
	
	            if(position >= 0)
	               movevec.index(position) = j;
	         }
	      }
	
	      LPColSetBase<R>::remove(j);
	   }
	
	   /// Internal helper method.
	   virtual void doRemoveCols(int perm[])
	   {
	      int nrows = nRows();
	
	      LPColSetBase<R>::remove(perm);
	
	      for(int i = 0; i < nrows; ++i)
	      {
	         SVectorBase<R>& vec = rowVector_w(i);
	
	         for(int k = vec.size() - 1; k >= 0; --k)
	         {
	            int idx = vec.index(k);
	
	            if(perm[idx] < 0)
	               vec.remove(k);
	            else
	               vec.index(k) = perm[idx];
	         }
	      }
	   }
	
	   /// Called after the last \p n rows have just been added.
	   virtual void addedRows(int newrows)
	   {}
	
	   /// Called after the last \p n columns have just been added.
	   virtual void addedCols(int newcols)
	   {}
	
	   ///
	   void added2Set(SVSetBase<R>& set, const SVSetBase<R>& addset, int n)
	   {
	
	      if(n == 0)
	         return;
	
	      DataArray<int> moreArray(set.num());
	      int* more = moreArray.get_ptr();
	
	      for(int i = set.num() - 1; i >= 0; --i)
	         more[i] = 0;
	
	      int tot = 0;
	      int end = addset.num();
	
	      for(int i = addset.num() - n; i < end; ++i)
	      {
	         const SVectorBase<R>& vec = addset[i];
	
	         tot += vec.size();
	
	         for(int j = vec.size() - 1; j >= 0; --j)
	            more[vec.index(j)]++;
	      }
	
	      if(set.memMax() < tot)
	         set.memRemax(tot);
	
	      for(int i = set.num() - 1; i >= 0; --i)
	      {
	         int j = set[i].size();
	         set.xtend(set[i], j + more[i]);
	         set[i].set_size(j + more[i]);
	         more[i] = j;
	      }
	
	      for(int i = addset.num() - n; i < addset.num(); ++i)
	      {
	         const SVectorBase<R>& vec = addset[i];
	
	         for(int j = vec.size() - 1; j >= 0; --j)
	         {
	            int k = vec.index(j);
	            int m = more[k]++;
	            SVectorBase<R>& l_xtend = set[k];
	            l_xtend.index(m) = i;
	            l_xtend.value(m) = vec.value(j);
	         }
	      }
	   }
	
	   ///@}
	
	
	private:
	
	   // ------------------------------------------------------------------------------------------------------------------
	   /**@name Private helpers */
	   ///@{
	
	   /// Returns the LP as an LPRowBase<R>Set.
	   SVectorBase<R>& colVector_w(int i)
	   {
	      return LPColSetBase<R>::colVector_w(i);
	   }
	
	   ///
	   SVectorBase<R>& rowVector_w(int i)
	   {
	      return LPRowSetBase<R>::rowVector_w(i);
	   }
	
	   ///
	   void doAddRow(const LPRowBase<R>& row, bool scale = false)
	   {
	      int idx = nRows();
	      int oldColNumber = nCols();
	      int newRowScaleExp = 0;
	
	      LPRowSetBase<R>::add(row);
	
	      SVectorBase<R>& vec = rowVector_w(idx);
	
	      DataArray <int>& colscaleExp = LPColSetBase<R>::scaleExp;
	
	      // compute new row scaling factor and apply it to the sides
	      if(scale && lp_scaler)
	      {
	         newRowScaleExp = lp_scaler->computeScaleExp(vec, colscaleExp);
	
	         if(rhs(idx) < R(infinity))
	            rhs_w(idx) = spxLdexp(rhs_w(idx), newRowScaleExp);
	
	         if(lhs(idx) > R(-infinity))
	            lhs_w(idx) = spxLdexp(lhs_w(idx), newRowScaleExp);
	
	         maxRowObj_w(idx) = spxLdexp(maxRowObj_w(idx), newRowScaleExp);
	
	         LPRowSetBase<R>::scaleExp[idx] = newRowScaleExp;
	      }
	
	      // now insert nonzeros to column file also
	      for(int j = vec.size() - 1; j >= 0; --j)
	      {
	         int i = vec.index(j);
	
	         // apply new row and existing column scaling factors to new values in RowSet
	         if(scale)
	            vec.value(j) = spxLdexp(vec.value(j), newRowScaleExp + colscaleExp[i]);
	
	         R val = vec.value(j);
	
	         // create new columns if required
	         if(i >= nCols())
	         {
	            LPColBase<R> empty;
	
	            for(int k = nCols(); k <= i; ++k)
	               LPColSetBase<R>::add(empty);
	         }
	
	         assert(i < nCols());
	         LPColSetBase<R>::add2(i, 1, &idx, &val);
	      }
	
	      addedRows(1);
	      addedCols(nCols() - oldColNumber);
	   }
	
	   ///
	   void doAddRow(const R& lhsValue, const SVectorBase<R>& rowVec, const R& rhsValue,
	                 bool scale = false)
	   {
	      int idx = nRows();
	      int oldColNumber = nCols();
	      int newRowScaleExp = 0;
	
	      LPRowSetBase<R>::add(lhsValue, rowVec, rhsValue);
	
	      DataArray <int>& colscaleExp = LPColSetBase<R>::scaleExp;
	
	      // compute new row scaling factor and apply it to the sides
	      if(scale)
	      {
	         newRowScaleExp = lp_scaler->computeScaleExp(rowVec, colscaleExp);
	
	         if(rhs(idx) < R(infinity))
	            rhs_w(idx) = spxLdexp(rhs_w(idx), newRowScaleExp);
	
	         if(lhs(idx) > R(-infinity))
	            lhs_w(idx) = spxLdexp(lhs_w(idx), newRowScaleExp);
	
	         maxRowObj_w(idx) = spxLdexp(maxRowObj_w(idx), newRowScaleExp);
	
	         LPRowSetBase<R>::scaleExp[idx] = newRowScaleExp;
	      }
	
	      SVectorBase<R>& vec = rowVector_w(idx);
	
	      // now insert nonzeros to column file also
	      for(int j = vec.size() - 1; j >= 0; --j)
	      {
	         int i = vec.index(j);
	
	         // apply new row and existing column scaling factors to new values in RowSet
	         if(scale)
	            vec.value(j) = spxLdexp(vec.value(j), newRowScaleExp + colscaleExp[i]);
	
	         R val = vec.value(j);
	
	         // create new columns if required
	         if(i >= nCols())
	         {
	            LPColBase<R> empty;
	
	            for(int k = nCols(); k <= i; ++k)
	               LPColSetBase<R>::add(empty);
	         }
	
	         assert(i < nCols());
	         LPColSetBase<R>::add2(i, 1, &idx, &val);
	      }
	
	      addedRows(1);
	      addedCols(nCols() - oldColNumber);
	   }
	
	   ///
	   void doAddRows(const LPRowSetBase<R>& set, bool scale = false)
	   {
	      int i, j, k, ii, idx;
	      SVectorBase<R>* col;
	      DataArray < int > newCols(nCols());
	      int oldRowNumber = nRows();
	      int oldColNumber = nCols();
	
	      if(&set != this)
	         LPRowSetBase<R>::add(set);
	
	      assert(LPRowSetBase<R>::isConsistent());
	      assert(LPColSetBase<R>::isConsistent());
	
	      // count additional nonzeros per column
	      for(i = nCols() - 1; i >= 0; --i)
	         newCols[i] = 0;
	
	      for(i = set.num() - 1; i >= 0; --i)
	      {
	         const SVectorBase<R>& vec = set.rowVector(i);
	
	         for(j = vec.size() - 1; j >= 0; --j)
	         {
	            // create new columns if required
	            ii = vec.index(j);
	
	            if(ii >= nCols())
	            {
	               LPColBase<R> empty;
	               newCols.reSize(ii + 1);
	
	               for(k = nCols(); k <= ii; ++k)
	               {
	                  newCols[k] = 0;
	                  LPColSetBase<R>::add(empty);
	               }
	            }
	
	            assert(ii < nCols());
	            newCols[ii]++;
	         }
	      }
	
	      // extend columns as required (backward because of memory efficiency reasons)
	      for(i = nCols() - 1; i >= 0; --i)
	      {
	         if(newCols[i] > 0)
	         {
	            int len = newCols[i] + colVector(i).size();
	            LPColSetBase<R>::xtend(i, len);
	
	            /* preset the sizes: beware that this can irritate a consistency check call from xtend(). We need to set the
	             * sizes here, because a possible garbage collection called from xtend might destroy the sizes again. */
	            colVector_w(i).set_size(len);
	         }
	      }
	
	      // compute new row scaling factor and insert new elements to column file
	      for(i = nRows() - 1; i >= oldRowNumber; --i)
	      {
	         SVectorBase<R>& vec = rowVector_w(i);
	         int newRowScaleExp = 0;
	
	         DataArray <int>& colscaleExp = LPColSetBase<R>::scaleExp;
	
	         // compute new row scaling factor and apply it to the sides
	         if(scale)
	         {
	            newRowScaleExp = lp_scaler->computeScaleExp(vec, colscaleExp);
	
	            if(rhs(i) < R(infinity))
	               rhs_w(i) = spxLdexp(rhs_w(i), newRowScaleExp);
	
	            if(lhs(i) > R(-infinity))
	               lhs_w(i) = spxLdexp(lhs_w(i), newRowScaleExp);
	
	            maxRowObj_w(i) = spxLdexp(maxRowObj_w(i), newRowScaleExp);
	
	            LPRowSetBase<R>::scaleExp[i] = newRowScaleExp;
	         }
	
	         for(j = vec.size() - 1; j >= 0; --j)
	         {
	            k = vec.index(j);
	            col = &colVector_w(k);
	            idx = col->size() - newCols[k];
	            assert(newCols[k] > 0);
	            assert(idx >= 0);
	            newCols[k]--;
	            col->index(idx) = i;
	
	            // apply new row and existing column scaling factors to both ColSet and RowSet
	            if(scale)
	               vec.value(j) = spxLdexp(vec.value(j), newRowScaleExp + colscaleExp[k]);
	
	            col->value(idx) = vec.value(j);
	         }
	      }
	
	#ifndef NDEBUG
	
	      for(i = 0; i < nCols(); ++i)
	         assert(newCols[i] == 0);
	
	#endif
	
	      assert(SPxLPBase<R>::isConsistent());
	
	      assert(set.num() == nRows() - oldRowNumber);
	      addedRows(nRows() - oldRowNumber);
	      addedCols(nCols() - oldColNumber);
	   }
	
	   ///
	   void doAddCol(const LPColBase<R>& col, bool scale = false)
	   {
	      int idx = nCols();
	      int oldRowNumber = nRows();
	      int newColScaleExp = 0;
	
	      LPColSetBase<R>::add(col);
	
	      if(thesense != MAXIMIZE)
	         LPColSetBase<R>::maxObj_w(idx) *= -1;
	
	      SVectorBase<R>& vec = colVector_w(idx);
	
	      DataArray <int>& rowscaleExp = LPRowSetBase<R>::scaleExp;
	
	      // compute new column scaling factor and apply it to the bounds
	      if(scale)
	      {
	         newColScaleExp = lp_scaler->computeScaleExp(vec, rowscaleExp);
	
	         if(upper(idx) < R(infinity))
	            upper_w(idx) = spxLdexp(upper_w(idx), - newColScaleExp);
	
	         if(lower(idx) > R(-infinity))
	            lower_w(idx) = spxLdexp(lower_w(idx), - newColScaleExp);
	
	         maxObj_w(idx) = spxLdexp(maxObj_w(idx), newColScaleExp);
	
	         LPColSetBase<R>::scaleExp[idx] = newColScaleExp;
	      }
	
	      // now insert nonzeros to row file also
	      for(int j = vec.size() - 1; j >= 0; --j)
	      {
	         int i = vec.index(j);
	
	         // apply new column and existing row scaling factors to new values in ColSet
	         if(scale)
	            vec.value(j) = spxLdexp(vec.value(j), newColScaleExp + rowscaleExp[i]);
	
	         R val = vec.value(j);
	
	         // create new rows if required
	         if(i >= nRows())
	         {
	            LPRowBase<R> empty;
	
	            for(int k = nRows(); k <= i; ++k)
	               LPRowSetBase<R>::add(empty);
	         }
	
	         assert(i < nRows());
	         LPRowSetBase<R>::add2(i, 1, &idx, &val);
	      }
	
	      addedCols(1);
	      addedRows(nRows() - oldRowNumber);
	   }
	
	   ///
	   void doAddCol(const R& objValue, const R& lowerValue, const SVectorBase<R>& colVec,
	                 const R& upperValue, bool scale = false)
	   {
	      int idx = nCols();
	      int oldRowNumber = nRows();
	      int newColScaleExp = 0;
	
	      LPColSetBase<R>::add(objValue, lowerValue, colVec, upperValue);
	
	      if(thesense != MAXIMIZE)
	         LPColSetBase<R>::maxObj_w(idx) *= -1;
	
	      DataArray <int>& rowscaleExp = LPRowSetBase<R>::scaleExp;
	
	      // compute new column scaling factor and apply it to the bounds
	      if(scale)
	      {
	         newColScaleExp = lp_scaler->computeScaleExp(colVec, rowscaleExp);
	
	         if(upper(idx) < R(infinity))
	            upper_w(idx) = spxLdexp(upper_w(idx), - newColScaleExp);
	
	         if(lower(idx) > R(-infinity))
	            lower_w(idx) = spxLdexp(lower_w(idx), - newColScaleExp);
	
	         maxObj_w(idx) = spxLdexp(maxObj_w(idx), newColScaleExp);
	
	         LPColSetBase<R>::scaleExp[idx] = newColScaleExp;
	      }
	
	      SVectorBase<R>& vec = colVector_w(idx);
	
	      // now insert nonzeros to row file also
	      for(int j = vec.size() - 1; j >= 0; --j)
	      {
	         int i = vec.index(j);
	
	         if(scale)
	            vec.value(j) = spxLdexp(vec.value(j), newColScaleExp + rowscaleExp[i]);
	
	         R val = vec.value(j);
	
	         // create new rows if required
	         if(i >= nRows())
	         {
	            LPRowBase<R> empty;
	
	            for(int k = nRows(); k <= i; ++k)
	               LPRowSetBase<R>::add(empty);
	         }
	
	         assert(i < nRows());
	         LPRowSetBase<R>::add2(i, 1, &idx, &val);
	      }
	
	      addedCols(1);
	      addedRows(nRows() - oldRowNumber);
	   }
	
	   ///
	   void doAddCols(const LPColSetBase<R>& set, bool scale = false)
	   {
	      int i, j;
	      int oldColNumber = nCols();
	      int oldRowNumber = nRows();
	      DataArray < int > newRows(nRows());
	
	      if(&set != this)
	         LPColSetBase<R>::add(set);
	
	      assert(LPColSetBase<R>::isConsistent());
	      assert(LPRowSetBase<R>::isConsistent());
	
	      // count additional nonzeros per row
	      for(i = nRows() - 1; i >= 0; --i)
	         newRows[i] = 0;
	
	      for(i = set.num() - 1; i >= 0; --i)
	      {
	         const SVectorBase<R>& vec = set.colVector(i);
	
	         for(j = vec.size() - 1; j >= 0; --j)
	         {
	            // create new rows if required
	            int l = vec.index(j);
	
	            if(l >= nRows())
	            {
	               LPRowBase<R> empty;
	               newRows.reSize(l + 1);
	
	               for(int k = nRows(); k <= l; ++k)
	               {
	                  newRows[k] = 0;
	                  LPRowSetBase<R>::add(empty);
	               }
	
	            }
	
	            assert(l < nRows());
	            newRows[l]++;
	         }
	      }
	
	      // extend rows as required
	      for(i = 0; i < nRows(); ++i)
	      {
	         if(newRows[i] > 0)
	         {
	            int len = newRows[i] + rowVector(i).size();
	            LPRowSetBase<R>::xtend(i, len);
	            rowVector_w(i).set_size(len);
	         }
	      }
	
	      // insert new elements to row file
	      for(i = oldColNumber; i < nCols(); ++i)
	      {
	         // @todo: Is there a better way to write the following if, else?
	         if(thesense == MAXIMIZE)
	         {
	            LPColSetBase<R>::maxObj_w(i) *= 1;
	         }
	         else                  // thesense is MINIMIZE = -1
	         {
	            LPColSetBase<R>::maxObj_w(i) *= -1;
	         }
	
	         SVectorBase<R>& vec = colVector_w(i);
	         int newColScaleExp = 0;
	
	         DataArray <int>& rowscaleExp = LPRowSetBase<R>::scaleExp;
	
	         // compute new column scaling factor and apply it to the bounds
	         if(scale)
	         {
	            newColScaleExp = lp_scaler->computeScaleExp(vec, rowscaleExp);
	
	            if(upper(i) < R(infinity))
	               upper_w(i) = spxLdexp(upper_w(i), - newColScaleExp);
	
	            if(lower(i) > R(-infinity))
	               lower_w(i) = spxLdexp(lower_w(i), - newColScaleExp);
	
	            maxObj_w(i) = spxLdexp(maxObj_w(i), newColScaleExp);
	
	            LPColSetBase<R>::scaleExp[i] = newColScaleExp;
	         }
	
	         for(j = vec.size() - 1; j >= 0; --j)
	         {
	            int k = vec.index(j);
	            SVectorBase<R>& row = rowVector_w(k);
	            int idx = row.size() - newRows[k];
	            assert(newRows[k] > 0);
	            newRows[k]--;
	            row.index(idx) = i;
	
	            // apply new column and existing row scaling factors to both ColSet and RowSet
	            if(scale)
	               vec.value(j) = spxLdexp(vec.value(j), newColScaleExp + rowscaleExp[k]);
	
	            row.value(idx) = vec.value(j);
	         }
	      }
	
	#ifndef NDEBUG
	
	      for(i = 0; i < nRows(); ++i)
	         assert(newRows[i] == 0);
	
	#endif
	
	      assert(SPxLPBase<R>::isConsistent());
	
	      assert(set.num() == nCols() - oldColNumber);
	      addedCols(nCols() - oldColNumber);
	      addedRows(nRows() - oldRowNumber);
	   }
	
	   ///@}
	
	public:
	
	   // ------------------------------------------------------------------------------------------------------------------
	   /**@name Constructors / Destructors */
	   ///@{
	
	   /// Default constructor.
	   SPxLPBase<R>()
	   {
	      SPxLPBase<R>::clear(); // clear is virtual.
	
	      assert(isConsistent());
	   }
	
	   /// Destructor.
	   virtual ~SPxLPBase<R>()
	   {}
	
	   /// Copy constructor.
	   SPxLPBase<R>(const SPxLPBase<R>& old)
	      : LPRowSetBase<R>(old)
	      , LPColSetBase<R>(old)
	      , thesense(old.thesense)
	      , offset(old.offset)
	      , _isScaled(old._isScaled)
	      , lp_scaler(old.lp_scaler)
	      , spxout(old.spxout)
	   {
	      assert(isConsistent());
	   }
	
	   /// Copy constructor.
	   template < class S >
	   SPxLPBase<R>(const SPxLPBase<S>& old)
	      : LPRowSetBase<R>(old)
	      , LPColSetBase<R>(old)
	      , thesense(old.thesense == SPxLPBase<S>::MINIMIZE ? SPxLPBase<R>::MINIMIZE : SPxLPBase<R>::MAXIMIZE)
	      , offset(old.offset)
	      , _isScaled(old._isScaled)
	      , spxout(old.spxout)
	   {
	      lp_scaler = nullptr;
	      assert(isConsistent());
	   }
	
	   /// Assignment operator.
	   SPxLPBase<R>& operator=(const SPxLPBase<R>& old)
	   {
	      if(this != &old)
	      {
	         LPRowSetBase<R>::operator=(old);
	         LPColSetBase<R>::operator=(old);
	         thesense = old.thesense;
	         offset = old.offset;
	         _isScaled = old._isScaled;
	         lp_scaler = old.lp_scaler;
	         spxout = old.spxout;
	
	         assert(isConsistent());
	      }
	
	      return *this;
	   }
	
	   /// Assignment operator.
	   template < class S >
	   SPxLPBase<R>& operator=(const SPxLPBase<S>& old)
	   {
	      if(this != (const SPxLPBase<R>*)(&old))
	      {
	         // The value of old.lp_scaler has to be nullptr
	         // Refer to issue #161 in soplex gitlab
	         assert(old.lp_scaler == nullptr);
	
	         LPRowSetBase<R>::operator=(old);
	         LPColSetBase<R>::operator=(old);
	         thesense = (old.thesense) == SPxLPBase<S>::MINIMIZE ? SPxLPBase<R>::MINIMIZE :
	                    SPxLPBase<R>::MAXIMIZE;
	         offset = R(old.offset);
	         _isScaled = old._isScaled;
	
	         // this may have un-intended consequences in the future
	         lp_scaler = nullptr;
	         spxout = old.spxout;
	
	         assert(isConsistent());
	      }
	
	      return *this;
	   }
	
	   ///@}
	};
	
	} // namespace soplex
	
	// For the general templated functions
	#include "spxlpbase_real.hpp"
	#include "spxlpbase_rational.hpp"
	
	/* reset the SOPLEX_DEBUG flag to its original value */
	#undef SOPLEX_DEBUG
	#ifdef SOPLEX_DEBUG_SPXLPBASE
	#define SOPLEX_DEBUG
	#undef SOPLEX_DEBUG_SPXLPBASE
	#endif
	
	#endif // _SPXLPBASE_H_