spxlpbase_real.cpp

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*                  This file is part of the class library                   */
/*       SoPlex --- the Sequential object-oriented simPlex.                  */
/*                                                                           */
/*    Copyright (C) 1996-2019 Konrad-Zuse-Zentrum                            */
/*                            fuer Informationstechnik Berlin                */
/*                                                                           */
/*  SoPlex is distributed under the terms of the ZIB Academic Licence.       */
/*                                                                           */
/*  You should have received a copy of the ZIB Academic License              */
/*  along with SoPlex; see the file COPYING. If not email to soplex@zib.de.  */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

/**@file  spxlpbase_real.cpp
 * @brief Saving LPs with Real values in a form suitable for SoPlex.
 */

#include <assert.h>
#include <stdio.h>
#include <ctype.h>
#include <iostream>

#include "soplex/spxdefines.h"
#include "soplex/spxlpbase.h"
#include "soplex/spxout.h"
#include "soplex/mpsinput.h"
#include "soplex/exceptions.h"
#include "soplex/spxscaler.h"

namespace soplex
{

template<>
void SPxLPBase<Real>::unscaleLP()
{
   MSG_INFO3((*spxout), (*spxout) << "remove persistent scaling of LP" << std::endl;)

   if(lp_scaler)
      lp_scaler->unscale(*this);
   else
      MSG_INFO3((*spxout), (*spxout) << "no LP scaler available" << std::endl;)
   }

template<>
void SPxLPBase<Real>::computePrimalActivity(const VectorBase<Real>& primal,
      VectorBase<Real>& activity, const bool unscaled) const
{
   if(primal.dim() != nCols())
      throw SPxInternalCodeException("XSPXLP01 Primal vector for computing row activity has wrong dimension");

   if(activity.dim() != nRows())
      throw SPxInternalCodeException("XSPXLP03 Activity vector computing row activity has wrong dimension");

   int c;

   for(c = 0; c < nCols() && primal[c] == 0; c++)
      ;

   if(c >= nCols())
   {
      activity.clear();
      return;
   }

   DSVector tmp(nRows());

   if(unscaled && _isScaled)
   {
      lp_scaler->getColUnscaled(*this, c, tmp);
      activity = tmp;
   }
   else
      activity = colVector(c);

   activity *= primal[c];
   c++;

   for(; c < nCols(); c++)
   {
      if(primal[c] != 0)
      {
         if(unscaled && _isScaled)
         {
            lp_scaler->getColUnscaled(*this, c, tmp);
            activity.multAdd(primal[c], tmp);
         }
         else
            activity.multAdd(primal[c], colVector(c));
      }
   }
}

template<>
void SPxLPBase<Real>::computeDualActivity(const VectorBase<Real>& dual, VectorBase<Real>& activity,
      const bool unscaled) const
{
   if(dual.dim() != nRows())
      throw SPxInternalCodeException("XSPXLP02 Dual vector for computing dual activity has wrong dimension");

   if(activity.dim() != nCols())
      throw SPxInternalCodeException("XSPXLP04 Activity vector computing dual activity has wrong dimension");

   int r;

   for(r = 0; r < nRows() && dual[r] == 0; r++)
      ;

   if(r >= nRows())
   {
      activity.clear();
      return;
   }

   DSVector tmp(nCols());

   if(unscaled && _isScaled)
   {
      lp_scaler->getRowUnscaled(*this, r, tmp);
      activity = tmp;
   }
   else
      activity = rowVector(r);

   activity *= dual[r];
   r++;

   for(; r < nRows(); r++)
   {
      if(dual[r] != 0)
      {
         if(unscaled && _isScaled)
         {
            lp_scaler->getRowUnscaled(*this, r, tmp);
            activity.multAdd(dual[r], tmp);
         }
         else
            activity.multAdd(dual[r], rowVector(r));
      }
   }
}

template<>
Real SPxLPBase<Real>::maxAbsNzo(bool unscaled) const
{
   Real maxi = 0.0;

   if(unscaled && _isScaled)
   {
      assert(lp_scaler != 0);

      for(int i = 0; i < nCols(); ++i)
      {
         Real m = lp_scaler->getColMaxAbsUnscaled(*this, i);

         if(m > maxi)
            maxi = m;
      }
   }
   else
   {
      for(int i = 0; i < nCols(); ++i)
      {
         Real m = colVector(i).maxAbs();

         if(m > maxi)
            maxi = m;
      }
   }

   assert(maxi >= 0.0);

   return maxi;
}

template<>
Real SPxLPBase<Real>::minAbsNzo(bool unscaled) const
{
   Real mini = infinity;

   if(unscaled && _isScaled)
   {
      assert(lp_scaler != 0);

      for(int i = 0; i < nCols(); ++i)
      {
         Real m = lp_scaler->getColMinAbsUnscaled(*this, i);

         if(m < mini)
            mini = m;
      }
   }
   else
   {
      for(int i = 0; i < nCols(); ++i)
      {
         Real m = colVector(i).minAbs();

         if(m < mini)
            mini = m;
      }
   }

   assert(mini >= 0.0);

   return mini;
}

/// Gets unscaled objective vector.
template <>
void SPxLPBase<Real>::getObjUnscaled(VectorBase<Real>& pobj) const
{
   if(_isScaled)
   {
      assert(lp_scaler);
      lp_scaler->getMaxObjUnscaled(*this, pobj);
   }
   else
   {
      pobj = LPColSetBase<Real>::maxObj();
   }

   if(spxSense() == MINIMIZE)
      pobj *= -1.0;
}

/// Gets unscaled row vector of row \p i.
template<>
void SPxLPBase<Real>::getRowVectorUnscaled(int i, DSVectorBase<Real>& vec) const
{
   assert(i >= 0 && i < nRows());

   if(_isScaled)
      lp_scaler->getRowUnscaled(*this, i, vec);
   else
      vec = DSVectorBase<Real>(LPRowSetBase<Real>::rowVector(i));
}

/// Gets unscaled right hand side vector.
template<>
void SPxLPBase<Real>::getRhsUnscaled(VectorBase<Real>& vec) const
{
   if(_isScaled)
      lp_scaler->getRhsUnscaled(*this, vec);
   else
      vec = LPRowSetBase<Real>::rhs();
}

/// Returns unscaled right hand side of row number \p i.
template<>
Real SPxLPBase<Real>::rhsUnscaled(int i) const
{
   assert(i >= 0 && i < nRows());

   if(_isScaled)
      return lp_scaler->rhsUnscaled(*this, i);
   else
      return LPRowSetBase<Real>::rhs(i);
}

/// Returns unscaled right hand side of row with identifier \p id.
template<>
Real SPxLPBase<Real>::rhsUnscaled(const SPxRowId& id) const
{
   assert(id.isValid());
   return rhsUnscaled(number(id));
}

/// Returns unscaled left hand side vector.
template<>
void SPxLPBase<Real>::getLhsUnscaled(VectorBase<Real>& vec) const
{
   if(_isScaled)
      lp_scaler->getLhsUnscaled(*this, vec);
   else
      vec = LPRowSetBase<Real>::lhs();
}

/// Returns unscaled left hand side of row number \p i.
template<>
Real SPxLPBase<Real>::lhsUnscaled(int i) const
{
   assert(i >= 0 && i < nRows());

   if(_isScaled)
      return lp_scaler->lhsUnscaled(*this, i);
   else
      return LPRowSetBase<Real>::lhs(i);
}

/// Returns left hand side of row with identifier \p id.
template<>
Real SPxLPBase<Real>::lhsUnscaled(const SPxRowId& id) const
{
   assert(id.isValid());
   return lhsUnscaled(number(id));
}

/// Gets column vector of column \p i.
template<>
void SPxLPBase<Real>::getColVectorUnscaled(int i, DSVectorBase<Real>& vec) const
{
   assert(i >= 0 && i < nCols());

   if(_isScaled)
      lp_scaler->getColUnscaled(*this, i, vec);
   else
      vec = LPColSetBase<Real>::colVector(i);
}

/// Gets column vector of column with identifier \p id.
template<>
void SPxLPBase<Real>::getColVectorUnscaled(const SPxColId& id, DSVectorBase<Real>& vec) const
{
   assert(id.isValid());
   getColVectorUnscaled(number(id), vec);
}

/// Returns unscaled objective value of column \p i.
template<>
Real SPxLPBase<Real>::objUnscaled(int i) const
{
   assert(i >= 0 && i < nCols());
   Real res;

   if(_isScaled)
   {
      res = lp_scaler->maxObjUnscaled(*this, i);
   }
   else
   {
      res = maxObj(i);
   }

   if(spxSense() == MINIMIZE)
      res *= -1;

   return res;
}

/// Returns unscaled objective value of column with identifier \p id.
template<>
Real SPxLPBase<Real>::objUnscaled(const SPxColId& id) const
{
   assert(id.isValid());
   return objUnscaled(number(id));
}

/// Returns unscaled objective vector for maximization problem.
template<>
void SPxLPBase<Real>::maxObjUnscaled(VectorBase<Real>& vec) const
{
   if(_isScaled)
      lp_scaler->getMaxObjUnscaled(*this, vec);
   else
      vec = LPColSetBase<Real>::maxObj();
}

/// Returns unscaled objective value of column \p i for maximization problem.
template<>
Real SPxLPBase<Real>::maxObjUnscaled(int i) const
{
   assert(i >= 0 && i < nCols());

   if(_isScaled)
      return lp_scaler->maxObjUnscaled(*this, i);
   else
      return LPColSetBase<Real>::maxObj(i);
}

/// Returns unscaled objective value of column with identifier \p id for maximization problem.
template<>
Real SPxLPBase<Real>::maxObjUnscaled(const SPxColId& id) const
{
   assert(id.isValid());
   return maxObjUnscaled(number(id));
}

/// Returns unscaled upper bound vector
template<>
void SPxLPBase<Real>::getUpperUnscaled(DVector& vec) const
{
   if(_isScaled)
      lp_scaler->getUpperUnscaled(*this, vec);
   else
      vec = DVector(LPColSetBase<Real>::upper());
}

/// Returns unscaled upper bound of column \p i.
template<>
Real SPxLPBase<Real>::upperUnscaled(int i) const
{
   assert(i >= 0 && i < nCols());

   if(_isScaled)
      return lp_scaler->upperUnscaled(*this, i);
   else
      return LPColSetBase<Real>::upper(i);
}

/// Returns unscaled upper bound of column with identifier \p id.
template<>
Real SPxLPBase<Real>::upperUnscaled(const SPxColId& id) const
{
   assert(id.isValid());
   return upperUnscaled(number(id));
}

/// Returns unscaled lower bound vector.
template<>
void SPxLPBase<Real>::getLowerUnscaled(DVector& vec) const
{
   if(_isScaled)
      lp_scaler->getLowerUnscaled(*this, vec);
   else
      vec = DVector(LPColSetBase<Real>::lower());
}

/// Returns unscaled lower bound of column \p i.
template<>
Real SPxLPBase<Real>::lowerUnscaled(int i) const
{
   assert(i >= 0 && i < nCols());

   if(_isScaled)
      return lp_scaler->lowerUnscaled(*this, i);
   else
      return LPColSetBase<Real>::lower(i);
}

/// Returns unscaled lower bound of column with identifier \p id.
template<>
Real SPxLPBase<Real>::lowerUnscaled(const SPxColId& id) const
{
   assert(id.isValid());
   return lowerUnscaled(number(id));
}

/// Changes objective vector to \p newObj.
template<>
void SPxLPBase<Real>::changeMaxObj(const VectorBase<Real>& newObj, bool scale)
{
   assert(maxObj().dim() == newObj.dim());

   if(scale)
   {
      assert(_isScaled);
      assert(lp_scaler);
      LPColSetBase<Real>::maxObj_w().scaleAssign(LPColSetBase<Real>::scaleExp.get_const_ptr(), newObj);
   }
   else
      LPColSetBase<Real>::maxObj_w() = newObj;

   assert(isConsistent());
}

/// Changes vector of lower bounds to \p newLower.
template<>
void SPxLPBase<Real>::changeLower(const VectorBase<Real>& newLower, bool scale)
{
   assert(lower().dim() == newLower.dim());

   if(scale)
   {
      assert(_isScaled);
      assert(lp_scaler);
      LPColSetBase<Real>::lower_w().scaleAssign(LPColSetBase<Real>::scaleExp.get_const_ptr(), newLower,
            true);
   }
   else
      LPColSetBase<Real>::lower_w() = newLower;

   assert(isConsistent());
}


/// Changes vector of upper bounds to \p newUpper.
template<>
void SPxLPBase<Real>::changeUpper(const VectorBase<Real>& newUpper, bool scale)
{
   assert(upper().dim() == newUpper.dim());

   if(scale)
   {
      assert(_isScaled);
      assert(lp_scaler);
      LPColSetBase<Real>::upper_w().scaleAssign(LPColSetBase<Real>::scaleExp.get_const_ptr(), newUpper,
            true);
   }
   else
      LPColSetBase<Real>::upper_w() = newUpper;

   assert(isConsistent());
}

/// Changes left hand side vector for constraints to \p newLhs.
template<>
void SPxLPBase<Real>::changeLhs(const VectorBase<Real>& newLhs, bool scale)
{
   assert(lhs().dim() == newLhs.dim());

   if(scale)
   {
      assert(_isScaled);
      assert(lp_scaler);
      LPRowSetBase<Real>::lhs_w().scaleAssign(LPRowSetBase<Real>::scaleExp.get_const_ptr(), newLhs);
   }
   else
      LPRowSetBase<Real>::lhs_w() = newLhs;

   assert(isConsistent());
}

/// Changes right hand side vector for constraints to \p newRhs.
template<>
void SPxLPBase<Real>::changeRhs(const VectorBase<Real>& newRhs, bool scale)
{
   assert(rhs().dim() == newRhs.dim());

   if(scale)
   {
      assert(_isScaled);
      assert(lp_scaler);
      LPRowSetBase<Real>::rhs_w().scaleAssign(LPRowSetBase<Real>::scaleExp.get_const_ptr(), newRhs);
   }
   else
      LPRowSetBase<Real>::rhs_w() = newRhs;

   assert(isConsistent());
}

/// Changes LP element (\p i, \p j) to \p val. \p scale determines whether the new data should be scaled
template <>
void SPxLPBase<Real>::changeElement(int i, int j, const Real& val, bool scale)
{
   if(i < 0 || j < 0)
      return;

   SVectorBase<Real>& row = rowVector_w(i);
   SVectorBase<Real>& col = colVector_w(j);

   if(isNotZero(val))
   {
      Real newVal;

      if(scale)
      {
         assert(_isScaled);
         assert(lp_scaler);
         newVal = lp_scaler->scaleElement(*this, i, j, val);
      }
      else
         newVal = val;

      if(row.pos(j) >= 0 && col.pos(i) >= 0)
      {
         row.value(row.pos(j)) = newVal;
         col.value(col.pos(i)) = newVal;
      }
      else
      {
         LPRowSetBase<Real>::add2(i, 1, &j, &newVal);
         LPColSetBase<Real>::add2(j, 1, &i, &newVal);
      }
   }
   else if(row.pos(j) >= 0 && col.pos(i) >= 0)
   {
      row.remove(row.pos(j));
      col.remove(col.pos(i));
   }

   assert(isConsistent());
}

// ---------------------------------------------------------------------------------------------------------------------
//  Specialization for reading LP format
// ---------------------------------------------------------------------------------------------------------------------

#define LPF_MAX_LINE_LEN  8192     ///< maximum length of a line (8190 + \\n + \\0)

/// Is \p c a \c space, \c tab, \c nl or \c cr ?
static inline bool LPFisSpace(int c)
{
   return (c == ' ') || (c == '\t') || (c == '\n') || (c == '\r');
}



/// Is there a number at the beginning of \p s ?
static bool LPFisValue(const char* s)
{
   return ((*s >= '0') && (*s <= '9')) || (*s == '+') || (*s == '-') || (*s == '.');
}



/// Is there a possible column name at the beginning of \p s ?
static bool LPFisColName(const char* s)
{
   // strchr() gives a true for the null char.
   if(*s == '\0')
      return false;

   return ((*s >= 'A') && (*s <= 'Z'))
          || ((*s >= 'a') && (*s <= 'z'))
          || (strchr("!\"#$%&()/,;?@_'`{}|~", *s) != 0);
}



/// Is there a comparison operator at the beginning of \p s ?
static bool LPFisSense(const char* s)
{
   return (*s == '<') || (*s == '>') || (*s == '=');
}



static bool LPFisInfinity(const char* s)
{
   return ((s[0] == '-') || (s[0] == '+'))
          && (tolower(s[1]) == 'i')
          && (tolower(s[2]) == 'n')
          && (tolower(s[3]) == 'f');
}



static bool LPFisFree(const char* s)
{
   return (tolower(s[0]) == 'f')
          && (tolower(s[1]) == 'r')
          && (tolower(s[2]) == 'e')
          && (tolower(s[3]) == 'e');
}



/// Read the next number and advance \p pos.
/** If only a sign is encountered, the number is assumed to be \c sign * 1.0.  This routine will not catch malformatted
 *  numbers like .e10 !
 */
static Real LPFreadValue(char*& pos, SPxOut* spxout)
{
   assert(LPFisValue(pos));

   char        tmp[LPF_MAX_LINE_LEN];
   const char* s = pos;
   char*       t;
   Real        value = 1.0;
   bool        has_digits = false;
   bool        has_emptyexponent = false;

   // 1. sign
   if((*s == '+') || (*s == '-'))
      s++;

   // 2. Digits before the decimal dot
   while((*s >= '0') && (*s <= '9'))
   {
      has_digits = true;
      s++;
   }

   // 3. Decimal dot
   if(*s == '.')
   {
      s++;

      // 4. If there was a dot, possible digit behind it
      while((*s >= '0') && (*s <= '9'))
      {
         has_digits = true;
         s++;
      }
   }

   // 5. Exponent
   if(tolower(*s) == 'e')
   {
      has_emptyexponent = true;
      s++;

      // 6. Exponent sign
      if((*s == '+') || (*s == '-'))
         s++;

      // 7. Exponent digits
      while((*s >= '0') && (*s <= '9'))
      {
         has_emptyexponent = false;
         s++;
      }
   }

   assert(s != pos);

   if(has_emptyexponent)
   {
      MSG_WARNING((*spxout), (*spxout) <<
                  "WLPFRD01 Warning: found empty exponent in LP file - check for forbidden variable names with initial 'e' or 'E'\n";
                 )
   }

   if(!has_digits)
      value = (*pos == '-') ? -1.0 : 1.0;
   else
   {
      for(t = tmp; pos != s; pos++)
         *t++ = *pos;

      *t = '\0';
      value = atof(tmp);
   }

   pos += s - pos;

   assert(pos == s);

   MSG_DEBUG(std::cout << "DLPFRD01 LPFreadValue = " << value << std::endl;)

   if(LPFisSpace(*pos))
      pos++;

   return value;
}



/// Read the next column name from the input.
/** The name read is looked up and if not found \p emptycol
 *  is added to \p colset. \p pos is advanced behind the name.
 *  @return The Index of the named column.
 */
static int LPFreadColName(char*& pos, NameSet* colnames, LPColSetBase<Real>& colset,
                          const LPColBase<Real>* emptycol, SPxOut* spxout)
{
   assert(LPFisColName(pos));
   assert(colnames != 0);

   char        name[LPF_MAX_LINE_LEN];
   const char* s = pos;
   int         i;
   int         colidx;

   // These are the characters that are not allowed in a column name.
   while((strchr("+-.<>= ", *s) == 0) && (*s != '\0'))
      s++;

   for(i = 0; pos != s; i++, pos++)
      name[i] = *pos;

   name[i] = '\0';

   if((colidx = colnames->number(name)) < 0)
   {
      // We only add the name if we got an empty column.
      if(emptycol == 0)
         MSG_WARNING((*spxout), (*spxout) << "WLPFRD02 Unknown variable \"" << name << "\" ";)
         else
         {
            colidx = colnames->num();
            colnames->add(name);
            colset.add(*emptycol);
         }
   }

   MSG_DEBUG(std::cout << "DLPFRD03 LPFreadColName [" << name << "] = " << colidx << std::endl;)

   if(LPFisSpace(*pos))
      pos++;

   return colidx;
}



/// Read the next <,>,=,==,<=,=<,>=,=> and advance \p pos.
static int LPFreadSense(char*& pos)
{
   assert(LPFisSense(pos));

   int sense = *pos++;

   if((*pos == '<') || (*pos == '>'))
      sense = *pos++;
   else if(*pos == '=')
      pos++;

   MSG_DEBUG(std::cout << "DLPFRD04 LPFreadSense = " << static_cast<char>(sense) << std::endl;)

   if(LPFisSpace(*pos))
      pos++;

   return sense;
}



/// Is the \p keyword present in \p buf ? If yes, advance \p pos.
/** \p keyword should be lower case. It can contain optional sections which are enclosed in '[' ']' like "min[imize]".
 */
static bool LPFhasKeyword(char*& pos, const char* keyword)
{
   int i;
   int k;

   assert(keyword != 0);

   for(i = 0, k = 0; keyword[i] != '\0'; i++, k++)
   {
      if(keyword[i] == '[')
      {
         i++;

         // Here we assumed that we have a ']' for the '['.
         while((tolower(pos[k]) == keyword[i]) && (pos[k] != '\0'))
         {
            k++;
            i++;
         }

         while(keyword[i] != ']')
            i++;

         --k;
      }
      else
      {
         if(keyword[i] != tolower(pos[k]))
            break;
      }
   }

   // we have to be at the end of the keyword and the word found on the line also has to end here.  Attention: The
   // LPFisSense is a kludge to allow LPFhasKeyword also to process Inf[inity] keywords in the bounds section.
   if(keyword[i] == '\0' && (pos[k] == '\0' || LPFisSpace(pos[k]) || LPFisSense(&pos[k])))
   {
      pos += k;

      MSG_DEBUG(std::cout << "DLPFRD05 LPFhasKeyword: " << keyword << std::endl;)

      return true;
   }

   return false;
}



/// If \p buf start with "name:" store the name in \p rownames and advance \p pos.
static bool LPFhasRowName(char*& pos, NameSet* rownames)
{
   const char* s = strchr(pos, ':');

   if(s == 0)
      return false;

   int dcolpos = int(s - pos);

   int end;
   int srt;

   // skip spaces between name and ":"
   for(end = dcolpos - 1; end >= 0; end--)
      if(pos[end] != ' ')
         break;

   // are there only spaces in front of the ":" ?
   if(end < 0)
   {
      pos = &(pos[dcolpos + 1]);
      return false;
   }

   // skip spaces in front of name
   for(srt = end - 1; srt >= 0; srt--)
      if(pos[srt] == ' ')
         break;

   // go back to the non-space character
   srt++;

   assert(srt <= end && pos[srt] != ' ');

   char name[LPF_MAX_LINE_LEN];
   int i;
   int k = 0;

   for(i = srt; i <= end; i++)
      name[k++] = pos[i];

   name[k] = '\0';

   if(rownames != 0)
      rownames->add(name);

   pos = &(pos[dcolpos + 1]);

   return true;
}



static Real LPFreadInfinity(char*& pos)
{
   assert(LPFisInfinity(pos));

   Real sense = (*pos == '-') ? -1.0 : 1.0;

   (void) LPFhasKeyword(++pos, "inf[inity]");

   return sense * infinity;
}



/// Read LP in "CPLEX LP File Format".
/** The specification is taken from the ILOG CPLEX 7.0 Reference Manual, Appendix E, Page 527.
 *
 *  This routine should read (most?) valid LP format files.  What it will not do, is find all cases where a file is ill
 *  formed.  If this happens it may complain and read nothing or read "something".
 *
 *  Problem: A line ending in '+' or '-' followed by a line starting with a number, will be regarded as an error.
 *
 *  The reader will accept the keyword INT[egers] as a synonym for GEN[erals] which is an undocumented feature in CPLEX.
 *
 *  A difference to the CPLEX reader, is that no name for the objective row is required.
 *
 * The manual says the maximum allowed line length is 255 characters, but CPLEX does not complain if the lines are
 * longer.
 *
 *  @return true if the file was read correctly
 */
template <>
bool SPxLPBase<Real>::readLPF(
   std::istream& p_input,                ///< input stream.
   NameSet*      p_rnames,               ///< row names.
   NameSet*      p_cnames,               ///< column names.
   DIdxSet*      p_intvars)              ///< integer variables.
{
   enum
   {
      START, OBJECTIVE, CONSTRAINTS, BOUNDS, INTEGERS, BINARIES
   } section = START;

   NameSet* rnames;                      ///< row names.
   NameSet* cnames;                      ///< column names.

   LPColSetBase<Real> cset;              ///< the set of columns read.
   LPRowSetBase<Real> rset;              ///< the set of rows read.
   LPColBase<Real> emptycol;             ///< reusable empty column.
   LPRowBase<Real> row;                  ///< last assembled row.
   DSVectorBase<Real> vec;               ///< last assembled vector (from row).

   Real val = 1.0;
   int colidx;
   int sense = 0;

   char buf[LPF_MAX_LINE_LEN];
   char tmp[LPF_MAX_LINE_LEN];
   char line[LPF_MAX_LINE_LEN];
   int lineno = 0;
   bool unnamed = true;
   bool finished = false;
   bool other;
   bool have_value = true;
   int i;
   int k;
   char* s;
   char* pos;
   char* pos_old = 0;

   if(p_cnames)
      cnames = p_cnames;
   else
   {
      cnames = 0;
      spx_alloc(cnames);
      cnames = new(cnames) NameSet();
   }

   cnames->clear();

   if(p_rnames)
      rnames = p_rnames;
   else
   {
      try
      {
         rnames = 0;
         spx_alloc(rnames);
         rnames = new(rnames) NameSet();
      }
      catch(const SPxMemoryException& x)
      {
         if(!p_cnames)
         {
            cnames->~NameSet();
            spx_free(cnames);
         }

         throw x;
      }
   }

   rnames->clear();

   SPxLPBase<Real>::clear(); // clear the LP.

   //--------------------------------------------------------------------------
   //--- Main Loop
   //--------------------------------------------------------------------------
   for(;;)
   {
      // 0. Read a line from the file.
      if(!p_input.getline(buf, sizeof(buf)))
      {
         if(strlen(buf) == LPF_MAX_LINE_LEN - 1)
         {
            MSG_ERROR(std::cerr << "ELPFRD06 Line exceeds " << LPF_MAX_LINE_LEN - 2
                      << " characters" << std::endl;)
         }
         else
         {
            MSG_ERROR(std::cerr << "ELPFRD07 No 'End' marker found" << std::endl;)
            finished = true;
         }

         break;
      }

      lineno++;
      i   = 0;
      pos = buf;

      MSG_DEBUG(std::cout << "DLPFRD08 Reading line " << lineno
                << " (pos=" << pos << ")" << std::endl;)

      // 1. Remove comments.
      if(0 != (s = strchr(buf, '\\')))
         * s = '\0';

      // 2. Look for keywords.
      if(section == START)
      {
         if(LPFhasKeyword(pos, "max[imize]"))
         {
            changeSense(SPxLPBase<Real>::MAXIMIZE);
            section = OBJECTIVE;
         }
         else if(LPFhasKeyword(pos, "min[imize]"))
         {
            changeSense(SPxLPBase<Real>::MINIMIZE);
            section = OBJECTIVE;
         }
      }
      else if(section == OBJECTIVE)
      {
         if(LPFhasKeyword(pos, "s[ubject][   ]t[o]")
               || LPFhasKeyword(pos, "s[uch][    ]t[hat]")
               || LPFhasKeyword(pos, "s[.][    ]t[.]")
               || LPFhasKeyword(pos, "lazy con[straints]"))
         {
            // store objective vector
            for(int j = vec.size() - 1; j >= 0; --j)
               cset.maxObj_w(vec.index(j)) = vec.value(j);

            // multiplication with -1 for minimization is done below
            vec.clear();
            have_value = true;
            val = 1.0;
            section = CONSTRAINTS;
         }
      }
      else if(section == CONSTRAINTS &&
              (LPFhasKeyword(pos, "s[ubject][   ]t[o]")
               || LPFhasKeyword(pos, "s[uch][    ]t[hat]")
               || LPFhasKeyword(pos, "s[.][    ]t[.]")))
      {
         have_value = true;
         val = 1.0;
      }
      else
      {
         if(LPFhasKeyword(pos, "lazy con[straints]"))
            ;
         else if(LPFhasKeyword(pos, "bound[s]"))
            section = BOUNDS;
         else if(LPFhasKeyword(pos, "bin[ary]"))
            section = BINARIES;
         else if(LPFhasKeyword(pos, "bin[aries]"))
            section = BINARIES;
         else if(LPFhasKeyword(pos, "gen[erals]"))
            section = INTEGERS;
         else if(LPFhasKeyword(pos, "int[egers]"))   // this is undocumented
            section = INTEGERS;
         else if(LPFhasKeyword(pos, "end"))
         {
            finished = true;
            break;
         }
         else if(LPFhasKeyword(pos, "s[ubject][   ]t[o]")  // second time
                 || LPFhasKeyword(pos, "s[uch][    ]t[hat]")
                 || LPFhasKeyword(pos, "s[.][    ]t[.]")
                 || LPFhasKeyword(pos, "lazy con[straints]"))
         {
            // In principle this has to checked for all keywords above,
            // otherwise we just ignore any half finished constraint
            if(have_value)
               goto syntax_error;

            have_value = true;
            val = 1.0;
         }
      }

      // 3a. Look for row names in objective and drop it.
      if(section == OBJECTIVE)
         LPFhasRowName(pos, 0);

      // 3b. Look for row name in constraint and store it.
      if(section == CONSTRAINTS)
         if(LPFhasRowName(pos, rnames))
            unnamed = false;

      // 4a. Remove initial spaces.
      while(LPFisSpace(pos[i]))
         i++;

      // 4b. remove spaces if they do not appear before the name of a vaiable.
      for(k = 0; pos[i] != '\0'; i++)
         if(!LPFisSpace(pos[i]) || LPFisColName(&pos[i + 1]))
            tmp[k++] = pos[i];

      tmp[k] = '\0';

      // 5. Is this an empty line ?
      if(tmp[0] == '\0')
         continue;

      // 6. Collapse sequences of '+' and '-'. e.g ++---+ => -
      for(i = 0, k = 0; tmp[i] != '\0'; i++)
      {
         while(((tmp[i] == '+') || (tmp[i] == '-')) && ((tmp[i + 1] == '+') || (tmp[i + 1] == '-')))
         {
            if(tmp[i++] == '-')
               tmp[i] = (tmp[i] == '-') ? '+' : '-';
         }

         line[k++] = tmp[i];
      }

      line[k] = '\0';

      //-----------------------------------------------------------------------
      //--- Line processing loop
      //-----------------------------------------------------------------------
      pos = line;

      MSG_DEBUG(std::cout << "DLPFRD09 pos=" << pos << std::endl;)

      // 7. We have something left to process.
      while((pos != 0) && (*pos != '\0'))
      {
         // remember our position, so we are sure we make progress.
         pos_old = pos;

         // now process the sections
         switch(section)
         {
         case OBJECTIVE:
            if(LPFisValue(pos))
            {
               Real pre_sign = 1.0;

               /* Already having here a value could only result from being the first number in a constraint, or a sign
                * '+' or '-' as last token on the previous line.
                */
               if(have_value)
               {
                  if(NE(spxAbs(val), 1.0))
                     goto syntax_error;

                  if(EQ(val, -1.0))
                     pre_sign = val;
               }

               /* non-finite coefficients are not allowed in the objective */
               if(LPFisInfinity(pos))
                  goto syntax_error;

               have_value = true;
               val = LPFreadValue(pos, spxout) * pre_sign;
            }

            if(*pos == '\0')
               continue;

            if(!have_value || !LPFisColName(pos))
               goto syntax_error;

            have_value = false;
            colidx = LPFreadColName(pos, cnames, cset, &emptycol, spxout);
            vec.add(colidx, val);
            break;

         case CONSTRAINTS:
            if(LPFisValue(pos))
            {
               Real pre_sign = 1.0;

               /* Already having here a value could only result from being the first number in a constraint, or a sign
                * '+' or '-' as last token on the previous line.
                */
               if(have_value)
               {
                  if(NE(spxAbs(val), 1.0))
                     goto syntax_error;

                  if(EQ(val, -1.0))
                     pre_sign = val;
               }

               if(LPFisInfinity(pos))
               {
                  /* non-finite coefficients are not allowed */
                  if(sense == 0)
                     goto syntax_error;

                  val = LPFreadInfinity(pos) * pre_sign;
               }
               else
                  val = LPFreadValue(pos, spxout) * pre_sign;

               have_value = true;

               if(sense != 0)
               {
                  if(sense == '<')
                  {
                     row.setLhs(-infinity);
                     row.setRhs(val);
                  }
                  else if(sense == '>')
                  {
                     row.setLhs(val);
                     row.setRhs(infinity);
                  }
                  else
                  {
                     assert(sense == '=');

                     row.setLhs(val);
                     row.setRhs(val);
                  }

                  row.setRowVector(vec);
                  rset.add(row);
                  vec.clear();

                  if(!unnamed)
                     unnamed = true;
                  else
                  {
                     char name[16];
                     spxSnprintf(name, 16, "C%d", rset.num());
                     rnames->add(name);
                  }

                  have_value = true;
                  val = 1.0;
                  sense = 0;
                  pos = 0;
                  // next line
                  continue;
               }
            }

            if(*pos == '\0')
               continue;

            if(have_value)
            {
               if(LPFisColName(pos))
               {
                  colidx = LPFreadColName(pos, cnames, cset, &emptycol, spxout);

                  if(val != 0.0)
                  {
                     // Do we have this index already in the row?
                     int n = vec.pos(colidx);

                     // if not, add it
                     if(n < 0)
                        vec.add(colidx, val);
                     // if yes, add them up and remove the element if it amounts to zero
                     else
                     {
                        assert(vec.index(n) == colidx);

                        val += vec.value(n);

                        if(val == 0.0)
                           vec.remove(n);
                        else
                           vec.value(n) = val;

                        assert(cnames->has(colidx));

                        MSG_WARNING((*spxout), (*spxout) << "WLPFRD10 Duplicate index "
                                    << (*cnames)[colidx]
                                    << " in line " << lineno
                                    << std::endl;)
                     }
                  }

                  have_value = false;
               }
               else
               {
                  // We have a row like c1: <= 5 with no variables. We can not handle 10 <= 5; issue a syntax error.
                  if(val != 1.0)
                     goto syntax_error;

                  // If the next thing is not the sense we give up also.
                  if(!LPFisSense(pos))
                     goto syntax_error;

                  have_value = false;
               }
            }

            assert(!have_value);

            if(LPFisSense(pos))
               sense = LPFreadSense(pos);

            break;

         case BOUNDS:
            other = false;
            sense = 0;

            if(LPFisValue(pos))
            {
               val = LPFisInfinity(pos) ? LPFreadInfinity(pos) : LPFreadValue(pos, spxout);

               if(!LPFisSense(pos))
                  goto syntax_error;

               sense = LPFreadSense(pos);
               other = true;
            }

            if(!LPFisColName(pos))
               goto syntax_error;

            if((colidx = LPFreadColName(pos, cnames, cset, 0, spxout)) < 0)
            {
               MSG_WARNING((*spxout), (*spxout) << "WLPFRD11 in Bounds section line "
                           << lineno << " ignored" << std::endl;)
               *pos = '\0';
               continue;
            }

            if(sense)
            {
               if(sense == '<')
                  cset.lower_w(colidx) = val;
               else if(sense == '>')
                  cset.upper_w(colidx) = val;
               else
               {
                  assert(sense == '=');
                  cset.lower_w(colidx) = val;
                  cset.upper_w(colidx) = val;
               }
            }

            if(LPFisFree(pos))
            {
               cset.lower_w(colidx) = -infinity;
               cset.upper_w(colidx) =  infinity;
               other = true;
               pos += 4;  // set position after the word "free"
            }
            else if(LPFisSense(pos))
            {
               sense = LPFreadSense(pos);
               other = true;

               if(!LPFisValue(pos))
                  goto syntax_error;

               val = LPFisInfinity(pos) ? LPFreadInfinity(pos) : LPFreadValue(pos, spxout);

               if(sense == '<')
                  cset.upper_w(colidx) = val;
               else if(sense == '>')
                  cset.lower_w(colidx) = val;
               else
               {
                  assert(sense == '=');
                  cset.lower_w(colidx) = val;
                  cset.upper_w(colidx) = val;
               }
            }

            /* Do we have only a single column name in the input line?  We could ignore this savely, but it is probably
             * a sign of some other error.
             */
            if(!other)
               goto syntax_error;

            break;

         case BINARIES:
         case INTEGERS:
            if((colidx = LPFreadColName(pos, cnames, cset, 0, spxout)) < 0)
            {
               MSG_WARNING((*spxout), (*spxout) << "WLPFRD12 in Binary/General section line " << lineno <<
                           " ignored" << std::endl;)
            }
            else
            {
               if(section == BINARIES)
               {
                  if(cset.lower(colidx) < 0.0)
                  {
                     cset.lower_w(colidx) = 0.0;
                  }

                  if(cset.upper(colidx) > 1.0)
                  {
                     cset.upper_w(colidx) = 1.0;
                  }
               }

               if(p_intvars != 0)
                  p_intvars->addIdx(colidx);
            }

            break;

         case START:
            MSG_ERROR(std::cerr << "ELPFRD13 This seems to be no LP format file" << std::endl;)
            goto syntax_error;

         default:
            throw SPxInternalCodeException("XLPFRD01 This should never happen.");
         }

         if(pos == pos_old)
            goto syntax_error;
      }
   }

   assert(isConsistent());

   addCols(cset);
   assert(isConsistent());

   addRows(rset);
   assert(isConsistent());

syntax_error:

   if(finished)
   {
      MSG_INFO2((*spxout), (*spxout) << "Finished reading " << lineno << " lines" << std::endl;)
   }
   else
      MSG_ERROR(std::cerr << "ELPFRD15 Syntax error in line " << lineno << std::endl;)

      if(p_cnames == 0)
         spx_free(cnames);

   if(p_rnames == 0)
      spx_free(rnames);

   return finished;
}



// ---------------------------------------------------------------------------------------------------------------------
// Specialization for reading MPS format
// ---------------------------------------------------------------------------------------------------------------------

/// Process NAME section.
static void MPSreadName(MPSInput& mps, SPxOut* spxout)
{
   do
   {
      // This has to be the Line with the NAME section.
      if(!mps.readLine() || (mps.field0() == 0) || strcmp(mps.field0(), "NAME"))
         break;

      // Sometimes the name is omitted.
      mps.setProbName((mps.field1() == 0) ? "_MPS_" : mps.field1());

      MSG_INFO2((*spxout), (*spxout) << "IMPSRD01 Problem name   : " << mps.probName() << std::endl;)

      // This has to be a new section
      if(!mps.readLine() || (mps.field0() == 0))
         break;

      if(!strcmp(mps.field0(), "ROWS"))
         mps.setSection(MPSInput::ROWS);
      else if(!strncmp(mps.field0(), "OBJSEN", 6))
         mps.setSection(MPSInput::OBJSEN);
      else if(!strcmp(mps.field0(), "OBJNAME"))
         mps.setSection(MPSInput::OBJNAME);
      else
         break;

      return;
   }
   while(false);

   mps.syntaxError();
}



/// Process OBJSEN section. This Section is an ILOG extension.
static void MPSreadObjsen(MPSInput& mps)
{
   do
   {
      // This has to be the Line with MIN or MAX.
      if(!mps.readLine() || (mps.field1() == 0))
         break;

      if(!strcmp(mps.field1(), "MIN"))
         mps.setObjSense(MPSInput::MINIMIZE);
      else if(!strcmp(mps.field1(), "MAX"))
         mps.setObjSense(MPSInput::MAXIMIZE);
      else
         break;

      // Look for ROWS or OBJNAME Section
      if(!mps.readLine() || (mps.field0() == 0))
         break;

      if(!strcmp(mps.field0(), "ROWS"))
         mps.setSection(MPSInput::ROWS);
      else if(!strcmp(mps.field0(), "OBJNAME"))
         mps.setSection(MPSInput::OBJNAME);
      else
         break;

      return;
   }
   while(false);

   mps.syntaxError();
}



/// Process OBJNAME section. This Section is an ILOG extension.
static void MPSreadObjname(MPSInput& mps)
{
   do
   {
      // This has to be the Line with the name.
      if(!mps.readLine() || (mps.field1() == 0))
         break;

      mps.setObjName(mps.field1());

      // Look for ROWS Section
      if(!mps.readLine() || (mps.field0() == 0))
         break;

      if(strcmp(mps.field0(), "ROWS"))
         break;

      mps.setSection(MPSInput::ROWS);

      return;
   }
   while(false);

   mps.syntaxError();
}



/// Process ROWS section.
static void MPSreadRows(MPSInput& mps, LPRowSetBase<Real>& rset, NameSet& rnames, SPxOut* spxout)
{
   LPRowBase<Real> row;

   while(mps.readLine())
   {
      if(mps.field0() != 0)
      {
         MSG_INFO2((*spxout), (*spxout) << "IMPSRD02 Objective name : " << mps.objName() << std::endl;)

         if(strcmp(mps.field0(), "COLUMNS"))
            break;

         mps.setSection(MPSInput::COLUMNS);

         return;
      }

      if((mps.field1() == 0) || (mps.field2() == 0))
         break;

      if(*mps.field1() == 'N')
      {
         if(*mps.objName() == '\0')
            mps.setObjName(mps.field2());
      }
      else
      {
         if(rnames.has(mps.field2()))
            break;

         rnames.add(mps.field2());

         switch(*mps.field1())
         {
         case 'G':
            row.setLhs(0.0);
            row.setRhs(infinity);
            break;

         case 'E':
            row.setLhs(0.0);
            row.setRhs(0.0);
            break;

         case 'L':
            row.setLhs(-infinity);
            row.setRhs(0.0);
            break;

         default:
            mps.syntaxError();
            return;
         }

         rset.add(row);
      }

      assert((*mps.field1() == 'N') || (rnames.number(mps.field2()) == rset.num() - 1));
   }

   mps.syntaxError();
}



/// Process COLUMNS section.
static void MPSreadCols(MPSInput& mps, const LPRowSetBase<Real>& rset, const NameSet&  rnames,
                        LPColSetBase<Real>& cset, NameSet& cnames, DIdxSet* intvars)
{
   Real val;
   int idx;
   char colname[MPSInput::MAX_LINE_LEN] = { '\0' };
   LPColBase<Real> col(rset.num());
   DSVectorBase<Real> vec;

   col.setObj(0.0);
   vec.clear();

   while(mps.readLine())
   {
      if(mps.field0() != 0)
      {
         if(strcmp(mps.field0(), "RHS"))
            break;

         if(colname[0] != '\0')
         {
            col.setColVector(vec);
            cset.add(col);
         }

         mps.setSection(MPSInput::RHS);

         return;
      }

      if((mps.field1() == 0) || (mps.field2() == 0) || (mps.field3() == 0))
         break;

      // new column?
      if(strcmp(colname, mps.field1()))
      {
         // first column?
         if(colname[0] != '\0')
         {
            col.setColVector(vec);
            cset.add(col);
         }

         // save copy of string (make sure string ends with \0)
         spxSnprintf(colname, MPSInput::MAX_LINE_LEN - 1, "%s", mps.field1());
         colname[MPSInput::MAX_LINE_LEN - 1] = '\0';

         int ncnames = cnames.size();
         cnames.add(colname);

         // check whether the new name is unique wrt previous column names
         if(cnames.size() <= ncnames)
         {
            MSG_ERROR(std::cerr << "ERROR in COLUMNS: duplicate column name or not column-wise ordering" <<
                      std::endl;)
            break;
         }

         vec.clear();
         col.setObj(0.0);
         col.setLower(0.0);
         col.setUpper(infinity);

         if(mps.isInteger())
         {
            assert(cnames.number(colname) == cset.num());

            if(intvars != 0)
               intvars->addIdx(cnames.number(colname));

            // for Integer variable the default bounds are 0/1
            col.setUpper(1.0);
         }
      }

      val = atof(mps.field3());

      if(!strcmp(mps.field2(), mps.objName()))
         col.setObj(val);
      else
      {
         if((idx = rnames.number(mps.field2())) < 0)
            mps.entryIgnored("Column", mps.field1(), "row", mps.field2());
         else if(val != 0.0)
            vec.add(idx, val);
      }

      if(mps.field5() != 0)
      {
         assert(mps.field4() != 0);

         val = atof(mps.field5());

         if(!strcmp(mps.field4(), mps.objName()))
            col.setObj(val);
         else
         {
            if((idx = rnames.number(mps.field4())) < 0)
               mps.entryIgnored("Column", mps.field1(), "row", mps.field4());
            else if(val != 0.0)
               vec.add(idx, val);
         }
      }
   }

   mps.syntaxError();
}



/// Process RHS section.
static void MPSreadRhs(MPSInput& mps, LPRowSetBase<Real>& rset, const NameSet& rnames,
                       SPxOut* spxout)
{
   char rhsname[MPSInput::MAX_LINE_LEN] = { '\0' };
   char addname[MPSInput::MAX_LINE_LEN] = { '\0' };
   int idx;
   Real val;

   while(mps.readLine())
   {
      if(mps.field0() != 0)
      {
         MSG_INFO2((*spxout), (*spxout) << "IMPSRD03 RHS name       : " << rhsname  << std::endl;);

         if(!strcmp(mps.field0(), "RANGES"))
            mps.setSection(MPSInput::RANGES);
         else if(!strcmp(mps.field0(), "BOUNDS"))
            mps.setSection(MPSInput::BOUNDS);
         else if(!strcmp(mps.field0(), "ENDATA"))
            mps.setSection(MPSInput::ENDATA);
         else
            break;

         return;
      }

      if(((mps.field2() != 0) && (mps.field3() == 0)) || ((mps.field4() != 0) && (mps.field5() == 0)))
         mps.insertName("_RHS_");

      if((mps.field1() == 0) || (mps.field2() == 0) || (mps.field3() == 0))
         break;

      if(*rhsname == '\0')
         spxSnprintf(rhsname, MPSInput::MAX_LINE_LEN, "%s", mps.field1());

      if(strcmp(rhsname, mps.field1()))
      {
         if(strcmp(addname, mps.field1()))
         {
            assert(strlen(mps.field1()) < MPSInput::MAX_LINE_LEN);
            spxSnprintf(addname, MPSInput::MAX_LINE_LEN, "%s", mps.field1());
            MSG_INFO3((*spxout), (*spxout) << "IMPSRD07 RHS ignored    : " << addname << std::endl);
         }
      }
      else
      {
         if((idx = rnames.number(mps.field2())) < 0)
            mps.entryIgnored("RHS", mps.field1(), "row", mps.field2());
         else
         {
            val = atof(mps.field3());

            // LE or EQ
            if(rset.rhs(idx) < infinity)
               rset.rhs_w(idx) = val;

            // GE or EQ
            if(rset.lhs(idx) > -infinity)
               rset.lhs_w(idx) = val;
         }

         if(mps.field5() != 0)
         {
            if((idx = rnames.number(mps.field4())) < 0)
               mps.entryIgnored("RHS", mps.field1(), "row", mps.field4());
            else
            {
               val = atof(mps.field5());

               // LE or EQ
               if(rset.rhs(idx) < infinity)
                  rset.rhs_w(idx) = val;

               // GE or EQ
               if(rset.lhs(idx) > -infinity)
                  rset.lhs_w(idx) = val;
            }
         }
      }
   }

   mps.syntaxError();
}



/// Process RANGES section.
static void MPSreadRanges(MPSInput& mps,  LPRowSetBase<Real>& rset, const NameSet& rnames,
                          SPxOut* spxout)
{
   char rngname[MPSInput::MAX_LINE_LEN] = { '\0' };
   int idx;
   Real val;

   while(mps.readLine())
   {
      if(mps.field0() != 0)
      {
         MSG_INFO2((*spxout), (*spxout) << "IMPSRD04 Range name     : " << rngname << std::endl;);

         if(!strcmp(mps.field0(), "BOUNDS"))
            mps.setSection(MPSInput::BOUNDS);
         else if(!strcmp(mps.field0(), "ENDATA"))
            mps.setSection(MPSInput::ENDATA);
         else
            break;

         return;
      }

      if(((mps.field2() != 0) && (mps.field3() == 0)) || ((mps.field4() != 0) && (mps.field5() == 0)))
         mps.insertName("_RNG_");

      if((mps.field1() == 0) || (mps.field2() == 0) || (mps.field3() == 0))
         break;

      if(*rngname == '\0')
      {
         assert(strlen(mps.field1()) < MPSInput::MAX_LINE_LEN);
         spxSnprintf(rngname, MPSInput::MAX_LINE_LEN, "%s", mps.field1());
      }

      /* The rules are:
       * Row Sign   LHS             RHS
       * ----------------------------------------
       *  G   +/-   rhs             rhs + |range|
       *  L   +/-   rhs - |range|   rhs
       *  E   +     rhs             rhs + range
       *  E   -     rhs + range     rhs
       * ----------------------------------------
       */
      if(!strcmp(rngname, mps.field1()))
      {
         if((idx = rnames.number(mps.field2())) < 0)
            mps.entryIgnored("Range", mps.field1(), "row", mps.field2());
         else
         {
            val = atof(mps.field3());

            // EQ
            if((rset.lhs(idx) > -infinity) && (rset.rhs_w(idx) <  infinity))
            {
               assert(rset.lhs(idx) == rset.rhs(idx));

               if(val >= 0)
                  rset.rhs_w(idx) += val;
               else
                  rset.lhs_w(idx) += val;
            }
            else
            {
               // GE
               if(rset.lhs(idx) > -infinity)
                  rset.rhs_w(idx)  = rset.lhs(idx) + spxAbs(val);
               // LE
               else
                  rset.lhs_w(idx)  = rset.rhs(idx) - spxAbs(val);
            }
         }

         if(mps.field5() != 0)
         {
            if((idx = rnames.number(mps.field4())) < 0)
               mps.entryIgnored("Range", mps.field1(), "row", mps.field4());
            else
            {
               val = atof(mps.field5());

               // EQ
               if((rset.lhs(idx) > -infinity) && (rset.rhs(idx) <  infinity))
               {
                  assert(rset.lhs(idx) == rset.rhs(idx));

                  if(val >= 0)
                     rset.rhs_w(idx) += val;
                  else
                     rset.lhs_w(idx) += val;
               }
               else
               {
                  // GE
                  if(rset.lhs(idx) > -infinity)
                     rset.rhs_w(idx)  = rset.lhs(idx) + spxAbs(val);
                  // LE
                  else
                     rset.lhs_w(idx)  = rset.rhs(idx) - spxAbs(val);
               }
            }
         }
      }
   }

   mps.syntaxError();
}



/// Process BOUNDS section.
static void MPSreadBounds(MPSInput& mps, LPColSetBase<Real>& cset, const NameSet& cnames,
                          DIdxSet* intvars, SPxOut* spxout)
{
   DIdxSet oldbinvars;
   char bndname[MPSInput::MAX_LINE_LEN] = { '\0' };
   int  idx;
   Real val;

   while(mps.readLine())
   {
      if(mps.field0() != 0)
      {
         MSG_INFO2((*spxout), (*spxout) << "IMPSRD05 Bound name     : " << bndname << std::endl;)

         if(strcmp(mps.field0(), "ENDATA"))
            break;

         mps.setSection(MPSInput::ENDATA);

         return;
      }

      // Is the value field used ?
      if((!strcmp(mps.field1(), "LO"))
            || (!strcmp(mps.field1(), "UP"))
            || (!strcmp(mps.field1(), "FX"))
            || (!strcmp(mps.field1(), "LI"))
            || (!strcmp(mps.field1(), "UI")))
      {
         if((mps.field3() != 0) && (mps.field4() == 0))
            mps.insertName("_BND_", true);
      }
      else
      {
         if((mps.field2() != 0) && (mps.field3() == 0))
            mps.insertName("_BND_", true);
      }

      if((mps.field1() == 0) || (mps.field2() == 0) || (mps.field3() == 0))
         break;

      if(*bndname == '\0')
      {
         assert(strlen(mps.field2()) < MPSInput::MAX_LINE_LEN);
         spxSnprintf(bndname, MPSInput::MAX_LINE_LEN, "%s", mps.field2());
      }

      // Only read the first Bound in section
      if(!strcmp(bndname, mps.field2()))
      {
         if((idx = cnames.number(mps.field3())) < 0)
            mps.entryIgnored("column", mps.field3(), "bound", bndname);
         else
         {
            if(mps.field4() == 0)
               val = 0.0;
            else if(!strcmp(mps.field4(), "-Inf") || !strcmp(mps.field4(), "-inf"))
               val = -infinity;
            else if(!strcmp(mps.field4(), "Inf") || !strcmp(mps.field4(), "inf")
                    || !strcmp(mps.field4(), "+Inf") || !strcmp(mps.field4(), "+inf"))
               val = infinity;
            else
               val = atof(mps.field4());

            // ILOG extension (Integer Bound)
            if(mps.field1()[1] == 'I')
            {
               if(intvars != 0)
                  intvars->addIdx(idx);

               // if the variable has appeared in the MARKER section of the COLUMNS section then its default bounds were
               // set to 0,1; the first time it is declared integer we need to change to default bounds 0,infinity
               if(oldbinvars.pos(idx) < 0)
               {
                  cset.upper_w(idx) = infinity;
                  oldbinvars.addIdx(idx);
               }
            }

            switch(*mps.field1())
            {
            case 'L':
               cset.lower_w(idx) = val;
               break;

            case 'U':
               cset.upper_w(idx) = val;
               break;

            case 'F':
               if(mps.field1()[1] == 'X')
               {
                  cset.lower_w(idx) = val;
                  cset.upper_w(idx) = val;
               }
               else
               {
                  cset.lower_w(idx) = -infinity;
                  cset.upper_w(idx) = infinity;
               }

               break;

            case 'M':
               cset.lower_w(idx) = -infinity;
               break;

            case 'P':
               cset.upper_w(idx) = infinity;
               break;

            // Ilog extension (Binary)
            case 'B':
               cset.lower_w(idx) = 0.0;
               cset.upper_w(idx) = 1.0;

               if(intvars != 0)
                  intvars->addIdx(idx);

               break;

            default:
               mps.syntaxError();
               return;
            }
         }
      }
   }

   mps.syntaxError();
}



/// Read LP in MPS File Format.
/**
 *  The specification is taken from the IBM Optimization Library Guide and Reference, online available at
 *  http://www.software.ibm.com/sos/features/libuser.htm and from the ILOG CPLEX 7.0 Reference Manual, Appendix E, Page
 *  531.
 *
 *  This routine should read all valid MPS format files.  What it will not do, is find all cases where a file is ill
 *  formed.  If this happens it may complain and read nothing or read "something".
 *
 *  @return true if the file was read correctly.
 */
#define INIT_COLS 10000 ///< initialy allocated columns.
#define INIT_NZOS 100000 ///< initialy allocated non zeros.
template <>
bool SPxLPBase<Real>::readMPS(
   std::istream& p_input,           ///< input stream.
   NameSet*      p_rnames,          ///< row names.
   NameSet*      p_cnames,          ///< column names.
   DIdxSet*      p_intvars)         ///< integer variables.
{
   LPRowSetBase<Real>& rset = *this;
   LPColSetBase<Real>& cset = *this;
   NameSet* rnames;
   NameSet* cnames;

   if(p_cnames)
      cnames = p_cnames;
   else
   {
      cnames = 0;
      spx_alloc(cnames);
      cnames = new(cnames) NameSet();
   }

   cnames->clear();

   if(p_rnames)
      rnames = p_rnames;
   else
   {
      try
      {
         rnames = 0;
         spx_alloc(rnames);
         rnames = new(rnames) NameSet();
      }
      catch(const SPxMemoryException& x)
      {
         if(!p_cnames)
         {
            cnames->~NameSet();
            spx_free(cnames);
         }

         throw x;
      }
   }

   rnames->clear();

   SPxLPBase<Real>::clear(); // clear the LP.

   cset.memRemax(INIT_NZOS);
   cset.reMax(INIT_COLS);

   MPSInput mps(p_input);

   MPSreadName(mps, spxout);

   if(mps.section() == MPSInput::OBJSEN)
      MPSreadObjsen(mps);

   if(mps.section() == MPSInput::OBJNAME)
      MPSreadObjname(mps);

   if(mps.section() == MPSInput::ROWS)
      MPSreadRows(mps, rset, *rnames, spxout);

   addedRows(rset.num());

   if(mps.section() == MPSInput::COLUMNS)
      MPSreadCols(mps, rset, *rnames, cset, *cnames, p_intvars);

   if(mps.section() == MPSInput::RHS)
      MPSreadRhs(mps, rset, *rnames, spxout);

   if(mps.section() == MPSInput::RANGES)
      MPSreadRanges(mps, rset, *rnames, spxout);

   if(mps.section() == MPSInput::BOUNDS)
      MPSreadBounds(mps, cset, *cnames, p_intvars, spxout);

   if(mps.section() != MPSInput::ENDATA)
      mps.syntaxError();

   if(mps.hasError())
      clear();
   else
   {
      changeSense(mps.objSense() == MPSInput::MINIMIZE ? SPxLPBase<Real>::MINIMIZE :
                  SPxLPBase<Real>::MAXIMIZE);

      MSG_INFO2((*spxout), (*spxout) << "IMPSRD06 Objective sense: " << ((mps.objSense() ==
                MPSInput::MINIMIZE) ? "Minimize\n" : "Maximize\n"));

      added2Set(
         *(reinterpret_cast<SVSetBase<Real>*>(static_cast<LPRowSetBase<Real>*>(this))),
         *(reinterpret_cast<SVSetBase<Real>*>(static_cast<LPColSetBase<Real>*>(this))),
         cset.num());
      addedCols(cset.num());

      assert(isConsistent());
   }

   if(p_cnames == 0)
   {
      cnames->~NameSet();
      spx_free(cnames);
   }

   if(p_rnames == 0)
   {
      rnames->~NameSet();
      spx_free(rnames);
   }

   return !mps.hasError();
}



// ---------------------------------------------------------------------------------------------------------------------
// Specialization for writing LP format
// ---------------------------------------------------------------------------------------------------------------------

// get the name of a row or construct one
static const char* LPFgetRowName(
   const SPxLPBase<Real>& p_lp,
   int                    p_idx,
   const NameSet*         p_rnames,
   char*                  p_buf,
   int                    p_num_written_rows
)
{
   assert(p_buf != 0);
   assert(p_idx >= 0);
   assert(p_idx <  p_lp.nRows());

   if(p_rnames != 0)
   {
      DataKey key = p_lp.rId(p_idx);

      if(p_rnames->has(key))
         return (*p_rnames)[key];
   }

   spxSnprintf(p_buf, 16, "C%d", p_num_written_rows);

   return p_buf;
}



// get the name of a column or construct one
static const char* getColName(
   const SPxLPBase<Real>& p_lp,
   int                    p_idx,
   const NameSet*         p_cnames,
   char*                  p_buf
)
{
   assert(p_buf != 0);
   assert(p_idx >= 0);
   assert(p_idx <  p_lp.nCols());

   if(p_cnames != 0)
   {
      DataKey key = p_lp.cId(p_idx);

      if(p_cnames->has(key))
         return (*p_cnames)[key];
   }

   spxSnprintf(p_buf, 16, "x%d", p_idx);

   return p_buf;
}



// write an SVector
#define NUM_ENTRIES_PER_LINE 5
static void LPFwriteSVector(
   const SPxLPBase<Real>&   p_lp,       ///< the LP
   std::ostream&            p_output,   ///< output stream
   const NameSet*           p_cnames,   ///< column names
   const SVectorBase<Real>& p_svec)     ///< vector to write
{

   char name[16];
   int num_coeffs = 0;

   for(int j = 0; j < p_lp.nCols(); ++j)
   {
      const Real coeff = p_svec[j];

      if(coeff == 0)
         continue;

      if(num_coeffs == 0)
         p_output << coeff << " " << getColName(p_lp, j, p_cnames, name);
      else
      {
         // insert a line break every NUM_ENTRIES_PER_LINE columns
         if(num_coeffs % NUM_ENTRIES_PER_LINE == 0)
            p_output << "\n\t";

         if(coeff < 0)
            p_output << " - " << -coeff;
         else
            p_output << " + " << coeff;

         p_output << " " << getColName(p_lp, j, p_cnames, name);
      }

      ++num_coeffs;
   }
}



// write the objective
static void LPFwriteObjective(
   const SPxLPBase<Real>& p_lp,       ///< the LP
   std::ostream&          p_output,   ///< output stream
   const NameSet*         p_cnames    ///< column names
)
{

   const int sense = p_lp.spxSense();

   p_output << ((sense == SPxLPBase<Real>::MINIMIZE) ? "Minimize\n" : "Maximize\n");
   p_output << "  obj: ";

   const VectorBase<Real>& obj = p_lp.maxObj();
   DSVectorBase<Real> svec(obj.dim());
   svec.operator = (obj);
   svec *= Real(sense);
   LPFwriteSVector(p_lp, p_output, p_cnames, svec);
   p_output << "\n";
}



// write non-ranged rows
static void LPFwriteRow(
   const SPxLPBase<Real>&   p_lp,       ///< the LP
   std::ostream&            p_output,   ///< output stream
   const NameSet*           p_cnames,   ///< column names
   const SVectorBase<Real>& p_svec,     ///< vector of the row
   const Real&              p_lhs,      ///< lhs of the row
   const Real&              p_rhs       ///< rhs of the row
)
{

   LPFwriteSVector(p_lp, p_output, p_cnames, p_svec);

   if(p_lhs == p_rhs)
      p_output << " = " << p_rhs;
   else if(p_lhs <= -infinity)
      p_output << " <= " << p_rhs;
   else
   {
      assert(p_rhs >= infinity);
      p_output << " >= " << p_lhs;
   }

   p_output << "\n";
}



// write all rows
static void LPFwriteRows(
   const SPxLPBase<Real>& p_lp,       ///< the LP
   std::ostream&          p_output,   ///< output stream
   const NameSet*         p_rnames,   ///< row names
   const NameSet*         p_cnames   ///< column names
)
{

   char name[16];

   p_output << "Subject To\n";

   for(int i = 0; i < p_lp.nRows(); ++i)
   {
      const Real lhs = p_lp.lhs(i);
      const Real rhs = p_lp.rhs(i);

      if(lhs > -infinity && rhs < infinity && lhs != rhs)
      {
         // ranged row -> write two non-ranged rows
         p_output << " " << LPFgetRowName(p_lp, i, p_rnames, name, i) << "_1 : ";
         LPFwriteRow(p_lp, p_output, p_cnames, p_lp.rowVector(i), lhs, infinity);

         p_output << " " << LPFgetRowName(p_lp, i, p_rnames, name, i) << "_2 : ";
         LPFwriteRow(p_lp, p_output, p_cnames, p_lp.rowVector(i), -infinity, rhs);
      }
      else
      {
         p_output << " " << LPFgetRowName(p_lp, i, p_rnames, name, i) << " : ";
         LPFwriteRow(p_lp, p_output, p_cnames, p_lp.rowVector(i), lhs, rhs);
      }
   }
}



// write the variable bounds
// (the default bounds 0 <= x <= infinity are not written)
static void LPFwriteBounds(
   const SPxLPBase<Real>&   p_lp,       ///< the LP to write
   std::ostream&            p_output,   ///< output stream
   const NameSet*           p_cnames    ///< column names
)
{

   char name[16];

   p_output << "Bounds\n";

   for(int j = 0; j < p_lp.nCols(); ++j)
   {
      const Real lower = p_lp.lower(j);
      const Real upper = p_lp.upper(j);

      if(lower == upper)
      {
         p_output << "  "   << getColName(p_lp, j, p_cnames, name) << " = "  << upper << '\n';
      }
      else if(lower > -infinity)
      {
         if(upper < infinity)
         {
            // range bound
            if(lower != 0)
               p_output << "  "   << lower << " <= "
                        << getColName(p_lp, j, p_cnames, name)
                        << " <= " << upper << '\n';
            else
               p_output << "  "   << getColName(p_lp, j, p_cnames, name)
                        << " <= " << upper << '\n';
         }
         else if(lower != 0)
            p_output << "  " << lower << " <= "
                     << getColName(p_lp, j, p_cnames, name)
                     << '\n';
      }
      else if(upper < infinity)
         p_output << "   -Inf <= "
                  << getColName(p_lp, j, p_cnames, name)
                  << " <= " << upper << '\n';
      else
         p_output << "  "   << getColName(p_lp, j, p_cnames, name)
                  << " free\n";
   }
}



// write the generals section
static void LPFwriteGenerals(
   const SPxLPBase<Real>&   p_lp,         ///< the LP to write
   std::ostream&            p_output,     ///< output stream
   const NameSet*           p_cnames,     ///< column names
   const DIdxSet*           p_intvars     ///< integer variables
)
{

   char name[16];

   if(p_intvars == NULL || p_intvars->size() <= 0)
      return;  // no integer variables

   p_output << "Generals\n";

   for(int j = 0; j < p_lp.nCols(); ++j)
      if(p_intvars->pos(j) >= 0)
         p_output << "  " << getColName(p_lp, j, p_cnames, name) << "\n";
}



/// Write LP in LP Format.
template <>
void SPxLPBase<Real>::writeLPF(
   std::ostream&  p_output,          ///< output stream
   const NameSet* p_rnames,          ///< row names
   const NameSet* p_cnames,          ///< column names
   const DIdxSet* p_intvars          ///< integer variables
) const
{
   SPxOut::setScientific(p_output, 16);

   LPFwriteObjective(*this, p_output, p_cnames);
   LPFwriteRows(*this, p_output, p_rnames, p_cnames);
   LPFwriteBounds(*this, p_output, p_cnames);
   LPFwriteGenerals(*this, p_output, p_cnames, p_intvars);

   p_output << "End" << std::endl;
}



// ---------------------------------------------------------------------------------------------------------------------
// Specialization for writing MPS format
// ---------------------------------------------------------------------------------------------------------------------

static void MPSwriteRecord(
   std::ostream&  os,
   const char*    indicator,
   const char*    name,
   const char*    name1  = 0,
   const Real     value1 = 0.0,
   const char*    name2  = 0,
   const Real     value2 = 0.0
)
{
   char buf[81];

   spxSnprintf(buf, sizeof(buf), " %-2.2s %-8.8s", (indicator == 0) ? "" : indicator,
               (name == 0)      ? "" : name);
   os << buf;

   if(name1 != 0)
   {
      spxSnprintf(buf, sizeof(buf), "%-8.8s  %.15" REAL_FORMAT, name1, value1);
      os << buf;

      if(name2 != 0)
      {
         spxSnprintf(buf, sizeof(buf), "   %-8.8s  %.15" REAL_FORMAT, name2, value2);
         os << buf;
      }
   }

   os << std::endl;
}



static Real MPSgetRHS(Real left, Real right)
{
   Real rhsval;

   if(left > -infinity)   /// This includes ranges
      rhsval = left;
   else if(right <  infinity)
      rhsval = right;
   else
      throw SPxInternalCodeException("XMPSWR01 This should never happen.");

   return rhsval;
}



static const char* MPSgetRowName(
   const SPxLPBase<Real>& lp,
   int                   idx,
   const NameSet*        rnames,
   char*                 buf
)
{
   assert(buf != 0);
   assert(idx >= 0);
   assert(idx <  lp.nRows());

   if(rnames != 0)
   {
      DataKey key = lp.rId(idx);

      if(rnames->has(key))
         return (*rnames)[key];
   }

   spxSnprintf(buf, 16, "C%d", idx);

   return buf;
}



/// Write LP in MPS format.
/** @note There will always be a BOUNDS section, even if there are no bounds.
 */
template <>
void SPxLPBase<Real>::writeMPS(
   std::ostream&  p_output,          ///< output stream.
   const NameSet* p_rnames,          ///< row names.
   const NameSet* p_cnames,          ///< column names.
   const DIdxSet* p_intvars          ///< integer variables.
) const
{

   const char*    indicator;
   char           name [16];
   char           name1[16];
   char           name2[16];
   bool           has_ranges = false;
   int            i;
   int            k;

   SPxOut::setScientific(p_output, 16);
   // --- NAME Section ---
   p_output << "NAME          MPSDATA" << std::endl;

   // --- ROWS Section ---
   p_output << "ROWS" << std::endl;

   for(i = 0; i < nRows(); i++)
   {
      if(lhs(i) == rhs(i))
         indicator = "E";
      else if((lhs(i) > -infinity) && (rhs(i) < infinity))
      {
         indicator = "E";
         has_ranges = true;
      }
      else if(lhs(i) > -infinity)
         indicator = "G";
      else if(rhs(i) <  infinity)
         indicator = "L";
      else
         throw SPxInternalCodeException("XMPSWR02 This should never happen.");

      MPSwriteRecord(p_output, indicator, MPSgetRowName(*this, i, p_rnames, name));
   }

   MPSwriteRecord(p_output, "N", "MINIMIZE");

   // --- COLUMNS Section ---
   p_output << "COLUMNS" << std::endl;

   bool has_intvars = (p_intvars != 0) && (p_intvars->size() > 0);

   for(int j = 0; j < (has_intvars ? 2 : 1); j++)
   {
      bool is_intrun = has_intvars && (j == 1);

      if(is_intrun)
         p_output << "    MARK0001  'MARKER'                 'INTORG'" << std::endl;

      for(i = 0; i < nCols(); i++)
      {
         bool is_intvar = has_intvars && (p_intvars->pos(i) >= 0);

         if((is_intrun && !is_intvar) || (!is_intrun &&  is_intvar))
            continue;

         const SVectorBase<Real>& col = colVector(i);
         int colsize2 = (col.size() / 2) * 2;

         assert(colsize2 % 2 == 0);

         for(k = 0; k < colsize2; k += 2)
            MPSwriteRecord(p_output, 0, getColName(*this, i, p_cnames, name),
                           MPSgetRowName(*this, col.index(k), p_rnames, name1), col.value(k),
                           MPSgetRowName(*this, col.index(k + 1), p_rnames, name2), col.value(k + 1));

         if(colsize2 != col.size())
            MPSwriteRecord(p_output, 0, getColName(*this, i, p_cnames, name),
                           MPSgetRowName(*this, col.index(k), p_rnames, name1), col.value(k));

         if(isNotZero(maxObj(i)))
            MPSwriteRecord(p_output, 0, getColName(*this, i, p_cnames, name), "MINIMIZE", -maxObj(i));
      }

      if(is_intrun)
         p_output << "    MARK0001  'MARKER'                 'INTEND'" << std::endl;
   }

   // --- RHS Section ---
   p_output << "RHS" << std::endl;

   i = 0;

   while(i < nRows())
   {
      Real rhsval1 = 0.0;
      Real rhsval2 = 0.0;

      for(; i < nRows(); i++)
         if((rhsval1 = MPSgetRHS(lhs(i), rhs(i))) != 0.0)
            break;

      if(i < nRows())
      {
         for(k = i + 1; k < nRows(); k++)
         {
            if((rhsval2 = MPSgetRHS(lhs(k), rhs(k))) != 0.0)
               break;
         }

         if(k < nRows())
         {
            MPSwriteRecord(p_output, 0, "RHS", MPSgetRowName(*this, i, p_rnames, name1), rhsval1,
                           MPSgetRowName(*this, k, p_rnames, name2), rhsval2);
         }
         else
            MPSwriteRecord(p_output, 0, "RHS", MPSgetRowName(*this, i, p_rnames, name1), rhsval1);

         i = k + 1;
      }
   }

   // --- RANGES Section ---
   if(has_ranges)
   {
      p_output << "RANGES" << std::endl;

      for(i = 0; i < nRows(); i++)
      {
         if((lhs(i) > -infinity) && (rhs(i) < infinity))
            MPSwriteRecord(p_output, "", "RANGE", MPSgetRowName(*this, i, p_rnames, name1), rhs(i) - lhs(i));
      }
   }

   // --- BOUNDS Section ---
   p_output << "BOUNDS" << std::endl;

   for(i = 0; i < nCols(); i++)
   {
      // skip variables that do not appear in the objective function or any constraint
      const SVectorBase<Real>& col = colVector(i);

      if(col.size() == 0 && isZero(maxObj(i)))
         continue;

      if(lower(i) == upper(i))
      {
         MPSwriteRecord(p_output, "FX", "BOUND", getColName(*this, i, p_cnames, name1), lower(i));
         continue;
      }

      if((lower(i) <= -infinity) && (upper(i) >= infinity))
      {
         MPSwriteRecord(p_output, "FR", "BOUND", getColName(*this, i, p_cnames, name1));
         continue;
      }

      if(lower(i) != 0.0)
      {
         if(lower(i) > -infinity)
            MPSwriteRecord(p_output, "LO", "BOUND", getColName(*this, i, p_cnames, name1), lower(i));
         else
            MPSwriteRecord(p_output, "MI", "BOUND", getColName(*this, i, p_cnames, name1));
      }

      if(has_intvars && (p_intvars->pos(i) >= 0))
      {
         // Integer variables have default upper bound 1.0, but we should write
         // it nevertheless since CPLEX seems to assume infinity otherwise.
         MPSwriteRecord(p_output, "UP", "BOUND", getColName(*this, i, p_cnames, name1), upper(i));
      }
      else
      {
         // Continous variables have default upper bound infinity
         if(upper(i) < infinity)
            MPSwriteRecord(p_output, "UP", "BOUND", getColName(*this, i, p_cnames, name1), upper(i));
      }
   }

   // --- ENDATA Section ---
   p_output << "ENDATA" << std::endl;

   // Output warning when writing a maximisation problem
   if(spxSense() == SPxLPBase<Real>::MAXIMIZE)
   {
      MSG_WARNING((*spxout), (*spxout) <<
                  "XMPSWR03 Warning: objective function inverted when writing maximization problem in MPS file format\n");
   }
}



/// Building the dual problem from a given LP
/// @note primalRows must be as large as the number of unranged primal rows + 2 * the number of ranged primal rows.
///       dualCols must have the identical size to the primal rows.
template < >
void SPxLPBase<Real>::buildDualProblem(SPxLPBase<Real>& dualLP, SPxRowId primalRowIds[],
                                       SPxColId primalColIds[],
                                       SPxRowId dualRowIds[], SPxColId dualColIds[], int* nprimalrows, int* nprimalcols, int* ndualrows,
                                       int* ndualcols)
{
   // Setting up the primalrowids and dualcolids arrays if not given as parameters
   if(primalRowIds == 0 || primalColIds == 0 || dualRowIds == 0 || dualColIds == 0)
   {
      DataArray < SPxRowId > primalrowids(2 * nRows());
      DataArray < SPxColId > primalcolids(2 * nCols());
      DataArray < SPxRowId > dualrowids(2 * nCols());
      DataArray < SPxColId > dualcolids(2 * nRows());
      int numprimalrows = 0;
      int numprimalcols = 0;
      int numdualrows = 0;
      int numdualcols = 0;

      buildDualProblem(dualLP, primalrowids.get_ptr(), primalcolids.get_ptr(), dualrowids.get_ptr(),
                       dualcolids.get_ptr(), &numprimalrows, &numprimalcols, &numdualrows, &numdualcols);

      if(primalRowIds != 0)
      {
         primalRowIds = primalrowids.get_ptr();
         (*nprimalrows) = numprimalrows;
      }

      if(primalColIds != 0)
      {
         primalColIds = primalcolids.get_ptr();
         (*nprimalcols) = numprimalcols;
      }

      if(dualRowIds != 0)
      {
         dualRowIds = dualrowids.get_ptr();
         (*ndualrows) = numdualrows;
      }

      if(dualColIds != 0)
      {
         dualColIds = dualcolids.get_ptr();
         (*ndualcols) = numdualcols;
      }

      return;
   }

   // setting the sense of the dual LP
   if(spxSense() == MINIMIZE)
      dualLP.changeSense(MAXIMIZE);
   else
      dualLP.changeSense(MINIMIZE);

   LPRowSetBase<Real> dualrows(nCols());
   LPColSetBase<Real> dualcols(nRows());
   DSVectorBase<Real> col(1);

   int numAddedRows = 0;
   int numVarBoundCols = 0;
   int primalrowsidx = 0;
   int primalcolsidx = 0;

   for(int i = 0; i < nCols(); ++i)
   {
      primalColIds[primalcolsidx] = cId(i);
      primalcolsidx++;

      if(lower(i) <= -infinity && upper(i) >= infinity)   // unrestricted variables
      {
         dualrows.create(0, obj(i), obj(i));
         numAddedRows++;
      }
      else if(lower(i) <= -infinity)   // no lower bound is set, indicating a <= 0 variable
      {
         if(isZero(upper(i)))   // standard bound variable
         {
            if(spxSense() == MINIMIZE)
               dualrows.create(0, obj(i), infinity);
            else
               dualrows.create(0, -infinity, obj(i));
         }
         else // additional upper bound on the variable
         {
            col.add(numAddedRows, 1.0);

            if(spxSense() == MINIMIZE)
            {
               dualrows.create(0, obj(i), obj(i));
               dualcols.add(upper(i), -infinity, col, 0.0);
            }
            else
            {
               dualrows.create(0, obj(i), obj(i));
               dualcols.add(upper(i), 0.0, col, infinity);
            }

            col.clear();

            numVarBoundCols++;
         }

         numAddedRows++;
      }
      else if(upper(i) >= infinity)   // no upper bound set, indicating a >= 0 variable
      {
         if(isZero(lower(i)))   // standard bound variable
         {
            if(spxSense() == MINIMIZE)
               dualrows.create(0, -infinity, obj(i));
            else
               dualrows.create(0, obj(i), infinity);
         }
         else // additional lower bound on the variable
         {
            col.add(numAddedRows, 1.0);

            if(spxSense() == MINIMIZE)
            {
               dualrows.create(0, obj(i), obj(i));
               dualcols.add(lower(i), 0.0, col, infinity);
            }
            else
            {
               dualrows.create(0, obj(i), obj(i));
               dualcols.add(lower(i), -infinity, col, 0.0);
            }

            col.clear();

            numVarBoundCols++;
         }

         numAddedRows++;
      }
      else if(NE(lower(i), upper(i)))   // a boxed variable
      {
         if(isZero(lower(i)))   // variable bounded between 0 and upper(i)
         {
            col.add(numAddedRows, 1.0);

            if(spxSense() == MINIMIZE)
            {
               dualrows.create(0, -infinity, obj(i));
               dualcols.add(upper(i), -infinity, col, 0.0);
            }
            else
            {
               dualrows.create(0, obj(i), infinity);
               dualcols.add(upper(i), 0.0, col, infinity);
            }

            col.clear();

            numVarBoundCols++;
         }
         else if(isZero(upper(i)))   // variable bounded between lower(i) and 0
         {
            col.add(numAddedRows, 1.0);

            if(spxSense() == MINIMIZE)
            {
               dualrows.create(0, obj(i), infinity);
               dualcols.add(lower(i), 0.0, col, infinity);
            }
            else
            {
               dualrows.create(0, -infinity, obj(i));
               dualcols.add(lower(i), -infinity, col, 0.0);
            }

            col.clear();

            numVarBoundCols++;
         }
         else // variable bounded between lower(i) and upper(i)
         {
            dualrows.create(0, obj(i), obj(i));

            col.add(numAddedRows, 1.0);

            if(spxSense() == MINIMIZE)
            {
               dualcols.add(lower(i), 0.0, col, infinity);
               dualcols.add(upper(i), -infinity, col, 0.0);
            }
            else
            {
               dualcols.add(lower(i), -infinity, col, 0.0);
               dualcols.add(upper(i), 0.0, col, infinity);
            }

            col.clear();

            numVarBoundCols += 2;
         }

         numAddedRows++;
      }
      else
      {
         assert(lower(i) == upper(i));

         dualrows.create(0, obj(i), obj(i));

         col.add(numAddedRows, 1.0);
         dualcols.add(lower(i), 0, col, infinity);
         dualcols.add(lower(i), -infinity, col, 0);
         col.clear();

         numVarBoundCols += 2;
         numAddedRows++;
      }
   }

   // adding the empty rows to the dual LP
   dualLP.addRows(dualrows);

   // setting the dual row ids for the related primal cols.
   // this assumes that the rows are added in sequential order.
   for(int i = 0; i < primalcolsidx; i++)
      dualRowIds[i] = dualLP.rId(i);

   (*nprimalcols) = primalcolsidx;
   (*ndualrows) = primalcolsidx;

   // iterating over each of the rows to create dual columns
   for(int i = 0; i < nRows(); ++i)
   {
      // checking the type of the row
      switch(rowType(i))
      {
      case LPRowBase<Real>::RANGE: // range constraint, requires the addition of two dual variables
         assert(lhs(i) > -infinity);
         assert(rhs(i) < infinity);

         if(spxSense() == MINIMIZE)
         {
            primalRowIds[primalrowsidx] = rId(i); // setting the rowid for the primal row
            primalrowsidx++;
            dualcols.add(lhs(i), 0.0, rowVector(i), infinity);

            primalRowIds[primalrowsidx] = rId(i); // setting the rowid for the primal row
            primalrowsidx++;
            dualcols.add(rhs(i), -infinity, rowVector(i), 0.0);
         }
         else
         {
            primalRowIds[primalrowsidx] = rId(i); // setting the rowid for the primal row
            primalrowsidx++;
            dualcols.add(lhs(i), -infinity, rowVector(i), 0.0);

            primalRowIds[primalrowsidx] = rId(i); // setting the rowid for the primal row
            primalrowsidx++;
            dualcols.add(rhs(i), 0.0, rowVector(i), infinity);
         }

         break;

      case LPRowBase<Real>::GREATER_EQUAL: // >= constraint
         assert(lhs(i) > -infinity);
         primalRowIds[primalrowsidx] = rId(i); // setting the rowid for the primal row
         primalrowsidx++;

         if(spxSense() == MINIMIZE)
            dualcols.add(lhs(i), 0.0, rowVector(i), infinity);
         else
            dualcols.add(lhs(i), -infinity, rowVector(i), 0.0);

         break;

      case LPRowBase<Real>::LESS_EQUAL: // <= constriant
         assert(rhs(i) < infinity);
         primalRowIds[primalrowsidx] = rId(i); // setting the rowid for the primal row
         primalrowsidx++;

         if(spxSense() == MINIMIZE)
            dualcols.add(rhs(i), -infinity, rowVector(i), 0.0);
         else
            dualcols.add(rhs(i), 0.0, rowVector(i), infinity);

         break;

      case LPRowBase<Real>::EQUAL: // Equality constraint
         assert(EQ(lhs(i), rhs(i)));
         primalRowIds[primalrowsidx] = rId(i); // setting the rowid for the primal row
         primalrowsidx++;
         dualcols.add(rhs(i), -infinity, rowVector(i), infinity);
         break;

      default:
         throw SPxInternalCodeException("XLPFRD01 This should never happen.");
      }
   }

   // adding the filled columns to the dual LP
   dualLP.addCols(dualcols);

   // setting the dual column ids for the related primal rows.
   // this assumes that the columns are added in sequential order.
   for(int i = 0; i < primalrowsidx; i++)
      dualColIds[i] = dualLP.cId(i + numVarBoundCols);

   (*nprimalrows) = primalrowsidx;
   (*ndualcols) = primalrowsidx;
}



// ---------------------------------------------------------------------------------------------------------------------
//  Explicit instantiation
// ---------------------------------------------------------------------------------------------------------------------

template class SPxLPBase < Real >;

} // namespace soplex