RhodeCode

    } else {

      switch (b) {

      case GLP_FR:

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

        break;

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

  }

  // GlpkLp members

  GlpkLp::GlpkLp()

    : LpBase(), GlpkBase(), LpSolver() {

    messageLevel(MESSAGE_NO_OUTPUT);

  }

  GlpkLp::GlpkLp(const GlpkLp& other)

    : LpBase(other), GlpkBase(other), LpSolver(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:

 * (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_GLPK_H

#define LEMON_GLPK_H

///\file

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

///\ingroup lp_group

#include <lemon/lp_base.h>

// forward declaration

#ifndef _GLP_PROB

#define _GLP_PROB

typedef struct { double _prob; } glp_prob;

/* LP/MIP problem object */

#endif

namespace lemon {

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

  ///

  /// This class implements the common interface of the GLPK LP and MIP solver.

  /// \ingroup lp_group

  class GlpkBase : virtual public LpBase {

  protected:

    typedef glp_prob LPX;

    glp_prob* lp;

    GlpkBase();

    GlpkBase(const GlpkBase&);

    virtual ~GlpkBase();

  protected:

    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;

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

    virtual Sense _getSense() const;

    virtual void _clear();

  public:

    ///Pointer to the underlying GLPK data structure.

    LPX *lpx() {return lp;}

    ///Const pointer to the underlying GLPK data structure.

    const LPX *lpx() const {return lp;}

    ///Returns the constraint identifier understood by GLPK.

    int lpxRow(Row r) const { return rows(id(r)); }

    ///Returns the variable identifier understood by GLPK.

    int lpxCol(Col c) const { return cols(id(c)); }

  };

  /// \brief Interface for the GLPK LP solver

  ///

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

  ///\ingroup lp_group

  class GlpkLp : public GlpkBase, public LpSolver {

  public:

    ///\e

    GlpkLp();

    ///\e

    GlpkLp(const GlpkLp&);

  private:

    mutable std::vector<double> _primal_ray;

    mutable std::vector<double> _dual_ray;

    void _clear_temporals();

  protected:

    virtual GlpkLp* _cloneSolver() const;

    virtual GlpkLp* _newSolver() const;

    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;

    ///\todo It should be clarified

    ///

    virtual ProblemType _getPrimalType() const;

    virtual ProblemType _getDualType() const;

  public:

    ///Solve with primal simplex

    SolveExitStatus solvePrimal();

    ///Solve with dual simplex

    SolveExitStatus solveDual();

    ///Turns on or off the presolver

    ///Turns on (\c b is \c true) or off (\c b is \c false) the presolver

    ///

    ///The presolver is off by default.

    void presolver(bool b);

    ///Enum for \c messageLevel() parameter

    enum MessageLevel {

      /// no output (default value)

      MESSAGE_NO_OUTPUT = 0,

      /// error messages only

      MESSAGE_ERROR_MESSAGE = 1,

      /// normal output

      MESSAGE_NORMAL_OUTPUT = 2,

      /// full output (includes informational messages)

      MESSAGE_FULL_OUTPUT = 3

    };

  private:

    MessageLevel _message_level;

  public:

    ///Set the verbosity of the messages

    ///Set the verbosity of the messages

    ///

    ///\param m is the level of the messages output by the solver routines.

    void messageLevel(MessageLevel m);

  };

  lp.addRow(x1-x2 >=-1);

  //Nonnegativity of the variables

  lp.colLowerBound(x1, 0);

  lp.colLowerBound(x2, 0);

  //Objective function

  lp.obj(x1+x2);

  lp.sense(lp.MAX);

  //Testing the problem retrieving routines

  check(lp.objCoeff(x1)==1,"First term should be 1 in the obj function!");

  check(lp.sense() == lp.MAX,"This is a maximization!");

  check(lp.coeff(upright,x1)==1,"The coefficient in question is 1!");

  check(lp.colLowerBound(x1)==0,

        "The lower bound for variable x1 should be 0.");

  check(lp.colUpperBound(x1)==LpSolver::INF,

        "The upper bound for variable x1 should be infty.");

  check(lp.rowLowerBound(upright) == -LpSolver::INF,

        "The lower bound for the first row should be -infty.");

  check(lp.rowUpperBound(upright)==1,

        "The upper bound for the first row should be 1.");

  LpSolver::Expr e = lp.row(upright);

  check(e[x1] == 1, "The first coefficient should 1.");

  check(e[x2] == 2, "The second coefficient should 1.");

  lp.row(upright, x1+x2 <=1);

  e = lp.row(upright);

  check(e[x1] == 1, "The first coefficient should 1.");

  check(e[x2] == 1, "The second coefficient should 1.");

  LpSolver::DualExpr de = lp.col(x1);

  check(  de[upright] == 1, "The first coefficient should 1.");

  LpSolver* clp = lp.cloneSolver();

  //Testing the problem retrieving routines

  check(clp->objCoeff(x1)==1,"First term should be 1 in the obj function!");

  check(clp->sense() == clp->MAX,"This is a maximization!");

  check(clp->coeff(upright,x1)==1,"The coefficient in question is 1!");

  //  std::cout<<lp.colLowerBound(x1)<<std::endl;

  check(clp->colLowerBound(x1)==0,

        "The lower bound for variable x1 should be 0.");

  check(clp->colUpperBound(x1)==LpSolver::INF,

        "The upper bound for variable x1 should be infty.");

  check(lp.rowLowerBound(upright)==-LpSolver::INF,

        "The lower bound for the first row should be -infty.");

  check(lp.rowUpperBound(upright)==1,

        "The upper bound for the first row should be 1.");

  e = clp->row(upright);

  check(e[x1] == 1, "The first coefficient should 1.");

  check(e[x2] == 1, "The second coefficient should 1.");

  de = clp->col(x1);

  check(de[upright] == 1, "The first coefficient should 1.");

  delete clp;

  //Maximization of x1+x2

  //over the triangle with vertices (0,0) (0,1) (1,0)

  double expected_opt=1;

  solveAndCheck(lp, LpSolver::OPTIMAL, expected_opt);

  //Minimization

  lp.sense(lp.MIN);

  expected_opt=0;

  solveAndCheck(lp, LpSolver::OPTIMAL, expected_opt);

  //Vertex (-1,0) instead of (0,0)

  lp.colLowerBound(x1, -LpSolver::INF);

  expected_opt=-1;

  solveAndCheck(lp, LpSolver::OPTIMAL, expected_opt);

  //Erase one constraint and return to maximization

  lp.erase(upright);

  lp.sense(lp.MAX);

  expected_opt=LpSolver::INF;

  solveAndCheck(lp, LpSolver::UNBOUNDED, expected_opt);

  //Infeasibilty

  lp.addRow(x1+x2 <=-2);

  solveAndCheck(lp, LpSolver::INFEASIBLE, expected_opt);

}

int main()

{

  LpSkeleton lp_skel;

  lpTest(lp_skel);

#ifdef HAVE_GLPK

  {

    GlpkLp lp_glpk1,lp_glpk2;

    lpTest(lp_glpk1);

    aTest(lp_glpk2);

  }

#endif

#ifdef HAVE_CPLEX

  try {

    CplexLp lp_cplex1,lp_cplex2;

    lpTest(lp_cplex1);

    aTest(lp_cplex2);

  } catch (CplexEnv::LicenseError& error) {

#ifdef LEMON_FORCE_CPLEX_CHECK

    check(false, error.what());

#else

    std::cerr << error.what() << std::endl;

    std::cerr << "Cplex license check failed, lp check skipped" << std::endl;

#endif

  }

#endif

#ifdef HAVE_SOPLEX

  {

    SoplexLp lp_soplex1,lp_soplex2;

    lpTest(lp_soplex1);

    aTest(lp_soplex2);

  }

#endif

#ifdef HAVE_CLP

  {

    ClpLp lp_clp1,lp_clp2;

    lpTest(lp_clp1);

    aTest(lp_clp2);

  }

#endif

  return 0;

}

#ifdef HAVE_CONFIG_H

#include <lemon/config.h>

#endif

#ifdef HAVE_CPLEX

#include <lemon/cplex.h>

#endif

#ifdef HAVE_GLPK

#include <lemon/glpk.h>

#endif

using namespace lemon;

void solveAndCheck(MipSolver& mip, MipSolver::ProblemType stat,

                   double exp_opt) {

  using std::string;

  mip.solve();

  //int decimal,sign;

  std::ostringstream buf;

  buf << "Type should be: " << int(stat)<<" and it is "<<int(mip.type());

  //  itoa(stat,buf1, 10);

  check(mip.type()==stat, buf.str());

  if (stat ==  MipSolver::OPTIMAL) {

    std::ostringstream sbuf;

    buf << "Wrong optimal value: the right optimum is " << exp_opt;

    check(std::abs(mip.solValue()-exp_opt) < 1e-3, sbuf.str());

    //+ecvt(exp_opt,2)

  }

}

void aTest(MipSolver& mip)

{

 //The following example is very simple

  typedef MipSolver::Row Row;

  typedef MipSolver::Col Col;

  Col x1 = mip.addCol();

  Col x2 = mip.addCol();

  //Objective function

  mip.obj(x1);

  mip.max();

  //Unconstrained optimization

  mip.solve();

  //Check it out!

  //Constraints

  mip.addRow(2*x1+x2 <=2);

  mip.addRow(x1-2*x2 <=0);

  //Nonnegativity of the variable x1

  mip.colLowerBound(x1, 0);

  //Maximization of x1

  //over the triangle with vertices (0,0),(4/5,2/5),(0,2)

  double expected_opt=4.0/5.0;

  solveAndCheck(mip, MipSolver::OPTIMAL, expected_opt);

  //Restrict x2 to integer

  mip.colType(x2,MipSolver::INTEGER);

  expected_opt=1.0/2.0;

  solveAndCheck(mip, MipSolver::OPTIMAL, expected_opt);

  //Restrict both to integer

  mip.colType(x1,MipSolver::INTEGER);

  expected_opt=0;

  solveAndCheck(mip, MipSolver::OPTIMAL, expected_opt);

}

int main()

{

#ifdef HAVE_GLPK

  {

    GlpkMip mip1;

    aTest(mip1);

  }

#endif

#ifdef HAVE_CPLEX

  try {

    CplexMip mip2;

    aTest(mip2);

  } catch (CplexEnv::LicenseError& error) {

#ifdef LEMON_FORCE_CPLEX_CHECK

    check(false, error.what());

#else

    std::cerr << error.what() << std::endl;

    std::cerr << "Cplex license check failed, lp check skipped" << std::endl;

#endif

  }

#endif

  return 0;

}