lemon/lp_glpk.cc
changeset 459 ed54c0d13df0
parent 458 7afc121e0689
     1.1 --- a/lemon/lp_glpk.cc	Tue Dec 02 21:40:33 2008 +0100
     1.2 +++ b/lemon/lp_glpk.cc	Tue Dec 02 22:48:28 2008 +0100
     1.3 @@ -17,628 +17,936 @@
     1.4   */
     1.5  
     1.6  ///\file
     1.7 -///\brief Implementation of the LEMON-GLPK lp solver interface.
     1.8 +///\brief Implementation of the LEMON GLPK LP and MIP solver interface.
     1.9  
    1.10  #include <lemon/lp_glpk.h>
    1.11 -//#include <iostream>
    1.12 +#include <glpk.h>
    1.13  
    1.14 -extern "C" {
    1.15 -#include <glpk.h>
    1.16 -}
    1.17 -
    1.18 -#if GLP_MAJOR_VERSION > 4 || (GLP_MAJOR_VERSION == 4 && GLP_MINOR_VERSION > 15)
    1.19 -#define LEMON_glp(func) (glp_##func)
    1.20 -#define LEMON_lpx(func) (lpx_##func)
    1.21 -
    1.22 -#define LEMON_GLP(def) (GLP_##def)
    1.23 -#define LEMON_LPX(def) (LPX_##def)
    1.24 -
    1.25 -#else
    1.26 -
    1.27 -#define LEMON_glp(func) (lpx_##func)
    1.28 -#define LEMON_lpx(func) (lpx_##func)
    1.29 -
    1.30 -#define LEMON_GLP(def) (LPX_##def)
    1.31 -#define LEMON_LPX(def) (LPX_##def)
    1.32 -
    1.33 -#endif
    1.34 +#include <lemon/assert.h>
    1.35  
    1.36  namespace lemon {
    1.37  
    1.38 -  LpGlpk::LpGlpk() : Parent() {
    1.39 -    solved = false;
    1.40 -    rows = _lp_bits::LpId(1);
    1.41 -    cols = _lp_bits::LpId(1);
    1.42 -    lp = LEMON_glp(create_prob)();
    1.43 -    LEMON_glp(create_index)(lp);
    1.44 -    messageLevel(0);
    1.45 +  // GlpkBase members
    1.46 +
    1.47 +  GlpkBase::GlpkBase() : LpBase() {
    1.48 +    lp = glp_create_prob();
    1.49 +    glp_create_index(lp);
    1.50    }
    1.51  
    1.52 -  LpGlpk::LpGlpk(const LpGlpk &glp) : Parent() {
    1.53 -    solved = false;
    1.54 -    rows = _lp_bits::LpId(1);
    1.55 -    cols = _lp_bits::LpId(1);
    1.56 -    lp = LEMON_glp(create_prob)();
    1.57 -    LEMON_glp(create_index)(lp);
    1.58 -    messageLevel(0);
    1.59 -    //Coefficient matrix, row bounds
    1.60 -    LEMON_glp(add_rows)(lp, LEMON_glp(get_num_rows)(glp.lp));
    1.61 -    LEMON_glp(add_cols)(lp, LEMON_glp(get_num_cols)(glp.lp));
    1.62 -    int len;
    1.63 -    std::vector<int> ind(1+LEMON_glp(get_num_cols)(glp.lp));
    1.64 -    std::vector<Value> val(1+LEMON_glp(get_num_cols)(glp.lp));
    1.65 -    for (int i=1;i<=LEMON_glp(get_num_rows)(glp.lp);++i)
    1.66 -      {
    1.67 -        len=LEMON_glp(get_mat_row)(glp.lp,i,&*ind.begin(),&*val.begin());
    1.68 -        LEMON_glp(set_mat_row)(lp, i,len,&*ind.begin(),&*val.begin());
    1.69 -        LEMON_glp(set_row_bnds)(lp,i,
    1.70 -                                LEMON_glp(get_row_type)(glp.lp,i),
    1.71 -                                LEMON_glp(get_row_lb)(glp.lp,i),
    1.72 -                                LEMON_glp(get_row_ub)(glp.lp,i));
    1.73 -      }
    1.74 -
    1.75 -    //Objective function, coloumn bounds
    1.76 -    LEMON_glp(set_obj_dir)(lp, LEMON_glp(get_obj_dir)(glp.lp));
    1.77 -    //Objectif function's constant term treated separately
    1.78 -    LEMON_glp(set_obj_coef)(lp,0,LEMON_glp(get_obj_coef)(glp.lp,0));
    1.79 -    for (int i=1;i<=LEMON_glp(get_num_cols)(glp.lp);++i)
    1.80 -      {
    1.81 -        LEMON_glp(set_obj_coef)(lp,i,
    1.82 -                                LEMON_glp(get_obj_coef)(glp.lp,i));
    1.83 -        LEMON_glp(set_col_bnds)(lp,i,
    1.84 -                                LEMON_glp(get_col_type)(glp.lp,i),
    1.85 -                                LEMON_glp(get_col_lb)(glp.lp,i),
    1.86 -                                LEMON_glp(get_col_ub)(glp.lp,i));
    1.87 -      }
    1.88 -    rows = glp.rows;
    1.89 -    cols = glp.cols;
    1.90 +  GlpkBase::GlpkBase(const GlpkBase &other) : LpBase() {
    1.91 +    lp = glp_create_prob();
    1.92 +    glp_copy_prob(lp, other.lp, GLP_ON);
    1.93 +    glp_create_index(lp);
    1.94 +    rows = other.rows;
    1.95 +    cols = other.cols;
    1.96    }
    1.97  
    1.98 -  LpGlpk::~LpGlpk() {
    1.99 -    LEMON_glp(delete_prob)(lp);
   1.100 +  GlpkBase::~GlpkBase() {
   1.101 +    glp_delete_prob(lp);
   1.102    }
   1.103  
   1.104 -  int LpGlpk::_addCol() {
   1.105 -    int i=LEMON_glp(add_cols)(lp, 1);
   1.106 -    LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(FR), 0.0, 0.0);
   1.107 -    solved = false;
   1.108 +  int GlpkBase::_addCol() {
   1.109 +    int i = glp_add_cols(lp, 1);
   1.110 +    glp_set_col_bnds(lp, i, GLP_FR, 0.0, 0.0);
   1.111      return i;
   1.112    }
   1.113  
   1.114 -  ///\e
   1.115 -
   1.116 -
   1.117 -  LpSolverBase* LpGlpk::_newLp()
   1.118 -  {
   1.119 -    LpGlpk* newlp = new LpGlpk;
   1.120 -    return newlp;
   1.121 -  }
   1.122 -
   1.123 -  ///\e
   1.124 -
   1.125 -  LpSolverBase* LpGlpk::_copyLp()
   1.126 -  {
   1.127 -    LpGlpk *newlp = new LpGlpk(*this);
   1.128 -    return newlp;
   1.129 -  }
   1.130 -
   1.131 -  int LpGlpk::_addRow() {
   1.132 -    int i=LEMON_glp(add_rows)(lp, 1);
   1.133 -    solved = false;
   1.134 +  int GlpkBase::_addRow() {
   1.135 +    int i = glp_add_rows(lp, 1);
   1.136 +    glp_set_row_bnds(lp, i, GLP_FR, 0.0, 0.0);
   1.137      return i;
   1.138    }
   1.139  
   1.140 -
   1.141 -  void LpGlpk::_eraseCol(int i) {
   1.142 +  void GlpkBase::_eraseCol(int i) {
   1.143      int ca[2];
   1.144 -    ca[1]=i;
   1.145 -    LEMON_glp(del_cols)(lp, 1, ca);
   1.146 -    solved = false;
   1.147 +    ca[1] = i;
   1.148 +    glp_del_cols(lp, 1, ca);
   1.149    }
   1.150  
   1.151 -  void LpGlpk::_eraseRow(int i) {
   1.152 +  void GlpkBase::_eraseRow(int i) {
   1.153      int ra[2];
   1.154 -    ra[1]=i;
   1.155 -    LEMON_glp(del_rows)(lp, 1, ra);
   1.156 -    solved = false;
   1.157 +    ra[1] = i;
   1.158 +    glp_del_rows(lp, 1, ra);
   1.159    }
   1.160  
   1.161 -  void LpGlpk::_getColName(int c, std::string & name) const
   1.162 -  {
   1.163 -
   1.164 -    const char *n = LEMON_glp(get_col_name)(lp,c);
   1.165 -    name = n?n:"";
   1.166 +  void GlpkBase::_eraseColId(int i) {
   1.167 +    cols.eraseIndex(i);
   1.168 +    cols.shiftIndices(i);
   1.169    }
   1.170  
   1.171 +  void GlpkBase::_eraseRowId(int i) {
   1.172 +    rows.eraseIndex(i);
   1.173 +    rows.shiftIndices(i);
   1.174 +  }
   1.175  
   1.176 -  void LpGlpk::_setColName(int c, const std::string & name)
   1.177 -  {
   1.178 -    LEMON_glp(set_col_name)(lp,c,const_cast<char*>(name.c_str()));
   1.179 +  void GlpkBase::_getColName(int c, std::string& name) const {
   1.180 +    const char *str = glp_get_col_name(lp, c);
   1.181 +    if (str) name = str;
   1.182 +    else name.clear();
   1.183 +  }
   1.184 +
   1.185 +  void GlpkBase::_setColName(int c, const std::string & name) {
   1.186 +    glp_set_col_name(lp, c, const_cast<char*>(name.c_str()));
   1.187  
   1.188    }
   1.189  
   1.190 -  int LpGlpk::_colByName(const std::string& name) const
   1.191 -  {
   1.192 -    int k = LEMON_glp(find_col)(lp, const_cast<char*>(name.c_str()));
   1.193 +  int GlpkBase::_colByName(const std::string& name) const {
   1.194 +    int k = glp_find_col(lp, const_cast<char*>(name.c_str()));
   1.195      return k > 0 ? k : -1;
   1.196    }
   1.197  
   1.198 +  void GlpkBase::_getRowName(int r, std::string& name) const {
   1.199 +    const char *str = glp_get_row_name(lp, r);
   1.200 +    if (str) name = str;
   1.201 +    else name.clear();
   1.202 +  }
   1.203  
   1.204 -  void LpGlpk::_setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e)
   1.205 -  {
   1.206 -    std::vector<int> indices;
   1.207 +  void GlpkBase::_setRowName(int r, const std::string & name) {
   1.208 +    glp_set_row_name(lp, r, const_cast<char*>(name.c_str()));
   1.209 +
   1.210 +  }
   1.211 +
   1.212 +  int GlpkBase::_rowByName(const std::string& name) const {
   1.213 +    int k = glp_find_row(lp, const_cast<char*>(name.c_str()));
   1.214 +    return k > 0 ? k : -1;
   1.215 +  }
   1.216 +
   1.217 +  void GlpkBase::_setRowCoeffs(int i, ExprIterator b, ExprIterator e) {
   1.218 +    std::vector<int> indexes;
   1.219      std::vector<Value> values;
   1.220  
   1.221 -    indices.push_back(0);
   1.222 +    indexes.push_back(0);
   1.223      values.push_back(0);
   1.224  
   1.225 -    for(ConstRowIterator it=b; it!=e; ++it) {
   1.226 -      indices.push_back(it->first);
   1.227 +    for(ExprIterator it = b; it != e; ++it) {
   1.228 +      indexes.push_back(it->first);
   1.229        values.push_back(it->second);
   1.230      }
   1.231  
   1.232 -    LEMON_glp(set_mat_row)(lp, i, values.size() - 1,
   1.233 -                                &indices[0], &values[0]);
   1.234 -
   1.235 -    solved = false;
   1.236 +    glp_set_mat_row(lp, i, values.size() - 1,
   1.237 +                    &indexes.front(), &values.front());
   1.238    }
   1.239  
   1.240 -  void LpGlpk::_getRowCoeffs(int ix, RowIterator b) const
   1.241 -  {
   1.242 -    int length = LEMON_glp(get_mat_row)(lp, ix, 0, 0);
   1.243 +  void GlpkBase::_getRowCoeffs(int ix, InsertIterator b) const {
   1.244 +    int length = glp_get_mat_row(lp, ix, 0, 0);
   1.245  
   1.246 -    std::vector<int> indices(length + 1);
   1.247 +    std::vector<int> indexes(length + 1);
   1.248      std::vector<Value> values(length + 1);
   1.249  
   1.250 -    LEMON_glp(get_mat_row)(lp, ix, &indices[0], &values[0]);
   1.251 +    glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
   1.252  
   1.253      for (int i = 1; i <= length; ++i) {
   1.254 -      *b = std::make_pair(indices[i], values[i]);
   1.255 +      *b = std::make_pair(indexes[i], values[i]);
   1.256        ++b;
   1.257      }
   1.258    }
   1.259  
   1.260 -  void LpGlpk::_setColCoeffs(int ix, ConstColIterator b, ConstColIterator e) {
   1.261 +  void GlpkBase::_setColCoeffs(int ix, ExprIterator b,
   1.262 +                                     ExprIterator e) {
   1.263  
   1.264 -    std::vector<int> indices;
   1.265 +    std::vector<int> indexes;
   1.266      std::vector<Value> values;
   1.267  
   1.268 -    indices.push_back(0);
   1.269 +    indexes.push_back(0);
   1.270      values.push_back(0);
   1.271  
   1.272 -    for(ConstColIterator it=b; it!=e; ++it) {
   1.273 -      indices.push_back(it->first);
   1.274 +    for(ExprIterator it = b; it != e; ++it) {
   1.275 +      indexes.push_back(it->first);
   1.276        values.push_back(it->second);
   1.277      }
   1.278  
   1.279 -    LEMON_glp(set_mat_col)(lp, ix, values.size() - 1,
   1.280 -                                &indices[0], &values[0]);
   1.281 -
   1.282 -    solved = false;
   1.283 +    glp_set_mat_col(lp, ix, values.size() - 1,
   1.284 +                    &indexes.front(), &values.front());
   1.285    }
   1.286  
   1.287 -  void LpGlpk::_getColCoeffs(int ix, ColIterator b) const
   1.288 -  {
   1.289 -    int length = LEMON_glp(get_mat_col)(lp, ix, 0, 0);
   1.290 +  void GlpkBase::_getColCoeffs(int ix, InsertIterator b) const {
   1.291 +    int length = glp_get_mat_col(lp, ix, 0, 0);
   1.292  
   1.293 -    std::vector<int> indices(length + 1);
   1.294 +    std::vector<int> indexes(length + 1);
   1.295      std::vector<Value> values(length + 1);
   1.296  
   1.297 -    LEMON_glp(get_mat_col)(lp, ix, &indices[0], &values[0]);
   1.298 +    glp_get_mat_col(lp, ix, &indexes.front(), &values.front());
   1.299  
   1.300 -    for (int i = 1; i <= length; ++i) {
   1.301 -      *b = std::make_pair(indices[i], values[i]);
   1.302 +    for (int i = 1; i  <= length; ++i) {
   1.303 +      *b = std::make_pair(indexes[i], values[i]);
   1.304        ++b;
   1.305      }
   1.306    }
   1.307  
   1.308 -  void LpGlpk::_setCoeff(int ix, int jx, Value value)
   1.309 -  {
   1.310 +  void GlpkBase::_setCoeff(int ix, int jx, Value value) {
   1.311  
   1.312 -    if (LEMON_glp(get_num_cols)(lp) < LEMON_glp(get_num_rows)(lp)) {
   1.313 +    if (glp_get_num_cols(lp) < glp_get_num_rows(lp)) {
   1.314  
   1.315 -      int length=LEMON_glp(get_mat_row)(lp, ix, 0, 0);
   1.316 +      int length = glp_get_mat_row(lp, ix, 0, 0);
   1.317  
   1.318 -      std::vector<int> indices(length + 2);
   1.319 +      std::vector<int> indexes(length + 2);
   1.320        std::vector<Value> values(length + 2);
   1.321  
   1.322 -      LEMON_glp(get_mat_row)(lp, ix, &indices[0], &values[0]);
   1.323 +      glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
   1.324  
   1.325        //The following code does not suppose that the elements of the
   1.326 -      //array indices are sorted
   1.327 -      bool found=false;
   1.328 -      for (int i = 1; i <= length; ++i) {
   1.329 -        if (indices[i]==jx){
   1.330 -          found=true;
   1.331 -          values[i]=value;
   1.332 +      //array indexes are sorted
   1.333 +      bool found = false;
   1.334 +      for (int i = 1; i  <= length; ++i) {
   1.335 +        if (indexes[i] == jx) {
   1.336 +          found = true;
   1.337 +          values[i] = value;
   1.338            break;
   1.339          }
   1.340        }
   1.341 -      if (!found){
   1.342 +      if (!found) {
   1.343          ++length;
   1.344 -        indices[length]=jx;
   1.345 -        values[length]=value;
   1.346 +        indexes[length] = jx;
   1.347 +        values[length] = value;
   1.348        }
   1.349  
   1.350 -      LEMON_glp(set_mat_row)(lp, ix, length, &indices[0], &values[0]);
   1.351 +      glp_set_mat_row(lp, ix, length, &indexes.front(), &values.front());
   1.352  
   1.353      } else {
   1.354  
   1.355 -      int length=LEMON_glp(get_mat_col)(lp, jx, 0, 0);
   1.356 +      int length = glp_get_mat_col(lp, jx, 0, 0);
   1.357  
   1.358 -      std::vector<int> indices(length + 2);
   1.359 +      std::vector<int> indexes(length + 2);
   1.360        std::vector<Value> values(length + 2);
   1.361  
   1.362 -      LEMON_glp(get_mat_col)(lp, jx, &indices[0], &values[0]);
   1.363 +      glp_get_mat_col(lp, jx, &indexes.front(), &values.front());
   1.364  
   1.365        //The following code does not suppose that the elements of the
   1.366 -      //array indices are sorted
   1.367 -      bool found=false;
   1.368 +      //array indexes are sorted
   1.369 +      bool found = false;
   1.370        for (int i = 1; i <= length; ++i) {
   1.371 -        if (indices[i]==ix){
   1.372 -          found=true;
   1.373 -          values[i]=value;
   1.374 +        if (indexes[i] == ix) {
   1.375 +          found = true;
   1.376 +          values[i] = value;
   1.377            break;
   1.378          }
   1.379        }
   1.380 -      if (!found){
   1.381 +      if (!found) {
   1.382          ++length;
   1.383 -        indices[length]=ix;
   1.384 -        values[length]=value;
   1.385 +        indexes[length] = ix;
   1.386 +        values[length] = value;
   1.387        }
   1.388  
   1.389 -      LEMON_glp(set_mat_col)(lp, jx, length, &indices[0], &values[0]);
   1.390 +      glp_set_mat_col(lp, jx, length, &indexes.front(), &values.front());
   1.391      }
   1.392  
   1.393 -    solved = false;
   1.394    }
   1.395  
   1.396 -  LpGlpk::Value LpGlpk::_getCoeff(int ix, int jx) const
   1.397 -  {
   1.398 +  GlpkBase::Value GlpkBase::_getCoeff(int ix, int jx) const {
   1.399  
   1.400 -    int length=LEMON_glp(get_mat_row)(lp, ix, 0, 0);
   1.401 +    int length = glp_get_mat_row(lp, ix, 0, 0);
   1.402  
   1.403 -    std::vector<int> indices(length + 1);
   1.404 +    std::vector<int> indexes(length + 1);
   1.405      std::vector<Value> values(length + 1);
   1.406  
   1.407 -    LEMON_glp(get_mat_row)(lp, ix, &indices[0], &values[0]);
   1.408 +    glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
   1.409  
   1.410 -    //The following code does not suppose that the elements of the
   1.411 -    //array indices are sorted
   1.412 -    for (int i = 1; i <= length; ++i) {
   1.413 -      if (indices[i]==jx){
   1.414 +    for (int i = 1; i  <= length; ++i) {
   1.415 +      if (indexes[i] == jx) {
   1.416          return values[i];
   1.417        }
   1.418      }
   1.419 +
   1.420      return 0;
   1.421 +  }
   1.422 +
   1.423 +  void GlpkBase::_setColLowerBound(int i, Value lo) {
   1.424 +    LEMON_ASSERT(lo != INF, "Invalid bound");
   1.425 +
   1.426 +    int b = glp_get_col_type(lp, i);
   1.427 +    double up = glp_get_col_ub(lp, i);
   1.428 +    if (lo == -INF) {
   1.429 +      switch (b) {
   1.430 +      case GLP_FR:
   1.431 +      case GLP_LO:
   1.432 +        glp_set_col_bnds(lp, i, GLP_FR, lo, up);
   1.433 +        break;
   1.434 +      case GLP_UP:
   1.435 +        break;
   1.436 +      case GLP_DB:
   1.437 +      case GLP_FX:
   1.438 +        glp_set_col_bnds(lp, i, GLP_UP, lo, up);
   1.439 +        break;
   1.440 +      default:
   1.441 +        break;
   1.442 +      }
   1.443 +    } else {
   1.444 +      switch (b) {
   1.445 +      case GLP_FR:
   1.446 +      case GLP_LO:
   1.447 +        glp_set_col_bnds(lp, i, GLP_LO, lo, up);
   1.448 +        break;
   1.449 +      case GLP_UP:
   1.450 +      case GLP_DB:
   1.451 +      case GLP_FX:
   1.452 +        if (lo == up)
   1.453 +          glp_set_col_bnds(lp, i, GLP_FX, lo, up);
   1.454 +        else
   1.455 +          glp_set_col_bnds(lp, i, GLP_DB, lo, up);
   1.456 +        break;
   1.457 +      default:
   1.458 +        break;
   1.459 +      }
   1.460 +    }
   1.461 +  }
   1.462 +
   1.463 +  GlpkBase::Value GlpkBase::_getColLowerBound(int i) const {
   1.464 +    int b = glp_get_col_type(lp, i);
   1.465 +    switch (b) {
   1.466 +    case GLP_LO:
   1.467 +    case GLP_DB:
   1.468 +    case GLP_FX:
   1.469 +      return glp_get_col_lb(lp, i);
   1.470 +    default:
   1.471 +      return -INF;
   1.472 +    }
   1.473 +  }
   1.474 +
   1.475 +  void GlpkBase::_setColUpperBound(int i, Value up) {
   1.476 +    LEMON_ASSERT(up != -INF, "Invalid bound");
   1.477 +
   1.478 +    int b = glp_get_col_type(lp, i);
   1.479 +    double lo = glp_get_col_lb(lp, i);
   1.480 +    if (up == INF) {
   1.481 +      switch (b) {
   1.482 +      case GLP_FR:
   1.483 +      case GLP_LO:
   1.484 +        break;
   1.485 +      case GLP_UP:
   1.486 +        glp_set_col_bnds(lp, i, GLP_FR, lo, up);
   1.487 +        break;
   1.488 +      case GLP_DB:
   1.489 +      case GLP_FX:
   1.490 +        glp_set_col_bnds(lp, i, GLP_LO, lo, up);
   1.491 +        break;
   1.492 +      default:
   1.493 +        break;
   1.494 +      }
   1.495 +    } else {
   1.496 +      switch (b) {
   1.497 +      case GLP_FR:
   1.498 +        glp_set_col_bnds(lp, i, GLP_UP, lo, up);
   1.499 +        break;
   1.500 +      case GLP_UP:
   1.501 +        glp_set_col_bnds(lp, i, GLP_UP, lo, up);
   1.502 +        break;
   1.503 +      case GLP_LO:
   1.504 +      case GLP_DB:
   1.505 +      case GLP_FX:
   1.506 +        if (lo == up)
   1.507 +          glp_set_col_bnds(lp, i, GLP_FX, lo, up);
   1.508 +        else
   1.509 +          glp_set_col_bnds(lp, i, GLP_DB, lo, up);
   1.510 +        break;
   1.511 +      default:
   1.512 +        break;
   1.513 +      }
   1.514 +    }
   1.515  
   1.516    }
   1.517  
   1.518 -
   1.519 -  void LpGlpk::_setColLowerBound(int i, Value lo)
   1.520 -  {
   1.521 -    if (lo==INF) {
   1.522 -      //FIXME error
   1.523 -    }
   1.524 -    int b=LEMON_glp(get_col_type)(lp, i);
   1.525 -    double up=LEMON_glp(get_col_ub)(lp, i);
   1.526 -    if (lo==-INF) {
   1.527 +  GlpkBase::Value GlpkBase::_getColUpperBound(int i) const {
   1.528 +    int b = glp_get_col_type(lp, i);
   1.529        switch (b) {
   1.530 -      case LEMON_GLP(FR):
   1.531 -      case LEMON_GLP(LO):
   1.532 -        LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(FR), lo, up);
   1.533 -        break;
   1.534 -      case LEMON_GLP(UP):
   1.535 -        break;
   1.536 -      case LEMON_GLP(DB):
   1.537 -      case LEMON_GLP(FX):
   1.538 -        LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(UP), lo, up);
   1.539 -        break;
   1.540 -      default: ;
   1.541 -        //FIXME error
   1.542 -      }
   1.543 -    } else {
   1.544 -      switch (b) {
   1.545 -      case LEMON_GLP(FR):
   1.546 -      case LEMON_GLP(LO):
   1.547 -        LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(LO), lo, up);
   1.548 -        break;
   1.549 -      case LEMON_GLP(UP):
   1.550 -      case LEMON_GLP(DB):
   1.551 -      case LEMON_GLP(FX):
   1.552 -        if (lo==up)
   1.553 -          LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(FX), lo, up);
   1.554 -        else
   1.555 -          LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(DB), lo, up);
   1.556 -        break;
   1.557 -      default: ;
   1.558 -        //FIXME error
   1.559 -      }
   1.560 -    }
   1.561 -
   1.562 -    solved = false;
   1.563 -  }
   1.564 -
   1.565 -  LpGlpk::Value LpGlpk::_getColLowerBound(int i) const
   1.566 -  {
   1.567 -    int b=LEMON_glp(get_col_type)(lp, i);
   1.568 -      switch (b) {
   1.569 -      case LEMON_GLP(LO):
   1.570 -      case LEMON_GLP(DB):
   1.571 -      case LEMON_GLP(FX):
   1.572 -        return LEMON_glp(get_col_lb)(lp, i);
   1.573 -      default: ;
   1.574 -        return -INF;
   1.575 -      }
   1.576 -  }
   1.577 -
   1.578 -  void LpGlpk::_setColUpperBound(int i, Value up)
   1.579 -  {
   1.580 -    if (up==-INF) {
   1.581 -      //FIXME error
   1.582 -    }
   1.583 -    int b=LEMON_glp(get_col_type)(lp, i);
   1.584 -    double lo=LEMON_glp(get_col_lb)(lp, i);
   1.585 -    if (up==INF) {
   1.586 -      switch (b) {
   1.587 -      case LEMON_GLP(FR):
   1.588 -      case LEMON_GLP(LO):
   1.589 -        break;
   1.590 -      case LEMON_GLP(UP):
   1.591 -        LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(FR), lo, up);
   1.592 -        break;
   1.593 -      case LEMON_GLP(DB):
   1.594 -      case LEMON_GLP(FX):
   1.595 -        LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(LO), lo, up);
   1.596 -        break;
   1.597 -      default: ;
   1.598 -        //FIXME error
   1.599 -      }
   1.600 -    } else {
   1.601 -      switch (b) {
   1.602 -      case LEMON_GLP(FR):
   1.603 -        LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(UP), lo, up);
   1.604 -        break;
   1.605 -      case LEMON_GLP(UP):
   1.606 -        LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(UP), lo, up);
   1.607 -        break;
   1.608 -      case LEMON_GLP(LO):
   1.609 -      case LEMON_GLP(DB):
   1.610 -      case LEMON_GLP(FX):
   1.611 -        if (lo==up)
   1.612 -          LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(FX), lo, up);
   1.613 -        else
   1.614 -          LEMON_glp(set_col_bnds)(lp, i, LEMON_GLP(DB), lo, up);
   1.615 -        break;
   1.616 -      default: ;
   1.617 -        //FIXME error
   1.618 -      }
   1.619 -    }
   1.620 -
   1.621 -    solved = false;
   1.622 -  }
   1.623 -
   1.624 -  LpGlpk::Value LpGlpk::_getColUpperBound(int i) const
   1.625 -  {
   1.626 -    int b=LEMON_glp(get_col_type)(lp, i);
   1.627 -      switch (b) {
   1.628 -      case LEMON_GLP(UP):
   1.629 -      case LEMON_GLP(DB):
   1.630 -      case LEMON_GLP(FX):
   1.631 -        return LEMON_glp(get_col_ub)(lp, i);
   1.632 -      default: ;
   1.633 +      case GLP_UP:
   1.634 +      case GLP_DB:
   1.635 +      case GLP_FX:
   1.636 +        return glp_get_col_ub(lp, i);
   1.637 +      default:
   1.638          return INF;
   1.639        }
   1.640    }
   1.641  
   1.642 -  void LpGlpk::_setRowBounds(int i, Value lb, Value ub)
   1.643 -  {
   1.644 -    //Bad parameter
   1.645 -    if (lb==INF || ub==-INF) {
   1.646 -      //FIXME error
   1.647 -    }
   1.648 +  void GlpkBase::_setRowLowerBound(int i, Value lo) {
   1.649 +    LEMON_ASSERT(lo != INF, "Invalid bound");
   1.650  
   1.651 -    if (lb == -INF){
   1.652 -      if (ub == INF){
   1.653 -        LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(FR), lb, ub);
   1.654 +    int b = glp_get_row_type(lp, i);
   1.655 +    double up = glp_get_row_ub(lp, i);
   1.656 +    if (lo == -INF) {
   1.657 +      switch (b) {
   1.658 +      case GLP_FR:
   1.659 +      case GLP_LO:
   1.660 +        glp_set_row_bnds(lp, i, GLP_FR, lo, up);
   1.661 +        break;
   1.662 +      case GLP_UP:
   1.663 +        break;
   1.664 +      case GLP_DB:
   1.665 +      case GLP_FX:
   1.666 +        glp_set_row_bnds(lp, i, GLP_UP, lo, up);
   1.667 +        break;
   1.668 +      default:
   1.669 +        break;
   1.670        }
   1.671 -      else{
   1.672 -        LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(UP), lb, ub);
   1.673 +    } else {
   1.674 +      switch (b) {
   1.675 +      case GLP_FR:
   1.676 +      case GLP_LO:
   1.677 +        glp_set_row_bnds(lp, i, GLP_LO, lo, up);
   1.678 +        break;
   1.679 +      case GLP_UP:
   1.680 +      case GLP_DB:
   1.681 +      case GLP_FX:
   1.682 +        if (lo == up)
   1.683 +          glp_set_row_bnds(lp, i, GLP_FX, lo, up);
   1.684 +        else
   1.685 +          glp_set_row_bnds(lp, i, GLP_DB, lo, up);
   1.686 +        break;
   1.687 +      default:
   1.688 +        break;
   1.689        }
   1.690      }
   1.691 -    else{
   1.692 -      if (ub==INF){
   1.693 -        LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(LO), lb, ub);
   1.694 -
   1.695 -      }
   1.696 -      else{
   1.697 -        if (lb == ub){
   1.698 -          LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(FX), lb, ub);
   1.699 -        }
   1.700 -        else{
   1.701 -          LEMON_glp(set_row_bnds)(lp, i, LEMON_GLP(DB), lb, ub);
   1.702 -        }
   1.703 -      }
   1.704 -    }
   1.705 -
   1.706 -    solved = false;
   1.707 -  }
   1.708 -
   1.709 -  void LpGlpk::_getRowBounds(int i, Value &lb, Value &ub) const
   1.710 -  {
   1.711 -
   1.712 -    int b=LEMON_glp(get_row_type)(lp, i);
   1.713 -    switch (b) {
   1.714 -    case LEMON_GLP(FR):
   1.715 -    case LEMON_GLP(UP):
   1.716 -      lb = -INF;
   1.717 -        break;
   1.718 -    default:
   1.719 -      lb=LEMON_glp(get_row_lb)(lp, i);
   1.720 -    }
   1.721 -
   1.722 -    switch (b) {
   1.723 -    case LEMON_GLP(FR):
   1.724 -    case LEMON_GLP(LO):
   1.725 -      ub = INF;
   1.726 -        break;
   1.727 -    default:
   1.728 -      ub=LEMON_glp(get_row_ub)(lp, i);
   1.729 -    }
   1.730  
   1.731    }
   1.732  
   1.733 -  void LpGlpk::_setObjCoeff(int i, Value obj_coef)
   1.734 -  {
   1.735 -    //i=0 means the constant term (shift)
   1.736 -    LEMON_glp(set_obj_coef)(lp, i, obj_coef);
   1.737 -
   1.738 -    solved = false;
   1.739 -  }
   1.740 -
   1.741 -  LpGlpk::Value LpGlpk::_getObjCoeff(int i) const {
   1.742 -    //i=0 means the constant term (shift)
   1.743 -    return LEMON_glp(get_obj_coef)(lp, i);
   1.744 -  }
   1.745 -
   1.746 -  void LpGlpk::_clearObj()
   1.747 -  {
   1.748 -    for (int i=0;i<=LEMON_glp(get_num_cols)(lp);++i){
   1.749 -      LEMON_glp(set_obj_coef)(lp, i, 0);
   1.750 -    }
   1.751 -
   1.752 -    solved = false;
   1.753 -  }
   1.754 -
   1.755 -  LpGlpk::SolveExitStatus LpGlpk::_solve()
   1.756 -  {
   1.757 -    // A way to check the problem to be solved
   1.758 -    //LEMON_glp(write_cpxlp(lp,"naittvan.cpx");
   1.759 -
   1.760 -    LEMON_lpx(std_basis)(lp);
   1.761 -    int i =  LEMON_lpx(simplex)(lp);
   1.762 -
   1.763 -    switch (i) {
   1.764 -    case LEMON_LPX(E_OK):
   1.765 -      solved = true;
   1.766 -      return SOLVED;
   1.767 +  GlpkBase::Value GlpkBase::_getRowLowerBound(int i) const {
   1.768 +    int b = glp_get_row_type(lp, i);
   1.769 +    switch (b) {
   1.770 +    case GLP_LO:
   1.771 +    case GLP_DB:
   1.772 +    case GLP_FX:
   1.773 +      return glp_get_row_lb(lp, i);
   1.774      default:
   1.775 -      return UNSOLVED;
   1.776 +      return -INF;
   1.777      }
   1.778    }
   1.779  
   1.780 -  LpGlpk::Value LpGlpk::_getPrimal(int i) const
   1.781 -  {
   1.782 -    return LEMON_glp(get_col_prim)(lp,i);
   1.783 -  }
   1.784 +  void GlpkBase::_setRowUpperBound(int i, Value up) {
   1.785 +    LEMON_ASSERT(up != -INF, "Invalid bound");
   1.786  
   1.787 -  LpGlpk::Value LpGlpk::_getDual(int i) const
   1.788 -  {
   1.789 -    return LEMON_glp(get_row_dual)(lp,i);
   1.790 -  }
   1.791 -
   1.792 -  LpGlpk::Value LpGlpk::_getPrimalValue() const
   1.793 -  {
   1.794 -    return LEMON_glp(get_obj_val)(lp);
   1.795 -  }
   1.796 -  bool LpGlpk::_isBasicCol(int i) const
   1.797 -  {
   1.798 -    return (LEMON_glp(get_col_stat)(lp, i)==LEMON_GLP(BS));
   1.799 -  }
   1.800 -
   1.801 -
   1.802 -  LpGlpk::SolutionStatus LpGlpk::_getPrimalStatus() const
   1.803 -  {
   1.804 -    if (!solved) return UNDEFINED;
   1.805 -    int stat=  LEMON_lpx(get_status)(lp);
   1.806 -    switch (stat) {
   1.807 -    case LEMON_LPX(UNDEF)://Undefined (no solve has been run yet)
   1.808 -      return UNDEFINED;
   1.809 -    case LEMON_LPX(NOFEAS)://There is no feasible solution (primal, I guess)
   1.810 -    case LEMON_LPX(INFEAS)://Infeasible
   1.811 -      return INFEASIBLE;
   1.812 -    case LEMON_LPX(UNBND)://Unbounded
   1.813 -      return INFINITE;
   1.814 -    case LEMON_LPX(FEAS)://Feasible
   1.815 -      return FEASIBLE;
   1.816 -    case LEMON_LPX(OPT)://Feasible
   1.817 -      return OPTIMAL;
   1.818 -    default:
   1.819 -      return UNDEFINED; //to avoid gcc warning
   1.820 -      //FIXME error
   1.821 +    int b = glp_get_row_type(lp, i);
   1.822 +    double lo = glp_get_row_lb(lp, i);
   1.823 +    if (up == INF) {
   1.824 +      switch (b) {
   1.825 +      case GLP_FR:
   1.826 +      case GLP_LO:
   1.827 +        break;
   1.828 +      case GLP_UP:
   1.829 +        glp_set_row_bnds(lp, i, GLP_FR, lo, up);
   1.830 +        break;
   1.831 +      case GLP_DB:
   1.832 +      case GLP_FX:
   1.833 +        glp_set_row_bnds(lp, i, GLP_LO, lo, up);
   1.834 +        break;
   1.835 +      default:
   1.836 +        break;
   1.837 +      }
   1.838 +    } else {
   1.839 +      switch (b) {
   1.840 +      case GLP_FR:
   1.841 +        glp_set_row_bnds(lp, i, GLP_UP, lo, up);
   1.842 +        break;
   1.843 +      case GLP_UP:
   1.844 +        glp_set_row_bnds(lp, i, GLP_UP, lo, up);
   1.845 +        break;
   1.846 +      case GLP_LO:
   1.847 +      case GLP_DB:
   1.848 +      case GLP_FX:
   1.849 +        if (lo == up)
   1.850 +          glp_set_row_bnds(lp, i, GLP_FX, lo, up);
   1.851 +        else
   1.852 +          glp_set_row_bnds(lp, i, GLP_DB, lo, up);
   1.853 +        break;
   1.854 +      default:
   1.855 +        break;
   1.856 +      }
   1.857      }
   1.858    }
   1.859  
   1.860 -  LpGlpk::SolutionStatus LpGlpk::_getDualStatus() const
   1.861 -  {
   1.862 -    if (!solved) return UNDEFINED;
   1.863 -    switch (LEMON_lpx(get_dual_stat)(lp)) {
   1.864 -    case LEMON_LPX(D_UNDEF)://Undefined (no solve has been run yet)
   1.865 -      return UNDEFINED;
   1.866 -    case LEMON_LPX(D_NOFEAS)://There is no dual feasible solution
   1.867 -//    case LEMON_LPX(D_INFEAS://Infeasible
   1.868 -      return INFEASIBLE;
   1.869 -    case LEMON_LPX(D_FEAS)://Feasible
   1.870 -      switch (LEMON_lpx(get_status)(lp)) {
   1.871 -      case LEMON_LPX(NOFEAS):
   1.872 -        return INFINITE;
   1.873 -      case LEMON_LPX(OPT):
   1.874 -        return OPTIMAL;
   1.875 -      default:
   1.876 -        return FEASIBLE;
   1.877 -      }
   1.878 +  GlpkBase::Value GlpkBase::_getRowUpperBound(int i) const {
   1.879 +    int b = glp_get_row_type(lp, i);
   1.880 +    switch (b) {
   1.881 +    case GLP_UP:
   1.882 +    case GLP_DB:
   1.883 +    case GLP_FX:
   1.884 +      return glp_get_row_ub(lp, i);
   1.885      default:
   1.886 -      return UNDEFINED; //to avoid gcc warning
   1.887 -      //FIXME error
   1.888 +      return INF;
   1.889      }
   1.890    }
   1.891  
   1.892 -  LpGlpk::ProblemTypes LpGlpk::_getProblemType() const
   1.893 -  {
   1.894 -    if (!solved) return UNKNOWN;
   1.895 -      //int stat=  LEMON_glp(get_status(lp);
   1.896 -    int statp=  LEMON_lpx(get_prim_stat)(lp);
   1.897 -    int statd=  LEMON_lpx(get_dual_stat)(lp);
   1.898 -    if (statp==LEMON_LPX(P_FEAS) && statd==LEMON_LPX(D_FEAS))
   1.899 -        return PRIMAL_DUAL_FEASIBLE;
   1.900 -    if (statp==LEMON_LPX(P_FEAS) && statd==LEMON_LPX(D_NOFEAS))
   1.901 -        return PRIMAL_FEASIBLE_DUAL_INFEASIBLE;
   1.902 -    if (statp==LEMON_LPX(P_NOFEAS) && statd==LEMON_LPX(D_FEAS))
   1.903 -        return PRIMAL_INFEASIBLE_DUAL_FEASIBLE;
   1.904 -    if (statp==LEMON_LPX(P_NOFEAS) && statd==LEMON_LPX(D_NOFEAS))
   1.905 -        return PRIMAL_DUAL_INFEASIBLE;
   1.906 -    //In all other cases
   1.907 -    return UNKNOWN;
   1.908 +  void GlpkBase::_setObjCoeffs(ExprIterator b, ExprIterator e) {
   1.909 +    for (int i = 1; i <= glp_get_num_cols(lp); ++i) {
   1.910 +      glp_set_obj_coef(lp, i, 0.0);
   1.911 +    }
   1.912 +    for (ExprIterator it = b; it != e; ++it) {
   1.913 +      glp_set_obj_coef(lp, it->first, it->second);
   1.914 +    }
   1.915    }
   1.916  
   1.917 -  void LpGlpk::_setMax()
   1.918 -  {
   1.919 -    solved = false;
   1.920 -    LEMON_glp(set_obj_dir)(lp, LEMON_GLP(MAX));
   1.921 +  void GlpkBase::_getObjCoeffs(InsertIterator b) const {
   1.922 +    for (int i = 1; i <= glp_get_num_cols(lp); ++i) {
   1.923 +      Value val = glp_get_obj_coef(lp, i);
   1.924 +      if (val != 0.0) {
   1.925 +        *b = std::make_pair(i, val);
   1.926 +        ++b;
   1.927 +      }
   1.928 +    }
   1.929    }
   1.930  
   1.931 -  void LpGlpk::_setMin()
   1.932 -  {
   1.933 -    solved = false;
   1.934 -    LEMON_glp(set_obj_dir)(lp, LEMON_GLP(MIN));
   1.935 +  void GlpkBase::_setObjCoeff(int i, Value obj_coef) {
   1.936 +    //i = 0 means the constant term (shift)
   1.937 +    glp_set_obj_coef(lp, i, obj_coef);
   1.938    }
   1.939  
   1.940 -  bool LpGlpk::_isMax() const
   1.941 -  {
   1.942 -    return (LEMON_glp(get_obj_dir)(lp)==LEMON_GLP(MAX));
   1.943 +  GlpkBase::Value GlpkBase::_getObjCoeff(int i) const {
   1.944 +    //i = 0 means the constant term (shift)
   1.945 +    return glp_get_obj_coef(lp, i);
   1.946    }
   1.947  
   1.948 -
   1.949 -
   1.950 -  void LpGlpk::messageLevel(int m)
   1.951 -  {
   1.952 -    LEMON_lpx(set_int_parm)(lp, LEMON_LPX(K_MSGLEV), m);
   1.953 +  void GlpkBase::_setSense(GlpkBase::Sense sense) {
   1.954 +    switch (sense) {
   1.955 +    case MIN:
   1.956 +      glp_set_obj_dir(lp, GLP_MIN);
   1.957 +      break;
   1.958 +    case MAX:
   1.959 +      glp_set_obj_dir(lp, GLP_MAX);
   1.960 +      break;
   1.961 +    }
   1.962    }
   1.963  
   1.964 -  void LpGlpk::presolver(bool b)
   1.965 -  {
   1.966 -    LEMON_lpx(set_int_parm)(lp, LEMON_LPX(K_PRESOL), b);
   1.967 +  GlpkBase::Sense GlpkBase::_getSense() const {
   1.968 +    switch(glp_get_obj_dir(lp)) {
   1.969 +    case GLP_MIN:
   1.970 +      return MIN;
   1.971 +    case GLP_MAX:
   1.972 +      return MAX;
   1.973 +    default:
   1.974 +      LEMON_ASSERT(false, "Wrong sense");
   1.975 +      return GlpkBase::Sense();
   1.976 +    }
   1.977    }
   1.978  
   1.979 +  void GlpkBase::_clear() {
   1.980 +    glp_erase_prob(lp);
   1.981 +    rows.clear();
   1.982 +    cols.clear();
   1.983 +  }
   1.984 +
   1.985 +  // LpGlpk members
   1.986 +
   1.987 +  LpGlpk::LpGlpk()
   1.988 +    : LpBase(), GlpkBase(), LpSolver() {
   1.989 +    messageLevel(MESSAGE_NO_OUTPUT);
   1.990 +  }
   1.991 +
   1.992 +  LpGlpk::LpGlpk(const LpGlpk& other)
   1.993 +    : LpBase(other), GlpkBase(other), LpSolver(other) {
   1.994 +    messageLevel(MESSAGE_NO_OUTPUT);
   1.995 +  }
   1.996 +
   1.997 +  LpGlpk* LpGlpk::_newSolver() const { return new LpGlpk; }
   1.998 +  LpGlpk* LpGlpk::_cloneSolver() const { return new LpGlpk(*this); }
   1.999 +
  1.1000 +  const char* LpGlpk::_solverName() const { return "LpGlpk"; }
  1.1001 +
  1.1002 +  void LpGlpk::_clear_temporals() {
  1.1003 +    _primal_ray.clear();
  1.1004 +    _dual_ray.clear();
  1.1005 +  }
  1.1006 +
  1.1007 +  LpGlpk::SolveExitStatus LpGlpk::_solve() {
  1.1008 +    return solvePrimal();
  1.1009 +  }
  1.1010 +
  1.1011 +  LpGlpk::SolveExitStatus LpGlpk::solvePrimal() {
  1.1012 +    _clear_temporals();
  1.1013 +
  1.1014 +    glp_smcp smcp;
  1.1015 +    glp_init_smcp(&smcp);
  1.1016 +
  1.1017 +    switch (_message_level) {
  1.1018 +    case MESSAGE_NO_OUTPUT:
  1.1019 +      smcp.msg_lev = GLP_MSG_OFF;
  1.1020 +      break;
  1.1021 +    case MESSAGE_ERROR_MESSAGE:
  1.1022 +      smcp.msg_lev = GLP_MSG_ERR;
  1.1023 +      break;
  1.1024 +    case MESSAGE_NORMAL_OUTPUT:
  1.1025 +      smcp.msg_lev = GLP_MSG_ON;
  1.1026 +      break;
  1.1027 +    case MESSAGE_FULL_OUTPUT:
  1.1028 +      smcp.msg_lev = GLP_MSG_ALL;
  1.1029 +      break;
  1.1030 +    }
  1.1031 +
  1.1032 +    if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
  1.1033 +    return SOLVED;
  1.1034 +  }
  1.1035 +
  1.1036 +  LpGlpk::SolveExitStatus LpGlpk::solveDual() {
  1.1037 +    _clear_temporals();
  1.1038 +
  1.1039 +    glp_smcp smcp;
  1.1040 +    glp_init_smcp(&smcp);
  1.1041 +
  1.1042 +    switch (_message_level) {
  1.1043 +    case MESSAGE_NO_OUTPUT:
  1.1044 +      smcp.msg_lev = GLP_MSG_OFF;
  1.1045 +      break;
  1.1046 +    case MESSAGE_ERROR_MESSAGE:
  1.1047 +      smcp.msg_lev = GLP_MSG_ERR;
  1.1048 +      break;
  1.1049 +    case MESSAGE_NORMAL_OUTPUT:
  1.1050 +      smcp.msg_lev = GLP_MSG_ON;
  1.1051 +      break;
  1.1052 +    case MESSAGE_FULL_OUTPUT:
  1.1053 +      smcp.msg_lev = GLP_MSG_ALL;
  1.1054 +      break;
  1.1055 +    }
  1.1056 +    smcp.meth = GLP_DUAL;
  1.1057 +
  1.1058 +    if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
  1.1059 +    return SOLVED;
  1.1060 +  }
  1.1061 +
  1.1062 +  LpGlpk::Value LpGlpk::_getPrimal(int i) const {
  1.1063 +    return glp_get_col_prim(lp, i);
  1.1064 +  }
  1.1065 +
  1.1066 +  LpGlpk::Value LpGlpk::_getDual(int i) const {
  1.1067 +    return glp_get_row_dual(lp, i);
  1.1068 +  }
  1.1069 +
  1.1070 +  LpGlpk::Value LpGlpk::_getPrimalValue() const {
  1.1071 +    return glp_get_obj_val(lp);
  1.1072 +  }
  1.1073 +
  1.1074 +  LpGlpk::VarStatus LpGlpk::_getColStatus(int i) const {
  1.1075 +    switch (glp_get_col_stat(lp, i)) {
  1.1076 +    case GLP_BS:
  1.1077 +      return BASIC;
  1.1078 +    case GLP_UP:
  1.1079 +      return UPPER;
  1.1080 +    case GLP_LO:
  1.1081 +      return LOWER;
  1.1082 +    case GLP_NF:
  1.1083 +      return FREE;
  1.1084 +    case GLP_NS:
  1.1085 +      return FIXED;
  1.1086 +    default:
  1.1087 +      LEMON_ASSERT(false, "Wrong column status");
  1.1088 +      return LpGlpk::VarStatus();
  1.1089 +    }
  1.1090 +  }
  1.1091 +
  1.1092 +  LpGlpk::VarStatus LpGlpk::_getRowStatus(int i) const {
  1.1093 +    switch (glp_get_row_stat(lp, i)) {
  1.1094 +    case GLP_BS:
  1.1095 +      return BASIC;
  1.1096 +    case GLP_UP:
  1.1097 +      return UPPER;
  1.1098 +    case GLP_LO:
  1.1099 +      return LOWER;
  1.1100 +    case GLP_NF:
  1.1101 +      return FREE;
  1.1102 +    case GLP_NS:
  1.1103 +      return FIXED;
  1.1104 +    default:
  1.1105 +      LEMON_ASSERT(false, "Wrong row status");
  1.1106 +      return LpGlpk::VarStatus();
  1.1107 +    }
  1.1108 +  }
  1.1109 +
  1.1110 +  LpGlpk::Value LpGlpk::_getPrimalRay(int i) const {
  1.1111 +    if (_primal_ray.empty()) {
  1.1112 +      int row_num = glp_get_num_rows(lp);
  1.1113 +      int col_num = glp_get_num_cols(lp);
  1.1114 +
  1.1115 +      _primal_ray.resize(col_num + 1, 0.0);
  1.1116 +
  1.1117 +      int index = glp_get_unbnd_ray(lp);
  1.1118 +      if (index != 0) {
  1.1119 +        // The primal ray is found in primal simplex second phase
  1.1120 +        LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :
  1.1121 +                      glp_get_col_stat(lp, index - row_num)) != GLP_BS,
  1.1122 +                     "Wrong primal ray");
  1.1123 +
  1.1124 +        bool negate = glp_get_obj_dir(lp) == GLP_MAX;
  1.1125 +
  1.1126 +        if (index > row_num) {
  1.1127 +          _primal_ray[index - row_num] = 1.0;
  1.1128 +          if (glp_get_col_dual(lp, index - row_num) > 0) {
  1.1129 +            negate = !negate;
  1.1130 +          }
  1.1131 +        } else {
  1.1132 +          if (glp_get_row_dual(lp, index) > 0) {
  1.1133 +            negate = !negate;
  1.1134 +          }
  1.1135 +        }
  1.1136 +
  1.1137 +        std::vector<int> ray_indexes(row_num + 1);
  1.1138 +        std::vector<Value> ray_values(row_num + 1);
  1.1139 +        int ray_length = glp_eval_tab_col(lp, index, &ray_indexes.front(),
  1.1140 +                                          &ray_values.front());
  1.1141 +
  1.1142 +        for (int i = 1; i <= ray_length; ++i) {
  1.1143 +          if (ray_indexes[i] > row_num) {
  1.1144 +            _primal_ray[ray_indexes[i] - row_num] = ray_values[i];
  1.1145 +          }
  1.1146 +        }
  1.1147 +
  1.1148 +        if (negate) {
  1.1149 +          for (int i = 1; i <= col_num; ++i) {
  1.1150 +            _primal_ray[i] = - _primal_ray[i];
  1.1151 +          }
  1.1152 +        }
  1.1153 +      } else {
  1.1154 +        for (int i = 1; i <= col_num; ++i) {
  1.1155 +          _primal_ray[i] = glp_get_col_prim(lp, i);
  1.1156 +        }
  1.1157 +      }
  1.1158 +    }
  1.1159 +    return _primal_ray[i];
  1.1160 +  }
  1.1161 +
  1.1162 +  LpGlpk::Value LpGlpk::_getDualRay(int i) const {
  1.1163 +    if (_dual_ray.empty()) {
  1.1164 +      int row_num = glp_get_num_rows(lp);
  1.1165 +
  1.1166 +      _dual_ray.resize(row_num + 1, 0.0);
  1.1167 +
  1.1168 +      int index = glp_get_unbnd_ray(lp);
  1.1169 +      if (index != 0) {
  1.1170 +        // The dual ray is found in dual simplex second phase
  1.1171 +        LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :
  1.1172 +                      glp_get_col_stat(lp, index - row_num)) == GLP_BS,
  1.1173 +
  1.1174 +                     "Wrong dual ray");
  1.1175 +
  1.1176 +        int idx;
  1.1177 +        bool negate = false;
  1.1178 +
  1.1179 +        if (index > row_num) {
  1.1180 +          idx = glp_get_col_bind(lp, index - row_num);
  1.1181 +          if (glp_get_col_prim(lp, index - row_num) >
  1.1182 +              glp_get_col_ub(lp, index - row_num)) {
  1.1183 +            negate = true;
  1.1184 +          }
  1.1185 +        } else {
  1.1186 +          idx = glp_get_row_bind(lp, index);
  1.1187 +          if (glp_get_row_prim(lp, index) > glp_get_row_ub(lp, index)) {
  1.1188 +            negate = true;
  1.1189 +          }
  1.1190 +        }
  1.1191 +
  1.1192 +        _dual_ray[idx] = negate ?  - 1.0 : 1.0;
  1.1193 +
  1.1194 +        glp_btran(lp, &_dual_ray.front());
  1.1195 +      } else {
  1.1196 +        double eps = 1e-7;
  1.1197 +        // The dual ray is found in primal simplex first phase
  1.1198 +        // We assume that the glpk minimizes the slack to get feasible solution
  1.1199 +        for (int i = 1; i <= row_num; ++i) {
  1.1200 +          int index = glp_get_bhead(lp, i);
  1.1201 +          if (index <= row_num) {
  1.1202 +            double res = glp_get_row_prim(lp, index);
  1.1203 +            if (res > glp_get_row_ub(lp, index) + eps) {
  1.1204 +              _dual_ray[i] = -1;
  1.1205 +            } else if (res < glp_get_row_lb(lp, index) - eps) {
  1.1206 +              _dual_ray[i] = 1;
  1.1207 +            } else {
  1.1208 +              _dual_ray[i] = 0;
  1.1209 +            }
  1.1210 +            _dual_ray[i] *= glp_get_rii(lp, index);
  1.1211 +          } else {
  1.1212 +            double res = glp_get_col_prim(lp, index - row_num);
  1.1213 +            if (res > glp_get_col_ub(lp, index - row_num) + eps) {
  1.1214 +              _dual_ray[i] = -1;
  1.1215 +            } else if (res < glp_get_col_lb(lp, index - row_num) - eps) {
  1.1216 +              _dual_ray[i] = 1;
  1.1217 +            } else {
  1.1218 +              _dual_ray[i] = 0;
  1.1219 +            }
  1.1220 +            _dual_ray[i] /= glp_get_sjj(lp, index - row_num);
  1.1221 +          }
  1.1222 +        }
  1.1223 +
  1.1224 +        glp_btran(lp, &_dual_ray.front());
  1.1225 +
  1.1226 +        for (int i = 1; i <= row_num; ++i) {
  1.1227 +          _dual_ray[i] /= glp_get_rii(lp, i);
  1.1228 +        }
  1.1229 +      }
  1.1230 +    }
  1.1231 +    return _dual_ray[i];
  1.1232 +  }
  1.1233 +
  1.1234 +  LpGlpk::ProblemType LpGlpk::_getPrimalType() const {
  1.1235 +    if (glp_get_status(lp) == GLP_OPT)
  1.1236 +      return OPTIMAL;
  1.1237 +    switch (glp_get_prim_stat(lp)) {
  1.1238 +    case GLP_UNDEF:
  1.1239 +      return UNDEFINED;
  1.1240 +    case GLP_FEAS:
  1.1241 +    case GLP_INFEAS:
  1.1242 +      if (glp_get_dual_stat(lp) == GLP_NOFEAS) {
  1.1243 +        return UNBOUNDED;
  1.1244 +      } else {
  1.1245 +        return UNDEFINED;
  1.1246 +      }
  1.1247 +    case GLP_NOFEAS:
  1.1248 +      return INFEASIBLE;
  1.1249 +    default:
  1.1250 +      LEMON_ASSERT(false, "Wrong primal type");
  1.1251 +      return  LpGlpk::ProblemType();
  1.1252 +    }
  1.1253 +  }
  1.1254 +
  1.1255 +  LpGlpk::ProblemType LpGlpk::_getDualType() const {
  1.1256 +    if (glp_get_status(lp) == GLP_OPT)
  1.1257 +      return OPTIMAL;
  1.1258 +    switch (glp_get_dual_stat(lp)) {
  1.1259 +    case GLP_UNDEF:
  1.1260 +      return UNDEFINED;
  1.1261 +    case GLP_FEAS:
  1.1262 +    case GLP_INFEAS:
  1.1263 +      if (glp_get_prim_stat(lp) == GLP_NOFEAS) {
  1.1264 +        return UNBOUNDED;
  1.1265 +      } else {
  1.1266 +        return UNDEFINED;
  1.1267 +      }
  1.1268 +    case GLP_NOFEAS:
  1.1269 +      return INFEASIBLE;
  1.1270 +    default:
  1.1271 +      LEMON_ASSERT(false, "Wrong primal type");
  1.1272 +      return  LpGlpk::ProblemType();
  1.1273 +    }
  1.1274 +  }
  1.1275 +
  1.1276 +  void LpGlpk::presolver(bool b) {
  1.1277 +    lpx_set_int_parm(lp, LPX_K_PRESOL, b ? 1 : 0);
  1.1278 +  }
  1.1279 +
  1.1280 +  void LpGlpk::messageLevel(MessageLevel m) {
  1.1281 +    _message_level = m;
  1.1282 +  }
  1.1283 +
  1.1284 +  // MipGlpk members
  1.1285 +
  1.1286 +  MipGlpk::MipGlpk()
  1.1287 +    : LpBase(), GlpkBase(), MipSolver() {
  1.1288 +    messageLevel(MESSAGE_NO_OUTPUT);
  1.1289 +  }
  1.1290 +
  1.1291 +  MipGlpk::MipGlpk(const MipGlpk& other)
  1.1292 +    : LpBase(), GlpkBase(other), MipSolver() {
  1.1293 +    messageLevel(MESSAGE_NO_OUTPUT);
  1.1294 +  }
  1.1295 +
  1.1296 +  void MipGlpk::_setColType(int i, MipGlpk::ColTypes col_type) {
  1.1297 +    switch (col_type) {
  1.1298 +    case INTEGER:
  1.1299 +      glp_set_col_kind(lp, i, GLP_IV);
  1.1300 +      break;
  1.1301 +    case REAL:
  1.1302 +      glp_set_col_kind(lp, i, GLP_CV);
  1.1303 +      break;
  1.1304 +    }
  1.1305 +  }
  1.1306 +
  1.1307 +  MipGlpk::ColTypes MipGlpk::_getColType(int i) const {
  1.1308 +    switch (glp_get_col_kind(lp, i)) {
  1.1309 +    case GLP_IV:
  1.1310 +    case GLP_BV:
  1.1311 +      return INTEGER;
  1.1312 +    default:
  1.1313 +      return REAL;
  1.1314 +    }
  1.1315 +
  1.1316 +  }
  1.1317 +
  1.1318 +  MipGlpk::SolveExitStatus MipGlpk::_solve() {
  1.1319 +    glp_smcp smcp;
  1.1320 +    glp_init_smcp(&smcp);
  1.1321 +
  1.1322 +    switch (_message_level) {
  1.1323 +    case MESSAGE_NO_OUTPUT:
  1.1324 +      smcp.msg_lev = GLP_MSG_OFF;
  1.1325 +      break;
  1.1326 +    case MESSAGE_ERROR_MESSAGE:
  1.1327 +      smcp.msg_lev = GLP_MSG_ERR;
  1.1328 +      break;
  1.1329 +    case MESSAGE_NORMAL_OUTPUT:
  1.1330 +      smcp.msg_lev = GLP_MSG_ON;
  1.1331 +      break;
  1.1332 +    case MESSAGE_FULL_OUTPUT:
  1.1333 +      smcp.msg_lev = GLP_MSG_ALL;
  1.1334 +      break;
  1.1335 +    }
  1.1336 +    smcp.meth = GLP_DUAL;
  1.1337 +
  1.1338 +    if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
  1.1339 +    if (glp_get_status(lp) != GLP_OPT) return SOLVED;
  1.1340 +
  1.1341 +    glp_iocp iocp;
  1.1342 +    glp_init_iocp(&iocp);
  1.1343 +
  1.1344 +    switch (_message_level) {
  1.1345 +    case MESSAGE_NO_OUTPUT:
  1.1346 +      iocp.msg_lev = GLP_MSG_OFF;
  1.1347 +      break;
  1.1348 +    case MESSAGE_ERROR_MESSAGE:
  1.1349 +      iocp.msg_lev = GLP_MSG_ERR;
  1.1350 +      break;
  1.1351 +    case MESSAGE_NORMAL_OUTPUT:
  1.1352 +      iocp.msg_lev = GLP_MSG_ON;
  1.1353 +      break;
  1.1354 +    case MESSAGE_FULL_OUTPUT:
  1.1355 +      iocp.msg_lev = GLP_MSG_ALL;
  1.1356 +      break;
  1.1357 +    }
  1.1358 +
  1.1359 +    if (glp_intopt(lp, &iocp) != 0) return UNSOLVED;
  1.1360 +    return SOLVED;
  1.1361 +  }
  1.1362 +
  1.1363 +
  1.1364 +  MipGlpk::ProblemType MipGlpk::_getType() const {
  1.1365 +    switch (glp_get_status(lp)) {
  1.1366 +    case GLP_OPT:
  1.1367 +      switch (glp_mip_status(lp)) {
  1.1368 +      case GLP_UNDEF:
  1.1369 +        return UNDEFINED;
  1.1370 +      case GLP_NOFEAS:
  1.1371 +        return INFEASIBLE;
  1.1372 +      case GLP_FEAS:
  1.1373 +        return FEASIBLE;
  1.1374 +      case GLP_OPT:
  1.1375 +        return OPTIMAL;
  1.1376 +      default:
  1.1377 +        LEMON_ASSERT(false, "Wrong problem type.");
  1.1378 +        return MipGlpk::ProblemType();
  1.1379 +      }
  1.1380 +    case GLP_NOFEAS:
  1.1381 +      return INFEASIBLE;
  1.1382 +    case GLP_INFEAS:
  1.1383 +    case GLP_FEAS:
  1.1384 +      if (glp_get_dual_stat(lp) == GLP_NOFEAS) {
  1.1385 +        return UNBOUNDED;
  1.1386 +      } else {
  1.1387 +        return UNDEFINED;
  1.1388 +      }
  1.1389 +    default:
  1.1390 +      LEMON_ASSERT(false, "Wrong problem type.");
  1.1391 +      return MipGlpk::ProblemType();
  1.1392 +    }
  1.1393 +  }
  1.1394 +
  1.1395 +  MipGlpk::Value MipGlpk::_getSol(int i) const {
  1.1396 +    return glp_mip_col_val(lp, i);
  1.1397 +  }
  1.1398 +
  1.1399 +  MipGlpk::Value MipGlpk::_getSolValue() const {
  1.1400 +    return glp_mip_obj_val(lp);
  1.1401 +  }
  1.1402 +
  1.1403 +  MipGlpk* MipGlpk::_newSolver() const { return new MipGlpk; }
  1.1404 +  MipGlpk* MipGlpk::_cloneSolver() const {return new MipGlpk(*this); }
  1.1405 +
  1.1406 +  const char* MipGlpk::_solverName() const { return "MipGlpk"; }
  1.1407 +
  1.1408 +  void MipGlpk::messageLevel(MessageLevel m) {
  1.1409 +    _message_level = m;
  1.1410 +  }
  1.1411  
  1.1412  } //END OF NAMESPACE LEMON