1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
3 * This file is a part of LEMON, a generic C++ optimization library.
5 * Copyright (C) 2003-2008
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
9 * Permission to use, modify and distribute this software is granted
10 * provided that this copyright notice appears in all copies. For
11 * precise terms see the accompanying LICENSE file.
13 * This software is provided "AS IS" with no warranty of any kind,
14 * express or implied, and with no claim as to its suitability for any
20 ///\brief Implementation of the LEMON GLPK LP and MIP solver interface.
22 #include <lemon/glpk.h>
25 #include <lemon/assert.h>
31 GlpkBase::GlpkBase() : LpBase() {
32 lp = glp_create_prob();
36 GlpkBase::GlpkBase(const GlpkBase &other) : LpBase() {
37 lp = glp_create_prob();
38 glp_copy_prob(lp, other.lp, GLP_ON);
44 GlpkBase::~GlpkBase() {
48 int GlpkBase::_addCol() {
49 int i = glp_add_cols(lp, 1);
50 glp_set_col_bnds(lp, i, GLP_FR, 0.0, 0.0);
54 int GlpkBase::_addRow() {
55 int i = glp_add_rows(lp, 1);
56 glp_set_row_bnds(lp, i, GLP_FR, 0.0, 0.0);
60 void GlpkBase::_eraseCol(int i) {
63 glp_del_cols(lp, 1, ca);
66 void GlpkBase::_eraseRow(int i) {
69 glp_del_rows(lp, 1, ra);
72 void GlpkBase::_eraseColId(int i) {
77 void GlpkBase::_eraseRowId(int i) {
82 void GlpkBase::_getColName(int c, std::string& name) const {
83 const char *str = glp_get_col_name(lp, c);
88 void GlpkBase::_setColName(int c, const std::string & name) {
89 glp_set_col_name(lp, c, const_cast<char*>(name.c_str()));
93 int GlpkBase::_colByName(const std::string& name) const {
94 int k = glp_find_col(lp, const_cast<char*>(name.c_str()));
95 return k > 0 ? k : -1;
98 void GlpkBase::_getRowName(int r, std::string& name) const {
99 const char *str = glp_get_row_name(lp, r);
104 void GlpkBase::_setRowName(int r, const std::string & name) {
105 glp_set_row_name(lp, r, const_cast<char*>(name.c_str()));
109 int GlpkBase::_rowByName(const std::string& name) const {
110 int k = glp_find_row(lp, const_cast<char*>(name.c_str()));
111 return k > 0 ? k : -1;
114 void GlpkBase::_setRowCoeffs(int i, ExprIterator b, ExprIterator e) {
115 std::vector<int> indexes;
116 std::vector<Value> values;
118 indexes.push_back(0);
121 for(ExprIterator it = b; it != e; ++it) {
122 indexes.push_back(it->first);
123 values.push_back(it->second);
126 glp_set_mat_row(lp, i, values.size() - 1,
127 &indexes.front(), &values.front());
130 void GlpkBase::_getRowCoeffs(int ix, InsertIterator b) const {
131 int length = glp_get_mat_row(lp, ix, 0, 0);
133 std::vector<int> indexes(length + 1);
134 std::vector<Value> values(length + 1);
136 glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
138 for (int i = 1; i <= length; ++i) {
139 *b = std::make_pair(indexes[i], values[i]);
144 void GlpkBase::_setColCoeffs(int ix, ExprIterator b,
147 std::vector<int> indexes;
148 std::vector<Value> values;
150 indexes.push_back(0);
153 for(ExprIterator it = b; it != e; ++it) {
154 indexes.push_back(it->first);
155 values.push_back(it->second);
158 glp_set_mat_col(lp, ix, values.size() - 1,
159 &indexes.front(), &values.front());
162 void GlpkBase::_getColCoeffs(int ix, InsertIterator b) const {
163 int length = glp_get_mat_col(lp, ix, 0, 0);
165 std::vector<int> indexes(length + 1);
166 std::vector<Value> values(length + 1);
168 glp_get_mat_col(lp, ix, &indexes.front(), &values.front());
170 for (int i = 1; i <= length; ++i) {
171 *b = std::make_pair(indexes[i], values[i]);
176 void GlpkBase::_setCoeff(int ix, int jx, Value value) {
178 if (glp_get_num_cols(lp) < glp_get_num_rows(lp)) {
180 int length = glp_get_mat_row(lp, ix, 0, 0);
182 std::vector<int> indexes(length + 2);
183 std::vector<Value> values(length + 2);
185 glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
187 //The following code does not suppose that the elements of the
188 //array indexes are sorted
190 for (int i = 1; i <= length; ++i) {
191 if (indexes[i] == jx) {
199 indexes[length] = jx;
200 values[length] = value;
203 glp_set_mat_row(lp, ix, length, &indexes.front(), &values.front());
207 int length = glp_get_mat_col(lp, jx, 0, 0);
209 std::vector<int> indexes(length + 2);
210 std::vector<Value> values(length + 2);
212 glp_get_mat_col(lp, jx, &indexes.front(), &values.front());
214 //The following code does not suppose that the elements of the
215 //array indexes are sorted
217 for (int i = 1; i <= length; ++i) {
218 if (indexes[i] == ix) {
226 indexes[length] = ix;
227 values[length] = value;
230 glp_set_mat_col(lp, jx, length, &indexes.front(), &values.front());
235 GlpkBase::Value GlpkBase::_getCoeff(int ix, int jx) const {
237 int length = glp_get_mat_row(lp, ix, 0, 0);
239 std::vector<int> indexes(length + 1);
240 std::vector<Value> values(length + 1);
242 glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
244 for (int i = 1; i <= length; ++i) {
245 if (indexes[i] == jx) {
253 void GlpkBase::_setColLowerBound(int i, Value lo) {
254 LEMON_ASSERT(lo != INF, "Invalid bound");
256 int b = glp_get_col_type(lp, i);
257 double up = glp_get_col_ub(lp, i);
262 glp_set_col_bnds(lp, i, GLP_FR, lo, up);
268 glp_set_col_bnds(lp, i, GLP_UP, lo, up);
277 glp_set_col_bnds(lp, i, GLP_LO, lo, up);
283 glp_set_col_bnds(lp, i, GLP_FX, lo, up);
285 glp_set_col_bnds(lp, i, GLP_DB, lo, up);
293 GlpkBase::Value GlpkBase::_getColLowerBound(int i) const {
294 int b = glp_get_col_type(lp, i);
299 return glp_get_col_lb(lp, i);
305 void GlpkBase::_setColUpperBound(int i, Value up) {
306 LEMON_ASSERT(up != -INF, "Invalid bound");
308 int b = glp_get_col_type(lp, i);
309 double lo = glp_get_col_lb(lp, i);
316 glp_set_col_bnds(lp, i, GLP_FR, lo, up);
320 glp_set_col_bnds(lp, i, GLP_LO, lo, up);
328 glp_set_col_bnds(lp, i, GLP_UP, lo, up);
331 glp_set_col_bnds(lp, i, GLP_UP, lo, up);
337 glp_set_col_bnds(lp, i, GLP_FX, lo, up);
339 glp_set_col_bnds(lp, i, GLP_DB, lo, up);
348 GlpkBase::Value GlpkBase::_getColUpperBound(int i) const {
349 int b = glp_get_col_type(lp, i);
354 return glp_get_col_ub(lp, i);
360 void GlpkBase::_setRowLowerBound(int i, Value lo) {
361 LEMON_ASSERT(lo != INF, "Invalid bound");
363 int b = glp_get_row_type(lp, i);
364 double up = glp_get_row_ub(lp, i);
369 glp_set_row_bnds(lp, i, GLP_FR, lo, up);
375 glp_set_row_bnds(lp, i, GLP_UP, lo, up);
384 glp_set_row_bnds(lp, i, GLP_LO, lo, up);
390 glp_set_row_bnds(lp, i, GLP_FX, lo, up);
392 glp_set_row_bnds(lp, i, GLP_DB, lo, up);
401 GlpkBase::Value GlpkBase::_getRowLowerBound(int i) const {
402 int b = glp_get_row_type(lp, i);
407 return glp_get_row_lb(lp, i);
413 void GlpkBase::_setRowUpperBound(int i, Value up) {
414 LEMON_ASSERT(up != -INF, "Invalid bound");
416 int b = glp_get_row_type(lp, i);
417 double lo = glp_get_row_lb(lp, i);
424 glp_set_row_bnds(lp, i, GLP_FR, lo, up);
428 glp_set_row_bnds(lp, i, GLP_LO, lo, up);
436 glp_set_row_bnds(lp, i, GLP_UP, lo, up);
439 glp_set_row_bnds(lp, i, GLP_UP, lo, up);
445 glp_set_row_bnds(lp, i, GLP_FX, lo, up);
447 glp_set_row_bnds(lp, i, GLP_DB, lo, up);
455 GlpkBase::Value GlpkBase::_getRowUpperBound(int i) const {
456 int b = glp_get_row_type(lp, i);
461 return glp_get_row_ub(lp, i);
467 void GlpkBase::_setObjCoeffs(ExprIterator b, ExprIterator e) {
468 for (int i = 1; i <= glp_get_num_cols(lp); ++i) {
469 glp_set_obj_coef(lp, i, 0.0);
471 for (ExprIterator it = b; it != e; ++it) {
472 glp_set_obj_coef(lp, it->first, it->second);
476 void GlpkBase::_getObjCoeffs(InsertIterator b) const {
477 for (int i = 1; i <= glp_get_num_cols(lp); ++i) {
478 Value val = glp_get_obj_coef(lp, i);
480 *b = std::make_pair(i, val);
486 void GlpkBase::_setObjCoeff(int i, Value obj_coef) {
487 //i = 0 means the constant term (shift)
488 glp_set_obj_coef(lp, i, obj_coef);
491 GlpkBase::Value GlpkBase::_getObjCoeff(int i) const {
492 //i = 0 means the constant term (shift)
493 return glp_get_obj_coef(lp, i);
496 void GlpkBase::_setSense(GlpkBase::Sense sense) {
499 glp_set_obj_dir(lp, GLP_MIN);
502 glp_set_obj_dir(lp, GLP_MAX);
507 GlpkBase::Sense GlpkBase::_getSense() const {
508 switch(glp_get_obj_dir(lp)) {
514 LEMON_ASSERT(false, "Wrong sense");
515 return GlpkBase::Sense();
519 void GlpkBase::_clear() {
525 void GlpkBase::freeEnv() {
529 GlpkBase::FreeEnvHelper GlpkBase::freeEnvHelper;
534 : LpBase(), GlpkBase(), LpSolver() {
535 messageLevel(MESSAGE_NO_OUTPUT);
538 GlpkLp::GlpkLp(const GlpkLp& other)
539 : LpBase(other), GlpkBase(other), LpSolver(other) {
540 messageLevel(MESSAGE_NO_OUTPUT);
543 GlpkLp* GlpkLp::newSolver() const { return new GlpkLp; }
544 GlpkLp* GlpkLp::cloneSolver() const { return new GlpkLp(*this); }
546 const char* GlpkLp::_solverName() const { return "GlpkLp"; }
548 void GlpkLp::_clear_temporals() {
553 GlpkLp::SolveExitStatus GlpkLp::_solve() {
554 return solvePrimal();
557 GlpkLp::SolveExitStatus GlpkLp::solvePrimal() {
561 glp_init_smcp(&smcp);
563 switch (_message_level) {
564 case MESSAGE_NO_OUTPUT:
565 smcp.msg_lev = GLP_MSG_OFF;
567 case MESSAGE_ERROR_MESSAGE:
568 smcp.msg_lev = GLP_MSG_ERR;
570 case MESSAGE_NORMAL_OUTPUT:
571 smcp.msg_lev = GLP_MSG_ON;
573 case MESSAGE_FULL_OUTPUT:
574 smcp.msg_lev = GLP_MSG_ALL;
578 if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
582 GlpkLp::SolveExitStatus GlpkLp::solveDual() {
586 glp_init_smcp(&smcp);
588 switch (_message_level) {
589 case MESSAGE_NO_OUTPUT:
590 smcp.msg_lev = GLP_MSG_OFF;
592 case MESSAGE_ERROR_MESSAGE:
593 smcp.msg_lev = GLP_MSG_ERR;
595 case MESSAGE_NORMAL_OUTPUT:
596 smcp.msg_lev = GLP_MSG_ON;
598 case MESSAGE_FULL_OUTPUT:
599 smcp.msg_lev = GLP_MSG_ALL;
602 smcp.meth = GLP_DUAL;
604 if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
608 GlpkLp::Value GlpkLp::_getPrimal(int i) const {
609 return glp_get_col_prim(lp, i);
612 GlpkLp::Value GlpkLp::_getDual(int i) const {
613 return glp_get_row_dual(lp, i);
616 GlpkLp::Value GlpkLp::_getPrimalValue() const {
617 return glp_get_obj_val(lp);
620 GlpkLp::VarStatus GlpkLp::_getColStatus(int i) const {
621 switch (glp_get_col_stat(lp, i)) {
633 LEMON_ASSERT(false, "Wrong column status");
634 return GlpkLp::VarStatus();
638 GlpkLp::VarStatus GlpkLp::_getRowStatus(int i) const {
639 switch (glp_get_row_stat(lp, i)) {
651 LEMON_ASSERT(false, "Wrong row status");
652 return GlpkLp::VarStatus();
656 GlpkLp::Value GlpkLp::_getPrimalRay(int i) const {
657 if (_primal_ray.empty()) {
658 int row_num = glp_get_num_rows(lp);
659 int col_num = glp_get_num_cols(lp);
661 _primal_ray.resize(col_num + 1, 0.0);
663 int index = glp_get_unbnd_ray(lp);
665 // The primal ray is found in primal simplex second phase
666 LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :
667 glp_get_col_stat(lp, index - row_num)) != GLP_BS,
670 bool negate = glp_get_obj_dir(lp) == GLP_MAX;
672 if (index > row_num) {
673 _primal_ray[index - row_num] = 1.0;
674 if (glp_get_col_dual(lp, index - row_num) > 0) {
678 if (glp_get_row_dual(lp, index) > 0) {
683 std::vector<int> ray_indexes(row_num + 1);
684 std::vector<Value> ray_values(row_num + 1);
685 int ray_length = glp_eval_tab_col(lp, index, &ray_indexes.front(),
686 &ray_values.front());
688 for (int i = 1; i <= ray_length; ++i) {
689 if (ray_indexes[i] > row_num) {
690 _primal_ray[ray_indexes[i] - row_num] = ray_values[i];
695 for (int i = 1; i <= col_num; ++i) {
696 _primal_ray[i] = - _primal_ray[i];
700 for (int i = 1; i <= col_num; ++i) {
701 _primal_ray[i] = glp_get_col_prim(lp, i);
705 return _primal_ray[i];
708 GlpkLp::Value GlpkLp::_getDualRay(int i) const {
709 if (_dual_ray.empty()) {
710 int row_num = glp_get_num_rows(lp);
712 _dual_ray.resize(row_num + 1, 0.0);
714 int index = glp_get_unbnd_ray(lp);
716 // The dual ray is found in dual simplex second phase
717 LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :
718 glp_get_col_stat(lp, index - row_num)) == GLP_BS,
725 if (index > row_num) {
726 idx = glp_get_col_bind(lp, index - row_num);
727 if (glp_get_col_prim(lp, index - row_num) >
728 glp_get_col_ub(lp, index - row_num)) {
732 idx = glp_get_row_bind(lp, index);
733 if (glp_get_row_prim(lp, index) > glp_get_row_ub(lp, index)) {
738 _dual_ray[idx] = negate ? - 1.0 : 1.0;
740 glp_btran(lp, &_dual_ray.front());
743 // The dual ray is found in primal simplex first phase
744 // We assume that the glpk minimizes the slack to get feasible solution
745 for (int i = 1; i <= row_num; ++i) {
746 int index = glp_get_bhead(lp, i);
747 if (index <= row_num) {
748 double res = glp_get_row_prim(lp, index);
749 if (res > glp_get_row_ub(lp, index) + eps) {
751 } else if (res < glp_get_row_lb(lp, index) - eps) {
756 _dual_ray[i] *= glp_get_rii(lp, index);
758 double res = glp_get_col_prim(lp, index - row_num);
759 if (res > glp_get_col_ub(lp, index - row_num) + eps) {
761 } else if (res < glp_get_col_lb(lp, index - row_num) - eps) {
766 _dual_ray[i] /= glp_get_sjj(lp, index - row_num);
770 glp_btran(lp, &_dual_ray.front());
772 for (int i = 1; i <= row_num; ++i) {
773 _dual_ray[i] /= glp_get_rii(lp, i);
780 GlpkLp::ProblemType GlpkLp::_getPrimalType() const {
781 if (glp_get_status(lp) == GLP_OPT)
783 switch (glp_get_prim_stat(lp)) {
788 if (glp_get_dual_stat(lp) == GLP_NOFEAS) {
796 LEMON_ASSERT(false, "Wrong primal type");
797 return GlpkLp::ProblemType();
801 GlpkLp::ProblemType GlpkLp::_getDualType() const {
802 if (glp_get_status(lp) == GLP_OPT)
804 switch (glp_get_dual_stat(lp)) {
809 if (glp_get_prim_stat(lp) == GLP_NOFEAS) {
817 LEMON_ASSERT(false, "Wrong primal type");
818 return GlpkLp::ProblemType();
822 void GlpkLp::presolver(bool b) {
823 lpx_set_int_parm(lp, LPX_K_PRESOL, b ? 1 : 0);
826 void GlpkLp::messageLevel(MessageLevel m) {
833 : LpBase(), GlpkBase(), MipSolver() {
834 messageLevel(MESSAGE_NO_OUTPUT);
837 GlpkMip::GlpkMip(const GlpkMip& other)
838 : LpBase(), GlpkBase(other), MipSolver() {
839 messageLevel(MESSAGE_NO_OUTPUT);
842 void GlpkMip::_setColType(int i, GlpkMip::ColTypes col_type) {
845 glp_set_col_kind(lp, i, GLP_IV);
848 glp_set_col_kind(lp, i, GLP_CV);
853 GlpkMip::ColTypes GlpkMip::_getColType(int i) const {
854 switch (glp_get_col_kind(lp, i)) {
864 GlpkMip::SolveExitStatus GlpkMip::_solve() {
866 glp_init_smcp(&smcp);
868 switch (_message_level) {
869 case MESSAGE_NO_OUTPUT:
870 smcp.msg_lev = GLP_MSG_OFF;
872 case MESSAGE_ERROR_MESSAGE:
873 smcp.msg_lev = GLP_MSG_ERR;
875 case MESSAGE_NORMAL_OUTPUT:
876 smcp.msg_lev = GLP_MSG_ON;
878 case MESSAGE_FULL_OUTPUT:
879 smcp.msg_lev = GLP_MSG_ALL;
882 smcp.meth = GLP_DUAL;
884 if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
885 if (glp_get_status(lp) != GLP_OPT) return SOLVED;
888 glp_init_iocp(&iocp);
890 switch (_message_level) {
891 case MESSAGE_NO_OUTPUT:
892 iocp.msg_lev = GLP_MSG_OFF;
894 case MESSAGE_ERROR_MESSAGE:
895 iocp.msg_lev = GLP_MSG_ERR;
897 case MESSAGE_NORMAL_OUTPUT:
898 iocp.msg_lev = GLP_MSG_ON;
900 case MESSAGE_FULL_OUTPUT:
901 iocp.msg_lev = GLP_MSG_ALL;
905 if (glp_intopt(lp, &iocp) != 0) return UNSOLVED;
910 GlpkMip::ProblemType GlpkMip::_getType() const {
911 switch (glp_get_status(lp)) {
913 switch (glp_mip_status(lp)) {
923 LEMON_ASSERT(false, "Wrong problem type.");
924 return GlpkMip::ProblemType();
930 if (glp_get_dual_stat(lp) == GLP_NOFEAS) {
936 LEMON_ASSERT(false, "Wrong problem type.");
937 return GlpkMip::ProblemType();
941 GlpkMip::Value GlpkMip::_getSol(int i) const {
942 return glp_mip_col_val(lp, i);
945 GlpkMip::Value GlpkMip::_getSolValue() const {
946 return glp_mip_obj_val(lp);
949 GlpkMip* GlpkMip::newSolver() const { return new GlpkMip; }
950 GlpkMip* GlpkMip::cloneSolver() const {return new GlpkMip(*this); }
952 const char* GlpkMip::_solverName() const { return "GlpkMip"; }
954 void GlpkMip::messageLevel(MessageLevel m) {
958 } //END OF NAMESPACE LEMON