# HG changeset patch
# User Balazs Dezso <deba@inf.elte.hu>
# Date 1228250433 -3600
# Node ID 7afc121e068927b0503e6c0cdcc58eb381ec50b1
# Parent  69928a704ffb1263a513409e81a4704fd2b7c754
Port LP and MIP solvers from SVN -r3509 (#44)

diff -r 69928a704ffb -r 7afc121e0689 configure.ac
--- a/configure.ac	Mon Jan 12 12:22:11 2009 +0000
+++ b/configure.ac	Tue Dec 02 21:40:33 2008 +0100
@@ -50,9 +50,9 @@
 AC_SUBST([WARNINGCXXFLAGS])
 
 dnl Checks for libraries.
-#LX_CHECK_GLPK
-#LX_CHECK_CPLEX
-#LX_CHECK_SOPLEX
+LX_CHECK_GLPK
+LX_CHECK_CPLEX
+LX_CHECK_SOPLEX
 
 AM_CONDITIONAL([HAVE_LP], [test x"$lx_lp_found" = x"yes"])
 AM_CONDITIONAL([HAVE_MIP], [test x"$lx_mip_found" = x"yes"])
@@ -117,10 +117,10 @@
 echo C++ compiler.................. : $CXX
 echo C++ compiles flags............ : $WARNINGCXXFLAGS $CXXFLAGS
 echo
-#echo GLPK support.................. : $lx_glpk_found
-#echo CPLEX support................. : $lx_cplex_found
-#echo SOPLEX support................ : $lx_soplex_found
-#echo
+echo GLPK support.................. : $lx_glpk_found
+echo CPLEX support................. : $lx_cplex_found
+echo SOPLEX support................ : $lx_soplex_found
+echo
 echo Build demo programs........... : $enable_demo
 echo Build additional tools........ : $enable_tools
 echo
diff -r 69928a704ffb -r 7afc121e0689 lemon/CMakeLists.txt
--- a/lemon/CMakeLists.txt	Mon Jan 12 12:22:11 2009 +0000
+++ b/lemon/CMakeLists.txt	Tue Dec 02 21:40:33 2008 +0100
@@ -4,6 +4,9 @@
   arg_parser.cc
   base.cc
   color.cc
+  lp_base.cc
+  lp_skeleton.cc
+  lp_utils.cc
   random.cc)
 
 INSTALL(
diff -r 69928a704ffb -r 7afc121e0689 lemon/Makefile.am
--- a/lemon/Makefile.am	Mon Jan 12 12:22:11 2009 +0000
+++ b/lemon/Makefile.am	Tue Dec 02 21:40:33 2008 +0100
@@ -10,10 +10,32 @@
 	lemon/arg_parser.cc \
 	lemon/base.cc \
 	lemon/color.cc \
+	lemon/lp_base.cc \
+	lemon/lp_skeleton.cc \
 	lemon/random.cc
 
-#lemon_libemon_la_CXXFLAGS = $(GLPK_CFLAGS) $(CPLEX_CFLAGS) $(SOPLEX_CXXFLAGS) $(AM_CXXFLAGS)
-#lemon_libemon_la_LDFLAGS = $(GLPK_LIBS) $(CPLEX_LIBS) $(SOPLEX_LIBS)
+
+lemon_libemon_la_CXXFLAGS = \
+	$(GLPK_CFLAGS) \
+	$(CPLEX_CFLAGS) \
+	$(SOPLEX_CXXFLAGS)
+
+lemon_libemon_la_LDFLAGS = \
+	$(GLPK_LIBS) \
+	$(CPLEX_LIBS) \
+	$(SOPLEX_LIBS)
+
+if HAVE_GLPK
+lemon_libemon_la_SOURCES += lemon/lp_glpk.cc lemon/mip_glpk.cc
+endif
+
+if HAVE_CPLEX
+lemon_libemon_la_SOURCES += lemon/lp_cplex.cc lemon/mip_cplex.cc
+endif
+
+if HAVE_SOPLEX
+lemon_libemon_la_SOURCES += lemon/lp_soplex.cc
+endif
 
 lemon_HEADERS += \
 	lemon/adaptors.h \
@@ -41,6 +63,14 @@
 	lemon/lgf_reader.h \
 	lemon/lgf_writer.h \
 	lemon/list_graph.h \
+	lemon/lp.h \
+	lemon/lp_base.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/maps.h \
 	lemon/math.h \
 	lemon/max_matching.h \
@@ -64,6 +94,7 @@
 	lemon/bits/enable_if.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/traits.h \
diff -r 69928a704ffb -r 7afc121e0689 lemon/bits/lp_id.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/bits/lp_id.h	Tue Dec 02 21:40:33 2008 +0100
@@ -0,0 +1,157 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_BITS_LP_SOLVER_ID_H
+#define LEMON_BITS_LP_SOLVER_ID_H
+
+namespace lemon {
+
+  namespace _lp_bits {
+
+    struct LpIdImpl {
+      std::vector<int> index;
+      std::vector<int> cross;
+      int first_index;
+      int first_free;
+    };
+
+    class LpId {
+    public:
+
+      class IdHandler {
+      public:
+        virtual ~IdHandler() {}
+        virtual int addId(LpIdImpl&) = 0;
+        virtual void eraseId(LpIdImpl&, int xn) = 0;
+      };
+
+      LpId(int min_index = 0) {
+        id_handler = 0;
+        impl.first_free = -1;
+        impl.first_index = min_index;
+        impl.cross.resize(impl.first_index);
+      }
+
+      LpId(const LpId& li) {
+        id_handler = 0;
+        impl = li.impl;
+      }
+
+      LpId& operator=(const LpId& li) {
+        id_handler = 0;
+        impl = li.impl;
+        return *this;
+      }
+
+      void setIdHandler(IdHandler& ih) {
+        id_handler = &ih;
+      }
+
+      int fixId(int fn) const {return impl.cross[fn];}
+      int floatingId(int xn) const {return impl.index[xn];}
+
+      int addId() {
+        if (id_handler == 0) {
+          int xn, fn = impl.cross.size();
+          if (impl.first_free == -1) {
+            xn = impl.index.size();
+            impl.index.push_back(fn);
+          } else {
+            xn = impl.first_free;
+            impl.first_free = impl.index[impl.first_free];
+            impl.index[xn] = fn;
+          }
+          impl.cross.push_back(xn);
+          return xn;
+        } else {
+          return id_handler->addId(impl);
+        }
+      }
+
+      void eraseId(int xn) {
+        if (id_handler == 0) {
+          int fn = impl.index[xn];
+          impl.index[xn] = impl.first_free;
+          impl.first_free = xn;
+          for(int i = fn + 1; i < int(impl.cross.size()); ++i) {
+            impl.cross[i - 1] = impl.cross[i];
+            impl.index[impl.cross[i]]--;
+          }
+          impl.cross.pop_back();
+        } else {
+          id_handler->eraseId(impl, xn);
+        }
+      }
+
+      void firstFloating(int& fn) const {
+        fn = impl.first_index;
+        if (fn == int(impl.cross.size())) fn = -1;
+      }
+
+      void nextFloating(int& fn) const {
+        ++fn;
+        if (fn == int(impl.cross.size())) fn = -1;
+      }
+
+      void firstFix(int& xn) const {
+        int fn;
+        firstFloating(fn);
+        xn = fn != -1 ? fixId(fn) : -1;
+      }
+
+      void nextFix(int& xn) const {
+        int fn = floatingId(xn);
+        nextFloating(fn);
+        xn = fn != -1 ? fixId(fn) : -1;
+      }
+
+    protected:
+      LpIdImpl impl;
+      IdHandler *id_handler;
+    };
+
+    class RelocateIdHandler : public LpId::IdHandler {
+    public:
+
+      virtual int addId(LpIdImpl& impl) {
+        int xn, fn = impl.cross.size();
+        if (impl.first_free == -1) {
+          xn = impl.index.size();
+          impl.index.push_back(fn);
+        } else {
+          xn = impl.first_free;
+          impl.first_free = impl.index[impl.first_free];
+          impl.index[xn] = fn;
+        }
+        impl.cross.push_back(xn);
+        return xn;
+      }
+
+      virtual void eraseId(LpIdImpl& impl, int xn) {
+        int fn = impl.index[xn];
+        impl.index[xn] = impl.first_free;
+        impl.first_free = xn;
+        impl.cross[fn] = impl.cross.back();
+        impl.index[impl.cross.back()] = fn;
+        impl.cross.pop_back();
+      }
+    };
+  }
+}
+
+#endif
diff -r 69928a704ffb -r 7afc121e0689 lemon/config.h.in
--- a/lemon/config.h.in	Mon Jan 12 12:22:11 2009 +0000
+++ b/lemon/config.h.in	Tue Dec 02 21:40:33 2008 +0100
@@ -9,3 +9,6 @@
 
 /* Define to 1 if you have GLPK. */
 #undef HAVE_GLPK
+
+/* Define to 1 if you have SOPLEX */
+#undef HAVE_SOPLEX
\ No newline at end of file
diff -r 69928a704ffb -r 7afc121e0689 lemon/lp.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/lp.h	Tue Dec 02 21:40:33 2008 +0100
@@ -0,0 +1,90 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_LP_H
+#define LEMON_LP_H
+
+#include<lemon/config.h>
+
+
+#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>
+#endif
+
+///\file
+///\brief Defines a default LP solver
+///\ingroup lp_group
+namespace lemon {
+
+#ifdef DOXYGEN
+  ///The default LP solver identifier
+
+  ///The default LP solver identifier.
+  ///\ingroup lp_group
+  ///
+  ///Currently, the possible values are \c GLPK or \c CPLEX
+#define DEFAULT_LP SOLVER
+  ///The default LP solver
+
+  ///The default LP solver.
+  ///\ingroup lp_group
+  ///
+  ///Currently, it is either \c LpGlpk or \c LpCplex
+  typedef LpGlpk Lp;
+  ///The default LP solver identifier string
+
+  ///The default LP solver identifier string.
+  ///\ingroup lp_group
+  ///
+  ///Currently, the possible values are "GLPK" or "CPLEX"
+  const char default_solver_name[]="SOLVER";
+
+  ///The default ILP solver.
+
+  ///The default ILP solver.
+  ///\ingroup lp_group
+  ///
+  ///Currently, it is either \c LpGlpk or \c LpCplex
+  typedef MipGlpk Mip;
+#else
+#ifdef HAVE_GLPK
+#define DEFAULT_LP GLPK
+  typedef LpGlpk Lp;
+  typedef MipGlpk Mip;
+  const char default_solver_name[]="GLPK";
+#elif HAVE_CPLEX
+#define DEFAULT_LP CPLEX
+  typedef LpCplex Lp;
+  typedef MipCplex Mip;
+  const char default_solver_name[]="CPLEX";
+#elif HAVE_SOPLEX
+#define DEFAULT_LP SOPLEX
+  typedef LpSoplex Lp;
+  const char default_solver_name[]="SOPLEX";
+#endif
+#endif
+
+} //namespace lemon
+
+#endif //LEMON_LP_H
diff -r 69928a704ffb -r 7afc121e0689 lemon/lp_base.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/lp_base.cc	Tue Dec 02 21:40:33 2008 +0100
@@ -0,0 +1,35 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+///\file
+///\brief The implementation of the LP solver interface.
+
+#include <lemon/lp_base.h>
+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();
+
+} //namespace lemon
diff -r 69928a704ffb -r 7afc121e0689 lemon/lp_base.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/lp_base.h	Tue Dec 02 21:40:33 2008 +0100
@@ -0,0 +1,1705 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_LP_BASE_H
+#define LEMON_LP_BASE_H
+
+#include<iostream>
+#include<vector>
+#include<map>
+#include<limits>
+#include<lemon/math.h>
+
+#include<lemon/core.h>
+#include<lemon/bits/lp_id.h>
+
+///\file
+///\brief The interface of the LP solver interface.
+///\ingroup lp_group
+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
+  ///\ingroup lp_group
+  class LpSolverBase {
+
+  protected:
+
+    _lp_bits::LpId rows;
+    _lp_bits::LpId cols;
+
+  public:
+
+    ///Possible outcomes of an LP solving procedure
+    enum SolveExitStatus {
+      ///This means that the problem has been successfully solved: either
+      ///an optimal solution has been found or infeasibility/unboundedness
+      ///has been proved.
+      SOLVED = 0,
+      ///Any other case (including the case when some user specified
+      ///limit has been exceeded)
+      UNSOLVED = 1
+    };
+
+      ///\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
+    };
+
+    ///The floating point type used by the solver
+    typedef double Value;
+    ///The infinity constant
+    static const Value INF;
+    ///The not a number constant
+    static const Value NaN;
+
+    static inline bool isNaN(const Value& v) { return v!=v; }
+
+    friend class Col;
+    friend class ColIt;
+    friend class Row;
+
+    ///Refer to a column of the LP.
+
+    ///This type is used to refer to a column of the LP.
+    ///
+    ///Its value remains valid and correct even after the addition or erase of
+    ///other columns.
+    ///
+    ///\todo Document what can one do with a Col (INVALID, comparing,
+    ///it is similar to Node/Edge)
+    class Col {
+    protected:
+      int id;
+      friend class LpSolverBase;
+      friend class MipSolverBase;
+      explicit Col(int _id) : id(_id) {}
+    public:
+      typedef Value ExprValue;
+      typedef True LpSolverCol;
+      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;}
+    };
+
+    class ColIt : public Col {
+      const LpSolverBase *_lp;
+    public:
+      ColIt() {}
+      ColIt(const LpSolverBase &lp) : _lp(&lp)
+      {
+        _lp->cols.firstFix(id);
+      }
+      ColIt(const Invalid&) : Col(INVALID) {}
+      ColIt &operator++()
+      {
+        _lp->cols.nextFix(id);
+        return *this;
+      }
+    };
+
+    static int id(const Col& col) { return col.id; }
+
+
+    ///Refer to a row of the LP.
+
+    ///This type is used to refer to a row of the LP.
+    ///
+    ///Its value remains valid and correct even after the addition or erase of
+    ///other rows.
+    ///
+    ///\todo Document what can one do with a Row (INVALID, comparing,
+    ///it is similar to Node/Edge)
+    class Row {
+    protected:
+      int id;
+      friend class LpSolverBase;
+      explicit Row(int _id) : id(_id) {}
+    public:
+      typedef Value ExprValue;
+      typedef True LpSolverRow;
+      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;}
+    };
+
+    class RowIt : public Row {
+      const LpSolverBase *_lp;
+    public:
+      RowIt() {}
+      RowIt(const LpSolverBase &lp) : _lp(&lp)
+      {
+        _lp->rows.firstFix(id);
+      }
+      RowIt(const Invalid&) : Row(INVALID) {}
+      RowIt &operator++()
+      {
+        _lp->rows.nextFix(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));
+    }
+
+
+  public:
+
+    ///Linear expression of variables and a constant component
+
+    ///This data structure stores a linear expression of the variables
+    ///(\ref Col "Col"s) and also has a constant component.
+    ///
+    ///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;
+    ///e[v]+=12;
+    ///e.erase(v);
+    ///\endcode
+    ///or you can also iterate through its elements.
+    ///\code
+    ///double s=0;
+    ///for(LpSolverBase::Expr::iterator i=e.begin();i!=e.end();++i)
+    ///  s+=i->second;
+    ///\endcode
+    ///(This code computes the sum of all coefficients).
+    ///- Numbers (<tt>double</tt>'s)
+    ///and variables (\ref Col "Col"s) directly convert to an
+    ///\ref Expr and the usual linear operations are defined, so
+    ///\code
+    ///v+w
+    ///2*v-3.12*(v-w/2)+2
+    ///v*2.1+(3*v+(v*12+w+6)*3)/2
+    ///\endcode
+    ///are valid \ref Expr "Expr"essions.
+    ///The usual assignment operations are also defined.
+    ///\code
+    ///e=v+w;
+    ///e+=2*v-3.12*(v-w/2)+2;
+    ///e*=3.4;
+    ///e/=5;
+    ///\endcode
+    ///- The constant member can be set and read by \ref constComp()
+    ///\code
+    ///e.constComp()=12;
+    ///double c=e.constComp();
+    ///\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>
+    {
+    public:
+      typedef LpSolverBase::Col Key;
+      typedef LpSolverBase::Value Value;
+
+    protected:
+      typedef std::map<Col,Value> Base;
+
+      Value const_comp;
+    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
+      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;
+        }
+      }
+
+      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;
+        }
+      }
+
+      ///Sets all coefficients and the constant component to 0.
+      void clear() {
+        Base::clear();
+        const_comp=0;
+      }
+
+      ///\e
+      Expr &operator+=(const Expr &e) {
+        for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
+          (*this)[j->first]+=j->second;
+        const_comp+=e.const_comp;
+        return *this;
+      }
+      ///\e
+      Expr &operator-=(const Expr &e) {
+        for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
+          (*this)[j->first]-=j->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;
+        return *this;
+      }
+      ///\e
+      Expr &operator/=(const Value &c) {
+        for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
+          j->second/=c;
+        const_comp/=c;
+        return *this;
+      }
+
+    };
+
+    ///Linear constraint
+
+    ///This data stucture represents a linear constraint in the LP.
+    ///Basically it is a linear expression with a lower or an upper bound
+    ///(or both). These parts of the constraint can be obtained by the member
+    ///functions \ref expr(), \ref lowerBound() and \ref upperBound(),
+    ///respectively.
+    ///There are two ways to construct a constraint.
+    ///- You can set the linear expression and the bounds directly
+    ///  by the functions above.
+    ///- The operators <tt>\<=</tt>, <tt>==</tt> and  <tt>\>=</tt>
+    ///  are defined between expressions, or even between constraints whenever
+    ///  it makes sense. Therefore if \c e and \c f are linear expressions and
+    ///  \c s and \c t are numbers, then the followings are valid expressions
+    ///  and thus they can be used directly e.g. in \ref addRow() whenever
+    ///  it makes sense.
+    ///\code
+    ///  e<=s
+    ///  e<=f
+    ///  e==f
+    ///  s<=e<=t
+    ///  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.
+    class Constr
+    {
+    public:
+      typedef LpSolverBase::Expr Expr;
+      typedef Expr::Key Key;
+      typedef Expr::Value Value;
+
+    protected:
+      Expr _expr;
+      Value _lb,_ub;
+    public:
+      ///\e
+      Constr() : _expr(), _lb(NaN), _ub(NaN) {}
+      ///\e
+      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) {}
+      ///\e
+      void clear()
+      {
+        _expr.clear();
+        _lb=_ub=NaN;
+      }
+
+      ///Reference to the linear expression
+      Expr &expr() { return _expr; }
+      ///Cont reference to the linear expression
+      const Expr &expr() const { return _expr; }
+      ///Reference to the lower bound.
+
+      ///\return
+      ///- \ref INF "INF": the constraint is lower unbounded.
+      ///- \ref NaN "NaN": lower bound has not been set.
+      ///- finite number: the lower bound
+      Value &lowerBound() { return _lb; }
+      ///The const version of \ref lowerBound()
+      const Value &lowerBound() const { return _lb; }
+      ///Reference to the upper bound.
+
+      ///\return
+      ///- \ref INF "INF": the constraint is upper unbounded.
+      ///- \ref NaN "NaN": upper bound has not been set.
+      ///- finite number: the upper bound
+      Value &upperBound() { return _ub; }
+      ///The const version of \ref upperBound()
+      const Value &upperBound() const { return _ub; }
+      ///Is the constraint lower bounded?
+      bool lowerBounded() const {
+        return isFinite(_lb);
+      }
+      ///Is the constraint upper bounded?
+      bool upperBounded() const {
+        return isFinite(_ub);
+      }
+
+    };
+
+    ///Linear expression of rows
+
+    ///This data structure represents a column of the matrix,
+    ///thas is it strores a linear expression of the dual variables
+    ///(\ref Row "Row"s).
+    ///
+    ///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;
+    ///e[v]+=12;
+    ///e.erase(v);
+    ///\endcode
+    ///or you can also iterate through its elements.
+    ///\code
+    ///double s=0;
+    ///for(LpSolverBase::DualExpr::iterator i=e.begin();i!=e.end();++i)
+    ///  s+=i->second;
+    ///\endcode
+    ///(This code computes the sum of all coefficients).
+    ///- Numbers (<tt>double</tt>'s)
+    ///and variables (\ref Row "Row"s) directly convert to an
+    ///\ref DualExpr and the usual linear operations are defined, so
+    ///\code
+    ///v+w
+    ///2*v-3.12*(v-w/2)
+    ///v*2.1+(3*v+(v*12+w)*3)/2
+    ///\endcode
+    ///are valid \ref DualExpr "DualExpr"essions.
+    ///The usual assignment operations are also defined.
+    ///\code
+    ///e=v+w;
+    ///e+=2*v-3.12*(v-w/2);
+    ///e*=3.4;
+    ///e/=5;
+    ///\endcode
+    ///
+    ///\sa Expr
+    ///
+    class DualExpr : public std::map<Row,Value>
+    {
+    public:
+      typedef LpSolverBase::Row Key;
+      typedef LpSolverBase::Value Value;
+
+    protected:
+      typedef std::map<Row,Value> Base;
+
+    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;
+        }
+      }
+
+      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;
+        }
+      }
+
+      ///Sets all coefficients to 0.
+      void clear() {
+        Base::clear();
+      }
+
+      ///\e
+      DualExpr &operator+=(const DualExpr &e) {
+        for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
+          (*this)[j->first]+=j->second;
+        return *this;
+      }
+      ///\e
+      DualExpr &operator-=(const DualExpr &e) {
+        for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
+          (*this)[j->first]-=j->second;
+        return *this;
+      }
+      ///\e
+      DualExpr &operator*=(const Value &c) {
+        for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
+          j->second*=c;
+        return *this;
+      }
+      ///\e
+      DualExpr &operator/=(const Value &c) {
+        for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
+          j->second/=c;
+        return *this;
+      }
+    };
+
+
+  private:
+
+    template <typename _Expr>
+    class MappedOutputIterator {
+    public:
+
+      typedef std::insert_iterator<_Expr> Base;
+
+      typedef std::output_iterator_tag iterator_category;
+      typedef void difference_type;
+      typedef void value_type;
+      typedef void reference;
+      typedef void pointer;
+
+      MappedOutputIterator(const Base& _base, const LpSolverBase& _lp)
+        : base(_base), lp(_lp) {}
+
+      MappedOutputIterator& operator*() {
+        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;
+      }
+
+    private:
+      Base base;
+      const LpSolverBase& lp;
+    };
+
+    template <typename Expr>
+    class MappedInputIterator {
+    public:
+
+      typedef typename Expr::const_iterator Base;
+
+      typedef typename Base::iterator_category iterator_category;
+      typedef typename Base::difference_type difference_type;
+      typedef const std::pair<int, Value> value_type;
+      typedef value_type reference;
+      class pointer {
+      public:
+        pointer(value_type& _value) : value(_value) {}
+        value_type* operator->() { return &value; }
+      private:
+        value_type value;
+      };
+
+      MappedInputIterator(const Base& _base, const LpSolverBase& _lp)
+        : base(_base), lp(_lp) {}
+
+      reference operator*() {
+        return std::make_pair(lp._lpId(base->first), base->second);
+      }
+
+      pointer operator->() {
+        return pointer(operator*());
+      }
+
+      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;
+    };
+
+  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 int _addCol() = 0;
+    virtual int _addRow() = 0;
+
+    virtual void _eraseCol(int col) = 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 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 _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 _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;
+
+    //Own protected stuff
+
+    //Constant component of the objective function
+    Value obj_const_comp;
+
+  public:
+
+    ///\e
+    LpSolverBase() : obj_const_comp(0) {}
+
+    ///\e
+    virtual ~LpSolverBase() {}
+
+    ///Creates a new LP problem
+    LpSolverBase* newLp() {return _newLp();}
+    ///Makes a copy of the LP problem
+    LpSolverBase* copyLp() {return _copyLp();}
+
+    ///\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)
+    ///\param t can be
+    ///- a standard STL compatible iterable container with
+    ///\ref Col as its \c values_type
+    ///like
+    ///\code
+    ///std::vector<LpSolverBase::Col>
+    ///std::list<LpSolverBase::Col>
+    ///\endcode
+    ///- a standard STL compatible iterable container with
+    ///\ref Col as its \c mapped_type
+    ///like
+    ///\code
+    ///std::map<AnyType,LpSolverBase::Col>
+    ///\endcode
+    ///- an iterable lemon \ref concepts::WriteMap "write map" like
+    ///\code
+    ///ListGraph::NodeMap<LpSolverBase::Col>
+    ///ListGraph::EdgeMap<LpSolverBase::Col>
+    ///\endcode
+    ///\return The number of the created column.
+#ifdef DOXYGEN
+    template<class T>
+    int addColSet(T &t) { return 0;}
+#else
+    template<class T>
+    typename enable_if<typename T::value_type::LpSolverCol,int>::type
+    addColSet(T &t,dummy<0> = 0) {
+      int s=0;
+      for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addCol();s++;}
+      return s;
+    }
+    template<class T>
+    typename enable_if<typename T::value_type::second_type::LpSolverCol,
+                       int>::type
+    addColSet(T &t,dummy<1> = 1) {
+      int s=0;
+      for(typename T::iterator i=t.begin();i!=t.end();++i) {
+        i->second=addCol();
+        s++;
+      }
+      return s;
+    }
+    template<class T>
+    typename enable_if<typename T::MapIt::Value::LpSolverCol,
+                       int>::type
+    addColSet(T &t,dummy<2> = 2) {
+      int s=0;
+      for(typename T::MapIt i(t); i!=INVALID; ++i)
+        {
+          i.set(addCol());
+          s++;
+        }
+      return s;
+    }
+#endif
+
+    ///Set a column (i.e a dual constraint) of the LP
+
+    ///\param c is the column to be modified
+    ///\param e is a dual linear expression (see \ref DualExpr)
+    ///a better one.
+    void col(Col c,const DualExpr &e) {
+      e.simplify();
+      _setColCoeffs(_lpId(c), ConstColIterator(e.begin(), *this),
+                    ConstColIterator(e.end(), *this));
+    }
+
+    ///Get a column (i.e a dual constraint) of the LP
+
+    ///\param r 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));
+      return e;
+    }
+
+    ///Add a new column to the LP
+
+    ///\param e is a dual linear expression (see \ref DualExpr)
+    ///\param obj is the corresponding component of the objective
+    ///function. It is 0 by default.
+    ///\return The created column.
+    Col addCol(const DualExpr &e, Value o = 0) {
+      Col c=addCol();
+      col(c,e);
+      objCoeff(c,o);
+      return c;
+    }
+
+    ///Add a new empty row (i.e a new constraint) to the LP
+
+    ///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)
+    ///\param t can be
+    ///- a standard STL compatible iterable container with
+    ///\ref Row as its \c values_type
+    ///like
+    ///\code
+    ///std::vector<LpSolverBase::Row>
+    ///std::list<LpSolverBase::Row>
+    ///\endcode
+    ///- a standard STL compatible iterable container with
+    ///\ref Row as its \c mapped_type
+    ///like
+    ///\code
+    ///std::map<AnyType,LpSolverBase::Row>
+    ///\endcode
+    ///- an iterable lemon \ref concepts::WriteMap "write map" like
+    ///\code
+    ///ListGraph::NodeMap<LpSolverBase::Row>
+    ///ListGraph::EdgeMap<LpSolverBase::Row>
+    ///\endcode
+    ///\return The number of rows created.
+#ifdef DOXYGEN
+    template<class T>
+    int addRowSet(T &t) { return 0;}
+#else
+    template<class T>
+    typename enable_if<typename T::value_type::LpSolverRow,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++;}
+      return s;
+    }
+    template<class T>
+    typename enable_if<typename T::value_type::second_type::LpSolverRow,
+                       int>::type
+    addRowSet(T &t,dummy<1> = 1) {
+      int s=0;
+      for(typename T::iterator i=t.begin();i!=t.end();++i) {
+        i->second=addRow();
+        s++;
+      }
+      return s;
+    }
+    template<class T>
+    typename enable_if<typename T::MapIt::Value::LpSolverRow,
+                       int>::type
+    addRowSet(T &t,dummy<2> = 2) {
+      int s=0;
+      for(typename T::MapIt i(t); i!=INVALID; ++i)
+        {
+          i.set(addRow());
+          s++;
+        }
+      return s;
+    }
+#endif
+
+    ///Set a row (i.e a constraint) of the LP
+
+    ///\param r is the row to be modified
+    ///\param l is lower bound (-\ref INF means no bound)
+    ///\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());
+    }
+
+    ///Set a row (i.e a constraint) of the LP
+
+    ///\param r is the row to be modified
+    ///\param c is a linear expression (see \ref Constr)
+    void row(Row r, const Constr &c) {
+      row(r, c.lowerBounded()?c.lowerBound():-INF,
+          c.expr(), c.upperBounded()?c.upperBound():INF);
+    }
+
+
+    ///Get a row (i.e a constraint) of the LP
+
+    ///\param r is the row to get
+    ///\return the expression associated to the row
+    Expr row(Row r) const {
+      Expr e;
+      _getRowCoeffs(_lpId(r), RowIterator(std::inserter(e, e.end()), *this));
+      return e;
+    }
+
+    ///Add a new row (i.e a new constraint) to the LP
+
+    ///\param l is the lower bound (-\ref INF means no bound)
+    ///\param e is a linear expression (see \ref Expr)
+    ///\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();
+      row(r,l,e,u);
+      return r;
+    }
+
+    ///Add a new row (i.e a new constraint) to the LP
+
+    ///\param c is a linear expression (see \ref Constr)
+    ///\return The created row.
+    Row addRow(const Constr &c) {
+      Row r=addRow();
+      row(r,c);
+      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
+
+    ///\param r is the row to be deleted
+    ///\todo Please check this
+    void eraseRow(Row r) {
+      _eraseRow(_lpId(r));
+      rows.eraseId(r.id);
+    }
+
+    /// Get the name of a column
+
+    ///\param c is the coresponding coloumn
+    ///\return The name of the colunm
+    std::string colName(Col c) const {
+      std::string name;
+      _getColName(_lpId(c), name);
+      return name;
+    }
+
+    /// Set the name of a column
+
+    ///\param c is the coresponding coloumn
+    ///\param name The name to be given
+    void colName(Col c, const std::string& name) {
+      _setColName(_lpId(c), name);
+    }
+
+    /// Get the column by its name
+
+    ///\param name The name of the column
+    ///\return the proper column or \c INVALID
+    Col colByName(const std::string& name) const {
+      int k = _colByName(name);
+      return k != -1 ? Col(cols.fixId(k)) : Col(INVALID);
+    }
+
+    /// Set an element of the coefficient matrix of the LP
+
+    ///\param r is the row of the element to be modified
+    ///\param c is the coloumn 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);
+    }
+
+    /// 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
+    ///\return the corresponding coefficient
+
+    Value coeff(Row r, Col c) const {
+      return _getCoeff(_lpId(r),_lpId(c));
+    }
+
+    /// Set the lower bound of a column (i.e a variable)
+
+    /// 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.
+    void colLowerBound(Col c, Value value) {
+      _setColLowerBound(_lpId(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 might be -\ref INF as well).
+    ///\return The lower bound for coloumn \t c
+    Value colLowerBound(Col c) const {
+      return _getColLowerBound(_lpId(c));
+    }
+
+    ///\brief Set the lower bound of  several columns
+    ///(i.e a 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.
+#ifdef DOXYGEN
+    template<class T>
+    void colLowerBound(T &t, Value value) { return 0;}
+#else
+    template<class T>
+    typename enable_if<typename T::value_type::LpSolverCol,void>::type
+    colLowerBound(T &t, Value value,dummy<0> = 0) {
+      for(typename T::iterator i=t.begin();i!=t.end();++i) {
+        colLowerBound(*i, value);
+      }
+    }
+    template<class T>
+    typename enable_if<typename T::value_type::second_type::LpSolverCol,
+                       void>::type
+    colLowerBound(T &t, Value value,dummy<1> = 1) {
+      for(typename T::iterator i=t.begin();i!=t.end();++i) {
+        colLowerBound(i->second, value);
+      }
+    }
+    template<class T>
+    typename enable_if<typename T::MapIt::Value::LpSolverCol,
+                       void>::type
+    colLowerBound(T &t, Value value,dummy<2> = 2) {
+      for(typename T::MapIt i(t); i!=INVALID; ++i){
+        colLowerBound(*i, value);
+      }
+    }
+#endif
+
+    /// Set the upper bound of a column (i.e a variable)
+
+    /// 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.
+    void colUpperBound(Col c, Value value) {
+      _setColUpperBound(_lpId(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 might be \ref INF as well).
+    ///\return The upper bound for coloumn \t c
+    Value colUpperBound(Col c) const {
+      return _getColUpperBound(_lpId(c));
+    }
+
+    ///\brief Set the upper bound of  several columns
+    ///(i.e a 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.
+#ifdef DOXYGEN
+    template<class T>
+    void colUpperBound(T &t, Value value) { return 0;}
+#else
+    template<class T>
+    typename enable_if<typename T::value_type::LpSolverCol,void>::type
+    colUpperBound(T &t, Value value,dummy<0> = 0) {
+      for(typename T::iterator i=t.begin();i!=t.end();++i) {
+        colUpperBound(*i, value);
+      }
+    }
+    template<class T>
+    typename enable_if<typename T::value_type::second_type::LpSolverCol,
+                       void>::type
+    colUpperBound(T &t, Value value,dummy<1> = 1) {
+      for(typename T::iterator i=t.begin();i!=t.end();++i) {
+        colUpperBound(i->second, value);
+      }
+    }
+    template<class T>
+    typename enable_if<typename T::MapIt::Value::LpSolverCol,
+                       void>::type
+    colUpperBound(T &t, Value value,dummy<2> = 2) {
+      for(typename T::MapIt i(t); i!=INVALID; ++i){
+        colUpperBound(*i, value);
+      }
+    }
+#endif
+
+    /// Set the lower and the upper bounds of a column (i.e a variable)
+
+    /// The lower and the upper bounds of
+    /// a variable (column) have to be given by an
+    /// extended number of type Value, i.e. a finite number of type
+    /// Value, -\ref INF or \ref INF.
+    void colBounds(Col c, Value lower, Value upper) {
+      _setColLowerBound(_lpId(c),lower);
+      _setColUpperBound(_lpId(c),upper);
+    }
+
+    ///\brief Set the lower and the upper bound of several columns
+    ///(i.e a variables) at once
+    ///
+    ///This magic function takes a container as its argument
+    ///and applies the function on all of its elements.
+    /// The lower and the upper bounds of
+    /// a variable (column) have to be given by an
+    /// extended number of type Value, i.e. a finite number of type
+    /// Value, -\ref INF or \ref INF.
+#ifdef DOXYGEN
+    template<class T>
+    void colBounds(T &t, Value lower, Value upper) { return 0;}
+#else
+    template<class T>
+    typename enable_if<typename T::value_type::LpSolverCol,void>::type
+    colBounds(T &t, Value lower, Value upper,dummy<0> = 0) {
+      for(typename T::iterator i=t.begin();i!=t.end();++i) {
+        colBounds(*i, lower, upper);
+      }
+    }
+    template<class T>
+    typename enable_if<typename T::value_type::second_type::LpSolverCol,
+                       void>::type
+    colBounds(T &t, Value lower, Value upper,dummy<1> = 1) {
+      for(typename T::iterator i=t.begin();i!=t.end();++i) {
+        colBounds(i->second, lower, upper);
+      }
+    }
+    template<class T>
+    typename enable_if<typename T::MapIt::Value::LpSolverCol,
+                       void>::type
+    colBounds(T &t, Value lower, Value upper,dummy<2> = 2) {
+      for(typename T::MapIt i(t); i!=INVALID; ++i){
+        colBounds(*i, lower, upper);
+      }
+    }
+#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 an element of the objective function
+    void objCoeff(Col c, Value v) {_setObjCoeff(_lpId(c),v); };
+
+    ///Get an element of the objective function
+    Value objCoeff(Col c) const { return _getObjCoeff(_lpId(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();
+    }
+
+    ///Get the objective function
+
+    ///\return the objective function as a linear expression of type \ref 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));
+        }
+      }
+      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(); }
+
+    ///@}
+
+
+    ///\name Solve the LP
+
+    ///@{
+
+    ///\e Solve the LP problem at hand
+    ///
+    ///\return The result of the optimization procedure. Possible
+    ///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(); }
+
+    ///@}
+
+    ///\name Obtain the solution
+
+    ///@{
+
+    /// 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
+    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;
+      }
+      return res;
+    }
+
+    ///\e
+    Value dual(Row r) const { return _getDual(_lpId(r)); }
+    ///\e
+    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;
+      }
+      return res;
+    }
+
+    ///\e
+    bool isBasicCol(Col c) const { return _isBasicCol(_lpId(c)); }
+
+    ///\e
+
+    ///\return
+    ///- \ref INF or -\ref INF means either infeasibility or unboundedness
+    /// of the primal 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 primalValue() const { return _getPrimalValue()+obj_const_comp;}
+    ///@}
+
+  };
+
+
+  /// \ingroup lp_group
+  ///
+  /// \brief Common base class for MIP solvers
+  /// \todo Much more docs
+  class MipSolverBase : virtual public LpSolverBase{
+  public:
+
+    ///Possible variable (coloumn) types (e.g. real, integer, binary etc.)
+    enum ColTypes {
+      ///Continuous variable
+      REAL = 0,
+      ///Integer variable
+
+      ///Unfortunately, cplex 7.5 somewhere writes something like
+      ///#define INTEGER 'I'
+      INT = 1
+      ///\todo No support for other types yet.
+    };
+
+    ///Sets the type of the given coloumn to the given type
+    ///
+    ///Sets the type of the given coloumn to the given type.
+    void colType(Col c, ColTypes col_type) {
+      _colType(_lpId(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();
+    }
+
+  protected:
+
+    virtual ColTypes _colType(int col) const = 0;
+    virtual void _colType(int col, ColTypes col_type) = 0;
+    virtual SolutionStatus _getMipStatus() 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;
+  }
+
+
+} //namespace lemon
+
+#endif //LEMON_LP_BASE_H
diff -r 69928a704ffb -r 7afc121e0689 lemon/lp_cplex.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/lp_cplex.cc	Tue Dec 02 21:40:33 2008 +0100
@@ -0,0 +1,699 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#include <iostream>
+#include <vector>
+#include <lemon/lp_cplex.h>
+
+extern "C" {
+#include <ilcplex/cplex.h>
+}
+
+
+///\file
+///\brief Implementation of the LEMON-CPLEX lp solver interface.
+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);
+    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);
+    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);
+    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) {
+      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
+  {
+    int index;
+    if (CPXgetcolindex(env, lp,
+                       const_cast<char*>(name.c_str()), &index) == 0) {
+      return index;
+    }
+    return -1;
+  }
+
+  ///\warning Data at index 0 is ignored in the arrays.
+  void LpCplex::_setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e)
+  {
+    std::vector<int> indices;
+    std::vector<int> rowlist;
+    std::vector<Value> values;
+
+    for(ConstRowIterator it=b; it!=e; ++it) {
+      indices.push_back(it->first);
+      values.push_back(it->second);
+      rowlist.push_back(i);
+    }
+
+    status = CPXchgcoeflist(env, lp, values.size(),
+                            &rowlist[0], &indices[0], &values[0]);
+  }
+
+  void LpCplex::_getRowCoeffs(int i, RowIterator b) const {
+    int tmp1, tmp2, tmp3, length;
+    CPXgetrows(env, lp, &tmp1, &tmp2, 0, 0, 0, &length, i, i);
+
+    length = -length;
+    std::vector<int> indices(length);
+    std::vector<double> values(length);
+
+    CPXgetrows(env, lp, &tmp1, &tmp2, &indices[0], &values[0],
+               length, &tmp3, i, i);
+
+    for (int i = 0; i < length; ++i) {
+      *b = std::make_pair(indices[i], values[i]);
+      ++b;
+    }
+
+    /// \todo implement
+  }
+
+  void LpCplex::_setColCoeffs(int i, ConstColIterator b, ConstColIterator e)
+  {
+    std::vector<int> indices;
+    std::vector<int> collist;
+    std::vector<Value> values;
+
+    for(ConstColIterator it=b; it!=e; ++it) {
+      indices.push_back(it->first);
+      values.push_back(it->second);
+      collist.push_back(i);
+    }
+
+    status = CPXchgcoeflist(env, lp, values.size(),
+                            &indices[0], &collist[0], &values[0]);
+  }
+
+  void LpCplex::_getColCoeffs(int i, ColIterator b) const {
+
+    int tmp1, tmp2, tmp3, length;
+    CPXgetcols(env, lp, &tmp1, &tmp2, 0, 0, 0, &length, i, i);
+
+    length = -length;
+    std::vector<int> indices(length);
+    std::vector<double> values(length);
+
+    CPXgetcols(env, lp, &tmp1, &tmp2, &indices[0], &values[0],
+               length, &tmp3, i, i);
+
+    for (int i = 0; i < length; ++i) {
+      *b = std::make_pair(indices[i], values[i]);
+      ++b;
+    }
+
+  }
+
+  void LpCplex::_setCoeff(int row, int col, Value value)
+  {
+    CPXchgcoef(env, lp, row, col, value);
+  }
+
+  LpCplex::Value LpCplex::_getCoeff(int row, int col) const
+  {
+    LpCplex::Value value;
+    CPXgetcoef(env, lp, row, col, &value);
+    return value;
+  }
+
+  void LpCplex::_setColLowerBound(int i, Value value)
+  {
+    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
+  {
+    LpCplex::Value x;
+    CPXgetlb (env, lp, &x, i, i);
+    if (x <= -CPX_INFBOUND) x = -INF;
+    return x;
+  }
+
+  void LpCplex::_setColUpperBound(int i, Value value)
+  {
+    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
+  {
+    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
+  {
+    Value x;
+    CPXgetcoef(env, lp, -1, i, &x);
+    return x;
+  }
+
+  void LpCplex::_clearObj()
+  {
+    for (int i=0;i< CPXgetnumcols(env, lp);++i){
+      CPXchgcoef(env, lp, -1, i, 0);
+    }
+
+  }
+  // The routine returns zero unless an error occurred during the
+  // optimization. Examples of errors include exhausting available
+  // memory (CPXERR_NO_MEMORY) or encountering invalid data in the
+  // CPLEX problem object (CPXERR_NO_PROBLEM). Exceeding a
+  // user-specified CPLEX limit, or proving the model infeasible or
+  // unbounded, are not considered errors. Note that a zero return
+  // value does not necessarily mean that a solution exists. Use query
+  // 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);
+#if CPX_VERSION >= 800
+    if (status)
+    {
+      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){
+      //We want to exclude some cases
+      switch (CPXgetstat(env, lp)){
+      case CPX_OBJ_LIM:
+      case CPX_IT_LIM_FEAS:
+      case CPX_IT_LIM_INFEAS:
+      case CPX_TIME_LIM_FEAS:
+      case CPX_TIME_LIM_INFEAS:
+        return UNSOLVED;
+      default:
+        return SOLVED;
+      }
+    }
+    else{
+      return UNSOLVED;
+    }
+#endif
+  }
+
+  LpCplex::Value LpCplex::_getPrimal(int i) const
+  {
+    Value x;
+    CPXgetx(env, lp, &x, i, i);
+    return x;
+  }
+
+  LpCplex::Value LpCplex::_getDual(int i) const
+  {
+    Value y;
+    CPXgetpi(env, lp, &y, i, i);
+    return y;
+  }
+
+  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);
+    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
+#if CPX_VERSION < 900
+  void statusSwitch(CPXENVptr env,int& stat){
+    int lpmethod;
+    CPXgetintparam (env,CPX_PARAM_LPMETHOD,&lpmethod);
+    if (lpmethod==2){
+      if (stat==CPX_UNBOUNDED){
+        stat=CPX_INFEASIBLE;
+      }
+      else{
+        if (stat==CPX_INFEASIBLE)
+          stat=CPX_UNBOUNDED;
+      }
+    }
+  }
+#else
+  void statusSwitch(CPXENVptr,int&){}
+#endif
+
+  LpCplex::SolutionStatus LpCplex::_getPrimalStatus() const
+  {
+    //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.
+
+    int stat = CPXgetstat(env, lp);
+#if CPX_VERSION >= 800
+    switch (stat)
+    {
+      case CPX_STAT_OPTIMAL:
+        return OPTIMAL;
+      case CPX_STAT_UNBOUNDED:
+        return INFINITE;
+      case CPX_STAT_INFEASIBLE:
+        return INFEASIBLE;
+      default:
+        return UNDEFINED;
+    }
+#else
+    statusSwitch(env,stat);
+    //CPXgetstat(env, lp);
+    //printf("A primal status: %d, CPX_OPTIMAL=%d \n",stat,CPX_OPTIMAL);
+    switch (stat) {
+    case 0:
+      return UNDEFINED; //Undefined
+    case CPX_OPTIMAL://Optimal
+      return OPTIMAL;
+    case CPX_UNBOUNDED://Unbounded
+      return INFEASIBLE;//In case of dual simplex
+      //return INFINITE;
+    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
+      //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;
+    default:
+      return UNDEFINED; //Everything else comes here
+      //FIXME error
+    }
+#endif
+  }
+
+//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);
+#if CPX_VERSION >= 800
+    switch (stat)
+    {
+      case CPX_STAT_OPTIMAL:
+        return OPTIMAL;
+      case CPX_STAT_UNBOUNDED:
+        return INFEASIBLE;
+      default:
+        return UNDEFINED;
+    }
+#else
+    statusSwitch(env,stat);
+    switch (stat) {
+    case 0:
+      return UNDEFINED; //Undefined
+    case CPX_OPTIMAL://Optimal
+      return OPTIMAL;
+    case CPX_UNBOUNDED:
+     return INFEASIBLE;
+    default:
+      return UNDEFINED; //Everything else comes here
+      //FIXME error
+    }
+#endif
+  }
+
+  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;
+    }
+#else
+    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;
+  }
+
+} //namespace lemon
+
diff -r 69928a704ffb -r 7afc121e0689 lemon/lp_cplex.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/lp_cplex.h	Tue Dec 02 21:40:33 2008 +0100
@@ -0,0 +1,113 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_LP_CPLEX_H
+#define LEMON_LP_CPLEX_H
+
+///\file
+///\brief Header of the LEMON-CPLEX lp solver interface.
+
+#include <lemon/lp_base.h>
+
+struct cpxenv;
+struct cpxlp;
+
+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();
+
+
+    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 _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;
+    virtual void _setObjCoeff(int i, Value obj_coef);
+    virtual Value _getObjCoeff(int i) const;
+    virtual void _clearObj();
+
+
+    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;
+
+    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; }
+
+  };
+} //END OF NAMESPACE LEMON
+
+#endif //LEMON_LP_CPLEX_H
+
diff -r 69928a704ffb -r 7afc121e0689 lemon/lp_glpk.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/lp_glpk.cc	Tue Dec 02 21:40:33 2008 +0100
@@ -0,0 +1,644 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+///\file
+///\brief Implementation of the LEMON-GLPK lp 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
+
+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;
+    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;
+    return i;
+  }
+
+
+  void LpGlpk::_eraseCol(int i) {
+    int ca[2];
+    ca[1]=i;
+    LEMON_glp(del_cols)(lp, 1, ca);
+    solved = false;
+  }
+
+  void LpGlpk::_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()));
+    return k > 0 ? k : -1;
+  }
+
+
+  void LpGlpk::_setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e)
+  {
+    std::vector<int> indices;
+    std::vector<Value> values;
+
+    indices.push_back(0);
+    values.push_back(0);
+
+    for(ConstRowIterator it=b; it!=e; ++it) {
+      indices.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);
+    std::vector<Value> values(length + 1);
+
+    LEMON_glp(get_mat_row)(lp, ix, &indices[0], &values[0]);
+
+    for (int i = 1; i <= length; ++i) {
+      *b = std::make_pair(indices[i], values[i]);
+      ++b;
+    }
+  }
+
+  void LpGlpk::_setColCoeffs(int ix, ConstColIterator b, ConstColIterator e) {
+
+    std::vector<int> indices;
+    std::vector<Value> values;
+
+    indices.push_back(0);
+    values.push_back(0);
+
+    for(ConstColIterator it=b; it!=e; ++it) {
+      indices.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);
+    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]);
+      ++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);
+      std::vector<Value> values(length + 2);
+
+      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
+      bool found=false;
+      for (int i = 1; i <= length; ++i) {
+        if (indices[i]==jx){
+          found=true;
+          values[i]=value;
+          break;
+        }
+      }
+      if (!found){
+        ++length;
+        indices[length]=jx;
+        values[length]=value;
+      }
+
+      LEMON_glp(set_mat_row)(lp, ix, length, &indices[0], &values[0]);
+
+    } else {
+
+      int length=LEMON_glp(get_mat_col)(lp, jx, 0, 0);
+
+      std::vector<int> indices(length + 2);
+      std::vector<Value> values(length + 2);
+
+      LEMON_glp(get_mat_col)(lp, jx, &indices[0], &values[0]);
+
+      //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]==ix){
+          found=true;
+          values[i]=value;
+          break;
+        }
+      }
+      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);
+    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){
+        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) {
+      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
+      }
+    } 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);
+        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);
+      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
+      }
+    } 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);
+        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);
+      switch (b) {
+      case LEMON_GLP(UP):
+      case LEMON_GLP(DB):
+      case LEMON_GLP(FX):
+        return LEMON_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);
+        }
+        else{
+          LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(DB), lb, ub);
+        }
+      }
+    }
+
+    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;
+    default:
+      lb=LEMON_glp(get_row_lb)(lp, i);
+    }
+
+    switch (b) {
+    case LEMON_GLP(FR):
+    case LEMON_GLP(LO):
+      ub = INF;
+        break;
+    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;
+    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 OPTIMAL;
+      default:
+        return FEASIBLE;
+      }
+    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);
+  }
+
+
+} //END OF NAMESPACE LEMON
diff -r 69928a704ffb -r 7afc121e0689 lemon/lp_glpk.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/lp_glpk.h	Tue Dec 02 21:40:33 2008 +0100
@@ -0,0 +1,139 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_LP_GLPK_H
+#define LEMON_LP_GLPK_H
+
+///\file
+///\brief Header of the LEMON-GLPK lp solver interface.
+///\ingroup lp_group
+
+#include <lemon/lp_base.h>
+
+// forward declaration
+#ifndef _GLP_PROB
+#define _GLP_PROB
+typedef struct { double _prob; } glp_prob;
+/* LP/MIP problem object */
+#endif
+
+namespace lemon {
+
+
+  /// \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;
+
+    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
+    ///
+    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
+
+    ///\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);
+    ///Turns on or off the presolver
+
+    ///Turns on (\c b is \c true) or off (\c b is \c false) the presolver
+    ///
+    ///The presolver is off by default.
+    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); }
+  };
+} //END OF NAMESPACE LEMON
+
+#endif //LEMON_LP_GLPK_H
+
diff -r 69928a704ffb -r 7afc121e0689 lemon/lp_skeleton.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/lp_skeleton.cc	Tue Dec 02 21:40:33 2008 +0100
@@ -0,0 +1,187 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#include <lemon/lp_skeleton.h>
+
+///\file
+///\brief A skeleton file to implement LP solver interfaces
+namespace lemon {
+
+  LpSolverBase* LpSkeleton::_newLp()
+  {
+    LpSolverBase *tmp=0;
+    return tmp;
+  }
+
+  LpSolverBase* LpSkeleton::_copyLp()
+  {
+    LpSolverBase *tmp=0;
+    return tmp;
+  }
+
+  int LpSkeleton::_addCol()
+  {
+    return ++col_num;
+  }
+
+  int LpSkeleton::_addRow()
+  {
+    return ++row_num;
+  }
+
+  void LpSkeleton::_eraseCol(int ) {
+  }
+
+  void LpSkeleton::_eraseRow(int) {
+  }
+
+  void LpSkeleton::_getColName(int, std::string &) const {
+  }
+
+
+  void LpSkeleton::_setColName(int, const std::string &) {
+  }
+
+  int LpSkeleton::_colByName(const std::string&) const { return -1; }
+
+
+  void LpSkeleton::_setRowCoeffs(int, ConstRowIterator, ConstRowIterator) {
+  }
+
+  void LpSkeleton::_getRowCoeffs(int, RowIterator) const {
+  }
+
+  void LpSkeleton::_setColCoeffs(int, ConstColIterator, ConstColIterator) {
+  }
+
+  void LpSkeleton::_getColCoeffs(int, ColIterator) const {
+  }
+
+  void LpSkeleton::_setCoeff(int, int, Value )
+  {
+  }
+
+  LpSkeleton::Value LpSkeleton::_getCoeff(int, int) const
+  {
+    return 0;
+  }
+
+
+  void LpSkeleton::_setColLowerBound(int, Value)
+  {
+  }
+
+  LpSkeleton::Value LpSkeleton::_getColLowerBound(int) const
+  {
+    return 0;
+  }
+
+  void LpSkeleton::_setColUpperBound(int, Value)
+  {
+  }
+
+  LpSkeleton::Value LpSkeleton::_getColUpperBound(int) const
+  {
+    return 0;
+  }
+
+//   void LpSkeleton::_setRowLowerBound(int, Value)
+//   {
+//   }
+
+//   void LpSkeleton::_setRowUpperBound(int, Value)
+//   {
+//   }
+
+  void LpSkeleton::_setRowBounds(int, Value, Value)
+  {
+  }
+
+  void LpSkeleton::_getRowBounds(int, Value&, Value&) const
+  {
+  }
+
+  void LpSkeleton::_setObjCoeff(int, Value)
+  {
+  }
+
+  LpSkeleton::Value LpSkeleton::_getObjCoeff(int) const
+  {
+    return 0;
+  }
+
+  void LpSkeleton::_setMax()
+  {
+  }
+
+  void LpSkeleton::_setMin()
+  {
+  }
+
+  bool LpSkeleton::_isMax() const
+  {
+    return true;
+  }
+
+
+  void LpSkeleton::_clearObj()
+  {
+  }
+
+  LpSkeleton::SolveExitStatus LpSkeleton::_solve()
+  {
+    return SOLVED;
+  }
+
+  LpSkeleton::Value LpSkeleton::_getPrimal(int) const
+  {
+    return 0;
+  }
+
+  LpSkeleton::Value LpSkeleton::_getDual(int) const
+  {
+    return 0;
+  }
+
+  LpSkeleton::Value LpSkeleton::_getPrimalValue() const
+  {
+    return 0;
+  }
+
+  LpSkeleton::SolutionStatus LpSkeleton::_getPrimalStatus() const
+  {
+    return UNDEFINED;
+  }
+
+  LpSkeleton::SolutionStatus LpSkeleton::_getDualStatus() const
+  {
+    return UNDEFINED;
+  }
+
+  LpSkeleton::ProblemTypes LpSkeleton::_getProblemType() const
+  {
+    return UNKNOWN;
+  }
+
+  bool LpSkeleton::_isBasicCol(int) const
+  {
+    return true;
+  }
+
+} //namespace lemon
+
diff -r 69928a704ffb -r 7afc121e0689 lemon/lp_skeleton.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/lp_skeleton.h	Tue Dec 02 21:40:33 2008 +0100
@@ -0,0 +1,183 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_LP_SKELETON
+#define LEMON_LP_SKELETON
+
+#include <lemon/lp_base.h>
+
+///\file
+///\brief A skeleton file to implement LP solver interfaces
+namespace lemon {
+
+  ///A skeleton class to implement LP solver interfaces
+  class LpSkeleton :public LpSolverBase {
+    int col_num,row_num;
+
+  protected:
+
+    ///\e
+    virtual LpSolverBase* _newLp();
+    ///\e
+    virtual LpSolverBase* _copyLp();
+    /// \e
+    virtual int _addCol();
+    /// \e
+    virtual int _addRow();
+    /// \e
+    virtual void _eraseCol(int i);
+    /// \e
+    virtual void _eraseRow(int i);
+    /// \e
+    virtual void _getColName(int col, std::string & name) const;
+    /// \e
+    virtual void _setColName(int col, const std::string & name);
+    /// \e
+    virtual int _colByName(const std::string& name) const;
+
+    /// \e
+    virtual void _setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e);
+    /// \e
+    virtual void _getRowCoeffs(int i, RowIterator b) const;
+    /// \e
+    virtual void _setColCoeffs(int i, ConstColIterator b, ConstColIterator e);
+    /// \e
+    virtual void _getColCoeffs(int i, ColIterator b) const;
+
+    /// Set one element of the coefficient matrix
+    virtual void _setCoeff(int row, int col, Value value);
+
+    /// Get one element of the coefficient matrix
+    virtual Value _getCoeff(int row, int col) const;
+
+    /// The lower bound of a variable (column) have to be given by an
+    /// extended number of type Value, i.e. a finite number of type
+    /// Value or -\ref INF.
+    virtual void _setColLowerBound(int i, Value value);
+    /// \e
+
+    /// The lower bound of a variable (column) is an
+    /// extended number of type Value, i.e. a finite number of type
+    /// Value or -\ref INF.
+    virtual Value _getColLowerBound(int i) const;
+
+    /// The upper bound of a variable (column) have to be given by an
+    /// extended number of type Value, i.e. a finite number of type
+    /// Value or \ref INF.
+    virtual void _setColUpperBound(int i, Value value);
+    /// \e
+
+    /// The upper bound of a variable (column) is an
+    /// extended number of type Value, i.e. a finite number of type
+    /// Value or \ref INF.
+    virtual Value _getColUpperBound(int i) const;
+
+//     /// The lower bound of a linear expression (row) have to be given by an
+//     /// extended number of type Value, i.e. a finite number of type
+//     /// Value or -\ref INF.
+//     virtual void _setRowLowerBound(int i, Value value);
+//     /// \e
+
+//     /// The upper bound of a linear expression (row) have to be given by an
+//     /// extended number of type Value, i.e. a finite number of type
+//     /// Value or \ref INF.
+//     virtual void _setRowUpperBound(int i, Value value);
+
+    /// The lower and upper bound of a linear expression (row) have to be
+    /// given by an
+    /// extended number of type Value, i.e. a finite number of type
+    /// Value or +/-\ref INF.
+    virtual void _setRowBounds(int i, Value lb, Value ub);
+    /// \e
+
+
+    /// 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.
+    virtual void _getRowBounds(int i, Value &lb, Value &ub) const;
+    /// \e
+
+
+    /// \e
+    virtual void _clearObj();
+    /// \e
+    virtual void _setObjCoeff(int i, Value obj_coef);
+
+    /// \e
+    virtual Value _getObjCoeff(int i) const;
+
+    ///\e
+
+    ///\bug Wrong interface
+    ///
+    virtual SolveExitStatus _solve();
+
+    ///\e
+
+    ///\bug Wrong interface
+    ///
+    virtual Value _getPrimal(int i) const;
+
+    ///\e
+
+    ///\bug Wrong interface
+    ///
+    virtual Value _getDual(int i) const;
+
+    ///\e
+
+    ///\bug Wrong interface
+    ///
+    virtual Value _getPrimalValue() const;
+
+    ///\e
+
+    ///\bug Wrong interface
+    ///
+    virtual SolutionStatus _getPrimalStatus() const;
+
+    ////e
+    virtual SolutionStatus _getDualStatus() const;
+
+
+    ///\e
+    virtual ProblemTypes _getProblemType() const;
+
+    ///\e
+    virtual void _setMax();
+    ///\e
+    virtual void _setMin();
+
+    ///\e
+    virtual bool _isMax() const;
+
+
+
+    ///\e
+    virtual bool _isBasicCol(int i) const;
+
+
+
+  public:
+    LpSkeleton() : LpSolverBase(), col_num(0), row_num(0) {}
+  };
+
+} //namespace lemon
+
+#endif // LEMON_LP_SKELETON
diff -r 69928a704ffb -r 7afc121e0689 lemon/lp_soplex.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/lp_soplex.cc	Tue Dec 02 21:40:33 2008 +0100
@@ -0,0 +1,316 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#include<iostream>
+#include<lemon/lp_soplex.h>
+
+#include <soplex/soplex.h>
+
+
+///\file
+///\brief Implementation of the LEMON-SOPLEX lp solver interface.
+namespace lemon {
+
+  LpSoplex::LpSoplex() : LpSolverBase() {
+    rows.setIdHandler(relocateIdHandler);
+    cols.setIdHandler(relocateIdHandler);
+    soplex = new soplex::SoPlex;
+    solved = false;
+  }
+
+  LpSoplex::~LpSoplex() {
+    delete soplex;
+  }
+
+  LpSoplex::LpSoplex(const LpSoplex& lp) : LpSolverBase() {
+    rows = lp.rows;
+    rows.setIdHandler(relocateIdHandler);
+
+    cols = lp.cols;
+    cols.setIdHandler(relocateIdHandler);
+
+    soplex = new soplex::SoPlex;
+    (*static_cast<soplex::SPxLP*>(soplex)) = *(lp.soplex);
+
+    colNames = lp.colNames;
+    invColNames = lp.invColNames;
+
+    primal_value = lp.primal_value;
+    dual_value = lp.dual_value;
+
+  }
+
+  LpSolverBase* LpSoplex::_newLp() {
+    LpSoplex* newlp = new LpSoplex();
+    return newlp;
+  }
+
+  LpSolverBase* LpSoplex::_copyLp() {
+    LpSoplex* newlp = new LpSoplex(*this);
+    return newlp;
+  }
+
+  int LpSoplex::_addCol() {
+    soplex::LPCol c;
+    c.setLower(-soplex::infinity);
+    c.setUpper(soplex::infinity);
+    soplex->addCol(c);
+
+    colNames.push_back(std::string());
+    primal_value.push_back(0.0);
+    solved = false;
+
+    return soplex->nCols() - 1;
+  }
+
+  int LpSoplex::_addRow() {
+    soplex::LPRow r;
+    r.setLhs(-soplex::infinity);
+    r.setRhs(soplex::infinity);
+    soplex->addRow(r);
+
+    dual_value.push_back(0.0);
+    solved = false;
+
+    return soplex->nRows() - 1;
+  }
+
+
+  void LpSoplex::_eraseCol(int i) {
+    soplex->removeCol(i);
+    invColNames.erase(colNames[i]);
+    colNames[i] = colNames.back();
+    invColNames[colNames.back()] = i;
+    colNames.pop_back();
+    primal_value[i] = primal_value.back();
+    primal_value.pop_back();
+    solved = false;
+  }
+
+  void LpSoplex::_eraseRow(int i) {
+    soplex->removeRow(i);
+    dual_value[i] = dual_value.back();
+    dual_value.pop_back();
+    solved = false;
+  }
+
+  void LpSoplex::_getColName(int c, std::string &name) const {
+    name = colNames[c];
+  }
+
+  void LpSoplex::_setColName(int c, const std::string &name) {
+    invColNames.erase(colNames[c]);
+    colNames[c] = name;
+    if (!name.empty()) {
+      invColNames.insert(std::make_pair(name, c));
+    }
+  }
+
+  int LpSoplex::_colByName(const std::string& name) const {
+    std::map<std::string, int>::const_iterator it =
+      invColNames.find(name);
+    if (it != invColNames.end()) {
+      return it->second;
+    } else {
+      return -1;
+    }
+  }
+
+
+  void LpSoplex::_setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e) {
+    for (int j = 0; j < soplex->nCols(); ++j) {
+      soplex->changeElement(i, j, 0.0);
+    }
+    for(ConstRowIterator it = b; it != e; ++it) {
+      soplex->changeElement(i, it->first, it->second);
+    }
+    solved = false;
+  }
+
+  void LpSoplex::_getRowCoeffs(int i, RowIterator b) const {
+    const soplex::SVector& vec = soplex->rowVector(i);
+    for (int k = 0; k < vec.size(); ++k) {
+      *b = std::make_pair(vec.index(k), vec.value(k));
+      ++b;
+    }
+  }
+
+  void LpSoplex::_setColCoeffs(int j, ConstColIterator b, ConstColIterator e) {
+    for (int i = 0; i < soplex->nRows(); ++i) {
+      soplex->changeElement(i, j, 0.0);
+    }
+    for(ConstColIterator it = b; it != e; ++it) {
+      soplex->changeElement(it->first, j, it->second);
+    }
+    solved = false;
+  }
+
+  void LpSoplex::_getColCoeffs(int i, ColIterator b) const {
+    const soplex::SVector& vec = soplex->colVector(i);
+    for (int k = 0; k < vec.size(); ++k) {
+      *b = std::make_pair(vec.index(k), vec.value(k));
+      ++b;
+    }
+  }
+
+  void LpSoplex::_setCoeff(int i, int j, Value value) {
+    soplex->changeElement(i, j, value);
+    solved = false;
+  }
+
+  LpSoplex::Value LpSoplex::_getCoeff(int i, int j) const {
+    return soplex->rowVector(i)[j];
+  }
+
+  void LpSoplex::_setColLowerBound(int i, Value value) {
+    soplex->changeLower(i, value != -INF ? value : -soplex::infinity);
+    solved = false;
+  }
+
+  LpSoplex::Value LpSoplex::_getColLowerBound(int i) const {
+    double value = soplex->lower(i);
+    return value != -soplex::infinity ? value : -INF;
+  }
+
+  void LpSoplex::_setColUpperBound(int i, Value value) {
+    soplex->changeUpper(i, value != INF ? value : soplex::infinity);
+    solved = false;
+  }
+
+  LpSoplex::Value LpSoplex::_getColUpperBound(int i) const {
+    double value = soplex->upper(i);
+    return value != soplex::infinity ? value : INF;
+  }
+
+  void LpSoplex::_setRowBounds(int i, Value lb, Value ub) {
+    soplex->changeRange(i, lb != -INF ? lb : -soplex::infinity,
+                        ub != INF ? ub : soplex::infinity);
+    solved = false;
+  }
+  void LpSoplex::_getRowBounds(int i, Value &lower, Value &upper) const {
+    lower = soplex->lhs(i);
+    if (lower == -soplex::infinity) lower = -INF;
+    upper = soplex->rhs(i);
+    if (upper == -soplex::infinity) upper = INF;
+  }
+
+  void LpSoplex::_setObjCoeff(int i, Value obj_coef) {
+    soplex->changeObj(i, obj_coef);
+    solved = false;
+  }
+
+  LpSoplex::Value LpSoplex::_getObjCoeff(int i) const {
+    return soplex->obj(i);
+  }
+
+  void LpSoplex::_clearObj() {
+    for (int i = 0; i < soplex->nCols(); ++i) {
+      soplex->changeObj(i, 0.0);
+    }
+    solved = false;
+  }
+
+  LpSoplex::SolveExitStatus LpSoplex::_solve() {
+    soplex::SPxSolver::Status status = soplex->solve();
+
+    soplex::Vector pv(primal_value.size(), &primal_value[0]);
+    soplex->getPrimal(pv);
+
+    soplex::Vector dv(dual_value.size(), &dual_value[0]);
+    soplex->getDual(dv);
+
+    switch (status) {
+    case soplex::SPxSolver::OPTIMAL:
+    case soplex::SPxSolver::INFEASIBLE:
+    case soplex::SPxSolver::UNBOUNDED:
+      solved = true;
+      return SOLVED;
+    default:
+      return UNSOLVED;
+    }
+  }
+
+  LpSoplex::Value LpSoplex::_getPrimal(int i) const {
+    return primal_value[i];
+  }
+
+  LpSoplex::Value LpSoplex::_getDual(int i) const {
+    return dual_value[i];
+  }
+
+  LpSoplex::Value LpSoplex::_getPrimalValue() const {
+    return soplex->objValue();
+  }
+
+  bool LpSoplex::_isBasicCol(int i) const {
+    return soplex->getBasisColStatus(i) == soplex::SPxSolver::BASIC;
+  }
+
+  LpSoplex::SolutionStatus LpSoplex::_getPrimalStatus() const {
+    if (!solved) return UNDEFINED;
+    switch (soplex->status()) {
+    case soplex::SPxSolver::OPTIMAL:
+      return OPTIMAL;
+    case soplex::SPxSolver::UNBOUNDED:
+      return INFINITE;
+    case soplex::SPxSolver::INFEASIBLE:
+      return INFEASIBLE;
+    default:
+      return UNDEFINED;
+    }
+  }
+
+  LpSoplex::SolutionStatus LpSoplex::_getDualStatus() const {
+    if (!solved) return UNDEFINED;
+    switch (soplex->status()) {
+    case soplex::SPxSolver::OPTIMAL:
+      return OPTIMAL;
+    case soplex::SPxSolver::UNBOUNDED:
+      return INFEASIBLE;
+    default:
+      return UNDEFINED;
+    }
+  }
+
+  LpSoplex::ProblemTypes LpSoplex::_getProblemType() const {
+    if (!solved) return UNKNOWN;
+    switch (soplex->status()) {
+    case soplex::SPxSolver::OPTIMAL:
+      return PRIMAL_DUAL_FEASIBLE;
+    case soplex::SPxSolver::UNBOUNDED:
+      return PRIMAL_FEASIBLE_DUAL_INFEASIBLE;
+    default:
+      return UNKNOWN;
+    }
+  }
+
+  void LpSoplex::_setMax() {
+    soplex->changeSense(soplex::SPxSolver::MAXIMIZE);
+    solved = false;
+  }
+  void LpSoplex::_setMin() {
+    soplex->changeSense(soplex::SPxSolver::MINIMIZE);
+    solved = false;
+  }
+  bool LpSoplex::_isMax() const {
+    return soplex->spxSense() == soplex::SPxSolver::MAXIMIZE;
+  }
+
+
+} //namespace lemon
+
diff -r 69928a704ffb -r 7afc121e0689 lemon/lp_soplex.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/lp_soplex.h	Tue Dec 02 21:40:33 2008 +0100
@@ -0,0 +1,120 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_LP_SOPLEX_H
+#define LEMON_LP_SOPLEX_H
+
+///\file
+///\brief Header of the LEMON-SOPLEX lp solver interface.
+
+#include <vector>
+#include <string>
+
+#include <lemon/lp_base.h>
+
+// Forward declaration
+namespace soplex {
+  class SoPlex;
+}
+
+namespace lemon {
+
+  /// \ingroup lp_group
+  ///
+  /// \brief Interface for the SOPLEX solver
+  ///
+  /// This class implements an interface for the SoPlex LP solver.
+  /// The SoPlex library is an object oriented lp solver library
+  /// developed at the Konrad-Zuse-Zentrum für Informationstechnik
+  /// Berlin (ZIB). You can find detailed information about it at the
+  /// <tt>http://soplex.zib.de</tt> address.
+  class LpSoplex :virtual public LpSolverBase {
+  protected:
+
+    _lp_bits::RelocateIdHandler relocateIdHandler;
+
+    soplex::SoPlex* soplex;
+    bool solved;
+
+    std::vector<std::string> colNames;
+    std::map<std::string, int> invColNames;
+
+    std::vector<Value> primal_value;
+    std::vector<Value> dual_value;
+
+
+  public:
+
+    typedef LpSolverBase Parent;
+
+
+    /// \e
+    LpSoplex();
+    /// \e
+    LpSoplex(const LpSoplex&);
+    /// \e
+    ~LpSoplex();
+
+  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 &lower, Value &upper) const;
+    virtual void _setObjCoeff(int i, Value obj_coef);
+    virtual Value _getObjCoeff(int i) const;
+    virtual void _clearObj();
+
+    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;
+
+    virtual SolutionStatus _getPrimalStatus() const;
+    virtual SolutionStatus _getDualStatus() const;
+    virtual ProblemTypes _getProblemType() const;
+
+
+    virtual void _setMax();
+    virtual void _setMin();
+    virtual bool _isMax() const;
+
+  };
+} //END OF NAMESPACE LEMON
+
+#endif //LEMON_LP_SOPLEX_H
+
diff -r 69928a704ffb -r 7afc121e0689 lemon/mip_cplex.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/mip_cplex.cc	Tue Dec 02 21:40:33 2008 +0100
@@ -0,0 +1,141 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+///\file
+///\brief Implementation of the LEMON-CPLEX mip solver interface.
+
+#include <lemon/mip_cplex.h>
+
+extern "C" {
+#include <ilcplex/cplex.h>
+}
+
+namespace lemon {
+
+  MipCplex::MipCplex() {
+    //This is unnecessary: setting integrality constraints on
+    //variables will set this, too
+
+    ///\todo The constant CPXPROB_MIP is
+    ///called CPXPROB_MILP in later versions
+#if CPX_VERSION < 800
+    CPXchgprobtype( env,  lp, CPXPROB_MIP);
+#else
+    CPXchgprobtype( env,  lp, CPXPROB_MILP);
+#endif
+
+  }
+
+  void MipCplex::_colType(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.
+
+    int indices[1];
+    indices[0]=i;
+    char ctype[1];
+    switch (col_type){
+      case INT:
+        ctype[0]=CPX_INTEGER;//'I'
+        break;
+      case REAL:
+        ctype[0]=CPX_CONTINUOUS        ;//'C'
+        break;
+    default:;
+        //FIXME problem
+    }
+    CPXchgctype (env, lp, 1, indices, ctype);
+  }
+
+  MipCplex::ColTypes MipCplex::_colType(int i) const {
+
+    char ctype[1];
+    CPXgetctype (env, lp, ctype, i, i);
+    switch (ctype[0]){
+
+    case CPX_INTEGER:
+      return INT;
+    case CPX_CONTINUOUS:
+      return REAL;
+    default:
+      return REAL;//Error!
+    }
+
+  }
+
+  LpCplex::SolveExitStatus MipCplex::_solve(){
+
+    status = CPXmipopt (env, lp);
+    if (status==0)
+      return SOLVED;
+    else
+      return UNSOLVED;
+
+  }
+
+
+  LpCplex::SolutionStatus MipCplex::_getMipStatus() const {
+
+    int stat = CPXgetstat(env, lp);
+
+    //Fortunately, MIP statuses did not change for cplex 8.0
+    switch (stat)
+    {
+      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 INFINITE;
+      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::_getPrimal(int i) const {
+    Value x;
+    CPXgetmipx(env, lp, &x, i, i);
+    return x;
+  }
+
+  MipCplex::Value MipCplex::_getPrimalValue() const {
+    Value objval;
+    CPXgetmipobjval(env, lp, &objval);
+    return objval;
+  }
+} //END OF NAMESPACE LEMON
diff -r 69928a704ffb -r 7afc121e0689 lemon/mip_cplex.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/mip_cplex.h	Tue Dec 02 21:40:33 2008 +0100
@@ -0,0 +1,61 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_MIP_CPLEX_H
+#define LEMON_MIP_CPLEX_H
+
+///\file
+///\brief Header of the LEMON-CPLEX mip solver interface.
+///\ingroup lp_group
+
+
+#include <lemon/lp_cplex.h>
+
+namespace lemon {
+
+  /// \brief Interface for the CPLEX MIP solver
+  ///
+  /// This class implements an interface for the CPLEX MIP solver.
+  ///\ingroup lp_group
+  class MipCplex : public MipSolverBase, public LpCplex{
+
+  public:
+
+    typedef MipSolverBase ParentMip;
+    typedef LpCplex ParentLp;
+
+    MipCplex();
+    //~MipCplex();
+
+
+
+
+  protected:
+
+    virtual ColTypes _colType(int col) const;
+    virtual void _colType(int col, ColTypes col_type);
+
+    virtual LpCplex::SolveExitStatus _solve();
+    virtual LpCplex::SolutionStatus _getMipStatus() const;
+    virtual ParentLp::Value _getPrimal(int i) const;
+    virtual ParentLp::Value _getPrimalValue() const;
+  };
+
+} //END OF NAMESPACE LEMON
+
+#endif // END OF LEMON_MIP_CPLEX_H
diff -r 69928a704ffb -r 7afc121e0689 lemon/mip_glpk.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/mip_glpk.cc	Tue Dec 02 21:40:33 2008 +0100
@@ -0,0 +1,154 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+///\file
+///\brief Implementation of the LEMON-GLPK mip solver interface.
+
+#include <lemon/mip_glpk.h>
+
+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
+
+namespace lemon {
+
+  MipGlpk::MipGlpk() {
+#if !(GLP_MAJOR_VERSION > 4 || \
+      (GLP_MAJOR_VERSION == 4 && GLP_MINOR_VERSION > 15))
+    LEMON_lpx(set_class)(lp,LEMON_GLP(MIP));
+#endif
+  }
+
+  void MipGlpk::_colType(int i, MipGlpk::ColTypes col_type){
+    switch (col_type){
+      case INT:
+        LEMON_glp(set_col_kind)(lp,i,LEMON_GLP(IV));
+        break;
+      case REAL:
+        LEMON_glp(set_col_kind)(lp,i,LEMON_GLP(CV));
+        break;
+    default:;
+        //FIXME problem
+    }
+  }
+
+  MipGlpk::ColTypes MipGlpk::_colType(int i) const {
+    switch (LEMON_glp(get_col_kind)(lp,i)){
+    case LEMON_GLP(IV):
+      return INT;//Or binary
+    case LEMON_GLP(CV):
+      return REAL;
+    default:
+      return REAL;//Error!
+    }
+
+  }
+
+  LpGlpk::SolveExitStatus MipGlpk::_solve() {
+    int result = LEMON_lpx(simplex)(lp);
+
+    // hack: mip does not contain integer variable
+#if GLP_MAJOR_VERSION == 4 && GLP_MINOR_VERSION == 16
+    int tmp = -1;
+    if (LEMON_glp(get_num_int(lp)) == 0) {
+      tmp = LEMON_lpx(add_cols)(lp, 1);
+      LEMON_glp(set_col_bnds)(lp, tmp, LEMON_GLP(FX), 0.0, 0.0);
+      LEMON_glp(set_col_kind)(lp, tmp, LEMON_GLP(IV));
+    }
+#endif
+
+    if (LEMON_lpx(get_status)(lp)==LEMON_LPX(OPT)) {
+      //Maybe we could try the routine lpx_intopt(lp), a revised
+      //version of lpx_integer
+
+      result = LEMON_lpx(integer)(lp);
+      switch (result){
+      case LEMON_LPX(E_OK):
+        solved = true;
+        break;
+      default:
+        solved = false;
+      }
+    } else {
+      solved = false;
+    }
+#if GLP_MAJOR_VERSION == 4 && GLP_MINOR_VERSION == 16
+    if (tmp != -1) {
+      int tmpa[2];
+      tmpa[1] = tmp;
+      LEMON_lpx(del_cols)(lp, 1, tmpa);
+    }
+#endif
+    return solved ? SOLVED : UNSOLVED;
+  }
+
+
+  LpGlpk::SolutionStatus MipGlpk::_getMipStatus() const {
+
+    if (LEMON_lpx(get_status)(lp)==LEMON_LPX(OPT)){
+      //Meg kell nezni: ha az LP is infinite, akkor ez is, ha az is
+      //infeasible, akkor ez is, de ez lehet maskepp is infeasible.
+      int stat= LEMON_lpx(mip_status)(lp);
+
+      switch (stat) {
+      case LEMON_LPX(I_UNDEF)://Undefined (no solve has been run yet)
+        return UNDEFINED;
+      case LEMON_LPX(I_NOFEAS)://There is no feasible integral solution
+        return INFEASIBLE;
+        //     case LEMON_LPX(UNBND)://Unbounded
+        //       return INFINITE;
+      case LEMON_LPX(I_FEAS)://Feasible
+        return FEASIBLE;
+      case LEMON_LPX(I_OPT)://Feasible
+        return OPTIMAL;
+      default:
+        return UNDEFINED; //to avoid gcc warning
+      //FIXME error
+      }
+    }
+    else
+      return UNDEFINED; //Maybe we could refine this: what does the LP
+                        //relaxation look like
+
+  }
+
+  MipGlpk::Value MipGlpk::_getPrimal(int i) const {
+    return LEMON_glp(mip_col_val)(lp,i);
+  }
+
+  MipGlpk::Value MipGlpk::_getPrimalValue() const {
+    return LEMON_glp(mip_obj_val)(lp);
+  }
+} //END OF NAMESPACE LEMON
diff -r 69928a704ffb -r 7afc121e0689 lemon/mip_glpk.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/mip_glpk.h	Tue Dec 02 21:40:33 2008 +0100
@@ -0,0 +1,59 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_MIP_GLPK_H
+#define LEMON_MIP_GLPK_H
+
+///\file
+///\brief Header of the LEMON-GLPK mip solver interface.
+///\ingroup lp_group
+
+
+#include <lemon/lp_glpk.h>
+
+namespace lemon {
+  /// \brief Interface for the GLPK MIP solver
+  ///
+  /// This class implements an interface for the GLPK MIP solver.
+  ///\ingroup lp_group
+  class MipGlpk : public MipSolverBase, public LpGlpk{
+
+  public:
+
+    typedef MipSolverBase ParentMip;
+    typedef LpGlpk ParentLp;
+
+    MipGlpk();
+    //~MipGlpk();
+
+
+
+  protected:
+
+    virtual ColTypes _colType(int col) const;
+    virtual void _colType(int col, ColTypes col_type);
+
+    virtual LpGlpk::SolveExitStatus _solve();
+    virtual LpGlpk::SolutionStatus _getMipStatus() const;
+    virtual ParentLp::Value _getPrimal(int i) const;
+    virtual ParentLp::Value _getPrimalValue() const;
+  };
+
+} //END OF NAMESPACE LEMON
+
+#endif // END OF LEMON_MIP_GLPK_H
diff -r 69928a704ffb -r 7afc121e0689 test/CMakeLists.txt
--- a/test/CMakeLists.txt	Mon Jan 12 12:22:11 2009 +0000
+++ b/test/CMakeLists.txt	Tue Dec 02 21:40:33 2008 +0100
@@ -18,6 +18,8 @@
   hao_orlin_test
   heap_test
   kruskal_test
+  lp_test
+  mip_test
   maps_test
   max_matching_test
   radix_sort_test
diff -r 69928a704ffb -r 7afc121e0689 test/Makefile.am
--- a/test/Makefile.am	Mon Jan 12 12:22:11 2009 +0000
+++ b/test/Makefile.am	Tue Dec 02 21:40:33 2008 +0100
@@ -33,6 +33,13 @@
 	test/time_measure_test \
 	test/unionfind_test
 
+if HAVE_LP
+check_PROGRAMS += test/lp_test
+endif HAVE_LP
+if HAVE_MIP
+check_PROGRAMS += test/mip_test
+endif HAVE_MIP
+
 TESTS += $(check_PROGRAMS)
 XFAIL_TESTS += test/test_tools_fail$(EXEEXT)
 
@@ -51,7 +58,9 @@
 test_heap_test_SOURCES = test/heap_test.cc
 test_kruskal_test_SOURCES = test/kruskal_test.cc
 test_hao_orlin_test_SOURCES = test/hao_orlin_test.cc
+test_lp_test_SOURCES = test/lp_test.cc
 test_maps_test_SOURCES = test/maps_test.cc
+test_mip_test_SOURCES = test/mip_test.cc
 test_max_matching_test_SOURCES = test/max_matching_test.cc
 test_path_test_SOURCES = test/path_test.cc
 test_preflow_test_SOURCES = test/preflow_test.cc
diff -r 69928a704ffb -r 7afc121e0689 test/lp_test.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/lp_test.cc	Tue Dec 02 21:40:33 2008 +0100
@@ -0,0 +1,423 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#include <sstream>
+#include <lemon/lp_skeleton.h>
+#include "test_tools.h"
+#include <lemon/tolerance.h>
+
+#ifdef HAVE_CONFIG_H
+#include <lemon/config.h>
+#endif
+
+#ifdef HAVE_GLPK
+#include <lemon/lp_glpk.h>
+#endif
+
+#ifdef HAVE_CPLEX
+#include <lemon/lp_cplex.h>
+#endif
+
+#ifdef HAVE_SOPLEX
+#include <lemon/lp_soplex.h>
+#endif
+
+using namespace lemon;
+
+void lpTest(LpSolverBase & lp)
+{
+
+
+
+  typedef LpSolverBase LP;
+
+  std::vector<LP::Col> x(10);
+  //  for(int i=0;i<10;i++) x.push_back(lp.addCol());
+  lp.addColSet(x);
+  lp.colLowerBound(x,1);
+  lp.colUpperBound(x,1);
+  lp.colBounds(x,1,2);
+#ifndef GYORSITAS
+
+  std::vector<LP::Col> y(10);
+  lp.addColSet(y);
+
+  lp.colLowerBound(y,1);
+  lp.colUpperBound(y,1);
+  lp.colBounds(y,1,2);
+
+  std::map<int,LP::Col> z;
+
+  z.insert(std::make_pair(12,INVALID));
+  z.insert(std::make_pair(2,INVALID));
+  z.insert(std::make_pair(7,INVALID));
+  z.insert(std::make_pair(5,INVALID));
+
+  lp.addColSet(z);
+
+  lp.colLowerBound(z,1);
+  lp.colUpperBound(z,1);
+  lp.colBounds(z,1,2);
+
+  {
+    LP::Expr e,f,g;
+    LP::Col p1,p2,p3,p4,p5;
+    LP::Constr c;
+
+    p1=lp.addCol();
+    p2=lp.addCol();
+    p3=lp.addCol();
+    p4=lp.addCol();
+    p5=lp.addCol();
+
+    e[p1]=2;
+    e.constComp()=12;
+    e[p1]+=2;
+    e.constComp()+=12;
+    e[p1]-=2;
+    e.constComp()-=12;
+
+    e=2;
+    e=2.2;
+    e=p1;
+    e=f;
+
+    e+=2;
+    e+=2.2;
+    e+=p1;
+    e+=f;
+
+    e-=2;
+    e-=2.2;
+    e-=p1;
+    e-=f;
+
+    e*=2;
+    e*=2.2;
+    e/=2;
+    e/=2.2;
+
+    e=((p1+p2)+(p1-p2)+(p1+12)+(12+p1)+(p1-12)+(12-p1)+
+       (f+12)+(12+f)+(p1+f)+(f+p1)+(f+g)+
+       (f-12)+(12-f)+(p1-f)+(f-p1)+(f-g)+
+       2.2*f+f*2.2+f/2.2+
+       2*f+f*2+f/2+
+       2.2*p1+p1*2.2+p1/2.2+
+       2*p1+p1*2+p1/2
+       );
+
+
+    c = (e  <= f  );
+    c = (e  <= 2.2);
+    c = (e  <= 2  );
+    c = (e  <= p1 );
+    c = (2.2<= f  );
+    c = (2  <= f  );
+    c = (p1 <= f  );
+    c = (p1 <= p2 );
+    c = (p1 <= 2.2);
+    c = (p1 <= 2  );
+    c = (2.2<= p2 );
+    c = (2  <= p2 );
+
+    c = (e  >= f  );
+    c = (e  >= 2.2);
+    c = (e  >= 2  );
+    c = (e  >= p1 );
+    c = (2.2>= f  );
+    c = (2  >= f  );
+    c = (p1 >= f  );
+    c = (p1 >= p2 );
+    c = (p1 >= 2.2);
+    c = (p1 >= 2  );
+    c = (2.2>= p2 );
+    c = (2  >= p2 );
+
+    c = (e  == f  );
+    c = (e  == 2.2);
+    c = (e  == 2  );
+    c = (e  == p1 );
+    c = (2.2== f  );
+    c = (2  == f  );
+    c = (p1 == f  );
+    //c = (p1 == p2 );
+    c = (p1 == 2.2);
+    c = (p1 == 2  );
+    c = (2.2== p2 );
+    c = (2  == p2 );
+
+    c = (2 <= e <= 3);
+    c = (2 <= p1<= 3);
+
+    c = (2 >= e >= 3);
+    c = (2 >= p1>= 3);
+
+    e[x[3]]=2;
+    e[x[3]]=4;
+    e[x[3]]=1;
+    e.constComp()=12;
+
+    lp.addRow(LP::INF,e,23);
+    lp.addRow(LP::INF,3.0*(x[1]+x[2]/2)-x[3],23);
+    lp.addRow(LP::INF,3.0*(x[1]+x[2]*2-5*x[3]+12-x[4]/3)+2*x[4]-4,23);
+
+    lp.addRow(x[1]+x[3]<=x[5]-3);
+    lp.addRow(-7<=x[1]+x[3]-12<=3);
+    lp.addRow(x[1]<=x[5]);
+
+    std::ostringstream buf;
+
+
+    //Checking the simplify function
+
+//     //How to check the simplify function? A map gives no information
+//     //on the question whether a given key is or is not stored in it, or
+//     //it does?
+//   Yes, it does, using the find() function.
+    e=((p1+p2)+(p1-p2));
+    e.simplify();
+    buf << "Coeff. of p2 should be 0";
+    //    std::cout<<e[p1]<<e[p2]<<e[p3]<<std::endl;
+    check(e.find(p2)==e.end(), buf.str());
+
+
+
+
+    e=((p1+p2)+(p1-0.99*p2));
+    //e.prettyPrint(std::cout);
+    //(e<=2).prettyPrint(std::cout);
+    double tolerance=0.001;
+    e.simplify(tolerance);
+    buf << "Coeff. of p2 should be 0.01";
+    check(e[p2]>0, buf.str());
+
+    tolerance=0.02;
+    e.simplify(tolerance);
+    buf << "Coeff. of p2 should be 0";
+    check(e.find(p2)==e.end(), buf.str());
+
+
+  }
+
+  {
+    LP::DualExpr e,f,g;
+    LP::Row p1 = INVALID, p2 = INVALID, p3 = INVALID,
+      p4 = INVALID, p5 = INVALID;
+
+    e[p1]=2;
+    e[p1]+=2;
+    e[p1]-=2;
+
+    e=p1;
+    e=f;
+
+    e+=p1;
+    e+=f;
+
+    e-=p1;
+    e-=f;
+
+    e*=2;
+    e*=2.2;
+    e/=2;
+    e/=2.2;
+
+    e=((p1+p2)+(p1-p2)+
+       (p1+f)+(f+p1)+(f+g)+
+       (p1-f)+(f-p1)+(f-g)+
+       2.2*f+f*2.2+f/2.2+
+       2*f+f*2+f/2+
+       2.2*p1+p1*2.2+p1/2.2+
+       2*p1+p1*2+p1/2
+       );
+  }
+
+#endif
+}
+
+void solveAndCheck(LpSolverBase& lp, LpSolverBase::SolutionStatus stat,
+                   double exp_opt) {
+  using std::string;
+  lp.solve();
+  //int decimal,sign;
+  std::ostringstream buf;
+  buf << "Primalstatus should be: " << int(stat);
+
+  //  itoa(stat,buf1, 10);
+  check(lp.primalStatus()==stat, buf.str());
+
+  if (stat ==  LpSolverBase::OPTIMAL) {
+    std::ostringstream sbuf;
+    sbuf << "Wrong optimal value: the right optimum is " << exp_opt;
+    check(std::abs(lp.primalValue()-exp_opt) < 1e-3, sbuf.str());
+    //+ecvt(exp_opt,2)
+  }
+}
+
+void aTest(LpSolverBase & lp)
+{
+  typedef LpSolverBase LP;
+
+ //The following example is very simple
+
+  typedef LpSolverBase::Row Row;
+  typedef LpSolverBase::Col Col;
+
+
+  Col x1 = lp.addCol();
+  Col x2 = lp.addCol();
+
+
+  //Constraints
+  Row upright=lp.addRow(x1+x2 <=1);
+  lp.addRow(x1+x2 >=-1);
+  lp.addRow(x1-x2 <=1);
+  lp.addRow(x1-x2 >=-1);
+  //Nonnegativity of the variables
+  lp.colLowerBound(x1, 0);
+  lp.colLowerBound(x2, 0);
+  //Objective function
+  lp.obj(x1+x2);
+
+  lp.max();
+
+  //Testing the problem retrieving routines
+  check(lp.objCoeff(x1)==1,"First term should be 1 in the obj function!");
+  check(lp.isMax(),"This is a maximization!");
+  check(lp.coeff(upright,x1)==1,"The coefficient in question is 1!");
+  //  std::cout<<lp.colLowerBound(x1)<<std::endl;
+  check(  lp.colLowerBound(x1)==0,
+          "The lower bound for variable x1 should be 0.");
+  check(  lp.colUpperBound(x1)==LpSolverBase::INF,
+          "The upper bound for variable x1 should be infty.");
+  LpSolverBase::Value lb,ub;
+  lp.getRowBounds(upright,lb,ub);
+  check(  lb==-LpSolverBase::INF,
+          "The lower bound for the first row should be -infty.");
+  check(  ub==1,"The upper bound for the first row should be 1.");
+  LpSolverBase::Expr e = lp.row(upright);
+  check(  e.size() == 2, "The row retrieval gives back wrong expression.");
+  check(  e[x1] == 1, "The first coefficient should 1.");
+  check(  e[x2] == 1, "The second coefficient should 1.");
+
+  LpSolverBase::DualExpr de = lp.col(x1);
+  check(  de.size() == 4, "The col retrieval gives back wrong expression.");
+  check(  de[upright] == 1, "The first coefficient should 1.");
+
+  LpSolverBase* clp = lp.copyLp();
+
+  //Testing the problem retrieving routines
+  check(clp->objCoeff(x1)==1,"First term should be 1 in the obj function!");
+  check(clp->isMax(),"This is a maximization!");
+  check(clp->coeff(upright,x1)==1,"The coefficient in question is 1!");
+  //  std::cout<<lp.colLowerBound(x1)<<std::endl;
+  check(  clp->colLowerBound(x1)==0,
+          "The lower bound for variable x1 should be 0.");
+  check(  clp->colUpperBound(x1)==LpSolverBase::INF,
+          "The upper bound for variable x1 should be infty.");
+
+  clp->getRowBounds(upright,lb,ub);
+  check(  lb==-LpSolverBase::INF,
+          "The lower bound for the first row should be -infty.");
+  check(  ub==1,"The upper bound for the first row should be 1.");
+  e = clp->row(upright);
+  check(  e.size() == 2, "The row retrieval gives back wrong expression.");
+  check(  e[x1] == 1, "The first coefficient should 1.");
+  check(  e[x2] == 1, "The second coefficient should 1.");
+
+  de = clp->col(x1);
+  check(  de.size() == 4, "The col retrieval gives back wrong expression.");
+  check(  de[upright] == 1, "The first coefficient should 1.");
+
+  delete clp;
+
+  //Maximization of x1+x2
+  //over the triangle with vertices (0,0) (0,1) (1,0)
+  double expected_opt=1;
+  solveAndCheck(lp, LpSolverBase::OPTIMAL, expected_opt);
+
+  //Minimization
+  lp.min();
+  expected_opt=0;
+  solveAndCheck(lp, LpSolverBase::OPTIMAL, expected_opt);
+
+  //Vertex (-1,0) instead of (0,0)
+  lp.colLowerBound(x1, -LpSolverBase::INF);
+  expected_opt=-1;
+  solveAndCheck(lp, LpSolverBase::OPTIMAL, expected_opt);
+
+  //Erase one constraint and return to maximization
+  lp.eraseRow(upright);
+  lp.max();
+  expected_opt=LpSolverBase::INF;
+  solveAndCheck(lp, LpSolverBase::INFINITE, expected_opt);
+
+  //Infeasibilty
+  lp.addRow(x1+x2 <=-2);
+  solveAndCheck(lp, LpSolverBase::INFEASIBLE, expected_opt);
+
+  //Change problem and forget to solve
+  lp.min();
+  check(lp.primalStatus()==LpSolverBase::UNDEFINED,
+        "Primalstatus should be UNDEFINED");
+
+
+//   lp.solve();
+//   if (lp.primalStatus()==LpSolverBase::OPTIMAL){
+//     std::cout<< "Z = "<<lp.primalValue()
+//              << " (error = " << lp.primalValue()-expected_opt
+//              << "); x1 = "<<lp.primal(x1)
+//              << "; x2 = "<<lp.primal(x2)
+//              <<std::endl;
+
+//   }
+//   else{
+//     std::cout<<lp.primalStatus()<<std::endl;
+//     std::cout<<"Optimal solution not found!"<<std::endl;
+//   }
+
+
+
+}
+
+
+int main()
+{
+  LpSkeleton lp_skel;
+  lpTest(lp_skel);
+
+#ifdef HAVE_GLPK
+  LpGlpk lp_glpk1,lp_glpk2;
+  lpTest(lp_glpk1);
+  aTest(lp_glpk2);
+#endif
+
+#ifdef HAVE_CPLEX
+  LpCplex lp_cplex1,lp_cplex2;
+  lpTest(lp_cplex1);
+  aTest(lp_cplex2);
+#endif
+
+#ifdef HAVE_SOPLEX
+  LpSoplex lp_soplex1,lp_soplex2;
+  lpTest(lp_soplex1);
+  aTest(lp_soplex2);
+#endif
+
+  return 0;
+}
diff -r 69928a704ffb -r 7afc121e0689 test/mip_test.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/mip_test.cc	Tue Dec 02 21:40:33 2008 +0100
@@ -0,0 +1,126 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#include "test_tools.h"
+
+
+#ifdef HAVE_CONFIG_H
+#include <lemon/config.h>
+#endif
+
+#ifdef HAVE_CPLEX
+#include <lemon/mip_cplex.h>
+#endif
+
+#ifdef HAVE_GLPK
+#include <lemon/mip_glpk.h>
+#endif
+
+
+using namespace lemon;
+
+void solveAndCheck(MipSolverBase& lp, MipSolverBase::SolutionStatus stat,
+                   double exp_opt) {
+  using std::string;
+
+  lp.solve();
+  //int decimal,sign;
+  std::ostringstream buf;
+  buf << "Primalstatus should be: " << int(stat)
+      <<" and it is "<<int(lp.mipStatus());
+
+
+  //  itoa(stat,buf1, 10);
+  check(lp.mipStatus()==stat, buf.str());
+
+  if (stat ==  MipSolverBase::OPTIMAL) {
+    std::ostringstream sbuf;
+    buf << "Wrong optimal value: the right optimum is " << exp_opt;
+    check(std::abs(lp.primalValue()-exp_opt) < 1e-3, sbuf.str());
+    //+ecvt(exp_opt,2)
+  }
+}
+
+void aTest(MipSolverBase& mip)
+{
+ //The following example is very simple
+
+
+  typedef MipSolverBase::Row Row;
+  typedef MipSolverBase::Col Col;
+
+
+
+  Col x1 = mip.addCol();
+  Col x2 = mip.addCol();
+
+
+  //Objective function
+  mip.obj(x1);
+
+  mip.max();
+
+
+  //Unconstrained optimization
+  mip.solve();
+  //Check it out!
+
+  //Constraints
+  mip.addRow(2*x1+x2 <=2);
+  mip.addRow(x1-2*x2 <=0);
+
+  //Nonnegativity of the variable x1
+  mip.colLowerBound(x1, 0);
+
+  //Maximization of x1
+  //over the triangle with vertices (0,0),(4/5,2/5),(0,2)
+  double expected_opt=4.0/5.0;
+  solveAndCheck(mip, MipSolverBase::OPTIMAL, expected_opt);
+
+  //Restrict x2 to integer
+  mip.colType(x2,MipSolverBase::INT);
+  expected_opt=1.0/2.0;
+  solveAndCheck(mip, MipSolverBase::OPTIMAL, expected_opt);
+
+
+  //Restrict both to integer
+  mip.colType(x1,MipSolverBase::INT);
+  expected_opt=0;
+  solveAndCheck(mip, MipSolverBase::OPTIMAL, expected_opt);
+
+
+
+}
+
+
+int main()
+{
+
+#ifdef HAVE_GLPK
+  MipGlpk mip1;
+  aTest(mip1);
+#endif
+
+#ifdef HAVE_CPLEX
+  MipCplex mip2;
+  aTest(mip2);
+#endif
+
+  return 0;
+
+}