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
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();