RhodeCode

/* -*- mode: C++; indent-tabs-mode: nil; -*-

 *

 * This file is a part of LEMON, a generic C++ optimization library.

 *

 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

 * (Egervary Research Group on Combinatorial Optimization, EGRES).

 *

 * Permission to use, modify and distribute this software is granted

 * provided that this copyright notice appears in all copies. For

 * precise terms see the accompanying LICENSE file.

 *

 * This software is provided "AS IS" with no warranty of any kind,

 * express or implied, and with no claim as to its suitability for any

 * purpose.

 *

 */

#include <iostream>

#include <vector>

#include <cstring>

#include <lemon/cplex.h>

extern "C" {

#include <ilcplex/cplex.h>

}

///\file

///\brief Implementation of the LEMON-CPLEX lp solver interface.

namespace lemon {

  CplexEnv::LicenseError::LicenseError(int status) {

    if (!CPXgeterrorstring(0, status, _message)) {

      std::strcpy(_message, "Cplex unknown error");

    }

  }

  CplexEnv::CplexEnv() {

    int status;

    _cnt = new int;

    _env = CPXopenCPLEX(&status);

    if (_env == 0) {

      delete _cnt;

      _cnt = 0;

      throw LicenseError(status);

    }

  }

  CplexEnv::CplexEnv(const CplexEnv& other) {

    _env = other._env;

    _cnt = other._cnt;

    ++(*_cnt);

  }

  CplexEnv& CplexEnv::operator=(const CplexEnv& other) {

    _env = other._env;

    _cnt = other._cnt;

    ++(*_cnt);

    return *this;

  }

  CplexEnv::~CplexEnv() {

    --(*_cnt);

    if (*_cnt == 0) {

      delete _cnt;

      CPXcloseCPLEX(&_env);

    }

  }

  CplexBase::CplexBase() : LpBase() {

    int status;

    _prob = CPXcreateprob(cplexEnv(), &status, "Cplex problem");

  }

  CplexBase::CplexBase(const CplexEnv& env)

    : LpBase(), _env(env) {

    int status;

    _prob = CPXcreateprob(cplexEnv(), &status, "Cplex problem");

  }

  CplexBase::CplexBase(const CplexBase& cplex)

    : LpBase() {

    int status;

    _prob = CPXcloneprob(cplexEnv(), cplex._prob, &status);

    rows = cplex.rows;

    cols = cplex.cols;

  }

  CplexBase::~CplexBase() {

    CPXfreeprob(cplexEnv(),&_prob);

  }

  int CplexBase::_addCol() {

    int i = CPXgetnumcols(cplexEnv(), _prob);

    double lb = -INF, ub = INF;

    CPXnewcols(cplexEnv(), _prob, 1, 0, &lb, &ub, 0, 0);

    return i;

  }

      CPXchgrhs(cplexEnv(), _prob, 1, &i, &lb);

    } else {

      const char s = 'R';

      CPXchgsense(cplexEnv(), _prob, 1, &i, &s);

      CPXchgrhs(cplexEnv(), _prob, 1, &i, &lb);

      double len = ub - lb;

      CPXchgrngval(cplexEnv(), _prob, 1, &i, &len);

    }

  }

  void CplexBase::_setRowLowerBound(int i, Value lb)

  {

    LEMON_ASSERT(lb != INF, "Invalid bound");

    _set_row_bounds(i, lb, CplexBase::_getRowUpperBound(i));

  }

  void CplexBase::_setRowUpperBound(int i, Value ub)

  {

    LEMON_ASSERT(ub != -INF, "Invalid bound");

    _set_row_bounds(i, CplexBase::_getRowLowerBound(i), ub);

  }

  void CplexBase::_setObjCoeffs(ExprIterator b, ExprIterator e)

  {

    std::vector<int> indices;

    std::vector<Value> values;

    for(ExprIterator it=b; it!=e; ++it) {

      indices.push_back(it->first);

      values.push_back(it->second);

    }

    CPXchgobj(cplexEnv(), _prob, values.size(),

              &indices.front(), &values.front());

  }

  void CplexBase::_getObjCoeffs(InsertIterator b) const

  {

    int num = CPXgetnumcols(cplexEnv(), _prob);

    std::vector<Value> x(num);

    CPXgetobj(cplexEnv(), _prob, &x.front(), 0, num - 1);

    for (int i = 0; i < num; ++i) {

      if (x[i] != 0.0) {

        *b = std::make_pair(i, x[i]);

        ++b;

      }

    }

  }

  void CplexBase::_setObjCoeff(int i, Value obj_coef)

  {

    CPXchgobj(cplexEnv(), _prob, 1, &i, &obj_coef);

  }

  CplexBase::Value CplexBase::_getObjCoeff(int i) const

  {

    Value x;

    CPXgetobj(cplexEnv(), _prob, &x, i, i);

    return x;

  }

  void CplexBase::_setSense(CplexBase::Sense sense) {

    switch (sense) {

    case MIN:

      CPXchgobjsen(cplexEnv(), _prob, CPX_MIN);

      break;

    case MAX:

      CPXchgobjsen(cplexEnv(), _prob, CPX_MAX);

      break;

    }

  }

  CplexBase::Sense CplexBase::_getSense() const {

    switch (CPXgetobjsen(cplexEnv(), _prob)) {

    case CPX_MIN:

      return MIN;

    case CPX_MAX:

      return MAX;

    default:

      LEMON_ASSERT(false, "Invalid sense");

      return CplexBase::Sense();

    }

  }

  void CplexBase::_clear() {

    CPXfreeprob(cplexEnv(),&_prob);

    int status;

    _prob = CPXcreateprob(cplexEnv(), &status, "Cplex problem");

    rows.clear();

    cols.clear();

  }

  // CplexLp members

  CplexLp::CplexLp()

  CplexLp::CplexLp(const CplexEnv& env)

  CplexLp::CplexLp(const CplexLp& other)

  CplexLp::~CplexLp() {}

  CplexLp* CplexLp::newSolver() const { return new CplexLp; }

  CplexLp* CplexLp::cloneSolver() const {return new CplexLp(*this); }

  const char* CplexLp::_solverName() const { return "CplexLp"; }

  void CplexLp::_clear_temporals() {

    _col_status.clear();

    _row_status.clear();

    _primal_ray.clear();

    _dual_ray.clear();

  }

  // The routine returns zero unless an error occurred during the

  // optimization. Examples of errors include exhausting available

  // memory (CPXERR_NO_MEMORY) or encountering invalid data in the

  // CPLEX problem object (CPXERR_NO_PROBLEM). Exceeding a

  // user-specified CPLEX limit, or proving the model infeasible or

  // unbounded, are not considered errors. Note that a zero return

  // value does not necessarily mean that a solution exists. Use query

  // routines CPXsolninfo, CPXgetstat, and CPXsolution to obtain

  // further information about the status of the optimization.

  CplexLp::SolveExitStatus CplexLp::convertStatus(int status) {

#if CPX_VERSION >= 800

    if (status == 0) {

      switch (CPXgetstat(cplexEnv(), _prob)) {

      case CPX_STAT_OPTIMAL:

      case CPX_STAT_INFEASIBLE:

      case CPX_STAT_UNBOUNDED:

        return SOLVED;

      default:

        return UNSOLVED;

      }

    } else {

      return UNSOLVED;

    }

#else

    if (status == 0) {

      //We want to exclude some cases

      switch (CPXgetstat(cplexEnv(), _prob)) {

      case CPX_OBJ_LIM:

      case CPX_IT_LIM_FEAS:

      case CPX_IT_LIM_INFEAS:

      case CPX_TIME_LIM_FEAS:

      case CPX_TIME_LIM_INFEAS:

        return UNSOLVED;

      default:

        return SOLVED;

      }

    } else {

      return UNSOLVED;

    }

#endif

  }

  CplexLp::SolveExitStatus CplexLp::_solve() {

    _clear_temporals();

    return convertStatus(CPXlpopt(cplexEnv(), _prob));

  }

  CplexLp::SolveExitStatus CplexLp::solvePrimal() {

    _clear_temporals();

    return convertStatus(CPXprimopt(cplexEnv(), _prob));

  }

  CplexLp::SolveExitStatus CplexLp::solveDual() {

    _clear_temporals();

    return convertStatus(CPXdualopt(cplexEnv(), _prob));

  }

  CplexLp::SolveExitStatus CplexLp::solveBarrier() {

    _clear_temporals();

    return convertStatus(CPXbaropt(cplexEnv(), _prob));

  }

  CplexLp::Value CplexLp::_getPrimal(int i) const {

    Value x;

    CPXgetx(cplexEnv(), _prob, &x, i, i);

    return x;

  }

  CplexLp::Value CplexLp::_getDual(int i) const {

    Value y;

    CPXgetpi(cplexEnv(), _prob, &y, i, i);

    return y;

  }

  CplexLp::Value CplexLp::_getPrimalValue() const {

    Value objval;

    CPXgetobjval(cplexEnv(), _prob, &objval);

    return objval;

  }

  CplexLp::VarStatus CplexLp::_getColStatus(int i) const {

    int stat = CPXgetstat(cplexEnv(), _prob);

#if CPX_VERSION >= 800

    switch (stat)

      {

      case CPX_STAT_OPTIMAL:

        return OPTIMAL;

      case CPX_STAT_UNBOUNDED:

        return UNBOUNDED;

      case CPX_STAT_INFEASIBLE:

        return INFEASIBLE;

      default:

        return UNDEFINED;

      }

#else

    statusSwitch(cplexEnv(),stat);

    //CPXgetstat(cplexEnv(), _prob);

    //printf("A primal status: %d, CPX_OPTIMAL=%d \n",stat,CPX_OPTIMAL);

    switch (stat) {

    case 0:

      return UNDEFINED; //Undefined

    case CPX_OPTIMAL://Optimal

      return OPTIMAL;

    case CPX_UNBOUNDED://Unbounded

      return INFEASIBLE;//In case of dual simplex

      //return UNBOUNDED;

    case CPX_INFEASIBLE://Infeasible

      //    case CPX_IT_LIM_INFEAS:

      //     case CPX_TIME_LIM_INFEAS:

      //     case CPX_NUM_BEST_INFEAS:

      //     case CPX_OPTIMAL_INFEAS:

      //     case CPX_ABORT_INFEAS:

      //     case CPX_ABORT_PRIM_INFEAS:

      //     case CPX_ABORT_PRIM_DUAL_INFEAS:

      return UNBOUNDED;//In case of dual simplex

      //return INFEASIBLE;

      //     case CPX_OBJ_LIM:

      //     case CPX_IT_LIM_FEAS:

      //     case CPX_TIME_LIM_FEAS:

      //     case CPX_NUM_BEST_FEAS:

      //     case CPX_ABORT_FEAS:

      //     case CPX_ABORT_PRIM_DUAL_FEAS:

      //       return FEASIBLE;

    default:

      return UNDEFINED; //Everything else comes here

      //FIXME error

    }

#endif

  }

  //9.0-as cplex verzio statusai

  // CPX_STAT_ABORT_DUAL_OBJ_LIM

  // CPX_STAT_ABORT_IT_LIM

  // CPX_STAT_ABORT_OBJ_LIM

  // CPX_STAT_ABORT_PRIM_OBJ_LIM

  // CPX_STAT_ABORT_TIME_LIM

  // CPX_STAT_ABORT_USER

  // CPX_STAT_FEASIBLE_RELAXED

  // CPX_STAT_INFEASIBLE

  // CPX_STAT_INForUNBD

  // CPX_STAT_NUM_BEST

  // CPX_STAT_OPTIMAL

  // CPX_STAT_OPTIMAL_FACE_UNBOUNDED

  // CPX_STAT_OPTIMAL_INFEAS

  // CPX_STAT_OPTIMAL_RELAXED

  // CPX_STAT_UNBOUNDED

  CplexLp::ProblemType CplexLp::_getDualType() const {

    int stat = CPXgetstat(cplexEnv(), _prob);

#if CPX_VERSION >= 800

    switch (stat) {

    case CPX_STAT_OPTIMAL:

      return OPTIMAL;

    case CPX_STAT_UNBOUNDED:

      return INFEASIBLE;

    default:

      return UNDEFINED;

    }

#else

    statusSwitch(cplexEnv(),stat);

    switch (stat) {

    case 0:

      return UNDEFINED; //Undefined

    case CPX_OPTIMAL://Optimal

      return OPTIMAL;

    case CPX_UNBOUNDED:

      return INFEASIBLE;

    default:

      return UNDEFINED; //Everything else comes here

      //FIXME error

    }

#endif

  }

  // CplexMip members

  CplexMip::CplexMip()

#if CPX_VERSION < 800

    CPXchgprobtype(cplexEnv(),  _prob, CPXPROB_MIP);

#else

    CPXchgprobtype(cplexEnv(),  _prob, CPXPROB_MILP);

#endif

  }

  CplexMip::CplexMip(const CplexEnv& env)

#if CPX_VERSION < 800

    CPXchgprobtype(cplexEnv(),  _prob, CPXPROB_MIP);

#else

    CPXchgprobtype(cplexEnv(),  _prob, CPXPROB_MILP);

#endif

  }

  CplexMip::CplexMip(const CplexMip& other)

  CplexMip::~CplexMip() {}

  CplexMip* CplexMip::newSolver() const { return new CplexMip; }

  CplexMip* CplexMip::cloneSolver() const {return new CplexMip(*this); }

  const char* CplexMip::_solverName() const { return "CplexMip"; }

  void CplexMip::_setColType(int i, CplexMip::ColTypes col_type) {

    // Note If a variable is to be changed to binary, a call to CPXchgbds

    // should also be made to change the bounds to 0 and 1.

    switch (col_type){

    case INTEGER: {

      const char t = 'I';

      CPXchgctype (cplexEnv(), _prob, 1, &i, &t);

    } break;

    case REAL: {

      const char t = 'C';

      CPXchgctype (cplexEnv(), _prob, 1, &i, &t);

    } break;

    default:

      break;

    }

  }

  CplexMip::ColTypes CplexMip::_getColType(int i) const {

    char t;

    CPXgetctype (cplexEnv(), _prob, &t, i, i);

    switch (t) {

    case 'I':

      return INTEGER;

    case 'C':

      return REAL;

    default:

      LEMON_ASSERT(false, "Invalid column type");

      return ColTypes();

    }

  }

  CplexMip::SolveExitStatus CplexMip::_solve() {

    int status;

    status = CPXmipopt (cplexEnv(), _prob);

    if (status==0)

      return SOLVED;

    else

      return UNSOLVED;

  }

  CplexMip::ProblemType CplexMip::_getType() const {

    int stat = CPXgetstat(cplexEnv(), _prob);

    //Fortunately, MIP statuses did not change for cplex 8.0

    switch (stat) {

    case CPXMIP_OPTIMAL:

      // Optimal integer solution has been found.

    case CPXMIP_OPTIMAL_TOL:

      // Optimal soluton with the tolerance defined by epgap or epagap has

      // been found.

      return OPTIMAL;

      //This also exists in later issues

      //    case CPXMIP_UNBOUNDED:

      //return UNBOUNDED;

      case CPXMIP_INFEASIBLE:

        return INFEASIBLE;

    default:

      return UNDEFINED;

    }

    //Unboundedness not treated well: the following is from cplex 9.0 doc

    // About Unboundedness

    // The treatment of models that are unbounded involves a few

    // subtleties. Specifically, a declaration of unboundedness means that

    // ILOG CPLEX has determined that the model has an unbounded

    // ray. Given any feasible solution x with objective z, a multiple of

    // the unbounded ray can be added to x to give a feasible solution

    // with objective z-1 (or z+1 for maximization models). Thus, if a

    // feasible solution exists, then the optimal objective is

    // unbounded. Note that ILOG CPLEX has not necessarily concluded that

    // a feasible solution exists. Users can call the routine CPXsolninfo

    // to determine whether ILOG CPLEX has also concluded that the model

    // has a feasible solution.

  }

  CplexMip::Value CplexMip::_getSol(int i) const {

    Value x;

    CPXgetmipx(cplexEnv(), _prob, &x, i, i);

    return x;

  }

  CplexMip::Value CplexMip::_getSolValue() const {

/* -*- mode: C++; indent-tabs-mode: nil; -*-

 *

 * This file is a part of LEMON, a generic C++ optimization library.

 *

 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

 * (Egervary Research Group on Combinatorial Optimization, EGRES).

 *

 * Permission to use, modify and distribute this software is granted

 * provided that this copyright notice appears in all copies. For

 * precise terms see the accompanying LICENSE file.

 *

 * This software is provided "AS IS" with no warranty of any kind,

 * express or implied, and with no claim as to its suitability for any

 * purpose.

 *

 */

#ifndef LEMON_CPLEX_H

#define LEMON_CPLEX_H

///\file

///\brief Header of the LEMON-CPLEX lp solver interface.

#include <lemon/lp_base.h>

struct cpxenv;

struct cpxlp;

namespace lemon {

  /// \brief Reference counted wrapper around cpxenv pointer

  ///

  /// The cplex uses environment object which is responsible for

  /// checking the proper license usage. This class provides a simple

  /// interface for share the environment object between different

  /// problems.

  class CplexEnv {

    friend class CplexBase;

  private:

    cpxenv* _env;

    mutable int* _cnt;

  public:

    /// \brief This exception is thrown when the license check is not

    /// sufficient

    class LicenseError : public Exception {

      friend class CplexEnv;

    private:

      LicenseError(int status);

      char _message[510];

    public:

      /// The short error message

      virtual const char* what() const throw() {

        return _message;

      }

    };

    /// Constructor

    CplexEnv();

    /// Shallow copy constructor

    CplexEnv(const CplexEnv&);

    /// Shallow assignement

    CplexEnv& operator=(const CplexEnv&);

    /// Destructor

    virtual ~CplexEnv();

  protected:

    cpxenv* cplexEnv() { return _env; }

    const cpxenv* cplexEnv() const { return _env; }

  };

  /// \brief Base interface for the CPLEX LP and MIP solver

  ///

  /// This class implements the common interface of the CPLEX LP and

  /// \ingroup lp_group

  class CplexBase : virtual public LpBase {

  protected:

    CplexEnv _env;

    cpxlp* _prob;

    CplexBase();

    CplexBase(const CplexEnv&);

    CplexBase(const CplexBase &);

    virtual ~CplexBase();

    virtual int _addCol();

    virtual int _addRow();

    virtual void _eraseCol(int i);

    virtual void _eraseRow(int i);

    virtual void _eraseColId(int i);

    virtual void _eraseRowId(int i);

    virtual void _getColName(int col, std::string& name) const;

    virtual void _setColName(int col, const std::string& name);

    virtual int _colByName(const std::string& name) const;

    virtual void _getRowName(int row, std::string& name) const;

    virtual void _setRowName(int row, const std::string& name);

    virtual int _rowByName(const std::string& name) const;

    virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e);

    virtual void _getRowCoeffs(int i, InsertIterator b) const;

    virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e);

    virtual void _getColCoeffs(int i, InsertIterator b) const;

    virtual void _setCoeff(int row, int col, Value value);

    virtual Value _getCoeff(int row, int col) const;

    virtual void _setColLowerBound(int i, Value value);

    virtual Value _getColLowerBound(int i) const;

    virtual void _setColUpperBound(int i, Value value);

    virtual Value _getColUpperBound(int i) const;

  private:

    void _set_row_bounds(int i, Value lb, Value ub);

  protected:

    virtual void _setRowLowerBound(int i, Value value);

    virtual Value _getRowLowerBound(int i) const;

    virtual void _setRowUpperBound(int i, Value value);

    virtual Value _getRowUpperBound(int i) const;

    virtual void _setObjCoeffs(ExprIterator b, ExprIterator e);

    virtual void _getObjCoeffs(InsertIterator b) const;

    virtual void _setObjCoeff(int i, Value obj_coef);

    virtual Value _getObjCoeff(int i) const;

    virtual void _setSense(Sense sense);

    virtual Sense _getSense() const;

    virtual void _clear();

  public:

    /// Returns the used \c CplexEnv instance

    const CplexEnv& env() const { return _env; }

    ///

    const cpxenv* cplexEnv() const { return _env.cplexEnv(); }

    cpxlp* cplexLp() { return _prob; }

    const cpxlp* cplexLp() const { return _prob; }

  };

  /// \brief Interface for the CPLEX LP solver

  ///

  /// This class implements an interface for the CPLEX LP solver.

  ///\ingroup lp_group

  class CplexLp : public LpSolver, public CplexBase {

  public:

    /// \e

    CplexLp();

    /// \e

    CplexLp(const CplexEnv&);

    /// \e

    CplexLp(const CplexLp&);

    /// \e

    virtual ~CplexLp();

    /// \e

    virtual CplexLp* cloneSolver() const;

    /// \e

    virtual CplexLp* newSolver() const;

  private:

    // these values cannot retrieved element by element

    mutable std::vector<int> _col_status;

    mutable std::vector<int> _row_status;

    mutable std::vector<Value> _primal_ray;

    mutable std::vector<Value> _dual_ray;

    void _clear_temporals();

    SolveExitStatus convertStatus(int status);

  protected:

    virtual const char* _solverName() const;

    virtual SolveExitStatus _solve();

    virtual Value _getPrimal(int i) const;

    virtual Value _getDual(int i) const;

    virtual Value _getPrimalValue() const;

    virtual VarStatus _getColStatus(int i) const;

    virtual VarStatus _getRowStatus(int i) const;

    virtual Value _getPrimalRay(int i) const;

    virtual Value _getDualRay(int i) const;

    virtual ProblemType _getPrimalType() const;

    virtual ProblemType _getDualType() const;

  public:

    /// Solve with primal simplex method

    SolveExitStatus solvePrimal();

    /// Solve with dual simplex method

    SolveExitStatus solveDual();

    /// Solve with barrier method

    SolveExitStatus solveBarrier();

  };

  /// \brief Interface for the CPLEX MIP solver

  ///

  /// This class implements an interface for the CPLEX MIP solver.

  ///\ingroup lp_group

  class CplexMip : public MipSolver, public CplexBase {

  public:

    /// \e

    CplexMip();

    /// \e

    CplexMip(const CplexEnv&);

    /// \e

    CplexMip(const CplexMip&);

    /// \e

    virtual ~CplexMip();

  protected:

    virtual const char* _solverName() const;

    virtual ColTypes _getColType(int col) const;

    virtual void _setColType(int col, ColTypes col_type);

    virtual SolveExitStatus _solve();

    virtual ProblemType _getType() const;

    virtual Value _getSol(int i) const;

    virtual Value _getSolValue() const;

  };

} //END OF NAMESPACE LEMON

#endif //LEMON_CPLEX_H

/* -*- mode: C++; indent-tabs-mode: nil; -*-

 *

 * This file is a part of LEMON, a generic C++ optimization library.

 *

 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

 * (Egervary Research Group on Combinatorial Optimization, EGRES).

 *

 * Permission to use, modify and distribute this software is granted

 * provided that this copyright notice appears in all copies. For

 * precise terms see the accompanying LICENSE file.

 *

 * This software is provided "AS IS" with no warranty of any kind,

 * express or implied, and with no claim as to its suitability for any

 * purpose.

 *

 */

///\file

///\brief Implementation of the LEMON GLPK LP and MIP solver interface.

#include <lemon/glpk.h>

#include <glpk.h>

#include <lemon/assert.h>

namespace lemon {

  // GlpkBase members

  GlpkBase::GlpkBase() : LpBase() {

    lp = glp_create_prob();

    glp_create_index(lp);

  }

  GlpkBase::GlpkBase(const GlpkBase &other) : LpBase() {

    lp = glp_create_prob();

    glp_copy_prob(lp, other.lp, GLP_ON);

    glp_create_index(lp);

    rows = other.rows;

    cols = other.cols;

  }

  GlpkBase::~GlpkBase() {

    glp_delete_prob(lp);

  }

  int GlpkBase::_addCol() {

    int i = glp_add_cols(lp, 1);

    glp_set_col_bnds(lp, i, GLP_FR, 0.0, 0.0);

    return i;

  }

  int GlpkBase::_addRow() {

    int i = glp_add_rows(lp, 1);

    glp_set_row_bnds(lp, i, GLP_FR, 0.0, 0.0);

    return i;

  }

  void GlpkBase::_eraseCol(int i) {

    int ca[2];

    ca[1] = i;

    glp_del_cols(lp, 1, ca);

  }

  void GlpkBase::_eraseRow(int i) {

    int ra[2];

    ra[1] = i;

    glp_del_rows(lp, 1, ra);

  }

  void GlpkBase::_eraseColId(int i) {

    cols.eraseIndex(i);

    cols.shiftIndices(i);

  }

  void GlpkBase::_eraseRowId(int i) {

    rows.eraseIndex(i);

    rows.shiftIndices(i);

  }

  void GlpkBase::_getColName(int c, std::string& name) const {

    const char *str = glp_get_col_name(lp, c);

    if (str) name = str;

    else name.clear();

  }

  void GlpkBase::_setColName(int c, const std::string & name) {

    glp_set_col_name(lp, c, const_cast<char*>(name.c_str()));

  }

  int GlpkBase::_colByName(const std::string& name) const {

    int k = glp_find_col(lp, const_cast<char*>(name.c_str()));

    return k > 0 ? k : -1;

  }

  void GlpkBase::_getRowName(int r, std::string& name) const {

    const char *str = glp_get_row_name(lp, r);

    if (str) name = str;

    else name.clear();

      case GLP_UP:

        glp_set_row_bnds(lp, i, GLP_UP, lo, up);

        break;

      case GLP_LO:

      case GLP_DB:

      case GLP_FX:

        if (lo == up)

          glp_set_row_bnds(lp, i, GLP_FX, lo, up);

        else

          glp_set_row_bnds(lp, i, GLP_DB, lo, up);

        break;

      default:

        break;

      }

    }

  }

  GlpkBase::Value GlpkBase::_getRowUpperBound(int i) const {

    int b = glp_get_row_type(lp, i);

    switch (b) {

    case GLP_UP:

    case GLP_DB:

    case GLP_FX:

      return glp_get_row_ub(lp, i);

    default:

      return INF;

    }

  }

  void GlpkBase::_setObjCoeffs(ExprIterator b, ExprIterator e) {

    for (int i = 1; i <= glp_get_num_cols(lp); ++i) {

      glp_set_obj_coef(lp, i, 0.0);

    }

    for (ExprIterator it = b; it != e; ++it) {

      glp_set_obj_coef(lp, it->first, it->second);

    }

  }

  void GlpkBase::_getObjCoeffs(InsertIterator b) const {

    for (int i = 1; i <= glp_get_num_cols(lp); ++i) {

      Value val = glp_get_obj_coef(lp, i);

      if (val != 0.0) {

        *b = std::make_pair(i, val);

        ++b;

      }

    }

  }

  void GlpkBase::_setObjCoeff(int i, Value obj_coef) {

    //i = 0 means the constant term (shift)

    glp_set_obj_coef(lp, i, obj_coef);

  }

  GlpkBase::Value GlpkBase::_getObjCoeff(int i) const {

    //i = 0 means the constant term (shift)

    return glp_get_obj_coef(lp, i);

  }

  void GlpkBase::_setSense(GlpkBase::Sense sense) {

    switch (sense) {

    case MIN:

      glp_set_obj_dir(lp, GLP_MIN);

      break;

    case MAX:

      glp_set_obj_dir(lp, GLP_MAX);

      break;

    }

  }

  GlpkBase::Sense GlpkBase::_getSense() const {

    switch(glp_get_obj_dir(lp)) {

    case GLP_MIN:

      return MIN;

    case GLP_MAX:

      return MAX;

    default:

      LEMON_ASSERT(false, "Wrong sense");

      return GlpkBase::Sense();

    }

  }

  void GlpkBase::_clear() {

    glp_erase_prob(lp);

    rows.clear();

    cols.clear();

  }

  void GlpkBase::freeEnv() {

    glp_free_env();

  }

  GlpkBase::FreeEnvHelper GlpkBase::freeEnvHelper;

  // GlpkLp members

  GlpkLp::GlpkLp()

    messageLevel(MESSAGE_NO_OUTPUT);

  }

  GlpkLp::GlpkLp(const GlpkLp& other)

    messageLevel(MESSAGE_NO_OUTPUT);

  }

  GlpkLp* GlpkLp::newSolver() const { return new GlpkLp; }

  GlpkLp* GlpkLp::cloneSolver() const { return new GlpkLp(*this); }

  const char* GlpkLp::_solverName() const { return "GlpkLp"; }

  void GlpkLp::_clear_temporals() {

    _primal_ray.clear();

    _dual_ray.clear();

  }

  GlpkLp::SolveExitStatus GlpkLp::_solve() {

    return solvePrimal();

  }

  GlpkLp::SolveExitStatus GlpkLp::solvePrimal() {

    _clear_temporals();

    glp_smcp smcp;

    glp_init_smcp(&smcp);

    switch (_message_level) {

    case MESSAGE_NO_OUTPUT:

      smcp.msg_lev = GLP_MSG_OFF;

      break;

    case MESSAGE_ERROR_MESSAGE:

      smcp.msg_lev = GLP_MSG_ERR;

      break;

    case MESSAGE_NORMAL_OUTPUT:

      smcp.msg_lev = GLP_MSG_ON;

      break;

    case MESSAGE_FULL_OUTPUT:

      smcp.msg_lev = GLP_MSG_ALL;

      break;

    }

    if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;

    return SOLVED;

  }

  GlpkLp::SolveExitStatus GlpkLp::solveDual() {

    _clear_temporals();

    glp_smcp smcp;

    glp_init_smcp(&smcp);

    switch (_message_level) {

    case MESSAGE_NO_OUTPUT:

      smcp.msg_lev = GLP_MSG_OFF;

      break;

    case MESSAGE_ERROR_MESSAGE:

      smcp.msg_lev = GLP_MSG_ERR;

      break;

    case MESSAGE_NORMAL_OUTPUT:

      smcp.msg_lev = GLP_MSG_ON;

      break;

    case MESSAGE_FULL_OUTPUT:

      smcp.msg_lev = GLP_MSG_ALL;

      break;

    }

    smcp.meth = GLP_DUAL;

    if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;

    return SOLVED;

  }

  GlpkLp::Value GlpkLp::_getPrimal(int i) const {

    return glp_get_col_prim(lp, i);

  }

  GlpkLp::Value GlpkLp::_getDual(int i) const {

    return glp_get_row_dual(lp, i);

  }

  GlpkLp::Value GlpkLp::_getPrimalValue() const {

    return glp_get_obj_val(lp);

  }

  GlpkLp::VarStatus GlpkLp::_getColStatus(int i) const {

    switch (glp_get_col_stat(lp, i)) {

    case GLP_BS:

      return BASIC;

    case GLP_UP:

      return UPPER;

    case GLP_LO:

      return LOWER;

    case GLP_NF:

      return FREE;

    case GLP_NS:

      return FIXED;

    default:

      LEMON_ASSERT(false, "Wrong column status");

      return GlpkLp::VarStatus();

    }

        _dual_ray[idx] = negate ?  - 1.0 : 1.0;

        glp_btran(lp, &_dual_ray.front());

      } else {

        double eps = 1e-7;

        // The dual ray is found in primal simplex first phase

        // We assume that the glpk minimizes the slack to get feasible solution

        for (int i = 1; i <= row_num; ++i) {