COIN-OR::LEMON - Graph Library

Changeset 482:ed54c0d13df0 in lemon


Ignore:
Timestamp:
12/02/08 22:48:28 (10 years ago)
Author:
Balazs Dezso <deba@…>
Branch:
default
Children:
483:76ec7bd57026, 547:17cabb114d52
Phase:
public
Message:

Thorough redesign of the LP/MIP interface (#44)

  • Redesigned class structure
  • Redesigned iterators
  • Some functions in the basic interface redesigned
  • More complete setting functions
  • Ray retrieving functions
  • Lot of improvements
  • Cplex common env
  • CLP macro definition to config.h.in
  • Update lp.h to also use soplex and clp
  • Remove default_solver_name
  • New solverName() function in solvers
  • Handle exceptions for MipCplex? test
  • Rename tolerance parameter to epsilon
  • Rename MapIt? to CoeffIt?
  • Lot of documentation improvements
  • Various bugfixes
Files:
4 added
5 deleted
18 edited

Legend:

Unmodified
Added
Removed
  • configure.ac

    r481 r482  
    5454LX_CHECK_CPLEX
    5555LX_CHECK_SOPLEX
     56LX_CHECK_CLP
    5657
    5758AM_CONDITIONAL([HAVE_LP], [test x"$lx_lp_found" = x"yes"])
     
    121122echo CPLEX support................. : $lx_cplex_found
    122123echo SOPLEX support................ : $lx_soplex_found
     124echo CLP support................... : $lx_clp_found
    123125echo
    124126echo Build demo programs........... : $enable_demo
  • lemon/CMakeLists.txt

    r481 r482  
    55  base.cc
    66  color.cc
    7   lp_base.cc
    8   lp_skeleton.cc
    9   lp_utils.cc
    107  random.cc)
    118
  • lemon/Makefile.am

    r481 r482  
    1919        $(GLPK_CFLAGS) \
    2020        $(CPLEX_CFLAGS) \
    21         $(SOPLEX_CXXFLAGS)
     21        $(SOPLEX_CXXFLAGS) \
     22        $(CLP_CXXFLAGS)
    2223
    2324lemon_libemon_la_LDFLAGS = \
    2425        $(GLPK_LIBS) \
    2526        $(CPLEX_LIBS) \
    26         $(SOPLEX_LIBS)
     27        $(SOPLEX_LIBS) \
     28        $(CLP_LIBS)
    2729
    2830if HAVE_GLPK
    29 lemon_libemon_la_SOURCES += lemon/lp_glpk.cc lemon/mip_glpk.cc
     31lemon_libemon_la_SOURCES += lemon/lp_glpk.cc
    3032endif
    3133
    3234if HAVE_CPLEX
    33 lemon_libemon_la_SOURCES += lemon/lp_cplex.cc lemon/mip_cplex.cc
     35lemon_libemon_la_SOURCES += lemon/lp_cplex.cc
    3436endif
    3537
    3638if HAVE_SOPLEX
    3739lemon_libemon_la_SOURCES += lemon/lp_soplex.cc
     40endif
     41
     42if HAVE_CLP
     43lemon_libemon_la_SOURCES += lemon/lp_clp.cc
    3844endif
    3945
     
    6672        lemon/lp.h \
    6773        lemon/lp_base.h \
     74        lemon/lp_clp.h \
    6875        lemon/lp_cplex.h \
    6976        lemon/lp_glpk.h \
    7077        lemon/lp_skeleton.h \
    7178        lemon/lp_soplex.h \
    72         lemon/mip_cplex.h \
    73         lemon/mip_glpk.h \
     79        lemon/list_graph.h \
    7480        lemon/maps.h \
    7581        lemon/math.h \
     
    95101        lemon/bits/graph_adaptor_extender.h \
    96102        lemon/bits/graph_extender.h \
    97         lemon/bits/lp_id.h \
    98103        lemon/bits/map_extender.h \
    99104        lemon/bits/path_dump.h \
     105        lemon/bits/solver_bits.h \
    100106        lemon/bits/traits.h \
    101107        lemon/bits/variant.h \
  • lemon/config.h.in

    r481 r482  
    1313/* Define to 1 if you have SOPLEX */
    1414#undef HAVE_SOPLEX
     15
     16/* Define to 1 if you have CLP */
     17#undef HAVE_CLP
  • lemon/lp.h

    r481 r482  
    2525#ifdef HAVE_GLPK
    2626#include <lemon/lp_glpk.h>
    27 #include <lemon/mip_glpk.h>
    2827#elif HAVE_CPLEX
    2928#include <lemon/lp_cplex.h>
    30 #include <lemon/mip_cplex.h>
    3129#elif HAVE_SOPLEX
    3230#include <lemon/lp_soplex.h>
     31#elif HAVE_CLP
     32#include <lemon/lp_clp.h>
    3333#endif
    3434
     
    4444  ///\ingroup lp_group
    4545  ///
    46   ///Currently, the possible values are \c GLPK or \c CPLEX
    47 #define DEFAULT_LP SOLVER
     46  ///Currently, the possible values are \c LP_GLPK, \c LP_CPLEX, \c
     47  ///LP_SOPLEX or \c LP_CLP
     48#define LEMON_DEFAULT_LP SOLVER
    4849  ///The default LP solver
    4950
     
    5152  ///\ingroup lp_group
    5253  ///
    53   ///Currently, it is either \c LpGlpk or \c LpCplex
     54  ///Currently, it is either \c LpGlpk, \c LpCplex, \c LpSoplex or \c LpClp
    5455  typedef LpGlpk Lp;
    55   ///The default LP solver identifier string
    5656
    57   ///The default LP solver identifier string.
     57  ///The default MIP solver identifier
     58
     59  ///The default MIP solver identifier.
    5860  ///\ingroup lp_group
    5961  ///
    60   ///Currently, the possible values are "GLPK" or "CPLEX"
    61   const char default_solver_name[]="SOLVER";
     62  ///Currently, the possible values are \c MIP_GLPK or \c MIP_CPLEX
     63#define LEMON_DEFAULT_MIP SOLVER
     64  ///The default MIP solver.
    6265
    63   ///The default ILP solver.
    64 
    65   ///The default ILP solver.
     66  ///The default MIP solver.
    6667  ///\ingroup lp_group
    6768  ///
    68   ///Currently, it is either \c LpGlpk or \c LpCplex
     69  ///Currently, it is either \c MipGlpk or \c MipCplex
    6970  typedef MipGlpk Mip;
    7071#else
    7172#ifdef HAVE_GLPK
    72 #define DEFAULT_LP GLPK
     73# define LEMON_DEFAULT_LP LP_GLPK
    7374  typedef LpGlpk Lp;
     75# define LEMON_DEFAULT_MIP MIP_GLPK
    7476  typedef MipGlpk Mip;
    75   const char default_solver_name[]="GLPK";
    7677#elif HAVE_CPLEX
    77 #define DEFAULT_LP CPLEX
     78# define LEMON_DEFAULT_LP LP_CPLEX
    7879  typedef LpCplex Lp;
     80# define LEMON_DEFAULT_MIP MIP_CPLEX
    7981  typedef MipCplex Mip;
    80   const char default_solver_name[]="CPLEX";
    8182#elif HAVE_SOPLEX
    82 #define DEFAULT_LP SOPLEX
     83# define DEFAULT_LP LP_SOPLEX
    8384  typedef LpSoplex Lp;
    84   const char default_solver_name[]="SOPLEX";
     85#elif HAVE_CLP
     86# define DEFAULT_LP LP_CLP
     87  typedef LpClp Lp; 
    8588#endif
    8689#endif
  • lemon/lp_base.cc

    r481 r482  
    2323namespace lemon {
    2424
    25   const LpSolverBase::Value
    26   LpSolverBase::INF = std::numeric_limits<Value>::infinity();
    27   const LpSolverBase::Value
    28   LpSolverBase::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();
     25  const LpBase::Value LpBase::INF = std::numeric_limits<Value>::infinity();
     26  const LpBase::Value LpBase::NaN = std::numeric_limits<Value>::quiet_NaN();
    3427
    3528} //namespace lemon
  • lemon/lp_base.h

    r481 r482  
    2626#include<lemon/math.h>
    2727
     28#include<lemon/error.h>
     29#include<lemon/assert.h>
     30
    2831#include<lemon/core.h>
    29 #include<lemon/bits/lp_id.h>
     32#include<lemon/bits/solver_bits.h>
    3033
    3134///\file
     
    3437namespace lemon {
    3538
    36   /// Function to decide whether a floating point value is finite or not.
    37 
    38   /// Retruns true if the argument is not infinity, minus infinity or NaN.
    39   /// It does the same as the isfinite() function defined by C99.
    40   template <typename T>
    41   bool isFinite(T value)
    42   {
    43     typedef std::numeric_limits<T> Lim;
    44     if ((Lim::has_infinity && (value == Lim::infinity() || value ==
    45                                -Lim::infinity())) ||
    46         ((Lim::has_quiet_NaN || Lim::has_signaling_NaN) && value != value))
    47     {
    48       return false;
    49     }
    50     return true;
    51   }
    52 
    53   ///Common base class for LP solvers
    54 
    55   ///\todo Much more docs
     39  ///Common base class for LP and MIP solvers
     40
     41  ///Usually this class is not used directly, please use one of the concrete
     42  ///implementations of the solver interface.
    5643  ///\ingroup lp_group
    57   class LpSolverBase {
     44  class LpBase {
    5845
    5946  protected:
    6047
    61     _lp_bits::LpId rows;
    62     _lp_bits::LpId cols;
     48    _solver_bits::VarIndex rows;
     49    _solver_bits::VarIndex cols;
    6350
    6451  public:
     
    7562    };
    7663
    77       ///\e
    78     enum SolutionStatus {
    79       ///Feasible solution hasn't been found (but may exist).
    80 
    81       ///\todo NOTFOUND might be a better name.
    82       ///
    83       UNDEFINED = 0,
    84       ///The problem has no feasible solution
    85       INFEASIBLE = 1,
    86       ///Feasible solution found
    87       FEASIBLE = 2,
    88       ///Optimal solution exists and found
    89       OPTIMAL = 3,
    90       ///The cost function is unbounded
    91 
    92       ///\todo Give a feasible solution and an infinite ray (and the
    93       ///corresponding bases)
    94       INFINITE = 4
    95     };
    96 
    97     ///\e The type of the investigated LP problem
    98     enum ProblemTypes {
    99       ///Primal-dual feasible
    100       PRIMAL_DUAL_FEASIBLE = 0,
    101       ///Primal feasible dual infeasible
    102       PRIMAL_FEASIBLE_DUAL_INFEASIBLE = 1,
    103       ///Primal infeasible dual feasible
    104       PRIMAL_INFEASIBLE_DUAL_FEASIBLE = 2,
    105       ///Primal-dual infeasible
    106       PRIMAL_DUAL_INFEASIBLE = 3,
    107       ///Could not determine so far
    108       UNKNOWN = 4
     64    ///Direction of the optimization
     65    enum Sense {
     66      /// Minimization
     67      MIN,
     68      /// Maximization
     69      MAX
    10970    };
    11071
     
    11677    static const Value NaN;
    11778
    118     static inline bool isNaN(const Value& v) { return v!=v; }
    119 
    12079    friend class Col;
    12180    friend class ColIt;
    12281    friend class Row;
     82    friend class RowIt;
    12383
    12484    ///Refer to a column of the LP.
     
    12989    ///other columns.
    13090    ///
    131     ///\todo Document what can one do with a Col (INVALID, comparing,
    132     ///it is similar to Node/Edge)
     91    ///\note This class is similar to other Item types in LEMON, like
     92    ///Node and Arc types in digraph.
    13393    class Col {
     94      friend class LpBase;
    13495    protected:
    135       int id;
    136       friend class LpSolverBase;
    137       friend class MipSolverBase;
    138       explicit Col(int _id) : id(_id) {}
     96      int _id;
     97      explicit Col(int id) : _id(id) {}
    13998    public:
    14099      typedef Value ExprValue;
    141       typedef True LpSolverCol;
     100      typedef True LpCol;
     101      /// Default constructor
     102     
     103      /// \warning The default constructor sets the Col to an
     104      /// undefined value.
    142105      Col() {}
    143       Col(const Invalid&) : id(-1) {}
    144       bool operator< (Col c) const  {return id< c.id;}
    145       bool operator> (Col c) const  {return id> c.id;}
    146       bool operator==(Col c) const  {return id==c.id;}
    147       bool operator!=(Col c) const  {return id!=c.id;}
     106      /// Invalid constructor \& conversion.
     107     
     108      /// This constructor initializes the Col to be invalid.
     109      /// \sa Invalid for more details.     
     110      Col(const Invalid&) : _id(-1) {}
     111      /// Equality operator
     112
     113      /// Two \ref Col "Col"s are equal if and only if they point to
     114      /// the same LP column or both are invalid.
     115      bool operator==(Col c) const  {return _id == c._id;}
     116      /// Inequality operator
     117
     118      /// \sa operator==(Col c)
     119      ///
     120      bool operator!=(Col c) const  {return _id != c._id;}
     121      /// Artificial ordering operator.
     122
     123      /// To allow the use of this object in std::map or similar
     124      /// associative container we require this.
     125      ///
     126      /// \note This operator only have to define some strict ordering of
     127      /// the items; this order has nothing to do with the iteration
     128      /// ordering of the items.
     129      bool operator<(Col c) const  {return _id < c._id;}
    148130    };
    149131
     132    ///Iterator for iterate over the columns of an LP problem
     133
     134    /// Its usage is quite simple, for example you can count the number
     135    /// of columns in an LP \c lp:
     136    ///\code
     137    /// int count=0;
     138    /// for (LpBase::ColIt c(lp); c!=INVALID; ++c) ++count;
     139    ///\endcode
    150140    class ColIt : public Col {
    151       const LpSolverBase *_lp;
     141      const LpBase *_solver;
    152142    public:
     143      /// Default constructor
     144     
     145      /// \warning The default constructor sets the iterator
     146      /// to an undefined value.
    153147      ColIt() {}
    154       ColIt(const LpSolverBase &lp) : _lp(&lp)
     148      /// Sets the iterator to the first Col
     149     
     150      /// Sets the iterator to the first Col.
     151      ///
     152      ColIt(const LpBase &solver) : _solver(&solver)
    155153      {
    156         _lp->cols.firstFix(id);
    157       }
     154        _solver->cols.firstItem(_id);
     155      }
     156      /// Invalid constructor \& conversion
     157     
     158      /// Initialize the iterator to be invalid.
     159      /// \sa Invalid for more details.
    158160      ColIt(const Invalid&) : Col(INVALID) {}
     161      /// Next column
     162     
     163      /// Assign the iterator to the next column.
     164      ///
    159165      ColIt &operator++()
    160166      {
    161         _lp->cols.nextFix(id);
     167        _solver->cols.nextItem(_id);
    162168        return *this;
    163169      }
    164170    };
    165171
    166     static int id(const Col& col) { return col.id; }
    167 
     172    /// \brief Returns the ID of the column.
     173    static int id(const Col& col) { return col._id; }
     174    /// \brief Returns the column with the given ID.
     175    ///
     176    /// \pre The argument should be a valid column ID in the LP problem.
     177    static Col colFromId(int id) { return Col(id); }
    168178
    169179    ///Refer to a row of the LP.
     
    174184    ///other rows.
    175185    ///
    176     ///\todo Document what can one do with a Row (INVALID, comparing,
    177     ///it is similar to Node/Edge)
     186    ///\note This class is similar to other Item types in LEMON, like
     187    ///Node and Arc types in digraph.
    178188    class Row {
     189      friend class LpBase;
    179190    protected:
    180       int id;
    181       friend class LpSolverBase;
    182       explicit Row(int _id) : id(_id) {}
     191      int _id;
     192      explicit Row(int id) : _id(id) {}
    183193    public:
    184194      typedef Value ExprValue;
    185       typedef True LpSolverRow;
     195      typedef True LpRow;
     196      /// Default constructor
     197     
     198      /// \warning The default constructor sets the Row to an
     199      /// undefined value.
    186200      Row() {}
    187       Row(const Invalid&) : id(-1) {}
    188 
    189       bool operator< (Row c) const  {return id< c.id;}
    190       bool operator> (Row c) const  {return id> c.id;}
    191       bool operator==(Row c) const  {return id==c.id;}
    192       bool operator!=(Row c) const  {return id!=c.id;}
     201      /// Invalid constructor \& conversion.
     202     
     203      /// This constructor initializes the Row to be invalid.
     204      /// \sa Invalid for more details.     
     205      Row(const Invalid&) : _id(-1) {}
     206      /// Equality operator
     207
     208      /// Two \ref Row "Row"s are equal if and only if they point to
     209      /// the same LP row or both are invalid.
     210      bool operator==(Row r) const  {return _id == r._id;}
     211      /// Inequality operator
     212     
     213      /// \sa operator==(Row r)
     214      ///
     215      bool operator!=(Row r) const  {return _id != r._id;}
     216      /// Artificial ordering operator.
     217
     218      /// To allow the use of this object in std::map or similar
     219      /// associative container we require this.
     220      ///
     221      /// \note This operator only have to define some strict ordering of
     222      /// the items; this order has nothing to do with the iteration
     223      /// ordering of the items.
     224      bool operator<(Row r) const  {return _id < r._id;}
    193225    };
    194226
     227    ///Iterator for iterate over the rows of an LP problem
     228
     229    /// Its usage is quite simple, for example you can count the number
     230    /// of rows in an LP \c lp:
     231    ///\code
     232    /// int count=0;
     233    /// for (LpBase::RowIt c(lp); c!=INVALID; ++c) ++count;
     234    ///\endcode
    195235    class RowIt : public Row {
    196       const LpSolverBase *_lp;
     236      const LpBase *_solver;
    197237    public:
     238      /// Default constructor
     239     
     240      /// \warning The default constructor sets the iterator
     241      /// to an undefined value.
    198242      RowIt() {}
    199       RowIt(const LpSolverBase &lp) : _lp(&lp)
     243      /// Sets the iterator to the first Row
     244     
     245      /// Sets the iterator to the first Row.
     246      ///
     247      RowIt(const LpBase &solver) : _solver(&solver)
    200248      {
    201         _lp->rows.firstFix(id);
    202       }
     249        _solver->rows.firstItem(_id);
     250      }
     251      /// Invalid constructor \& conversion
     252     
     253      /// Initialize the iterator to be invalid.
     254      /// \sa Invalid for more details.
    203255      RowIt(const Invalid&) : Row(INVALID) {}
     256      /// Next row
     257     
     258      /// Assign the iterator to the next row.
     259      ///
    204260      RowIt &operator++()
    205261      {
    206         _lp->rows.nextFix(id);
     262        _solver->rows.nextItem(_id);
    207263        return *this;
    208264      }
    209265    };
    210266
    211     static int id(const Row& row) { return row.id; }
    212 
    213   protected:
    214 
    215     int _lpId(const Col& c) const {
    216       return cols.floatingId(id(c));
    217     }
    218 
    219     int _lpId(const Row& r) const {
    220       return rows.floatingId(id(r));
    221     }
    222 
    223     Col _item(int i, Col) const {
    224       return Col(cols.fixId(i));
    225     }
    226 
    227     Row _item(int i, Row) const {
    228       return Row(rows.fixId(i));
    229     }
    230 
     267    /// \brief Returns the ID of the row.
     268    static int id(const Row& row) { return row._id; }
     269    /// \brief Returns the row with the given ID.
     270    ///
     271    /// \pre The argument should be a valid row ID in the LP problem.
     272    static Row rowFromId(int id) { return Row(id); }
    231273
    232274  public:
     
    239281    ///There are several ways to access and modify the contents of this
    240282    ///container.
    241     ///- Its it fully compatible with \c std::map<Col,double>, so for expamle
    242     ///if \c e is an Expr and \c v and \c w are of type \ref Col, then you can
    243     ///read and modify the coefficients like
    244     ///these.
    245283    ///\code
    246284    ///e[v]=5;
     
    251289    ///\code
    252290    ///double s=0;
    253     ///for(LpSolverBase::Expr::iterator i=e.begin();i!=e.end();++i)
    254     ///  s+=i->second;
     291    ///for(LpBase::Expr::ConstCoeffIt i(e);i!=INVALID;++i)
     292    ///  s+=*i * primal(i);
    255293    ///\endcode
    256     ///(This code computes the sum of all coefficients).
     294    ///(This code computes the primal value of the expression).
    257295    ///- Numbers (<tt>double</tt>'s)
    258296    ///and variables (\ref Col "Col"s) directly convert to an
     
    263301    ///v*2.1+(3*v+(v*12+w+6)*3)/2
    264302    ///\endcode
    265     ///are valid \ref Expr "Expr"essions.
     303    ///are valid expressions.
    266304    ///The usual assignment operations are also defined.
    267305    ///\code
     
    271309    ///e/=5;
    272310    ///\endcode
    273     ///- The constant member can be set and read by \ref constComp()
     311    ///- The constant member can be set and read by dereference
     312    ///  operator (unary *)
     313    ///
    274314    ///\code
    275     ///e.constComp()=12;
    276     ///double c=e.constComp();
     315    ///*e=12;
     316    ///double c=*e;
    277317    ///\endcode
    278318    ///
    279     ///\note \ref clear() not only sets all coefficients to 0 but also
    280     ///clears the constant components.
    281     ///
    282319    ///\sa Constr
    283     ///
    284     class Expr : public std::map<Col,Value>
    285     {
     320    class Expr {
     321      friend class LpBase;
    286322    public:
    287       typedef LpSolverBase::Col Key;
    288       typedef LpSolverBase::Value Value;
     323      /// The key type of the expression
     324      typedef LpBase::Col Key;
     325      /// The value type of the expression
     326      typedef LpBase::Value Value;
    289327
    290328    protected:
    291       typedef std::map<Col,Value> Base;
    292 
    293329      Value const_comp;
     330      std::map<int, Value> comps;
     331
    294332    public:
    295       typedef True IsLinExpression;
    296       ///\e
    297       Expr() : Base(), const_comp(0) { }
    298       ///\e
    299       Expr(const Key &v) : const_comp(0) {
    300         Base::insert(std::make_pair(v, 1));
    301       }
    302       ///\e
     333      typedef True SolverExpr;
     334      /// Default constructor
     335     
     336      /// Construct an empty expression, the coefficients and
     337      /// the constant component are initialized to zero.
     338      Expr() : const_comp(0) {}
     339      /// Construct an expression from a column
     340
     341      /// Construct an expression, which has a term with \c c variable
     342      /// and 1.0 coefficient.
     343      Expr(const Col &c) : const_comp(0) {
     344        typedef std::map<int, Value>::value_type pair_type;
     345        comps.insert(pair_type(id(c), 1));
     346      }
     347      /// Construct an expression from a constant
     348
     349      /// Construct an expression, which's constant component is \c v.
     350      ///
    303351      Expr(const Value &v) : const_comp(v) {}
    304       ///\e
    305       void set(const Key &v,const Value &c) {
    306         Base::insert(std::make_pair(v, c));
    307       }
    308       ///\e
    309       Value &constComp() { return const_comp; }
    310       ///\e
    311       const Value &constComp() const { return const_comp; }
    312 
    313       ///Removes the components with zero coefficient.
    314       void simplify() {
    315         for (Base::iterator i=Base::begin(); i!=Base::end();) {
    316           Base::iterator j=i;
    317           ++j;
    318           if ((*i).second==0) Base::erase(i);
    319           i=j;
     352      /// Returns the coefficient of the column
     353      Value operator[](const Col& c) const {
     354        std::map<int, Value>::const_iterator it=comps.find(id(c));
     355        if (it != comps.end()) {
     356          return it->second;
     357        } else {
     358          return 0;
    320359        }
    321360      }
    322 
    323       void simplify() const {
    324         const_cast<Expr*>(this)->simplify();
    325       }
    326 
    327       ///Removes the coefficients closer to zero than \c tolerance.
    328       void simplify(double &tolerance) {
    329         for (Base::iterator i=Base::begin(); i!=Base::end();) {
    330           Base::iterator j=i;
    331           ++j;
    332           if (std::fabs((*i).second)<tolerance) Base::erase(i);
    333           i=j;
     361      /// Returns the coefficient of the column
     362      Value& operator[](const Col& c) {
     363        return comps[id(c)];
     364      }
     365      /// Sets the coefficient of the column
     366      void set(const Col &c, const Value &v) {
     367        if (v != 0.0) {
     368          typedef std::map<int, Value>::value_type pair_type;
     369          comps.insert(pair_type(id(c), v));
     370        } else {
     371          comps.erase(id(c));
    334372        }
     373      }
     374      /// Returns the constant component of the expression
     375      Value& operator*() { return const_comp; }
     376      /// Returns the constant component of the expression
     377      const Value& operator*() const { return const_comp; }
     378      /// \brief Removes the coefficients which's absolute value does
     379      /// not exceed \c epsilon. It also sets to zero the constant
     380      /// component, if it does not exceed epsilon in absolute value.
     381      void simplify(Value epsilon = 0.0) {
     382        std::map<int, Value>::iterator it=comps.begin();
     383        while (it != comps.end()) {
     384          std::map<int, Value>::iterator jt=it;
     385          ++jt;
     386          if (std::fabs((*it).second) <= epsilon) comps.erase(it);
     387          it=jt;
     388        }
     389        if (std::fabs(const_comp) <= epsilon) const_comp = 0;
     390      }
     391
     392      void simplify(Value epsilon = 0.0) const {
     393        const_cast<Expr*>(this)->simplify(epsilon);
    335394      }
    336395
    337396      ///Sets all coefficients and the constant component to 0.
    338397      void clear() {
    339         Base::clear();
     398        comps.clear();
    340399        const_comp=0;
    341400      }
    342401
    343       ///\e
     402      ///Compound assignment
    344403      Expr &operator+=(const Expr &e) {
    345         for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
    346           (*this)[j->first]+=j->second;
     404        for (std::map<int, Value>::const_iterator it=e.comps.begin();
     405             it!=e.comps.end(); ++it)
     406          comps[it->first]+=it->second;
    347407        const_comp+=e.const_comp;
    348408        return *this;
    349409      }
    350       ///\e
     410      ///Compound assignment
    351411      Expr &operator-=(const Expr &e) {
    352         for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
    353           (*this)[j->first]-=j->second;
     412        for (std::map<int, Value>::const_iterator it=e.comps.begin();
     413             it!=e.comps.end(); ++it)
     414          comps[it->first]-=it->second;
    354415        const_comp-=e.const_comp;
    355416        return *this;
    356417      }
    357       ///\e
    358       Expr &operator*=(const Value &c) {
    359         for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
    360           j->second*=c;
    361         const_comp*=c;
     418      ///Multiply with a constant
     419      Expr &operator*=(const Value &v) {
     420        for (std::map<int, Value>::iterator it=comps.begin();
     421             it!=comps.end(); ++it)
     422          it->second*=v;
     423        const_comp*=v;
    362424        return *this;
    363425      }
    364       ///\e
     426      ///Division with a constant
    365427      Expr &operator/=(const Value &c) {
    366         for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
    367           j->second/=c;
     428        for (std::map<int, Value>::iterator it=comps.begin();
     429             it!=comps.end(); ++it)
     430          it->second/=c;
    368431        const_comp/=c;
    369432        return *this;
    370433      }
     434
     435      ///Iterator over the expression
     436     
     437      ///The iterator iterates over the terms of the expression.
     438      ///
     439      ///\code
     440      ///double s=0;
     441      ///for(LpBase::Expr::CoeffIt i(e);i!=INVALID;++i)
     442      ///  s+= *i * primal(i);
     443      ///\endcode
     444      class CoeffIt {
     445      private:
     446
     447        std::map<int, Value>::iterator _it, _end;
     448
     449      public:
     450
     451        /// Sets the iterator to the first term
     452       
     453        /// Sets the iterator to the first term of the expression.
     454        ///
     455        CoeffIt(Expr& e)
     456          : _it(e.comps.begin()), _end(e.comps.end()){}
     457
     458        /// Convert the iterator to the column of the term
     459        operator Col() const {
     460          return colFromId(_it->first);
     461        }
     462
     463        /// Returns the coefficient of the term
     464        Value& operator*() { return _it->second; }
     465
     466        /// Returns the coefficient of the term
     467        const Value& operator*() const { return _it->second; }
     468        /// Next term
     469       
     470        /// Assign the iterator to the next term.
     471        ///
     472        CoeffIt& operator++() { ++_it; return *this; }
     473
     474        /// Equality operator
     475        bool operator==(Invalid) const { return _it == _end; }
     476        /// Inequality operator
     477        bool operator!=(Invalid) const { return _it != _end; }
     478      };
     479
     480      /// Const iterator over the expression
     481     
     482      ///The iterator iterates over the terms of the expression.
     483      ///
     484      ///\code
     485      ///double s=0;
     486      ///for(LpBase::Expr::ConstCoeffIt i(e);i!=INVALID;++i)
     487      ///  s+=*i * primal(i);
     488      ///\endcode
     489      class ConstCoeffIt {
     490      private:
     491
     492        std::map<int, Value>::const_iterator _it, _end;
     493
     494      public:
     495
     496        /// Sets the iterator to the first term
     497       
     498        /// Sets the iterator to the first term of the expression.
     499        ///
     500        ConstCoeffIt(const Expr& e)
     501          : _it(e.comps.begin()), _end(e.comps.end()){}
     502
     503        /// Convert the iterator to the column of the term
     504        operator Col() const {
     505          return colFromId(_it->first);
     506        }
     507
     508        /// Returns the coefficient of the term
     509        const Value& operator*() const { return _it->second; }
     510
     511        /// Next term
     512       
     513        /// Assign the iterator to the next term.
     514        ///
     515        ConstCoeffIt& operator++() { ++_it; return *this; }
     516
     517        /// Equality operator
     518        bool operator==(Invalid) const { return _it == _end; }
     519        /// Inequality operator
     520        bool operator!=(Invalid) const { return _it != _end; }
     521      };
    371522
    372523    };
     
    395546    ///  e>=t
    396547    ///\endcode
    397     ///\warning The validity of a constraint is checked only at run time, so
    398     ///e.g. \ref addRow(<tt>x[1]\<=x[2]<=5</tt>) will compile, but will throw
    399     ///an assertion.
     548    ///\warning The validity of a constraint is checked only at run
     549    ///time, so e.g. \ref addRow(<tt>x[1]\<=x[2]<=5</tt>) will
     550    ///compile, but will fail an assertion.
    400551    class Constr
    401552    {
    402553    public:
    403       typedef LpSolverBase::Expr Expr;
     554      typedef LpBase::Expr Expr;
    404555      typedef Expr::Key Key;
    405556      typedef Expr::Value Value;
     
    412563      Constr() : _expr(), _lb(NaN), _ub(NaN) {}
    413564      ///\e
    414       Constr(Value lb,const Expr &e,Value ub) :
     565      Constr(Value lb, const Expr &e, Value ub) :
    415566        _expr(e), _lb(lb), _ub(ub) {}
    416       ///\e
    417       Constr(const Expr &e,Value ub) :
    418         _expr(e), _lb(NaN), _ub(ub) {}
    419       ///\e
    420       Constr(Value lb,const Expr &e) :
    421         _expr(e), _lb(lb), _ub(NaN) {}
    422       ///\e
    423567      Constr(const Expr &e) :
    424568        _expr(e), _lb(NaN), _ub(NaN) {}
     
    454598      ///Is the constraint lower bounded?
    455599      bool lowerBounded() const {
    456         return isFinite(_lb);
     600        return _lb != -INF && !std::isnan(_lb);
    457601      }
    458602      ///Is the constraint upper bounded?
    459603      bool upperBounded() const {
    460         return isFinite(_ub);
     604        return _ub != INF && !std::isnan(_ub);
    461605      }
    462606
     
    471615    ///There are several ways to access and modify the contents of this
    472616    ///container.
    473     ///- Its it fully compatible with \c std::map<Row,double>, so for expamle
    474     ///if \c e is an DualExpr and \c v
    475     ///and \c w are of type \ref Row, then you can
    476     ///read and modify the coefficients like
    477     ///these.
    478617    ///\code
    479618    ///e[v]=5;
     
    484623    ///\code
    485624    ///double s=0;
    486     ///for(LpSolverBase::DualExpr::iterator i=e.begin();i!=e.end();++i)
    487     ///  s+=i->second;
     625    ///for(LpBase::DualExpr::ConstCoeffIt i(e);i!=INVALID;++i)
     626    ///  s+=*i;
    488627    ///\endcode
    489628    ///(This code computes the sum of all coefficients).
     
    496635    ///v*2.1+(3*v+(v*12+w)*3)/2
    497636    ///\endcode
    498     ///are valid \ref DualExpr "DualExpr"essions.
     637    ///are valid \ref DualExpr dual expressions.
    499638    ///The usual assignment operations are also defined.
    500639    ///\code
     
    506645    ///
    507646    ///\sa Expr
    508     ///
    509     class DualExpr : public std::map<Row,Value>
    510     {
     647    class DualExpr {
     648      friend class LpBase;
    511649    public:
    512       typedef LpSolverBase::Row Key;
    513       typedef LpSolverBase::Value Value;
     650      /// The key type of the expression
     651      typedef LpBase::Row Key;
     652      /// The value type of the expression
     653      typedef LpBase::Value Value;
    514654
    515655    protected:
    516       typedef std::map<Row,Value> Base;
     656      std::map<int, Value> comps;
    517657
    518658    public:
    519       typedef True IsLinExpression;
    520       ///\e
    521       DualExpr() : Base() { }
    522       ///\e
    523       DualExpr(const Key &v) {
    524         Base::insert(std::make_pair(v, 1));
    525       }
    526       ///\e
    527       void set(const Key &v,const Value &c) {
    528         Base::insert(std::make_pair(v, c));
    529       }
    530 
    531       ///Removes the components with zero coefficient.
    532       void simplify() {
    533         for (Base::iterator i=Base::begin(); i!=Base::end();) {
    534           Base::iterator j=i;
    535           ++j;
    536           if ((*i).second==0) Base::erase(i);
    537           i=j;
     659      typedef True SolverExpr;
     660      /// Default constructor
     661     
     662      /// Construct an empty expression, the coefficients are
     663      /// initialized to zero.
     664      DualExpr() {}
     665      /// Construct an expression from a row
     666
     667      /// Construct an expression, which has a term with \c r dual
     668      /// variable and 1.0 coefficient.
     669      DualExpr(const Row &r) {
     670        typedef std::map<int, Value>::value_type pair_type;
     671        comps.insert(pair_type(id(r), 1));
     672      }
     673      /// Returns the coefficient of the row
     674      Value operator[](const Row& r) const {
     675        std::map<int, Value>::const_iterator it = comps.find(id(r));
     676        if (it != comps.end()) {
     677          return it->second;
     678        } else {
     679          return 0;
    538680        }
    539681      }
    540 
    541       void simplify() const {
    542         const_cast<DualExpr*>(this)->simplify();
    543       }
    544 
    545       ///Removes the coefficients closer to zero than \c tolerance.
    546       void simplify(double &tolerance) {
    547         for (Base::iterator i=Base::begin(); i!=Base::end();) {
    548           Base::iterator j=i;
    549           ++j;
    550           if (std::fabs((*i).second)<tolerance) Base::erase(i);
    551           i=j;
     682      /// Returns the coefficient of the row
     683      Value& operator[](const Row& r) {
     684        return comps[id(r)];
     685      }
     686      /// Sets the coefficient of the row
     687      void set(const Row &r, const Value &v) {
     688        if (v != 0.0) {
     689          typedef std::map<int, Value>::value_type pair_type;
     690          comps.insert(pair_type(id(r), v));
     691        } else {
     692          comps.erase(id(r));
    552693        }
     694      }
     695      /// \brief Removes the coefficients which's absolute value does
     696      /// not exceed \c epsilon.
     697      void simplify(Value epsilon = 0.0) {
     698        std::map<int, Value>::iterator it=comps.begin();
     699        while (it != comps.end()) {
     700          std::map<int, Value>::iterator jt=it;
     701          ++jt;
     702          if (std::fabs((*it).second) <= epsilon) comps.erase(it);
     703          it=jt;
     704        }
     705      }
     706
     707      void simplify(Value epsilon = 0.0) const {
     708        const_cast<DualExpr*>(this)->simplify(epsilon);
    553709      }
    554710
    555711      ///Sets all coefficients to 0.
    556712      void clear() {
    557         Base::clear();
    558       }
    559 
    560       ///\e
     713        comps.clear();
     714      }
     715      ///Compound assignment
    561716      DualExpr &operator+=(const DualExpr &e) {
    562         for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
    563           (*this)[j->first]+=j->second;
     717        for (std::map<int, Value>::const_iterator it=e.comps.begin();
     718             it!=e.comps.end(); ++it)
     719          comps[it->first]+=it->second;
    564720        return *this;
    565721      }
    566       ///\e
     722      ///Compound assignment
    567723      DualExpr &operator-=(const DualExpr &e) {
    568         for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
    569           (*this)[j->first]-=j->second;
     724        for (std::map<int, Value>::const_iterator it=e.comps.begin();
     725             it!=e.comps.end(); ++it)
     726          comps[it->first]-=it->second;
    570727        return *this;
    571728      }
    572       ///\e
    573       DualExpr &operator*=(const Value &c) {
    574         for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
    575           j->second*=c;
     729      ///Multiply with a constant
     730      DualExpr &operator*=(const Value &v) {
     731        for (std::map<int, Value>::iterator it=comps.begin();
     732             it!=comps.end(); ++it)
     733          it->second*=v;
    576734        return *this;
    577735      }
    578       ///\e
    579       DualExpr &operator/=(const Value &c) {
    580         for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
    581           j->second/=c;
     736      ///Division with a constant
     737      DualExpr &operator/=(const Value &v) {
     738        for (std::map<int, Value>::iterator it=comps.begin();
     739             it!=comps.end(); ++it)
     740          it->second/=v;
    582741        return *this;
    583742      }
     743
     744      ///Iterator over the expression
     745     
     746      ///The iterator iterates over the terms of the expression.
     747      ///
     748      ///\code
     749      ///double s=0;
     750      ///for(LpBase::DualExpr::CoeffIt i(e);i!=INVALID;++i)
     751      ///  s+= *i * dual(i);
     752      ///\endcode
     753      class CoeffIt {
     754      private:
     755
     756        std::map<int, Value>::iterator _it, _end;
     757
     758      public:
     759
     760        /// Sets the iterator to the first term
     761       
     762        /// Sets the iterator to the first term of the expression.
     763        ///
     764        CoeffIt(DualExpr& e)
     765          : _it(e.comps.begin()), _end(e.comps.end()){}
     766
     767        /// Convert the iterator to the row of the term
     768        operator Row() const {
     769          return rowFromId(_it->first);
     770        }
     771
     772        /// Returns the coefficient of the term
     773        Value& operator*() { return _it->second; }
     774
     775        /// Returns the coefficient of the term
     776        const Value& operator*() const { return _it->second; }
     777
     778        /// Next term
     779       
     780        /// Assign the iterator to the next term.
     781        ///
     782        CoeffIt& operator++() { ++_it; return *this; }
     783
     784        /// Equality operator
     785        bool operator==(Invalid) const { return _it == _end; }
     786        /// Inequality operator
     787        bool operator!=(Invalid) const { return _it != _end; }
     788      };
     789
     790      ///Iterator over the expression
     791     
     792      ///The iterator iterates over the terms of the expression.
     793      ///
     794      ///\code
     795      ///double s=0;
     796      ///for(LpBase::DualExpr::ConstCoeffIt i(e);i!=INVALID;++i)
     797      ///  s+= *i * dual(i);
     798      ///\endcode
     799      class ConstCoeffIt {
     800      private:
     801
     802        std::map<int, Value>::const_iterator _it, _end;
     803
     804      public:
     805
     806        /// Sets the iterator to the first term
     807       
     808        /// Sets the iterator to the first term of the expression.
     809        ///
     810        ConstCoeffIt(const DualExpr& e)
     811          : _it(e.comps.begin()), _end(e.comps.end()){}
     812
     813        /// Convert the iterator to the row of the term
     814        operator Row() const {
     815          return rowFromId(_it->first);
     816        }
     817
     818        /// Returns the coefficient of the term
     819        const Value& operator*() const { return _it->second; }
     820
     821        /// Next term
     822       
     823        /// Assign the iterator to the next term.
     824        ///
     825        ConstCoeffIt& operator++() { ++_it; return *this; }
     826
     827        /// Equality operator
     828        bool operator==(Invalid) const { return _it == _end; }
     829        /// Inequality operator
     830        bool operator!=(Invalid) const { return _it != _end; }
     831      };
    584832    };
    585833
    586834
    587   private:
    588 
    589     template <typename _Expr>
    590     class MappedOutputIterator {
     835  protected:
     836
     837    class InsertIterator {
     838    private:
     839
     840      std::map<int, Value>& _host;
     841      const _solver_bits::VarIndex& _index;
     842
    591843    public:
    592 
    593       typedef std::insert_iterator<_Expr> Base;
    594844
    595845      typedef std::output_iterator_tag iterator_category;
     
    599849      typedef void pointer;
    600850
    601       MappedOutputIterator(const Base& _base, const LpSolverBase& _lp)
    602         : base(_base), lp(_lp) {}
    603 
    604       MappedOutputIterator& operator*() {
     851      InsertIterator(std::map<int, Value>& host,
     852                   const _solver_bits::VarIndex& index)
     853        : _host(host), _index(index) {}
     854
     855      InsertIterator& operator=(const std::pair<int, Value>& value) {
     856        typedef std::map<int, Value>::value_type pair_type;
     857        _host.insert(pair_type(_index[value.first], value.second));
    605858        return *this;
    606859      }
    607860
    608       MappedOutputIterator& operator=(const std::pair<int, Value>& value) {
    609         *base = std::make_pair(lp._item(value.first, typename _Expr::Key()),
    610                                value.second);
    611         return *this;
    612       }
    613 
    614       MappedOutputIterator& operator++() {
    615         ++base;
    616         return *this;
    617       }
    618 
    619       MappedOutputIterator operator++(int) {
    620         MappedOutputIterator tmp(*this);
    621         ++base;
    622         return tmp;
    623       }
    624 
    625       bool operator==(const MappedOutputIterator& it) const {
    626         return base == it.base;
    627       }
    628 
    629       bool operator!=(const MappedOutputIterator& it) const {
    630         return base != it.base;
    631       }
    632 
     861      InsertIterator& operator*() { return *this; }
     862      InsertIterator& operator++() { return *this; }
     863      InsertIterator operator++(int) { return *this; }
     864
     865    };
     866
     867    class ExprIterator {
    633868    private:
    634       Base base;
    635       const LpSolverBase& lp;
    636     };
    637 
    638     template <typename Expr>
    639     class MappedInputIterator {
     869      std::map<int, Value>::const_iterator _host_it;
     870      const _solver_bits::VarIndex& _index;
    640871    public:
    641872
    642       typedef typename Expr::const_iterator Base;
    643 
    644       typedef typename Base::iterator_category iterator_category;
    645       typedef typename Base::difference_type difference_type;
     873      typedef std::bidirectional_iterator_tag iterator_category;
     874      typedef std::ptrdiff_t difference_type;
    646875      typedef const std::pair<int, Value> value_type;
    647876      typedef value_type reference;
     877
    648878      class pointer {
    649879      public:
     
    654884      };
    655885
    656       MappedInputIterator(const Base& _base, const LpSolverBase& _lp)
    657         : base(_base), lp(_lp) {}
     886      ExprIterator(const std::map<int, Value>::const_iterator& host_it,
     887                   const _solver_bits::VarIndex& index)
     888        : _host_it(host_it), _index(index) {}
    658889
    659890      reference operator*() {
    660         return std::make_pair(lp._lpId(base->first), base->second);
     891        return std::make_pair(_index(_host_it->first), _host_it->second);
    661892      }
    662893
     
    665896      }
    666897
    667       MappedInputIterator& operator++() {
    668         ++base;
    669         return *this;
    670       }
    671 
    672       MappedInputIterator operator++(int) {
    673         MappedInputIterator tmp(*this);
    674         ++base;
    675         return tmp;
    676       }
    677 
    678       bool operator==(const MappedInputIterator& it) const {
    679         return base == it.base;
    680       }
    681 
    682       bool operator!=(const MappedInputIterator& it) const {
    683         return base != it.base;
    684       }
    685 
    686     private:
    687       Base base;
    688       const LpSolverBase& lp;
     898      ExprIterator& operator++() { ++_host_it; return *this; }
     899      ExprIterator operator++(int) {
     900        ExprIterator tmp(*this); ++_host_it; return tmp;
     901      }
     902
     903      ExprIterator& operator--() { --_host_it; return *this; }
     904      ExprIterator operator--(int) {
     905        ExprIterator tmp(*this); --_host_it; return tmp;
     906      }
     907
     908      bool operator==(const ExprIterator& it) const {
     909        return _host_it == it._host_it;
     910      }
     911
     912      bool operator!=(const ExprIterator& it) const {
     913        return _host_it != it._host_it;
     914      }
     915
    689916    };
    690917
    691918  protected:
    692919
    693     /// STL compatible iterator for lp col
    694     typedef MappedInputIterator<Expr> ConstRowIterator;
    695     /// STL compatible iterator for lp row
    696     typedef MappedInputIterator<DualExpr> ConstColIterator;
    697 
    698     /// STL compatible iterator for lp col
    699     typedef MappedOutputIterator<Expr> RowIterator;
    700     /// STL compatible iterator for lp row
    701     typedef MappedOutputIterator<DualExpr> ColIterator;
    702 
    703920    //Abstract virtual functions
    704     virtual LpSolverBase* _newLp() = 0;
    705     virtual LpSolverBase* _copyLp(){
    706       LpSolverBase* newlp = _newLp();
    707 
    708       std::map<Col, Col> ref;
    709       for (LpSolverBase::ColIt it(*this); it != INVALID; ++it) {
    710         Col ccol = newlp->addCol();
    711         ref[it] = ccol;
    712         newlp->colName(ccol, colName(it));
    713         newlp->colLowerBound(ccol, colLowerBound(it));
    714         newlp->colUpperBound(ccol, colUpperBound(it));
    715       }
    716 
    717       for (LpSolverBase::RowIt it(*this); it != INVALID; ++it) {
    718         Expr e = row(it), ce;
    719         for (Expr::iterator jt = e.begin(); jt != e.end(); ++jt) {
    720           ce[ref[jt->first]] = jt->second;
    721         }
    722         ce += e.constComp();
    723         Row r = newlp->addRow(ce);
    724 
    725         double lower, upper;
    726         getRowBounds(it, lower, upper);
    727         newlp->rowBounds(r, lower, upper);
    728       }
    729 
    730       return newlp;
    731     };
     921    virtual LpBase* _newSolver() const = 0;
     922    virtual LpBase* _cloneSolver() const = 0;
     923
     924    virtual int _addColId(int col) { return cols.addIndex(col); }
     925    virtual int _addRowId(int row) { return rows.addIndex(row); }
     926
     927    virtual void _eraseColId(int col) { cols.eraseIndex(col); }
     928    virtual void _eraseRowId(int row) { rows.eraseIndex(row); }
    732929
    733930    virtual int _addCol() = 0;
     
    737934    virtual void _eraseRow(int row) = 0;
    738935
    739     virtual void _getColName(int col, std::string & name) const = 0;
    740     virtual void _setColName(int col, const std::string & name) = 0;
     936    virtual void _getColName(int col, std::string& name) const = 0;
     937    virtual void _setColName(int col, const std::string& name) = 0;
    741938    virtual int _colByName(const std::string& name) const = 0;
    742939
    743     virtual void _setRowCoeffs(int i, ConstRowIterator b,
    744                                ConstRowIterator e) = 0;
    745     virtual void _getRowCoeffs(int i, RowIterator b) const = 0;
    746     virtual void _setColCoeffs(int i, ConstColIterator b,
    747                                ConstColIterator e) = 0;
    748     virtual void _getColCoeffs(int i, ColIterator b) const = 0;
     940    virtual void _getRowName(int row, std::string& name) const = 0;
     941    virtual void _setRowName(int row, const std::string& name) = 0;
     942    virtual int _rowByName(const std::string& name) const = 0;
     943
     944    virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e) = 0;
     945    virtual void _getRowCoeffs(int i, InsertIterator b) const = 0;
     946
     947    virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e) = 0;
     948    virtual void _getColCoeffs(int i, InsertIterator b) const = 0;
     949
    749950    virtual void _setCoeff(int row, int col, Value value) = 0;
    750951    virtual Value _getCoeff(int row, int col) const = 0;
     952
    751953    virtual void _setColLowerBound(int i, Value value) = 0;
    752954    virtual Value _getColLowerBound(int i) const = 0;
     955
    753956    virtual void _setColUpperBound(int i, Value value) = 0;
    754957    virtual Value _getColUpperBound(int i) const = 0;
    755     virtual void _setRowBounds(int i, Value lower, Value upper) = 0;
    756     virtual void _getRowBounds(int i, Value &lower, Value &upper) const = 0;
     958
     959    virtual void _setRowLowerBound(int i, Value value) = 0;
     960    virtual Value _getRowLowerBound(int i) const = 0;
     961
     962    virtual void _setRowUpperBound(int i, Value value) = 0;
     963    virtual Value _getRowUpperBound(int i) const = 0;
     964
     965    virtual void _setObjCoeffs(ExprIterator b, ExprIterator e) = 0;
     966    virtual void _getObjCoeffs(InsertIterator b) const = 0;
    757967
    758968    virtual void _setObjCoeff(int i, Value obj_coef) = 0;
    759969    virtual Value _getObjCoeff(int i) const = 0;
    760     virtual void _clearObj()=0;
    761 
    762     virtual SolveExitStatus _solve() = 0;
    763     virtual Value _getPrimal(int i) const = 0;
    764     virtual Value _getDual(int i) const = 0;
    765     virtual Value _getPrimalValue() const = 0;
    766     virtual bool _isBasicCol(int i) const = 0;
    767     virtual SolutionStatus _getPrimalStatus() const = 0;
    768     virtual SolutionStatus _getDualStatus() const = 0;
    769     virtual ProblemTypes _getProblemType() const = 0;
    770 
    771     virtual void _setMax() = 0;
    772     virtual void _setMin() = 0;
    773 
    774 
    775     virtual bool _isMax() const = 0;
     970
     971    virtual void _setSense(Sense) = 0;
     972    virtual Sense _getSense() const = 0;
     973
     974    virtual void _clear() = 0;
     975
     976    virtual const char* _solverName() const = 0;
    776977
    777978    //Own protected stuff
     
    780981    Value obj_const_comp;
    781982
     983    LpBase() : rows(), cols(), obj_const_comp(0) {}
     984
    782985  public:
    783986
    784     ///\e
    785     LpSolverBase() : obj_const_comp(0) {}
    786 
    787     ///\e
    788     virtual ~LpSolverBase() {}
     987    /// Virtual destructor
     988    virtual ~LpBase() {}
    789989
    790990    ///Creates a new LP problem
    791     LpSolverBase* newLp() {return _newLp();}
     991    LpBase* newSolver() {return _newSolver();}
    792992    ///Makes a copy of the LP problem
    793     LpSolverBase* copyLp() {return _copyLp();}
     993    LpBase* cloneSolver() {return _cloneSolver();}
     994
     995    ///Gives back the name of the solver.
     996    const char* solverName() const {return _solverName();}
    794997
    795998    ///\name Build up and modify the LP
     
    7981001
    7991002    ///Add a new empty column (i.e a new variable) to the LP
    800     Col addCol() { Col c; _addCol(); c.id = cols.addId(); return c;}
    801 
    802     ///\brief Adds several new columns
    803     ///(i.e a variables) at once
    804     ///
    805     ///This magic function takes a container as its argument
    806     ///and fills its elements
    807     ///with new columns (i.e. variables)
     1003    Col addCol() { Col c; c._id = _addColId(_addCol()); return c;}
     1004
     1005    ///\brief Adds several new columns (i.e variables) at once
     1006    ///
     1007    ///This magic function takes a container as its argument and fills
     1008    ///its elements with new columns (i.e. variables)
    8081009    ///\param t can be
    8091010    ///- a standard STL compatible iterable container with
    810     ///\ref Col as its \c values_type
    811     ///like
     1011    ///\ref Col as its \c values_type like
    8121012    ///\code
    813     ///std::vector<LpSolverBase::Col>
    814     ///std::list<LpSolverBase::Col>
     1013    ///std::vector<LpBase::Col>
     1014    ///std::list<LpBase::Col>
    8151015    ///\endcode
    8161016    ///- a standard STL compatible iterable container with
    817     ///\ref Col as its \c mapped_type
    818     ///like
     1017    ///\ref Col as its \c mapped_type like
    8191018    ///\code
    820     ///std::map<AnyType,LpSolverBase::Col>
     1019    ///std::map<AnyType,LpBase::Col>
    8211020    ///\endcode
    8221021    ///- an iterable lemon \ref concepts::WriteMap "write map" like
    8231022    ///\code
    824     ///ListGraph::NodeMap<LpSolverBase::Col>
    825     ///ListGraph::EdgeMap<LpSolverBase::Col>
     1023    ///ListGraph::NodeMap<LpBase::Col>
     1024    ///ListGraph::ArcMap<LpBase::Col>
    8261025    ///\endcode
    8271026    ///\return The number of the created column.
     
    8311030#else
    8321031    template<class T>
    833     typename enable_if<typename T::value_type::LpSolverCol,int>::type
     1032    typename enable_if<typename T::value_type::LpCol,int>::type
    8341033    addColSet(T &t,dummy<0> = 0) {
    8351034      int s=0;
     
    8381037    }
    8391038    template<class T>
    840     typename enable_if<typename T::value_type::second_type::LpSolverCol,
     1039    typename enable_if<typename T::value_type::second_type::LpCol,
    8411040                       int>::type
    8421041    addColSet(T &t,dummy<1> = 1) {
     
    8491048    }
    8501049    template<class T>
    851     typename enable_if<typename T::MapIt::Value::LpSolverCol,
     1050    typename enable_if<typename T::MapIt::Value::LpCol,
    8521051                       int>::type
    8531052    addColSet(T &t,dummy<2> = 2) {
     
    8671066    ///\param e is a dual linear expression (see \ref DualExpr)
    8681067    ///a better one.
    869     void col(Col c,const DualExpr &e) {
     1068    void col(Col c, const DualExpr &e) {
    8701069      e.simplify();
    871       _setColCoeffs(_lpId(c), ConstColIterator(e.begin(), *this),
    872                     ConstColIterator(e.end(), *this));
     1070      _setColCoeffs(cols(id(c)), ExprIterator(e.comps.begin(), cols),
     1071                    ExprIterator(e.comps.end(), cols));
    8731072    }
    8741073
    8751074    ///Get a column (i.e a dual constraint) of the LP
    8761075
    877     ///\param r is the column to get
     1076    ///\param c is the column to get
    8781077    ///\return the dual expression associated to the column
    8791078    DualExpr col(Col c) const {
    8801079      DualExpr e;
    881       _getColCoeffs(_lpId(c), ColIterator(std::inserter(e, e.end()), *this));
     1080      _getColCoeffs(cols(id(c)), InsertIterator(e.comps, rows));
    8821081      return e;
    8831082    }
     
    8861085
    8871086    ///\param e is a dual linear expression (see \ref DualExpr)
    888     ///\param obj is the corresponding component of the objective
     1087    ///\param o is the corresponding component of the objective
    8891088    ///function. It is 0 by default.
    8901089    ///\return The created column.
     
    9001099    ///This function adds a new empty row (i.e a new constraint) to the LP.
    9011100    ///\return The created row
    902     Row addRow() { Row r; _addRow(); r.id = rows.addId(); return r;}
    903 
    904     ///\brief Add several new rows
    905     ///(i.e a constraints) at once
    906     ///
    907     ///This magic function takes a container as its argument
    908     ///and fills its elements
    909     ///with new row (i.e. variables)
     1101    Row addRow() { Row r; r._id = _addRowId(_addRow()); return r;}
     1102
     1103    ///\brief Add several new rows (i.e constraints) at once
     1104    ///
     1105    ///This magic function takes a container as its argument and fills
     1106    ///its elements with new row (i.e. variables)
    9101107    ///\param t can be
    9111108    ///- a standard STL compatible iterable container with
    912     ///\ref Row as its \c values_type
    913     ///like
     1109    ///\ref Row as its \c values_type like
    9141110    ///\code
    915     ///std::vector<LpSolverBase::Row>
    916     ///std::list<LpSolverBase::Row>
     1111    ///std::vector<LpBase::Row>
     1112    ///std::list<LpBase::Row>
    9171113    ///\endcode
    9181114    ///- a standard STL compatible iterable container with
    919     ///\ref Row as its \c mapped_type
    920     ///like
     1115    ///\ref Row as its \c mapped_type like
    9211116    ///\code
    922     ///std::map<AnyType,LpSolverBase::Row>
     1117    ///std::map<AnyType,LpBase::Row>
    9231118    ///\endcode
    9241119    ///- an iterable lemon \ref concepts::WriteMap "write map" like
    9251120    ///\code
    926     ///ListGraph::NodeMap<LpSolverBase::Row>
    927     ///ListGraph::EdgeMap<LpSolverBase::Row>
     1121    ///ListGraph::NodeMap<LpBase::Row>
     1122    ///ListGraph::ArcMap<LpBase::Row>
    9281123    ///\endcode
    9291124    ///\return The number of rows created.
     
    9331128#else
    9341129    template<class T>
    935     typename enable_if<typename T::value_type::LpSolverRow,int>::type
    936     addRowSet(T &t,dummy<0> = 0) {
     1130    typename enable_if<typename T::value_type::LpRow,int>::type
     1131    addRowSet(T &t, dummy<0> = 0) {
    9371132      int s=0;
    9381133      for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addRow();s++;}
     
    9401135    }
    9411136    template<class T>
    942     typename enable_if<typename T::value_type::second_type::LpSolverRow,
    943                        int>::type
    944     addRowSet(T &t,dummy<1> = 1) {
     1137    typename enable_if<typename T::value_type::second_type::LpRow, int>::type
     1138    addRowSet(T &t, dummy<1> = 1) {
    9451139      int s=0;
    9461140      for(typename T::iterator i=t.begin();i!=t.end();++i) {
     
    9511145    }
    9521146    template<class T>
    953     typename enable_if<typename T::MapIt::Value::LpSolverRow,
    954                        int>::type
    955     addRowSet(T &t,dummy<2> = 2) {
     1147    typename enable_if<typename T::MapIt::Value::LpRow, int>::type
     1148    addRowSet(T &t, dummy<2> = 2) {
    9561149      int s=0;
    9571150      for(typename T::MapIt i(t); i!=INVALID; ++i)
     
    9701163    ///\param e is a linear expression (see \ref Expr)
    9711164    ///\param u is the upper bound (\ref INF means no bound)
    972     ///\bug This is a temporary function. The interface will change to
    973     ///a better one.
    974     ///\todo Option to control whether a constraint with a single variable is
    975     ///added or not.
    9761165    void row(Row r, Value l, const Expr &e, Value u) {
    9771166      e.simplify();
    978       _setRowCoeffs(_lpId(r), ConstRowIterator(e.begin(), *this),
    979                     ConstRowIterator(e.end(), *this));
    980       _setRowBounds(_lpId(r),l-e.constComp(),u-e.constComp());
     1167      _setRowCoeffs(rows(id(r)), ExprIterator(e.comps.begin(), cols),
     1168                    ExprIterator(e.comps.end(), cols));
     1169      _setRowLowerBound(rows(id(r)),l - *e);
     1170      _setRowUpperBound(rows(id(r)),u - *e);
    9811171    }
    9821172
     
    9971187    Expr row(Row r) const {
    9981188      Expr e;
    999       _getRowCoeffs(_lpId(r), RowIterator(std::inserter(e, e.end()), *this));
     1189      _getRowCoeffs(rows(id(r)), InsertIterator(e.comps, cols));
    10001190      return e;
    10011191    }
     
    10071197    ///\param u is the upper bound (\ref INF means no bound)
    10081198    ///\return The created row.
    1009     ///\bug This is a temporary function. The interface will change to
    1010     ///a better one.
    10111199    Row addRow(Value l,const Expr &e, Value u) {
    10121200      Row r=addRow();
     
    10241212      return r;
    10251213    }
    1026     ///Erase a coloumn (i.e a variable) from the LP
    1027 
    1028     ///\param c is the coloumn to be deleted
    1029     ///\todo Please check this
    1030     void eraseCol(Col c) {
    1031       _eraseCol(_lpId(c));
    1032       cols.eraseId(c.id);
    1033     }
    1034     ///Erase a  row (i.e a constraint) from the LP
     1214    ///Erase a column (i.e a variable) from the LP
     1215
     1216    ///\param c is the column to be deleted
     1217    void erase(Col c) {
     1218      _eraseCol(cols(id(c)));
     1219      _eraseColId(cols(id(c)));
     1220    }
     1221    ///Erase a row (i.e a constraint) from the LP
    10351222
    10361223    ///\param r is the row to be deleted
    1037     ///\todo Please check this
    1038     void eraseRow(Row r) {
    1039       _eraseRow(_lpId(r));
    1040       rows.eraseId(r.id);
     1224    void erase(Row r) {
     1225      _eraseRow(rows(id(r)));
     1226      _eraseRowId(rows(id(r)));
    10411227    }
    10421228
    10431229    /// Get the name of a column
    10441230
    1045     ///\param c is the coresponding coloumn
     1231    ///\param c is the coresponding column
    10461232    ///\return The name of the colunm
    10471233    std::string colName(Col c) const {
    10481234      std::string name;
    1049       _getColName(_lpId(c), name);
     1235      _getColName(cols(id(c)), name);
    10501236      return name;
    10511237    }
     
    10531239    /// Set the name of a column
    10541240
    1055     ///\param c is the coresponding coloumn
     1241    ///\param c is the coresponding column
    10561242    ///\param name The name to be given
    10571243    void colName(Col c, const std::string& name) {
    1058       _setColName(_lpId(c), name);
     1244      _setColName(cols(id(c)), name);
    10591245    }
    10601246
     
    10651251    Col colByName(const std::string& name) const {
    10661252      int k = _colByName(name);
    1067       return k != -1 ? Col(cols.fixId(k)) : Col(INVALID);
     1253      return k != -1 ? Col(cols[k]) : Col(INVALID);
     1254    }
     1255
     1256    /// Get the name of a row
     1257
     1258    ///\param r is the coresponding row
     1259    ///\return The name of the row
     1260    std::string rowName(Row r) const {
     1261      std::string name;
     1262      _getRowName(rows(id(r)), name);
     1263      return name;
     1264    }
     1265
     1266    /// Set the name of a row
     1267
     1268    ///\param r is the coresponding row
     1269    ///\param name The name to be given
     1270    void rowName(Row r, const std::string& name) {
     1271      _setRowName(rows(id(r)), name);
     1272    }
     1273
     1274    /// Get the row by its name
     1275
     1276    ///\param name The name of the row
     1277    ///\return the proper row or \c INVALID
     1278    Row rowByName(const std::string& name) const {
     1279      int k = _rowByName(name);
     1280      return k != -1 ? Row(rows[k]) : Row(INVALID);
    10681281    }
    10691282
     
    10711284
    10721285    ///\param r is the row of the element to be modified
    1073     ///\param c is the coloumn of the element to be modified
     1286    ///\param c is the column of the element to be modified
    10741287    ///\param val is the new value of the coefficient
    1075 
    10761288    void coeff(Row r, Col c, Value val) {
    1077       _setCoeff(_lpId(r),_lpId(c), val);
     1289      _setCoeff(rows(id(r)),cols(id(c)), val);
    10781290    }
    10791291
    10801292    /// Get an element of the coefficient matrix of the LP
    10811293
    1082     ///\param r is the row of the element in question
    1083     ///\param c is the coloumn of the element in question
     1294    ///\param r is the row of the element
     1295    ///\param c is the column of the element
    10841296    ///\return the corresponding coefficient
    1085 
    10861297    Value coeff(Row r, Col c) const {
    1087       return _getCoeff(_lpId(r),_lpId(c));
     1298      return _getCoeff(rows(id(r)),cols(id(c)));
    10881299    }
    10891300
     
    10941305    /// Value or -\ref INF.
    10951306    void colLowerBound(Col c, Value value) {
    1096       _setColLowerBound(_lpId(c),value);
     1307      _setColLowerBound(cols(id(c)),value);
    10971308    }
    10981309
    10991310    /// Get the lower bound of a column (i.e a variable)
    11001311
    1101     /// This function returns the lower bound for column (variable) \t c
     1312    /// This function returns the lower bound for column (variable) \c c
    11021313    /// (this might be -\ref INF as well).
    1103     ///\return The lower bound for coloumn \t c
     1314    ///\return The lower bound for column \c c
    11041315    Value colLowerBound(Col c) const {
    1105       return _getColLowerBound(_lpId(c));
     1316      return _getColLowerBound(cols(id(c)));
    11061317    }
    11071318
    11081319    ///\brief Set the lower bound of  several columns
    1109     ///(i.e a variables) at once
     1320    ///(i.e variables) at once
    11101321    ///
    11111322    ///This magic function takes a container as its argument
    11121323    ///and applies the function on all of its elements.
    1113     /// The lower bound of a variable (column) has to be given by an
    1114     /// extended number of type Value, i.e. a finite number of type
    1115     /// Value or -\ref INF.
     1324    ///The lower bound of a variable (column) has to be given by an
     1325    ///extended number of type Value, i.e. a finite number of type
     1326    ///Value or -\ref INF.
    11161327#ifdef DOXYGEN
    11171328    template<class T>
     
    11191330#else
    11201331    template<class T>
    1121     typename enable_if<typename T::value_type::LpSolverCol,void>::type
     1332    typename enable_if<typename T::value_type::LpCol,void>::type
    11221333    colLowerBound(T &t, Value value,dummy<0> = 0) {
    11231334      for(typename T::iterator i=t.begin();i!=t.end();++i) {
     
    11261337    }
    11271338    template<class T>
    1128     typename enable_if<typename T::value_type::second_type::LpSolverCol,
     1339    typename enable_if<typename T::value_type::second_type::LpCol,
    11291340                       void>::type
    11301341    colLowerBound(T &t, Value value,dummy<1> = 1) {
     
    11341345    }
    11351346    template<class T>
    1136     typename enable_if<typename T::MapIt::Value::LpSolverCol,
     1347    typename enable_if<typename T::MapIt::Value::LpCol,
    11371348                       void>::type
    11381349    colLowerBound(T &t, Value value,dummy<2> = 2) {
     
    11491360    /// Value or \ref INF.
    11501361    void colUpperBound(Col c, Value value) {
    1151       _setColUpperBound(_lpId(c),value);
     1362      _setColUpperBound(cols(id(c)),value);
    11521363    };
    11531364
    11541365    /// Get the upper bound of a column (i.e a variable)
    11551366
    1156     /// This function returns the upper bound for column (variable) \t c
     1367    /// This function returns the upper bound for column (variable) \c c
    11571368    /// (this might be \ref INF as well).
    1158     ///\return The upper bound for coloumn \t c
     1369    /// \return The upper bound for column \c c
    11591370    Value colUpperBound(Col c) const {
    1160       return _getColUpperBound(_lpId(c));
     1371      return _getColUpperBound(cols(id(c)));
    11611372    }
    11621373
    11631374    ///\brief Set the upper bound of  several columns
    1164     ///(i.e a variables) at once
     1375    ///(i.e variables) at once
    11651376    ///
    11661377    ///This magic function takes a container as its argument
    11671378    ///and applies the function on all of its elements.
    1168     /// The upper bound of a variable (column) has to be given by an
    1169     /// extended number of type Value, i.e. a finite number of type
    1170     /// Value or \ref INF.
     1379    ///The upper bound of a variable (column) has to be given by an
     1380    ///extended number of type Value, i.e. a finite number of type
     1381    ///Value or \ref INF.
    11711382#ifdef DOXYGEN
    11721383    template<class T>
     
    11741385#else
    11751386    template<class T>
    1176     typename enable_if<typename T::value_type::LpSolverCol,void>::type
     1387    typename enable_if<typename T::value_type::LpCol,void>::type
    11771388    colUpperBound(T &t, Value value,dummy<0> = 0) {
    11781389      for(typename T::iterator i=t.begin();i!=t.end();++i) {
     
    11811392    }
    11821393    template<class T>
    1183     typename enable_if<typename T::value_type::second_type::LpSolverCol,
     1394    typename enable_if<typename T::value_type::second_type::LpCol,
    11841395                       void>::type
    11851396    colUpperBound(T &t, Value value,dummy<1> = 1) {
     
    11891400    }
    11901401    template<class T>
    1191     typename enable_if<typename T::MapIt::Value::LpSolverCol,
     1402    typename enable_if<typename T::MapIt::Value::LpCol,
    11921403                       void>::type
    11931404    colUpperBound(T &t, Value value,dummy<2> = 2) {
     
    12051416    /// Value, -\ref INF or \ref INF.
    12061417    void colBounds(Col c, Value lower, Value upper) {
    1207       _setColLowerBound(_lpId(c),lower);
    1208       _setColUpperBound(_lpId(c),upper);
     1418      _setColLowerBound(cols(id(c)),lower);
     1419      _setColUpperBound(cols(id(c)),upper);
    12091420    }
    12101421
    12111422    ///\brief Set the lower and the upper bound of several columns
    1212     ///(i.e a variables) at once
     1423    ///(i.e variables) at once
    12131424    ///
    12141425    ///This magic function takes a container as its argument
     
    12231434#else
    12241435    template<class T>
    1225     typename enable_if<typename T::value_type::LpSolverCol,void>::type
     1436    typename enable_if<typename T::value_type::LpCol,void>::type
    12261437    colBounds(T &t, Value lower, Value upper,dummy<0> = 0) {
    12271438      for(typename T::iterator i=t.begin();i!=t.end();++i) {
     
    12301441    }
    12311442    template<class T>
    1232     typename enable_if<typename T::value_type::second_type::LpSolverCol,
    1233                        void>::type
     1443    typename enable_if<typename T::value_type::second_type::LpCol, void>::type
    12341444    colBounds(T &t, Value lower, Value upper,dummy<1> = 1) {
    12351445      for(typename T::iterator i=t.begin();i!=t.end();++i) {
     
    12381448    }
    12391449    template<class T>
    1240     typename enable_if<typename T::MapIt::Value::LpSolverCol,
    1241                        void>::type
     1450    typename enable_if<typename T::MapIt::Value::LpCol, void>::type
    12421451    colBounds(T &t, Value lower, Value upper,dummy<2> = 2) {
    12431452      for(typename T::MapIt i(t); i!=INVALID; ++i){
     
    12471456#endif
    12481457
    1249 
    1250     /// Set the lower and the upper bounds of a row (i.e a constraint)
    1251 
    1252     /// The lower and the upper bound of a constraint (row) have to be
    1253     /// given by an extended number of type Value, i.e. a finite
    1254     /// number of type Value, -\ref INF or \ref INF. There is no
    1255     /// separate function for the lower and the upper bound because
    1256     /// that would have been hard to implement for CPLEX.
    1257     void rowBounds(Row c, Value lower, Value upper) {
    1258       _setRowBounds(_lpId(c),lower, upper);
    1259     }
    1260 
    1261     /// Get the lower and the upper bounds of a row (i.e a constraint)
    1262 
    1263     /// The lower and the upper bound of
    1264     /// a constraint (row) are
    1265     /// extended numbers of type Value, i.e.  finite numbers of type
    1266     /// Value, -\ref INF or \ref INF.
    1267     /// \todo There is no separate function for the
    1268     /// lower and the upper bound because we had problems with the
    1269     /// implementation of the setting functions for CPLEX:
    1270     /// check out whether this can be done for these functions.
    1271     void getRowBounds(Row c, Value &lower, Value &upper) const {
    1272       _getRowBounds(_lpId(c),lower, upper);
     1458    /// Set the lower bound of a row (i.e a constraint)
     1459
     1460    /// The lower bound of a constraint (row) has to be given by an
     1461    /// extended number of type Value, i.e. a finite number of type
     1462    /// Value or -\ref INF.
     1463    void rowLowerBound(Row r, Value value) {
     1464      _setRowLowerBound(rows(id(r)),value);
     1465    }
     1466
     1467    /// Get the lower bound of a row (i.e a constraint)
     1468
     1469    /// This function returns the lower bound for row (constraint) \c c
     1470    /// (this might be -\ref INF as well).
     1471    ///\return The lower bound for row \c r
     1472    Value rowLowerBound(Row r) const {
     1473      return _getRowLowerBound(rows(id(r)));
     1474    }
     1475
     1476    /// Set the upper bound of a row (i.e a constraint)
     1477
     1478    /// The upper bound of a constraint (row) has to be given by an
     1479    /// extended number of type Value, i.e. a finite number of type
     1480    /// Value or -\ref INF.
     1481    void rowUpperBound(Row r, Value value) {
     1482      _setRowUpperBound(rows(id(r)),value);
     1483    }
     1484
     1485    /// Get the upper bound of a row (i.e a constraint)
     1486
     1487    /// This function returns the upper bound for row (constraint) \c c
     1488    /// (this might be -\ref INF as well).
     1489    ///\return The upper bound for row \c r
     1490    Value rowUpperBound(Row r) const {
     1491      return _getRowUpperBound(rows(id(r)));
    12731492    }
    12741493
    12751494    ///Set an element of the objective function
    1276     void objCoeff(Col c, Value v) {_setObjCoeff(_lpId(c),v); };
     1495    void objCoeff(Col c, Value v) {_setObjCoeff(cols(id(c)),v); };
    12771496
    12781497    ///Get an element of the objective function
    1279     Value objCoeff(Col c) const { return _getObjCoeff(_lpId(c)); };
     1498    Value objCoeff(Col c) const { return _getObjCoeff(cols(id(c))); };
    12801499
    12811500    ///Set the objective function
    12821501
    12831502    ///\param e is a linear expression of type \ref Expr.
    1284     void obj(Expr e) {
    1285       _clearObj();
    1286       for (Expr::iterator i=e.begin(); i!=e.end(); ++i)
    1287         objCoeff((*i).first,(*i).second);
    1288       obj_const_comp=e.constComp();
     1503    ///
     1504    void obj(const Expr& e) {
     1505      _setObjCoeffs(ExprIterator(e.comps.begin(), cols),
     1506                    ExprIterator(e.comps.end(), cols));
     1507      obj_const_comp = *e;
    12891508    }
    12901509
    12911510    ///Get the objective function
    12921511
    1293     ///\return the objective function as a linear expression of type \ref Expr.
     1512    ///\return the objective function as a linear expression of type
     1513    ///Expr.
    12941514    Expr obj() const {
    12951515      Expr e;
    1296       for (ColIt it(*this); it != INVALID; ++it) {
    1297         double c = objCoeff(it);
    1298         if (c != 0.0) {
    1299           e.insert(std::make_pair(it, c));
    1300         }
    1301       }
     1516      _getObjCoeffs(InsertIterator(e.comps, cols));
     1517      *e = obj_const_comp;
    13021518      return e;
    13031519    }
    13041520
    13051521
    1306     ///Maximize
    1307     void max() { _setMax(); }
    1308     ///Minimize
    1309     void min() { _setMin(); }
    1310 
    1311     ///Query function: is this a maximization problem?
    1312     bool isMax() const {return _isMax(); }
    1313 
    1314     ///Query function: is this a minimization problem?
    1315     bool isMin() const {return !isMax(); }
     1522    ///Set the direction of optimization
     1523    void sense(Sense sense) { _setSense(sense); }
     1524
     1525    ///Query the direction of the optimization
     1526    Sense sense() const {return _getSense(); }
     1527
     1528    ///Set the sense to maximization
     1529    void max() { _setSense(MAX); }
     1530
     1531    ///Set the sense to maximization
     1532    void min() { _setSense(MIN); }
     1533
     1534    ///Clears the problem
     1535    void clear() { _clear(); }
    13161536
    13171537    ///@}
    13181538
     1539  };
     1540
     1541  /// Addition
     1542
     1543  ///\relates LpBase::Expr
     1544  ///
     1545  inline LpBase::Expr operator+(const LpBase::Expr &a, const LpBase::Expr &b) {
     1546    LpBase::Expr tmp(a);
     1547    tmp+=b;
     1548    return tmp;
     1549  }
     1550  ///Substraction
     1551
     1552  ///\relates LpBase::Expr
     1553  ///
     1554  inline LpBase::Expr operator-(const LpBase::Expr &a, const LpBase::Expr &b) {
     1555    LpBase::Expr tmp(a);
     1556    tmp-=b;
     1557    return tmp;
     1558  }
     1559  ///Multiply with constant
     1560
     1561  ///\relates LpBase::Expr
     1562  ///
     1563  inline LpBase::Expr operator*(const LpBase::Expr &a, const LpBase::Value &b) {
     1564    LpBase::Expr tmp(a);
     1565    tmp*=b;
     1566    return tmp;
     1567  }
     1568
     1569  ///Multiply with constant
     1570
     1571  ///\relates LpBase::Expr
     1572  ///
     1573  inline LpBase::Expr operator*(const LpBase::Value &a, const LpBase::Expr &b) {
     1574    LpBase::Expr tmp(b);
     1575    tmp*=a;
     1576    return tmp;
     1577  }
     1578  ///Divide with constant
     1579
     1580  ///\relates LpBase::Expr
     1581  ///
     1582  inline LpBase::Expr operator/(const LpBase::Expr &a, const LpBase::Value &b) {
     1583    LpBase::Expr tmp(a);
     1584    tmp/=b;
     1585    return tmp;
     1586  }
     1587
     1588  ///Create constraint
     1589
     1590  ///\relates LpBase::Constr
     1591  ///
     1592  inline LpBase::Constr operator<=(const LpBase::Expr &e,
     1593                                   const LpBase::Expr &f) {
     1594    return LpBase::Constr(0, f - e, LpBase::INF);
     1595  }
     1596
     1597  ///Create constraint
     1598
     1599  ///\relates LpBase::Constr
     1600  ///
     1601  inline LpBase::Constr operator<=(const LpBase::Value &e,
     1602                                   const LpBase::Expr &f) {
     1603    return LpBase::Constr(e, f, LpBase::NaN);
     1604  }
     1605
     1606  ///Create constraint
     1607
     1608  ///\relates LpBase::Constr
     1609  ///
     1610  inline LpBase::Constr operator<=(const LpBase::Expr &e,
     1611                                   const LpBase::Value &f) {
     1612    return LpBase::Constr(- LpBase::INF, e, f);
     1613  }
     1614
     1615  ///Create constraint
     1616
     1617  ///\relates LpBase::Constr
     1618  ///
     1619  inline LpBase::Constr operator>=(const LpBase::Expr &e,
     1620                                   const LpBase::Expr &f) {
     1621    return LpBase::Constr(0, e - f, LpBase::INF);
     1622  }
     1623
     1624
     1625  ///Create constraint
     1626
     1627  ///\relates LpBase::Constr
     1628  ///
     1629  inline LpBase::Constr operator>=(const LpBase::Value &e,
     1630                                   const LpBase::Expr &f) {
     1631    return LpBase::Constr(LpBase::NaN, f, e);
     1632  }
     1633
     1634
     1635  ///Create constraint
     1636
     1637  ///\relates LpBase::Constr
     1638  ///
     1639  inline LpBase::Constr operator>=(const LpBase::Expr &e,
     1640                                   const LpBase::Value &f) {
     1641    return LpBase::Constr(f, e, LpBase::INF);
     1642  }
     1643
     1644  ///Create constraint
     1645
     1646  ///\relates LpBase::Constr
     1647  ///
     1648  inline LpBase::Constr operator==(const LpBase::Expr &e,
     1649                                   const LpBase::Value &f) {
     1650    return LpBase::Constr(f, e, f);
     1651  }
     1652
     1653  ///Create constraint
     1654
     1655  ///\relates LpBase::Constr
     1656  ///
     1657  inline LpBase::Constr operator==(const LpBase::Expr &e,
     1658                                   const LpBase::Expr &f) {
     1659    return LpBase::Constr(0, f - e, 0);
     1660  }
     1661
     1662  ///Create constraint
     1663
     1664  ///\relates LpBase::Constr
     1665  ///
     1666  inline LpBase::Constr operator<=(const LpBase::Value &n,
     1667                                   const LpBase::Constr &c) {
     1668    LpBase::Constr tmp(c);
     1669    LEMON_ASSERT(std::isnan(tmp.lowerBound()), "Wrong LP constraint");
     1670    tmp.lowerBound()=n;
     1671    return tmp;
     1672  }
     1673  ///Create constraint
     1674
     1675  ///\relates LpBase::Constr
     1676  ///
     1677  inline LpBase::Constr operator<=(const LpBase::Constr &c,
     1678                                   const LpBase::Value &n)
     1679  {
     1680    LpBase::Constr tmp(c);
     1681    LEMON_ASSERT(std::isnan(tmp.upperBound()), "Wrong LP constraint");
     1682    tmp.upperBound()=n;
     1683    return tmp;
     1684  }
     1685
     1686  ///Create constraint
     1687
     1688  ///\relates LpBase::Constr
     1689  ///
     1690  inline LpBase::Constr operator>=(const LpBase::Value &n,
     1691                                   const LpBase::Constr &c) {
     1692    LpBase::Constr tmp(c);
     1693    LEMON_ASSERT(std::isnan(tmp.upperBound()), "Wrong LP constraint");
     1694    tmp.upperBound()=n;
     1695    return tmp;
     1696  }
     1697  ///Create constraint
     1698
     1699  ///\relates LpBase::Constr
     1700  ///
     1701  inline LpBase::Constr operator>=(const LpBase::Constr &c,
     1702                                   const LpBase::Value &n)
     1703  {
     1704    LpBase::Constr tmp(c);
     1705    LEMON_ASSERT(std::isnan(tmp.lowerBound()), "Wrong LP constraint");
     1706    tmp.lowerBound()=n;
     1707    return tmp;
     1708  }
     1709
     1710  ///Addition
     1711
     1712  ///\relates LpBase::DualExpr
     1713  ///
     1714  inline LpBase::DualExpr operator+(const LpBase::DualExpr &a,
     1715                                    const LpBase::DualExpr &b) {
     1716    LpBase::DualExpr tmp(a);
     1717    tmp+=b;
     1718    return tmp;
     1719  }
     1720  ///Substraction
     1721
     1722  ///\relates LpBase::DualExpr
     1723  ///
     1724  inline LpBase::DualExpr operator-(const LpBase::DualExpr &a,
     1725                                    const LpBase::DualExpr &b) {
     1726    LpBase::DualExpr tmp(a);
     1727    tmp-=b;
     1728    return tmp;
     1729  }
     1730  ///Multiply with constant
     1731
     1732  ///\relates LpBase::DualExpr
     1733  ///
     1734  inline LpBase::DualExpr operator*(const LpBase::DualExpr &a,
     1735                                    const LpBase::Value &b) {
     1736    LpBase::DualExpr tmp(a);
     1737    tmp*=b;
     1738    return tmp;
     1739  }
     1740
     1741  ///Multiply with constant
     1742
     1743  ///\relates LpBase::DualExpr
     1744  ///
     1745  inline LpBase::DualExpr operator*(const LpBase::Value &a,
     1746                                    const LpBase::DualExpr &b) {
     1747    LpBase::DualExpr tmp(b);
     1748    tmp*=a;
     1749    return tmp;
     1750  }
     1751  ///Divide with constant
     1752
     1753  ///\relates LpBase::DualExpr
     1754  ///
     1755  inline LpBase::DualExpr operator/(const LpBase::DualExpr &a,
     1756                                    const LpBase::Value &b) {
     1757    LpBase::DualExpr tmp(a);
     1758    tmp/=b;
     1759    return tmp;
     1760  }
     1761
     1762  /// \ingroup lp_group
     1763  ///
     1764  /// \brief Common base class for LP solvers
     1765  ///
     1766  /// This class is an abstract base class for LP solvers. This class
     1767  /// provides a full interface for set and modify an LP problem,
     1768  /// solve it and retrieve the solution. You can use one of the
     1769  /// descendants as a concrete implementation, or the \c Lp
     1770  /// default LP solver. However, if you would like to handle LP
     1771  /// solvers as reference or pointer in a generic way, you can use
     1772  /// this class directly.
     1773  class LpSolver : virtual public LpBase {
     1774  public:
     1775
     1776    /// The problem types for primal and dual problems
     1777    enum ProblemType {
     1778      ///Feasible solution hasn't been found (but may exist).
     1779      UNDEFINED = 0,
     1780      ///The problem has no feasible solution
     1781      INFEASIBLE = 1,
     1782      ///Feasible solution found
     1783      FEASIBLE = 2,
     1784      ///Optimal solution exists and found
     1785      OPTIMAL = 3,
     1786      ///The cost function is unbounded
     1787      UNBOUNDED = 4
     1788    };
     1789
     1790    ///The basis status of variables
     1791    enum VarStatus {
     1792      /// The variable is in the basis
     1793      BASIC,
     1794      /// The variable is free, but not basic
     1795      FREE,
     1796      /// The variable has active lower bound
     1797      LOWER,
     1798      /// The variable has active upper bound
     1799      UPPER,
     1800      /// The variable is non-basic and fixed
     1801      FIXED
     1802    };
     1803
     1804  protected:
     1805
     1806    virtual SolveExitStatus _solve() = 0;
     1807
     1808    virtual Value _getPrimal(int i) const = 0;
     1809    virtual Value _getDual(int i) const = 0;
     1810
     1811    virtual Value _getPrimalRay(int i) const = 0;
     1812    virtual Value _getDualRay(int i) const = 0;
     1813
     1814    virtual Value _getPrimalValue() const = 0;
     1815
     1816    virtual VarStatus _getColStatus(int i) const = 0;
     1817    virtual VarStatus _getRowStatus(int i) const = 0;
     1818
     1819    virtual ProblemType _getPrimalType() const = 0;
     1820    virtual ProblemType _getDualType() const = 0;
     1821
     1822  public:
    13191823
    13201824    ///\name Solve the LP
     
    13271831    ///values and their meanings can be found in the documentation of
    13281832    ///\ref SolveExitStatus.
    1329     ///
    1330     ///\todo Which method is used to solve the problem
    13311833    SolveExitStatus solve() { return _solve(); }
    13321834
     
    13371839    ///@{
    13381840
    1339     /// The status of the primal problem (the original LP problem)
    1340     SolutionStatus primalStatus() const {
    1341       return _getPrimalStatus();
    1342     }
    1343 
    1344     /// The status of the dual (of the original LP) problem
    1345     SolutionStatus dualStatus() const {
    1346       return _getDualStatus();
    1347     }
    1348 
    1349     ///The type of the original LP problem
    1350     ProblemTypes problemType() const {
    1351       return _getProblemType();
    1352     }
    1353 
    1354     ///\e
    1355     Value primal(Col c) const { return _getPrimal(_lpId(c)); }
    1356     ///\e
     1841    /// The type of the primal problem
     1842    ProblemType primalType() const {
     1843      return _getPrimalType();
     1844    }
     1845
     1846    /// The type of the dual problem
     1847    ProblemType dualType() const {
     1848      return _getDualType();
     1849    }
     1850
     1851    /// Return the primal value of the column
     1852
     1853    /// Return the primal value of the column.
     1854    /// \pre The problem is solved.
     1855    Value primal(Col c) const { return _getPrimal(cols(id(c))); }
     1856
     1857    /// Return the primal value of the expression
     1858
     1859    /// Return the primal value of the expression, i.e. the dot
     1860    /// product of the primal solution and the expression.
     1861    /// \pre The problem is solved.
    13571862    Value primal(const Expr& e) const {
    1358       double res = e.constComp();
    1359       for (std::map<Col, double>::const_iterator it = e.begin();
    1360            it != e.end(); ++it) {
    1361         res += _getPrimal(_lpId(it->first)) * it->second;
     1863      double res = *e;
     1864      for (Expr::ConstCoeffIt c(e); c != INVALID; ++c) {
     1865        res += *c * primal(c);
    13621866      }
    13631867      return res;
    13641868    }
    1365 
    1366     ///\e
    1367     Value dual(Row r) const { return _getDual(_lpId(r)); }
    1368     ///\e
     1869    /// Returns a component of the primal ray
     1870   
     1871    /// The primal ray is solution of the modified primal problem,
     1872    /// where we change each finite bound to 0, and we looking for a
     1873    /// negative objective value in case of minimization, and positive
     1874    /// objective value for maximization. If there is such solution,
     1875    /// that proofs the unsolvability of the dual problem, and if a
     1876    /// feasible primal solution exists, then the unboundness of
     1877    /// primal problem.
     1878    ///
     1879    /// \pre The problem is solved and the dual problem is infeasible.
     1880    /// \note Some solvers does not provide primal ray calculation
     1881    /// functions.
     1882    Value primalRay(Col c) const { return _getPrimalRay(cols(id(c))); }
     1883
     1884    /// Return the dual value of the row
     1885
     1886    /// Return the dual value of the row.
     1887    /// \pre The problem is solved.
     1888    Value dual(Row r) const { return _getDual(rows(id(r))); }
     1889
     1890    /// Return the dual value of the dual expression
     1891
     1892    /// Return the dual value of the dual expression, i.e. the dot
     1893    /// product of the dual solution and the dual expression.
     1894    /// \pre The problem is solved.
    13691895    Value dual(const DualExpr& e) const {
    13701896      double res = 0.0;
    1371       for (std::map<Row, double>::const_iterator it = e.begin();
    1372            it != e.end(); ++it) {
    1373         res += _getPrimal(_lpId(it->first)) * it->second;
     1897      for (DualExpr::ConstCoeffIt r(e); r != INVALID; ++r) {
     1898        res += *r * dual(r);
    13741899      }
    13751900      return res;
    13761901    }
    13771902
    1378     ///\e
    1379     bool isBasicCol(Col c) const { return _isBasicCol(_lpId(c)); }
    1380 
    1381     ///\e
     1903    /// Returns a component of the dual ray
     1904   
     1905    /// The dual ray is solution of the modified primal problem, where
     1906    /// we change each finite bound to 0 (i.e. the objective function
     1907    /// coefficients in the primal problem), and we looking for a
     1908    /// ositive objective value. If there is such solution, that
     1909    /// proofs the unsolvability of the primal problem, and if a
     1910    /// feasible dual solution exists, then the unboundness of
     1911    /// dual problem.
     1912    ///
     1913    /// \pre The problem is solved and the primal problem is infeasible.
     1914    /// \note Some solvers does not provide dual ray calculation
     1915    /// functions.
     1916    Value dualRay(Row r) const { return _getDualRay(rows(id(r))); }
     1917
     1918    /// Return the basis status of the column
     1919
     1920    /// \see VarStatus
     1921    VarStatus colStatus(Col c) const { return _getColStatus(cols(id(c))); }
     1922
     1923    /// Return the basis status of the row
     1924
     1925    /// \see VarStatus
     1926    VarStatus rowStatus(Row r) const { return _getRowStatus(rows(id(r))); }
     1927
     1928    ///The value of the objective function
    13821929
    13831930    ///\return
     
    13861933    ///- \ref NaN if no primal solution is found.
    13871934    ///- The (finite) objective value if an optimal solution is found.
    1388     Value primalValue() const { return _getPrimalValue()+obj_const_comp;}
     1935    Value primal() const { return _getPrimalValue()+obj_const_comp;}
    13891936    ///@}
    13901937
     1938    LpSolver* newSolver() {return _newSolver();}
     1939    LpSolver* cloneSolver() {return _cloneSolver();}
     1940
     1941  protected:
     1942
     1943    virtual LpSolver* _newSolver() const = 0;
     1944    virtual LpSolver* _cloneSolver() const = 0;
    13911945  };
    13921946
     
    13951949  ///
    13961950  /// \brief Common base class for MIP solvers
    1397   /// \todo Much more docs
    1398   class MipSolverBase : virtual public LpSolverBase{
     1951  ///
     1952  /// This class is an abstract base class for MIP solvers. This class
     1953  /// provides a full interface for set and modify an MIP problem,
     1954  /// solve it and retrieve the solution. You can use one of the
     1955  /// descendants as a concrete implementation, or the \c Lp
     1956  /// default MIP solver. However, if you would like to handle MIP
     1957  /// solvers as reference or pointer in a generic way, you can use
     1958  /// this class directly.
     1959  class MipSolver : virtual public LpBase {
    13991960  public:
    14001961
    1401     ///Possible variable (coloumn) types (e.g. real, integer, binary etc.)
     1962    /// The problem types for MIP problems
     1963    enum ProblemType {
     1964      ///Feasible solution hasn't been found (but may exist).
     1965      UNDEFINED = 0,
     1966      ///The problem has no feasible solution
     1967      INFEASIBLE = 1,
     1968      ///Feasible solution found
     1969      FEASIBLE = 2,
     1970      ///Optimal solution exists and found
     1971      OPTIMAL = 3,
     1972      ///The cost function is unbounded
     1973      ///
     1974      ///The Mip or at least the relaxed problem is unbounded
     1975      UNBOUNDED = 4
     1976    };
     1977
     1978    ///\name Solve the MIP
     1979
     1980    ///@{
     1981
     1982    /// Solve the MIP problem at hand
     1983    ///
     1984    ///\return The result of the optimization procedure. Possible
     1985    ///values and their meanings can be found in the documentation of
     1986    ///\ref SolveExitStatus.
     1987    SolveExitStatus solve() { return _solve(); }
     1988
     1989    ///@}
     1990
     1991    ///\name Setting column type
     1992    ///@{
     1993
     1994    ///Possible variable (column) types (e.g. real, integer, binary etc.)
    14021995    enum ColTypes {
    1403       ///Continuous variable
     1996      ///Continuous variable (default)
    14041997      REAL = 0,
    14051998      ///Integer variable
    1406 
    1407       ///Unfortunately, cplex 7.5 somewhere writes something like
    1408       ///#define INTEGER 'I'
    1409       INT = 1
    1410       ///\todo No support for other types yet.
     1999      INTEGER = 1
    14112000    };
    14122001
    1413     ///Sets the type of the given coloumn to the given type
    1414     ///
    1415     ///Sets the type of the given coloumn to the given type.
     2002    ///Sets the type of the given column to the given type
     2003
     2004    ///Sets the type of the given column to the given type.
     2005    ///
    14162006    void colType(Col c, ColTypes col_type) {
    1417       _colType(_lpId(c),col_type);
     2007      _setColType(cols(id(c)),col_type);
    14182008    }
    14192009
    14202010    ///Gives back the type of the column.
    1421     ///
     2011
    14222012    ///Gives back the type of the column.
     2013    ///
    14232014    ColTypes colType(Col c) const {
    1424       return _colType(_lpId(c));
    1425     }
    1426 
    1427     ///Sets the type of the given Col to integer or remove that property.
    1428     ///
    1429     ///Sets the type of the given Col to integer or remove that property.
    1430     void integer(Col c, bool enable) {
    1431       if (enable)
    1432         colType(c,INT);
    1433       else
    1434         colType(c,REAL);
    1435     }
    1436 
    1437     ///Gives back whether the type of the column is integer or not.
    1438     ///
    1439     ///Gives back the type of the column.
    1440     ///\return true if the column has integer type and false if not.
    1441     bool integer(Col c) const {
    1442       return (colType(c)==INT);
    1443     }
    1444 
    1445     /// The status of the MIP problem
    1446     SolutionStatus mipStatus() const {
    1447       return _getMipStatus();
    1448     }
     2015      return _getColType(cols(id(c)));
     2016    }
     2017    ///@}
     2018
     2019    ///\name Obtain the solution
     2020
     2021    ///@{
     2022
     2023    /// The type of the MIP problem
     2024    ProblemType type() const {
     2025      return _getType();
     2026    }
     2027
     2028    /// Return the value of the row in the solution
     2029
     2030    ///  Return the value of the row in the solution.
     2031    /// \pre The problem is solved.
     2032    Value sol(Col c) const { return _getSol(cols(id(c))); }
     2033
     2034    /// Return the value of the expression in the solution
     2035
     2036    /// Return the value of the expression in the solution, i.e. the
     2037    /// dot product of the solution and the expression.
     2038    /// \pre The problem is solved.
     2039    Value sol(const Expr& e) const {
     2040      double res = *e;
     2041      for (Expr::ConstCoeffIt c(e); c != INVALID; ++c) {
     2042        res += *c * sol(c);
     2043      }
     2044      return res;
     2045    }
     2046    ///The value of the objective function
     2047   
     2048    ///\return
     2049    ///- \ref INF or -\ref INF means either infeasibility or unboundedness
     2050    /// of the problem, depending on whether we minimize or maximize.
     2051    ///- \ref NaN if no primal solution is found.
     2052    ///- The (finite) objective value if an optimal solution is found.
     2053    Value solValue() const { return _getSolValue()+obj_const_comp;}
     2054    ///@}
    14492055
    14502056  protected:
    14512057
    1452     virtual ColTypes _colType(int col) const = 0;
    1453     virtual void _colType(int col, ColTypes col_type) = 0;
    1454     virtual SolutionStatus _getMipStatus() const = 0;
    1455 
     2058    virtual SolveExitStatus _solve() = 0;
     2059    virtual ColTypes _getColType(int col) const = 0;
     2060    virtual void _setColType(int col, ColTypes col_type) = 0;
     2061    virtual ProblemType _getType() const = 0;
     2062    virtual Value _getSol(int i) const = 0;
     2063    virtual Value _getSolValue() const = 0;
     2064
     2065  public:
     2066
     2067    MipSolver* newSolver() {return _newSolver();}
     2068    MipSolver* cloneSolver() {return _cloneSolver();}
     2069
     2070  protected:
     2071
     2072    virtual MipSolver* _newSolver() const = 0;
     2073    virtual MipSolver* _cloneSolver() const = 0;
    14562074  };
    14572075
    1458   ///\relates LpSolverBase::Expr
    1459   ///
    1460   inline LpSolverBase::Expr operator+(const LpSolverBase::Expr &a,
    1461                                       const LpSolverBase::Expr &b)
    1462   {
    1463     LpSolverBase::Expr tmp(a);
    1464     tmp+=b;
    1465     return tmp;
    1466   }
    1467   ///\e
    1468 
    1469   ///\relates LpSolverBase::Expr
    1470   ///
    1471   inline LpSolverBase::Expr operator-(const LpSolverBase::Expr &a,
    1472                                       const LpSolverBase::Expr &b)
    1473   {
    1474     LpSolverBase::Expr tmp(a);
    1475     tmp-=b;
    1476     return tmp;
    1477   }
    1478   ///\e
    1479 
    1480   ///\relates LpSolverBase::Expr
    1481   ///
    1482   inline LpSolverBase::Expr operator*(const LpSolverBase::Expr &a,
    1483                                       const LpSolverBase::Value &b)
    1484   {
    1485     LpSolverBase::Expr tmp(a);
    1486     tmp*=b;
    1487     return tmp;
    1488   }
    1489 
    1490   ///\e
    1491 
    1492   ///\relates LpSolverBase::Expr
    1493   ///
    1494   inline LpSolverBase::Expr operator*(const LpSolverBase::Value &a,
    1495                                       const LpSolverBase::Expr &b)
    1496   {
    1497     LpSolverBase::Expr tmp(b);
    1498     tmp*=a;
    1499     return tmp;
    1500   }
    1501   ///\e
    1502 
    1503   ///\relates LpSolverBase::Expr
    1504   ///
    1505   inline LpSolverBase::Expr operator/(const LpSolverBase::Expr &a,
    1506                                       const LpSolverBase::Value &b)
    1507   {
    1508     LpSolverBase::Expr tmp(a);
    1509     tmp/=b;
    1510     return tmp;
    1511   }
    1512 
    1513   ///\e
    1514 
    1515   ///\relates LpSolverBase::Constr
    1516   ///
    1517   inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,
    1518                                          const LpSolverBase::Expr &f)
    1519   {
    1520     return LpSolverBase::Constr(-LpSolverBase::INF,e-f,0);
    1521   }
    1522 
    1523   ///\e
    1524 
    1525   ///\relates LpSolverBase::Constr
    1526   ///
    1527   inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &e,
    1528                                          const LpSolverBase::Expr &f)
    1529   {
    1530     return LpSolverBase::Constr(e,f);
    1531   }
    1532 
    1533   ///\e
    1534 
    1535   ///\relates LpSolverBase::Constr
    1536   ///
    1537   inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,
    1538                                          const LpSolverBase::Value &f)
    1539   {
    1540     return LpSolverBase::Constr(-LpSolverBase::INF,e,f);
    1541   }
    1542 
    1543   ///\e
    1544 
    1545   ///\relates LpSolverBase::Constr
    1546   ///
    1547   inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,
    1548                                          const LpSolverBase::Expr &f)
    1549   {
    1550     return LpSolverBase::Constr(-LpSolverBase::INF,f-e,0);
    1551   }
    1552 
    1553 
    1554   ///\e
    1555 
    1556   ///\relates LpSolverBase::Constr
    1557   ///
    1558   inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &e,
    1559                                          const LpSolverBase::Expr &f)
    1560   {
    1561     return LpSolverBase::Constr(f,e);
    1562   }
    1563 
    1564 
    1565   ///\e
    1566 
    1567   ///\relates LpSolverBase::Constr
    1568   ///
    1569   inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,
    1570                                          const LpSolverBase::Value &f)
    1571   {
    1572     return LpSolverBase::Constr(f,e,LpSolverBase::INF);
    1573   }
    1574 
    1575   ///\e
    1576 
    1577   ///\relates LpSolverBase::Constr
    1578   ///
    1579   inline LpSolverBase::Constr operator==(const LpSolverBase::Expr &e,
    1580                                          const LpSolverBase::Value &f)
    1581   {
    1582     return LpSolverBase::Constr(f,e,f);
    1583   }
    1584 
    1585   ///\e
    1586 
    1587   ///\relates LpSolverBase::Constr
    1588   ///
    1589   inline LpSolverBase::Constr operator==(const LpSolverBase::Expr &e,
    1590                                          const LpSolverBase::Expr &f)
    1591   {
    1592     return LpSolverBase::Constr(0,e-f,0);
    1593   }
    1594 
    1595   ///\e
    1596 
    1597   ///\relates LpSolverBase::Constr
    1598   ///
    1599   inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &n,
    1600                                          const LpSolverBase::Constr&c)
    1601   {
    1602     LpSolverBase::Constr tmp(c);
    1603     LEMON_ASSERT(LpSolverBase::isNaN(tmp.lowerBound()), "Wrong LP constraint");
    1604     tmp.lowerBound()=n;
    1605     return tmp;
    1606   }
    1607   ///\e
    1608 
    1609   ///\relates LpSolverBase::Constr
    1610   ///
    1611   inline LpSolverBase::Constr operator<=(const LpSolverBase::Constr& c,
    1612                                          const LpSolverBase::Value &n)
    1613   {
    1614     LpSolverBase::Constr tmp(c);
    1615     LEMON_ASSERT(LpSolverBase::isNaN(tmp.upperBound()), "Wrong LP constraint");
    1616     tmp.upperBound()=n;
    1617     return tmp;
    1618   }
    1619 
    1620   ///\e
    1621 
    1622   ///\relates LpSolverBase::Constr
    1623   ///
    1624   inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &n,
    1625                                          const LpSolverBase::Constr&c)
    1626   {
    1627     LpSolverBase::Constr tmp(c);
    1628     LEMON_ASSERT(LpSolverBase::isNaN(tmp.upperBound()), "Wrong LP constraint");
    1629     tmp.upperBound()=n;
    1630     return tmp;
    1631   }
    1632   ///\e
    1633 
    1634   ///\relates LpSolverBase::Constr
    1635   ///
    1636   inline LpSolverBase::Constr operator>=(const LpSolverBase::Constr& c,
    1637                                          const LpSolverBase::Value &n)
    1638   {
    1639     LpSolverBase::Constr tmp(c);
    1640     LEMON_ASSERT(LpSolverBase::isNaN(tmp.lowerBound()), "Wrong LP constraint");
    1641     tmp.lowerBound()=n;
    1642     return tmp;
    1643   }
    1644 
    1645   ///\e
    1646 
    1647   ///\relates LpSolverBase::DualExpr
    1648   ///
    1649   inline LpSolverBase::DualExpr operator+(const LpSolverBase::DualExpr &a,
    1650                                           const LpSolverBase::DualExpr &b)
    1651   {
    1652     LpSolverBase::DualExpr tmp(a);
    1653     tmp+=b;
    1654     return tmp;
    1655   }
    1656   ///\e
    1657 
    1658   ///\relates LpSolverBase::DualExpr
    1659   ///
    1660   inline LpSolverBase::DualExpr operator-(const LpSolverBase::DualExpr &a,
    1661                                           const LpSolverBase::DualExpr &b)
    1662   {
    1663     LpSolverBase::DualExpr tmp(a);
    1664     tmp-=b;
    1665     return tmp;
    1666   }
    1667   ///\e
    1668 
    1669   ///\relates LpSolverBase::DualExpr
    1670   ///
    1671   inline LpSolverBase::DualExpr operator*(const LpSolverBase::DualExpr &a,
    1672                                           const LpSolverBase::Value &b)
    1673   {
    1674     LpSolverBase::DualExpr tmp(a);
    1675     tmp*=b;
    1676     return tmp;
    1677   }
    1678 
    1679   ///\e
    1680 
    1681   ///\relates LpSolverBase::DualExpr
    1682   ///
    1683   inline LpSolverBase::DualExpr operator*(const LpSolverBase::Value &a,
    1684                                           const LpSolverBase::DualExpr &b)
    1685   {
    1686     LpSolverBase::DualExpr tmp(b);
    1687     tmp*=a;
    1688     return tmp;
    1689   }
    1690   ///\e
    1691 
    1692   ///\relates LpSolverBase::DualExpr
    1693   ///
    1694   inline LpSolverBase::DualExpr operator/(const LpSolverBase::DualExpr &a,
    1695                                           const LpSolverBase::Value &b)
    1696   {
    1697     LpSolverBase::DualExpr tmp(a);
    1698     tmp/=b;
    1699     return tmp;
    1700   }
    17012076
    17022077
  • lemon/lp_cplex.cc

    r481 r482  
    1919#include <iostream>
    2020#include <vector>
     21#include <cstring>
     22
    2123#include <lemon/lp_cplex.h>
    2224
     
    3032namespace lemon {
    3133
    32   LpCplex::LpCplex() {
    33     //    env = CPXopenCPLEXdevelop(&status);
    34     env = CPXopenCPLEX(&status);
    35     lp = CPXcreateprob(env, &status, "LP problem");
    36   }
    37 
    38   LpCplex::LpCplex(const LpCplex& cplex) : LpSolverBase() {
    39     env = CPXopenCPLEX(&status);
    40     lp = CPXcloneprob(env, cplex.lp, &status);
     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  CplexBase::CplexBase() : LpBase() {
     73    int status;
     74    _prob = CPXcreateprob(cplexEnv(), &status, "Cplex problem");
     75  }
     76
     77  CplexBase::CplexBase(const CplexEnv& env)
     78    : LpBase(), _env(env) {
     79    int status;
     80    _prob = CPXcreateprob(cplexEnv(), &status, "Cplex problem");
     81  }
     82
     83  CplexBase::CplexBase(const CplexBase& cplex)
     84    : LpBase() {
     85    int status;
     86    _prob = CPXcloneprob(cplexEnv(), cplex._prob, &status);
    4187    rows = cplex.rows;
    4288    cols = cplex.cols;
    4389  }
    4490
    45   LpCplex::~LpCplex() {
    46     CPXfreeprob(env,&lp);
    47     CPXcloseCPLEX(&env);
    48   }
    49 
    50   LpSolverBase* LpCplex::_newLp()
    51   {
    52     //The first approach opens a new environment
    53     return new LpCplex();
    54   }
    55 
    56   LpSolverBase* LpCplex::_copyLp() {
    57     return new LpCplex(*this);
    58   }
    59 
    60   int LpCplex::_addCol()
    61   {
    62     int i = CPXgetnumcols(env, lp);
    63     Value lb[1],ub[1];
    64     lb[0]=-INF;
    65     ub[0]=INF;
    66     status = CPXnewcols(env, lp, 1, NULL, lb, ub, NULL, NULL);
     91  CplexBase::~CplexBase() {
     92    CPXfreeprob(cplexEnv(),&_prob);
     93  }
     94
     95  int CplexBase::_addCol() {
     96    int i = CPXgetnumcols(cplexEnv(), _prob);
     97    double lb = -INF, ub = INF;
     98    CPXnewcols(cplexEnv(), _prob, 1, 0, &lb, &ub, 0, 0);
    6799    return i;
    68100  }
    69101
    70102
    71   int LpCplex::_addRow()
    72   {
    73     //We want a row that is not constrained
    74     char sense[1];
    75     sense[0]='L';//<= constraint
    76     Value rhs[1];
    77     rhs[0]=INF;
    78     int i = CPXgetnumrows(env, lp);
    79     status = CPXnewrows(env, lp, 1, rhs, sense, NULL, NULL);
     103  int CplexBase::_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);
    80108    return i;
    81109  }
    82110
    83111
    84   void LpCplex::_eraseCol(int i) {
    85     CPXdelcols(env, lp, i, i);
    86   }
    87 
    88   void LpCplex::_eraseRow(int i) {
    89     CPXdelrows(env, lp, i, i);
    90   }
    91 
    92   void LpCplex::_getColName(int col, std::string &name) const
    93   {
    94     ///\bug Untested
    95     int storespace;
    96     CPXgetcolname(env, lp, 0, 0, 0, &storespace, col, col);
    97     if (storespace == 0) {
     112  void CplexBase::_eraseCol(int i) {
     113    CPXdelcols(cplexEnv(), _prob, i, i);
     114  }
     115
     116  void CplexBase::_eraseRow(int i) {
     117    CPXdelrows(cplexEnv(), _prob, i, i);
     118  }
     119
     120  void CplexBase::_eraseColId(int i) {
     121    cols.eraseIndex(i);
     122    cols.shiftIndices(i);
     123  }
     124  void CplexBase::_eraseRowId(int i) {
     125    rows.eraseIndex(i);
     126    rows.shiftIndices(i);
     127  }
     128
     129  void CplexBase::_getColName(int col, std::string &name) const {
     130    int size;
     131    CPXgetcolname(cplexEnv(), _prob, 0, 0, 0, &size, col, col);
     132    if (size == 0) {
    98133      name.clear();
    99134      return;
    100135    }
    101136
    102     storespace *= -1;
    103     std::vector<char> buf(storespace);
    104     char *names[1];
    105     int dontcare;
    106     ///\bug return code unchecked for error
    107     CPXgetcolname(env, lp, names, &*buf.begin(), storespace,
    108                   &dontcare, col, col);
    109     name = names[0];
    110   }
    111 
    112   void LpCplex::_setColName(int col, const std::string &name)
    113   {
    114     ///\bug Untested
    115     char *names[1];
    116     names[0] = const_cast<char*>(name.c_str());
    117     ///\bug return code unchecked for error
    118     CPXchgcolname(env, lp, 1, &col, names);
    119   }
    120 
    121   int LpCplex::_colByName(const std::string& name) const
    122   {
     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 CplexBase::_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 CplexBase::_colByName(const std::string& name) const {
    123153    int index;
    124     if (CPXgetcolindex(env, lp,
     154    if (CPXgetcolindex(cplexEnv(), _prob,
    125155                       const_cast<char*>(name.c_str()), &index) == 0) {
    126156      return index;
     
    129159  }
    130160
    131   ///\warning Data at index 0 is ignored in the arrays.
    132   void LpCplex::_setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e)
     161  void CplexBase::_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 CplexBase::_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 CplexBase::_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 CplexBase::_setRowCoeffs(int i, ExprIterator b,
     194                                      ExprIterator e)
    133195  {
    134196    std::vector<int> indices;
     
    136198    std::vector<Value> values;
    137199
    138     for(ConstRowIterator it=b; it!=e; ++it) {
     200    for(ExprIterator it=b; it!=e; ++it) {
    139201      indices.push_back(it->first);
    140202      values.push_back(it->second);
     
    142204    }
    143205
    144     status = CPXchgcoeflist(env, lp, values.size(),
    145                             &rowlist[0], &indices[0], &values[0]);
    146   }
    147 
    148   void LpCplex::_getRowCoeffs(int i, RowIterator b) const {
     206    CPXchgcoeflist(cplexEnv(), _prob, values.size(),
     207                   &rowlist.front(), &indices.front(), &values.front());
     208  }
     209
     210  void CplexBase::_getRowCoeffs(int i, InsertIterator b) const {
    149211    int tmp1, tmp2, tmp3, length;
    150     CPXgetrows(env, lp, &tmp1, &tmp2, 0, 0, 0, &length, i, i);
     212    CPXgetrows(cplexEnv(), _prob, &tmp1, &tmp2, 0, 0, 0, &length, i, i);
    151213
    152214    length = -length;
     
    154216    std::vector<double> values(length);
    155217
    156     CPXgetrows(env, lp, &tmp1, &tmp2, &indices[0], &values[0],
     218    CPXgetrows(cplexEnv(), _prob, &tmp1, &tmp2,
     219               &indices.front(), &values.front(),
    157220               length, &tmp3, i, i);
    158221
     
    161224      ++b;
    162225    }
    163 
    164     /// \todo implement
    165   }
    166 
    167   void LpCplex::_setColCoeffs(int i, ConstColIterator b, ConstColIterator e)
    168   {
     226  }
     227
     228  void CplexBase::_setColCoeffs(int i, ExprIterator b, ExprIterator e) {
    169229    std::vector<int> indices;
    170230    std::vector<int> collist;
    171231    std::vector<Value> values;
    172232
    173     for(ConstColIterator it=b; it!=e; ++it) {
     233    for(ExprIterator it=b; it!=e; ++it) {
    174234      indices.push_back(it->first);
    175235      values.push_back(it->second);
     
    177237    }
    178238
    179     status = CPXchgcoeflist(env, lp, values.size(),
    180                             &indices[0], &collist[0], &values[0]);
    181   }
    182 
    183   void LpCplex::_getColCoeffs(int i, ColIterator b) const {
     239    CPXchgcoeflist(cplexEnv(), _prob, values.size(),
     240                   &indices.front(), &collist.front(), &values.front());
     241  }
     242
     243  void CplexBase::_getColCoeffs(int i, InsertIterator b) const {
    184244
    185245    int tmp1, tmp2, tmp3, length;
    186     CPXgetcols(env, lp, &tmp1, &tmp2, 0, 0, 0, &length, i, i);
     246    CPXgetcols(cplexEnv(), _prob, &tmp1, &tmp2, 0, 0, 0, &length, i, i);
    187247
    188248    length = -length;
     
    190250    std::vector<double> values(length);
    191251
    192     CPXgetcols(env, lp, &tmp1, &tmp2, &indices[0], &values[0],
     252    CPXgetcols(cplexEnv(), _prob, &tmp1, &tmp2,
     253               &indices.front(), &values.front(),
    193254               length, &tmp3, i, i);
    194255
     
    200261  }
    201262
    202   void LpCplex::_setCoeff(int row, int col, Value value)
     263  void CplexBase::_setCoeff(int row, int col, Value value) {
     264    CPXchgcoef(cplexEnv(), _prob, row, col, value);
     265  }
     266
     267  CplexBase::Value CplexBase::_getCoeff(int row, int col) const {
     268    CplexBase::Value value;
     269    CPXgetcoef(cplexEnv(), _prob, row, col, &value);
     270    return value;
     271  }
     272
     273  void CplexBase::_setColLowerBound(int i, Value value) {
     274    const char s = 'L';
     275    CPXchgbds(cplexEnv(), _prob, 1, &i, &s, &value);
     276  }
     277
     278  CplexBase::Value CplexBase::_getColLowerBound(int i) const {
     279    CplexBase::Value res;
     280    CPXgetlb(cplexEnv(), _prob, &res, i, i);
     281    return res <= -CPX_INFBOUND ? -INF : res;
     282  }
     283
     284  void CplexBase::_setColUpperBound(int i, Value value)
    203285  {
    204     CPXchgcoef(env, lp, row, col, value);
    205   }
    206 
    207   LpCplex::Value LpCplex::_getCoeff(int row, int col) const
     286    const char s = 'U';
     287    CPXchgbds(cplexEnv(), _prob, 1, &i, &s, &value);
     288  }
     289
     290  CplexBase::Value CplexBase::_getColUpperBound(int i) const {
     291    CplexBase::Value res;
     292    CPXgetub(cplexEnv(), _prob, &res, i, i);
     293    return res >= CPX_INFBOUND ? INF : res;
     294  }
     295
     296  CplexBase::Value CplexBase::_getRowLowerBound(int i) const {
     297    char s;
     298    CPXgetsense(cplexEnv(), _prob, &s, i, i);
     299    CplexBase::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  CplexBase::Value CplexBase::_getRowUpperBound(int i) const {
     313    char s;
     314    CPXgetsense(cplexEnv(), _prob, &s, i, i);
     315    CplexBase::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 CplexBase::_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 CplexBase::_setRowLowerBound(int i, Value lb)
    208359  {
    209     LpCplex::Value value;
    210     CPXgetcoef(env, lp, row, col, &value);
    211     return value;
    212   }
    213 
    214   void LpCplex::_setColLowerBound(int i, Value value)
     360    LEMON_ASSERT(lb != INF, "Invalid bound");
     361    _set_row_bounds(i, lb, CplexBase::_getRowUpperBound(i));
     362  }
     363
     364  void CplexBase::_setRowUpperBound(int i, Value ub)
    215365  {
    216     int indices[1];
    217     indices[0]=i;
    218     char lu[1];
    219     lu[0]='L';
    220     Value bd[1];
    221     bd[0]=value;
    222     status = CPXchgbds(env, lp, 1, indices, lu, bd);
    223 
    224   }
    225 
    226   LpCplex::Value LpCplex::_getColLowerBound(int i) const
     366
     367    LEMON_ASSERT(ub != -INF, "Invalid bound");
     368    _set_row_bounds(i, CplexBase::_getRowLowerBound(i), ub);
     369  }
     370
     371  void CplexBase::_setObjCoeffs(ExprIterator b, ExprIterator e)
    227372  {
    228     LpCplex::Value x;
    229     CPXgetlb (env, lp, &x, i, i);
    230     if (x <= -CPX_INFBOUND) x = -INF;
    231     return x;
    232   }
    233 
    234   void LpCplex::_setColUpperBound(int i, Value value)
     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 CplexBase::_getObjCoeffs(InsertIterator b) const
    235385  {
    236     int indices[1];
    237     indices[0]=i;
    238     char lu[1];
    239     lu[0]='U';
    240     Value bd[1];
    241     bd[0]=value;
    242     status = CPXchgbds(env, lp, 1, indices, lu, bd);
    243   }
    244 
    245   LpCplex::Value LpCplex::_getColUpperBound(int i) const
     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 CplexBase::_setObjCoeff(int i, Value obj_coef)
    246399  {
    247     LpCplex::Value x;
    248     CPXgetub (env, lp, &x, i, i);
    249     if (x >= CPX_INFBOUND) x = INF;
    250     return x;
    251   }
    252 
    253   //This will be easier to implement
    254   void LpCplex::_setRowBounds(int i, Value lb, Value ub)
    255   {
    256     //Bad parameter
    257     if (lb==INF || ub==-INF) {
    258       //FIXME error
    259     }
    260 
    261     int cnt=1;
    262     int indices[1];
    263     indices[0]=i;
    264     char sense[1];
    265 
    266     if (lb==-INF){
    267       sense[0]='L';
    268       CPXchgsense(env, lp, cnt, indices, sense);
    269       CPXchgcoef(env, lp, i, -1, ub);
    270 
    271     }
    272     else{
    273       if (ub==INF){
    274         sense[0]='G';
    275         CPXchgsense(env, lp, cnt, indices, sense);
    276         CPXchgcoef(env, lp, i, -1, lb);
    277       }
    278       else{
    279         if (lb == ub){
    280           sense[0]='E';
    281           CPXchgsense(env, lp, cnt, indices, sense);
    282           CPXchgcoef(env, lp, i, -1, lb);
    283         }
    284         else{
    285           sense[0]='R';
    286           CPXchgsense(env, lp, cnt, indices, sense);
    287           CPXchgcoef(env, lp, i, -1, lb);
    288           CPXchgcoef(env, lp, i, -2, ub-lb);
    289         }
    290       }
    291     }
    292   }
    293 
    294 //   void LpCplex::_setRowLowerBound(int i, Value value)
    295 //   {
    296 //     //Not implemented, obsolete
    297 //   }
    298 
    299 //   void LpCplex::_setRowUpperBound(int i, Value value)
    300 //   {
    301 //     //Not implemented, obsolete
    302 // //     //TODO Ezt kell meg megirni
    303 // //     //type of the problem
    304 // //     char sense[1];
    305 // //     status = CPXgetsense(env, lp, sense, i, i);
    306 // //     Value rhs[1];
    307 // //     status = CPXgetrhs(env, lp, rhs, i, i);
    308 
    309 // //     switch (sense[0]) {
    310 // //     case 'L'://<= constraint
    311 // //       break;
    312 // //     case 'E'://= constraint
    313 // //       break;
    314 // //     case 'G'://>= constraint
    315 // //       break;
    316 // //     case 'R'://ranged constraint
    317 // //       break;
    318 // //     default: ;
    319 // //       //FIXME error
    320 // //     }
    321 
    322 // //     status = CPXchgcoef(env, lp, i, -2, value_rng);
    323 //   }
    324 
    325   void LpCplex::_getRowBounds(int i, Value &lb, Value &ub) const
    326   {
    327     char sense;
    328     CPXgetsense(env, lp, &sense,i,i);
    329     lb=-INF;
    330     ub=INF;
    331     switch (sense)
    332       {
    333       case 'L':
    334         CPXgetcoef(env, lp, i, -1, &ub);
    335         break;
    336       case 'G':
    337         CPXgetcoef(env, lp, i, -1, &lb);
    338         break;
    339       case 'E':
    340         CPXgetcoef(env, lp, i, -1, &lb);
    341         ub=lb;
    342         break;
    343       case 'R':
    344         CPXgetcoef(env, lp, i, -1, &lb);
    345         Value x;
    346         CPXgetcoef(env, lp, i, -2, &x);
    347         ub=lb+x;
    348         break;
    349       }
    350   }
    351 
    352   void LpCplex::_setObjCoeff(int i, Value obj_coef)
    353   {
    354     CPXchgcoef(env, lp, -1, i, obj_coef);
    355   }
    356 
    357   LpCplex::Value LpCplex::_getObjCoeff(int i) const
     400    CPXchgobj(cplexEnv(), _prob, 1, &i, &obj_coef);
     401  }
     402
     403  CplexBase::Value CplexBase::_getObjCoeff(int i) const
    358404  {
    359405    Value x;
    360     CPXgetcoef(env, lp, -1, i, &x);
     406    CPXgetobj(cplexEnv(), _prob, &x, i, i);
    361407    return x;
    362408  }
    363409
    364   void LpCplex::_clearObj()
    365   {
    366     for (int i=0;i< CPXgetnumcols(env, lp);++i){
    367       CPXchgcoef(env, lp, -1, i, 0);
    368     }
    369 
    370   }
     410  void CplexBase::_setSense(CplexBase::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  CplexBase::Sense CplexBase::_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 CplexBase::Sense();
     430    }
     431  }
     432
     433  void CplexBase::_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    : LpBase(), CplexBase(), LpSolver() {}
     445
     446  LpCplex::LpCplex(const CplexEnv& env)
     447    : LpBase(), CplexBase(env), LpSolver() {}
     448
     449  LpCplex::LpCplex(const LpCplex& other)
     450    : LpBase(), CplexBase(other), LpSolver() {}
     451
     452  LpCplex::~LpCplex() {}
     453
     454  LpCplex* LpCplex::_newSolver() const { return new LpCplex; }
     455  LpCplex* LpCplex::_cloneSolver() const {return new LpCplex(*this); }
     456
     457  const char* LpCplex::_solverName() const { return "LpCplex"; }
     458
     459  void LpCplex::_clear_temporals() {
     460    _col_status.clear();
     461    _row_status.clear();
     462    _primal_ray.clear();
     463    _dual_ray.clear();
     464  }
     465
    371466  // The routine returns zero unless an error occurred during the
    372467  // optimization. Examples of errors include exhausting available
     
    378473  // routines CPXsolninfo, CPXgetstat, and CPXsolution to obtain
    379474  // further information about the status of the optimization.
    380   LpCplex::SolveExitStatus LpCplex::_solve()
    381   {
    382     //CPX_PARAM_LPMETHOD
    383     status = CPXlpopt(env, lp);
    384     //status = CPXprimopt(env, lp);
     475  LpCplex::SolveExitStatus LpCplex::convertStatus(int status) {
    385476#if CPX_VERSION >= 800
    386     if (status)
    387     {
     477    if (status == 0) {
     478      switch (CPXgetstat(cplexEnv(), _prob)) {
     479      case CPX_STAT_OPTIMAL:
     480      case CPX_STAT_INFEASIBLE:
     481      case CPX_STAT_UNBOUNDED:
     482        return SOLVED;
     483      default:
     484        return UNSOLVED;
     485      }
     486    } else {
    388487      return UNSOLVED;
    389488    }
    390     else
    391     {
    392       switch (CPXgetstat(env, lp))
    393       {
    394         case CPX_STAT_OPTIMAL:
    395         case CPX_STAT_INFEASIBLE:
    396         case CPX_STAT_UNBOUNDED:
    397           return SOLVED;
    398         default:
    399           return UNSOLVED;
    400       }
    401     }
    402489#else
    403     if (status == 0){
     490    if (status == 0) {
    404491      //We want to exclude some cases
    405       switch (CPXgetstat(env, lp)){
     492      switch (CPXgetstat(cplexEnv(), _prob)) {
    406493      case CPX_OBJ_LIM:
    407494      case CPX_IT_LIM_FEAS:
     
    413500        return SOLVED;
    414501      }
    415     }
    416     else{
     502    } else {
    417503      return UNSOLVED;
    418504    }
     
    420506  }
    421507
    422   LpCplex::Value LpCplex::_getPrimal(int i) const
    423   {
     508  LpCplex::SolveExitStatus LpCplex::_solve() {
     509    _clear_temporals();
     510    return convertStatus(CPXlpopt(cplexEnv(), _prob));
     511  }
     512
     513  LpCplex::SolveExitStatus LpCplex::solvePrimal() {
     514    _clear_temporals();
     515    return convertStatus(CPXprimopt(cplexEnv(), _prob));
     516  }
     517
     518  LpCplex::SolveExitStatus LpCplex::solveDual() {
     519    _clear_temporals();
     520    return convertStatus(CPXdualopt(cplexEnv(), _prob));
     521  }
     522
     523  LpCplex::SolveExitStatus LpCplex::solveBarrier() {
     524    _clear_temporals();
     525    return convertStatus(CPXbaropt(cplexEnv(), _prob));
     526  }
     527
     528  LpCplex::Value LpCplex::_getPrimal(int i) const {
    424529    Value x;
    425     CPXgetx(env, lp, &x, i, i);
     530    CPXgetx(cplexEnv(), _prob, &x, i, i);
    426531    return x;
    427532  }
    428533
    429   LpCplex::Value LpCplex::_getDual(int i) const
    430   {
     534  LpCplex::Value LpCplex::_getDual(int i) const {
    431535    Value y;
    432     CPXgetpi(env, lp, &y, i, i);
     536    CPXgetpi(cplexEnv(), _prob, &y, i, i);
    433537    return y;
    434538  }
    435539
    436   LpCplex::Value LpCplex::_getPrimalValue() const
    437   {
     540  LpCplex::Value LpCplex::_getPrimalValue() const {
    438541    Value objval;
    439     //method = CPXgetmethod (env, lp);
    440     //printf("CPXgetprobtype %d \n",CPXgetprobtype(env,lp));
    441     CPXgetobjval(env, lp, &objval);
    442     //printf("Objective value: %g \n",objval);
     542    CPXgetobjval(cplexEnv(), _prob, &objval);
    443543    return objval;
    444544  }
    445   bool LpCplex::_isBasicCol(int i) const
    446   {
    447     std::vector<int> cstat(CPXgetnumcols(env, lp));
    448     CPXgetbase(env, lp, &*cstat.begin(), NULL);
    449     return (cstat[i]==CPX_BASIC);
    450   }
    451 
    452 //7.5-os cplex statusai (Vigyazat: a 9.0-asei masok!)
    453 // This table lists the statuses, returned by the CPXgetstat()
    454 // routine, for solutions to LP problems or mixed integer problems. If
    455 // no solution exists, the return value is zero.
    456 
    457 // For Simplex, Barrier
    458 // 1          CPX_OPTIMAL
    459 //          Optimal solution found
    460 // 2          CPX_INFEASIBLE
    461 //          Problem infeasible
    462 // 3    CPX_UNBOUNDED
    463 //          Problem unbounded
    464 // 4          CPX_OBJ_LIM
    465 //          Objective limit exceeded in Phase II
    466 // 5          CPX_IT_LIM_FEAS
    467 //          Iteration limit exceeded in Phase II
    468 // 6          CPX_IT_LIM_INFEAS
    469 //          Iteration limit exceeded in Phase I
    470 // 7          CPX_TIME_LIM_FEAS
    471 //          Time limit exceeded in Phase II
    472 // 8          CPX_TIME_LIM_INFEAS
    473 //          Time limit exceeded in Phase I
    474 // 9          CPX_NUM_BEST_FEAS
    475 //          Problem non-optimal, singularities in Phase II
    476 // 10         CPX_NUM_BEST_INFEAS
    477 //          Problem non-optimal, singularities in Phase I
    478 // 11         CPX_OPTIMAL_INFEAS
    479 //          Optimal solution found, unscaled infeasibilities
    480 // 12         CPX_ABORT_FEAS
    481 //          Aborted in Phase II
    482 // 13         CPX_ABORT_INFEAS
    483 //          Aborted in Phase I
    484 // 14          CPX_ABORT_DUAL_INFEAS
    485 //          Aborted in barrier, dual infeasible
    486 // 15          CPX_ABORT_PRIM_INFEAS
    487 //          Aborted in barrier, primal infeasible
    488 // 16          CPX_ABORT_PRIM_DUAL_INFEAS
    489 //          Aborted in barrier, primal and dual infeasible
    490 // 17          CPX_ABORT_PRIM_DUAL_FEAS
    491 //          Aborted in barrier, primal and dual feasible
    492 // 18          CPX_ABORT_CROSSOVER
    493 //          Aborted in crossover
    494 // 19          CPX_INForUNBD
    495 //          Infeasible or unbounded
    496 // 20   CPX_PIVOT
    497 //       User pivot used
    498 //
    499 //     Ezeket hova tegyem:
    500 // ??case CPX_ABORT_DUAL_INFEAS
    501 // ??case CPX_ABORT_CROSSOVER
    502 // ??case CPX_INForUNBD
    503 // ??case CPX_PIVOT
    504 
    505 //Some more interesting stuff:
    506 
    507 // CPX_PARAM_LPMETHOD  1062  int  LPMETHOD
    508 // 0 Automatic
    509 // 1 Primal Simplex
    510 // 2 Dual Simplex
    511 // 3 Network Simplex
    512 // 4 Standard Barrier
    513 // Default: 0
    514 // Description: Method for linear optimization.
    515 // Determines which algorithm is used when CPXlpopt() (or "optimize"
    516 // in the Interactive Optimizer) is called. Currently the behavior of
    517 // the "Automatic" setting is that CPLEX simply invokes the dual
    518 // simplex method, but this capability may be expanded in the future
    519 // so that CPLEX chooses the method based on problem characteristics
     545
     546  LpCplex::VarStatus LpCplex::_getColStatus(int i) const {
     547    if (_col_status.empty()) {
     548      _col_status.resize(CPXgetnumcols(cplexEnv(), _prob));
     549      CPXgetbase(cplexEnv(), _prob, &_col_status.front(), 0);
     550    }
     551    switch (_col_status[i]) {
     552    case CPX_BASIC:
     553      return BASIC;
     554    case CPX_FREE_SUPER:
     555      return FREE;
     556    case CPX_AT_LOWER:
     557      return LOWER;
     558    case CPX_AT_UPPER:
     559      return UPPER;
     560    default:
     561      LEMON_ASSERT(false, "Wrong column status");
     562      return LpCplex::VarStatus();
     563    }
     564  }
     565
     566  LpCplex::VarStatus LpCplex::_getRowStatus(int i) const {
     567    if (_row_status.empty()) {
     568      _row_status.resize(CPXgetnumrows(cplexEnv(), _prob));
     569      CPXgetbase(cplexEnv(), _prob, 0, &_row_status.front());
     570    }
     571    switch (_row_status[i]) {
     572    case CPX_BASIC:
     573      return BASIC;
     574    case CPX_AT_LOWER:
     575      {
     576        char s;
     577        CPXgetsense(cplexEnv(), _prob, &s, i, i);
     578        return s != 'L' ? LOWER : UPPER;
     579      }
     580    case CPX_AT_UPPER:
     581      return UPPER;
     582    default:
     583      LEMON_ASSERT(false, "Wrong row status");
     584      return LpCplex::VarStatus();
     585    }
     586  }
     587
     588  LpCplex::Value LpCplex::_getPrimalRay(int i) const {
     589    if (_primal_ray.empty()) {
     590      _primal_ray.resize(CPXgetnumcols(cplexEnv(), _prob));
     591      CPXgetray(cplexEnv(), _prob, &_primal_ray.front());
     592    }
     593    return _primal_ray[i];
     594  }
     595
     596  LpCplex::Value LpCplex::_getDualRay(int i) const {
     597    if (_dual_ray.empty()) {
     598
     599    }
     600    return _dual_ray[i];
     601  }
     602
     603  //7.5-os cplex statusai (Vigyazat: a 9.0-asei masok!)
     604  // This table lists the statuses, returned by the CPXgetstat()
     605  // routine, for solutions to LP problems or mixed integer problems. If
     606  // no solution exists, the return value is zero.
     607
     608  // For Simplex, Barrier
     609  // 1          CPX_OPTIMAL
     610  //          Optimal solution found
     611  // 2          CPX_INFEASIBLE
     612  //          Problem infeasible
     613  // 3    CPX_UNBOUNDED
     614  //          Problem unbounded
     615  // 4          CPX_OBJ_LIM
     616  //          Objective limit exceeded in Phase II
     617  // 5          CPX_IT_LIM_FEAS
     618  //          Iteration limit exceeded in Phase II
     619  // 6          CPX_IT_LIM_INFEAS
     620  //          Iteration limit exceeded in Phase I
     621  // 7          CPX_TIME_LIM_FEAS
     622  //          Time limit exceeded in Phase II
     623  // 8          CPX_TIME_LIM_INFEAS
     624  //          Time limit exceeded in Phase I
     625  // 9          CPX_NUM_BEST_FEAS
     626  //          Problem non-optimal, singularities in Phase II
     627  // 10         CPX_NUM_BEST_INFEAS
     628  //          Problem non-optimal, singularities in Phase I
     629  // 11         CPX_OPTIMAL_INFEAS
     630  //          Optimal solution found, unscaled infeasibilities
     631  // 12         CPX_ABORT_FEAS
     632  //          Aborted in Phase II
     633  // 13         CPX_ABORT_INFEAS
     634  //          Aborted in Phase I
     635  // 14          CPX_ABORT_DUAL_INFEAS
     636  //          Aborted in barrier, dual infeasible
     637  // 15          CPX_ABORT_PRIM_INFEAS
     638  //          Aborted in barrier, primal infeasible
     639  // 16          CPX_ABORT_PRIM_DUAL_INFEAS
     640  //          Aborted in barrier, primal and dual infeasible
     641  // 17          CPX_ABORT_PRIM_DUAL_FEAS
     642  //          Aborted in barrier, primal and dual feasible
     643  // 18          CPX_ABORT_CROSSOVER
     644  //          Aborted in crossover
     645  // 19          CPX_INForUNBD
     646  //          Infeasible or unbounded
     647  // 20   CPX_PIVOT
     648  //       User pivot used
     649  //
     650  //     Ezeket hova tegyem:
     651  // ??case CPX_ABORT_DUAL_INFEAS
     652  // ??case CPX_ABORT_CROSSOVER
     653  // ??case CPX_INForUNBD
     654  // ??case CPX_PIVOT
     655
     656  //Some more interesting stuff:
     657
     658  // CPX_PARAM_PROBMETHOD  1062  int  LPMETHOD
     659  // 0 Automatic
     660  // 1 Primal Simplex
     661  // 2 Dual Simplex
     662  // 3 Network Simplex
     663  // 4 Standard Barrier
     664  // Default: 0
     665  // Description: Method for linear optimization.
     666  // Determines which algorithm is used when CPXlpopt() (or "optimize"
     667  // in the Interactive Optimizer) is called. Currently the behavior of
     668  // the "Automatic" setting is that CPLEX simply invokes the dual
     669  // simplex method, but this capability may be expanded in the future
     670  // so that CPLEX chooses the method based on problem characteristics
    520671#if CPX_VERSION < 900
    521   void statusSwitch(CPXENVptr env,int& stat){
     672  void statusSwitch(CPXENVptr cplexEnv(),int& stat){
    522673    int lpmethod;
    523     CPXgetintparam (env,CPX_PARAM_LPMETHOD,&lpmethod);
     674    CPXgetintparam (cplexEnv(),CPX_PARAM_PROBMETHOD,&lpmethod);
    524675    if (lpmethod==2){
    525676      if (stat==CPX_UNBOUNDED){
     
    536687#endif
    537688
    538   LpCplex::SolutionStatus LpCplex::_getPrimalStatus() const
    539   {
    540     //Unboundedness not treated well: the following is from cplex 9.0 doc
     689  LpCplex::ProblemType LpCplex::_getPrimalType() const {
     690    // Unboundedness not treated well: the following is from cplex 9.0 doc
    541691    // About Unboundedness
    542692
     
    553703    // has a feasible solution.
    554704
    555     int stat = CPXgetstat(env, lp);
     705    int stat = CPXgetstat(cplexEnv(), _prob);
    556706#if CPX_VERSION >= 800
    557707    switch (stat)
    558     {
     708      {
    559709      case CPX_STAT_OPTIMAL:
    560710        return OPTIMAL;
    561711      case CPX_STAT_UNBOUNDED:
    562         return INFINITE;
     712        return UNBOUNDED;
    563713      case CPX_STAT_INFEASIBLE:
    564714        return INFEASIBLE;
    565715      default:
    566716        return UNDEFINED;
    567     }
     717      }
    568718#else
    569     statusSwitch(env,stat);
    570     //CPXgetstat(env, lp);
     719    statusSwitch(cplexEnv(),stat);
     720    //CPXgetstat(cplexEnv(), _prob);
    571721    //printf("A primal status: %d, CPX_OPTIMAL=%d \n",stat,CPX_OPTIMAL);
    572722    switch (stat) {
     
    577727    case CPX_UNBOUNDED://Unbounded
    578728      return INFEASIBLE;//In case of dual simplex
    579       //return INFINITE;
     729      //return UNBOUNDED;
    580730    case CPX_INFEASIBLE://Infeasible
    581  //    case CPX_IT_LIM_INFEAS:
    582 //     case CPX_TIME_LIM_INFEAS:
    583 //     case CPX_NUM_BEST_INFEAS:
    584 //     case CPX_OPTIMAL_INFEAS:
    585 //     case CPX_ABORT_INFEAS:
    586 //     case CPX_ABORT_PRIM_INFEAS:
    587 //     case CPX_ABORT_PRIM_DUAL_INFEAS:
    588       return INFINITE;//In case of dual simplex
     731      //    case CPX_IT_LIM_INFEAS:
     732      //     case CPX_TIME_LIM_INFEAS:
     733      //     case CPX_NUM_BEST_INFEAS:
     734      //     case CPX_OPTIMAL_INFEAS:
     735      //     case CPX_ABORT_INFEAS:
     736      //     case CPX_ABORT_PRIM_INFEAS:
     737      //     case CPX_ABORT_PRIM_DUAL_INFEAS:
     738      return UNBOUNDED;//In case of dual simplex
    589739      //return INFEASIBLE;
    590 //     case CPX_OBJ_LIM:
    591 //     case CPX_IT_LIM_FEAS:
    592 //     case CPX_TIME_LIM_FEAS:
    593 //     case CPX_NUM_BEST_FEAS:
    594 //     case CPX_ABORT_FEAS:
    595 //     case CPX_ABORT_PRIM_DUAL_FEAS:
    596 //       return FEASIBLE;
     740      //     case CPX_OBJ_LIM:
     741      //     case CPX_IT_LIM_FEAS:
     742      //     case CPX_TIME_LIM_FEAS:
     743      //     case CPX_NUM_BEST_FEAS:
     744      //     case CPX_ABORT_FEAS:
     745      //     case CPX_ABORT_PRIM_DUAL_FEAS:
     746      //       return FEASIBLE;
    597747    default:
    598748      return UNDEFINED; //Everything else comes here
     
    602752  }
    603753
    604 //9.0-as cplex verzio statusai
    605 // CPX_STAT_ABORT_DUAL_OBJ_LIM
    606 // CPX_STAT_ABORT_IT_LIM
    607 // CPX_STAT_ABORT_OBJ_LIM
    608 // CPX_STAT_ABORT_PRIM_OBJ_LIM
    609 // CPX_STAT_ABORT_TIME_LIM
    610 // CPX_STAT_ABORT_USER
    611 // CPX_STAT_FEASIBLE_RELAXED
    612 // CPX_STAT_INFEASIBLE
    613 // CPX_STAT_INForUNBD
    614 // CPX_STAT_NUM_BEST
    615 // CPX_STAT_OPTIMAL
    616 // CPX_STAT_OPTIMAL_FACE_UNBOUNDED
    617 // CPX_STAT_OPTIMAL_INFEAS
    618 // CPX_STAT_OPTIMAL_RELAXED
    619 // CPX_STAT_UNBOUNDED
    620 
    621   LpCplex::SolutionStatus LpCplex::_getDualStatus() const
    622   {
    623     int stat = CPXgetstat(env, lp);
     754  //9.0-as cplex verzio statusai
     755  // CPX_STAT_ABORT_DUAL_OBJ_LIM
     756  // CPX_STAT_ABORT_IT_LIM
     757  // CPX_STAT_ABORT_OBJ_LIM
     758  // CPX_STAT_ABORT_PRIM_OBJ_LIM
     759  // CPX_STAT_ABORT_TIME_LIM
     760  // CPX_STAT_ABORT_USER
     761  // CPX_STAT_FEASIBLE_RELAXED
     762  // CPX_STAT_INFEASIBLE
     763  // CPX_STAT_INForUNBD
     764  // CPX_STAT_NUM_BEST
     765  // CPX_STAT_OPTIMAL
     766  // CPX_STAT_OPTIMAL_FACE_UNBOUNDED
     767  // CPX_STAT_OPTIMAL_INFEAS
     768  // CPX_STAT_OPTIMAL_RELAXED
     769  // CPX_STAT_UNBOUNDED
     770
     771  LpCplex::ProblemType LpCplex::_getDualType() const {
     772    int stat = CPXgetstat(cplexEnv(), _prob);
    624773#if CPX_VERSION >= 800
    625     switch (stat)
    626     {
    627       case CPX_STAT_OPTIMAL:
    628         return OPTIMAL;
    629       case CPX_STAT_UNBOUNDED:
    630         return INFEASIBLE;
    631       default:
    632         return UNDEFINED;
     774    switch (stat) {
     775    case CPX_STAT_OPTIMAL:
     776      return OPTIMAL;
     777    case CPX_STAT_UNBOUNDED:
     778      return INFEASIBLE;
     779    default:
     780      return UNDEFINED;
    633781    }
    634782#else
    635     statusSwitch(env,stat);
     783    statusSwitch(cplexEnv(),stat);
    636784    switch (stat) {
    637785    case 0:
     
    640788      return OPTIMAL;
    641789    case CPX_UNBOUNDED:
    642      return INFEASIBLE;
     790      return INFEASIBLE;
    643791    default:
    644792      return UNDEFINED; //Everything else comes here
     
    648796  }
    649797
    650   LpCplex::ProblemTypes LpCplex::_getProblemType() const
    651   {
    652     int stat = CPXgetstat(env, lp);
    653 #if CPX_VERSION >= 800
    654     switch (stat)
    655     {
    656       case CPX_STAT_OPTIMAL:
    657         return PRIMAL_DUAL_FEASIBLE;
    658       case CPX_STAT_UNBOUNDED:
    659          return PRIMAL_FEASIBLE_DUAL_INFEASIBLE;
    660       default:
    661         return UNKNOWN;
    662     }
     798  // MipCplex members
     799
     800  MipCplex::MipCplex()
     801    : LpBase(), CplexBase(), MipSolver() {
     802
     803#if CPX_VERSION < 800
     804    CPXchgprobtype(cplexEnv(),  _prob, CPXPROB_MIP);
    663805#else
     806    CPXchgprobtype(cplexEnv(),  _prob, CPXPROB_MILP);
     807#endif
     808  }
     809
     810  MipCplex::MipCplex(const CplexEnv& env)
     811    : LpBase(), CplexBase(env), MipSolver() {
     812
     813#if CPX_VERSION < 800
     814    CPXchgprobtype(cplexEnv(),  _prob, CPXPROB_MIP);
     815#else
     816    CPXchgprobtype(cplexEnv(),  _prob, CPXPROB_MILP);
     817#endif
     818
     819  }
     820
     821  MipCplex::MipCplex(const MipCplex& other)
     822    : LpBase(), CplexBase(other), MipSolver() {}
     823
     824  MipCplex::~MipCplex() {}
     825
     826  MipCplex* MipCplex::_newSolver() const { return new MipCplex; }
     827  MipCplex* MipCplex::_cloneSolver() const {return new MipCplex(*this); }
     828
     829  const char* MipCplex::_solverName() const { return "MipCplex"; }
     830
     831  void MipCplex::_setColType(int i, MipCplex::ColTypes col_type) {
     832
     833    // Note If a variable is to be changed to binary, a call to CPXchgbds
     834    // should also be made to change the bounds to 0 and 1.
     835
     836    switch (col_type){
     837    case INTEGER: {
     838      const char t = 'I';
     839      CPXchgctype (cplexEnv(), _prob, 1, &i, &t);
     840    } break;
     841    case REAL: {
     842      const char t = 'C';
     843      CPXchgctype (cplexEnv(), _prob, 1, &i, &t);
     844    } break;
     845    default:
     846      break;
     847    }
     848  }
     849
     850  MipCplex::ColTypes MipCplex::_getColType(int i) const {
     851    char t;
     852    CPXgetctype (cplexEnv(), _prob, &t, i, i);
     853    switch (t) {
     854    case 'I':
     855      return INTEGER;
     856    case 'C':
     857      return REAL;
     858    default:
     859      LEMON_ASSERT(false, "Invalid column type");
     860      return ColTypes();
     861    }
     862
     863  }
     864
     865  MipCplex::SolveExitStatus MipCplex::_solve() {
     866    int status;
     867    status = CPXmipopt (cplexEnv(), _prob);
     868    if (status==0)
     869      return SOLVED;
     870    else
     871      return UNSOLVED;
     872
     873  }
     874
     875
     876  MipCplex::ProblemType MipCplex::_getType() const {
     877
     878    int stat = CPXgetstat(cplexEnv(), _prob);
     879
     880    //Fortunately, MIP statuses did not change for cplex 8.0
    664881    switch (stat) {
    665     case CPX_OPTIMAL://Optimal
    666         return PRIMAL_DUAL_FEASIBLE;
    667     case CPX_UNBOUNDED:
    668          return PRIMAL_FEASIBLE_DUAL_INFEASIBLE;
    669 //         return PRIMAL_INFEASIBLE_DUAL_FEASIBLE;
    670 //         return PRIMAL_DUAL_INFEASIBLE;
    671 
    672 //Seems to be that this is all we can say for sure
    673     default:
    674         //In all other cases
    675         return UNKNOWN;
    676       //FIXME error
    677     }
    678 #endif
    679   }
    680 
    681   void LpCplex::_setMax()
    682   {
    683     CPXchgobjsen(env, lp, CPX_MAX);
    684    }
    685   void LpCplex::_setMin()
    686   {
    687     CPXchgobjsen(env, lp, CPX_MIN);
    688    }
    689 
    690   bool LpCplex::_isMax() const
    691   {
    692     if (CPXgetobjsen(env, lp)==CPX_MAX)
    693       return true;
    694     else
    695       return false;
     882    case CPXMIP_OPTIMAL:
     883      // Optimal integer solution has been found.
     884    case CPXMIP_OPTIMAL_TOL:
     885      // Optimal soluton with the tolerance defined by epgap or epagap has
     886      // been found.
     887      return OPTIMAL;
     888      //This also exists in later issues
     889      //    case CPXMIP_UNBOUNDED:
     890      //return UNBOUNDED;
     891      case CPXMIP_INFEASIBLE:
     892        return INFEASIBLE;
     893    default:
     894      return UNDEFINED;
     895    }
     896    //Unboundedness not treated well: the following is from cplex 9.0 doc
     897    // About Unboundedness
     898
     899    // The treatment of models that are unbounded involves a few
     900    // subtleties. Specifically, a declaration of unboundedness means that
     901    // ILOG CPLEX has determined that the model has an unbounded
     902    // ray. Given any feasible solution x with objective z, a multiple of
     903    // the unbounded ray can be added to x to give a feasible solution
     904    // with objective z-1 (or z+1 for maximization models). Thus, if a
     905    // feasible solution exists, then the optimal objective is
     906    // unbounded. Note that ILOG CPLEX has not necessarily concluded that
     907    // a feasible solution exists. Users can call the routine CPXsolninfo
     908    // to determine whether ILOG CPLEX has also concluded that the model
     909    // has a feasible solution.
     910  }
     911
     912  MipCplex::Value MipCplex::_getSol(int i) const {
     913    Value x;
     914    CPXgetmipx(cplexEnv(), _prob, &x, i, i);
     915    return x;
     916  }
     917
     918  MipCplex::Value MipCplex::_getSolValue() const {
     919    Value objval;
     920    CPXgetmipobjval(cplexEnv(), _prob, &objval);
     921    return objval;
    696922  }
    697923
  • lemon/lp_cplex.h

    r481 r482  
    3030namespace lemon {
    3131
    32 
    33   /// \brief Interface for the CPLEX solver
    34   ///
    35   /// This class implements an interface for the CPLEX LP solver.
    36   class LpCplex :virtual public LpSolverBase {
    37 
    38   public:
    39 
    40     typedef LpSolverBase Parent;
    41 
    42     /// \e
    43     int status;
    44     cpxenv* env;
    45     cpxlp* lp;
    46 
    47 
    48     /// \e
    49     LpCplex();
    50     /// \e
    51     LpCplex(const LpCplex&);
    52     /// \e
    53     ~LpCplex();
    54 
    55   protected:
    56     virtual LpSolverBase* _newLp();
    57     virtual LpSolverBase* _copyLp();
    58 
     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 CplexBase;
     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
     81  /// MIP solvers. 
     82  /// \ingroup lp_group
     83  class CplexBase : virtual public LpBase {
     84  protected:
     85
     86    CplexEnv _env;
     87    cpxlp* _prob;
     88
     89    CplexBase();
     90    CplexBase(const CplexEnv&);
     91    CplexBase(const CplexBase &);
     92    virtual ~CplexBase();
    5993
    6094    virtual int _addCol();
    6195    virtual int _addRow();
     96
    6297    virtual void _eraseCol(int i);
    6398    virtual void _eraseRow(int i);
    64     virtual void _getColName(int col, std::string & name) const;
    65     virtual void _setColName(int col, const std::string & name);
     99
     100    virtual void _eraseColId(int i);
     101    virtual void _eraseRowId(int i);
     102
     103    virtual void _getColName(int col, std::string& name) const;
     104    virtual void _setColName(int col, const std::string& name);
    66105    virtual int _colByName(const std::string& name) const;
    67     virtual void _setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e);
    68     virtual void _getRowCoeffs(int i, RowIterator b) const;
    69     virtual void _setColCoeffs(int i, ConstColIterator b, ConstColIterator e);
    70     virtual void _getColCoeffs(int i, ColIterator b) const;
     106
     107    virtual void _getRowName(int row, std::string& name) const;
     108    virtual void _setRowName(int row, const std::string& name);
     109    virtual int _rowByName(const std::string& name) const;
     110
     111    virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e);
     112    virtual void _getRowCoeffs(int i, InsertIterator b) const;
     113
     114    virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e);
     115    virtual void _getColCoeffs(int i, InsertIterator b) const;
     116
    71117    virtual void _setCoeff(int row, int col, Value value);
    72118    virtual Value _getCoeff(int row, int col) const;
     
    74120    virtual void _setColLowerBound(int i, Value value);
    75121    virtual Value _getColLowerBound(int i) const;
     122
    76123    virtual void _setColUpperBound(int i, Value value);
    77124    virtual Value _getColUpperBound(int i) const;
    78125
    79 //     virtual void _setRowLowerBound(int i, Value value);
    80 //     virtual void _setRowUpperBound(int i, Value value);
    81     virtual void _setRowBounds(int i, Value lower, Value upper);
    82     virtual void _getRowBounds(int i, Value &lb, Value &ub) const;
     126  private:
     127    void _set_row_bounds(int i, Value lb, Value ub);
     128  protected:
     129
     130    virtual void _setRowLowerBound(int i, Value value);
     131    virtual Value _getRowLowerBound(int i) const;
     132
     133    virtual void _setRowUpperBound(int i, Value value);
     134    virtual Value _getRowUpperBound(int i) const;
     135
     136    virtual void _setObjCoeffs(ExprIterator b, ExprIterator e);
     137    virtual void _getObjCoeffs(InsertIterator b) const;
     138
    83139    virtual void _setObjCoeff(int i, Value obj_coef);
    84140    virtual Value _getObjCoeff(int i) const;
    85     virtual void _clearObj();
    86 
     141
     142    virtual void _setSense(Sense sense);
     143    virtual Sense _getSense() const;
     144
     145    virtual void _clear();
     146
     147  public:
     148
     149    /// Returns the used \c CplexEnv instance
     150    const CplexEnv& env() const { return _env; }
     151    ///
     152    const cpxenv* cplexEnv() const { return _env.cplexEnv(); }
     153
     154    cpxlp* cplexLp() { return _prob; }
     155    const cpxlp* cplexLp() const { return _prob; }
     156
     157  };
     158
     159  /// \brief Interface for the CPLEX LP solver
     160  ///
     161  /// This class implements an interface for the CPLEX LP solver.
     162  ///\ingroup lp_group
     163  class LpCplex : public CplexBase, public LpSolver {
     164  public:
     165    /// \e
     166    LpCplex();
     167    /// \e
     168    LpCplex(const CplexEnv&);
     169    /// \e
     170    LpCplex(const LpCplex&);
     171    /// \e
     172    virtual ~LpCplex();
     173
     174  private:
     175
     176    // these values cannot retrieved element by element
     177    mutable std::vector<int> _col_status;
     178    mutable std::vector<int> _row_status;
     179
     180    mutable std::vector<Value> _primal_ray;
     181    mutable std::vector<Value> _dual_ray;
     182
     183    void _clear_temporals();
     184
     185    SolveExitStatus convertStatus(int status);
     186
     187  protected:
     188
     189    virtual LpCplex* _cloneSolver() const;
     190    virtual LpCplex* _newSolver() const;
     191
     192    virtual const char* _solverName() const;
    87193
    88194    virtual SolveExitStatus _solve();
     
    90196    virtual Value _getDual(int i) const;
    91197    virtual Value _getPrimalValue() const;
    92     virtual bool _isBasicCol(int i) const;
    93 
    94     virtual SolutionStatus _getPrimalStatus() const;
    95     virtual SolutionStatus _getDualStatus() const;
    96     virtual ProblemTypes _getProblemType() const;
    97 
    98 
    99     virtual void _setMax();
    100     virtual void _setMin();
    101 
    102     virtual bool _isMax() const;
    103 
    104   public:
    105 
    106     cpxenv* cplexEnv() { return env; }
    107     cpxlp* cplexLp() { return lp; }
    108 
    109   };
     198
     199    virtual VarStatus _getColStatus(int i) const;
     200    virtual VarStatus _getRowStatus(int i) const;
     201
     202    virtual Value _getPrimalRay(int i) const;
     203    virtual Value _getDualRay(int i) const;
     204
     205    virtual ProblemType _getPrimalType() const;
     206    virtual ProblemType _getDualType() const;
     207
     208  public:
     209
     210    /// Solve with primal simplex method
     211    SolveExitStatus solvePrimal();
     212
     213    /// Solve with dual simplex method
     214    SolveExitStatus solveDual();
     215
     216    /// Solve with barrier method
     217    SolveExitStatus solveBarrier();
     218
     219  };
     220
     221  /// \brief Interface for the CPLEX MIP solver
     222  ///
     223  /// This class implements an interface for the CPLEX MIP solver.
     224  ///\ingroup lp_group
     225  class MipCplex : public CplexBase, public MipSolver {
     226  public:
     227    /// \e
     228    MipCplex();
     229    /// \e
     230    MipCplex(const CplexEnv&);
     231    /// \e
     232    MipCplex(const MipCplex&);
     233    /// \e
     234    virtual ~MipCplex();
     235
     236  protected:
     237
     238    virtual MipCplex* _cloneSolver() const;
     239    virtual MipCplex* _newSolver() const;
     240
     241    virtual const char* _solverName() const;
     242
     243    virtual ColTypes _getColType(int col) const;
     244    virtual void _setColType(int col, ColTypes col_type);
     245
     246    virtual SolveExitStatus _solve();
     247    virtual ProblemType _getType() const;
     248    virtual Value _getSol(int i) const;
     249    virtual Value _getSolValue() const;
     250
     251  };
     252
    110253} //END OF NAMESPACE LEMON
    111254
  • lemon/lp_glpk.cc

    r481 r482  
    1818
    1919///\file
    20 ///\brief Implementation of the LEMON-GLPK lp solver interface.
     20///\brief Implementation of the LEMON GLPK LP and MIP solver interface.
    2121
    2222#include <lemon/lp_glpk.h>
    23 //#include <iostream>
    24 
    25 extern "C" {
    2623#include <glpk.h>
    27 }
    28 
    29 #if GLP_MAJOR_VERSION > 4 || (GLP_MAJOR_VERSION == 4 && GLP_MINOR_VERSION > 15)
    30 #define LEMON_glp(func) (glp_##func)
    31 #define LEMON_lpx(func) (lpx_##func)
    32 
    33 #define LEMON_GLP(def) (GLP_##def)
    34 #define LEMON_LPX(def) (LPX_##def)
    35 
    36 #else
    37 
    38 #define LEMON_glp(func) (lpx_##func)
    39 #define LEMON_lpx(func) (lpx_##func)
    40 
    41 #define LEMON_GLP(def) (LPX_##def)
    42 #define LEMON_LPX(def) (LPX_##def)
    43 
    44 #endif
     24
     25#include <lemon/assert.h>
    4526
    4627namespace lemon {
    4728
    48   LpGlpk::LpGlpk() : Parent() {
    49     solved = false;
    50     rows = _lp_bits::LpId(1);
    51     cols = _lp_bits::LpId(1);
    52     lp = LEMON_glp(create_prob)();
    53     LEMON_glp(create_index)(lp);
    54     messageLevel(0);
    55   }
    56 
    57   LpGlpk::LpGlpk(const LpGlpk &glp) : Parent() {
    58     solved = false;
    59     rows = _lp_bits::LpId(1);
    60     cols = _lp_bits::LpId(1);
    61     lp = LEMON_glp(create_prob)();
    62     LEMON_glp(create_index)(lp);
    63     messageLevel(0);
    64     //Coefficient matrix, row bounds
    65     LEMON_glp(add_rows)(lp, LEMON_glp(get_num_rows)(glp.lp));
    66     LEMON_glp(add_cols)(lp, LEMON_glp(get_num_cols)(glp.lp));
    67     int len;
    68     std::vector<int> ind(1+LEMON_glp(get_num_cols)(glp.lp));
    69     std::vector<Value> val(1+LEMON_glp(get_num_cols)(glp.lp));
    70     for (int i=1;i<=LEMON_glp(get_num_rows)(glp.lp);++i)
    71       {
    72         len=LEMON_glp(get_mat_row)(glp.lp,i,&*ind.begin(),&*val.begin());
    73         LEMON_glp(set_mat_row)(lp, i,len,&*ind.begin(),&*val.begin());
    74         LEMON_glp(set_row_bnds)(lp,i,
    75                                 LEMON_glp(get_row_type)(glp.lp,i),
    76                                 LEMON_glp(get_row_lb)(glp.lp,i),
    77                                 LEMON_glp(get_row_ub)(glp.lp,i));
    78       }
    79 
    80     //Objective function, coloumn bounds
    81     LEMON_glp(set_obj_dir)(lp, LEMON_glp(get_obj_dir)(glp.lp));
    82     //Objectif function's constant term treated separately
    83     LEMON_glp(set_obj_coef)(lp,0,LEMON_glp(get_obj_coef)(glp.lp,0));
    84     for (int i=1;i<=LEMON_glp(get_num_cols)(glp.lp);++i)
    85       {
    86         LEMON_glp(set_obj_coef)(lp,i,
    87                                 LEMON_glp(get_obj_coef)(glp.lp,i));
    88         LEMON_glp(set_col_bnds)(lp,i,
    89                                 LEMON_glp(get_col_type)(glp.lp,i),
    90                                 LEMON_glp(get_col_lb)(glp.lp,i),
    91                                 LEMON_glp(get_col_ub)(glp.lp,i));
    92       }
    93     rows = glp.rows;
    94     cols = glp.cols;
    95   }
    96 
    97   LpGlpk::~LpGlpk() {
    98     LEMON_glp(delete_prob)(lp);
    99   }
    100 
    101   int LpGlpk::_addCol() {
    102     int i=LEMON_glp(add_cols)(lp, 1);
    103     LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(FR), 0.0, 0.0);
    104     solved = false;
     29  // GlpkBase members
     30
     31  GlpkBase::GlpkBase() : LpBase() {
     32    lp = glp_create_prob();
     33    glp_create_index(lp);
     34  }
     35
     36  GlpkBase::GlpkBase(const GlpkBase &other) : LpBase() {
     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  GlpkBase::~GlpkBase() {
     45    glp_delete_prob(lp);
     46  }
     47
     48  int GlpkBase::_addCol() {
     49    int i = glp_add_cols(lp, 1);
     50    glp_set_col_bnds(lp, i, GLP_FR, 0.0, 0.0);
    10551    return i;
    10652  }
    10753
    108   ///\e
    109 
    110 
    111   LpSolverBase* LpGlpk::_newLp()
    112   {
    113     LpGlpk* newlp = new LpGlpk;
    114     return newlp;
    115   }
    116 
    117   ///\e
    118 
    119   LpSolverBase* LpGlpk::_copyLp()
    120   {
    121     LpGlpk *newlp = new LpGlpk(*this);
    122     return newlp;
    123   }
    124 
    125   int LpGlpk::_addRow() {
    126     int i=LEMON_glp(add_rows)(lp, 1);
    127     solved = false;
     54  int GlpkBase::_addRow() {
     55    int i = glp_add_rows(lp, 1);
     56    glp_set_row_bnds(lp, i, GLP_FR, 0.0, 0.0);
    12857    return i;
    12958  }
    13059
    131 
    132   void LpGlpk::_eraseCol(int i) {
     60  void GlpkBase::_eraseCol(int i) {
    13361    int ca[2];
    134     ca[1]=i;
    135     LEMON_glp(del_cols)(lp, 1, ca);
    136     solved = false;
    137   }
    138 
    139   void LpGlpk::_eraseRow(int i) {
     62    ca[1] = i;
     63    glp_del_cols(lp, 1, ca);
     64  }
     65
     66  void GlpkBase::_eraseRow(int i) {
    14067    int ra[2];
    141     ra[1]=i;
    142     LEMON_glp(del_rows)(lp, 1, ra);
    143     solved = false;
    144   }
    145 
    146   void LpGlpk::_getColName(int c, std::string & name) const
    147   {
    148 
    149     const char *n = LEMON_glp(get_col_name)(lp,c);
    150     name = n?n:"";
    151   }
    152 
    153 
    154   void LpGlpk::_setColName(int c, const std::string & name)
    155   {
    156     LEMON_glp(set_col_name)(lp,c,const_cast<char*>(name.c_str()));
    157 
    158   }
    159 
    160   int LpGlpk::_colByName(const std::string& name) const
    161   {
    162     int k = LEMON_glp(find_col)(lp, const_cast<char*>(name.c_str()));
     68    ra[1] = i;
     69    glp_del_rows(lp, 1, ra);
     70  }
     71
     72  void GlpkBase::_eraseColId(int i) {
     73    cols.eraseIndex(i);
     74    cols.shiftIndices(i);
     75  }
     76
     77  void GlpkBase::_eraseRowId(int i) {
     78    rows.eraseIndex(i);
     79    rows.shiftIndices(i);
     80  }
     81
     82  void GlpkBase::_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 GlpkBase::_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 GlpkBase::_colByName(const std::string& name) const {
     94    int k = glp_find_col(lp, const_cast<char*>(name.c_str()));
    16395    return k > 0 ? k : -1;
    16496  }
    16597
    166 
    167   void LpGlpk::_setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e)
    168   {
    169     std::vector<int> indices;
     98  void GlpkBase::_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 GlpkBase::_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 GlpkBase::_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 GlpkBase::_setRowCoeffs(int i, ExprIterator b, ExprIterator e) {
     115    std::vector<int> indexes;
    170116    std::vector<Value> values;
    171117
    172     indices.push_back(0);
     118    indexes.push_back(0);
    173119    values.push_back(0);
    174120
    175     for(ConstRowIterator it=b; it!=e; ++it) {
    176       indices.push_back(it->first);
     121    for(ExprIterator it = b; it != e; ++it) {
     122      indexes.push_back(it->first);
    177123      values.push_back(it->second);
    178124    }
    179125
    180     LEMON_glp(set_mat_row)(lp, i, values.size() - 1,
    181                                 &indices[0], &values[0]);
    182 
    183     solved = false;
    184   }
    185 
    186   void LpGlpk::_getRowCoeffs(int ix, RowIterator b) const
    187   {
    188     int length = LEMON_glp(get_mat_row)(lp, ix, 0, 0);
    189 
    190     std::vector<int> indices(length + 1);
     126    glp_set_mat_row(lp, i, values.size() - 1,
     127                    &indexes.front(), &values.front());
     128  }
     129
     130  void GlpkBase::_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);
    191134    std::vector<Value> values(length + 1);
    192135
    193     LEMON_glp(get_mat_row)(lp, ix, &indices[0], &values[0]);
     136    glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
    194137
    195138    for (int i = 1; i <= length; ++i) {
    196       *b = std::make_pair(indices[i], values[i]);
     139      *b = std::make_pair(indexes[i], values[i]);
    197140      ++b;
    198141    }
    199142  }
    200143
    201   void LpGlpk::_setColCoeffs(int ix, ConstColIterator b, ConstColIterator e) {
    202 
    203     std::vector<int> indices;
     144  void GlpkBase::_setColCoeffs(int ix, ExprIterator b,
     145                                     ExprIterator e) {
     146
     147    std::vector<int> indexes;
    204148    std::vector<Value> values;
    205149
    206     indices.push_back(0);
     150    indexes.push_back(0);
    207151    values.push_back(0);
    208152
    209     for(ConstColIterator it=b; it!=e; ++it) {
    210       indices.push_back(it->first);
     153    for(ExprIterator it = b; it != e; ++it) {
     154      indexes.push_back(it->first);
    211155      values.push_back(it->second);
    212156    }
    213157
    214     LEMON_glp(set_mat_col)(lp, ix, values.size() - 1,
    215                                 &indices[0], &values[0]);
    216 
    217     solved = false;
    218   }
    219 
    220   void LpGlpk::_getColCoeffs(int ix, ColIterator b) const
    221   {
    222     int length = LEMON_glp(get_mat_col)(lp, ix, 0, 0);
    223 
    224     std::vector<int> indices(length + 1);
     158    glp_set_mat_col(lp, ix, values.size() - 1,
     159                    &indexes.front(), &values.front());
     160  }
     161
     162  void GlpkBase::_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);
    225166    std::vector<Value> values(length + 1);
    226167
    227     LEMON_glp(get_mat_col)(lp, ix, &indices[0], &values[0]);
    228 
    229     for (int i = 1; i <= length; ++i) {
    230       *b = std::make_pair(indices[i], values[i]);
     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]);
    231172      ++b;
    232173    }
    233174  }
    234175
    235   void LpGlpk::_setCoeff(int ix, int jx, Value value)
    236   {
    237 
    238     if (LEMON_glp(get_num_cols)(lp) < LEMON_glp(get_num_rows)(lp)) {
    239 
    240       int length=LEMON_glp(get_mat_row)(lp, ix, 0, 0);
    241 
    242       std::vector<int> indices(length + 2);
     176  void GlpkBase::_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);
    243183      std::vector<Value> values(length + 2);
    244184
    245       LEMON_glp(get_mat_row)(lp, ix, &indices[0], &values[0]);
     185      glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
    246186
    247187      //The following code does not suppose that the elements of the
    248       //array indices are sorted
    249       bool found=false;
    250       for (int i = 1; i <= length; ++i) {
    251         if (indices[i]==jx){
    252           found=true;
    253           values[i]=value;
     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;
    254194          break;
    255195        }
    256196      }
    257       if (!found){
     197      if (!found) {
    258198        ++length;
    259         indices[length]=jx;
    260         values[length]=value;
    261       }
    262 
    263       LEMON_glp(set_mat_row)(lp, ix, length, &indices[0], &values[0]);
     199        indexes[length] = jx;
     200        values[length] = value;
     201      }
     202
     203      glp_set_mat_row(lp, ix, length, &indexes.front(), &values.front());
    264204
    265205    } else {
    266206
    267       int length=LEMON_glp(get_mat_col)(lp, jx, 0, 0);
    268 
    269       std::vector<int> indices(length + 2);
     207      int length = glp_get_mat_col(lp, jx, 0, 0);
     208
     209      std::vector<int> indexes(length + 2);
    270210      std::vector<Value> values(length + 2);
    271211
    272       LEMON_glp(get_mat_col)(lp, jx, &indices[0], &values[0]);
     212      glp_get_mat_col(lp, jx, &indexes.front(), &values.front());
    273213
    274214      //The following code does not suppose that the elements of the
    275       //array indices are sorted
    276       bool found=false;
     215      //array indexes are sorted
     216      bool found = false;
    277217      for (int i = 1; i <= length; ++i) {
    278         if (indices[i]==ix){
    279           found=true;
    280           values[i]=value;
     218        if (indexes[i] == ix) {
     219          found = true;
     220          values[i] = value;
    281221          break;
    282222        }
    283223      }
    284       if (!found){
     224      if (!found) {
    285225        ++length;
    286         indices[length]=ix;
    287         values[length]=value;
    288       }
    289 
    290       LEMON_glp(set_mat_col)(lp, jx, length, &indices[0], &values[0]);
    291     }
    292 
    293     solved = false;
    294   }
    295 
    296   LpGlpk::Value LpGlpk::_getCoeff(int ix, int jx) const
    297   {
    298 
    299     int length=LEMON_glp(get_mat_row)(lp, ix, 0, 0);
    300 
    301     std::vector<int> indices(length + 1);
     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  GlpkBase::Value GlpkBase::_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);
    302240    std::vector<Value> values(length + 1);
    303241
    304     LEMON_glp(get_mat_row)(lp, ix, &indices[0], &values[0]);
    305 
    306     //The following code does not suppose that the elements of the
    307     //array indices are sorted
    308     for (int i = 1; i <= length; ++i) {
    309       if (indices[i]==jx){
     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) {
    310246        return values[i];
    311247      }
    312248    }
     249
    313250    return 0;
    314 
    315   }
    316 
    317 
    318   void LpGlpk::_setColLowerBound(int i, Value lo)
    319   {
    320     if (lo==INF) {
    321       //FIXME error
    322     }
    323     int b=LEMON_glp(get_col_type)(lp, i);
    324     double up=LEMON_glp(get_col_ub)(lp, i);
    325     if (lo==-INF) {
     251  }
     252
     253  void GlpkBase::_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) {
    326259      switch (b) {
    327       case LEMON_GLP(FR):
    328       case LEMON_GLP(LO):
    329         LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(FR), lo, up);
    330         break;
    331       case LEMON_GLP(UP):
    332         break;
    333       case LEMON_GLP(DB):
    334       case LEMON_GLP(FX):
    335         LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(UP), lo, up);
    336         break;
    337       default: ;
    338         //FIXME error
     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;
    339272      }
    340273    } else {
    341274      switch (b) {
    342       case LEMON_GLP(FR):
    343       case LEMON_GLP(LO):
    344         LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(LO), lo, up);
    345         break;
    346       case LEMON_GLP(UP):
    347       case LEMON_GLP(DB):
    348       case LEMON_GLP(FX):
    349         if (lo==up)
    350           LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(FX), lo, up);
     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);
    351284        else
    352           LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(DB), lo, up);
    353         break;
    354       default: ;
    355         //FIXME error
    356       }
    357     }
    358 
    359     solved = false;
    360   }
    361 
    362   LpGlpk::Value LpGlpk::_getColLowerBound(int i) const
    363   {
    364     int b=LEMON_glp(get_col_type)(lp, i);
     285          glp_set_col_bnds(lp, i, GLP_DB, lo, up);
     286        break;
     287      default:
     288        break;
     289      }
     290    }
     291  }
     292
     293  GlpkBase::Value GlpkBase::_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 GlpkBase::_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) {
    365311      switch (b) {
    366       case LEMON_GLP(LO):
    367       case LEMON_GLP(DB):
    368       case LEMON_GLP(FX):
    369         return LEMON_glp(get_col_lb)(lp, i);
    370       default: ;
    371         return -INF;
    372       }
    373   }
    374 
    375   void LpGlpk::_setColUpperBound(int i, Value up)
    376   {
    377     if (up==-INF) {
    378       //FIXME error
    379     }
    380     int b=LEMON_glp(get_col_type)(lp, i);
    381     double lo=LEMON_glp(get_col_lb)(lp, i);
    382     if (up==INF) {
    383       switch (b) {
    384       case LEMON_GLP(FR):
    385       case LEMON_GLP(LO):
    386         break;
    387       case LEMON_GLP(UP):
    388         LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(FR), lo, up);
    389         break;
    390       case LEMON_GLP(DB):
    391       case LEMON_GLP(FX):
    392         LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(LO), lo, up);
    393         break;
    394       default: ;
    395         //FIXME error
     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;
    396324      }
    397325    } else {
    398326      switch (b) {
    399       case LEMON_GLP(FR):
    400         LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(UP), lo, up);
    401         break;
    402       case LEMON_GLP(UP):
    403         LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(UP), lo, up);
    404         break;
    405       case LEMON_GLP(LO):
    406       case LEMON_GLP(DB):
    407       case LEMON_GLP(FX):
    408         if (lo==up)
    409           LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(FX), lo, up);
     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);
    410338        else
    411           LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(DB), lo, up);
    412         break;
    413       default: ;
    414         //FIXME error
    415       }
    416     }
    417 
    418     solved = false;
    419   }
    420 
    421   LpGlpk::Value LpGlpk::_getColUpperBound(int i) const
    422   {
    423     int b=LEMON_glp(get_col_type)(lp, i);
     339          glp_set_col_bnds(lp, i, GLP_DB, lo, up);
     340        break;
     341      default:
     342        break;
     343      }
     344    }
     345
     346  }
     347
     348  GlpkBase::Value GlpkBase::_getColUpperBound(int i) const {
     349    int b = glp_get_col_type(lp, i);
    424350      switch (b) {
    425       case LEMON_GLP(UP):
    426       case LEMON_GLP(DB):
    427       case LEMON_GLP(FX):
    428         return LEMON_glp(get_col_ub)(lp, i);
    429       default: ;
     351      case GLP_UP:
     352      case GLP_DB:
     353      case GLP_FX:
     354        return glp_get_col_ub(lp, i);
     355      default:
    430356        return INF;
    431357      }
    432358  }
    433359
    434   void LpGlpk::_setRowBounds(int i, Value lb, Value ub)
    435   {
    436     //Bad parameter
    437     if (lb==INF || ub==-INF) {
    438       //FIXME error
    439     }
    440 
    441     if (lb == -INF){
    442       if (ub == INF){
    443         LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(FR), lb, ub);
    444       }
    445       else{
    446         LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(UP), lb, ub);
    447       }
    448     }
    449     else{
    450       if (ub==INF){
    451         LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(LO), lb, ub);
    452 
    453       }
    454       else{
    455         if (lb == ub){
    456           LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(FX), lb, ub);
     360  void GlpkBase::_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  GlpkBase::Value GlpkBase::_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 GlpkBase::_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  GlpkBase::Value GlpkBase::_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 GlpkBase::_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 GlpkBase::_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 GlpkBase::_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  GlpkBase::Value GlpkBase::_getObjCoeff(int i) const {
     492    //i = 0 means the constant term (shift)
     493    return glp_get_obj_coef(lp, i);
     494  }
     495
     496  void GlpkBase::_setSense(GlpkBase::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  GlpkBase::Sense GlpkBase::_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 GlpkBase::Sense();
     516    }
     517  }
     518
     519  void GlpkBase::_clear() {
     520    glp_erase_prob(lp);
     521    rows.clear();
     522    cols.clear();
     523  }
     524
     525  // LpGlpk members
     526
     527  LpGlpk::LpGlpk()
     528    : LpBase(), GlpkBase(), LpSolver() {
     529    messageLevel(MESSAGE_NO_OUTPUT);
     530  }
     531
     532  LpGlpk::LpGlpk(const LpGlpk& other)
     533    : LpBase(other), GlpkBase(other), LpSolver(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  const char* LpGlpk::_solverName() const { return "LpGlpk"; }
     541
     542  void LpGlpk::_clear_temporals() {
     543    _primal_ray.clear();
     544    _dual_ray.clear();
     545  }
     546
     547  LpGlpk::SolveExitStatus LpGlpk::_solve() {
     548    return solvePrimal();
     549  }
     550
     551  LpGlpk::SolveExitStatus LpGlpk::solvePrimal() {
     552    _clear_temporals();
     553
     554    glp_smcp smcp;
     555    glp_init_smcp(&smcp);
     556
     557    switch (_message_level) {
     558    case MESSAGE_NO_OUTPUT:
     559      smcp.msg_lev = GLP_MSG_OFF;
     560      break;
     561    case MESSAGE_ERROR_MESSAGE:
     562      smcp.msg_lev = GLP_MSG_ERR;
     563      break;
     564    case MESSAGE_NORMAL_OUTPUT:
     565      smcp.msg_lev = GLP_MSG_ON;
     566      break;
     567    case MESSAGE_FULL_OUTPUT:
     568      smcp.msg_lev = GLP_MSG_ALL;
     569      break;
     570    }
     571
     572    if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
     573    return SOLVED;
     574  }
     575
     576  LpGlpk::SolveExitStatus LpGlpk::solveDual() {
     577    _clear_temporals();
     578
     579    glp_smcp smcp;
     580    glp_init_smcp(&smcp);
     581
     582    switch (_message_level) {
     583    case MESSAGE_NO_OUTPUT:
     584      smcp.msg_lev = GLP_MSG_OFF;
     585      break;
     586    case MESSAGE_ERROR_MESSAGE:
     587      smcp.msg_lev = GLP_MSG_ERR;
     588      break;
     589    case MESSAGE_NORMAL_OUTPUT:
     590      smcp.msg_lev = GLP_MSG_ON;
     591      break;
     592    case MESSAGE_FULL_OUTPUT:
     593      smcp.msg_lev = GLP_MSG_ALL;
     594      break;
     595    }
     596    smcp.meth = GLP_DUAL;
     597
     598    if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
     599    return SOLVED;
     600  }
     601
     602  LpGlpk::Value LpGlpk::_getPrimal(int i) const {
     603    return glp_get_col_prim(lp, i);
     604  }
     605
     606  LpGlpk::Value LpGlpk::_getDual(int i) const {
     607    return glp_get_row_dual(lp, i);
     608  }
     609
     610  LpGlpk::Value LpGlpk::_getPrimalValue() const {
     611    return glp_get_obj_val(lp);
     612  }
     613
     614  LpGlpk::VarStatus LpGlpk::_getColStatus(int i) const {
     615    switch (glp_get_col_stat(lp, i)) {
     616    case GLP_BS:
     617      return BASIC;
     618    case GLP_UP:
     619      return UPPER;
     620    case GLP_LO:
     621      return LOWER;
     622    case GLP_NF:
     623      return FREE;
     624    case GLP_NS:
     625      return FIXED;
     626    default:
     627      LEMON_ASSERT(false, "Wrong column status");
     628      return LpGlpk::VarStatus();
     629    }
     630  }
     631
     632  LpGlpk::VarStatus LpGlpk::_getRowStatus(int i) const {
     633    switch (glp_get_row_stat(lp, i)) {
     634    case GLP_BS:
     635      return BASIC;
     636    case GLP_UP:
     637      return UPPER;
     638    case GLP_LO:
     639      return LOWER;
     640    case GLP_NF:
     641      return FREE;
     642    case GLP_NS:
     643      return FIXED;
     644    default:
     645      LEMON_ASSERT(false, "Wrong row status");
     646      return LpGlpk::VarStatus();
     647    }
     648  }
     649
     650  LpGlpk::Value LpGlpk::_getPrimalRay(int i) const {
     651    if (_primal_ray.empty()) {
     652      int row_num = glp_get_num_rows(lp);
     653      int col_num = glp_get_num_cols(lp);
     654
     655      _primal_ray.resize(col_num + 1, 0.0);
     656
     657      int index = glp_get_unbnd_ray(lp);
     658      if (index != 0) {
     659        // The primal ray is found in primal simplex second phase
     660        LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :
     661                      glp_get_col_stat(lp, index - row_num)) != GLP_BS,
     662                     "Wrong primal ray");
     663
     664        bool negate = glp_get_obj_dir(lp) == GLP_MAX;
     665
     666        if (index > row_num) {
     667          _primal_ray[index - row_num] = 1.0;
     668          if (glp_get_col_dual(lp, index - row_num) > 0) {
     669            negate = !negate;
     670          }
     671        } else {
     672          if (glp_get_row_dual(lp, index) > 0) {
     673            negate = !negate;
     674          }
    457675        }
    458         else{
    459           LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(DB), lb, ub);
     676
     677        std::vector<int> ray_indexes(row_num + 1);
     678        std::vector<Value> ray_values(row_num + 1);
     679        int ray_length = glp_eval_tab_col(lp, index, &ray_indexes.front(),
     680                                          &ray_values.front());
     681
     682        for (int i = 1; i <= ray_length; ++i) {
     683          if (ray_indexes[i] > row_num) {
     684            _primal_ray[ray_indexes[i] - row_num] = ray_values[i];
     685          }
    460686        }
    461       }
    462     }
    463 
    464     solved = false;
    465   }
    466 
    467   void LpGlpk::_getRowBounds(int i, Value &lb, Value &ub) const
    468   {
    469 
    470     int b=LEMON_glp(get_row_type)(lp, i);
    471     switch (b) {
    472     case LEMON_GLP(FR):
    473     case LEMON_GLP(UP):
    474       lb = -INF;
    475         break;
     687
     688        if (negate) {
     689          for (int i = 1; i <= col_num; ++i) {
     690            _primal_ray[i] = - _primal_ray[i];
     691          }
     692        }
     693      } else {
     694        for (int i = 1; i <= col_num; ++i) {
     695          _primal_ray[i] = glp_get_col_prim(lp, i);
     696        }
     697      }
     698    }
     699    return _primal_ray[i];
     700  }
     701
     702  LpGlpk::Value LpGlpk::_getDualRay(int i) const {
     703    if (_dual_ray.empty()) {
     704      int row_num = glp_get_num_rows(lp);
     705
     706      _dual_ray.resize(row_num + 1, 0.0);
     707
     708      int index = glp_get_unbnd_ray(lp);
     709      if (index != 0) {
     710        // The dual ray is found in dual simplex second phase
     711        LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :
     712                      glp_get_col_stat(lp, index - row_num)) == GLP_BS,
     713
     714                     "Wrong dual ray");
     715
     716        int idx;
     717        bool negate = false;
     718
     719        if (index > row_num) {
     720          idx = glp_get_col_bind(lp, index - row_num);
     721          if (glp_get_col_prim(lp, index - row_num) >
     722              glp_get_col_ub(lp, index - row_num)) {
     723            negate = true;
     724          }
     725        } else {
     726          idx = glp_get_row_bind(lp, index);
     727          if (glp_get_row_prim(lp, index) > glp_get_row_ub(lp, index)) {
     728            negate = true;
     729          }
     730        }
     731
     732        _dual_ray[idx] = negate ?  - 1.0 : 1.0;
     733
     734        glp_btran(lp, &_dual_ray.front());
     735      } else {
     736        double eps = 1e-7;
     737        // The dual ray is found in primal simplex first phase
     738        // We assume that the glpk minimizes the slack to get feasible solution
     739        for (int i = 1; i <= row_num; ++i) {
     740          int index = glp_get_bhead(lp, i);
     741          if (index <= row_num) {
     742            double res = glp_get_row_prim(lp, index);
     743            if (res > glp_get_row_ub(lp, index) + eps) {
     744              _dual_ray[i] = -1;
     745            } else if (res < glp_get_row_lb(lp, index) - eps) {
     746              _dual_ray[i] = 1;
     747            } else {
     748              _dual_ray[i] = 0;
     749            }
     750            _dual_ray[i] *= glp_get_rii(lp, index);
     751          } else {
     752            double res = glp_get_col_prim(lp, index - row_num);
     753            if (res > glp_get_col_ub(lp, index - row_num) + eps) {
     754              _dual_ray[i] = -1;
     755            } else if (res < glp_get_col_lb(lp, index - row_num) - eps) {
     756              _dual_ray[i] = 1;
     757            } else {
     758              _dual_ray[i] = 0;
     759            }
     760            _dual_ray[i] /= glp_get_sjj(lp, index - row_num);
     761          }
     762        }
     763
     764        glp_btran(lp, &_dual_ray.front());
     765
     766        for (int i = 1; i <= row_num; ++i) {
     767          _dual_ray[i] /= glp_get_rii(lp, i);
     768        }
     769      }
     770    }
     771    return _dual_ray[i];
     772  }
     773
     774  LpGlpk::ProblemType LpGlpk::_getPrimalType() const {
     775    if (glp_get_status(lp) == GLP_OPT)
     776      return OPTIMAL;
     777    switch (glp_get_prim_stat(lp)) {
     778    case GLP_UNDEF:
     779      return UNDEFINED;
     780    case GLP_FEAS:
     781    case GLP_INFEAS:
     782      if (glp_get_dual_stat(lp) == GLP_NOFEAS) {
     783        return UNBOUNDED;
     784      } else {
     785        return UNDEFINED;
     786      }
     787    case GLP_NOFEAS:
     788      return INFEASIBLE;
    476789    default:
    477       lb=LEMON_glp(get_row_lb)(lp, i);
    478     }
    479 
    480     switch (b) {
    481     case LEMON_GLP(FR):
    482     case LEMON_GLP(LO):
    483       ub = INF;
    484         break;
     790      LEMON_ASSERT(false, "Wrong primal type");
     791      return  LpGlpk::ProblemType();
     792    }
     793  }
     794
     795  LpGlpk::ProblemType LpGlpk::_getDualType() const {
     796    if (glp_get_status(lp) == GLP_OPT)
     797      return OPTIMAL;
     798    switch (glp_get_dual_stat(lp)) {
     799    case GLP_UNDEF:
     800      return UNDEFINED;
     801    case GLP_FEAS:
     802    case GLP_INFEAS:
     803      if (glp_get_prim_stat(lp) == GLP_NOFEAS) {
     804        return UNBOUNDED;
     805      } else {
     806        return UNDEFINED;
     807      }
     808    case GLP_NOFEAS:
     809      return INFEASIBLE;
    485810    default:
    486       ub=LEMON_glp(get_row_ub)(lp, i);
    487     }
    488 
    489   }
    490 
    491   void LpGlpk::_setObjCoeff(int i, Value obj_coef)
    492   {
    493     //i=0 means the constant term (shift)
    494     LEMON_glp(set_obj_coef)(lp, i, obj_coef);
    495 
    496     solved = false;
    497   }
    498 
    499   LpGlpk::Value LpGlpk::_getObjCoeff(int i) const {
    500     //i=0 means the constant term (shift)
    501     return LEMON_glp(get_obj_coef)(lp, i);
    502   }
    503 
    504   void LpGlpk::_clearObj()
    505   {
    506     for (int i=0;i<=LEMON_glp(get_num_cols)(lp);++i){
    507       LEMON_glp(set_obj_coef)(lp, i, 0);
    508     }
    509 
    510     solved = false;
    511   }
    512 
    513   LpGlpk::SolveExitStatus LpGlpk::_solve()
    514   {
    515     // A way to check the problem to be solved
    516     //LEMON_glp(write_cpxlp(lp,"naittvan.cpx");
    517 
    518     LEMON_lpx(std_basis)(lp);
    519     int i =  LEMON_lpx(simplex)(lp);
    520 
    521     switch (i) {
    522     case LEMON_LPX(E_OK):
    523       solved = true;
    524       return SOLVED;
     811      LEMON_ASSERT(false, "Wrong primal type");
     812      return  LpGlpk::ProblemType();
     813    }
     814  }
     815
     816  void LpGlpk::presolver(bool b) {
     817    lpx_set_int_parm(lp, LPX_K_PRESOL, b ? 1 : 0);
     818  }
     819
     820  void LpGlpk::messageLevel(MessageLevel m) {
     821    _message_level = m;
     822  }
     823
     824  // MipGlpk members
     825
     826  MipGlpk::MipGlpk()
     827    : LpBase(), GlpkBase(), MipSolver() {
     828    messageLevel(MESSAGE_NO_OUTPUT);
     829  }
     830
     831  MipGlpk::MipGlpk(const MipGlpk& other)
     832    : LpBase(), GlpkBase(other), MipSolver() {
     833    messageLevel(MESSAGE_NO_OUTPUT);
     834  }
     835
     836  void MipGlpk::_setColType(int i, MipGlpk::ColTypes col_type) {
     837    switch (col_type) {
     838    case INTEGER:
     839      glp_set_col_kind(lp, i, GLP_IV);
     840      break;
     841    case REAL:
     842      glp_set_col_kind(lp, i, GLP_CV);
     843      break;
     844    }
     845  }
     846
     847  MipGlpk::ColTypes MipGlpk::_getColType(int i) const {
     848    switch (glp_get_col_kind(lp, i)) {
     849    case GLP_IV:
     850    case GLP_BV:
     851      return INTEGER;
    525852    default:
    526       return UNSOLVED;
    527     }
    528   }
    529 
    530   LpGlpk::Value LpGlpk::_getPrimal(int i) const
    531   {
    532     return LEMON_glp(get_col_prim)(lp,i);
    533   }
    534 
    535   LpGlpk::Value LpGlpk::_getDual(int i) const
    536   {
    537     return LEMON_glp(get_row_dual)(lp,i);
    538   }
    539 
    540   LpGlpk::Value LpGlpk::_getPrimalValue() const
    541   {
    542     return LEMON_glp(get_obj_val)(lp);
    543   }
    544   bool LpGlpk::_isBasicCol(int i) const
    545   {
    546     return (LEMON_glp(get_col_stat)(lp, i)==LEMON_GLP(BS));
    547   }
    548 
    549 
    550   LpGlpk::SolutionStatus LpGlpk::_getPrimalStatus() const
    551   {
    552     if (!solved) return UNDEFINED;
    553     int stat=  LEMON_lpx(get_status)(lp);
    554     switch (stat) {
    555     case LEMON_LPX(UNDEF)://Undefined (no solve has been run yet)
    556       return UNDEFINED;
    557     case LEMON_LPX(NOFEAS)://There is no feasible solution (primal, I guess)
    558     case LEMON_LPX(INFEAS)://Infeasible
    559       return INFEASIBLE;
    560     case LEMON_LPX(UNBND)://Unbounded
    561       return INFINITE;
    562     case LEMON_LPX(FEAS)://Feasible
    563       return FEASIBLE;
    564     case LEMON_LPX(OPT)://Feasible
    565       return OPTIMAL;
    566     default:
    567       return UNDEFINED; //to avoid gcc warning
    568       //FIXME error
    569     }
    570   }
    571 
    572   LpGlpk::SolutionStatus LpGlpk::_getDualStatus() const
    573   {
    574     if (!solved) return UNDEFINED;
    575     switch (LEMON_lpx(get_dual_stat)(lp)) {
    576     case LEMON_LPX(D_UNDEF)://Undefined (no solve has been run yet)
    577       return UNDEFINED;
    578     case LEMON_LPX(D_NOFEAS)://There is no dual feasible solution
    579 //    case LEMON_LPX(D_INFEAS://Infeasible
    580       return INFEASIBLE;
    581     case LEMON_LPX(D_FEAS)://Feasible
    582       switch (LEMON_lpx(get_status)(lp)) {
    583       case LEMON_LPX(NOFEAS):
    584         return INFINITE;
    585       case LEMON_LPX(OPT):
     853      return REAL;
     854    }
     855
     856  }
     857
     858  MipGlpk::SolveExitStatus MipGlpk::_solve() {
     859    glp_smcp smcp;
     860    glp_init_smcp(&smcp);
     861
     862    switch (_message_level) {
     863    case MESSAGE_NO_OUTPUT:
     864      smcp.msg_lev = GLP_MSG_OFF;
     865      break;
     866    case MESSAGE_ERROR_MESSAGE:
     867      smcp.msg_lev = GLP_MSG_ERR;
     868      break;
     869    case MESSAGE_NORMAL_OUTPUT:
     870      smcp.msg_lev = GLP_MSG_ON;
     871      break;
     872    case MESSAGE_FULL_OUTPUT:
     873      smcp.msg_lev = GLP_MSG_ALL;
     874      break;
     875    }
     876    smcp.meth = GLP_DUAL;
     877
     878    if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
     879    if (glp_get_status(lp) != GLP_OPT) return SOLVED;
     880
     881    glp_iocp iocp;
     882    glp_init_iocp(&iocp);
     883
     884    switch (_message_level) {
     885    case MESSAGE_NO_OUTPUT:
     886      iocp.msg_lev = GLP_MSG_OFF;
     887      break;
     888    case MESSAGE_ERROR_MESSAGE:
     889      iocp.msg_lev = GLP_MSG_ERR;
     890      break;
     891    case MESSAGE_NORMAL_OUTPUT:
     892      iocp.msg_lev = GLP_MSG_ON;
     893      break;
     894    case MESSAGE_FULL_OUTPUT:
     895      iocp.msg_lev = GLP_MSG_ALL;
     896      break;
     897    }
     898
     899    if (glp_intopt(lp, &iocp) != 0) return UNSOLVED;
     900    return SOLVED;
     901  }
     902
     903
     904  MipGlpk::ProblemType MipGlpk::_getType() const {
     905    switch (glp_get_status(lp)) {
     906    case GLP_OPT:
     907      switch (glp_mip_status(lp)) {
     908      case GLP_UNDEF:
     909        return UNDEFINED;
     910      case GLP_NOFEAS:
     911        return INFEASIBLE;
     912      case GLP_FEAS:
     913        return FEASIBLE;
     914      case GLP_OPT:
    586915        return OPTIMAL;
    587916      default:
    588         return FEASIBLE;
     917        LEMON_ASSERT(false, "Wrong problem type.");
     918        return MipGlpk::ProblemType();
     919      }
     920    case GLP_NOFEAS:
     921      return INFEASIBLE;
     922    case GLP_INFEAS:
     923    case GLP_FEAS:
     924      if (glp_get_dual_stat(lp) == GLP_NOFEAS) {
     925        return UNBOUNDED;
     926      } else {
     927        return UNDEFINED;
    589928      }
    590929    default:
    591       return UNDEFINED; //to avoid gcc warning
    592       //FIXME error
    593     }
    594   }
    595 
    596   LpGlpk::ProblemTypes LpGlpk::_getProblemType() const
    597   {
    598     if (!solved) return UNKNOWN;
    599       //int stat=  LEMON_glp(get_status(lp);
    600     int statp=  LEMON_lpx(get_prim_stat)(lp);
    601     int statd=  LEMON_lpx(get_dual_stat)(lp);
    602     if (statp==LEMON_LPX(P_FEAS) && statd==LEMON_LPX(D_FEAS))
    603         return PRIMAL_DUAL_FEASIBLE;
    604     if (statp==LEMON_LPX(P_FEAS) && statd==LEMON_LPX(D_NOFEAS))
    605         return PRIMAL_FEASIBLE_DUAL_INFEASIBLE;
    606     if (statp==LEMON_LPX(P_NOFEAS) && statd==LEMON_LPX(D_FEAS))
    607         return PRIMAL_INFEASIBLE_DUAL_FEASIBLE;
    608     if (statp==LEMON_LPX(P_NOFEAS) && statd==LEMON_LPX(D_NOFEAS))
    609         return PRIMAL_DUAL_INFEASIBLE;
    610     //In all other cases
    611     return UNKNOWN;
    612   }
    613 
    614   void LpGlpk::_setMax()
    615   {
    616     solved = false;
    617     LEMON_glp(set_obj_dir)(lp, LEMON_GLP(MAX));
    618   }
    619 
    620   void LpGlpk::_setMin()
    621   {
    622     solved = false;
    623     LEMON_glp(set_obj_dir)(lp, LEMON_GLP(MIN));
    624   }
    625 
    626   bool LpGlpk::_isMax() const
    627   {
    628     return (LEMON_glp(get_obj_dir)(lp)==LEMON_GLP(MAX));
    629   }
    630 
    631 
    632 
    633   void LpGlpk::messageLevel(int m)
    634   {
    635     LEMON_lpx(set_int_parm)(lp, LEMON_LPX(K_MSGLEV), m);
    636   }
    637 
    638   void LpGlpk::presolver(bool b)
    639   {
    640     LEMON_lpx(set_int_parm)(lp, LEMON_LPX(K_PRESOL), b);
    641   }
    642 
     930      LEMON_ASSERT(false, "Wrong problem type.");
     931      return MipGlpk::ProblemType();
     932    }
     933  }
     934
     935  MipGlpk::Value MipGlpk::_getSol(int i) const {
     936    return glp_mip_col_val(lp, i);
     937  }
     938
     939  MipGlpk::Value MipGlpk::_getSolValue() const {
     940    return glp_mip_obj_val(lp);
     941  }
     942
     943  MipGlpk* MipGlpk::_newSolver() const { return new MipGlpk; }
     944  MipGlpk* MipGlpk::_cloneSolver() const {return new MipGlpk(*this); }
     945
     946  const char* MipGlpk::_solverName() const { return "MipGlpk"; }
     947
     948  void MipGlpk::messageLevel(MessageLevel m) {
     949    _message_level = m;
     950  }
    643951
    644952} //END OF NAMESPACE LEMON
  • lemon/lp_glpk.h

    r481 r482  
    3636
    3737
     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 GlpkBase : virtual public LpBase {
     43  protected:
     44
     45    typedef glp_prob LPX;
     46    glp_prob* lp;
     47
     48    GlpkBase();
     49    GlpkBase(const GlpkBase&);
     50    virtual ~GlpkBase();
     51
     52  protected:
     53
     54    virtual int _addCol();
     55    virtual int _addRow();
     56
     57    virtual void _eraseCol(int i);
     58    virtual void _eraseRow(int i);
     59
     60    virtual void _eraseColId(int i);
     61    virtual void _eraseRowId(int i);
     62
     63    virtual void _getColName(int col, std::string& name) const;
     64    virtual void _setColName(int col, const std::string& name);
     65    virtual int _colByName(const std::string& name) const;
     66
     67    virtual void _getRowName(int row, std::string& name) const;
     68    virtual void _setRowName(int row, const std::string& name);
     69    virtual int _rowByName(const std::string& name) const;
     70
     71    virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e);
     72    virtual void _getRowCoeffs(int i, InsertIterator b) const;
     73
     74    virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e);
     75    virtual void _getColCoeffs(int i, InsertIterator b) const;
     76
     77    virtual void _setCoeff(int row, int col, Value value);
     78    virtual Value _getCoeff(int row, int col) const;
     79
     80    virtual void _setColLowerBound(int i, Value value);
     81    virtual Value _getColLowerBound(int i) const;
     82
     83    virtual void _setColUpperBound(int i, Value value);
     84    virtual Value _getColUpperBound(int i) const;
     85
     86    virtual void _setRowLowerBound(int i, Value value);
     87    virtual Value _getRowLowerBound(int i) const;
     88
     89    virtual void _setRowUpperBound(int i, Value value);
     90    virtual Value _getRowUpperBound(int i) const;
     91
     92    virtual void _setObjCoeffs(ExprIterator b, ExprIterator e);
     93    virtual void _getObjCoeffs(InsertIterator b) const;
     94
     95    virtual void _setObjCoeff(int i, Value obj_coef);
     96    virtual Value _getObjCoeff(int i) const;
     97
     98    virtual void _setSense(Sense);
     99    virtual Sense _getSense() const;
     100
     101    virtual void _clear();
     102
     103  public:
     104
     105    ///Pointer to the underlying GLPK data structure.
     106    LPX *lpx() {return lp;}
     107    ///Const pointer to the underlying GLPK data structure.
     108    const LPX *lpx() const {return lp;}
     109
     110    ///Returns the constraint identifier understood by GLPK.
     111    int lpxRow(Row r) const { return rows(id(r)); }
     112
     113    ///Returns the variable identifier understood by GLPK.
     114    int lpxCol(Col c) const { return cols(id(c)); }
     115
     116  };
     117
    38118  /// \brief Interface for the GLPK LP solver
    39119  ///
    40120  /// This class implements an interface for the GLPK LP solver.
    41121  ///\ingroup lp_group
    42   class LpGlpk : virtual public LpSolverBase {
    43   protected:
    44 
    45     typedef glp_prob LPX;
    46     glp_prob* lp;
    47     bool solved;
    48 
    49   public:
    50 
    51     typedef LpSolverBase Parent;
    52 
     122  class LpGlpk : public GlpkBase, public LpSolver {
     123  public:
     124
     125    ///\e
    53126    LpGlpk();
    54     LpGlpk(const LpGlpk &);
    55     ~LpGlpk();
    56 
    57   protected:
    58     virtual LpSolverBase* _newLp();
    59     virtual LpSolverBase* _copyLp();
    60 
    61     virtual int _addCol();
    62     virtual int _addRow();
    63     virtual void _eraseCol(int i);
    64     virtual void _eraseRow(int i);
    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     virtual void _setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e);
    69     virtual void _getRowCoeffs(int i, RowIterator b) const;
    70     virtual void _setColCoeffs(int i, ConstColIterator b, ConstColIterator e);
    71     virtual void _getColCoeffs(int i, ColIterator b) const;
    72     virtual void _setCoeff(int row, int col, Value value);
    73     virtual Value _getCoeff(int row, int col) const;
    74 
    75     virtual void _setColLowerBound(int i, Value value);
    76     virtual Value _getColLowerBound(int i) const;
    77     virtual void _setColUpperBound(int i, Value value);
    78     virtual Value _getColUpperBound(int i) const;
    79 
    80     virtual void _setRowBounds(int i, Value lower, Value upper);
    81     virtual void _getRowBounds(int i, Value &lb, Value &ub) const;
    82     virtual void _setObjCoeff(int i, Value obj_coef);
    83     virtual Value _getObjCoeff(int i) const;
    84     virtual void _clearObj();
    85 
    86     ///\e
    87 
    88     ///\todo It should be clarified
    89     ///
     127    ///\e
     128    LpGlpk(const LpGlpk&);
     129
     130  private:
     131
     132    mutable std::vector<double> _primal_ray;
     133    mutable std::vector<double> _dual_ray;
     134
     135    void _clear_temporals();
     136
     137  protected:
     138
     139    virtual LpGlpk* _cloneSolver() const;
     140    virtual LpGlpk* _newSolver() const;
     141
     142    virtual const char* _solverName() const;
     143
    90144    virtual SolveExitStatus _solve();
    91145    virtual Value _getPrimal(int i) const;
    92146    virtual Value _getDual(int i) const;
     147
    93148    virtual Value _getPrimalValue() const;
    94     virtual bool _isBasicCol(int i) const;
    95     ///\e
     149
     150    virtual VarStatus _getColStatus(int i) const;
     151    virtual VarStatus _getRowStatus(int i) const;
     152
     153    virtual Value _getPrimalRay(int i) const;
     154    virtual Value _getDualRay(int i) const;
    96155
    97156    ///\todo It should be clarified
    98157    ///
    99     virtual SolutionStatus _getPrimalStatus() const;
    100     virtual SolutionStatus _getDualStatus() const;
    101     virtual ProblemTypes _getProblemType() const;
    102 
    103     virtual void _setMax();
    104     virtual void _setMin();
    105 
    106     virtual bool _isMax() const;
    107 
    108   public:
    109     ///Set the verbosity of the messages
    110 
    111     ///Set the verbosity of the messages
    112     ///
    113     ///\param m is the level of the messages output by the solver routines.
    114     ///The possible values are:
    115     ///- 0 --- no output (default value)
    116     ///- 1 --- error messages only
    117     ///- 2 --- normal output
    118     ///- 3 --- full output (includes informational messages)
    119     void messageLevel(int m);
     158    virtual ProblemType _getPrimalType() const;
     159    virtual ProblemType _getDualType() const;
     160
     161  public:
     162
     163    ///Solve with primal simplex
     164    SolveExitStatus solvePrimal();
     165
     166    ///Solve with dual simplex
     167    SolveExitStatus solveDual();
     168
    120169    ///Turns on or off the presolver
    121170
     
    125174    void presolver(bool b);
    126175
    127     ///Pointer to the underlying GLPK data structure.
    128     LPX *lpx() {return lp;}
    129 
    130     ///Returns the constraint identifier understood by GLPK.
    131     int lpxRow(Row r) { return _lpId(r); }
    132 
    133     ///Returns the variable identifier understood by GLPK.
    134     int lpxCol(Col c) { return _lpId(c); }
     176    ///Enum for \c messageLevel() parameter
     177    enum MessageLevel {
     178      /// no output (default value)
     179      MESSAGE_NO_OUTPUT = 0,
     180      /// error messages only
     181      MESSAGE_ERROR_MESSAGE = 1,
     182      /// normal output
     183      MESSAGE_NORMAL_OUTPUT = 2,
     184      /// full output (includes informational messages)
     185      MESSAGE_FULL_OUTPUT = 3
     186    };
     187
     188  private:
     189
     190    MessageLevel _message_level;
     191
     192  public:
     193
     194    ///Set the verbosity of the messages
     195
     196    ///Set the verbosity of the messages
     197    ///
     198    ///\param m is the level of the messages output by the solver routines.
     199    void messageLevel(MessageLevel m);
    135200  };
     201
     202  /// \brief Interface for the GLPK MIP solver
     203  ///
     204  /// This class implements an interface for the GLPK MIP solver.
     205  ///\ingroup lp_group
     206  class MipGlpk : public GlpkBase, public MipSolver {
     207  public:
     208
     209    ///\e
     210    MipGlpk();
     211    ///\e
     212    MipGlpk(const MipGlpk&);
     213
     214  protected:
     215
     216    virtual MipGlpk* _cloneSolver() const;
     217    virtual MipGlpk* _newSolver() const;
     218
     219    virtual const char* _solverName() 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