COIN-OR::LEMON - Graph Library

Ticket #44: aa6dff56b2e1.patch

File aa6dff56b2e1.patch, 186.4 KB (added by Balazs Dezso, 10 years ago)

Under development version of lp port

  • Makefile.am

    # HG changeset patch
    # User Balazs Dezso <deba@inf.elte.hu>
    # Date 1226695561 -3600
    # Node ID aa6dff56b2e1e658d712cd8685023928673d73d4
    # Parent  cbe3ec2d59d2d7c539ab62f2ca6796c960795237
    Development version of lp interface
    
     - Redesigned class structure
     - Redesigned iterattors
     - Some functions in the basic interface redesigned
     - More complete setting functions
     - Ray retrieving functions
     - And lot of improvements
     - Cplex common env
    
    diff -r cbe3ec2d59d2 -r aa6dff56b2e1 Makefile.am
    a b  
    55
    66EXTRA_DIST = \
    77        LICENSE \
     8        m4/lx_check_clp.m4 \
    89        m4/lx_check_cplex.m4 \
    910        m4/lx_check_glpk.m4 \
    1011        m4/lx_check_soplex.m4 \
  • configure.ac

    diff -r cbe3ec2d59d2 -r aa6dff56b2e1 configure.ac
    a b  
    3737fi
    3838
    3939dnl Checks for libraries.
    40 #LX_CHECK_GLPK
    41 #LX_CHECK_CPLEX
    42 #LX_CHECK_SOPLEX
     40LX_CHECK_GLPK
     41LX_CHECK_CPLEX
     42LX_CHECK_SOPLEX
     43LX_CHECK_CLP
    4344
    4445dnl Disable/enable building the demo programs.
    4546AC_ARG_ENABLE([demo],
     
    114115echo C++ compiler.................. : $CXX
    115116echo C++ compiles flags............ : $CXXFLAGS
    116117echo
    117 #echo GLPK support.................. : $lx_glpk_found
    118 #echo CPLEX support................. : $lx_cplex_found
    119 #echo SOPLEX support................ : $lx_soplex_found
    120 #echo
     118echo GLPK support.................. : $lx_glpk_found
     119echo CPLEX support................. : $lx_cplex_found
     120echo SOPLEX support................ : $lx_soplex_found
     121echo CLP support................... : $lx_clp_found
     122echo
    121123echo Build benchmarks.............. : $enable_benchmark
    122124echo Build demo programs........... : $enable_demo
    123125echo Build additional tools........ : $enable_tools
  • lemon/Makefile.am

    diff -r cbe3ec2d59d2 -r aa6dff56b2e1 lemon/Makefile.am
    a b  
    1010        lemon/arg_parser.cc \
    1111        lemon/base.cc \
    1212        lemon/color.cc \
     13        lemon/lp_base.cc \
     14        lemon/lp_skeleton.cc \
    1315        lemon/random.cc
    1416
    15 #lemon_libemon_la_CXXFLAGS = $(GLPK_CFLAGS) $(CPLEX_CFLAGS) $(SOPLEX_CXXFLAGS)
    16 #lemon_libemon_la_LDFLAGS = $(GLPK_LIBS) $(CPLEX_LIBS) $(SOPLEX_LIBS)
     17
     18lemon_libemon_la_CXXFLAGS = \
     19        $(GLPK_CFLAGS) \
     20        $(CPLEX_CFLAGS) \
     21        $(SOPLEX_CXXFLAGS) \
     22        $(CLP_CXXFLAGS)
     23
     24lemon_libemon_la_LDFLAGS = \
     25        $(GLPK_LIBS) \
     26        $(CPLEX_LIBS) \
     27        $(SOPLEX_LIBS) \
     28        $(CLP_LIBS)
     29
     30if HAVE_GLPK
     31lemon_libemon_la_SOURCES += lemon/lp_glpk.cc
     32endif
     33
     34if HAVE_CPLEX
     35lemon_libemon_la_SOURCES += lemon/lp_cplex.cc
     36endif
     37
     38if HAVE_SOPLEX
     39lemon_libemon_la_SOURCES += lemon/lp_soplex.cc
     40endif
     41
     42if HAVE_CLP
     43lemon_libemon_la_SOURCES += lemon/lp_clp.cc
     44endif
    1745
    1846lemon_HEADERS += \
    1947        lemon/arg_parser.h \
     
    3260        lemon/kruskal.h \
    3361        lemon/lgf_reader.h \
    3462        lemon/lgf_writer.h \
     63        lemon/lp_base.h \
     64        lemon/lp_cplex.h \
     65        lemon/lp_glpk.h \
    3566        lemon/list_graph.h \
    3667        lemon/maps.h \
    3768        lemon/math.h \
     
    5283        lemon/bits/graph_extender.h \
    5384        lemon/bits/map_extender.h \
    5485        lemon/bits/path_dump.h \
     86        lemon/bits/solver_bits.h \
    5587        lemon/bits/traits.h \
    5688        lemon/bits/vector_map.h
    5789
  • new file lemon/bits/solver_bits.h

    diff -r cbe3ec2d59d2 -r aa6dff56b2e1 lemon/bits/solver_bits.h
    - +  
     1/* -*- C++ -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library
     4 *
     5 * Copyright (C) 2003-2008
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#ifndef LEMON_BITS_SOLVER_BITS_H
     20#define LEMON_BITS_SOLVER_BITS_H
     21
     22namespace lemon {
     23 
     24  namespace _solver_bits {
     25
     26    class VarIndex {
     27    private:
     28      struct ItemT {
     29        int prev, next;
     30        int index;
     31      };
     32      std::vector<ItemT> items;
     33      int first_item, last_item, first_free_item;
     34     
     35      std::vector<int> cross;
     36
     37    public:
     38
     39      VarIndex()
     40        : first_item(-1), last_item(-1), first_free_item(-1) {
     41      }
     42
     43      void clear() {
     44        first_item = -1;
     45        first_free_item = -1;
     46        items.clear();
     47        cross.clear();
     48      }
     49
     50      int addIndex(int idx) {
     51        int n;
     52        if (first_free_item == -1) {
     53          n = items.size();
     54          items.push_back(ItemT());
     55        } else {
     56          n = first_free_item;
     57          first_free_item = items[n].next;
     58          if (first_free_item != -1) {
     59            items[first_free_item].prev = -1;
     60          }
     61        }
     62        items[n].index = idx;
     63        if (static_cast<int>(cross.size()) <= idx) {
     64          cross.resize(idx + 1, -1);
     65        }
     66        cross[idx] = n;
     67
     68        items[n].prev = last_item;
     69        items[n].next = -1;
     70        if (last_item != -1) {
     71          items[last_item].next = n;
     72        } else {
     73          first_item = n;
     74        }
     75        last_item = n;
     76       
     77        return n;
     78      }
     79
     80      int addIndex(int idx, int n) {
     81        while (n >= static_cast<int>(items.size())) {
     82          items.push_back(ItemT());
     83          items.back().prev = -1;
     84          items.back().next = first_free_item;
     85          if (first_free_item != -1) {
     86            items[first_free_item].prev = items.size() - 1;
     87          }
     88          first_free_item = items.size() - 1;
     89        }
     90        if (items[n].next != -1) {
     91          items[items[n].next].prev = items[n].prev;
     92        }
     93        if (items[n].prev != -1) {
     94          items[items[n].prev].next = items[n].next;
     95        } else {
     96          first_free_item = items[n].next;
     97        }
     98
     99        items[n].index = idx;
     100        if (static_cast<int>(cross.size()) <= idx) {
     101          cross.resize(idx + 1, -1);
     102        }
     103        cross[idx] = n;
     104
     105        items[n].prev = last_item;
     106        items[n].next = -1;
     107        if (last_item != -1) {
     108          items[last_item].next = n;
     109        } else {
     110          first_item = n;
     111        }
     112        last_item = n;
     113       
     114        return n;
     115      }
     116
     117      void eraseIndex(int idx) {
     118        int n = cross[idx];
     119       
     120        if (items[n].prev != -1) {
     121          items[items[n].prev].next = items[n].next;
     122        } else {
     123          first_item = items[n].next;
     124        }
     125        if (items[n].next != -1) {
     126          items[items[n].next].prev = items[n].prev;
     127        } else {
     128          last_item = items[n].prev;
     129        }
     130
     131        if (first_free_item != -1) {
     132          items[first_free_item].prev = n;
     133        }
     134        items[n].next = first_free_item;
     135        items[n].prev = -1;
     136        first_free_item = n;
     137
     138        while (!cross.empty() && cross.back() == -1) {
     139          cross.pop_back();
     140        }
     141      }
     142
     143      int maxIndex() const {
     144        return cross.size() - 1;
     145      }
     146     
     147      void shiftIndexes(int idx) {
     148        for (int i = idx + 1; i < static_cast<int>(cross.size()); ++i) {
     149          cross[i - 1] = cross[i];
     150          if (cross[i] != -1) {
     151            --items[cross[i]].index;
     152          }
     153        }
     154        cross.back() = -1;
     155        cross.pop_back();
     156        while (!cross.empty() && cross.back() == -1) {
     157          cross.pop_back();
     158        }
     159      }
     160
     161      void relocateIndex(int idx, int jdx) {
     162        cross[idx] = cross[jdx];
     163        items[cross[jdx]].index = idx;
     164        cross[jdx] = -1;
     165
     166        while (!cross.empty() && cross.back() == -1) {
     167          cross.pop_back();
     168        }
     169      }
     170
     171      int operator[](int idx) const {
     172        return cross[idx];
     173      }
     174
     175      int operator()(int fdx) const {
     176        return items[fdx].index;
     177      }
     178
     179      void firstItem(int& fdx) const {
     180        fdx = first_item;
     181      }
     182
     183      void nextItem(int& fdx) const {
     184        fdx = items[fdx].next;
     185      }
     186
     187    };
     188  } 
     189}
     190
     191#endif
  • new file lemon/lp.h

    diff -r cbe3ec2d59d2 -r aa6dff56b2e1 lemon/lp.h
    - +  
     1/* -*- C++ -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library
     4 *
     5 * Copyright (C) 2003-2008
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#ifndef LEMON_LP_H
     20#define LEMON_LP_H
     21
     22#include<lemon/config.h>
     23
     24
     25#ifdef HAVE_GLPK
     26#include <lemon/lp_glpk.h>
     27#elif HAVE_CPLEX
     28#include <lemon/lp_cplex.h>
     29#elif HAVE_SOPLEX
     30#include <lemon/lp_soplex.h>
     31#endif
     32
     33///\file
     34///\brief Defines a default LP solver
     35///\ingroup lp_group
     36namespace lemon {
     37
     38#ifdef DOXYGEN
     39  ///The default LP solver identifier
     40
     41  ///The default LP solver identifier.
     42  ///\ingroup lp_group
     43  ///
     44  ///Currently, the possible values are \c GLPK or \c CPLEX
     45#define DEFAULT_LP SOLVER
     46  ///The default LP solver
     47
     48  ///The default LP solver.
     49  ///\ingroup lp_group
     50  ///
     51  ///Currently, it is either \c LpGlpk or \c LpCplex
     52  typedef LpGlpk Lp;
     53  ///The default LP solver identifier string
     54
     55  ///The default LP solver identifier string.
     56  ///\ingroup lp_group
     57  ///
     58  ///Currently, the possible values are "GLPK" or "CPLEX"
     59  const char default_solver_name[]="SOLVER"; 
     60
     61  ///The default ILP solver.
     62
     63  ///The default ILP solver.
     64  ///\ingroup lp_group
     65  ///
     66  ///Currently, it is either \c LpGlpk or \c LpCplex
     67  typedef MipGlpk Mip;
     68#else
     69#ifdef HAVE_GLPK
     70#define DEFAULT_LP GLPK
     71  typedef LpGlpk Lp;
     72  typedef MipGlpk Mip;
     73  const char default_solver_name[]="GLPK";
     74#elif HAVE_CPLEX
     75#define DEFAULT_LP CPLEX
     76  typedef LpCplex Lp;
     77  typedef MipCplex Mip;
     78  const char default_solver_name[]="CPLEX";
     79#elif HAVE_SOPLEX
     80#define DEFAULT_LP SOPLEX
     81  typedef LpSoplex Lp;
     82  const char default_solver_name[]="SOPLEX";
     83#endif
     84#endif
     85 
     86} //namespace lemon
     87
     88#endif //LEMON_LP_H
  • new file lemon/lp_base.cc

    diff -r cbe3ec2d59d2 -r aa6dff56b2e1 lemon/lp_base.cc
    - +  
     1/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library.
     4 *
     5 * Copyright (C) 2003-2008
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19///\file
     20///\brief The implementation of the LP solver interface.
     21
     22#include <lemon/lp_base.h>
     23namespace lemon {
     24
     25  const LpSolverBase::Value
     26  SolverBase::INF = std::numeric_limits<Value>::infinity();
     27  const LpSolverBase::Value
     28  SolverBase::NaN = std::numeric_limits<Value>::quiet_NaN();
     29
     30//   const LpSolverBase::Constr::Value
     31//   LpSolverBase::Constr::INF = std::numeric_limits<Value>::infinity();
     32//   const LpSolverBase::Constr::Value
     33//   LpSolverBase::Constr::NaN = std::numeric_limits<Value>::quiet_NaN();
     34
     35} //namespace lemon
  • new file lemon/lp_base.h

    diff -r cbe3ec2d59d2 -r aa6dff56b2e1 lemon/lp_base.h
    - +  
     1/* -*- C++ -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library
     4 *
     5 * Copyright (C) 2003-2008
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#ifndef LEMON_LP_BASE_H
     20#define LEMON_LP_BASE_H
     21
     22#include<iostream>
     23#include<vector>
     24#include<map>
     25#include<limits>
     26#include<lemon/math.h>
     27
     28#include<lemon/error.h>
     29#include<lemon/assert.h>
     30
     31#include<lemon/core.h>
     32#include<lemon/bits/solver_bits.h>
     33
     34///\file
     35///\brief The interface of the LP solver interface.
     36///\ingroup lp_group
     37namespace lemon {
     38
     39  ///Common base class for LP solvers
     40 
     41  ///\todo Much more docs
     42  ///\ingroup lp_group
     43  class SolverBase {
     44
     45  protected:
     46
     47    _solver_bits::VarIndex rows;
     48    _solver_bits::VarIndex cols;
     49   
     50  public:
     51   
     52    ///Possible outcomes of an LP solving procedure
     53    enum SolveExitStatus {
     54      ///This means that the problem has been successfully solved: either
     55      ///an optimal solution has been found or infeasibility/unboundedness
     56      ///has been proved.
     57      SOLVED = 0,
     58      ///Any other case (including the case when some user specified
     59      ///limit has been exceeded)
     60      UNSOLVED = 1
     61    };
     62
     63    ///Direction of the optimization
     64    enum Sense {
     65      /// Minimization
     66      MIN,
     67      /// Maximization
     68      MAX
     69    };
     70     
     71    ///The floating point type used by the solver
     72    typedef double Value;
     73    ///The infinity constant
     74    static const Value INF;
     75    ///The not a number constant
     76    static const Value NaN;
     77   
     78    friend class Col;
     79    friend class ColIt;
     80    friend class Row;
     81    friend class RowIt;
     82   
     83    ///Refer to a column of the LP.
     84
     85    ///This type is used to refer to a column of the LP.
     86    ///
     87    ///Its value remains valid and correct even after the addition or erase of
     88    ///other columns.
     89    ///
     90    ///\todo Document what can one do with a Col (INVALID, comparing,
     91    ///it is similar to Node/Edge)
     92    class Col {
     93      friend class SolverBase;
     94    protected:
     95      int _id;
     96      explicit Col(int id) : _id(id) {}
     97    public:
     98      typedef Value ExprValue;
     99      typedef True LpSolverCol;
     100      Col() {}
     101      Col(const Invalid&) : _id(-1) {}
     102      bool operator<(Col c) const  {return _id < c._id;}
     103      bool operator==(Col c) const  {return _id == c._id;}
     104      bool operator!=(Col c) const  {return _id != c._id;}
     105    };
     106
     107    class ColIt : public Col {
     108      const SolverBase *_solver;
     109    public:
     110      ColIt() {}
     111      ColIt(const SolverBase &solver) : _solver(&solver)
     112      {
     113        _solver->cols.firstItem(_id);
     114      }
     115      ColIt(const Invalid&) : Col(INVALID) {}
     116      ColIt &operator++()
     117      {
     118        _solver->cols.nextItem(_id);
     119        return *this;
     120      }
     121    };
     122
     123    static int id(const Col& col) { return col._id; }
     124    static Col colFromId(int id) { return Col(id); }
     125     
     126    ///Refer to a row of the LP.
     127
     128    ///This type is used to refer to a row of the LP.
     129    ///
     130    ///Its value remains valid and correct even after the addition or erase of
     131    ///other rows.
     132    ///
     133    ///\todo Document what can one do with a Row (INVALID, comparing,
     134    ///it is similar to Node/Edge)
     135    class Row {
     136      friend class SolverBase;
     137    protected:
     138      int _id;
     139      explicit Row(int id) : _id(id) {}
     140    public:
     141      typedef Value ExprValue;
     142      typedef True LpSolverRow;
     143      Row() {}
     144      Row(const Invalid&) : _id(-1) {}
     145
     146      bool operator<(Row c) const  {return _id < c._id;}
     147      bool operator==(Row c) const  {return _id == c._id;}
     148      bool operator!=(Row c) const  {return _id != c._id;}
     149    };
     150
     151    class RowIt : public Row {
     152      const SolverBase *_solver;
     153    public:
     154      RowIt() {}
     155      RowIt(const SolverBase &solver) : _solver(&solver)
     156      {
     157        _solver->rows.firstItem(_id);
     158      }
     159      RowIt(const Invalid&) : Row(INVALID) {}
     160      RowIt &operator++()
     161      {
     162        _solver->rows.nextItem(_id);
     163        return *this;
     164      }
     165    };
     166
     167    static int id(const Row& row) { return row._id; }
     168    static Row rowFromId(int id) { return Row(id); }
     169
     170  public:
     171   
     172    ///Linear expression of variables and a constant component
     173   
     174    ///This data structure stores a linear expression of the variables
     175    ///(\ref Col "Col"s) and also has a constant component.
     176    ///
     177    ///There are several ways to access and modify the contents of this
     178    ///container.
     179    ///\code
     180    ///e[v]=5;
     181    ///e[v]+=12;
     182    ///e.erase(v);
     183    ///\endcode
     184    ///or you can also iterate through its elements.
     185    ///\code
     186    ///double s=0;
     187    ///for(SolverBase::Expr::ConstMapIt i(e);i!=INVALID;++i)
     188    ///  s+=*i * primal(i);
     189    ///\endcode
     190    ///(This code computes the primal value of the expression).
     191    ///- Numbers (<tt>double</tt>'s)
     192    ///and variables (\ref Col "Col"s) directly convert to an
     193    ///\ref Expr and the usual linear operations are defined, so 
     194    ///\code
     195    ///v+w
     196    ///2*v-3.12*(v-w/2)+2
     197    ///v*2.1+(3*v+(v*12+w+6)*3)/2
     198    ///\endcode
     199    ///are valid \ref Expr "Expr"essions.
     200    ///The usual assignment operations are also defined.
     201    ///\code
     202    ///e=v+w;
     203    ///e+=2*v-3.12*(v-w/2)+2;
     204    ///e*=3.4;
     205    ///e/=5;
     206    ///\endcode
     207    ///- The constant member can be set and read by \ref operator*()
     208    ///\code
     209    ///*e=12;
     210    ///double c=*e;
     211    ///\endcode
     212    ///
     213    ///\note \ref clear() not only sets all coefficients to 0 but also
     214    ///clears the constant components.
     215    ///
     216    ///\sa Constr
     217    ///
     218    class Expr {
     219      friend class SolverBase;
     220    public:
     221      typedef SolverBase::Col Key;
     222      typedef SolverBase::Value Value;
     223     
     224    protected:
     225      Value const_comp;
     226      std::map<int, Value> comps;
     227
     228    public:
     229      typedef True SolverExpr;
     230      ///\e
     231      Expr() : const_comp(0) {}
     232      ///\e
     233      Expr(const Key &k) : const_comp(0) {
     234        typedef std::map<int, Value>::value_type pair_type;
     235        comps.insert(pair_type(id(k), 1));
     236      }
     237      ///\e
     238      Expr(const Value &v) : const_comp(v) {}
     239      ///\e
     240      Value operator[](const Key& k) const {
     241        std::map<int, Value>::const_iterator it=comps.find(id(k));
     242        if (it != comps.end()) {
     243          return it->second;
     244        } else {
     245          return 0;
     246        }
     247      }
     248      ///\e
     249      Value& operator[](const Key& k) {
     250        return comps[id(k)];
     251      }
     252      ///\e
     253      void set(const Key &k, const Value &v) {
     254        if (v != 0.0) {
     255          typedef std::map<int, Value>::value_type pair_type;
     256          comps.insert(pair_type(id(k), v));
     257        } else {
     258          comps.erase(id(k));
     259        }
     260      }
     261
     262      ///\e
     263      Value& operator*() { return const_comp; }
     264      ///\e
     265      const Value& operator*() const { return const_comp; }
     266     
     267      ///Removes the coefficients closer to zero than \c tolerance.
     268      void simplify(Value tolerance = 0.0) {
     269        std::map<int, Value>::iterator it=comps.begin();
     270        while (it != comps.end()) {
     271          std::map<int, Value>::iterator jt=it;
     272          ++jt;
     273          if (std::fabs((*it).second) <= tolerance) comps.erase(it);
     274          it=jt;
     275        }
     276        if (std::fabs(const_comp) <= tolerance) const_comp = 0;
     277      }
     278
     279      void simplify(Value tolerance = 0.0) const {
     280        const_cast<Expr*>(this)->simplify(tolerance);
     281      }
     282
     283      ///Sets all coefficients and the constant component to 0.
     284      void clear() {
     285        comps.clear();
     286        const_comp=0;
     287      }
     288
     289      ///\e
     290      Expr &operator+=(const Expr &e) {
     291        for (std::map<int, Value>::const_iterator it=e.comps.begin();
     292             it!=e.comps.end(); ++it)
     293          comps[it->first]+=it->second;
     294        const_comp+=e.const_comp;
     295        return *this;
     296      }
     297      ///\e
     298      Expr &operator-=(const Expr &e) {
     299        for (std::map<int, Value>::const_iterator it=e.comps.begin();
     300             it!=e.comps.end(); ++it)
     301          comps[it->first]-=it->second;
     302        const_comp-=e.const_comp;
     303        return *this;
     304      }
     305      ///\e
     306      Expr &operator*=(const Value &c) {
     307        for (std::map<int, Value>::iterator it=comps.begin();
     308             it!=comps.end(); ++it)
     309          it->second*=c;
     310        const_comp*=c;
     311        return *this;
     312      }
     313      ///\e
     314      Expr &operator/=(const Value &c) {
     315        for (std::map<int, Value>::iterator it=comps.begin();
     316             it!=comps.end(); ++it)
     317          it->second/=c;
     318        const_comp/=c;
     319        return *this;
     320      }
     321
     322      class MapIt {
     323      private:
     324       
     325        std::map<int, Value>::iterator _it, _end;
     326       
     327      public:
     328
     329        MapIt(Expr& e)
     330          : _it(e.comps.begin()), _end(e.comps.end()){}
     331
     332        operator Col() const {
     333          return colFromId(_it->first);
     334        }
     335       
     336        Value& operator*() { return _it->second; }
     337        const Value& operator*() const { return _it->second; }
     338
     339        MapIt& operator++() { ++_it; return *this; }
     340
     341        bool operator==(Invalid) const { return _it == _end; }
     342        bool operator!=(Invalid) const { return _it != _end; }
     343      };
     344
     345      class ConstMapIt {
     346      private:
     347       
     348        std::map<int, Value>::const_iterator _it, _end;
     349       
     350      public:
     351
     352        ConstMapIt(const Expr& e)
     353          : _it(e.comps.begin()), _end(e.comps.end()){}
     354
     355        operator Col() const {
     356          return colFromId(_it->first);
     357        }
     358       
     359        const Value& operator*() const { return _it->second; }
     360
     361        ConstMapIt& operator++() { ++_it; return *this; }
     362
     363        bool operator==(Invalid) const { return _it == _end; }
     364        bool operator!=(Invalid) const { return _it != _end; }
     365      };
     366
     367    };
     368   
     369    ///Linear constraint
     370
     371    ///This data stucture represents a linear constraint in the LP.
     372    ///Basically it is a linear expression with a lower or an upper bound
     373    ///(or both). These parts of the constraint can be obtained by the member
     374    ///functions \ref expr(), \ref lowerBound() and \ref upperBound(),
     375    ///respectively.
     376    ///There are two ways to construct a constraint.
     377    ///- You can set the linear expression and the bounds directly
     378    ///  by the functions above.
     379    ///- The operators <tt>\<=</tt>, <tt>==</tt> and  <tt>\>=</tt>
     380    ///  are defined between expressions, or even between constraints whenever
     381    ///  it makes sense. Therefore if \c e and \c f are linear expressions and
     382    ///  \c s and \c t are numbers, then the followings are valid expressions
     383    ///  and thus they can be used directly e.g. in \ref addRow() whenever
     384    ///  it makes sense.
     385    ///\code
     386    ///  e<=s
     387    ///  e<=f
     388    ///  e==f
     389    ///  s<=e<=t
     390    ///  e>=t
     391    ///\endcode
     392    ///\warning The validity of a constraint is checked only at run
     393    ///time, so e.g. \ref addRow(<tt>x[1]\<=x[2]<=5</tt>) will
     394    ///compile, but will fail an assertion.
     395    class Constr
     396    {
     397    public:
     398      typedef SolverBase::Expr Expr;
     399      typedef Expr::Key Key;
     400      typedef Expr::Value Value;
     401     
     402    protected:
     403      Expr _expr;
     404      Value _lb,_ub;
     405    public:
     406      ///\e
     407      Constr() : _expr(), _lb(NaN), _ub(NaN) {}
     408      ///\e
     409      Constr(Value lb,const Expr &e,Value ub) :
     410        _expr(e), _lb(lb), _ub(ub) {}
     411      Constr(const Expr &e) :
     412        _expr(e), _lb(NaN), _ub(NaN) {}
     413      ///\e
     414      void clear()
     415      {
     416        _expr.clear();
     417        _lb=_ub=NaN;
     418      }
     419
     420      ///Reference to the linear expression
     421      Expr &expr() { return _expr; }
     422      ///Cont reference to the linear expression
     423      const Expr &expr() const { return _expr; }
     424      ///Reference to the lower bound.
     425
     426      ///\return
     427      ///- \ref INF "INF": the constraint is lower unbounded.
     428      ///- \ref NaN "NaN": lower bound has not been set.
     429      ///- finite number: the lower bound
     430      Value &lowerBound() { return _lb; }
     431      ///The const version of \ref lowerBound()
     432      const Value &lowerBound() const { return _lb; }
     433      ///Reference to the upper bound.
     434
     435      ///\return
     436      ///- \ref INF "INF": the constraint is upper unbounded.
     437      ///- \ref NaN "NaN": upper bound has not been set.
     438      ///- finite number: the upper bound
     439      Value &upperBound() { return _ub; }
     440      ///The const version of \ref upperBound()
     441      const Value &upperBound() const { return _ub; }
     442      ///Is the constraint lower bounded?
     443      bool lowerBounded() const {
     444        return _lb != -INF && !std::isnan(_lb);
     445      }
     446      ///Is the constraint upper bounded?
     447      bool upperBounded() const {
     448        return _ub != INF && !std::isnan(_ub);
     449      }
     450
     451    };
     452   
     453    ///Linear expression of rows
     454   
     455    ///This data structure represents a column of the matrix,
     456    ///thas is it strores a linear expression of the dual variables
     457    ///(\ref Row "Row"s).
     458    ///
     459    ///There are several ways to access and modify the contents of this
     460    ///container.
     461    ///\code
     462    ///e[v]=5;
     463    ///e[v]+=12;
     464    ///e.erase(v);
     465    ///\endcode
     466    ///or you can also iterate through its elements.
     467    ///\code
     468    ///double s=0;
     469    ///for(SolverBase::DualExpr::ConstMapIt i(e);i!=INVALID;++i)
     470    ///  s+=*i;
     471    ///\endcode
     472    ///(This code computes the sum of all coefficients).
     473    ///- Numbers (<tt>double</tt>'s)
     474    ///and variables (\ref Row "Row"s) directly convert to an
     475    ///\ref DualExpr and the usual linear operations are defined, so
     476    ///\code
     477    ///v+w
     478    ///2*v-3.12*(v-w/2)
     479    ///v*2.1+(3*v+(v*12+w)*3)/2
     480    ///\endcode
     481    ///are valid \ref DualExpr "DualExpr"essions.
     482    ///The usual assignment operations are also defined.
     483    ///\code
     484    ///e=v+w;
     485    ///e+=2*v-3.12*(v-w/2);
     486    ///e*=3.4;
     487    ///e/=5;
     488    ///\endcode
     489    ///
     490    ///\sa Expr
     491    ///
     492    class DualExpr {
     493      friend class SolverBase;
     494    public:
     495      typedef SolverBase::Row Key;
     496      typedef SolverBase::Value Value;
     497     
     498    protected:
     499      std::map<int, Value> comps;
     500     
     501    public:
     502      typedef True SolverExpr;
     503      ///\e
     504      DualExpr() {}
     505      ///\e
     506      DualExpr(const Key &k) {
     507        typedef std::map<int, Value>::value_type pair_type;
     508        comps.insert(pair_type(id(k), 1));
     509      }
     510      ///\e
     511      Value operator[](const Key& k) const {
     512        std::map<int, Value>::const_iterator it = comps.find(id(k));
     513        if (it != comps.end()) {
     514          return it->second;
     515        } else {
     516          return 0;
     517        }
     518      }
     519      ///\e
     520      Value& operator[](const Key& k) {
     521        return comps[id(k)];
     522      }
     523      ///\e
     524      void set(const Key &k, const Value &v) {
     525        if (v != 0.0) {
     526          typedef std::map<int, Value>::value_type pair_type;
     527          comps.insert(pair_type(id(k), v));
     528        } else {
     529          comps.erase(id(k));
     530        }
     531      }
     532     
     533      ///Removes the coefficients closer to zero than \c tolerance.
     534      void simplify(Value tolerance = 0.0) {
     535        std::map<int, Value>::iterator it=comps.begin();
     536        while (it != comps.end()) {
     537          std::map<int, Value>::iterator jt=it;
     538          ++jt;
     539          if (std::fabs((*it).second) <= tolerance) comps.erase(it);
     540          it=jt;
     541        }
     542      }
     543
     544      void simplify(Value tolerance = 0.0) const {
     545        const_cast<DualExpr*>(this)->simplify(tolerance);
     546      }
     547
     548      ///Sets all coefficients to 0.
     549      void clear() {
     550        comps.clear();
     551      }
     552
     553      ///\e
     554      DualExpr &operator+=(const DualExpr &e) {
     555        for (std::map<int, Value>::const_iterator it=e.comps.begin();
     556             it!=e.comps.end(); ++it)
     557          comps[it->first]+=it->second;
     558        return *this;
     559      }
     560      ///\e
     561      DualExpr &operator-=(const DualExpr &e) {
     562        for (std::map<int, Value>::const_iterator it=e.comps.begin();
     563             it!=e.comps.end(); ++it)
     564          comps[it->first]-=it->second;
     565        return *this;
     566      }
     567      ///\e
     568      DualExpr &operator*=(const Value &c) {
     569        for (std::map<int, Value>::iterator it=comps.begin();
     570             it!=comps.end(); ++it)
     571          it->second*=c;
     572        return *this;
     573      }
     574      ///\e
     575      DualExpr &operator/=(const Value &c) {
     576        for (std::map<int, Value>::iterator it=comps.begin();
     577             it!=comps.end(); ++it)
     578          it->second/=c;
     579        return *this;
     580      }
     581
     582      class MapIt {
     583      private:
     584       
     585        std::map<int, Value>::iterator _it, _end;
     586       
     587      public:
     588
     589        MapIt(DualExpr& e)
     590          : _it(e.comps.begin()), _end(e.comps.end()){}
     591
     592        operator Row() const {
     593          return rowFromId(_it->first);
     594        }
     595       
     596        Value& operator*() { return _it->second; }
     597        const Value& operator*() const { return _it->second; }
     598
     599        MapIt& operator++() { ++_it; return *this; }
     600
     601        bool operator==(Invalid) const { return _it == _end; }
     602        bool operator!=(Invalid) const { return _it != _end; }
     603      };
     604
     605      class ConstMapIt {
     606      private:
     607       
     608        std::map<int, Value>::const_iterator _it, _end;
     609       
     610      public:
     611
     612        ConstMapIt(const DualExpr& e)
     613          : _it(e.comps.begin()), _end(e.comps.end()){}
     614
     615        operator Row() const {
     616          return rowFromId(_it->first);
     617        }
     618       
     619        const Value& operator*() const { return _it->second; }
     620
     621        ConstMapIt& operator++() { ++_it; return *this; }
     622
     623        bool operator==(Invalid) const { return _it == _end; }
     624        bool operator!=(Invalid) const { return _it != _end; }
     625      };
     626    };
     627   
     628
     629  protected:
     630
     631    class InsertIterator {
     632    private:
     633     
     634      std::map<int, Value>& _host;
     635      const _solver_bits::VarIndex& _index;
     636
     637    public:
     638
     639      typedef std::output_iterator_tag iterator_category;
     640      typedef void difference_type;
     641      typedef void value_type;
     642      typedef void reference;
     643      typedef void pointer;
     644     
     645      InsertIterator(std::map<int, Value>& host,
     646                   const _solver_bits::VarIndex& index)
     647        : _host(host), _index(index) {}
     648
     649      InsertIterator& operator=(const std::pair<int, Value>& value) {
     650        typedef std::map<int, Value>::value_type pair_type;
     651        _host.insert(pair_type(_index[value.first], value.second));
     652        return *this;
     653      }
     654
     655      InsertIterator& operator*() { return *this; }
     656      InsertIterator& operator++() { return *this; }
     657      InsertIterator operator++(int) { return *this; }
     658
     659    };
     660
     661    class ExprIterator {
     662    private:
     663      std::map<int, Value>::const_iterator _host_it;
     664      const _solver_bits::VarIndex& _index;
     665    public:
     666
     667      typedef std::bidirectional_iterator_tag iterator_category;
     668      typedef std::ptrdiff_t difference_type;
     669      typedef const std::pair<int, Value> value_type;
     670      typedef value_type reference;
     671
     672      class pointer {
     673      public:
     674        pointer(value_type& _value) : value(_value) {}
     675        value_type* operator->() { return &value; }
     676      private:
     677        value_type value;
     678      };
     679
     680      ExprIterator(const std::map<int, Value>::const_iterator& host_it,
     681                   const _solver_bits::VarIndex& index)
     682        : _host_it(host_it), _index(index) {}
     683
     684      reference operator*() {
     685        return std::make_pair(_index(_host_it->first), _host_it->second);
     686      }
     687
     688      pointer operator->() {
     689        return pointer(operator*());
     690      }
     691
     692      ExprIterator& operator++() { ++_host_it; return *this; }
     693      ExprIterator operator++(int) {
     694        ExprIterator tmp(*this); ++_host_it; return tmp;
     695      }
     696
     697      ExprIterator& operator--() { --_host_it; return *this; }
     698      ExprIterator operator--(int) {
     699        ExprIterator tmp(*this); --_host_it; return tmp;
     700      }
     701
     702      bool operator==(const ExprIterator& it) const {
     703        return _host_it == it._host_it;
     704      }
     705
     706      bool operator!=(const ExprIterator& it) const {
     707        return _host_it != it._host_it;
     708      }
     709
     710    };
     711
     712  protected:
     713
     714    //Abstract virtual functions
     715    virtual SolverBase* _newSolver() const = 0;
     716    virtual SolverBase* _cloneSolver() const = 0;
     717
     718    virtual int _addColId(int col) { return cols.addIndex(col); }
     719    virtual int _addRowId(int row) { return rows.addIndex(row); }
     720
     721    virtual void _eraseColId(int col) { cols.eraseIndex(col); }
     722    virtual void _eraseRowId(int row) { rows.eraseIndex(row); }
     723
     724    virtual int _addCol() = 0;
     725    virtual int _addRow() = 0;
     726
     727    virtual void _eraseCol(int col) = 0;
     728    virtual void _eraseRow(int row) = 0;
     729
     730    virtual void _getColName(int col, std::string& name) const = 0;
     731    virtual void _setColName(int col, const std::string& name) = 0;
     732    virtual int _colByName(const std::string& name) const = 0;
     733
     734    virtual void _getRowName(int row, std::string& name) const = 0;
     735    virtual void _setRowName(int row, const std::string& name) = 0;
     736    virtual int _rowByName(const std::string& name) const = 0;
     737
     738    virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e) = 0;
     739    virtual void _getRowCoeffs(int i, InsertIterator b) const = 0;
     740
     741    virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e) = 0;
     742    virtual void _getColCoeffs(int i, InsertIterator b) const = 0;
     743
     744    virtual void _setCoeff(int row, int col, Value value) = 0;
     745    virtual Value _getCoeff(int row, int col) const = 0;
     746
     747    virtual void _setColLowerBound(int i, Value value) = 0;
     748    virtual Value _getColLowerBound(int i) const = 0;
     749
     750    virtual void _setColUpperBound(int i, Value value) = 0;
     751    virtual Value _getColUpperBound(int i) const = 0;
     752
     753    virtual void _setRowLowerBound(int i, Value value) = 0;
     754    virtual Value _getRowLowerBound(int i) const = 0;
     755
     756    virtual void _setRowUpperBound(int i, Value value) = 0;
     757    virtual Value _getRowUpperBound(int i) const = 0;
     758
     759    virtual void _setObjCoeffs(ExprIterator b, ExprIterator e) = 0;
     760    virtual void _getObjCoeffs(InsertIterator b) const = 0;
     761
     762    virtual void _setObjCoeff(int i, Value obj_coef) = 0;
     763    virtual Value _getObjCoeff(int i) const = 0;
     764
     765    virtual void _setSense(Sense) = 0;
     766    virtual Sense _getSense() const = 0;
     767
     768    virtual void _clear() = 0;
     769
     770    //Own protected stuff
     771   
     772    //Constant component of the objective function
     773    Value obj_const_comp;
     774       
     775  public:
     776
     777    ///\e
     778    SolverBase() : rows(), cols(), obj_const_comp(0) {
     779    }
     780
     781    ///\e
     782    virtual ~SolverBase() {}
     783
     784    ///Creates a new LP problem
     785    SolverBase* newSolver() {return _newSolver();}
     786    ///Makes a copy of the LP problem
     787    SolverBase* cloneSolver() {return _cloneSolver();}
     788   
     789    ///\name Build up and modify the LP
     790
     791    ///@{
     792
     793    ///Add a new empty column (i.e a new variable) to the LP
     794    Col addCol() { Col c; c._id = _addColId(_addCol()); return c;}
     795
     796    ///\brief Adds several new columns
     797    ///(i.e a variables) at once
     798    ///
     799    ///This magic function takes a container as its argument
     800    ///and fills its elements
     801    ///with new columns (i.e. variables)
     802    ///\param t can be
     803    ///- a standard STL compatible iterable container with
     804    ///\ref Col as its \c values_type
     805    ///like
     806    ///\code
     807    ///std::vector<SolverBase::Col>
     808    ///std::list<SolverBase::Col>
     809    ///\endcode
     810    ///- a standard STL compatible iterable container with
     811    ///\ref Col as its \c mapped_type
     812    ///like
     813    ///\code
     814    ///std::map<AnyType,SolverBase::Col>
     815    ///\endcode
     816    ///- an iterable lemon \ref concepts::WriteMap "write map" like
     817    ///\code
     818    ///ListGraph::NodeMap<SolverBase::Col>
     819    ///ListGraph::ArcMap<SolverBase::Col>
     820    ///\endcode
     821    ///\return The number of the created column.
     822#ifdef DOXYGEN
     823    template<class T>
     824    int addColSet(T &t) { return 0;}
     825#else
     826    template<class T>
     827    typename enable_if<typename T::value_type::LpSolverCol,int>::type
     828    addColSet(T &t,dummy<0> = 0) {
     829      int s=0;
     830      for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addCol();s++;}
     831      return s;
     832    }
     833    template<class T>
     834    typename enable_if<typename T::value_type::second_type::LpSolverCol,
     835                       int>::type
     836    addColSet(T &t,dummy<1> = 1) {
     837      int s=0;
     838      for(typename T::iterator i=t.begin();i!=t.end();++i) {
     839        i->second=addCol();
     840        s++;
     841      }
     842      return s;
     843    }
     844    template<class T>
     845    typename enable_if<typename T::MapIt::Value::LpSolverCol,
     846                       int>::type
     847    addColSet(T &t,dummy<2> = 2) {
     848      int s=0;
     849      for(typename T::MapIt i(t); i!=INVALID; ++i)
     850        {
     851          i.set(addCol());
     852          s++;
     853        }
     854      return s;
     855    }
     856#endif
     857
     858    ///Set a column (i.e a dual constraint) of the LP
     859
     860    ///\param c is the column to be modified
     861    ///\param e is a dual linear expression (see \ref DualExpr)
     862    ///a better one.
     863    void col(Col c, const DualExpr &e) {
     864      e.simplify();
     865      _setColCoeffs(cols(id(c)), ExprIterator(e.comps.begin(), cols),
     866                    ExprIterator(e.comps.end(), cols));
     867    }
     868
     869    ///Get a column (i.e a dual constraint) of the LP
     870
     871    ///\param r is the column to get
     872    ///\return the dual expression associated to the column
     873    DualExpr col(Col c) const {
     874      DualExpr e;
     875      _getColCoeffs(cols(id(c)), InsertIterator(e.comps, rows));
     876      return e;
     877    }
     878
     879    ///Add a new column to the LP
     880
     881    ///\param e is a dual linear expression (see \ref DualExpr)
     882    ///\param obj is the corresponding component of the objective
     883    ///function. It is 0 by default.
     884    ///\return The created column.
     885    Col addCol(const DualExpr &e, Value o = 0) {
     886      Col c=addCol();
     887      col(c,e);
     888      objCoeff(c,o);
     889      return c;
     890    }
     891
     892    ///Add a new empty row (i.e a new constraint) to the LP
     893
     894    ///This function adds a new empty row (i.e a new constraint) to the LP.
     895    ///\return The created row
     896    Row addRow() { Row r; r._id = _addRowId(_addRow()); return r;}
     897
     898    ///\brief Add several new rows
     899    ///(i.e a constraints) at once
     900    ///
     901    ///This magic function takes a container as its argument
     902    ///and fills its elements
     903    ///with new row (i.e. variables)
     904    ///\param t can be
     905    ///- a standard STL compatible iterable container with
     906    ///\ref Row as its \c values_type
     907    ///like
     908    ///\code
     909    ///std::vector<SolverBase::Row>
     910    ///std::list<SolverBase::Row>
     911    ///\endcode
     912    ///- a standard STL compatible iterable container with
     913    ///\ref Row as its \c mapped_type
     914    ///like
     915    ///\code
     916    ///std::map<AnyType,SolverBase::Row>
     917    ///\endcode
     918    ///- an iterable lemon \ref concepts::WriteMap "write map" like
     919    ///\code
     920    ///ListGraph::NodeMap<SolverBase::Row>
     921    ///ListGraph::ArcMap<SolverBase::Row>
     922    ///\endcode
     923    ///\return The number of rows created.
     924#ifdef DOXYGEN
     925    template<class T>
     926    int addRowSet(T &t) { return 0;}
     927#else
     928    template<class T>
     929    typename enable_if<typename T::value_type::LpSolverRow,int>::type
     930    addRowSet(T &t, dummy<0> = 0) {
     931      int s=0;
     932      for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addRow();s++;}
     933      return s;
     934    }
     935    template<class T>
     936    typename enable_if<typename T::value_type::second_type::LpSolverRow,
     937                       int>::type
     938    addRowSet(T &t, dummy<1> = 1) {
     939      int s=0;
     940      for(typename T::iterator i=t.begin();i!=t.end();++i) {
     941        i->second=addRow();
     942        s++;
     943      }
     944      return s;
     945    }
     946    template<class T>
     947    typename enable_if<typename T::MapIt::Value::LpSolverRow,
     948                       int>::type
     949    addRowSet(T &t, dummy<2> = 2) {
     950      int s=0;
     951      for(typename T::MapIt i(t); i!=INVALID; ++i)
     952        {
     953          i.set(addRow());
     954          s++;
     955        }
     956      return s;
     957    }
     958#endif
     959
     960    ///Set a row (i.e a constraint) of the LP
     961
     962    ///\param r is the row to be modified
     963    ///\param l is lower bound (-\ref INF means no bound)
     964    ///\param e is a linear expression (see \ref Expr)
     965    ///\param u is the upper bound (\ref INF means no bound)
     966    void row(Row r, Value l, const Expr &e, Value u) {
     967      e.simplify();
     968      _setRowCoeffs(rows(id(r)), ExprIterator(e.comps.begin(), cols),
     969                    ExprIterator(e.comps.end(), cols));
     970      _setRowLowerBound(rows(id(r)),l - *e);
     971      _setRowUpperBound(rows(id(r)),u - *e);
     972    }
     973
     974    ///Set a row (i.e a constraint) of the LP
     975
     976    ///\param r is the row to be modified
     977    ///\param c is a linear expression (see \ref Constr)
     978    void row(Row r, const Constr &c) {
     979      row(r, c.lowerBounded()?c.lowerBound():-INF,
     980          c.expr(), c.upperBounded()?c.upperBound():INF);
     981    }
     982
     983   
     984    ///Get a row (i.e a constraint) of the LP
     985
     986    ///\param r is the row to get
     987    ///\return the expression associated to the row
     988    Expr row(Row r) const {
     989      Expr e;
     990      _getRowCoeffs(rows(id(r)), InsertIterator(e.comps, cols));
     991      return e;
     992    }
     993
     994    ///Add a new row (i.e a new constraint) to the LP
     995
     996    ///\param l is the lower bound (-\ref INF means no bound)
     997    ///\param e is a linear expression (see \ref Expr)
     998    ///\param u is the upper bound (\ref INF means no bound)
     999    ///\return The created row.
     1000    Row addRow(Value l,const Expr &e, Value u) {
     1001      Row r=addRow();
     1002      row(r,l,e,u);
     1003      return r;
     1004    }
     1005
     1006    ///Add a new row (i.e a new constraint) to the LP
     1007
     1008    ///\param c is a linear expression (see \ref Constr)
     1009    ///\return The created row.
     1010    Row addRow(const Constr &c) {
     1011      Row r=addRow();
     1012      row(r,c);
     1013      return r;
     1014    }
     1015    ///Erase a column (i.e a variable) from the LP
     1016
     1017    ///\param c is the column to be deleted
     1018    ///\todo Please check this
     1019    void erase(Col c) {
     1020      _eraseCol(cols(id(c)));
     1021      _eraseColId(cols(id(c)));
     1022    }
     1023    ///Erase a row (i.e a constraint) from the LP
     1024
     1025    ///\param r is the row to be deleted
     1026    ///\todo Please check this
     1027    void erase(Row r) {
     1028      _eraseRow(rows(id(r)));
     1029      _eraseRowId(rows(id(r)));
     1030    }
     1031
     1032    /// Get the name of a column
     1033   
     1034    ///\param c is the coresponding column
     1035    ///\return The name of the colunm
     1036    std::string colName(Col c) const {
     1037      std::string name;
     1038      _getColName(cols(id(c)), name);
     1039      return name;
     1040    }
     1041   
     1042    /// Set the name of a column
     1043   
     1044    ///\param c is the coresponding column
     1045    ///\param name The name to be given
     1046    void colName(Col c, const std::string& name) {
     1047      _setColName(cols(id(c)), name);
     1048    }
     1049
     1050    /// Get the column by its name
     1051   
     1052    ///\param name The name of the column
     1053    ///\return the proper column or \c INVALID
     1054    Col colByName(const std::string& name) const {
     1055      int k = _colByName(name);
     1056      return k != -1 ? Col(cols[k]) : Col(INVALID);
     1057    }
     1058
     1059    /// Get the name of a row
     1060   
     1061    ///\param r is the coresponding row
     1062    ///\return The name of the row
     1063    std::string rowName(Row r) const {
     1064      std::string name;
     1065      _getRowName(rows(id(r)), name);
     1066      return name;
     1067    }
     1068   
     1069    /// Set the name of a row
     1070   
     1071    ///\param r is the coresponding row
     1072    ///\param name The name to be given
     1073    void rowName(Row r, const std::string& name) {
     1074      _setRowName(rows(id(r)), name);
     1075    }
     1076
     1077    /// Get the row by its name
     1078   
     1079    ///\param name The name of the row
     1080    ///\return the proper row or \c INVALID
     1081    Row rowByName(const std::string& name) const {
     1082      int k = _rowByName(name);
     1083      return k != -1 ? Row(rows[k]) : Row(INVALID);
     1084    }
     1085   
     1086    /// Set an element of the coefficient matrix of the LP
     1087
     1088    ///\param r is the row of the element to be modified
     1089    ///\param c is the column of the element to be modified
     1090    ///\param val is the new value of the coefficient
     1091    void coeff(Row r, Col c, Value val) {
     1092      _setCoeff(rows(id(r)),cols(id(c)), val);
     1093    }
     1094
     1095    /// Get an element of the coefficient matrix of the LP
     1096
     1097    ///\param r is the row of the element in question
     1098    ///\param c is the column of the element in question
     1099    ///\return the corresponding coefficient
     1100    Value coeff(Row r, Col c) const {
     1101      return _getCoeff(rows(id(r)),cols(id(c)));
     1102    }
     1103
     1104    /// Set the lower bound of a column (i.e a variable)
     1105
     1106    /// The lower bound of a variable (column) has to be given by an
     1107    /// extended number of type Value, i.e. a finite number of type
     1108    /// Value or -\ref INF.
     1109    void colLowerBound(Col c, Value value) {
     1110      _setColLowerBound(cols(id(c)),value);
     1111    }
     1112
     1113    /// Get the lower bound of a column (i.e a variable)
     1114
     1115    /// This function returns the lower bound for column (variable) \t c
     1116    /// (this might be -\ref INF as well). 
     1117    ///\return The lower bound for column \t c
     1118    Value colLowerBound(Col c) const {
     1119      return _getColLowerBound(cols(id(c)));
     1120    }
     1121   
     1122    ///\brief Set the lower bound of  several columns
     1123    ///(i.e a variables) at once
     1124    ///
     1125    ///This magic function takes a container as its argument
     1126    ///and applies the function on all of its elements.
     1127    /// The lower bound of a variable (column) has to be given by an
     1128    /// extended number of type Value, i.e. a finite number of type
     1129    /// Value or -\ref INF.
     1130#ifdef DOXYGEN
     1131    template<class T>
     1132    void colLowerBound(T &t, Value value) { return 0;}
     1133#else
     1134    template<class T>
     1135    typename enable_if<typename T::value_type::LpSolverCol,void>::type
     1136    colLowerBound(T &t, Value value,dummy<0> = 0) {
     1137      for(typename T::iterator i=t.begin();i!=t.end();++i) {
     1138        colLowerBound(*i, value);
     1139      }
     1140    }
     1141    template<class T>
     1142    typename enable_if<typename T::value_type::second_type::LpSolverCol,
     1143                       void>::type
     1144    colLowerBound(T &t, Value value,dummy<1> = 1) {
     1145      for(typename T::iterator i=t.begin();i!=t.end();++i) {
     1146        colLowerBound(i->second, value);
     1147      }
     1148    }
     1149    template<class T>
     1150    typename enable_if<typename T::MapIt::Value::LpSolverCol,
     1151                       void>::type
     1152    colLowerBound(T &t, Value value,dummy<2> = 2) {
     1153      for(typename T::MapIt i(t); i!=INVALID; ++i){
     1154        colLowerBound(*i, value);
     1155      }
     1156    }
     1157#endif
     1158   
     1159    /// Set the upper bound of a column (i.e a variable)
     1160
     1161    /// The upper bound of a variable (column) has to be given by an
     1162    /// extended number of type Value, i.e. a finite number of type
     1163    /// Value or \ref INF.
     1164    void colUpperBound(Col c, Value value) {
     1165      _setColUpperBound(cols(id(c)),value);
     1166    };
     1167
     1168    /// Get the upper bound of a column (i.e a variable)
     1169
     1170    /// This function returns the upper bound for column (variable) \t c
     1171    /// (this might be \ref INF as well). 
     1172    ///\return The upper bound for column \t c
     1173    Value colUpperBound(Col c) const {
     1174      return _getColUpperBound(cols(id(c)));
     1175    }
     1176
     1177    ///\brief Set the upper bound of  several columns
     1178    ///(i.e a variables) at once
     1179    ///
     1180    ///This magic function takes a container as its argument
     1181    ///and applies the function on all of its elements.
     1182    /// The upper bound of a variable (column) has to be given by an
     1183    /// extended number of type Value, i.e. a finite number of type
     1184    /// Value or \ref INF.
     1185#ifdef DOXYGEN
     1186    template<class T>
     1187    void colUpperBound(T &t, Value value) { return 0;}
     1188#else
     1189    template<class T>
     1190    typename enable_if<typename T::value_type::LpSolverCol,void>::type
     1191    colUpperBound(T &t, Value value,dummy<0> = 0) {
     1192      for(typename T::iterator i=t.begin();i!=t.end();++i) {
     1193        colUpperBound(*i, value);
     1194      }
     1195    }
     1196    template<class T>
     1197    typename enable_if<typename T::value_type::second_type::LpSolverCol,
     1198                       void>::type
     1199    colUpperBound(T &t, Value value,dummy<1> = 1) {
     1200      for(typename T::iterator i=t.begin();i!=t.end();++i) {
     1201        colUpperBound(i->second, value);
     1202      }
     1203    }
     1204    template<class T>
     1205    typename enable_if<typename T::MapIt::Value::LpSolverCol,
     1206                       void>::type
     1207    colUpperBound(T &t, Value value,dummy<2> = 2) {
     1208      for(typename T::MapIt i(t); i!=INVALID; ++i){
     1209        colUpperBound(*i, value);
     1210      }
     1211    }
     1212#endif
     1213
     1214    /// Set the lower and the upper bounds of a column (i.e a variable)
     1215
     1216    /// The lower and the upper bounds of
     1217    /// a variable (column) have to be given by an
     1218    /// extended number of type Value, i.e. a finite number of type
     1219    /// Value, -\ref INF or \ref INF.
     1220    void colBounds(Col c, Value lower, Value upper) {
     1221      _setColLowerBound(cols(id(c)),lower);
     1222      _setColUpperBound(cols(id(c)),upper);
     1223    }
     1224   
     1225    ///\brief Set the lower and the upper bound of several columns
     1226    ///(i.e a variables) at once
     1227    ///
     1228    ///This magic function takes a container as its argument
     1229    ///and applies the function on all of its elements.
     1230    /// The lower and the upper bounds of
     1231    /// a variable (column) have to be given by an
     1232    /// extended number of type Value, i.e. a finite number of type
     1233    /// Value, -\ref INF or \ref INF.
     1234#ifdef DOXYGEN
     1235    template<class T>
     1236    void colBounds(T &t, Value lower, Value upper) { return 0;}
     1237#else
     1238    template<class T>
     1239    typename enable_if<typename T::value_type::LpSolverCol,void>::type
     1240    colBounds(T &t, Value lower, Value upper,dummy<0> = 0) {
     1241      for(typename T::iterator i=t.begin();i!=t.end();++i) {
     1242        colBounds(*i, lower, upper);
     1243      }
     1244    }
     1245    template<class T>
     1246    typename enable_if<typename T::value_type::second_type::LpSolverCol,
     1247                       void>::type
     1248    colBounds(T &t, Value lower, Value upper,dummy<1> = 1) {
     1249      for(typename T::iterator i=t.begin();i!=t.end();++i) {
     1250        colBounds(i->second, lower, upper);
     1251      }
     1252    }
     1253    template<class T>
     1254    typename enable_if<typename T::MapIt::Value::LpSolverCol,
     1255                       void>::type
     1256    colBounds(T &t, Value lower, Value upper,dummy<2> = 2) {
     1257      for(typename T::MapIt i(t); i!=INVALID; ++i){
     1258        colBounds(*i, lower, upper);
     1259      }
     1260    }
     1261#endif
     1262   
     1263    /// Set the lower bound of a row (i.e a constraint)
     1264
     1265    /// The lower bound of a constraint (row) has to be given by an
     1266    /// extended number of type Value, i.e. a finite number of type
     1267    /// Value or -\ref INF.
     1268    void rowLowerBound(Row r, Value value) {
     1269      _setRowLowerBound(rows(id(r)),value);
     1270    }
     1271
     1272    /// Get the lower bound of a row (i.e a constraint)
     1273
     1274    /// This function returns the lower bound for row (constraint) \t c
     1275    /// (this might be -\ref INF as well). 
     1276    ///\return The lower bound for row \t r
     1277    Value rowLowerBound(Row r) const {
     1278      return _getRowLowerBound(rows(id(r)));
     1279    }
     1280
     1281    /// Set the upper bound of a row (i.e a constraint)
     1282
     1283    /// The upper bound of a constraint (row) has to be given by an
     1284    /// extended number of type Value, i.e. a finite number of type
     1285    /// Value or -\ref INF.
     1286    void rowUpperBound(Row r, Value value) {
     1287      _setRowUpperBound(rows(id(r)),value);
     1288    }
     1289
     1290    /// Get the upper bound of a row (i.e a constraint)
     1291
     1292    /// This function returns the upper bound for row (constraint) \t c
     1293    /// (this might be -\ref INF as well). 
     1294    ///\return The upper bound for row \t r
     1295    Value rowUpperBound(Row r) const {
     1296      return _getRowUpperBound(rows(id(r)));
     1297    }
     1298
     1299    ///Set an element of the objective function
     1300    void objCoeff(Col c, Value v) {_setObjCoeff(cols(id(c)),v); };
     1301
     1302    ///Get an element of the objective function
     1303    Value objCoeff(Col c) const { return _getObjCoeff(cols(id(c))); };
     1304
     1305    ///Set the objective function
     1306
     1307    ///\param e is a linear expression of type \ref Expr.
     1308    void obj(const Expr& e) {
     1309      _setObjCoeffs(ExprIterator(e.comps.begin(), cols),
     1310                    ExprIterator(e.comps.end(), cols));
     1311      obj_const_comp = *e;
     1312    }
     1313
     1314    ///Get the objective function
     1315
     1316    ///\return the objective function as a linear expression of type \ref Expr.
     1317    Expr obj() const {
     1318      Expr e;
     1319      _getObjCoeffs(InsertIterator(e.comps, cols));
     1320      *e = obj_const_comp;
     1321      return e;
     1322    }
     1323   
     1324
     1325    ///Set the direction of optimization
     1326    void sense(Sense sense) { _setSense(sense); }
     1327
     1328    ///Query the direction of the optimization
     1329    Sense sense() const {return _getSense(); }
     1330
     1331    ///Set the sense to maximization
     1332    void max() { _setSense(MAX); }
     1333
     1334    ///Set the sense to maximization
     1335    void min() { _setSense(MIN); }
     1336
     1337    ///Clears the problem
     1338    void clear() { _clear(); }
     1339   
     1340    ///@}
     1341
     1342  };
     1343
     1344  ///\relates SolverBase::Expr
     1345  ///
     1346  inline SolverBase::Expr operator+(const SolverBase::Expr &a,
     1347                                    const SolverBase::Expr &b) {
     1348    SolverBase::Expr tmp(a);
     1349    tmp+=b;
     1350    return tmp;
     1351  }
     1352  ///\e
     1353 
     1354  ///\relates SolverBase::Expr
     1355  ///
     1356  inline SolverBase::Expr operator-(const SolverBase::Expr &a,
     1357                                    const SolverBase::Expr &b) {
     1358    SolverBase::Expr tmp(a);
     1359    tmp-=b;
     1360    return tmp;
     1361  }
     1362  ///\e
     1363 
     1364  ///\relates SolverBase::Expr
     1365  ///
     1366  inline SolverBase::Expr operator*(const SolverBase::Expr &a,
     1367                                    const SolverBase::Value &b) {
     1368    SolverBase::Expr tmp(a);
     1369    tmp*=b;
     1370    return tmp;
     1371  }
     1372 
     1373  ///\e
     1374 
     1375  ///\relates SolverBase::Expr
     1376  ///
     1377  inline SolverBase::Expr operator*(const SolverBase::Value &a,
     1378                                    const SolverBase::Expr &b) {
     1379    SolverBase::Expr tmp(b);
     1380    tmp*=a;
     1381    return tmp;
     1382  }
     1383  ///\e
     1384 
     1385  ///\relates SolverBase::Expr
     1386  ///
     1387  inline SolverBase::Expr operator/(const SolverBase::Expr &a,
     1388                                    const SolverBase::Value &b) {
     1389    SolverBase::Expr tmp(a);
     1390    tmp/=b;
     1391    return tmp;
     1392  }
     1393 
     1394  ///\e
     1395 
     1396  ///\relates SolverBase::Constr
     1397  ///
     1398  inline SolverBase::Constr operator<=(const SolverBase::Expr &e,
     1399                                       const SolverBase::Expr &f) {
     1400    return SolverBase::Constr(0, f - e, SolverBase::INF);
     1401  }
     1402
     1403  ///\e
     1404 
     1405  ///\relates SolverBase::Constr
     1406  ///
     1407  inline SolverBase::Constr operator<=(const SolverBase::Value &e,
     1408                                       const SolverBase::Expr &f) {
     1409    return SolverBase::Constr(e, f, SolverBase::NaN);
     1410  }
     1411
     1412  ///\e
     1413 
     1414  ///\relates SolverBase::Constr
     1415  ///
     1416  inline SolverBase::Constr operator<=(const SolverBase::Expr &e,
     1417                                       const SolverBase::Value &f) {
     1418    return SolverBase::Constr(- SolverBase::INF, e, f);
     1419  }
     1420
     1421  ///\e
     1422 
     1423  ///\relates SolverBase::Constr
     1424  ///
     1425  inline SolverBase::Constr operator>=(const SolverBase::Expr &e,
     1426                                       const SolverBase::Expr &f) {
     1427    return SolverBase::Constr(0, e - f, SolverBase::INF);
     1428  }
     1429
     1430
     1431  ///\e
     1432 
     1433  ///\relates SolverBase::Constr
     1434  ///
     1435  inline SolverBase::Constr operator>=(const SolverBase::Value &e,
     1436                                       const SolverBase::Expr &f) {
     1437    return SolverBase::Constr(SolverBase::NaN, f, e);
     1438  }
     1439
     1440
     1441  ///\e
     1442 
     1443  ///\relates SolverBase::Constr
     1444  ///
     1445  inline SolverBase::Constr operator>=(const SolverBase::Expr &e,
     1446                                       const SolverBase::Value &f) {
     1447    return SolverBase::Constr(f, e, SolverBase::INF);
     1448  }
     1449
     1450  ///\e
     1451
     1452  ///\relates SolverBase::Constr
     1453  ///
     1454  inline SolverBase::Constr operator==(const SolverBase::Expr &e,
     1455                                       const SolverBase::Value &f) {
     1456    return SolverBase::Constr(f, e, f);
     1457  }
     1458
     1459  ///\e
     1460 
     1461  ///\relates SolverBase::Constr
     1462  ///
     1463  inline SolverBase::Constr operator==(const SolverBase::Expr &e,
     1464                                       const SolverBase::Expr &f) {
     1465    return SolverBase::Constr(0, f - e, 0);
     1466  }
     1467
     1468  ///\e
     1469 
     1470  ///\relates SolverBase::Constr
     1471  ///
     1472  inline SolverBase::Constr operator<=(const SolverBase::Value &n,
     1473                                       const SolverBase::Constr &c) {
     1474    SolverBase::Constr tmp(c);
     1475    LEMON_ASSERT(std::isnan(tmp.lowerBound()), "Wrong LP constraint");
     1476    tmp.lowerBound()=n;
     1477    return tmp;
     1478  }
     1479  ///\e
     1480 
     1481  ///\relates SolverBase::Constr
     1482  ///
     1483  inline SolverBase::Constr operator<=(const SolverBase::Constr &c,
     1484                                       const SolverBase::Value &n)
     1485  {
     1486    SolverBase::Constr tmp(c);
     1487    LEMON_ASSERT(std::isnan(tmp.upperBound()), "Wrong LP constraint");
     1488    tmp.upperBound()=n;
     1489    return tmp;
     1490  }
     1491
     1492  ///\e
     1493 
     1494  ///\relates SolverBase::Constr
     1495  ///
     1496  inline SolverBase::Constr operator>=(const SolverBase::Value &n,
     1497                                       const SolverBase::Constr &c) {
     1498    SolverBase::Constr tmp(c);
     1499    LEMON_ASSERT(std::isnan(tmp.upperBound()), "Wrong LP constraint");
     1500    tmp.upperBound()=n;
     1501    return tmp;
     1502  }
     1503  ///\e
     1504 
     1505  ///\relates SolverBase::Constr
     1506  ///
     1507  inline SolverBase::Constr operator>=(const SolverBase::Constr &c,
     1508                                       const SolverBase::Value &n)
     1509  {
     1510    SolverBase::Constr tmp(c);
     1511    LEMON_ASSERT(std::isnan(tmp.lowerBound()), "Wrong LP constraint");
     1512    tmp.lowerBound()=n;
     1513    return tmp;
     1514  }
     1515
     1516  ///\e
     1517 
     1518  ///\relates SolverBase::DualExpr
     1519  ///
     1520  inline SolverBase::DualExpr operator+(const SolverBase::DualExpr &a,
     1521                                        const SolverBase::DualExpr &b) {
     1522    SolverBase::DualExpr tmp(a);
     1523    tmp+=b;
     1524    return tmp;
     1525  }
     1526  ///\e
     1527 
     1528  ///\relates SolverBase::DualExpr
     1529  ///
     1530  inline SolverBase::DualExpr operator-(const SolverBase::DualExpr &a,
     1531                                        const SolverBase::DualExpr &b) {
     1532    SolverBase::DualExpr tmp(a);
     1533    tmp-=b;
     1534    return tmp;
     1535  }
     1536  ///\e
     1537 
     1538  ///\relates SolverBase::DualExpr
     1539  ///
     1540  inline SolverBase::DualExpr operator*(const SolverBase::DualExpr &a,
     1541                                        const SolverBase::Value &b) {
     1542    SolverBase::DualExpr tmp(a);
     1543    tmp*=b;
     1544    return tmp;
     1545  }
     1546 
     1547  ///\e
     1548 
     1549  ///\relates SolverBase::DualExpr
     1550  ///
     1551  inline SolverBase::DualExpr operator*(const SolverBase::Value &a,
     1552                                        const SolverBase::DualExpr &b) {
     1553    SolverBase::DualExpr tmp(b);
     1554    tmp*=a;
     1555    return tmp;
     1556  }
     1557  ///\e
     1558 
     1559  ///\relates SolverBase::DualExpr
     1560  ///
     1561  inline SolverBase::DualExpr operator/(const SolverBase::DualExpr &a,
     1562                                        const SolverBase::Value &b) {
     1563    SolverBase::DualExpr tmp(a);
     1564    tmp/=b;
     1565    return tmp;
     1566  }
     1567
     1568  class LpSolverBase : virtual public SolverBase {
     1569  public:
     1570
     1571    ///\e
     1572    enum ProblemType {
     1573      ///Feasible solution hasn't been found (but may exist).
     1574      UNDEFINED = 0,
     1575      ///The problem has no feasible solution
     1576      INFEASIBLE = 1,
     1577      ///Feasible solution found
     1578      FEASIBLE = 2,
     1579      ///Optimal solution exists and found
     1580      OPTIMAL = 3,
     1581      ///The cost function is unbounded
     1582      UNBOUNDED = 4
     1583    };
     1584
     1585    enum VarStatus {
     1586      BASIC, FREE, LOWER, UPPER, FIXED
     1587    };
     1588
     1589  protected:
     1590
     1591    virtual SolveExitStatus _solve() = 0;
     1592
     1593    virtual Value _getPrimal(int i) const = 0;
     1594    virtual Value _getDual(int i) const = 0;
     1595
     1596    virtual Value _getPrimalRay(int i) const = 0;
     1597    virtual Value _getDualRay(int i) const = 0;
     1598
     1599    virtual Value _getPrimalValue() const = 0;
     1600
     1601    virtual VarStatus _getColStatus(int i) const = 0;
     1602    virtual VarStatus _getRowStatus(int i) const = 0;
     1603
     1604    virtual ProblemType _getPrimalType() const = 0;
     1605    virtual ProblemType _getDualType() const = 0;
     1606
     1607  public:
     1608
     1609    ///\name Solve the LP
     1610
     1611    ///@{
     1612
     1613    ///\e Solve the LP problem at hand
     1614    ///
     1615    ///\return The result of the optimization procedure. Possible
     1616    ///values and their meanings can be found in the documentation of
     1617    ///\ref SolveExitStatus.
     1618    ///
     1619    ///\todo Which method is used to solve the problem
     1620    SolveExitStatus solve() { return _solve(); }
     1621   
     1622    ///@}
     1623   
     1624    ///\name Obtain the solution
     1625
     1626    ///@{
     1627
     1628    /// The type of the primal problem
     1629    ProblemType primalType() const {
     1630      return _getPrimalType();
     1631    }
     1632
     1633    /// The type of the dual problem
     1634    ProblemType dualType() const {
     1635      return _getDualType();
     1636    }
     1637
     1638    ///\e
     1639    Value primal(Col c) const { return _getPrimal(cols(id(c))); }
     1640    ///\e
     1641    Value primal(const Expr& e) const {
     1642      double res = *e;
     1643      for (Expr::ConstMapIt c(e); c != INVALID; ++c) {
     1644        res += *c * primal(c);
     1645      }
     1646      return res;
     1647    }
     1648    ///\e
     1649    Value primalRay(Col c) const { return _getPrimalRay(cols(id(c))); }
     1650
     1651    ///\e
     1652    Value dual(Row r) const { return _getDual(rows(id(r))); }
     1653    ///\e
     1654    Value dual(const DualExpr& e) const {
     1655      double res = 0.0;
     1656      for (DualExpr::ConstMapIt r(e); r != INVALID; ++r) {
     1657        res += *r * dual(r);
     1658      }
     1659      return res;
     1660    }
     1661    ///\e
     1662    Value dualRay(Row r) const { return _getDualRay(rows(id(r))); }
     1663
     1664    ///\e
     1665    VarStatus colStatus(Col c) const { return _getColStatus(cols(id(c))); }
     1666
     1667    ///\e
     1668    VarStatus rowStatus(Col c) const { return _getRowStatus(cols(id(c))); }
     1669
     1670    ///\e
     1671
     1672    ///\return
     1673    ///- \ref INF or -\ref INF means either infeasibility or unboundedness
     1674    /// of the primal problem, depending on whether we minimize or maximize.
     1675    ///- \ref NaN if no primal solution is found.
     1676    ///- The (finite) objective value if an optimal solution is found.
     1677    Value primal() const { return _getPrimalValue()+obj_const_comp;}
     1678    ///@}
     1679   
     1680    LpSolverBase* newSolver() {return _newSolver();}
     1681    LpSolverBase* cloneSolver() {return _cloneSolver();}
     1682
     1683  protected:
     1684
     1685    virtual LpSolverBase* _newSolver() const = 0;
     1686    virtual LpSolverBase* _cloneSolver() const = 0;
     1687  }; 
     1688
     1689
     1690  /// \ingroup lp_group
     1691  ///
     1692  /// \brief Common base class for MIP solvers
     1693  /// \todo Much more docs
     1694  class MipSolverBase : virtual public SolverBase {
     1695  public:
     1696
     1697    ///\e
     1698    enum ProblemType {
     1699      ///Feasible solution hasn't been found (but may exist).
     1700      UNDEFINED = 0,
     1701      ///The problem has no feasible solution
     1702      INFEASIBLE = 1,
     1703      ///Feasible solution found
     1704      FEASIBLE = 2,
     1705      ///Optimal solution exists and found
     1706      OPTIMAL = 3,
     1707      ///The cost function is unbounded
     1708      ///
     1709      ///The Mip or at least the relaxed problem is unbounded
     1710      UNBOUNDED = 4
     1711    };
     1712
     1713    ///\name Solve the MIP
     1714
     1715    ///@{
     1716
     1717    ///\e Solve the MIP problem at hand
     1718    ///
     1719    ///\return The result of the optimization procedure. Possible
     1720    ///values and their meanings can be found in the documentation of
     1721    ///\ref SolveExitStatus.
     1722    ///
     1723    ///\todo Which method is used to solve the problem
     1724    SolveExitStatus solve() { return _solve(); }
     1725   
     1726    ///@}
     1727
     1728    ///\name Setting column type
     1729    ///@{
     1730
     1731    ///Possible variable (column) types (e.g. real, integer, binary etc.)
     1732    enum ColTypes {
     1733      ///Continuous variable (default)
     1734      REAL = 0,
     1735      ///Integer variable
     1736      INTEGER = 1
     1737    };
     1738
     1739    ///Sets the type of the given column to the given type
     1740    ///
     1741    ///Sets the type of the given column to the given type.
     1742    void colType(Col c, ColTypes col_type) {
     1743      _setColType(cols(id(c)),col_type);
     1744    }
     1745
     1746    ///Gives back the type of the column.
     1747    ///
     1748    ///Gives back the type of the column.
     1749    ColTypes colType(Col c) const {
     1750      return _getColType(cols(id(c)));
     1751    }
     1752
     1753    ///Sets the type of the given Col to integer.
     1754    ///
     1755    ///Sets the type of the given Col to integer.
     1756    void integer(Col c) {
     1757      colType(c, INTEGER);
     1758    }
     1759
     1760    ///Gives back whether the type of the column is integer or not.
     1761    ///
     1762    ///Gives back the type of the column.
     1763    ///\return true if the column has integer type and false if not.
     1764    bool integer(Col c) const {
     1765      return (colType(c) == INTEGER);
     1766    }
     1767
     1768    ///Sets the type of the given Col to continuous.
     1769    ///
     1770    ///Sets the type of the given Col to continuous.
     1771    void real(Col c) {
     1772      colType(c, REAL);
     1773    }
     1774
     1775    ///Gives back whether the type of the column is continuous or not.
     1776    ///
     1777    ///Gives back the type of the column.
     1778    ///\return true if the column has continuous type and false if not.
     1779    bool real(Col c) const {
     1780      return (colType(c) == REAL);
     1781    }
     1782    ///@}
     1783
     1784    ///\name Obtain the solution
     1785
     1786    ///@{
     1787
     1788    /// The type of the MIP problem
     1789    ProblemType type() const {
     1790      return _getType();
     1791    }
     1792
     1793    ///\e
     1794    Value sol(Col c) const { return _getSol(cols(id(c))); }
     1795    ///\e
     1796    Value sol(const Expr& e) const {
     1797      double res = *e;
     1798      for (Expr::ConstMapIt c(e); c != INVALID; ++c) {
     1799        res += *c * sol(c);
     1800      }
     1801      return res;
     1802    }
     1803    ///\return
     1804    ///- \ref INF or -\ref INF means either infeasibility or unboundedness
     1805    /// of the problem, depending on whether we minimize or maximize.
     1806    ///- \ref NaN if no primal solution is found.
     1807    ///- The (finite) objective value if an optimal solution is found.
     1808    Value solValue() const { return _getSolValue()+obj_const_comp;}
     1809    ///@}
     1810
     1811  protected:
     1812
     1813    virtual SolveExitStatus _solve() = 0;
     1814    virtual ColTypes _getColType(int col) const = 0;
     1815    virtual void _setColType(int col, ColTypes col_type) = 0;
     1816    virtual ProblemType _getType() const = 0;
     1817    virtual Value _getSol(int i) const = 0;
     1818    virtual Value _getSolValue() const = 0;
     1819
     1820  public:
     1821
     1822    MipSolverBase* newSolver() {return _newSolver();}
     1823    MipSolverBase* cloneSolver() {return _cloneSolver();}
     1824
     1825  protected:
     1826
     1827    virtual MipSolverBase* _newSolver() const = 0;
     1828    virtual MipSolverBase* _cloneSolver() const = 0;
     1829  };
     1830 
     1831 
     1832
     1833} //namespace lemon
     1834
     1835#endif //LEMON_LP_BASE_H
  • new file lemon/lp_clp.cc

    diff -r cbe3ec2d59d2 -r aa6dff56b2e1 lemon/lp_clp.cc
    - +  
     1/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library.
     4 *
     5 * Copyright (C) 2003-2008
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#include <lemon/lp_clp.h>
     20#include <coin/ClpSimplex.hpp>
     21
     22namespace lemon {
     23
     24  LpClp::LpClp() {
     25    _prob = new ClpSimplex();
     26    _init_temporals();
     27    messageLevel(MESSAGE_NO_OUTPUT);
     28  }
     29
     30  LpClp::LpClp(const LpClp& other) {
     31    _prob = new ClpSimplex(*other._prob);
     32    rows = other.rows;
     33    cols = other.cols;
     34    _init_temporals();
     35    messageLevel(MESSAGE_NO_OUTPUT);
     36  }
     37
     38  LpClp::~LpClp() {
     39    delete _prob;
     40    _clear_temporals();
     41  }
     42
     43  void LpClp::_init_temporals() {
     44    _primal_ray = 0;
     45    _dual_ray = 0;
     46  }
     47
     48  void LpClp::_clear_temporals() {
     49    if (_primal_ray) {
     50      delete[] _primal_ray;
     51      _primal_ray = 0;
     52    }
     53    if (_dual_ray) {
     54      delete[] _dual_ray;
     55      _dual_ray = 0;
     56    }
     57  }
     58
     59  LpClp* LpClp::_newSolver() const {
     60    LpClp* newlp = new LpClp;
     61    return newlp;
     62  }
     63
     64  LpClp* LpClp::_cloneSolver() const {
     65    LpClp* copylp = new LpClp(*this);
     66    return copylp;
     67  }
     68
     69  int LpClp::_addCol() {
     70    _prob->addColumn(0, 0, 0, -COIN_DBL_MAX, COIN_DBL_MAX, 0.0);
     71    return _prob->numberColumns() - 1;
     72  }
     73
     74  int LpClp::_addRow() {
     75    _prob->addRow(0, 0, 0, -COIN_DBL_MAX, COIN_DBL_MAX);
     76    return _prob->numberRows() - 1;
     77  }
     78
     79
     80  void LpClp::_eraseCol(int c) {
     81    _col_names_ref.erase(_prob->getColumnName(c));
     82    _prob->deleteColumns(1, &c);
     83  }
     84
     85  void LpClp::_eraseRow(int r) {
     86    _row_names_ref.erase(_prob->getRowName(r));
     87    _prob->deleteRows(1, &r);
     88  }
     89
     90  void LpClp::_eraseColId(int i) {
     91    cols.eraseIndex(i);
     92    cols.shiftIndexes(i);
     93  }
     94
     95  void LpClp::_eraseRowId(int i) {
     96    rows.eraseIndex(i);
     97    rows.shiftIndexes(i);
     98  }
     99
     100  void LpClp::_getColName(int c, std::string& name) const {
     101    name = _prob->getColumnName(c);
     102  }
     103
     104  void LpClp::_setColName(int c, const std::string& name) {
     105    _prob->setColumnName(c, const_cast<std::string&>(name));
     106    _col_names_ref[name] = c;
     107  }
     108
     109  int LpClp::_colByName(const std::string& name) const {
     110    std::map<std::string, int>::const_iterator it = _col_names_ref.find(name);
     111    return it != _col_names_ref.end() ? it->second : -1;
     112  }
     113
     114  void LpClp::_getRowName(int r, std::string& name) const {
     115    name = _prob->getRowName(r);
     116  }
     117
     118  void LpClp::_setRowName(int r, const std::string& name) {
     119    _prob->setRowName(r, const_cast<std::string&>(name));
     120    _row_names_ref[name] = r;
     121  }
     122
     123  int LpClp::_rowByName(const std::string& name) const {
     124    std::map<std::string, int>::const_iterator it = _row_names_ref.find(name);
     125    return it != _row_names_ref.end() ? it->second : -1;
     126  }
     127
     128
     129  void LpClp::_setRowCoeffs(int ix, ExprIterator b, ExprIterator e) {
     130    std::map<int, Value> coeffs;
     131
     132    int n = _prob->clpMatrix()->getNumCols();
     133
     134    const int* indices = _prob->clpMatrix()->getIndices();
     135    const double* elements = _prob->clpMatrix()->getElements();
     136
     137    for (int i = 0; i < n; ++i) {
     138      CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[i];
     139      CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[i];
     140
     141      const int* it = std::lower_bound(indices + begin, indices + end, ix);
     142      if (it != indices + end && *it == ix && elements[it - indices] != 0.0) {
     143        coeffs[i] = 0.0;
     144      }
     145    }
     146   
     147    for (ExprIterator it = b; it != e; ++it) {
     148      coeffs[it->first] = it->second;
     149    }
     150   
     151    for (std::map<int, Value>::iterator it = coeffs.begin();
     152         it != coeffs.end(); ++it) {
     153      _prob->modifyCoefficient(ix, it->first, it->second);
     154    }   
     155  }
     156
     157  void LpClp::_getRowCoeffs(int ix, InsertIterator b) const {
     158    int n = _prob->clpMatrix()->getNumCols();
     159
     160    const int* indices = _prob->clpMatrix()->getIndices();
     161    const double* elements = _prob->clpMatrix()->getElements();
     162
     163    for (int i = 0; i < n; ++i) {
     164      CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[i];
     165      CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[i];
     166
     167      const int* it = std::lower_bound(indices + begin, indices + end, ix);
     168      if (it != indices + end && *it == ix) {
     169        *b = std::make_pair(i, elements[it - indices]);
     170      }
     171    }
     172  }
     173
     174  void LpClp::_setColCoeffs(int ix, ExprIterator b, ExprIterator e) {
     175    std::map<int, Value> coeffs;
     176
     177    CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[ix];
     178    CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[ix];
     179   
     180    const int* indices = _prob->clpMatrix()->getIndices();
     181    const double* elements = _prob->clpMatrix()->getElements();
     182
     183    for (CoinBigIndex i = begin; i != end; ++i) {
     184      if (elements[i] != 0.0) {
     185        coeffs[indices[i]] = 0.0;
     186      }
     187    }
     188    for (ExprIterator it = b; it != e; ++it) {     
     189      coeffs[it->first] = it->second;
     190    }
     191    for (std::map<int, Value>::iterator it = coeffs.begin();
     192         it != coeffs.end(); ++it) {
     193      _prob->modifyCoefficient(it->first, ix, it->second);
     194    }
     195  }
     196
     197  void LpClp::_getColCoeffs(int ix, InsertIterator b) const {
     198    CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[ix];
     199    CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[ix];
     200   
     201    const int* indices = _prob->clpMatrix()->getIndices();
     202    const double* elements = _prob->clpMatrix()->getElements();
     203
     204    for (CoinBigIndex i = begin; i != end; ++i) {
     205      *b = std::make_pair(indices[i], elements[i]);
     206      ++b;
     207    }
     208  }
     209
     210  void LpClp::_setCoeff(int ix, int jx, Value value) {
     211    _prob->modifyCoefficient(ix, jx, value);
     212  }
     213
     214  LpClp::Value LpClp::_getCoeff(int ix, int jx) const {
     215    CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[ix];
     216    CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[ix];
     217   
     218    const int* indices = _prob->clpMatrix()->getIndices();
     219    const double* elements = _prob->clpMatrix()->getElements();
     220
     221    const int* it = std::lower_bound(indices + begin, indices + end, jx);
     222    if (it != indices + end && *it == jx) {
     223      return elements[it - indices];
     224    } else {
     225      return 0.0;
     226    }
     227  }
     228
     229  void LpClp::_setColLowerBound(int i, Value lo) {
     230    _prob->setColumnLower(i, lo == - INF ? - COIN_DBL_MAX : lo);
     231  }
     232
     233  LpClp::Value LpClp::_getColLowerBound(int i) const {
     234    double val = _prob->getColLower()[i];
     235    return val == - COIN_DBL_MAX ? - INF : val;
     236  }
     237
     238  void LpClp::_setColUpperBound(int i, Value up) {
     239    _prob->setColumnUpper(i, up == INF ? COIN_DBL_MAX : up);
     240  }
     241
     242  LpClp::Value LpClp::_getColUpperBound(int i) const {
     243    double val = _prob->getColUpper()[i];
     244    return val == COIN_DBL_MAX ? INF : val;
     245  }
     246
     247  void LpClp::_setRowLowerBound(int i, Value lo) {
     248    _prob->setRowLower(i, lo == - INF ? - COIN_DBL_MAX : lo);
     249  }
     250
     251  LpClp::Value LpClp::_getRowLowerBound(int i) const {
     252    double val = _prob->getRowLower()[i];
     253    return val == - COIN_DBL_MAX ? - INF : val;
     254  }
     255
     256  void LpClp::_setRowUpperBound(int i, Value up) {
     257    _prob->setRowUpper(i, up == INF ? COIN_DBL_MAX : up);
     258  }
     259
     260  LpClp::Value LpClp::_getRowUpperBound(int i) const {
     261    double val = _prob->getRowUpper()[i];
     262    return val == COIN_DBL_MAX ? INF : val;
     263  }
     264
     265  void LpClp::_setObjCoeffs(ExprIterator b, ExprIterator e) {
     266    int num = _prob->clpMatrix()->getNumCols();
     267    for (int i = 0; i < num; ++i) {
     268      _prob->setObjectiveCoefficient(i, 0.0);
     269    }
     270    for (ExprIterator it = b; it != e; ++it) {
     271      _prob->setObjectiveCoefficient(it->first, it->second);
     272    }
     273  }
     274
     275  void LpClp::_getObjCoeffs(InsertIterator b) const {
     276    int num = _prob->clpMatrix()->getNumCols();
     277    for (int i = 0; i < num; ++i) {
     278      Value coef = _prob->getObjCoefficients()[i];
     279      if (coef != 0.0) {
     280        *b = std::make_pair(i, coef);
     281        ++b;
     282      }
     283    }
     284  }
     285
     286  void LpClp::_setObjCoeff(int i, Value obj_coef) {
     287    _prob->setObjectiveCoefficient(i, obj_coef);
     288  }
     289
     290  LpClp::Value LpClp::_getObjCoeff(int i) const {
     291    return _prob->getObjCoefficients()[i];
     292  }
     293
     294  LpClp::SolveExitStatus LpClp::_solve() {
     295    return _prob->primal() >= 0 ? SOLVED : UNSOLVED;
     296  }
     297 
     298  LpClp::SolveExitStatus LpClp::solvePrimal() {
     299    return _prob->primal() >= 0 ? SOLVED : UNSOLVED;
     300  }
     301
     302  LpClp::SolveExitStatus LpClp::solveDual() {
     303    return _prob->dual() >= 0 ? SOLVED : UNSOLVED;
     304  }
     305
     306  LpClp::SolveExitStatus LpClp::solveBarrier() {
     307    return _prob->barrier() >= 0 ? SOLVED : UNSOLVED;
     308  }
     309
     310  LpClp::Value LpClp::_getPrimal(int i) const {
     311    return _prob->primalColumnSolution()[i];
     312  }
     313  LpClp::Value LpClp::_getPrimalValue() const {
     314    return _prob->rawObjectiveValue();
     315  }
     316
     317  LpClp::Value LpClp::_getDual(int i) const {
     318    return _prob->dualRowSolution()[i];
     319  }
     320
     321  LpClp::Value LpClp::_getPrimalRay(int i) const {
     322    if (!_primal_ray) {
     323      _primal_ray = _prob->unboundedRay();
     324      LEMON_ASSERT(_primal_ray != 0, "Primal ray is not provided");
     325    }
     326    return _primal_ray[i];
     327  }
     328
     329  LpClp::Value LpClp::_getDualRay(int i) const {
     330    if (!_dual_ray) {
     331      _dual_ray = _prob->infeasibilityRay();
     332      LEMON_ASSERT(_dual_ray != 0, "Dual ray is not provided");
     333    }
     334    return _dual_ray[i];
     335  }
     336
     337  LpClp::VarStatus LpClp::_getColStatus(int i) const {
     338    switch (_prob->getColumnStatus(i)) {
     339    case ClpSimplex::basic:
     340      return BASIC;
     341    case ClpSimplex::isFree:
     342      return FREE;
     343    case ClpSimplex::atUpperBound:
     344      return UPPER;
     345    case ClpSimplex::atLowerBound:
     346      return LOWER;
     347    case ClpSimplex::isFixed:
     348      return FIXED;
     349    case ClpSimplex::superBasic:
     350      return FREE;
     351    default:
     352      LEMON_ASSERT(false, "Wrong column status");
     353      return VarStatus();
     354    }
     355  }
     356
     357  LpClp::VarStatus LpClp::_getRowStatus(int i) const {
     358    switch (_prob->getColumnStatus(i)) {
     359    case ClpSimplex::basic:
     360      return BASIC;
     361    case ClpSimplex::isFree:
     362      return FREE;
     363    case ClpSimplex::atUpperBound:
     364      return UPPER;
     365    case ClpSimplex::atLowerBound:
     366      return LOWER;     
     367    case ClpSimplex::isFixed:
     368      return FIXED;
     369    case ClpSimplex::superBasic:
     370      return FREE;
     371    default:
     372      LEMON_ASSERT(false, "Wrong row status");
     373      return VarStatus();
     374    }
     375  }
     376
     377
     378  LpClp::ProblemType LpClp::_getPrimalType() const {
     379    if (_prob->isProvenOptimal()) {
     380      return OPTIMAL;
     381    } else if (_prob->isProvenPrimalInfeasible()) {
     382      return INFEASIBLE;
     383    } else if (_prob->isProvenDualInfeasible()) {
     384      return UNBOUNDED;
     385    } else {
     386      return UNDEFINED;
     387    }
     388  }
     389
     390  LpClp::ProblemType LpClp::_getDualType() const {
     391    if (_prob->isProvenOptimal()) {
     392      return OPTIMAL;
     393    } else if (_prob->isProvenDualInfeasible()) {
     394      return INFEASIBLE;     
     395    } else if (_prob->isProvenPrimalInfeasible()) {
     396      return INFEASIBLE;
     397    } else {
     398      return UNDEFINED;
     399    }
     400  }
     401
     402  void LpClp::_setSense(SolverBase::Sense sense) {
     403    switch (sense) {
     404    case MIN:
     405      _prob->setOptimizationDirection(1);
     406      break;
     407    case MAX:
     408      _prob->setOptimizationDirection(-1);
     409      break;
     410    }
     411  }
     412
     413  SolverBase::Sense LpClp::_getSense() const {
     414    double dir = _prob->optimizationDirection();
     415    if (dir > 0.0) {
     416      return MIN;
     417    } else {
     418      return MAX;
     419    }
     420  }
     421
     422  void LpClp::_clear() {
     423    delete _prob;
     424    _prob = new ClpSimplex();
     425    rows.clear();
     426    cols.clear();
     427    _col_names_ref.clear();
     428    _clear_temporals();
     429  }
     430
     431  void LpClp::messageLevel(MessageLevel m) {
     432    _prob->setLogLevel(static_cast<int>(m));
     433  }
     434
     435} //END OF NAMESPACE LEMON
  • new file lemon/lp_clp.h

    diff -r cbe3ec2d59d2 -r aa6dff56b2e1 lemon/lp_clp.h
    - +  
     1/* -*- C++ -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library
     4 *
     5 * Copyright (C) 2003-2008
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#ifndef LEMON_LP_CLP_H
     20#define LEMON_LP_CLP_H
     21
     22///\file
     23///\brief Header of the LEMON-CLP lp solver interface.
     24
     25#include <vector>
     26#include <string>
     27
     28#include <lemon/lp_base.h>
     29
     30class ClpSimplex;
     31
     32namespace lemon {
     33
     34  /// \ingroup lp_group
     35  ///
     36  /// \brief Interface for the CLP solver
     37  ///
     38  /// This class implements an interface for the Clp LP solver.  The
     39  /// Clp library is an object oriented lp solver library developed at
     40  /// the IBM. The CLP is part of the COIN-OR package and it can be
     41  /// used with Common Public License.
     42  class LpClp : public LpSolverBase {
     43  protected:
     44
     45    ClpSimplex* _prob;
     46
     47    std::map<std::string, int> _col_names_ref;
     48    std::map<std::string, int> _row_names_ref;
     49
     50  public:
     51
     52    typedef LpSolverBase Parent;
     53   
     54    /// \e
     55    LpClp();
     56    /// \e
     57    LpClp(const LpClp&);
     58    /// \e
     59    ~LpClp();
     60
     61  protected:
     62
     63    mutable double* _primal_ray;
     64    mutable double* _dual_ray;
     65
     66    void _init_temporals();
     67    void _clear_temporals();
     68
     69  protected:
     70
     71    virtual LpClp* _newSolver() const;
     72    virtual LpClp* _cloneSolver() const ;
     73
     74    virtual int _addCol();
     75    virtual int _addRow();
     76
     77    virtual void _eraseCol(int i);
     78    virtual void _eraseRow(int i);
     79
     80    virtual void _eraseColId(int i);
     81    virtual void _eraseRowId(int i);
     82
     83    virtual void _getColName(int col, std::string& name) const;
     84    virtual void _setColName(int col, const std::string& name);
     85    virtual int _colByName(const std::string& name) const;
     86
     87    virtual void _getRowName(int row, std::string& name) const;
     88    virtual void _setRowName(int row, const std::string& name);
     89    virtual int _rowByName(const std::string& name) const;
     90
     91    virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e);
     92    virtual void _getRowCoeffs(int i, InsertIterator b) const;
     93
     94    virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e);
     95    virtual void _getColCoeffs(int i, InsertIterator b) const;
     96
     97    virtual void _setCoeff(int row, int col, Value value);
     98    virtual Value _getCoeff(int row, int col) const;
     99
     100    virtual void _setColLowerBound(int i, Value value);
     101    virtual Value _getColLowerBound(int i) const;
     102    virtual void _setColUpperBound(int i, Value value);
     103    virtual Value _getColUpperBound(int i) const;
     104
     105    virtual void _setRowLowerBound(int i, Value value);
     106    virtual Value _getRowLowerBound(int i) const;
     107    virtual void _setRowUpperBound(int i, Value value);
     108    virtual Value _getRowUpperBound(int i) const;
     109
     110    virtual void _setObjCoeffs(ExprIterator, ExprIterator);
     111    virtual void _getObjCoeffs(InsertIterator) const;
     112
     113    virtual void _setObjCoeff(int i, Value obj_coef);
     114    virtual Value _getObjCoeff(int i) const;
     115   
     116    virtual void _setSense(Sense sense);
     117    virtual Sense _getSense() const;
     118
     119    virtual SolveExitStatus _solve();
     120
     121    virtual Value _getPrimal(int i) const;
     122    virtual Value _getDual(int i) const;
     123
     124    virtual Value _getPrimalValue() const;
     125
     126    virtual Value _getPrimalRay(int i) const;
     127    virtual Value _getDualRay(int i) const;
     128
     129    virtual VarStatus _getColStatus(int i) const;
     130    virtual VarStatus _getRowStatus(int i) const;
     131   
     132    virtual ProblemType _getPrimalType() const;
     133    virtual ProblemType _getDualType() const;
     134
     135    virtual void _clear();
     136
     137  public:
     138
     139    ///Solves LP with primal simplex method.
     140    SolveExitStatus solvePrimal();
     141
     142    ///Solves LP with dual simplex method.
     143    SolveExitStatus solveDual();
     144
     145    ///Solves LP with barrier method.
     146    SolveExitStatus solveBarrier();
     147   
     148    ///Returns the constraint identifier understood by CLP.
     149    int clpRow(Row r) const { return rows(id(r)); }
     150
     151    ///Returns the variable identifier understood by CLP.
     152    int clpCol(Col c) const { return cols(id(c)); }
     153
     154    ///Enum for \c messageLevel() parameter
     155    enum MessageLevel {
     156      /// no output (default value)
     157      MESSAGE_NO_OUTPUT = 0,
     158      /// print final solution
     159      MESSAGE_FINAL_SOLUTION = 1,
     160      /// print factorization
     161      MESSAGE_FACTORIZATION = 2,
     162      /// normal output
     163      MESSAGE_NORMAL_OUTPUT = 3,
     164      /// verbose output
     165      MESSAGE_VERBOSE_OUTPUT = 4
     166    };
     167    ///Set the verbosity of the messages
     168
     169    ///Set the verbosity of the messages
     170    ///
     171    ///\param m is the level of the messages output by the solver routines.
     172    void messageLevel(MessageLevel m);
     173
     174  };
     175
     176} //END OF NAMESPACE LEMON
     177
     178#endif //LEMON_LP_CLP_H
     179
  • new file lemon/lp_cplex.cc

    diff -r cbe3ec2d59d2 -r aa6dff56b2e1 lemon/lp_cplex.cc
    - +  
     1/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library.
     4 *
     5 * Copyright (C) 2003-2008
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#include <iostream>
     20#include <vector>
     21#include <cstring>
     22
     23#include <lemon/lp_cplex.h>
     24
     25extern "C" {
     26#include <ilcplex/cplex.h>
     27}
     28
     29
     30///\file
     31///\brief Implementation of the LEMON-CPLEX lp solver interface.
     32namespace lemon {
     33
     34  CplexEnv::LicenseError::LicenseError(int status) {
     35    if (!CPXgeterrorstring(0, status, _message)) {
     36      std::strcpy(_message, "Cplex unknown error");
     37    }
     38  }
     39
     40  CplexEnv::CplexEnv() {
     41    int status;
     42    _cnt = new int;
     43    _env = CPXopenCPLEX(&status);
     44    if (_env == 0) {
     45      delete _cnt;
     46      _cnt = 0;
     47      throw LicenseError(status);
     48    }
     49  }
     50
     51  CplexEnv::CplexEnv(const CplexEnv& other) {
     52    _env = other._env;
     53    _cnt = other._cnt;
     54    ++(*_cnt);
     55  }
     56
     57  CplexEnv& CplexEnv::operator=(const CplexEnv& other) {
     58    _env = other._env;
     59    _cnt = other._cnt;
     60    ++(*_cnt);
     61    return *this;
     62  }
     63
     64  CplexEnv::~CplexEnv() {
     65    --(*_cnt);
     66    if (*_cnt == 0) {
     67      delete _cnt;
     68      CPXcloseCPLEX(&_env);
     69    }
     70  }
     71
     72  CplexSolverBase::CplexSolverBase() : SolverBase() {
     73    int status;
     74    _prob = CPXcreateprob(cplexEnv(), &status, "Cplex problem");
     75  }
     76
     77  CplexSolverBase::CplexSolverBase(const CplexEnv& env)
     78    : SolverBase(), _env(env) {
     79    int status;
     80    _prob = CPXcreateprob(cplexEnv(), &status, "Cplex problem");
     81  }
     82
     83  CplexSolverBase::CplexSolverBase(const CplexSolverBase& cplex)
     84    : SolverBase() {
     85    int status;
     86    _prob = CPXcloneprob(cplexEnv(), cplex._prob, &status);
     87    rows = cplex.rows;
     88    cols = cplex.cols;
     89  }
     90
     91  CplexSolverBase::~CplexSolverBase() {
     92    CPXfreeprob(cplexEnv(),&_prob);
     93  }
     94
     95  int CplexSolverBase::_addCol() {
     96    int i = CPXgetnumcols(cplexEnv(), _prob);
     97    double lb = -INF, ub = INF;
     98    CPXnewcols(cplexEnv(), _prob, 1, 0, &lb, &ub, 0, 0);
     99    return i;
     100  }
     101
     102
     103  int CplexSolverBase::_addRow() {
     104    int i = CPXgetnumrows(cplexEnv(), _prob);
     105    const double ub = INF;
     106    const char s = 'L';
     107    CPXnewrows(cplexEnv(), _prob, 1, &ub, &s, 0, 0);
     108    return i;
     109  }
     110
     111
     112  void CplexSolverBase::_eraseCol(int i) {
     113    CPXdelcols(cplexEnv(), _prob, i, i);
     114  }
     115
     116  void CplexSolverBase::_eraseRow(int i) {
     117    CPXdelrows(cplexEnv(), _prob, i, i);
     118  }
     119
     120  void CplexSolverBase::_eraseColId(int i) {
     121    cols.eraseIndex(i);
     122    cols.shiftIndexes(i);
     123  }
     124  void CplexSolverBase::_eraseRowId(int i) {
     125    rows.eraseIndex(i);
     126    rows.shiftIndexes(i);
     127  }
     128
     129  void CplexSolverBase::_getColName(int col, std::string &name) const {
     130    int size;
     131    CPXgetcolname(cplexEnv(), _prob, 0, 0, 0, &size, col, col);
     132    if (size == 0) {
     133      name.clear();
     134      return;
     135    }
     136
     137    size *= -1;
     138    std::vector<char> buf(size);
     139    char *cname;
     140    int tmp;
     141    CPXgetcolname(cplexEnv(), _prob, &cname, &buf.front(), size,
     142                  &tmp, col, col);
     143    name = cname;
     144  }
     145
     146  void CplexSolverBase::_setColName(int col, const std::string &name) {
     147    char *cname;
     148    cname = const_cast<char*>(name.c_str());
     149    CPXchgcolname(cplexEnv(), _prob, 1, &col, &cname);
     150  }
     151
     152  int CplexSolverBase::_colByName(const std::string& name) const {
     153    int index;
     154    if (CPXgetcolindex(cplexEnv(), _prob,
     155                       const_cast<char*>(name.c_str()), &index) == 0) {
     156      return index;
     157    }
     158    return -1;
     159  }
     160
     161  void CplexSolverBase::_getRowName(int row, std::string &name) const {
     162    int size;
     163    CPXgetrowname(cplexEnv(), _prob, 0, 0, 0, &size, row, row);
     164    if (size == 0) {
     165      name.clear();
     166      return;
     167    }
     168
     169    size *= -1;
     170    std::vector<char> buf(size);
     171    char *cname;
     172    int tmp;
     173    CPXgetrowname(cplexEnv(), _prob, &cname, &buf.front(), size,
     174                  &tmp, row, row);
     175    name = cname;
     176  }
     177
     178  void CplexSolverBase::_setRowName(int row, const std::string &name) {
     179    char *cname;
     180    cname = const_cast<char*>(name.c_str());
     181    CPXchgrowname(cplexEnv(), _prob, 1, &row, &cname);
     182  }
     183
     184  int CplexSolverBase::_rowByName(const std::string& name) const {
     185    int index;
     186    if (CPXgetrowindex(cplexEnv(), _prob,
     187                       const_cast<char*>(name.c_str()), &index) == 0) {
     188      return index;
     189    }
     190    return -1;
     191  }
     192
     193  void CplexSolverBase::_setRowCoeffs(int i, ExprIterator b,
     194                                      ExprIterator e)
     195  {
     196    std::vector<int> indices;
     197    std::vector<int> rowlist;
     198    std::vector<Value> values;
     199
     200    for(ExprIterator it=b; it!=e; ++it) {
     201      indices.push_back(it->first);
     202      values.push_back(it->second);
     203      rowlist.push_back(i);
     204    }
     205
     206    CPXchgcoeflist(cplexEnv(), _prob, values.size(),
     207                   &rowlist.front(), &indices.front(), &values.front());
     208  }
     209
     210  void CplexSolverBase::_getRowCoeffs(int i, InsertIterator b) const {
     211    int tmp1, tmp2, tmp3, length;
     212    CPXgetrows(cplexEnv(), _prob, &tmp1, &tmp2, 0, 0, 0, &length, i, i);
     213
     214    length = -length;
     215    std::vector<int> indices(length);
     216    std::vector<double> values(length);
     217
     218    CPXgetrows(cplexEnv(), _prob, &tmp1, &tmp2,
     219               &indices.front(), &values.front(),
     220               length, &tmp3, i, i);
     221
     222    for (int i = 0; i < length; ++i) {
     223      *b = std::make_pair(indices[i], values[i]);
     224      ++b;
     225    }
     226  }
     227
     228  void CplexSolverBase::_setColCoeffs(int i, ExprIterator b, ExprIterator e) {
     229    std::vector<int> indices;
     230    std::vector<int> collist;
     231    std::vector<Value> values;
     232
     233    for(ExprIterator it=b; it!=e; ++it) {
     234      indices.push_back(it->first);
     235      values.push_back(it->second);
     236      collist.push_back(i);
     237    }
     238
     239    CPXchgcoeflist(cplexEnv(), _prob, values.size(),
     240                   &indices.front(), &collist.front(), &values.front());
     241  }
     242
     243  void CplexSolverBase::_getColCoeffs(int i, InsertIterator b) const {
     244
     245    int tmp1, tmp2, tmp3, length;
     246    CPXgetcols(cplexEnv(), _prob, &tmp1, &tmp2, 0, 0, 0, &length, i, i);
     247
     248    length = -length;
     249    std::vector<int> indices(length);
     250    std::vector<double> values(length);
     251
     252    CPXgetcols(cplexEnv(), _prob, &tmp1, &tmp2,
     253               &indices.front(), &values.front(),
     254               length, &tmp3, i, i);
     255
     256    for (int i = 0; i < length; ++i) {
     257      *b = std::make_pair(indices[i], values[i]);
     258      ++b;
     259    }
     260
     261  }
     262
     263  void CplexSolverBase::_setCoeff(int row, int col, Value value) {
     264    CPXchgcoef(cplexEnv(), _prob, row, col, value);
     265  }
     266
     267  CplexSolverBase::Value CplexSolverBase::_getCoeff(int row, int col) const {
     268    CplexSolverBase::Value value;
     269    CPXgetcoef(cplexEnv(), _prob, row, col, &value);
     270    return value;
     271  }
     272
     273  void CplexSolverBase::_setColLowerBound(int i, Value value) {
     274    const char s = 'L';
     275    CPXchgbds(cplexEnv(), _prob, 1, &i, &s, &value);
     276  }
     277
     278  CplexSolverBase::Value CplexSolverBase::_getColLowerBound(int i) const {
     279    CplexSolverBase::Value res;
     280    CPXgetlb(cplexEnv(), _prob, &res, i, i);
     281    return res <= -CPX_INFBOUND ? -INF : res;
     282  }
     283
     284  void CplexSolverBase::_setColUpperBound(int i, Value value)
     285  {
     286    const char s = 'U';
     287    CPXchgbds(cplexEnv(), _prob, 1, &i, &s, &value);
     288  }
     289
     290  CplexSolverBase::Value CplexSolverBase::_getColUpperBound(int i) const {
     291    CplexSolverBase::Value res;
     292    CPXgetub(cplexEnv(), _prob, &res, i, i);
     293    return res >= CPX_INFBOUND ? INF : res;
     294  }
     295
     296  CplexSolverBase::Value CplexSolverBase::_getRowLowerBound(int i) const {
     297    char s;
     298    CPXgetsense(cplexEnv(), _prob, &s, i, i);
     299    CplexSolverBase::Value res;
     300
     301    switch (s) {
     302    case 'G':
     303    case 'R':
     304    case 'E':
     305      CPXgetrhs(cplexEnv(), _prob, &res, i, i);
     306      return res <= -CPX_INFBOUND ? -INF : res;
     307    default:
     308      return -INF;
     309    }
     310  }
     311
     312  CplexSolverBase::Value CplexSolverBase::_getRowUpperBound(int i) const {
     313    char s;
     314    CPXgetsense(cplexEnv(), _prob, &s, i, i);
     315    CplexSolverBase::Value res;
     316
     317    switch (s) {
     318    case 'L':
     319    case 'E':
     320      CPXgetrhs(cplexEnv(), _prob, &res, i, i);
     321      return res >= CPX_INFBOUND ? INF : res;
     322    case 'R':
     323      CPXgetrhs(cplexEnv(), _prob, &res, i, i);
     324      {
     325        double rng;
     326        CPXgetrngval(cplexEnv(), _prob, &rng, i, i);
     327        res += rng;
     328      }
     329      return res >= CPX_INFBOUND ? INF : res;
     330    default:
     331      return INF;
     332    }
     333  }
     334
     335  //This is easier to implement
     336  void CplexSolverBase::_set_row_bounds(int i, Value lb, Value ub) {
     337    if (lb == -INF) {
     338      const char s = 'L';
     339      CPXchgsense(cplexEnv(), _prob, 1, &i, &s);
     340      CPXchgrhs(cplexEnv(), _prob, 1, &i, &ub);
     341    } else if (ub == INF) {
     342      const char s = 'G';
     343      CPXchgsense(cplexEnv(), _prob, 1, &i, &s);
     344      CPXchgrhs(cplexEnv(), _prob, 1, &i, &lb);
     345    } else if (lb == ub){
     346      const char s = 'E';
     347      CPXchgsense(cplexEnv(), _prob, 1, &i, &s);
     348      CPXchgrhs(cplexEnv(), _prob, 1, &i, &lb);
     349    } else {
     350      const char s = 'R';
     351      CPXchgsense(cplexEnv(), _prob, 1, &i, &s);
     352      CPXchgrhs(cplexEnv(), _prob, 1, &i, &lb);
     353      double len = ub - lb;
     354      CPXchgrngval(cplexEnv(), _prob, 1, &i, &len);
     355    }
     356  }
     357
     358  void CplexSolverBase::_setRowLowerBound(int i, Value lb)
     359  {
     360    LEMON_ASSERT(lb != INF, "Invalid bound");
     361    _set_row_bounds(i, lb, CplexSolverBase::_getRowUpperBound(i));
     362  }
     363
     364  void CplexSolverBase::_setRowUpperBound(int i, Value ub)
     365  {
     366
     367    LEMON_ASSERT(ub != -INF, "Invalid bound");
     368    _set_row_bounds(i, CplexSolverBase::_getRowLowerBound(i), ub);
     369  }
     370
     371  void CplexSolverBase::_setObjCoeffs(ExprIterator b, ExprIterator e)
     372  {
     373    std::vector<int> indices;
     374    std::vector<Value> values;
     375    for(ExprIterator it=b; it!=e; ++it) {
     376      indices.push_back(it->first);
     377      values.push_back(it->second);
     378    }
     379    CPXchgobj(cplexEnv(), _prob, values.size(),
     380              &indices.front(), &values.front());
     381   
     382  }
     383
     384  void CplexSolverBase::_getObjCoeffs(InsertIterator b) const
     385  {
     386    int num = CPXgetnumcols(cplexEnv(), _prob);
     387    std::vector<Value> x(num);
     388   
     389    CPXgetobj(cplexEnv(), _prob, &x.front(), 0, num - 1);
     390    for (int i = 0; i < num; ++i) {
     391      if (x[i] != 0.0) {
     392        *b = std::make_pair(i, x[i]);
     393        ++b;
     394      }
     395    }
     396  }
     397
     398  void CplexSolverBase::_setObjCoeff(int i, SolverBase::Value obj_coef)
     399  {
     400    CPXchgobj(cplexEnv(), _prob, 1, &i, &obj_coef);
     401  }
     402
     403  CplexSolverBase::Value CplexSolverBase::_getObjCoeff(int i) const
     404  {
     405    Value x;
     406    CPXgetobj(cplexEnv(), _prob, &x, i, i);
     407    return x;
     408  }
     409
     410  void CplexSolverBase::_setSense(CplexSolverBase::Sense sense) {
     411    switch (sense) {
     412    case MIN:
     413      CPXchgobjsen(cplexEnv(), _prob, CPX_MIN);
     414      break;
     415    case MAX:
     416      CPXchgobjsen(cplexEnv(), _prob, CPX_MAX);
     417      break;
     418    }
     419  }
     420
     421  CplexSolverBase::Sense CplexSolverBase::_getSense() const {
     422    switch (CPXgetobjsen(cplexEnv(), _prob)) {
     423    case CPX_MIN:
     424      return MIN;
     425    case CPX_MAX:
     426      return MAX;
     427    default:
     428      LEMON_ASSERT(false, "Invalid sense");
     429      return CplexSolverBase::Sense();
     430    }
     431  }
     432
     433  void CplexSolverBase::_clear() {
     434    CPXfreeprob(cplexEnv(),&_prob);
     435    int status;
     436    _prob = CPXcreateprob(cplexEnv(), &status, "Cplex problem");
     437    rows.clear();
     438    cols.clear();
     439  }
     440
     441  // LpCplex members
     442
     443  LpCplex::LpCplex()
     444    : SolverBase(), CplexSolverBase(), LpSolverBase() {}
     445
     446  LpCplex::LpCplex(const CplexEnv& env)
     447    : SolverBase(), CplexSolverBase(env), LpSolverBase() {}
     448
     449  LpCplex::LpCplex(const LpCplex& other)
     450    : SolverBase(), CplexSolverBase(other), LpSolverBase() {}
     451
     452  LpCplex::~LpCplex() {}
     453
     454  LpCplex* LpCplex::_newSolver() const { return new LpCplex; }
     455  LpCplex* LpCplex::_cloneSolver() const {return new LpCplex(*this); }
     456
     457  void LpCplex::_clear_temporals() {
     458    _col_status.clear();
     459    _row_status.clear();
     460    _primal_ray.clear();
     461    _dual_ray.clear();
     462  }
     463
     464  // The routine returns zero unless an error occurred during the
     465  // optimization. Examples of errors include exhausting available
     466  // memory (CPXERR_NO_MEMORY) or encountering invalid data in the
     467  // CPLEX problem object (CPXERR_NO_PROBLEM). Exceeding a
     468  // user-specified CPLEX limit, or proving the model infeasible or
     469  // unbounded, are not considered errors. Note that a zero return
     470  // value does not necessarily mean that a solution exists. Use query
     471  // routines CPXsolninfo, CPXgetstat, and CPXsolution to obtain
     472  // further information about the status of the optimization.
     473  LpCplex::SolveExitStatus LpCplex::convertStatus(int status) {
     474#if CPX_VERSION >= 800
     475    if (status == 0) {
     476      switch (CPXgetstat(cplexEnv(), _prob)) {
     477      case CPX_STAT_OPTIMAL:
     478      case CPX_STAT_INFEASIBLE:
     479      case CPX_STAT_UNBOUNDED:
     480        return SOLVED;
     481      default:
     482        return UNSOLVED;
     483      }
     484    } else {
     485      return UNSOLVED;
     486    }
     487#else
     488    if (status == 0) {
     489      //We want to exclude some cases
     490      switch (CPXgetstat(cplexEnv(), _prob)) {
     491      case CPX_OBJ_LIM:
     492      case CPX_IT_LIM_FEAS:
     493      case CPX_IT_LIM_INFEAS:
     494      case CPX_TIME_LIM_FEAS:
     495      case CPX_TIME_LIM_INFEAS:
     496        return UNSOLVED;
     497      default:
     498        return SOLVED;
     499      }
     500    } else {
     501      return UNSOLVED;
     502    }
     503#endif
     504  }
     505 
     506  LpCplex::SolveExitStatus LpCplex::_solve() {
     507    _clear_temporals();
     508    return convertStatus(CPXlpopt(cplexEnv(), _prob));
     509  }
     510
     511  LpCplex::SolveExitStatus LpCplex::solvePrimal() {
     512    _clear_temporals();
     513    return convertStatus(CPXprimopt(cplexEnv(), _prob));
     514  }
     515
     516  LpCplex::SolveExitStatus LpCplex::solveDual() {
     517    _clear_temporals();
     518    return convertStatus(CPXdualopt(cplexEnv(), _prob));
     519  }
     520
     521  LpCplex::SolveExitStatus LpCplex::solveBarrier() {
     522    _clear_temporals();
     523    return convertStatus(CPXbaropt(cplexEnv(), _prob));
     524  }
     525
     526  LpCplex::Value LpCplex::_getPrimal(int i) const {
     527    Value x;
     528    CPXgetx(cplexEnv(), _prob, &x, i, i);
     529    return x;
     530  }
     531
     532  LpCplex::Value LpCplex::_getDual(int i) const {
     533    Value y;
     534    CPXgetpi(cplexEnv(), _prob, &y, i, i);
     535    return y;
     536  }
     537
     538  LpCplex::Value LpCplex::_getPrimalValue() const {
     539    Value objval;
     540    CPXgetobjval(cplexEnv(), _prob, &objval);
     541    return objval;
     542  }
     543
     544  LpCplex::VarStatus LpCplex::_getColStatus(int i) const {
     545    if (_col_status.empty()) {
     546      _col_status.resize(CPXgetnumcols(cplexEnv(), _prob));
     547      CPXgetbase(cplexEnv(), _prob, &_col_status.front(), 0);
     548    }
     549    switch (_col_status[i]) {
     550    case CPX_BASIC:
     551      return BASIC;
     552    case CPX_FREE_SUPER:
     553      return FREE;
     554    case CPX_AT_LOWER:
     555      return LOWER;
     556    case CPX_AT_UPPER:
     557      return UPPER;
     558    default:
     559      LEMON_ASSERT(false, "Wrong column status");
     560      return LpCplex::VarStatus();
     561    }
     562  }
     563
     564  LpCplex::VarStatus LpCplex::_getRowStatus(int i) const {
     565    if (_row_status.empty()) {
     566      _row_status.resize(CPXgetnumrows(cplexEnv(), _prob));
     567      CPXgetbase(cplexEnv(), _prob, 0, &_row_status.front());
     568    }
     569    switch (_row_status[i]) {
     570    case CPX_BASIC:
     571      return BASIC;
     572    case CPX_AT_LOWER:
     573      {
     574        char s;
     575        CPXgetsense(cplexEnv(), _prob, &s, i, i);
     576        return s != 'L' ? LOWER : UPPER;
     577      }
     578    case CPX_AT_UPPER:
     579      return UPPER;
     580    default:
     581      LEMON_ASSERT(false, "Wrong row status");
     582      return LpCplex::VarStatus();
     583    }
     584  }
     585
     586  LpCplex::Value LpCplex::_getPrimalRay(int i) const {
     587    if (_primal_ray.empty()) {
     588      _primal_ray.resize(CPXgetnumcols(cplexEnv(), _prob));
     589      CPXgetray(cplexEnv(), _prob, &_primal_ray.front());
     590    }
     591    return _primal_ray[i];
     592  }
     593
     594  LpCplex::Value LpCplex::_getDualRay(int i) const {
     595    if (_dual_ray.empty()) {
     596     
     597    }
     598    return _dual_ray[i];
     599  }
     600
     601  //7.5-os cplex statusai (Vigyazat: a 9.0-asei masok!)
     602  // This table lists the statuses, returned by the CPXgetstat()
     603  // routine, for solutions to LP problems or mixed integer problems. If
     604  // no solution exists, the return value is zero.
     605
     606  // For Simplex, Barrier
     607  // 1          CPX_OPTIMAL
     608  //          Optimal solution found
     609  // 2          CPX_INFEASIBLE
     610  //          Problem infeasible
     611  // 3    CPX_UNBOUNDED
     612  //          Problem unbounded
     613  // 4          CPX_OBJ_LIM
     614  //          Objective limit exceeded in Phase II
     615  // 5          CPX_IT_LIM_FEAS
     616  //          Iteration limit exceeded in Phase II
     617  // 6          CPX_IT_LIM_INFEAS
     618  //          Iteration limit exceeded in Phase I
     619  // 7          CPX_TIME_LIM_FEAS
     620  //          Time limit exceeded in Phase II
     621  // 8          CPX_TIME_LIM_INFEAS
     622  //          Time limit exceeded in Phase I
     623  // 9          CPX_NUM_BEST_FEAS
     624  //          Problem non-optimal, singularities in Phase II
     625  // 10         CPX_NUM_BEST_INFEAS
     626  //          Problem non-optimal, singularities in Phase I
     627  // 11         CPX_OPTIMAL_INFEAS
     628  //          Optimal solution found, unscaled infeasibilities
     629  // 12         CPX_ABORT_FEAS
     630  //          Aborted in Phase II
     631  // 13         CPX_ABORT_INFEAS
     632  //          Aborted in Phase I
     633  // 14          CPX_ABORT_DUAL_INFEAS
     634  //          Aborted in barrier, dual infeasible
     635  // 15          CPX_ABORT_PRIM_INFEAS
     636  //          Aborted in barrier, primal infeasible
     637  // 16          CPX_ABORT_PRIM_DUAL_INFEAS
     638  //          Aborted in barrier, primal and dual infeasible
     639  // 17          CPX_ABORT_PRIM_DUAL_FEAS
     640  //          Aborted in barrier, primal and dual feasible
     641  // 18          CPX_ABORT_CROSSOVER
     642  //          Aborted in crossover
     643  // 19          CPX_INForUNBD
     644  //          Infeasible or unbounded
     645  // 20   CPX_PIVOT
     646  //       User pivot used
     647  //
     648  //     Ezeket hova tegyem:
     649  // ??case CPX_ABORT_DUAL_INFEAS
     650  // ??case CPX_ABORT_CROSSOVER
     651  // ??case CPX_INForUNBD
     652  // ??case CPX_PIVOT
     653
     654  //Some more interesting stuff:
     655
     656  // CPX_PARAM_PROBMETHOD  1062  int  LPMETHOD
     657  // 0 Automatic
     658  // 1 Primal Simplex
     659  // 2 Dual Simplex
     660  // 3 Network Simplex
     661  // 4 Standard Barrier
     662  // Default: 0
     663  // Description: Method for linear optimization.
     664  // Determines which algorithm is used when CPXlpopt() (or "optimize"
     665  // in the Interactive Optimizer) is called. Currently the behavior of
     666  // the "Automatic" setting is that CPLEX simply invokes the dual
     667  // simplex method, but this capability may be expanded in the future
     668  // so that CPLEX chooses the method based on problem characteristics
     669#if CPX_VERSION < 900
     670  void statusSwitch(CPXENVptr cplexEnv(),int& stat){
     671    int lpmethod;
     672    CPXgetintparam (cplexEnv(),CPX_PARAM_PROBMETHOD,&lpmethod);
     673    if (lpmethod==2){
     674      if (stat==CPX_UNBOUNDED){
     675        stat=CPX_INFEASIBLE;
     676      }
     677      else{
     678        if (stat==CPX_INFEASIBLE)
     679          stat=CPX_UNBOUNDED;
     680      }
     681    }
     682  }
     683#else
     684  void statusSwitch(CPXENVptr,int&){}
     685#endif
     686
     687  LpCplex::ProblemType LpCplex::_getPrimalType() const {
     688    // Unboundedness not treated well: the following is from cplex 9.0 doc
     689    // About Unboundedness
     690
     691    // The treatment of models that are unbounded involves a few
     692    // subtleties. Specifically, a declaration of unboundedness means that
     693    // ILOG CPLEX has determined that the model has an unbounded
     694    // ray. Given any feasible solution x with objective z, a multiple of
     695    // the unbounded ray can be added to x to give a feasible solution
     696    // with objective z-1 (or z+1 for maximization models). Thus, if a
     697    // feasible solution exists, then the optimal objective is
     698    // unbounded. Note that ILOG CPLEX has not necessarily concluded that
     699    // a feasible solution exists. Users can call the routine CPXsolninfo
     700    // to determine whether ILOG CPLEX has also concluded that the model
     701    // has a feasible solution.
     702
     703    int stat = CPXgetstat(cplexEnv(), _prob);
     704#if CPX_VERSION >= 800
     705    switch (stat)
     706      {
     707      case CPX_STAT_OPTIMAL:
     708        return OPTIMAL;
     709      case CPX_STAT_UNBOUNDED:
     710        return UNBOUNDED;
     711      case CPX_STAT_INFEASIBLE:
     712        return INFEASIBLE;
     713      default:
     714        return UNDEFINED;
     715      }
     716#else
     717    statusSwitch(cplexEnv(),stat);
     718    //CPXgetstat(cplexEnv(), _prob);
     719    //printf("A primal status: %d, CPX_OPTIMAL=%d \n",stat,CPX_OPTIMAL);
     720    switch (stat) {
     721    case 0:
     722      return UNDEFINED; //Undefined
     723    case CPX_OPTIMAL://Optimal
     724      return OPTIMAL;
     725    case CPX_UNBOUNDED://Unbounded
     726      return INFEASIBLE;//In case of dual simplex
     727      //return UNBOUNDED;
     728    case CPX_INFEASIBLE://Infeasible
     729      //    case CPX_IT_LIM_INFEAS:
     730      //     case CPX_TIME_LIM_INFEAS:
     731      //     case CPX_NUM_BEST_INFEAS:
     732      //     case CPX_OPTIMAL_INFEAS:
     733      //     case CPX_ABORT_INFEAS:
     734      //     case CPX_ABORT_PRIM_INFEAS:
     735      //     case CPX_ABORT_PRIM_DUAL_INFEAS:
     736      return UNBOUNDED;//In case of dual simplex
     737      //return INFEASIBLE;
     738      //     case CPX_OBJ_LIM:
     739      //     case CPX_IT_LIM_FEAS:
     740      //     case CPX_TIME_LIM_FEAS:
     741      //     case CPX_NUM_BEST_FEAS:
     742      //     case CPX_ABORT_FEAS:
     743      //     case CPX_ABORT_PRIM_DUAL_FEAS:
     744      //       return FEASIBLE;
     745    default:
     746      return UNDEFINED; //Everything else comes here
     747      //FIXME error
     748    }
     749#endif
     750  }
     751
     752  //9.0-as cplex verzio statusai
     753  // CPX_STAT_ABORT_DUAL_OBJ_LIM
     754  // CPX_STAT_ABORT_IT_LIM
     755  // CPX_STAT_ABORT_OBJ_LIM
     756  // CPX_STAT_ABORT_PRIM_OBJ_LIM
     757  // CPX_STAT_ABORT_TIME_LIM
     758  // CPX_STAT_ABORT_USER
     759  // CPX_STAT_FEASIBLE_RELAXED
     760  // CPX_STAT_INFEASIBLE
     761  // CPX_STAT_INForUNBD
     762  // CPX_STAT_NUM_BEST
     763  // CPX_STAT_OPTIMAL
     764  // CPX_STAT_OPTIMAL_FACE_UNBOUNDED
     765  // CPX_STAT_OPTIMAL_INFEAS
     766  // CPX_STAT_OPTIMAL_RELAXED
     767  // CPX_STAT_UNBOUNDED
     768
     769  LpCplex::ProblemType LpCplex::_getDualType() const {
     770    int stat = CPXgetstat(cplexEnv(), _prob);
     771#if CPX_VERSION >= 800
     772    switch (stat) {
     773    case CPX_STAT_OPTIMAL:
     774      return OPTIMAL;
     775    case CPX_STAT_UNBOUNDED:
     776      return INFEASIBLE;
     777    default:
     778      return UNDEFINED;
     779    }
     780#else
     781    statusSwitch(cplexEnv(),stat);
     782    switch (stat) {
     783    case 0:
     784      return UNDEFINED; //Undefined
     785    case CPX_OPTIMAL://Optimal
     786      return OPTIMAL;
     787    case CPX_UNBOUNDED:
     788      return INFEASIBLE;
     789    default:
     790      return UNDEFINED; //Everything else comes here
     791      //FIXME error
     792    }
     793#endif
     794  }
     795
     796  // MipCplex members
     797
     798  MipCplex::MipCplex()
     799    : SolverBase(), CplexSolverBase(), MipSolverBase() {
     800
     801#if CPX_VERSION < 800
     802    CPXchgprobtype(cplexEnv(),  _prob, CPXPROB_MIP);
     803#else
     804    CPXchgprobtype(cplexEnv(),  _prob, CPXPROB_MILP);
     805#endif
     806  }
     807
     808  MipCplex::MipCplex(const CplexEnv& env)
     809    : SolverBase(), CplexSolverBase(env), MipSolverBase() {
     810
     811#if CPX_VERSION < 800
     812    CPXchgprobtype(cplexEnv(),  _prob, CPXPROB_MIP);
     813#else
     814    CPXchgprobtype(cplexEnv(),  _prob, CPXPROB_MILP);
     815#endif
     816
     817  }
     818
     819  MipCplex::MipCplex(const MipCplex& other)
     820    : SolverBase(), CplexSolverBase(other), MipSolverBase() {}
     821
     822  MipCplex::~MipCplex() {}
     823
     824  MipCplex* MipCplex::_newSolver() const { return new MipCplex; }
     825  MipCplex* MipCplex::_cloneSolver() const {return new MipCplex(*this); }
     826
     827  void MipCplex::_setColType(int i, MipCplex::ColTypes col_type) {
     828
     829    // Note If a variable is to be changed to binary, a call to CPXchgbds
     830    // should also be made to change the bounds to 0 and 1.
     831
     832    switch (col_type){
     833    case INTEGER: {
     834      const char t = 'I';
     835      CPXchgctype (cplexEnv(), _prob, 1, &i, &t);
     836    } break;
     837    case REAL: {
     838      const char t = 'C';
     839      CPXchgctype (cplexEnv(), _prob, 1, &i, &t);
     840    } break;
     841    default:
     842      break;
     843    }
     844  }
     845
     846  MipSolverBase::ColTypes MipCplex::_getColType(int i) const {
     847    char t;
     848    CPXgetctype (cplexEnv(), _prob, &t, i, i);
     849    switch (t) {
     850    case 'I':
     851      return INTEGER;
     852    case 'C':
     853      return REAL;
     854    default:
     855      LEMON_ASSERT(false, "Invalid column type");
     856      return ColTypes();
     857    }
     858
     859  }
     860
     861  MipCplex::SolveExitStatus MipCplex::_solve() {
     862    int status;
     863    status = CPXmipopt (cplexEnv(), _prob);
     864    if (status==0)
     865      return SOLVED;
     866    else
     867      return UNSOLVED;
     868
     869  }
     870
     871
     872  MipCplex::ProblemType MipCplex::_getType() const {
     873
     874    int stat = CPXgetstat(cplexEnv(), _prob);
     875
     876    //Fortunately, MIP statuses did not change for cplex 8.0
     877    switch (stat) {
     878    case CPXMIP_OPTIMAL:
     879      // Optimal integer solution has been found.
     880    case CPXMIP_OPTIMAL_TOL:
     881      // Optimal soluton with the tolerance defined by epgap or epagap has
     882      // been found.
     883      return OPTIMAL;
     884      //This also exists in later issues
     885      //    case CPXMIP_UNBOUNDED:
     886      //return UNBOUNDED;
     887      case CPXMIP_INFEASIBLE:
     888        return INFEASIBLE;
     889    default:
     890      return UNDEFINED;
     891    }
     892    //Unboundedness not treated well: the following is from cplex 9.0 doc
     893    // About Unboundedness
     894
     895    // The treatment of models that are unbounded involves a few
     896    // subtleties. Specifically, a declaration of unboundedness means that
     897    // ILOG CPLEX has determined that the model has an unbounded
     898    // ray. Given any feasible solution x with objective z, a multiple of
     899    // the unbounded ray can be added to x to give a feasible solution
     900    // with objective z-1 (or z+1 for maximization models). Thus, if a
     901    // feasible solution exists, then the optimal objective is
     902    // unbounded. Note that ILOG CPLEX has not necessarily concluded that
     903    // a feasible solution exists. Users can call the routine CPXsolninfo
     904    // to determine whether ILOG CPLEX has also concluded that the model
     905    // has a feasible solution.
     906  }
     907
     908  MipCplex::Value MipCplex::_getSol(int i) const {
     909    Value x;
     910    CPXgetmipx(cplexEnv(), _prob, &x, i, i);
     911    return x;
     912  }
     913
     914  MipCplex::Value MipCplex::_getSolValue() const {
     915    Value objval;
     916    CPXgetmipobjval(cplexEnv(), _prob, &objval);
     917    return objval;
     918  }
     919
     920} //namespace lemon
     921
  • new file lemon/lp_cplex.h

    diff -r cbe3ec2d59d2 -r aa6dff56b2e1 lemon/lp_cplex.h
    - +  
     1/* -*- C++ -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library
     4 *
     5 * Copyright (C) 2003-2008
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#ifndef LEMON_LP_CPLEX_H
     20#define LEMON_LP_CPLEX_H
     21
     22///\file
     23///\brief Header of the LEMON-CPLEX lp solver interface.
     24
     25#include <lemon/lp_base.h>
     26
     27struct cpxenv;
     28struct cpxlp;
     29
     30namespace lemon {
     31
     32  /// \brief Reference counted wrapper around cpxenv pointer
     33  ///
     34  /// The cplex uses environment object which is responsible for
     35  /// checking the proper license usage. This class provides a simple
     36  /// interface for share the environment object between different
     37  /// problems.
     38  class CplexEnv {
     39    friend class CplexSolverBase;
     40  private:
     41    cpxenv* _env;
     42    mutable int* _cnt;
     43
     44  public:
     45
     46    /// \brief This exception is thrown when the license check is not
     47    /// sufficient
     48    class LicenseError : public Exception {
     49      friend class CplexEnv;
     50    private:
     51
     52      LicenseError(int status);
     53      char _message[510];     
     54
     55    public:
     56     
     57      /// The short error message
     58      virtual const char* what() const throw() {
     59        return _message;
     60      }
     61    };
     62
     63    /// Constructor
     64    CplexEnv();
     65    /// Shallow copy constructor
     66    CplexEnv(const CplexEnv&);
     67    /// Shallow assignement
     68    CplexEnv& operator=(const CplexEnv&);
     69    /// Destructor
     70    virtual ~CplexEnv();
     71
     72  protected:
     73   
     74    cpxenv* cplexEnv() { return _env; }
     75    const cpxenv* cplexEnv() const { return _env; }
     76  };
     77
     78  /// \brief Base interface for the CPLEX LP and MIP solver
     79  ///
     80  /// This class implements the common interface of the CPLEX LP and MIP solver.
     81  ///\ingroup lp_group
     82  class CplexSolverBase : virtual public SolverBase {
     83  protected:
     84
     85    typedef SolverBase Parent;
     86   
     87    CplexEnv _env;
     88    cpxlp* _prob;
     89
     90    CplexSolverBase();
     91    CplexSolverBase(const CplexEnv&);
     92    CplexSolverBase(const CplexSolverBase &);
     93    virtual ~CplexSolverBase();
     94
     95    virtual int _addCol();
     96    virtual int _addRow();
     97
     98    virtual void _eraseCol(int i);
     99    virtual void _eraseRow(int i);
     100
     101    virtual void _eraseColId(int i);
     102    virtual void _eraseRowId(int i);
     103
     104    virtual void _getColName(int col, std::string& name) const;
     105    virtual void _setColName(int col, const std::string& name);
     106    virtual int _colByName(const std::string& name) const;
     107
     108    virtual void _getRowName(int row, std::string& name) const;
     109    virtual void _setRowName(int row, const std::string& name);
     110    virtual int _rowByName(const std::string& name) const;
     111
     112    virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e);
     113    virtual void _getRowCoeffs(int i, InsertIterator b) const;
     114
     115    virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e);
     116    virtual void _getColCoeffs(int i, InsertIterator b) const;
     117
     118    virtual void _setCoeff(int row, int col, Value value);
     119    virtual Value _getCoeff(int row, int col) const;
     120
     121    virtual void _setColLowerBound(int i, Value value);
     122    virtual Value _getColLowerBound(int i) const;
     123
     124    virtual void _setColUpperBound(int i, Value value);
     125    virtual Value _getColUpperBound(int i) const;
     126
     127  private:
     128    void _set_row_bounds(int i, Value lb, Value ub);
     129  protected:
     130
     131    virtual void _setRowLowerBound(int i, Value value);
     132    virtual Value _getRowLowerBound(int i) const;
     133
     134    virtual void _setRowUpperBound(int i, Value value);
     135    virtual Value _getRowUpperBound(int i) const;
     136
     137    virtual void _setObjCoeffs(ExprIterator b, ExprIterator e);
     138    virtual void _getObjCoeffs(InsertIterator b) const;
     139
     140    virtual void _setObjCoeff(int i, Value obj_coef);
     141    virtual Value _getObjCoeff(int i) const;
     142   
     143    virtual void _setSense(Sense sense);
     144    virtual Sense _getSense() const;
     145
     146    virtual void _clear();
     147
     148  public:
     149
     150    /// Returns the used \c CplexEnv instance
     151    const CplexEnv& env() const { return _env; }
     152    ///
     153    const cpxenv* cplexEnv() const { return _env.cplexEnv(); }
     154
     155    cpxlp* cplexLp() { return _prob; }
     156    const cpxlp* cplexLp() const { return _prob; }
     157
     158  };
     159
     160  /// \brief Interface for the CPLEX LP solver
     161  ///
     162  /// This class implements an interface for the CPLEX LP solver.
     163  ///\ingroup lp_group
     164  class LpCplex : public CplexSolverBase, public LpSolverBase {
     165  public:
     166    /// \e
     167    LpCplex();
     168    /// \e
     169    LpCplex(const CplexEnv&);
     170    /// \e
     171    LpCplex(const LpCplex&);
     172    /// \e
     173    virtual ~LpCplex();
     174
     175  private:
     176
     177    // these values cannot retrieved element by element
     178    mutable std::vector<int> _col_status;
     179    mutable std::vector<int> _row_status;
     180
     181    mutable std::vector<Value> _primal_ray;
     182    mutable std::vector<Value> _dual_ray;
     183
     184    void _clear_temporals();
     185
     186    SolveExitStatus convertStatus(int status);
     187
     188  protected:
     189
     190    virtual LpCplex* _cloneSolver() const;
     191    virtual LpCplex* _newSolver() const;
     192
     193    virtual SolveExitStatus _solve();
     194    virtual Value _getPrimal(int i) const;
     195    virtual Value _getDual(int i) const;
     196    virtual Value _getPrimalValue() const;
     197
     198    virtual VarStatus _getColStatus(int i) const;
     199    virtual VarStatus _getRowStatus(int i) const;
     200
     201    virtual Value _getPrimalRay(int i) const;
     202    virtual Value _getDualRay(int i) const;
     203   
     204    virtual ProblemType _getPrimalType() const;
     205    virtual ProblemType _getDualType() const;
     206
     207  public:
     208   
     209    /// Solve with primal simplex method
     210    SolveExitStatus solvePrimal();
     211
     212    /// Solve with dual simplex method
     213    SolveExitStatus solveDual();
     214
     215    /// Solve with barrier method
     216    SolveExitStatus solveBarrier();
     217
     218  };
     219
     220  /// \brief Interface for the CPLEX MIP solver
     221  ///
     222  /// This class implements an interface for the CPLEX MIP solver.
     223  ///\ingroup lp_group
     224  class MipCplex : public CplexSolverBase, public MipSolverBase {
     225  public:
     226    /// \e
     227    MipCplex();
     228    /// \e
     229    MipCplex(const CplexEnv&);
     230    /// \e
     231    MipCplex(const MipCplex&);
     232    /// \e
     233    virtual ~MipCplex();
     234
     235  protected:
     236
     237    virtual MipCplex* _cloneSolver() const;
     238    virtual MipCplex* _newSolver() const;
     239
     240    virtual ColTypes _getColType(int col) const;
     241    virtual void _setColType(int col, ColTypes col_type);
     242   
     243    virtual SolveExitStatus _solve();
     244    virtual ProblemType _getType() const;
     245    virtual Value _getSol(int i) const;
     246    virtual Value _getSolValue() const;
     247   
     248  };
     249
     250} //END OF NAMESPACE LEMON
     251
     252#endif //LEMON_LP_CPLEX_H
     253
  • new file lemon/lp_glpk.cc

    diff -r cbe3ec2d59d2 -r aa6dff56b2e1 lemon/lp_glpk.cc
    - +  
     1/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library.
     4 *
     5 * Copyright (C) 2003-2008
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19///\file
     20///\brief Implementation of the LEMON GLPK LP and MIP solver interface.
     21
     22#include <lemon/lp_glpk.h>
     23#include <glpk.h>
     24
     25#include <lemon/assert.h>
     26
     27namespace lemon {
     28
     29  // GlpkSolverBase members
     30
     31  GlpkSolverBase::GlpkSolverBase() : SolverBase() {
     32    lp = glp_create_prob();
     33    glp_create_index(lp);
     34  }
     35
     36  GlpkSolverBase::GlpkSolverBase(const GlpkSolverBase &other) : SolverBase() {
     37    lp = glp_create_prob();
     38    glp_copy_prob(lp, other.lp, GLP_ON);
     39    glp_create_index(lp);
     40    rows = other.rows;
     41    cols = other.cols;
     42  }
     43
     44  GlpkSolverBase::~GlpkSolverBase() {
     45    glp_delete_prob(lp);
     46  }
     47
     48  int GlpkSolverBase::_addCol() {
     49    int i = glp_add_cols(lp, 1);
     50    glp_set_col_bnds(lp, i, GLP_FR, 0.0, 0.0);
     51    return i;
     52  }
     53
     54  int GlpkSolverBase::_addRow() {
     55    int i = glp_add_rows(lp, 1);
     56    glp_set_row_bnds(lp, i, GLP_FR, 0.0, 0.0);
     57    return i;
     58  }
     59
     60  void GlpkSolverBase::_eraseCol(int i) {
     61    int ca[2];
     62    ca[1] = i;
     63    glp_del_cols(lp, 1, ca);
     64  }
     65
     66  void GlpkSolverBase::_eraseRow(int i) {
     67    int ra[2];
     68    ra[1] = i;
     69    glp_del_rows(lp, 1, ra);
     70  }
     71
     72  void GlpkSolverBase::_eraseColId(int i) {
     73    cols.eraseIndex(i);
     74    cols.shiftIndexes(i);
     75  }
     76
     77  void GlpkSolverBase::_eraseRowId(int i) {
     78    rows.eraseIndex(i);
     79    rows.shiftIndexes(i);
     80  }
     81
     82  void GlpkSolverBase::_getColName(int c, std::string& name) const {
     83    const char *str = glp_get_col_name(lp, c);
     84    if (str) name = str;
     85    else name.clear();
     86  }
     87
     88  void GlpkSolverBase::_setColName(int c, const std::string & name) {
     89    glp_set_col_name(lp, c, const_cast<char*>(name.c_str()));
     90
     91  }
     92
     93  int GlpkSolverBase::_colByName(const std::string& name) const {
     94    int k = glp_find_col(lp, const_cast<char*>(name.c_str()));
     95    return k > 0 ? k : -1;
     96  }
     97
     98  void GlpkSolverBase::_getRowName(int r, std::string& name) const {
     99    const char *str = glp_get_row_name(lp, r);
     100    if (str) name = str;
     101    else name.clear();
     102  }
     103
     104  void GlpkSolverBase::_setRowName(int r, const std::string & name) {
     105    glp_set_row_name(lp, r, const_cast<char*>(name.c_str()));
     106
     107  }
     108
     109  int GlpkSolverBase::_rowByName(const std::string& name) const {
     110    int k = glp_find_row(lp, const_cast<char*>(name.c_str()));
     111    return k > 0 ? k : -1;
     112  }
     113
     114  void GlpkSolverBase::_setRowCoeffs(int i, ExprIterator b, ExprIterator e) {
     115    std::vector<int> indexes;
     116    std::vector<Value> values;
     117
     118    indexes.push_back(0);
     119    values.push_back(0);
     120
     121    for(ExprIterator it = b; it != e; ++it) {
     122      indexes.push_back(it->first);
     123      values.push_back(it->second);
     124    }
     125
     126    glp_set_mat_row(lp, i, values.size() - 1,
     127                    &indexes.front(), &values.front());
     128  }
     129
     130  void GlpkSolverBase::_getRowCoeffs(int ix, InsertIterator b) const {
     131    int length = glp_get_mat_row(lp, ix, 0, 0);
     132
     133    std::vector<int> indexes(length + 1);
     134    std::vector<Value> values(length + 1);
     135
     136    glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
     137
     138    for (int i = 1; i <= length; ++i) {
     139      *b = std::make_pair(indexes[i], values[i]);
     140      ++b;
     141    }
     142  }
     143
     144  void GlpkSolverBase::_setColCoeffs(int ix, ExprIterator b,
     145                                     ExprIterator e) {
     146
     147    std::vector<int> indexes;
     148    std::vector<Value> values;
     149
     150    indexes.push_back(0);
     151    values.push_back(0);
     152
     153    for(ExprIterator it = b; it != e; ++it) {
     154      indexes.push_back(it->first);
     155      values.push_back(it->second);
     156    }
     157
     158    glp_set_mat_col(lp, ix, values.size() - 1,
     159                    &indexes.front(), &values.front());
     160  }
     161
     162  void GlpkSolverBase::_getColCoeffs(int ix, InsertIterator b) const {
     163    int length = glp_get_mat_col(lp, ix, 0, 0);
     164
     165    std::vector<int> indexes(length + 1);
     166    std::vector<Value> values(length + 1);
     167
     168    glp_get_mat_col(lp, ix, &indexes.front(), &values.front());
     169
     170    for (int i = 1; i  <= length; ++i) {
     171      *b = std::make_pair(indexes[i], values[i]);
     172      ++b;
     173    }
     174  }
     175
     176  void GlpkSolverBase::_setCoeff(int ix, int jx, Value value) {
     177
     178    if (glp_get_num_cols(lp) < glp_get_num_rows(lp)) {
     179
     180      int length = glp_get_mat_row(lp, ix, 0, 0);
     181
     182      std::vector<int> indexes(length + 2);
     183      std::vector<Value> values(length + 2);
     184
     185      glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
     186
     187      //The following code does not suppose that the elements of the
     188      //array indexes are sorted
     189      bool found = false;
     190      for (int i = 1; i  <= length; ++i) {
     191        if (indexes[i] == jx) {
     192          found = true;
     193          values[i] = value;
     194          break;
     195        }
     196      }
     197      if (!found) {
     198        ++length;
     199        indexes[length] = jx;
     200        values[length] = value;
     201      }
     202
     203      glp_set_mat_row(lp, ix, length, &indexes.front(), &values.front());
     204
     205    } else {
     206
     207      int length = glp_get_mat_col(lp, jx, 0, 0);
     208
     209      std::vector<int> indexes(length + 2);
     210      std::vector<Value> values(length + 2);
     211
     212      glp_get_mat_col(lp, jx, &indexes.front(), &values.front());
     213
     214      //The following code does not suppose that the elements of the
     215      //array indexes are sorted
     216      bool found = false;
     217      for (int i = 1; i <= length; ++i) {
     218        if (indexes[i] == ix) {
     219          found = true;
     220          values[i] = value;
     221          break;
     222        }
     223      }
     224      if (!found) {
     225        ++length;
     226        indexes[length] = ix;
     227        values[length] = value;
     228      }
     229
     230      glp_set_mat_col(lp, jx, length, &indexes.front(), &values.front());
     231    }
     232
     233  }
     234
     235  SolverBase::Value GlpkSolverBase::_getCoeff(int ix, int jx) const {
     236
     237    int length = glp_get_mat_row(lp, ix, 0, 0);
     238
     239    std::vector<int> indexes(length + 1);
     240    std::vector<Value> values(length + 1);
     241
     242    glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
     243
     244    for (int i = 1; i  <= length; ++i) {
     245      if (indexes[i] == jx) {
     246        return values[i];
     247      }
     248    }
     249
     250    return 0;
     251  }
     252
     253  void GlpkSolverBase::_setColLowerBound(int i, Value lo) {
     254    LEMON_ASSERT(lo != INF, "Invalid bound");
     255
     256    int b = glp_get_col_type(lp, i);
     257    double up = glp_get_col_ub(lp, i);
     258    if (lo == -INF) {
     259      switch (b) {
     260      case GLP_FR:
     261      case GLP_LO:
     262        glp_set_col_bnds(lp, i, GLP_FR, lo, up);
     263        break;
     264      case GLP_UP:
     265        break;
     266      case GLP_DB:
     267      case GLP_FX:
     268        glp_set_col_bnds(lp, i, GLP_UP, lo, up);
     269        break;
     270      default:
     271        break;
     272      }
     273    } else {
     274      switch (b) {
     275      case GLP_FR:
     276      case GLP_LO:
     277        glp_set_col_bnds(lp, i, GLP_LO, lo, up);
     278        break;
     279      case GLP_UP:
     280      case GLP_DB:
     281      case GLP_FX:
     282        if (lo == up)
     283          glp_set_col_bnds(lp, i, GLP_FX, lo, up);
     284        else
     285          glp_set_col_bnds(lp, i, GLP_DB, lo, up);
     286        break;
     287      default:
     288        break;
     289      }
     290    }
     291  }
     292
     293  SolverBase::Value GlpkSolverBase::_getColLowerBound(int i) const {
     294    int b = glp_get_col_type(lp, i);
     295    switch (b) {
     296    case GLP_LO:
     297    case GLP_DB:
     298    case GLP_FX:
     299      return glp_get_col_lb(lp, i);
     300    default:
     301      return -INF;
     302    }
     303  }
     304
     305  void GlpkSolverBase::_setColUpperBound(int i, Value up) {
     306    LEMON_ASSERT(up != -INF, "Invalid bound");
     307
     308    int b = glp_get_col_type(lp, i);
     309    double lo = glp_get_col_lb(lp, i);
     310    if (up == INF) {
     311      switch (b) {
     312      case GLP_FR:
     313      case GLP_LO:
     314        break;
     315      case GLP_UP:
     316        glp_set_col_bnds(lp, i, GLP_FR, lo, up);
     317        break;
     318      case GLP_DB:
     319      case GLP_FX:
     320        glp_set_col_bnds(lp, i, GLP_LO, lo, up);
     321        break;
     322      default:
     323        break;
     324      }
     325    } else {
     326      switch (b) {
     327      case GLP_FR:
     328        glp_set_col_bnds(lp, i, GLP_UP, lo, up);
     329        break;
     330      case GLP_UP:
     331        glp_set_col_bnds(lp, i, GLP_UP, lo, up);
     332        break;
     333      case GLP_LO:
     334      case GLP_DB:
     335      case GLP_FX:
     336        if (lo == up)
     337          glp_set_col_bnds(lp, i, GLP_FX, lo, up);
     338        else
     339          glp_set_col_bnds(lp, i, GLP_DB, lo, up);
     340        break;
     341      default:
     342        break;
     343      }
     344    }
     345
     346  }
     347
     348  SolverBase::Value GlpkSolverBase::_getColUpperBound(int i) const {
     349    int b = glp_get_col_type(lp, i);
     350      switch (b) {
     351      case GLP_UP:
     352      case GLP_DB:
     353      case GLP_FX:
     354        return glp_get_col_ub(lp, i);
     355      default:
     356        return INF;
     357      }
     358  }
     359
     360  void GlpkSolverBase::_setRowLowerBound(int i, Value lo) {
     361    LEMON_ASSERT(lo != INF, "Invalid bound");
     362
     363    int b = glp_get_row_type(lp, i);
     364    double up = glp_get_row_ub(lp, i);
     365    if (lo == -INF) {
     366      switch (b) {
     367      case GLP_FR:
     368      case GLP_LO:
     369        glp_set_row_bnds(lp, i, GLP_FR, lo, up);
     370        break;
     371      case GLP_UP:
     372        break;
     373      case GLP_DB:
     374      case GLP_FX:
     375        glp_set_row_bnds(lp, i, GLP_UP, lo, up);
     376        break;
     377      default:
     378        break;
     379      }
     380    } else {
     381      switch (b) {
     382      case GLP_FR:
     383      case GLP_LO:
     384        glp_set_row_bnds(lp, i, GLP_LO, lo, up);
     385        break;
     386      case GLP_UP:
     387      case GLP_DB:
     388      case GLP_FX:
     389        if (lo == up)
     390          glp_set_row_bnds(lp, i, GLP_FX, lo, up);
     391        else
     392          glp_set_row_bnds(lp, i, GLP_DB, lo, up);
     393        break;
     394      default:
     395        break;
     396      }
     397    }
     398
     399  }
     400
     401  SolverBase::Value GlpkSolverBase::_getRowLowerBound(int i) const {
     402    int b = glp_get_row_type(lp, i);
     403    switch (b) {
     404    case GLP_LO:
     405    case GLP_DB:
     406    case GLP_FX:
     407      return glp_get_row_lb(lp, i);
     408    default:
     409      return -INF;
     410    }
     411  }
     412
     413  void GlpkSolverBase::_setRowUpperBound(int i, Value up) {
     414    LEMON_ASSERT(up != -INF, "Invalid bound");
     415
     416    int b = glp_get_row_type(lp, i);
     417    double lo = glp_get_row_lb(lp, i);
     418    if (up == INF) {
     419      switch (b) {
     420      case GLP_FR:
     421      case GLP_LO:
     422        break;
     423      case GLP_UP:
     424        glp_set_row_bnds(lp, i, GLP_FR, lo, up);
     425        break;
     426      case GLP_DB:
     427      case GLP_FX:
     428        glp_set_row_bnds(lp, i, GLP_LO, lo, up);
     429        break;
     430      default:
     431        break;
     432      }
     433    } else {
     434      switch (b) {
     435      case GLP_FR:
     436        glp_set_row_bnds(lp, i, GLP_UP, lo, up);
     437        break;
     438      case GLP_UP:
     439        glp_set_row_bnds(lp, i, GLP_UP, lo, up);
     440        break;
     441      case GLP_LO:
     442      case GLP_DB:
     443      case GLP_FX:
     444        if (lo == up)
     445          glp_set_row_bnds(lp, i, GLP_FX, lo, up);
     446        else
     447          glp_set_row_bnds(lp, i, GLP_DB, lo, up);
     448        break;
     449      default:
     450        break;
     451      }
     452    }
     453  }
     454
     455  SolverBase::Value GlpkSolverBase::_getRowUpperBound(int i) const {
     456    int b = glp_get_row_type(lp, i);
     457    switch (b) {
     458    case GLP_UP:
     459    case GLP_DB:
     460    case GLP_FX:
     461      return glp_get_row_ub(lp, i);
     462    default:
     463      return INF;
     464    }
     465  }
     466
     467  void GlpkSolverBase::_setObjCoeffs(ExprIterator b, ExprIterator e) {
     468    for (int i = 1; i <= glp_get_num_cols(lp); ++i) {
     469      glp_set_obj_coef(lp, i, 0.0);
     470    }
     471    for (ExprIterator it = b; it != e; ++it) {
     472      glp_set_obj_coef(lp, it->first, it->second);
     473    }
     474  }
     475
     476  void GlpkSolverBase::_getObjCoeffs(InsertIterator b) const {
     477    for (int i = 1; i <= glp_get_num_cols(lp); ++i) {
     478      Value val = glp_get_obj_coef(lp, i);
     479      if (val != 0.0) {
     480        *b = std::make_pair(i, val);
     481        ++b;
     482      }
     483    }
     484  }
     485
     486  void GlpkSolverBase::_setObjCoeff(int i, Value obj_coef) {
     487    //i = 0 means the constant term (shift)
     488    glp_set_obj_coef(lp, i, obj_coef);
     489  }
     490
     491  SolverBase::Value GlpkSolverBase::_getObjCoeff(int i) const {
     492    //i = 0 means the constant term (shift)
     493    return glp_get_obj_coef(lp, i);
     494  }
     495
     496  void GlpkSolverBase::_setSense(SolverBase::Sense sense) {
     497    switch (sense) {
     498    case MIN:
     499      glp_set_obj_dir(lp, GLP_MIN);
     500      break;
     501    case MAX:
     502      glp_set_obj_dir(lp, GLP_MAX);
     503      break;
     504    }
     505  }
     506
     507  SolverBase::Sense GlpkSolverBase::_getSense() const {
     508    switch(glp_get_obj_dir(lp)) {
     509    case GLP_MIN:
     510      return MIN;
     511    case GLP_MAX:
     512      return MAX;
     513    default:
     514      LEMON_ASSERT(false, "Wrong sense");
     515      return SolverBase::Sense();
     516    }
     517  }
     518
     519  void GlpkSolverBase::_clear() {
     520    glp_erase_prob(lp);
     521    rows.clear();
     522    cols.clear();
     523  }
     524
     525  // LpGlpk members
     526
     527  LpGlpk::LpGlpk()
     528    : SolverBase(), GlpkSolverBase(), LpSolverBase() {
     529    messageLevel(MESSAGE_NO_OUTPUT);
     530  }
     531
     532  LpGlpk::LpGlpk(const LpGlpk& other)
     533    : SolverBase(other), GlpkSolverBase(other), LpSolverBase(other) {
     534    messageLevel(MESSAGE_NO_OUTPUT);
     535  }
     536
     537  LpGlpk* LpGlpk::_newSolver() const { return new LpGlpk; }
     538  LpGlpk* LpGlpk::_cloneSolver() const {return new LpGlpk(*this); }
     539
     540  void LpGlpk::_clear_temporals() {
     541    _primal_ray.clear();
     542    _dual_ray.clear();
     543  }
     544
     545  LpGlpk::SolveExitStatus LpGlpk::_solve() {
     546    return solvePrimal();
     547  }
     548
     549  LpGlpk::SolveExitStatus LpGlpk::solvePrimal() {
     550    _clear_temporals();
     551
     552    glp_smcp smcp;
     553    glp_init_smcp(&smcp);
     554
     555    switch (_message_level) {
     556    case MESSAGE_NO_OUTPUT:
     557      smcp.msg_lev = GLP_MSG_OFF;
     558      break;
     559    case MESSAGE_ERROR_MESSAGE:
     560      smcp.msg_lev = GLP_MSG_ERR;
     561      break;
     562    case MESSAGE_NORMAL_OUTPUT:
     563      smcp.msg_lev = GLP_MSG_ON;     
     564      break;
     565    case MESSAGE_FULL_OUTPUT:
     566      smcp.msg_lev = GLP_MSG_ALL;
     567      break;
     568    }
     569
     570    if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
     571    return SOLVED;
     572  }
     573
     574  LpGlpk::SolveExitStatus LpGlpk::solveDual() {
     575    _clear_temporals();
     576
     577    glp_smcp smcp;
     578    glp_init_smcp(&smcp);
     579
     580    switch (_message_level) {
     581    case MESSAGE_NO_OUTPUT:
     582      smcp.msg_lev = GLP_MSG_OFF;
     583      break;
     584    case MESSAGE_ERROR_MESSAGE:
     585      smcp.msg_lev = GLP_MSG_ERR;
     586      break;
     587    case MESSAGE_NORMAL_OUTPUT:
     588      smcp.msg_lev = GLP_MSG_ON;     
     589      break;
     590    case MESSAGE_FULL_OUTPUT:
     591      smcp.msg_lev = GLP_MSG_ALL;
     592      break;
     593    }
     594    smcp.meth = GLP_DUAL;
     595   
     596    if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
     597    return SOLVED;
     598  }
     599
     600  LpGlpk::Value LpGlpk::_getPrimal(int i) const {
     601    return glp_get_col_prim(lp, i);
     602  }
     603
     604  LpGlpk::Value LpGlpk::_getDual(int i) const {
     605    return glp_get_row_dual(lp, i);
     606  }
     607
     608  LpGlpk::Value LpGlpk::_getPrimalValue() const {
     609    return glp_get_obj_val(lp);
     610  }
     611
     612  LpGlpk::VarStatus LpGlpk::_getColStatus(int i) const {
     613    switch (glp_get_col_stat(lp, i)) {
     614    case GLP_BS:
     615      return BASIC;
     616    case GLP_UP:
     617      return UPPER;
     618    case GLP_LO:
     619      return LOWER;
     620    case GLP_NF:
     621      return FREE;
     622    case GLP_NS:
     623      return FIXED;
     624    default:
     625      LEMON_ASSERT(false, "Wrong column status");
     626      return LpGlpk::VarStatus();
     627    }
     628  }
     629
     630  LpGlpk::VarStatus LpGlpk::_getRowStatus(int i) const {
     631    switch (glp_get_row_stat(lp, i)) {
     632    case GLP_BS:
     633      return BASIC;
     634    case GLP_UP:
     635      return UPPER;
     636    case GLP_LO:
     637      return LOWER;
     638    case GLP_NF:
     639      return FREE;
     640    case GLP_NS:
     641      return FIXED;
     642    default:
     643      LEMON_ASSERT(false, "Wrong row status");
     644      return LpGlpk::VarStatus();
     645    }
     646  }
     647
     648  LpGlpk::Value LpGlpk::_getPrimalRay(int i) const {
     649    if (_primal_ray.empty()) {
     650      int row_num = glp_get_num_rows(lp);
     651      int col_num = glp_get_num_cols(lp);
     652
     653      _primal_ray.resize(col_num + 1, 0.0);
     654     
     655      int index = glp_get_unbnd_ray(lp);
     656      if (index != 0) {
     657        // The primal ray is found in primal simplex second phase       
     658        LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :
     659                      glp_get_col_stat(lp, index - row_num)) != GLP_BS,
     660                     "Wrong primal ray");
     661
     662        bool negate = glp_get_obj_dir(lp) == GLP_MAX;
     663       
     664        if (index > row_num) {
     665          _primal_ray[index - row_num] = 1.0;
     666          if (glp_get_col_dual(lp, index - row_num) > 0) {
     667            negate = !negate;
     668          }
     669        } else {
     670          if (glp_get_row_dual(lp, index) > 0) {
     671            negate = !negate;
     672          }
     673        }
     674       
     675        std::vector<int> ray_indexes(row_num + 1);
     676        std::vector<Value> ray_values(row_num + 1);
     677        int ray_length = glp_eval_tab_col(lp, index, &ray_indexes.front(),
     678                                          &ray_values.front());
     679       
     680        for (int i = 1; i <= ray_length; ++i) {
     681          if (ray_indexes[i] > row_num) {
     682            _primal_ray[ray_indexes[i] - row_num] = ray_values[i];
     683          }
     684        }
     685       
     686        if (negate) {
     687          for (int i = 1; i <= col_num; ++i) {
     688            _primal_ray[i] = - _primal_ray[i];
     689          }
     690        }
     691      } else {
     692        for (int i = 1; i <= col_num; ++i) {
     693          _primal_ray[i] = glp_get_col_prim(lp, i);
     694        }
     695      }
     696    }
     697    return _primal_ray[i];
     698  }
     699
     700  LpGlpk::Value LpGlpk::_getDualRay(int i) const {
     701    if (_dual_ray.empty()) {
     702      int row_num = glp_get_num_rows(lp);
     703
     704      _dual_ray.resize(row_num + 1, 0.0);
     705     
     706      int index = glp_get_unbnd_ray(lp);
     707      if (index != 0) {
     708        // The dual ray is found in dual simplex second phase
     709        LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :
     710                      glp_get_col_stat(lp, index - row_num)) == GLP_BS,
     711
     712                     "Wrong dual ray");
     713
     714        int idx;
     715        bool negate = false;
     716       
     717        if (index > row_num) {
     718          idx = glp_get_col_bind(lp, index - row_num);
     719          if (glp_get_col_prim(lp, index - row_num) >
     720              glp_get_col_ub(lp, index - row_num)) {
     721            negate = true;
     722          }
     723        } else {
     724          idx = glp_get_row_bind(lp, index);
     725          if (glp_get_row_prim(lp, index) > glp_get_row_ub(lp, index)) {
     726            negate = true;
     727          }
     728        }
     729
     730        _dual_ray[idx] = negate ?  - 1.0 : 1.0;
     731
     732        glp_btran(lp, &_dual_ray.front());         
     733      } else {
     734        double eps = 1e-7;
     735        // The dual ray is found in primal simplex first phase
     736        // We assume that the glpk minimizes the slack to get feasible solution
     737        for (int i = 1; i <= row_num; ++i) {
     738          int index = glp_get_bhead(lp, i);
     739          if (index <= row_num) {
     740            double res = glp_get_row_prim(lp, index);
     741            if (res > glp_get_row_ub(lp, index) + eps) {
     742              _dual_ray[i] = -1;
     743            } else if (res < glp_get_row_lb(lp, index) - eps) {
     744              _dual_ray[i] = 1;
     745            } else {
     746              _dual_ray[i] = 0;
     747            }
     748            _dual_ray[i] *= glp_get_rii(lp, index);
     749          } else {
     750            double res = glp_get_col_prim(lp, index - row_num);
     751            if (res > glp_get_col_ub(lp, index - row_num) + eps) {
     752              _dual_ray[i] = -1;
     753            } else if (res < glp_get_col_lb(lp, index - row_num) - eps) {
     754              _dual_ray[i] = 1;
     755            } else {
     756              _dual_ray[i] = 0;
     757            }
     758            _dual_ray[i] /= glp_get_sjj(lp, index - row_num);
     759          }         
     760        }
     761
     762        glp_btran(lp, &_dual_ray.front());
     763
     764        for (int i = 1; i <= row_num; ++i) {
     765          _dual_ray[i] /= glp_get_rii(lp, i);
     766        }
     767      }
     768    }
     769    return _dual_ray[i];
     770  }
     771
     772  LpGlpk::ProblemType LpGlpk::_getPrimalType() const {
     773    if (glp_get_status(lp) == GLP_OPT)
     774      return OPTIMAL;
     775    switch (glp_get_prim_stat(lp)) {
     776    case GLP_UNDEF:
     777      return UNDEFINED;
     778    case GLP_FEAS:
     779    case GLP_INFEAS:
     780      if (glp_get_dual_stat(lp) == GLP_NOFEAS) {
     781        return UNBOUNDED;
     782      } else {
     783        return UNDEFINED;
     784      }
     785    case GLP_NOFEAS:
     786      return INFEASIBLE;
     787    default:
     788      LEMON_ASSERT(false, "Wrong primal type");
     789      return  LpGlpk::ProblemType();
     790    }
     791  }
     792
     793  LpGlpk::ProblemType LpGlpk::_getDualType() const {
     794    if (glp_get_status(lp) == GLP_OPT)
     795      return OPTIMAL;
     796    switch (glp_get_dual_stat(lp)) {
     797    case GLP_UNDEF:
     798      return UNDEFINED;
     799    case GLP_FEAS:
     800    case GLP_INFEAS:
     801      if (glp_get_prim_stat(lp) == GLP_NOFEAS) {
     802        return UNBOUNDED;
     803      } else {
     804        return UNDEFINED;
     805      }
     806    case GLP_NOFEAS:
     807      return INFEASIBLE;
     808    default:
     809      LEMON_ASSERT(false, "Wrong primal type");
     810      return  LpGlpk::ProblemType();
     811    }
     812  }
     813
     814  void LpGlpk::presolver(bool b) {
     815    lpx_set_int_parm(lp, LPX_K_PRESOL, b ? 1 : 0);
     816  }
     817
     818  void LpGlpk::messageLevel(MessageLevel m) {
     819    _message_level = m;
     820  }
     821
     822  // MipGlpk members
     823
     824  MipGlpk::MipGlpk()
     825    : SolverBase(), GlpkSolverBase(), MipSolverBase() {
     826    messageLevel(MESSAGE_NO_OUTPUT);
     827  }
     828
     829  MipGlpk::MipGlpk(const MipGlpk& other)
     830    : SolverBase(), GlpkSolverBase(other), MipSolverBase() {
     831    messageLevel(MESSAGE_NO_OUTPUT);
     832  }
     833
     834  void MipGlpk::_setColType(int i, MipGlpk::ColTypes col_type) {
     835    switch (col_type) {
     836    case INTEGER:
     837      glp_set_col_kind(lp, i, GLP_IV);
     838      break;
     839    case REAL:
     840      glp_set_col_kind(lp, i, GLP_CV);
     841      break;
     842    }
     843  }
     844
     845  MipGlpk::ColTypes MipGlpk::_getColType(int i) const {
     846    switch (glp_get_col_kind(lp, i)) {
     847    case GLP_IV:
     848    case GLP_BV:
     849      return INTEGER;
     850    default:
     851      return REAL;
     852    }
     853
     854  }
     855
     856  MipGlpk::SolveExitStatus MipGlpk::_solve() {
     857    glp_smcp smcp;
     858    glp_init_smcp(&smcp);
     859
     860    switch (_message_level) {
     861    case MESSAGE_NO_OUTPUT:
     862      smcp.msg_lev = GLP_MSG_OFF;
     863      break;
     864    case MESSAGE_ERROR_MESSAGE:
     865      smcp.msg_lev = GLP_MSG_ERR;
     866      break;
     867    case MESSAGE_NORMAL_OUTPUT:
     868      smcp.msg_lev = GLP_MSG_ON;     
     869      break;
     870    case MESSAGE_FULL_OUTPUT:
     871      smcp.msg_lev = GLP_MSG_ALL;
     872      break;
     873    }
     874    smcp.meth = GLP_DUAL;
     875
     876    if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
     877    if (glp_get_status(lp) != GLP_OPT) return SOLVED;
     878
     879    glp_iocp iocp;
     880    glp_init_iocp(&iocp);
     881
     882    switch (_message_level) {
     883    case MESSAGE_NO_OUTPUT:
     884      iocp.msg_lev = GLP_MSG_OFF;
     885      break;
     886    case MESSAGE_ERROR_MESSAGE:
     887      iocp.msg_lev = GLP_MSG_ERR;
     888      break;
     889    case MESSAGE_NORMAL_OUTPUT:
     890      iocp.msg_lev = GLP_MSG_ON;     
     891      break;
     892    case MESSAGE_FULL_OUTPUT:
     893      iocp.msg_lev = GLP_MSG_ALL;
     894      break;
     895    }
     896
     897    if (glp_intopt(lp, &iocp) != 0) return UNSOLVED;
     898    return SOLVED;
     899  }
     900
     901
     902  MipGlpk::ProblemType MipGlpk::_getType() const {
     903    switch (glp_get_status(lp)) {
     904    case GLP_OPT:
     905      switch (glp_mip_status(lp)) {
     906      case GLP_UNDEF:
     907        return UNDEFINED;
     908      case GLP_NOFEAS:
     909        return INFEASIBLE;
     910      case GLP_FEAS:
     911        return FEASIBLE;
     912      case GLP_OPT:
     913        return OPTIMAL;
     914      default:
     915        LEMON_ASSERT(false, "Wrong problem type.");
     916        return MipGlpk::ProblemType();
     917      }
     918    case GLP_NOFEAS:
     919      return INFEASIBLE;
     920    case GLP_INFEAS:
     921    case GLP_FEAS:
     922      if (glp_get_dual_stat(lp) == GLP_NOFEAS) {
     923        return UNBOUNDED;
     924      } else {
     925        return UNDEFINED;
     926      }
     927    default:
     928      LEMON_ASSERT(false, "Wrong problem type.");
     929      return MipGlpk::ProblemType();
     930    }
     931  }
     932
     933  MipGlpk::Value MipGlpk::_getSol(int i) const {
     934    return glp_mip_col_val(lp, i);
     935  }
     936
     937  MipGlpk::Value MipGlpk::_getSolValue() const {
     938    return glp_mip_obj_val(lp);
     939  }
     940
     941  MipGlpk* MipGlpk::_newSolver() const { return new MipGlpk; }
     942  MipGlpk* MipGlpk::_cloneSolver() const {return new MipGlpk(*this); }
     943
     944  void MipGlpk::messageLevel(MessageLevel m) {
     945    _message_level = m;
     946  }
     947
     948} //END OF NAMESPACE LEMON
  • new file lemon/lp_glpk.h

    diff -r cbe3ec2d59d2 -r aa6dff56b2e1 lemon/lp_glpk.h
    - +  
     1/* -*- C++ -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library
     4 *
     5 * Copyright (C) 2003-2008
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#ifndef LEMON_LP_GLPK_H
     20#define LEMON_LP_GLPK_H
     21
     22///\file
     23///\brief Header of the LEMON-GLPK lp solver interface.
     24///\ingroup lp_group
     25
     26#include <lemon/lp_base.h>
     27
     28// forward declaration
     29#ifndef _GLP_PROB
     30#define _GLP_PROB
     31typedef struct { double _prob; } glp_prob;
     32/* LP/MIP problem object */
     33#endif
     34
     35namespace lemon {
     36
     37 
     38  /// \brief Base interface for the GLPK LP and MIP solver
     39  ///
     40  /// This class implements the common interface of the GLPK LP and MIP solver.
     41  /// \ingroup lp_group
     42  class GlpkSolverBase : virtual public SolverBase {
     43  protected:
     44   
     45    typedef glp_prob LPX;
     46    glp_prob* lp;
     47   
     48    typedef LpSolverBase Parent;
     49   
     50    GlpkSolverBase();
     51    GlpkSolverBase(const GlpkSolverBase&);
     52    virtual ~GlpkSolverBase();
     53
     54  protected:
     55
     56    virtual int _addCol();
     57    virtual int _addRow();
     58
     59    virtual void _eraseCol(int i);
     60    virtual void _eraseRow(int i);
     61
     62    virtual void _eraseColId(int i);
     63    virtual void _eraseRowId(int i);
     64
     65    virtual void _getColName(int col, std::string& name) const;
     66    virtual void _setColName(int col, const std::string& name);
     67    virtual int _colByName(const std::string& name) const;
     68
     69    virtual void _getRowName(int row, std::string& name) const;
     70    virtual void _setRowName(int row, const std::string& name);
     71    virtual int _rowByName(const std::string& name) const;
     72
     73    virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e);
     74    virtual void _getRowCoeffs(int i, InsertIterator b) const;
     75
     76    virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e);
     77    virtual void _getColCoeffs(int i, InsertIterator b) const;
     78
     79    virtual void _setCoeff(int row, int col, Value value);
     80    virtual Value _getCoeff(int row, int col) const;
     81
     82    virtual void _setColLowerBound(int i, Value value);
     83    virtual Value _getColLowerBound(int i) const;
     84
     85    virtual void _setColUpperBound(int i, Value value);
     86    virtual Value _getColUpperBound(int i) const;
     87
     88    virtual void _setRowLowerBound(int i, Value value);
     89    virtual Value _getRowLowerBound(int i) const;
     90
     91    virtual void _setRowUpperBound(int i, Value value);
     92    virtual Value _getRowUpperBound(int i) const;
     93
     94    virtual void _setObjCoeffs(ExprIterator b, ExprIterator e);
     95    virtual void _getObjCoeffs(InsertIterator b) const;
     96   
     97    virtual void _setObjCoeff(int i, Value obj_coef);
     98    virtual Value _getObjCoeff(int i) const;
     99   
     100    virtual void _setSense(Sense);
     101    virtual Sense _getSense() const;
     102
     103    virtual void _clear();
     104
     105  public:
     106
     107    ///Pointer to the underlying GLPK data structure.
     108    LPX *lpx() {return lp;}
     109    ///Const pointer to the underlying GLPK data structure.
     110    const LPX *lpx() const {return lp;}
     111
     112    ///Returns the constraint identifier understood by GLPK.
     113    int lpxRow(Row r) const { return rows(id(r)); }
     114
     115    ///Returns the variable identifier understood by GLPK.
     116    int lpxCol(Col c) const { return cols(id(c)); }
     117
     118  };
     119 
     120  /// \brief Interface for the GLPK LP solver
     121  ///
     122  /// This class implements an interface for the GLPK LP solver.
     123  ///\ingroup lp_group
     124  class LpGlpk : public GlpkSolverBase, public LpSolverBase {
     125  public:
     126   
     127    ///\e
     128    LpGlpk();
     129    ///\e
     130    LpGlpk(const LpGlpk&);
     131
     132  private:
     133   
     134    mutable std::vector<double> _primal_ray;
     135    mutable std::vector<double> _dual_ray;
     136
     137    void _clear_temporals();
     138   
     139  protected:
     140
     141    virtual LpGlpk* _cloneSolver() const;
     142    virtual LpGlpk* _newSolver() const;
     143
     144    ///\todo It should be clarified
     145    ///
     146    virtual SolveExitStatus _solve();
     147    virtual Value _getPrimal(int i) const;
     148    virtual Value _getDual(int i) const;
     149
     150    virtual Value _getPrimalValue() const;
     151
     152    virtual VarStatus _getColStatus(int i) const;
     153    virtual VarStatus _getRowStatus(int i) const;
     154
     155    virtual Value _getPrimalRay(int i) const;
     156    virtual Value _getDualRay(int i) const;
     157   
     158    ///\todo It should be clarified
     159    ///
     160    virtual ProblemType _getPrimalType() const;
     161    virtual ProblemType _getDualType() const;
     162
     163  public:
     164
     165    ///Solve with primal simplex
     166    SolveExitStatus solvePrimal();
     167
     168    ///Solve with dual simplex
     169    SolveExitStatus solveDual();
     170
     171    ///Turns on or off the presolver
     172
     173    ///Turns on (\c b is \c true) or off (\c b is \c false) the presolver
     174    ///
     175    ///The presolver is off by default.
     176    void presolver(bool b);
     177
     178    ///Enum for \c messageLevel() parameter
     179    enum MessageLevel {
     180      /// no output (default value)
     181      MESSAGE_NO_OUTPUT = 0,
     182      /// error messages only
     183      MESSAGE_ERROR_MESSAGE = 1,
     184      /// normal output
     185      MESSAGE_NORMAL_OUTPUT = 2,
     186      /// full output (includes informational messages)
     187      MESSAGE_FULL_OUTPUT = 3
     188    };
     189
     190  private:
     191   
     192    MessageLevel _message_level;
     193
     194  public:
     195
     196    ///Set the verbosity of the messages
     197
     198    ///Set the verbosity of the messages
     199    ///
     200    ///\param m is the level of the messages output by the solver routines.
     201    void messageLevel(MessageLevel m);
     202  };
     203
     204  /// \brief Interface for the GLPK MIP solver
     205  ///
     206  /// This class implements an interface for the GLPK MIP solver.
     207  ///\ingroup lp_group
     208  class MipGlpk : public GlpkSolverBase, public MipSolverBase { 
     209  public:
     210 
     211    ///\e
     212    MipGlpk();
     213    ///\e
     214    MipGlpk(const MipGlpk&);
     215   
     216  protected:
     217
     218    virtual MipGlpk* _cloneSolver() const;
     219    virtual MipGlpk* _newSolver() const;
     220 
     221    virtual ColTypes _getColType(int col) const;
     222    virtual void _setColType(int col, ColTypes col_type);
     223   
     224    virtual SolveExitStatus _solve();
     225    virtual ProblemType _getType() const;
     226    virtual Value _getSol(int i) const;
     227    virtual Value _getSolValue() const;
     228
     229    ///Enum for \c messageLevel() parameter
     230    enum MessageLevel {
     231      /// no output (default value)
     232      MESSAGE_NO_OUTPUT = 0,
     233      /// error messages only
     234      MESSAGE_ERROR_MESSAGE = 1,
     235      /// normal output
     236      MESSAGE_NORMAL_OUTPUT = 2,
     237      /// full output (includes informational messages)
     238      MESSAGE_FULL_OUTPUT = 3
     239    };
     240
     241  private:
     242   
     243    MessageLevel _message_level;
     244
     245  public:
     246
     247    ///Set the verbosity of the messages
     248
     249    ///Set the verbosity of the messages
     250    ///
     251    ///\param m is the level of the messages output by the solver routines.
     252    void messageLevel(MessageLevel m);
     253  };
     254
     255
     256} //END OF NAMESPACE LEMON
     257
     258#endif //LEMON_LP_GLPK_H
     259
  • new file lemon/lp_skeleton.cc

    diff -r cbe3ec2d59d2 -r aa6dff56b2e1 lemon/lp_skeleton.cc
    - +  
     1/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library.
     4 *
     5 * Copyright (C) 2003-2008
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#include <lemon/lp_skeleton.h>
     20
     21///\file
     22///\brief A skeleton file to implement LP solver interfaces
     23namespace lemon {
     24
     25  int SkeletonSolverBase::_addCol()
     26  {
     27    return ++col_num;
     28  }
     29
     30  int SkeletonSolverBase::_addRow()
     31  {
     32    return ++row_num;
     33  }
     34
     35  void SkeletonSolverBase::_eraseCol(int) {}
     36  void SkeletonSolverBase::_eraseRow(int) {}
     37
     38  void SkeletonSolverBase::_getColName(int, std::string &) const {}
     39  void SkeletonSolverBase::_setColName(int, const std::string &) {}
     40  int SkeletonSolverBase::_colByName(const std::string&) const { return -1; }
     41
     42  void SkeletonSolverBase::_getRowName(int, std::string &) const {}
     43  void SkeletonSolverBase::_setRowName(int, const std::string &) {}
     44  int SkeletonSolverBase::_rowByName(const std::string&) const { return -1; }
     45
     46  void SkeletonSolverBase::_setRowCoeffs(int, ExprIterator, ExprIterator) {}
     47  void SkeletonSolverBase::_getRowCoeffs(int, InsertIterator) const {}
     48
     49  void SkeletonSolverBase::_setColCoeffs(int, ExprIterator, ExprIterator) {}
     50  void SkeletonSolverBase::_getColCoeffs(int, InsertIterator) const {}
     51
     52  void SkeletonSolverBase::_setCoeff(int, int, Value) {}
     53  SkeletonSolverBase::Value SkeletonSolverBase::_getCoeff(int, int) const
     54  { return 0; }
     55
     56  void SkeletonSolverBase::_setColLowerBound(int, Value) {}
     57  SkeletonSolverBase::Value SkeletonSolverBase::_getColLowerBound(int) const
     58  {  return 0; }
     59
     60  void SkeletonSolverBase::_setColUpperBound(int, Value) {}
     61  SkeletonSolverBase::Value SkeletonSolverBase::_getColUpperBound(int) const
     62  {  return 0; }
     63
     64  void SkeletonSolverBase::_setRowLowerBound(int, Value) {}
     65  SkeletonSolverBase::Value SkeletonSolverBase::_getRowLowerBound(int) const
     66  {  return 0; }
     67
     68  void SkeletonSolverBase::_setRowUpperBound(int, Value) {}
     69  SkeletonSolverBase::Value SkeletonSolverBase::_getRowUpperBound(int) const
     70  {  return 0; }
     71
     72  void SkeletonSolverBase::_setObjCoeffs(ExprIterator, ExprIterator) {}
     73  void SkeletonSolverBase::_getObjCoeffs(InsertIterator) const {};
     74
     75  void SkeletonSolverBase::_setObjCoeff(int, Value) {}
     76  SkeletonSolverBase::Value SkeletonSolverBase::_getObjCoeff(int) const
     77  {  return 0; }
     78
     79  void SkeletonSolverBase::_setSense(Sense) {}
     80  SkeletonSolverBase::Sense SkeletonSolverBase::_getSense() const
     81  { return MIN; }
     82
     83  void SkeletonSolverBase::_clear() {
     84    row_num = col_num = 0;
     85  }
     86
     87  LpSkeleton::SolveExitStatus LpSkeleton::_solve() { return SOLVED; }
     88
     89  LpSkeleton::Value LpSkeleton::_getPrimal(int) const { return 0; }
     90  LpSkeleton::Value LpSkeleton::_getDual(int) const { return 0; }
     91  LpSkeleton::Value LpSkeleton::_getPrimalValue() const { return 0; }
     92
     93  LpSkeleton::Value LpSkeleton::_getPrimalRay(int) const { return 0; }
     94  LpSkeleton::Value LpSkeleton::_getDualRay(int) const { return 0; }
     95
     96  LpSkeleton::ProblemType LpSkeleton::_getPrimalType() const
     97  { return UNDEFINED; }
     98
     99  LpSkeleton::ProblemType LpSkeleton::_getDualType() const
     100  { return UNDEFINED; }
     101
     102  LpSkeleton::VarStatus LpSkeleton::_getColStatus(int) const
     103  { return BASIC; }
     104
     105  LpSkeleton::VarStatus LpSkeleton::_getRowStatus(int) const
     106  { return BASIC; }
     107
     108  LpSkeleton* LpSkeleton::_newSolver() const
     109  { return static_cast<LpSkeleton*>(0); }
     110
     111  LpSkeleton* LpSkeleton::_cloneSolver() const
     112  { return static_cast<LpSkeleton*>(0); }
     113
     114  MipSkeleton::SolveExitStatus MipSkeleton::_solve()
     115  { return SOLVED; }
     116
     117  MipSkeleton::Value MipSkeleton::_getSol(int) const { return 0; }
     118  MipSkeleton::Value MipSkeleton::_getSolValue() const { return 0; }
     119
     120  MipSkeleton::ProblemType MipSkeleton::_getType() const
     121  { return UNDEFINED; }
     122
     123  MipSkeleton* MipSkeleton::_newSolver() const
     124  { return static_cast<MipSkeleton*>(0); }
     125
     126  MipSkeleton* MipSkeleton::_cloneSolver() const
     127  { return static_cast<MipSkeleton*>(0); }
     128
     129} //namespace lemon
     130
  • new file lemon/lp_skeleton.h

    diff -r cbe3ec2d59d2 -r aa6dff56b2e1 lemon/lp_skeleton.h
    - +  
     1/* -*- C++ -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library
     4 *
     5 * Copyright (C) 2003-2008
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#ifndef LEMON_LP_SKELETON
     20#define LEMON_LP_SKELETON
     21
     22#include <lemon/lp_base.h>
     23
     24///\file
     25///\brief A skeleton file to implement LP solver interfaces
     26namespace lemon {
     27
     28  ///A skeleton class to implement LP solver interfaces
     29  class SkeletonSolverBase : public virtual SolverBase {
     30    int col_num,row_num;
     31   
     32  protected:
     33
     34    SkeletonSolverBase()
     35      : col_num(-1), row_num(-1) {}
     36
     37    /// \e
     38    virtual int _addCol();
     39    /// \e
     40    virtual int _addRow();
     41    /// \e
     42    virtual void _eraseCol(int i);
     43    /// \e
     44    virtual void _eraseRow(int i);
     45
     46    /// \e
     47    virtual void _getColName(int col, std::string& name) const;
     48    /// \e
     49    virtual void _setColName(int col, const std::string& name);
     50    /// \e
     51    virtual int _colByName(const std::string& name) const;
     52
     53    /// \e
     54    virtual void _getRowName(int row, std::string& name) const;
     55    /// \e
     56    virtual void _setRowName(int row, const std::string& name);
     57    /// \e
     58    virtual int _rowByName(const std::string& name) const;
     59
     60    /// \e
     61    virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e);
     62    /// \e
     63    virtual void _getRowCoeffs(int i, InsertIterator b) const;
     64    /// \e
     65    virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e);
     66    /// \e
     67    virtual void _getColCoeffs(int i, InsertIterator b) const;
     68   
     69    /// Set one element of the coefficient matrix
     70    virtual void _setCoeff(int row, int col, Value value);
     71
     72    /// Get one element of the coefficient matrix
     73    virtual Value _getCoeff(int row, int col) const;
     74
     75    /// The lower bound of a variable (column) have to be given by an
     76    /// extended number of type Value, i.e. a finite number of type
     77    /// Value or -\ref INF.
     78    virtual void _setColLowerBound(int i, Value value);
     79    /// \e
     80
     81    /// The lower bound of a variable (column) is an
     82    /// extended number of type Value, i.e. a finite number of type
     83    /// Value or -\ref INF.
     84    virtual Value _getColLowerBound(int i) const;
     85
     86    /// The upper bound of a variable (column) have to be given by an
     87    /// extended number of type Value, i.e. a finite number of type
     88    /// Value or \ref INF.
     89    virtual void _setColUpperBound(int i, Value value);
     90    /// \e
     91
     92    /// The upper bound of a variable (column) is an
     93    /// extended number of type Value, i.e. a finite number of type
     94    /// Value or \ref INF.
     95    virtual Value _getColUpperBound(int i) const;
     96
     97    /// The lower bound of a constraint (row) have to be given by an
     98    /// extended number of type Value, i.e. a finite number of type
     99    /// Value or -\ref INF.
     100    virtual void _setRowLowerBound(int i, Value value);
     101    /// \e
     102
     103    /// The lower bound of a constraint (row) is an
     104    /// extended number of type Value, i.e. a finite number of type
     105    /// Value or -\ref INF.
     106    virtual Value _getRowLowerBound(int i) const;
     107
     108    /// The upper bound of a constraint (row) have to be given by an
     109    /// extended number of type Value, i.e. a finite number of type
     110    /// Value or \ref INF.
     111    virtual void _setRowUpperBound(int i, Value value);
     112    /// \e
     113
     114    /// The upper bound of a constraint (row) is an
     115    /// extended number of type Value, i.e. a finite number of type
     116    /// Value or \ref INF.
     117    virtual Value _getRowUpperBound(int i) const;
     118
     119    /// \e
     120    virtual void _setObjCoeffs(ExprIterator b, ExprIterator e);
     121    /// \e
     122    virtual void _getObjCoeffs(InsertIterator b) const;
     123
     124    /// \e
     125    virtual void _setObjCoeff(int i, Value obj_coef);
     126    /// \e
     127    virtual Value _getObjCoeff(int i) const;
     128
     129    ///\e
     130    virtual void _setSense(Sense);
     131    ///\e
     132    virtual Sense _getSense() const;
     133
     134    ///\e
     135    virtual void _clear();
     136
     137  };
     138
     139  /// \brief Interface for a skeleton LP solver
     140  ///
     141  /// This class implements an interface for a skeleton LP solver.
     142  ///\ingroup lp_group
     143  class LpSkeleton : public SkeletonSolverBase, public LpSolverBase {
     144  public:
     145    LpSkeleton() : SkeletonSolverBase(), LpSolverBase() {}
     146
     147  protected:
     148
     149    ///\e
     150    virtual SolveExitStatus _solve();
     151
     152    ///\e
     153    virtual Value _getPrimal(int i) const;
     154    ///\e
     155    virtual Value _getDual(int i) const;
     156
     157    ///\e
     158    virtual Value _getPrimalValue() const;
     159
     160    ///\e
     161    virtual Value _getPrimalRay(int i) const;
     162    ///\e
     163    virtual Value _getDualRay(int i) const;
     164
     165    ///\e
     166    virtual ProblemType _getPrimalType() const;
     167    ///\e
     168    virtual ProblemType _getDualType() const;
     169
     170    ///\e
     171    virtual VarStatus _getColStatus(int i) const;
     172    ///\e
     173    virtual VarStatus _getRowStatus(int i) const;
     174
     175    ///\e
     176    virtual LpSkeleton* _newSolver() const;
     177    ///\e
     178    virtual LpSkeleton* _cloneSolver() const;
     179
     180  };
     181
     182  /// \brief Interface for a skeleton MIP solver
     183  ///
     184  /// This class implements an interface for a skeleton MIP solver.
     185  ///\ingroup lp_group
     186  class MipSkeleton : public SkeletonSolverBase, public MipSolverBase {
     187  public:
     188    MipSkeleton() : SkeletonSolverBase(), MipSolverBase() {}
     189
     190  protected:
     191    ///\e
     192   
     193    ///\bug Wrong interface
     194    ///
     195    virtual SolveExitStatus _solve();
     196
     197    ///\e
     198
     199    ///\bug Wrong interface
     200    ///
     201    virtual Value _getSol(int i) const;
     202
     203    ///\e
     204
     205    ///\bug Wrong interface
     206    ///
     207    virtual Value _getSolValue() const;
     208
     209    ///\e
     210
     211    ///\bug Wrong interface
     212    ///
     213    virtual ProblemType _getType() const;
     214
     215    ///\e
     216    virtual MipSkeleton* _newSolver() const;
     217
     218    ///\e
     219    virtual MipSkeleton* _cloneSolver() const;
     220
     221  };
     222
     223} //namespace lemon
     224
     225#endif // LEMON_LP_SKELETON
  • new file lemon/lp_soplex.cc

    diff -r cbe3ec2d59d2 -r aa6dff56b2e1 lemon/lp_soplex.cc
    - +  
     1/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library.
     4 *
     5 * Copyright (C) 2003-2008
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#include <iostream>
     20#include <lemon/lp_soplex.h>
     21
     22#include <soplex/soplex.h>
     23
     24
     25///\file
     26///\brief Implementation of the LEMON-SOPLEX lp solver interface.
     27namespace lemon {
     28
     29  LpSoplex::LpSoplex()
     30    : SolverBase(), LpSolverBase() {
     31    soplex = new soplex::SoPlex;
     32  }
     33
     34  LpSoplex::~LpSoplex() {
     35    delete soplex;
     36  }
     37
     38  LpSoplex::LpSoplex(const LpSoplex& lp)
     39    : SolverBase(), LpSolverBase() {
     40    rows = lp.rows;
     41    cols = lp.cols;
     42
     43    soplex = new soplex::SoPlex;
     44    (*static_cast<soplex::SPxLP*>(soplex)) = *(lp.soplex);
     45
     46    _col_names = lp._col_names;
     47    _col_names_ref = lp._col_names_ref;
     48
     49    _row_names = lp._row_names;
     50    _row_names_ref = lp._row_names_ref;
     51
     52  }
     53 
     54  void LpSoplex::_clear_temporals() {
     55    _primal_values.clear();
     56    _dual_values.clear();
     57  }
     58
     59  LpSoplex* LpSoplex::_newSolver() const {
     60    LpSoplex* newlp = new LpSoplex();
     61    return newlp;
     62  }
     63
     64  LpSoplex* LpSoplex::_cloneSolver() const {
     65    LpSoplex* newlp = new LpSoplex(*this);
     66    return newlp;
     67  }
     68
     69  int LpSoplex::_addCol() {
     70    soplex::LPCol c;
     71    c.setLower(-soplex::infinity);
     72    c.setUpper(soplex::infinity);
     73    soplex->addCol(c);
     74
     75    _col_names.push_back(std::string());
     76
     77    return soplex->nCols() - 1;
     78  }
     79
     80  int LpSoplex::_addRow() {
     81    soplex::LPRow r;
     82    r.setLhs(-soplex::infinity);
     83    r.setRhs(soplex::infinity);
     84    soplex->addRow(r);
     85
     86    _row_names.push_back(std::string());
     87
     88    return soplex->nRows() - 1;
     89  }
     90
     91
     92  void LpSoplex::_eraseCol(int i) {
     93    soplex->removeCol(i);
     94    _col_names_ref.erase(_col_names[i]);
     95    _col_names[i] = _col_names.back();
     96    _col_names_ref[_col_names.back()] = i;
     97    _col_names.pop_back();
     98  }
     99
     100  void LpSoplex::_eraseRow(int i) {
     101    soplex->removeRow(i);
     102    _row_names_ref.erase(_row_names[i]);
     103    _row_names[i] = _row_names.back();
     104    _row_names_ref[_row_names.back()] = i;
     105    _row_names.pop_back();
     106  }
     107
     108  void LpSoplex::_eraseColId(int i) {
     109    cols.eraseIndex(i);
     110    cols.relocateIndex(i, cols.maxIndex());
     111  }
     112  void LpSoplex::_eraseRowId(int i) {
     113    rows.eraseIndex(i);
     114    rows.relocateIndex(i, rows.maxIndex());
     115  }
     116
     117  void LpSoplex::_getColName(int c, std::string &name) const {
     118    name = _col_names[c];
     119  }
     120
     121  void LpSoplex::_setColName(int c, const std::string &name) {
     122    _col_names_ref.erase(_col_names[c]);
     123    _col_names[c] = name;
     124    if (!name.empty()) {
     125      _col_names_ref.insert(std::make_pair(name, c));
     126    }
     127  }
     128
     129  int LpSoplex::_colByName(const std::string& name) const {
     130    std::map<std::string, int>::const_iterator it =
     131      _col_names_ref.find(name);
     132    if (it != _col_names_ref.end()) {
     133      return it->second;
     134    } else {
     135      return -1;
     136    }
     137  }
     138
     139  void LpSoplex::_getRowName(int r, std::string &name) const {
     140    name = _row_names[r];
     141  }
     142
     143  void LpSoplex::_setRowName(int r, const std::string &name) {
     144    _row_names_ref.erase(_row_names[r]);
     145    _row_names[r] = name;
     146    if (!name.empty()) {
     147      _row_names_ref.insert(std::make_pair(name, r));
     148    }
     149  }
     150
     151  int LpSoplex::_rowByName(const std::string& name) const {
     152    std::map<std::string, int>::const_iterator it =
     153      _row_names_ref.find(name);
     154    if (it != _row_names_ref.end()) {
     155      return it->second;
     156    } else {
     157      return -1;
     158    }
     159  }
     160
     161
     162  void LpSoplex::_setRowCoeffs(int i, ExprIterator b, ExprIterator e) {
     163    for (int j = 0; j < soplex->nCols(); ++j) {
     164      soplex->changeElement(i, j, 0.0);
     165    }
     166    for(ExprIterator it = b; it != e; ++it) {
     167      soplex->changeElement(i, it->first, it->second);
     168    }
     169  }
     170
     171  void LpSoplex::_getRowCoeffs(int i, InsertIterator b) const {
     172    const soplex::SVector& vec = soplex->rowVector(i);
     173    for (int k = 0; k < vec.size(); ++k) {
     174      *b = std::make_pair(vec.index(k), vec.value(k));
     175      ++b;
     176    }
     177  }
     178
     179  void LpSoplex::_setColCoeffs(int j, ExprIterator b, ExprIterator e) {
     180    for (int i = 0; i < soplex->nRows(); ++i) {
     181      soplex->changeElement(i, j, 0.0);
     182    }
     183    for(ExprIterator it = b; it != e; ++it) {
     184      soplex->changeElement(it->first, j, it->second);
     185    }
     186  }
     187
     188  void LpSoplex::_getColCoeffs(int i, InsertIterator b) const {
     189    const soplex::SVector& vec = soplex->colVector(i);
     190    for (int k = 0; k < vec.size(); ++k) {
     191      *b = std::make_pair(vec.index(k), vec.value(k));
     192      ++b;
     193    }
     194  }
     195
     196  void LpSoplex::_setCoeff(int i, int j, Value value) {
     197    soplex->changeElement(i, j, value);
     198  }
     199
     200  LpSoplex::Value LpSoplex::_getCoeff(int i, int j) const {
     201    return soplex->rowVector(i)[j];
     202  }
     203
     204  void LpSoplex::_setColLowerBound(int i, Value value) {
     205    LEMON_ASSERT(value != INF, "Invalid bound");
     206    soplex->changeLower(i, value != -INF ? value : -soplex::infinity);
     207  }
     208
     209  LpSoplex::Value LpSoplex::_getColLowerBound(int i) const {
     210    double value = soplex->lower(i);
     211    return value != -soplex::infinity ? value : -INF;
     212  }
     213
     214  void LpSoplex::_setColUpperBound(int i, Value value) {
     215    LEMON_ASSERT(value != -INF, "Invalid bound");
     216    soplex->changeUpper(i, value != INF ? value : soplex::infinity);
     217  }
     218
     219  LpSoplex::Value LpSoplex::_getColUpperBound(int i) const {
     220    double value = soplex->upper(i);
     221    return value != soplex::infinity ? value : INF;
     222  }
     223
     224  void LpSoplex::_setRowLowerBound(int i, Value lb) {
     225    LEMON_ASSERT(lb != INF, "Invalid bound");
     226    soplex->changeRange(i, lb != -INF ? lb : -soplex::infinity, soplex->rhs(i));
     227  }
     228
     229  LpSoplex::Value LpSoplex::_getRowLowerBound(int i) const {
     230    double res = soplex->lhs(i);
     231    return res == -soplex::infinity ? -INF : res;
     232  }
     233
     234  void LpSoplex::_setRowUpperBound(int i, Value ub) {
     235    LEMON_ASSERT(ub != -INF, "Invalid bound");
     236    soplex->changeRange(i, soplex->lhs(i), ub != INF ? ub : soplex::infinity);
     237  }
     238
     239  LpSoplex::Value LpSoplex::_getRowUpperBound(int i) const {
     240    double res = soplex->rhs(i);
     241    return res == soplex::infinity ? INF : res;
     242  }
     243
     244  void LpSoplex::_setObjCoeffs(ExprIterator b, ExprIterator e) {
     245    for (int j = 0; j < soplex->nCols(); ++j) {
     246      soplex->changeObj(j, 0.0);
     247    }
     248    for (ExprIterator it = b; it != e; ++it) {
     249      soplex->changeObj(it->first, it->second);
     250    }
     251  }
     252
     253  void LpSoplex::_getObjCoeffs(InsertIterator b) const {
     254    for (int j = 0; j < soplex->nCols(); ++j) {
     255      Value coef = soplex->obj(j);
     256      if (coef != 0.0) {
     257        *b = std::make_pair(j, coef);
     258        ++b;
     259      }
     260    }
     261  }
     262
     263  void LpSoplex::_setObjCoeff(int i, Value obj_coef) {
     264    soplex->changeObj(i, obj_coef);
     265  }
     266
     267  LpSoplex::Value LpSoplex::_getObjCoeff(int i) const {
     268    return soplex->obj(i);
     269  }
     270
     271  LpSoplex::SolveExitStatus LpSoplex::_solve() {
     272   
     273    _clear_temporals();
     274
     275    soplex::SPxSolver::Status status = soplex->solve();
     276
     277    switch (status) {
     278    case soplex::SPxSolver::OPTIMAL:
     279    case soplex::SPxSolver::INFEASIBLE:
     280    case soplex::SPxSolver::UNBOUNDED:
     281      return SOLVED;
     282    default:
     283      return UNSOLVED;
     284    }
     285  }
     286
     287  LpSoplex::Value LpSoplex::_getPrimal(int i) const {
     288    if (_primal_values.empty()) {
     289      _primal_values.resize(soplex->nCols());
     290      soplex::Vector pv(_primal_values.size(), &_primal_values.front());
     291      soplex->getPrimal(pv);
     292    }
     293    return _primal_values[i];
     294  }
     295
     296  LpSoplex::Value LpSoplex::_getDual(int i) const {
     297    if (_dual_values.empty()) {
     298      _dual_values.resize(soplex->nRows());
     299      soplex::Vector dv(_dual_values.size(), &_dual_values.front());
     300      soplex->getDual(dv);
     301    }
     302    return _dual_values[i];
     303  }
     304
     305  LpSoplex::Value LpSoplex::_getPrimalValue() const {
     306    return soplex->objValue();
     307  }
     308
     309  LpSoplex::VarStatus LpSoplex::_getColStatus(int i) const {
     310    switch (soplex->getBasisColStatus(i)) {
     311    case soplex::SPxSolver::BASIC:
     312      return BASIC;
     313    case soplex::SPxSolver::ON_UPPER:
     314      return UPPER;
     315    case soplex::SPxSolver::ON_LOWER:
     316      return LOWER;
     317    case soplex::SPxSolver::FIXED:
     318      return FIXED;
     319    case soplex::SPxSolver::ZERO:
     320      return FREE;
     321    default:
     322      LEMON_ASSERT(false, "Wrong column status");
     323      return VarStatus();
     324    }
     325  }
     326
     327  LpSoplex::VarStatus LpSoplex::_getRowStatus(int i) const {
     328    switch (soplex->getBasisRowStatus(i)) {
     329    case soplex::SPxSolver::BASIC:
     330      return BASIC;
     331    case soplex::SPxSolver::ON_UPPER:
     332      return UPPER;
     333    case soplex::SPxSolver::ON_LOWER:
     334      return LOWER;
     335    case soplex::SPxSolver::FIXED:
     336      return FIXED;
     337    case soplex::SPxSolver::ZERO:
     338      return FREE;
     339    default:
     340      LEMON_ASSERT(false, "Wrong row status");
     341      return VarStatus();
     342    }
     343  }
     344
     345  LpSoplex::Value LpSoplex::_getPrimalRay(int i) const {
     346    if (_primal_ray.empty()) {
     347      _primal_ray.resize(soplex->nCols());
     348      soplex::Vector pv(_primal_ray.size(), &_primal_ray.front());
     349      soplex->getDualfarkas(pv);
     350    }
     351    return _primal_ray[i];
     352  }
     353
     354  LpSoplex::Value LpSoplex::_getDualRay(int i) const {
     355    if (_dual_ray.empty()) {
     356      _dual_ray.resize(soplex->nRows());
     357      soplex::Vector dv(_dual_ray.size(), &_dual_ray.front());
     358      soplex->getDualfarkas(dv);
     359    }
     360    return _dual_ray[i];
     361  }
     362
     363  LpSoplex::ProblemType LpSoplex::_getPrimalType() const {
     364    switch (soplex->status()) {
     365    case soplex::SPxSolver::OPTIMAL:
     366      return OPTIMAL;
     367    case soplex::SPxSolver::UNBOUNDED:
     368      return UNBOUNDED;
     369    case soplex::SPxSolver::INFEASIBLE:
     370      return INFEASIBLE;
     371    default:
     372      return UNDEFINED;
     373    }
     374  }
     375
     376  LpSoplex::ProblemType LpSoplex::_getDualType() const {
     377    switch (soplex->status()) {
     378    case soplex::SPxSolver::OPTIMAL:
     379      return OPTIMAL;
     380    case soplex::SPxSolver::UNBOUNDED:
     381      return UNBOUNDED;
     382    case soplex::SPxSolver::INFEASIBLE:
     383      return INFEASIBLE;
     384    default:
     385      return UNDEFINED;
     386    }
     387  }
     388
     389  void LpSoplex::_setSense(Sense sense) {
     390    switch (sense) {
     391    case MIN:
     392      soplex->changeSense(soplex::SPxSolver::MINIMIZE);
     393      break;
     394    case MAX:
     395      soplex->changeSense(soplex::SPxSolver::MAXIMIZE);
     396    }
     397  }
     398
     399  LpSoplex::Sense LpSoplex::_getSense() const {
     400    switch (soplex->spxSense()) {
     401    case soplex::SPxSolver::MAXIMIZE:
     402      return MAX;
     403    case soplex::SPxSolver::MINIMIZE:
     404      return MIN;
     405    default:
     406      LEMON_ASSERT(false, "Wrong sense.");
     407      return LpSoplex::Sense();
     408    }
     409  }
     410
     411  void LpSoplex::_clear() {
     412    soplex->clear();
     413    _col_names.clear();
     414    _col_names_ref.clear();
     415    _row_names.clear();
     416    _row_names_ref.clear();
     417    cols.clear();
     418    rows.clear();
     419    _clear_temporals();
     420  }
     421
     422} //namespace lemon
     423
  • new file lemon/lp_soplex.h

    diff -r cbe3ec2d59d2 -r aa6dff56b2e1 lemon/lp_soplex.h
    - +  
     1/* -*- C++ -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library
     4 *
     5 * Copyright (C) 2003-2008
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#ifndef LEMON_LP_SOPLEX_H
     20#define LEMON_LP_SOPLEX_H
     21
     22///\file
     23///\brief Header of the LEMON-SOPLEX lp solver interface.
     24
     25#include <vector>
     26#include <string>
     27
     28#include <lemon/lp_base.h>
     29
     30// Forward declaration
     31namespace soplex {
     32  class SoPlex;
     33}
     34
     35namespace lemon {
     36
     37  /// \ingroup lp_group
     38  ///
     39  /// \brief Interface for the SOPLEX solver
     40  ///
     41  /// This class implements an interface for the SoPlex LP solver.
     42  /// The SoPlex library is an object oriented lp solver library
     43  /// developed at the Konrad-Zuse-Zentrum für Informationstechnik
     44  /// Berlin (ZIB). You can find detailed information about it at the
     45  /// <tt>http://soplex.zib.de</tt> address.
     46  class LpSoplex : public LpSolverBase {
     47  private:
     48
     49    soplex::SoPlex* soplex;
     50
     51    std::vector<std::string> _col_names;
     52    std::map<std::string, int> _col_names_ref;
     53
     54    std::vector<std::string> _row_names;
     55    std::map<std::string, int> _row_names_ref;
     56
     57  private:
     58
     59    // these values cannot be retrieved element by element
     60    mutable std::vector<Value> _primal_values;
     61    mutable std::vector<Value> _dual_values;
     62
     63    mutable std::vector<Value> _primal_ray;
     64    mutable std::vector<Value> _dual_ray;
     65
     66    void _clear_temporals();
     67
     68  public:
     69
     70    typedef LpSolverBase Parent;
     71   
     72
     73    /// \e
     74    LpSoplex();
     75    /// \e
     76    LpSoplex(const LpSoplex&);
     77    /// \e
     78    ~LpSoplex();
     79
     80  protected:
     81
     82    virtual LpSoplex* _newSolver() const;
     83    virtual LpSoplex* _cloneSolver() const ;
     84
     85    virtual int _addCol();
     86    virtual int _addRow();
     87
     88    virtual void _eraseCol(int i);
     89    virtual void _eraseRow(int i);
     90
     91    virtual void _eraseColId(int i);
     92    virtual void _eraseRowId(int i);
     93
     94    virtual void _getColName(int col, std::string& name) const;
     95    virtual void _setColName(int col, const std::string& name);
     96    virtual int _colByName(const std::string& name) const;
     97
     98    virtual void _getRowName(int row, std::string& name) const;
     99    virtual void _setRowName(int row, const std::string& name);
     100    virtual int _rowByName(const std::string& name) const;
     101
     102    virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e);
     103    virtual void _getRowCoeffs(int i, InsertIterator b) const;
     104
     105    virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e);
     106    virtual void _getColCoeffs(int i, InsertIterator b) const;
     107
     108    virtual void _setCoeff(int row, int col, Value value);
     109    virtual Value _getCoeff(int row, int col) const;
     110
     111    virtual void _setColLowerBound(int i, Value value);
     112    virtual Value _getColLowerBound(int i) const;
     113    virtual void _setColUpperBound(int i, Value value);
     114    virtual Value _getColUpperBound(int i) const;
     115
     116    virtual void _setRowLowerBound(int i, Value value);
     117    virtual Value _getRowLowerBound(int i) const;
     118    virtual void _setRowUpperBound(int i, Value value);
     119    virtual Value _getRowUpperBound(int i) const;
     120
     121    virtual void _setObjCoeffs(ExprIterator b, ExprIterator e);
     122    virtual void _getObjCoeffs(InsertIterator b) const;
     123
     124    virtual void _setObjCoeff(int i, Value obj_coef);
     125    virtual Value _getObjCoeff(int i) const;
     126   
     127    virtual void _setSense(Sense sense);
     128    virtual Sense _getSense() const;
     129
     130    virtual SolveExitStatus _solve();
     131    virtual Value _getPrimal(int i) const;
     132    virtual Value _getDual(int i) const;
     133
     134    virtual Value _getPrimalValue() const;
     135
     136    virtual Value _getPrimalRay(int i) const;
     137    virtual Value _getDualRay(int i) const;
     138
     139    virtual VarStatus _getColStatus(int i) const;
     140    virtual VarStatus _getRowStatus(int i) const;
     141   
     142    virtual ProblemType _getPrimalType() const;
     143    virtual ProblemType _getDualType() const;
     144   
     145    virtual void _clear();
     146
     147  };
     148
     149} //END OF NAMESPACE LEMON
     150
     151#endif //LEMON_LP_SOPLEX_H
     152
  • new file m4/lx_check_clp.m4

    diff -r cbe3ec2d59d2 -r aa6dff56b2e1 m4/lx_check_clp.m4
    - +  
     1AC_DEFUN([LX_CHECK_CLP],
     2[
     3  AC_ARG_WITH([clp],
     4AS_HELP_STRING([--with-clp@<:@=PREFIX@:>@], [search for CLP under PREFIX or under the default search paths if PREFIX is not given @<:@default@:>@])
     5AS_HELP_STRING([--without-clp], [disable checking for CLP]),
     6              [], [with_clp=yes])
     7
     8  AC_ARG_WITH([clp-includedir],
     9AS_HELP_STRING([--with-clp-includedir=DIR], [search for CLP headers in DIR]),
     10              [], [with_clp_includedir=no])
     11
     12  AC_ARG_WITH([clp-libdir],
     13AS_HELP_STRING([--with-clp-libdir=DIR], [search for CLP libraries in DIR]),
     14              [], [with_clp_libdir=no])
     15
     16  lx_clp_found=no
     17  if test x"$with_clp" != x"no"; then
     18    AC_MSG_CHECKING([for CLP])
     19
     20    if test x"$with_clp_includedir" != x"no"; then
     21      CLP_CXXFLAGS="-I$with_clp_includedir"
     22    elif test x"$with_clp" != x"yes"; then
     23      CLP_CXXFLAGS="-I$with_clp/include"
     24    fi
     25
     26    if test x"$with_clp_libdir" != x"no"; then
     27      CLP_LDFLAGS="-L$with_clp_libdir"
     28    elif test x"$with_clp" != x"yes"; then
     29      CLP_LDFLAGS="-L$with_clp/lib"
     30    fi
     31    CLP_LIBS="-lClp -lCoinUtils -lm"
     32
     33    lx_save_cxxflags="$CXXFLAGS"
     34    lx_save_ldflags="$LDFLAGS"
     35    lx_save_libs="$LIBS"
     36    CXXFLAGS="$CLP_CXXFLAGS"
     37    LDFLAGS="$CLP_LDFLAGS"
     38    LIBS="$CLP_LIBS"
     39
     40    lx_clp_test_prog='
     41      #include <coin/ClpModel.hpp>
     42
     43      int main(int argc, char** argv)
     44      {
     45        ClpModel clp;
     46        return 0;
     47      }'
     48
     49    AC_LANG_PUSH(C++)
     50    AC_LINK_IFELSE([$lx_clp_test_prog], [lx_clp_found=yes], [lx_clp_found=no])
     51    AC_LANG_POP(C++)
     52
     53    CXXFLAGS="$lx_save_cxxflags"
     54    LDFLAGS="$lx_save_ldflags"
     55    LIBS="$lx_save_libs"
     56
     57    if test x"$lx_clp_found" = x"yes"; then
     58      AC_DEFINE([HAVE_CLP], [1], [Define to 1 if you have CLP.])
     59      AC_MSG_RESULT([yes])
     60    else
     61      CLP_CXXFLAGS=""
     62      CLP_LDFLAGS=""
     63      CLP_LIBS=""
     64      AC_MSG_RESULT([no])
     65    fi
     66  fi
     67  CLP_LIBS="$CLP_LDFLAGS $CLP_LIBS"
     68  AC_SUBST(CLP_CXXFLAGS)
     69  AC_SUBST(CLP_LIBS)
     70  AM_CONDITIONAL([HAVE_CLP], [test x"$lx_clp_found" = x"yes"])
     71])
  • test/CMakeLists.txt

    diff -r cbe3ec2d59d2 -r aa6dff56b2e1 test/CMakeLists.txt
    a b  
    1515  graph_utils_test
    1616  heap_test
    1717  kruskal_test
     18  lp_test
    1819  maps_test
     20  mip_test
    1921  random_test
    2022  path_test
    2123  time_measure_test
  • test/Makefile.am

    diff -r cbe3ec2d59d2 -r aa6dff56b2e1 test/Makefile.am
    a b  
    1818        test/graph_utils_test \
    1919        test/heap_test \
    2020        test/kruskal_test \
     21        test/lp_test \
    2122        test/maps_test \
     23        test/mip_test \
    2224        test/random_test \
    2325        test/path_test \
    2426        test/test_tools_fail \
     
    4143test_graph_utils_test_SOURCES = test/graph_utils_test.cc
    4244test_heap_test_SOURCES = test/heap_test.cc
    4345test_kruskal_test_SOURCES = test/kruskal_test.cc
     46test_lp_test_SOURCES = test/lp_test.cc
    4447test_maps_test_SOURCES = test/maps_test.cc
     48test_mip_test_SOURCES = test/mip_test.cc
    4549test_path_test_SOURCES = test/path_test.cc
    4650test_random_test_SOURCES = test/random_test.cc
    4751test_test_tools_fail_SOURCES = test/test_tools_fail.cc
  • new file test/lp_test.cc

    diff -r cbe3ec2d59d2 -r aa6dff56b2e1 test/lp_test.cc
    - +  
     1/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library.
     4 *
     5 * Copyright (C) 2003-2008
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#include <sstream>
     20#include <lemon/lp_skeleton.h>
     21#include "test_tools.h"
     22#include <lemon/tolerance.h>
     23#include <lemon/lp_utils.h>
     24#include <lemon/lgf_writer.h>
     25
     26#ifdef HAVE_CONFIG_H
     27#include <lemon/config.h>
     28#endif
     29
     30#ifdef HAVE_GLPK
     31#include <lemon/lp_glpk.h>
     32#endif
     33
     34#ifdef HAVE_CPLEX
     35#include <lemon/lp_cplex.h>
     36#endif
     37
     38#ifdef HAVE_SOPLEX
     39#include <lemon/lp_soplex.h>
     40#endif
     41
     42#ifdef HAVE_CLP
     43#include <lemon/lp_clp.h>
     44#endif
     45
     46using namespace lemon;
     47
     48void lpTest(LpSolverBase & lp)
     49{
     50
     51  typedef LpSolverBase LP;
     52
     53  std::vector<LP::Col> x(10);
     54  //  for(int i=0;i<10;i++) x.push_back(lp.addCol());
     55  lp.addColSet(x);
     56  lp.colLowerBound(x,1);
     57  lp.colUpperBound(x,1);
     58  lp.colBounds(x,1,2);
     59
     60  std::vector<LP::Col> y(10);
     61  lp.addColSet(y);
     62
     63  lp.colLowerBound(y,1);
     64  lp.colUpperBound(y,1);
     65  lp.colBounds(y,1,2);
     66
     67  std::map<int,LP::Col> z;
     68
     69  z.insert(std::make_pair(12,INVALID));
     70  z.insert(std::make_pair(2,INVALID));
     71  z.insert(std::make_pair(7,INVALID));
     72  z.insert(std::make_pair(5,INVALID));
     73
     74  lp.addColSet(z);
     75
     76  lp.colLowerBound(z,1);
     77  lp.colUpperBound(z,1);
     78  lp.colBounds(z,1,2);
     79
     80  {
     81    LP::Expr e,f,g;
     82    LP::Col p1,p2,p3,p4,p5;
     83    LP::Constr c;
     84
     85    p1=lp.addCol();
     86    p2=lp.addCol();
     87    p3=lp.addCol();
     88    p4=lp.addCol();
     89    p5=lp.addCol();
     90
     91    e[p1]=2;
     92    *e=12;
     93    e[p1]+=2;
     94    *e+=12;
     95    e[p1]-=2;
     96    *e-=12;
     97
     98    e=2;
     99    e=2.2;
     100    e=p1;
     101    e=f;
     102
     103    e+=2;
     104    e+=2.2;
     105    e+=p1;
     106    e+=f;
     107
     108    e-=2;
     109    e-=2.2;
     110    e-=p1;
     111    e-=f;
     112
     113    e*=2;
     114    e*=2.2;
     115    e/=2;
     116    e/=2.2;
     117
     118    e=((p1+p2)+(p1-p2)+(p1+12)+(12+p1)+(p1-12)+(12-p1)+
     119       (f+12)+(12+f)+(p1+f)+(f+p1)+(f+g)+
     120       (f-12)+(12-f)+(p1-f)+(f-p1)+(f-g)+
     121       2.2*f+f*2.2+f/2.2+
     122       2*f+f*2+f/2+
     123       2.2*p1+p1*2.2+p1/2.2+
     124       2*p1+p1*2+p1/2
     125       );
     126
     127
     128    c = (e  <= f  );
     129    c = (e  <= 2.2);
     130    c = (e  <= 2  );
     131    c = (e  <= p1 );
     132    c = (2.2<= f  );
     133    c = (2  <= f  );
     134    c = (p1 <= f  );
     135    c = (p1 <= p2 );
     136    c = (p1 <= 2.2);
     137    c = (p1 <= 2  );
     138    c = (2.2<= p2 );
     139    c = (2  <= p2 );
     140
     141    c = (e  >= f  );
     142    c = (e  >= 2.2);
     143    c = (e  >= 2  );
     144    c = (e  >= p1 );
     145    c = (2.2>= f  );
     146    c = (2  >= f  );
     147    c = (p1 >= f  );
     148    c = (p1 >= p2 );
     149    c = (p1 >= 2.2);
     150    c = (p1 >= 2  );
     151    c = (2.2>= p2 );
     152    c = (2  >= p2 );
     153
     154    c = (e  == f  );
     155    c = (e  == 2.2);
     156    c = (e  == 2  );
     157    c = (e  == p1 );
     158    c = (2.2== f  );
     159    c = (2  == f  );
     160    c = (p1 == f  );
     161    //c = (p1 == p2 );
     162    c = (p1 == 2.2);
     163    c = (p1 == 2  );
     164    c = (2.2== p2 );
     165    c = (2  == p2 );
     166
     167    c = (2 <= e <= 3);
     168    c = (2 <= p1<= 3);
     169
     170    c = (2 >= e >= 3);
     171    c = (2 >= p1>= 3);
     172
     173    e[x[3]]=2;
     174    e[x[3]]=4;
     175    e[x[3]]=1;
     176    *e=12;
     177
     178    lp.addRow(-LP::INF,e,23);
     179    lp.addRow(-LP::INF,3.0*(x[1]+x[2]/2)-x[3],23);
     180    lp.addRow(-LP::INF,3.0*(x[1]+x[2]*2-5*x[3]+12-x[4]/3)+2*x[4]-4,23);
     181
     182    lp.addRow(x[1]+x[3]<=x[5]-3);
     183    lp.addRow(-7<=x[1]+x[3]-12<=3);
     184    lp.addRow(x[1]<=x[5]);
     185
     186    std::ostringstream buf;
     187
     188
     189    e=((p1+p2)+(p1-0.99*p2));
     190    //e.prettyPrint(std::cout);
     191    //(e<=2).prettyPrint(std::cout);
     192    double tolerance=0.001;
     193    e.simplify(tolerance);
     194    buf << "Coeff. of p2 should be 0.01";
     195    check(e[p2]>0, buf.str());
     196
     197    tolerance=0.02;
     198    e.simplify(tolerance);
     199    buf << "Coeff. of p2 should be 0";
     200    check(const_cast<const LpSolverBase::Expr&>(e)[p2]==0, buf.str());
     201
     202
     203  }
     204
     205  {
     206    LP::DualExpr e,f,g;
     207    LP::Row p1 = INVALID, p2 = INVALID, p3 = INVALID,
     208      p4 = INVALID, p5 = INVALID;
     209
     210    e[p1]=2;
     211    e[p1]+=2;
     212    e[p1]-=2;
     213
     214    e=p1;
     215    e=f;
     216
     217    e+=p1;
     218    e+=f;
     219
     220    e-=p1;
     221    e-=f;
     222
     223    e*=2;
     224    e*=2.2;
     225    e/=2;
     226    e/=2.2;
     227
     228    e=((p1+p2)+(p1-p2)+
     229       (p1+f)+(f+p1)+(f+g)+
     230       (p1-f)+(f-p1)+(f-g)+
     231       2.2*f+f*2.2+f/2.2+
     232       2*f+f*2+f/2+
     233       2.2*p1+p1*2.2+p1/2.2+
     234       2*p1+p1*2+p1/2
     235       );
     236  }
     237
     238}
     239
     240void solveAndCheck(LpSolverBase& lp, LpSolverBase::ProblemType stat,
     241                   double exp_opt) {
     242  using std::string;
     243  lp.solve();
     244 
     245  std::ostringstream buf;
     246  buf << "PrimalType should be: " << int(stat) << int(lp.primalType());
     247
     248  check(lp.primalType()==stat, buf.str());
     249
     250  if (stat ==  LpSolverBase::OPTIMAL) {
     251    std::ostringstream sbuf;
     252    sbuf << "Wrong optimal value: the right optimum is " << exp_opt;
     253    check(std::abs(lp.primal()-exp_opt) < 1e-3, sbuf.str());
     254  }
     255}
     256
     257void aTest(LpSolverBase & lp)
     258{
     259  typedef LpSolverBase LP;
     260
     261 //The following example is very simple
     262
     263  typedef LpSolverBase::Row Row;
     264  typedef LpSolverBase::Col Col;
     265
     266
     267  Col x1 = lp.addCol();
     268  Col x2 = lp.addCol();
     269
     270
     271  //Constraints
     272  Row upright=lp.addRow(x1+2*x2 <=1);
     273  lp.addRow(x1+x2 >=-1);
     274  lp.addRow(x1-x2 <=1);
     275  lp.addRow(x1-x2 >=-1);
     276  //Nonnegativity of the variables
     277  lp.colLowerBound(x1, 0);
     278  lp.colLowerBound(x2, 0);
     279  //Objective function
     280  lp.obj(x1+x2);
     281
     282  lp.sense(lp.MAX);
     283
     284  //Testing the problem retrieving routines
     285  check(lp.objCoeff(x1)==1,"First term should be 1 in the obj function!");
     286  check(lp.sense() == lp.MAX,"This is a maximization!");
     287  check(lp.coeff(upright,x1)==1,"The coefficient in question is 1!");
     288  check(lp.colLowerBound(x1)==0,
     289        "The lower bound for variable x1 should be 0.");
     290  check(lp.colUpperBound(x1)==LpSolverBase::INF,
     291        "The upper bound for variable x1 should be infty.");
     292  check(lp.rowLowerBound(upright) == -LpSolverBase::INF,
     293        "The lower bound for the first row should be -infty.");
     294  check(lp.rowUpperBound(upright)==1,
     295        "The upper bound for the first row should be 1.");
     296  LpSolverBase::Expr e = lp.row(upright);
     297  check(e[x1] == 1, "The first coefficient should 1.");
     298  check(e[x2] == 2, "The second coefficient should 1.");
     299
     300  lp.row(upright, x1+x2 <=1);
     301  e = lp.row(upright);
     302  check(e[x1] == 1, "The first coefficient should 1.");
     303  check(e[x2] == 1, "The second coefficient should 1.");
     304
     305  LpSolverBase::DualExpr de = lp.col(x1);
     306  check(  de[upright] == 1, "The first coefficient should 1.");
     307
     308  LpSolverBase* clp = lp.cloneSolver();
     309
     310  //Testing the problem retrieving routines
     311  check(clp->objCoeff(x1)==1,"First term should be 1 in the obj function!");
     312  check(clp->sense() == clp->MAX,"This is a maximization!");
     313  check(clp->coeff(upright,x1)==1,"The coefficient in question is 1!");
     314  //  std::cout<<lp.colLowerBound(x1)<<std::endl;
     315  check(clp->colLowerBound(x1)==0,
     316        "The lower bound for variable x1 should be 0.");
     317  check(clp->colUpperBound(x1)==LpSolverBase::INF,
     318        "The upper bound for variable x1 should be infty.");
     319
     320  check(lp.rowLowerBound(upright)==-LpSolverBase::INF,
     321        "The lower bound for the first row should be -infty.");
     322  check(lp.rowUpperBound(upright)==1,
     323        "The upper bound for the first row should be 1.");
     324  e = clp->row(upright);
     325  check(e[x1] == 1, "The first coefficient should 1.");
     326  check(e[x2] == 1, "The second coefficient should 1.");
     327
     328  de = clp->col(x1);
     329  check(de[upright] == 1, "The first coefficient should 1.");
     330
     331  delete clp;
     332
     333  //Maximization of x1+x2
     334  //over the triangle with vertices (0,0) (0,1) (1,0)
     335  double expected_opt=1;
     336  solveAndCheck(lp, LpSolverBase::OPTIMAL, expected_opt);
     337
     338  //Minimization
     339  lp.sense(lp.MIN);
     340  expected_opt=0;
     341  solveAndCheck(lp, LpSolverBase::OPTIMAL, expected_opt);
     342
     343  //Vertex (-1,0) instead of (0,0)
     344  lp.colLowerBound(x1, -LpSolverBase::INF);
     345  expected_opt=-1;
     346  solveAndCheck(lp, LpSolverBase::OPTIMAL, expected_opt);
     347
     348  //Erase one constraint and return to maximization
     349  lp.erase(upright);
     350  lp.sense(lp.MAX);
     351  expected_opt=LpSolverBase::INF;
     352  solveAndCheck(lp, LpSolverBase::UNBOUNDED, expected_opt);
     353
     354  //Infeasibilty
     355  lp.addRow(x1+x2 <=-2);
     356  solveAndCheck(lp, LpSolverBase::INFEASIBLE, expected_opt);
     357
     358}
     359
     360int main()
     361{
     362  LpSkeleton lp_skel;
     363  lpTest(lp_skel);
     364
     365#ifdef HAVE_GLPK
     366  {
     367    std::cerr << "glpk" << std::endl;
     368    LpGlpk lp_glpk1,lp_glpk2;
     369    lpTest(lp_glpk1);
     370    aTest(lp_glpk2);
     371  }
     372#endif
     373
     374#ifdef HAVE_CPLEX
     375  try {
     376    std::cerr << "cplex" << std::endl;
     377    LpCplex lp_cplex1,lp_cplex2;
     378    lpTest(lp_cplex1);
     379    aTest(lp_cplex2);
     380  } catch (CplexEnv::LicenseError& error) {
     381#ifdef LEMON_FORCE_CPLEX_CHECK
     382    check(false, error.what());
     383#else
     384    std::cerr << error.what() << std::endl;
     385    std::cerr << "Cplex license check failed, lp check skipped" << std::endl;
     386#endif
     387  }
     388#endif
     389
     390#ifdef HAVE_SOPLEX
     391  {
     392    std::cerr << "soplex" << std::endl;
     393    LpSoplex lp_soplex1,lp_soplex2;
     394    lpTest(lp_soplex1);
     395    aTest(lp_soplex2);
     396  }
     397#endif
     398
     399#ifdef HAVE_CLP
     400  {
     401    std::cerr << "clp" << std::endl;
     402    LpClp lp_clp1,lp_clp2;
     403    lpTest(lp_clp1);
     404    aTest(lp_clp2);
     405  }
     406#endif
     407
     408  return 0;
     409}
  • new file test/mip_test.cc

    diff -r cbe3ec2d59d2 -r aa6dff56b2e1 test/mip_test.cc
    - +  
     1/* -*- C++ -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library
     4 *
     5 * Copyright (C) 2003-2008
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#include "test_tools.h"
     20
     21
     22#ifdef HAVE_CONFIG_H
     23#include <lemon/config.h>
     24#endif
     25
     26#ifdef HAVE_CPLEX
     27#include <lemon/lp_cplex.h>
     28#endif
     29
     30#ifdef HAVE_GLPK
     31#include <lemon/lp_glpk.h>
     32#endif
     33
     34
     35using namespace lemon;
     36
     37void solveAndCheck(MipSolverBase& mip, MipSolverBase::ProblemType stat,
     38                   double exp_opt) {
     39  using std::string;
     40
     41  mip.solve();
     42  //int decimal,sign;
     43  std::ostringstream buf;
     44  buf << "Type should be: " << int(stat)<<" and it is "<<int(mip.type());
     45
     46
     47  //  itoa(stat,buf1, 10);
     48  check(mip.type()==stat, buf.str());
     49
     50  if (stat ==  MipSolverBase::OPTIMAL) {
     51    std::ostringstream sbuf;
     52    buf << "Wrong optimal value: the right optimum is " << exp_opt;
     53    check(std::abs(mip.solValue()-exp_opt) < 1e-3, sbuf.str());
     54    //+ecvt(exp_opt,2)
     55  }
     56}
     57
     58void aTest(MipSolverBase& mip)
     59{
     60 //The following example is very simple
     61
     62
     63  typedef MipSolverBase::Row Row;
     64  typedef MipSolverBase::Col Col;
     65
     66
     67
     68  Col x1 = mip.addCol();
     69  Col x2 = mip.addCol();
     70
     71
     72  //Objective function
     73  mip.obj(x1);
     74
     75  mip.max();
     76
     77
     78  //Unconstrained optimization
     79  mip.solve();
     80  //Check it out!
     81
     82  //Constraints
     83  mip.addRow(2*x1+x2 <=2); 
     84  mip.addRow(x1-2*x2 <=0);
     85
     86  //Nonnegativity of the variable x1
     87  mip.colLowerBound(x1, 0);
     88
     89  //Maximization of x1
     90  //over the triangle with vertices (0,0),(4/5,2/5),(0,2)
     91  double expected_opt=4.0/5.0;
     92  solveAndCheck(mip, MipSolverBase::OPTIMAL, expected_opt);
     93
     94  //Restrict x2 to integer
     95  mip.colType(x2,MipSolverBase::INTEGER); 
     96  expected_opt=1.0/2.0;
     97  solveAndCheck(mip, MipSolverBase::OPTIMAL, expected_opt);
     98
     99
     100  //Restrict both to integer
     101  mip.colType(x1,MipSolverBase::INTEGER); 
     102  expected_opt=0;
     103  solveAndCheck(mip, MipSolverBase::OPTIMAL, expected_opt);
     104
     105 
     106
     107}
     108
     109
     110int main()
     111{
     112
     113#ifdef HAVE_GLPK
     114  MipGlpk mip1;
     115  aTest(mip1);
     116#endif
     117
     118#ifdef HAVE_CPLEX
     119  MipCplex mip2;
     120  aTest(mip2);
     121#endif
     122
     123  return 0;
     124
     125}