Fix wrong iteration in ListGraph snapshot, part II. (#598)
1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
3 * This file is a part of LEMON, a generic C++ optimization library.
5 * Copyright (C) 2003-2013
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();
34 messageLevel(MESSAGE_NOTHING);
37 GlpkBase::GlpkBase(const GlpkBase &other) : LpBase() {
38 lp = glp_create_prob();
39 glp_copy_prob(lp, other.lp, GLP_ON);
43 messageLevel(MESSAGE_NOTHING);
46 GlpkBase::~GlpkBase() {
50 int GlpkBase::_addCol() {
51 int i = glp_add_cols(lp, 1);
52 glp_set_col_bnds(lp, i, GLP_FR, 0.0, 0.0);
56 int GlpkBase::_addRow() {
57 int i = glp_add_rows(lp, 1);
58 glp_set_row_bnds(lp, i, GLP_FR, 0.0, 0.0);
62 int GlpkBase::_addRow(Value lo, ExprIterator b,
63 ExprIterator e, Value up) {
64 int i = glp_add_rows(lp, 1);
68 glp_set_row_bnds(lp, i, GLP_FR, lo, up);
70 glp_set_row_bnds(lp, i, GLP_UP, lo, up);
74 glp_set_row_bnds(lp, i, GLP_LO, lo, up);
75 } else if (lo != up) {
76 glp_set_row_bnds(lp, i, GLP_DB, lo, up);
78 glp_set_row_bnds(lp, i, GLP_FX, lo, up);
82 std::vector<int> indexes;
83 std::vector<Value> values;
88 for(ExprIterator it = b; it != e; ++it) {
89 indexes.push_back(it->first);
90 values.push_back(it->second);
93 glp_set_mat_row(lp, i, values.size() - 1,
94 &indexes.front(), &values.front());
98 void GlpkBase::_eraseCol(int i) {
101 glp_del_cols(lp, 1, ca);
104 void GlpkBase::_eraseRow(int i) {
107 glp_del_rows(lp, 1, ra);
110 void GlpkBase::_eraseColId(int i) {
112 cols.shiftIndices(i);
115 void GlpkBase::_eraseRowId(int i) {
117 rows.shiftIndices(i);
120 void GlpkBase::_getColName(int c, std::string& name) const {
121 const char *str = glp_get_col_name(lp, c);
126 void GlpkBase::_setColName(int c, const std::string & name) {
127 glp_set_col_name(lp, c, const_cast<char*>(name.c_str()));
131 int GlpkBase::_colByName(const std::string& name) const {
132 int k = glp_find_col(lp, const_cast<char*>(name.c_str()));
133 return k > 0 ? k : -1;
136 void GlpkBase::_getRowName(int r, std::string& name) const {
137 const char *str = glp_get_row_name(lp, r);
142 void GlpkBase::_setRowName(int r, const std::string & name) {
143 glp_set_row_name(lp, r, const_cast<char*>(name.c_str()));
147 int GlpkBase::_rowByName(const std::string& name) const {
148 int k = glp_find_row(lp, const_cast<char*>(name.c_str()));
149 return k > 0 ? k : -1;
152 void GlpkBase::_setRowCoeffs(int i, ExprIterator b, ExprIterator e) {
153 std::vector<int> indexes;
154 std::vector<Value> values;
156 indexes.push_back(0);
159 for(ExprIterator it = b; it != e; ++it) {
160 indexes.push_back(it->first);
161 values.push_back(it->second);
164 glp_set_mat_row(lp, i, values.size() - 1,
165 &indexes.front(), &values.front());
168 void GlpkBase::_getRowCoeffs(int ix, InsertIterator b) const {
169 int length = glp_get_mat_row(lp, ix, 0, 0);
171 std::vector<int> indexes(length + 1);
172 std::vector<Value> values(length + 1);
174 glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
176 for (int i = 1; i <= length; ++i) {
177 *b = std::make_pair(indexes[i], values[i]);
182 void GlpkBase::_setColCoeffs(int ix, ExprIterator b,
185 std::vector<int> indexes;
186 std::vector<Value> values;
188 indexes.push_back(0);
191 for(ExprIterator it = b; it != e; ++it) {
192 indexes.push_back(it->first);
193 values.push_back(it->second);
196 glp_set_mat_col(lp, ix, values.size() - 1,
197 &indexes.front(), &values.front());
200 void GlpkBase::_getColCoeffs(int ix, InsertIterator b) const {
201 int length = glp_get_mat_col(lp, ix, 0, 0);
203 std::vector<int> indexes(length + 1);
204 std::vector<Value> values(length + 1);
206 glp_get_mat_col(lp, ix, &indexes.front(), &values.front());
208 for (int i = 1; i <= length; ++i) {
209 *b = std::make_pair(indexes[i], values[i]);
214 void GlpkBase::_setCoeff(int ix, int jx, Value value) {
216 if (glp_get_num_cols(lp) < glp_get_num_rows(lp)) {
218 int length = glp_get_mat_row(lp, ix, 0, 0);
220 std::vector<int> indexes(length + 2);
221 std::vector<Value> values(length + 2);
223 glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
225 //The following code does not suppose that the elements of the
226 //array indexes are sorted
228 for (int i = 1; i <= length; ++i) {
229 if (indexes[i] == jx) {
237 indexes[length] = jx;
238 values[length] = value;
241 glp_set_mat_row(lp, ix, length, &indexes.front(), &values.front());
245 int length = glp_get_mat_col(lp, jx, 0, 0);
247 std::vector<int> indexes(length + 2);
248 std::vector<Value> values(length + 2);
250 glp_get_mat_col(lp, jx, &indexes.front(), &values.front());
252 //The following code does not suppose that the elements of the
253 //array indexes are sorted
255 for (int i = 1; i <= length; ++i) {
256 if (indexes[i] == ix) {
264 indexes[length] = ix;
265 values[length] = value;
268 glp_set_mat_col(lp, jx, length, &indexes.front(), &values.front());
273 GlpkBase::Value GlpkBase::_getCoeff(int ix, int jx) const {
275 int length = glp_get_mat_row(lp, ix, 0, 0);
277 std::vector<int> indexes(length + 1);
278 std::vector<Value> values(length + 1);
280 glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
282 for (int i = 1; i <= length; ++i) {
283 if (indexes[i] == jx) {
291 void GlpkBase::_setColLowerBound(int i, Value lo) {
292 LEMON_ASSERT(lo != INF, "Invalid bound");
294 int b = glp_get_col_type(lp, i);
295 double up = glp_get_col_ub(lp, i);
300 glp_set_col_bnds(lp, i, GLP_FR, lo, up);
306 glp_set_col_bnds(lp, i, GLP_UP, lo, up);
315 glp_set_col_bnds(lp, i, GLP_LO, lo, up);
321 glp_set_col_bnds(lp, i, GLP_FX, lo, up);
323 glp_set_col_bnds(lp, i, GLP_DB, lo, up);
331 GlpkBase::Value GlpkBase::_getColLowerBound(int i) const {
332 int b = glp_get_col_type(lp, i);
337 return glp_get_col_lb(lp, i);
343 void GlpkBase::_setColUpperBound(int i, Value up) {
344 LEMON_ASSERT(up != -INF, "Invalid bound");
346 int b = glp_get_col_type(lp, i);
347 double lo = glp_get_col_lb(lp, i);
354 glp_set_col_bnds(lp, i, GLP_FR, lo, up);
358 glp_set_col_bnds(lp, i, GLP_LO, lo, up);
366 glp_set_col_bnds(lp, i, GLP_UP, lo, up);
369 glp_set_col_bnds(lp, i, GLP_UP, lo, up);
375 glp_set_col_bnds(lp, i, GLP_FX, lo, up);
377 glp_set_col_bnds(lp, i, GLP_DB, lo, up);
386 GlpkBase::Value GlpkBase::_getColUpperBound(int i) const {
387 int b = glp_get_col_type(lp, i);
392 return glp_get_col_ub(lp, i);
398 void GlpkBase::_setRowLowerBound(int i, Value lo) {
399 LEMON_ASSERT(lo != INF, "Invalid bound");
401 int b = glp_get_row_type(lp, i);
402 double up = glp_get_row_ub(lp, i);
407 glp_set_row_bnds(lp, i, GLP_FR, lo, up);
413 glp_set_row_bnds(lp, i, GLP_UP, lo, up);
422 glp_set_row_bnds(lp, i, GLP_LO, lo, up);
428 glp_set_row_bnds(lp, i, GLP_FX, lo, up);
430 glp_set_row_bnds(lp, i, GLP_DB, lo, up);
439 GlpkBase::Value GlpkBase::_getRowLowerBound(int i) const {
440 int b = glp_get_row_type(lp, i);
445 return glp_get_row_lb(lp, i);
451 void GlpkBase::_setRowUpperBound(int i, Value up) {
452 LEMON_ASSERT(up != -INF, "Invalid bound");
454 int b = glp_get_row_type(lp, i);
455 double lo = glp_get_row_lb(lp, i);
462 glp_set_row_bnds(lp, i, GLP_FR, lo, up);
466 glp_set_row_bnds(lp, i, GLP_LO, lo, up);
474 glp_set_row_bnds(lp, i, GLP_UP, lo, up);
477 glp_set_row_bnds(lp, i, GLP_UP, lo, up);
483 glp_set_row_bnds(lp, i, GLP_FX, lo, up);
485 glp_set_row_bnds(lp, i, GLP_DB, lo, up);
493 GlpkBase::Value GlpkBase::_getRowUpperBound(int i) const {
494 int b = glp_get_row_type(lp, i);
499 return glp_get_row_ub(lp, i);
505 void GlpkBase::_setObjCoeffs(ExprIterator b, ExprIterator e) {
506 for (int i = 1; i <= glp_get_num_cols(lp); ++i) {
507 glp_set_obj_coef(lp, i, 0.0);
509 for (ExprIterator it = b; it != e; ++it) {
510 glp_set_obj_coef(lp, it->first, it->second);
514 void GlpkBase::_getObjCoeffs(InsertIterator b) const {
515 for (int i = 1; i <= glp_get_num_cols(lp); ++i) {
516 Value val = glp_get_obj_coef(lp, i);
518 *b = std::make_pair(i, val);
524 void GlpkBase::_setObjCoeff(int i, Value obj_coef) {
525 //i = 0 means the constant term (shift)
526 glp_set_obj_coef(lp, i, obj_coef);
529 GlpkBase::Value GlpkBase::_getObjCoeff(int i) const {
530 //i = 0 means the constant term (shift)
531 return glp_get_obj_coef(lp, i);
534 void GlpkBase::_setSense(GlpkBase::Sense sense) {
537 glp_set_obj_dir(lp, GLP_MIN);
540 glp_set_obj_dir(lp, GLP_MAX);
545 GlpkBase::Sense GlpkBase::_getSense() const {
546 switch(glp_get_obj_dir(lp)) {
552 LEMON_ASSERT(false, "Wrong sense");
553 return GlpkBase::Sense();
557 void GlpkBase::_clear() {
561 void GlpkBase::freeEnv() {
565 void GlpkBase::_messageLevel(MessageLevel level) {
567 case MESSAGE_NOTHING:
568 _message_level = GLP_MSG_OFF;
571 _message_level = GLP_MSG_ERR;
573 case MESSAGE_WARNING:
574 _message_level = GLP_MSG_ERR;
577 _message_level = GLP_MSG_ON;
579 case MESSAGE_VERBOSE:
580 _message_level = GLP_MSG_ALL;
585 void GlpkBase::_write(std::string file, std::string format) const
588 glp_write_mps(lp, GLP_MPS_FILE, 0, file.c_str());
589 else if(format == "LP")
590 glp_write_lp(lp, 0, file.c_str());
591 else throw UnsupportedFormatError(format);
594 GlpkBase::FreeEnvHelper GlpkBase::freeEnvHelper;
599 : LpBase(), LpSolver(), GlpkBase() {
603 GlpkLp::GlpkLp(const GlpkLp& other)
604 : LpBase(other), LpSolver(other), GlpkBase(other) {
608 GlpkLp* GlpkLp::newSolver() const { return new GlpkLp; }
609 GlpkLp* GlpkLp::cloneSolver() const { return new GlpkLp(*this); }
611 const char* GlpkLp::_solverName() const { return "GlpkLp"; }
613 void GlpkLp::_clear_temporals() {
618 GlpkLp::SolveExitStatus GlpkLp::_solve() {
619 return solvePrimal();
622 GlpkLp::SolveExitStatus GlpkLp::solvePrimal() {
626 glp_init_smcp(&smcp);
628 smcp.msg_lev = _message_level;
629 smcp.presolve = _presolve;
631 // If the basis is not valid we get an error return value.
632 // In this case we can try to create a new basis.
633 switch (glp_simplex(lp, &smcp)) {
640 glp_adv_basis(lp, 0);
642 if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
651 GlpkLp::SolveExitStatus GlpkLp::solveDual() {
655 glp_init_smcp(&smcp);
657 smcp.msg_lev = _message_level;
658 smcp.meth = GLP_DUAL;
659 smcp.presolve = _presolve;
661 // If the basis is not valid we get an error return value.
662 // In this case we can try to create a new basis.
663 switch (glp_simplex(lp, &smcp)) {
670 glp_adv_basis(lp, 0);
672 if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
680 GlpkLp::Value GlpkLp::_getPrimal(int i) const {
681 return glp_get_col_prim(lp, i);
684 GlpkLp::Value GlpkLp::_getDual(int i) const {
685 return glp_get_row_dual(lp, i);
688 GlpkLp::Value GlpkLp::_getPrimalValue() const {
689 return glp_get_obj_val(lp);
692 GlpkLp::VarStatus GlpkLp::_getColStatus(int i) const {
693 switch (glp_get_col_stat(lp, i)) {
705 LEMON_ASSERT(false, "Wrong column status");
706 return GlpkLp::VarStatus();
710 GlpkLp::VarStatus GlpkLp::_getRowStatus(int i) const {
711 switch (glp_get_row_stat(lp, i)) {
723 LEMON_ASSERT(false, "Wrong row status");
724 return GlpkLp::VarStatus();
728 GlpkLp::Value GlpkLp::_getPrimalRay(int i) const {
729 if (_primal_ray.empty()) {
730 int row_num = glp_get_num_rows(lp);
731 int col_num = glp_get_num_cols(lp);
733 _primal_ray.resize(col_num + 1, 0.0);
735 int index = glp_get_unbnd_ray(lp);
737 // The primal ray is found in primal simplex second phase
738 LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :
739 glp_get_col_stat(lp, index - row_num)) != GLP_BS,
742 bool negate = glp_get_obj_dir(lp) == GLP_MAX;
744 if (index > row_num) {
745 _primal_ray[index - row_num] = 1.0;
746 if (glp_get_col_dual(lp, index - row_num) > 0) {
750 if (glp_get_row_dual(lp, index) > 0) {
755 std::vector<int> ray_indexes(row_num + 1);
756 std::vector<Value> ray_values(row_num + 1);
757 int ray_length = glp_eval_tab_col(lp, index, &ray_indexes.front(),
758 &ray_values.front());
760 for (int i = 1; i <= ray_length; ++i) {
761 if (ray_indexes[i] > row_num) {
762 _primal_ray[ray_indexes[i] - row_num] = ray_values[i];
767 for (int i = 1; i <= col_num; ++i) {
768 _primal_ray[i] = - _primal_ray[i];
772 for (int i = 1; i <= col_num; ++i) {
773 _primal_ray[i] = glp_get_col_prim(lp, i);
777 return _primal_ray[i];
780 GlpkLp::Value GlpkLp::_getDualRay(int i) const {
781 if (_dual_ray.empty()) {
782 int row_num = glp_get_num_rows(lp);
784 _dual_ray.resize(row_num + 1, 0.0);
786 int index = glp_get_unbnd_ray(lp);
788 // The dual ray is found in dual simplex second phase
789 LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :
790 glp_get_col_stat(lp, index - row_num)) == GLP_BS,
797 if (index > row_num) {
798 idx = glp_get_col_bind(lp, index - row_num);
799 if (glp_get_col_prim(lp, index - row_num) >
800 glp_get_col_ub(lp, index - row_num)) {
804 idx = glp_get_row_bind(lp, index);
805 if (glp_get_row_prim(lp, index) > glp_get_row_ub(lp, index)) {
810 _dual_ray[idx] = negate ? - 1.0 : 1.0;
812 glp_btran(lp, &_dual_ray.front());
815 // The dual ray is found in primal simplex first phase
816 // We assume that the glpk minimizes the slack to get feasible solution
817 for (int i = 1; i <= row_num; ++i) {
818 int index = glp_get_bhead(lp, i);
819 if (index <= row_num) {
820 double res = glp_get_row_prim(lp, index);
821 if (res > glp_get_row_ub(lp, index) + eps) {
823 } else if (res < glp_get_row_lb(lp, index) - eps) {
828 _dual_ray[i] *= glp_get_rii(lp, index);
830 double res = glp_get_col_prim(lp, index - row_num);
831 if (res > glp_get_col_ub(lp, index - row_num) + eps) {
833 } else if (res < glp_get_col_lb(lp, index - row_num) - eps) {
838 _dual_ray[i] /= glp_get_sjj(lp, index - row_num);
842 glp_btran(lp, &_dual_ray.front());
844 for (int i = 1; i <= row_num; ++i) {
845 _dual_ray[i] /= glp_get_rii(lp, i);
852 GlpkLp::ProblemType GlpkLp::_getPrimalType() const {
853 if (glp_get_status(lp) == GLP_OPT)
855 switch (glp_get_prim_stat(lp)) {
860 if (glp_get_dual_stat(lp) == GLP_NOFEAS) {
868 LEMON_ASSERT(false, "Wrong primal type");
869 return GlpkLp::ProblemType();
873 GlpkLp::ProblemType GlpkLp::_getDualType() const {
874 if (glp_get_status(lp) == GLP_OPT)
876 switch (glp_get_dual_stat(lp)) {
881 if (glp_get_prim_stat(lp) == GLP_NOFEAS) {
889 LEMON_ASSERT(false, "Wrong primal type");
890 return GlpkLp::ProblemType();
894 void GlpkLp::presolver(bool presolve) {
895 _presolve = presolve;
901 : LpBase(), MipSolver(), GlpkBase() {
904 GlpkMip::GlpkMip(const GlpkMip& other)
905 : LpBase(), MipSolver(), GlpkBase(other) {
908 void GlpkMip::_setColType(int i, GlpkMip::ColTypes col_type) {
911 glp_set_col_kind(lp, i, GLP_IV);
914 glp_set_col_kind(lp, i, GLP_CV);
919 GlpkMip::ColTypes GlpkMip::_getColType(int i) const {
920 switch (glp_get_col_kind(lp, i)) {
930 GlpkMip::SolveExitStatus GlpkMip::_solve() {
932 glp_init_smcp(&smcp);
934 smcp.msg_lev = _message_level;
935 smcp.meth = GLP_DUAL;
937 // If the basis is not valid we get an error return value.
938 // In this case we can try to create a new basis.
939 switch (glp_simplex(lp, &smcp)) {
946 glp_adv_basis(lp, 0);
948 if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
954 if (glp_get_status(lp) != GLP_OPT) return SOLVED;
957 glp_init_iocp(&iocp);
959 iocp.msg_lev = _message_level;
961 if (glp_intopt(lp, &iocp) != 0) return UNSOLVED;
966 GlpkMip::ProblemType GlpkMip::_getType() const {
967 switch (glp_get_status(lp)) {
969 switch (glp_mip_status(lp)) {
979 LEMON_ASSERT(false, "Wrong problem type.");
980 return GlpkMip::ProblemType();
986 if (glp_get_dual_stat(lp) == GLP_NOFEAS) {
992 LEMON_ASSERT(false, "Wrong problem type.");
993 return GlpkMip::ProblemType();
997 GlpkMip::Value GlpkMip::_getSol(int i) const {
998 return glp_mip_col_val(lp, i);
1001 GlpkMip::Value GlpkMip::_getSolValue() const {
1002 return glp_mip_obj_val(lp);
1005 GlpkMip* GlpkMip::newSolver() const { return new GlpkMip; }
1006 GlpkMip* GlpkMip::cloneSolver() const {return new GlpkMip(*this); }
1008 const char* GlpkMip::_solverName() const { return "GlpkMip"; }
1012 } //END OF NAMESPACE LEMON