Ticket #44: aa6dff56b2e1.patch

Makefile.am

@@ -5 +5 @@
 
 EXTRA_DIST = \
         LICENSE \
+        m4/lx_check_clp.m4 \
         m4/lx_check_cplex.m4 \
         m4/lx_check_glpk.m4 \
         m4/lx_check_soplex.m4 \

configure.ac

@@ -37 +37 @@
 fi
 
 dnl Checks for libraries.
-#LX_CHECK_GLPK
-#LX_CHECK_CPLEX
-#LX_CHECK_SOPLEX
+LX_CHECK_GLPK
+LX_CHECK_CPLEX
+LX_CHECK_SOPLEX
+LX_CHECK_CLP
 
 dnl Disable/enable building the demo programs.
 AC_ARG_ENABLE([demo],
@@ -114 +115 @@
 echo C++ compiler.................. : $CXX
 echo C++ compiles flags............ : $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 CLP support................... : $lx_clp_found
+echo
 echo Build benchmarks.............. : $enable_benchmark
 echo Build demo programs........... : $enable_demo
 echo Build additional tools........ : $enable_tools

lemon/Makefile.am

@@ -10 +10 @@
         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)
-#lemon_libemon_la_LDFLAGS = $(GLPK_LIBS) $(CPLEX_LIBS) $(SOPLEX_LIBS)
+
+lemon_libemon_la_CXXFLAGS = \
+        $(GLPK_CFLAGS) \
+        $(CPLEX_CFLAGS) \
+        $(SOPLEX_CXXFLAGS) \
+        $(CLP_CXXFLAGS)
+
+lemon_libemon_la_LDFLAGS = \
+        $(GLPK_LIBS) \
+        $(CPLEX_LIBS) \
+        $(SOPLEX_LIBS) \
+        $(CLP_LIBS)
+
+if HAVE_GLPK
+lemon_libemon_la_SOURCES += lemon/lp_glpk.cc
+endif
+
+if HAVE_CPLEX
+lemon_libemon_la_SOURCES += lemon/lp_cplex.cc
+endif
+
+if HAVE_SOPLEX
+lemon_libemon_la_SOURCES += lemon/lp_soplex.cc
+endif
+
+if HAVE_CLP
+lemon_libemon_la_SOURCES += lemon/lp_clp.cc
+endif
 
 lemon_HEADERS += \
         lemon/arg_parser.h \
@@ -32 +60 @@
         lemon/kruskal.h \
         lemon/lgf_reader.h \
         lemon/lgf_writer.h \
+        lemon/lp_base.h \
+        lemon/lp_cplex.h \
+        lemon/lp_glpk.h \
         lemon/list_graph.h \
         lemon/maps.h \
         lemon/math.h \
@@ -52 +83 @@
         lemon/bits/graph_extender.h \
         lemon/bits/map_extender.h \
         lemon/bits/path_dump.h \
+        lemon/bits/solver_bits.h \
         lemon/bits/traits.h \
         lemon/bits/vector_map.h
 

new file lemon/bits/solver_bits.h

@@ +1 @@
+/* -*- C++ -*-
+ *
+ * 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_SOLVER_BITS_H
+#define LEMON_BITS_SOLVER_BITS_H
+
+namespace lemon {
+  
+  namespace _solver_bits {
+
+    class VarIndex {
+    private:
+      struct ItemT {
+        int prev, next;
+        int index;
+      };
+      std::vector<ItemT> items;
+      int first_item, last_item, first_free_item;
+      
+      std::vector<int> cross;
+
+    public:
+
+      VarIndex() 
+        : first_item(-1), last_item(-1), first_free_item(-1) {
+      }
+
+      void clear() {
+        first_item = -1;
+        first_free_item = -1;
+        items.clear();
+        cross.clear();
+      }
+
+      int addIndex(int idx) {
+        int n;
+        if (first_free_item == -1) {
+          n = items.size();
+          items.push_back(ItemT());
+        } else {
+          n = first_free_item;
+          first_free_item = items[n].next;
+          if (first_free_item != -1) {
+            items[first_free_item].prev = -1;
+          }
+        }
+        items[n].index = idx;
+        if (static_cast<int>(cross.size()) <= idx) {
+          cross.resize(idx + 1, -1);
+        }
+        cross[idx] = n;
+
+        items[n].prev = last_item;
+        items[n].next = -1;
+        if (last_item != -1) {
+          items[last_item].next = n;
+        } else {
+          first_item = n;
+        }
+        last_item = n;
+        
+        return n;
+      }
+
+      int addIndex(int idx, int n) {
+        while (n >= static_cast<int>(items.size())) {
+          items.push_back(ItemT());
+          items.back().prev = -1;
+          items.back().next = first_free_item;
+          if (first_free_item != -1) {
+            items[first_free_item].prev = items.size() - 1;
+          }
+          first_free_item = items.size() - 1;
+        }
+        if (items[n].next != -1) {
+          items[items[n].next].prev = items[n].prev;
+        }
+        if (items[n].prev != -1) {
+          items[items[n].prev].next = items[n].next;
+        } else {
+          first_free_item = items[n].next;
+        }
+
+        items[n].index = idx;
+        if (static_cast<int>(cross.size()) <= idx) {
+          cross.resize(idx + 1, -1);
+        }
+        cross[idx] = n;
+
+        items[n].prev = last_item;
+        items[n].next = -1;
+        if (last_item != -1) {
+          items[last_item].next = n;
+        } else {
+          first_item = n;
+        }
+        last_item = n;
+        
+        return n;
+      }
+
+      void eraseIndex(int idx) {
+        int n = cross[idx];
+        
+        if (items[n].prev != -1) {
+          items[items[n].prev].next = items[n].next;
+        } else {
+          first_item = items[n].next;
+        }
+        if (items[n].next != -1) {
+          items[items[n].next].prev = items[n].prev;
+        } else {
+          last_item = items[n].prev;
+        }
+
+        if (first_free_item != -1) {
+          items[first_free_item].prev = n;
+        }
+        items[n].next = first_free_item;
+        items[n].prev = -1;
+        first_free_item = n;
+
+        while (!cross.empty() && cross.back() == -1) {
+          cross.pop_back();
+        }
+      }
+
+      int maxIndex() const {
+        return cross.size() - 1;
+      }
+      
+      void shiftIndexes(int idx) {
+        for (int i = idx + 1; i < static_cast<int>(cross.size()); ++i) {
+          cross[i - 1] = cross[i];
+          if (cross[i] != -1) {
+            --items[cross[i]].index;
+          }
+        }
+        cross.back() = -1;
+        cross.pop_back();
+        while (!cross.empty() && cross.back() == -1) {
+          cross.pop_back();
+        }
+      }
+
+      void relocateIndex(int idx, int jdx) {
+        cross[idx] = cross[jdx];
+        items[cross[jdx]].index = idx;
+        cross[jdx] = -1;
+
+        while (!cross.empty() && cross.back() == -1) {
+          cross.pop_back();
+        }
+      }
+
+      int operator[](int idx) const {
+        return cross[idx];
+      }
+
+      int operator()(int fdx) const {
+        return items[fdx].index;
+      }
+
+      void firstItem(int& fdx) const {
+        fdx = first_item;
+      }
+
+      void nextItem(int& fdx) const {
+        fdx = items[fdx].next;
+      }
+
+    };
+  }  
+}
+
+#endif

new file lemon/lp.h

@@ +1 @@
+/* -*- C++ -*-
+ *
+ * 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>
+#elif HAVE_CPLEX
+#include <lemon/lp_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

new file lemon/lp_base.cc

@@ +1 @@
+/* -*- 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
+  SolverBase::INF = std::numeric_limits<Value>::infinity();
+  const LpSolverBase::Value
+  SolverBase::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

new file lemon/lp_base.h

@@ +1 @@
+/* -*- C++ -*-
+ *
+ * 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/error.h>
+#include<lemon/assert.h>
+
+#include<lemon/core.h>
+#include<lemon/bits/solver_bits.h>
+
+///\file
+///\brief The interface of the LP solver interface.
+///\ingroup lp_group
+namespace lemon {
+
+  ///Common base class for LP solvers
+  
+  ///\todo Much more docs
+  ///\ingroup lp_group
+  class SolverBase {
+
+  protected:
+
+    _solver_bits::VarIndex rows;
+    _solver_bits::VarIndex 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
+    };
+
+    ///Direction of the optimization
+    enum Sense {
+      /// Minimization
+      MIN,
+      /// Maximization
+      MAX
+    };
+      
+    ///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;
+    
+    friend class Col;
+    friend class ColIt;
+    friend class Row;
+    friend class RowIt;
+    
+    ///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 {
+      friend class SolverBase;
+    protected:
+      int _id;
+      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;}
+    };
+
+    class ColIt : public Col {
+      const SolverBase *_solver;
+    public:
+      ColIt() {}
+      ColIt(const SolverBase &solver) : _solver(&solver)
+      {
+        _solver->cols.firstItem(_id);
+      }
+      ColIt(const Invalid&) : Col(INVALID) {}
+      ColIt &operator++() 
+      {
+        _solver->cols.nextItem(_id);
+        return *this;
+      }
+    };
+
+    static int id(const Col& col) { return col._id; }
+    static Col colFromId(int id) { 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 {
+      friend class SolverBase;
+    protected:
+      int _id;
+      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;} 
+    };
+
+    class RowIt : public Row {
+      const SolverBase *_solver;
+    public:
+      RowIt() {}
+      RowIt(const SolverBase &solver) : _solver(&solver)
+      {
+        _solver->rows.firstItem(_id);
+      }
+      RowIt(const Invalid&) : Row(INVALID) {}
+      RowIt &operator++() 
+      {
+        _solver->rows.nextItem(_id);
+        return *this;
+      }
+    };
+
+    static int id(const Row& row) { return row._id; }
+    static Row rowFromId(int id) { return Row(id); }
+
+  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.
+    ///\code
+    ///e[v]=5;
+    ///e[v]+=12;
+    ///e.erase(v);
+    ///\endcode
+    ///or you can also iterate through its elements.
+    ///\code
+    ///double s=0;
+    ///for(SolverBase::Expr::ConstMapIt i(e);i!=INVALID;++i)
+    ///  s+=*i * primal(i);
+    ///\endcode
+    ///(This code computes the primal value of the expression).
+    ///- 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 operator*()
+    ///\code
+    ///*e=12;
+    ///double c=*e;
+    ///\endcode
+    ///
+    ///\note \ref clear() not only sets all coefficients to 0 but also
+    ///clears the constant components.
+    ///
+    ///\sa Constr
+    ///
+    class Expr {
+      friend class SolverBase;
+    public:
+      typedef SolverBase::Col Key; 
+      typedef SolverBase::Value Value;
+      
+    protected:
+      Value const_comp;
+      std::map<int, Value> comps;
+
+    public:
+      typedef True SolverExpr;
+      ///\e
+      Expr() : const_comp(0) {}
+      ///\e
+      Expr(const Key &k) : const_comp(0) {
+        typedef std::map<int, Value>::value_type pair_type;
+        comps.insert(pair_type(id(k), 1));
+      }
+      ///\e
+      Expr(const Value &v) : const_comp(v) {}
+      ///\e
+      Value operator[](const Key& k) const {
+        std::map<int, Value>::const_iterator it=comps.find(id(k));
+        if (it != comps.end()) {
+          return it->second;
+        } else {
+          return 0;
+        }
+      }
+      ///\e
+      Value& operator[](const Key& k) {
+        return comps[id(k)];
+      }
+      ///\e
+      void set(const Key &k, const Value &v) {
+        if (v != 0.0) {
+          typedef std::map<int, Value>::value_type pair_type;
+          comps.insert(pair_type(id(k), v));
+        } else {
+          comps.erase(id(k));
+        }
+      }
+
+      ///\e
+      Value& operator*() { return const_comp; }
+      ///\e
+      const Value& operator*() const { return const_comp; }
+      
+      ///Removes the coefficients closer to zero than \c tolerance.
+      void simplify(Value tolerance = 0.0) {
+        std::map<int, Value>::iterator it=comps.begin();
+        while (it != comps.end()) { 
+          std::map<int, Value>::iterator jt=it;
+          ++jt;
+          if (std::fabs((*it).second) <= tolerance) comps.erase(it);
+          it=jt;
+        }
+        if (std::fabs(const_comp) <= tolerance) const_comp = 0;
+      }
+
+      void simplify(Value tolerance = 0.0) const {
+        const_cast<Expr*>(this)->simplify(tolerance);
+      }
+
+      ///Sets all coefficients and the constant component to 0.
+      void clear() {
+        comps.clear();
+        const_comp=0;
+      }
+
+      ///\e
+      Expr &operator+=(const Expr &e) {
+        for (std::map<int, Value>::const_iterator it=e.comps.begin(); 
+             it!=e.comps.end(); ++it)
+          comps[it->first]+=it->second;
+        const_comp+=e.const_comp;
+        return *this;
+      }
+      ///\e
+      Expr &operator-=(const Expr &e) {
+        for (std::map<int, Value>::const_iterator it=e.comps.begin(); 
+             it!=e.comps.end(); ++it)
+          comps[it->first]-=it->second;
+        const_comp-=e.const_comp;
+        return *this;
+      }
+      ///\e
+      Expr &operator*=(const Value &c) {
+        for (std::map<int, Value>::iterator it=comps.begin(); 
+             it!=comps.end(); ++it)
+          it->second*=c;
+        const_comp*=c;
+        return *this;
+      }
+      ///\e
+      Expr &operator/=(const Value &c) {
+        for (std::map<int, Value>::iterator it=comps.begin(); 
+             it!=comps.end(); ++it)
+          it->second/=c;
+        const_comp/=c;
+        return *this;
+      }
+
+      class MapIt {
+      private:
+        
+        std::map<int, Value>::iterator _it, _end;
+        
+      public:
+
+        MapIt(Expr& e) 
+          : _it(e.comps.begin()), _end(e.comps.end()){}
+
+        operator Col() const { 
+          return colFromId(_it->first); 
+        }
+        
+        Value& operator*() { return _it->second; }
+        const Value& operator*() const { return _it->second; }
+
+        MapIt& operator++() { ++_it; return *this; }
+
+        bool operator==(Invalid) const { return _it == _end; }
+        bool operator!=(Invalid) const { return _it != _end; }
+      };
+
+      class ConstMapIt {
+      private:
+        
+        std::map<int, Value>::const_iterator _it, _end;
+        
+      public:
+
+        ConstMapIt(const Expr& e) 
+          : _it(e.comps.begin()), _end(e.comps.end()){}
+
+        operator Col() const { 
+          return colFromId(_it->first); 
+        }
+        
+        const Value& operator*() const { return _it->second; }
+
+        ConstMapIt& operator++() { ++_it; return *this; }
+
+        bool operator==(Invalid) const { return _it == _end; }
+        bool operator!=(Invalid) const { return _it != _end; }
+      };
+
+    };
+    
+    ///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 fail an assertion.
+    class Constr
+    {
+    public:
+      typedef SolverBase::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) {}
+      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 _lb != -INF && !std::isnan(_lb);
+      }
+      ///Is the constraint upper bounded?
+      bool upperBounded() const {
+        return _ub != INF && !std::isnan(_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.
+    ///\code
+    ///e[v]=5;
+    ///e[v]+=12;
+    ///e.erase(v);
+    ///\endcode
+    ///or you can also iterate through its elements.
+    ///\code
+    ///double s=0;
+    ///for(SolverBase::DualExpr::ConstMapIt i(e);i!=INVALID;++i)
+    ///  s+=*i;
+    ///\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 {
+      friend class SolverBase;
+    public:
+      typedef SolverBase::Row Key; 
+      typedef SolverBase::Value Value;
+      
+    protected:
+      std::map<int, Value> comps;
+      
+    public:
+      typedef True SolverExpr;
+      ///\e
+      DualExpr() {}
+      ///\e
+      DualExpr(const Key &k) {
+        typedef std::map<int, Value>::value_type pair_type;
+        comps.insert(pair_type(id(k), 1));
+      }
+      ///\e
+      Value operator[](const Key& k) const {
+        std::map<int, Value>::const_iterator it = comps.find(id(k));
+        if (it != comps.end()) {
+          return it->second;
+        } else {
+          return 0;
+        }
+      }
+      ///\e
+      Value& operator[](const Key& k) {
+        return comps[id(k)];
+      }
+      ///\e
+      void set(const Key &k, const Value &v) {
+        if (v != 0.0) {
+          typedef std::map<int, Value>::value_type pair_type;
+          comps.insert(pair_type(id(k), v));
+        } else {
+          comps.erase(id(k));
+        }
+      }
+      
+      ///Removes the coefficients closer to zero than \c tolerance.
+      void simplify(Value tolerance = 0.0) {
+        std::map<int, Value>::iterator it=comps.begin();
+        while (it != comps.end()) { 
+          std::map<int, Value>::iterator jt=it;
+          ++jt;
+          if (std::fabs((*it).second) <= tolerance) comps.erase(it);
+          it=jt;
+        }
+      }
+
+      void simplify(Value tolerance = 0.0) const {
+        const_cast<DualExpr*>(this)->simplify(tolerance);
+      }
+
+      ///Sets all coefficients to 0.
+      void clear() {
+        comps.clear();
+      }
+
+      ///\e
+      DualExpr &operator+=(const DualExpr &e) {
+        for (std::map<int, Value>::const_iterator it=e.comps.begin(); 
+             it!=e.comps.end(); ++it)
+          comps[it->first]+=it->second;
+        return *this;
+      }
+      ///\e
+      DualExpr &operator-=(const DualExpr &e) {
+        for (std::map<int, Value>::const_iterator it=e.comps.begin(); 
+             it!=e.comps.end(); ++it)
+          comps[it->first]-=it->second;
+        return *this;
+      }
+      ///\e
+      DualExpr &operator*=(const Value &c) {
+        for (std::map<int, Value>::iterator it=comps.begin(); 
+             it!=comps.end(); ++it)
+          it->second*=c;
+        return *this;
+      }
+      ///\e
+      DualExpr &operator/=(const Value &c) {
+        for (std::map<int, Value>::iterator it=comps.begin(); 
+             it!=comps.end(); ++it)
+          it->second/=c;
+        return *this;
+      }
+
+      class MapIt {
+      private:
+        
+        std::map<int, Value>::iterator _it, _end;
+        
+      public:
+
+        MapIt(DualExpr& e) 
+          : _it(e.comps.begin()), _end(e.comps.end()){}
+
+        operator Row() const { 
+          return rowFromId(_it->first); 
+        }
+        
+        Value& operator*() { return _it->second; }
+        const Value& operator*() const { return _it->second; }
+
+        MapIt& operator++() { ++_it; return *this; }
+
+        bool operator==(Invalid) const { return _it == _end; }
+        bool operator!=(Invalid) const { return _it != _end; }
+      };
+
+      class ConstMapIt {
+      private:
+        
+        std::map<int, Value>::const_iterator _it, _end;
+        
+      public:
+
+        ConstMapIt(const DualExpr& e) 
+          : _it(e.comps.begin()), _end(e.comps.end()){}
+
+        operator Row() const { 
+          return rowFromId(_it->first); 
+        }
+        
+        const Value& operator*() const { return _it->second; }
+
+        ConstMapIt& operator++() { ++_it; return *this; }
+
+        bool operator==(Invalid) const { return _it == _end; }
+        bool operator!=(Invalid) const { return _it != _end; }
+      };
+    };
+    
+
+  protected:
+
+    class InsertIterator {
+    private:
+      
+      std::map<int, Value>& _host;
+      const _solver_bits::VarIndex& _index;
+
+    public:
+
+      typedef std::output_iterator_tag iterator_category;
+      typedef void difference_type;
+      typedef void value_type;
+      typedef void reference;
+      typedef void pointer;
+      
+      InsertIterator(std::map<int, Value>& host, 
+                   const _solver_bits::VarIndex& index) 
+        : _host(host), _index(index) {}
+
+      InsertIterator& operator=(const std::pair<int, Value>& value) {
+        typedef std::map<int, Value>::value_type pair_type;
+        _host.insert(pair_type(_index[value.first], value.second));
+        return *this;
+      }
+
+      InsertIterator& operator*() { return *this; }
+      InsertIterator& operator++() { return *this; }
+      InsertIterator operator++(int) { return *this; }
+
+    };
+
+    class ExprIterator {
+    private:
+      std::map<int, Value>::const_iterator _host_it;
+      const _solver_bits::VarIndex& _index;
+    public:
+
+      typedef std::bidirectional_iterator_tag iterator_category;
+      typedef std::ptrdiff_t difference_type;
+      typedef const std::pair<int, Value> value_type;
+      typedef value_type reference;
+
+      class pointer {
+      public:
+        pointer(value_type& _value) : value(_value) {}
+        value_type* operator->() { return &value; }
+      private:
+        value_type value;
+      };
+
+      ExprIterator(const std::map<int, Value>::const_iterator& host_it,
+                   const _solver_bits::VarIndex& index) 
+        : _host_it(host_it), _index(index) {}
+
+      reference operator*() {
+        return std::make_pair(_index(_host_it->first), _host_it->second);
+      }
+
+      pointer operator->() {
+        return pointer(operator*());
+      }
+
+      ExprIterator& operator++() { ++_host_it; return *this; }
+      ExprIterator operator++(int) {
+        ExprIterator tmp(*this); ++_host_it; return tmp;
+      }
+
+      ExprIterator& operator--() { --_host_it; return *this; }
+      ExprIterator operator--(int) {
+        ExprIterator tmp(*this); --_host_it; return tmp;
+      }
+
+      bool operator==(const ExprIterator& it) const {
+        return _host_it == it._host_it;
+      }
+
+      bool operator!=(const ExprIterator& it) const {
+        return _host_it != it._host_it;
+      }
+
+    };
+
+  protected:
+
+    //Abstract virtual functions
+    virtual SolverBase* _newSolver() const = 0;
+    virtual SolverBase* _cloneSolver() const = 0;
+
+    virtual int _addColId(int col) { return cols.addIndex(col); }
+    virtual int _addRowId(int row) { return rows.addIndex(row); }
+
+    virtual void _eraseColId(int col) { cols.eraseIndex(col); }
+    virtual void _eraseRowId(int row) { rows.eraseIndex(row); }
+
+    virtual int _addCol() = 0;
+    virtual 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 _getRowName(int row, std::string& name) const = 0;
+    virtual void _setRowName(int row, const std::string& name) = 0;
+    virtual int _rowByName(const std::string& name) const = 0;
+
+    virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e) = 0;
+    virtual void _getRowCoeffs(int i, InsertIterator b) const = 0;
+
+    virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e) = 0;
+    virtual void _getColCoeffs(int i, InsertIterator b) const = 0;
+
+    virtual void _setCoeff(int row, int col, Value value) = 0;
+    virtual Value _getCoeff(int row, int col) const = 0;
+
+    virtual void _setColLowerBound(int i, Value value) = 0;
+    virtual Value _getColLowerBound(int i) const = 0;
+
+    virtual void _setColUpperBound(int i, Value value) = 0;
+    virtual Value _getColUpperBound(int i) const = 0;
+
+    virtual void _setRowLowerBound(int i, Value value) = 0;
+    virtual Value _getRowLowerBound(int i) const = 0;
+
+    virtual void _setRowUpperBound(int i, Value value) = 0;
+    virtual Value _getRowUpperBound(int i) const = 0;
+
+    virtual void _setObjCoeffs(ExprIterator b, ExprIterator e) = 0;
+    virtual void _getObjCoeffs(InsertIterator b) const = 0;
+
+    virtual void _setObjCoeff(int i, Value obj_coef) = 0;
+    virtual Value _getObjCoeff(int i) const = 0;
+
+    virtual void _setSense(Sense) = 0;
+    virtual Sense _getSense() const = 0;
+
+    virtual void _clear() = 0;
+
+    //Own protected stuff
+    
+    //Constant component of the objective function
+    Value obj_const_comp;
+        
+  public:
+
+    ///\e
+    SolverBase() : rows(), cols(), obj_const_comp(0) {
+    }
+
+    ///\e
+    virtual ~SolverBase() {}
+
+    ///Creates a new LP problem
+    SolverBase* newSolver() {return _newSolver();}
+    ///Makes a copy of the LP problem
+    SolverBase* cloneSolver() {return _cloneSolver();}
+    
+    ///\name Build up and modify the LP
+
+    ///@{
+
+    ///Add a new empty column (i.e a new variable) to the LP
+    Col addCol() { Col c; c._id = _addColId(_addCol()); 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<SolverBase::Col>
+    ///std::list<SolverBase::Col>
+    ///\endcode
+    ///- a standard STL compatible iterable container with
+    ///\ref Col as its \c mapped_type
+    ///like
+    ///\code
+    ///std::map<AnyType,SolverBase::Col>
+    ///\endcode
+    ///- an iterable lemon \ref concepts::WriteMap "write map" like 
+    ///\code
+    ///ListGraph::NodeMap<SolverBase::Col>
+    ///ListGraph::ArcMap<SolverBase::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(cols(id(c)), ExprIterator(e.comps.begin(), cols), 
+                    ExprIterator(e.comps.end(), cols));
+    }
+
+    ///Get a column (i.e a dual constraint) of the LP
+
+    ///\param r is the column to get
+    ///\return the dual expression associated to the column
+    DualExpr col(Col c) const {
+      DualExpr e;
+      _getColCoeffs(cols(id(c)), InsertIterator(e.comps, rows));
+      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; r._id = _addRowId(_addRow()); 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<SolverBase::Row>
+    ///std::list<SolverBase::Row>
+    ///\endcode
+    ///- a standard STL compatible iterable container with
+    ///\ref Row as its \c mapped_type
+    ///like
+    ///\code
+    ///std::map<AnyType,SolverBase::Row>
+    ///\endcode
+    ///- an iterable lemon \ref concepts::WriteMap "write map" like 
+    ///\code
+    ///ListGraph::NodeMap<SolverBase::Row>
+    ///ListGraph::ArcMap<SolverBase::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)
+    void row(Row r, Value l, const Expr &e, Value u) {
+      e.simplify();
+      _setRowCoeffs(rows(id(r)), ExprIterator(e.comps.begin(), cols),
+                    ExprIterator(e.comps.end(), cols));
+      _setRowLowerBound(rows(id(r)),l - *e);
+      _setRowUpperBound(rows(id(r)),u - *e);
+    }
+
+    ///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(rows(id(r)), InsertIterator(e.comps, cols));
+      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.
+    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 column (i.e a variable) from the LP
+
+    ///\param c is the column to be deleted
+    ///\todo Please check this
+    void erase(Col c) {
+      _eraseCol(cols(id(c)));
+      _eraseColId(cols(id(c)));
+    }
+    ///Erase a row (i.e a constraint) from the LP
+
+    ///\param r is the row to be deleted
+    ///\todo Please check this
+    void erase(Row r) {
+      _eraseRow(rows(id(r)));
+      _eraseRowId(rows(id(r)));
+    }
+
+    /// Get the name of a column
+    
+    ///\param c is the coresponding column 
+    ///\return The name of the colunm
+    std::string colName(Col c) const {
+      std::string name;
+      _getColName(cols(id(c)), name);
+      return name;
+    }
+    
+    /// Set the name of a column
+    
+    ///\param c is the coresponding column 
+    ///\param name The name to be given
+    void colName(Col c, const std::string& name) {
+      _setColName(cols(id(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[k]) : Col(INVALID);
+    }
+
+    /// Get the name of a row
+    
+    ///\param r is the coresponding row 
+    ///\return The name of the row
+    std::string rowName(Row r) const {
+      std::string name;
+      _getRowName(rows(id(r)), name);
+      return name;
+    }
+    
+    /// Set the name of a row
+    
+    ///\param r is the coresponding row 
+    ///\param name The name to be given
+    void rowName(Row r, const std::string& name) {
+      _setRowName(rows(id(r)), name);
+    }
+
+    /// Get the row by its name
+    
+    ///\param name The name of the row
+    ///\return the proper row or \c INVALID
+    Row rowByName(const std::string& name) const {
+      int k = _rowByName(name);
+      return k != -1 ? Row(rows[k]) : Row(INVALID);
+    }
+    
+    /// 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 column of the element to be modified
+    ///\param val is the new value of the coefficient
+    void coeff(Row r, Col c, Value val) {
+      _setCoeff(rows(id(r)),cols(id(c)), val);
+    }
+
+    /// Get an element of the coefficient matrix of the LP
+
+    ///\param r is the row of the element in question
+    ///\param c is the column of the element in question
+    ///\return the corresponding coefficient
+    Value coeff(Row r, Col c) const {
+      return _getCoeff(rows(id(r)),cols(id(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(cols(id(c)),value);
+    }
+
+    /// Get the lower bound of a column (i.e a variable)
+
+    /// This function returns the lower bound for column (variable) \t c
+    /// (this might be -\ref INF as well).  
+    ///\return The lower bound for column \t c
+    Value colLowerBound(Col c) const {
+      return _getColLowerBound(cols(id(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(cols(id(c)),value);
+    };
+
+    /// Get the upper bound of a column (i.e a variable)
+
+    /// This function returns the upper bound for column (variable) \t c
+    /// (this might be \ref INF as well).  
+    ///\return The upper bound for column \t c
+    Value colUpperBound(Col c) const {
+      return _getColUpperBound(cols(id(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(cols(id(c)),lower);
+      _setColUpperBound(cols(id(c)),upper);
+    }
+    
+    ///\brief Set the lower and the upper bound of several columns
+    ///(i.e a variables) at once
+    ///
+    ///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 bound of a row (i.e a constraint)
+
+    /// The lower bound of a constraint (row) has to be given by an 
+    /// extended number of type Value, i.e. a finite number of type 
+    /// Value or -\ref INF.
+    void rowLowerBound(Row r, Value value) {
+      _setRowLowerBound(rows(id(r)),value);
+    }
+
+    /// Get the lower bound of a row (i.e a constraint)
+
+    /// This function returns the lower bound for row (constraint) \t c
+    /// (this might be -\ref INF as well).  
+    ///\return The lower bound for row \t r
+    Value rowLowerBound(Row r) const {
+      return _getRowLowerBound(rows(id(r)));
+    }
+
+    /// Set the upper bound of a row (i.e a constraint)
+
+    /// The upper bound of a constraint (row) has to be given by an 
+    /// extended number of type Value, i.e. a finite number of type 
+    /// Value or -\ref INF.
+    void rowUpperBound(Row r, Value value) {
+      _setRowUpperBound(rows(id(r)),value);
+    }
+
+    /// Get the upper bound of a row (i.e a constraint)
+
+    /// This function returns the upper bound for row (constraint) \t c
+    /// (this might be -\ref INF as well).  
+    ///\return The upper bound for row \t r
+    Value rowUpperBound(Row r) const {
+      return _getRowUpperBound(rows(id(r)));
+    }
+
+    ///Set an element of the objective function
+    void objCoeff(Col c, Value v) {_setObjCoeff(cols(id(c)),v); };
+
+    ///Get an element of the objective function
+    Value objCoeff(Col c) const { return _getObjCoeff(cols(id(c))); };
+
+    ///Set the objective function
+
+    ///\param e is a linear expression of type \ref Expr.
+    void obj(const Expr& e) {
+      _setObjCoeffs(ExprIterator(e.comps.begin(), cols), 
+                    ExprIterator(e.comps.end(), cols));
+      obj_const_comp = *e;
+    }
+
+    ///Get the objective function
+
+    ///\return the objective function as a linear expression of type \ref Expr.
+    Expr obj() const {
+      Expr e;
+      _getObjCoeffs(InsertIterator(e.comps, cols));
+      *e = obj_const_comp;
+      return e;
+    }
+    
+
+    ///Set the direction of optimization
+    void sense(Sense sense) { _setSense(sense); }
+
+    ///Query the direction of the optimization
+    Sense sense() const {return _getSense(); }
+
+    ///Set the sense to maximization
+    void max() { _setSense(MAX); }
+
+    ///Set the sense to maximization
+    void min() { _setSense(MIN); }
+
+    ///Clears the problem
+    void clear() { _clear(); }
+    
+    ///@}
+
+  };
+
+  ///\relates SolverBase::Expr
+  ///
+  inline SolverBase::Expr operator+(const SolverBase::Expr &a,
+                                    const SolverBase::Expr &b) {
+    SolverBase::Expr tmp(a);
+    tmp+=b;
+    return tmp;
+  }
+  ///\e
+  
+  ///\relates SolverBase::Expr
+  ///
+  inline SolverBase::Expr operator-(const SolverBase::Expr &a,
+                                    const SolverBase::Expr &b) {
+    SolverBase::Expr tmp(a);
+    tmp-=b;
+    return tmp;
+  }
+  ///\e
+  
+  ///\relates SolverBase::Expr
+  ///
+  inline SolverBase::Expr operator*(const SolverBase::Expr &a,
+                                    const SolverBase::Value &b) {
+    SolverBase::Expr tmp(a);
+    tmp*=b;
+    return tmp;
+  }
+  
+  ///\e
+  
+  ///\relates SolverBase::Expr
+  ///
+  inline SolverBase::Expr operator*(const SolverBase::Value &a,
+                                    const SolverBase::Expr &b) {
+    SolverBase::Expr tmp(b);
+    tmp*=a;
+    return tmp;
+  }
+  ///\e
+  
+  ///\relates SolverBase::Expr
+  ///
+  inline SolverBase::Expr operator/(const SolverBase::Expr &a,
+                                    const SolverBase::Value &b) {
+    SolverBase::Expr tmp(a);
+    tmp/=b;
+    return tmp;
+  }
+  
+  ///\e
+  
+  ///\relates SolverBase::Constr
+  ///
+  inline SolverBase::Constr operator<=(const SolverBase::Expr &e,
+                                       const SolverBase::Expr &f) {
+    return SolverBase::Constr(0, f - e, SolverBase::INF);
+  }
+
+  ///\e
+  
+  ///\relates SolverBase::Constr
+  ///
+  inline SolverBase::Constr operator<=(const SolverBase::Value &e,
+                                       const SolverBase::Expr &f) {
+    return SolverBase::Constr(e, f, SolverBase::NaN);
+  }
+
+  ///\e
+  
+  ///\relates SolverBase::Constr
+  ///
+  inline SolverBase::Constr operator<=(const SolverBase::Expr &e,
+                                       const SolverBase::Value &f) {
+    return SolverBase::Constr(- SolverBase::INF, e, f);
+  }
+
+  ///\e
+  
+  ///\relates SolverBase::Constr
+  ///
+  inline SolverBase::Constr operator>=(const SolverBase::Expr &e,
+                                       const SolverBase::Expr &f) {
+    return SolverBase::Constr(0, e - f, SolverBase::INF);
+  }
+
+
+  ///\e
+  
+  ///\relates SolverBase::Constr
+  ///
+  inline SolverBase::Constr operator>=(const SolverBase::Value &e,
+                                       const SolverBase::Expr &f) {
+    return SolverBase::Constr(SolverBase::NaN, f, e);
+  }
+
+
+  ///\e
+  
+  ///\relates SolverBase::Constr
+  ///
+  inline SolverBase::Constr operator>=(const SolverBase::Expr &e,
+                                       const SolverBase::Value &f) {
+    return SolverBase::Constr(f, e, SolverBase::INF);
+  }
+
+  ///\e
+
+  ///\relates SolverBase::Constr
+  ///
+  inline SolverBase::Constr operator==(const SolverBase::Expr &e,
+                                       const SolverBase::Value &f) {
+    return SolverBase::Constr(f, e, f);
+  }
+
+  ///\e
+  
+  ///\relates SolverBase::Constr
+  ///
+  inline SolverBase::Constr operator==(const SolverBase::Expr &e,
+                                       const SolverBase::Expr &f) {
+    return SolverBase::Constr(0, f - e, 0);
+  }
+
+  ///\e
+  
+  ///\relates SolverBase::Constr
+  ///
+  inline SolverBase::Constr operator<=(const SolverBase::Value &n,
+                                       const SolverBase::Constr &c) {
+    SolverBase::Constr tmp(c);
+    LEMON_ASSERT(std::isnan(tmp.lowerBound()), "Wrong LP constraint");
+    tmp.lowerBound()=n;
+    return tmp;
+  }
+  ///\e
+  
+  ///\relates SolverBase::Constr
+  ///
+  inline SolverBase::Constr operator<=(const SolverBase::Constr &c,
+                                       const SolverBase::Value &n)
+  {
+    SolverBase::Constr tmp(c);
+    LEMON_ASSERT(std::isnan(tmp.upperBound()), "Wrong LP constraint");
+    tmp.upperBound()=n;
+    return tmp;
+  }
+
+  ///\e
+  
+  ///\relates SolverBase::Constr
+  ///
+  inline SolverBase::Constr operator>=(const SolverBase::Value &n,
+                                       const SolverBase::Constr &c) {
+    SolverBase::Constr tmp(c);
+    LEMON_ASSERT(std::isnan(tmp.upperBound()), "Wrong LP constraint");
+    tmp.upperBound()=n;
+    return tmp;
+  }
+  ///\e
+  
+  ///\relates SolverBase::Constr
+  ///
+  inline SolverBase::Constr operator>=(const SolverBase::Constr &c,
+                                       const SolverBase::Value &n)
+  {
+    SolverBase::Constr tmp(c);
+    LEMON_ASSERT(std::isnan(tmp.lowerBound()), "Wrong LP constraint");
+    tmp.lowerBound()=n;
+    return tmp;
+  }
+
+  ///\e
+  
+  ///\relates SolverBase::DualExpr
+  ///
+  inline SolverBase::DualExpr operator+(const SolverBase::DualExpr &a,
+                                        const SolverBase::DualExpr &b) {
+    SolverBase::DualExpr tmp(a);
+    tmp+=b;
+    return tmp;
+  }
+  ///\e
+  
+  ///\relates SolverBase::DualExpr
+  ///
+  inline SolverBase::DualExpr operator-(const SolverBase::DualExpr &a,
+                                        const SolverBase::DualExpr &b) {
+    SolverBase::DualExpr tmp(a);
+    tmp-=b;
+    return tmp;
+  }
+  ///\e
+  
+  ///\relates SolverBase::DualExpr
+  ///
+  inline SolverBase::DualExpr operator*(const SolverBase::DualExpr &a,
+                                        const SolverBase::Value &b) {
+    SolverBase::DualExpr tmp(a);
+    tmp*=b;
+    return tmp;
+  }
+  
+  ///\e
+  
+  ///\relates SolverBase::DualExpr
+  ///
+  inline SolverBase::DualExpr operator*(const SolverBase::Value &a,
+                                        const SolverBase::DualExpr &b) {
+    SolverBase::DualExpr tmp(b);
+    tmp*=a;
+    return tmp;
+  }
+  ///\e
+  
+  ///\relates SolverBase::DualExpr
+  ///
+  inline SolverBase::DualExpr operator/(const SolverBase::DualExpr &a,
+                                        const SolverBase::Value &b) {
+    SolverBase::DualExpr tmp(a);
+    tmp/=b;
+    return tmp;
+  }
+
+  class LpSolverBase : virtual public SolverBase {
+  public:
+
+    ///\e
+    enum ProblemType {
+      ///Feasible solution hasn't been found (but may exist).
+      UNDEFINED = 0,
+      ///The problem has no feasible solution
+      INFEASIBLE = 1,
+      ///Feasible solution found
+      FEASIBLE = 2,
+      ///Optimal solution exists and found
+      OPTIMAL = 3,
+      ///The cost function is unbounded
+      UNBOUNDED = 4
+    };
+
+    enum VarStatus {
+      BASIC, FREE, LOWER, UPPER, FIXED
+    };
+
+  protected:
+
+    virtual SolveExitStatus _solve() = 0;
+
+    virtual Value _getPrimal(int i) const = 0;
+    virtual Value _getDual(int i) const = 0;
+
+    virtual Value _getPrimalRay(int i) const = 0;
+    virtual Value _getDualRay(int i) const = 0;
+
+    virtual Value _getPrimalValue() const = 0;
+
+    virtual VarStatus _getColStatus(int i) const = 0;
+    virtual VarStatus _getRowStatus(int i) const = 0;
+
+    virtual ProblemType _getPrimalType() const = 0;
+    virtual ProblemType _getDualType() const = 0;
+
+  public:
+
+    ///\name Solve the LP
+
+    ///@{
+
+    ///\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 type of the primal problem
+    ProblemType primalType() const {
+      return _getPrimalType();
+    }
+
+    /// The type of the dual problem 
+    ProblemType dualType() const {
+      return _getDualType();
+    }
+
+    ///\e
+    Value primal(Col c) const { return _getPrimal(cols(id(c))); }
+    ///\e
+    Value primal(const Expr& e) const {
+      double res = *e;
+      for (Expr::ConstMapIt c(e); c != INVALID; ++c) {
+        res += *c * primal(c);
+      }
+      return res; 
+    }
+    ///\e
+    Value primalRay(Col c) const { return _getPrimalRay(cols(id(c))); }
+
+    ///\e
+    Value dual(Row r) const { return _getDual(rows(id(r))); }
+    ///\e
+    Value dual(const DualExpr& e) const {
+      double res = 0.0;
+      for (DualExpr::ConstMapIt r(e); r != INVALID; ++r) {
+        res += *r * dual(r);
+      }
+      return res; 
+    }
+    ///\e
+    Value dualRay(Row r) const { return _getDualRay(rows(id(r))); }
+
+    ///\e
+    VarStatus colStatus(Col c) const { return _getColStatus(cols(id(c))); }
+
+    ///\e
+    VarStatus rowStatus(Col c) const { return _getRowStatus(cols(id(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 primal() const { return _getPrimalValue()+obj_const_comp;}
+    ///@}
+    
+    LpSolverBase* newSolver() {return _newSolver();}
+    LpSolverBase* cloneSolver() {return _cloneSolver();}
+
+  protected:
+
+    virtual LpSolverBase* _newSolver() const = 0;
+    virtual LpSolverBase* _cloneSolver() const = 0;
+  };  
+
+
+  /// \ingroup lp_group
+  ///
+  /// \brief Common base class for MIP solvers
+  /// \todo Much more docs
+  class MipSolverBase : virtual public SolverBase {
+  public:
+
+    ///\e
+    enum ProblemType {
+      ///Feasible solution hasn't been found (but may exist).
+      UNDEFINED = 0,
+      ///The problem has no feasible solution
+      INFEASIBLE = 1,
+      ///Feasible solution found
+      FEASIBLE = 2,
+      ///Optimal solution exists and found
+      OPTIMAL = 3,
+      ///The cost function is unbounded
+      ///
+      ///The Mip or at least the relaxed problem is unbounded
+      UNBOUNDED = 4
+    };
+
+    ///\name Solve the MIP
+
+    ///@{
+
+    ///\e Solve the MIP problem at hand
+    ///
+    ///\return The result of the optimization procedure. Possible 
+    ///values and their meanings can be found in the documentation of 
+    ///\ref SolveExitStatus.
+    ///
+    ///\todo Which method is used to solve the problem
+    SolveExitStatus solve() { return _solve(); }
+    
+    ///@}
+
+    ///\name Setting column type
+    ///@{
+
+    ///Possible variable (column) types (e.g. real, integer, binary etc.)
+    enum ColTypes {
+      ///Continuous variable (default)
+      REAL = 0,
+      ///Integer variable
+      INTEGER = 1
+    };
+
+    ///Sets the type of the given column to the given type
+    ///
+    ///Sets the type of the given column to the given type.
+    void colType(Col c, ColTypes col_type) {
+      _setColType(cols(id(c)),col_type);
+    }
+
+    ///Gives back the type of the column.
+    ///
+    ///Gives back the type of the column.
+    ColTypes colType(Col c) const {
+      return _getColType(cols(id(c)));
+    }
+
+    ///Sets the type of the given Col to integer.
+    ///
+    ///Sets the type of the given Col to integer.
+    void integer(Col c) {
+      colType(c, INTEGER);
+    }
+
+    ///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) == INTEGER);
+    }
+
+    ///Sets the type of the given Col to continuous.
+    ///
+    ///Sets the type of the given Col to continuous.
+    void real(Col c) {
+      colType(c, REAL);
+    }
+
+    ///Gives back whether the type of the column is continuous or not.
+    ///
+    ///Gives back the type of the column.
+    ///\return true if the column has continuous type and false if not.
+    bool real(Col c) const {
+      return (colType(c) == REAL);
+    }
+    ///@}
+
+    ///\name Obtain the solution
+
+    ///@{
+
+    /// The type of the MIP problem
+    ProblemType type() const {
+      return _getType();
+    }
+
+    ///\e
+    Value sol(Col c) const { return _getSol(cols(id(c))); }
+    ///\e
+    Value sol(const Expr& e) const {
+      double res = *e;
+      for (Expr::ConstMapIt c(e); c != INVALID; ++c) {
+        res += *c * sol(c);
+      }
+      return res; 
+    }
+    ///\return
+    ///- \ref INF or -\ref INF means either infeasibility or unboundedness
+    /// of the problem, depending on whether we minimize or maximize.
+    ///- \ref NaN if no primal solution is found.
+    ///- The (finite) objective value if an optimal solution is found.
+    Value solValue() const { return _getSolValue()+obj_const_comp;}
+    ///@}
+
+  protected:
+
+    virtual SolveExitStatus _solve() = 0;
+    virtual ColTypes _getColType(int col) const = 0;
+    virtual void _setColType(int col, ColTypes col_type) = 0;
+    virtual ProblemType _getType() const = 0;
+    virtual Value _getSol(int i) const = 0;
+    virtual Value _getSolValue() const = 0;
+
+  public:
+
+    MipSolverBase* newSolver() {return _newSolver();}
+    MipSolverBase* cloneSolver() {return _cloneSolver();}
+
+  protected:
+
+    virtual MipSolverBase* _newSolver() const = 0;
+    virtual MipSolverBase* _cloneSolver() const = 0;
+  };
+  
+  
+
+} //namespace lemon
+
+#endif //LEMON_LP_BASE_H

new file lemon/lp_clp.cc

@@ +1 @@
+/* -*- 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_clp.h>
+#include <coin/ClpSimplex.hpp>
+
+namespace lemon {
+
+  LpClp::LpClp() {
+    _prob = new ClpSimplex();
+    _init_temporals();
+    messageLevel(MESSAGE_NO_OUTPUT);
+  }
+
+  LpClp::LpClp(const LpClp& other) {
+    _prob = new ClpSimplex(*other._prob);
+    rows = other.rows;
+    cols = other.cols;
+    _init_temporals();
+    messageLevel(MESSAGE_NO_OUTPUT);
+  }
+
+  LpClp::~LpClp() {
+    delete _prob;
+    _clear_temporals();
+  }
+
+  void LpClp::_init_temporals() {
+    _primal_ray = 0;
+    _dual_ray = 0;
+  }
+
+  void LpClp::_clear_temporals() {
+    if (_primal_ray) {
+      delete[] _primal_ray;
+      _primal_ray = 0;
+    }
+    if (_dual_ray) {
+      delete[] _dual_ray;
+      _dual_ray = 0;
+    }
+  }
+
+  LpClp* LpClp::_newSolver() const {
+    LpClp* newlp = new LpClp;
+    return newlp;
+  }
+
+  LpClp* LpClp::_cloneSolver() const {
+    LpClp* copylp = new LpClp(*this);
+    return copylp;
+  }
+
+  int LpClp::_addCol() {
+    _prob->addColumn(0, 0, 0, -COIN_DBL_MAX, COIN_DBL_MAX, 0.0);
+    return _prob->numberColumns() - 1;
+  }
+
+  int LpClp::_addRow() {
+    _prob->addRow(0, 0, 0, -COIN_DBL_MAX, COIN_DBL_MAX);
+    return _prob->numberRows() - 1;
+  }
+
+
+  void LpClp::_eraseCol(int c) {
+    _col_names_ref.erase(_prob->getColumnName(c));
+    _prob->deleteColumns(1, &c);
+  }
+
+  void LpClp::_eraseRow(int r) {
+    _row_names_ref.erase(_prob->getRowName(r));
+    _prob->deleteRows(1, &r);
+  }
+
+  void LpClp::_eraseColId(int i) {
+    cols.eraseIndex(i);
+    cols.shiftIndexes(i);
+  }
+
+  void LpClp::_eraseRowId(int i) {
+    rows.eraseIndex(i);
+    rows.shiftIndexes(i);
+  }
+
+  void LpClp::_getColName(int c, std::string& name) const {
+    name = _prob->getColumnName(c);
+  }
+
+  void LpClp::_setColName(int c, const std::string& name) {
+    _prob->setColumnName(c, const_cast<std::string&>(name));
+    _col_names_ref[name] = c;
+  }
+
+  int LpClp::_colByName(const std::string& name) const {
+    std::map<std::string, int>::const_iterator it = _col_names_ref.find(name);
+    return it != _col_names_ref.end() ? it->second : -1;
+  }
+
+  void LpClp::_getRowName(int r, std::string& name) const {
+    name = _prob->getRowName(r);
+  }
+
+  void LpClp::_setRowName(int r, const std::string& name) {
+    _prob->setRowName(r, const_cast<std::string&>(name));
+    _row_names_ref[name] = r;
+  }
+
+  int LpClp::_rowByName(const std::string& name) const {
+    std::map<std::string, int>::const_iterator it = _row_names_ref.find(name);
+    return it != _row_names_ref.end() ? it->second : -1;
+  }
+
+
+  void LpClp::_setRowCoeffs(int ix, ExprIterator b, ExprIterator e) {
+    std::map<int, Value> coeffs;
+
+    int n = _prob->clpMatrix()->getNumCols();
+
+    const int* indices = _prob->clpMatrix()->getIndices();
+    const double* elements = _prob->clpMatrix()->getElements();
+
+    for (int i = 0; i < n; ++i) {
+      CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[i];
+      CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[i];
+
+      const int* it = std::lower_bound(indices + begin, indices + end, ix);
+      if (it != indices + end && *it == ix && elements[it - indices] != 0.0) {
+        coeffs[i] = 0.0;
+      }
+    }
+    
+    for (ExprIterator it = b; it != e; ++it) {
+      coeffs[it->first] = it->second;
+    }
+    
+    for (std::map<int, Value>::iterator it = coeffs.begin(); 
+         it != coeffs.end(); ++it) {
+      _prob->modifyCoefficient(ix, it->first, it->second);
+    }   
+  }
+
+  void LpClp::_getRowCoeffs(int ix, InsertIterator b) const {
+    int n = _prob->clpMatrix()->getNumCols();
+
+    const int* indices = _prob->clpMatrix()->getIndices();
+    const double* elements = _prob->clpMatrix()->getElements();
+
+    for (int i = 0; i < n; ++i) {
+      CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[i];
+      CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[i];
+
+      const int* it = std::lower_bound(indices + begin, indices + end, ix);
+      if (it != indices + end && *it == ix) {
+        *b = std::make_pair(i, elements[it - indices]);
+      }
+    }
+  }
+
+  void LpClp::_setColCoeffs(int ix, ExprIterator b, ExprIterator e) {
+    std::map<int, Value> coeffs;
+
+    CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[ix];
+    CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[ix];
+    
+    const int* indices = _prob->clpMatrix()->getIndices();
+    const double* elements = _prob->clpMatrix()->getElements();
+
+    for (CoinBigIndex i = begin; i != end; ++i) {
+      if (elements[i] != 0.0) {
+        coeffs[indices[i]] = 0.0;
+      }
+    }
+    for (ExprIterator it = b; it != e; ++it) {      
+      coeffs[it->first] = it->second;
+    }
+    for (std::map<int, Value>::iterator it = coeffs.begin(); 
+         it != coeffs.end(); ++it) {
+      _prob->modifyCoefficient(it->first, ix, it->second);
+    }
+  }
+
+  void LpClp::_getColCoeffs(int ix, InsertIterator b) const {
+    CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[ix];
+    CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[ix];
+    
+    const int* indices = _prob->clpMatrix()->getIndices();
+    const double* elements = _prob->clpMatrix()->getElements();
+
+    for (CoinBigIndex i = begin; i != end; ++i) {
+      *b = std::make_pair(indices[i], elements[i]);
+      ++b;
+    }
+  }
+
+  void LpClp::_setCoeff(int ix, int jx, Value value) {
+    _prob->modifyCoefficient(ix, jx, value);
+  }
+
+  LpClp::Value LpClp::_getCoeff(int ix, int jx) const {
+    CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[ix];
+    CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[ix];
+    
+    const int* indices = _prob->clpMatrix()->getIndices();
+    const double* elements = _prob->clpMatrix()->getElements();
+
+    const int* it = std::lower_bound(indices + begin, indices + end, jx);
+    if (it != indices + end && *it == jx) {
+      return elements[it - indices];
+    } else {
+      return 0.0;
+    }
+  }
+
+  void LpClp::_setColLowerBound(int i, Value lo) {
+    _prob->setColumnLower(i, lo == - INF ? - COIN_DBL_MAX : lo);
+  }
+
+  LpClp::Value LpClp::_getColLowerBound(int i) const {
+    double val = _prob->getColLower()[i];
+    return val == - COIN_DBL_MAX ? - INF : val;
+  }
+
+  void LpClp::_setColUpperBound(int i, Value up) {
+    _prob->setColumnUpper(i, up == INF ? COIN_DBL_MAX : up);
+  }
+
+  LpClp::Value LpClp::_getColUpperBound(int i) const {
+    double val = _prob->getColUpper()[i];
+    return val == COIN_DBL_MAX ? INF : val;
+  }
+
+  void LpClp::_setRowLowerBound(int i, Value lo) {
+    _prob->setRowLower(i, lo == - INF ? - COIN_DBL_MAX : lo);
+  }
+
+  LpClp::Value LpClp::_getRowLowerBound(int i) const {
+    double val = _prob->getRowLower()[i];
+    return val == - COIN_DBL_MAX ? - INF : val;
+  }
+
+  void LpClp::_setRowUpperBound(int i, Value up) {
+    _prob->setRowUpper(i, up == INF ? COIN_DBL_MAX : up);
+  }
+
+  LpClp::Value LpClp::_getRowUpperBound(int i) const {
+    double val = _prob->getRowUpper()[i];
+    return val == COIN_DBL_MAX ? INF : val;
+  }
+
+  void LpClp::_setObjCoeffs(ExprIterator b, ExprIterator e) {
+    int num = _prob->clpMatrix()->getNumCols();
+    for (int i = 0; i < num; ++i) {
+      _prob->setObjectiveCoefficient(i, 0.0);
+    }
+    for (ExprIterator it = b; it != e; ++it) {
+      _prob->setObjectiveCoefficient(it->first, it->second);
+    }
+  }
+
+  void LpClp::_getObjCoeffs(InsertIterator b) const {
+    int num = _prob->clpMatrix()->getNumCols();
+    for (int i = 0; i < num; ++i) {
+      Value coef = _prob->getObjCoefficients()[i];
+      if (coef != 0.0) {
+        *b = std::make_pair(i, coef);
+        ++b;
+      }
+    }
+  }
+
+  void LpClp::_setObjCoeff(int i, Value obj_coef) {
+    _prob->setObjectiveCoefficient(i, obj_coef);
+  }
+
+  LpClp::Value LpClp::_getObjCoeff(int i) const {
+    return _prob->getObjCoefficients()[i];
+  }
+
+  LpClp::SolveExitStatus LpClp::_solve() {
+    return _prob->primal() >= 0 ? SOLVED : UNSOLVED;
+  }
+ 
+  LpClp::SolveExitStatus LpClp::solvePrimal() {
+    return _prob->primal() >= 0 ? SOLVED : UNSOLVED;
+  }
+
+  LpClp::SolveExitStatus LpClp::solveDual() {
+    return _prob->dual() >= 0 ? SOLVED : UNSOLVED;
+  }
+
+  LpClp::SolveExitStatus LpClp::solveBarrier() {
+    return _prob->barrier() >= 0 ? SOLVED : UNSOLVED;
+  }
+
+  LpClp::Value LpClp::_getPrimal(int i) const {
+    return _prob->primalColumnSolution()[i];
+  }
+  LpClp::Value LpClp::_getPrimalValue() const {
+    return _prob->rawObjectiveValue();
+  }
+
+  LpClp::Value LpClp::_getDual(int i) const {
+    return _prob->dualRowSolution()[i];
+  }
+
+  LpClp::Value LpClp::_getPrimalRay(int i) const {
+    if (!_primal_ray) {
+      _primal_ray = _prob->unboundedRay();
+      LEMON_ASSERT(_primal_ray != 0, "Primal ray is not provided");
+    }
+    return _primal_ray[i];
+  }
+
+  LpClp::Value LpClp::_getDualRay(int i) const {
+    if (!_dual_ray) {
+      _dual_ray = _prob->infeasibilityRay();
+      LEMON_ASSERT(_dual_ray != 0, "Dual ray is not provided");
+    }
+    return _dual_ray[i]; 
+  }
+
+  LpClp::VarStatus LpClp::_getColStatus(int i) const {
+    switch (_prob->getColumnStatus(i)) {
+    case ClpSimplex::basic:
+      return BASIC;
+    case ClpSimplex::isFree:
+      return FREE;
+    case ClpSimplex::atUpperBound:
+      return UPPER; 
+    case ClpSimplex::atLowerBound:
+      return LOWER;
+    case ClpSimplex::isFixed:
+      return FIXED;
+    case ClpSimplex::superBasic:
+      return FREE;
+    default:
+      LEMON_ASSERT(false, "Wrong column status");
+      return VarStatus();
+    }
+  }
+
+  LpClp::VarStatus LpClp::_getRowStatus(int i) const {
+    switch (_prob->getColumnStatus(i)) {
+    case ClpSimplex::basic:
+      return BASIC;
+    case ClpSimplex::isFree:
+      return FREE;
+    case ClpSimplex::atUpperBound:
+      return UPPER; 
+    case ClpSimplex::atLowerBound:
+      return LOWER;     
+    case ClpSimplex::isFixed:
+      return FIXED;
+    case ClpSimplex::superBasic:
+      return FREE;
+    default:
+      LEMON_ASSERT(false, "Wrong row status");
+      return VarStatus();
+    }
+  }
+
+
+  LpClp::ProblemType LpClp::_getPrimalType() const {
+    if (_prob->isProvenOptimal()) {
+      return OPTIMAL;
+    } else if (_prob->isProvenPrimalInfeasible()) {
+      return INFEASIBLE;
+    } else if (_prob->isProvenDualInfeasible()) {
+      return UNBOUNDED;
+    } else {
+      return UNDEFINED;
+    }
+  }
+
+  LpClp::ProblemType LpClp::_getDualType() const {
+    if (_prob->isProvenOptimal()) {
+      return OPTIMAL;
+    } else if (_prob->isProvenDualInfeasible()) {
+      return INFEASIBLE;      
+    } else if (_prob->isProvenPrimalInfeasible()) {
+      return INFEASIBLE;
+    } else {
+      return UNDEFINED;
+    }
+  }
+
+  void LpClp::_setSense(SolverBase::Sense sense) {
+    switch (sense) {
+    case MIN:
+      _prob->setOptimizationDirection(1);
+      break;
+    case MAX:
+      _prob->setOptimizationDirection(-1);
+      break;
+    }
+  }
+
+  SolverBase::Sense LpClp::_getSense() const {
+    double dir = _prob->optimizationDirection();
+    if (dir > 0.0) {
+      return MIN;
+    } else {
+      return MAX;
+    }
+  }
+
+  void LpClp::_clear() {
+    delete _prob;
+    _prob = new ClpSimplex();
+    rows.clear();
+    cols.clear();
+    _col_names_ref.clear();
+    _clear_temporals();
+  }
+
+  void LpClp::messageLevel(MessageLevel m) {
+    _prob->setLogLevel(static_cast<int>(m));
+  }
+
+} //END OF NAMESPACE LEMON

new file lemon/lp_clp.h

@@ +1 @@
+/* -*- C++ -*-
+ *
+ * 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_CLP_H
+#define LEMON_LP_CLP_H
+
+///\file
+///\brief Header of the LEMON-CLP lp solver interface.
+
+#include <vector>
+#include <string>
+
+#include <lemon/lp_base.h>
+
+class ClpSimplex;
+
+namespace lemon {
+
+  /// \ingroup lp_group
+  ///
+  /// \brief Interface for the CLP solver
+  /// 
+  /// This class implements an interface for the Clp LP solver.  The
+  /// Clp library is an object oriented lp solver library developed at
+  /// the IBM. The CLP is part of the COIN-OR package and it can be
+  /// used with Common Public License.
+  class LpClp : public LpSolverBase {
+  protected:
+
+    ClpSimplex* _prob;
+
+    std::map<std::string, int> _col_names_ref;
+    std::map<std::string, int> _row_names_ref;
+
+  public:
+
+    typedef LpSolverBase Parent;
+    
+    /// \e
+    LpClp();
+    /// \e
+    LpClp(const LpClp&);
+    /// \e
+    ~LpClp();
+
+  protected:
+
+    mutable double* _primal_ray;
+    mutable double* _dual_ray;
+
+    void _init_temporals();
+    void _clear_temporals();
+
+  protected:
+
+    virtual LpClp* _newSolver() const;
+    virtual LpClp* _cloneSolver() const ;
+
+    virtual int _addCol();
+    virtual int _addRow();
+
+    virtual void _eraseCol(int i);
+    virtual void _eraseRow(int i);
+
+    virtual void _eraseColId(int i);
+    virtual void _eraseRowId(int i);
+
+    virtual void _getColName(int col, std::string& name) const;
+    virtual void _setColName(int col, const std::string& name);
+    virtual int _colByName(const std::string& name) const;
+
+    virtual void _getRowName(int row, std::string& name) const;
+    virtual void _setRowName(int row, const std::string& name);
+    virtual int _rowByName(const std::string& name) const;
+
+    virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e);
+    virtual void _getRowCoeffs(int i, InsertIterator b) const;
+
+    virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e);
+    virtual void _getColCoeffs(int i, InsertIterator b) const;
+
+    virtual void _setCoeff(int row, int col, Value value);
+    virtual Value _getCoeff(int row, int col) const;
+
+    virtual void _setColLowerBound(int i, Value value);
+    virtual Value _getColLowerBound(int i) const;
+    virtual void _setColUpperBound(int i, Value value);
+    virtual Value _getColUpperBound(int i) const;
+
+    virtual void _setRowLowerBound(int i, Value value);
+    virtual Value _getRowLowerBound(int i) const;
+    virtual void _setRowUpperBound(int i, Value value);
+    virtual Value _getRowUpperBound(int i) const;
+
+    virtual void _setObjCoeffs(ExprIterator, ExprIterator);
+    virtual void _getObjCoeffs(InsertIterator) const;
+
+    virtual void _setObjCoeff(int i, Value obj_coef);
+    virtual Value _getObjCoeff(int i) const;
+    
+    virtual void _setSense(Sense sense);
+    virtual Sense _getSense() const;
+
+    virtual SolveExitStatus _solve();
+
+    virtual Value _getPrimal(int i) const;
+    virtual Value _getDual(int i) const;
+
+    virtual Value _getPrimalValue() const;
+
+    virtual Value _getPrimalRay(int i) const;
+    virtual Value _getDualRay(int i) const;
+
+    virtual VarStatus _getColStatus(int i) const;
+    virtual VarStatus _getRowStatus(int i) const;
+    
+    virtual ProblemType _getPrimalType() const;
+    virtual ProblemType _getDualType() const;
+
+    virtual void _clear();
+
+  public:
+
+    ///Solves LP with primal simplex method.
+    SolveExitStatus solvePrimal();
+
+    ///Solves LP with dual simplex method.
+    SolveExitStatus solveDual();
+
+    ///Solves LP with barrier method.
+    SolveExitStatus solveBarrier();
+    
+    ///Returns the constraint identifier understood by CLP.
+    int clpRow(Row r) const { return rows(id(r)); }
+
+    ///Returns the variable identifier understood by CLP.
+    int clpCol(Col c) const { return cols(id(c)); }
+
+    ///Enum for \c messageLevel() parameter
+    enum MessageLevel {
+      /// no output (default value)
+      MESSAGE_NO_OUTPUT = 0,
+      /// print final solution
+      MESSAGE_FINAL_SOLUTION = 1,
+      /// print factorization
+      MESSAGE_FACTORIZATION = 2,
+      /// normal output
+      MESSAGE_NORMAL_OUTPUT = 3,
+      /// verbose output
+      MESSAGE_VERBOSE_OUTPUT = 4
+    };
+    ///Set the verbosity of the messages
+
+    ///Set the verbosity of the messages
+    ///
+    ///\param m is the level of the messages output by the solver routines.
+    void messageLevel(MessageLevel m);
+
+  };
+
+} //END OF NAMESPACE LEMON
+
+#endif //LEMON_LP_CLP_H
+

new file lemon/lp_cplex.cc

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

new file lemon/lp_cplex.h

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

new file lemon/lp_glpk.cc

@@ +1 @@
+/* -*- 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 and MIP solver interface.
+
+#include <lemon/lp_glpk.h>
+#include <glpk.h>
+
+#include <lemon/assert.h>
+
+namespace lemon {
+
+  // GlpkSolverBase members
+
+  GlpkSolverBase::GlpkSolverBase() : SolverBase() {
+    lp = glp_create_prob();
+    glp_create_index(lp);
+  }
+
+  GlpkSolverBase::GlpkSolverBase(const GlpkSolverBase &other) : SolverBase() {
+    lp = glp_create_prob();
+    glp_copy_prob(lp, other.lp, GLP_ON);
+    glp_create_index(lp);
+    rows = other.rows;
+    cols = other.cols;
+  }
+
+  GlpkSolverBase::~GlpkSolverBase() {
+    glp_delete_prob(lp);
+  }
+
+  int GlpkSolverBase::_addCol() {
+    int i = glp_add_cols(lp, 1);
+    glp_set_col_bnds(lp, i, GLP_FR, 0.0, 0.0);
+    return i;
+  }
+
+  int GlpkSolverBase::_addRow() {
+    int i = glp_add_rows(lp, 1);
+    glp_set_row_bnds(lp, i, GLP_FR, 0.0, 0.0);
+    return i;
+  }
+
+  void GlpkSolverBase::_eraseCol(int i) {
+    int ca[2];
+    ca[1] = i;
+    glp_del_cols(lp, 1, ca);
+  }
+
+  void GlpkSolverBase::_eraseRow(int i) {
+    int ra[2];
+    ra[1] = i;
+    glp_del_rows(lp, 1, ra);
+  }
+
+  void GlpkSolverBase::_eraseColId(int i) {
+    cols.eraseIndex(i);
+    cols.shiftIndexes(i);
+  }
+
+  void GlpkSolverBase::_eraseRowId(int i) {
+    rows.eraseIndex(i);
+    rows.shiftIndexes(i);
+  }
+
+  void GlpkSolverBase::_getColName(int c, std::string& name) const {
+    const char *str = glp_get_col_name(lp, c);
+    if (str) name = str;
+    else name.clear();
+  }
+
+  void GlpkSolverBase::_setColName(int c, const std::string & name) {
+    glp_set_col_name(lp, c, const_cast<char*>(name.c_str()));
+
+  }
+
+  int GlpkSolverBase::_colByName(const std::string& name) const {
+    int k = glp_find_col(lp, const_cast<char*>(name.c_str()));
+    return k > 0 ? k : -1;
+  }
+
+  void GlpkSolverBase::_getRowName(int r, std::string& name) const {
+    const char *str = glp_get_row_name(lp, r);
+    if (str) name = str;
+    else name.clear();
+  }
+
+  void GlpkSolverBase::_setRowName(int r, const std::string & name) {
+    glp_set_row_name(lp, r, const_cast<char*>(name.c_str()));
+
+  }
+
+  int GlpkSolverBase::_rowByName(const std::string& name) const {
+    int k = glp_find_row(lp, const_cast<char*>(name.c_str()));
+    return k > 0 ? k : -1;
+  }
+
+  void GlpkSolverBase::_setRowCoeffs(int i, ExprIterator b, ExprIterator e) {
+    std::vector<int> indexes;
+    std::vector<Value> values;
+
+    indexes.push_back(0);
+    values.push_back(0);
+
+    for(ExprIterator it = b; it != e; ++it) {
+      indexes.push_back(it->first);
+      values.push_back(it->second);
+    }
+
+    glp_set_mat_row(lp, i, values.size() - 1,
+                    &indexes.front(), &values.front());
+  }
+
+  void GlpkSolverBase::_getRowCoeffs(int ix, InsertIterator b) const {
+    int length = glp_get_mat_row(lp, ix, 0, 0);
+
+    std::vector<int> indexes(length + 1);
+    std::vector<Value> values(length + 1);
+
+    glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
+
+    for (int i = 1; i <= length; ++i) {
+      *b = std::make_pair(indexes[i], values[i]);
+      ++b;
+    }
+  }
+
+  void GlpkSolverBase::_setColCoeffs(int ix, ExprIterator b,
+                                     ExprIterator e) {
+
+    std::vector<int> indexes;
+    std::vector<Value> values;
+
+    indexes.push_back(0);
+    values.push_back(0);
+
+    for(ExprIterator it = b; it != e; ++it) {
+      indexes.push_back(it->first);
+      values.push_back(it->second);
+    }
+
+    glp_set_mat_col(lp, ix, values.size() - 1,
+                    &indexes.front(), &values.front());
+  }
+
+  void GlpkSolverBase::_getColCoeffs(int ix, InsertIterator b) const {
+    int length = glp_get_mat_col(lp, ix, 0, 0);
+
+    std::vector<int> indexes(length + 1);
+    std::vector<Value> values(length + 1);
+
+    glp_get_mat_col(lp, ix, &indexes.front(), &values.front());
+
+    for (int i = 1; i  <= length; ++i) {
+      *b = std::make_pair(indexes[i], values[i]);
+      ++b;
+    }
+  }
+
+  void GlpkSolverBase::_setCoeff(int ix, int jx, Value value) {
+
+    if (glp_get_num_cols(lp) < glp_get_num_rows(lp)) {
+
+      int length = glp_get_mat_row(lp, ix, 0, 0);
+
+      std::vector<int> indexes(length + 2);
+      std::vector<Value> values(length + 2);
+
+      glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
+
+      //The following code does not suppose that the elements of the
+      //array indexes are sorted
+      bool found = false;
+      for (int i = 1; i  <= length; ++i) {
+        if (indexes[i] == jx) {
+          found = true;
+          values[i] = value;
+          break;
+        }
+      }
+      if (!found) {
+        ++length;
+        indexes[length] = jx;
+        values[length] = value;
+      }
+
+      glp_set_mat_row(lp, ix, length, &indexes.front(), &values.front());
+
+    } else {
+
+      int length = glp_get_mat_col(lp, jx, 0, 0);
+
+      std::vector<int> indexes(length + 2);
+      std::vector<Value> values(length + 2);
+
+      glp_get_mat_col(lp, jx, &indexes.front(), &values.front());
+
+      //The following code does not suppose that the elements of the
+      //array indexes are sorted
+      bool found = false;
+      for (int i = 1; i <= length; ++i) {
+        if (indexes[i] == ix) {
+          found = true;
+          values[i] = value;
+          break;
+        }
+      }
+      if (!found) {
+        ++length;
+        indexes[length] = ix;
+        values[length] = value;
+      }
+
+      glp_set_mat_col(lp, jx, length, &indexes.front(), &values.front());
+    }
+
+  }
+
+  SolverBase::Value GlpkSolverBase::_getCoeff(int ix, int jx) const {
+
+    int length = glp_get_mat_row(lp, ix, 0, 0);
+
+    std::vector<int> indexes(length + 1);
+    std::vector<Value> values(length + 1);
+
+    glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
+
+    for (int i = 1; i  <= length; ++i) {
+      if (indexes[i] == jx) {
+        return values[i];
+      }
+    }
+
+    return 0;
+  }
+
+  void GlpkSolverBase::_setColLowerBound(int i, Value lo) {
+    LEMON_ASSERT(lo != INF, "Invalid bound");
+
+    int b = glp_get_col_type(lp, i);
+    double up = glp_get_col_ub(lp, i);
+    if (lo == -INF) {
+      switch (b) {
+      case GLP_FR:
+      case GLP_LO:
+        glp_set_col_bnds(lp, i, GLP_FR, lo, up);
+        break;
+      case GLP_UP:
+        break;
+      case GLP_DB:
+      case GLP_FX:
+        glp_set_col_bnds(lp, i, GLP_UP, lo, up);
+        break;
+      default:
+        break;
+      }
+    } else {
+      switch (b) {
+      case GLP_FR:
+      case GLP_LO:
+        glp_set_col_bnds(lp, i, GLP_LO, lo, up);
+        break;
+      case GLP_UP:
+      case GLP_DB:
+      case GLP_FX:
+        if (lo == up)
+          glp_set_col_bnds(lp, i, GLP_FX, lo, up);
+        else
+          glp_set_col_bnds(lp, i, GLP_DB, lo, up);
+        break;
+      default:
+        break;
+      }
+    }
+  }
+
+  SolverBase::Value GlpkSolverBase::_getColLowerBound(int i) const {
+    int b = glp_get_col_type(lp, i);
+    switch (b) {
+    case GLP_LO:
+    case GLP_DB:
+    case GLP_FX:
+      return glp_get_col_lb(lp, i);
+    default:
+      return -INF;
+    }
+  }
+
+  void GlpkSolverBase::_setColUpperBound(int i, Value up) {
+    LEMON_ASSERT(up != -INF, "Invalid bound");
+
+    int b = glp_get_col_type(lp, i);
+    double lo = glp_get_col_lb(lp, i);
+    if (up == INF) {
+      switch (b) {
+      case GLP_FR:
+      case GLP_LO:
+        break;
+      case GLP_UP:
+        glp_set_col_bnds(lp, i, GLP_FR, lo, up);
+        break;
+      case GLP_DB:
+      case GLP_FX:
+        glp_set_col_bnds(lp, i, GLP_LO, lo, up);
+        break;
+      default:
+        break;
+      }
+    } else {
+      switch (b) {
+      case GLP_FR:
+        glp_set_col_bnds(lp, i, GLP_UP, lo, up);
+        break;
+      case GLP_UP:
+        glp_set_col_bnds(lp, i, GLP_UP, lo, up);
+        break;
+      case GLP_LO:
+      case GLP_DB:
+      case GLP_FX:
+        if (lo == up)
+          glp_set_col_bnds(lp, i, GLP_FX, lo, up);
+        else
+          glp_set_col_bnds(lp, i, GLP_DB, lo, up);
+        break;
+      default:
+        break;
+      }
+    }
+
+  }
+
+  SolverBase::Value GlpkSolverBase::_getColUpperBound(int i) const {
+    int b = glp_get_col_type(lp, i);
+      switch (b) {
+      case GLP_UP:
+      case GLP_DB:
+      case GLP_FX:
+        return glp_get_col_ub(lp, i);
+      default:
+        return INF;
+      }
+  }
+
+  void GlpkSolverBase::_setRowLowerBound(int i, Value lo) {
+    LEMON_ASSERT(lo != INF, "Invalid bound");
+
+    int b = glp_get_row_type(lp, i);
+    double up = glp_get_row_ub(lp, i);
+    if (lo == -INF) {
+      switch (b) {
+      case GLP_FR:
+      case GLP_LO:
+        glp_set_row_bnds(lp, i, GLP_FR, lo, up);
+        break;
+      case GLP_UP:
+        break;
+      case GLP_DB:
+      case GLP_FX:
+        glp_set_row_bnds(lp, i, GLP_UP, lo, up);
+        break;
+      default:
+        break;
+      }
+    } else {
+      switch (b) {
+      case GLP_FR:
+      case GLP_LO:
+        glp_set_row_bnds(lp, i, GLP_LO, lo, up);
+        break;
+      case GLP_UP:
+      case GLP_DB:
+      case GLP_FX:
+        if (lo == up)
+          glp_set_row_bnds(lp, i, GLP_FX, lo, up);
+        else
+          glp_set_row_bnds(lp, i, GLP_DB, lo, up);
+        break;
+      default:
+        break;
+      }
+    }
+
+  }
+
+  SolverBase::Value GlpkSolverBase::_getRowLowerBound(int i) const {
+    int b = glp_get_row_type(lp, i);
+    switch (b) {
+    case GLP_LO:
+    case GLP_DB:
+    case GLP_FX:
+      return glp_get_row_lb(lp, i);
+    default:
+      return -INF;
+    }
+  }
+
+  void GlpkSolverBase::_setRowUpperBound(int i, Value up) {
+    LEMON_ASSERT(up != -INF, "Invalid bound");
+
+    int b = glp_get_row_type(lp, i);
+    double lo = glp_get_row_lb(lp, i);
+    if (up == INF) {
+      switch (b) {
+      case GLP_FR:
+      case GLP_LO:
+        break;
+      case GLP_UP:
+        glp_set_row_bnds(lp, i, GLP_FR, lo, up);
+        break;
+      case GLP_DB:
+      case GLP_FX:
+        glp_set_row_bnds(lp, i, GLP_LO, lo, up);
+        break;
+      default:
+        break;
+      }
+    } else {
+      switch (b) {
+      case GLP_FR:
+        glp_set_row_bnds(lp, i, GLP_UP, lo, up);
+        break;
+      case GLP_UP:
+        glp_set_row_bnds(lp, i, GLP_UP, lo, up);
+        break;
+      case GLP_LO:
+      case GLP_DB:
+      case GLP_FX:
+        if (lo == up)
+          glp_set_row_bnds(lp, i, GLP_FX, lo, up);
+        else
+          glp_set_row_bnds(lp, i, GLP_DB, lo, up);
+        break;
+      default:
+        break;
+      }
+    }
+  }
+
+  SolverBase::Value GlpkSolverBase::_getRowUpperBound(int i) const {
+    int b = glp_get_row_type(lp, i);
+    switch (b) {
+    case GLP_UP:
+    case GLP_DB:
+    case GLP_FX:
+      return glp_get_row_ub(lp, i);
+    default:
+      return INF;
+    }
+  }
+
+  void GlpkSolverBase::_setObjCoeffs(ExprIterator b, ExprIterator e) {
+    for (int i = 1; i <= glp_get_num_cols(lp); ++i) {
+      glp_set_obj_coef(lp, i, 0.0);
+    }
+    for (ExprIterator it = b; it != e; ++it) {
+      glp_set_obj_coef(lp, it->first, it->second);
+    }
+  }
+
+  void GlpkSolverBase::_getObjCoeffs(InsertIterator b) const {
+    for (int i = 1; i <= glp_get_num_cols(lp); ++i) {
+      Value val = glp_get_obj_coef(lp, i);
+      if (val != 0.0) {
+        *b = std::make_pair(i, val);
+        ++b;
+      }
+    }
+  }
+
+  void GlpkSolverBase::_setObjCoeff(int i, Value obj_coef) {
+    //i = 0 means the constant term (shift)
+    glp_set_obj_coef(lp, i, obj_coef);
+  }
+
+  SolverBase::Value GlpkSolverBase::_getObjCoeff(int i) const {
+    //i = 0 means the constant term (shift)
+    return glp_get_obj_coef(lp, i);
+  }
+
+  void GlpkSolverBase::_setSense(SolverBase::Sense sense) {
+    switch (sense) {
+    case MIN:
+      glp_set_obj_dir(lp, GLP_MIN);
+      break;
+    case MAX:
+      glp_set_obj_dir(lp, GLP_MAX);
+      break;
+    }
+  }
+
+  SolverBase::Sense GlpkSolverBase::_getSense() const {
+    switch(glp_get_obj_dir(lp)) {
+    case GLP_MIN:
+      return MIN;
+    case GLP_MAX:
+      return MAX;
+    default:
+      LEMON_ASSERT(false, "Wrong sense");
+      return SolverBase::Sense();
+    }
+  }
+
+  void GlpkSolverBase::_clear() {
+    glp_erase_prob(lp);
+    rows.clear();
+    cols.clear();
+  }
+
+  // LpGlpk members
+
+  LpGlpk::LpGlpk()
+    : SolverBase(), GlpkSolverBase(), LpSolverBase() {
+    messageLevel(MESSAGE_NO_OUTPUT);
+  }
+
+  LpGlpk::LpGlpk(const LpGlpk& other)
+    : SolverBase(other), GlpkSolverBase(other), LpSolverBase(other) {
+    messageLevel(MESSAGE_NO_OUTPUT);
+  }
+
+  LpGlpk* LpGlpk::_newSolver() const { return new LpGlpk; }
+  LpGlpk* LpGlpk::_cloneSolver() const {return new LpGlpk(*this); }
+
+  void LpGlpk::_clear_temporals() {
+    _primal_ray.clear();
+    _dual_ray.clear();
+  }
+
+  LpGlpk::SolveExitStatus LpGlpk::_solve() {
+    return solvePrimal();
+  }
+
+  LpGlpk::SolveExitStatus LpGlpk::solvePrimal() {
+    _clear_temporals();
+
+    glp_smcp smcp;
+    glp_init_smcp(&smcp);
+
+    switch (_message_level) {
+    case MESSAGE_NO_OUTPUT:
+      smcp.msg_lev = GLP_MSG_OFF;
+      break;
+    case MESSAGE_ERROR_MESSAGE:
+      smcp.msg_lev = GLP_MSG_ERR;
+      break;
+    case MESSAGE_NORMAL_OUTPUT:
+      smcp.msg_lev = GLP_MSG_ON;      
+      break;
+    case MESSAGE_FULL_OUTPUT:
+      smcp.msg_lev = GLP_MSG_ALL;
+      break;
+    }
+
+    if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
+    return SOLVED;
+  }
+
+  LpGlpk::SolveExitStatus LpGlpk::solveDual() {
+    _clear_temporals();
+
+    glp_smcp smcp;
+    glp_init_smcp(&smcp);
+
+    switch (_message_level) {
+    case MESSAGE_NO_OUTPUT:
+      smcp.msg_lev = GLP_MSG_OFF;
+      break;
+    case MESSAGE_ERROR_MESSAGE:
+      smcp.msg_lev = GLP_MSG_ERR;
+      break;
+    case MESSAGE_NORMAL_OUTPUT:
+      smcp.msg_lev = GLP_MSG_ON;      
+      break;
+    case MESSAGE_FULL_OUTPUT:
+      smcp.msg_lev = GLP_MSG_ALL;
+      break;
+    }
+    smcp.meth = GLP_DUAL;
+    
+    if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
+    return SOLVED;
+  }
+
+  LpGlpk::Value LpGlpk::_getPrimal(int i) const {
+    return glp_get_col_prim(lp, i);
+  }
+
+  LpGlpk::Value LpGlpk::_getDual(int i) const {
+    return glp_get_row_dual(lp, i);
+  }
+
+  LpGlpk::Value LpGlpk::_getPrimalValue() const {
+    return glp_get_obj_val(lp);
+  }
+
+  LpGlpk::VarStatus LpGlpk::_getColStatus(int i) const {
+    switch (glp_get_col_stat(lp, i)) {
+    case GLP_BS:
+      return BASIC;
+    case GLP_UP:
+      return UPPER;
+    case GLP_LO:
+      return LOWER;
+    case GLP_NF:
+      return FREE;
+    case GLP_NS:
+      return FIXED;
+    default:
+      LEMON_ASSERT(false, "Wrong column status");
+      return LpGlpk::VarStatus();
+    }
+  }
+
+  LpGlpk::VarStatus LpGlpk::_getRowStatus(int i) const {
+    switch (glp_get_row_stat(lp, i)) {
+    case GLP_BS:
+      return BASIC;
+    case GLP_UP:
+      return UPPER;
+    case GLP_LO:
+      return LOWER;
+    case GLP_NF:
+      return FREE;
+    case GLP_NS:
+      return FIXED;
+    default:
+      LEMON_ASSERT(false, "Wrong row status");
+      return LpGlpk::VarStatus();
+    }
+  }
+
+  LpGlpk::Value LpGlpk::_getPrimalRay(int i) const {
+    if (_primal_ray.empty()) {
+      int row_num = glp_get_num_rows(lp);
+      int col_num = glp_get_num_cols(lp);
+
+      _primal_ray.resize(col_num + 1, 0.0);
+      
+      int index = glp_get_unbnd_ray(lp);
+      if (index != 0) {
+        // The primal ray is found in primal simplex second phase        
+        LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :
+                      glp_get_col_stat(lp, index - row_num)) != GLP_BS, 
+                     "Wrong primal ray");
+
+        bool negate = glp_get_obj_dir(lp) == GLP_MAX;
+        
+        if (index > row_num) {
+          _primal_ray[index - row_num] = 1.0;
+          if (glp_get_col_dual(lp, index - row_num) > 0) {
+            negate = !negate;
+          }
+        } else {
+          if (glp_get_row_dual(lp, index) > 0) {
+            negate = !negate;
+          }
+        }
+        
+        std::vector<int> ray_indexes(row_num + 1);
+        std::vector<Value> ray_values(row_num + 1);
+        int ray_length = glp_eval_tab_col(lp, index, &ray_indexes.front(), 
+                                          &ray_values.front());
+        
+        for (int i = 1; i <= ray_length; ++i) {
+          if (ray_indexes[i] > row_num) {
+            _primal_ray[ray_indexes[i] - row_num] = ray_values[i];
+          }
+        }
+        
+        if (negate) {
+          for (int i = 1; i <= col_num; ++i) {
+            _primal_ray[i] = - _primal_ray[i];
+          }
+        }
+      } else {
+        for (int i = 1; i <= col_num; ++i) {
+          _primal_ray[i] = glp_get_col_prim(lp, i);
+        }
+      }
+    }
+    return _primal_ray[i];
+  }
+
+  LpGlpk::Value LpGlpk::_getDualRay(int i) const {
+    if (_dual_ray.empty()) {
+      int row_num = glp_get_num_rows(lp);
+
+      _dual_ray.resize(row_num + 1, 0.0);
+      
+      int index = glp_get_unbnd_ray(lp);
+      if (index != 0) {
+        // The dual ray is found in dual simplex second phase
+        LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :
+                      glp_get_col_stat(lp, index - row_num)) == GLP_BS, 
+
+                     "Wrong dual ray");
+
+        int idx;
+        bool negate = false;
+        
+        if (index > row_num) {
+          idx = glp_get_col_bind(lp, index - row_num);
+          if (glp_get_col_prim(lp, index - row_num) > 
+              glp_get_col_ub(lp, index - row_num)) { 
+            negate = true;
+          }
+        } else {
+          idx = glp_get_row_bind(lp, index);
+          if (glp_get_row_prim(lp, index) > glp_get_row_ub(lp, index)) {
+            negate = true;
+          }
+        }
+
+        _dual_ray[idx] = negate ?  - 1.0 : 1.0;
+
+        glp_btran(lp, &_dual_ray.front());          
+      } else {
+        double eps = 1e-7;
+        // The dual ray is found in primal simplex first phase
+        // We assume that the glpk minimizes the slack to get feasible solution
+        for (int i = 1; i <= row_num; ++i) {
+          int index = glp_get_bhead(lp, i);
+          if (index <= row_num) {
+            double res = glp_get_row_prim(lp, index);
+            if (res > glp_get_row_ub(lp, index) + eps) {
+              _dual_ray[i] = -1;
+            } else if (res < glp_get_row_lb(lp, index) - eps) {
+              _dual_ray[i] = 1;
+            } else {
+              _dual_ray[i] = 0;
+            }
+            _dual_ray[i] *= glp_get_rii(lp, index);
+          } else {
+            double res = glp_get_col_prim(lp, index - row_num);
+            if (res > glp_get_col_ub(lp, index - row_num) + eps) {
+              _dual_ray[i] = -1;
+            } else if (res < glp_get_col_lb(lp, index - row_num) - eps) {
+              _dual_ray[i] = 1;
+            } else {
+              _dual_ray[i] = 0;
+            }
+            _dual_ray[i] /= glp_get_sjj(lp, index - row_num);
+          }          
+        }
+
+        glp_btran(lp, &_dual_ray.front());
+
+        for (int i = 1; i <= row_num; ++i) {
+          _dual_ray[i] /= glp_get_rii(lp, i);
+        }
+      }
+    }
+    return _dual_ray[i];
+  }
+
+  LpGlpk::ProblemType LpGlpk::_getPrimalType() const {
+    if (glp_get_status(lp) == GLP_OPT)
+      return OPTIMAL;
+    switch (glp_get_prim_stat(lp)) {
+    case GLP_UNDEF:
+      return UNDEFINED;
+    case GLP_FEAS:
+    case GLP_INFEAS:
+      if (glp_get_dual_stat(lp) == GLP_NOFEAS) { 
+        return UNBOUNDED;
+      } else {
+        return UNDEFINED;
+      }
+    case GLP_NOFEAS:
+      return INFEASIBLE;
+    default:
+      LEMON_ASSERT(false, "Wrong primal type");
+      return  LpGlpk::ProblemType();
+    }
+  }
+
+  LpGlpk::ProblemType LpGlpk::_getDualType() const {
+    if (glp_get_status(lp) == GLP_OPT)
+      return OPTIMAL;
+    switch (glp_get_dual_stat(lp)) {
+    case GLP_UNDEF:
+      return UNDEFINED;
+    case GLP_FEAS:
+    case GLP_INFEAS:
+      if (glp_get_prim_stat(lp) == GLP_NOFEAS) { 
+        return UNBOUNDED;
+      } else {
+        return UNDEFINED;
+      }
+    case GLP_NOFEAS:
+      return INFEASIBLE;
+    default:
+      LEMON_ASSERT(false, "Wrong primal type");
+      return  LpGlpk::ProblemType();
+    }
+  }
+
+  void LpGlpk::presolver(bool b) {
+    lpx_set_int_parm(lp, LPX_K_PRESOL, b ? 1 : 0);
+  }
+
+  void LpGlpk::messageLevel(MessageLevel m) {
+    _message_level = m;
+  }
+
+  // MipGlpk members
+
+  MipGlpk::MipGlpk()
+    : SolverBase(), GlpkSolverBase(), MipSolverBase() {
+    messageLevel(MESSAGE_NO_OUTPUT);
+  }
+
+  MipGlpk::MipGlpk(const MipGlpk& other)
+    : SolverBase(), GlpkSolverBase(other), MipSolverBase() {
+    messageLevel(MESSAGE_NO_OUTPUT);
+  }
+
+  void MipGlpk::_setColType(int i, MipGlpk::ColTypes col_type) {
+    switch (col_type) {
+    case INTEGER:
+      glp_set_col_kind(lp, i, GLP_IV);
+      break;
+    case REAL:
+      glp_set_col_kind(lp, i, GLP_CV);
+      break;
+    }
+  }
+
+  MipGlpk::ColTypes MipGlpk::_getColType(int i) const {
+    switch (glp_get_col_kind(lp, i)) {
+    case GLP_IV:
+    case GLP_BV:
+      return INTEGER;
+    default:
+      return REAL;
+    }
+
+  }
+
+  MipGlpk::SolveExitStatus MipGlpk::_solve() {
+    glp_smcp smcp;
+    glp_init_smcp(&smcp);
+
+    switch (_message_level) {
+    case MESSAGE_NO_OUTPUT:
+      smcp.msg_lev = GLP_MSG_OFF;
+      break;
+    case MESSAGE_ERROR_MESSAGE:
+      smcp.msg_lev = GLP_MSG_ERR;
+      break;
+    case MESSAGE_NORMAL_OUTPUT:
+      smcp.msg_lev = GLP_MSG_ON;      
+      break;
+    case MESSAGE_FULL_OUTPUT:
+      smcp.msg_lev = GLP_MSG_ALL;
+      break;
+    }
+    smcp.meth = GLP_DUAL;
+
+    if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
+    if (glp_get_status(lp) != GLP_OPT) return SOLVED;
+
+    glp_iocp iocp;
+    glp_init_iocp(&iocp);
+
+    switch (_message_level) {
+    case MESSAGE_NO_OUTPUT:
+      iocp.msg_lev = GLP_MSG_OFF;
+      break;
+    case MESSAGE_ERROR_MESSAGE:
+      iocp.msg_lev = GLP_MSG_ERR;
+      break;
+    case MESSAGE_NORMAL_OUTPUT:
+      iocp.msg_lev = GLP_MSG_ON;      
+      break;
+    case MESSAGE_FULL_OUTPUT:
+      iocp.msg_lev = GLP_MSG_ALL;
+      break;
+    }
+
+    if (glp_intopt(lp, &iocp) != 0) return UNSOLVED;
+    return SOLVED;
+  }
+
+
+  MipGlpk::ProblemType MipGlpk::_getType() const {
+    switch (glp_get_status(lp)) {
+    case GLP_OPT:
+      switch (glp_mip_status(lp)) {
+      case GLP_UNDEF:
+        return UNDEFINED;
+      case GLP_NOFEAS:
+        return INFEASIBLE;
+      case GLP_FEAS:
+        return FEASIBLE;
+      case GLP_OPT:
+        return OPTIMAL;
+      default:
+        LEMON_ASSERT(false, "Wrong problem type.");
+        return MipGlpk::ProblemType();
+      }
+    case GLP_NOFEAS:
+      return INFEASIBLE;
+    case GLP_INFEAS:
+    case GLP_FEAS:
+      if (glp_get_dual_stat(lp) == GLP_NOFEAS) { 
+        return UNBOUNDED;
+      } else {
+        return UNDEFINED;
+      }
+    default:
+      LEMON_ASSERT(false, "Wrong problem type.");
+      return MipGlpk::ProblemType();
+    }
+  }
+
+  MipGlpk::Value MipGlpk::_getSol(int i) const {
+    return glp_mip_col_val(lp, i);
+  }
+
+  MipGlpk::Value MipGlpk::_getSolValue() const {
+    return glp_mip_obj_val(lp);
+  }
+
+  MipGlpk* MipGlpk::_newSolver() const { return new MipGlpk; }
+  MipGlpk* MipGlpk::_cloneSolver() const {return new MipGlpk(*this); }
+
+  void MipGlpk::messageLevel(MessageLevel m) {
+    _message_level = m;
+  }
+
+} //END OF NAMESPACE LEMON

new file lemon/lp_glpk.h

@@ +1 @@
+/* -*- C++ -*-
+ *
+ * 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 Base interface for the GLPK LP and MIP solver
+  /// 
+  /// This class implements the common interface of the GLPK LP and MIP solver.
+  /// \ingroup lp_group
+  class GlpkSolverBase : virtual public SolverBase {
+  protected:
+    
+    typedef glp_prob LPX;
+    glp_prob* lp;
+    
+    typedef LpSolverBase Parent;
+    
+    GlpkSolverBase();
+    GlpkSolverBase(const GlpkSolverBase&);
+    virtual ~GlpkSolverBase();
+
+  protected:
+
+    virtual int _addCol();
+    virtual int _addRow();
+
+    virtual void _eraseCol(int i);
+    virtual void _eraseRow(int i);
+
+    virtual void _eraseColId(int i);
+    virtual void _eraseRowId(int i);
+
+    virtual void _getColName(int col, std::string& name) const;
+    virtual void _setColName(int col, const std::string& name);
+    virtual int _colByName(const std::string& name) const;
+
+    virtual void _getRowName(int row, std::string& name) const;
+    virtual void _setRowName(int row, const std::string& name);
+    virtual int _rowByName(const std::string& name) const;
+
+    virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e);
+    virtual void _getRowCoeffs(int i, InsertIterator b) const;
+
+    virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e);
+    virtual void _getColCoeffs(int i, InsertIterator b) const;
+
+    virtual void _setCoeff(int row, int col, Value value);
+    virtual Value _getCoeff(int row, int col) const;
+
+    virtual void _setColLowerBound(int i, Value value);
+    virtual Value _getColLowerBound(int i) const;
+
+    virtual void _setColUpperBound(int i, Value value);
+    virtual Value _getColUpperBound(int i) const;
+
+    virtual void _setRowLowerBound(int i, Value value);
+    virtual Value _getRowLowerBound(int i) const;
+
+    virtual void _setRowUpperBound(int i, Value value);
+    virtual Value _getRowUpperBound(int i) const;
+
+    virtual void _setObjCoeffs(ExprIterator b, ExprIterator e);
+    virtual void _getObjCoeffs(InsertIterator b) const;
+    
+    virtual void _setObjCoeff(int i, Value obj_coef);
+    virtual Value _getObjCoeff(int i) const;
+    
+    virtual void _setSense(Sense);
+    virtual Sense _getSense() const;
+
+    virtual void _clear();
+
+  public:
+
+    ///Pointer to the underlying GLPK data structure.
+    LPX *lpx() {return lp;}
+    ///Const pointer to the underlying GLPK data structure.
+    const LPX *lpx() const {return lp;}
+
+    ///Returns the constraint identifier understood by GLPK.
+    int lpxRow(Row r) const { return rows(id(r)); }
+
+    ///Returns the variable identifier understood by GLPK.
+    int lpxCol(Col c) const { return cols(id(c)); }
+
+  };
+  
+  /// \brief Interface for the GLPK LP solver
+  /// 
+  /// This class implements an interface for the GLPK LP solver.
+  ///\ingroup lp_group
+  class LpGlpk : public GlpkSolverBase, public LpSolverBase {
+  public:
+    
+    ///\e
+    LpGlpk();
+    ///\e
+    LpGlpk(const LpGlpk&);
+
+  private:
+    
+    mutable std::vector<double> _primal_ray;
+    mutable std::vector<double> _dual_ray;
+
+    void _clear_temporals();
+    
+  protected:
+
+    virtual LpGlpk* _cloneSolver() const;
+    virtual LpGlpk* _newSolver() const;
+
+    ///\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 VarStatus _getColStatus(int i) const;
+    virtual VarStatus _getRowStatus(int i) const;
+
+    virtual Value _getPrimalRay(int i) const;
+    virtual Value _getDualRay(int i) const;
+    
+    ///\todo It should be clarified
+    ///
+    virtual ProblemType _getPrimalType() const;
+    virtual ProblemType _getDualType() const;
+
+  public:
+
+    ///Solve with primal simplex
+    SolveExitStatus solvePrimal();
+
+    ///Solve with dual simplex
+    SolveExitStatus solveDual();
+
+    ///Turns on or off the presolver
+
+    ///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);
+
+    ///Enum for \c messageLevel() parameter
+    enum MessageLevel {
+      /// no output (default value)
+      MESSAGE_NO_OUTPUT = 0,
+      /// error messages only
+      MESSAGE_ERROR_MESSAGE = 1,
+      /// normal output
+      MESSAGE_NORMAL_OUTPUT = 2,
+      /// full output (includes informational messages)
+      MESSAGE_FULL_OUTPUT = 3
+    };
+
+  private:
+    
+    MessageLevel _message_level;
+
+  public:
+
+    ///Set the verbosity of the messages
+
+    ///Set the verbosity of the messages
+    ///
+    ///\param m is the level of the messages output by the solver routines.
+    void messageLevel(MessageLevel m);
+  };
+
+  /// \brief Interface for the GLPK MIP solver
+  /// 
+  /// This class implements an interface for the GLPK MIP solver.
+  ///\ingroup lp_group
+  class MipGlpk : public GlpkSolverBase, public MipSolverBase {  
+  public:
+  
+    ///\e
+    MipGlpk();
+    ///\e
+    MipGlpk(const MipGlpk&);
+    
+  protected:
+
+    virtual MipGlpk* _cloneSolver() const;
+    virtual MipGlpk* _newSolver() const;
+  
+    virtual ColTypes _getColType(int col) const;
+    virtual void _setColType(int col, ColTypes col_type);
+    
+    virtual SolveExitStatus _solve();
+    virtual ProblemType _getType() const;
+    virtual Value _getSol(int i) const;
+    virtual Value _getSolValue() const;
+
+    ///Enum for \c messageLevel() parameter
+    enum MessageLevel {
+      /// no output (default value)
+      MESSAGE_NO_OUTPUT = 0,
+      /// error messages only
+      MESSAGE_ERROR_MESSAGE = 1,
+      /// normal output
+      MESSAGE_NORMAL_OUTPUT = 2,
+      /// full output (includes informational messages)
+      MESSAGE_FULL_OUTPUT = 3
+    };
+
+  private:
+    
+    MessageLevel _message_level;
+
+  public:
+
+    ///Set the verbosity of the messages
+
+    ///Set the verbosity of the messages
+    ///
+    ///\param m is the level of the messages output by the solver routines.
+    void messageLevel(MessageLevel m);
+  };
+
+
+} //END OF NAMESPACE LEMON
+
+#endif //LEMON_LP_GLPK_H
+

new file lemon/lp_skeleton.cc

@@ +1 @@
+/* -*- 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 {
+
+  int SkeletonSolverBase::_addCol()
+  {
+    return ++col_num;
+  }
+
+  int SkeletonSolverBase::_addRow()
+  {
+    return ++row_num;
+  }
+
+  void SkeletonSolverBase::_eraseCol(int) {}
+  void SkeletonSolverBase::_eraseRow(int) {}
+
+  void SkeletonSolverBase::_getColName(int, std::string &) const {}
+  void SkeletonSolverBase::_setColName(int, const std::string &) {}
+  int SkeletonSolverBase::_colByName(const std::string&) const { return -1; }
+
+  void SkeletonSolverBase::_getRowName(int, std::string &) const {}
+  void SkeletonSolverBase::_setRowName(int, const std::string &) {}
+  int SkeletonSolverBase::_rowByName(const std::string&) const { return -1; }
+
+  void SkeletonSolverBase::_setRowCoeffs(int, ExprIterator, ExprIterator) {}
+  void SkeletonSolverBase::_getRowCoeffs(int, InsertIterator) const {}
+
+  void SkeletonSolverBase::_setColCoeffs(int, ExprIterator, ExprIterator) {}
+  void SkeletonSolverBase::_getColCoeffs(int, InsertIterator) const {}
+
+  void SkeletonSolverBase::_setCoeff(int, int, Value) {}
+  SkeletonSolverBase::Value SkeletonSolverBase::_getCoeff(int, int) const 
+  { return 0; }
+
+  void SkeletonSolverBase::_setColLowerBound(int, Value) {}
+  SkeletonSolverBase::Value SkeletonSolverBase::_getColLowerBound(int) const
+  {  return 0; }
+
+  void SkeletonSolverBase::_setColUpperBound(int, Value) {}
+  SkeletonSolverBase::Value SkeletonSolverBase::_getColUpperBound(int) const
+  {  return 0; }
+
+  void SkeletonSolverBase::_setRowLowerBound(int, Value) {}
+  SkeletonSolverBase::Value SkeletonSolverBase::_getRowLowerBound(int) const
+  {  return 0; }
+
+  void SkeletonSolverBase::_setRowUpperBound(int, Value) {}
+  SkeletonSolverBase::Value SkeletonSolverBase::_getRowUpperBound(int) const
+  {  return 0; }
+
+  void SkeletonSolverBase::_setObjCoeffs(ExprIterator, ExprIterator) {}
+  void SkeletonSolverBase::_getObjCoeffs(InsertIterator) const {};
+
+  void SkeletonSolverBase::_setObjCoeff(int, Value) {}
+  SkeletonSolverBase::Value SkeletonSolverBase::_getObjCoeff(int) const
+  {  return 0; }
+
+  void SkeletonSolverBase::_setSense(Sense) {}
+  SkeletonSolverBase::Sense SkeletonSolverBase::_getSense() const
+  { return MIN; }
+
+  void SkeletonSolverBase::_clear() {
+    row_num = col_num = 0;
+  }
+
+  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::Value LpSkeleton::_getPrimalRay(int) const { return 0; }
+  LpSkeleton::Value LpSkeleton::_getDualRay(int) const { return 0; }
+
+  LpSkeleton::ProblemType LpSkeleton::_getPrimalType() const
+  { return UNDEFINED; }
+
+  LpSkeleton::ProblemType LpSkeleton::_getDualType() const
+  { return UNDEFINED; }
+
+  LpSkeleton::VarStatus LpSkeleton::_getColStatus(int) const
+  { return BASIC; }
+
+  LpSkeleton::VarStatus LpSkeleton::_getRowStatus(int) const
+  { return BASIC; }
+
+  LpSkeleton* LpSkeleton::_newSolver() const
+  { return static_cast<LpSkeleton*>(0); }
+
+  LpSkeleton* LpSkeleton::_cloneSolver() const
+  { return static_cast<LpSkeleton*>(0); }
+
+  MipSkeleton::SolveExitStatus MipSkeleton::_solve()
+  { return SOLVED; }
+
+  MipSkeleton::Value MipSkeleton::_getSol(int) const { return 0; }
+  MipSkeleton::Value MipSkeleton::_getSolValue() const { return 0; }
+
+  MipSkeleton::ProblemType MipSkeleton::_getType() const
+  { return UNDEFINED; }
+
+  MipSkeleton* MipSkeleton::_newSolver() const
+  { return static_cast<MipSkeleton*>(0); }
+
+  MipSkeleton* MipSkeleton::_cloneSolver() const
+  { return static_cast<MipSkeleton*>(0); }
+
+} //namespace lemon
+

new file lemon/lp_skeleton.h

@@ +1 @@
+/* -*- C++ -*-
+ *
+ * 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 SkeletonSolverBase : public virtual SolverBase {
+    int col_num,row_num;
+    
+  protected:
+
+    SkeletonSolverBase() 
+      : col_num(-1), row_num(-1) {} 
+
+    /// \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 _getRowName(int row, std::string& name) const;
+    /// \e
+    virtual void _setRowName(int row, const std::string& name);
+    /// \e
+    virtual int _rowByName(const std::string& name) const;
+
+    /// \e
+    virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e);
+    /// \e
+    virtual void _getRowCoeffs(int i, InsertIterator b) const;
+    /// \e
+    virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e);
+    /// \e
+    virtual void _getColCoeffs(int i, InsertIterator 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 constraint (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 lower bound of a constraint (row) is an 
+    /// extended number of type Value, i.e. a finite number of type 
+    /// Value or -\ref INF.
+    virtual Value _getRowLowerBound(int i) const;
+
+    /// 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 or \ref INF.
+    virtual void _setRowUpperBound(int i, Value value);
+    /// \e
+
+    /// The upper bound of a constraint (row) is an 
+    /// extended number of type Value, i.e. a finite number of type 
+    /// Value or \ref INF.
+    virtual Value _getRowUpperBound(int i) const;
+
+    /// \e
+    virtual void _setObjCoeffs(ExprIterator b, ExprIterator e);
+    /// \e
+    virtual void _getObjCoeffs(InsertIterator b) const;
+
+    /// \e
+    virtual void _setObjCoeff(int i, Value obj_coef);
+    /// \e
+    virtual Value _getObjCoeff(int i) const;
+
+    ///\e
+    virtual void _setSense(Sense);
+    ///\e
+    virtual Sense _getSense() const;
+
+    ///\e
+    virtual void _clear();
+
+  };
+
+  /// \brief Interface for a skeleton LP solver
+  /// 
+  /// This class implements an interface for a skeleton LP solver.
+  ///\ingroup lp_group
+  class LpSkeleton : public SkeletonSolverBase, public LpSolverBase {
+  public:
+    LpSkeleton() : SkeletonSolverBase(), LpSolverBase() {}
+
+  protected:
+
+    ///\e
+    virtual SolveExitStatus _solve();
+
+    ///\e
+    virtual Value _getPrimal(int i) const;
+    ///\e
+    virtual Value _getDual(int i) const;
+
+    ///\e
+    virtual Value _getPrimalValue() const;
+
+    ///\e
+    virtual Value _getPrimalRay(int i) const;
+    ///\e
+    virtual Value _getDualRay(int i) const;
+
+    ///\e
+    virtual ProblemType _getPrimalType() const;
+    ///\e
+    virtual ProblemType _getDualType() const;
+
+    ///\e
+    virtual VarStatus _getColStatus(int i) const;
+    ///\e
+    virtual VarStatus _getRowStatus(int i) const;
+
+    ///\e
+    virtual LpSkeleton* _newSolver() const;
+    ///\e
+    virtual LpSkeleton* _cloneSolver() const;
+
+  };
+
+  /// \brief Interface for a skeleton MIP solver
+  /// 
+  /// This class implements an interface for a skeleton MIP solver.
+  ///\ingroup lp_group
+  class MipSkeleton : public SkeletonSolverBase, public MipSolverBase {
+  public:
+    MipSkeleton() : SkeletonSolverBase(), MipSolverBase() {}
+
+  protected:
+    ///\e
+    
+    ///\bug Wrong interface
+    ///
+    virtual SolveExitStatus _solve();
+
+    ///\e
+
+    ///\bug Wrong interface
+    ///
+    virtual Value _getSol(int i) const;
+
+    ///\e
+
+    ///\bug Wrong interface
+    ///
+    virtual Value _getSolValue() const;
+
+    ///\e
+
+    ///\bug Wrong interface
+    ///
+    virtual ProblemType _getType() const;
+
+    ///\e
+    virtual MipSkeleton* _newSolver() const;
+
+    ///\e
+    virtual MipSkeleton* _cloneSolver() const;
+
+  };
+
+} //namespace lemon
+
+#endif // LEMON_LP_SKELETON

new file lemon/lp_soplex.cc

@@ +1 @@
+/* -*- 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() 
+    : SolverBase(), LpSolverBase() {
+    soplex = new soplex::SoPlex;
+  }
+
+  LpSoplex::~LpSoplex() {
+    delete soplex;
+  }
+
+  LpSoplex::LpSoplex(const LpSoplex& lp) 
+    : SolverBase(), LpSolverBase() {
+    rows = lp.rows;
+    cols = lp.cols;
+
+    soplex = new soplex::SoPlex;
+    (*static_cast<soplex::SPxLP*>(soplex)) = *(lp.soplex);
+
+    _col_names = lp._col_names;
+    _col_names_ref = lp._col_names_ref;
+
+    _row_names = lp._row_names;
+    _row_names_ref = lp._row_names_ref;
+
+  }
+  
+  void LpSoplex::_clear_temporals() {
+    _primal_values.clear();
+    _dual_values.clear();
+  }
+
+  LpSoplex* LpSoplex::_newSolver() const {
+    LpSoplex* newlp = new LpSoplex();
+    return newlp;
+  }
+
+  LpSoplex* LpSoplex::_cloneSolver() const {
+    LpSoplex* newlp = new LpSoplex(*this);
+    return newlp;
+  }
+
+  int LpSoplex::_addCol() {
+    soplex::LPCol c;
+    c.setLower(-soplex::infinity);
+    c.setUpper(soplex::infinity);
+    soplex->addCol(c);
+
+    _col_names.push_back(std::string());
+
+    return soplex->nCols() - 1;
+  }
+
+  int LpSoplex::_addRow() {
+    soplex::LPRow r;
+    r.setLhs(-soplex::infinity);
+    r.setRhs(soplex::infinity);
+    soplex->addRow(r);
+
+    _row_names.push_back(std::string());
+
+    return soplex->nRows() - 1;
+  }
+
+
+  void LpSoplex::_eraseCol(int i) {
+    soplex->removeCol(i);
+    _col_names_ref.erase(_col_names[i]);
+    _col_names[i] = _col_names.back();
+    _col_names_ref[_col_names.back()] = i;
+    _col_names.pop_back();
+  }
+
+  void LpSoplex::_eraseRow(int i) {
+    soplex->removeRow(i);
+    _row_names_ref.erase(_row_names[i]);
+    _row_names[i] = _row_names.back();
+    _row_names_ref[_row_names.back()] = i;
+    _row_names.pop_back();
+  }
+
+  void LpSoplex::_eraseColId(int i) {
+    cols.eraseIndex(i);
+    cols.relocateIndex(i, cols.maxIndex());
+  }
+  void LpSoplex::_eraseRowId(int i) {
+    rows.eraseIndex(i);
+    rows.relocateIndex(i, rows.maxIndex());
+  }
+
+  void LpSoplex::_getColName(int c, std::string &name) const {
+    name = _col_names[c];
+  }
+
+  void LpSoplex::_setColName(int c, const std::string &name) {
+    _col_names_ref.erase(_col_names[c]);
+    _col_names[c] = name;
+    if (!name.empty()) {
+      _col_names_ref.insert(std::make_pair(name, c));
+    }
+  }
+
+  int LpSoplex::_colByName(const std::string& name) const {
+    std::map<std::string, int>::const_iterator it =
+      _col_names_ref.find(name);
+    if (it != _col_names_ref.end()) {
+      return it->second;
+    } else {
+      return -1;
+    }
+  }
+
+  void LpSoplex::_getRowName(int r, std::string &name) const {
+    name = _row_names[r];
+  }
+
+  void LpSoplex::_setRowName(int r, const std::string &name) {
+    _row_names_ref.erase(_row_names[r]);
+    _row_names[r] = name;
+    if (!name.empty()) {
+      _row_names_ref.insert(std::make_pair(name, r));
+    }
+  }
+
+  int LpSoplex::_rowByName(const std::string& name) const {
+    std::map<std::string, int>::const_iterator it =
+      _row_names_ref.find(name);
+    if (it != _row_names_ref.end()) {
+      return it->second;
+    } else {
+      return -1;
+    }
+  }
+
+
+  void LpSoplex::_setRowCoeffs(int i, ExprIterator b, ExprIterator e) {
+    for (int j = 0; j < soplex->nCols(); ++j) {
+      soplex->changeElement(i, j, 0.0);
+    }
+    for(ExprIterator it = b; it != e; ++it) {
+      soplex->changeElement(i, it->first, it->second);
+    }
+  }
+
+  void LpSoplex::_getRowCoeffs(int i, InsertIterator 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, ExprIterator b, ExprIterator e) {
+    for (int i = 0; i < soplex->nRows(); ++i) {
+      soplex->changeElement(i, j, 0.0);
+    }
+    for(ExprIterator it = b; it != e; ++it) {
+      soplex->changeElement(it->first, j, it->second);
+    }
+  }
+
+  void LpSoplex::_getColCoeffs(int i, InsertIterator 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);
+  }
+
+  LpSoplex::Value LpSoplex::_getCoeff(int i, int j) const {
+    return soplex->rowVector(i)[j];
+  }
+
+  void LpSoplex::_setColLowerBound(int i, Value value) {
+    LEMON_ASSERT(value != INF, "Invalid bound");
+    soplex->changeLower(i, value != -INF ? value : -soplex::infinity);
+  }
+
+  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) {
+    LEMON_ASSERT(value != -INF, "Invalid bound");
+    soplex->changeUpper(i, value != INF ? value : soplex::infinity);
+  }
+
+  LpSoplex::Value LpSoplex::_getColUpperBound(int i) const {
+    double value = soplex->upper(i);
+    return value != soplex::infinity ? value : INF;
+  }
+
+  void LpSoplex::_setRowLowerBound(int i, Value lb) {
+    LEMON_ASSERT(lb != INF, "Invalid bound");
+    soplex->changeRange(i, lb != -INF ? lb : -soplex::infinity, soplex->rhs(i));
+  }
+
+  LpSoplex::Value LpSoplex::_getRowLowerBound(int i) const {
+    double res = soplex->lhs(i);
+    return res == -soplex::infinity ? -INF : res;
+  }
+
+  void LpSoplex::_setRowUpperBound(int i, Value ub) {
+    LEMON_ASSERT(ub != -INF, "Invalid bound");
+    soplex->changeRange(i, soplex->lhs(i), ub != INF ? ub : soplex::infinity);
+  }
+
+  LpSoplex::Value LpSoplex::_getRowUpperBound(int i) const {
+    double res = soplex->rhs(i);
+    return res == soplex::infinity ? INF : res;
+  }
+
+  void LpSoplex::_setObjCoeffs(ExprIterator b, ExprIterator e) {
+    for (int j = 0; j < soplex->nCols(); ++j) {
+      soplex->changeObj(j, 0.0);
+    }
+    for (ExprIterator it = b; it != e; ++it) {
+      soplex->changeObj(it->first, it->second);
+    }
+  }
+
+  void LpSoplex::_getObjCoeffs(InsertIterator b) const {
+    for (int j = 0; j < soplex->nCols(); ++j) {
+      Value coef = soplex->obj(j);
+      if (coef != 0.0) {
+        *b = std::make_pair(j, coef);
+        ++b;
+      }
+    }
+  }
+
+  void LpSoplex::_setObjCoeff(int i, Value obj_coef) {
+    soplex->changeObj(i, obj_coef);
+  }
+
+  LpSoplex::Value LpSoplex::_getObjCoeff(int i) const {
+    return soplex->obj(i);
+  }
+
+  LpSoplex::SolveExitStatus LpSoplex::_solve() {
+    
+    _clear_temporals();
+
+    soplex::SPxSolver::Status status = soplex->solve();
+
+    switch (status) {
+    case soplex::SPxSolver::OPTIMAL:
+    case soplex::SPxSolver::INFEASIBLE:
+    case soplex::SPxSolver::UNBOUNDED:
+      return SOLVED;
+    default:
+      return UNSOLVED;
+    }
+  }
+
+  LpSoplex::Value LpSoplex::_getPrimal(int i) const {
+    if (_primal_values.empty()) {
+      _primal_values.resize(soplex->nCols());
+      soplex::Vector pv(_primal_values.size(), &_primal_values.front());
+      soplex->getPrimal(pv);
+    }
+    return _primal_values[i];
+  }
+
+  LpSoplex::Value LpSoplex::_getDual(int i) const {
+    if (_dual_values.empty()) {
+      _dual_values.resize(soplex->nRows());
+      soplex::Vector dv(_dual_values.size(), &_dual_values.front());
+      soplex->getDual(dv);
+    }
+    return _dual_values[i];
+  }
+
+  LpSoplex::Value LpSoplex::_getPrimalValue() const {
+    return soplex->objValue();
+  }
+
+  LpSoplex::VarStatus LpSoplex::_getColStatus(int i) const {
+    switch (soplex->getBasisColStatus(i)) {
+    case soplex::SPxSolver::BASIC:
+      return BASIC;
+    case soplex::SPxSolver::ON_UPPER:
+      return UPPER;
+    case soplex::SPxSolver::ON_LOWER:
+      return LOWER;
+    case soplex::SPxSolver::FIXED:
+      return FIXED;
+    case soplex::SPxSolver::ZERO:
+      return FREE;
+    default:
+      LEMON_ASSERT(false, "Wrong column status");
+      return VarStatus();
+    }
+  }
+
+  LpSoplex::VarStatus LpSoplex::_getRowStatus(int i) const {
+    switch (soplex->getBasisRowStatus(i)) {
+    case soplex::SPxSolver::BASIC:
+      return BASIC;
+    case soplex::SPxSolver::ON_UPPER:
+      return UPPER;
+    case soplex::SPxSolver::ON_LOWER:
+      return LOWER;
+    case soplex::SPxSolver::FIXED:
+      return FIXED;
+    case soplex::SPxSolver::ZERO:
+      return FREE;
+    default:
+      LEMON_ASSERT(false, "Wrong row status");
+      return VarStatus();
+    }
+  }
+
+  LpSoplex::Value LpSoplex::_getPrimalRay(int i) const {
+    if (_primal_ray.empty()) {
+      _primal_ray.resize(soplex->nCols());
+      soplex::Vector pv(_primal_ray.size(), &_primal_ray.front());
+      soplex->getDualfarkas(pv);
+    }
+    return _primal_ray[i];
+  }
+
+  LpSoplex::Value LpSoplex::_getDualRay(int i) const {
+    if (_dual_ray.empty()) {
+      _dual_ray.resize(soplex->nRows());
+      soplex::Vector dv(_dual_ray.size(), &_dual_ray.front());
+      soplex->getDualfarkas(dv);
+    }
+    return _dual_ray[i];
+  }
+
+  LpSoplex::ProblemType LpSoplex::_getPrimalType() const {
+    switch (soplex->status()) {
+    case soplex::SPxSolver::OPTIMAL:
+      return OPTIMAL;
+    case soplex::SPxSolver::UNBOUNDED:
+      return UNBOUNDED;
+    case soplex::SPxSolver::INFEASIBLE:
+      return INFEASIBLE;
+    default:
+      return UNDEFINED;
+    }
+  }
+
+  LpSoplex::ProblemType LpSoplex::_getDualType() const {
+    switch (soplex->status()) {
+    case soplex::SPxSolver::OPTIMAL:
+      return OPTIMAL;
+    case soplex::SPxSolver::UNBOUNDED:
+      return UNBOUNDED;
+    case soplex::SPxSolver::INFEASIBLE:
+      return INFEASIBLE;
+    default:
+      return UNDEFINED;
+    }
+  }
+
+  void LpSoplex::_setSense(Sense sense) {
+    switch (sense) {
+    case MIN:
+      soplex->changeSense(soplex::SPxSolver::MINIMIZE);
+      break;
+    case MAX:
+      soplex->changeSense(soplex::SPxSolver::MAXIMIZE);
+    }
+  }
+
+  LpSoplex::Sense LpSoplex::_getSense() const {
+    switch (soplex->spxSense()) {
+    case soplex::SPxSolver::MAXIMIZE:
+      return MAX;
+    case soplex::SPxSolver::MINIMIZE:
+      return MIN;
+    default:
+      LEMON_ASSERT(false, "Wrong sense.");
+      return LpSoplex::Sense();
+    }
+  }
+
+  void LpSoplex::_clear() {
+    soplex->clear();
+    _col_names.clear();
+    _col_names_ref.clear();
+    _row_names.clear();
+    _row_names_ref.clear();
+    cols.clear();
+    rows.clear();
+    _clear_temporals();
+  }
+
+} //namespace lemon
+

new file lemon/lp_soplex.h

@@ +1 @@
+/* -*- C++ -*-
+ *
+ * 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 : public LpSolverBase {
+  private:
+
+    soplex::SoPlex* soplex;
+
+    std::vector<std::string> _col_names;
+    std::map<std::string, int> _col_names_ref;
+
+    std::vector<std::string> _row_names;
+    std::map<std::string, int> _row_names_ref;
+
+  private:
+
+    // these values cannot be retrieved element by element
+    mutable std::vector<Value> _primal_values;
+    mutable std::vector<Value> _dual_values;
+
+    mutable std::vector<Value> _primal_ray;
+    mutable std::vector<Value> _dual_ray;
+
+    void _clear_temporals();
+
+  public:
+
+    typedef LpSolverBase Parent;
+    
+
+    /// \e
+    LpSoplex();
+    /// \e
+    LpSoplex(const LpSoplex&);
+    /// \e
+    ~LpSoplex();
+
+  protected:
+
+    virtual LpSoplex* _newSolver() const;
+    virtual LpSoplex* _cloneSolver() const ;
+
+    virtual int _addCol();
+    virtual int _addRow();
+
+    virtual void _eraseCol(int i);
+    virtual void _eraseRow(int i);
+
+    virtual void _eraseColId(int i);
+    virtual void _eraseRowId(int i);
+
+    virtual void _getColName(int col, std::string& name) const;
+    virtual void _setColName(int col, const std::string& name);
+    virtual int _colByName(const std::string& name) const;
+
+    virtual void _getRowName(int row, std::string& name) const;
+    virtual void _setRowName(int row, const std::string& name);
+    virtual int _rowByName(const std::string& name) const;
+
+    virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e);
+    virtual void _getRowCoeffs(int i, InsertIterator b) const;
+
+    virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e);
+    virtual void _getColCoeffs(int i, InsertIterator b) const;
+
+    virtual void _setCoeff(int row, int col, Value value);
+    virtual Value _getCoeff(int row, int col) const;
+
+    virtual void _setColLowerBound(int i, Value value);
+    virtual Value _getColLowerBound(int i) const;
+    virtual void _setColUpperBound(int i, Value value);
+    virtual Value _getColUpperBound(int i) const;
+
+    virtual void _setRowLowerBound(int i, Value value);
+    virtual Value _getRowLowerBound(int i) const;
+    virtual void _setRowUpperBound(int i, Value value);
+    virtual Value _getRowUpperBound(int i) const;
+
+    virtual void _setObjCoeffs(ExprIterator b, ExprIterator e);
+    virtual void _getObjCoeffs(InsertIterator b) const;
+
+    virtual void _setObjCoeff(int i, Value obj_coef);
+    virtual Value _getObjCoeff(int i) const;
+    
+    virtual void _setSense(Sense sense);
+    virtual Sense _getSense() const;
+
+    virtual SolveExitStatus _solve();
+    virtual Value _getPrimal(int i) const;
+    virtual Value _getDual(int i) const;
+
+    virtual Value _getPrimalValue() const;
+
+    virtual Value _getPrimalRay(int i) const;
+    virtual Value _getDualRay(int i) const;
+
+    virtual VarStatus _getColStatus(int i) const;
+    virtual VarStatus _getRowStatus(int i) const;
+    
+    virtual ProblemType _getPrimalType() const;
+    virtual ProblemType _getDualType() const;
+    
+    virtual void _clear();
+
+  };
+
+} //END OF NAMESPACE LEMON
+
+#endif //LEMON_LP_SOPLEX_H
+

new file m4/lx_check_clp.m4

@@ +1 @@
+AC_DEFUN([LX_CHECK_CLP],
+[
+  AC_ARG_WITH([clp],
+AS_HELP_STRING([--with-clp@<:@=PREFIX@:>@], [search for CLP under PREFIX or under the default search paths if PREFIX is not given @<:@default@:>@])
+AS_HELP_STRING([--without-clp], [disable checking for CLP]),
+              [], [with_clp=yes])
+
+  AC_ARG_WITH([clp-includedir],
+AS_HELP_STRING([--with-clp-includedir=DIR], [search for CLP headers in DIR]),
+              [], [with_clp_includedir=no])
+
+  AC_ARG_WITH([clp-libdir],
+AS_HELP_STRING([--with-clp-libdir=DIR], [search for CLP libraries in DIR]),
+              [], [with_clp_libdir=no])
+
+  lx_clp_found=no
+  if test x"$with_clp" != x"no"; then
+    AC_MSG_CHECKING([for CLP])
+
+    if test x"$with_clp_includedir" != x"no"; then
+      CLP_CXXFLAGS="-I$with_clp_includedir"
+    elif test x"$with_clp" != x"yes"; then
+      CLP_CXXFLAGS="-I$with_clp/include"
+    fi
+
+    if test x"$with_clp_libdir" != x"no"; then
+      CLP_LDFLAGS="-L$with_clp_libdir"
+    elif test x"$with_clp" != x"yes"; then
+      CLP_LDFLAGS="-L$with_clp/lib"
+    fi
+    CLP_LIBS="-lClp -lCoinUtils -lm"
+
+    lx_save_cxxflags="$CXXFLAGS"
+    lx_save_ldflags="$LDFLAGS"
+    lx_save_libs="$LIBS"
+    CXXFLAGS="$CLP_CXXFLAGS"
+    LDFLAGS="$CLP_LDFLAGS"
+    LIBS="$CLP_LIBS"
+
+    lx_clp_test_prog='
+      #include <coin/ClpModel.hpp>
+
+      int main(int argc, char** argv)
+      {
+        ClpModel clp;
+        return 0;
+      }'
+
+    AC_LANG_PUSH(C++)
+    AC_LINK_IFELSE([$lx_clp_test_prog], [lx_clp_found=yes], [lx_clp_found=no])
+    AC_LANG_POP(C++)
+
+    CXXFLAGS="$lx_save_cxxflags"
+    LDFLAGS="$lx_save_ldflags"
+    LIBS="$lx_save_libs"
+
+    if test x"$lx_clp_found" = x"yes"; then
+      AC_DEFINE([HAVE_CLP], [1], [Define to 1 if you have CLP.])
+      AC_MSG_RESULT([yes])
+    else
+      CLP_CXXFLAGS=""
+      CLP_LDFLAGS=""
+      CLP_LIBS=""
+      AC_MSG_RESULT([no])
+    fi
+  fi
+  CLP_LIBS="$CLP_LDFLAGS $CLP_LIBS"
+  AC_SUBST(CLP_CXXFLAGS)
+  AC_SUBST(CLP_LIBS)
+  AM_CONDITIONAL([HAVE_CLP], [test x"$lx_clp_found" = x"yes"])
+])

test/CMakeLists.txt

@@ -15 +15 @@
   graph_utils_test
   heap_test
   kruskal_test
+  lp_test
   maps_test
+  mip_test
   random_test
   path_test
   time_measure_test

test/Makefile.am

@@ -18 +18 @@
         test/graph_utils_test \
         test/heap_test \
         test/kruskal_test \
+        test/lp_test \
         test/maps_test \
+        test/mip_test \
         test/random_test \
         test/path_test \
         test/test_tools_fail \
@@ -41 +43 @@
 test_graph_utils_test_SOURCES = test/graph_utils_test.cc
 test_heap_test_SOURCES = test/heap_test.cc
 test_kruskal_test_SOURCES = test/kruskal_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_path_test_SOURCES = test/path_test.cc
 test_random_test_SOURCES = test/random_test.cc
 test_test_tools_fail_SOURCES = test/test_tools_fail.cc

new file test/lp_test.cc

@@ +1 @@
+/* -*- 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>
+#include <lemon/lp_utils.h>
+#include <lemon/lgf_writer.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
+
+#ifdef HAVE_CLP
+#include <lemon/lp_clp.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);
+
+  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=12;
+    e[p1]+=2;
+    *e+=12;
+    e[p1]-=2;
+    *e-=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=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;
+
+
+    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(const_cast<const LpSolverBase::Expr&>(e)[p2]==0, 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
+       );
+  }
+
+}
+
+void solveAndCheck(LpSolverBase& lp, LpSolverBase::ProblemType stat,
+                   double exp_opt) {
+  using std::string;
+  lp.solve();
+  
+  std::ostringstream buf;
+  buf << "PrimalType should be: " << int(stat) << int(lp.primalType());
+
+  check(lp.primalType()==stat, buf.str());
+
+  if (stat ==  LpSolverBase::OPTIMAL) {
+    std::ostringstream sbuf;
+    sbuf << "Wrong optimal value: the right optimum is " << exp_opt;
+    check(std::abs(lp.primal()-exp_opt) < 1e-3, sbuf.str());
+  }
+}
+
+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+2*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.sense(lp.MAX);
+
+  //Testing the problem retrieving routines
+  check(lp.objCoeff(x1)==1,"First term should be 1 in the obj function!");
+  check(lp.sense() == lp.MAX,"This is a maximization!");
+  check(lp.coeff(upright,x1)==1,"The coefficient in question is 1!");
+  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.");
+  check(lp.rowLowerBound(upright) == -LpSolverBase::INF,
+        "The lower bound for the first row should be -infty.");
+  check(lp.rowUpperBound(upright)==1,
+        "The upper bound for the first row should be 1.");
+  LpSolverBase::Expr e = lp.row(upright);
+  check(e[x1] == 1, "The first coefficient should 1.");
+  check(e[x2] == 2, "The second coefficient should 1.");
+
+  lp.row(upright, x1+x2 <=1);
+  e = lp.row(upright);
+  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[upright] == 1, "The first coefficient should 1.");
+
+  LpSolverBase* clp = lp.cloneSolver();
+
+  //Testing the problem retrieving routines
+  check(clp->objCoeff(x1)==1,"First term should be 1 in the obj function!");
+  check(clp->sense() == clp->MAX,"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.");
+
+  check(lp.rowLowerBound(upright)==-LpSolverBase::INF,
+        "The lower bound for the first row should be -infty.");
+  check(lp.rowUpperBound(upright)==1,
+        "The upper bound for the first row should be 1.");
+  e = clp->row(upright);
+  check(e[x1] == 1, "The first coefficient should 1.");
+  check(e[x2] == 1, "The second coefficient should 1.");
+
+  de = clp->col(x1);
+  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.sense(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.erase(upright);
+  lp.sense(lp.MAX);
+  expected_opt=LpSolverBase::INF;
+  solveAndCheck(lp, LpSolverBase::UNBOUNDED, expected_opt);
+
+  //Infeasibilty
+  lp.addRow(x1+x2 <=-2);
+  solveAndCheck(lp, LpSolverBase::INFEASIBLE, expected_opt);
+
+}
+
+int main()
+{
+  LpSkeleton lp_skel;
+  lpTest(lp_skel);
+
+#ifdef HAVE_GLPK
+  {
+    std::cerr << "glpk" << std::endl;
+    LpGlpk lp_glpk1,lp_glpk2;
+    lpTest(lp_glpk1);
+    aTest(lp_glpk2);
+  }
+#endif
+
+#ifdef HAVE_CPLEX
+  try {
+    std::cerr << "cplex" << std::endl;
+    LpCplex lp_cplex1,lp_cplex2;
+    lpTest(lp_cplex1);
+    aTest(lp_cplex2);
+  } catch (CplexEnv::LicenseError& error) {
+#ifdef LEMON_FORCE_CPLEX_CHECK
+    check(false, error.what());
+#else
+    std::cerr << error.what() << std::endl;
+    std::cerr << "Cplex license check failed, lp check skipped" << std::endl;
+#endif
+  }
+#endif
+
+#ifdef HAVE_SOPLEX
+  {
+    std::cerr << "soplex" << std::endl;
+    LpSoplex lp_soplex1,lp_soplex2;
+    lpTest(lp_soplex1);
+    aTest(lp_soplex2);
+  }
+#endif
+
+#ifdef HAVE_CLP
+  {
+    std::cerr << "clp" << std::endl;
+    LpClp lp_clp1,lp_clp2;
+    lpTest(lp_clp1);
+    aTest(lp_clp2);
+  }
+#endif
+
+  return 0;
+}

new file test/mip_test.cc

@@ +1 @@
+/* -*- C++ -*-
+ *
+ * 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/lp_cplex.h>
+#endif
+
+#ifdef HAVE_GLPK
+#include <lemon/lp_glpk.h>
+#endif
+
+
+using namespace lemon;
+
+void solveAndCheck(MipSolverBase& mip, MipSolverBase::ProblemType stat, 
+                   double exp_opt) {
+  using std::string;
+
+  mip.solve();
+  //int decimal,sign;
+  std::ostringstream buf;
+  buf << "Type should be: " << int(stat)<<" and it is "<<int(mip.type());
+
+
+  //  itoa(stat,buf1, 10);
+  check(mip.type()==stat, buf.str());
+
+  if (stat ==  MipSolverBase::OPTIMAL) {
+    std::ostringstream sbuf;
+    buf << "Wrong optimal value: the right optimum is " << exp_opt; 
+    check(std::abs(mip.solValue()-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::INTEGER);  
+  expected_opt=1.0/2.0;
+  solveAndCheck(mip, MipSolverBase::OPTIMAL, expected_opt);
+
+
+  //Restrict both to integer
+  mip.colType(x1,MipSolverBase::INTEGER);  
+  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;
+
+}