/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
	/*                                                                           */
	/*                  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.  */
	/*                                                                           */
	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
	
	#include <iostream>
	

	#include "soplex/spxmainsm.h"
	#include "soplex/array.h"
	#include "soplex/dataarray.h"
	#include "soplex/sorter.h"
	#include "soplex/spxout.h"
	#include <sstream>
	#include <iostream>
	#include <fstream>
	#include <memory>
	
	
	//rows
	#define FREE_LHS_RHS            1
	#define FREE_CONSTRAINT         1
	#define EMPTY_CONSTRAINT        1
	#define ROW_SINGLETON           1
	#define AGGREGATE_VARS          1
	#define FORCE_CONSTRAINT        1
	//cols
	#define FREE_BOUNDS             1
	#define EMPTY_COLUMN            1
	#define FIX_VARIABLE            1
	#define FREE_ZERO_OBJ_VARIABLE  1
	#define ZERO_OBJ_COL_SINGLETON  1
	#define DOUBLETON_EQUATION      1
	#define FREE_COL_SINGLETON      1
	//dual
	#define DOMINATED_COLUMN        1
	#define WEAKLY_DOMINATED_COLUMN 1
	#define MULTI_AGGREGATE         1
	//other
	#define TRIVIAL_HEURISTICS      1
	#define PSEUDOOBJ               1
	
	
	#define EXTREMES                1
	#define ROWS_SPXMAINSM                    1
	#define COLS_SPXMAINSM                    1
	#define DUAL_SPXMAINSM                    1
	///@todo check: with this simplification step, the unsimplified basis seems to be slightly suboptimal for some instances
	#define DUPLICATE_ROWS          1
	#define DUPLICATE_COLS          1
	
	
	#ifndef NDEBUG
	#define CHECK_BASIC_DIM
	#endif  // NDEBUG
	
	namespace soplex
	{
	
	template <class R>
	bool SPxMainSM<R>::PostStep::checkBasisDim(DataArray<typename SPxSolverBase<R>::VarStatus> rows,
	      DataArray<typename SPxSolverBase<R>::VarStatus> cols) const
	{
	   int numBasis = 0;
	
	   for(int rs = 0; rs < nRows; ++rs)
	   {
	      if(rows[rs] == SPxSolverBase<R>::BASIC)
	         numBasis++;
	   }
	
	   for(int cs = 0; cs < nCols; ++cs)
	   {
	      if(cols[cs] == SPxSolverBase<R>::BASIC)
	         numBasis++;
	   }
	
	   assert(numBasis == nRows);
	   return numBasis == nRows;
	}
	
	template <class R>
	void SPxMainSM<R>::RowObjPS::execute(VectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& s,
	                                     VectorBase<R>&,
	                                     DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
	                                     DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
	{
	   s[m_i] = s[m_i] - x[m_j];
	
	   assert(rStatus[m_i] != SPxSolverBase<R>::UNDEFINED);
	   assert(cStatus[m_j] != SPxSolverBase<R>::UNDEFINED);
	   assert(rStatus[m_i] != SPxSolverBase<R>::BASIC || cStatus[m_j] != SPxSolverBase<R>::BASIC);
	
	   MSG_DEBUG(std::cout << "RowObjPS: removing slack column " << m_j << " (" << cStatus[m_j] <<
	             ") for row " << m_i << " (" << rStatus[m_i] << ").\n");
	
	   if(rStatus[m_i] != SPxSolverBase<R>::BASIC)
	   {
	      switch(cStatus[m_j])
	      {
	      case SPxSolverBase<R>::ON_UPPER:
	         rStatus[m_i] = SPxSolverBase<R>::ON_LOWER;
	         break;
	
	      case SPxSolverBase<R>::ON_LOWER:
	         rStatus[m_i] = SPxSolverBase<R>::ON_UPPER;
	         break;
	
	      default:
	         rStatus[m_i] = cStatus[m_j];
	      }
	
	      // otherwise checkBasisDim() may fail
	      cStatus[m_j] = SPxSolverBase<R>::ZERO;
	   }
	
	#ifdef CHECK_BASIC_DIM
	
	   if(!this->checkBasisDim(rStatus, cStatus))
	   {
	      assert(false);
	      throw SPxInternalCodeException("XMAISM15 Dimension doesn't match after this step.");
	   }
	
	#endif
	}
	
	template <class R>
	void SPxMainSM<R>::FreeConstraintPS::execute(VectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& s,
	      VectorBase<R>&,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
	{
	   // correcting the change of idx by deletion of the row:
	   if(m_i != m_old_i)
	   {
	      s[m_old_i] = s[m_i];
	      y[m_old_i] = y[m_i];
	      rStatus[m_old_i] = rStatus[m_i];
	   }
	
	   // primal:
	   R slack = 0.0;
	
	   for(int k = 0; k < m_row.size(); ++k)
	      slack += m_row.value(k) * x[m_row.index(k)];
	
	   s[m_i] = slack;
	
	   // dual:
	   y[m_i] = m_row_obj;
	
	   // basis:
	   rStatus[m_i] = SPxSolverBase<R>::BASIC;
	
	#ifdef CHECK_BASIC_DIM
	
	   if(!this->checkBasisDim(rStatus, cStatus))
	   {
	      throw SPxInternalCodeException("XMAISM15 Dimension doesn't match after this step.");
	   }
	
	#endif
	}
	
	template <class R>
	void SPxMainSM<R>::EmptyConstraintPS::execute(VectorBase<R>&, VectorBase<R>& y, VectorBase<R>& s,
	      VectorBase<R>&,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
	{
	   // correcting the change of idx by deletion of the row:
	   if(m_i != m_old_i)
	   {
	      s[m_old_i] = s[m_i];
	      y[m_old_i] = y[m_i];
	      rStatus[m_old_i] = rStatus[m_i];
	   }
	
	   // primal:
	   s[m_i] = 0.0;
	
	   // dual:
	   y[m_i] = m_row_obj;
	
	   // basis:
	   rStatus[m_i] = SPxSolverBase<R>::BASIC;
	
	#ifdef CHECK_BASIC_DIM
	
	   if(!this->checkBasisDim(rStatus, cStatus))
	   {
	      throw SPxInternalCodeException("XMAISM16 Dimension doesn't match after this step.");
	   }
	
	#endif
	}
	
	template <class R>
	void SPxMainSM<R>::RowSingletonPS::execute(VectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& s,
	      VectorBase<R>& r,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
	{
	   // correcting the change of idx by deletion of the row:
	   if(m_i != m_old_i)
	   {
	      y[m_old_i] = y[m_i];
	      s[m_old_i] = s[m_i];
	      rStatus[m_old_i] = rStatus[m_i];
	   }
	
	   R aij = m_col[m_i];
	   assert(aij != 0.0);
	
	   // primal:
	   s[m_i] = aij * x[m_j];
	
	   // dual & basis:
	   R val = m_obj;
	
	   for(int k = 0; k < m_col.size(); ++k)
	   {
	      if(m_col.index(k) != m_i)
	         val -= m_col.value(k) * y[m_col.index(k)];
	   }
	
	   R newLo = (aij > 0) ? m_lhs / aij : m_rhs / aij; // implicit lhs
	   R newUp = (aij > 0) ? m_rhs / aij : m_lhs / aij; // implicit rhs
	
	   switch(cStatus[m_j])
	   {
	   case SPxSolverBase<R>::FIXED:
	      if(newLo <= m_oldLo && newUp >= m_oldUp)
	      {
	         // this row is totally redundant, has not changed bound of xj
	         rStatus[m_i] = SPxSolverBase<R>::BASIC;
	         y[m_i] = m_row_obj;
	      }
	      else if(EQrel(newLo, newUp, this->eps()))
	      {
	         // row is in the type  aij * xj = b
	         assert(EQrel(newLo, x[m_j], this->eps()));
	
	         if(EQrel(m_oldLo, m_oldUp, this->eps()))
	         {
	            // xj has been fixed in other row
	            rStatus[m_i] = SPxSolverBase<R>::BASIC;
	            y[m_i] = m_row_obj;
	         }
	         else if((EQrel(m_oldLo, x[m_j], this->eps()) && r[m_j] <= -this->eps())
	                 || (EQrel(m_oldUp, x[m_j], this->eps()) && r[m_j] >= this->eps())
	                 || (!EQrel(m_oldLo, x[m_j], this->eps()) && !(EQrel(m_oldUp, x[m_j], this->eps()))))
	         {
	            // if x_j on lower but reduced cost is negative, or x_j on upper but reduced cost is positive, or x_j not on bound: basic
	            rStatus[m_i] = (EQrel(m_lhs, x[m_j] * aij,
	                                  this->eps())) ? SPxSolverBase<R>::ON_LOWER : SPxSolverBase<R>::ON_UPPER;
	            cStatus[m_j] = SPxSolverBase<R>::BASIC;
	            y[m_i] = val / aij;
	            r[m_j] = 0.0;
	         }
	         else
	         {
	            // set x_j on one of the bound
	            assert(EQrel(m_oldLo, x[m_j], this->eps()) || EQrel(m_oldUp, x[m_j], this->eps()));
	
	            cStatus[m_j] = EQrel(m_oldLo, x[m_j],
	                                 this->eps()) ? SPxSolverBase<R>::ON_LOWER : SPxSolverBase<R>::ON_UPPER;
	            rStatus[m_i] = SPxSolverBase<R>::BASIC;
	            y[m_i] = m_row_obj;
	            r[m_j] = val;
	         }
	      }
	      else if(EQrel(newLo, m_oldUp, this->eps()))
	      {
	         // row is in the type  xj >= b/aij, try to set xj on upper
	         if(r[m_j] >= this->eps())
	         {
	            // the reduced cost is positive, xj should in the basic
	            assert(EQrel(m_rhs / aij, x[m_j], this->eps()) || EQrel(m_lhs / aij, x[m_j], this->eps()));
	
	            rStatus[m_i] = (EQrel(m_lhs / aij, x[m_j],
	                                  this->eps())) ? SPxSolverBase<R>::ON_LOWER : SPxSolverBase<R>::ON_UPPER;
	            cStatus[m_j] = SPxSolverBase<R>::BASIC;
	            y[m_i] = val / aij;
	            r[m_j] = 0.0;
	         }
	         else
	         {
	            assert(EQrel(m_oldUp, x[m_j], this->eps()));
	
	            cStatus[m_j] = SPxSolverBase<R>::ON_UPPER;
	            rStatus[m_i] = SPxSolverBase<R>::BASIC;
	            y[m_i] = m_row_obj;
	            r[m_j] = val;
	         }
	      }
	      else if(EQrel(newUp, m_oldLo, this->eps()))
	      {
	         // row is in the type  xj <= b/aij, try to set xj on lower
	         if(r[m_j] <= -this->eps())
	         {
	            // the reduced cost is negative, xj should in the basic
	            assert(EQrel(m_rhs / aij, x[m_j], this->eps()) || EQrel(m_lhs / aij, x[m_j], this->eps()));
	
	            rStatus[m_i] = (EQrel(m_lhs / aij, x[m_j],
	                                  this->eps())) ? SPxSolverBase<R>::ON_LOWER : SPxSolverBase<R>::ON_UPPER;
	            cStatus[m_j] = SPxSolverBase<R>::BASIC;
	            y[m_i] = val / aij;
	            r[m_j] = 0.0;
	         }
	         else
	         {
	            assert(EQrel(m_oldLo, x[m_j], this->eps()));
	
	            cStatus[m_j] = SPxSolverBase<R>::ON_LOWER;
	            rStatus[m_i] = SPxSolverBase<R>::BASIC;
	            y[m_i] = m_row_obj;
	            r[m_j] = val;
	         }
	      }
	      else
	      {
	         // the variable is set to FIXED by other constraints, i.e., this singleton row is redundant
	         rStatus[m_i] = SPxSolverBase<R>::BASIC;
	         y[m_i] = m_row_obj;
	      }
	
	      break;
	
	   case SPxSolverBase<R>::BASIC:
	      rStatus[m_i] = SPxSolverBase<R>::BASIC;
	      y[m_i] = m_row_obj;
	      r[m_j] = 0.0;
	      break;
	
	   case SPxSolverBase<R>::ON_LOWER:
	      if(EQrel(m_oldLo, x[m_j], this->eps())) // xj may stay on lower
	      {
	         rStatus[m_i] = SPxSolverBase<R>::BASIC;
	         y[m_i] = m_row_obj;
	         r[m_j] = val;
	      }
	      else // if reduced costs are negative or old lower bound not equal to xj, we need to change xj into the basis
	      {
	         assert(!isOptimal || EQrel(m_rhs / aij, x[m_j], this->eps())
	                || EQrel(m_lhs / aij, x[m_j], this->eps()));
	
	         cStatus[m_j] = SPxSolverBase<R>::BASIC;
	         rStatus[m_i] = (EQrel(m_lhs / aij, x[m_j],
	                               this->eps())) ? SPxSolverBase<R>::ON_LOWER : SPxSolverBase<R>::ON_UPPER;
	         y[m_i] = val / aij;
	         r[m_j] = 0.0;
	      }
	
	      break;
	
	   case SPxSolverBase<R>::ON_UPPER:
	      if(EQrel(m_oldUp, x[m_j], this->eps())) // xj may stay on upper
	      {
	         rStatus[m_i] = SPxSolverBase<R>::BASIC;
	         y[m_i] = m_row_obj;
	         r[m_j] = val;
	      }
	      else // if reduced costs are positive or old upper bound not equal to xj, we need to change xj into the basis
	      {
	         assert(!isOptimal || EQrel(m_rhs / aij, x[m_j], this->eps())
	                || EQrel(m_lhs / aij, x[m_j], this->eps()));
	
	         cStatus[m_j] = SPxSolverBase<R>::BASIC;
	         rStatus[m_i] = (EQrel(m_lhs / aij, x[m_j],
	                               this->eps())) ? SPxSolverBase<R>::ON_LOWER : SPxSolverBase<R>::ON_UPPER;
	         y[m_i] = val / aij;
	         r[m_j] = 0.0;
	      }
	
	      break;
	
	   case SPxSolverBase<R>::ZERO:
	      rStatus[m_i] = SPxSolverBase<R>::BASIC;
	      y[m_i] = m_row_obj;
	      r[m_j] = val;
	      break;
	
	   default:
	      break;
	   }
	
	#ifdef CHECK_BASIC_DIM
	
	   if(!this->checkBasisDim(rStatus, cStatus))
	   {
	      throw SPxInternalCodeException("XMAISM17 Dimension doesn't match after this step.");
	   }
	
	#endif
	}
	
	template <class R>
	void SPxMainSM<R>::ForceConstraintPS::execute(VectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& s,
	      VectorBase<R>& r,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
	{
	   // correcting the change of idx by deletion of the row:
	   if(m_i != m_old_i)
	   {
	      s[m_old_i] = s[m_i];
	      y[m_old_i] = y[m_i];
	      rStatus[m_old_i] = rStatus[m_i];
	   }
	
	   // primal:
	   s[m_i] = m_lRhs;
	
	   // basis:
	   int cBasisCandidate = -1;
	   R maxViolation = -1.0;
	   int bas_k = -1;
	
	   for(int k = 0; k < m_row.size(); ++k)
	   {
	      int  cIdx  = m_row.index(k);
	      R aij   = m_row.value(k);
	      R oldLo = m_oldLowers[k];
	      R oldUp = m_oldUppers[k];
	
	      switch(cStatus[cIdx])
	      {
	      case SPxSolverBase<R>::FIXED:
	         if(m_fixed[k])
	         {
	            assert(EQrel(oldLo, x[cIdx], this->eps()) || EQrel(oldUp, x[cIdx], this->eps()));
	
	            R violation = spxAbs(r[cIdx] / aij);
	
	            cStatus[cIdx] = EQrel(oldLo, x[cIdx],
	                                  this->eps()) ? SPxSolverBase<R>::ON_LOWER : SPxSolverBase<R>::ON_UPPER;
	
	            if(violation > maxViolation && ((EQrel(oldLo, x[cIdx], this->eps()) && r[cIdx] < -this->eps())
	                                            || (EQrel(oldUp, x[cIdx], this->eps()) && r[cIdx] > this->eps())))
	            {
	               maxViolation = violation;
	               cBasisCandidate = cIdx;
	               bas_k = k;
	            }
	         } // do nothing, if the old bounds are equal, i.e. variable has been not fixed in this row
	
	         break;
	
	      case SPxSolverBase<R>::ON_LOWER:
	      case SPxSolverBase<R>::ON_UPPER:
	      case SPxSolverBase<R>::BASIC:
	         break;
	
	      default:
	         break;
	      }
	   }
	
	   // dual and basis :
	   if(cBasisCandidate >= 0)  // one of the variable in the row should in the basis
	   {
	      assert(EQrel(m_lRhs, m_rhs, this->eps()) || EQrel(m_lRhs, m_lhs, this->eps()));
	      assert(bas_k >= 0);
	      assert(cBasisCandidate == m_row.index(bas_k));
	
	      cStatus[cBasisCandidate] = SPxSolverBase<R>::BASIC;
	      rStatus[m_i] = (EQrel(m_lRhs, m_lhs,
	                            this->eps())) ? SPxSolverBase<R>::ON_LOWER : SPxSolverBase<R>::ON_UPPER;
	
	      R aij = m_row.value(bas_k);
	      R multiplier = r[cBasisCandidate] / aij;
	      r[cBasisCandidate] = 0.0;
	
	      for(int k = 0; k < m_row.size(); ++k)  // update the reduced cost
	      {
	         if(k == bas_k)
	         {
	            continue;
	         }
	
	         r[m_row.index(k)] -= m_row.value(k) * multiplier;
	      }
	
	      // compute the value of new dual variable (because we have a new row)
	      R val = m_objs[bas_k];
	      DSVectorBase<R> basis_col = m_cols[bas_k];
	
	      for(int k = 0; k < basis_col.size(); ++k)
	      {
	         if(basis_col.index(k) != m_i)
	            val -= basis_col.value(k) * y[basis_col.index(k)];
	      }
	
	      y[m_i] = val / aij;
	   }
	   else // slack in the basis
	   {
	      rStatus[m_i] = SPxSolverBase<R>::BASIC;
	      y[m_i] = m_rowobj;
	   }
	
	#ifdef CHECK_BASIC_DIM
	
	   if(!this->checkBasisDim(rStatus, cStatus))
	   {
	      throw SPxInternalCodeException("XMAISM18 Dimension doesn't match after this step.");
	   }
	
	#endif
	}
	
	template <class R>
	void SPxMainSM<R>::FixVariablePS::execute(VectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& s,
	      VectorBase<R>& r,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
	{
	   // update the index mapping; if m_correctIdx is false, we assume that this has happened already
	   if(m_correctIdx)
	   {
	      x[m_old_j] = x[m_j];
	      r[m_old_j] = r[m_j];
	      cStatus[m_old_j] = cStatus[m_j];
	   }
	
	   // primal:
	   x[m_j] = m_val;
	
	   for(int k = 0; k < m_col.size(); ++k)
	      s[m_col.index(k)] += m_col.value(k) * x[m_j];
	
	   // dual:
	   R val = m_obj;
	
	   for(int k = 0; k < m_col.size(); ++k)
	      val -= m_col.value(k) * y[m_col.index(k)];
	
	   r[m_j] = val;
	
	   // basis:
	   if(m_lower == m_upper)
	   {
	      assert(EQrel(m_lower, m_val));
	
	      cStatus[m_j] = SPxSolverBase<R>::FIXED;
	   }
	   else
	   {
	      assert(EQrel(m_val, m_lower) || EQrel(m_val, m_upper) || m_val == 0.0);
	
	      cStatus[m_j] = EQrel(m_val, m_lower) ? SPxSolverBase<R>::ON_LOWER : (EQrel(m_val,
	                     m_upper) ? SPxSolverBase<R>::ON_UPPER : SPxSolverBase<R>::ZERO);
	   }
	
	#ifdef CHECK_BASIC_DIM
	
	   if(m_correctIdx)
	   {
	      if(!this->checkBasisDim(rStatus, cStatus))
	      {
	         throw SPxInternalCodeException("XMAISM19 Dimension doesn't match after this step.");
	      }
	   }
	
	#endif
	}
	
	template <class R>
	void SPxMainSM<R>::FixBoundsPS::execute(VectorBase<R>&, VectorBase<R>&, VectorBase<R>&,
	                                        VectorBase<R>&,
	                                        DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
	                                        DataArray<typename SPxSolverBase<R>::VarStatus>&, bool isOptimal) const
	{
	   // basis:
	   cStatus[m_j] = m_status;
	}
	
	template <class R>
	void SPxMainSM<R>::FreeZeroObjVariablePS::execute(VectorBase<R>& x, VectorBase<R>& y,
	      VectorBase<R>& s, VectorBase<R>& r,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
	{
	   // correcting the change of idx by deletion of the column and corresponding rows:
	   if(m_j != m_old_j)
	   {
	      x[m_old_j] = x[m_j];
	      r[m_old_j] = r[m_j];
	      cStatus[m_old_j] = cStatus[m_j];
	   }
	
	   int rIdx = m_old_i - m_col.size() + 1;
	
	   for(int k = 0; k < m_col.size(); ++k)
	   {
	      int rIdx_new = m_col.index(k);
	      s[rIdx] = s[rIdx_new];
	      y[rIdx] = y[rIdx_new];
	      rStatus[rIdx] = rStatus[rIdx_new];
	      rIdx++;
	   }
	
	   // primal:
	   int      domIdx = -1;
	   DSVectorBase<R> slack(m_col.size());
	
	   if(m_loFree)
	   {
	      R minRowUp = R(infinity);
	
	      for(int k = 0; k < m_rows.size(); ++k)
	      {
	         R           val = 0.0;
	         const SVectorBase<R>& row = m_rows[k];
	
	         for(int l = 0; l < row.size(); ++l)
	         {
	            if(row.index(l) != m_j)
	               val += row.value(l) * x[row.index(l)];
	         }
	
	         R scale = maxAbs(m_lRhs[k], val);
	
	         if(scale < 1.0)
	            scale = 1.0;
	
	         R z = (m_lRhs[k] / scale) - (val / scale);
	
	         if(isZero(z))
	            z = 0.0;
	
	         R up = z * scale / row[m_j];
	         slack.add(k, val);
	
	         if(up < minRowUp)
	         {
	            minRowUp = up;
	            domIdx   = k;
	         }
	      }
	
	      if(m_bnd < minRowUp)
	      {
	         x[m_j] = m_bnd;
	         domIdx = -1;
	      }
	      else
	         x[m_j] = minRowUp;
	   }
	   else
	   {
	      R maxRowLo = R(-infinity);
	
	      for(int k = 0; k < m_rows.size(); ++k)
	      {
	         R val = 0.0;
	         const SVectorBase<R>& row = m_rows[k];
	
	         for(int l = 0; l < row.size(); ++l)
	         {
	            if(row.index(l) != m_j)
	               val += row.value(l) * x[row.index(l)];
	         }
	
	         R scale = maxAbs(m_lRhs[k], val);
	
	         if(scale < 1.0)
	            scale = 1.0;
	
	         R z = (m_lRhs[k] / scale) - (val / scale);
	
	         if(isZero(z))
	            z = 0.0;
	
	         R lo = z * scale / row[m_j];
	         slack.add(k, val);
	
	         if(lo > maxRowLo)
	         {
	            maxRowLo = lo;
	            domIdx   = k;
	         }
	      }
	
	      if(m_bnd > maxRowLo)
	      {
	         x[m_j] = m_bnd;
	         domIdx = -1;
	      }
	      else
	         x[m_j] = maxRowLo;
	   }
	
	   for(int k = 0; k < m_col.size(); ++k)
	      s[m_col.index(k)] = slack[k] + m_col.value(k) * x[m_j];
	
	   // dual:
	   r[m_j] = 0.0;
	
	   for(int k = 0; k < m_col.size(); ++k)
	   {
	      int idx = m_col.index(k);
	      y[idx] = m_rowObj[idx];
	   }
	
	   // basis:
	   for(int k = 0; k < m_col.size(); ++k)
	   {
	      if(k != domIdx)
	         rStatus[m_col.index(k)] = SPxSolverBase<R>::BASIC;
	
	      else
	      {
	         cStatus[m_j] = SPxSolverBase<R>::BASIC;
	
	         if(m_loFree)
	            rStatus[m_col.index(k)] = (m_col.value(k) > 0) ? SPxSolverBase<R>::ON_UPPER :
	                                      SPxSolverBase<R>::ON_LOWER;
	         else
	            rStatus[m_col.index(k)] = (m_col.value(k) > 0) ? SPxSolverBase<R>::ON_LOWER :
	                                      SPxSolverBase<R>::ON_UPPER;
	      }
	   }
	
	   if(domIdx == -1)
	   {
	      if(m_loFree)
	         cStatus[m_j] = SPxSolverBase<R>::ON_UPPER;
	      else
	         cStatus[m_j] = SPxSolverBase<R>::ON_LOWER;
	   }
	
	#ifdef CHECK_BASIC_DIM
	
	   if(!this->checkBasisDim(rStatus, cStatus))
	   {
	      throw SPxInternalCodeException("XMAISM20 Dimension doesn't match after this step.");
	   }
	
	#endif
	}
	
	template <class R>
	void SPxMainSM<R>::ZeroObjColSingletonPS::execute(VectorBase<R>& x, VectorBase<R>& y,
	      VectorBase<R>& s, VectorBase<R>& r,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
	{
	   // correcting the change of idx by deletion of the column and corresponding rows:
	   if(m_j != m_old_j)
	   {
	      x[m_old_j] = x[m_j];
	      r[m_old_j] = r[m_j];
	      cStatus[m_old_j] = cStatus[m_j];
	   }
	
	   // primal & basis:
	   R aij = m_row[m_j];
	
	   if(isZero(s[m_i], R(1e-6)))
	      s[m_i] = 0.0;
	   else if(s[m_i] >= R(infinity))
	      // this is a fix for a highly ill conditioned instance that is "solved" in presolving (ilaser0 from MINLP, mittelmann)
	      throw SPxException("Simplifier: infinite activities - aborting unsimplification");
	
	   R scale1 = maxAbs(m_lhs, s[m_i]);
	   R scale2 = maxAbs(m_rhs, s[m_i]);
	
	   if(scale1 < 1.0)
	      scale1 = 1.0;
	
	   if(scale2 < 1.0)
	      scale2 = 1.0;
	
	   R z1 = (m_lhs / scale1) - (s[m_i] / scale1);
	   R z2 = (m_rhs / scale2) - (s[m_i] / scale2);
	
	   if(isZero(z1))
	      z1 = 0.0;
	
	   if(isZero(z2))
	      z2 = 0.0;
	
	   R lo = (aij > 0) ? z1 * scale1 / aij : z2 * scale2 / aij;
	   R up = (aij > 0) ? z2 * scale2 / aij : z1 * scale1 / aij;
	
	   if(isZero(lo, this->eps()))
	      lo = 0.0;
	
	   if(isZero(up, this->eps()))
	      up = 0.0;
	
	   assert(LErel(lo, up));
	   ASSERT_WARN("WMAISM01", isNotZero(aij, R(1.0 / R(infinity))));
	
	   if(rStatus[m_i] == SPxSolverBase<R>::ON_LOWER)
	   {
	      if(m_lower <= R(-infinity) && m_upper >= R(infinity))
	      {
	         x[m_j] = 0.0;
	         cStatus[m_j] = SPxSolverBase<R>::ZERO;
	      }
	      else if(m_lower == m_upper)
	      {
	         x[m_j]       = m_lower;
	         cStatus[m_j] = SPxSolverBase<R>::FIXED;
	      }
	      else if(aij > 0)
	      {
	         x[m_j]       = m_upper;
	         cStatus[m_j] = SPxSolverBase<R>::ON_UPPER;
	      }
	      else if(aij < 0)
	      {
	         x[m_j]       = m_lower;
	         cStatus[m_j] = SPxSolverBase<R>::ON_LOWER;
	      }
	      else
	         throw SPxInternalCodeException("XMAISM01 This should never happen.");
	   }
	   else if(rStatus[m_i] == SPxSolverBase<R>::ON_UPPER)
	   {
	      if(m_lower <= R(-infinity) && m_upper >= R(infinity))
	      {
	         x[m_j] = 0.0;
	         cStatus[m_j] = SPxSolverBase<R>::ZERO;
	      }
	      else if(m_lower == m_upper)
	      {
	         x[m_j]       = m_lower;
	         cStatus[m_j] = SPxSolverBase<R>::FIXED;
	      }
	      else if(aij > 0)
	      {
	         x[m_j]       = m_lower;
	         cStatus[m_j] = SPxSolverBase<R>::ON_LOWER;
	      }
	      else if(aij < 0)
	      {
	         x[m_j]       = m_upper;
	         cStatus[m_j] = SPxSolverBase<R>::ON_UPPER;
	      }
	      else
	         throw SPxInternalCodeException("XMAISM02 This should never happen.");
	   }
	   else if(rStatus[m_i] == SPxSolverBase<R>::FIXED)
	   {
	      if(m_lower <= R(-infinity) && m_upper >= R(infinity))
	      {
	         x[m_j] = 0.0;
	         cStatus[m_j] = SPxSolverBase<R>::ZERO;
	      }
	      else
	      {
	         assert(EQrel(m_lower, m_upper, this->eps()));
	
	         x[m_j]        = (m_lower + m_upper) / 2.0;
	         cStatus[m_j]  = SPxSolverBase<R>::FIXED;
	      }
	   }
	   else if(rStatus[m_i] == SPxSolverBase<R>::BASIC)
	   {
	      if(GErel(m_lower, lo, this->eps()) && m_lower > R(-infinity))
	      {
	         x[m_j]       = m_lower;
	         cStatus[m_j] = (m_lower == m_upper) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_LOWER;
	      }
	      else if(LErel(m_upper, up, this->eps()) && m_upper < R(infinity))
	      {
	         x[m_j]       = m_upper;
	         cStatus[m_j] = (m_lower == m_upper) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_UPPER;
	      }
	      else if(lo > R(-infinity))
	      {
	         // make m_i non-basic and m_j basic
	         x[m_j]       = lo;
	         cStatus[m_j] = SPxSolverBase<R>::BASIC;
	         rStatus[m_i] = (aij > 0 ? SPxSolverBase<R>::ON_LOWER : SPxSolverBase<R>::ON_UPPER);
	      }
	      else if(up < R(infinity))
	      {
	         // make m_i non-basic and m_j basic
	         x[m_j]       = up;
	         cStatus[m_j] = SPxSolverBase<R>::BASIC;
	         rStatus[m_i] = (aij > 0 ? SPxSolverBase<R>::ON_UPPER : SPxSolverBase<R>::ON_LOWER);
	      }
	      else
	         throw SPxInternalCodeException("XMAISM03 This should never happen.");
	   }
	   else
	      throw SPxInternalCodeException("XMAISM04 This should never happen.");
	
	   s[m_i] += aij * x[m_j];
	
	   // dual:
	   r[m_j] = -1.0 * aij * y[m_i];
	
	   assert(!isOptimal || (cStatus[m_j] != SPxSolverBase<R>::BASIC || isZero(r[m_j], this->eps())));
	
	#ifdef CHECK_BASIC_DIM
	
	   if(!this->checkBasisDim(rStatus, cStatus))
	   {
	      throw SPxInternalCodeException("XMAISM21 Dimension doesn't match after this step.");
	   }
	
	#endif
	}
	
	template <class R>
	void SPxMainSM<R>::FreeColSingletonPS::execute(VectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& s,
	      VectorBase<R>& r,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
	{
	
	   // correcting the change of idx by deletion of the row:
	   if(m_i != m_old_i)
	   {
	      s[m_old_i] = s[m_i];
	      y[m_old_i] = y[m_i];
	      rStatus[m_old_i] = rStatus[m_i];
	   }
	
	   // correcting the change of idx by deletion of the column:
	   if(m_j != m_old_j)
	   {
	      x[m_old_j] = x[m_j];
	      r[m_old_j] = r[m_j];
	      cStatus[m_old_j] = cStatus[m_j];
	   }
	
	   // primal:
	   R val = 0.0;
	   R aij = m_row[m_j];
	
	   for(int k = 0; k < m_row.size(); ++k)
	   {
	      if(m_row.index(k) != m_j)
	         val += m_row.value(k) * x[m_row.index(k)];
	   }
	
	   R scale = maxAbs(m_lRhs, val);
	
	   if(scale < 1.0)
	      scale = 1.0;
	
	   R z = (m_lRhs / scale) - (val / scale);
	
	   if(isZero(z))
	      z = 0.0;
	
	   x[m_j] = z * scale / aij;
	   s[m_i] = m_lRhs;
	
	   // dual:
	   y[m_i] = m_obj / aij;
	   r[m_j] = 0.0;
	
	   // basis:
	   cStatus[m_j] = SPxSolverBase<R>::BASIC;
	
	   if(m_eqCons)
	      rStatus[m_i] = SPxSolverBase<R>::FIXED;
	   else if(m_onLhs)
	      rStatus[m_i] = SPxSolverBase<R>::ON_LOWER;
	   else
	      rStatus[m_i] = SPxSolverBase<R>::ON_UPPER;
	
	#ifdef CHECK_BASIC_DIM
	
	   if(!this->checkBasisDim(rStatus, cStatus))
	   {
	      throw SPxInternalCodeException("XMAISM22 Dimension doesn't match after this step.");
	   }
	
	#endif
	}
	
	template <class R>
	void SPxMainSM<R>::DoubletonEquationPS::execute(VectorBase<R>& x, VectorBase<R>& y, VectorBase<R>&,
	      VectorBase<R>& r,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
	{
	   // dual:
	   if((cStatus[m_k]  != SPxSolverBase<R>::BASIC) &&
	         ((cStatus[m_k] == SPxSolverBase<R>::ON_LOWER && m_strictLo) ||
	          (cStatus[m_k] == SPxSolverBase<R>::ON_UPPER && m_strictUp) ||
	          (cStatus[m_k] == SPxSolverBase<R>::FIXED    &&
	           ((m_maxSense && ((r[m_j] > 0 && m_strictUp) || (r[m_j] < 0 && m_strictLo))) ||
	            (!m_maxSense && ((r[m_j] > 0 && m_strictLo) || (r[m_j] < 0 && m_strictUp)))))))
	   {
	      R val  = m_kObj;
	      R aik  = m_col[m_i];
	
	      for(int _k = 0; _k < m_col.size(); ++_k)
	      {
	         if(m_col.index(_k) != m_i)
	            val -= m_col.value(_k) * y[m_col.index(_k)];
	      }
	
	      y[m_i] = val / aik;
	      r[m_k] = 0.0;
	
	      r[m_j] = m_jObj - val * m_aij / aik;
	
	      ASSERT_WARN("WMAISM73", isNotZero(m_aij * aik));
	
	      // basis:
	      if(m_jFixed)
	         cStatus[m_j] = SPxSolverBase<R>::FIXED;
	      else
	      {
	         if(GT(r[m_j], (R) 0) || (isZero(r[m_j]) && EQ(x[m_j], m_Lo_j)))
	            cStatus[m_j] = SPxSolverBase<R>::ON_LOWER;
	         else
	            cStatus[m_j] = SPxSolverBase<R>::ON_UPPER;
	      }
	
	      cStatus[m_k] = SPxSolverBase<R>::BASIC;
	   }
	
	#ifdef CHECK_BASIC_DIM
	
	   if(!this->checkBasisDim(rStatus, cStatus))
	   {
	      throw SPxInternalCodeException("XMAISM23 Dimension doesn't match after this step.");
	   }
	
	#endif
	}
	
	template <class R>
	void SPxMainSM<R>::DuplicateRowsPS::execute(VectorBase<R>&, VectorBase<R>& y, VectorBase<R>& s,
	      VectorBase<R>&,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
	{
	   // correcting the change of idx by deletion of the duplicated rows:
	   if(m_isLast)
	   {
	      for(int i = m_perm.size() - 1; i >= 0; --i)
	      {
	         if(m_perm[i] >= 0)
	         {
	            int rIdx_new = m_perm[i];
	            int rIdx = i;
	            s[rIdx] = s[rIdx_new];
	            y[rIdx] = y[rIdx_new];
	            rStatus[rIdx] = rStatus[rIdx_new];
	         }
	      }
	   }
	
	   // primal:
	   for(int k = 0; k < m_scale.size(); ++k)
	   {
	      if(m_scale.index(k) != m_i)
	         s[m_scale.index(k)] = s[m_i] / m_scale.value(k);
	   }
	
	   // dual & basis:
	   bool haveSetBasis = false;
	
	   for(int k = 0; k < m_scale.size(); ++k)
	   {
	      int i = m_scale.index(k);
	
	      if(rStatus[m_i] == SPxSolverBase<R>::BASIC || (haveSetBasis && i != m_i))
	         // if the row with tightest lower and upper bound in the basic, every duplicate row should in basic
	         // or basis status of row m_i has been set, this row should be in basis
	      {
	         y[i]       = m_rowObj.value(k);
	         rStatus[i] = SPxSolverBase<R>::BASIC;
	         continue;
	      }
	
	      ASSERT_WARN("WMAISM02", isNotZero(m_scale.value(k)));
	
	      if(rStatus[m_i] == SPxSolverBase<R>::FIXED && (i == m_maxLhsIdx || i == m_minRhsIdx))
	      {
	         // this row leads to the tightest lower or upper bound, slack should not be in the basis
	         y[i]   = y[m_i] * m_scale.value(k);
	         y[m_i] = m_i_rowObj;
	
	         if(m_isLhsEqualRhs[k])
	         {
	            rStatus[i] = SPxSolverBase<R>::FIXED;
	         }
	         else if(i == m_maxLhsIdx)
	         {
	            rStatus[i] = m_scale.value(k) * m_scale.value(0) > 0 ? SPxSolverBase<R>::ON_LOWER :
	                         SPxSolverBase<R>::ON_UPPER;
	         }
	         else
	         {
	            assert(i == m_minRhsIdx);
	
	            rStatus[i] = m_scale.value(k) * m_scale.value(0) > 0 ? SPxSolverBase<R>::ON_UPPER :
	                         SPxSolverBase<R>::ON_LOWER;
	         }
	
	         if(i != m_i)
	            rStatus[m_i] = SPxSolverBase<R>::BASIC;
	
	         haveSetBasis = true;
	      }
	      else if(i == m_maxLhsIdx && rStatus[m_i] == SPxSolverBase<R>::ON_LOWER)
	      {
	         // this row leads to the tightest lower bound, slack should not be in the basis
	         y[i]   = y[m_i] * m_scale.value(k);
	         y[m_i] = m_i_rowObj;
	
	         rStatus[i] = m_scale.value(k) * m_scale.value(0) > 0 ? SPxSolverBase<R>::ON_LOWER :
	                      SPxSolverBase<R>::ON_UPPER;
	
	         if(i != m_i)
	            rStatus[m_i] = SPxSolverBase<R>::BASIC;
	
	         haveSetBasis = true;
	      }
	      else if(i == m_minRhsIdx && rStatus[m_i] == SPxSolverBase<R>::ON_UPPER)
	      {
	         // this row leads to the tightest upper bound, slack should not be in the basis
	         y[i]   = y[m_i] * m_scale.value(k);
	         y[m_i] = m_i_rowObj;
	
	         rStatus[i] = m_scale.value(k) * m_scale.value(0) > 0 ? SPxSolverBase<R>::ON_UPPER :
	                      SPxSolverBase<R>::ON_LOWER;
	
	         if(i != m_i)
	            rStatus[m_i] = SPxSolverBase<R>::BASIC;
	
	         haveSetBasis = true;
	      }
	      else if(i != m_i)
	      {
	         // this row does not lead to the tightest lower or upper bound, slack should be in the basis
	         y[i]       = m_rowObj.value(k);
	         rStatus[i] = SPxSolverBase<R>::BASIC;
	      }
	   }
	
	#ifdef CHECK_BASIC_DIM
	
	   if(m_isFirst && !this->checkBasisDim(rStatus, cStatus))
	   {
	      throw SPxInternalCodeException("XMAISM24 Dimension doesn't match after this step.");
	   }
	
	#endif
	
	   // nothing to do for the reduced cost values
	}
	
	template <class R>
	void SPxMainSM<R>::DuplicateColsPS::execute(VectorBase<R>& x,
	      VectorBase<R>&,
	      VectorBase<R>&,
	      VectorBase<R>& r,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
	{
	
	   if(m_isFirst)
	   {
	#ifdef CHECK_BASIC_DIM
	
	      if(!this->checkBasisDim(rStatus, cStatus))
	      {
	         throw SPxInternalCodeException("XMAISM25 Dimension doesn't match after this step.");
	      }
	
	#endif
	      return;
	   }
	
	
	   // correcting the change of idx by deletion of the columns:
	   if(m_isLast)
	   {
	      for(int i = m_perm.size() - 1; i >= 0; --i)
	      {
	         if(m_perm[i] >= 0)
	         {
	            int cIdx_new = m_perm[i];
	            int cIdx = i;
	            x[cIdx] = x[cIdx_new];
	            r[cIdx] = r[cIdx_new];
	            cStatus[cIdx] = cStatus[cIdx_new];
	         }
	      }
	
	      return;
	   }
	
	   // primal & basis:
	   ASSERT_WARN("WMAISM03", isNotZero(m_scale));
	
	   if(cStatus[m_k] == SPxSolverBase<R>::ON_LOWER)
	   {
	      x[m_k] = m_loK;
	
	      if(m_scale > 0)
	      {
	         x[m_j]       = m_loJ;
	         cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_LOWER;
	      }
	      else
	      {
	         x[m_j]       = m_upJ;
	         cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_UPPER;
	      }
	   }
	   else if(cStatus[m_k] == SPxSolverBase<R>::ON_UPPER)
	   {
	      x[m_k] = m_upK;
	
	      if(m_scale > 0)
	      {
	         x[m_j]       = m_upJ;
	         cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_UPPER;
	      }
	      else
	      {
	         x[m_j]       = m_loJ;
	         cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_LOWER;
	      }
	   }
	   else if(cStatus[m_k] == SPxSolverBase<R>::FIXED)
	   {
	      // => x[m_k] and x[m_j] are also fixed before the corresponding preprocessing step
	      x[m_j]       = m_loJ;
	      cStatus[m_j] = SPxSolverBase<R>::FIXED;
	   }
	   else if(cStatus[m_k] == SPxSolverBase<R>::ZERO)
	   {
	      /* we only aggregate duplicate columns if 0 is contained in their bounds, so we can handle this case properly */
	      assert(isZero(x[m_k]));
	      assert(LErel(m_loJ, R(0.0)));
	      assert(GErel(m_upJ, R(0.0)));
	      assert(LErel(m_loK, R(0.0)));
	      assert(GErel(m_upK, R(0.0)));
	
	      if(isZero(m_loK) && isZero(m_upK) && m_loK == m_upK)
	         cStatus[m_k] = SPxSolverBase<R>::FIXED;
	      else if(isZero(m_loK))
	         cStatus[m_k] = SPxSolverBase<R>::ON_LOWER;
	      else if(isZero(m_upK))
	         cStatus[m_k] = SPxSolverBase<R>::ON_UPPER;
	      else if(LErel(m_loK, R(0.0)) && GErel(m_upK, R(0.0)))
	         cStatus[m_k] = SPxSolverBase<R>::ZERO;
	      else
	         throw SPxInternalCodeException("XMAISM05 This should never happen.");
	
	      x[m_j] = 0.0;
	
	      if(isZero(m_loJ) && isZero(m_upJ) && m_loJ == m_upJ)
	         cStatus[m_j] = SPxSolverBase<R>::FIXED;
	      else if(isZero(m_loJ))
	         cStatus[m_j] = SPxSolverBase<R>::ON_LOWER;
	      else if(isZero(m_upJ))
	         cStatus[m_j] = SPxSolverBase<R>::ON_UPPER;
	      else if(LErel(m_loJ, R(0.0)) && GErel(m_upJ, R(0.0)))
	         cStatus[m_j] = SPxSolverBase<R>::ZERO;
	      else
	         throw SPxInternalCodeException("XMAISM06 This should never happen.");
	   }
	   else if(cStatus[m_k] == SPxSolverBase<R>::BASIC)
	   {
	      R scale1 = maxAbs(x[m_k], m_loK);
	      R scale2 = maxAbs(x[m_k], m_upK);
	
	      if(scale1 < 1.0)
	         scale1 = 1.0;
	
	      if(scale2 < 1.0)
	         scale2 = 1.0;
	
	      R z1 = (x[m_k] / scale1) - (m_loK / scale1);
	      R z2 = (x[m_k] / scale2) - (m_upK / scale2);
	
	      if(isZero(z1))
	         z1 = 0.0;
	
	      if(isZero(z2))
	         z2 = 0.0;
	
	      if(m_loJ <= R(-infinity) && m_upJ >= R(infinity) && m_loK <= R(-infinity) && m_upK >= R(infinity))
	      {
	         cStatus[m_j] = SPxSolverBase<R>::ZERO;
	         x[m_j] = 0.0;
	      }
	      else if(m_scale > 0.0)
	      {
	         if(GErel(x[m_k], m_upK + m_scale * m_upJ))
	         {
	            assert(m_upJ < R(infinity));
	            cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_UPPER;
	            x[m_j] = m_upJ;
	            x[m_k] -= m_scale * x[m_j];
	         }
	         else if(GErel(x[m_k], m_loK + m_scale * m_upJ) && m_upJ < R(infinity))
	         {
	            cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_UPPER;
	            x[m_j] = m_upJ;
	            x[m_k] -= m_scale * x[m_j];
	         }
	         else if(GErel(x[m_k], m_upK + m_scale * m_loJ) && m_upK < R(infinity))
	         {
	            cStatus[m_k] = (m_loK == m_upK) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_UPPER;
	            x[m_k] = m_upK;
	            cStatus[m_j] = SPxSolverBase<R>::BASIC;
	            x[m_j] = z2 * scale2 / m_scale;
	         }
	         else if(GErel(x[m_k], m_loK + m_scale * m_loJ) && m_loJ > R(-infinity))
	         {
	            cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_LOWER;
	            x[m_j] = m_loJ;
	            x[m_k] -= m_scale * x[m_j];
	         }
	         else if(GErel(x[m_k], m_loK + m_scale * m_loJ) && m_loK > R(-infinity))
	         {
	            cStatus[m_k] = (m_loK == m_upK) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_LOWER;
	            x[m_k] = m_loK;
	            cStatus[m_j] = SPxSolverBase<R>::BASIC;
	            x[m_j] = z1 * scale1 / m_scale;
	         }
	         else if(LTrel(x[m_k], m_loK + m_scale * m_loJ))
	         {
	            assert(m_loJ > R(-infinity));
	            cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_LOWER;
	            x[m_j] = m_loJ;
	            x[m_k] -= m_scale * x[m_j];
	         }
	         else
	         {
	            throw SPxInternalCodeException("XMAISM08 This should never happen.");
	         }
	      }
	      else
	      {
	         assert(m_scale < 0.0);
	
	         if(GErel(x[m_k], m_upK + m_scale * m_loJ))
	         {
	            assert(m_loJ > R(-infinity));
	            cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_LOWER;
	            x[m_j] = m_loJ;
	            x[m_k] -= m_scale * x[m_j];
	         }
	         else if(GErel(x[m_k], m_loK + m_scale * m_loJ) && m_loJ > R(-infinity))
	         {
	            cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_LOWER;
	            x[m_j] = m_loJ;
	            x[m_k] -= m_scale * x[m_j];
	         }
	         else if(GErel(x[m_k], m_upK + m_scale * m_upJ) && m_upK < R(infinity))
	         {
	            cStatus[m_k] = (m_loK == m_upK) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_UPPER;
	            x[m_k] = m_upK;
	            cStatus[m_j] = SPxSolverBase<R>::BASIC;
	            x[m_j] = z2 * scale2 / m_scale;
	         }
	         else if(GErel(x[m_k], m_loK + m_scale * m_upJ) && m_upJ < R(infinity))
	         {
	            cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_UPPER;
	            x[m_j] = m_upJ;
	            x[m_k] -= m_scale * x[m_j];
	         }
	         else if(GErel(x[m_k], m_loK + m_scale * m_upJ) && m_loK > R(-infinity))
	         {
	            cStatus[m_k] = (m_loK == m_upK) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_LOWER;
	            x[m_k] = m_loK;
	            cStatus[m_j] = SPxSolverBase<R>::BASIC;
	            x[m_j] = z1 * scale1 / m_scale;
	         }
	         else if(LTrel(x[m_k], m_loK + m_scale * m_upJ))
	         {
	            assert(m_upJ < R(infinity));
	            cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_UPPER;
	            x[m_j] = m_upJ;
	            x[m_k] -= m_scale * x[m_j];
	         }
	         else
	         {
	            throw SPxInternalCodeException("XMAISM09 This should never happen.");
	         }
	      }
	   }
	
	   // dual:
	   r[m_j] = m_scale * r[m_k];
	}
	
	template <class R>
	void SPxMainSM<R>::AggregationPS::execute(VectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& s,
	      VectorBase<R>& r,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
	{
	   // correcting the change of idx by deletion of the row:
	   if(m_i != m_old_i)
	   {
	      s[m_old_i] = s[m_i];
	      y[m_old_i] = y[m_i];
	      rStatus[m_old_i] = rStatus[m_i];
	   }
	
	   // correcting the change of idx by deletion of the column:
	   if(m_j != m_old_j)
	   {
	      x[m_old_j] = x[m_j];
	      r[m_old_j] = r[m_j];
	      cStatus[m_old_j] = cStatus[m_j];
	   }
	
	   // primal:
	   R val = 0.0;
	   R aij = m_row[m_j];
	   int active_idx = -1;
	
	   assert(m_row.size() == 2);
	
	   for(int k = 0; k < 2; ++k)
	   {
	      if(m_row.index(k) != m_j)
	      {
	         active_idx = m_row.index(k);
	         val = m_row.value(k) * x[active_idx];
	      }
	   }
	
	   assert(active_idx >= 0);
	
	   R scale = maxAbs(m_rhs, val);
	
	   if(scale < 1.0)
	      scale = 1.0;
	
	   R z = (m_rhs / scale) - (val / scale);
	
	   if(isZero(z))
	      z = 0.0;
	
	   x[m_j] = z * scale / aij;
	   s[m_i] = m_rhs;
	
	   if(isOptimal && (LT(x[m_j], m_lower, this->eps()) || GT(x[m_j], m_upper, this->eps())))
	   {
	      MSG_ERROR(std::cerr << "EMAISM: numerical violation after disaggregating variable" << std::endl;)
	   }
	
	   // dual:
	   R dualVal = 0.0;
	
	   for(int k = 0; k < m_col.size(); ++k)
	   {
	      if(m_col.index(k) != m_i)
	         dualVal += m_col.value(k) * y[m_col.index(k)];
	   }
	
	   z = m_obj - dualVal;
	
	   y[m_i] = z / aij;
	   r[m_j] = 0.0;
	
	   // basis:
	   if(((cStatus[active_idx] == SPxSolverBase<R>::ON_UPPER
	         || cStatus[active_idx] == SPxSolverBase<R>::FIXED)
	         && NE(x[active_idx], m_oldupper, this->eps())) ||
	         ((cStatus[active_idx] == SPxSolverBase<R>::ON_LOWER
	           || cStatus[active_idx] == SPxSolverBase<R>::FIXED)
	          && NE(x[active_idx], m_oldlower, this->eps())))
	   {
	      cStatus[active_idx] = SPxSolverBase<R>::BASIC;
	      r[active_idx] = 0.0;
	      assert(NE(m_upper, m_lower));
	
	      if(EQ(x[m_j], m_upper, this->eps()))
	         cStatus[m_j] = SPxSolverBase<R>::ON_UPPER;
	      else if(EQ(x[m_j], m_lower, this->eps()))
	         cStatus[m_j] = SPxSolverBase<R>::ON_LOWER;
	      else if(m_upper >= R(infinity) && m_lower <= R(-infinity))
	         cStatus[m_j] = SPxSolverBase<R>::ZERO;
	      else
	         throw SPxInternalCodeException("XMAISM unexpected basis status in aggregation unsimplifier.");
	   }
	   else
	   {
	      cStatus[m_j] = SPxSolverBase<R>::BASIC;
	   }
	
	   // sides may not be equal and we always only consider the rhs during aggregation, so set ON_UPPER
	   // (in theory and with exact arithmetic setting it to FIXED would be correct)
	   rStatus[m_i] = SPxSolverBase<R>::ON_UPPER;
	
	#ifdef CHECK_BASIC_DIM
	
	   if(!this->checkBasisDim(rStatus, cStatus))
	   {
	      throw SPxInternalCodeException("XMAISM22 Dimension doesn't match after this step.");
	   }
	
	#endif
	}
	
	template <class R>
	void SPxMainSM<R>::MultiAggregationPS::execute(VectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& s,
	      VectorBase<R>& r,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
	{
	
	   // correcting the change of idx by deletion of the row:
	   if(m_i != m_old_i)
	   {
	      s[m_old_i] = s[m_i];
	      y[m_old_i] = y[m_i];
	      rStatus[m_old_i] = rStatus[m_i];
	   }
	
	   // correcting the change of idx by deletion of the column:
	   if(m_j != m_old_j)
	   {
	      x[m_old_j] = x[m_j];
	      r[m_old_j] = r[m_j];
	      cStatus[m_old_j] = cStatus[m_j];
	   }
	
	   // primal:
	   R val = 0.0;
	   R aij = m_row[m_j];
	
	   for(int k = 0; k < m_row.size(); ++k)
	   {
	      if(m_row.index(k) != m_j)
	         val += m_row.value(k) * x[m_row.index(k)];
	   }
	
	   R scale = maxAbs(m_const, val);
	
	   if(scale < 1.0)
	      scale = 1.0;
	
	   R z = (m_const / scale) - (val / scale);
	
	   if(isZero(z))
	      z = 0.0;
	
	   x[m_j] = z * scale / aij;
	   s[m_i] = 0.0;
	
	#ifndef NDEBUG
	
	   if(isOptimal && (LT(x[m_j], m_lower, this->eps()) || GT(x[m_j], m_upper, this->eps())))
	      MSG_ERROR(std::cerr << "numerical violation in original space due to MultiAggregation\n";)
	#endif
	
	      // dual:
	      R dualVal = 0.0;
	
	   for(int k = 0; k < m_col.size(); ++k)
	   {
	      if(m_col.index(k) != m_i)
	         dualVal += m_col.value(k) * y[m_col.index(k)];
	   }
	
	   z = m_obj - dualVal;
	
	   y[m_i] = z / aij;
	   r[m_j] = 0.0;
	
	   // basis:
	   cStatus[m_j] = SPxSolverBase<R>::BASIC;
	
	   if(m_eqCons)
	      rStatus[m_i] = SPxSolverBase<R>::FIXED;
	   else if(m_onLhs)
	      rStatus[m_i] = SPxSolverBase<R>::ON_LOWER;
	   else
	      rStatus[m_i] = SPxSolverBase<R>::ON_UPPER;
	
	#ifdef CHECK_BASIC_DIM
	
	   if(!this->checkBasisDim(rStatus, cStatus))
	   {
	      throw SPxInternalCodeException("XMAISM22 Dimension doesn't match after this step.");
	   }
	
	#endif
	}
	
	template <class R>
	void SPxMainSM<R>::TightenBoundsPS::execute(VectorBase<R>& x, VectorBase<R>&, VectorBase<R>&,
	      VectorBase<R>&,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus,
	      DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const
	{
	   // basis:
	   switch(cStatus[m_j])
	   {
	   case SPxSolverBase<R>::FIXED:
	      if(LT(x[m_j], m_origupper, this->eps()) && GT(x[m_j], m_origlower, this->eps()))
	         cStatus[m_j] = SPxSolverBase<R>::BASIC;
	      else if(LT(x[m_j], m_origupper, this->eps()))
	         cStatus[m_j] = SPxSolverBase<R>::ON_LOWER;
	      else if(GT(x[m_j], m_origlower, this->eps()))
	         cStatus[m_j] = SPxSolverBase<R>::ON_UPPER;
	
	      break;
	
	   case SPxSolverBase<R>::ON_LOWER:
	      if(GT(x[m_j], m_origlower, this->eps()))
	         cStatus[m_j] = SPxSolverBase<R>::BASIC;
	
	      break;
	
	   case SPxSolverBase<R>::ON_UPPER:
	      if(LT(x[m_j], m_origupper, this->eps()))
	         cStatus[m_j] = SPxSolverBase<R>::BASIC;
	
	      break;
	
	   default:
	      break;
	   }
	
	#ifdef CHECK_BASIC_DIM
	
	   if(!this->checkBasisDim(rStatus, cStatus))
	   {
	      throw SPxInternalCodeException("XMAISM22 Dimension doesn't match after this step.");
	   }
	
	#endif
	}
	
	template <class R>
	void SPxMainSM<R>::handleRowObjectives(SPxLPBase<R>& lp)
	{
	   for(int i = lp.nRows() - 1; i >= 0; --i)
	   {
	      if(lp.maxRowObj(i) != 0.0)
	      {
	         std::shared_ptr<PostStep> ptr(new RowObjPS(lp, i, lp.nCols()));
	         m_hist.append(ptr);
	         lp.addCol(lp.rowObj(i), -lp.rhs(i), UnitVectorBase<R>(i), -lp.lhs(i));
	         lp.changeRange(i, R(0.0), R(0.0));
	         lp.changeRowObj(i, R(0.0));
	         m_addedcols++;
	      }
	   }
	}
	
	template <class R>
	void SPxMainSM<R>::handleExtremes(SPxLPBase<R>& lp)
	{
	
	   // This method handles extreme value of the given LP by
	   //
	   // 1. setting numbers of very small absolute values to zero and
	   // 2. setting numbers of very large absolute values to R(-infinity) or +R(infinity), respectively.
	
	   R maxVal  = R(infinity) / 5.0;
	   R tol = feastol() * 1e-2;
	   tol = (tol < this->epsZero()) ? this->epsZero() : tol;
	   int  remRows = 0;
	   int  remNzos = 0;
	   int  chgBnds = 0;
	   int  chgLRhs = 0;
	   int  objCnt  = 0;
	
	   for(int i = lp.nRows() - 1; i >= 0; --i)
	   {
	      // lhs
	      R lhs = lp.lhs(i);
	
	      if(lhs != 0.0 && isZero(lhs, this->epsZero()))
	      {
	         lp.changeLhs(i, R(0.0));
	         ++chgLRhs;
	      }
	      else if(lhs > R(-infinity) && lhs < -maxVal)
	      {
	         lp.changeLhs(i, R(-infinity));
	         ++chgLRhs;
	      }
	      else if(lhs <  R(infinity) && lhs >  maxVal)
	      {
	         lp.changeLhs(i,  R(infinity));
	         ++chgLRhs;
	      }
	
	      // rhs
	      R rhs = lp.rhs(i);
	
	      if(rhs != 0.0 && isZero(rhs, this->epsZero()))
	      {
	         lp.changeRhs(i, R(0.0));
	         ++chgLRhs;
	      }
	      else if(rhs > R(-infinity) && rhs < -maxVal)
	      {
	         lp.changeRhs(i, R(-infinity));
	         ++chgLRhs;
	      }
	      else if(rhs <  R(infinity) && rhs >  maxVal)
	      {
	         lp.changeRhs(i,  R(infinity));
	         ++chgLRhs;
	      }
	
	      if(lp.lhs(i) <= R(-infinity) && lp.rhs(i) >= R(infinity))
	      {
	         std::shared_ptr<PostStep> ptr(new FreeConstraintPS(lp, i));
	         m_hist.append(ptr);
	
	         removeRow(lp, i);
	         ++remRows;
	
	         ++m_stat[FREE_ROW];
	      }
	   }
	
	   for(int j = 0; j < lp.nCols(); ++j)
	   {
	      // lower bound
	      R lo = lp.lower(j);
	
	      if(lo != 0.0 && isZero(lo, this->epsZero()))
	      {
	         lp.changeLower(j, R(0.0));
	         ++chgBnds;
	      }
	      else if(lo > R(-infinity) && lo < -maxVal)
	      {
	         lp.changeLower(j, R(-infinity));
	         ++chgBnds;
	      }
	      else if(lo <  R(infinity) && lo >  maxVal)
	      {
	         lp.changeLower(j,  R(infinity));
	         ++chgBnds;
	      }
	
	      // upper bound
	      R up = lp.upper(j);
	
	      if(up != 0.0 && isZero(up, this->epsZero()))
	      {
	         lp.changeUpper(j, R(0.0));
	         ++chgBnds;
	      }
	      else if(up > R(-infinity) && up < -maxVal)
	      {
	         lp.changeUpper(j, R(-infinity));
	         ++chgBnds;
	      }
	      else if(up <  R(infinity) && up >  maxVal)
	      {
	         lp.changeUpper(j,  R(infinity));
	         ++chgBnds;
	      }
	
	      // fixed columns will be eliminated later
	      if(NE(lo, up))
	      {
	         lo = spxAbs(lo);
	         up = spxAbs(up);
	
	         R absBnd = (lo > up) ? lo : up;
	
	         if(absBnd < 1.0)
	            absBnd = 1.0;
	
	         // non-zeros
	         SVectorBase<R>& col = lp.colVector_w(j);
	         int        i = 0;
	
	         while(i < col.size())
	         {
	            R aij = spxAbs(col.value(i));
	
	            if(isZero(aij * absBnd, tol))
	            {
	               SVectorBase<R>& row = lp.rowVector_w(col.index(i));
	               int row_j = row.pos(j);
	
	               // this changes col.size()
	               if(row_j >= 0)
	                  row.remove(row_j);
	
	               col.remove(i);
	
	               MSG_DEBUG((*this->spxout) << "IMAISM04 aij=" << aij
	                         << " removed, absBnd=" << absBnd
	                         << std::endl;)
	               ++remNzos;
	            }
	            else
	            {
	               if(aij > maxVal)
	               {
	                  MSG_WARNING((*this->spxout), (*this->spxout) << "WMAISM05 Warning! Big matrix coefficient " << aij
	                              << std::endl);
	               }
	               else if(isZero(aij, tol))
	               {
	                  MSG_WARNING((*this->spxout), (*this->spxout) << "WMAISM06 Warning! Tiny matrix coefficient " << aij
	                              << std::endl);
	               }
	
	               ++i;
	            }
	         }
	      }
	
	      // objective
	      R obj = lp.obj(j);
	
	      if(obj != 0.0 && isZero(obj, this->epsZero()))
	      {
	         lp.changeObj(j, R(0.0));
	         ++objCnt;
	      }
	      else if(obj > R(-infinity) && obj < -maxVal)
	      {
	         lp.changeObj(j, R(-infinity));
	         ++objCnt;
	      }
	      else if(obj <  R(infinity) && obj >  maxVal)
	      {
	         lp.changeObj(j,  R(infinity));
	         ++objCnt;
	      }
	   }
	
	   if(remRows + remNzos + chgLRhs + chgBnds + objCnt > 0)
	   {
	      this->m_remRows += remRows;
	      this->m_remNzos += remNzos;
	      this->m_chgLRhs += chgLRhs;
	      this->m_chgBnds += chgBnds;
	
	      MSG_INFO2((*this->spxout), (*this->spxout) << "Simplifier (extremes) removed "
	                << remRows << " rows, "
	                << remNzos << " non-zeros, "
	                << chgBnds << " col bounds, "
	                << chgLRhs << " row bounds, "
	                << objCnt  << " objective coefficients" << std::endl;)
	   }
	
	   assert(lp.isConsistent());
	}
	
	/// computes the minimum and maximum residual activity for a given variable
	template <class R>
	void SPxMainSM<R>::computeMinMaxResidualActivity(SPxLPBase<R>& lp, int rowNumber, int colNumber,
	      R& minAct, R& maxAct)
	{
	   const SVectorBase<R>& row = lp.rowVector(rowNumber);
	   bool minNegInfinite = false;
	   bool maxInfinite = false;
	
	   minAct = 0;   // this is the minimum value that the aggregation can attain
	   maxAct = 0;   // this is the maximum value that the aggregation can attain
	
	   for(int l = 0; l < row.size(); ++l)
	   {
	      if(colNumber < 0 || row.index(l) != colNumber)
	      {
	         // computing the minimum activity of the aggregated variables
	         if(GT(row.value(l), R(0.0)))
	         {
	            if(LE(lp.lower(row.index(l)), R(-infinity)))
	               minNegInfinite = true;
	            else
	               minAct += row.value(l) * lp.lower(row.index(l));
	         }
	         else if(LT(row.value(l), R(0.0)))
	         {
	            if(GE(lp.upper(row.index(l)), R(infinity)))
	               minNegInfinite = true;
	            else
	               minAct += row.value(l) * lp.upper(row.index(l));
	         }
	
	         // computing the maximum activity of the aggregated variables
	         if(GT(row.value(l), R(0.0)))
	         {
	            if(GE(lp.upper(row.index(l)), R(infinity)))
	               maxInfinite = true;
	            else
	               maxAct += row.value(l) * lp.upper(row.index(l));
	         }
	         else if(LT(row.value(l), R(0.0)))
	         {
	            if(LE(lp.lower(row.index(l)), R(-infinity)))
	               maxInfinite = true;
	            else
	               maxAct += row.value(l) * lp.lower(row.index(l));
	         }
	      }
	   }
	
	   // if an infinite value exists for the minimum activity, then that it taken
	   if(minNegInfinite)
	      minAct = R(-infinity);
	
	   // if an -infinite value exists for the maximum activity, then that value is taken
	   if(maxInfinite)
	      maxAct = R(infinity);
	}
	
	
	/// calculate min/max value for the multi aggregated variables
	template <class R>
	void SPxMainSM<R>::computeMinMaxValues(SPxLPBase<R>& lp, R side, R val, R minRes, R maxRes,
	                                       R& minVal, R& maxVal)
	{
	   minVal = 0;
	   maxVal = 0;
	
	   if(LT(val, R(0.0)))
	   {
	      if(LE(minRes, R(-infinity)))
	         minVal = R(-infinity);
	      else
	         minVal = (side - minRes) / val;
	
	      if(GE(maxRes, R(infinity)))
	         maxVal = R(infinity);
	      else
	         maxVal = (side - maxRes) / val;
	   }
	   else if(GT(val, R(0.0)))
	   {
	      if(GE(maxRes, R(infinity)))
	         minVal = R(-infinity);
	      else
	         minVal = (side - maxRes) / val;
	
	      if(LE(minRes, R(-infinity)))
	         maxVal = R(infinity);
	      else
	         maxVal = (side - minRes) / val;
	   }
	}
	
	
	/// tries to find good lower bound solutions by applying some trivial heuristics
	template <class R>
	void SPxMainSM<R>::trivialHeuristic(SPxLPBase<R>& lp)
	{
	   VectorBase<R>         zerosol(lp.nCols());  // the zero solution VectorBase<R>
	   VectorBase<R>         lowersol(lp.nCols()); // the lower bound solution VectorBase<R>
	   VectorBase<R>         uppersol(lp.nCols()); // the upper bound solution VectorBase<R>
	   VectorBase<R>         locksol(lp.nCols());  // the locks solution VectorBase<R>
	
	   VectorBase<R>         upLocks(lp.nCols());
	   VectorBase<R>         downLocks(lp.nCols());
	
	   R            zeroObj = this->m_objoffset;
	   R            lowerObj = this->m_objoffset;
	   R            upperObj = this->m_objoffset;
	   R            lockObj = this->m_objoffset;
	
	   bool            zerovalid = true;
	
	   R largeValue = R(infinity);
	
	   if(LT(R(1.0 / feastol()), R(infinity)))
	      largeValue = 1.0 / feastol();
	
	
	
	   for(int j = lp.nCols() - 1; j >= 0; --j)
	   {
	      upLocks[j] = 0;
	      downLocks[j] = 0;
	
	      // computing the locks on the variables
	      const SVectorBase<R>& col = lp.colVector(j);
	
	      for(int k = 0; k < col.size(); ++k)
	      {
	         R aij = col.value(k);
	
	         ASSERT_WARN("WMAISM45", isNotZero(aij, R(1.0 / R(infinity))));
	
	         if(GT(lp.lhs(col.index(k)), R(-infinity)) && LT(lp.rhs(col.index(k)), R(infinity)))
	         {
	            upLocks[j]++;
	            downLocks[j]++;
	         }
	         else if(GT(lp.lhs(col.index(k)), R(-infinity)))
	         {
	            if(aij > 0)
	               downLocks[j]++;
	            else if(aij < 0)
	               upLocks[j]++;
	         }
	         else if(LT(lp.rhs(col.index(k)), R(infinity)))
	         {
	            if(aij > 0)
	               upLocks[j]++;
	            else if(aij < 0)
	               downLocks[j]++;
	         }
	      }
	
	      R lower = lp.lower(j);
	      R upper = lp.upper(j);
	
	      if(LE(lower, R(-infinity)))
	         lower = MINIMUM(-largeValue, upper);
	
	      if(GE(upper, R(infinity)))
	         upper = MAXIMUM(lp.lower(j), largeValue);
	
	      if(zerovalid)
	      {
	         if(LE(lower, R(0.0), feastol()) && GE(upper, R(0.0), feastol()))
	            zerosol[j] = 0.0;
	         else
	            zerovalid = false;
	      }
	
	      lowersol[j] = lower;
	      uppersol[j] = upper;
	
	      if(downLocks[j] > upLocks[j])
	         locksol[j] = upper;
	      else if(downLocks[j] < upLocks[j])
	         locksol[j] = lower;
	      else
	         locksol[j] = (lower + upper) / 2.0;
	
	      lowerObj += lp.maxObj(j) * lowersol[j];
	      upperObj += lp.maxObj(j) * uppersol[j];
	      lockObj += lp.maxObj(j) * locksol[j];
	   }
	
	   // trying the lower bound solution
	   if(checkSolution(lp, lowersol))
	   {
	      if(lowerObj > m_cutoffbound)
	         m_cutoffbound = lowerObj;
	   }
	
	   // trying the upper bound solution
	   if(checkSolution(lp, uppersol))
	   {
	      if(upperObj > m_cutoffbound)
	         m_cutoffbound = upperObj;
	   }
	
	   // trying the zero solution
	   if(zerovalid && checkSolution(lp, zerosol))
	   {
	      if(zeroObj > m_cutoffbound)
	         m_cutoffbound = zeroObj;
	   }
	
	   // trying the lock solution
	   if(checkSolution(lp, locksol))
	   {
	      if(lockObj > m_cutoffbound)
	         m_cutoffbound = lockObj;
	   }
	}
	
	
	
	/// checks a solution for feasibility
	template <class R>
	bool SPxMainSM<R>::checkSolution(SPxLPBase<R>& lp, VectorBase<R> sol)
	{
	   for(int i = lp.nRows() - 1; i >= 0; --i)
	   {
	      const SVectorBase<R>& row = lp.rowVector(i);
	      R activity = 0;
	
	      for(int k = 0; k < row.size(); k++)
	         activity += row.value(k) * sol[row.index(k)];
	
	      if(!GE(activity, lp.lhs(i), feastol()) || !LE(activity, lp.rhs(i), feastol()))
	         return false;
	   }
	
	   return true;
	}
	
	
	/// tightens variable bounds by propagating the pseudo objective function value.
	template <class R>
	void SPxMainSM<R>::propagatePseudoobj(SPxLPBase<R>& lp)
	{
	   R pseudoObj = this->m_objoffset;
	
	   for(int j = lp.nCols() - 1; j >= 0; --j)
	   {
	      R val = lp.maxObj(j);
	
	      if(val < 0)
	      {
	         if(lp.lower(j) <= R(-infinity))
	            return;
	
	         pseudoObj += val * lp.lower(j);
	      }
	      else if(val > 0)
	      {
	         if(lp.upper(j) >= R(-infinity))
	            return;
	
	         pseudoObj += val * lp.upper(j);
	      }
	   }
	
	   if(GT(m_cutoffbound, R(-infinity)) && LT(m_cutoffbound, R(infinity)))
	   {
	      if(pseudoObj > m_pseudoobj)
	         m_pseudoobj = pseudoObj;
	
	      for(int j = lp.nCols() - 1; j >= 0; --j)
	      {
	         R objval = lp.maxObj(j);
	
	         if(EQ(objval, R(0.0)))
	            continue;
	
	         if(objval < 0.0)
	         {
	            R newbound = lp.lower(j) + (m_cutoffbound - m_pseudoobj) / objval;
	
	            if(LT(newbound, lp.upper(j)))
	            {
	               std::shared_ptr<PostStep> ptr(new TightenBoundsPS(lp, j, lp.upper(j), lp.lower(j)));
	               m_hist.append(ptr);
	               lp.changeUpper(j, newbound);
	            }
	         }
	         else if(objval > 0.0)
	         {
	            R newbound = lp.upper(j) + (m_cutoffbound - m_pseudoobj) / objval;
	
	            if(GT(newbound, lp.lower(j)))
	            {
	               std::shared_ptr<PostStep> ptr(new TightenBoundsPS(lp, j, lp.upper(j), lp.lower(j)));
	               m_hist.append(ptr);
	               lp.changeLower(j, newbound);
	            }
	         }
	      }
	   }
	}
	
	
	
	template <class R>
	typename SPxSimplifier<R>::Result SPxMainSM<R>::removeEmpty(SPxLPBase<R>& lp)
	{
	
	   // This method removes empty rows and columns from the LP.
	
	   int remRows = 0;
	   int remCols = 0;
	
	   for(int i = lp.nRows() - 1; i >= 0; --i)
	   {
	      const SVectorBase<R>& row = lp.rowVector(i);
	
	      if(row.size() == 0)
	      {
	         MSG_DEBUG((*this->spxout) << "IMAISM07 row " << i
	                   << ": empty ->";)
	
	         if(LT(lp.rhs(i), R(0.0), feastol()) || GT(lp.lhs(i), R(0.0), feastol()))
	         {
	            MSG_DEBUG((*this->spxout) << " infeasible lhs=" << lp.lhs(i)
	                      << " rhs=" << lp.rhs(i) << std::endl;)
	            return this->INFEASIBLE;
	         }
	
	         MSG_DEBUG((*this->spxout) << " removed" << std::endl;)
	
	         std::shared_ptr<PostStep> ptr(new EmptyConstraintPS(lp, i));
	         m_hist.append(ptr);
	
	         ++remRows;
	         removeRow(lp, i);
	
	         ++m_stat[EMPTY_ROW];
	      }
	   }
	
	   for(int j = lp.nCols() - 1; j >= 0; --j)
	   {
	      const SVectorBase<R>& col = lp.colVector(j);
	
	      if(col.size() == 0)
	      {
	         MSG_DEBUG((*this->spxout) << "IMAISM08 col " << j
	                   << ": empty -> maxObj=" << lp.maxObj(j)
	                   << " lower=" << lp.lower(j)
	                   << " upper=" << lp.upper(j);)
	
	         R val;
	
	         if(GT(lp.maxObj(j), R(0.0), this->epsZero()))
	         {
	            if(lp.upper(j) >= R(infinity))
	            {
	               MSG_DEBUG((*this->spxout) << " unbounded" << std::endl;)
	               return this->UNBOUNDED;
	            }
	
	            val = lp.upper(j);
	         }
	         else if(LT(lp.maxObj(j), R(0.0), this->epsZero()))
	         {
	            if(lp.lower(j) <= R(-infinity))
	            {
	               MSG_DEBUG((*this->spxout) << " unbounded" << std::endl;)
	               return this->UNBOUNDED;
	            }
	
	            val = lp.lower(j);
	         }
	         else
	         {
	            ASSERT_WARN("WMAISM09", isZero(lp.maxObj(j), this->epsZero()));
	
	            // any value within the bounds is ok
	            if(lp.lower(j) > R(-infinity))
	               val = lp.lower(j);
	            else if(lp.upper(j) < R(infinity))
	               val = lp.upper(j);
	            else
	               val = 0.0;
	         }
	
	         MSG_DEBUG((*this->spxout) << " removed" << std::endl;)
	
	         std::shared_ptr<PostStep> ptr1(new FixBoundsPS(lp, j, val));
	         std::shared_ptr<PostStep> ptr2(new FixVariablePS(lp, *this, j, val));
	         m_hist.append(ptr1);
	         m_hist.append(ptr2);
	
	         ++remCols;
	         removeCol(lp, j);
	
	         ++m_stat[EMPTY_COL];
	      }
	   }
	
	   if(remRows + remCols > 0)
	   {
	      this->m_remRows += remRows;
	      this->m_remCols += remCols;
	
	      MSG_INFO2((*this->spxout), (*this->spxout) << "Simplifier (empty rows/colums) removed "
	                << remRows << " rows, "
	                << remCols << " cols"
	                << std::endl;)
	
	   }
	
	   return this->OKAY;
	}
	
	template <class R>
	typename SPxSimplifier<R>::Result SPxMainSM<R>::removeRowSingleton(SPxLPBase<R>& lp,
	      const SVectorBase<R>& row, int& i)
	{
	   assert(row.size() == 1);
	
	   R aij = row.value(0);
	   int  j   = row.index(0);
	   R lo  = R(-infinity);
	   R up  =  R(infinity);
	
	   MSG_DEBUG((*this->spxout) << "IMAISM22 row " << i
	             << ": singleton -> val=" << aij
	             << " lhs=" << lp.lhs(i)
	             << " rhs=" << lp.rhs(i);)
	
	   if(GT(aij, R(0.0), this->epsZero()))            // aij > 0
	   {
	      lo = (lp.lhs(i) <= R(-infinity)) ? R(-infinity) : lp.lhs(i) / aij;
	      up = (lp.rhs(i) >=  R(infinity)) ?  R(infinity) : lp.rhs(i) / aij;
	   }
	   else if(LT(aij, R(0.0), this->epsZero()))       // aij < 0
	   {
	      lo = (lp.rhs(i) >=  R(infinity)) ? R(-infinity) : lp.rhs(i) / aij;
	      up = (lp.lhs(i) <= R(-infinity)) ?  R(infinity) : lp.lhs(i) / aij;
	   }
	   else if(LT(lp.rhs(i), R(0.0), feastol()) || GT(lp.lhs(i), R(0.0), feastol()))
	   {
	      // aij == 0, rhs < 0 or lhs > 0
	      MSG_DEBUG((*this->spxout) << " infeasible" << std::endl;)
	      return this->INFEASIBLE;
	   }
	
	   if(isZero(lo, this->epsZero()))
	      lo = 0.0;
	
	   if(isZero(up, this->epsZero()))
	      up = 0.0;
	
	   MSG_DEBUG((*this->spxout) << " removed, lower=" << lo
	             << " (" << lp.lower(j)
	             << ") upper=" << up
	             << " (" << lp.upper(j)
	             << ")" << std::endl;)
	
	   bool stricterUp = false;
	   bool stricterLo = false;
	
	   R oldLo = lp.lower(j);
	   R oldUp = lp.upper(j);
	
	   if(LTrel(up, lp.upper(j), feastol()))
	   {
	      lp.changeUpper(j, up);
	      stricterUp = true;
	   }
	
	   if(GTrel(lo, lp.lower(j), feastol()))
	   {
	      lp.changeLower(j, lo);
	      stricterLo = true;
	   }
	
	   std::shared_ptr<PostStep> ptr(new RowSingletonPS(lp, i, j, stricterLo, stricterUp, lp.lower(j),
	                                 lp.upper(j), oldLo, oldUp));
	   m_hist.append(ptr);
	
	   removeRow(lp, i);
	
	   this->m_remRows++;
	   this->m_remNzos++;
	   ++m_stat[SINGLETON_ROW];
	
	   return this->OKAY;
	}
	
	/// aggregate variable x_j to x_j = (rhs - aik * x_k) / aij from row i: aij * x_j + aik * x_k = rhs
	template <class R>
	typename SPxSimplifier<R>::Result SPxMainSM<R>::aggregateVars(SPxLPBase<R>& lp,
	      const SVectorBase<R>& row, int& i)
	{
	   assert(row.size() == 2);
	   assert(EQrel(lp.lhs(i), lp.rhs(i), feastol()));
	
	   R rhs = lp.rhs(i);
	   assert(rhs < R(infinity) && rhs > R(-infinity));
	
	   int j = row.index(0);
	   int k = row.index(1);
	   R aij = row.value(0);
	   R aik = row.value(1);
	   R lower_j = lp.lower(j);
	   R upper_j = lp.upper(j);
	   R lower_k = lp.lower(k);
	   R upper_k = lp.upper(k);
	
	   // fixed variables should be removed by simplifyCols()
	   if(EQrel(lower_j, upper_j, feastol()) || EQrel(lower_k, upper_k, feastol()))
	      return this->OKAY;
	
	   assert(isNotZero(aij, this->epsZero()) && isNotZero(aik, this->epsZero()));
	
	   MSG_DEBUG((*this->spxout) << "IMAISM22 row " << i << ": doubleton equation -> "
	             << aij << " x_" << j << " + " << aik << " x_" << k << " = " << rhs;)
	
	   // determine which variable can be aggregated without requiring bound tightening of the other variable
	   R new_lo_j;
	   R new_up_j;
	   R new_lo_k;
	   R new_up_k;
	
	   if(aij * aik < 0.0)
	   {
	      // orientation persists
	      new_lo_j = (upper_k >=  R(infinity)) ? R(-infinity) : (rhs - aik * upper_k) / aij;
	      new_up_j = (lower_k <= R(-infinity)) ?  R(infinity) : (rhs - aik * lower_k) / aij;
	      new_lo_k = (upper_j >=  R(infinity)) ? R(-infinity) : (rhs - aij * upper_j) / aik;
	      new_up_k = (lower_j <= R(-infinity)) ?  R(infinity) : (rhs - aij * lower_j) / aik;
	   }
	   else if(aij * aik > 0.0)
	   {
	      // orientation is reversed
	      new_lo_j = (lower_k <= R(-infinity)) ? R(-infinity) : (rhs - aik * lower_k) / aij;
	      new_up_j = (upper_k >=  R(infinity)) ?  R(infinity) : (rhs - aik * upper_k) / aij;
	      new_lo_k = (lower_j <= R(-infinity)) ? R(-infinity) : (rhs - aij * lower_j) / aik;
	      new_up_k = (upper_j >=  R(infinity)) ?  R(infinity) : (rhs - aij * upper_j) / aik;
	   }
	   else
	      throw SPxInternalCodeException("XMAISM12 This should never happen.");
	
	   bool flip_jk = false;
	
	   if(new_lo_j <= R(-infinity) && new_up_j >= R(infinity))
	   {
	      // no bound tightening on x_j when x_k is aggregated
	      flip_jk = true;
	   }
	   else if(new_lo_k <= R(-infinity) && new_up_k >= R(infinity))
	   {
	      // no bound tightening on x_k when x_j is aggregated
	      flip_jk = false;
	   }
	   else if(LE(new_lo_j, lower_j) && GE(new_up_j, upper_j))
	   {
	      if(LE(new_lo_k, lower_k) && GE(new_up_k, upper_k))
	      {
	         // both variables' bounds are not affected by aggregation; choose the better aggregation coeff (aik/aij)
	         if(spxAbs(aij) > spxAbs(aik))
	            flip_jk = false;
	         else
	            flip_jk = true;
	      }
	      else
	         flip_jk = false;
	   }
	   else if(LE(new_lo_k, lower_k) && GE(new_up_k, upper_k))
	   {
	      flip_jk = true;
	   }
	   else
	   {
	      if(spxAbs(aij) > spxAbs(aik))
	         flip_jk = false;
	      else
	         flip_jk = true;
	   }
	
	   if(flip_jk)
	   {
	      int _j = j;
	      R _aij = aij;
	      R _lower_j = lower_j;
	      R _upper_j = upper_j;
	      j = k;
	      k = _j;
	      aij = aik;
	      aik = _aij;
	      lower_j = lower_k;
	      lower_k = _lower_j;
	      upper_j = upper_k;
	      upper_k = _upper_j;
	   }
	
	   const SVectorBase<R>& col_j = lp.colVector(j);
	   const SVectorBase<R>& col_k = lp.colVector(k);
	
	   // aggregation coefficients (x_j = aggr_coef * x_k + aggr_const)
	   R aggr_coef = - (aik / aij);
	   R aggr_const = rhs / aij;
	
	   MSG_DEBUG((*this->spxout) << " removed, replacing x_" << j << " with "
	             << aggr_const << " + " << aggr_coef << " * x_" << k << std::endl;)
	
	   // replace all occurrences of x_j
	   for(int r = 0; r < col_j.size(); ++r)
	   {
	      int row_r = col_j.index(r);
	      R arj = col_j.value(r);
	
	      // skip row i
	      if(row_r == i)
	         continue;
	
	      // adapt sides of row r
	      R lhs_r = lp.lhs(row_r);
	      R rhs_r = lp.rhs(row_r);
	
	      if(lhs_r > R(-infinity))
	      {
	         lp.changeLhs(row_r, lhs_r - aggr_const * arj);
	         this->m_chgLRhs++;
	      }
	
	      if(rhs_r < R(infinity))
	      {
	         lp.changeRhs(row_r, rhs_r - aggr_const * arj);
	         this->m_chgLRhs++;
	      }
	
	      R newcoef = aggr_coef * arj;
	      int pos_rk = col_k.pos(row_r);
	
	      // check whether x_k is also present in row r and get its coefficient
	      if(pos_rk >= 0)
	      {
	         R ark = col_k.value(pos_rk);
	         newcoef += ark;
	         this->m_remNzos++;
	      }
	
	      // add new column k to row r or adapt the coefficient a_rk
	      lp.changeElement(row_r, k, newcoef);
	   }
	
	   // adapt objective function
	   R obj_j = lp.obj(j);
	
	   if(isNotZero(obj_j, this->epsZero()))
	   {
	      this->addObjoffset(aggr_const * obj_j);
	      R obj_k = lp.obj(k);
	      lp.changeObj(k, obj_k + aggr_coef * obj_j);
	   }
	
	   // adapt bounds of x_k
	   R scale1 = maxAbs(rhs, aij * upper_j);
	   R scale2 = maxAbs(rhs, aij * lower_j);
	
	   if(scale1 < 1.0)
	      scale1 = 1.0;
	
	   if(scale2 < 1.0)
	      scale2 = 1.0;
	
	   R z1 = (rhs / scale1) - (aij * upper_j / scale1);
	   R z2 = (rhs / scale2) - (aij * lower_j / scale2);
	
	   // just some rounding
	   if(isZero(z1, this->epsZero()))
	      z1 = 0.0;
	
	   if(isZero(z2, this->epsZero()))
	      z2 = 0.0;
	
	   // determine which side has to be used for the bounds comparison below
	   if(aik * aij > 0.0)
	   {
	      new_lo_k = (upper_j >=  R(infinity)) ? R(-infinity) : z1 * scale1 / aik;
	      new_up_k = (lower_j <= R(-infinity)) ?  R(infinity) : z2 * scale2 / aik;
	   }
	   else if(aik * aij < 0.0)
	   {
	      new_lo_k = (lower_j <= R(-infinity)) ? R(-infinity) : z2 * scale2 / aik;
	      new_up_k = (upper_j >=  R(infinity)) ?  R(infinity) : z1 * scale1 / aik;
	   }
	   else
	      throw SPxInternalCodeException("XMAISM12 This should never happen.");
	
	   // change bounds of x_k if the new ones are tighter
	   R oldlower_k = lower_k;
	   R oldupper_k = upper_k;
	
	   if(GT(new_lo_k, lower_k, this->epsZero()))
	   {
	      lp.changeLower(k, new_lo_k);
	      this->m_chgBnds++;
	   }
	
	   if(LT(new_up_k, upper_k, this->epsZero()))
	   {
	      lp.changeUpper(k, new_up_k);
	      this->m_chgBnds++;
	   }
	
	   std::shared_ptr<PostStep> ptr(new AggregationPS(lp, i, j, rhs, oldupper_k, oldlower_k));
	   m_hist.append(ptr);
	
	   removeRow(lp, i);
	   removeCol(lp, j);
	
	   this->m_remRows++;
	   this->m_remCols++;
	   this->m_remNzos += 2;
	
	   ++m_stat[AGGREGATION];
	
	   return this->OKAY;
	}
	
	template <class R>
	typename SPxSimplifier<R>::Result SPxMainSM<R>::simplifyRows(SPxLPBase<R>& lp, bool& again)
	{
	
	   // This method simplifies the rows of the LP.
	   //
	   // The following operations are done:
	   // 1. detect implied free variables
	   // 2. detect implied free constraints
	   // 3. detect infeasible constraints
	   // 4. remove unconstrained constraints
	   // 5. remove empty constraints
	   // 6. remove row singletons and tighten the corresponding variable bounds if necessary
	   // 7. remove doubleton equation, aka aggregation
	   // 8. detect forcing rows and fix the corresponding variables
	
	   int remRows = 0;
	   int remNzos = 0;
	   int chgLRhs = 0;
	   int chgBnds = 0;
	   int keptBnds = 0;
	   int keptLRhs = 0;
	
	   int oldRows = lp.nRows();
	
	   bool redundantLower;
	   bool redundantUpper;
	   bool redundantLhs;
	   bool redundantRhs;
	
	   for(int i = lp.nRows() - 1; i >= 0; --i)
	   {
	      const SVectorBase<R>& row = lp.rowVector(i);
	
	
	      // compute bounds on constraint value
	      R lhsBnd = 0.0; // minimal activity (finite summands)
	      R rhsBnd = 0.0; // maximal activity (finite summands)
	      int  lhsCnt = 0; // number of R(-infinity) summands in minimal activity
	      int  rhsCnt = 0; // number of +R(infinity) summands in maximal activity
	
	      for(int k = 0; k < row.size(); ++k)
	      {
	         R aij = row.value(k);
	         int  j   = row.index(k);
	
	         if(!isNotZero(aij, R(1.0 / R(infinity))))
	         {
	            MSG_WARNING((*this->spxout), (*this->spxout) << "Warning: tiny nonzero coefficient " << aij <<
	                        " in row " << i << "\n");
	         }
	
	         if(aij > 0.0)
	         {
	            if(lp.lower(j) <= R(-infinity))
	               ++lhsCnt;
	            else
	               lhsBnd += aij * lp.lower(j);
	
	            if(lp.upper(j) >= R(infinity))
	               ++rhsCnt;
	            else
	               rhsBnd += aij * lp.upper(j);
	         }
	         else if(aij < 0.0)
	         {
	            if(lp.lower(j) <= R(-infinity))
	               ++rhsCnt;
	            else
	               rhsBnd += aij * lp.lower(j);
	
	            if(lp.upper(j) >= R(infinity))
	               ++lhsCnt;
	            else
	               lhsBnd += aij * lp.upper(j);
	         }
	      }
	
	#if FREE_BOUNDS
	
	      // 1. detect implied free variables
	      if(rhsCnt <= 1 || lhsCnt <= 1)
	      {
	         for(int k = 0; k < row.size(); ++k)
	         {
	            R aij = row.value(k);
	            int  j   = row.index(k);
	
	            redundantLower = false;
	            redundantUpper = false;
	
	            ASSERT_WARN("WMAISM12", isNotZero(aij, R(1.0 / R(infinity))));
	
	            if(aij > 0.0)
	            {
	               if(lp.lhs(i) > R(-infinity) && lp.lower(j) > R(-infinity) && rhsCnt <= 1
	                     && NErel(lp.lhs(i), rhsBnd, feastol())
	                     // do not perform if strongly different orders of magnitude occur
	                     && spxAbs(lp.lhs(i) / maxAbs(rhsBnd, R(1.0))) > Param::epsilon())
	               {
	                  R lo    = R(-infinity);
	                  R scale = maxAbs(lp.lhs(i), rhsBnd);
	
	                  if(scale < 1.0)
	                     scale = 1.0;
	
	                  R z = (lp.lhs(i) / scale) - (rhsBnd / scale);
	
	                  if(isZero(z, this->epsZero()))
	                     z = 0.0;
	
	                  assert(rhsCnt > 0 || lp.upper(j) < R(infinity));
	
	                  if(rhsCnt == 0)
	                     lo = lp.upper(j) + z * scale / aij;
	                  else if(lp.upper(j) >= R(infinity))
	                     lo = z * scale / aij;
	
	                  if(isZero(lo, this->epsZero()))
	                     lo = 0.0;
	
	                  if(GErel(lo, lp.lower(j), feastol()))
	                  {
	                     MSG_DEBUG((*this->spxout) << "IMAISM13 row " << i
	                               << ": redundant lower bound on x" << j
	                               << " -> lower=" << lo
	                               << " (" << lp.lower(j)
	                               << ")" << std::endl;)
	
	                     redundantLower = true;
	                  }
	
	               }
	
	               if(lp.rhs(i) < R(infinity) && lp.upper(j) < R(infinity) && lhsCnt <= 1
	                     && NErel(lp.rhs(i), lhsBnd, feastol())
	                     // do not perform if strongly different orders of magnitude occur
	                     && spxAbs(lp.rhs(i) / maxAbs(lhsBnd, R(1.0))) > Param::epsilon())
	               {
	                  R up    = R(infinity);
	                  R scale = maxAbs(lp.rhs(i), lhsBnd);
	
	                  if(scale < 1.0)
	                     scale = 1.0;
	
	                  R z = (lp.rhs(i) / scale) - (lhsBnd / scale);
	
	                  if(isZero(z, this->epsZero()))
	                     z = 0.0;
	
	                  assert(lhsCnt > 0 || lp.lower(j) > R(-infinity));
	
	                  if(lhsCnt == 0)
	                     up = lp.lower(j) + z * scale / aij;
	                  else if(lp.lower(j) <= R(-infinity))
	                     up = z * scale / aij;
	
	                  if(isZero(up, this->epsZero()))
	                     up = 0.0;
	
	                  if(LErel(up, lp.upper(j), feastol()))
	                  {
	                     MSG_DEBUG((*this->spxout) << "IMAISM14 row " << i
	                               << ": redundant upper bound on x" << j
	                               << " -> upper=" << up
	                               << " (" << lp.upper(j)
	                               << ")" << std::endl;)
	
	                     redundantUpper = true;
	                  }
	               }
	
	               if(redundantLower)
	               {
	                  // no upper bound on x_j OR redundant upper bound
	                  if((lp.upper(j) >= R(infinity)) || redundantUpper || (!m_keepbounds))
	                  {
	                     ++lhsCnt;
	                     lhsBnd -= aij * lp.lower(j);
	
	                     lp.changeLower(j, R(-infinity));
	                     ++chgBnds;
	                  }
	                  else
	                     ++keptBnds;
	               }
	
	               if(redundantUpper)
	               {
	                  // no lower bound on x_j OR redundant lower bound
	                  if((lp.lower(j) <= R(-infinity)) || redundantLower || (!m_keepbounds))
	                  {
	                     ++rhsCnt;
	                     rhsBnd -= aij * lp.upper(j);
	
	                     lp.changeUpper(j, R(infinity));
	                     ++chgBnds;
	                  }
	                  else
	                     ++keptBnds;
	               }
	            }
	            else if(aij < 0.0)
	            {
	               if(lp.lhs(i) > R(-infinity) && lp.upper(j) < R(infinity) && rhsCnt <= 1
	                     && NErel(lp.lhs(i), rhsBnd, feastol())
	                     // do not perform if strongly different orders of magnitude occur
	                     && spxAbs(lp.lhs(i) / maxAbs(rhsBnd, R(1.0))) > Param::epsilon())
	               {
	                  R up    = R(infinity);
	                  R scale = maxAbs(lp.lhs(i), rhsBnd);
	
	                  if(scale < 1.0)
	                     scale = 1.0;
	
	                  R z = (lp.lhs(i) / scale) - (rhsBnd / scale);
	
	                  if(isZero(z, this->epsZero()))
	                     z = 0.0;
	
	                  assert(rhsCnt > 0 || lp.lower(j) > R(-infinity));
	
	                  if(rhsCnt == 0)
	                     up = lp.lower(j) + z * scale / aij;
	                  else if(lp.lower(j) <= R(-infinity))
	                     up = z * scale / aij;
	
	                  if(isZero(up, this->epsZero()))
	                     up = 0.0;
	
	                  if(LErel(up, lp.upper(j), feastol()))
	                  {
	                     MSG_DEBUG((*this->spxout) << "IMAISM15 row " << i
	                               << ": redundant upper bound on x" << j
	                               << " -> upper=" << up
	                               << " (" << lp.upper(j)
	                               << ")" << std::endl;)
	
	                     redundantUpper = true;
	                  }
	               }
	
	               if(lp.rhs(i) < R(infinity) && lp.lower(j) > R(-infinity) && lhsCnt <= 1
	                     && NErel(lp.rhs(i), lhsBnd, feastol())
	                     // do not perform if strongly different orders of magnitude occur
	                     && spxAbs(lp.rhs(i) / maxAbs(lhsBnd, R(1.0))) > Param::epsilon())
	               {
	                  R lo    = R(-infinity);
	                  R scale = maxAbs(lp.rhs(i), lhsBnd);
	
	                  if(scale < 1.0)
	                     scale = 1.0;
	
	                  R z = (lp.rhs(i) / scale) - (lhsBnd / scale);
	
	                  if(isZero(z, this->epsZero()))
	                     z = 0.0;
	
	                  assert(lhsCnt > 0 || lp.upper(j) < R(infinity));
	
	                  if(lhsCnt == 0)
	                     lo = lp.upper(j) + z * scale / aij;
	                  else if(lp.upper(j) >= R(infinity))
	                     lo = z * scale / aij;
	
	                  if(isZero(lo, this->epsZero()))
	                     lo = 0.0;
	
	                  if(GErel(lo, lp.lower(j)))
	                  {
	                     MSG_DEBUG((*this->spxout) << "IMAISM16 row " << i
	                               << ": redundant lower bound on x" << j
	                               << " -> lower=" << lo
	                               << " (" << lp.lower(j)
	                               << ")" << std::endl;)
	
	                     redundantLower = true;
	                  }
	               }
	
	               if(redundantUpper)
	               {
	                  // no lower bound on x_j OR redundant lower bound
	                  if((lp.lower(j) <= R(-infinity)) || redundantLower || (!m_keepbounds))
	                  {
	                     ++lhsCnt;
	                     lhsBnd -= aij * lp.upper(j);
	
	                     lp.changeUpper(j, R(infinity));
	                     ++chgBnds;
	                  }
	                  else
	                     ++keptBnds;
	               }
	
	               if(redundantLower)
	               {
	                  // no upper bound on x_j OR redundant upper bound
	                  if((lp.upper(j) >= R(infinity)) || redundantUpper || (!m_keepbounds))
	                  {
	                     ++rhsCnt;
	                     rhsBnd -= aij * lp.lower(j);
	
	                     lp.changeLower(j, R(-infinity));
	                     ++chgBnds;
	                  }
	                  else
	                     ++keptBnds;
	               }
	            }
	         }
	      }
	
	#endif
	
	#if FREE_LHS_RHS
	
	      redundantLhs = false;
	      redundantRhs = false;
	
	      // 2. detect implied free constraints
	      if(lp.lhs(i) > R(-infinity) && lhsCnt == 0 && GErel(lhsBnd, lp.lhs(i), feastol()))
	      {
	         MSG_DEBUG((*this->spxout) << "IMAISM17 row " << i
	                   << ": redundant lhs -> lhsBnd=" << lhsBnd
	                   << " lhs=" << lp.lhs(i)
	                   << std::endl;)
	
	         redundantLhs = true;
	      }
	
	      if(lp.rhs(i) <  R(infinity) && rhsCnt == 0 && LErel(rhsBnd, lp.rhs(i), feastol()))
	      {
	         MSG_DEBUG((*this->spxout) << "IMAISM18 row " << i
	                   << ": redundant rhs -> rhsBnd=" << rhsBnd
	                   << " rhs=" << lp.rhs(i)
	                   << std::endl;)
	
	         redundantRhs = true;
	      }
	
	      if(redundantLhs)
	      {
	         // no rhs for constraint i OR redundant rhs
	         if((lp.rhs(i) >= R(infinity)) || redundantRhs || (!m_keepbounds))
	         {
	            lp.changeLhs(i, R(-infinity));
	            ++chgLRhs;
	         }
	         else
	            ++keptLRhs;
	      }
	
	      if(redundantRhs)
	      {
	         // no lhs for constraint i OR redundant lhs
	         if((lp.lhs(i) <= R(-infinity)) || redundantLhs || (!m_keepbounds))
	         {
	            lp.changeRhs(i, R(infinity));
	            ++chgLRhs;
	         }
	         else
	            ++keptLRhs;
	      }
	
	#endif
	
	      // 3. infeasible constraint
	      if(LTrel(lp.rhs(i), lp.lhs(i), feastol())                 ||
	            (LTrel(rhsBnd,   lp.lhs(i), feastol()) && rhsCnt == 0) ||
	            (GTrel(lhsBnd,   lp.rhs(i), feastol()) && lhsCnt == 0))
	      {
	         MSG_DEBUG((*this->spxout) << "IMAISM19 row " << std::setprecision(20) << i
	                   << ": infeasible -> lhs=" << lp.lhs(i)
	                   << " rhs=" << lp.rhs(i)
	                   << " lhsBnd=" << lhsBnd
	                   << " rhsBnd=" << rhsBnd
	                   << std::endl;)
	         return this->INFEASIBLE;
	      }
	
	#if FREE_CONSTRAINT
	
	      // 4. unconstrained constraint
	      if(lp.lhs(i) <= R(-infinity) && lp.rhs(i) >= R(infinity))
	      {
	         MSG_DEBUG((*this->spxout) << "IMAISM20 row " << i
	                   << ": unconstrained -> removed" << std::endl;)
	
	         std::shared_ptr<PostStep> ptr(new FreeConstraintPS(lp, i));
	         m_hist.append(ptr);
	
	         ++remRows;
	         remNzos += row.size();
	         removeRow(lp, i);
	
	         ++m_stat[FREE_ROW];
	
	         continue;
	      }
	
	#endif
	
	#if EMPTY_CONSTRAINT
	
	      // 5. empty constraint
	      if(row.size() == 0)
	      {
	         MSG_DEBUG((*this->spxout) << "IMAISM21 row " << i
	                   << ": empty ->";)
	
	         if(LT(lp.rhs(i), R(0.0), feastol()) || GT(lp.lhs(i), R(0.0), feastol()))
	         {
	            MSG_DEBUG((*this->spxout) << " infeasible lhs=" << lp.lhs(i)
	                      << " rhs=" << lp.rhs(i) << std::endl;)
	            return this->INFEASIBLE;
	         }
	
	         MSG_DEBUG((*this->spxout) << " removed" << std::endl;)
	
	         std::shared_ptr<PostStep> ptr(new EmptyConstraintPS(lp, i));
	         m_hist.append(ptr);
	
	         ++remRows;
	         removeRow(lp, i);
	
	         ++m_stat[EMPTY_ROW];
	
	         continue;
	      }
	
	#endif
	
	#if ROW_SINGLETON
	
	      // 6. row singleton
	      if(row.size() == 1)
	      {
	         removeRowSingleton(lp, row, i);
	         continue;
	      }
	
	#endif
	
	#if AGGREGATE_VARS
	
	      // 7. row doubleton, aka. simple aggregation of two variables in an equation
	      if(row.size() == 2 && EQrel(lp.lhs(i), lp.rhs(i), feastol()))
	      {
	         aggregateVars(lp, row, i);
	         continue;
	      }
	
	#endif
	
	#if FORCE_CONSTRAINT
	
	      // 8. forcing constraint (postsolving)
	      // fix variables to obtain the upper bound on constraint value
	      if(rhsCnt == 0 && EQrel(rhsBnd, lp.lhs(i), feastol()))
	      {
	         MSG_DEBUG((*this->spxout) << "IMAISM24 row " << i
	                   << ": forcing constraint fix on lhs ->"
	                   << " lhs=" << lp.lhs(i)
	                   << " rhsBnd=" << rhsBnd
	                   << std::endl;)
	
	         DataArray<bool> fixedCol(row.size());
	         Array<R> lowers(row.size());
	         Array<R> uppers(row.size());
	
	         for(int k = 0; k < row.size(); ++k)
	         {
	            R aij = row.value(k);
	            int  j   = row.index(k);
	
	            fixedCol[k] = !(EQrel(lp.upper(j), lp.lower(j), m_epsilon));
	
	            lowers[k] = lp.lower(j);
	            uppers[k] = lp.upper(j);
	
	            ASSERT_WARN("WMAISM25", isNotZero(aij, R(1.0 / R(infinity))));
	
	            if(aij > 0.0)
	               lp.changeLower(j, lp.upper(j));
	            else
	               lp.changeUpper(j, lp.lower(j));
	         }
	
	         std::shared_ptr<PostStep> ptr(new ForceConstraintPS(lp, i, true, fixedCol, lowers, uppers));
	         m_hist.append(ptr);
	
	         ++remRows;
	         remNzos += row.size();
	         removeRow(lp, i);
	
	         ++m_stat[FORCE_ROW];
	
	         continue;
	      }
	
	      // fix variables to obtain the lower bound on constraint value
	      if(lhsCnt == 0 && EQrel(lhsBnd, lp.rhs(i), feastol()))
	      {
	         MSG_DEBUG((*this->spxout) << "IMAISM26 row " << i
	                   << ": forcing constraint fix on rhs ->"
	                   << " rhs=" << lp.rhs(i)
	                   << " lhsBnd=" << lhsBnd
	                   << std::endl;)
	
	         DataArray<bool> fixedCol(row.size());
	         Array<R> lowers(row.size());
	         Array<R> uppers(row.size());
	
	         for(int k = 0; k < row.size(); ++k)
	         {
	            R aij   = row.value(k);
	            int  j     = row.index(k);
	
	            fixedCol[k] = !(EQrel(lp.upper(j), lp.lower(j), m_epsilon));
	
	            lowers[k] = lp.lower(j);
	            uppers[k] = lp.upper(j);
	
	            ASSERT_WARN("WMAISM27", isNotZero(aij, R(1.0 / R(infinity))));
	
	            if(aij > 0.0)
	               lp.changeUpper(j, lp.lower(j));
	            else
	               lp.changeLower(j, lp.upper(j));
	         }
	
	         std::shared_ptr<PostStep> ptr(new ForceConstraintPS(lp, i, false, fixedCol, lowers, uppers));
	         m_hist.append(ptr);
	
	         ++remRows;
	         remNzos += row.size();
	         removeRow(lp, i);
	
	         ++m_stat[FORCE_ROW];
	
	         continue;
	      }
	
	#endif
	   }
	
	   assert(remRows > 0 || remNzos == 0);
	
	   if(remRows + chgLRhs + chgBnds > 0)
	   {
	      this->m_remRows += remRows;
	      this->m_remNzos += remNzos;
	      this->m_chgLRhs += chgLRhs;
	      this->m_chgBnds += chgBnds;
	      this->m_keptBnds += keptBnds;
	      this->m_keptLRhs += keptLRhs;
	
	      MSG_INFO2((*this->spxout), (*this->spxout) << "Simplifier (rows) removed "
	                << remRows << " rows, "
	                << remNzos << " non-zeros, "
	                << chgBnds << " col bounds, "
	                << chgLRhs << " row bounds; kept "
	                << keptBnds << " column bounds, "
	                << keptLRhs << " row bounds"
	                << std::endl;)
	
	      if(remRows > this->m_minReduction * oldRows)
	         again = true;
	   }
	
	   return this->OKAY;
	}
	
	template <class R>
	typename SPxSimplifier<R>::Result SPxMainSM<R>::simplifyCols(SPxLPBase<R>& lp, bool& again)
	{
	
	   // This method simplifies the columns of the LP.
	   //
	   // The following operations are done:
	   // 1. detect empty columns and fix corresponding variables
	   // 2. detect variables that are unconstrained from below or above
	   //    and fix corresponding variables or remove involved constraints
	   // 3. fix variables
	   // 4. use column singleton variables with zero objective to adjust constraint bounds
	   // 5. (not free) column singleton combined with doubleton equation are
	   //    used to make the column singleton variable free
	   // 6. substitute (implied) free column singletons
	
	   int remRows = 0;
	   int remCols = 0;
	   int remNzos = 0;
	   int chgBnds = 0;
	
	   int oldCols = lp.nCols();
	   int oldRows = lp.nRows();
	
	   for(int j = lp.nCols() - 1; j >= 0; --j)
	   {
	      const SVectorBase<R>& col = lp.colVector(j);
	
	      // infeasible bounds
	      if(GTrel(lp.lower(j), lp.upper(j), feastol()))
	      {
	         MSG_DEBUG((*this->spxout) << "IMAISM29 col " << j
	                   << ": infeasible bounds on x" << j
	                   << " -> lower=" << lp.lower(j)
	                   << " upper=" << lp.upper(j)
	                   << std::endl;)
	         return this->INFEASIBLE;
	      }
	
	      // 1. empty column
	      if(col.size() == 0)
	      {
	#if EMPTY_COLUMN
	         MSG_DEBUG((*this->spxout) << "IMAISM30 col " << j
	                   << ": empty -> maxObj=" << lp.maxObj(j)
	                   << " lower=" << lp.lower(j)
	                   << " upper=" << lp.upper(j);)
	
	         R val;
	
	         if(GT(lp.maxObj(j), R(0.0), this->epsZero()))
	         {
	            if(lp.upper(j) >= R(infinity))
	            {
	               MSG_DEBUG((*this->spxout) << " unbounded" << std::endl;)
	               return this->UNBOUNDED;
	            }
	
	            val = lp.upper(j);
	         }
	         else if(LT(lp.maxObj(j), R(0.0), this->epsZero()))
	         {
	            if(lp.lower(j) <= R(-infinity))
	            {
	               MSG_DEBUG((*this->spxout) << " unbounded" << std::endl;)
	               return this->UNBOUNDED;
	            }
	
	            val = lp.lower(j);
	         }
	         else
	         {
	            assert(isZero(lp.maxObj(j), this->epsZero()));
	
	            // any value within the bounds is ok
	            if(lp.lower(j) > R(-infinity))
	               val = lp.lower(j);
	            else if(lp.upper(j) < R(infinity))
	               val = lp.upper(j);
	            else
	               val = 0.0;
	         }
	
	         MSG_DEBUG((*this->spxout) << " removed" << std::endl;)
	
	         std::shared_ptr<PostStep> ptr1(new FixBoundsPS(lp, j, val));
	         std::shared_ptr<PostStep> ptr2(new FixVariablePS(lp, *this, j, val));
	         m_hist.append(ptr1);
	         m_hist.append(ptr2);
	
	         ++remCols;
	         removeCol(lp, j);
	
	         ++m_stat[EMPTY_COL];
	
	         continue;
	#endif
	      }
	
	      if(NErel(lp.lower(j), lp.upper(j), feastol()))
	      {
	         // will be set to false if any constraint implies a bound on the variable
	         bool loFree = true;
	         bool upFree = true;
	
	         // 1. fix and remove variables
	         for(int k = 0; k < col.size(); ++k)
	         {
	            if(!loFree && !upFree)
	               break;
	
	            int i = col.index(k);
	
	            // warn since this unhandled case may slip through unnoticed otherwise
	            ASSERT_WARN("WMAISM31", isNotZero(col.value(k), R(1.0 / R(infinity))));
	
	            if(col.value(k) > 0.0)
	            {
	               if(lp.rhs(i) <  R(infinity))
	                  upFree = false;
	
	               if(lp.lhs(i) > R(-infinity))
	                  loFree = false;
	            }
	            else if(col.value(k) < 0.0)
	            {
	               if(lp.rhs(i) <  R(infinity))
	                  loFree = false;
	
	               if(lp.lhs(i) > R(-infinity))
	                  upFree = false;
	            }
	         }
	
	         // 2. detect variables that are unconstrained from below or above
	         // max  3 x
	         // s.t. 5 x >= 8
	         if(GT(lp.maxObj(j), R(0.0), this->epsZero()) && upFree)
	         {
	#if FIX_VARIABLE
	            MSG_DEBUG((*this->spxout) << "IMAISM32 col " << j
	                      << ": x" << j
	                      << " unconstrained above ->";)
	
	            if(lp.upper(j) >= R(infinity))
	            {
	               MSG_DEBUG((*this->spxout) << " unbounded" << std::endl;)
	
	               return this->UNBOUNDED;
	            }
	
	            MSG_DEBUG((*this->spxout) << " fixed at upper=" << lp.upper(j) << std::endl;)
	
	            std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j, lp.upper(j)));
	            m_hist.append(ptr);
	            lp.changeLower(j, lp.upper(j));
	         }
	         // max -3 x
	         // s.t. 5 x <= 8
	         else if(LT(lp.maxObj(j), R(0.0), this->epsZero()) && loFree)
	         {
	            MSG_DEBUG((*this->spxout) << "IMAISM33 col " << j
	                      << ": x" << j
	                      << " unconstrained below ->";)
	
	            if(lp.lower(j) <= R(-infinity))
	            {
	               MSG_DEBUG((*this->spxout) << " unbounded" << std::endl;)
	
	               return this->UNBOUNDED;
	            }
	
	            MSG_DEBUG((*this->spxout) << " fixed at lower=" << lp.lower(j) << std::endl;)
	
	            std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j, lp.lower(j)));
	            m_hist.append(ptr);
	            lp.changeUpper(j, lp.lower(j));
	#endif
	         }
	         else if(isZero(lp.maxObj(j), this->epsZero()))
	         {
	#if FREE_ZERO_OBJ_VARIABLE
	            bool unconstrained_below = loFree && lp.lower(j) <= R(-infinity);
	            bool unconstrained_above = upFree && lp.upper(j) >= R(infinity);
	
	            if(unconstrained_below || unconstrained_above)
	            {
	               MSG_DEBUG((*this->spxout) << "IMAISM34 col " << j
	                         << ": x" << j
	                         << " unconstrained "
	                         << (unconstrained_below ? "below" : "above")
	                         << " with zero objective (" << lp.maxObj(j)
	                         << ")" << std::endl;)
	
	               DSVectorBase<R> col_idx_sorted(col);
	
	               // sort col elements by increasing idx
	               IdxCompare compare;
	               SPxQuicksort(col_idx_sorted.mem(), col_idx_sorted.size(), compare);
	
	               std::shared_ptr<PostStep> ptr(new FreeZeroObjVariablePS(lp, j, unconstrained_below,
	                                             col_idx_sorted));
	               m_hist.append(ptr);
	
	               // we have to remove the rows with larger idx first, because otherwise the rows are reorder and indices
	               // are out-of-date
	               remRows += col.size();
	
	               for(int k = col_idx_sorted.size() - 1; k >= 0; --k)
	                  removeRow(lp, col_idx_sorted.index(k));
	
	               // remove column
	               removeCol(lp, j);
	
	               // statistics
	               for(int k = 0; k < col.size(); ++k)
	               {
	                  int l   =  col.index(k);
	                  remNzos += lp.rowVector(l).size();
	               }
	
	               ++m_stat[FREE_ZOBJ_COL];
	               ++remCols;
	
	               continue;
	            }
	
	#endif
	         }
	      }
	
	#if FIX_VARIABLE
	
	      // 3. fix variable
	      if(EQrel(lp.lower(j), lp.upper(j), feastol()))
	      {
	         MSG_DEBUG((*this->spxout) << "IMAISM36 col " << j
	                   << ": x" << j
	                   << " fixed -> lower=" << lp.lower(j)
	                   << " upper=" << lp.upper(j) << std::endl;)
	
	         fixColumn(lp, j);
	
	         ++remCols;
	         remNzos += col.size();
	         removeCol(lp, j);
	
	         ++m_stat[FIX_COL];
	
	         continue;
	      }
	
	#endif
	
	      // handle column singletons
	      if(col.size() == 1)
	      {
	         R aij = col.value(0);
	         int  i   = col.index(0);
	
	         // 4. column singleton with zero objective
	         if(isZero(lp.maxObj(j), this->epsZero()))
	         {
	#if ZERO_OBJ_COL_SINGLETON
	            MSG_DEBUG((*this->spxout) << "IMAISM37 col " << j
	                      << ": singleton in row " << i
	                      << " with zero objective";)
	
	            R lhs = R(-infinity);
	            R rhs = +R(infinity);
	
	            if(GT(aij, R(0.0), this->epsZero()))
	            {
	               if(lp.lhs(i) > R(-infinity) && lp.upper(j) <  R(infinity))
	                  lhs = lp.lhs(i) - aij * lp.upper(j);
	
	               if(lp.rhs(i) <  R(infinity) && lp.lower(j) > R(-infinity))
	                  rhs = lp.rhs(i) - aij * lp.lower(j);
	            }
	            else if(LT(aij, R(0.0), this->epsZero()))
	            {
	               if(lp.lhs(i) > R(-infinity) && lp.lower(j) > R(-infinity))
	                  lhs = lp.lhs(i) - aij * lp.lower(j);
	
	               if(lp.rhs(i) <  R(infinity) && lp.upper(j) <  R(infinity))
	                  rhs = lp.rhs(i) - aij * lp.upper(j);
	            }
	            else
	            {
	               lhs = lp.lhs(i);
	               rhs = lp.rhs(i);
	            }
	
	            if(isZero(lhs, this->epsZero()))
	               lhs = 0.0;
	
	            if(isZero(rhs, this->epsZero()))
	               rhs = 0.0;
	
	            MSG_DEBUG((*this->spxout) << " removed -> lhs=" << lhs
	                      << " (" << lp.lhs(i)
	                      << ") rhs=" << rhs
	                      << " (" << lp.rhs(i)
	                      << ")" << std::endl;)
	
	            std::shared_ptr<PostStep> ptr(new ZeroObjColSingletonPS(lp, *this, j, i));
	            m_hist.append(ptr);
	
	            lp.changeRange(i, lhs, rhs);
	
	            ++remCols;
	            ++remNzos;
	            removeCol(lp, j);
	
	            ++m_stat[ZOBJ_SINGLETON_COL];
	
	            if(lp.lhs(i) <= R(-infinity) && lp.rhs(i) >= R(infinity))
	            {
	               std::shared_ptr<PostStep> ptr2(new FreeConstraintPS(lp, i));
	               m_hist.append(ptr2);
	
	               ++remRows;
	               removeRow(lp, i);
	
	               ++m_stat[FREE_ROW];
	            }
	
	            continue;
	#endif
	         }
	
	         // 5. not free column singleton combined with doubleton equation
	         else if(EQrel(lp.lhs(i), lp.rhs(i), feastol())             &&
	                 lp.rowVector(i).size() == 2                         &&
	                 (lp.lower(j) > R(-infinity) || lp.upper(j) < R(infinity)))
	         {
	#if DOUBLETON_EQUATION
	            MSG_DEBUG((*this->spxout) << "IMAISM38 col " << j
	                      << ": singleton in row " << i
	                      << " with doubleton equation ->";)
	
	            R lhs = lp.lhs(i);
	
	            const SVectorBase<R>& row = lp.rowVector(i);
	
	            R aik;
	            int  k;
	
	            if(row.index(0) == j)
	            {
	               aik = row.value(1);
	               k   = row.index(1);
	            }
	            else if(row.index(1) == j)
	            {
	               aik = row.value(0);
	               k   = row.index(0);
	            }
	            else
	               throw SPxInternalCodeException("XMAISM11 This should never happen.");
	
	            ASSERT_WARN("WMAISM39", isNotZero(aik, R(1.0 / R(infinity))));
	
	            R lo, up;
	            R oldLower = lp.lower(k);
	            R oldUpper = lp.upper(k);
	
	            R scale1 = maxAbs(lhs, aij * lp.upper(j));
	            R scale2 = maxAbs(lhs, aij * lp.lower(j));
	
	            if(scale1 < 1.0)
	               scale1 = 1.0;
	
	            if(scale2 < 1.0)
	               scale2 = 1.0;
	
	            R z1 = (lhs / scale1) - (aij * lp.upper(j) / scale1);
	            R z2 = (lhs / scale2) - (aij * lp.lower(j) / scale2);
	
	            if(isZero(z1, this->epsZero()))
	               z1 = 0.0;
	
	            if(isZero(z2, this->epsZero()))
	               z2 = 0.0;
	
	            if(aij * aik > 0.0)
	            {
	               lo = (lp.upper(j) >=  R(infinity)) ? R(-infinity) : z1 * scale1 / aik;
	               up = (lp.lower(j) <= R(-infinity)) ?  R(infinity) : z2 * scale2 / aik;
	            }
	            else if(aij * aik < 0.0)
	            {
	               lo = (lp.lower(j) <= R(-infinity)) ? R(-infinity) : z2 * scale2 / aik;
	               up = (lp.upper(j) >=  R(infinity)) ?  R(infinity) : z1 * scale1 / aik;
	            }
	            else
	               throw SPxInternalCodeException("XMAISM12 This should never happen.");
	
	            if(GTrel(lo, lp.lower(k), this->epsZero()))
	               lp.changeLower(k, lo);
	
	            if(LTrel(up, lp.upper(k), this->epsZero()))
	               lp.changeUpper(k, up);
	
	            MSG_DEBUG((*this->spxout) << " made free, bounds on x" << k
	                      << ": lower=" << lp.lower(k)
	                      << " (" << oldLower
	                      << ") upper=" << lp.upper(k)
	                      << " (" << oldUpper
	                      << ")" << std::endl;)
	
	            // infeasible bounds
	            if(GTrel(lp.lower(k), lp.upper(k), feastol()))
	            {
	               MSG_DEBUG((*this->spxout) << "new bounds are infeasible"
	                         << std::endl;)
	               return this->INFEASIBLE;
	            }
	
	            std::shared_ptr<PostStep> ptr(new DoubletonEquationPS(lp, j, k, i, oldLower, oldUpper));
	            m_hist.append(ptr);
	
	            if(lp.lower(j) > R(-infinity) && lp.upper(j) < R(infinity))
	               chgBnds += 2;
	            else
	               ++chgBnds;
	
	            lp.changeBounds(j, R(-infinity), R(infinity));
	
	            ++m_stat[DOUBLETON_ROW];
	#endif
	         }
	
	         // 6. (implied) free column singleton
	         if(lp.lower(j) <= R(-infinity) && lp.upper(j) >= R(infinity))
	         {
	#if FREE_COL_SINGLETON
	            R slackVal = lp.lhs(i);
	
	            // constraint i is an inequality constraint -> transform into equation type
	            if(NErel(lp.lhs(i), lp.rhs(i), feastol()))
	            {
	               MSG_DEBUG((*this->spxout) << "IMAISM40 col " << j
	                         << ": free singleton in inequality constraint" << std::endl;)
	
	               // do nothing if constraint i is unconstrained
	               if(lp.lhs(i) <= R(-infinity) && lp.rhs(i) >= R(infinity))
	                  continue;
	
	               // introduce slack variable to obtain equality constraint
	               R sMaxObj = lp.maxObj(j) / aij; // after substituting variable j in objective
	               R sLo     = lp.lhs(i);
	               R sUp     = lp.rhs(i);
	
	               if(GT(sMaxObj, R(0.0), this->epsZero()))
	               {
	                  if(sUp >= R(infinity))
	                  {
	                     MSG_DEBUG((*this->spxout) << " -> problem unbounded" << std::endl;)
	                     return this->UNBOUNDED;
	                  }
	
	                  slackVal = sUp;
	               }
	               else if(LT(sMaxObj, R(0.0), this->epsZero()))
	               {
	                  if(sLo <= R(-infinity))
	                  {
	                     MSG_DEBUG((*this->spxout) << " -> problem unbounded" << std::endl;)
	                     return this->UNBOUNDED;
	                  }
	
	                  slackVal = sLo;
	               }
	               else
	               {
	                  assert(isZero(sMaxObj, this->epsZero()));
	
	                  // any value within the bounds is ok
	                  if(sLo > R(-infinity))
	                     slackVal = sLo;
	                  else if(sUp < R(infinity))
	                     slackVal = sUp;
	                  else
	                     throw SPxInternalCodeException("XMAISM13 This should never happen.");
	               }
	            }
	
	            std::shared_ptr<PostStep> ptr(new FreeColSingletonPS(lp, *this, j, i, slackVal));
	            m_hist.append(ptr);
	
	            MSG_DEBUG((*this->spxout) << "IMAISM41 col " << j
	                      << ": free singleton removed" << std::endl;)
	
	            const SVectorBase<R>& row = lp.rowVector(i);
	
	            for(int h = 0; h < row.size(); ++h)
	            {
	               int k = row.index(h);
	
	               if(k != j)
	               {
	                  R new_obj = lp.obj(k) - (lp.obj(j) * row.value(h) / aij);
	                  lp.changeObj(k, new_obj);
	               }
	            }
	
	            ++remRows;
	            ++remCols;
	            remNzos += row.size();
	            removeRow(lp, i);
	            removeCol(lp, j);
	
	            ++m_stat[FREE_SINGLETON_COL];
	
	            continue;
	#endif
	         }
	      }
	   }
	
	   if(remCols + remRows > 0)
	   {
	      this->m_remRows += remRows;
	      this->m_remCols += remCols;
	      this->m_remNzos += remNzos;
	      this->m_chgBnds += chgBnds;
	
	      MSG_INFO2((*this->spxout), (*this->spxout) << "Simplifier (columns) removed "
	                << remRows << " rows, "
	                << remCols << " cols, "
	                << remNzos << " non-zeros, "
	                << chgBnds << " col bounds"
	                << std::endl;)
	
	      if(remCols + remRows > this->m_minReduction * (oldCols + oldRows))
	         again = true;
	   }
	
	   return this->OKAY;
	}
	
	template <class R>
	typename SPxSimplifier<R>::Result SPxMainSM<R>::simplifyDual(SPxLPBase<R>& lp, bool& again)
	{
	
	   // This method simplifies LP using the following dual structures:
	   //
	   // 1. dominated columns
	   // 2. weakly dominated columns
	   //
	   // For constructing the dual variables, it is assumed that the objective sense is max
	
	   int remRows = 0;
	   int remCols = 0;
	   int remNzos = 0;
	
	   int oldRows = lp.nRows();
	   int oldCols = lp.nCols();
	
	   DataArray<bool> colSingleton(lp.nCols());
	   VectorBase<R>         dualVarLo(lp.nRows());
	   VectorBase<R>         dualVarUp(lp.nRows());
	   VectorBase<R>         dualConsLo(lp.nCols());
	   VectorBase<R>         dualConsUp(lp.nCols());
	
	   // init
	   for(int i = lp.nRows() - 1; i >= 0; --i)
	   {
	      // check for unconstrained constraints
	      if(lp.lhs(i) <= R(-infinity) && lp.rhs(i) >= R(infinity))
	      {
	         MSG_DEBUG((*this->spxout) << "IMAISM43 row " << i
	                   << ": unconstrained" << std::endl;)
	
	         std::shared_ptr<PostStep> ptr(new FreeConstraintPS(lp, i));
	         m_hist.append(ptr);
	
	         ++remRows;
	         remNzos += lp.rowVector(i).size();
	         removeRow(lp, i);
	
	         ++m_stat[FREE_ROW];
	
	         continue;
	      }
	
	      // corresponds to maximization sense
	      dualVarLo[i] = (lp.lhs(i) <= R(-infinity)) ? 0.0 : R(-infinity);
	      dualVarUp[i] = (lp.rhs(i) >=  R(infinity)) ? 0.0 :  R(infinity);
	   }
	
	   // compute bounds on the dual variables using column singletons
	   for(int j = 0; j < lp.nCols(); ++j)
	   {
	      if(lp.colVector(j).size() == 1)
	      {
	         int  i   = lp.colVector(j).index(0);
	         R aij = lp.colVector(j).value(0);
	
	         ASSERT_WARN("WMAISM44", isNotZero(aij, R(1.0 / R(infinity))));
	
	         R bound = lp.maxObj(j) / aij;
	
	         if(aij > 0)
	         {
	            if(lp.lower(j) <= R(-infinity) && bound < dualVarUp[i])
	               dualVarUp[i] = bound;
	
	            if(lp.upper(j) >=  R(infinity) && bound > dualVarLo[i])
	               dualVarLo[i] = bound;
	         }
	         else if(aij < 0)
	         {
	            if(lp.lower(j) <= R(-infinity) && bound > dualVarLo[i])
	               dualVarLo[i] = bound;
	
	            if(lp.upper(j) >=  R(infinity) && bound < dualVarUp[i])
	               dualVarUp[i] = bound;
	         }
	      }
	
	   }
	
	   // compute bounds on the dual constraints
	   for(int j = 0; j < lp.nCols(); ++j)
	   {
	      dualConsLo[j] = dualConsUp[j] = 0.0;
	
	      const SVectorBase<R>& col = lp.colVector(j);
	
	      for(int k = 0; k < col.size(); ++k)
	      {
	         if(dualConsLo[j] <= R(-infinity) && dualConsUp[j] >= R(infinity))
	            break;
	
	         R aij = col.value(k);
	         int  i   = col.index(k);
	
	         ASSERT_WARN("WMAISM45", isNotZero(aij, R(1.0 / R(infinity))));
	
	         if(aij > 0)
	         {
	            if(dualVarLo[i] <= R(-infinity))
	               dualConsLo[j] = R(-infinity);
	            else
	               dualConsLo[j] += aij * dualVarLo[i];
	
	            if(dualVarUp[i] >= R(infinity))
	               dualConsUp[j] = R(infinity);
	            else
	               dualConsUp[j] += aij * dualVarUp[i];
	         }
	         else if(aij < 0)
	         {
	            if(dualVarLo[i] <= R(-infinity))
	               dualConsUp[j] = R(infinity);
	            else
	               dualConsUp[j] += aij * dualVarLo[i];
	
	            if(dualVarUp[i] >= R(infinity))
	               dualConsLo[j] = R(-infinity);
	            else
	               dualConsLo[j] += aij * dualVarUp[i];
	         }
	      }
	   }
	
	   for(int j = lp.nCols() - 1; j >= 0; --j)
	   {
	      if(lp.colVector(j).size() <= 1)
	         continue;
	
	      // dual infeasibility checks
	      if(LTrel(dualConsUp[j], dualConsLo[j], opttol()))
	      {
	         MSG_DEBUG((*this->spxout) << "IMAISM46 col " << j
	                   << ": dual infeasible -> dual lhs bound=" << dualConsLo[j]
	                   << " dual rhs bound=" << dualConsUp[j] << std::endl;)
	         return this->DUAL_INFEASIBLE;
	      }
	
	      R obj = lp.maxObj(j);
	
	      // 1. dominated column
	      // Is the problem really unbounded in the cases below ??? Or is only dual infeasibility be shown
	      if(GTrel(obj, dualConsUp[j], opttol()))
	      {
	#if DOMINATED_COLUMN
	         MSG_DEBUG((*this->spxout) << "IMAISM47 col " << j
	                   << ": dominated -> maxObj=" << obj
	                   << " dual rhs bound=" << dualConsUp[j] << std::endl;)
	
	         if(lp.upper(j) >= R(infinity))
	         {
	            MSG_INFO2((*this->spxout), (*this->spxout) << " unbounded" << std::endl;)
	            return this->UNBOUNDED;
	         }
	
	         MSG_DEBUG((*this->spxout) << " fixed at upper=" << lp.upper(j) << std::endl;)
	
	         std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j, lp.upper(j)));
	         m_hist.append(ptr);
	         lp.changeLower(j, lp.upper(j));
	
	         ++m_stat[DOMINATED_COL];
	#endif
	      }
	      else if(LTrel(obj, dualConsLo[j], opttol()))
	      {
	#if DOMINATED_COLUMN
	         MSG_DEBUG((*this->spxout) << "IMAISM48 col " << j
	                   << ": dominated -> maxObj=" << obj
	                   << " dual lhs bound=" << dualConsLo[j] << std::endl;)
	
	         if(lp.lower(j) <= R(-infinity))
	         {
	            MSG_INFO2((*this->spxout), (*this->spxout) << " unbounded" << std::endl;)
	            return this->UNBOUNDED;
	         }
	
	         MSG_DEBUG((*this->spxout) << " fixed at lower=" << lp.lower(j) << std::endl;)
	
	         std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j, lp.lower(j)));
	         m_hist.append(ptr);
	         lp.changeUpper(j, lp.lower(j));
	
	         ++m_stat[DOMINATED_COL];
	#endif
	      }
	
	      // 2. weakly dominated column (no postsolving)
	      else if(lp.upper(j) < R(infinity) && EQrel(obj, dualConsUp[j], opttol()))
	      {
	#if WEAKLY_DOMINATED_COLUMN
	         MSG_DEBUG((*this->spxout) << "IMAISM49 col " << j
	                   << ": weakly dominated -> maxObj=" << obj
	                   << " dual rhs bound=" << dualConsUp[j] << std::endl;)
	
	         std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j, lp.upper(j)));
	         m_hist.append(ptr);
	         lp.changeLower(j, lp.upper(j));
	
	         ++m_stat[WEAKLY_DOMINATED_COL];
	#endif
	      }
	      else if(lp.lower(j) > R(-infinity) && EQrel(obj, dualConsLo[j], opttol()))
	      {
	#if WEAKLY_DOMINATED_COLUMN
	         MSG_DEBUG((*this->spxout) << "IMAISM50 col " << j
	                   << ": weakly dominated -> maxObj=" << obj
	                   << " dual lhs bound=" << dualConsLo[j] << std::endl;)
	
	         std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j, lp.lower(j)));
	         m_hist.append(ptr);
	         lp.changeUpper(j, lp.lower(j));
	
	         ++m_stat[WEAKLY_DOMINATED_COL];
	#endif
	      }
	
	      // fix column
	      if(EQrel(lp.lower(j), lp.upper(j), feastol()))
	      {
	#if FIX_VARIABLE
	         fixColumn(lp, j);
	
	         ++remCols;
	         remNzos += lp.colVector(j).size();
	         removeCol(lp, j);
	
	         ++m_stat[FIX_COL];
	#endif
	      }
	   }
	
	
	   assert(remRows > 0 || remCols > 0 || remNzos == 0);
	
	   if(remCols + remRows > 0)
	   {
	      this->m_remRows += remRows;
	      this->m_remCols += remCols;
	      this->m_remNzos += remNzos;
	
	      MSG_INFO2((*this->spxout), (*this->spxout) << "Simplifier (dual) removed "
	                << remRows << " rows, "
	                << remCols << " cols, "
	                << remNzos << " non-zeros"
	                << std::endl;)
	
	      if(remCols + remRows > this->m_minReduction * (oldCols + oldRows))
	         again = true;
	   }
	
	   return this->OKAY;
	}
	
	
	
	template <class R>
	typename SPxSimplifier<R>::Result SPxMainSM<R>::multiaggregation(SPxLPBase<R>& lp, bool& again)
	{
	   // this simplifier eliminates rows and columns by performing multi aggregations as identified by the constraint
	   // activities.
	   int remRows = 0;
	   int remCols = 0;
	   int remNzos = 0;
	
	   int oldRows = lp.nRows();
	   int oldCols = lp.nCols();
	
	   VectorBase<R> upLocks(lp.nCols());
	   VectorBase<R> downLocks(lp.nCols());
	
	   for(int j = lp.nCols() - 1; j >= 0; --j)
	   {
	      // setting the locks on the variables
	      upLocks[j] = 0;
	      downLocks[j] = 0;
	
	      if(lp.colVector(j).size() <= 1)
	         continue;
	
	      const SVectorBase<R>& col = lp.colVector(j);
	
	      for(int k = 0; k < col.size(); ++k)
	      {
	         R aij = col.value(k);
	
	         ASSERT_WARN("WMAISM45", isNotZero(aij, R(1.0 / R(infinity))));
	
	         if(GT(lp.lhs(col.index(k)), R(-infinity)) && LT(lp.rhs(col.index(k)), R(infinity)))
	         {
	            upLocks[j]++;
	            downLocks[j]++;
	         }
	         else if(GT(lp.lhs(col.index(k)), R(-infinity)))
	         {
	            if(aij > 0)
	               downLocks[j]++;
	            else if(aij < 0)
	               upLocks[j]++;
	         }
	         else if(LT(lp.rhs(col.index(k)), R(infinity)))
	         {
	            if(aij > 0)
	               upLocks[j]++;
	            else if(aij < 0)
	               downLocks[j]++;
	         }
	      }
	
	      // multi-aggregate column
	      if(upLocks[j] == 1 || downLocks[j] == 1)
	      {
	         R lower = lp.lower(j);
	         R upper = lp.upper(j);
	         int maxOtherLocks;
	         int bestpos = -1;
	         bool bestislhs = true;
	
	         for(int k = 0; k < col.size(); ++k)
	         {
	            int rowNumber;
	            R lhs;
	            R rhs;
	            bool lhsExists;
	            bool rhsExists;
	            bool aggLhs;
	            bool aggRhs;
	
	            R val = col.value(k);
	
	            rowNumber = col.index(k);
	            lhs = lp.lhs(rowNumber);
	            rhs = lp.rhs(rowNumber);
	
	            if(EQ(lhs, rhs, feastol()))
	               continue;
	
	            lhsExists = GT(lhs, R(-infinity));
	            rhsExists = LT(rhs, R(infinity));
	
	            if(lp.rowVector(rowNumber).size() <= 2)
	               maxOtherLocks = INT_MAX;
	            else if(lp.rowVector(rowNumber).size() == 3)
	               maxOtherLocks = 3;
	            else if(lp.rowVector(rowNumber).size() == 4)
	               maxOtherLocks = 2;
	            else
	               maxOtherLocks = 1;
	
	            aggLhs = lhsExists
	                     && ((col.value(k) > 0.0 && lp.maxObj(j) <= 0.0 && downLocks[j] == 1 && upLocks[j] <= maxOtherLocks)
	                         || (col.value(k) < 0.0 && lp.maxObj(j) >= 0.0 && upLocks[j] == 1 && downLocks[j] <= maxOtherLocks));
	            aggRhs = rhsExists
	                     && ((col.value(k) > 0.0 && lp.maxObj(j) >= 0.0 && upLocks[j] == 1 && downLocks[j] <= maxOtherLocks)
	                         || (col.value(k) < 0.0 && lp.maxObj(j) <= 0.0 && downLocks[j] == 1 && upLocks[j] <= maxOtherLocks));
	
	            if(aggLhs || aggRhs)
	            {
	               R minRes = 0;   // this is the minimum value that the aggregation can attain
	               R maxRes = 0;   // this is the maximum value that the aggregation can attain
	
	               // computing the minimum and maximum residuals if variable j is set to zero.
	               computeMinMaxResidualActivity(lp, rowNumber, j, minRes, maxRes);
	
	               // we will try to aggregate to the lhs
	               if(aggLhs)
	               {
	                  R minVal;
	                  R maxVal;
	
	                  // computing the values of the upper and lower bounds for the aggregated variables
	                  computeMinMaxValues(lp, lhs, val, minRes, maxRes, minVal, maxVal);
	
	                  assert(LE(minVal, maxVal));
	
	                  // if the bounds of the aggregation and the original variable are equivalent, then we can reduce
	                  if((minVal > R(-infinity) && GT(minVal, lower, feastol()))
	                        && (maxVal < R(infinity) && LT(maxVal, upper, feastol())))
	                  {
	                     bestpos = col.index(k);
	                     bestislhs = true;
	                     break;
	                  }
	               }
	
	               // we will try to aggregate to the rhs
	               if(aggRhs)
	               {
	                  R minVal;
	                  R maxVal;
	
	                  // computing the values of the upper and lower bounds for the aggregated variables
	                  computeMinMaxValues(lp, rhs, val, minRes, maxRes, minVal, maxVal);
	
	                  assert(LE(minVal, maxVal));
	
	                  if((minVal > R(-infinity) && GT(minVal, lower, feastol()))
	                        && (maxVal < R(infinity) && LT(maxVal, upper, feastol())))
	                  {
	                     bestpos = col.index(k);
	                     bestislhs = false;
	                     break;
	                  }
	               }
	            }
	         }
	
	         // it is only possible to aggregate if a best position has been found
	         if(bestpos >= 0)
	         {
	            const SVectorBase<R>& bestRow = lp.rowVector(bestpos);
	            // aggregating the variable and applying the fixings to the all other constraints
	            R aggConstant = (bestislhs ? lp.lhs(bestpos) : lp.rhs(
	                                bestpos));   // this is the lhs or rhs of the aggregated row
	            R aggAij =
	               bestRow[j];                                   // this is the coefficient of the deleted col
	
	            MSG_DEBUG(
	               (*this->spxout) << "IMAISM51 col " << j
	               << ": Aggregating row: " << bestpos
	               << " Aggregation Constant=" << aggConstant
	               << " Coefficient of aggregated col=" << aggAij << std::endl;
	            )
	
	            std::shared_ptr<PostStep> ptr(new MultiAggregationPS(lp, *this, bestpos, j, aggConstant));
	            m_hist.append(ptr);
	
	            for(int k = 0; k < col.size(); ++k)
	            {
	               if(col.index(k) != bestpos)
	               {
	                  int rowNumber = col.index(k);
	                  VectorBase<R> updateRow(lp.nCols());
	                  R updateRhs = lp.rhs(col.index(k));
	                  R updateLhs = lp.lhs(col.index(k));
	
	                  updateRow = lp.rowVector(col.index(k));
	
	                  // updating the row with the best row
	                  for(int l = 0; l < bestRow.size(); l++)
	                  {
	                     if(bestRow.index(l) != j)
	                     {
	                        if(lp.rowVector(rowNumber).pos(bestRow.index(l)) >= 0)
	                           lp.changeElement(rowNumber, bestRow.index(l), updateRow[bestRow.index(l)]
	                                            - updateRow[j]*bestRow.value(l) / aggAij);
	                        else
	                           lp.changeElement(rowNumber, bestRow.index(l), -1.0 * updateRow[j]*bestRow.value(l) / aggAij);
	                     }
	                  }
	
	                  // NOTE: I don't know whether we should change the LHS and RHS if they are currently at R(infinity)
	                  if(LT(lp.rhs(rowNumber), R(infinity)))
	                     lp.changeRhs(rowNumber, updateRhs - updateRow[j]*aggConstant / aggAij);
	
	                  if(GT(lp.lhs(rowNumber), R(-infinity)))
	                     lp.changeLhs(rowNumber, updateLhs - updateRow[j]*aggConstant / aggAij);
	
	                  assert(LE(lp.lhs(rowNumber), lp.rhs(rowNumber)));
	               }
	            }
	
	            for(int l = 0; l < bestRow.size(); l++)
	            {
	               if(bestRow.index(l) != j)
	                  lp.changeMaxObj(bestRow.index(l),
	                                  lp.maxObj(bestRow.index(l)) - lp.maxObj(j)*bestRow.value(l) / aggAij);
	            }
	
	            ++remCols;
	            remNzos += lp.colVector(j).size();
	            removeCol(lp, j);
	            ++remRows;
	            remNzos += lp.rowVector(bestpos).size();
	            removeRow(lp, bestpos);
	
	            ++m_stat[MULTI_AGG];
	         }
	      }
	   }
	
	
	   assert(remRows > 0 || remCols > 0 || remNzos == 0);
	
	   if(remCols + remRows > 0)
	   {
	      this->m_remRows += remRows;
	      this->m_remCols += remCols;
	      this->m_remNzos += remNzos;
	
	      MSG_INFO2((*this->spxout), (*this->spxout) << "Simplifier (multi-aggregation) removed "
	                << remRows << " rows, "
	                << remCols << " cols, "
	                << remNzos << " non-zeros"
	                << std::endl;)
	
	      if(remCols + remRows > this->m_minReduction * (oldCols + oldRows))
	         again = true;
	   }
	
	   return this->OKAY;
	}
	
	
	
	template <class R>
	typename SPxSimplifier<R>::Result SPxMainSM<R>::duplicateRows(SPxLPBase<R>& lp, bool& again)
	{
	
	   // This method simplifies the LP by removing duplicate rows
	   // Duplicates are detected using the algorithm of Bixby and Wagner [1987]
	
	   // Possible extension: use generalized definition of duplicate rows according to Andersen and Andersen
	   // However: the resulting sparsification is often very small since the involved rows are usually very sparse
	
	   int remRows = 0;
	   int remNzos = 0;
	
	   int oldRows = lp.nRows();
	
	   // remove empty rows and columns
	   typename SPxSimplifier<R>::Result ret = removeEmpty(lp);
	
	   if(ret != this->OKAY)
	      return ret;
	
	#if ROW_SINGLETON
	   int rs_remRows = 0;
	
	   for(int i = 0; i < lp.nRows(); ++i)
	   {
	      const SVectorBase<R>& row = lp.rowVector(i);
	
	      if(row.size() == 1)
	      {
	         removeRowSingleton(lp, row, i);
	         rs_remRows++;
	      }
	   }
	
	   if(rs_remRows > 0)
	   {
	      MSG_INFO2((*this->spxout), (*this->spxout) <<
	                "Simplifier duplicate rows (row singleton stage) removed "
	                << rs_remRows << " rows, "
	                << rs_remRows << " non-zeros"
	                << std::endl;)
	   }
	
	#endif
	
	   if(lp.nRows() < 2)
	      return this->OKAY;
	
	   DataArray<int>    pClass(lp.nRows());           // class of parallel rows
	   DataArray<int>    classSize(lp.nRows());        // size of each class
	   Array<R>   scale(lp.nRows());            // scaling factor for each row
	   int*              idxMem = 0;
	
	   try
	   {
	      spx_alloc(idxMem, lp.nRows());
	   }
	   catch(const SPxMemoryException& x)
	   {
	      spx_free(idxMem);
	      throw x;
	   }
	
	   IdxSet idxSet(lp.nRows(), idxMem);           // set of feasible indices for new pClass
	
	   // init
	   pClass[0]    = 0;
	   scale[0]     = 0.0;
	   classSize[0] = lp.nRows();
	
	   for(int i = 1; i < lp.nRows(); ++i)
	   {
	      pClass[i] = 0;
	      scale[i]  = 0.0;
	      classSize[i] = 0;
	      idxSet.addIdx(i);
	   }
	
	   R oldVal = 0.0;
	
	   // main loop
	   for(int j = 0; j < lp.nCols(); ++j)
	   {
	      const SVectorBase<R>& col = lp.colVector(j);
	
	      for(int k = 0; k < col.size(); ++k)
	      {
	         R aij = col.value(k);
	         int  i   = col.index(k);
	
	         if(scale[i] == 0.0)
	            scale[i] = aij;
	
	         m_classSetRows[pClass[i]].add(i, aij / scale[i]);
	
	         if(--classSize[pClass[i]] == 0)
	            idxSet.addIdx(pClass[i]);
	      }
	
	      // update each parallel class with non-zero column entry
	      for(int m = 0; m < col.size(); ++m)
	      {
	         int k = pClass[col.index(m)];
	
	         if(m_classSetRows[k].size() > 0)
	         {
	            // sort classSet[k] w.r.t. scaled column values
	            ElementCompare compare;
	
	            if(m_classSetRows[k].size() > 1)
	               SPxQuicksort(m_classSetRows[k].mem(), m_classSetRows[k].size(), compare);
	
	            // use new index first
	            int classIdx = idxSet.index(0);
	            idxSet.remove(0);
	
	            for(int l = 0; l < m_classSetRows[k].size(); ++l)
	            {
	               if(l != 0 && NErel(m_classSetRows[k].value(l), oldVal, this->epsZero()))
	               {
	                  classIdx = idxSet.index(0);
	                  idxSet.remove(0);
	               }
	
	               pClass[m_classSetRows[k].index(l)] = classIdx;
	               ++classSize[classIdx];
	
	               oldVal = m_classSetRows[k].value(l);
	            }
	
	            m_classSetRows[k].clear();
	         }
	      }
	   }
	
	   spx_free(idxMem);
	
	   DataArray<bool> remRow(lp.nRows());
	
	   for(int k = 0; k < lp.nRows(); ++k)
	      m_dupRows[k].clear();
	
	   for(int k = 0; k < lp.nRows(); ++k)
	   {
	      remRow[k] = false;
	      m_dupRows[pClass[k]].add(k, 0.0);
	   }
	
	   const int nRowsOld_tmp = lp.nRows();
	   int* perm_tmp = 0;
	   spx_alloc(perm_tmp, nRowsOld_tmp);
	
	   for(int j = 0; j < nRowsOld_tmp; ++j)
	   {
	      perm_tmp[j] = 0;
	   }
	
	   int idxFirstDupRows = -1;
	   int idxLastDupRows = -1;
	   int numDelRows = 0;
	
	   for(int k = 0; k < lp.nRows(); ++k)
	   {
	      if(m_dupRows[k].size() > 1 && !(lp.rowVector(m_dupRows[k].index(0)).size() == 1))
	      {
	         idxLastDupRows = k;
	
	         if(idxFirstDupRows < 0)
	         {
	            idxFirstDupRows = k;
	         }
	
	         for(int l = 1; l < m_dupRows[k].size(); ++l)
	         {
	            int i = m_dupRows[k].index(l);
	            perm_tmp[i] = -1;
	         }
	
	         numDelRows += (m_dupRows[k].size() - 1);
	      }
	   }
	
	   {
	      int k_tmp, j_tmp = -1;
	
	      for(k_tmp = j_tmp = 0; k_tmp < nRowsOld_tmp; ++k_tmp)
	      {
	         if(perm_tmp[k_tmp] >= 0)
	            perm_tmp[k_tmp] = j_tmp++;
	      }
	   }
	
	   // store rhs and lhs changes for combined update
	   bool doChangeRanges = false;
	   VectorBase<R> newLhsVec(lp.lhs());
	   VectorBase<R> newRhsVec(lp.rhs());
	
	   for(int k = 0; k < lp.nRows(); ++k)
	   {
	      if(m_dupRows[k].size() > 1 && !(lp.rowVector(m_dupRows[k].index(0)).size() == 1))
	      {
	         MSG_DEBUG((*this->spxout) << "IMAISM53 " << m_dupRows[k].size()
	                   << " duplicate rows found" << std::endl;)
	
	         m_stat[DUPLICATE_ROW] += m_dupRows[k].size() - 1;
	
	         // index of one non-column singleton row in dupRows[k]
	         int  rowIdx    = -1;
	         int  maxLhsIdx = -1;
	         int  minRhsIdx = -1;
	         R maxLhs    = R(-infinity);
	         R minRhs    = +R(infinity);
	
	         DataArray<bool> isLhsEqualRhs(m_dupRows[k].size());
	
	         // determine strictest bounds on constraint
	         for(int l = 0; l < m_dupRows[k].size(); ++l)
	         {
	            int i = m_dupRows[k].index(l);
	            isLhsEqualRhs[l] = (lp.lhs(i) == lp.rhs(i));
	
	            ASSERT_WARN("WMAISM54", isNotZero(scale[i], R(1.0 / R(infinity))));
	
	            if(rowIdx == -1)
	            {
	               rowIdx = i;
	               maxLhs = lp.lhs(rowIdx);
	               minRhs = lp.rhs(rowIdx);
	            }
	            else
	            {
	               R scaledLhs, scaledRhs;
	               R factor = scale[rowIdx] / scale[i];
	
	               if(factor > 0)
	               {
	                  scaledLhs = (lp.lhs(i) <= R(-infinity)) ? R(-infinity) : lp.lhs(i) * factor;
	                  scaledRhs = (lp.rhs(i) >=  R(infinity)) ?  R(infinity) : lp.rhs(i) * factor;
	               }
	               else
	               {
	                  scaledLhs = (lp.rhs(i) >=  R(infinity)) ? R(-infinity) : lp.rhs(i) * factor;
	                  scaledRhs = (lp.lhs(i) <= R(-infinity)) ?  R(infinity) : lp.lhs(i) * factor;
	               }
	
	               if(scaledLhs > maxLhs)
	               {
	                  maxLhs    = scaledLhs;
	                  maxLhsIdx = i;
	               }
	
	               if(scaledRhs < minRhs)
	               {
	                  minRhs    = scaledRhs;
	                  minRhsIdx = i;
	               }
	
	               remRow[i] = true;
	            }
	         }
	
	         if(rowIdx != -1)
	         {
	            R newLhs = (maxLhs > lp.lhs(rowIdx)) ? maxLhs : lp.lhs(rowIdx);
	            R newRhs = (minRhs < lp.rhs(rowIdx)) ? minRhs : lp.rhs(rowIdx);
	
	            if(k == idxLastDupRows)
	            {
	               DataArray<int> da_perm(nRowsOld_tmp);
	
	               for(int j = 0; j < nRowsOld_tmp; ++j)
	               {
	                  da_perm[j] = perm_tmp[j];
	               }
	
	               std::shared_ptr<PostStep> ptr(new DuplicateRowsPS(lp, rowIdx, maxLhsIdx, minRhsIdx,
	                                             m_dupRows[k], scale, da_perm, isLhsEqualRhs, true,
	                                             EQrel(newLhs, newRhs), k == idxFirstDupRows));
	               m_hist.append(ptr);
	            }
	            else
	            {
	               DataArray<int> da_perm_empty(0);
	               std::shared_ptr<PostStep> ptr(new DuplicateRowsPS(lp, rowIdx, maxLhsIdx, minRhsIdx,
	                                             m_dupRows[k], scale, da_perm_empty, isLhsEqualRhs, false, EQrel(newLhs, newRhs),
	                                             k == idxFirstDupRows));
	               m_hist.append(ptr);
	            }
	
	            if(maxLhs > lp.lhs(rowIdx) || minRhs < lp.rhs(rowIdx))
	            {
	               // modify lhs and rhs of constraint rowIdx
	               doChangeRanges = true;
	
	               if(LTrel(newRhs, newLhs, feastol()))
	               {
	                  MSG_DEBUG((*this->spxout) << "IMAISM55 duplicate rows yield infeasible bounds:"
	                            << " lhs=" << newLhs
	                            << " rhs=" << newRhs << std::endl;)
	                  spx_free(perm_tmp);
	                  return this->INFEASIBLE;
	               }
	
	               // if we accept the infeasibility we should clean up the values to avoid problems later
	               if(newRhs < newLhs)
	                  newRhs = newLhs;
	
	               newLhsVec[rowIdx] = newLhs;
	               newRhsVec[rowIdx] = newRhs;
	            }
	         }
	      }
	   }
	
	   // change ranges for all modified constraints by one single call (more efficient)
	   if(doChangeRanges)
	   {
	      lp.changeRange(newLhsVec, newRhsVec);
	   }
	
	   // remove all rows by one single method call (more efficient)
	   const int nRowsOld = lp.nRows();
	   int* perm = 0;
	   spx_alloc(perm, nRowsOld);
	
	   for(int i = 0; i < nRowsOld; ++i)
	   {
	      if(remRow[i])
	      {
	         perm[i] = -1;
	         ++remRows;
	         remNzos += lp.rowVector(i).size();
	      }
	      else
	         perm[i] = 0;
	   }
	
	   lp.removeRows(perm);
	
	   for(int i = 0; i < nRowsOld; ++i)
	   {
	      // assert that the pre-computed permutation was correct
	      assert(perm[i] == perm_tmp[i]);
	
	      // update the global index mapping
	      if(perm[i] >= 0)
	         m_rIdx[perm[i]] = m_rIdx[i];
	   }
	
	   spx_free(perm);
	   spx_free(perm_tmp);
	
	   if(remRows + remNzos > 0)
	   {
	      this->m_remRows += remRows;
	      this->m_remNzos += remNzos;
	
	      MSG_INFO2((*this->spxout), (*this->spxout) << "Simplifier (duplicate rows) removed "
	                << remRows << " rows, "
	                << remNzos << " non-zeros"
	                << std::endl;)
	
	      if(remRows > this->m_minReduction * oldRows)
	         again = true;
	
	   }
	
	   return this->OKAY;
	}
	
	template <class R>
	typename SPxSimplifier<R>::Result SPxMainSM<R>::duplicateCols(SPxLPBase<R>& lp, bool& again)
	{
	
	   // This method simplifies the LP by removing duplicate columns
	   // Duplicates are detected using the algorithm of Bixby and Wagner [1987]
	
	   int remCols = 0;
	   int remNzos = 0;
	
	   // remove empty rows and columns
	   typename SPxSimplifier<R>::Result ret = removeEmpty(lp);
	
	   if(ret != this->OKAY)
	      return ret;
	
	   if(lp.nCols() < 2)
	      return this->OKAY;
	
	   DataArray<int>    pClass(lp.nCols());          // class of parallel columns
	   DataArray<int>    classSize(lp.nCols());       // size of each class
	   Array<R>   scale(lp.nCols());           // scaling factor for each column
	   int*              idxMem = 0;
	
	   try
	   {
	      spx_alloc(idxMem, lp.nCols());
	   }
	   catch(const SPxMemoryException& x)
	   {
	      spx_free(idxMem);
	      throw x;
	   }
	
	   IdxSet idxSet(lp.nCols(), idxMem);  // set of feasible indices for new pClass
	
	   // init
	   pClass[0]    = 0;
	   scale[0]     = 0.0;
	   classSize[0] = lp.nCols();
	
	   for(int j = 1; j < lp.nCols(); ++j)
	   {
	      pClass[j] = 0;
	      scale[j]  = 0.0;
	      classSize[j] = 0;
	      idxSet.addIdx(j);
	   }
	
	   R oldVal = 0.0;
	
	   // main loop
	   for(int i = 0; i < lp.nRows(); ++i)
	   {
	      const SVectorBase<R>& row = lp.rowVector(i);
	
	      for(int k = 0; k < row.size(); ++k)
	      {
	         R aij = row.value(k);
	         int  j   = row.index(k);
	
	         if(scale[j] == 0.0)
	            scale[j] = aij;
	
	         m_classSetCols[pClass[j]].add(j, aij / scale[j]);
	
	         if(--classSize[pClass[j]] == 0)
	            idxSet.addIdx(pClass[j]);
	      }
	
	      // update each parallel class with non-zero row entry
	      for(int m = 0; m < row.size(); ++m)
	      {
	         int k = pClass[row.index(m)];
	
	         if(m_classSetCols[k].size() > 0)
	         {
	            // sort classSet[k] w.r.t. scaled row values
	            ElementCompare compare;
	
	            if(m_classSetCols[k].size() > 1)
	               SPxQuicksort(m_classSetCols[k].mem(), m_classSetCols[k].size(), compare);
	
	            // use new index first
	            int classIdx = idxSet.index(0);
	            idxSet.remove(0);
	
	            for(int l = 0; l < m_classSetCols[k].size(); ++l)
	            {
	               if(l != 0 && NErel(m_classSetCols[k].value(l), oldVal, this->epsZero()))
	               {
	                  // start new parallel class
	                  classIdx = idxSet.index(0);
	                  idxSet.remove(0);
	               }
	
	               pClass[m_classSetCols[k].index(l)] = classIdx;
	               ++classSize[classIdx];
	
	               oldVal = m_classSetCols[k].value(l);
	            }
	
	            m_classSetCols[k].clear();
	         }
	      }
	   }
	
	   spx_free(idxMem);
	
	   DataArray<bool> remCol(lp.nCols());
	   DataArray<bool> fixAndRemCol(lp.nCols());
	
	   for(int k = 0; k < lp.nCols(); ++k)
	      m_dupCols[k].clear();
	
	   for(int k = 0; k < lp.nCols(); ++k)
	   {
	      remCol[k] = false;
	      fixAndRemCol[k] = false;
	      m_dupCols[pClass[k]].add(k, 0.0);
	   }
	
	   bool hasDuplicateCol = false;
	   DataArray<int>  m_perm_empty(0);
	
	   for(int k = 0; k < lp.nCols(); ++k)
	   {
	      if(m_dupCols[k].size() > 1 && !(lp.colVector(m_dupCols[k].index(0)).size() == 1))
	      {
	         MSG_DEBUG((*this->spxout) << "IMAISM58 " << m_dupCols[k].size()
	                   << " duplicate columns found" << std::endl;)
	
	         if(!hasDuplicateCol)
	         {
	            std::shared_ptr<PostStep> ptr(new DuplicateColsPS(lp, 0, 0, 1.0, m_perm_empty, true));
	            m_hist.append(ptr);
	            hasDuplicateCol = true;
	         }
	
	         for(int l = 0; l < m_dupCols[k].size(); ++l)
	         {
	            for(int m = 0; m < m_dupCols[k].size(); ++m)
	            {
	               int j1  = m_dupCols[k].index(l);
	               int j2  = m_dupCols[k].index(m);
	
	               if(l != m && !remCol[j1] && !remCol[j2])
	               {
	                  R cj1 = lp.maxObj(j1);
	                  R cj2 = lp.maxObj(j2);
	
	                  // A.j1 = factor * A.j2
	                  R factor = scale[j1] / scale[j2];
	                  R objDif = cj1 - cj2 * scale[j1] / scale[j2];
	
	                  ASSERT_WARN("WMAISM59", isNotZero(factor, this->epsZero()));
	
	                  if(isZero(objDif, this->epsZero()))
	                  {
	                     // case 1: objectives also duplicate
	
	                     // if 0 is not within the column bounds, we are not able to postsolve if the aggregated column has
	                     // status ZERO, hence we skip this case
	                     if(LErel(lp.lower(j1), R(0.0)) && GErel(lp.upper(j1), R(0.0))
	                           && LErel(lp.lower(j2), R(0.0)) && GErel(lp.upper(j2), R(0.0)))
	                     {
	                        std::shared_ptr<PostStep> ptr(new DuplicateColsPS(lp, j1, j2, factor, m_perm_empty));
	                        // variable substitution xj2' := xj2 + factor * xj1 <=> xj2 = -factor * xj1 + xj2'
	                        m_hist.append(ptr);
	
	                        // update bounds of remaining column j2 (new column j2')
	                        if(factor > 0)
	                        {
	                           if(lp.lower(j2) <= R(-infinity) || lp.lower(j1) <= R(-infinity))
	                              lp.changeLower(j2, R(-infinity));
	                           else
	                              lp.changeLower(j2, lp.lower(j2) + factor * lp.lower(j1));
	
	                           if(lp.upper(j2) >= R(infinity) || lp.upper(j1) >= R(infinity))
	                              lp.changeUpper(j2, R(infinity));
	                           else
	                              lp.changeUpper(j2, lp.upper(j2) + factor * lp.upper(j1));
	                        }
	                        else if(factor < 0)
	                        {
	                           if(lp.lower(j2) <= R(-infinity) || lp.upper(j1) >= R(infinity))
	                              lp.changeLower(j2, R(-infinity));
	                           else
	                              lp.changeLower(j2, lp.lower(j2) + factor * lp.upper(j1));
	
	                           if(lp.upper(j2) >= R(infinity) || lp.lower(j1) <= R(-infinity))
	                              lp.changeUpper(j2, R(infinity));
	                           else
	                              lp.changeUpper(j2, lp.upper(j2) + factor * lp.lower(j1));
	                        }
	
	                        MSG_DEBUG((*this->spxout) << "IMAISM60 two duplicate columns " << j1
	                                  << ", " << j2
	                                  << " replaced by one" << std::endl;)
	
	                        remCol[j1] = true;
	
	                        ++m_stat[SUB_DUPLICATE_COL];
	                     }
	                     else
	                     {
	                        MSG_DEBUG((*this->spxout) << "IMAISM80 not removing two duplicate columns " << j1
	                                  << ", " << j2
	                                  << " because zero not contained in their bounds" << std::endl;)
	                     }
	                  }
	                  else
	                  {
	                     // case 2: objectives not duplicate
	                     // considered for maximization sense
	                     if(lp.lower(j2) <= R(-infinity))
	                     {
	                        if(factor > 0 && objDif > 0)
	                        {
	                           if(lp.upper(j1) >= R(infinity))
	                           {
	                              MSG_DEBUG((*this->spxout) << "IMAISM75 LP unbounded" << std::endl;)
	                              return this->UNBOUNDED;
	                           }
	
	                           // fix j1 at upper bound
	                           MSG_DEBUG((*this->spxout) << "IMAISM61 two duplicate columns " << j1
	                                     << ", " << j2
	                                     << " first one fixed at upper bound=" << lp.upper(j1) << std::endl;)
	
	                           std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j1, lp.upper(j1)));
	                           m_hist.append(ptr);
	                           lp.changeLower(j1, lp.upper(j1));
	                        }
	                        else if(factor < 0 && objDif < 0)
	                        {
	                           if(lp.lower(j1) <= R(-infinity))
	                           {
	                              MSG_DEBUG((*this->spxout) << "IMAISM76 LP unbounded" << std::endl;)
	                              return this->UNBOUNDED;
	                           }
	
	                           // fix j1 at lower bound
	                           MSG_DEBUG((*this->spxout) << "IMAISM62 two duplicate columns " << j1
	                                     << ", " << j2
	                                     << " first one fixed at lower bound=" << lp.lower(j1) << std::endl;)
	
	                           std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j1, lp.lower(j1)));
	                           m_hist.append(ptr);
	                           lp.changeUpper(j1, lp.lower(j1));
	                        }
	                     }
	                     else if(lp.upper(j2) >= R(infinity))
	                     {
	                        // fix j1 at upper bound
	                        if(factor < 0 && objDif > 0)
	                        {
	                           if(lp.upper(j1) >= R(infinity))
	                           {
	                              MSG_DEBUG((*this->spxout) << "IMAISM77 LP unbounded" << std::endl;)
	                              return this->UNBOUNDED;
	                           }
	
	                           // fix j1 at upper bound
	                           MSG_DEBUG((*this->spxout) << "IMAISM63 two duplicate columns " << j1
	                                     << ", " << j2
	                                     << " first one fixed at upper bound=" << lp.upper(j1) << std::endl;)
	
	                           std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j1, lp.upper(j1)));
	                           m_hist.append(ptr);
	                           lp.changeLower(j1, lp.upper(j1));
	                        }
	
	                        // fix j1 at lower bound
	                        else if(factor > 0 && objDif < 0)
	                        {
	                           if(lp.lower(j1) <= R(-infinity))
	                           {
	                              MSG_DEBUG((*this->spxout) << "IMAISM78 LP unbounded" << std::endl;)
	                              return this->UNBOUNDED;
	                           }
	
	                           // fix j1 at lower bound
	                           MSG_DEBUG((*this->spxout) << "IMAISM64 two duplicate columns " << j1
	                                     << ", " << j2
	                                     << " first one fixed at lower bound=" << lp.lower(j1) << std::endl;)
	
	                           std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j1, lp.lower(j1)));
	                           m_hist.append(ptr);
	                           lp.changeUpper(j1, lp.lower(j1));
	                        }
	                     }
	
	                     if(EQrel(lp.lower(j1), lp.upper(j1), feastol()))
	                     {
	                        remCol[j1] = true;
	                        fixAndRemCol[j1] = true;
	
	                        ++m_stat[FIX_DUPLICATE_COL];
	                     }
	                  }
	               }
	            }
	         }
	      }
	   }
	
	   for(int j = 0; j < lp.nCols(); ++j)
	   {
	      if(fixAndRemCol[j])
	      {
	         assert(remCol[j]);
	
	         // correctIdx == false, because the index mapping will be handled by the postsolving in DuplicateColsPS
	         fixColumn(lp, j, false);
	      }
	   }
	
	   // remove all columns by one single method call (more efficient)
	   const int nColsOld = lp.nCols();
	   int* perm = 0;
	   spx_alloc(perm, nColsOld);
	
	   for(int j = 0; j < nColsOld; ++j)
	   {
	      if(remCol[j])
	      {
	         perm[j] = -1;
	         ++remCols;
	         remNzos += lp.colVector(j).size();
	      }
	      else
	         perm[j] = 0;
	   }
	
	   lp.removeCols(perm);
	
	   for(int j = 0; j < nColsOld; ++j)
	   {
	      if(perm[j] >= 0)
	         m_cIdx[perm[j]] = m_cIdx[j];
	   }
	
	   DataArray<int> da_perm(nColsOld);
	
	   for(int j = 0; j < nColsOld; ++j)
	   {
	      da_perm[j] = perm[j];
	   }
	
	   if(hasDuplicateCol)
	   {
	      std::shared_ptr<PostStep> ptr(new DuplicateColsPS(lp, 0, 0, 1.0, da_perm, false, true));
	      m_hist.append(ptr);
	   }
	
	   spx_free(perm);
	
	   assert(remCols > 0 || remNzos == 0);
	
	   if(remCols > 0)
	   {
	      this->m_remCols += remCols;
	      this->m_remNzos += remNzos;
	
	      MSG_INFO2((*this->spxout), (*this->spxout) << "Simplifier (duplicate columns) removed "
	                << remCols << " cols, "
	                << remNzos << " non-zeros"
	                << std::endl;)
	
	      if(remCols > this->m_minReduction * nColsOld)
	         again = true;
	   }
	
	   return this->OKAY;
	}
	
	template <class R>
	void SPxMainSM<R>::fixColumn(SPxLPBase<R>& lp, int j, bool correctIdx)
	{
	
	   assert(EQrel(lp.lower(j), lp.upper(j), feastol()));
	
	   R lo            = lp.lower(j);
	   R up            = lp.upper(j);
	   const SVectorBase<R>& col = lp.colVector(j);
	   R mid           = lo;
	
	   // use the center value between slightly different bounds to improve numerics
	   if(NE(lo, up))
	      mid = (up + lo) / 2.0;
	
	   assert(LT(lo, R(infinity)) && GT(lo, R(-infinity)));
	   assert(LT(up, R(infinity)) && GT(up, R(-infinity)));
	
	   MSG_DEBUG((*this->spxout) << "IMAISM66 fix variable x" << j
	             << ": lower=" << lo
	             << " upper=" << up
	             << "to new value: " << mid
	             << std::endl;)
	
	   if(isNotZero(lo, this->epsZero()))
	   {
	      for(int k = 0; k < col.size(); ++k)
	      {
	         int i = col.index(k);
	
	         if(lp.rhs(i) < R(infinity))
	         {
	            R y     = mid * col.value(k);
	            R scale = maxAbs(lp.rhs(i), y);
	
	            if(scale < 1.0)
	               scale = 1.0;
	
	            R rhs = (lp.rhs(i) / scale) - (y / scale);
	
	            if(isZero(rhs, this->epsZero()))
	               rhs = 0.0;
	            else
	               rhs *= scale;
	
	            MSG_DEBUG((*this->spxout) << "IMAISM67 row " << i
	                      << ": rhs=" << rhs
	                      << " (" << lp.rhs(i)
	                      << ") aij=" << col.value(k)
	                      << std::endl;)
	
	            lp.changeRhs(i, rhs);
	         }
	
	         if(lp.lhs(i) > R(-infinity))
	         {
	            R y     = mid * col.value(k);
	            R scale = maxAbs(lp.lhs(i), y);
	
	            if(scale < 1.0)
	               scale = 1.0;
	
	            R lhs = (lp.lhs(i) / scale) - (y / scale);
	
	            if(isZero(lhs, this->epsZero()))
	               lhs = 0.0;
	            else
	               lhs *= scale;
	
	            MSG_DEBUG((*this->spxout) << "IMAISM68 row " << i
	                      << ": lhs=" << lhs
	                      << " (" << lp.lhs(i)
	                      << ") aij=" << col.value(k)
	                      << std::endl;)
	
	            lp.changeLhs(i, lhs);
	         }
	
	         assert(lp.lhs(i) <= lp.rhs(i) + feastol());
	      }
	   }
	
	   std::shared_ptr<PostStep> ptr(new FixVariablePS(lp, *this, j, lp.lower(j), correctIdx));
	   m_hist.append(ptr);
	}
	
	template <class R>
	typename SPxSimplifier<R>::Result SPxMainSM<R>::simplify(SPxLPBase<R>& lp, R eps, R ftol, R otol,
	      Real remainingTime,
	      bool keepbounds, uint32_t seed)
	{
	   // transfer message handler
	   this->spxout = lp.spxout;
	   assert(this->spxout != 0);
	
	   m_thesense = lp.spxSense();
	   this->m_timeUsed->reset();
	   this->m_timeUsed->start();
	
	   this->m_objoffset = 0.0;
	   m_cutoffbound = R(-infinity);
	   m_pseudoobj = R(-infinity);
	
	   this->m_remRows = 0;
	   this->m_remCols = 0;
	   this->m_remNzos = 0;
	   this->m_chgBnds = 0;
	   this->m_chgLRhs = 0;
	   this->m_keptBnds = 0;
	   this->m_keptLRhs = 0;
	
	   m_result     = this->OKAY;
	   bool   again = true;
	   int nrounds = 0;
	
	   if(m_hist.size() > 0)
	   {
	      m_hist.clear();
	   }
	
	   m_hist.reSize(0);
	   m_postsolved = false;
	
	   if(eps < 0.0)
	      throw SPxInterfaceException("XMAISM30 Cannot use negative this->epsilon in simplify().");
	
	   if(ftol < 0.0)
	      throw SPxInterfaceException("XMAISM31 Cannot use negative feastol in simplify().");
	
	   if(otol < 0.0)
	      throw SPxInterfaceException("XMAISM32 Cannot use negative opttol in simplify().");
	
	   m_epsilon = eps;
	   m_feastol = ftol;
	   m_opttol = otol;
	
	
	   MSG_INFO2((*this->spxout),
	             int numRangedRows = 0;
	             int numBoxedCols = 0;
	             int numEqualities = 0;
	
	             for(int i = 0; i < lp.nRows(); ++i)
	{
	   if(lp.lhs(i) > R(-infinity) && lp.rhs(i) < R(infinity))
	      {
	         if(EQ(lp.lhs(i), lp.rhs(i)))
	            ++numEqualities;
	         else
	            ++numRangedRows;
	      }
	   }
	   for(int j = 0; j < lp.nCols(); ++j)
	   if(lp.lower(j) > R(-infinity) && lp.upper(j) < R(infinity))
	      ++numBoxedCols;
	
	      (*this->spxout) << "LP has "
	                      << numEqualities << " equations, "
	                      << numRangedRows << " ranged rows, "
	                      << numBoxedCols << " boxed columns"
	                      << std::endl;
	               )
	
	         m_stat.reSize(17);
	
	   for(int k = 0; k < m_stat.size(); ++k)
	      m_stat[k] = 0;
	
	   m_addedcols = 0;
	   handleRowObjectives(lp);
	
	   m_prim.reDim(lp.nCols());
	   m_slack.reDim(lp.nRows());
	   m_dual.reDim(lp.nRows());
	   m_redCost.reDim(lp.nCols());
	   m_cBasisStat.reSize(lp.nCols());
	   m_rBasisStat.reSize(lp.nRows());
	   m_cIdx.reSize(lp.nCols());
	   m_rIdx.reSize(lp.nRows());
	
	   m_classSetRows.reSize(lp.nRows());
	   m_classSetCols.reSize(lp.nCols());
	   m_dupRows.reSize(lp.nRows());
	   m_dupCols.reSize(lp.nCols());
	
	   m_keepbounds = keepbounds;
	
	   for(int i = 0; i < lp.nRows(); ++i)
	      m_rIdx[i] = i;
	
	   for(int j = 0; j < lp.nCols(); ++j)
	      m_cIdx[j] = j;
	
	   // round extreme values (set all values smaller than this->eps to zero and all values bigger than R(infinity)/5 to R(infinity))
	#if EXTREMES
	   handleExtremes(lp);
	#endif
	
	   // main presolving loop
	   while(again && m_result == this->OKAY)
	   {
	      nrounds++;
	      MSG_INFO3((*this->spxout), (*this->spxout) << "Round " << nrounds << ":" << std::endl;)
	      again = false;
	
	#if ROWS_SPXMAINSM
	
	      if(m_result == this->OKAY)
	         m_result = simplifyRows(lp, again);
	
	#endif
	
	#if COLS_SPXMAINSM
	
	      if(m_result == this->OKAY)
	         m_result = simplifyCols(lp, again);
	
	#endif
	
	#if DUAL_SPXMAINSM
	
	      if(m_result == this->OKAY)
	         m_result = simplifyDual(lp, again);
	
	#endif
	
	#if DUPLICATE_ROWS
	
	      if(m_result == this->OKAY)
	         m_result = duplicateRows(lp, again);
	
	#endif
	
	#if DUPLICATE_COLS
	
	      if(m_result == this->OKAY)
	         m_result = duplicateCols(lp, again);
	
	#endif
	
	      if(!again)
	      {
	#if TRIVIAL_HEURISTICS
	         trivialHeuristic(lp);
	#endif
	
	#if PSEUDOOBJ
	         propagatePseudoobj(lp);
	#endif
	
	#if MULTI_AGGREGATE
	
	         if(m_result == this->OKAY)
	            m_result = multiaggregation(lp, again);
	
	#endif
	      }
	
	   }
	
	   // preprocessing detected infeasibility or unboundedness
	   if(m_result != this->OKAY)
	   {
	      MSG_INFO1((*this->spxout), (*this->spxout) << "Simplifier result: " << static_cast<int>
	                (m_result) << std::endl;)
	      return m_result;
	   }
	
	   this->m_remCols -= m_addedcols;
	   this->m_remNzos -= m_addedcols;
	   MSG_INFO1((*this->spxout), (*this->spxout) << "Simplifier removed "
	             << this->m_remRows << " rows, "
	             << this->m_remCols << " columns, "
	             << this->m_remNzos << " nonzeros, "
	             << this->m_chgBnds << " col bounds, "
	             << this->m_chgLRhs << " row bounds"
	             << std::endl;)
	
	   if(keepbounds)
	      MSG_INFO2((*this->spxout), (*this->spxout) << "Simplifier kept "
	                << this->m_keptBnds << " column bounds, "
	                << this->m_keptLRhs << " row bounds"
	                << std::endl;)
	
	      MSG_INFO1((*this->spxout), (*this->spxout) << "Reduced LP has "
	                << lp.nRows() << " rows "
	                << lp.nCols() << " columns "
	                << lp.nNzos() << " nonzeros"
	                << std::endl;)
	
	      MSG_INFO2((*this->spxout),
	                int numRangedRows = 0;
	                int numBoxedCols  = 0;
	                int numEqualities = 0;
	
	                for(int i = 0; i < lp.nRows(); ++i)
	   {
	      if(lp.lhs(i) > R(-infinity) && lp.rhs(i) < R(infinity))
	         {
	            if(EQ(lp.lhs(i), lp.rhs(i)))
	               ++numEqualities;
	            else
	               ++numRangedRows;
	         }
	      }
	   for(int j = 0; j < lp.nCols(); ++j)
	   if(lp.lower(j) > R(-infinity) && lp.upper(j) < R(infinity))
	      ++numBoxedCols;
	
	      (*this->spxout) << "Reduced LP has "
	                      << numEqualities << " equations, "
	                      << numRangedRows << " ranged rows, "
	                      << numBoxedCols << " boxed columns"
	                      << std::endl;
	               )
	
	         if(lp.nCols() == 0 && lp.nRows() == 0)
	         {
	            MSG_INFO1((*this->spxout), (*this->spxout) << "Simplifier removed all rows and columns" <<
	                      std::endl;)
	            m_result = this->VANISHED;
	         }
	
	   MSG_INFO2((*this->spxout), (*this->spxout) << "\nSimplifier performed:\n"
	             << m_stat[EMPTY_ROW]            << " empty rows\n"
	             << m_stat[FREE_ROW]             << " free rows\n"
	             << m_stat[SINGLETON_ROW]        << " singleton rows\n"
	             << m_stat[FORCE_ROW]            << " forcing rows\n"
	             << m_stat[EMPTY_COL]            << " empty columns\n"
	             << m_stat[FIX_COL]              << " fixed columns\n"
	             << m_stat[FREE_ZOBJ_COL]        << " free columns with zero objective\n"
	             << m_stat[ZOBJ_SINGLETON_COL]   << " singleton columns with zero objective\n"
	             << m_stat[DOUBLETON_ROW]        << " singleton columns combined with a doubleton equation\n"
	             << m_stat[FREE_SINGLETON_COL]   << " free singleton columns\n"
	             << m_stat[DOMINATED_COL]        << " dominated columns\n"
	             << m_stat[WEAKLY_DOMINATED_COL] << " weakly dominated columns\n"
	             << m_stat[DUPLICATE_ROW]        << " duplicate rows\n"
	             << m_stat[FIX_DUPLICATE_COL]    << " duplicate columns (fixed)\n"
	             << m_stat[SUB_DUPLICATE_COL]    << " duplicate columns (substituted)\n"
	             << m_stat[AGGREGATION]          << " variable aggregations\n"
	             << m_stat[MULTI_AGG]            << " multi aggregations\n"
	             << std::endl;);
	
	   this->m_timeUsed->stop();
	
	   return m_result;
	}
	
	template <class R>
	void SPxMainSM<R>::unsimplify(const VectorBase<R>& x, const VectorBase<R>& y,
	                              const VectorBase<R>& s, const VectorBase<R>& r,
	                              const typename SPxSolverBase<R>::VarStatus rows[],
	                              const typename SPxSolverBase<R>::VarStatus cols[], bool isOptimal)
	{
	   MSG_INFO1((*this->spxout), (*this->spxout) << " --- unsimplifying solution and basis" << std::endl;)
	   assert(x.dim() <= m_prim.dim());
	   assert(y.dim() <= m_dual.dim());
	   assert(x.dim() == r.dim());
	   assert(y.dim() == s.dim());
	
	   // assign values of variables in reduced LP
	   // NOTE: for maximization problems, we have to switch signs of dual and reduced cost values,
	   // since simplifier assumes minimization problem
	   for(int j = 0; j < x.dim(); ++j)
	   {
	      m_prim[j] = isZero(x[j], this->epsZero()) ? 0.0 : x[j];
	      m_redCost[j] = isZero(r[j], this->epsZero()) ? 0.0 : (m_thesense == SPxLPBase<R>::MAXIMIZE ? -r[j] :
	                     r[j]);
	      m_cBasisStat[j] = cols[j];
	   }
	
	   for(int i = 0; i < y.dim(); ++i)
	   {
	      m_dual[i] = isZero(y[i], this->epsZero()) ? 0.0 : (m_thesense == SPxLPBase<R>::MAXIMIZE ? -y[i] :
	                  y[i]);
	      m_slack[i] = isZero(s[i], this->epsZero()) ? 0.0 : s[i];
	      m_rBasisStat[i] = rows[i];
	   }
	
	   // undo preprocessing
	   for(int k = m_hist.size() - 1; k >= 0; --k)
	   {
	      MSG_DEBUG(std::cout << "unsimplifying " << m_hist[k]->getName() << "\n");
	
	      try
	      {
	         m_hist[k]->execute(m_prim, m_dual, m_slack, m_redCost, m_cBasisStat, m_rBasisStat, isOptimal);
	      }
	      catch(const SPxException& ex)
	      {
	         MSG_INFO1((*this->spxout), (*this->spxout) << "Exception thrown while unsimplifying " <<
	                   m_hist[k]->getName() << ":\n" << ex.what() << "\n");
	         throw SPxInternalCodeException("XMAISM00 Exception thrown during unsimply().");
	      }
	
	      m_hist.reSize(k);
	   }
	
	   // for maximization problems, we have to switch signs of dual and reduced cost values back
	   if(m_thesense == SPxLPBase<R>::MAXIMIZE)
	   {
	      for(int j = 0; j < m_redCost.dim(); ++j)
	         m_redCost[j] = -m_redCost[j];
	
	      for(int i = 0; i < m_dual.dim(); ++i)
	         m_dual[i] = -m_dual[i];
	   }
	
	   if(m_addedcols > 0)
	   {
	      assert(m_prim.dim() >= m_addedcols);
	      m_prim.reDim(m_prim.dim() - m_addedcols);
	      m_redCost.reDim(m_redCost.dim() - m_addedcols);
	      m_cBasisStat.reSize(m_cBasisStat.size() - m_addedcols);
	      m_cIdx.reSize(m_cIdx.size() - m_addedcols);
	   }
	
	#ifdef CHECK_BASIC_DIM
	   int numBasis = 0;
	
	   for(int rs = 0; rs < m_rBasisStat.size(); ++rs)
	   {
	      if(m_rBasisStat[rs] == SPxSolverBase<R>::BASIC)
	         numBasis ++;
	   }
	
	   for(int cs = 0; cs < m_cBasisStat.size(); ++cs)
	   {
	      if(m_cBasisStat[cs] == SPxSolverBase<R>::BASIC)
	         numBasis ++;
	   }
	
	   if(numBasis != m_rBasisStat.size())
	   {
	      throw SPxInternalCodeException("XMAISM26 Dimension doesn't match after this step.");
	   }
	
	#endif
	
	   m_hist.clear();
	   m_postsolved = true;
	}
	
	// Pretty-printing of solver status.
	template <class R>
	std::ostream& operator<<(std::ostream& os, const typename SPxSimplifier<R>::Result& status)
	{
	   switch(status)
	   {
	   case SPxSimplifier<R>::OKAY:
	      os << "SUCCESS";
	      break;
	
	   case SPxSimplifier<R>::INFEASIBLE:
	      os << "INFEASIBLE";
	      break;
	
	   case SPxSimplifier<R>::DUAL_INFEASIBLE:
	      os << "DUAL_INFEASIBLE";
	      break;
	
	   case SPxSimplifier<R>::UNBOUNDED:
	      os << "UNBOUNDED";
	      break;
	
	   case SPxSimplifier<R>::VANISHED:
	      os << "VANISHED";
	      break;
	
	   default:
	      os << "UNKNOWN";
	      break;
	   }
	
	   return os;
	}
	
	} //namespace soplex