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