RhodeCode

        break;

      default:

        break;

      }

    } else {

      switch (b) {

      case GLP_FR:

      case GLP_LO:

        glp_set_col_bnds(lp, i, GLP_LO, lo, up);

        break;

      case GLP_UP:

      case GLP_DB:

      case GLP_FX:

        if (lo == up)

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

        else

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

        break;

      default:

        break;

      }

    }

  }

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

    int b = glp_get_col_type(lp, i);

    switch (b) {

    case GLP_LO:

    case GLP_DB:

    case GLP_FX:

      return glp_get_col_lb(lp, i);

    default:

      return -INF;

    }

  }

  void GlpkBase::_setColUpperBound(int i, Value up) {

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

    int b = glp_get_col_type(lp, i);

    double lo = glp_get_col_lb(lp, i);

    if (up == INF) {

      switch (b) {

      case GLP_FR:

      case GLP_LO:

        break;

      case GLP_UP:

        glp_set_col_bnds(lp, i, GLP_FR, lo, up);

        break;

      case GLP_DB:

      case GLP_FX:

        glp_set_col_bnds(lp, i, GLP_LO, lo, up);

        break;

      default:

        break;

      }

    } else {

      switch (b) {

      case GLP_FR:

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

        break;

      case GLP_UP:

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

        break;

      case GLP_LO:

      case GLP_DB:

      case GLP_FX:

        if (lo == up)

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

        else

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

        break;

      default:

        break;

      }

    }

  }

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

    int b = glp_get_col_type(lp, i);

      switch (b) {

      case GLP_UP:

      case GLP_DB:

      case GLP_FX:

        return glp_get_col_ub(lp, i);

      default:

        return INF;

      }

  }

  void GlpkBase::_setRowLowerBound(int i, Value lo) {

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

    int b = glp_get_row_type(lp, i);

    double up = glp_get_row_ub(lp, i);

    if (lo == -INF) {

      switch (b) {

      case GLP_FR:

      case GLP_LO:

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

        break;

      case GLP_UP:

        break;

      case GLP_DB:

      case GLP_FX:

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

        break;

      default:

        break;

      }

    } else {

      switch (b) {

      case GLP_FR:

      case GLP_LO:

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

        break;

      case GLP_UP:

      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::_getRowLowerBound(int i) const {

    int b = glp_get_row_type(lp, i);

    switch (b) {

    case GLP_LO:

    case GLP_DB:

    case GLP_FX:

      return glp_get_row_lb(lp, i);

    default:

      return -INF;

    }

  }

  void GlpkBase::_setRowUpperBound(int i, Value up) {

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

    int b = glp_get_row_type(lp, i);

    double lo = glp_get_row_lb(lp, i);

    if (up == INF) {

      switch (b) {

      case GLP_FR:

      case GLP_LO:

        break;

      case GLP_UP:

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

        break;

      case GLP_DB:

      case GLP_FX:

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

        break;

      default:

        break;

      }

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

  }

  // 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:

      return LOWER;

    case GLP_NF:

      return FREE;

    case GLP_NS:

      return FIXED;

    default:

      LEMON_ASSERT(false, "Wrong column status");

      return GlpkLp::VarStatus();

    }

  }

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

    switch (glp_get_row_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 row status");

      return GlpkLp::VarStatus();

    }

  }

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

    if (_primal_ray.empty()) {

      int row_num = glp_get_num_rows(lp);

      int col_num = glp_get_num_cols(lp);

      _primal_ray.resize(col_num + 1, 0.0);

      int index = glp_get_unbnd_ray(lp);

      if (index != 0) {

        // The primal ray is found in primal simplex second phase

        LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :

                      glp_get_col_stat(lp, index - row_num)) != GLP_BS,

                     "Wrong primal ray");

        bool negate = glp_get_obj_dir(lp) == GLP_MAX;

        if (index > row_num) {

          _primal_ray[index - row_num] = 1.0;

          if (glp_get_col_dual(lp, index - row_num) > 0) {

            negate = !negate;

          }

        } else {

          if (glp_get_row_dual(lp, index) > 0) {

            negate = !negate;

          }

        }

        std::vector<int> ray_indexes(row_num + 1);

        std::vector<Value> ray_values(row_num + 1);

        int ray_length = glp_eval_tab_col(lp, index, &ray_indexes.front(),

                                          &ray_values.front());

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

          if (ray_indexes[i] > row_num) {

            _primal_ray[ray_indexes[i] - row_num] = ray_values[i];

          }

        }

        if (negate) {

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

            _primal_ray[i] = - _primal_ray[i];

          }

        }

      } else {

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

          _primal_ray[i] = glp_get_col_prim(lp, i);

        }

      }

    }

    return _primal_ray[i];

  }

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

    if (_dual_ray.empty()) {

      int row_num = glp_get_num_rows(lp);

      _dual_ray.resize(row_num + 1, 0.0);

      int index = glp_get_unbnd_ray(lp);

      if (index != 0) {

        // The dual ray is found in dual simplex second phase

        LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :

                      glp_get_col_stat(lp, index - row_num)) == GLP_BS,

                     "Wrong dual ray");

        int idx;

        bool negate = false;

        if (index > row_num) {

          idx = glp_get_col_bind(lp, index - row_num);

          if (glp_get_col_prim(lp, index - row_num) >

              glp_get_col_ub(lp, index - row_num)) {

            negate = true;

          }

        } else {

          idx = glp_get_row_bind(lp, index);

          if (glp_get_row_prim(lp, index) > glp_get_row_ub(lp, index)) {

            negate = true;

          }

        }

        _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) {

          int index = glp_get_bhead(lp, i);

          if (index <= row_num) {

            double res = glp_get_row_prim(lp, index);

            if (res > glp_get_row_ub(lp, index) + eps) {

              _dual_ray[i] = -1;

            } else if (res < glp_get_row_lb(lp, index) - eps) {

              _dual_ray[i] = 1;

            } else {

              _dual_ray[i] = 0;

            }

            _dual_ray[i] *= glp_get_rii(lp, index);

          } else {

            double res = glp_get_col_prim(lp, index - row_num);

            if (res > glp_get_col_ub(lp, index - row_num) + eps) {

              _dual_ray[i] = -1;

            } else if (res < glp_get_col_lb(lp, index - row_num) - eps) {

              _dual_ray[i] = 1;

            } else {

              _dual_ray[i] = 0;

            }

            _dual_ray[i] /= glp_get_sjj(lp, index - row_num);

          }

        }

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

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

          _dual_ray[i] /= glp_get_rii(lp, i);

        }

      }

    }

    return _dual_ray[i];

  }

  GlpkLp::ProblemType GlpkLp::_getPrimalType() const {

    if (glp_get_status(lp) == GLP_OPT)

      return OPTIMAL;

    switch (glp_get_prim_stat(lp)) {

    case GLP_UNDEF:

      return UNDEFINED;

    case GLP_FEAS:

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

 *

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

 *

 * Copyright (C) 2003-2008

 * 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_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 LpSolver, public GlpkBase {

  public:

    ///\e

    GlpkLp();

    ///\e

    GlpkLp(const GlpkLp&);

    ///\e

    virtual GlpkLp* cloneSolver() const;

    ///\e

    virtual GlpkLp* newSolver() const;

  private:

    mutable std::vector<double> _primal_ray;

    mutable std::vector<double> _dual_ray;

    void _clear_temporals();

  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

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

  };

  /// \brief Interface for the GLPK MIP solver

  ///

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

  ///\ingroup lp_group

  class GlpkMip : public MipSolver, public GlpkBase {

  public:

    ///\e

    GlpkMip();

    ///\e

    GlpkMip(const GlpkMip&);

    virtual GlpkMip* cloneSolver() const;

    virtual GlpkMip* newSolver() const;

  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;

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

  };

} //END OF NAMESPACE LEMON

#endif //LEMON_GLPK_H