spxlpbase.h

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*                  This file is part of the class library                   */
/*       SoPlex --- the Sequential object-oriented simPlex.                  */
/*                                                                           */
/*  Copyright (c) 1996-2023 Zuse Institute Berlin (ZIB)                      */
/*                                                                           */
/*  Licensed under the Apache License, Version 2.0 (the "License");          */
/*  you may not use this file except in compliance with the License.         */
/*  You may obtain a copy of the License at                                  */
/*                                                                           */
/*      http://www.apache.org/licenses/LICENSE-2.0                           */
/*                                                                           */
/*  Unless required by applicable law or agreed to in writing, software      */
/*  distributed under the License is distributed on an "AS IS" BASIS,        */
/*  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
/*  See the License for the specific language governing permissions and      */
/*  limitations under the License.                                           */
/*                                                                           */
/*  You should have received a copy of the Apache-2.0 license                */
/*  along with SoPlex; see the file LICENSE. 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 current tolerances
   const std::shared_ptr<Tolerances> tolerances() const
   {
      return _tolerances;
   }

   /// set tolerances
   virtual void setTolerances(std::shared_ptr<Tolerances> tolerances)
   {
      this->_tolerances = tolerances;
   }

   /// 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), this->tolerances()->epsilon()))
               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, this->tolerances()->epsilon()))
      {
         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
   {
      if(this->_tolerances == nullptr && nCols() != 0)
         return SPX_MSG_INCONSISTENT("SPxLPBase");

#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 SPX_MSG_INCONSISTENT("SPxLPBase");

            if(v.value(j) != w.value(n))
               return SPX_MSG_INCONSISTENT("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 SPX_MSG_INCONSISTENT("SPxLPBase");

            if(v.value(j) != w.value(n))
               return SPX_MSG_INCONSISTENT("SPxLPBase");
         }
      }

      return LPRowSetBase<R>::isConsistent() && LPColSetBase<R>::isConsistent();
#else
      return true;
#endif
   }

   ///@}

protected:

   std::shared_ptr<Tolerances> _tolerances;

   // ------------------------------------------------------------------------------------------------------------------
   /**@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()
   {
      SPxLPBase<R>::clear(); // clear is virtual.

      assert(isConsistent());
   }

   /// Destructor.
   virtual ~SPxLPBase()
   {}

   /// Copy constructor.
   SPxLPBase(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)
   {
      _tolerances = old._tolerances;
      assert(isConsistent());
   }

   /// Copy constructor.
   template < class S >
   SPxLPBase(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;
      _tolerances = old._tolerances;
      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;
         _tolerances = old._tolerances;

         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;
         _tolerances = old._tolerances;

         // 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_