lemon/lp_glpk.cc
changeset 461 08d495d48089
parent 458 7afc121e0689
equal deleted inserted replaced
1:00264c0fde1e -1:000000000000
     1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
       
     2  *
       
     3  * This file is a part of LEMON, a generic C++ optimization library.
       
     4  *
       
     5  * Copyright (C) 2003-2008
       
     6  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
       
     7  * (Egervary Research Group on Combinatorial Optimization, EGRES).
       
     8  *
       
     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.
       
    12  *
       
    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
       
    15  * purpose.
       
    16  *
       
    17  */
       
    18 
       
    19 ///\file
       
    20 ///\brief Implementation of the LEMON GLPK LP and MIP solver interface.
       
    21 
       
    22 #include <lemon/lp_glpk.h>
       
    23 #include <glpk.h>
       
    24 
       
    25 #include <lemon/assert.h>
       
    26 
       
    27 namespace lemon {
       
    28 
       
    29   // GlpkBase members
       
    30 
       
    31   GlpkBase::GlpkBase() : LpBase() {
       
    32     lp = glp_create_prob();
       
    33     glp_create_index(lp);
       
    34   }
       
    35 
       
    36   GlpkBase::GlpkBase(const GlpkBase &other) : LpBase() {
       
    37     lp = glp_create_prob();
       
    38     glp_copy_prob(lp, other.lp, GLP_ON);
       
    39     glp_create_index(lp);
       
    40     rows = other.rows;
       
    41     cols = other.cols;
       
    42   }
       
    43 
       
    44   GlpkBase::~GlpkBase() {
       
    45     glp_delete_prob(lp);
       
    46   }
       
    47 
       
    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);
       
    51     return i;
       
    52   }
       
    53 
       
    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);
       
    57     return i;
       
    58   }
       
    59 
       
    60   void GlpkBase::_eraseCol(int i) {
       
    61     int ca[2];
       
    62     ca[1] = i;
       
    63     glp_del_cols(lp, 1, ca);
       
    64   }
       
    65 
       
    66   void GlpkBase::_eraseRow(int i) {
       
    67     int ra[2];
       
    68     ra[1] = i;
       
    69     glp_del_rows(lp, 1, ra);
       
    70   }
       
    71 
       
    72   void GlpkBase::_eraseColId(int i) {
       
    73     cols.eraseIndex(i);
       
    74     cols.shiftIndices(i);
       
    75   }
       
    76 
       
    77   void GlpkBase::_eraseRowId(int i) {
       
    78     rows.eraseIndex(i);
       
    79     rows.shiftIndices(i);
       
    80   }
       
    81 
       
    82   void GlpkBase::_getColName(int c, std::string& name) const {
       
    83     const char *str = glp_get_col_name(lp, c);
       
    84     if (str) name = str;
       
    85     else name.clear();
       
    86   }
       
    87 
       
    88   void GlpkBase::_setColName(int c, const std::string & name) {
       
    89     glp_set_col_name(lp, c, const_cast<char*>(name.c_str()));
       
    90 
       
    91   }
       
    92 
       
    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;
       
    96   }
       
    97 
       
    98   void GlpkBase::_getRowName(int r, std::string& name) const {
       
    99     const char *str = glp_get_row_name(lp, r);
       
   100     if (str) name = str;
       
   101     else name.clear();
       
   102   }
       
   103 
       
   104   void GlpkBase::_setRowName(int r, const std::string & name) {
       
   105     glp_set_row_name(lp, r, const_cast<char*>(name.c_str()));
       
   106 
       
   107   }
       
   108 
       
   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;
       
   112   }
       
   113 
       
   114   void GlpkBase::_setRowCoeffs(int i, ExprIterator b, ExprIterator e) {
       
   115     std::vector<int> indexes;
       
   116     std::vector<Value> values;
       
   117 
       
   118     indexes.push_back(0);
       
   119     values.push_back(0);
       
   120 
       
   121     for(ExprIterator it = b; it != e; ++it) {
       
   122       indexes.push_back(it->first);
       
   123       values.push_back(it->second);
       
   124     }
       
   125 
       
   126     glp_set_mat_row(lp, i, values.size() - 1,
       
   127                     &indexes.front(), &values.front());
       
   128   }
       
   129 
       
   130   void GlpkBase::_getRowCoeffs(int ix, InsertIterator b) const {
       
   131     int length = glp_get_mat_row(lp, ix, 0, 0);
       
   132 
       
   133     std::vector<int> indexes(length + 1);
       
   134     std::vector<Value> values(length + 1);
       
   135 
       
   136     glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
       
   137 
       
   138     for (int i = 1; i <= length; ++i) {
       
   139       *b = std::make_pair(indexes[i], values[i]);
       
   140       ++b;
       
   141     }
       
   142   }
       
   143 
       
   144   void GlpkBase::_setColCoeffs(int ix, ExprIterator b,
       
   145                                      ExprIterator e) {
       
   146 
       
   147     std::vector<int> indexes;
       
   148     std::vector<Value> values;
       
   149 
       
   150     indexes.push_back(0);
       
   151     values.push_back(0);
       
   152 
       
   153     for(ExprIterator it = b; it != e; ++it) {
       
   154       indexes.push_back(it->first);
       
   155       values.push_back(it->second);
       
   156     }
       
   157 
       
   158     glp_set_mat_col(lp, ix, values.size() - 1,
       
   159                     &indexes.front(), &values.front());
       
   160   }
       
   161 
       
   162   void GlpkBase::_getColCoeffs(int ix, InsertIterator b) const {
       
   163     int length = glp_get_mat_col(lp, ix, 0, 0);
       
   164 
       
   165     std::vector<int> indexes(length + 1);
       
   166     std::vector<Value> values(length + 1);
       
   167 
       
   168     glp_get_mat_col(lp, ix, &indexes.front(), &values.front());
       
   169 
       
   170     for (int i = 1; i  <= length; ++i) {
       
   171       *b = std::make_pair(indexes[i], values[i]);
       
   172       ++b;
       
   173     }
       
   174   }
       
   175 
       
   176   void GlpkBase::_setCoeff(int ix, int jx, Value value) {
       
   177 
       
   178     if (glp_get_num_cols(lp) < glp_get_num_rows(lp)) {
       
   179 
       
   180       int length = glp_get_mat_row(lp, ix, 0, 0);
       
   181 
       
   182       std::vector<int> indexes(length + 2);
       
   183       std::vector<Value> values(length + 2);
       
   184 
       
   185       glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
       
   186 
       
   187       //The following code does not suppose that the elements of the
       
   188       //array indexes are sorted
       
   189       bool found = false;
       
   190       for (int i = 1; i  <= length; ++i) {
       
   191         if (indexes[i] == jx) {
       
   192           found = true;
       
   193           values[i] = value;
       
   194           break;
       
   195         }
       
   196       }
       
   197       if (!found) {
       
   198         ++length;
       
   199         indexes[length] = jx;
       
   200         values[length] = value;
       
   201       }
       
   202 
       
   203       glp_set_mat_row(lp, ix, length, &indexes.front(), &values.front());
       
   204 
       
   205     } else {
       
   206 
       
   207       int length = glp_get_mat_col(lp, jx, 0, 0);
       
   208 
       
   209       std::vector<int> indexes(length + 2);
       
   210       std::vector<Value> values(length + 2);
       
   211 
       
   212       glp_get_mat_col(lp, jx, &indexes.front(), &values.front());
       
   213 
       
   214       //The following code does not suppose that the elements of the
       
   215       //array indexes are sorted
       
   216       bool found = false;
       
   217       for (int i = 1; i <= length; ++i) {
       
   218         if (indexes[i] == ix) {
       
   219           found = true;
       
   220           values[i] = value;
       
   221           break;
       
   222         }
       
   223       }
       
   224       if (!found) {
       
   225         ++length;
       
   226         indexes[length] = ix;
       
   227         values[length] = value;
       
   228       }
       
   229 
       
   230       glp_set_mat_col(lp, jx, length, &indexes.front(), &values.front());
       
   231     }
       
   232 
       
   233   }
       
   234 
       
   235   GlpkBase::Value GlpkBase::_getCoeff(int ix, int jx) const {
       
   236 
       
   237     int length = glp_get_mat_row(lp, ix, 0, 0);
       
   238 
       
   239     std::vector<int> indexes(length + 1);
       
   240     std::vector<Value> values(length + 1);
       
   241 
       
   242     glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
       
   243 
       
   244     for (int i = 1; i  <= length; ++i) {
       
   245       if (indexes[i] == jx) {
       
   246         return values[i];
       
   247       }
       
   248     }
       
   249 
       
   250     return 0;
       
   251   }
       
   252 
       
   253   void GlpkBase::_setColLowerBound(int i, Value lo) {
       
   254     LEMON_ASSERT(lo != INF, "Invalid bound");
       
   255 
       
   256     int b = glp_get_col_type(lp, i);
       
   257     double up = glp_get_col_ub(lp, i);
       
   258     if (lo == -INF) {
       
   259       switch (b) {
       
   260       case GLP_FR:
       
   261       case GLP_LO:
       
   262         glp_set_col_bnds(lp, i, GLP_FR, lo, up);
       
   263         break;
       
   264       case GLP_UP:
       
   265         break;
       
   266       case GLP_DB:
       
   267       case GLP_FX:
       
   268         glp_set_col_bnds(lp, i, GLP_UP, lo, up);
       
   269         break;
       
   270       default:
       
   271         break;
       
   272       }
       
   273     } else {
       
   274       switch (b) {
       
   275       case GLP_FR:
       
   276       case GLP_LO:
       
   277         glp_set_col_bnds(lp, i, GLP_LO, lo, up);
       
   278         break;
       
   279       case GLP_UP:
       
   280       case GLP_DB:
       
   281       case GLP_FX:
       
   282         if (lo == up)
       
   283           glp_set_col_bnds(lp, i, GLP_FX, lo, up);
       
   284         else
       
   285           glp_set_col_bnds(lp, i, GLP_DB, lo, up);
       
   286         break;
       
   287       default:
       
   288         break;
       
   289       }
       
   290     }
       
   291   }
       
   292 
       
   293   GlpkBase::Value GlpkBase::_getColLowerBound(int i) const {
       
   294     int b = glp_get_col_type(lp, i);
       
   295     switch (b) {
       
   296     case GLP_LO:
       
   297     case GLP_DB:
       
   298     case GLP_FX:
       
   299       return glp_get_col_lb(lp, i);
       
   300     default:
       
   301       return -INF;
       
   302     }
       
   303   }
       
   304 
       
   305   void GlpkBase::_setColUpperBound(int i, Value up) {
       
   306     LEMON_ASSERT(up != -INF, "Invalid bound");
       
   307 
       
   308     int b = glp_get_col_type(lp, i);
       
   309     double lo = glp_get_col_lb(lp, i);
       
   310     if (up == INF) {
       
   311       switch (b) {
       
   312       case GLP_FR:
       
   313       case GLP_LO:
       
   314         break;
       
   315       case GLP_UP:
       
   316         glp_set_col_bnds(lp, i, GLP_FR, lo, up);
       
   317         break;
       
   318       case GLP_DB:
       
   319       case GLP_FX:
       
   320         glp_set_col_bnds(lp, i, GLP_LO, lo, up);
       
   321         break;
       
   322       default:
       
   323         break;
       
   324       }
       
   325     } else {
       
   326       switch (b) {
       
   327       case GLP_FR:
       
   328         glp_set_col_bnds(lp, i, GLP_UP, lo, up);
       
   329         break;
       
   330       case GLP_UP:
       
   331         glp_set_col_bnds(lp, i, GLP_UP, lo, up);
       
   332         break;
       
   333       case GLP_LO:
       
   334       case GLP_DB:
       
   335       case GLP_FX:
       
   336         if (lo == up)
       
   337           glp_set_col_bnds(lp, i, GLP_FX, lo, up);
       
   338         else
       
   339           glp_set_col_bnds(lp, i, GLP_DB, lo, up);
       
   340         break;
       
   341       default:
       
   342         break;
       
   343       }
       
   344     }
       
   345 
       
   346   }
       
   347 
       
   348   GlpkBase::Value GlpkBase::_getColUpperBound(int i) const {
       
   349     int b = glp_get_col_type(lp, i);
       
   350       switch (b) {
       
   351       case GLP_UP:
       
   352       case GLP_DB:
       
   353       case GLP_FX:
       
   354         return glp_get_col_ub(lp, i);
       
   355       default:
       
   356         return INF;
       
   357       }
       
   358   }
       
   359 
       
   360   void GlpkBase::_setRowLowerBound(int i, Value lo) {
       
   361     LEMON_ASSERT(lo != INF, "Invalid bound");
       
   362 
       
   363     int b = glp_get_row_type(lp, i);
       
   364     double up = glp_get_row_ub(lp, i);
       
   365     if (lo == -INF) {
       
   366       switch (b) {
       
   367       case GLP_FR:
       
   368       case GLP_LO:
       
   369         glp_set_row_bnds(lp, i, GLP_FR, lo, up);
       
   370         break;
       
   371       case GLP_UP:
       
   372         break;
       
   373       case GLP_DB:
       
   374       case GLP_FX:
       
   375         glp_set_row_bnds(lp, i, GLP_UP, lo, up);
       
   376         break;
       
   377       default:
       
   378         break;
       
   379       }
       
   380     } else {
       
   381       switch (b) {
       
   382       case GLP_FR:
       
   383       case GLP_LO:
       
   384         glp_set_row_bnds(lp, i, GLP_LO, lo, up);
       
   385         break;
       
   386       case GLP_UP:
       
   387       case GLP_DB:
       
   388       case GLP_FX:
       
   389         if (lo == up)
       
   390           glp_set_row_bnds(lp, i, GLP_FX, lo, up);
       
   391         else
       
   392           glp_set_row_bnds(lp, i, GLP_DB, lo, up);
       
   393         break;
       
   394       default:
       
   395         break;
       
   396       }
       
   397     }
       
   398 
       
   399   }
       
   400 
       
   401   GlpkBase::Value GlpkBase::_getRowLowerBound(int i) const {
       
   402     int b = glp_get_row_type(lp, i);
       
   403     switch (b) {
       
   404     case GLP_LO:
       
   405     case GLP_DB:
       
   406     case GLP_FX:
       
   407       return glp_get_row_lb(lp, i);
       
   408     default:
       
   409       return -INF;
       
   410     }
       
   411   }
       
   412 
       
   413   void GlpkBase::_setRowUpperBound(int i, Value up) {
       
   414     LEMON_ASSERT(up != -INF, "Invalid bound");
       
   415 
       
   416     int b = glp_get_row_type(lp, i);
       
   417     double lo = glp_get_row_lb(lp, i);
       
   418     if (up == INF) {
       
   419       switch (b) {
       
   420       case GLP_FR:
       
   421       case GLP_LO:
       
   422         break;
       
   423       case GLP_UP:
       
   424         glp_set_row_bnds(lp, i, GLP_FR, lo, up);
       
   425         break;
       
   426       case GLP_DB:
       
   427       case GLP_FX:
       
   428         glp_set_row_bnds(lp, i, GLP_LO, lo, up);
       
   429         break;
       
   430       default:
       
   431         break;
       
   432       }
       
   433     } else {
       
   434       switch (b) {
       
   435       case GLP_FR:
       
   436         glp_set_row_bnds(lp, i, GLP_UP, lo, up);
       
   437         break;
       
   438       case GLP_UP:
       
   439         glp_set_row_bnds(lp, i, GLP_UP, lo, up);
       
   440         break;
       
   441       case GLP_LO:
       
   442       case GLP_DB:
       
   443       case GLP_FX:
       
   444         if (lo == up)
       
   445           glp_set_row_bnds(lp, i, GLP_FX, lo, up);
       
   446         else
       
   447           glp_set_row_bnds(lp, i, GLP_DB, lo, up);
       
   448         break;
       
   449       default:
       
   450         break;
       
   451       }
       
   452     }
       
   453   }
       
   454 
       
   455   GlpkBase::Value GlpkBase::_getRowUpperBound(int i) const {
       
   456     int b = glp_get_row_type(lp, i);
       
   457     switch (b) {
       
   458     case GLP_UP:
       
   459     case GLP_DB:
       
   460     case GLP_FX:
       
   461       return glp_get_row_ub(lp, i);
       
   462     default:
       
   463       return INF;
       
   464     }
       
   465   }
       
   466 
       
   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);
       
   470     }
       
   471     for (ExprIterator it = b; it != e; ++it) {
       
   472       glp_set_obj_coef(lp, it->first, it->second);
       
   473     }
       
   474   }
       
   475 
       
   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);
       
   479       if (val != 0.0) {
       
   480         *b = std::make_pair(i, val);
       
   481         ++b;
       
   482       }
       
   483     }
       
   484   }
       
   485 
       
   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);
       
   489   }
       
   490 
       
   491   GlpkBase::Value GlpkBase::_getObjCoeff(int i) const {
       
   492     //i = 0 means the constant term (shift)
       
   493     return glp_get_obj_coef(lp, i);
       
   494   }
       
   495 
       
   496   void GlpkBase::_setSense(GlpkBase::Sense sense) {
       
   497     switch (sense) {
       
   498     case MIN:
       
   499       glp_set_obj_dir(lp, GLP_MIN);
       
   500       break;
       
   501     case MAX:
       
   502       glp_set_obj_dir(lp, GLP_MAX);
       
   503       break;
       
   504     }
       
   505   }
       
   506 
       
   507   GlpkBase::Sense GlpkBase::_getSense() const {
       
   508     switch(glp_get_obj_dir(lp)) {
       
   509     case GLP_MIN:
       
   510       return MIN;
       
   511     case GLP_MAX:
       
   512       return MAX;
       
   513     default:
       
   514       LEMON_ASSERT(false, "Wrong sense");
       
   515       return GlpkBase::Sense();
       
   516     }
       
   517   }
       
   518 
       
   519   void GlpkBase::_clear() {
       
   520     glp_erase_prob(lp);
       
   521     rows.clear();
       
   522     cols.clear();
       
   523   }
       
   524 
       
   525   // LpGlpk members
       
   526 
       
   527   LpGlpk::LpGlpk()
       
   528     : LpBase(), GlpkBase(), LpSolver() {
       
   529     messageLevel(MESSAGE_NO_OUTPUT);
       
   530   }
       
   531 
       
   532   LpGlpk::LpGlpk(const LpGlpk& other)
       
   533     : LpBase(other), GlpkBase(other), LpSolver(other) {
       
   534     messageLevel(MESSAGE_NO_OUTPUT);
       
   535   }
       
   536 
       
   537   LpGlpk* LpGlpk::_newSolver() const { return new LpGlpk; }
       
   538   LpGlpk* LpGlpk::_cloneSolver() const { return new LpGlpk(*this); }
       
   539 
       
   540   const char* LpGlpk::_solverName() const { return "LpGlpk"; }
       
   541 
       
   542   void LpGlpk::_clear_temporals() {
       
   543     _primal_ray.clear();
       
   544     _dual_ray.clear();
       
   545   }
       
   546 
       
   547   LpGlpk::SolveExitStatus LpGlpk::_solve() {
       
   548     return solvePrimal();
       
   549   }
       
   550 
       
   551   LpGlpk::SolveExitStatus LpGlpk::solvePrimal() {
       
   552     _clear_temporals();
       
   553 
       
   554     glp_smcp smcp;
       
   555     glp_init_smcp(&smcp);
       
   556 
       
   557     switch (_message_level) {
       
   558     case MESSAGE_NO_OUTPUT:
       
   559       smcp.msg_lev = GLP_MSG_OFF;
       
   560       break;
       
   561     case MESSAGE_ERROR_MESSAGE:
       
   562       smcp.msg_lev = GLP_MSG_ERR;
       
   563       break;
       
   564     case MESSAGE_NORMAL_OUTPUT:
       
   565       smcp.msg_lev = GLP_MSG_ON;
       
   566       break;
       
   567     case MESSAGE_FULL_OUTPUT:
       
   568       smcp.msg_lev = GLP_MSG_ALL;
       
   569       break;
       
   570     }
       
   571 
       
   572     if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
       
   573     return SOLVED;
       
   574   }
       
   575 
       
   576   LpGlpk::SolveExitStatus LpGlpk::solveDual() {
       
   577     _clear_temporals();
       
   578 
       
   579     glp_smcp smcp;
       
   580     glp_init_smcp(&smcp);
       
   581 
       
   582     switch (_message_level) {
       
   583     case MESSAGE_NO_OUTPUT:
       
   584       smcp.msg_lev = GLP_MSG_OFF;
       
   585       break;
       
   586     case MESSAGE_ERROR_MESSAGE:
       
   587       smcp.msg_lev = GLP_MSG_ERR;
       
   588       break;
       
   589     case MESSAGE_NORMAL_OUTPUT:
       
   590       smcp.msg_lev = GLP_MSG_ON;
       
   591       break;
       
   592     case MESSAGE_FULL_OUTPUT:
       
   593       smcp.msg_lev = GLP_MSG_ALL;
       
   594       break;
       
   595     }
       
   596     smcp.meth = GLP_DUAL;
       
   597 
       
   598     if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
       
   599     return SOLVED;
       
   600   }
       
   601 
       
   602   LpGlpk::Value LpGlpk::_getPrimal(int i) const {
       
   603     return glp_get_col_prim(lp, i);
       
   604   }
       
   605 
       
   606   LpGlpk::Value LpGlpk::_getDual(int i) const {
       
   607     return glp_get_row_dual(lp, i);
       
   608   }
       
   609 
       
   610   LpGlpk::Value LpGlpk::_getPrimalValue() const {
       
   611     return glp_get_obj_val(lp);
       
   612   }
       
   613 
       
   614   LpGlpk::VarStatus LpGlpk::_getColStatus(int i) const {
       
   615     switch (glp_get_col_stat(lp, i)) {
       
   616     case GLP_BS:
       
   617       return BASIC;
       
   618     case GLP_UP:
       
   619       return UPPER;
       
   620     case GLP_LO:
       
   621       return LOWER;
       
   622     case GLP_NF:
       
   623       return FREE;
       
   624     case GLP_NS:
       
   625       return FIXED;
       
   626     default:
       
   627       LEMON_ASSERT(false, "Wrong column status");
       
   628       return LpGlpk::VarStatus();
       
   629     }
       
   630   }
       
   631 
       
   632   LpGlpk::VarStatus LpGlpk::_getRowStatus(int i) const {
       
   633     switch (glp_get_row_stat(lp, i)) {
       
   634     case GLP_BS:
       
   635       return BASIC;
       
   636     case GLP_UP:
       
   637       return UPPER;
       
   638     case GLP_LO:
       
   639       return LOWER;
       
   640     case GLP_NF:
       
   641       return FREE;
       
   642     case GLP_NS:
       
   643       return FIXED;
       
   644     default:
       
   645       LEMON_ASSERT(false, "Wrong row status");
       
   646       return LpGlpk::VarStatus();
       
   647     }
       
   648   }
       
   649 
       
   650   LpGlpk::Value LpGlpk::_getPrimalRay(int i) const {
       
   651     if (_primal_ray.empty()) {
       
   652       int row_num = glp_get_num_rows(lp);
       
   653       int col_num = glp_get_num_cols(lp);
       
   654 
       
   655       _primal_ray.resize(col_num + 1, 0.0);
       
   656 
       
   657       int index = glp_get_unbnd_ray(lp);
       
   658       if (index != 0) {
       
   659         // The primal ray is found in primal simplex second phase
       
   660         LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :
       
   661                       glp_get_col_stat(lp, index - row_num)) != GLP_BS,
       
   662                      "Wrong primal ray");
       
   663 
       
   664         bool negate = glp_get_obj_dir(lp) == GLP_MAX;
       
   665 
       
   666         if (index > row_num) {
       
   667           _primal_ray[index - row_num] = 1.0;
       
   668           if (glp_get_col_dual(lp, index - row_num) > 0) {
       
   669             negate = !negate;
       
   670           }
       
   671         } else {
       
   672           if (glp_get_row_dual(lp, index) > 0) {
       
   673             negate = !negate;
       
   674           }
       
   675         }
       
   676 
       
   677         std::vector<int> ray_indexes(row_num + 1);
       
   678         std::vector<Value> ray_values(row_num + 1);
       
   679         int ray_length = glp_eval_tab_col(lp, index, &ray_indexes.front(),
       
   680                                           &ray_values.front());
       
   681 
       
   682         for (int i = 1; i <= ray_length; ++i) {
       
   683           if (ray_indexes[i] > row_num) {
       
   684             _primal_ray[ray_indexes[i] - row_num] = ray_values[i];
       
   685           }
       
   686         }
       
   687 
       
   688         if (negate) {
       
   689           for (int i = 1; i <= col_num; ++i) {
       
   690             _primal_ray[i] = - _primal_ray[i];
       
   691           }
       
   692         }
       
   693       } else {
       
   694         for (int i = 1; i <= col_num; ++i) {
       
   695           _primal_ray[i] = glp_get_col_prim(lp, i);
       
   696         }
       
   697       }
       
   698     }
       
   699     return _primal_ray[i];
       
   700   }
       
   701 
       
   702   LpGlpk::Value LpGlpk::_getDualRay(int i) const {
       
   703     if (_dual_ray.empty()) {
       
   704       int row_num = glp_get_num_rows(lp);
       
   705 
       
   706       _dual_ray.resize(row_num + 1, 0.0);
       
   707 
       
   708       int index = glp_get_unbnd_ray(lp);
       
   709       if (index != 0) {
       
   710         // The dual ray is found in dual simplex second phase
       
   711         LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :
       
   712                       glp_get_col_stat(lp, index - row_num)) == GLP_BS,
       
   713 
       
   714                      "Wrong dual ray");
       
   715 
       
   716         int idx;
       
   717         bool negate = false;
       
   718 
       
   719         if (index > row_num) {
       
   720           idx = glp_get_col_bind(lp, index - row_num);
       
   721           if (glp_get_col_prim(lp, index - row_num) >
       
   722               glp_get_col_ub(lp, index - row_num)) {
       
   723             negate = true;
       
   724           }
       
   725         } else {
       
   726           idx = glp_get_row_bind(lp, index);
       
   727           if (glp_get_row_prim(lp, index) > glp_get_row_ub(lp, index)) {
       
   728             negate = true;
       
   729           }
       
   730         }
       
   731 
       
   732         _dual_ray[idx] = negate ?  - 1.0 : 1.0;
       
   733 
       
   734         glp_btran(lp, &_dual_ray.front());
       
   735       } else {
       
   736         double eps = 1e-7;
       
   737         // The dual ray is found in primal simplex first phase
       
   738         // We assume that the glpk minimizes the slack to get feasible solution
       
   739         for (int i = 1; i <= row_num; ++i) {
       
   740           int index = glp_get_bhead(lp, i);
       
   741           if (index <= row_num) {
       
   742             double res = glp_get_row_prim(lp, index);
       
   743             if (res > glp_get_row_ub(lp, index) + eps) {
       
   744               _dual_ray[i] = -1;
       
   745             } else if (res < glp_get_row_lb(lp, index) - eps) {
       
   746               _dual_ray[i] = 1;
       
   747             } else {
       
   748               _dual_ray[i] = 0;
       
   749             }
       
   750             _dual_ray[i] *= glp_get_rii(lp, index);
       
   751           } else {
       
   752             double res = glp_get_col_prim(lp, index - row_num);
       
   753             if (res > glp_get_col_ub(lp, index - row_num) + eps) {
       
   754               _dual_ray[i] = -1;
       
   755             } else if (res < glp_get_col_lb(lp, index - row_num) - eps) {
       
   756               _dual_ray[i] = 1;
       
   757             } else {
       
   758               _dual_ray[i] = 0;
       
   759             }
       
   760             _dual_ray[i] /= glp_get_sjj(lp, index - row_num);
       
   761           }
       
   762         }
       
   763 
       
   764         glp_btran(lp, &_dual_ray.front());
       
   765 
       
   766         for (int i = 1; i <= row_num; ++i) {
       
   767           _dual_ray[i] /= glp_get_rii(lp, i);
       
   768         }
       
   769       }
       
   770     }
       
   771     return _dual_ray[i];
       
   772   }
       
   773 
       
   774   LpGlpk::ProblemType LpGlpk::_getPrimalType() const {
       
   775     if (glp_get_status(lp) == GLP_OPT)
       
   776       return OPTIMAL;
       
   777     switch (glp_get_prim_stat(lp)) {
       
   778     case GLP_UNDEF:
       
   779       return UNDEFINED;
       
   780     case GLP_FEAS:
       
   781     case GLP_INFEAS:
       
   782       if (glp_get_dual_stat(lp) == GLP_NOFEAS) {
       
   783         return UNBOUNDED;
       
   784       } else {
       
   785         return UNDEFINED;
       
   786       }
       
   787     case GLP_NOFEAS:
       
   788       return INFEASIBLE;
       
   789     default:
       
   790       LEMON_ASSERT(false, "Wrong primal type");
       
   791       return  LpGlpk::ProblemType();
       
   792     }
       
   793   }
       
   794 
       
   795   LpGlpk::ProblemType LpGlpk::_getDualType() const {
       
   796     if (glp_get_status(lp) == GLP_OPT)
       
   797       return OPTIMAL;
       
   798     switch (glp_get_dual_stat(lp)) {
       
   799     case GLP_UNDEF:
       
   800       return UNDEFINED;
       
   801     case GLP_FEAS:
       
   802     case GLP_INFEAS:
       
   803       if (glp_get_prim_stat(lp) == GLP_NOFEAS) {
       
   804         return UNBOUNDED;
       
   805       } else {
       
   806         return UNDEFINED;
       
   807       }
       
   808     case GLP_NOFEAS:
       
   809       return INFEASIBLE;
       
   810     default:
       
   811       LEMON_ASSERT(false, "Wrong primal type");
       
   812       return  LpGlpk::ProblemType();
       
   813     }
       
   814   }
       
   815 
       
   816   void LpGlpk::presolver(bool b) {
       
   817     lpx_set_int_parm(lp, LPX_K_PRESOL, b ? 1 : 0);
       
   818   }
       
   819 
       
   820   void LpGlpk::messageLevel(MessageLevel m) {
       
   821     _message_level = m;
       
   822   }
       
   823 
       
   824   // MipGlpk members
       
   825 
       
   826   MipGlpk::MipGlpk()
       
   827     : LpBase(), GlpkBase(), MipSolver() {
       
   828     messageLevel(MESSAGE_NO_OUTPUT);
       
   829   }
       
   830 
       
   831   MipGlpk::MipGlpk(const MipGlpk& other)
       
   832     : LpBase(), GlpkBase(other), MipSolver() {
       
   833     messageLevel(MESSAGE_NO_OUTPUT);
       
   834   }
       
   835 
       
   836   void MipGlpk::_setColType(int i, MipGlpk::ColTypes col_type) {
       
   837     switch (col_type) {
       
   838     case INTEGER:
       
   839       glp_set_col_kind(lp, i, GLP_IV);
       
   840       break;
       
   841     case REAL:
       
   842       glp_set_col_kind(lp, i, GLP_CV);
       
   843       break;
       
   844     }
       
   845   }
       
   846 
       
   847   MipGlpk::ColTypes MipGlpk::_getColType(int i) const {
       
   848     switch (glp_get_col_kind(lp, i)) {
       
   849     case GLP_IV:
       
   850     case GLP_BV:
       
   851       return INTEGER;
       
   852     default:
       
   853       return REAL;
       
   854     }
       
   855 
       
   856   }
       
   857 
       
   858   MipGlpk::SolveExitStatus MipGlpk::_solve() {
       
   859     glp_smcp smcp;
       
   860     glp_init_smcp(&smcp);
       
   861 
       
   862     switch (_message_level) {
       
   863     case MESSAGE_NO_OUTPUT:
       
   864       smcp.msg_lev = GLP_MSG_OFF;
       
   865       break;
       
   866     case MESSAGE_ERROR_MESSAGE:
       
   867       smcp.msg_lev = GLP_MSG_ERR;
       
   868       break;
       
   869     case MESSAGE_NORMAL_OUTPUT:
       
   870       smcp.msg_lev = GLP_MSG_ON;
       
   871       break;
       
   872     case MESSAGE_FULL_OUTPUT:
       
   873       smcp.msg_lev = GLP_MSG_ALL;
       
   874       break;
       
   875     }
       
   876     smcp.meth = GLP_DUAL;
       
   877 
       
   878     if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
       
   879     if (glp_get_status(lp) != GLP_OPT) return SOLVED;
       
   880 
       
   881     glp_iocp iocp;
       
   882     glp_init_iocp(&iocp);
       
   883 
       
   884     switch (_message_level) {
       
   885     case MESSAGE_NO_OUTPUT:
       
   886       iocp.msg_lev = GLP_MSG_OFF;
       
   887       break;
       
   888     case MESSAGE_ERROR_MESSAGE:
       
   889       iocp.msg_lev = GLP_MSG_ERR;
       
   890       break;
       
   891     case MESSAGE_NORMAL_OUTPUT:
       
   892       iocp.msg_lev = GLP_MSG_ON;
       
   893       break;
       
   894     case MESSAGE_FULL_OUTPUT:
       
   895       iocp.msg_lev = GLP_MSG_ALL;
       
   896       break;
       
   897     }
       
   898 
       
   899     if (glp_intopt(lp, &iocp) != 0) return UNSOLVED;
       
   900     return SOLVED;
       
   901   }
       
   902 
       
   903 
       
   904   MipGlpk::ProblemType MipGlpk::_getType() const {
       
   905     switch (glp_get_status(lp)) {
       
   906     case GLP_OPT:
       
   907       switch (glp_mip_status(lp)) {
       
   908       case GLP_UNDEF:
       
   909         return UNDEFINED;
       
   910       case GLP_NOFEAS:
       
   911         return INFEASIBLE;
       
   912       case GLP_FEAS:
       
   913         return FEASIBLE;
       
   914       case GLP_OPT:
       
   915         return OPTIMAL;
       
   916       default:
       
   917         LEMON_ASSERT(false, "Wrong problem type.");
       
   918         return MipGlpk::ProblemType();
       
   919       }
       
   920     case GLP_NOFEAS:
       
   921       return INFEASIBLE;
       
   922     case GLP_INFEAS:
       
   923     case GLP_FEAS:
       
   924       if (glp_get_dual_stat(lp) == GLP_NOFEAS) {
       
   925         return UNBOUNDED;
       
   926       } else {
       
   927         return UNDEFINED;
       
   928       }
       
   929     default:
       
   930       LEMON_ASSERT(false, "Wrong problem type.");
       
   931       return MipGlpk::ProblemType();
       
   932     }
       
   933   }
       
   934 
       
   935   MipGlpk::Value MipGlpk::_getSol(int i) const {
       
   936     return glp_mip_col_val(lp, i);
       
   937   }
       
   938 
       
   939   MipGlpk::Value MipGlpk::_getSolValue() const {
       
   940     return glp_mip_obj_val(lp);
       
   941   }
       
   942 
       
   943   MipGlpk* MipGlpk::_newSolver() const { return new MipGlpk; }
       
   944   MipGlpk* MipGlpk::_cloneSolver() const {return new MipGlpk(*this); }
       
   945 
       
   946   const char* MipGlpk::_solverName() const { return "MipGlpk"; }
       
   947 
       
   948   void MipGlpk::messageLevel(MessageLevel m) {
       
   949     _message_level = m;
       
   950   }
       
   951 
       
   952 } //END OF NAMESPACE LEMON