1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lemon/glpk.cc Mon Jan 12 12:26:01 2009 +0000
1.3 @@ -0,0 +1,952 @@
1.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
1.5 + *
1.6 + * This file is a part of LEMON, a generic C++ optimization library.
1.7 + *
1.8 + * Copyright (C) 2003-2008
1.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
1.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
1.11 + *
1.12 + * Permission to use, modify and distribute this software is granted
1.13 + * provided that this copyright notice appears in all copies. For
1.14 + * precise terms see the accompanying LICENSE file.
1.15 + *
1.16 + * This software is provided "AS IS" with no warranty of any kind,
1.17 + * express or implied, and with no claim as to its suitability for any
1.18 + * purpose.
1.19 + *
1.20 + */
1.21 +
1.22 +///\file
1.23 +///\brief Implementation of the LEMON GLPK LP and MIP solver interface.
1.24 +
1.25 +#include <lemon/glpk.h>
1.26 +#include <glpk.h>
1.27 +
1.28 +#include <lemon/assert.h>
1.29 +
1.30 +namespace lemon {
1.31 +
1.32 + // GlpkBase members
1.33 +
1.34 + GlpkBase::GlpkBase() : LpBase() {
1.35 + lp = glp_create_prob();
1.36 + glp_create_index(lp);
1.37 + }
1.38 +
1.39 + GlpkBase::GlpkBase(const GlpkBase &other) : LpBase() {
1.40 + lp = glp_create_prob();
1.41 + glp_copy_prob(lp, other.lp, GLP_ON);
1.42 + glp_create_index(lp);
1.43 + rows = other.rows;
1.44 + cols = other.cols;
1.45 + }
1.46 +
1.47 + GlpkBase::~GlpkBase() {
1.48 + glp_delete_prob(lp);
1.49 + }
1.50 +
1.51 + int GlpkBase::_addCol() {
1.52 + int i = glp_add_cols(lp, 1);
1.53 + glp_set_col_bnds(lp, i, GLP_FR, 0.0, 0.0);
1.54 + return i;
1.55 + }
1.56 +
1.57 + int GlpkBase::_addRow() {
1.58 + int i = glp_add_rows(lp, 1);
1.59 + glp_set_row_bnds(lp, i, GLP_FR, 0.0, 0.0);
1.60 + return i;
1.61 + }
1.62 +
1.63 + void GlpkBase::_eraseCol(int i) {
1.64 + int ca[2];
1.65 + ca[1] = i;
1.66 + glp_del_cols(lp, 1, ca);
1.67 + }
1.68 +
1.69 + void GlpkBase::_eraseRow(int i) {
1.70 + int ra[2];
1.71 + ra[1] = i;
1.72 + glp_del_rows(lp, 1, ra);
1.73 + }
1.74 +
1.75 + void GlpkBase::_eraseColId(int i) {
1.76 + cols.eraseIndex(i);
1.77 + cols.shiftIndices(i);
1.78 + }
1.79 +
1.80 + void GlpkBase::_eraseRowId(int i) {
1.81 + rows.eraseIndex(i);
1.82 + rows.shiftIndices(i);
1.83 + }
1.84 +
1.85 + void GlpkBase::_getColName(int c, std::string& name) const {
1.86 + const char *str = glp_get_col_name(lp, c);
1.87 + if (str) name = str;
1.88 + else name.clear();
1.89 + }
1.90 +
1.91 + void GlpkBase::_setColName(int c, const std::string & name) {
1.92 + glp_set_col_name(lp, c, const_cast<char*>(name.c_str()));
1.93 +
1.94 + }
1.95 +
1.96 + int GlpkBase::_colByName(const std::string& name) const {
1.97 + int k = glp_find_col(lp, const_cast<char*>(name.c_str()));
1.98 + return k > 0 ? k : -1;
1.99 + }
1.100 +
1.101 + void GlpkBase::_getRowName(int r, std::string& name) const {
1.102 + const char *str = glp_get_row_name(lp, r);
1.103 + if (str) name = str;
1.104 + else name.clear();
1.105 + }
1.106 +
1.107 + void GlpkBase::_setRowName(int r, const std::string & name) {
1.108 + glp_set_row_name(lp, r, const_cast<char*>(name.c_str()));
1.109 +
1.110 + }
1.111 +
1.112 + int GlpkBase::_rowByName(const std::string& name) const {
1.113 + int k = glp_find_row(lp, const_cast<char*>(name.c_str()));
1.114 + return k > 0 ? k : -1;
1.115 + }
1.116 +
1.117 + void GlpkBase::_setRowCoeffs(int i, ExprIterator b, ExprIterator e) {
1.118 + std::vector<int> indexes;
1.119 + std::vector<Value> values;
1.120 +
1.121 + indexes.push_back(0);
1.122 + values.push_back(0);
1.123 +
1.124 + for(ExprIterator it = b; it != e; ++it) {
1.125 + indexes.push_back(it->first);
1.126 + values.push_back(it->second);
1.127 + }
1.128 +
1.129 + glp_set_mat_row(lp, i, values.size() - 1,
1.130 + &indexes.front(), &values.front());
1.131 + }
1.132 +
1.133 + void GlpkBase::_getRowCoeffs(int ix, InsertIterator b) const {
1.134 + int length = glp_get_mat_row(lp, ix, 0, 0);
1.135 +
1.136 + std::vector<int> indexes(length + 1);
1.137 + std::vector<Value> values(length + 1);
1.138 +
1.139 + glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
1.140 +
1.141 + for (int i = 1; i <= length; ++i) {
1.142 + *b = std::make_pair(indexes[i], values[i]);
1.143 + ++b;
1.144 + }
1.145 + }
1.146 +
1.147 + void GlpkBase::_setColCoeffs(int ix, ExprIterator b,
1.148 + ExprIterator e) {
1.149 +
1.150 + std::vector<int> indexes;
1.151 + std::vector<Value> values;
1.152 +
1.153 + indexes.push_back(0);
1.154 + values.push_back(0);
1.155 +
1.156 + for(ExprIterator it = b; it != e; ++it) {
1.157 + indexes.push_back(it->first);
1.158 + values.push_back(it->second);
1.159 + }
1.160 +
1.161 + glp_set_mat_col(lp, ix, values.size() - 1,
1.162 + &indexes.front(), &values.front());
1.163 + }
1.164 +
1.165 + void GlpkBase::_getColCoeffs(int ix, InsertIterator b) const {
1.166 + int length = glp_get_mat_col(lp, ix, 0, 0);
1.167 +
1.168 + std::vector<int> indexes(length + 1);
1.169 + std::vector<Value> values(length + 1);
1.170 +
1.171 + glp_get_mat_col(lp, ix, &indexes.front(), &values.front());
1.172 +
1.173 + for (int i = 1; i <= length; ++i) {
1.174 + *b = std::make_pair(indexes[i], values[i]);
1.175 + ++b;
1.176 + }
1.177 + }
1.178 +
1.179 + void GlpkBase::_setCoeff(int ix, int jx, Value value) {
1.180 +
1.181 + if (glp_get_num_cols(lp) < glp_get_num_rows(lp)) {
1.182 +
1.183 + int length = glp_get_mat_row(lp, ix, 0, 0);
1.184 +
1.185 + std::vector<int> indexes(length + 2);
1.186 + std::vector<Value> values(length + 2);
1.187 +
1.188 + glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
1.189 +
1.190 + //The following code does not suppose that the elements of the
1.191 + //array indexes are sorted
1.192 + bool found = false;
1.193 + for (int i = 1; i <= length; ++i) {
1.194 + if (indexes[i] == jx) {
1.195 + found = true;
1.196 + values[i] = value;
1.197 + break;
1.198 + }
1.199 + }
1.200 + if (!found) {
1.201 + ++length;
1.202 + indexes[length] = jx;
1.203 + values[length] = value;
1.204 + }
1.205 +
1.206 + glp_set_mat_row(lp, ix, length, &indexes.front(), &values.front());
1.207 +
1.208 + } else {
1.209 +
1.210 + int length = glp_get_mat_col(lp, jx, 0, 0);
1.211 +
1.212 + std::vector<int> indexes(length + 2);
1.213 + std::vector<Value> values(length + 2);
1.214 +
1.215 + glp_get_mat_col(lp, jx, &indexes.front(), &values.front());
1.216 +
1.217 + //The following code does not suppose that the elements of the
1.218 + //array indexes are sorted
1.219 + bool found = false;
1.220 + for (int i = 1; i <= length; ++i) {
1.221 + if (indexes[i] == ix) {
1.222 + found = true;
1.223 + values[i] = value;
1.224 + break;
1.225 + }
1.226 + }
1.227 + if (!found) {
1.228 + ++length;
1.229 + indexes[length] = ix;
1.230 + values[length] = value;
1.231 + }
1.232 +
1.233 + glp_set_mat_col(lp, jx, length, &indexes.front(), &values.front());
1.234 + }
1.235 +
1.236 + }
1.237 +
1.238 + GlpkBase::Value GlpkBase::_getCoeff(int ix, int jx) const {
1.239 +
1.240 + int length = glp_get_mat_row(lp, ix, 0, 0);
1.241 +
1.242 + std::vector<int> indexes(length + 1);
1.243 + std::vector<Value> values(length + 1);
1.244 +
1.245 + glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
1.246 +
1.247 + for (int i = 1; i <= length; ++i) {
1.248 + if (indexes[i] == jx) {
1.249 + return values[i];
1.250 + }
1.251 + }
1.252 +
1.253 + return 0;
1.254 + }
1.255 +
1.256 + void GlpkBase::_setColLowerBound(int i, Value lo) {
1.257 + LEMON_ASSERT(lo != INF, "Invalid bound");
1.258 +
1.259 + int b = glp_get_col_type(lp, i);
1.260 + double up = glp_get_col_ub(lp, i);
1.261 + if (lo == -INF) {
1.262 + switch (b) {
1.263 + case GLP_FR:
1.264 + case GLP_LO:
1.265 + glp_set_col_bnds(lp, i, GLP_FR, lo, up);
1.266 + break;
1.267 + case GLP_UP:
1.268 + break;
1.269 + case GLP_DB:
1.270 + case GLP_FX:
1.271 + glp_set_col_bnds(lp, i, GLP_UP, lo, up);
1.272 + break;
1.273 + default:
1.274 + break;
1.275 + }
1.276 + } else {
1.277 + switch (b) {
1.278 + case GLP_FR:
1.279 + case GLP_LO:
1.280 + glp_set_col_bnds(lp, i, GLP_LO, lo, up);
1.281 + break;
1.282 + case GLP_UP:
1.283 + case GLP_DB:
1.284 + case GLP_FX:
1.285 + if (lo == up)
1.286 + glp_set_col_bnds(lp, i, GLP_FX, lo, up);
1.287 + else
1.288 + glp_set_col_bnds(lp, i, GLP_DB, lo, up);
1.289 + break;
1.290 + default:
1.291 + break;
1.292 + }
1.293 + }
1.294 + }
1.295 +
1.296 + GlpkBase::Value GlpkBase::_getColLowerBound(int i) const {
1.297 + int b = glp_get_col_type(lp, i);
1.298 + switch (b) {
1.299 + case GLP_LO:
1.300 + case GLP_DB:
1.301 + case GLP_FX:
1.302 + return glp_get_col_lb(lp, i);
1.303 + default:
1.304 + return -INF;
1.305 + }
1.306 + }
1.307 +
1.308 + void GlpkBase::_setColUpperBound(int i, Value up) {
1.309 + LEMON_ASSERT(up != -INF, "Invalid bound");
1.310 +
1.311 + int b = glp_get_col_type(lp, i);
1.312 + double lo = glp_get_col_lb(lp, i);
1.313 + if (up == INF) {
1.314 + switch (b) {
1.315 + case GLP_FR:
1.316 + case GLP_LO:
1.317 + break;
1.318 + case GLP_UP:
1.319 + glp_set_col_bnds(lp, i, GLP_FR, lo, up);
1.320 + break;
1.321 + case GLP_DB:
1.322 + case GLP_FX:
1.323 + glp_set_col_bnds(lp, i, GLP_LO, lo, up);
1.324 + break;
1.325 + default:
1.326 + break;
1.327 + }
1.328 + } else {
1.329 + switch (b) {
1.330 + case GLP_FR:
1.331 + glp_set_col_bnds(lp, i, GLP_UP, lo, up);
1.332 + break;
1.333 + case GLP_UP:
1.334 + glp_set_col_bnds(lp, i, GLP_UP, lo, up);
1.335 + break;
1.336 + case GLP_LO:
1.337 + case GLP_DB:
1.338 + case GLP_FX:
1.339 + if (lo == up)
1.340 + glp_set_col_bnds(lp, i, GLP_FX, lo, up);
1.341 + else
1.342 + glp_set_col_bnds(lp, i, GLP_DB, lo, up);
1.343 + break;
1.344 + default:
1.345 + break;
1.346 + }
1.347 + }
1.348 +
1.349 + }
1.350 +
1.351 + GlpkBase::Value GlpkBase::_getColUpperBound(int i) const {
1.352 + int b = glp_get_col_type(lp, i);
1.353 + switch (b) {
1.354 + case GLP_UP:
1.355 + case GLP_DB:
1.356 + case GLP_FX:
1.357 + return glp_get_col_ub(lp, i);
1.358 + default:
1.359 + return INF;
1.360 + }
1.361 + }
1.362 +
1.363 + void GlpkBase::_setRowLowerBound(int i, Value lo) {
1.364 + LEMON_ASSERT(lo != INF, "Invalid bound");
1.365 +
1.366 + int b = glp_get_row_type(lp, i);
1.367 + double up = glp_get_row_ub(lp, i);
1.368 + if (lo == -INF) {
1.369 + switch (b) {
1.370 + case GLP_FR:
1.371 + case GLP_LO:
1.372 + glp_set_row_bnds(lp, i, GLP_FR, lo, up);
1.373 + break;
1.374 + case GLP_UP:
1.375 + break;
1.376 + case GLP_DB:
1.377 + case GLP_FX:
1.378 + glp_set_row_bnds(lp, i, GLP_UP, lo, up);
1.379 + break;
1.380 + default:
1.381 + break;
1.382 + }
1.383 + } else {
1.384 + switch (b) {
1.385 + case GLP_FR:
1.386 + case GLP_LO:
1.387 + glp_set_row_bnds(lp, i, GLP_LO, lo, up);
1.388 + break;
1.389 + case GLP_UP:
1.390 + case GLP_DB:
1.391 + case GLP_FX:
1.392 + if (lo == up)
1.393 + glp_set_row_bnds(lp, i, GLP_FX, lo, up);
1.394 + else
1.395 + glp_set_row_bnds(lp, i, GLP_DB, lo, up);
1.396 + break;
1.397 + default:
1.398 + break;
1.399 + }
1.400 + }
1.401 +
1.402 + }
1.403 +
1.404 + GlpkBase::Value GlpkBase::_getRowLowerBound(int i) const {
1.405 + int b = glp_get_row_type(lp, i);
1.406 + switch (b) {
1.407 + case GLP_LO:
1.408 + case GLP_DB:
1.409 + case GLP_FX:
1.410 + return glp_get_row_lb(lp, i);
1.411 + default:
1.412 + return -INF;
1.413 + }
1.414 + }
1.415 +
1.416 + void GlpkBase::_setRowUpperBound(int i, Value up) {
1.417 + LEMON_ASSERT(up != -INF, "Invalid bound");
1.418 +
1.419 + int b = glp_get_row_type(lp, i);
1.420 + double lo = glp_get_row_lb(lp, i);
1.421 + if (up == INF) {
1.422 + switch (b) {
1.423 + case GLP_FR:
1.424 + case GLP_LO:
1.425 + break;
1.426 + case GLP_UP:
1.427 + glp_set_row_bnds(lp, i, GLP_FR, lo, up);
1.428 + break;
1.429 + case GLP_DB:
1.430 + case GLP_FX:
1.431 + glp_set_row_bnds(lp, i, GLP_LO, lo, up);
1.432 + break;
1.433 + default:
1.434 + break;
1.435 + }
1.436 + } else {
1.437 + switch (b) {
1.438 + case GLP_FR:
1.439 + glp_set_row_bnds(lp, i, GLP_UP, lo, up);
1.440 + break;
1.441 + case GLP_UP:
1.442 + glp_set_row_bnds(lp, i, GLP_UP, lo, up);
1.443 + break;
1.444 + case GLP_LO:
1.445 + case GLP_DB:
1.446 + case GLP_FX:
1.447 + if (lo == up)
1.448 + glp_set_row_bnds(lp, i, GLP_FX, lo, up);
1.449 + else
1.450 + glp_set_row_bnds(lp, i, GLP_DB, lo, up);
1.451 + break;
1.452 + default:
1.453 + break;
1.454 + }
1.455 + }
1.456 + }
1.457 +
1.458 + GlpkBase::Value GlpkBase::_getRowUpperBound(int i) const {
1.459 + int b = glp_get_row_type(lp, i);
1.460 + switch (b) {
1.461 + case GLP_UP:
1.462 + case GLP_DB:
1.463 + case GLP_FX:
1.464 + return glp_get_row_ub(lp, i);
1.465 + default:
1.466 + return INF;
1.467 + }
1.468 + }
1.469 +
1.470 + void GlpkBase::_setObjCoeffs(ExprIterator b, ExprIterator e) {
1.471 + for (int i = 1; i <= glp_get_num_cols(lp); ++i) {
1.472 + glp_set_obj_coef(lp, i, 0.0);
1.473 + }
1.474 + for (ExprIterator it = b; it != e; ++it) {
1.475 + glp_set_obj_coef(lp, it->first, it->second);
1.476 + }
1.477 + }
1.478 +
1.479 + void GlpkBase::_getObjCoeffs(InsertIterator b) const {
1.480 + for (int i = 1; i <= glp_get_num_cols(lp); ++i) {
1.481 + Value val = glp_get_obj_coef(lp, i);
1.482 + if (val != 0.0) {
1.483 + *b = std::make_pair(i, val);
1.484 + ++b;
1.485 + }
1.486 + }
1.487 + }
1.488 +
1.489 + void GlpkBase::_setObjCoeff(int i, Value obj_coef) {
1.490 + //i = 0 means the constant term (shift)
1.491 + glp_set_obj_coef(lp, i, obj_coef);
1.492 + }
1.493 +
1.494 + GlpkBase::Value GlpkBase::_getObjCoeff(int i) const {
1.495 + //i = 0 means the constant term (shift)
1.496 + return glp_get_obj_coef(lp, i);
1.497 + }
1.498 +
1.499 + void GlpkBase::_setSense(GlpkBase::Sense sense) {
1.500 + switch (sense) {
1.501 + case MIN:
1.502 + glp_set_obj_dir(lp, GLP_MIN);
1.503 + break;
1.504 + case MAX:
1.505 + glp_set_obj_dir(lp, GLP_MAX);
1.506 + break;
1.507 + }
1.508 + }
1.509 +
1.510 + GlpkBase::Sense GlpkBase::_getSense() const {
1.511 + switch(glp_get_obj_dir(lp)) {
1.512 + case GLP_MIN:
1.513 + return MIN;
1.514 + case GLP_MAX:
1.515 + return MAX;
1.516 + default:
1.517 + LEMON_ASSERT(false, "Wrong sense");
1.518 + return GlpkBase::Sense();
1.519 + }
1.520 + }
1.521 +
1.522 + void GlpkBase::_clear() {
1.523 + glp_erase_prob(lp);
1.524 + rows.clear();
1.525 + cols.clear();
1.526 + }
1.527 +
1.528 + // LpGlpk members
1.529 +
1.530 + LpGlpk::LpGlpk()
1.531 + : LpBase(), GlpkBase(), LpSolver() {
1.532 + messageLevel(MESSAGE_NO_OUTPUT);
1.533 + }
1.534 +
1.535 + LpGlpk::LpGlpk(const LpGlpk& other)
1.536 + : LpBase(other), GlpkBase(other), LpSolver(other) {
1.537 + messageLevel(MESSAGE_NO_OUTPUT);
1.538 + }
1.539 +
1.540 + LpGlpk* LpGlpk::_newSolver() const { return new LpGlpk; }
1.541 + LpGlpk* LpGlpk::_cloneSolver() const { return new LpGlpk(*this); }
1.542 +
1.543 + const char* LpGlpk::_solverName() const { return "LpGlpk"; }
1.544 +
1.545 + void LpGlpk::_clear_temporals() {
1.546 + _primal_ray.clear();
1.547 + _dual_ray.clear();
1.548 + }
1.549 +
1.550 + LpGlpk::SolveExitStatus LpGlpk::_solve() {
1.551 + return solvePrimal();
1.552 + }
1.553 +
1.554 + LpGlpk::SolveExitStatus LpGlpk::solvePrimal() {
1.555 + _clear_temporals();
1.556 +
1.557 + glp_smcp smcp;
1.558 + glp_init_smcp(&smcp);
1.559 +
1.560 + switch (_message_level) {
1.561 + case MESSAGE_NO_OUTPUT:
1.562 + smcp.msg_lev = GLP_MSG_OFF;
1.563 + break;
1.564 + case MESSAGE_ERROR_MESSAGE:
1.565 + smcp.msg_lev = GLP_MSG_ERR;
1.566 + break;
1.567 + case MESSAGE_NORMAL_OUTPUT:
1.568 + smcp.msg_lev = GLP_MSG_ON;
1.569 + break;
1.570 + case MESSAGE_FULL_OUTPUT:
1.571 + smcp.msg_lev = GLP_MSG_ALL;
1.572 + break;
1.573 + }
1.574 +
1.575 + if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
1.576 + return SOLVED;
1.577 + }
1.578 +
1.579 + LpGlpk::SolveExitStatus LpGlpk::solveDual() {
1.580 + _clear_temporals();
1.581 +
1.582 + glp_smcp smcp;
1.583 + glp_init_smcp(&smcp);
1.584 +
1.585 + switch (_message_level) {
1.586 + case MESSAGE_NO_OUTPUT:
1.587 + smcp.msg_lev = GLP_MSG_OFF;
1.588 + break;
1.589 + case MESSAGE_ERROR_MESSAGE:
1.590 + smcp.msg_lev = GLP_MSG_ERR;
1.591 + break;
1.592 + case MESSAGE_NORMAL_OUTPUT:
1.593 + smcp.msg_lev = GLP_MSG_ON;
1.594 + break;
1.595 + case MESSAGE_FULL_OUTPUT:
1.596 + smcp.msg_lev = GLP_MSG_ALL;
1.597 + break;
1.598 + }
1.599 + smcp.meth = GLP_DUAL;
1.600 +
1.601 + if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
1.602 + return SOLVED;
1.603 + }
1.604 +
1.605 + LpGlpk::Value LpGlpk::_getPrimal(int i) const {
1.606 + return glp_get_col_prim(lp, i);
1.607 + }
1.608 +
1.609 + LpGlpk::Value LpGlpk::_getDual(int i) const {
1.610 + return glp_get_row_dual(lp, i);
1.611 + }
1.612 +
1.613 + LpGlpk::Value LpGlpk::_getPrimalValue() const {
1.614 + return glp_get_obj_val(lp);
1.615 + }
1.616 +
1.617 + LpGlpk::VarStatus LpGlpk::_getColStatus(int i) const {
1.618 + switch (glp_get_col_stat(lp, i)) {
1.619 + case GLP_BS:
1.620 + return BASIC;
1.621 + case GLP_UP:
1.622 + return UPPER;
1.623 + case GLP_LO:
1.624 + return LOWER;
1.625 + case GLP_NF:
1.626 + return FREE;
1.627 + case GLP_NS:
1.628 + return FIXED;
1.629 + default:
1.630 + LEMON_ASSERT(false, "Wrong column status");
1.631 + return LpGlpk::VarStatus();
1.632 + }
1.633 + }
1.634 +
1.635 + LpGlpk::VarStatus LpGlpk::_getRowStatus(int i) const {
1.636 + switch (glp_get_row_stat(lp, i)) {
1.637 + case GLP_BS:
1.638 + return BASIC;
1.639 + case GLP_UP:
1.640 + return UPPER;
1.641 + case GLP_LO:
1.642 + return LOWER;
1.643 + case GLP_NF:
1.644 + return FREE;
1.645 + case GLP_NS:
1.646 + return FIXED;
1.647 + default:
1.648 + LEMON_ASSERT(false, "Wrong row status");
1.649 + return LpGlpk::VarStatus();
1.650 + }
1.651 + }
1.652 +
1.653 + LpGlpk::Value LpGlpk::_getPrimalRay(int i) const {
1.654 + if (_primal_ray.empty()) {
1.655 + int row_num = glp_get_num_rows(lp);
1.656 + int col_num = glp_get_num_cols(lp);
1.657 +
1.658 + _primal_ray.resize(col_num + 1, 0.0);
1.659 +
1.660 + int index = glp_get_unbnd_ray(lp);
1.661 + if (index != 0) {
1.662 + // The primal ray is found in primal simplex second phase
1.663 + LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :
1.664 + glp_get_col_stat(lp, index - row_num)) != GLP_BS,
1.665 + "Wrong primal ray");
1.666 +
1.667 + bool negate = glp_get_obj_dir(lp) == GLP_MAX;
1.668 +
1.669 + if (index > row_num) {
1.670 + _primal_ray[index - row_num] = 1.0;
1.671 + if (glp_get_col_dual(lp, index - row_num) > 0) {
1.672 + negate = !negate;
1.673 + }
1.674 + } else {
1.675 + if (glp_get_row_dual(lp, index) > 0) {
1.676 + negate = !negate;
1.677 + }
1.678 + }
1.679 +
1.680 + std::vector<int> ray_indexes(row_num + 1);
1.681 + std::vector<Value> ray_values(row_num + 1);
1.682 + int ray_length = glp_eval_tab_col(lp, index, &ray_indexes.front(),
1.683 + &ray_values.front());
1.684 +
1.685 + for (int i = 1; i <= ray_length; ++i) {
1.686 + if (ray_indexes[i] > row_num) {
1.687 + _primal_ray[ray_indexes[i] - row_num] = ray_values[i];
1.688 + }
1.689 + }
1.690 +
1.691 + if (negate) {
1.692 + for (int i = 1; i <= col_num; ++i) {
1.693 + _primal_ray[i] = - _primal_ray[i];
1.694 + }
1.695 + }
1.696 + } else {
1.697 + for (int i = 1; i <= col_num; ++i) {
1.698 + _primal_ray[i] = glp_get_col_prim(lp, i);
1.699 + }
1.700 + }
1.701 + }
1.702 + return _primal_ray[i];
1.703 + }
1.704 +
1.705 + LpGlpk::Value LpGlpk::_getDualRay(int i) const {
1.706 + if (_dual_ray.empty()) {
1.707 + int row_num = glp_get_num_rows(lp);
1.708 +
1.709 + _dual_ray.resize(row_num + 1, 0.0);
1.710 +
1.711 + int index = glp_get_unbnd_ray(lp);
1.712 + if (index != 0) {
1.713 + // The dual ray is found in dual simplex second phase
1.714 + LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :
1.715 + glp_get_col_stat(lp, index - row_num)) == GLP_BS,
1.716 +
1.717 + "Wrong dual ray");
1.718 +
1.719 + int idx;
1.720 + bool negate = false;
1.721 +
1.722 + if (index > row_num) {
1.723 + idx = glp_get_col_bind(lp, index - row_num);
1.724 + if (glp_get_col_prim(lp, index - row_num) >
1.725 + glp_get_col_ub(lp, index - row_num)) {
1.726 + negate = true;
1.727 + }
1.728 + } else {
1.729 + idx = glp_get_row_bind(lp, index);
1.730 + if (glp_get_row_prim(lp, index) > glp_get_row_ub(lp, index)) {
1.731 + negate = true;
1.732 + }
1.733 + }
1.734 +
1.735 + _dual_ray[idx] = negate ? - 1.0 : 1.0;
1.736 +
1.737 + glp_btran(lp, &_dual_ray.front());
1.738 + } else {
1.739 + double eps = 1e-7;
1.740 + // The dual ray is found in primal simplex first phase
1.741 + // We assume that the glpk minimizes the slack to get feasible solution
1.742 + for (int i = 1; i <= row_num; ++i) {
1.743 + int index = glp_get_bhead(lp, i);
1.744 + if (index <= row_num) {
1.745 + double res = glp_get_row_prim(lp, index);
1.746 + if (res > glp_get_row_ub(lp, index) + eps) {
1.747 + _dual_ray[i] = -1;
1.748 + } else if (res < glp_get_row_lb(lp, index) - eps) {
1.749 + _dual_ray[i] = 1;
1.750 + } else {
1.751 + _dual_ray[i] = 0;
1.752 + }
1.753 + _dual_ray[i] *= glp_get_rii(lp, index);
1.754 + } else {
1.755 + double res = glp_get_col_prim(lp, index - row_num);
1.756 + if (res > glp_get_col_ub(lp, index - row_num) + eps) {
1.757 + _dual_ray[i] = -1;
1.758 + } else if (res < glp_get_col_lb(lp, index - row_num) - eps) {
1.759 + _dual_ray[i] = 1;
1.760 + } else {
1.761 + _dual_ray[i] = 0;
1.762 + }
1.763 + _dual_ray[i] /= glp_get_sjj(lp, index - row_num);
1.764 + }
1.765 + }
1.766 +
1.767 + glp_btran(lp, &_dual_ray.front());
1.768 +
1.769 + for (int i = 1; i <= row_num; ++i) {
1.770 + _dual_ray[i] /= glp_get_rii(lp, i);
1.771 + }
1.772 + }
1.773 + }
1.774 + return _dual_ray[i];
1.775 + }
1.776 +
1.777 + LpGlpk::ProblemType LpGlpk::_getPrimalType() const {
1.778 + if (glp_get_status(lp) == GLP_OPT)
1.779 + return OPTIMAL;
1.780 + switch (glp_get_prim_stat(lp)) {
1.781 + case GLP_UNDEF:
1.782 + return UNDEFINED;
1.783 + case GLP_FEAS:
1.784 + case GLP_INFEAS:
1.785 + if (glp_get_dual_stat(lp) == GLP_NOFEAS) {
1.786 + return UNBOUNDED;
1.787 + } else {
1.788 + return UNDEFINED;
1.789 + }
1.790 + case GLP_NOFEAS:
1.791 + return INFEASIBLE;
1.792 + default:
1.793 + LEMON_ASSERT(false, "Wrong primal type");
1.794 + return LpGlpk::ProblemType();
1.795 + }
1.796 + }
1.797 +
1.798 + LpGlpk::ProblemType LpGlpk::_getDualType() const {
1.799 + if (glp_get_status(lp) == GLP_OPT)
1.800 + return OPTIMAL;
1.801 + switch (glp_get_dual_stat(lp)) {
1.802 + case GLP_UNDEF:
1.803 + return UNDEFINED;
1.804 + case GLP_FEAS:
1.805 + case GLP_INFEAS:
1.806 + if (glp_get_prim_stat(lp) == GLP_NOFEAS) {
1.807 + return UNBOUNDED;
1.808 + } else {
1.809 + return UNDEFINED;
1.810 + }
1.811 + case GLP_NOFEAS:
1.812 + return INFEASIBLE;
1.813 + default:
1.814 + LEMON_ASSERT(false, "Wrong primal type");
1.815 + return LpGlpk::ProblemType();
1.816 + }
1.817 + }
1.818 +
1.819 + void LpGlpk::presolver(bool b) {
1.820 + lpx_set_int_parm(lp, LPX_K_PRESOL, b ? 1 : 0);
1.821 + }
1.822 +
1.823 + void LpGlpk::messageLevel(MessageLevel m) {
1.824 + _message_level = m;
1.825 + }
1.826 +
1.827 + // MipGlpk members
1.828 +
1.829 + MipGlpk::MipGlpk()
1.830 + : LpBase(), GlpkBase(), MipSolver() {
1.831 + messageLevel(MESSAGE_NO_OUTPUT);
1.832 + }
1.833 +
1.834 + MipGlpk::MipGlpk(const MipGlpk& other)
1.835 + : LpBase(), GlpkBase(other), MipSolver() {
1.836 + messageLevel(MESSAGE_NO_OUTPUT);
1.837 + }
1.838 +
1.839 + void MipGlpk::_setColType(int i, MipGlpk::ColTypes col_type) {
1.840 + switch (col_type) {
1.841 + case INTEGER:
1.842 + glp_set_col_kind(lp, i, GLP_IV);
1.843 + break;
1.844 + case REAL:
1.845 + glp_set_col_kind(lp, i, GLP_CV);
1.846 + break;
1.847 + }
1.848 + }
1.849 +
1.850 + MipGlpk::ColTypes MipGlpk::_getColType(int i) const {
1.851 + switch (glp_get_col_kind(lp, i)) {
1.852 + case GLP_IV:
1.853 + case GLP_BV:
1.854 + return INTEGER;
1.855 + default:
1.856 + return REAL;
1.857 + }
1.858 +
1.859 + }
1.860 +
1.861 + MipGlpk::SolveExitStatus MipGlpk::_solve() {
1.862 + glp_smcp smcp;
1.863 + glp_init_smcp(&smcp);
1.864 +
1.865 + switch (_message_level) {
1.866 + case MESSAGE_NO_OUTPUT:
1.867 + smcp.msg_lev = GLP_MSG_OFF;
1.868 + break;
1.869 + case MESSAGE_ERROR_MESSAGE:
1.870 + smcp.msg_lev = GLP_MSG_ERR;
1.871 + break;
1.872 + case MESSAGE_NORMAL_OUTPUT:
1.873 + smcp.msg_lev = GLP_MSG_ON;
1.874 + break;
1.875 + case MESSAGE_FULL_OUTPUT:
1.876 + smcp.msg_lev = GLP_MSG_ALL;
1.877 + break;
1.878 + }
1.879 + smcp.meth = GLP_DUAL;
1.880 +
1.881 + if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
1.882 + if (glp_get_status(lp) != GLP_OPT) return SOLVED;
1.883 +
1.884 + glp_iocp iocp;
1.885 + glp_init_iocp(&iocp);
1.886 +
1.887 + switch (_message_level) {
1.888 + case MESSAGE_NO_OUTPUT:
1.889 + iocp.msg_lev = GLP_MSG_OFF;
1.890 + break;
1.891 + case MESSAGE_ERROR_MESSAGE:
1.892 + iocp.msg_lev = GLP_MSG_ERR;
1.893 + break;
1.894 + case MESSAGE_NORMAL_OUTPUT:
1.895 + iocp.msg_lev = GLP_MSG_ON;
1.896 + break;
1.897 + case MESSAGE_FULL_OUTPUT:
1.898 + iocp.msg_lev = GLP_MSG_ALL;
1.899 + break;
1.900 + }
1.901 +
1.902 + if (glp_intopt(lp, &iocp) != 0) return UNSOLVED;
1.903 + return SOLVED;
1.904 + }
1.905 +
1.906 +
1.907 + MipGlpk::ProblemType MipGlpk::_getType() const {
1.908 + switch (glp_get_status(lp)) {
1.909 + case GLP_OPT:
1.910 + switch (glp_mip_status(lp)) {
1.911 + case GLP_UNDEF:
1.912 + return UNDEFINED;
1.913 + case GLP_NOFEAS:
1.914 + return INFEASIBLE;
1.915 + case GLP_FEAS:
1.916 + return FEASIBLE;
1.917 + case GLP_OPT:
1.918 + return OPTIMAL;
1.919 + default:
1.920 + LEMON_ASSERT(false, "Wrong problem type.");
1.921 + return MipGlpk::ProblemType();
1.922 + }
1.923 + case GLP_NOFEAS:
1.924 + return INFEASIBLE;
1.925 + case GLP_INFEAS:
1.926 + case GLP_FEAS:
1.927 + if (glp_get_dual_stat(lp) == GLP_NOFEAS) {
1.928 + return UNBOUNDED;
1.929 + } else {
1.930 + return UNDEFINED;
1.931 + }
1.932 + default:
1.933 + LEMON_ASSERT(false, "Wrong problem type.");
1.934 + return MipGlpk::ProblemType();
1.935 + }
1.936 + }
1.937 +
1.938 + MipGlpk::Value MipGlpk::_getSol(int i) const {
1.939 + return glp_mip_col_val(lp, i);
1.940 + }
1.941 +
1.942 + MipGlpk::Value MipGlpk::_getSolValue() const {
1.943 + return glp_mip_obj_val(lp);
1.944 + }
1.945 +
1.946 + MipGlpk* MipGlpk::_newSolver() const { return new MipGlpk; }
1.947 + MipGlpk* MipGlpk::_cloneSolver() const {return new MipGlpk(*this); }
1.948 +
1.949 + const char* MipGlpk::_solverName() const { return "MipGlpk"; }
1.950 +
1.951 + void MipGlpk::messageLevel(MessageLevel m) {
1.952 + _message_level = m;
1.953 + }
1.954 +
1.955 +} //END OF NAMESPACE LEMON