-                      r458
+                      r459
         $(GLPK_CFLAGS) \
         $(CPLEX_CFLAGS) \
+        $(SOPLEX_CXXFLAGS)
+        $(SOPLEX_CXXFLAGS) \
+        $(CLP_CXXFLAGS)
 lemon_libemon_la_LDFLAGS = \
         $(GLPK_LIBS) \
         $(CPLEX_LIBS) \
+        $(SOPLEX_LIBS)
+        $(SOPLEX_LIBS) \
+        $(CLP_LIBS)
 if HAVE_GLPK
 lemon_libemon_la_SOURCES += lemon/lp_glpk.cc lemon/mip_glpk.cc
+lemon_libemon_la_SOURCES += lemon/lp_glpk.cc
 endif
 if HAVE_CPLEX
 lemon_libemon_la_SOURCES += lemon/lp_cplex.cc lemon/mip_cplex.cc
+lemon_libemon_la_SOURCES += lemon/lp_cplex.cc
 endif
 if HAVE_SOPLEX
 lemon_libemon_la_SOURCES += lemon/lp_soplex.cc
+endif
+if HAVE_CLP
+lemon_libemon_la_SOURCES += lemon/lp_clp.cc
 endif
 …
         lemon/lp.h \
         lemon/lp_base.h \
+        lemon/lp_clp.h \
         lemon/lp_cplex.h \
         lemon/lp_glpk.h \
         lemon/lp_skeleton.h \
         lemon/lp_soplex.h \
+        lemon/mip_cplex.h \
+        lemon/mip_glpk.h \
+        lemon/list_graph.h \
         lemon/maps.h \
         lemon/math.h \
 …
         lemon/bits/graph_adaptor_extender.h \
         lemon/bits/graph_extender.h \
-        lemon/bits/lp_id.h \
         lemon/bits/map_extender.h \
         lemon/bits/path_dump.h \
+        lemon/bits/solver_bits.h \
         lemon/bits/traits.h \
         lemon/bits/variant.h \

lemon/config.h.in

r458	r459
13	13	/* Define to 1 if you have SOPLEX */
14	14	#undef HAVE_SOPLEX
	15
	16	/* Define to 1 if you have CLP */
	17	#undef HAVE_CLP

lemon/lp.h

-                      r458
+                      r459
 #ifdef HAVE_GLPK
 #include <lemon/lp_glpk.h>
-#include <lemon/mip_glpk.h>
 #elif HAVE_CPLEX
 #include <lemon/lp_cplex.h>
-#include <lemon/mip_cplex.h>
 #elif HAVE_SOPLEX
 #include <lemon/lp_soplex.h>
+#elif HAVE_CLP
+#include <lemon/lp_clp.h>
 #endif
 …
   ///\ingroup lp_group
   ///
+  ///Currently, the possible values are \c GLPK or \c CPLEX
+#define DEFAULT_LP SOLVER
+  ///Currently, the possible values are \c LP_GLPK, \c LP_CPLEX, \c
+  ///LP_SOPLEX or \c LP_CLP
+#define LEMON_DEFAULT_LP SOLVER
   ///The default LP solver
 …
   ///\ingroup lp_group
   ///
   ///Currently, it is either \c LpGlpk or \c LpCplex
+  ///Currently, it is either \c LpGlpk, \c LpCplex, \c LpSoplex or \c LpClp
   typedef LpGlpk Lp;
-  ///The default LP solver identifier string
+  ///The default LP solver identifier string.
+  ///The default MIP solver identifier
+  ///The default MIP solver identifier.
   ///\ingroup lp_group
   ///
+  ///Currently, the possible values are "GLPK" or "CPLEX"
+  const char default_solver_name[]="SOLVER";
+  ///Currently, the possible values are \c MIP_GLPK or \c MIP_CPLEX
+#define LEMON_DEFAULT_MIP SOLVER
+  ///The default MIP solver.
+  ///The default ILP solver.
+  ///The default ILP solver.
+  ///The default MIP solver.
   ///\ingroup lp_group
   ///
   ///Currently, it is either \c LpGlpk or \c LpCplex
+  ///Currently, it is either \c MipGlpk or \c MipCplex
   typedef MipGlpk Mip;
 #else
 #ifdef HAVE_GLPK
 #define DEFAULT_LP GLPK
+# define LEMON_DEFAULT_LP LP_GLPK
   typedef LpGlpk Lp;
+# define LEMON_DEFAULT_MIP MIP_GLPK
   typedef MipGlpk Mip;
-  const char default_solver_name[]="GLPK";
 #elif HAVE_CPLEX
 #define DEFAULT_LP CPLEX
+# define LEMON_DEFAULT_LP LP_CPLEX
   typedef LpCplex Lp;
+# define LEMON_DEFAULT_MIP MIP_CPLEX
   typedef MipCplex Mip;
-  const char default_solver_name[]="CPLEX";
 #elif HAVE_SOPLEX
 #define DEFAULT_LP SOPLEX
+# define DEFAULT_LP LP_SOPLEX
   typedef LpSoplex Lp;
+  const char default_solver_name[]="SOPLEX";
+#elif HAVE_CLP
+# define DEFAULT_LP LP_CLP
+  typedef LpClp Lp;
 #endif
 #endif

lemon/lp_base.cc

-                      r458
+                      r459
 namespace lemon {
+  const LpSolverBase::Value
+  LpSolverBase::INF = std::numeric_limits<Value>::infinity();
+  const LpSolverBase::Value
+  LpSolverBase::NaN = std::numeric_limits<Value>::quiet_NaN();
+//   const LpSolverBase::Constr::Value
+//   LpSolverBase::Constr::INF = std::numeric_limits<Value>::infinity();
+//   const LpSolverBase::Constr::Value
+//   LpSolverBase::Constr::NaN = std::numeric_limits<Value>::quiet_NaN();
+  const LpBase::Value LpBase::INF = std::numeric_limits<Value>::infinity();
+  const LpBase::Value LpBase::NaN = std::numeric_limits<Value>::quiet_NaN();
 } //namespace lemon

lemon/lp_base.h

-                      r458
+                      r459
 #include<lemon/math.h>
+#include<lemon/error.h>
+#include<lemon/assert.h>
 #include<lemon/core.h>
 #include<lemon/bits/lp_id.h>
+#include<lemon/bits/solver_bits.h>
 ///\file
 …
 namespace lemon {
+  /// Function to decide whether a floating point value is finite or not.
+  /// Retruns true if the argument is not infinity, minus infinity or NaN.
+  /// It does the same as the isfinite() function defined by C99.
+  template <typename T>
+  bool isFinite(T value)
+  {
+    typedef std::numeric_limits<T> Lim;
+    if ((Lim::has_infinity && (value == Lim::infinity() || value ==
+                               -Lim::infinity())) ||
+        ((Lim::has_quiet_NaN || Lim::has_signaling_NaN) && value != value))
+    {
+      return false;
+    }
+    return true;
+  }
+  ///Common base class for LP solvers
+  ///\todo Much more docs
+  ///Common base class for LP and MIP solvers
+  ///Usually this class is not used directly, please use one of the concrete
+  ///implementations of the solver interface.
   ///\ingroup lp_group
   class LpSolverBase {
+  class LpBase {
   protected:
     _lp_bits::LpId rows;
     _lp_bits::LpId cols;
+    _solver_bits::VarIndex rows;
+    _solver_bits::VarIndex cols;
   public:
 …
     };
+      ///\e
+    enum SolutionStatus {
+      ///Feasible solution hasn't been found (but may exist).
+      ///\todo NOTFOUND might be a better name.
+      ///
+      UNDEFINED = 0,
+      ///The problem has no feasible solution
+      INFEASIBLE = 1,
+      ///Feasible solution found
+      FEASIBLE = 2,
+      ///Optimal solution exists and found
+      OPTIMAL = 3,
+      ///The cost function is unbounded
+      ///\todo Give a feasible solution and an infinite ray (and the
+      ///corresponding bases)
+      INFINITE = 4
+    };
+    ///\e The type of the investigated LP problem
+    enum ProblemTypes {
+      ///Primal-dual feasible
+      PRIMAL_DUAL_FEASIBLE = 0,
+      ///Primal feasible dual infeasible
+      PRIMAL_FEASIBLE_DUAL_INFEASIBLE = 1,
+      ///Primal infeasible dual feasible
+      PRIMAL_INFEASIBLE_DUAL_FEASIBLE = 2,
+      ///Primal-dual infeasible
+      PRIMAL_DUAL_INFEASIBLE = 3,
+      ///Could not determine so far
+      UNKNOWN = 4
+    ///Direction of the optimization
+    enum Sense {
+      /// Minimization
+      MIN,
+      /// Maximization
+      MAX
     };
 …
     static const Value NaN;
-    static inline bool isNaN(const Value& v) { return v!=v; }
     friend class Col;
     friend class ColIt;
     friend class Row;
+    friend class RowIt;
     ///Refer to a column of the LP.
 …
     ///other columns.
     ///
     ///\todo Document what can one do with a Col (INVALID, comparing,
     ///it is similar to Node/Edge)
+    ///\note This class is similar to other Item types in LEMON, like
+    ///Node and Arc types in digraph.
     class Col {
+      friend class LpBase;
     protected:
+      int id;
+      friend class LpSolverBase;
+      friend class MipSolverBase;
+      explicit Col(int _id) : id(_id) {}
+      int _id;
+      explicit Col(int id) : _id(id) {}
     public:
       typedef Value ExprValue;
+      typedef True LpSolverCol;
+      typedef True LpCol;
+      /// Default constructor
+      /// \warning The default constructor sets the Col to an
+      /// undefined value.
       Col() {}
+      Col(const Invalid&) : id(-1) {}
+      bool operator< (Col c) const  {return id< c.id;}
+      bool operator> (Col c) const  {return id> c.id;}
+      bool operator==(Col c) const  {return id==c.id;}
+      bool operator!=(Col c) const  {return id!=c.id;}
+      /// Invalid constructor \& conversion.
+      /// This constructor initializes the Col to be invalid.
+      /// \sa Invalid for more details.
+      Col(const Invalid&) : _id(-1) {}
+      /// Equality operator
+      /// Two \ref Col "Col"s are equal if and only if they point to
+      /// the same LP column or both are invalid.
+      bool operator==(Col c) const  {return _id == c._id;}
+      /// Inequality operator
+      /// \sa operator==(Col c)
+      ///
+      bool operator!=(Col c) const  {return _id != c._id;}
+      /// Artificial ordering operator.
+      /// To allow the use of this object in std::map or similar
+      /// associative container we require this.
+      ///
+      /// \note This operator only have to define some strict ordering of
+      /// the items; this order has nothing to do with the iteration
+      /// ordering of the items.
+      bool operator<(Col c) const  {return _id < c._id;}
     };
+    ///Iterator for iterate over the columns of an LP problem
+    /// Its usage is quite simple, for example you can count the number
+    /// of columns in an LP \c lp:
+    ///\code
+    /// int count=0;
+    /// for (LpBase::ColIt c(lp); c!=INVALID; ++c) ++count;
+    ///\endcode
     class ColIt : public Col {
       const LpSolverBase *_lp;
+      const LpBase *_solver;
     public:
+      /// Default constructor
+      /// \warning The default constructor sets the iterator
+      /// to an undefined value.
       ColIt() {}
+      ColIt(const LpSolverBase &lp) : _lp(&lp)
+      /// Sets the iterator to the first Col
+      /// Sets the iterator to the first Col.
+      ///
+      ColIt(const LpBase &solver) : _solver(&solver)
+      {
+        _lp->cols.firstFix(id);
+      }
+        _solver->cols.firstItem(_id);
+      }
+      /// Invalid constructor \& conversion
+      /// Initialize the iterator to be invalid.
+      /// \sa Invalid for more details.
       ColIt(const Invalid&) : Col(INVALID) {}
+      /// Next column
+      /// Assign the iterator to the next column.
+      ///
       ColIt &operator++()
+      {
         _lp->cols.nextFix(id);
+        _solver->cols.nextItem(_id);
         return *this;
+      }
     };
+    static int id(const Col& col) { return col.id; }
+    /// \brief Returns the ID of the column.
+    static int id(const Col& col) { return col._id; }
+    /// \brief Returns the column with the given ID.
+    ///
+    /// \pre The argument should be a valid column ID in the LP problem.
+    static Col colFromId(int id) { return Col(id); }
     ///Refer to a row of the LP.
 …
     ///other rows.
     ///
     ///\todo Document what can one do with a Row (INVALID, comparing,
     ///it is similar to Node/Edge)
+    ///\note This class is similar to other Item types in LEMON, like
+    ///Node and Arc types in digraph.
     class Row {
+      friend class LpBase;
     protected:
+      int id;
+      friend class LpSolverBase;
+      explicit Row(int _id) : id(_id) {}
+      int _id;
+      explicit Row(int id) : _id(id) {}
     public:
       typedef Value ExprValue;
+      typedef True LpSolverRow;
+      typedef True LpRow;
+      /// Default constructor
+      /// \warning The default constructor sets the Row to an
+      /// undefined value.
       Row() {}
+      Row(const Invalid&) : id(-1) {}
+      bool operator< (Row c) const  {return id< c.id;}
+      bool operator> (Row c) const  {return id> c.id;}
+      bool operator==(Row c) const  {return id==c.id;}
+      bool operator!=(Row c) const  {return id!=c.id;}
+      /// Invalid constructor \& conversion.
+      /// This constructor initializes the Row to be invalid.
+      /// \sa Invalid for more details.
+      Row(const Invalid&) : _id(-1) {}
+      /// Equality operator
+      /// Two \ref Row "Row"s are equal if and only if they point to
+      /// the same LP row or both are invalid.
+      bool operator==(Row r) const  {return _id == r._id;}
+      /// Inequality operator
+      /// \sa operator==(Row r)
+      ///
+      bool operator!=(Row r) const  {return _id != r._id;}
+      /// Artificial ordering operator.
+      /// To allow the use of this object in std::map or similar
+      /// associative container we require this.
+      ///
+      /// \note This operator only have to define some strict ordering of
+      /// the items; this order has nothing to do with the iteration
+      /// ordering of the items.
+      bool operator<(Row r) const  {return _id < r._id;}
     };
+    ///Iterator for iterate over the rows of an LP problem
+    /// Its usage is quite simple, for example you can count the number
+    /// of rows in an LP \c lp:
+    ///\code
+    /// int count=0;
+    /// for (LpBase::RowIt c(lp); c!=INVALID; ++c) ++count;
+    ///\endcode
     class RowIt : public Row {
       const LpSolverBase *_lp;
+      const LpBase *_solver;
     public:
+      /// Default constructor
+      /// \warning The default constructor sets the iterator
+      /// to an undefined value.
       RowIt() {}
+      RowIt(const LpSolverBase &lp) : _lp(&lp)
+      /// Sets the iterator to the first Row
+      /// Sets the iterator to the first Row.
+      ///
+      RowIt(const LpBase &solver) : _solver(&solver)
+      {
+        _lp->rows.firstFix(id);
+      }
+        _solver->rows.firstItem(_id);
+      }
+      /// Invalid constructor \& conversion
+      /// Initialize the iterator to be invalid.
+      /// \sa Invalid for more details.
       RowIt(const Invalid&) : Row(INVALID) {}
+      /// Next row
+      /// Assign the iterator to the next row.
+      ///
       RowIt &operator++()
+      {
         _lp->rows.nextFix(id);
+        _solver->rows.nextItem(_id);
         return *this;
+      }
     };
+    static int id(const Row& row) { return row.id; }
+  protected:
+    int _lpId(const Col& c) const {
+      return cols.floatingId(id(c));
+    }
+    int _lpId(const Row& r) const {
+      return rows.floatingId(id(r));
+    }
+    Col _item(int i, Col) const {
+      return Col(cols.fixId(i));
+    }
+    Row _item(int i, Row) const {
+      return Row(rows.fixId(i));
+    }
+    /// \brief Returns the ID of the row.
+    static int id(const Row& row) { return row._id; }
+    /// \brief Returns the row with the given ID.
+    ///
+    /// \pre The argument should be a valid row ID in the LP problem.
+    static Row rowFromId(int id) { return Row(id); }
   public:
 …
     ///There are several ways to access and modify the contents of this
     ///container.
-    ///- Its it fully compatible with \c std::map<Col,double>, so for expamle
-    ///if \c e is an Expr and \c v and \c w are of type \ref Col, then you can
-    ///read and modify the coefficients like
-    ///these.
     ///\code
     ///e[v]=5;
 …
     ///\code
     ///double s=0;
     ///for(LpSolverBase::Expr::iterator i=e.begin();i!=e.end();++i)
     ///  s+=i->second;
+    ///for(LpBase::Expr::ConstCoeffIt i(e);i!=INVALID;++i)
+    ///  s+=*i * primal(i);
     ///\endcode
     ///(This code computes the sum of all coefficients).
+    ///(This code computes the primal value of the expression).
     ///- Numbers (<tt>double</tt>'s)
     ///and variables (\ref Col "Col"s) directly convert to an
 …
     ///v*2.1+(3*v+(v*12+w+6)*3)/2
     ///\endcode
     ///are valid \ref Expr "Expr"essions.
+    ///are valid expressions.
     ///The usual assignment operations are also defined.
     ///\code
 …
     ///e/=5;
     ///\endcode
+    ///- The constant member can be set and read by \ref constComp()
+    ///- The constant member can be set and read by dereference
+    ///  operator (unary *)
+    ///
     ///\code
     ///e.constComp()=12;
     ///double c=e.constComp();
+    ///*e=12;
+    ///double c=*e;
     ///\endcode
     ///
-    ///\note \ref clear() not only sets all coefficients to 0 but also
-    ///clears the constant components.
-    ///
     ///\sa Constr
+    ///
+    class Expr : public std::map<Col,Value>
+    {
+    class Expr {
+      friend class LpBase;
     public:
+      typedef LpSolverBase::Col Key;
+      typedef LpSolverBase::Value Value;
+      /// The key type of the expression
+      typedef LpBase::Col Key;
+      /// The value type of the expression
+      typedef LpBase::Value Value;
     protected:
-      typedef std::map<Col,Value> Base;
       Value const_comp;
+      std::map<int, Value> comps;
     public:
+      typedef True IsLinExpression;
+      ///\e
+      Expr() : Base(), const_comp(0) { }
+      ///\e
+      Expr(const Key &v) : const_comp(0) {
+        Base::insert(std::make_pair(v, 1));
+      }
+      ///\e
+      typedef True SolverExpr;
+      /// Default constructor
+      /// Construct an empty expression, the coefficients and
+      /// the constant component are initialized to zero.
+      Expr() : const_comp(0) {}
+      /// Construct an expression from a column
+      /// Construct an expression, which has a term with \c c variable
+      /// and 1.0 coefficient.
+      Expr(const Col &c) : const_comp(0) {
+        typedef std::map<int, Value>::value_type pair_type;
+        comps.insert(pair_type(id(c), 1));
+      }
+      /// Construct an expression from a constant
+      /// Construct an expression, which's constant component is \c v.
+      ///
       Expr(const Value &v) : const_comp(v) {}
+      ///\e
+      void set(const Key &v,const Value &c) {
+        Base::insert(std::make_pair(v, c));
+      }
+      ///\e
+      Value &constComp() { return const_comp; }
+      ///\e
+      const Value &constComp() const { return const_comp; }
+      ///Removes the components with zero coefficient.
+      void simplify() {
+        for (Base::iterator i=Base::begin(); i!=Base::end();) {
+          Base::iterator j=i;
+          ++j;
+          if ((*i).second==0) Base::erase(i);
+          i=j;
+      /// Returns the coefficient of the column
+      Value operator[](const Col& c) const {
+        std::map<int, Value>::const_iterator it=comps.find(id(c));
+        if (it != comps.end()) {
+          return it->second;
+        } else {
+          return 0;
+        }
+      }
+      void simplify() const {
+        const_cast<Expr*>(this)->simplify();
+      }
+      ///Removes the coefficients closer to zero than \c tolerance.
+      void simplify(double &tolerance) {
+        for (Base::iterator i=Base::begin(); i!=Base::end();) {
+          Base::iterator j=i;
+          ++j;
+          if (std::fabs((*i).second)<tolerance) Base::erase(i);
+          i=j;
+      /// Returns the coefficient of the column
+      Value& operator[](const Col& c) {
+        return comps[id(c)];
+      }
+      /// Sets the coefficient of the column
+      void set(const Col &c, const Value &v) {
+        if (v != 0.0) {
+          typedef std::map<int, Value>::value_type pair_type;
+          comps.insert(pair_type(id(c), v));
+        } else {
+          comps.erase(id(c));
+        }
+      }
+      /// Returns the constant component of the expression
+      Value& operator*() { return const_comp; }
+      /// Returns the constant component of the expression
+      const Value& operator*() const { return const_comp; }
+      /// \brief Removes the coefficients which's absolute value does
+      /// not exceed \c epsilon. It also sets to zero the constant
+      /// component, if it does not exceed epsilon in absolute value.
+      void simplify(Value epsilon = 0.0) {
+        std::map<int, Value>::iterator it=comps.begin();
+        while (it != comps.end()) {
+          std::map<int, Value>::iterator jt=it;
+          ++jt;
+          if (std::fabs((*it).second) <= epsilon) comps.erase(it);
+          it=jt;
+        }
+        if (std::fabs(const_comp) <= epsilon) const_comp = 0;
+      }
+      void simplify(Value epsilon = 0.0) const {
+        const_cast<Expr*>(this)->simplify(epsilon);
+      }
       ///Sets all coefficients and the constant component to 0.
       void clear() {
         Base::clear();
+        comps.clear();
         const_comp=0;
+      }
       ///\e
+      ///Compound assignment
       Expr &operator+=(const Expr &e) {
+        for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
+          (*this)[j->first]+=j->second;
+        for (std::map<int, Value>::const_iterator it=e.comps.begin();
+             it!=e.comps.end(); ++it)
+          comps[it->first]+=it->second;
         const_comp+=e.const_comp;
         return *this;
+      }
       ///\e
+      ///Compound assignment
       Expr &operator-=(const Expr &e) {
+        for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
+          (*this)[j->first]-=j->second;
+        for (std::map<int, Value>::const_iterator it=e.comps.begin();
+             it!=e.comps.end(); ++it)
+          comps[it->first]-=it->second;
         const_comp-=e.const_comp;
         return *this;
+      }
+      ///\e
+      Expr &operator*=(const Value &c) {
+        for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
+          j->second*=c;
+        const_comp*=c;
+      ///Multiply with a constant
+      Expr &operator*=(const Value &v) {
+        for (std::map<int, Value>::iterator it=comps.begin();
+             it!=comps.end(); ++it)
+          it->second*=v;
+        const_comp*=v;
         return *this;
+      }
       ///\e
+      ///Division with a constant
       Expr &operator/=(const Value &c) {
+        for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
+          j->second/=c;
+        for (std::map<int, Value>::iterator it=comps.begin();
+             it!=comps.end(); ++it)
+          it->second/=c;
         const_comp/=c;
         return *this;
+      }
+      ///Iterator over the expression
+      ///The iterator iterates over the terms of the expression.
+      ///
+      ///\code
+      ///double s=0;
+      ///for(LpBase::Expr::CoeffIt i(e);i!=INVALID;++i)
+      ///  s+= *i * primal(i);
+      ///\endcode
+      class CoeffIt {
+      private:
+        std::map<int, Value>::iterator _it, _end;
+      public:
+        /// Sets the iterator to the first term
+        /// Sets the iterator to the first term of the expression.
+        ///
+        CoeffIt(Expr& e)
+          : _it(e.comps.begin()), _end(e.comps.end()){}
+        /// Convert the iterator to the column of the term
+        operator Col() const {
+          return colFromId(_it->first);
+        }
+        /// Returns the coefficient of the term
+        Value& operator*() { return _it->second; }
+        /// Returns the coefficient of the term
+        const Value& operator*() const { return _it->second; }
+        /// Next term
+        /// Assign the iterator to the next term.
+        ///
+        CoeffIt& operator++() { ++_it; return *this; }
+        /// Equality operator
+        bool operator==(Invalid) const { return _it == _end; }
+        /// Inequality operator
+        bool operator!=(Invalid) const { return _it != _end; }
+      };
+      /// Const iterator over the expression
+      ///The iterator iterates over the terms of the expression.
+      ///
+      ///\code
+      ///double s=0;
+      ///for(LpBase::Expr::ConstCoeffIt i(e);i!=INVALID;++i)
+      ///  s+=*i * primal(i);
+      ///\endcode
+      class ConstCoeffIt {
+      private:
+        std::map<int, Value>::const_iterator _it, _end;
+      public:
+        /// Sets the iterator to the first term
+        /// Sets the iterator to the first term of the expression.
+        ///
+        ConstCoeffIt(const Expr& e)
+          : _it(e.comps.begin()), _end(e.comps.end()){}
+        /// Convert the iterator to the column of the term
+        operator Col() const {
+          return colFromId(_it->first);
+        }
+        /// Returns the coefficient of the term
+        const Value& operator*() const { return _it->second; }
+        /// Next term
+        /// Assign the iterator to the next term.
+        ///
+        ConstCoeffIt& operator++() { ++_it; return *this; }
+        /// Equality operator
+        bool operator==(Invalid) const { return _it == _end; }
+        /// Inequality operator
+        bool operator!=(Invalid) const { return _it != _end; }
+      };
     };
 …
     ///  e>=t
     ///\endcode
     ///\warning The validity of a constraint is checked only at run time, so
     ///e.g. \ref addRow(<tt>x[1]\<=x[2]<=5</tt>) will compile, but will throw
     ///an assertion.
+    ///\warning The validity of a constraint is checked only at run
+    ///time, so e.g. \ref addRow(<tt>x[1]\<=x[2]<=5</tt>) will
+    ///compile, but will fail an assertion.
     class Constr
+    {
     public:
       typedef LpSolverBase::Expr Expr;
+      typedef LpBase::Expr Expr;
       typedef Expr::Key Key;
       typedef Expr::Value Value;
 …
       Constr() : _expr(), _lb(NaN), _ub(NaN) {}
       ///\e
       Constr(Value lb,const Expr &e,Value ub) :
+      Constr(Value lb, const Expr &e, Value ub) :
         _expr(e), _lb(lb), _ub(ub) {}
-      ///\e
-      Constr(const Expr &e,Value ub) :
-        _expr(e), _lb(NaN), _ub(ub) {}
-      ///\e
-      Constr(Value lb,const Expr &e) :
-        _expr(e), _lb(lb), _ub(NaN) {}
-      ///\e
       Constr(const Expr &e) :
         _expr(e), _lb(NaN), _ub(NaN) {}
 …
       ///Is the constraint lower bounded?
       bool lowerBounded() const {
         return isFinite(_lb);
+        return _lb != -INF && !std::isnan(_lb);
+      }
       ///Is the constraint upper bounded?
       bool upperBounded() const {
         return isFinite(_ub);
+        return _ub != INF && !std::isnan(_ub);
+      }
 …
     ///There are several ways to access and modify the contents of this
     ///container.
-    ///- Its it fully compatible with \c std::map<Row,double>, so for expamle
-    ///if \c e is an DualExpr and \c v
-    ///and \c w are of type \ref Row, then you can
-    ///read and modify the coefficients like
-    ///these.
     ///\code
     ///e[v]=5;
 …
     ///\code
     ///double s=0;
     ///for(LpSolverBase::DualExpr::iterator i=e.begin();i!=e.end();++i)
     ///  s+=i->second;
+    ///for(LpBase::DualExpr::ConstCoeffIt i(e);i!=INVALID;++i)
+    ///  s+=*i;
     ///\endcode
     ///(This code computes the sum of all coefficients).
 …
     ///v*2.1+(3*v+(v*12+w)*3)/2
     ///\endcode
     ///are valid \ref DualExpr "DualExpr"essions.
+    ///are valid \ref DualExpr dual expressions.
     ///The usual assignment operations are also defined.
     ///\code
 …
     ///
     ///\sa Expr
+    ///
+    class DualExpr : public std::map<Row,Value>
+    {
+    class DualExpr {
+      friend class LpBase;
     public:
+      typedef LpSolverBase::Row Key;
+      typedef LpSolverBase::Value Value;
+      /// The key type of the expression
+      typedef LpBase::Row Key;
+      /// The value type of the expression
+      typedef LpBase::Value Value;
     protected:
       typedef std::map<Row,Value> Base;
+      std::map<int, Value> comps;
     public:
+      typedef True IsLinExpression;
+      ///\e
+      DualExpr() : Base() { }
+      ///\e
+      DualExpr(const Key &v) {
+        Base::insert(std::make_pair(v, 1));
+      }
+      ///\e
+      void set(const Key &v,const Value &c) {
+        Base::insert(std::make_pair(v, c));
+      }
+      ///Removes the components with zero coefficient.
+      void simplify() {
+        for (Base::iterator i=Base::begin(); i!=Base::end();) {
+          Base::iterator j=i;
+          ++j;
+          if ((*i).second==0) Base::erase(i);
+          i=j;
+      typedef True SolverExpr;
+      /// Default constructor
+      /// Construct an empty expression, the coefficients are
+      /// initialized to zero.
+      DualExpr() {}
+      /// Construct an expression from a row
+      /// Construct an expression, which has a term with \c r dual
+      /// variable and 1.0 coefficient.
+      DualExpr(const Row &r) {
+        typedef std::map<int, Value>::value_type pair_type;
+        comps.insert(pair_type(id(r), 1));
+      }
+      /// Returns the coefficient of the row
+      Value operator[](const Row& r) const {
+        std::map<int, Value>::const_iterator it = comps.find(id(r));
+        if (it != comps.end()) {
+          return it->second;
+        } else {
+          return 0;
+        }
+      }
+      void simplify() const {
+        const_cast<DualExpr*>(this)->simplify();
+      }
+      ///Removes the coefficients closer to zero than \c tolerance.
+      void simplify(double &tolerance) {
+        for (Base::iterator i=Base::begin(); i!=Base::end();) {
+          Base::iterator j=i;
+          ++j;
+          if (std::fabs((*i).second)<tolerance) Base::erase(i);
+          i=j;
+      /// Returns the coefficient of the row
+      Value& operator[](const Row& r) {
+        return comps[id(r)];
+      }
+      /// Sets the coefficient of the row
+      void set(const Row &r, const Value &v) {
+        if (v != 0.0) {
+          typedef std::map<int, Value>::value_type pair_type;
+          comps.insert(pair_type(id(r), v));
+        } else {
+          comps.erase(id(r));
+        }
+      }
+      /// \brief Removes the coefficients which's absolute value does
+      /// not exceed \c epsilon.
+      void simplify(Value epsilon = 0.0) {
+        std::map<int, Value>::iterator it=comps.begin();
+        while (it != comps.end()) {
+          std::map<int, Value>::iterator jt=it;
+          ++jt;
+          if (std::fabs((*it).second) <= epsilon) comps.erase(it);
+          it=jt;
+        }
+      }
+      void simplify(Value epsilon = 0.0) const {
+        const_cast<DualExpr*>(this)->simplify(epsilon);
+      }
       ///Sets all coefficients to 0.
       void clear() {
+        Base::clear();
+      }
+      ///\e
+        comps.clear();
+      }
+      ///Compound assignment
       DualExpr &operator+=(const DualExpr &e) {
+        for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
+          (*this)[j->first]+=j->second;
+        for (std::map<int, Value>::const_iterator it=e.comps.begin();
+             it!=e.comps.end(); ++it)
+          comps[it->first]+=it->second;
         return *this;
+      }
       ///\e
+      ///Compound assignment
       DualExpr &operator-=(const DualExpr &e) {
+        for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
+          (*this)[j->first]-=j->second;
+        for (std::map<int, Value>::const_iterator it=e.comps.begin();
+             it!=e.comps.end(); ++it)
+          comps[it->first]-=it->second;
         return *this;
+      }
+      ///\e
+      DualExpr &operator*=(const Value &c) {
+        for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
+          j->second*=c;
+      ///Multiply with a constant
+      DualExpr &operator*=(const Value &v) {
+        for (std::map<int, Value>::iterator it=comps.begin();
+             it!=comps.end(); ++it)
+          it->second*=v;
         return *this;
+      }
+      ///\e
+      DualExpr &operator/=(const Value &c) {
+        for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
+          j->second/=c;
+      ///Division with a constant
+      DualExpr &operator/=(const Value &v) {
+        for (std::map<int, Value>::iterator it=comps.begin();
+             it!=comps.end(); ++it)
+          it->second/=v;
         return *this;
+      }
+      ///Iterator over the expression
+      ///The iterator iterates over the terms of the expression.
+      ///
+      ///\code
+      ///double s=0;
+      ///for(LpBase::DualExpr::CoeffIt i(e);i!=INVALID;++i)
+      ///  s+= *i * dual(i);
+      ///\endcode
+      class CoeffIt {
+      private:
+        std::map<int, Value>::iterator _it, _end;
+      public:
+        /// Sets the iterator to the first term
+        /// Sets the iterator to the first term of the expression.
+        ///
+        CoeffIt(DualExpr& e)
+          : _it(e.comps.begin()), _end(e.comps.end()){}
+        /// Convert the iterator to the row of the term
+        operator Row() const {
+          return rowFromId(_it->first);
+        }
+        /// Returns the coefficient of the term
+        Value& operator*() { return _it->second; }
+        /// Returns the coefficient of the term
+        const Value& operator*() const { return _it->second; }
+        /// Next term
+        /// Assign the iterator to the next term.
+        ///
+        CoeffIt& operator++() { ++_it; return *this; }
+        /// Equality operator
+        bool operator==(Invalid) const { return _it == _end; }
+        /// Inequality operator
+        bool operator!=(Invalid) const { return _it != _end; }
+      };
+      ///Iterator over the expression
+      ///The iterator iterates over the terms of the expression.
+      ///
+      ///\code
+      ///double s=0;
+      ///for(LpBase::DualExpr::ConstCoeffIt i(e);i!=INVALID;++i)
+      ///  s+= *i * dual(i);
+      ///\endcode
+      class ConstCoeffIt {
+      private:
+        std::map<int, Value>::const_iterator _it, _end;
+      public:
+        /// Sets the iterator to the first term
+        /// Sets the iterator to the first term of the expression.
+        ///
+        ConstCoeffIt(const DualExpr& e)
+          : _it(e.comps.begin()), _end(e.comps.end()){}
+        /// Convert the iterator to the row of the term
+        operator Row() const {
+          return rowFromId(_it->first);
+        }
+        /// Returns the coefficient of the term
+        const Value& operator*() const { return _it->second; }
+        /// Next term
+        /// Assign the iterator to the next term.
+        ///
+        ConstCoeffIt& operator++() { ++_it; return *this; }
+        /// Equality operator
+        bool operator==(Invalid) const { return _it == _end; }
+        /// Inequality operator
+        bool operator!=(Invalid) const { return _it != _end; }
+      };
     };
+  private:
+    template <typename _Expr>
+    class MappedOutputIterator {
+  protected:
+    class InsertIterator {
+    private:
+      std::map<int, Value>& _host;
+      const _solver_bits::VarIndex& _index;
     public:
-      typedef std::insert_iterator<_Expr> Base;
       typedef std::output_iterator_tag iterator_category;
 …
       typedef void pointer;
+      MappedOutputIterator(const Base& _base, const LpSolverBase& _lp)
+        : base(_base), lp(_lp) {}
+      MappedOutputIterator& operator*() {
+      InsertIterator(std::map<int, Value>& host,
+                   const _solver_bits::VarIndex& index)
+        : _host(host), _index(index) {}
+      InsertIterator& operator=(const std::pair<int, Value>& value) {
+        typedef std::map<int, Value>::value_type pair_type;
+        _host.insert(pair_type(_index[value.first], value.second));
         return *this;
+      }
+      MappedOutputIterator& operator=(const std::pair<int, Value>& value) {
+        *base = std::make_pair(lp._item(value.first, typename _Expr::Key()),
+                               value.second);
+        return *this;
+      }
+      MappedOutputIterator& operator++() {
+        ++base;
+        return *this;
+      }
+      MappedOutputIterator operator++(int) {
+        MappedOutputIterator tmp(*this);
+        ++base;
+        return tmp;
+      }
+      bool operator==(const MappedOutputIterator& it) const {
+        return base == it.base;
+      }
+      bool operator!=(const MappedOutputIterator& it) const {
+        return base != it.base;
+      }
+      InsertIterator& operator*() { return *this; }
+      InsertIterator& operator++() { return *this; }
+      InsertIterator operator++(int) { return *this; }
+    };
+    class ExprIterator {
     private:
+      Base base;
+      const LpSolverBase& lp;
+    };
+    template <typename Expr>
+    class MappedInputIterator {
+      std::map<int, Value>::const_iterator _host_it;
+      const _solver_bits::VarIndex& _index;
     public:
+      typedef typename Expr::const_iterator Base;
+      typedef typename Base::iterator_category iterator_category;
+      typedef typename Base::difference_type difference_type;
+      typedef std::bidirectional_iterator_tag iterator_category;
+      typedef std::ptrdiff_t difference_type;
       typedef const std::pair<int, Value> value_type;
       typedef value_type reference;
       class pointer {
       public:
 …
       };
+      MappedInputIterator(const Base& _base, const LpSolverBase& _lp)
+        : base(_base), lp(_lp) {}
+      ExprIterator(const std::map<int, Value>::const_iterator& host_it,
+                   const _solver_bits::VarIndex& index)
+        : _host_it(host_it), _index(index) {}
       reference operator*() {
         return std::make_pair(lp._lpId(base->first), base->second);
+        return std::make_pair(_index(_host_it->first), _host_it->second);
+      }
 …
+      }
+      MappedInputIterator& operator++() {
+        ++base;
+        return *this;
+      }
+      MappedInputIterator operator++(int) {
+        MappedInputIterator tmp(*this);
+        ++base;
+        return tmp;
+      }
+      bool operator==(const MappedInputIterator& it) const {
+        return base == it.base;
+      }
+      bool operator!=(const MappedInputIterator& it) const {
+        return base != it.base;
+      }
+    private:
+      Base base;
+      const LpSolverBase& lp;
+      ExprIterator& operator++() { ++_host_it; return *this; }
+      ExprIterator operator++(int) {
+        ExprIterator tmp(*this); ++_host_it; return tmp;
+      }
+      ExprIterator& operator--() { --_host_it; return *this; }
+      ExprIterator operator--(int) {
+        ExprIterator tmp(*this); --_host_it; return tmp;
+      }
+      bool operator==(const ExprIterator& it) const {
+        return _host_it == it._host_it;
+      }
+      bool operator!=(const ExprIterator& it) const {
+        return _host_it != it._host_it;
+      }
     };
   protected:
-    /// STL compatible iterator for lp col
-    typedef MappedInputIterator<Expr> ConstRowIterator;
-    /// STL compatible iterator for lp row
-    typedef MappedInputIterator<DualExpr> ConstColIterator;
-    /// STL compatible iterator for lp col
-    typedef MappedOutputIterator<Expr> RowIterator;
-    /// STL compatible iterator for lp row
-    typedef MappedOutputIterator<DualExpr> ColIterator;
     //Abstract virtual functions
+    virtual LpSolverBase* _newLp() = 0;
+    virtual LpSolverBase* _copyLp(){
+      LpSolverBase* newlp = _newLp();
+      std::map<Col, Col> ref;
+      for (LpSolverBase::ColIt it(*this); it != INVALID; ++it) {
+        Col ccol = newlp->addCol();
+        ref[it] = ccol;
+        newlp->colName(ccol, colName(it));
+        newlp->colLowerBound(ccol, colLowerBound(it));
+        newlp->colUpperBound(ccol, colUpperBound(it));
+      }
+      for (LpSolverBase::RowIt it(*this); it != INVALID; ++it) {
+        Expr e = row(it), ce;
+        for (Expr::iterator jt = e.begin(); jt != e.end(); ++jt) {
+          ce[ref[jt->first]] = jt->second;
+        }
+        ce += e.constComp();
+        Row r = newlp->addRow(ce);
+        double lower, upper;
+        getRowBounds(it, lower, upper);
+        newlp->rowBounds(r, lower, upper);
+      }
+      return newlp;
+    };
+    virtual LpBase* _newSolver() const = 0;
+    virtual LpBase* _cloneSolver() const = 0;
+    virtual int _addColId(int col) { return cols.addIndex(col); }
+    virtual int _addRowId(int row) { return rows.addIndex(row); }
+    virtual void _eraseColId(int col) { cols.eraseIndex(col); }
+    virtual void _eraseRowId(int row) { rows.eraseIndex(row); }
     virtual int _addCol() = 0;
 …
     virtual void _eraseRow(int row) = 0;
     virtual void _getColName(int col, std::string & name) const = 0;
     virtual void _setColName(int col, const std::string & name) = 0;
+    virtual void _getColName(int col, std::string& name) const = 0;
+    virtual void _setColName(int col, const std::string& name) = 0;
     virtual int _colByName(const std::string& name) const = 0;
+    virtual void _setRowCoeffs(int i, ConstRowIterator b,
+                               ConstRowIterator e) = 0;
+    virtual void _getRowCoeffs(int i, RowIterator b) const = 0;
+    virtual void _setColCoeffs(int i, ConstColIterator b,
+                               ConstColIterator e) = 0;
+    virtual void _getColCoeffs(int i, ColIterator b) const = 0;
+    virtual void _getRowName(int row, std::string& name) const = 0;
+    virtual void _setRowName(int row, const std::string& name) = 0;
+    virtual int _rowByName(const std::string& name) const = 0;
+    virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e) = 0;
+    virtual void _getRowCoeffs(int i, InsertIterator b) const = 0;
+    virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e) = 0;
+    virtual void _getColCoeffs(int i, InsertIterator b) const = 0;
     virtual void _setCoeff(int row, int col, Value value) = 0;
     virtual Value _getCoeff(int row, int col) const = 0;
     virtual void _setColLowerBound(int i, Value value) = 0;
     virtual Value _getColLowerBound(int i) const = 0;
     virtual void _setColUpperBound(int i, Value value) = 0;
     virtual Value _getColUpperBound(int i) const = 0;
+    virtual void _setRowBounds(int i, Value lower, Value upper) = 0;
+    virtual void _getRowBounds(int i, Value &lower, Value &upper) const = 0;
+    virtual void _setRowLowerBound(int i, Value value) = 0;
+    virtual Value _getRowLowerBound(int i) const = 0;
+    virtual void _setRowUpperBound(int i, Value value) = 0;
+    virtual Value _getRowUpperBound(int i) const = 0;
+    virtual void _setObjCoeffs(ExprIterator b, ExprIterator e) = 0;
+    virtual void _getObjCoeffs(InsertIterator b) const = 0;
     virtual void _setObjCoeff(int i, Value obj_coef) = 0;
     virtual Value _getObjCoeff(int i) const = 0;
+    virtual void _clearObj()=0;
+    virtual SolveExitStatus _solve() = 0;
+    virtual Value _getPrimal(int i) const = 0;
+    virtual Value _getDual(int i) const = 0;
+    virtual Value _getPrimalValue() const = 0;
+    virtual bool _isBasicCol(int i) const = 0;
+    virtual SolutionStatus _getPrimalStatus() const = 0;
+    virtual SolutionStatus _getDualStatus() const = 0;
+    virtual ProblemTypes _getProblemType() const = 0;
+    virtual void _setMax() = 0;
+    virtual void _setMin() = 0;
+    virtual bool _isMax() const = 0;
+    virtual void _setSense(Sense) = 0;
+    virtual Sense _getSense() const = 0;
+    virtual void _clear() = 0;
+    virtual const char* _solverName() const = 0;
     //Own protected stuff
 …
     Value obj_const_comp;
+    LpBase() : rows(), cols(), obj_const_comp(0) {}
   public:
+    ///\e
+    LpSolverBase() : obj_const_comp(0) {}
+    ///\e
+    virtual ~LpSolverBase() {}
+    /// Virtual destructor
+    virtual ~LpBase() {}
     ///Creates a new LP problem
     LpSolverBase* newLp() {return _newLp();}
+    LpBase* newSolver() {return _newSolver();}
     ///Makes a copy of the LP problem
+    LpSolverBase* copyLp() {return _copyLp();}
+    LpBase* cloneSolver() {return _cloneSolver();}
+    ///Gives back the name of the solver.
+    const char* solverName() const {return _solverName();}
     ///\name Build up and modify the LP
 …
     ///Add a new empty column (i.e a new variable) to the LP
+    Col addCol() { Col c; _addCol(); c.id = cols.addId(); return c;}
+    ///\brief Adds several new columns
+    ///(i.e a variables) at once
+    ///
+    ///This magic function takes a container as its argument
+    ///and fills its elements
+    ///with new columns (i.e. variables)
+    Col addCol() { Col c; c._id = _addColId(_addCol()); return c;}
+    ///\brief Adds several new columns (i.e variables) at once
+    ///
+    ///This magic function takes a container as its argument and fills
+    ///its elements with new columns (i.e. variables)
     ///\param t can be
     ///- a standard STL compatible iterable container with
+    ///\ref Col as its \c values_type
+    ///like
+    ///\ref Col as its \c values_type like
     ///\code
     ///std::vector<LpSolverBase::Col>
     ///std::list<LpSolverBase::Col>
+    ///std::vector<LpBase::Col>
+    ///std::list<LpBase::Col>
     ///\endcode
     ///- a standard STL compatible iterable container with
+    ///\ref Col as its \c mapped_type
+    ///like
+    ///\ref Col as its \c mapped_type like
     ///\code
     ///std::map<AnyType,LpSolverBase::Col>
+    ///std::map<AnyType,LpBase::Col>
     ///\endcode
     ///- an iterable lemon \ref concepts::WriteMap "write map" like
     ///\code
     ///ListGraph::NodeMap<LpSolverBase::Col>
     ///ListGraph::EdgeMap<LpSolverBase::Col>
+    ///ListGraph::NodeMap<LpBase::Col>
+    ///ListGraph::ArcMap<LpBase::Col>
     ///\endcode
     ///\return The number of the created column.
 …
 #else
     template<class T>
     typename enable_if<typename T::value_type::LpSolverCol,int>::type
+    typename enable_if<typename T::value_type::LpCol,int>::type
     addColSet(T &t,dummy<0> = 0) {
       int s=0;
 …
+    }
     template<class T>
     typename enable_if<typename T::value_type::second_type::LpSolverCol,
+    typename enable_if<typename T::value_type::second_type::LpCol,
                        int>::type
     addColSet(T &t,dummy<1> = 1) {
 …
+    }
     template<class T>
     typename enable_if<typename T::MapIt::Value::LpSolverCol,
+    typename enable_if<typename T::MapIt::Value::LpCol,
                        int>::type
     addColSet(T &t,dummy<2> = 2) {
 …
     ///\param e is a dual linear expression (see \ref DualExpr)
     ///a better one.
     void col(Col c,const DualExpr &e) {
+    void col(Col c, const DualExpr &e) {
       e.simplify();
       _setColCoeffs(_lpId(c), ConstColIterator(e.begin(), *this),
                     ConstColIterator(e.end(), *this));
+      _setColCoeffs(cols(id(c)), ExprIterator(e.comps.begin(), cols),
+                    ExprIterator(e.comps.end(), cols));
+    }
     ///Get a column (i.e a dual constraint) of the LP
     ///\param r is the column to get
+    ///\param c is the column to get
     ///\return the dual expression associated to the column
     DualExpr col(Col c) const {
       DualExpr e;
       _getColCoeffs(_lpId(c), ColIterator(std::inserter(e, e.end()), *this));
+      _getColCoeffs(cols(id(c)), InsertIterator(e.comps, rows));
       return e;
+    }
 …
     ///\param e is a dual linear expression (see \ref DualExpr)
     ///\param obj is the corresponding component of the objective
+    ///\param o is the corresponding component of the objective
     ///function. It is 0 by default.
     ///\return The created column.
 …
     ///This function adds a new empty row (i.e a new constraint) to the LP.
     ///\return The created row
+    Row addRow() { Row r; _addRow(); r.id = rows.addId(); return r;}
+    ///\brief Add several new rows
+    ///(i.e a constraints) at once
+    ///
+    ///This magic function takes a container as its argument
+    ///and fills its elements
+    ///with new row (i.e. variables)
+    Row addRow() { Row r; r._id = _addRowId(_addRow()); return r;}
+    ///\brief Add several new rows (i.e constraints) at once
+    ///
+    ///This magic function takes a container as its argument and fills
+    ///its elements with new row (i.e. variables)
     ///\param t can be
     ///- a standard STL compatible iterable container with
+    ///\ref Row as its \c values_type
+    ///like
+    ///\ref Row as its \c values_type like
     ///\code
     ///std::vector<LpSolverBase::Row>
     ///std::list<LpSolverBase::Row>
+    ///std::vector<LpBase::Row>
+    ///std::list<LpBase::Row>
     ///\endcode
     ///- a standard STL compatible iterable container with
+    ///\ref Row as its \c mapped_type
+    ///like
+    ///\ref Row as its \c mapped_type like
     ///\code
     ///std::map<AnyType,LpSolverBase::Row>
+    ///std::map<AnyType,LpBase::Row>
     ///\endcode
     ///- an iterable lemon \ref concepts::WriteMap "write map" like
     ///\code
     ///ListGraph::NodeMap<LpSolverBase::Row>
     ///ListGraph::EdgeMap<LpSolverBase::Row>
+    ///ListGraph::NodeMap<LpBase::Row>
+    ///ListGraph::ArcMap<LpBase::Row>
     ///\endcode
     ///\return The number of rows created.
 …
 #else
     template<class T>
     typename enable_if<typename T::value_type::LpSolverRow,int>::type
     addRowSet(T &t,dummy<0> = 0) {
+    typename enable_if<typename T::value_type::LpRow,int>::type
+    addRowSet(T &t, dummy<0> = 0) {
       int s=0;
       for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addRow();s++;}
 …
+    }
     template<class T>
+    typename enable_if<typename T::value_type::second_type::LpSolverRow,
+                       int>::type
+    addRowSet(T &t,dummy<1> = 1) {
+    typename enable_if<typename T::value_type::second_type::LpRow, int>::type
+    addRowSet(T &t, dummy<1> = 1) {
       int s=0;
       for(typename T::iterator i=t.begin();i!=t.end();++i) {
 …
+    }
     template<class T>
+    typename enable_if<typename T::MapIt::Value::LpSolverRow,
+                       int>::type
+    addRowSet(T &t,dummy<2> = 2) {
+    typename enable_if<typename T::MapIt::Value::LpRow, int>::type
+    addRowSet(T &t, dummy<2> = 2) {
       int s=0;
       for(typename T::MapIt i(t); i!=INVALID; ++i)
 …
     ///\param e is a linear expression (see \ref Expr)
     ///\param u is the upper bound (\ref INF means no bound)
-    ///\bug This is a temporary function. The interface will change to
-    ///a better one.
-    ///\todo Option to control whether a constraint with a single variable is
-    ///added or not.
     void row(Row r, Value l, const Expr &e, Value u) {
       e.simplify();
+      _setRowCoeffs(_lpId(r), ConstRowIterator(e.begin(), *this),
+                    ConstRowIterator(e.end(), *this));
+      _setRowBounds(_lpId(r),l-e.constComp(),u-e.constComp());
+      _setRowCoeffs(rows(id(r)), ExprIterator(e.comps.begin(), cols),
+                    ExprIterator(e.comps.end(), cols));
+      _setRowLowerBound(rows(id(r)),l - *e);
+      _setRowUpperBound(rows(id(r)),u - *e);
+    }
 …
     Expr row(Row r) const {
       Expr e;
       _getRowCoeffs(_lpId(r), RowIterator(std::inserter(e, e.end()), *this));
+      _getRowCoeffs(rows(id(r)), InsertIterator(e.comps, cols));
       return e;
+    }
 …
     ///\param u is the upper bound (\ref INF means no bound)
     ///\return The created row.
-    ///\bug This is a temporary function. The interface will change to
-    ///a better one.
     Row addRow(Value l,const Expr &e, Value u) {
       Row r=addRow();
 …
       return r;
+    }
+    ///Erase a coloumn (i.e a variable) from the LP
+    ///\param c is the coloumn to be deleted
+    ///\todo Please check this
+    void eraseCol(Col c) {
+      _eraseCol(_lpId(c));
+      cols.eraseId(c.id);
+    }
+    ///Erase a  row (i.e a constraint) from the LP
+    ///Erase a column (i.e a variable) from the LP
+    ///\param c is the column to be deleted
+    void erase(Col c) {
+      _eraseCol(cols(id(c)));
+      _eraseColId(cols(id(c)));
+    }
+    ///Erase a row (i.e a constraint) from the LP
     ///\param r is the row to be deleted
+    ///\todo Please check this
+    void eraseRow(Row r) {
+      _eraseRow(_lpId(r));
+      rows.eraseId(r.id);
+    void erase(Row r) {
+      _eraseRow(rows(id(r)));
+      _eraseRowId(rows(id(r)));
+    }
     /// Get the name of a column
     ///\param c is the coresponding coloumn
+    ///\param c is the coresponding column
     ///\return The name of the colunm
     std::string colName(Col c) const {
       std::string name;
       _getColName(_lpId(c), name);
+      _getColName(cols(id(c)), name);
       return name;
+    }
 …
     /// Set the name of a column
     ///\param c is the coresponding coloumn
+    ///\param c is the coresponding column
     ///\param name The name to be given
     void colName(Col c, const std::string& name) {
       _setColName(_lpId(c), name);
+      _setColName(cols(id(c)), name);
+    }
 …
     Col colByName(const std::string& name) const {
       int k = _colByName(name);
+      return k != -1 ? Col(cols.fixId(k)) : Col(INVALID);
+      return k != -1 ? Col(cols[k]) : Col(INVALID);
+    }
+    /// Get the name of a row
+    ///\param r is the coresponding row
+    ///\return The name of the row
+    std::string rowName(Row r) const {
+      std::string name;
+      _getRowName(rows(id(r)), name);
+      return name;
+    }
+    /// Set the name of a row
+    ///\param r is the coresponding row
+    ///\param name The name to be given
+    void rowName(Row r, const std::string& name) {
+      _setRowName(rows(id(r)), name);
+    }
+    /// Get the row by its name
+    ///\param name The name of the row
+    ///\return the proper row or \c INVALID
+    Row rowByName(const std::string& name) const {
+      int k = _rowByName(name);
+      return k != -1 ? Row(rows[k]) : Row(INVALID);
+    }
 …
     ///\param r is the row of the element to be modified
     ///\param c is the coloumn of the element to be modified
+    ///\param c is the column of the element to be modified
     ///\param val is the new value of the coefficient
     void coeff(Row r, Col c, Value val) {
       _setCoeff(_lpId(r),_lpId(c), val);
+      _setCoeff(rows(id(r)),cols(id(c)), val);
+    }
     /// Get an element of the coefficient matrix of the LP
     ///\param r is the row of the element in question
     ///\param c is the coloumn of the element in question
+    ///\param r is the row of the element
+    ///\param c is the column of the element
     ///\return the corresponding coefficient
     Value coeff(Row r, Col c) const {
       return _getCoeff(_lpId(r),_lpId(c));
+      return _getCoeff(rows(id(r)),cols(id(c)));
+    }
 …
     /// Value or -\ref INF.
     void colLowerBound(Col c, Value value) {
       _setColLowerBound(_lpId(c),value);
+      _setColLowerBound(cols(id(c)),value);
+    }
     /// Get the lower bound of a column (i.e a variable)
     /// This function returns the lower bound for column (variable) \t c
+    /// This function returns the lower bound for column (variable) \c c
     /// (this might be -\ref INF as well).
     ///\return The lower bound for coloumn \t c
+    ///\return The lower bound for column \c c
     Value colLowerBound(Col c) const {
       return _getColLowerBound(_lpId(c));
+      return _getColLowerBound(cols(id(c)));
+    }
     ///\brief Set the lower bound of  several columns
     ///(i.e a variables) at once
+    ///(i.e variables) at once
     ///
     ///This magic function takes a container as its argument
     ///and applies the function on all of its elements.
     /// The lower bound of a variable (column) has to be given by an
     /// extended number of type Value, i.e. a finite number of type
     /// Value or -\ref INF.
+    ///The lower bound of a variable (column) has to be given by an
+    ///extended number of type Value, i.e. a finite number of type
+    ///Value or -\ref INF.
 #ifdef DOXYGEN
     template<class T>
 …
 #else
     template<class T>
     typename enable_if<typename T::value_type::LpSolverCol,void>::type
+    typename enable_if<typename T::value_type::LpCol,void>::type
     colLowerBound(T &t, Value value,dummy<0> = 0) {
       for(typename T::iterator i=t.begin();i!=t.end();++i) {
 …
+    }
     template<class T>
     typename enable_if<typename T::value_type::second_type::LpSolverCol,
+    typename enable_if<typename T::value_type::second_type::LpCol,
                        void>::type
     colLowerBound(T &t, Value value,dummy<1> = 1) {
 …
+    }
     template<class T>
     typename enable_if<typename T::MapIt::Value::LpSolverCol,
+    typename enable_if<typename T::MapIt::Value::LpCol,
                        void>::type
     colLowerBound(T &t, Value value,dummy<2> = 2) {
 …
     /// Value or \ref INF.
     void colUpperBound(Col c, Value value) {
       _setColUpperBound(_lpId(c),value);
+      _setColUpperBound(cols(id(c)),value);
     };
     /// Get the upper bound of a column (i.e a variable)
     /// This function returns the upper bound for column (variable) \t c
+    /// This function returns the upper bound for column (variable) \c c
     /// (this might be \ref INF as well).
     ///\return The upper bound for coloumn \t c
+    /// \return The upper bound for column \c c
     Value colUpperBound(Col c) const {
       return _getColUpperBound(_lpId(c));
+      return _getColUpperBound(cols(id(c)));
+    }
     ///\brief Set the upper bound of  several columns
     ///(i.e a variables) at once
+    ///(i.e variables) at once
     ///
     ///This magic function takes a container as its argument
     ///and applies the function on all of its elements.
     /// The upper bound of a variable (column) has to be given by an
     /// extended number of type Value, i.e. a finite number of type
     /// Value or \ref INF.
+    ///The upper bound of a variable (column) has to be given by an
+    ///extended number of type Value, i.e. a finite number of type
+    ///Value or \ref INF.
 #ifdef DOXYGEN
     template<class T>
 …
 #else
     template<class T>
     typename enable_if<typename T::value_type::LpSolverCol,void>::type
+    typename enable_if<typename T::value_type::LpCol,void>::type
     colUpperBound(T &t, Value value,dummy<0> = 0) {
       for(typename T::iterator i=t.begin();i!=t.end();++i) {
 …
+    }
     template<class T>
     typename enable_if<typename T::value_type::second_type::LpSolverCol,
+    typename enable_if<typename T::value_type::second_type::LpCol,
                        void>::type
     colUpperBound(T &t, Value value,dummy<1> = 1) {
 …
+    }
     template<class T>
     typename enable_if<typename T::MapIt::Value::LpSolverCol,
+    typename enable_if<typename T::MapIt::Value::LpCol,
                        void>::type
     colUpperBound(T &t, Value value,dummy<2> = 2) {
 …
     /// Value, -\ref INF or \ref INF.
     void colBounds(Col c, Value lower, Value upper) {
       _setColLowerBound(_lpId(c),lower);
       _setColUpperBound(_lpId(c),upper);
+      _setColLowerBound(cols(id(c)),lower);
+      _setColUpperBound(cols(id(c)),upper);
+    }
     ///\brief Set the lower and the upper bound of several columns
     ///(i.e a variables) at once
+    ///(i.e variables) at once
     ///
     ///This magic function takes a container as its argument
 …
 #else
     template<class T>
     typename enable_if<typename T::value_type::LpSolverCol,void>::type
+    typename enable_if<typename T::value_type::LpCol,void>::type
     colBounds(T &t, Value lower, Value upper,dummy<0> = 0) {
       for(typename T::iterator i=t.begin();i!=t.end();++i) {
 …
+    }
     template<class T>
+    typename enable_if<typename T::value_type::second_type::LpSolverCol,
+                       void>::type
+    typename enable_if<typename T::value_type::second_type::LpCol, void>::type
     colBounds(T &t, Value lower, Value upper,dummy<1> = 1) {
       for(typename T::iterator i=t.begin();i!=t.end();++i) {
 …
+    }
     template<class T>
+    typename enable_if<typename T::MapIt::Value::LpSolverCol,
+                       void>::type
+    typename enable_if<typename T::MapIt::Value::LpCol, void>::type
     colBounds(T &t, Value lower, Value upper,dummy<2> = 2) {
       for(typename T::MapIt i(t); i!=INVALID; ++i){
 …
 #endif
+    /// Set the lower and the upper bounds of a row (i.e a constraint)
+    /// The lower and the upper bound of a constraint (row) have to be
+    /// given by an extended number of type Value, i.e. a finite
+    /// number of type Value, -\ref INF or \ref INF. There is no
+    /// separate function for the lower and the upper bound because
+    /// that would have been hard to implement for CPLEX.
+    void rowBounds(Row c, Value lower, Value upper) {
+      _setRowBounds(_lpId(c),lower, upper);
+    }
+    /// Get the lower and the upper bounds of a row (i.e a constraint)
+    /// The lower and the upper bound of
+    /// a constraint (row) are
+    /// extended numbers of type Value, i.e.  finite numbers of type
+    /// Value, -\ref INF or \ref INF.
+    /// \todo There is no separate function for the
+    /// lower and the upper bound because we had problems with the
+    /// implementation of the setting functions for CPLEX:
+    /// check out whether this can be done for these functions.
+    void getRowBounds(Row c, Value &lower, Value &upper) const {
+      _getRowBounds(_lpId(c),lower, upper);
+    /// Set the lower bound of a row (i.e a constraint)
+    /// The lower bound of a constraint (row) has to be given by an
+    /// extended number of type Value, i.e. a finite number of type
+    /// Value or -\ref INF.
+    void rowLowerBound(Row r, Value value) {
+      _setRowLowerBound(rows(id(r)),value);
+    }
+    /// Get the lower bound of a row (i.e a constraint)
+    /// This function returns the lower bound for row (constraint) \c c
+    /// (this might be -\ref INF as well).
+    ///\return The lower bound for row \c r
+    Value rowLowerBound(Row r) const {
+      return _getRowLowerBound(rows(id(r)));
+    }
+    /// Set the upper bound of a row (i.e a constraint)
+    /// The upper bound of a constraint (row) has to be given by an
+    /// extended number of type Value, i.e. a finite number of type
+    /// Value or -\ref INF.
+    void rowUpperBound(Row r, Value value) {
+      _setRowUpperBound(rows(id(r)),value);
+    }
+    /// Get the upper bound of a row (i.e a constraint)
+    /// This function returns the upper bound for row (constraint) \c c
+    /// (this might be -\ref INF as well).
+    ///\return The upper bound for row \c r
+    Value rowUpperBound(Row r) const {
+      return _getRowUpperBound(rows(id(r)));
+    }
     ///Set an element of the objective function
     void objCoeff(Col c, Value v) {_setObjCoeff(_lpId(c),v); };
+    void objCoeff(Col c, Value v) {_setObjCoeff(cols(id(c)),v); };
     ///Get an element of the objective function
     Value objCoeff(Col c) const { return _getObjCoeff(_lpId(c)); };
+    Value objCoeff(Col c) const { return _getObjCoeff(cols(id(c))); };
     ///Set the objective function
     ///\param e is a linear expression of type \ref Expr.
     void obj(Expr e) {
       _clearObj();
       for (Expr::iterator i=e.begin(); i!=e.end(); ++i)
         objCoeff((*i).first,(*i).second);
       obj_const_comp=e.constComp();
+    ///
+    void obj(const Expr& e) {
+      _setObjCoeffs(ExprIterator(e.comps.begin(), cols),
+                    ExprIterator(e.comps.end(), cols));
+      obj_const_comp = *e;
+    }
     ///Get the objective function
+    ///\return the objective function as a linear expression of type \ref Expr.
+    ///\return the objective function as a linear expression of type
+    ///Expr.
     Expr obj() const {
       Expr e;
+      for (ColIt it(*this); it != INVALID; ++it) {
+        double c = objCoeff(it);
+        if (c != 0.0) {
+          e.insert(std::make_pair(it, c));
+        }
+      }
+      _getObjCoeffs(InsertIterator(e.comps, cols));
+      *e = obj_const_comp;
       return e;
+    }
+    ///Maximize
+    void max() { _setMax(); }
+    ///Minimize
+    void min() { _setMin(); }
+    ///Query function: is this a maximization problem?
+    bool isMax() const {return _isMax(); }
+    ///Query function: is this a minimization problem?
+    bool isMin() const {return !isMax(); }
+    ///Set the direction of optimization
+    void sense(Sense sense) { _setSense(sense); }
+    ///Query the direction of the optimization
+    Sense sense() const {return _getSense(); }
+    ///Set the sense to maximization
+    void max() { _setSense(MAX); }
+    ///Set the sense to maximization
+    void min() { _setSense(MIN); }
+    ///Clears the problem
+    void clear() { _clear(); }
     ///@}
+  };
+  /// Addition
+  ///\relates LpBase::Expr
+  ///
+  inline LpBase::Expr operator+(const LpBase::Expr &a, const LpBase::Expr &b) {
+    LpBase::Expr tmp(a);
+    tmp+=b;
+    return tmp;
+  }
+  ///Substraction
+  ///\relates LpBase::Expr
+  ///
+  inline LpBase::Expr operator-(const LpBase::Expr &a, const LpBase::Expr &b) {
+    LpBase::Expr tmp(a);
+    tmp-=b;
+    return tmp;
+  }
+  ///Multiply with constant
+  ///\relates LpBase::Expr
+  ///
+  inline LpBase::Expr operator*(const LpBase::Expr &a, const LpBase::Value &b) {
+    LpBase::Expr tmp(a);
+    tmp*=b;
+    return tmp;
+  }
+  ///Multiply with constant
+  ///\relates LpBase::Expr
+  ///
+  inline LpBase::Expr operator*(const LpBase::Value &a, const LpBase::Expr &b) {
+    LpBase::Expr tmp(b);
+    tmp*=a;
+    return tmp;
+  }
+  ///Divide with constant
+  ///\relates LpBase::Expr
+  ///
+  inline LpBase::Expr operator/(const LpBase::Expr &a, const LpBase::Value &b) {
+    LpBase::Expr tmp(a);
+    tmp/=b;
+    return tmp;
+  }
+  ///Create constraint
+  ///\relates LpBase::Constr
+  ///
+  inline LpBase::Constr operator<=(const LpBase::Expr &e,
+                                   const LpBase::Expr &f) {
+    return LpBase::Constr(0, f - e, LpBase::INF);
+  }
+  ///Create constraint
+  ///\relates LpBase::Constr
+  ///
+  inline LpBase::Constr operator<=(const LpBase::Value &e,
+                                   const LpBase::Expr &f) {
+    return LpBase::Constr(e, f, LpBase::NaN);
+  }
+  ///Create constraint
+  ///\relates LpBase::Constr
+  ///
+  inline LpBase::Constr operator<=(const LpBase::Expr &e,
+                                   const LpBase::Value &f) {
+    return LpBase::Constr(- LpBase::INF, e, f);
+  }
+  ///Create constraint
+  ///\relates LpBase::Constr
+  ///
+  inline LpBase::Constr operator>=(const LpBase::Expr &e,
+                                   const LpBase::Expr &f) {
+    return LpBase::Constr(0, e - f, LpBase::INF);
+  }
+  ///Create constraint
+  ///\relates LpBase::Constr
+  ///
+  inline LpBase::Constr operator>=(const LpBase::Value &e,
+                                   const LpBase::Expr &f) {
+    return LpBase::Constr(LpBase::NaN, f, e);
+  }
+  ///Create constraint
+  ///\relates LpBase::Constr
+  ///
+  inline LpBase::Constr operator>=(const LpBase::Expr &e,
+                                   const LpBase::Value &f) {
+    return LpBase::Constr(f, e, LpBase::INF);
+  }
+  ///Create constraint
+  ///\relates LpBase::Constr
+  ///
+  inline LpBase::Constr operator==(const LpBase::Expr &e,
+                                   const LpBase::Value &f) {
+    return LpBase::Constr(f, e, f);
+  }
+  ///Create constraint
+  ///\relates LpBase::Constr
+  ///
+  inline LpBase::Constr operator==(const LpBase::Expr &e,
+                                   const LpBase::Expr &f) {
+    return LpBase::Constr(0, f - e, 0);
+  }
+  ///Create constraint
+  ///\relates LpBase::Constr
+  ///
+  inline LpBase::Constr operator<=(const LpBase::Value &n,
+                                   const LpBase::Constr &c) {
+    LpBase::Constr tmp(c);
+    LEMON_ASSERT(std::isnan(tmp.lowerBound()), "Wrong LP constraint");
+    tmp.lowerBound()=n;
+    return tmp;
+  }
+  ///Create constraint
+  ///\relates LpBase::Constr
+  ///
+  inline LpBase::Constr operator<=(const LpBase::Constr &c,
+                                   const LpBase::Value &n)
+  {
+    LpBase::Constr tmp(c);
+    LEMON_ASSERT(std::isnan(tmp.upperBound()), "Wrong LP constraint");
+    tmp.upperBound()=n;
+    return tmp;
+  }
+  ///Create constraint
+  ///\relates LpBase::Constr
+  ///
+  inline LpBase::Constr operator>=(const LpBase::Value &n,
+                                   const LpBase::Constr &c) {
+    LpBase::Constr tmp(c);
+    LEMON_ASSERT(std::isnan(tmp.upperBound()), "Wrong LP constraint");
+    tmp.upperBound()=n;
+    return tmp;
+  }
+  ///Create constraint
+  ///\relates LpBase::Constr
+  ///
+  inline LpBase::Constr operator>=(const LpBase::Constr &c,
+                                   const LpBase::Value &n)
+  {
+    LpBase::Constr tmp(c);
+    LEMON_ASSERT(std::isnan(tmp.lowerBound()), "Wrong LP constraint");
+    tmp.lowerBound()=n;
+    return tmp;
+  }
+  ///Addition
+  ///\relates LpBase::DualExpr
+  ///
+  inline LpBase::DualExpr operator+(const LpBase::DualExpr &a,
+                                    const LpBase::DualExpr &b) {
+    LpBase::DualExpr tmp(a);
+    tmp+=b;
+    return tmp;
+  }
+  ///Substraction
+  ///\relates LpBase::DualExpr
+  ///
+  inline LpBase::DualExpr operator-(const LpBase::DualExpr &a,
+                                    const LpBase::DualExpr &b) {
+    LpBase::DualExpr tmp(a);
+    tmp-=b;
+    return tmp;
+  }
+  ///Multiply with constant
+  ///\relates LpBase::DualExpr
+  ///
+  inline LpBase::DualExpr operator*(const LpBase::DualExpr &a,
+                                    const LpBase::Value &b) {
+    LpBase::DualExpr tmp(a);
+    tmp*=b;
+    return tmp;
+  }
+  ///Multiply with constant
+  ///\relates LpBase::DualExpr
+  ///
+  inline LpBase::DualExpr operator*(const LpBase::Value &a,
+                                    const LpBase::DualExpr &b) {
+    LpBase::DualExpr tmp(b);
+    tmp*=a;
+    return tmp;
+  }
+  ///Divide with constant
+  ///\relates LpBase::DualExpr
+  ///
+  inline LpBase::DualExpr operator/(const LpBase::DualExpr &a,
+                                    const LpBase::Value &b) {
+    LpBase::DualExpr tmp(a);
+    tmp/=b;
+    return tmp;
+  }
+  /// \ingroup lp_group
+  ///
+  /// \brief Common base class for LP solvers
+  ///
+  /// This class is an abstract base class for LP solvers. This class
+  /// provides a full interface for set and modify an LP problem,
+  /// solve it and retrieve the solution. You can use one of the
+  /// descendants as a concrete implementation, or the \c Lp
+  /// default LP solver. However, if you would like to handle LP
+  /// solvers as reference or pointer in a generic way, you can use
+  /// this class directly.
+  class LpSolver : virtual public LpBase {
+  public:
+    /// The problem types for primal and dual problems
+    enum ProblemType {
+      ///Feasible solution hasn't been found (but may exist).
+      UNDEFINED = 0,
+      ///The problem has no feasible solution
+      INFEASIBLE = 1,
+      ///Feasible solution found
+      FEASIBLE = 2,
+      ///Optimal solution exists and found
+      OPTIMAL = 3,
+      ///The cost function is unbounded
+      UNBOUNDED = 4
+    };
+    ///The basis status of variables
+    enum VarStatus {
+      /// The variable is in the basis
+      BASIC,
+      /// The variable is free, but not basic
+      FREE,
+      /// The variable has active lower bound
+      LOWER,
+      /// The variable has active upper bound
+      UPPER,
+      /// The variable is non-basic and fixed
+      FIXED
+    };
+  protected:
+    virtual SolveExitStatus _solve() = 0;
+    virtual Value _getPrimal(int i) const = 0;
+    virtual Value _getDual(int i) const = 0;
+    virtual Value _getPrimalRay(int i) const = 0;
+    virtual Value _getDualRay(int i) const = 0;
+    virtual Value _getPrimalValue() const = 0;
+    virtual VarStatus _getColStatus(int i) const = 0;
+    virtual VarStatus _getRowStatus(int i) const = 0;
+    virtual ProblemType _getPrimalType() const = 0;
+    virtual ProblemType _getDualType() const = 0;
+  public:
     ///\name Solve the LP
 …
     ///values and their meanings can be found in the documentation of
     ///\ref SolveExitStatus.
-    ///
-    ///\todo Which method is used to solve the problem
     SolveExitStatus solve() { return _solve(); }
 …
     ///@{
+    /// The status of the primal problem (the original LP problem)
+    SolutionStatus primalStatus() const {
+      return _getPrimalStatus();
+    }
+    /// The status of the dual (of the original LP) problem
+    SolutionStatus dualStatus() const {
+      return _getDualStatus();
+    }
+    ///The type of the original LP problem
+    ProblemTypes problemType() const {
+      return _getProblemType();
+    }
+    ///\e
+    Value primal(Col c) const { return _getPrimal(_lpId(c)); }
+    ///\e
+    /// The type of the primal problem
+    ProblemType primalType() const {
+      return _getPrimalType();
+    }
+    /// The type of the dual problem
+    ProblemType dualType() const {
+      return _getDualType();
+    }
+    /// Return the primal value of the column
+    /// Return the primal value of the column.
+    /// \pre The problem is solved.
+    Value primal(Col c) const { return _getPrimal(cols(id(c))); }
+    /// Return the primal value of the expression
+    /// Return the primal value of the expression, i.e. the dot
+    /// product of the primal solution and the expression.
+    /// \pre The problem is solved.
     Value primal(const Expr& e) const {
+      double res = e.constComp();
+      for (std::map<Col, double>::const_iterator it = e.begin();
+           it != e.end(); ++it) {
+        res += _getPrimal(_lpId(it->first)) * it->second;
+      double res = *e;
+      for (Expr::ConstCoeffIt c(e); c != INVALID; ++c) {
+        res += *c * primal(c);
+      }
       return res;
+    }
+    ///\e
+    Value dual(Row r) const { return _getDual(_lpId(r)); }
+    ///\e
+    /// Returns a component of the primal ray
+    /// The primal ray is solution of the modified primal problem,
+    /// where we change each finite bound to 0, and we looking for a
+    /// negative objective value in case of minimization, and positive
+    /// objective value for maximization. If there is such solution,
+    /// that proofs the unsolvability of the dual problem, and if a
+    /// feasible primal solution exists, then the unboundness of
+    /// primal problem.
+    ///
+    /// \pre The problem is solved and the dual problem is infeasible.
+    /// \note Some solvers does not provide primal ray calculation
+    /// functions.
+    Value primalRay(Col c) const { return _getPrimalRay(cols(id(c))); }
+    /// Return the dual value of the row
+    /// Return the dual value of the row.
+    /// \pre The problem is solved.
+    Value dual(Row r) const { return _getDual(rows(id(r))); }
+    /// Return the dual value of the dual expression
+    /// Return the dual value of the dual expression, i.e. the dot
+    /// product of the dual solution and the dual expression.
+    /// \pre The problem is solved.
     Value dual(const DualExpr& e) const {
       double res = 0.0;
+      for (std::map<Row, double>::const_iterator it = e.begin();
+           it != e.end(); ++it) {
+        res += _getPrimal(_lpId(it->first)) * it->second;
+      for (DualExpr::ConstCoeffIt r(e); r != INVALID; ++r) {
+        res += *r * dual(r);
+      }
       return res;
+    }
+    ///\e
+    bool isBasicCol(Col c) const { return _isBasicCol(_lpId(c)); }
+    ///\e
+    /// Returns a component of the dual ray
+    /// The dual ray is solution of the modified primal problem, where
+    /// we change each finite bound to 0 (i.e. the objective function
+    /// coefficients in the primal problem), and we looking for a
+    /// ositive objective value. If there is such solution, that
+    /// proofs the unsolvability of the primal problem, and if a
+    /// feasible dual solution exists, then the unboundness of
+    /// dual problem.
+    ///
+    /// \pre The problem is solved and the primal problem is infeasible.
+    /// \note Some solvers does not provide dual ray calculation
+    /// functions.
+    Value dualRay(Row r) const { return _getDualRay(rows(id(r))); }
+    /// Return the basis status of the column
+    /// \see VarStatus
+    VarStatus colStatus(Col c) const { return _getColStatus(cols(id(c))); }
+    /// Return the basis status of the row
+    /// \see VarStatus
+    VarStatus rowStatus(Row r) const { return _getRowStatus(rows(id(r))); }
+    ///The value of the objective function
     ///\return
 …
     ///- \ref NaN if no primal solution is found.
     ///- The (finite) objective value if an optimal solution is found.
     Value primalValue() const { return _getPrimalValue()+obj_const_comp;}
+    Value primal() const { return _getPrimalValue()+obj_const_comp;}
     ///@}
+    LpSolver* newSolver() {return _newSolver();}
+    LpSolver* cloneSolver() {return _cloneSolver();}
+  protected:
+    virtual LpSolver* _newSolver() const = 0;
+    virtual LpSolver* _cloneSolver() const = 0;
   };
 …
   ///
   /// \brief Common base class for MIP solvers
+  /// \todo Much more docs
+  class MipSolverBase : virtual public LpSolverBase{
+  ///
+  /// This class is an abstract base class for MIP solvers. This class
+  /// provides a full interface for set and modify an MIP problem,
+  /// solve it and retrieve the solution. You can use one of the
+  /// descendants as a concrete implementation, or the \c Lp
+  /// default MIP solver. However, if you would like to handle MIP
+  /// solvers as reference or pointer in a generic way, you can use
+  /// this class directly.
+  class MipSolver : virtual public LpBase {
   public:
+    ///Possible variable (coloumn) types (e.g. real, integer, binary etc.)
+    /// The problem types for MIP problems
+    enum ProblemType {
+      ///Feasible solution hasn't been found (but may exist).
+      UNDEFINED = 0,
+      ///The problem has no feasible solution
+      INFEASIBLE = 1,
+      ///Feasible solution found
+      FEASIBLE = 2,
+      ///Optimal solution exists and found
+      OPTIMAL = 3,
+      ///The cost function is unbounded
+      ///
+      ///The Mip or at least the relaxed problem is unbounded
+      UNBOUNDED = 4
+    };
+    ///\name Solve the MIP
+    ///@{
+    /// Solve the MIP problem at hand
+    ///
+    ///\return The result of the optimization procedure. Possible
+    ///values and their meanings can be found in the documentation of
+    ///\ref SolveExitStatus.
+    SolveExitStatus solve() { return _solve(); }
+    ///@}
+    ///\name Setting column type
+    ///@{
+    ///Possible variable (column) types (e.g. real, integer, binary etc.)
     enum ColTypes {
       ///Continuous variable
+      ///Continuous variable (default)
       REAL = 0,
       ///Integer variable
+      ///Unfortunately, cplex 7.5 somewhere writes something like
+      ///#define INTEGER 'I'
+      INT = 1
+      ///\todo No support for other types yet.
+      INTEGER = 1
     };
+    ///Sets the type of the given coloumn to the given type
+    ///
+    ///Sets the type of the given coloumn to the given type.
+    ///Sets the type of the given column to the given type
+    ///Sets the type of the given column to the given type.
+    ///
     void colType(Col c, ColTypes col_type) {
       _colType(_lpId(c),col_type);
+      _setColType(cols(id(c)),col_type);
+    }
     ///Gives back the type of the column.
+    ///
     ///Gives back the type of the column.
+    ///
     ColTypes colType(Col c) const {
+      return _colType(_lpId(c));
+    }
+    ///Sets the type of the given Col to integer or remove that property.
+    ///
+    ///Sets the type of the given Col to integer or remove that property.
+    void integer(Col c, bool enable) {
+      if (enable)
+        colType(c,INT);
+      else
+        colType(c,REAL);
+    }
+    ///Gives back whether the type of the column is integer or not.
+    ///
+    ///Gives back the type of the column.
+    ///\return true if the column has integer type and false if not.
+    bool integer(Col c) const {
+      return (colType(c)==INT);
+    }
+    /// The status of the MIP problem
+    SolutionStatus mipStatus() const {
+      return _getMipStatus();
+    }
+      return _getColType(cols(id(c)));
+    }
+    ///@}
+    ///\name Obtain the solution
+    ///@{
+    /// The type of the MIP problem
+    ProblemType type() const {
+      return _getType();
+    }
+    /// Return the value of the row in the solution
+    ///  Return the value of the row in the solution.
+    /// \pre The problem is solved.
+    Value sol(Col c) const { return _getSol(cols(id(c))); }
+    /// Return the value of the expression in the solution
+    /// Return the value of the expression in the solution, i.e. the
+    /// dot product of the solution and the expression.
+    /// \pre The problem is solved.
+    Value sol(const Expr& e) const {
+      double res = *e;
+      for (Expr::ConstCoeffIt c(e); c != INVALID; ++c) {
+        res += *c * sol(c);
+      }
+      return res;
+    }
+    ///The value of the objective function
+    ///\return
+    ///- \ref INF or -\ref INF means either infeasibility or unboundedness
+    /// of the problem, depending on whether we minimize or maximize.
+    ///- \ref NaN if no primal solution is found.
+    ///- The (finite) objective value if an optimal solution is found.
+    Value solValue() const { return _getSolValue()+obj_const_comp;}
+    ///@}
   protected:
+    virtual ColTypes _colType(int col) const = 0;
+    virtual void _colType(int col, ColTypes col_type) = 0;
+    virtual SolutionStatus _getMipStatus() const = 0;
+    virtual SolveExitStatus _solve() = 0;
+    virtual ColTypes _getColType(int col) const = 0;
+    virtual void _setColType(int col, ColTypes col_type) = 0;
+    virtual ProblemType _getType() const = 0;
+    virtual Value _getSol(int i) const = 0;
+    virtual Value _getSolValue() const = 0;
+  public:
+    MipSolver* newSolver() {return _newSolver();}
+    MipSolver* cloneSolver() {return _cloneSolver();}
+  protected:
+    virtual MipSolver* _newSolver() const = 0;
+    virtual MipSolver* _cloneSolver() const = 0;
   };
-  ///\relates LpSolverBase::Expr
-  ///
-  inline LpSolverBase::Expr operator+(const LpSolverBase::Expr &a,
-                                      const LpSolverBase::Expr &b)
+  {
-    LpSolverBase::Expr tmp(a);
-    tmp+=b;
-    return tmp;
+  }
-  ///\e
-  ///\relates LpSolverBase::Expr
-  ///
-  inline LpSolverBase::Expr operator-(const LpSolverBase::Expr &a,
-                                      const LpSolverBase::Expr &b)
+  {
-    LpSolverBase::Expr tmp(a);
-    tmp-=b;
-    return tmp;
+  }
-  ///\e
-  ///\relates LpSolverBase::Expr
-  ///
-  inline LpSolverBase::Expr operator*(const LpSolverBase::Expr &a,
-                                      const LpSolverBase::Value &b)
+  {
-    LpSolverBase::Expr tmp(a);
-    tmp*=b;
-    return tmp;
+  }
-  ///\e
-  ///\relates LpSolverBase::Expr
-  ///
-  inline LpSolverBase::Expr operator*(const LpSolverBase::Value &a,
-                                      const LpSolverBase::Expr &b)
+  {
-    LpSolverBase::Expr tmp(b);
-    tmp*=a;
-    return tmp;
+  }
-  ///\e
-  ///\relates LpSolverBase::Expr
-  ///
-  inline LpSolverBase::Expr operator/(const LpSolverBase::Expr &a,
-                                      const LpSolverBase::Value &b)
+  {
-    LpSolverBase::Expr tmp(a);
-    tmp/=b;
-    return tmp;
+  }
-  ///\e
-  ///\relates LpSolverBase::Constr
-  ///
-  inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,
-                                         const LpSolverBase::Expr &f)
+  {
-    return LpSolverBase::Constr(-LpSolverBase::INF,e-f,0);
+  }
-  ///\e
-  ///\relates LpSolverBase::Constr
-  ///
-  inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &e,
-                                         const LpSolverBase::Expr &f)
+  {
-    return LpSolverBase::Constr(e,f);
+  }
-  ///\e
-  ///\relates LpSolverBase::Constr
-  ///
-  inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,
-                                         const LpSolverBase::Value &f)
+  {
-    return LpSolverBase::Constr(-LpSolverBase::INF,e,f);
+  }
-  ///\e
-  ///\relates LpSolverBase::Constr
-  ///
-  inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,
-                                         const LpSolverBase::Expr &f)
+  {
-    return LpSolverBase::Constr(-LpSolverBase::INF,f-e,0);
+  }
-  ///\e
-  ///\relates LpSolverBase::Constr
-  ///
-  inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &e,
-                                         const LpSolverBase::Expr &f)
+  {
-    return LpSolverBase::Constr(f,e);
+  }
-  ///\e
-  ///\relates LpSolverBase::Constr
-  ///
-  inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,
-                                         const LpSolverBase::Value &f)
+  {
-    return LpSolverBase::Constr(f,e,LpSolverBase::INF);
+  }
-  ///\e
-  ///\relates LpSolverBase::Constr
-  ///
-  inline LpSolverBase::Constr operator==(const LpSolverBase::Expr &e,
-                                         const LpSolverBase::Value &f)
+  {
-    return LpSolverBase::Constr(f,e,f);
+  }
-  ///\e
-  ///\relates LpSolverBase::Constr
-  ///
-  inline LpSolverBase::Constr operator==(const LpSolverBase::Expr &e,
-                                         const LpSolverBase::Expr &f)
+  {
-    return LpSolverBase::Constr(0,e-f,0);
+  }
-  ///\e
-  ///\relates LpSolverBase::Constr
-  ///
-  inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &n,
-                                         const LpSolverBase::Constr&c)
+  {
-    LpSolverBase::Constr tmp(c);
-    LEMON_ASSERT(LpSolverBase::isNaN(tmp.lowerBound()), "Wrong LP constraint");
-    tmp.lowerBound()=n;
-    return tmp;
+  }
-  ///\e
-  ///\relates LpSolverBase::Constr
-  ///
-  inline LpSolverBase::Constr operator<=(const LpSolverBase::Constr& c,
-                                         const LpSolverBase::Value &n)
+  {
-    LpSolverBase::Constr tmp(c);
-    LEMON_ASSERT(LpSolverBase::isNaN(tmp.upperBound()), "Wrong LP constraint");
-    tmp.upperBound()=n;
-    return tmp;
+  }
-  ///\e
-  ///\relates LpSolverBase::Constr
-  ///
-  inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &n,
-                                         const LpSolverBase::Constr&c)
+  {
-    LpSolverBase::Constr tmp(c);
-    LEMON_ASSERT(LpSolverBase::isNaN(tmp.upperBound()), "Wrong LP constraint");
-    tmp.upperBound()=n;
-    return tmp;
+  }
-  ///\e
-  ///\relates LpSolverBase::Constr
-  ///
-  inline LpSolverBase::Constr operator>=(const LpSolverBase::Constr& c,
-                                         const LpSolverBase::Value &n)
+  {
-    LpSolverBase::Constr tmp(c);
-    LEMON_ASSERT(LpSolverBase::isNaN(tmp.lowerBound()), "Wrong LP constraint");
-    tmp.lowerBound()=n;
-    return tmp;
+  }
-  ///\e
-  ///\relates LpSolverBase::DualExpr
-  ///
-  inline LpSolverBase::DualExpr operator+(const LpSolverBase::DualExpr &a,
-                                          const LpSolverBase::DualExpr &b)
+  {
-    LpSolverBase::DualExpr tmp(a);
-    tmp+=b;
-    return tmp;
+  }
-  ///\e
-  ///\relates LpSolverBase::DualExpr
-  ///
-  inline LpSolverBase::DualExpr operator-(const LpSolverBase::DualExpr &a,
-                                          const LpSolverBase::DualExpr &b)
+  {
-    LpSolverBase::DualExpr tmp(a);
-    tmp-=b;
-    return tmp;
+  }
-  ///\e
-  ///\relates LpSolverBase::DualExpr
-  ///
-  inline LpSolverBase::DualExpr operator*(const LpSolverBase::DualExpr &a,
-                                          const LpSolverBase::Value &b)
+  {
-    LpSolverBase::DualExpr tmp(a);
-    tmp*=b;
-    return tmp;
+  }
-  ///\e
-  ///\relates LpSolverBase::DualExpr
-  ///
-  inline LpSolverBase::DualExpr operator*(const LpSolverBase::Value &a,
-                                          const LpSolverBase::DualExpr &b)
+  {
-    LpSolverBase::DualExpr tmp(b);
-    tmp*=a;
-    return tmp;
+  }
-  ///\e
-  ///\relates LpSolverBase::DualExpr
-  ///
-  inline LpSolverBase::DualExpr operator/(const LpSolverBase::DualExpr &a,
-                                          const LpSolverBase::Value &b)
+  {
-    LpSolverBase::DualExpr tmp(a);
-    tmp/=b;
-    return tmp;
+  }

lemon/lp_cplex.cc

-                      r458
+                      r459
 #include <iostream>
 #include <vector>
+#include <cstring>
 #include <lemon/lp_cplex.h>
 …
 namespace lemon {
+  LpCplex::LpCplex() {
+    //    env = CPXopenCPLEXdevelop(&status);
+    env = CPXopenCPLEX(&status);
+    lp = CPXcreateprob(env, &status, "LP problem");
+  }
+  LpCplex::LpCplex(const LpCplex& cplex) : LpSolverBase() {
+    env = CPXopenCPLEX(&status);
+    lp = CPXcloneprob(env, cplex.lp, &status);
+  CplexEnv::LicenseError::LicenseError(int status) {
+    if (!CPXgeterrorstring(0, status, _message)) {
+      std::strcpy(_message, "Cplex unknown error");
+    }
+  }
+  CplexEnv::CplexEnv() {
+    int status;
+    _cnt = new int;
+    _env = CPXopenCPLEX(&status);
+    if (_env == 0) {
+      delete _cnt;
+      _cnt = 0;
+      throw LicenseError(status);
+    }
+  }
+  CplexEnv::CplexEnv(const CplexEnv& other) {
+    _env = other._env;
+    _cnt = other._cnt;
+    ++(*_cnt);
+  }
+  CplexEnv& CplexEnv::operator=(const CplexEnv& other) {
+    _env = other._env;
+    _cnt = other._cnt;
+    ++(*_cnt);
+    return *this;
+  }
+  CplexEnv::~CplexEnv() {
+    --(*_cnt);
+    if (*_cnt == 0) {
+      delete _cnt;
+      CPXcloseCPLEX(&_env);
+    }
+  }
+  CplexBase::CplexBase() : LpBase() {
+    int status;
+    _prob = CPXcreateprob(cplexEnv(), &status, "Cplex problem");
+  }
+  CplexBase::CplexBase(const CplexEnv& env)
+    : LpBase(), _env(env) {
+    int status;
+    _prob = CPXcreateprob(cplexEnv(), &status, "Cplex problem");
+  }
+  CplexBase::CplexBase(const CplexBase& cplex)
+    : LpBase() {
+    int status;
+    _prob = CPXcloneprob(cplexEnv(), cplex._prob, &status);
     rows = cplex.rows;
     cols = cplex.cols;
+  }
+  LpCplex::~LpCplex() {
+    CPXfreeprob(env,&lp);
+    CPXcloseCPLEX(&env);
+  }
+  LpSolverBase* LpCplex::_newLp()
+  {
+    //The first approach opens a new environment
+    return new LpCplex();
+  }
+  LpSolverBase* LpCplex::_copyLp() {
+    return new LpCplex(*this);
+  }
+  int LpCplex::_addCol()
+  {
+    int i = CPXgetnumcols(env, lp);
+    Value lb[1],ub[1];
+    lb[0]=-INF;
+    ub[0]=INF;
+    status = CPXnewcols(env, lp, 1, NULL, lb, ub, NULL, NULL);
+  CplexBase::~CplexBase() {
+    CPXfreeprob(cplexEnv(),&_prob);
+  }
+  int CplexBase::_addCol() {
+    int i = CPXgetnumcols(cplexEnv(), _prob);
+    double lb = -INF, ub = INF;
+    CPXnewcols(cplexEnv(), _prob, 1, 0, &lb, &ub, 0, 0);
     return i;
+  }
+  int LpCplex::_addRow()
+  {
+    //We want a row that is not constrained
+    char sense[1];
+    sense[0]='L';//<= constraint
+    Value rhs[1];
+    rhs[0]=INF;
+    int i = CPXgetnumrows(env, lp);
+    status = CPXnewrows(env, lp, 1, rhs, sense, NULL, NULL);
+  int CplexBase::_addRow() {
+    int i = CPXgetnumrows(cplexEnv(), _prob);
+    const double ub = INF;
+    const char s = 'L';
+    CPXnewrows(cplexEnv(), _prob, 1, &ub, &s, 0, 0);
     return i;
+  }
+  void LpCplex::_eraseCol(int i) {
+    CPXdelcols(env, lp, i, i);
+  }
+  void LpCplex::_eraseRow(int i) {
+    CPXdelrows(env, lp, i, i);
+  }
+  void LpCplex::_getColName(int col, std::string &name) const
+  {
+    ///\bug Untested
+    int storespace;
+    CPXgetcolname(env, lp, 0, 0, 0, &storespace, col, col);
+    if (storespace == 0) {
+  void CplexBase::_eraseCol(int i) {
+    CPXdelcols(cplexEnv(), _prob, i, i);
+  }
+  void CplexBase::_eraseRow(int i) {
+    CPXdelrows(cplexEnv(), _prob, i, i);
+  }
+  void CplexBase::_eraseColId(int i) {
+    cols.eraseIndex(i);
+    cols.shiftIndices(i);
+  }
+  void CplexBase::_eraseRowId(int i) {
+    rows.eraseIndex(i);
+    rows.shiftIndices(i);
+  }
+  void CplexBase::_getColName(int col, std::string &name) const {
+    int size;
+    CPXgetcolname(cplexEnv(), _prob, 0, 0, 0, &size, col, col);
+    if (size == 0) {
       name.clear();
       return;
+    }
+    storespace *= -1;
+    std::vector<char> buf(storespace);
+    char *names[1];
+    int dontcare;
+    ///\bug return code unchecked for error
+    CPXgetcolname(env, lp, names, &*buf.begin(), storespace,
+                  &dontcare, col, col);
+    name = names[0];
+  }
+  void LpCplex::_setColName(int col, const std::string &name)
+  {
+    ///\bug Untested
+    char *names[1];
+    names[0] = const_cast<char*>(name.c_str());
+    ///\bug return code unchecked for error
+    CPXchgcolname(env, lp, 1, &col, names);
+  }
+  int LpCplex::_colByName(const std::string& name) const
+  {
+    size *= -1;
+    std::vector<char> buf(size);
+    char *cname;
+    int tmp;
+    CPXgetcolname(cplexEnv(), _prob, &cname, &buf.front(), size,
+                  &tmp, col, col);
+    name = cname;
+  }
+  void CplexBase::_setColName(int col, const std::string &name) {
+    char *cname;
+    cname = const_cast<char*>(name.c_str());
+    CPXchgcolname(cplexEnv(), _prob, 1, &col, &cname);
+  }
+  int CplexBase::_colByName(const std::string& name) const {
     int index;
     if (CPXgetcolindex(env, lp,
+    if (CPXgetcolindex(cplexEnv(), _prob,
                        const_cast<char*>(name.c_str()), &index) == 0) {
       return index;
 …
+  }
+  ///\warning Data at index 0 is ignored in the arrays.
+  void LpCplex::_setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e)
+  void CplexBase::_getRowName(int row, std::string &name) const {
+    int size;
+    CPXgetrowname(cplexEnv(), _prob, 0, 0, 0, &size, row, row);
+    if (size == 0) {
+      name.clear();
+      return;
+    }
+    size *= -1;
+    std::vector<char> buf(size);
+    char *cname;
+    int tmp;
+    CPXgetrowname(cplexEnv(), _prob, &cname, &buf.front(), size,
+                  &tmp, row, row);
+    name = cname;
+  }
+  void CplexBase::_setRowName(int row, const std::string &name) {
+    char *cname;
+    cname = const_cast<char*>(name.c_str());
+    CPXchgrowname(cplexEnv(), _prob, 1, &row, &cname);
+  }
+  int CplexBase::_rowByName(const std::string& name) const {
+    int index;
+    if (CPXgetrowindex(cplexEnv(), _prob,
+                       const_cast<char*>(name.c_str()), &index) == 0) {
+      return index;
+    }
+    return -1;
+  }
+  void CplexBase::_setRowCoeffs(int i, ExprIterator b,
+                                      ExprIterator e)
+  {
     std::vector<int> indices;
 …
     std::vector<Value> values;
     for(ConstRowIterator it=b; it!=e; ++it) {
+    for(ExprIterator it=b; it!=e; ++it) {
       indices.push_back(it->first);
       values.push_back(it->second);
 …
+    }
     status = CPXchgcoeflist(env, lp, values.size(),
                             &rowlist[0], &indices[0], &values[0]);
+  }
   void LpCplex::_getRowCoeffs(int i, RowIterator b) const {
+    CPXchgcoeflist(cplexEnv(), _prob, values.size(),
+                   &rowlist.front(), &indices.front(), &values.front());
+  }
+  void CplexBase::_getRowCoeffs(int i, InsertIterator b) const {
     int tmp1, tmp2, tmp3, length;
     CPXgetrows(env, lp, &tmp1, &tmp2, 0, 0, 0, &length, i, i);
+    CPXgetrows(cplexEnv(), _prob, &tmp1, &tmp2, 0, 0, 0, &length, i, i);
     length = -length;
 …
     std::vector<double> values(length);
+    CPXgetrows(env, lp, &tmp1, &tmp2, &indices[0], &values[0],
+    CPXgetrows(cplexEnv(), _prob, &tmp1, &tmp2,
+               &indices.front(), &values.front(),
                length, &tmp3, i, i);
 …
       ++b;
+    }
+    /// \todo implement
+  }
+  void LpCplex::_setColCoeffs(int i, ConstColIterator b, ConstColIterator e)
+  {
+  }
+  void CplexBase::_setColCoeffs(int i, ExprIterator b, ExprIterator e) {
     std::vector<int> indices;
     std::vector<int> collist;
     std::vector<Value> values;
     for(ConstColIterator it=b; it!=e; ++it) {
+    for(ExprIterator it=b; it!=e; ++it) {
       indices.push_back(it->first);
       values.push_back(it->second);
 …
+    }
     status = CPXchgcoeflist(env, lp, values.size(),
                             &indices[0], &collist[0], &values[0]);
+  }
   void LpCplex::_getColCoeffs(int i, ColIterator b) const {
+    CPXchgcoeflist(cplexEnv(), _prob, values.size(),
+                   &indices.front(), &collist.front(), &values.front());
+  }
+  void CplexBase::_getColCoeffs(int i, InsertIterator b) const {
     int tmp1, tmp2, tmp3, length;
     CPXgetcols(env, lp, &tmp1, &tmp2, 0, 0, 0, &length, i, i);
+    CPXgetcols(cplexEnv(), _prob, &tmp1, &tmp2, 0, 0, 0, &length, i, i);
     length = -length;
 …
     std::vector<double> values(length);
+    CPXgetcols(env, lp, &tmp1, &tmp2, &indices[0], &values[0],
+    CPXgetcols(cplexEnv(), _prob, &tmp1, &tmp2,
+               &indices.front(), &values.front(),
                length, &tmp3, i, i);
 …
+  }
+  void LpCplex::_setCoeff(int row, int col, Value value)
+  void CplexBase::_setCoeff(int row, int col, Value value) {
+    CPXchgcoef(cplexEnv(), _prob, row, col, value);
+  }
+  CplexBase::Value CplexBase::_getCoeff(int row, int col) const {
+    CplexBase::Value value;
+    CPXgetcoef(cplexEnv(), _prob, row, col, &value);
+    return value;
+  }
+  void CplexBase::_setColLowerBound(int i, Value value) {
+    const char s = 'L';
+    CPXchgbds(cplexEnv(), _prob, 1, &i, &s, &value);
+  }
+  CplexBase::Value CplexBase::_getColLowerBound(int i) const {
+    CplexBase::Value res;
+    CPXgetlb(cplexEnv(), _prob, &res, i, i);
+    return res <= -CPX_INFBOUND ? -INF : res;
+  }
+  void CplexBase::_setColUpperBound(int i, Value value)
+  {
+    CPXchgcoef(env, lp, row, col, value);
+  }
+  LpCplex::Value LpCplex::_getCoeff(int row, int col) const
+    const char s = 'U';
+    CPXchgbds(cplexEnv(), _prob, 1, &i, &s, &value);
+  }
+  CplexBase::Value CplexBase::_getColUpperBound(int i) const {
+    CplexBase::Value res;
+    CPXgetub(cplexEnv(), _prob, &res, i, i);
+    return res >= CPX_INFBOUND ? INF : res;
+  }
+  CplexBase::Value CplexBase::_getRowLowerBound(int i) const {
+    char s;
+    CPXgetsense(cplexEnv(), _prob, &s, i, i);
+    CplexBase::Value res;
+    switch (s) {
+    case 'G':
+    case 'R':
+    case 'E':
+      CPXgetrhs(cplexEnv(), _prob, &res, i, i);
+      return res <= -CPX_INFBOUND ? -INF : res;
+    default:
+      return -INF;
+    }
+  }
+  CplexBase::Value CplexBase::_getRowUpperBound(int i) const {
+    char s;
+    CPXgetsense(cplexEnv(), _prob, &s, i, i);
+    CplexBase::Value res;
+    switch (s) {
+    case 'L':
+    case 'E':
+      CPXgetrhs(cplexEnv(), _prob, &res, i, i);
+      return res >= CPX_INFBOUND ? INF : res;
+    case 'R':
+      CPXgetrhs(cplexEnv(), _prob, &res, i, i);
+      {
+        double rng;
+        CPXgetrngval(cplexEnv(), _prob, &rng, i, i);
+        res += rng;
+      }
+      return res >= CPX_INFBOUND ? INF : res;
+    default:
+      return INF;
+    }
+  }
+  //This is easier to implement
+  void CplexBase::_set_row_bounds(int i, Value lb, Value ub) {
+    if (lb == -INF) {
+      const char s = 'L';
+      CPXchgsense(cplexEnv(), _prob, 1, &i, &s);
+      CPXchgrhs(cplexEnv(), _prob, 1, &i, &ub);
+    } else if (ub == INF) {
+      const char s = 'G';
+      CPXchgsense(cplexEnv(), _prob, 1, &i, &s);
+      CPXchgrhs(cplexEnv(), _prob, 1, &i, &lb);
+    } else if (lb == ub){
+      const char s = 'E';
+      CPXchgsense(cplexEnv(), _prob, 1, &i, &s);
+      CPXchgrhs(cplexEnv(), _prob, 1, &i, &lb);
+    } else {
+      const char s = 'R';
+      CPXchgsense(cplexEnv(), _prob, 1, &i, &s);
+      CPXchgrhs(cplexEnv(), _prob, 1, &i, &lb);
+      double len = ub - lb;
+      CPXchgrngval(cplexEnv(), _prob, 1, &i, &len);
+    }
+  }
+  void CplexBase::_setRowLowerBound(int i, Value lb)
+  {
+    LpCplex::Value value;
+    CPXgetcoef(env, lp, row, col, &value);
+    return value;
+  }
+  void LpCplex::_setColLowerBound(int i, Value value)
+    LEMON_ASSERT(lb != INF, "Invalid bound");
+    _set_row_bounds(i, lb, CplexBase::_getRowUpperBound(i));
+  }
+  void CplexBase::_setRowUpperBound(int i, Value ub)
+  {
+    int indices[1];
+    indices[0]=i;
+    char lu[1];
+    lu[0]='L';
+    Value bd[1];
+    bd[0]=value;
+    status = CPXchgbds(env, lp, 1, indices, lu, bd);
+  }
+  LpCplex::Value LpCplex::_getColLowerBound(int i) const
+    LEMON_ASSERT(ub != -INF, "Invalid bound");
+    _set_row_bounds(i, CplexBase::_getRowLowerBound(i), ub);
+  }
+  void CplexBase::_setObjCoeffs(ExprIterator b, ExprIterator e)
+  {
+    LpCplex::Value x;
+    CPXgetlb (env, lp, &x, i, i);
+    if (x <= -CPX_INFBOUND) x = -INF;
+    return x;
+  }
+  void LpCplex::_setColUpperBound(int i, Value value)
+    std::vector<int> indices;
+    std::vector<Value> values;
+    for(ExprIterator it=b; it!=e; ++it) {
+      indices.push_back(it->first);
+      values.push_back(it->second);
+    }
+    CPXchgobj(cplexEnv(), _prob, values.size(),
+              &indices.front(), &values.front());
+  }
+  void CplexBase::_getObjCoeffs(InsertIterator b) const
+  {
+    int indices[1];
+    indices[0]=i;
+    char lu[1];
+    lu[0]='U';
+    Value bd[1];
+    bd[0]=value;
+    status = CPXchgbds(env, lp, 1, indices, lu, bd);
+  }
+  LpCplex::Value LpCplex::_getColUpperBound(int i) const
+    int num = CPXgetnumcols(cplexEnv(), _prob);
+    std::vector<Value> x(num);
+    CPXgetobj(cplexEnv(), _prob, &x.front(), 0, num - 1);
+    for (int i = 0; i < num; ++i) {
+      if (x[i] != 0.0) {
+        *b = std::make_pair(i, x[i]);
+        ++b;
+      }
+    }
+  }
+  void CplexBase::_setObjCoeff(int i, Value obj_coef)
+  {
+    LpCplex::Value x;
+    CPXgetub (env, lp, &x, i, i);
+    if (x >= CPX_INFBOUND) x = INF;
+    return x;
+  }
+  //This will be easier to implement
+  void LpCplex::_setRowBounds(int i, Value lb, Value ub)
+  {
+    //Bad parameter
+    if (lb==INF || ub==-INF) {
+      //FIXME error
+    }
+    int cnt=1;
+    int indices[1];
+    indices[0]=i;
+    char sense[1];
+    if (lb==-INF){
+      sense[0]='L';
+      CPXchgsense(env, lp, cnt, indices, sense);
+      CPXchgcoef(env, lp, i, -1, ub);
+    }
+    else{
+      if (ub==INF){
+        sense[0]='G';
+        CPXchgsense(env, lp, cnt, indices, sense);
+        CPXchgcoef(env, lp, i, -1, lb);
+      }
+      else{
+        if (lb == ub){
+          sense[0]='E';
+          CPXchgsense(env, lp, cnt, indices, sense);
+          CPXchgcoef(env, lp, i, -1, lb);
+        }
+        else{
+          sense[0]='R';
+          CPXchgsense(env, lp, cnt, indices, sense);
+          CPXchgcoef(env, lp, i, -1, lb);
+          CPXchgcoef(env, lp, i, -2, ub-lb);
+        }
+      }
+    }
+  }
+//   void LpCplex::_setRowLowerBound(int i, Value value)
+//   {
+//     //Not implemented, obsolete
+//   }
+//   void LpCplex::_setRowUpperBound(int i, Value value)
+//   {
+//     //Not implemented, obsolete
+// //     //TODO Ezt kell meg megirni
+// //     //type of the problem
+// //     char sense[1];
+// //     status = CPXgetsense(env, lp, sense, i, i);
+// //     Value rhs[1];
+// //     status = CPXgetrhs(env, lp, rhs, i, i);
+// //     switch (sense[0]) {
+// //     case 'L'://<= constraint
+// //       break;
+// //     case 'E'://= constraint
+// //       break;
+// //     case 'G'://>= constraint
+// //       break;
+// //     case 'R'://ranged constraint
+// //       break;
+// //     default: ;
+// //       //FIXME error
+// //     }
+// //     status = CPXchgcoef(env, lp, i, -2, value_rng);
+//   }
+  void LpCplex::_getRowBounds(int i, Value &lb, Value &ub) const
+  {
+    char sense;
+    CPXgetsense(env, lp, &sense,i,i);
+    lb=-INF;
+    ub=INF;
+    switch (sense)
+      {
+      case 'L':
+        CPXgetcoef(env, lp, i, -1, &ub);
+        break;
+      case 'G':
+        CPXgetcoef(env, lp, i, -1, &lb);
+        break;
+      case 'E':
+        CPXgetcoef(env, lp, i, -1, &lb);
+        ub=lb;
+        break;
+      case 'R':
+        CPXgetcoef(env, lp, i, -1, &lb);
+        Value x;
+        CPXgetcoef(env, lp, i, -2, &x);
+        ub=lb+x;
+        break;
+      }
+  }
+  void LpCplex::_setObjCoeff(int i, Value obj_coef)
+  {
+    CPXchgcoef(env, lp, -1, i, obj_coef);
+  }
+  LpCplex::Value LpCplex::_getObjCoeff(int i) const
+    CPXchgobj(cplexEnv(), _prob, 1, &i, &obj_coef);
+  }
+  CplexBase::Value CplexBase::_getObjCoeff(int i) const
+  {
     Value x;
     CPXgetcoef(env, lp, -1, i, &x);
+    CPXgetobj(cplexEnv(), _prob, &x, i, i);
     return x;
+  }
+  void LpCplex::_clearObj()
+  {
+    for (int i=0;i< CPXgetnumcols(env, lp);++i){
+      CPXchgcoef(env, lp, -1, i, 0);
+    }
+  }
+  void CplexBase::_setSense(CplexBase::Sense sense) {
+    switch (sense) {
+    case MIN:
+      CPXchgobjsen(cplexEnv(), _prob, CPX_MIN);
+      break;
+    case MAX:
+      CPXchgobjsen(cplexEnv(), _prob, CPX_MAX);
+      break;
+    }
+  }
+  CplexBase::Sense CplexBase::_getSense() const {
+    switch (CPXgetobjsen(cplexEnv(), _prob)) {
+    case CPX_MIN:
+      return MIN;
+    case CPX_MAX:
+      return MAX;
+    default:
+      LEMON_ASSERT(false, "Invalid sense");
+      return CplexBase::Sense();
+    }
+  }
+  void CplexBase::_clear() {
+    CPXfreeprob(cplexEnv(),&_prob);
+    int status;
+    _prob = CPXcreateprob(cplexEnv(), &status, "Cplex problem");
+    rows.clear();
+    cols.clear();
+  }
+  // LpCplex members
+  LpCplex::LpCplex()
+    : LpBase(), CplexBase(), LpSolver() {}
+  LpCplex::LpCplex(const CplexEnv& env)
+    : LpBase(), CplexBase(env), LpSolver() {}
+  LpCplex::LpCplex(const LpCplex& other)
+    : LpBase(), CplexBase(other), LpSolver() {}
+  LpCplex::~LpCplex() {}
+  LpCplex* LpCplex::_newSolver() const { return new LpCplex; }
+  LpCplex* LpCplex::_cloneSolver() const {return new LpCplex(*this); }
+  const char* LpCplex::_solverName() const { return "LpCplex"; }
+  void LpCplex::_clear_temporals() {
+    _col_status.clear();
+    _row_status.clear();
+    _primal_ray.clear();
+    _dual_ray.clear();
+  }
   // The routine returns zero unless an error occurred during the
   // optimization. Examples of errors include exhausting available
 …
   // routines CPXsolninfo, CPXgetstat, and CPXsolution to obtain
   // further information about the status of the optimization.
+  LpCplex::SolveExitStatus LpCplex::_solve()
+  {
+    //CPX_PARAM_LPMETHOD
+    status = CPXlpopt(env, lp);
+    //status = CPXprimopt(env, lp);
+  LpCplex::SolveExitStatus LpCplex::convertStatus(int status) {
 #if CPX_VERSION >= 800
+    if (status)
+    {
+    if (status == 0) {
+      switch (CPXgetstat(cplexEnv(), _prob)) {
+      case CPX_STAT_OPTIMAL:
+      case CPX_STAT_INFEASIBLE:
+      case CPX_STAT_UNBOUNDED:
+        return SOLVED;
+      default:
+        return UNSOLVED;
+      }
+    } else {
       return UNSOLVED;
+    }
-    else
+    {
-      switch (CPXgetstat(env, lp))
+      {
-        case CPX_STAT_OPTIMAL:
-        case CPX_STAT_INFEASIBLE:
-        case CPX_STAT_UNBOUNDED:
-          return SOLVED;
-        default:
-          return UNSOLVED;
+      }
+    }
 #else
     if (status == 0){
+    if (status == 0) {
       //We want to exclude some cases
       switch (CPXgetstat(env, lp)){
+      switch (CPXgetstat(cplexEnv(), _prob)) {
       case CPX_OBJ_LIM:
       case CPX_IT_LIM_FEAS:
 …
         return SOLVED;
+      }
+    }
+    else{
+    } else {
       return UNSOLVED;
+    }
 …
+  }
+  LpCplex::Value LpCplex::_getPrimal(int i) const
+  {
+  LpCplex::SolveExitStatus LpCplex::_solve() {
+    _clear_temporals();
+    return convertStatus(CPXlpopt(cplexEnv(), _prob));
+  }
+  LpCplex::SolveExitStatus LpCplex::solvePrimal() {
+    _clear_temporals();
+    return convertStatus(CPXprimopt(cplexEnv(), _prob));
+  }
+  LpCplex::SolveExitStatus LpCplex::solveDual() {
+    _clear_temporals();
+    return convertStatus(CPXdualopt(cplexEnv(), _prob));
+  }
+  LpCplex::SolveExitStatus LpCplex::solveBarrier() {
+    _clear_temporals();
+    return convertStatus(CPXbaropt(cplexEnv(), _prob));
+  }
+  LpCplex::Value LpCplex::_getPrimal(int i) const {
     Value x;
     CPXgetx(env, lp, &x, i, i);
+    CPXgetx(cplexEnv(), _prob, &x, i, i);
     return x;
+  }
+  LpCplex::Value LpCplex::_getDual(int i) const
+  {
+  LpCplex::Value LpCplex::_getDual(int i) const {
     Value y;
     CPXgetpi(env, lp, &y, i, i);
+    CPXgetpi(cplexEnv(), _prob, &y, i, i);
     return y;
+  }
+  LpCplex::Value LpCplex::_getPrimalValue() const
+  {
+  LpCplex::Value LpCplex::_getPrimalValue() const {
     Value objval;
+    //method = CPXgetmethod (env, lp);
+    //printf("CPXgetprobtype %d \n",CPXgetprobtype(env,lp));
+    CPXgetobjval(env, lp, &objval);
+    //printf("Objective value: %g \n",objval);
+    CPXgetobjval(cplexEnv(), _prob, &objval);
     return objval;
+  }
+  bool LpCplex::_isBasicCol(int i) const
+  {
+    std::vector<int> cstat(CPXgetnumcols(env, lp));
+    CPXgetbase(env, lp, &*cstat.begin(), NULL);
+    return (cstat[i]==CPX_BASIC);
+  }
+//7.5-os cplex statusai (Vigyazat: a 9.0-asei masok!)
+// This table lists the statuses, returned by the CPXgetstat()
+// routine, for solutions to LP problems or mixed integer problems. If
+// no solution exists, the return value is zero.
+// For Simplex, Barrier
+// 1          CPX_OPTIMAL
+//          Optimal solution found
+// 2          CPX_INFEASIBLE
+//          Problem infeasible
+// 3    CPX_UNBOUNDED
+//          Problem unbounded
+// 4          CPX_OBJ_LIM
+//          Objective limit exceeded in Phase II
+// 5          CPX_IT_LIM_FEAS
+//          Iteration limit exceeded in Phase II
+// 6          CPX_IT_LIM_INFEAS
+//          Iteration limit exceeded in Phase I
+// 7          CPX_TIME_LIM_FEAS
+//          Time limit exceeded in Phase II
+// 8          CPX_TIME_LIM_INFEAS
+//          Time limit exceeded in Phase I
+// 9          CPX_NUM_BEST_FEAS
+//          Problem non-optimal, singularities in Phase II
+// 10         CPX_NUM_BEST_INFEAS
+//          Problem non-optimal, singularities in Phase I
+// 11         CPX_OPTIMAL_INFEAS
+//          Optimal solution found, unscaled infeasibilities
+// 12         CPX_ABORT_FEAS
+//          Aborted in Phase II
+// 13         CPX_ABORT_INFEAS
+//          Aborted in Phase I
+// 14          CPX_ABORT_DUAL_INFEAS
+//          Aborted in barrier, dual infeasible
+// 15          CPX_ABORT_PRIM_INFEAS
+//          Aborted in barrier, primal infeasible
+// 16          CPX_ABORT_PRIM_DUAL_INFEAS
+//          Aborted in barrier, primal and dual infeasible
+// 17          CPX_ABORT_PRIM_DUAL_FEAS
+//          Aborted in barrier, primal and dual feasible
+// 18          CPX_ABORT_CROSSOVER
+//          Aborted in crossover
+// 19          CPX_INForUNBD
+//          Infeasible or unbounded
+// 20   CPX_PIVOT
+//       User pivot used
+//
+//     Ezeket hova tegyem:
+// ??case CPX_ABORT_DUAL_INFEAS
+// ??case CPX_ABORT_CROSSOVER
+// ??case CPX_INForUNBD
+// ??case CPX_PIVOT
+//Some more interesting stuff:
+// CPX_PARAM_LPMETHOD  1062  int  LPMETHOD
+// 0 Automatic
+// 1 Primal Simplex
+// 2 Dual Simplex
+// 3 Network Simplex
+// 4 Standard Barrier
+// Default: 0
+// Description: Method for linear optimization.
+// Determines which algorithm is used when CPXlpopt() (or "optimize"
+// in the Interactive Optimizer) is called. Currently the behavior of
+// the "Automatic" setting is that CPLEX simply invokes the dual
+// simplex method, but this capability may be expanded in the future
+// so that CPLEX chooses the method based on problem characteristics
+  LpCplex::VarStatus LpCplex::_getColStatus(int i) const {
+    if (_col_status.empty()) {
+      _col_status.resize(CPXgetnumcols(cplexEnv(), _prob));
+      CPXgetbase(cplexEnv(), _prob, &_col_status.front(), 0);
+    }
+    switch (_col_status[i]) {
+    case CPX_BASIC:
+      return BASIC;
+    case CPX_FREE_SUPER:
+      return FREE;
+    case CPX_AT_LOWER:
+      return LOWER;
+    case CPX_AT_UPPER:
+      return UPPER;
+    default:
+      LEMON_ASSERT(false, "Wrong column status");
+      return LpCplex::VarStatus();
+    }
+  }
+  LpCplex::VarStatus LpCplex::_getRowStatus(int i) const {
+    if (_row_status.empty()) {
+      _row_status.resize(CPXgetnumrows(cplexEnv(), _prob));
+      CPXgetbase(cplexEnv(), _prob, 0, &_row_status.front());
+    }
+    switch (_row_status[i]) {
+    case CPX_BASIC:
+      return BASIC;
+    case CPX_AT_LOWER:
+      {
+        char s;
+        CPXgetsense(cplexEnv(), _prob, &s, i, i);
+        return s != 'L' ? LOWER : UPPER;
+      }
+    case CPX_AT_UPPER:
+      return UPPER;
+    default:
+      LEMON_ASSERT(false, "Wrong row status");
+      return LpCplex::VarStatus();
+    }
+  }
+  LpCplex::Value LpCplex::_getPrimalRay(int i) const {
+    if (_primal_ray.empty()) {
+      _primal_ray.resize(CPXgetnumcols(cplexEnv(), _prob));
+      CPXgetray(cplexEnv(), _prob, &_primal_ray.front());
+    }
+    return _primal_ray[i];
+  }
+  LpCplex::Value LpCplex::_getDualRay(int i) const {
+    if (_dual_ray.empty()) {
+    }
+    return _dual_ray[i];
+  }
+  //7.5-os cplex statusai (Vigyazat: a 9.0-asei masok!)
+  // This table lists the statuses, returned by the CPXgetstat()
+  // routine, for solutions to LP problems or mixed integer problems. If
+  // no solution exists, the return value is zero.
+  // For Simplex, Barrier
+  // 1          CPX_OPTIMAL
+  //          Optimal solution found
+  // 2          CPX_INFEASIBLE
+  //          Problem infeasible
+  // 3    CPX_UNBOUNDED
+  //          Problem unbounded
+  // 4          CPX_OBJ_LIM
+  //          Objective limit exceeded in Phase II
+  // 5          CPX_IT_LIM_FEAS
+  //          Iteration limit exceeded in Phase II
+  // 6          CPX_IT_LIM_INFEAS
+  //          Iteration limit exceeded in Phase I
+  // 7          CPX_TIME_LIM_FEAS
+  //          Time limit exceeded in Phase II
+  // 8          CPX_TIME_LIM_INFEAS
+  //          Time limit exceeded in Phase I
+  // 9          CPX_NUM_BEST_FEAS
+  //          Problem non-optimal, singularities in Phase II
+  // 10         CPX_NUM_BEST_INFEAS
+  //          Problem non-optimal, singularities in Phase I
+  // 11         CPX_OPTIMAL_INFEAS
+  //          Optimal solution found, unscaled infeasibilities
+  // 12         CPX_ABORT_FEAS
+  //          Aborted in Phase II
+  // 13         CPX_ABORT_INFEAS
+  //          Aborted in Phase I
+  // 14          CPX_ABORT_DUAL_INFEAS
+  //          Aborted in barrier, dual infeasible
+  // 15          CPX_ABORT_PRIM_INFEAS
+  //          Aborted in barrier, primal infeasible
+  // 16          CPX_ABORT_PRIM_DUAL_INFEAS
+  //          Aborted in barrier, primal and dual infeasible
+  // 17          CPX_ABORT_PRIM_DUAL_FEAS
+  //          Aborted in barrier, primal and dual feasible
+  // 18          CPX_ABORT_CROSSOVER
+  //          Aborted in crossover
+  // 19          CPX_INForUNBD
+  //          Infeasible or unbounded
+  // 20   CPX_PIVOT
+  //       User pivot used
+  //
+  //     Ezeket hova tegyem:
+  // ??case CPX_ABORT_DUAL_INFEAS
+  // ??case CPX_ABORT_CROSSOVER
+  // ??case CPX_INForUNBD
+  // ??case CPX_PIVOT
+  //Some more interesting stuff:
+  // CPX_PARAM_PROBMETHOD  1062  int  LPMETHOD
+  // 0 Automatic
+  // 1 Primal Simplex
+  // 2 Dual Simplex
+  // 3 Network Simplex
+  // 4 Standard Barrier
+  // Default: 0
+  // Description: Method for linear optimization.
+  // Determines which algorithm is used when CPXlpopt() (or "optimize"
+  // in the Interactive Optimizer) is called. Currently the behavior of
+  // the "Automatic" setting is that CPLEX simply invokes the dual
+  // simplex method, but this capability may be expanded in the future
+  // so that CPLEX chooses the method based on problem characteristics
 #if CPX_VERSION < 900
   void statusSwitch(CPXENVptr env,int& stat){
+  void statusSwitch(CPXENVptr cplexEnv(),int& stat){
     int lpmethod;
     CPXgetintparam (env,CPX_PARAM_LPMETHOD,&lpmethod);
+    CPXgetintparam (cplexEnv(),CPX_PARAM_PROBMETHOD,&lpmethod);
     if (lpmethod==2){
       if (stat==CPX_UNBOUNDED){
 …
 #endif
+  LpCplex::SolutionStatus LpCplex::_getPrimalStatus() const
+  {
+    //Unboundedness not treated well: the following is from cplex 9.0 doc
+  LpCplex::ProblemType LpCplex::_getPrimalType() const {
+    // Unboundedness not treated well: the following is from cplex 9.0 doc
     // About Unboundedness
 …
     // has a feasible solution.
     int stat = CPXgetstat(env, lp);
+    int stat = CPXgetstat(cplexEnv(), _prob);
 #if CPX_VERSION >= 800
     switch (stat)
+    {
+      {
       case CPX_STAT_OPTIMAL:
         return OPTIMAL;
       case CPX_STAT_UNBOUNDED:
         return INFINITE;
+        return UNBOUNDED;
       case CPX_STAT_INFEASIBLE:
         return INFEASIBLE;
       default:
         return UNDEFINED;
+    }
+      }
 #else
     statusSwitch(env,stat);
     //CPXgetstat(env, lp);
+    statusSwitch(cplexEnv(),stat);
+    //CPXgetstat(cplexEnv(), _prob);
     //printf("A primal status: %d, CPX_OPTIMAL=%d \n",stat,CPX_OPTIMAL);
     switch (stat) {
 …
     case CPX_UNBOUNDED://Unbounded
       return INFEASIBLE;//In case of dual simplex
       //return INFINITE;
+      //return UNBOUNDED;
     case CPX_INFEASIBLE://Infeasible
  //    case CPX_IT_LIM_INFEAS:
 //     case CPX_TIME_LIM_INFEAS:
 //     case CPX_NUM_BEST_INFEAS:
 //     case CPX_OPTIMAL_INFEAS:
 //     case CPX_ABORT_INFEAS:
 //     case CPX_ABORT_PRIM_INFEAS:
 //     case CPX_ABORT_PRIM_DUAL_INFEAS:
       return INFINITE;//In case of dual simplex
+      //    case CPX_IT_LIM_INFEAS:
+      //     case CPX_TIME_LIM_INFEAS:
+      //     case CPX_NUM_BEST_INFEAS:
+      //     case CPX_OPTIMAL_INFEAS:
+      //     case CPX_ABORT_INFEAS:
+      //     case CPX_ABORT_PRIM_INFEAS:
+      //     case CPX_ABORT_PRIM_DUAL_INFEAS:
+      return UNBOUNDED;//In case of dual simplex
       //return INFEASIBLE;
 //     case CPX_OBJ_LIM:
 //     case CPX_IT_LIM_FEAS:
 //     case CPX_TIME_LIM_FEAS:
 //     case CPX_NUM_BEST_FEAS:
 //     case CPX_ABORT_FEAS:
 //     case CPX_ABORT_PRIM_DUAL_FEAS:
 //       return FEASIBLE;
+      //     case CPX_OBJ_LIM:
+      //     case CPX_IT_LIM_FEAS:
+      //     case CPX_TIME_LIM_FEAS:
+      //     case CPX_NUM_BEST_FEAS:
+      //     case CPX_ABORT_FEAS:
+      //     case CPX_ABORT_PRIM_DUAL_FEAS:
+      //       return FEASIBLE;
     default:
       return UNDEFINED; //Everything else comes here
 …
+  }
+//9.0-as cplex verzio statusai
+// CPX_STAT_ABORT_DUAL_OBJ_LIM
+// CPX_STAT_ABORT_IT_LIM
+// CPX_STAT_ABORT_OBJ_LIM
+// CPX_STAT_ABORT_PRIM_OBJ_LIM
+// CPX_STAT_ABORT_TIME_LIM
+// CPX_STAT_ABORT_USER
+// CPX_STAT_FEASIBLE_RELAXED
+// CPX_STAT_INFEASIBLE
+// CPX_STAT_INForUNBD
+// CPX_STAT_NUM_BEST
+// CPX_STAT_OPTIMAL
+// CPX_STAT_OPTIMAL_FACE_UNBOUNDED
+// CPX_STAT_OPTIMAL_INFEAS
+// CPX_STAT_OPTIMAL_RELAXED
+// CPX_STAT_UNBOUNDED
+  LpCplex::SolutionStatus LpCplex::_getDualStatus() const
+  {
+    int stat = CPXgetstat(env, lp);
+  //9.0-as cplex verzio statusai
+  // CPX_STAT_ABORT_DUAL_OBJ_LIM
+  // CPX_STAT_ABORT_IT_LIM
+  // CPX_STAT_ABORT_OBJ_LIM
+  // CPX_STAT_ABORT_PRIM_OBJ_LIM
+  // CPX_STAT_ABORT_TIME_LIM
+  // CPX_STAT_ABORT_USER
+  // CPX_STAT_FEASIBLE_RELAXED
+  // CPX_STAT_INFEASIBLE
+  // CPX_STAT_INForUNBD
+  // CPX_STAT_NUM_BEST
+  // CPX_STAT_OPTIMAL
+  // CPX_STAT_OPTIMAL_FACE_UNBOUNDED
+  // CPX_STAT_OPTIMAL_INFEAS
+  // CPX_STAT_OPTIMAL_RELAXED
+  // CPX_STAT_UNBOUNDED
+  LpCplex::ProblemType LpCplex::_getDualType() const {
+    int stat = CPXgetstat(cplexEnv(), _prob);
 #if CPX_VERSION >= 800
+    switch (stat)
+    {
+      case CPX_STAT_OPTIMAL:
+        return OPTIMAL;
+      case CPX_STAT_UNBOUNDED:
+        return INFEASIBLE;
+      default:
+        return UNDEFINED;
+    switch (stat) {
+    case CPX_STAT_OPTIMAL:
+      return OPTIMAL;
+    case CPX_STAT_UNBOUNDED:
+      return INFEASIBLE;
+    default:
+      return UNDEFINED;
+    }
 #else
     statusSwitch(env,stat);
+    statusSwitch(cplexEnv(),stat);
     switch (stat) {
     case 0:
 …
       return OPTIMAL;
     case CPX_UNBOUNDED:
      return INFEASIBLE;
+      return INFEASIBLE;
     default:
       return UNDEFINED; //Everything else comes here
 …
+  }
+  LpCplex::ProblemTypes LpCplex::_getProblemType() const
+  {
+    int stat = CPXgetstat(env, lp);
+#if CPX_VERSION >= 800
+    switch (stat)
+    {
+      case CPX_STAT_OPTIMAL:
+        return PRIMAL_DUAL_FEASIBLE;
+      case CPX_STAT_UNBOUNDED:
+         return PRIMAL_FEASIBLE_DUAL_INFEASIBLE;
+      default:
+        return UNKNOWN;
+    }
+  // MipCplex members
+  MipCplex::MipCplex()
+    : LpBase(), CplexBase(), MipSolver() {
+#if CPX_VERSION < 800
+    CPXchgprobtype(cplexEnv(),  _prob, CPXPROB_MIP);
 #else
+    CPXchgprobtype(cplexEnv(),  _prob, CPXPROB_MILP);
+#endif
+  }
+  MipCplex::MipCplex(const CplexEnv& env)
+    : LpBase(), CplexBase(env), MipSolver() {
+#if CPX_VERSION < 800
+    CPXchgprobtype(cplexEnv(),  _prob, CPXPROB_MIP);
+#else
+    CPXchgprobtype(cplexEnv(),  _prob, CPXPROB_MILP);
+#endif
+  }
+  MipCplex::MipCplex(const MipCplex& other)
+    : LpBase(), CplexBase(other), MipSolver() {}
+  MipCplex::~MipCplex() {}
+  MipCplex* MipCplex::_newSolver() const { return new MipCplex; }
+  MipCplex* MipCplex::_cloneSolver() const {return new MipCplex(*this); }
+  const char* MipCplex::_solverName() const { return "MipCplex"; }
+  void MipCplex::_setColType(int i, MipCplex::ColTypes col_type) {
+    // Note If a variable is to be changed to binary, a call to CPXchgbds
+    // should also be made to change the bounds to 0 and 1.
+    switch (col_type){
+    case INTEGER: {
+      const char t = 'I';
+      CPXchgctype (cplexEnv(), _prob, 1, &i, &t);
+    } break;
+    case REAL: {
+      const char t = 'C';
+      CPXchgctype (cplexEnv(), _prob, 1, &i, &t);
+    } break;
+    default:
+      break;
+    }
+  }
+  MipCplex::ColTypes MipCplex::_getColType(int i) const {
+    char t;
+    CPXgetctype (cplexEnv(), _prob, &t, i, i);
+    switch (t) {
+    case 'I':
+      return INTEGER;
+    case 'C':
+      return REAL;
+    default:
+      LEMON_ASSERT(false, "Invalid column type");
+      return ColTypes();
+    }
+  }
+  MipCplex::SolveExitStatus MipCplex::_solve() {
+    int status;
+    status = CPXmipopt (cplexEnv(), _prob);
+    if (status==0)
+      return SOLVED;
+    else
+      return UNSOLVED;
+  }
+  MipCplex::ProblemType MipCplex::_getType() const {
+    int stat = CPXgetstat(cplexEnv(), _prob);
+    //Fortunately, MIP statuses did not change for cplex 8.0
     switch (stat) {
+    case CPX_OPTIMAL://Optimal
+        return PRIMAL_DUAL_FEASIBLE;
+    case CPX_UNBOUNDED:
+         return PRIMAL_FEASIBLE_DUAL_INFEASIBLE;
+//         return PRIMAL_INFEASIBLE_DUAL_FEASIBLE;
+//         return PRIMAL_DUAL_INFEASIBLE;
+//Seems to be that this is all we can say for sure
+    default:
+        //In all other cases
+        return UNKNOWN;
+      //FIXME error
+    }
+#endif
+  }
+  void LpCplex::_setMax()
+  {
+    CPXchgobjsen(env, lp, CPX_MAX);
+   }
+  void LpCplex::_setMin()
+  {
+    CPXchgobjsen(env, lp, CPX_MIN);
+   }
+  bool LpCplex::_isMax() const
+  {
+    if (CPXgetobjsen(env, lp)==CPX_MAX)
+      return true;
+    else
+      return false;
+    case CPXMIP_OPTIMAL:
+      // Optimal integer solution has been found.
+    case CPXMIP_OPTIMAL_TOL:
+      // Optimal soluton with the tolerance defined by epgap or epagap has
+      // been found.
+      return OPTIMAL;
+      //This also exists in later issues
+      //    case CPXMIP_UNBOUNDED:
+      //return UNBOUNDED;
+      case CPXMIP_INFEASIBLE:
+        return INFEASIBLE;
+    default:
+      return UNDEFINED;
+    }
+    //Unboundedness not treated well: the following is from cplex 9.0 doc
+    // About Unboundedness
+    // The treatment of models that are unbounded involves a few
+    // subtleties. Specifically, a declaration of unboundedness means that
+    // ILOG CPLEX has determined that the model has an unbounded
+    // ray. Given any feasible solution x with objective z, a multiple of
+    // the unbounded ray can be added to x to give a feasible solution
+    // with objective z-1 (or z+1 for maximization models). Thus, if a
+    // feasible solution exists, then the optimal objective is
+    // unbounded. Note that ILOG CPLEX has not necessarily concluded that
+    // a feasible solution exists. Users can call the routine CPXsolninfo
+    // to determine whether ILOG CPLEX has also concluded that the model
+    // has a feasible solution.
+  }
+  MipCplex::Value MipCplex::_getSol(int i) const {
+    Value x;
+    CPXgetmipx(cplexEnv(), _prob, &x, i, i);
+    return x;
+  }
+  MipCplex::Value MipCplex::_getSolValue() const {
+    Value objval;
+    CPXgetmipobjval(cplexEnv(), _prob, &objval);
+    return objval;
+  }

lemon/lp_cplex.h

-                      r458
+                      r459
 namespace lemon {
+  /// \brief Interface for the CPLEX solver
+  ///
+  /// This class implements an interface for the CPLEX LP solver.
+  class LpCplex :virtual public LpSolverBase {
+  public:
+    typedef LpSolverBase Parent;
+    /// \e
+    int status;
+    cpxenv* env;
+    cpxlp* lp;
+    /// \e
+    LpCplex();
+    /// \e
+    LpCplex(const LpCplex&);
+    /// \e
+    ~LpCplex();
+  protected:
+    virtual LpSolverBase* _newLp();
+    virtual LpSolverBase* _copyLp();
+  /// \brief Reference counted wrapper around cpxenv pointer
+  ///
+  /// The cplex uses environment object which is responsible for
+  /// checking the proper license usage. This class provides a simple
+  /// interface for share the environment object between different
+  /// problems.
+  class CplexEnv {
+    friend class CplexBase;
+  private:
+    cpxenv* _env;
+    mutable int* _cnt;
+  public:
+    /// \brief This exception is thrown when the license check is not
+    /// sufficient
+    class LicenseError : public Exception {
+      friend class CplexEnv;
+    private:
+      LicenseError(int status);
+      char _message[510];
+    public:
+      /// The short error message
+      virtual const char* what() const throw() {
+        return _message;
+      }
+    };
+    /// Constructor
+    CplexEnv();
+    /// Shallow copy constructor
+    CplexEnv(const CplexEnv&);
+    /// Shallow assignement
+    CplexEnv& operator=(const CplexEnv&);
+    /// Destructor
+    virtual ~CplexEnv();
+  protected:
+    cpxenv* cplexEnv() { return _env; }
+    const cpxenv* cplexEnv() const { return _env; }
+  };
+  /// \brief Base interface for the CPLEX LP and MIP solver
+  ///
+  /// This class implements the common interface of the CPLEX LP and
+  /// MIP solvers.
+  /// \ingroup lp_group
+  class CplexBase : virtual public LpBase {
+  protected:
+    CplexEnv _env;
+    cpxlp* _prob;
+    CplexBase();
+    CplexBase(const CplexEnv&);
+    CplexBase(const CplexBase &);
+    virtual ~CplexBase();
     virtual int _addCol();
     virtual int _addRow();
     virtual void _eraseCol(int i);
     virtual void _eraseRow(int i);
+    virtual void _getColName(int col, std::string & name) const;
+    virtual void _setColName(int col, const std::string & name);
+    virtual void _eraseColId(int i);
+    virtual void _eraseRowId(int i);
+    virtual void _getColName(int col, std::string& name) const;
+    virtual void _setColName(int col, const std::string& name);
     virtual int _colByName(const std::string& name) const;
+    virtual void _setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e);
+    virtual void _getRowCoeffs(int i, RowIterator b) const;
+    virtual void _setColCoeffs(int i, ConstColIterator b, ConstColIterator e);
+    virtual void _getColCoeffs(int i, ColIterator b) const;
+    virtual void _getRowName(int row, std::string& name) const;
+    virtual void _setRowName(int row, const std::string& name);
+    virtual int _rowByName(const std::string& name) const;
+    virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e);
+    virtual void _getRowCoeffs(int i, InsertIterator b) const;
+    virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e);
+    virtual void _getColCoeffs(int i, InsertIterator b) const;
     virtual void _setCoeff(int row, int col, Value value);
     virtual Value _getCoeff(int row, int col) const;
 …
     virtual void _setColLowerBound(int i, Value value);
     virtual Value _getColLowerBound(int i) const;
     virtual void _setColUpperBound(int i, Value value);
     virtual Value _getColUpperBound(int i) const;
+//     virtual void _setRowLowerBound(int i, Value value);
+//     virtual void _setRowUpperBound(int i, Value value);
+    virtual void _setRowBounds(int i, Value lower, Value upper);
+    virtual void _getRowBounds(int i, Value &lb, Value &ub) const;
+  private:
+    void _set_row_bounds(int i, Value lb, Value ub);
+  protected:
+    virtual void _setRowLowerBound(int i, Value value);
+    virtual Value _getRowLowerBound(int i) const;
+    virtual void _setRowUpperBound(int i, Value value);
+    virtual Value _getRowUpperBound(int i) const;
+    virtual void _setObjCoeffs(ExprIterator b, ExprIterator e);
+    virtual void _getObjCoeffs(InsertIterator b) const;
     virtual void _setObjCoeff(int i, Value obj_coef);
     virtual Value _getObjCoeff(int i) const;
+    virtual void _clearObj();
+    virtual void _setSense(Sense sense);
+    virtual Sense _getSense() const;
+    virtual void _clear();
+  public:
+    /// Returns the used \c CplexEnv instance
+    const CplexEnv& env() const { return _env; }
+    ///
+    const cpxenv* cplexEnv() const { return _env.cplexEnv(); }
+    cpxlp* cplexLp() { return _prob; }
+    const cpxlp* cplexLp() const { return _prob; }
+  };
+  /// \brief Interface for the CPLEX LP solver
+  ///
+  /// This class implements an interface for the CPLEX LP solver.
+  ///\ingroup lp_group
+  class LpCplex : public CplexBase, public LpSolver {
+  public:
+    /// \e
+    LpCplex();
+    /// \e
+    LpCplex(const CplexEnv&);
+    /// \e
+    LpCplex(const LpCplex&);
+    /// \e
+    virtual ~LpCplex();
+  private:
+    // these values cannot retrieved element by element
+    mutable std::vector<int> _col_status;
+    mutable std::vector<int> _row_status;
+    mutable std::vector<Value> _primal_ray;
+    mutable std::vector<Value> _dual_ray;
+    void _clear_temporals();
+    SolveExitStatus convertStatus(int status);
+  protected:
+    virtual LpCplex* _cloneSolver() const;
+    virtual LpCplex* _newSolver() const;
+    virtual const char* _solverName() const;
     virtual SolveExitStatus _solve();
 …
     virtual Value _getDual(int i) const;
     virtual Value _getPrimalValue() const;
+    virtual bool _isBasicCol(int i) const;
+    virtual SolutionStatus _getPrimalStatus() const;
+    virtual SolutionStatus _getDualStatus() const;
+    virtual ProblemTypes _getProblemType() const;
+    virtual void _setMax();
+    virtual void _setMin();
+    virtual bool _isMax() const;
+  public:
+    cpxenv* cplexEnv() { return env; }
+    cpxlp* cplexLp() { return lp; }
+  };
+    virtual VarStatus _getColStatus(int i) const;
+    virtual VarStatus _getRowStatus(int i) const;
+    virtual Value _getPrimalRay(int i) const;
+    virtual Value _getDualRay(int i) const;
+    virtual ProblemType _getPrimalType() const;
+    virtual ProblemType _getDualType() const;
+  public:
+    /// Solve with primal simplex method
+    SolveExitStatus solvePrimal();
+    /// Solve with dual simplex method
+    SolveExitStatus solveDual();
+    /// Solve with barrier method
+    SolveExitStatus solveBarrier();
+  };
+  /// \brief Interface for the CPLEX MIP solver
+  ///
+  /// This class implements an interface for the CPLEX MIP solver.
+  ///\ingroup lp_group
+  class MipCplex : public CplexBase, public MipSolver {
+  public:
+    /// \e
+    MipCplex();
+    /// \e
+    MipCplex(const CplexEnv&);
+    /// \e
+    MipCplex(const MipCplex&);
+    /// \e
+    virtual ~MipCplex();
+  protected:
+    virtual MipCplex* _cloneSolver() const;
+    virtual MipCplex* _newSolver() const;
+    virtual const char* _solverName() const;
+    virtual ColTypes _getColType(int col) const;
+    virtual void _setColType(int col, ColTypes col_type);
+    virtual SolveExitStatus _solve();
+    virtual ProblemType _getType() const;
+    virtual Value _getSol(int i) const;
+    virtual Value _getSolValue() const;
+  };
 } //END OF NAMESPACE LEMON

lemon/lp_glpk.cc

-                      r458
+                      r459
 ///\file
 ///\brief Implementation of the LEMON-GLPK lp solver interface.
+///\brief Implementation of the LEMON GLPK LP and MIP solver interface.
 #include <lemon/lp_glpk.h>
-//#include <iostream>
-extern "C" {
 #include <glpk.h>
+}
+#if GLP_MAJOR_VERSION > 4 || (GLP_MAJOR_VERSION == 4 && GLP_MINOR_VERSION > 15)
+#define LEMON_glp(func) (glp_##func)
+#define LEMON_lpx(func) (lpx_##func)
+#define LEMON_GLP(def) (GLP_##def)
+#define LEMON_LPX(def) (LPX_##def)
+#else
+#define LEMON_glp(func) (lpx_##func)
+#define LEMON_lpx(func) (lpx_##func)
+#define LEMON_GLP(def) (LPX_##def)
+#define LEMON_LPX(def) (LPX_##def)
+#endif
+#include <lemon/assert.h>
 namespace lemon {
+  LpGlpk::LpGlpk() : Parent() {
+    solved = false;
+    rows = _lp_bits::LpId(1);
+    cols = _lp_bits::LpId(1);
+    lp = LEMON_glp(create_prob)();
+    LEMON_glp(create_index)(lp);
+    messageLevel(0);
+  }
+  LpGlpk::LpGlpk(const LpGlpk &glp) : Parent() {
+    solved = false;
+    rows = _lp_bits::LpId(1);
+    cols = _lp_bits::LpId(1);
+    lp = LEMON_glp(create_prob)();
+    LEMON_glp(create_index)(lp);
+    messageLevel(0);
+    //Coefficient matrix, row bounds
+    LEMON_glp(add_rows)(lp, LEMON_glp(get_num_rows)(glp.lp));
+    LEMON_glp(add_cols)(lp, LEMON_glp(get_num_cols)(glp.lp));
+    int len;
+    std::vector<int> ind(1+LEMON_glp(get_num_cols)(glp.lp));
+    std::vector<Value> val(1+LEMON_glp(get_num_cols)(glp.lp));
+    for (int i=1;i<=LEMON_glp(get_num_rows)(glp.lp);++i)
+      {
+        len=LEMON_glp(get_mat_row)(glp.lp,i,&*ind.begin(),&*val.begin());
+        LEMON_glp(set_mat_row)(lp, i,len,&*ind.begin(),&*val.begin());
+        LEMON_glp(set_row_bnds)(lp,i,
+                                LEMON_glp(get_row_type)(glp.lp,i),
+                                LEMON_glp(get_row_lb)(glp.lp,i),
+                                LEMON_glp(get_row_ub)(glp.lp,i));
+      }
+    //Objective function, coloumn bounds
+    LEMON_glp(set_obj_dir)(lp, LEMON_glp(get_obj_dir)(glp.lp));
+    //Objectif function's constant term treated separately
+    LEMON_glp(set_obj_coef)(lp,0,LEMON_glp(get_obj_coef)(glp.lp,0));
+    for (int i=1;i<=LEMON_glp(get_num_cols)(glp.lp);++i)
+      {
+        LEMON_glp(set_obj_coef)(lp,i,
+                                LEMON_glp(get_obj_coef)(glp.lp,i));
+        LEMON_glp(set_col_bnds)(lp,i,
+                                LEMON_glp(get_col_type)(glp.lp,i),
+                                LEMON_glp(get_col_lb)(glp.lp,i),
+                                LEMON_glp(get_col_ub)(glp.lp,i));
+      }
+    rows = glp.rows;
+    cols = glp.cols;
+  }
+  LpGlpk::~LpGlpk() {
+    LEMON_glp(delete_prob)(lp);
+  }
+  int LpGlpk::_addCol() {
+    int i=LEMON_glp(add_cols)(lp, 1);
+    LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(FR), 0.0, 0.0);
+    solved = false;
+  // GlpkBase members
+  GlpkBase::GlpkBase() : LpBase() {
+    lp = glp_create_prob();
+    glp_create_index(lp);
+  }
+  GlpkBase::GlpkBase(const GlpkBase &other) : LpBase() {
+    lp = glp_create_prob();
+    glp_copy_prob(lp, other.lp, GLP_ON);
+    glp_create_index(lp);
+    rows = other.rows;
+    cols = other.cols;
+  }
+  GlpkBase::~GlpkBase() {
+    glp_delete_prob(lp);
+  }
+  int GlpkBase::_addCol() {
+    int i = glp_add_cols(lp, 1);
+    glp_set_col_bnds(lp, i, GLP_FR, 0.0, 0.0);
     return i;
+  }
+  ///\e
+  LpSolverBase* LpGlpk::_newLp()
+  {
+    LpGlpk* newlp = new LpGlpk;
+    return newlp;
+  }
+  ///\e
+  LpSolverBase* LpGlpk::_copyLp()
+  {
+    LpGlpk *newlp = new LpGlpk(*this);
+    return newlp;
+  }
+  int LpGlpk::_addRow() {
+    int i=LEMON_glp(add_rows)(lp, 1);
+    solved = false;
+  int GlpkBase::_addRow() {
+    int i = glp_add_rows(lp, 1);
+    glp_set_row_bnds(lp, i, GLP_FR, 0.0, 0.0);
     return i;
+  }
+  void LpGlpk::_eraseCol(int i) {
+  void GlpkBase::_eraseCol(int i) {
     int ca[2];
+    ca[1]=i;
+    LEMON_glp(del_cols)(lp, 1, ca);
+    solved = false;
+  }
+  void LpGlpk::_eraseRow(int i) {
+    ca[1] = i;
+    glp_del_cols(lp, 1, ca);
+  }
+  void GlpkBase::_eraseRow(int i) {
     int ra[2];
+    ra[1]=i;
+    LEMON_glp(del_rows)(lp, 1, ra);
+    solved = false;
+  }
+  void LpGlpk::_getColName(int c, std::string & name) const
+  {
+    const char *n = LEMON_glp(get_col_name)(lp,c);
+    name = n?n:"";
+  }
+  void LpGlpk::_setColName(int c, const std::string & name)
+  {
+    LEMON_glp(set_col_name)(lp,c,const_cast<char*>(name.c_str()));
+  }
+  int LpGlpk::_colByName(const std::string& name) const
+  {
+    int k = LEMON_glp(find_col)(lp, const_cast<char*>(name.c_str()));
+    ra[1] = i;
+    glp_del_rows(lp, 1, ra);
+  }
+  void GlpkBase::_eraseColId(int i) {
+    cols.eraseIndex(i);
+    cols.shiftIndices(i);
+  }
+  void GlpkBase::_eraseRowId(int i) {
+    rows.eraseIndex(i);
+    rows.shiftIndices(i);
+  }
+  void GlpkBase::_getColName(int c, std::string& name) const {
+    const char *str = glp_get_col_name(lp, c);
+    if (str) name = str;
+    else name.clear();
+  }
+  void GlpkBase::_setColName(int c, const std::string & name) {
+    glp_set_col_name(lp, c, const_cast<char*>(name.c_str()));
+  }
+  int GlpkBase::_colByName(const std::string& name) const {
+    int k = glp_find_col(lp, const_cast<char*>(name.c_str()));
     return k > 0 ? k : -1;
+  }
+  void LpGlpk::_setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e)
+  {
+    std::vector<int> indices;
+  void GlpkBase::_getRowName(int r, std::string& name) const {
+    const char *str = glp_get_row_name(lp, r);
+    if (str) name = str;
+    else name.clear();
+  }
+  void GlpkBase::_setRowName(int r, const std::string & name) {
+    glp_set_row_name(lp, r, const_cast<char*>(name.c_str()));
+  }
+  int GlpkBase::_rowByName(const std::string& name) const {
+    int k = glp_find_row(lp, const_cast<char*>(name.c_str()));
+    return k > 0 ? k : -1;
+  }
+  void GlpkBase::_setRowCoeffs(int i, ExprIterator b, ExprIterator e) {
+    std::vector<int> indexes;
     std::vector<Value> values;
     indices.push_back(0);
+    indexes.push_back(0);
     values.push_back(0);
     for(ConstRowIterator it=b; it!=e; ++it) {
       indices.push_back(it->first);
+    for(ExprIterator it = b; it != e; ++it) {
+      indexes.push_back(it->first);
       values.push_back(it->second);
+    }
+    LEMON_glp(set_mat_row)(lp, i, values.size() - 1,
+                                &indices[0], &values[0]);
+    solved = false;
+  }
+  void LpGlpk::_getRowCoeffs(int ix, RowIterator b) const
+  {
+    int length = LEMON_glp(get_mat_row)(lp, ix, 0, 0);
+    std::vector<int> indices(length + 1);
+    glp_set_mat_row(lp, i, values.size() - 1,
+                    &indexes.front(), &values.front());
+  }
+  void GlpkBase::_getRowCoeffs(int ix, InsertIterator b) const {
+    int length = glp_get_mat_row(lp, ix, 0, 0);
+    std::vector<int> indexes(length + 1);
     std::vector<Value> values(length + 1);
     LEMON_glp(get_mat_row)(lp, ix, &indices[0], &values[0]);
+    glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
     for (int i = 1; i <= length; ++i) {
       *b = std::make_pair(indices[i], values[i]);
+      *b = std::make_pair(indexes[i], values[i]);
       ++b;
+    }
+  }
+  void LpGlpk::_setColCoeffs(int ix, ConstColIterator b, ConstColIterator e) {
+    std::vector<int> indices;
+  void GlpkBase::_setColCoeffs(int ix, ExprIterator b,
+                                     ExprIterator e) {
+    std::vector<int> indexes;
     std::vector<Value> values;
     indices.push_back(0);
+    indexes.push_back(0);
     values.push_back(0);
     for(ConstColIterator it=b; it!=e; ++it) {
       indices.push_back(it->first);
+    for(ExprIterator it = b; it != e; ++it) {
+      indexes.push_back(it->first);
       values.push_back(it->second);
+    }
+    LEMON_glp(set_mat_col)(lp, ix, values.size() - 1,
+                                &indices[0], &values[0]);
+    solved = false;
+  }
+  void LpGlpk::_getColCoeffs(int ix, ColIterator b) const
+  {
+    int length = LEMON_glp(get_mat_col)(lp, ix, 0, 0);
+    std::vector<int> indices(length + 1);
+    glp_set_mat_col(lp, ix, values.size() - 1,
+                    &indexes.front(), &values.front());
+  }
+  void GlpkBase::_getColCoeffs(int ix, InsertIterator b) const {
+    int length = glp_get_mat_col(lp, ix, 0, 0);
+    std::vector<int> indexes(length + 1);
     std::vector<Value> values(length + 1);
     LEMON_glp(get_mat_col)(lp, ix, &indices[0], &values[0]);
     for (int i = 1; i <= length; ++i) {
       *b = std::make_pair(indices[i], values[i]);
+    glp_get_mat_col(lp, ix, &indexes.front(), &values.front());
+    for (int i = 1; i  <= length; ++i) {
+      *b = std::make_pair(indexes[i], values[i]);
       ++b;
+    }
+  }
+  void LpGlpk::_setCoeff(int ix, int jx, Value value)
+  {
+    if (LEMON_glp(get_num_cols)(lp) < LEMON_glp(get_num_rows)(lp)) {
+      int length=LEMON_glp(get_mat_row)(lp, ix, 0, 0);
+      std::vector<int> indices(length + 2);
+  void GlpkBase::_setCoeff(int ix, int jx, Value value) {
+    if (glp_get_num_cols(lp) < glp_get_num_rows(lp)) {
+      int length = glp_get_mat_row(lp, ix, 0, 0);
+      std::vector<int> indexes(length + 2);
       std::vector<Value> values(length + 2);
       LEMON_glp(get_mat_row)(lp, ix, &indices[0], &values[0]);
+      glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
       //The following code does not suppose that the elements of the
       //array indices are sorted
       bool found=false;
       for (int i = 1; i <= length; ++i) {
         if (indices[i]==jx){
           found=true;
           values[i]=value;
+      //array indexes are sorted
+      bool found = false;
+      for (int i = 1; i  <= length; ++i) {
+        if (indexes[i] == jx) {
+          found = true;
+          values[i] = value;
           break;
+        }
+      }
       if (!found){
+      if (!found) {
         ++length;
         indices[length]=jx;
         values[length]=value;
+      }
       LEMON_glp(set_mat_row)(lp, ix, length, &indices[0], &values[0]);
+        indexes[length] = jx;
+        values[length] = value;
+      }
+      glp_set_mat_row(lp, ix, length, &indexes.front(), &values.front());
     } else {
       int length=LEMON_glp(get_mat_col)(lp, jx, 0, 0);
       std::vector<int> indices(length + 2);
+      int length = glp_get_mat_col(lp, jx, 0, 0);
+      std::vector<int> indexes(length + 2);
       std::vector<Value> values(length + 2);
       LEMON_glp(get_mat_col)(lp, jx, &indices[0], &values[0]);
+      glp_get_mat_col(lp, jx, &indexes.front(), &values.front());
       //The following code does not suppose that the elements of the
       //array indices are sorted
       bool found=false;
+      //array indexes are sorted
+      bool found = false;
       for (int i = 1; i <= length; ++i) {
         if (indices[i]==ix){
           found=true;
           values[i]=value;
+        if (indexes[i] == ix) {
+          found = true;
+          values[i] = value;
           break;
+        }
+      }
       if (!found){
+      if (!found) {
         ++length;
+        indices[length]=ix;
+        values[length]=value;
+      }
+      LEMON_glp(set_mat_col)(lp, jx, length, &indices[0], &values[0]);
+    }
+    solved = false;
+  }
+  LpGlpk::Value LpGlpk::_getCoeff(int ix, int jx) const
+  {
+    int length=LEMON_glp(get_mat_row)(lp, ix, 0, 0);
+    std::vector<int> indices(length + 1);
+        indexes[length] = ix;
+        values[length] = value;
+      }
+      glp_set_mat_col(lp, jx, length, &indexes.front(), &values.front());
+    }
+  }
+  GlpkBase::Value GlpkBase::_getCoeff(int ix, int jx) const {
+    int length = glp_get_mat_row(lp, ix, 0, 0);
+    std::vector<int> indexes(length + 1);
     std::vector<Value> values(length + 1);
+    LEMON_glp(get_mat_row)(lp, ix, &indices[0], &values[0]);
+    //The following code does not suppose that the elements of the
+    //array indices are sorted
+    for (int i = 1; i <= length; ++i) {
+      if (indices[i]==jx){
+    glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
+    for (int i = 1; i  <= length; ++i) {
+      if (indexes[i] == jx) {
         return values[i];
+      }
+    }
     return 0;
+  }
+  void LpGlpk::_setColLowerBound(int i, Value lo)
+  {
+    if (lo==INF) {
+      //FIXME error
+    }
+    int b=LEMON_glp(get_col_type)(lp, i);
+    double up=LEMON_glp(get_col_ub)(lp, i);
+    if (lo==-INF) {
+  }
+  void GlpkBase::_setColLowerBound(int i, Value lo) {
+    LEMON_ASSERT(lo != INF, "Invalid bound");
+    int b = glp_get_col_type(lp, i);
+    double up = glp_get_col_ub(lp, i);
+    if (lo == -INF) {
       switch (b) {
       case LEMON_GLP(FR):
       case LEMON_GLP(LO):
         LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(FR), lo, up);
         break;
       case LEMON_GLP(UP):
         break;
       case LEMON_GLP(DB):
       case LEMON_GLP(FX):
         LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(UP), lo, up);
         break;
       default: ;
         //FIXME error
+      case GLP_FR:
+      case GLP_LO:
+        glp_set_col_bnds(lp, i, GLP_FR, lo, up);
+        break;
+      case GLP_UP:
+        break;
+      case GLP_DB:
+      case GLP_FX:
+        glp_set_col_bnds(lp, i, GLP_UP, lo, up);
+        break;
+      default:
+        break;
+      }
     } else {
       switch (b) {
       case LEMON_GLP(FR):
       case LEMON_GLP(LO):
         LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(LO), lo, up);
         break;
       case LEMON_GLP(UP):
       case LEMON_GLP(DB):
       case LEMON_GLP(FX):
         if (lo==up)
           LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(FX), lo, up);
+      case GLP_FR:
+      case GLP_LO:
+        glp_set_col_bnds(lp, i, GLP_LO, lo, up);
+        break;
+      case GLP_UP:
+      case GLP_DB:
+      case GLP_FX:
+        if (lo == up)
+          glp_set_col_bnds(lp, i, GLP_FX, lo, up);
         else
+          LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(DB), lo, up);
+        break;
+      default: ;
+        //FIXME error
+      }
+    }
+    solved = false;
+  }
+  LpGlpk::Value LpGlpk::_getColLowerBound(int i) const
+  {
+    int b=LEMON_glp(get_col_type)(lp, i);
+          glp_set_col_bnds(lp, i, GLP_DB, lo, up);
+        break;
+      default:
+        break;
+      }
+    }
+  }
+  GlpkBase::Value GlpkBase::_getColLowerBound(int i) const {
+    int b = glp_get_col_type(lp, i);
+    switch (b) {
+    case GLP_LO:
+    case GLP_DB:
+    case GLP_FX:
+      return glp_get_col_lb(lp, i);
+    default:
+      return -INF;
+    }
+  }
+  void GlpkBase::_setColUpperBound(int i, Value up) {
+    LEMON_ASSERT(up != -INF, "Invalid bound");
+    int b = glp_get_col_type(lp, i);
+    double lo = glp_get_col_lb(lp, i);
+    if (up == INF) {
       switch (b) {
+      case LEMON_GLP(LO):
+      case LEMON_GLP(DB):
+      case LEMON_GLP(FX):
+        return LEMON_glp(get_col_lb)(lp, i);
+      default: ;
+        return -INF;
+      }
+  }
+  void LpGlpk::_setColUpperBound(int i, Value up)
+  {
+    if (up==-INF) {
+      //FIXME error
+    }
+    int b=LEMON_glp(get_col_type)(lp, i);
+    double lo=LEMON_glp(get_col_lb)(lp, i);
+    if (up==INF) {
+      switch (b) {
+      case LEMON_GLP(FR):
+      case LEMON_GLP(LO):
+        break;
+      case LEMON_GLP(UP):
+        LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(FR), lo, up);
+        break;
+      case LEMON_GLP(DB):
+      case LEMON_GLP(FX):
+        LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(LO), lo, up);
+        break;
+      default: ;
+        //FIXME error
+      case GLP_FR:
+      case GLP_LO:
+        break;
+      case GLP_UP:
+        glp_set_col_bnds(lp, i, GLP_FR, lo, up);
+        break;
+      case GLP_DB:
+      case GLP_FX:
+        glp_set_col_bnds(lp, i, GLP_LO, lo, up);
+        break;
+      default:
+        break;
+      }
     } else {
       switch (b) {
       case LEMON_GLP(FR):
         LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(UP), lo, up);
         break;
       case LEMON_GLP(UP):
         LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(UP), lo, up);
         break;
       case LEMON_GLP(LO):
       case LEMON_GLP(DB):
       case LEMON_GLP(FX):
         if (lo==up)
           LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(FX), lo, up);
+      case GLP_FR:
+        glp_set_col_bnds(lp, i, GLP_UP, lo, up);
+        break;
+      case GLP_UP:
+        glp_set_col_bnds(lp, i, GLP_UP, lo, up);
+        break;
+      case GLP_LO:
+      case GLP_DB:
+      case GLP_FX:
+        if (lo == up)
+          glp_set_col_bnds(lp, i, GLP_FX, lo, up);
         else
+          LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(DB), lo, up);
+        break;
+      default: ;
+        //FIXME error
+      }
+    }
+    solved = false;
+  }
+  LpGlpk::Value LpGlpk::_getColUpperBound(int i) const
+  {
+    int b=LEMON_glp(get_col_type)(lp, i);
+          glp_set_col_bnds(lp, i, GLP_DB, lo, up);
+        break;
+      default:
+        break;
+      }
+    }
+  }
+  GlpkBase::Value GlpkBase::_getColUpperBound(int i) const {
+    int b = glp_get_col_type(lp, i);
       switch (b) {
       case LEMON_GLP(UP):
       case LEMON_GLP(DB):
       case LEMON_GLP(FX):
         return LEMON_glp(get_col_ub)(lp, i);
       default: ;
+      case GLP_UP:
+      case GLP_DB:
+      case GLP_FX:
+        return glp_get_col_ub(lp, i);
+      default:
         return INF;
+      }
+  }
+  void LpGlpk::_setRowBounds(int i, Value lb, Value ub)
+  {
+    //Bad parameter
+    if (lb==INF || ub==-INF) {
+      //FIXME error
+    }
+    if (lb == -INF){
+      if (ub == INF){
+        LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(FR), lb, ub);
+      }
+      else{
+        LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(UP), lb, ub);
+      }
+    }
+    else{
+      if (ub==INF){
+        LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(LO), lb, ub);
+      }
+      else{
+        if (lb == ub){
+          LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(FX), lb, ub);
+  void GlpkBase::_setRowLowerBound(int i, Value lo) {
+    LEMON_ASSERT(lo != INF, "Invalid bound");
+    int b = glp_get_row_type(lp, i);
+    double up = glp_get_row_ub(lp, i);
+    if (lo == -INF) {
+      switch (b) {
+      case GLP_FR:
+      case GLP_LO:
+        glp_set_row_bnds(lp, i, GLP_FR, lo, up);
+        break;
+      case GLP_UP:
+        break;
+      case GLP_DB:
+      case GLP_FX:
+        glp_set_row_bnds(lp, i, GLP_UP, lo, up);
+        break;
+      default:
+        break;
+      }
+    } else {
+      switch (b) {
+      case GLP_FR:
+      case GLP_LO:
+        glp_set_row_bnds(lp, i, GLP_LO, lo, up);
+        break;
+      case GLP_UP:
+      case GLP_DB:
+      case GLP_FX:
+        if (lo == up)
+          glp_set_row_bnds(lp, i, GLP_FX, lo, up);
+        else
+          glp_set_row_bnds(lp, i, GLP_DB, lo, up);
+        break;
+      default:
+        break;
+      }
+    }
+  }
+  GlpkBase::Value GlpkBase::_getRowLowerBound(int i) const {
+    int b = glp_get_row_type(lp, i);
+    switch (b) {
+    case GLP_LO:
+    case GLP_DB:
+    case GLP_FX:
+      return glp_get_row_lb(lp, i);
+    default:
+      return -INF;
+    }
+  }
+  void GlpkBase::_setRowUpperBound(int i, Value up) {
+    LEMON_ASSERT(up != -INF, "Invalid bound");
+    int b = glp_get_row_type(lp, i);
+    double lo = glp_get_row_lb(lp, i);
+    if (up == INF) {
+      switch (b) {
+      case GLP_FR:
+      case GLP_LO:
+        break;
+      case GLP_UP:
+        glp_set_row_bnds(lp, i, GLP_FR, lo, up);
+        break;
+      case GLP_DB:
+      case GLP_FX:
+        glp_set_row_bnds(lp, i, GLP_LO, lo, up);
+        break;
+      default:
+        break;
+      }
+    } else {
+      switch (b) {
+      case GLP_FR:
+        glp_set_row_bnds(lp, i, GLP_UP, lo, up);
+        break;
+      case GLP_UP:
+        glp_set_row_bnds(lp, i, GLP_UP, lo, up);
+        break;
+      case GLP_LO:
+      case GLP_DB:
+      case GLP_FX:
+        if (lo == up)
+          glp_set_row_bnds(lp, i, GLP_FX, lo, up);
+        else
+          glp_set_row_bnds(lp, i, GLP_DB, lo, up);
+        break;
+      default:
+        break;
+      }
+    }
+  }
+  GlpkBase::Value GlpkBase::_getRowUpperBound(int i) const {
+    int b = glp_get_row_type(lp, i);
+    switch (b) {
+    case GLP_UP:
+    case GLP_DB:
+    case GLP_FX:
+      return glp_get_row_ub(lp, i);
+    default:
+      return INF;
+    }
+  }
+  void GlpkBase::_setObjCoeffs(ExprIterator b, ExprIterator e) {
+    for (int i = 1; i <= glp_get_num_cols(lp); ++i) {
+      glp_set_obj_coef(lp, i, 0.0);
+    }
+    for (ExprIterator it = b; it != e; ++it) {
+      glp_set_obj_coef(lp, it->first, it->second);
+    }
+  }
+  void GlpkBase::_getObjCoeffs(InsertIterator b) const {
+    for (int i = 1; i <= glp_get_num_cols(lp); ++i) {
+      Value val = glp_get_obj_coef(lp, i);
+      if (val != 0.0) {
+        *b = std::make_pair(i, val);
+        ++b;
+      }
+    }
+  }
+  void GlpkBase::_setObjCoeff(int i, Value obj_coef) {
+    //i = 0 means the constant term (shift)
+    glp_set_obj_coef(lp, i, obj_coef);
+  }
+  GlpkBase::Value GlpkBase::_getObjCoeff(int i) const {
+    //i = 0 means the constant term (shift)
+    return glp_get_obj_coef(lp, i);
+  }
+  void GlpkBase::_setSense(GlpkBase::Sense sense) {
+    switch (sense) {
+    case MIN:
+      glp_set_obj_dir(lp, GLP_MIN);
+      break;
+    case MAX:
+      glp_set_obj_dir(lp, GLP_MAX);
+      break;
+    }
+  }
+  GlpkBase::Sense GlpkBase::_getSense() const {
+    switch(glp_get_obj_dir(lp)) {
+    case GLP_MIN:
+      return MIN;
+    case GLP_MAX:
+      return MAX;
+    default:
+      LEMON_ASSERT(false, "Wrong sense");
+      return GlpkBase::Sense();
+    }
+  }
+  void GlpkBase::_clear() {
+    glp_erase_prob(lp);
+    rows.clear();
+    cols.clear();
+  }
+  // LpGlpk members
+  LpGlpk::LpGlpk()
+    : LpBase(), GlpkBase(), LpSolver() {
+    messageLevel(MESSAGE_NO_OUTPUT);
+  }
+  LpGlpk::LpGlpk(const LpGlpk& other)
+    : LpBase(other), GlpkBase(other), LpSolver(other) {
+    messageLevel(MESSAGE_NO_OUTPUT);
+  }
+  LpGlpk* LpGlpk::_newSolver() const { return new LpGlpk; }
+  LpGlpk* LpGlpk::_cloneSolver() const { return new LpGlpk(*this); }
+  const char* LpGlpk::_solverName() const { return "LpGlpk"; }
+  void LpGlpk::_clear_temporals() {
+    _primal_ray.clear();
+    _dual_ray.clear();
+  }
+  LpGlpk::SolveExitStatus LpGlpk::_solve() {
+    return solvePrimal();
+  }
+  LpGlpk::SolveExitStatus LpGlpk::solvePrimal() {
+    _clear_temporals();
+    glp_smcp smcp;
+    glp_init_smcp(&smcp);
+    switch (_message_level) {
+    case MESSAGE_NO_OUTPUT:
+      smcp.msg_lev = GLP_MSG_OFF;
+      break;
+    case MESSAGE_ERROR_MESSAGE:
+      smcp.msg_lev = GLP_MSG_ERR;
+      break;
+    case MESSAGE_NORMAL_OUTPUT:
+      smcp.msg_lev = GLP_MSG_ON;
+      break;
+    case MESSAGE_FULL_OUTPUT:
+      smcp.msg_lev = GLP_MSG_ALL;
+      break;
+    }
+    if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
+    return SOLVED;
+  }
+  LpGlpk::SolveExitStatus LpGlpk::solveDual() {
+    _clear_temporals();
+    glp_smcp smcp;
+    glp_init_smcp(&smcp);
+    switch (_message_level) {
+    case MESSAGE_NO_OUTPUT:
+      smcp.msg_lev = GLP_MSG_OFF;
+      break;
+    case MESSAGE_ERROR_MESSAGE:
+      smcp.msg_lev = GLP_MSG_ERR;
+      break;
+    case MESSAGE_NORMAL_OUTPUT:
+      smcp.msg_lev = GLP_MSG_ON;
+      break;
+    case MESSAGE_FULL_OUTPUT:
+      smcp.msg_lev = GLP_MSG_ALL;
+      break;
+    }
+    smcp.meth = GLP_DUAL;
+    if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
+    return SOLVED;
+  }
+  LpGlpk::Value LpGlpk::_getPrimal(int i) const {
+    return glp_get_col_prim(lp, i);
+  }
+  LpGlpk::Value LpGlpk::_getDual(int i) const {
+    return glp_get_row_dual(lp, i);
+  }
+  LpGlpk::Value LpGlpk::_getPrimalValue() const {
+    return glp_get_obj_val(lp);
+  }
+  LpGlpk::VarStatus LpGlpk::_getColStatus(int i) const {
+    switch (glp_get_col_stat(lp, i)) {
+    case GLP_BS:
+      return BASIC;
+    case GLP_UP:
+      return UPPER;
+    case GLP_LO:
+      return LOWER;
+    case GLP_NF:
+      return FREE;
+    case GLP_NS:
+      return FIXED;
+    default:
+      LEMON_ASSERT(false, "Wrong column status");
+      return LpGlpk::VarStatus();
+    }
+  }
+  LpGlpk::VarStatus LpGlpk::_getRowStatus(int i) const {
+    switch (glp_get_row_stat(lp, i)) {
+    case GLP_BS:
+      return BASIC;
+    case GLP_UP:
+      return UPPER;
+    case GLP_LO:
+      return LOWER;
+    case GLP_NF:
+      return FREE;
+    case GLP_NS:
+      return FIXED;
+    default:
+      LEMON_ASSERT(false, "Wrong row status");
+      return LpGlpk::VarStatus();
+    }
+  }
+  LpGlpk::Value LpGlpk::_getPrimalRay(int i) const {
+    if (_primal_ray.empty()) {
+      int row_num = glp_get_num_rows(lp);
+      int col_num = glp_get_num_cols(lp);
+      _primal_ray.resize(col_num + 1, 0.0);
+      int index = glp_get_unbnd_ray(lp);
+      if (index != 0) {
+        // The primal ray is found in primal simplex second phase
+        LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :
+                      glp_get_col_stat(lp, index - row_num)) != GLP_BS,
+                     "Wrong primal ray");
+        bool negate = glp_get_obj_dir(lp) == GLP_MAX;
+        if (index > row_num) {
+          _primal_ray[index - row_num] = 1.0;
+          if (glp_get_col_dual(lp, index - row_num) > 0) {
+            negate = !negate;
+          }
+        } else {
+          if (glp_get_row_dual(lp, index) > 0) {
+            negate = !negate;
+          }
+        }
+        else{
+          LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(DB), lb, ub);
+        std::vector<int> ray_indexes(row_num + 1);
+        std::vector<Value> ray_values(row_num + 1);
+        int ray_length = glp_eval_tab_col(lp, index, &ray_indexes.front(),
+                                          &ray_values.front());
+        for (int i = 1; i <= ray_length; ++i) {
+          if (ray_indexes[i] > row_num) {
+            _primal_ray[ray_indexes[i] - row_num] = ray_values[i];
+          }
+        }
+      }
+    }
+    solved = false;
+  }
+  void LpGlpk::_getRowBounds(int i, Value &lb, Value &ub) const
+  {
+    int b=LEMON_glp(get_row_type)(lp, i);
+    switch (b) {
+    case LEMON_GLP(FR):
+    case LEMON_GLP(UP):
+      lb = -INF;
+        break;
+        if (negate) {
+          for (int i = 1; i <= col_num; ++i) {
+            _primal_ray[i] = - _primal_ray[i];
+          }
+        }
+      } else {
+        for (int i = 1; i <= col_num; ++i) {
+          _primal_ray[i] = glp_get_col_prim(lp, i);
+        }
+      }
+    }
+    return _primal_ray[i];
+  }
+  LpGlpk::Value LpGlpk::_getDualRay(int i) const {
+    if (_dual_ray.empty()) {
+      int row_num = glp_get_num_rows(lp);
+      _dual_ray.resize(row_num + 1, 0.0);
+      int index = glp_get_unbnd_ray(lp);
+      if (index != 0) {
+        // The dual ray is found in dual simplex second phase
+        LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :
+                      glp_get_col_stat(lp, index - row_num)) == GLP_BS,
+                     "Wrong dual ray");
+        int idx;
+        bool negate = false;
+        if (index > row_num) {
+          idx = glp_get_col_bind(lp, index - row_num);
+          if (glp_get_col_prim(lp, index - row_num) >
+              glp_get_col_ub(lp, index - row_num)) {
+            negate = true;
+          }
+        } else {
+          idx = glp_get_row_bind(lp, index);
+          if (glp_get_row_prim(lp, index) > glp_get_row_ub(lp, index)) {
+            negate = true;
+          }
+        }
+        _dual_ray[idx] = negate ?  - 1.0 : 1.0;
+        glp_btran(lp, &_dual_ray.front());
+      } else {
+        double eps = 1e-7;
+        // The dual ray is found in primal simplex first phase
+        // We assume that the glpk minimizes the slack to get feasible solution
+        for (int i = 1; i <= row_num; ++i) {
+          int index = glp_get_bhead(lp, i);
+          if (index <= row_num) {
+            double res = glp_get_row_prim(lp, index);
+            if (res > glp_get_row_ub(lp, index) + eps) {
+              _dual_ray[i] = -1;
+            } else if (res < glp_get_row_lb(lp, index) - eps) {
+              _dual_ray[i] = 1;
+            } else {
+              _dual_ray[i] = 0;
+            }
+            _dual_ray[i] *= glp_get_rii(lp, index);
+          } else {
+            double res = glp_get_col_prim(lp, index - row_num);
+            if (res > glp_get_col_ub(lp, index - row_num) + eps) {
+              _dual_ray[i] = -1;
+            } else if (res < glp_get_col_lb(lp, index - row_num) - eps) {
+              _dual_ray[i] = 1;
+            } else {
+              _dual_ray[i] = 0;
+            }
+            _dual_ray[i] /= glp_get_sjj(lp, index - row_num);
+          }
+        }
+        glp_btran(lp, &_dual_ray.front());
+        for (int i = 1; i <= row_num; ++i) {
+          _dual_ray[i] /= glp_get_rii(lp, i);
+        }
+      }
+    }
+    return _dual_ray[i];
+  }
+  LpGlpk::ProblemType LpGlpk::_getPrimalType() const {
+    if (glp_get_status(lp) == GLP_OPT)
+      return OPTIMAL;
+    switch (glp_get_prim_stat(lp)) {
+    case GLP_UNDEF:
+      return UNDEFINED;
+    case GLP_FEAS:
+    case GLP_INFEAS:
+      if (glp_get_dual_stat(lp) == GLP_NOFEAS) {
+        return UNBOUNDED;
+      } else {
+        return UNDEFINED;
+      }
+    case GLP_NOFEAS:
+      return INFEASIBLE;
     default:
+      lb=LEMON_glp(get_row_lb)(lp, i);
+    }
+    switch (b) {
+    case LEMON_GLP(FR):
+    case LEMON_GLP(LO):
+      ub = INF;
+        break;
+      LEMON_ASSERT(false, "Wrong primal type");
+      return  LpGlpk::ProblemType();
+    }
+  }
+  LpGlpk::ProblemType LpGlpk::_getDualType() const {
+    if (glp_get_status(lp) == GLP_OPT)
+      return OPTIMAL;
+    switch (glp_get_dual_stat(lp)) {
+    case GLP_UNDEF:
+      return UNDEFINED;
+    case GLP_FEAS:
+    case GLP_INFEAS:
+      if (glp_get_prim_stat(lp) == GLP_NOFEAS) {
+        return UNBOUNDED;
+      } else {
+        return UNDEFINED;
+      }
+    case GLP_NOFEAS:
+      return INFEASIBLE;
     default:
+      ub=LEMON_glp(get_row_ub)(lp, i);
+    }
+  }
+  void LpGlpk::_setObjCoeff(int i, Value obj_coef)
+  {
+    //i=0 means the constant term (shift)
+    LEMON_glp(set_obj_coef)(lp, i, obj_coef);
+    solved = false;
+  }
+  LpGlpk::Value LpGlpk::_getObjCoeff(int i) const {
+    //i=0 means the constant term (shift)
+    return LEMON_glp(get_obj_coef)(lp, i);
+  }
+  void LpGlpk::_clearObj()
+  {
+    for (int i=0;i<=LEMON_glp(get_num_cols)(lp);++i){
+      LEMON_glp(set_obj_coef)(lp, i, 0);
+    }
+    solved = false;
+  }
+  LpGlpk::SolveExitStatus LpGlpk::_solve()
+  {
+    // A way to check the problem to be solved
+    //LEMON_glp(write_cpxlp(lp,"naittvan.cpx");
+    LEMON_lpx(std_basis)(lp);
+    int i =  LEMON_lpx(simplex)(lp);
+    switch (i) {
+    case LEMON_LPX(E_OK):
+      solved = true;
+      return SOLVED;
+      LEMON_ASSERT(false, "Wrong primal type");
+      return  LpGlpk::ProblemType();
+    }
+  }
+  void LpGlpk::presolver(bool b) {
+    lpx_set_int_parm(lp, LPX_K_PRESOL, b ? 1 : 0);
+  }
+  void LpGlpk::messageLevel(MessageLevel m) {
+    _message_level = m;
+  }
+  // MipGlpk members
+  MipGlpk::MipGlpk()
+    : LpBase(), GlpkBase(), MipSolver() {
+    messageLevel(MESSAGE_NO_OUTPUT);
+  }
+  MipGlpk::MipGlpk(const MipGlpk& other)
+    : LpBase(), GlpkBase(other), MipSolver() {
+    messageLevel(MESSAGE_NO_OUTPUT);
+  }
+  void MipGlpk::_setColType(int i, MipGlpk::ColTypes col_type) {
+    switch (col_type) {
+    case INTEGER:
+      glp_set_col_kind(lp, i, GLP_IV);
+      break;
+    case REAL:
+      glp_set_col_kind(lp, i, GLP_CV);
+      break;
+    }
+  }
+  MipGlpk::ColTypes MipGlpk::_getColType(int i) const {
+    switch (glp_get_col_kind(lp, i)) {
+    case GLP_IV:
+    case GLP_BV:
+      return INTEGER;
     default:
+      return UNSOLVED;
+    }
+  }
+  LpGlpk::Value LpGlpk::_getPrimal(int i) const
+  {
+    return LEMON_glp(get_col_prim)(lp,i);
+  }
+  LpGlpk::Value LpGlpk::_getDual(int i) const
+  {
+    return LEMON_glp(get_row_dual)(lp,i);
+  }
+  LpGlpk::Value LpGlpk::_getPrimalValue() const
+  {
+    return LEMON_glp(get_obj_val)(lp);
+  }
+  bool LpGlpk::_isBasicCol(int i) const
+  {
+    return (LEMON_glp(get_col_stat)(lp, i)==LEMON_GLP(BS));
+  }
+  LpGlpk::SolutionStatus LpGlpk::_getPrimalStatus() const
+  {
+    if (!solved) return UNDEFINED;
+    int stat=  LEMON_lpx(get_status)(lp);
+    switch (stat) {
+    case LEMON_LPX(UNDEF)://Undefined (no solve has been run yet)
+      return UNDEFINED;
+    case LEMON_LPX(NOFEAS)://There is no feasible solution (primal, I guess)
+    case LEMON_LPX(INFEAS)://Infeasible
+      return INFEASIBLE;
+    case LEMON_LPX(UNBND)://Unbounded
+      return INFINITE;
+    case LEMON_LPX(FEAS)://Feasible
+      return FEASIBLE;
+    case LEMON_LPX(OPT)://Feasible
+      return OPTIMAL;
+    default:
+      return UNDEFINED; //to avoid gcc warning
+      //FIXME error
+    }
+  }
+  LpGlpk::SolutionStatus LpGlpk::_getDualStatus() const
+  {
+    if (!solved) return UNDEFINED;
+    switch (LEMON_lpx(get_dual_stat)(lp)) {
+    case LEMON_LPX(D_UNDEF)://Undefined (no solve has been run yet)
+      return UNDEFINED;
+    case LEMON_LPX(D_NOFEAS)://There is no dual feasible solution
+//    case LEMON_LPX(D_INFEAS://Infeasible
+      return INFEASIBLE;
+    case LEMON_LPX(D_FEAS)://Feasible
+      switch (LEMON_lpx(get_status)(lp)) {
+      case LEMON_LPX(NOFEAS):
+        return INFINITE;
+      case LEMON_LPX(OPT):
+      return REAL;
+    }
+  }
+  MipGlpk::SolveExitStatus MipGlpk::_solve() {
+    glp_smcp smcp;
+    glp_init_smcp(&smcp);
+    switch (_message_level) {
+    case MESSAGE_NO_OUTPUT:
+      smcp.msg_lev = GLP_MSG_OFF;
+      break;
+    case MESSAGE_ERROR_MESSAGE:
+      smcp.msg_lev = GLP_MSG_ERR;
+      break;
+    case MESSAGE_NORMAL_OUTPUT:
+      smcp.msg_lev = GLP_MSG_ON;
+      break;
+    case MESSAGE_FULL_OUTPUT:
+      smcp.msg_lev = GLP_MSG_ALL;
+      break;
+    }
+    smcp.meth = GLP_DUAL;
+    if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
+    if (glp_get_status(lp) != GLP_OPT) return SOLVED;
+    glp_iocp iocp;
+    glp_init_iocp(&iocp);
+    switch (_message_level) {
+    case MESSAGE_NO_OUTPUT:
+      iocp.msg_lev = GLP_MSG_OFF;
+      break;
+    case MESSAGE_ERROR_MESSAGE:
+      iocp.msg_lev = GLP_MSG_ERR;
+      break;
+    case MESSAGE_NORMAL_OUTPUT:
+      iocp.msg_lev = GLP_MSG_ON;
+      break;
+    case MESSAGE_FULL_OUTPUT:
+      iocp.msg_lev = GLP_MSG_ALL;
+      break;
+    }
+    if (glp_intopt(lp, &iocp) != 0) return UNSOLVED;
+    return SOLVED;
+  }
+  MipGlpk::ProblemType MipGlpk::_getType() const {
+    switch (glp_get_status(lp)) {
+    case GLP_OPT:
+      switch (glp_mip_status(lp)) {
+      case GLP_UNDEF:
+        return UNDEFINED;
+      case GLP_NOFEAS:
+        return INFEASIBLE;
+      case GLP_FEAS:
+        return FEASIBLE;
+      case GLP_OPT:
         return OPTIMAL;
       default:
+        return FEASIBLE;
+        LEMON_ASSERT(false, "Wrong problem type.");
+        return MipGlpk::ProblemType();
+      }
+    case GLP_NOFEAS:
+      return INFEASIBLE;
+    case GLP_INFEAS:
+    case GLP_FEAS:
+      if (glp_get_dual_stat(lp) == GLP_NOFEAS) {
+        return UNBOUNDED;
+      } else {
+        return UNDEFINED;
+      }
     default:
+      return UNDEFINED; //to avoid gcc warning
+      //FIXME error
+    }
+  }
+  LpGlpk::ProblemTypes LpGlpk::_getProblemType() const
+  {
+    if (!solved) return UNKNOWN;
+      //int stat=  LEMON_glp(get_status(lp);
+    int statp=  LEMON_lpx(get_prim_stat)(lp);
+    int statd=  LEMON_lpx(get_dual_stat)(lp);
+    if (statp==LEMON_LPX(P_FEAS) && statd==LEMON_LPX(D_FEAS))
+        return PRIMAL_DUAL_FEASIBLE;
+    if (statp==LEMON_LPX(P_FEAS) && statd==LEMON_LPX(D_NOFEAS))
+        return PRIMAL_FEASIBLE_DUAL_INFEASIBLE;
+    if (statp==LEMON_LPX(P_NOFEAS) && statd==LEMON_LPX(D_FEAS))
+        return PRIMAL_INFEASIBLE_DUAL_FEASIBLE;
+    if (statp==LEMON_LPX(P_NOFEAS) && statd==LEMON_LPX(D_NOFEAS))
+        return PRIMAL_DUAL_INFEASIBLE;
+    //In all other cases
+    return UNKNOWN;
+  }
+  void LpGlpk::_setMax()
+  {
+    solved = false;
+    LEMON_glp(set_obj_dir)(lp, LEMON_GLP(MAX));
+  }
+  void LpGlpk::_setMin()
+  {
+    solved = false;
+    LEMON_glp(set_obj_dir)(lp, LEMON_GLP(MIN));
+  }
+  bool LpGlpk::_isMax() const
+  {
+    return (LEMON_glp(get_obj_dir)(lp)==LEMON_GLP(MAX));
+  }
+  void LpGlpk::messageLevel(int m)
+  {
+    LEMON_lpx(set_int_parm)(lp, LEMON_LPX(K_MSGLEV), m);
+  }
+  void LpGlpk::presolver(bool b)
+  {
+    LEMON_lpx(set_int_parm)(lp, LEMON_LPX(K_PRESOL), b);
+  }
+      LEMON_ASSERT(false, "Wrong problem type.");
+      return MipGlpk::ProblemType();
+    }
+  }
+  MipGlpk::Value MipGlpk::_getSol(int i) const {
+    return glp_mip_col_val(lp, i);
+  }
+  MipGlpk::Value MipGlpk::_getSolValue() const {
+    return glp_mip_obj_val(lp);
+  }
+  MipGlpk* MipGlpk::_newSolver() const { return new MipGlpk; }
+  MipGlpk* MipGlpk::_cloneSolver() const {return new MipGlpk(*this); }
+  const char* MipGlpk::_solverName() const { return "MipGlpk"; }
+  void MipGlpk::messageLevel(MessageLevel m) {
+    _message_level = m;
+  }
 } //END OF NAMESPACE LEMON

lemon/lp_glpk.h

-                      r458
+                      r459
+  /// \brief Base interface for the GLPK LP and MIP solver
+  ///
+  /// This class implements the common interface of the GLPK LP and MIP solver.
+  /// \ingroup lp_group
+  class GlpkBase : virtual public LpBase {
+  protected:
+    typedef glp_prob LPX;
+    glp_prob* lp;
+    GlpkBase();
+    GlpkBase(const GlpkBase&);
+    virtual ~GlpkBase();
+  protected:
+    virtual int _addCol();
+    virtual int _addRow();
+    virtual void _eraseCol(int i);
+    virtual void _eraseRow(int i);
+    virtual void _eraseColId(int i);
+    virtual void _eraseRowId(int i);
+    virtual void _getColName(int col, std::string& name) const;
+    virtual void _setColName(int col, const std::string& name);
+    virtual int _colByName(const std::string& name) const;
+    virtual void _getRowName(int row, std::string& name) const;
+    virtual void _setRowName(int row, const std::string& name);
+    virtual int _rowByName(const std::string& name) const;
+    virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e);
+    virtual void _getRowCoeffs(int i, InsertIterator b) const;
+    virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e);
+    virtual void _getColCoeffs(int i, InsertIterator b) const;
+    virtual void _setCoeff(int row, int col, Value value);
+    virtual Value _getCoeff(int row, int col) const;
+    virtual void _setColLowerBound(int i, Value value);
+    virtual Value _getColLowerBound(int i) const;
+    virtual void _setColUpperBound(int i, Value value);
+    virtual Value _getColUpperBound(int i) const;
+    virtual void _setRowLowerBound(int i, Value value);
+    virtual Value _getRowLowerBound(int i) const;
+    virtual void _setRowUpperBound(int i, Value value);
+    virtual Value _getRowUpperBound(int i) const;
+    virtual void _setObjCoeffs(ExprIterator b, ExprIterator e);
+    virtual void _getObjCoeffs(InsertIterator b) const;
+    virtual void _setObjCoeff(int i, Value obj_coef);
+    virtual Value _getObjCoeff(int i) const;
+    virtual void _setSense(Sense);
+    virtual Sense _getSense() const;
+    virtual void _clear();
+  public:
+    ///Pointer to the underlying GLPK data structure.
+    LPX *lpx() {return lp;}
+    ///Const pointer to the underlying GLPK data structure.
+    const LPX *lpx() const {return lp;}
+    ///Returns the constraint identifier understood by GLPK.
+    int lpxRow(Row r) const { return rows(id(r)); }
+    ///Returns the variable identifier understood by GLPK.
+    int lpxCol(Col c) const { return cols(id(c)); }
+  };
   /// \brief Interface for the GLPK LP solver
   ///
   /// This class implements an interface for the GLPK LP solver.
   ///\ingroup lp_group
+  class LpGlpk : virtual public LpSolverBase {
+  protected:
+    typedef glp_prob LPX;
+    glp_prob* lp;
+    bool solved;
+  public:
+    typedef LpSolverBase Parent;
+  class LpGlpk : public GlpkBase, public LpSolver {
+  public:
+    ///\e
     LpGlpk();
+    LpGlpk(const LpGlpk &);
+    ~LpGlpk();
+  protected:
+    virtual LpSolverBase* _newLp();
+    virtual LpSolverBase* _copyLp();
+    virtual int _addCol();
+    virtual int _addRow();
+    virtual void _eraseCol(int i);
+    virtual void _eraseRow(int i);
+    virtual void _getColName(int col, std::string & name) const;
+    virtual void _setColName(int col, const std::string & name);
+    virtual int _colByName(const std::string& name) const;
+    virtual void _setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e);
+    virtual void _getRowCoeffs(int i, RowIterator b) const;
+    virtual void _setColCoeffs(int i, ConstColIterator b, ConstColIterator e);
+    virtual void _getColCoeffs(int i, ColIterator b) const;
+    virtual void _setCoeff(int row, int col, Value value);
+    virtual Value _getCoeff(int row, int col) const;
+    virtual void _setColLowerBound(int i, Value value);
+    virtual Value _getColLowerBound(int i) const;
+    virtual void _setColUpperBound(int i, Value value);
+    virtual Value _getColUpperBound(int i) const;
+    virtual void _setRowBounds(int i, Value lower, Value upper);
+    virtual void _getRowBounds(int i, Value &lb, Value &ub) const;
+    virtual void _setObjCoeff(int i, Value obj_coef);
+    virtual Value _getObjCoeff(int i) const;
+    virtual void _clearObj();
+    ///\e
+    ///\todo It should be clarified
+    ///
+    ///\e
+    LpGlpk(const LpGlpk&);
+  private:
+    mutable std::vector<double> _primal_ray;
+    mutable std::vector<double> _dual_ray;
+    void _clear_temporals();
+  protected:
+    virtual LpGlpk* _cloneSolver() const;
+    virtual LpGlpk* _newSolver() const;
+    virtual const char* _solverName() const;
     virtual SolveExitStatus _solve();
     virtual Value _getPrimal(int i) const;
     virtual Value _getDual(int i) const;
     virtual Value _getPrimalValue() const;
+    virtual bool _isBasicCol(int i) const;
+    ///\e
+    virtual VarStatus _getColStatus(int i) const;
+    virtual VarStatus _getRowStatus(int i) const;
+    virtual Value _getPrimalRay(int i) const;
+    virtual Value _getDualRay(int i) const;
     ///\todo It should be clarified
     ///
+    virtual SolutionStatus _getPrimalStatus() const;
+    virtual SolutionStatus _getDualStatus() const;
+    virtual ProblemTypes _getProblemType() const;
+    virtual void _setMax();
+    virtual void _setMin();
+    virtual bool _isMax() const;
+  public:
+    ///Set the verbosity of the messages
+    ///Set the verbosity of the messages
+    ///
+    ///\param m is the level of the messages output by the solver routines.
+    ///The possible values are:
+    ///- 0 --- no output (default value)
+    ///- 1 --- error messages only
+    ///- 2 --- normal output
+    ///- 3 --- full output (includes informational messages)
+    void messageLevel(int m);
+    virtual ProblemType _getPrimalType() const;
+    virtual ProblemType _getDualType() const;
+  public:
+    ///Solve with primal simplex
+    SolveExitStatus solvePrimal();
+    ///Solve with dual simplex
+    SolveExitStatus solveDual();
     ///Turns on or off the presolver
 …
     void presolver(bool b);
+    ///Pointer to the underlying GLPK data structure.
+    LPX *lpx() {return lp;}
+    ///Returns the constraint identifier understood by GLPK.
+    int lpxRow(Row r) { return _lpId(r); }
+    ///Returns the variable identifier understood by GLPK.
+    int lpxCol(Col c) { return _lpId(c); }
+    ///Enum for \c messageLevel() parameter
+    enum MessageLevel {
+      /// no output (default value)
+      MESSAGE_NO_OUTPUT = 0,
+      /// error messages only
+      MESSAGE_ERROR_MESSAGE = 1,
+      /// normal output
+      MESSAGE_NORMAL_OUTPUT = 2,
+      /// full output (includes informational messages)
+      MESSAGE_FULL_OUTPUT = 3
+    };
+  private:
+    MessageLevel _message_level;
+  public:
+    ///Set the verbosity of the messages
+    ///Set the verbosity of the messages
+    ///
+    ///\param m is the level of the messages output by the solver routines.
+    void messageLevel(MessageLevel m);
   };
+  /// \brief Interface for the GLPK MIP solver
+  ///
+  /// This class implements an interface for the GLPK MIP solver.
+  ///\ingroup lp_group
+  class MipGlpk : public GlpkBase, public MipSolver {
+  public:
+    ///\e
+    MipGlpk();
+    ///\e
+    MipGlpk(const MipGlpk&);
+  protected:
+    virtual MipGlpk* _cloneSolver() const;
+    virtual MipGlpk* _newSolver() const;
+    virtual const char* _solverName() const;
+    virtual ColTypes _getColType(int col) const;
+    virtual void _setColType(int col, ColTypes col_type);
+    virtual SolveExitStatus _solve();
+    virtual ProblemType _getType() const;
+    virtual Value _getSol(int i) const;
+    virtual Value _getSolValue() const;
+    ///Enum for \c messageLevel() parameter
+    enum MessageLevel {
+      /// no output (default value)
+      MESSAGE_NO_OUTPUT = 0,
+      /// error messages only
+      MESSAGE_ERROR_MESSAGE = 1,
+      /// normal output
+      MESSAGE_NORMAL_OUTPUT = 2,
+      /// full output (includes informational messages)
+      MESSAGE_FULL_OUTPUT = 3
+    };
+  private:
+    MessageLevel _message_level;
+  public:
+    ///Set the verbosity of the messages
+    ///Set the verbosity of the messages
+    ///
+    ///\param m is the level of the messages output by the solver routines.
+    void messageLevel(MessageLevel m);
+  };
 } //END OF NAMESPACE LEMON