alpar@461: /* -*- mode: C++; indent-tabs-mode: nil; -*-
alpar@461:  *
alpar@461:  * This file is a part of LEMON, a generic C++ optimization library.
alpar@461:  *
alpar@877:  * Copyright (C) 2003-2010
alpar@461:  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
alpar@461:  * (Egervary Research Group on Combinatorial Optimization, EGRES).
alpar@461:  *
alpar@461:  * Permission to use, modify and distribute this software is granted
alpar@461:  * provided that this copyright notice appears in all copies. For
alpar@461:  * precise terms see the accompanying LICENSE file.
alpar@461:  *
alpar@461:  * This software is provided "AS IS" with no warranty of any kind,
alpar@461:  * express or implied, and with no claim as to its suitability for any
alpar@461:  * purpose.
alpar@461:  *
alpar@461:  */
alpar@461: 
alpar@461: ///\file
alpar@461: ///\brief Implementation of the LEMON GLPK LP and MIP solver interface.
alpar@461: 
alpar@461: #include <lemon/glpk.h>
alpar@461: #include <glpk.h>
alpar@461: 
alpar@461: #include <lemon/assert.h>
alpar@461: 
alpar@461: namespace lemon {
alpar@461: 
alpar@461:   // GlpkBase members
alpar@461: 
alpar@461:   GlpkBase::GlpkBase() : LpBase() {
alpar@461:     lp = glp_create_prob();
alpar@461:     glp_create_index(lp);
deba@576:     messageLevel(MESSAGE_NOTHING);
alpar@461:   }
alpar@461: 
alpar@461:   GlpkBase::GlpkBase(const GlpkBase &other) : LpBase() {
alpar@461:     lp = glp_create_prob();
alpar@461:     glp_copy_prob(lp, other.lp, GLP_ON);
alpar@461:     glp_create_index(lp);
alpar@461:     rows = other.rows;
alpar@461:     cols = other.cols;
deba@576:     messageLevel(MESSAGE_NOTHING);
alpar@461:   }
alpar@461: 
alpar@461:   GlpkBase::~GlpkBase() {
alpar@461:     glp_delete_prob(lp);
alpar@461:   }
alpar@461: 
alpar@461:   int GlpkBase::_addCol() {
alpar@461:     int i = glp_add_cols(lp, 1);
alpar@461:     glp_set_col_bnds(lp, i, GLP_FR, 0.0, 0.0);
alpar@461:     return i;
alpar@461:   }
alpar@461: 
alpar@461:   int GlpkBase::_addRow() {
alpar@461:     int i = glp_add_rows(lp, 1);
alpar@461:     glp_set_row_bnds(lp, i, GLP_FR, 0.0, 0.0);
alpar@461:     return i;
alpar@461:   }
alpar@461: 
alpar@877:   int GlpkBase::_addRow(Value lo, ExprIterator b,
deba@746:                         ExprIterator e, Value up) {
deba@746:     int i = glp_add_rows(lp, 1);
deba@746: 
deba@746:     if (lo == -INF) {
deba@746:       if (up == INF) {
deba@746:         glp_set_row_bnds(lp, i, GLP_FR, lo, up);
deba@746:       } else {
deba@746:         glp_set_row_bnds(lp, i, GLP_UP, lo, up);
alpar@877:       }
deba@746:     } else {
deba@746:       if (up == INF) {
deba@746:         glp_set_row_bnds(lp, i, GLP_LO, lo, up);
alpar@877:       } else if (lo != up) {
deba@746:         glp_set_row_bnds(lp, i, GLP_DB, lo, up);
deba@746:       } else {
deba@746:         glp_set_row_bnds(lp, i, GLP_FX, lo, up);
deba@746:       }
deba@746:     }
deba@746: 
deba@746:     std::vector<int> indexes;
deba@746:     std::vector<Value> values;
deba@746: 
deba@746:     indexes.push_back(0);
deba@746:     values.push_back(0);
deba@746: 
deba@746:     for(ExprIterator it = b; it != e; ++it) {
deba@746:       indexes.push_back(it->first);
deba@746:       values.push_back(it->second);
deba@746:     }
deba@746: 
deba@746:     glp_set_mat_row(lp, i, values.size() - 1,
deba@746:                     &indexes.front(), &values.front());
deba@746:     return i;
deba@746:   }
deba@746: 
alpar@461:   void GlpkBase::_eraseCol(int i) {
alpar@461:     int ca[2];
alpar@461:     ca[1] = i;
alpar@461:     glp_del_cols(lp, 1, ca);
alpar@461:   }
alpar@461: 
alpar@461:   void GlpkBase::_eraseRow(int i) {
alpar@461:     int ra[2];
alpar@461:     ra[1] = i;
alpar@461:     glp_del_rows(lp, 1, ra);
alpar@461:   }
alpar@461: 
alpar@461:   void GlpkBase::_eraseColId(int i) {
alpar@461:     cols.eraseIndex(i);
alpar@461:     cols.shiftIndices(i);
alpar@461:   }
alpar@461: 
alpar@461:   void GlpkBase::_eraseRowId(int i) {
alpar@461:     rows.eraseIndex(i);
alpar@461:     rows.shiftIndices(i);
alpar@461:   }
alpar@461: 
alpar@461:   void GlpkBase::_getColName(int c, std::string& name) const {
alpar@461:     const char *str = glp_get_col_name(lp, c);
alpar@461:     if (str) name = str;
alpar@461:     else name.clear();
alpar@461:   }
alpar@461: 
alpar@461:   void GlpkBase::_setColName(int c, const std::string & name) {
alpar@461:     glp_set_col_name(lp, c, const_cast<char*>(name.c_str()));
alpar@461: 
alpar@461:   }
alpar@461: 
alpar@461:   int GlpkBase::_colByName(const std::string& name) const {
alpar@461:     int k = glp_find_col(lp, const_cast<char*>(name.c_str()));
alpar@461:     return k > 0 ? k : -1;
alpar@461:   }
alpar@461: 
alpar@461:   void GlpkBase::_getRowName(int r, std::string& name) const {
alpar@461:     const char *str = glp_get_row_name(lp, r);
alpar@461:     if (str) name = str;
alpar@461:     else name.clear();
alpar@461:   }
alpar@461: 
alpar@461:   void GlpkBase::_setRowName(int r, const std::string & name) {
alpar@461:     glp_set_row_name(lp, r, const_cast<char*>(name.c_str()));
alpar@461: 
alpar@461:   }
alpar@461: 
alpar@461:   int GlpkBase::_rowByName(const std::string& name) const {
alpar@461:     int k = glp_find_row(lp, const_cast<char*>(name.c_str()));
alpar@461:     return k > 0 ? k : -1;
alpar@461:   }
alpar@461: 
alpar@461:   void GlpkBase::_setRowCoeffs(int i, ExprIterator b, ExprIterator e) {
alpar@461:     std::vector<int> indexes;
alpar@461:     std::vector<Value> values;
alpar@461: 
alpar@461:     indexes.push_back(0);
alpar@461:     values.push_back(0);
alpar@461: 
alpar@461:     for(ExprIterator it = b; it != e; ++it) {
alpar@461:       indexes.push_back(it->first);
alpar@461:       values.push_back(it->second);
alpar@461:     }
alpar@461: 
alpar@461:     glp_set_mat_row(lp, i, values.size() - 1,
alpar@461:                     &indexes.front(), &values.front());
alpar@461:   }
alpar@461: 
alpar@461:   void GlpkBase::_getRowCoeffs(int ix, InsertIterator b) const {
alpar@461:     int length = glp_get_mat_row(lp, ix, 0, 0);
alpar@461: 
alpar@461:     std::vector<int> indexes(length + 1);
alpar@461:     std::vector<Value> values(length + 1);
alpar@461: 
alpar@461:     glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
alpar@461: 
alpar@461:     for (int i = 1; i <= length; ++i) {
alpar@461:       *b = std::make_pair(indexes[i], values[i]);
alpar@461:       ++b;
alpar@461:     }
alpar@461:   }
alpar@461: 
alpar@461:   void GlpkBase::_setColCoeffs(int ix, ExprIterator b,
alpar@461:                                      ExprIterator e) {
alpar@461: 
alpar@461:     std::vector<int> indexes;
alpar@461:     std::vector<Value> values;
alpar@461: 
alpar@461:     indexes.push_back(0);
alpar@461:     values.push_back(0);
alpar@461: 
alpar@461:     for(ExprIterator it = b; it != e; ++it) {
alpar@461:       indexes.push_back(it->first);
alpar@461:       values.push_back(it->second);
alpar@461:     }
alpar@461: 
alpar@461:     glp_set_mat_col(lp, ix, values.size() - 1,
alpar@461:                     &indexes.front(), &values.front());
alpar@461:   }
alpar@461: 
alpar@461:   void GlpkBase::_getColCoeffs(int ix, InsertIterator b) const {
alpar@461:     int length = glp_get_mat_col(lp, ix, 0, 0);
alpar@461: 
alpar@461:     std::vector<int> indexes(length + 1);
alpar@461:     std::vector<Value> values(length + 1);
alpar@461: 
alpar@461:     glp_get_mat_col(lp, ix, &indexes.front(), &values.front());
alpar@461: 
alpar@461:     for (int i = 1; i  <= length; ++i) {
alpar@461:       *b = std::make_pair(indexes[i], values[i]);
alpar@461:       ++b;
alpar@461:     }
alpar@461:   }
alpar@461: 
alpar@461:   void GlpkBase::_setCoeff(int ix, int jx, Value value) {
alpar@461: 
alpar@461:     if (glp_get_num_cols(lp) < glp_get_num_rows(lp)) {
alpar@461: 
alpar@461:       int length = glp_get_mat_row(lp, ix, 0, 0);
alpar@461: 
alpar@461:       std::vector<int> indexes(length + 2);
alpar@461:       std::vector<Value> values(length + 2);
alpar@461: 
alpar@461:       glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
alpar@461: 
alpar@461:       //The following code does not suppose that the elements of the
alpar@461:       //array indexes are sorted
alpar@461:       bool found = false;
alpar@461:       for (int i = 1; i  <= length; ++i) {
alpar@461:         if (indexes[i] == jx) {
alpar@461:           found = true;
alpar@461:           values[i] = value;
alpar@461:           break;
alpar@461:         }
alpar@461:       }
alpar@461:       if (!found) {
alpar@461:         ++length;
alpar@461:         indexes[length] = jx;
alpar@461:         values[length] = value;
alpar@461:       }
alpar@461: 
alpar@461:       glp_set_mat_row(lp, ix, length, &indexes.front(), &values.front());
alpar@461: 
alpar@461:     } else {
alpar@461: 
alpar@461:       int length = glp_get_mat_col(lp, jx, 0, 0);
alpar@461: 
alpar@461:       std::vector<int> indexes(length + 2);
alpar@461:       std::vector<Value> values(length + 2);
alpar@461: 
alpar@461:       glp_get_mat_col(lp, jx, &indexes.front(), &values.front());
alpar@461: 
alpar@461:       //The following code does not suppose that the elements of the
alpar@461:       //array indexes are sorted
alpar@461:       bool found = false;
alpar@461:       for (int i = 1; i <= length; ++i) {
alpar@461:         if (indexes[i] == ix) {
alpar@461:           found = true;
alpar@461:           values[i] = value;
alpar@461:           break;
alpar@461:         }
alpar@461:       }
alpar@461:       if (!found) {
alpar@461:         ++length;
alpar@461:         indexes[length] = ix;
alpar@461:         values[length] = value;
alpar@461:       }
alpar@461: 
alpar@461:       glp_set_mat_col(lp, jx, length, &indexes.front(), &values.front());
alpar@461:     }
alpar@461: 
alpar@461:   }
alpar@461: 
alpar@461:   GlpkBase::Value GlpkBase::_getCoeff(int ix, int jx) const {
alpar@461: 
alpar@461:     int length = glp_get_mat_row(lp, ix, 0, 0);
alpar@461: 
alpar@461:     std::vector<int> indexes(length + 1);
alpar@461:     std::vector<Value> values(length + 1);
alpar@461: 
alpar@461:     glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
alpar@461: 
alpar@461:     for (int i = 1; i  <= length; ++i) {
alpar@461:       if (indexes[i] == jx) {
alpar@461:         return values[i];
alpar@461:       }
alpar@461:     }
alpar@461: 
alpar@461:     return 0;
alpar@461:   }
alpar@461: 
alpar@461:   void GlpkBase::_setColLowerBound(int i, Value lo) {
alpar@461:     LEMON_ASSERT(lo != INF, "Invalid bound");
alpar@461: 
alpar@461:     int b = glp_get_col_type(lp, i);
alpar@461:     double up = glp_get_col_ub(lp, i);
alpar@461:     if (lo == -INF) {
alpar@461:       switch (b) {
alpar@461:       case GLP_FR:
alpar@461:       case GLP_LO:
alpar@461:         glp_set_col_bnds(lp, i, GLP_FR, lo, up);
alpar@461:         break;
alpar@461:       case GLP_UP:
alpar@461:         break;
alpar@461:       case GLP_DB:
alpar@461:       case GLP_FX:
alpar@461:         glp_set_col_bnds(lp, i, GLP_UP, lo, up);
alpar@461:         break;
alpar@461:       default:
alpar@461:         break;
alpar@461:       }
alpar@461:     } else {
alpar@461:       switch (b) {
alpar@461:       case GLP_FR:
alpar@461:       case GLP_LO:
alpar@461:         glp_set_col_bnds(lp, i, GLP_LO, lo, up);
alpar@461:         break;
alpar@461:       case GLP_UP:
alpar@461:       case GLP_DB:
alpar@461:       case GLP_FX:
alpar@461:         if (lo == up)
alpar@461:           glp_set_col_bnds(lp, i, GLP_FX, lo, up);
alpar@461:         else
alpar@461:           glp_set_col_bnds(lp, i, GLP_DB, lo, up);
alpar@461:         break;
alpar@461:       default:
alpar@461:         break;
alpar@461:       }
alpar@461:     }
alpar@461:   }
alpar@461: 
alpar@461:   GlpkBase::Value GlpkBase::_getColLowerBound(int i) const {
alpar@461:     int b = glp_get_col_type(lp, i);
alpar@461:     switch (b) {
alpar@461:     case GLP_LO:
alpar@461:     case GLP_DB:
alpar@461:     case GLP_FX:
alpar@461:       return glp_get_col_lb(lp, i);
alpar@461:     default:
alpar@461:       return -INF;
alpar@461:     }
alpar@461:   }
alpar@461: 
alpar@461:   void GlpkBase::_setColUpperBound(int i, Value up) {
alpar@461:     LEMON_ASSERT(up != -INF, "Invalid bound");
alpar@461: 
alpar@461:     int b = glp_get_col_type(lp, i);
alpar@461:     double lo = glp_get_col_lb(lp, i);
alpar@461:     if (up == INF) {
alpar@461:       switch (b) {
alpar@461:       case GLP_FR:
alpar@461:       case GLP_LO:
alpar@461:         break;
alpar@461:       case GLP_UP:
alpar@461:         glp_set_col_bnds(lp, i, GLP_FR, lo, up);
alpar@461:         break;
alpar@461:       case GLP_DB:
alpar@461:       case GLP_FX:
alpar@461:         glp_set_col_bnds(lp, i, GLP_LO, lo, up);
alpar@461:         break;
alpar@461:       default:
alpar@461:         break;
alpar@461:       }
alpar@461:     } else {
alpar@461:       switch (b) {
alpar@461:       case GLP_FR:
alpar@461:         glp_set_col_bnds(lp, i, GLP_UP, lo, up);
alpar@461:         break;
alpar@461:       case GLP_UP:
alpar@461:         glp_set_col_bnds(lp, i, GLP_UP, lo, up);
alpar@461:         break;
alpar@461:       case GLP_LO:
alpar@461:       case GLP_DB:
alpar@461:       case GLP_FX:
alpar@461:         if (lo == up)
alpar@461:           glp_set_col_bnds(lp, i, GLP_FX, lo, up);
alpar@461:         else
alpar@461:           glp_set_col_bnds(lp, i, GLP_DB, lo, up);
alpar@461:         break;
alpar@461:       default:
alpar@461:         break;
alpar@461:       }
alpar@461:     }
alpar@461: 
alpar@461:   }
alpar@461: 
alpar@461:   GlpkBase::Value GlpkBase::_getColUpperBound(int i) const {
alpar@461:     int b = glp_get_col_type(lp, i);
alpar@461:       switch (b) {
alpar@461:       case GLP_UP:
alpar@461:       case GLP_DB:
alpar@461:       case GLP_FX:
alpar@461:         return glp_get_col_ub(lp, i);
alpar@461:       default:
alpar@461:         return INF;
alpar@461:       }
alpar@461:   }
alpar@461: 
alpar@461:   void GlpkBase::_setRowLowerBound(int i, Value lo) {
alpar@461:     LEMON_ASSERT(lo != INF, "Invalid bound");
alpar@461: 
alpar@461:     int b = glp_get_row_type(lp, i);
alpar@461:     double up = glp_get_row_ub(lp, i);
alpar@461:     if (lo == -INF) {
alpar@461:       switch (b) {
alpar@461:       case GLP_FR:
alpar@461:       case GLP_LO:
alpar@461:         glp_set_row_bnds(lp, i, GLP_FR, lo, up);
alpar@461:         break;
alpar@461:       case GLP_UP:
alpar@461:         break;
alpar@461:       case GLP_DB:
alpar@461:       case GLP_FX:
alpar@461:         glp_set_row_bnds(lp, i, GLP_UP, lo, up);
alpar@461:         break;
alpar@461:       default:
alpar@461:         break;
alpar@461:       }
alpar@461:     } else {
alpar@461:       switch (b) {
alpar@461:       case GLP_FR:
alpar@461:       case GLP_LO:
alpar@461:         glp_set_row_bnds(lp, i, GLP_LO, lo, up);
alpar@461:         break;
alpar@461:       case GLP_UP:
alpar@461:       case GLP_DB:
alpar@461:       case GLP_FX:
alpar@461:         if (lo == up)
alpar@461:           glp_set_row_bnds(lp, i, GLP_FX, lo, up);
alpar@461:         else
alpar@461:           glp_set_row_bnds(lp, i, GLP_DB, lo, up);
alpar@461:         break;
alpar@461:       default:
alpar@461:         break;
alpar@461:       }
alpar@461:     }
alpar@461: 
alpar@461:   }
alpar@461: 
alpar@461:   GlpkBase::Value GlpkBase::_getRowLowerBound(int i) const {
alpar@461:     int b = glp_get_row_type(lp, i);
alpar@461:     switch (b) {
alpar@461:     case GLP_LO:
alpar@461:     case GLP_DB:
alpar@461:     case GLP_FX:
alpar@461:       return glp_get_row_lb(lp, i);
alpar@461:     default:
alpar@461:       return -INF;
alpar@461:     }
alpar@461:   }
alpar@461: 
alpar@461:   void GlpkBase::_setRowUpperBound(int i, Value up) {
alpar@461:     LEMON_ASSERT(up != -INF, "Invalid bound");
alpar@461: 
alpar@461:     int b = glp_get_row_type(lp, i);
alpar@461:     double lo = glp_get_row_lb(lp, i);
alpar@461:     if (up == INF) {
alpar@461:       switch (b) {
alpar@461:       case GLP_FR:
alpar@461:       case GLP_LO:
alpar@461:         break;
alpar@461:       case GLP_UP:
alpar@461:         glp_set_row_bnds(lp, i, GLP_FR, lo, up);
alpar@461:         break;
alpar@461:       case GLP_DB:
alpar@461:       case GLP_FX:
alpar@461:         glp_set_row_bnds(lp, i, GLP_LO, lo, up);
alpar@461:         break;
alpar@461:       default:
alpar@461:         break;
alpar@461:       }
alpar@461:     } else {
alpar@461:       switch (b) {
alpar@461:       case GLP_FR:
alpar@461:         glp_set_row_bnds(lp, i, GLP_UP, lo, up);
alpar@461:         break;
alpar@461:       case GLP_UP:
alpar@461:         glp_set_row_bnds(lp, i, GLP_UP, lo, up);
alpar@461:         break;
alpar@461:       case GLP_LO:
alpar@461:       case GLP_DB:
alpar@461:       case GLP_FX:
alpar@461:         if (lo == up)
alpar@461:           glp_set_row_bnds(lp, i, GLP_FX, lo, up);
alpar@461:         else
alpar@461:           glp_set_row_bnds(lp, i, GLP_DB, lo, up);
alpar@461:         break;
alpar@461:       default:
alpar@461:         break;
alpar@461:       }
alpar@461:     }
alpar@461:   }
alpar@461: 
alpar@461:   GlpkBase::Value GlpkBase::_getRowUpperBound(int i) const {
alpar@461:     int b = glp_get_row_type(lp, i);
alpar@461:     switch (b) {
alpar@461:     case GLP_UP:
alpar@461:     case GLP_DB:
alpar@461:     case GLP_FX:
alpar@461:       return glp_get_row_ub(lp, i);
alpar@461:     default:
alpar@461:       return INF;
alpar@461:     }
alpar@461:   }
alpar@461: 
alpar@461:   void GlpkBase::_setObjCoeffs(ExprIterator b, ExprIterator e) {
alpar@461:     for (int i = 1; i <= glp_get_num_cols(lp); ++i) {
alpar@461:       glp_set_obj_coef(lp, i, 0.0);
alpar@461:     }
alpar@461:     for (ExprIterator it = b; it != e; ++it) {
alpar@461:       glp_set_obj_coef(lp, it->first, it->second);
alpar@461:     }
alpar@461:   }
alpar@461: 
alpar@461:   void GlpkBase::_getObjCoeffs(InsertIterator b) const {
alpar@461:     for (int i = 1; i <= glp_get_num_cols(lp); ++i) {
alpar@461:       Value val = glp_get_obj_coef(lp, i);
alpar@461:       if (val != 0.0) {
alpar@461:         *b = std::make_pair(i, val);
alpar@461:         ++b;
alpar@461:       }
alpar@461:     }
alpar@461:   }
alpar@461: 
alpar@461:   void GlpkBase::_setObjCoeff(int i, Value obj_coef) {
alpar@461:     //i = 0 means the constant term (shift)
alpar@461:     glp_set_obj_coef(lp, i, obj_coef);
alpar@461:   }
alpar@461: 
alpar@461:   GlpkBase::Value GlpkBase::_getObjCoeff(int i) const {
alpar@461:     //i = 0 means the constant term (shift)
alpar@461:     return glp_get_obj_coef(lp, i);
alpar@461:   }
alpar@461: 
alpar@461:   void GlpkBase::_setSense(GlpkBase::Sense sense) {
alpar@461:     switch (sense) {
alpar@461:     case MIN:
alpar@461:       glp_set_obj_dir(lp, GLP_MIN);
alpar@461:       break;
alpar@461:     case MAX:
alpar@461:       glp_set_obj_dir(lp, GLP_MAX);
alpar@461:       break;
alpar@461:     }
alpar@461:   }
alpar@461: 
alpar@461:   GlpkBase::Sense GlpkBase::_getSense() const {
alpar@461:     switch(glp_get_obj_dir(lp)) {
alpar@461:     case GLP_MIN:
alpar@461:       return MIN;
alpar@461:     case GLP_MAX:
alpar@461:       return MAX;
alpar@461:     default:
alpar@461:       LEMON_ASSERT(false, "Wrong sense");
alpar@461:       return GlpkBase::Sense();
alpar@461:     }
alpar@461:   }
alpar@461: 
alpar@461:   void GlpkBase::_clear() {
alpar@461:     glp_erase_prob(lp);
alpar@461:   }
alpar@461: 
deba@537:   void GlpkBase::freeEnv() {
deba@537:     glp_free_env();
deba@537:   }
deba@537: 
deba@576:   void GlpkBase::_messageLevel(MessageLevel level) {
deba@576:     switch (level) {
deba@576:     case MESSAGE_NOTHING:
deba@576:       _message_level = GLP_MSG_OFF;
deba@576:       break;
deba@576:     case MESSAGE_ERROR:
deba@576:       _message_level = GLP_MSG_ERR;
deba@576:       break;
deba@576:     case MESSAGE_WARNING:
deba@576:       _message_level = GLP_MSG_ERR;
deba@576:       break;
deba@576:     case MESSAGE_NORMAL:
deba@576:       _message_level = GLP_MSG_ON;
deba@576:       break;
deba@576:     case MESSAGE_VERBOSE:
deba@576:       _message_level = GLP_MSG_ALL;
deba@576:       break;
deba@576:     }
deba@576:   }
deba@576: 
deba@538:   GlpkBase::FreeEnvHelper GlpkBase::freeEnvHelper;
deba@538: 
alpar@462:   // GlpkLp members
alpar@461: 
alpar@462:   GlpkLp::GlpkLp()
deba@551:     : LpBase(), LpSolver(), GlpkBase() {
deba@565:     presolver(false);
alpar@461:   }
alpar@461: 
alpar@462:   GlpkLp::GlpkLp(const GlpkLp& other)
deba@551:     : LpBase(other), LpSolver(other), GlpkBase(other) {
deba@565:     presolver(false);
alpar@461:   }
alpar@461: 
alpar@540:   GlpkLp* GlpkLp::newSolver() const { return new GlpkLp; }
alpar@540:   GlpkLp* GlpkLp::cloneSolver() const { return new GlpkLp(*this); }
alpar@461: 
alpar@462:   const char* GlpkLp::_solverName() const { return "GlpkLp"; }
alpar@461: 
alpar@462:   void GlpkLp::_clear_temporals() {
alpar@461:     _primal_ray.clear();
alpar@461:     _dual_ray.clear();
alpar@461:   }
alpar@461: 
alpar@462:   GlpkLp::SolveExitStatus GlpkLp::_solve() {
alpar@461:     return solvePrimal();
alpar@461:   }
alpar@461: 
alpar@462:   GlpkLp::SolveExitStatus GlpkLp::solvePrimal() {
alpar@461:     _clear_temporals();
alpar@461: 
alpar@461:     glp_smcp smcp;
alpar@461:     glp_init_smcp(&smcp);
alpar@461: 
deba@576:     smcp.msg_lev = _message_level;
deba@565:     smcp.presolve = _presolve;
alpar@461: 
deba@565:     // If the basis is not valid we get an error return value.
deba@565:     // In this case we can try to create a new basis.
deba@565:     switch (glp_simplex(lp, &smcp)) {
deba@565:     case 0:
deba@565:       break;
deba@565:     case GLP_EBADB:
deba@565:     case GLP_ESING:
deba@565:     case GLP_ECOND:
deba@566:       glp_term_out(false);
deba@565:       glp_adv_basis(lp, 0);
deba@566:       glp_term_out(true);
deba@565:       if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
deba@565:       break;
deba@565:     default:
deba@565:       return UNSOLVED;
deba@565:     }
deba@565: 
alpar@461:     return SOLVED;
alpar@461:   }
alpar@461: 
alpar@462:   GlpkLp::SolveExitStatus GlpkLp::solveDual() {
alpar@461:     _clear_temporals();
alpar@461: 
alpar@461:     glp_smcp smcp;
alpar@461:     glp_init_smcp(&smcp);
alpar@461: 
deba@576:     smcp.msg_lev = _message_level;
alpar@461:     smcp.meth = GLP_DUAL;
deba@565:     smcp.presolve = _presolve;
alpar@461: 
deba@565:     // If the basis is not valid we get an error return value.
deba@565:     // In this case we can try to create a new basis.
deba@565:     switch (glp_simplex(lp, &smcp)) {
deba@565:     case 0:
deba@565:       break;
deba@565:     case GLP_EBADB:
deba@565:     case GLP_ESING:
deba@565:     case GLP_ECOND:
deba@566:       glp_term_out(false);
deba@565:       glp_adv_basis(lp, 0);
deba@566:       glp_term_out(true);
deba@565:       if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
deba@565:       break;
deba@565:     default:
deba@565:       return UNSOLVED;
deba@565:     }
alpar@461:     return SOLVED;
alpar@461:   }
alpar@461: 
alpar@462:   GlpkLp::Value GlpkLp::_getPrimal(int i) const {
alpar@461:     return glp_get_col_prim(lp, i);
alpar@461:   }
alpar@461: 
alpar@462:   GlpkLp::Value GlpkLp::_getDual(int i) const {
alpar@461:     return glp_get_row_dual(lp, i);
alpar@461:   }
alpar@461: 
alpar@462:   GlpkLp::Value GlpkLp::_getPrimalValue() const {
alpar@461:     return glp_get_obj_val(lp);
alpar@461:   }
alpar@461: 
alpar@462:   GlpkLp::VarStatus GlpkLp::_getColStatus(int i) const {
alpar@461:     switch (glp_get_col_stat(lp, i)) {
alpar@461:     case GLP_BS:
alpar@461:       return BASIC;
alpar@461:     case GLP_UP:
alpar@461:       return UPPER;
alpar@461:     case GLP_LO:
alpar@461:       return LOWER;
alpar@461:     case GLP_NF:
alpar@461:       return FREE;
alpar@461:     case GLP_NS:
alpar@461:       return FIXED;
alpar@461:     default:
alpar@461:       LEMON_ASSERT(false, "Wrong column status");
alpar@462:       return GlpkLp::VarStatus();
alpar@461:     }
alpar@461:   }
alpar@461: 
alpar@462:   GlpkLp::VarStatus GlpkLp::_getRowStatus(int i) const {
alpar@461:     switch (glp_get_row_stat(lp, i)) {
alpar@461:     case GLP_BS:
alpar@461:       return BASIC;
alpar@461:     case GLP_UP:
alpar@461:       return UPPER;
alpar@461:     case GLP_LO:
alpar@461:       return LOWER;
alpar@461:     case GLP_NF:
alpar@461:       return FREE;
alpar@461:     case GLP_NS:
alpar@461:       return FIXED;
alpar@461:     default:
alpar@461:       LEMON_ASSERT(false, "Wrong row status");
alpar@462:       return GlpkLp::VarStatus();
alpar@461:     }
alpar@461:   }
alpar@461: 
alpar@462:   GlpkLp::Value GlpkLp::_getPrimalRay(int i) const {
alpar@461:     if (_primal_ray.empty()) {
alpar@461:       int row_num = glp_get_num_rows(lp);
alpar@461:       int col_num = glp_get_num_cols(lp);
alpar@461: 
alpar@461:       _primal_ray.resize(col_num + 1, 0.0);
alpar@461: 
alpar@461:       int index = glp_get_unbnd_ray(lp);
alpar@461:       if (index != 0) {
alpar@461:         // The primal ray is found in primal simplex second phase
alpar@461:         LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :
alpar@461:                       glp_get_col_stat(lp, index - row_num)) != GLP_BS,
alpar@461:                      "Wrong primal ray");
alpar@461: 
alpar@461:         bool negate = glp_get_obj_dir(lp) == GLP_MAX;
alpar@461: 
alpar@461:         if (index > row_num) {
alpar@461:           _primal_ray[index - row_num] = 1.0;
alpar@461:           if (glp_get_col_dual(lp, index - row_num) > 0) {
alpar@461:             negate = !negate;
alpar@461:           }
alpar@461:         } else {
alpar@461:           if (glp_get_row_dual(lp, index) > 0) {
alpar@461:             negate = !negate;
alpar@461:           }
alpar@461:         }
alpar@461: 
alpar@461:         std::vector<int> ray_indexes(row_num + 1);
alpar@461:         std::vector<Value> ray_values(row_num + 1);
alpar@461:         int ray_length = glp_eval_tab_col(lp, index, &ray_indexes.front(),
alpar@461:                                           &ray_values.front());
alpar@461: 
alpar@461:         for (int i = 1; i <= ray_length; ++i) {
alpar@461:           if (ray_indexes[i] > row_num) {
alpar@461:             _primal_ray[ray_indexes[i] - row_num] = ray_values[i];
alpar@461:           }
alpar@461:         }
alpar@461: 
alpar@461:         if (negate) {
alpar@461:           for (int i = 1; i <= col_num; ++i) {
alpar@461:             _primal_ray[i] = - _primal_ray[i];
alpar@461:           }
alpar@461:         }
alpar@461:       } else {
alpar@461:         for (int i = 1; i <= col_num; ++i) {
alpar@461:           _primal_ray[i] = glp_get_col_prim(lp, i);
alpar@461:         }
alpar@461:       }
alpar@461:     }
alpar@461:     return _primal_ray[i];
alpar@461:   }
alpar@461: 
alpar@462:   GlpkLp::Value GlpkLp::_getDualRay(int i) const {
alpar@461:     if (_dual_ray.empty()) {
alpar@461:       int row_num = glp_get_num_rows(lp);
alpar@461: 
alpar@461:       _dual_ray.resize(row_num + 1, 0.0);
alpar@461: 
alpar@461:       int index = glp_get_unbnd_ray(lp);
alpar@461:       if (index != 0) {
alpar@461:         // The dual ray is found in dual simplex second phase
alpar@461:         LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :
alpar@461:                       glp_get_col_stat(lp, index - row_num)) == GLP_BS,
alpar@461: 
alpar@461:                      "Wrong dual ray");
alpar@461: 
alpar@461:         int idx;
alpar@461:         bool negate = false;
alpar@461: 
alpar@461:         if (index > row_num) {
alpar@461:           idx = glp_get_col_bind(lp, index - row_num);
alpar@461:           if (glp_get_col_prim(lp, index - row_num) >
alpar@461:               glp_get_col_ub(lp, index - row_num)) {
alpar@461:             negate = true;
alpar@461:           }
alpar@461:         } else {
alpar@461:           idx = glp_get_row_bind(lp, index);
alpar@461:           if (glp_get_row_prim(lp, index) > glp_get_row_ub(lp, index)) {
alpar@461:             negate = true;
alpar@461:           }
alpar@461:         }
alpar@461: 
alpar@461:         _dual_ray[idx] = negate ?  - 1.0 : 1.0;
alpar@461: 
alpar@461:         glp_btran(lp, &_dual_ray.front());
alpar@461:       } else {
alpar@461:         double eps = 1e-7;
alpar@461:         // The dual ray is found in primal simplex first phase
alpar@461:         // We assume that the glpk minimizes the slack to get feasible solution
alpar@461:         for (int i = 1; i <= row_num; ++i) {
alpar@461:           int index = glp_get_bhead(lp, i);
alpar@461:           if (index <= row_num) {
alpar@461:             double res = glp_get_row_prim(lp, index);
alpar@461:             if (res > glp_get_row_ub(lp, index) + eps) {
alpar@461:               _dual_ray[i] = -1;
alpar@461:             } else if (res < glp_get_row_lb(lp, index) - eps) {
alpar@461:               _dual_ray[i] = 1;
alpar@461:             } else {
alpar@461:               _dual_ray[i] = 0;
alpar@461:             }
alpar@461:             _dual_ray[i] *= glp_get_rii(lp, index);
alpar@461:           } else {
alpar@461:             double res = glp_get_col_prim(lp, index - row_num);
alpar@461:             if (res > glp_get_col_ub(lp, index - row_num) + eps) {
alpar@461:               _dual_ray[i] = -1;
alpar@461:             } else if (res < glp_get_col_lb(lp, index - row_num) - eps) {
alpar@461:               _dual_ray[i] = 1;
alpar@461:             } else {
alpar@461:               _dual_ray[i] = 0;
alpar@461:             }
alpar@461:             _dual_ray[i] /= glp_get_sjj(lp, index - row_num);
alpar@461:           }
alpar@461:         }
alpar@461: 
alpar@461:         glp_btran(lp, &_dual_ray.front());
alpar@461: 
alpar@461:         for (int i = 1; i <= row_num; ++i) {
alpar@461:           _dual_ray[i] /= glp_get_rii(lp, i);
alpar@461:         }
alpar@461:       }
alpar@461:     }
alpar@461:     return _dual_ray[i];
alpar@461:   }
alpar@461: 
alpar@462:   GlpkLp::ProblemType GlpkLp::_getPrimalType() const {
alpar@461:     if (glp_get_status(lp) == GLP_OPT)
alpar@461:       return OPTIMAL;
alpar@461:     switch (glp_get_prim_stat(lp)) {
alpar@461:     case GLP_UNDEF:
alpar@461:       return UNDEFINED;
alpar@461:     case GLP_FEAS:
alpar@461:     case GLP_INFEAS:
alpar@461:       if (glp_get_dual_stat(lp) == GLP_NOFEAS) {
alpar@461:         return UNBOUNDED;
alpar@461:       } else {
alpar@461:         return UNDEFINED;
alpar@461:       }
alpar@461:     case GLP_NOFEAS:
alpar@461:       return INFEASIBLE;
alpar@461:     default:
alpar@461:       LEMON_ASSERT(false, "Wrong primal type");
alpar@462:       return  GlpkLp::ProblemType();
alpar@461:     }
alpar@461:   }
alpar@461: 
alpar@462:   GlpkLp::ProblemType GlpkLp::_getDualType() const {
alpar@461:     if (glp_get_status(lp) == GLP_OPT)
alpar@461:       return OPTIMAL;
alpar@461:     switch (glp_get_dual_stat(lp)) {
alpar@461:     case GLP_UNDEF:
alpar@461:       return UNDEFINED;
alpar@461:     case GLP_FEAS:
alpar@461:     case GLP_INFEAS:
alpar@461:       if (glp_get_prim_stat(lp) == GLP_NOFEAS) {
alpar@461:         return UNBOUNDED;
alpar@461:       } else {
alpar@461:         return UNDEFINED;
alpar@461:       }
alpar@461:     case GLP_NOFEAS:
alpar@461:       return INFEASIBLE;
alpar@461:     default:
alpar@461:       LEMON_ASSERT(false, "Wrong primal type");
alpar@462:       return  GlpkLp::ProblemType();
alpar@461:     }
alpar@461:   }
alpar@461: 
deba@565:   void GlpkLp::presolver(bool presolve) {
deba@565:     _presolve = presolve;
alpar@461:   }
alpar@461: 
alpar@462:   // GlpkMip members
alpar@461: 
alpar@462:   GlpkMip::GlpkMip()
deba@551:     : LpBase(), MipSolver(), GlpkBase() {
alpar@461:   }
alpar@461: 
alpar@462:   GlpkMip::GlpkMip(const GlpkMip& other)
deba@551:     : LpBase(), MipSolver(), GlpkBase(other) {
alpar@461:   }
alpar@461: 
alpar@462:   void GlpkMip::_setColType(int i, GlpkMip::ColTypes col_type) {
alpar@461:     switch (col_type) {
alpar@461:     case INTEGER:
alpar@461:       glp_set_col_kind(lp, i, GLP_IV);
alpar@461:       break;
alpar@461:     case REAL:
alpar@461:       glp_set_col_kind(lp, i, GLP_CV);
alpar@461:       break;
alpar@461:     }
alpar@461:   }
alpar@461: 
alpar@462:   GlpkMip::ColTypes GlpkMip::_getColType(int i) const {
alpar@461:     switch (glp_get_col_kind(lp, i)) {
alpar@461:     case GLP_IV:
alpar@461:     case GLP_BV:
alpar@461:       return INTEGER;
alpar@461:     default:
alpar@461:       return REAL;
alpar@461:     }
alpar@461: 
alpar@461:   }
alpar@461: 
alpar@462:   GlpkMip::SolveExitStatus GlpkMip::_solve() {
alpar@461:     glp_smcp smcp;
alpar@461:     glp_init_smcp(&smcp);
alpar@461: 
deba@576:     smcp.msg_lev = _message_level;
alpar@461:     smcp.meth = GLP_DUAL;
alpar@461: 
deba@565:     // If the basis is not valid we get an error return value.
deba@565:     // In this case we can try to create a new basis.
deba@565:     switch (glp_simplex(lp, &smcp)) {
deba@565:     case 0:
deba@565:       break;
deba@565:     case GLP_EBADB:
deba@565:     case GLP_ESING:
deba@565:     case GLP_ECOND:
deba@566:       glp_term_out(false);
deba@565:       glp_adv_basis(lp, 0);
deba@566:       glp_term_out(true);
deba@565:       if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
deba@565:       break;
deba@565:     default:
deba@565:       return UNSOLVED;
deba@565:     }
deba@565: 
alpar@461:     if (glp_get_status(lp) != GLP_OPT) return SOLVED;
alpar@461: 
alpar@461:     glp_iocp iocp;
alpar@461:     glp_init_iocp(&iocp);
alpar@461: 
deba@576:     iocp.msg_lev = _message_level;
alpar@461: 
alpar@461:     if (glp_intopt(lp, &iocp) != 0) return UNSOLVED;
alpar@461:     return SOLVED;
alpar@461:   }
alpar@461: 
alpar@461: 
alpar@462:   GlpkMip::ProblemType GlpkMip::_getType() const {
alpar@461:     switch (glp_get_status(lp)) {
alpar@461:     case GLP_OPT:
alpar@461:       switch (glp_mip_status(lp)) {
alpar@461:       case GLP_UNDEF:
alpar@461:         return UNDEFINED;
alpar@461:       case GLP_NOFEAS:
alpar@461:         return INFEASIBLE;
alpar@461:       case GLP_FEAS:
alpar@461:         return FEASIBLE;
alpar@461:       case GLP_OPT:
alpar@461:         return OPTIMAL;
alpar@461:       default:
alpar@461:         LEMON_ASSERT(false, "Wrong problem type.");
alpar@462:         return GlpkMip::ProblemType();
alpar@461:       }
alpar@461:     case GLP_NOFEAS:
alpar@461:       return INFEASIBLE;
alpar@461:     case GLP_INFEAS:
alpar@461:     case GLP_FEAS:
alpar@461:       if (glp_get_dual_stat(lp) == GLP_NOFEAS) {
alpar@461:         return UNBOUNDED;
alpar@461:       } else {
alpar@461:         return UNDEFINED;
alpar@461:       }
alpar@461:     default:
alpar@461:       LEMON_ASSERT(false, "Wrong problem type.");
alpar@462:       return GlpkMip::ProblemType();
alpar@461:     }
alpar@461:   }
alpar@461: 
alpar@462:   GlpkMip::Value GlpkMip::_getSol(int i) const {
alpar@461:     return glp_mip_col_val(lp, i);
alpar@461:   }
alpar@461: 
alpar@462:   GlpkMip::Value GlpkMip::_getSolValue() const {
alpar@461:     return glp_mip_obj_val(lp);
alpar@461:   }
alpar@461: 
alpar@540:   GlpkMip* GlpkMip::newSolver() const { return new GlpkMip; }
alpar@540:   GlpkMip* GlpkMip::cloneSolver() const {return new GlpkMip(*this); }
alpar@461: 
alpar@462:   const char* GlpkMip::_solverName() const { return "GlpkMip"; }
alpar@461: 
alpar@461: } //END OF NAMESPACE LEMON