lemon/glpk.cc
author Alpar Juttner <alpar@cs.elte.hu>
Wed, 07 Aug 2013 06:55:05 +0200
changeset 983 8b2d4e5d96e4
parent 877 141f9c0db4a3
parent 955 8d281761dea4
permissions -rw-r--r--
Merge #294 to branches >=1.2
     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-2010
     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   }
   560 
   561   void GlpkBase::freeEnv() {
   562     glp_free_env();
   563   }
   564 
   565   void GlpkBase::_messageLevel(MessageLevel level) {
   566     switch (level) {
   567     case MESSAGE_NOTHING:
   568       _message_level = GLP_MSG_OFF;
   569       break;
   570     case MESSAGE_ERROR:
   571       _message_level = GLP_MSG_ERR;
   572       break;
   573     case MESSAGE_WARNING:
   574       _message_level = GLP_MSG_ERR;
   575       break;
   576     case MESSAGE_NORMAL:
   577       _message_level = GLP_MSG_ON;
   578       break;
   579     case MESSAGE_VERBOSE:
   580       _message_level = GLP_MSG_ALL;
   581       break;
   582     }
   583   }
   584 
   585   GlpkBase::FreeEnvHelper GlpkBase::freeEnvHelper;
   586 
   587   // GlpkLp members
   588 
   589   GlpkLp::GlpkLp()
   590     : LpBase(), LpSolver(), GlpkBase() {
   591     presolver(false);
   592   }
   593 
   594   GlpkLp::GlpkLp(const GlpkLp& other)
   595     : LpBase(other), LpSolver(other), GlpkBase(other) {
   596     presolver(false);
   597   }
   598 
   599   GlpkLp* GlpkLp::newSolver() const { return new GlpkLp; }
   600   GlpkLp* GlpkLp::cloneSolver() const { return new GlpkLp(*this); }
   601 
   602   const char* GlpkLp::_solverName() const { return "GlpkLp"; }
   603 
   604   void GlpkLp::_clear_temporals() {
   605     _primal_ray.clear();
   606     _dual_ray.clear();
   607   }
   608 
   609   GlpkLp::SolveExitStatus GlpkLp::_solve() {
   610     return solvePrimal();
   611   }
   612 
   613   GlpkLp::SolveExitStatus GlpkLp::solvePrimal() {
   614     _clear_temporals();
   615 
   616     glp_smcp smcp;
   617     glp_init_smcp(&smcp);
   618 
   619     smcp.msg_lev = _message_level;
   620     smcp.presolve = _presolve;
   621 
   622     // If the basis is not valid we get an error return value.
   623     // In this case we can try to create a new basis.
   624     switch (glp_simplex(lp, &smcp)) {
   625     case 0:
   626       break;
   627     case GLP_EBADB:
   628     case GLP_ESING:
   629     case GLP_ECOND:
   630       glp_term_out(false);
   631       glp_adv_basis(lp, 0);
   632       glp_term_out(true);
   633       if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
   634       break;
   635     default:
   636       return UNSOLVED;
   637     }
   638 
   639     return SOLVED;
   640   }
   641 
   642   GlpkLp::SolveExitStatus GlpkLp::solveDual() {
   643     _clear_temporals();
   644 
   645     glp_smcp smcp;
   646     glp_init_smcp(&smcp);
   647 
   648     smcp.msg_lev = _message_level;
   649     smcp.meth = GLP_DUAL;
   650     smcp.presolve = _presolve;
   651 
   652     // If the basis is not valid we get an error return value.
   653     // In this case we can try to create a new basis.
   654     switch (glp_simplex(lp, &smcp)) {
   655     case 0:
   656       break;
   657     case GLP_EBADB:
   658     case GLP_ESING:
   659     case GLP_ECOND:
   660       glp_term_out(false);
   661       glp_adv_basis(lp, 0);
   662       glp_term_out(true);
   663       if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
   664       break;
   665     default:
   666       return UNSOLVED;
   667     }
   668     return SOLVED;
   669   }
   670 
   671   GlpkLp::Value GlpkLp::_getPrimal(int i) const {
   672     return glp_get_col_prim(lp, i);
   673   }
   674 
   675   GlpkLp::Value GlpkLp::_getDual(int i) const {
   676     return glp_get_row_dual(lp, i);
   677   }
   678 
   679   GlpkLp::Value GlpkLp::_getPrimalValue() const {
   680     return glp_get_obj_val(lp);
   681   }
   682 
   683   GlpkLp::VarStatus GlpkLp::_getColStatus(int i) const {
   684     switch (glp_get_col_stat(lp, i)) {
   685     case GLP_BS:
   686       return BASIC;
   687     case GLP_UP:
   688       return UPPER;
   689     case GLP_LO:
   690       return LOWER;
   691     case GLP_NF:
   692       return FREE;
   693     case GLP_NS:
   694       return FIXED;
   695     default:
   696       LEMON_ASSERT(false, "Wrong column status");
   697       return GlpkLp::VarStatus();
   698     }
   699   }
   700 
   701   GlpkLp::VarStatus GlpkLp::_getRowStatus(int i) const {
   702     switch (glp_get_row_stat(lp, i)) {
   703     case GLP_BS:
   704       return BASIC;
   705     case GLP_UP:
   706       return UPPER;
   707     case GLP_LO:
   708       return LOWER;
   709     case GLP_NF:
   710       return FREE;
   711     case GLP_NS:
   712       return FIXED;
   713     default:
   714       LEMON_ASSERT(false, "Wrong row status");
   715       return GlpkLp::VarStatus();
   716     }
   717   }
   718 
   719   GlpkLp::Value GlpkLp::_getPrimalRay(int i) const {
   720     if (_primal_ray.empty()) {
   721       int row_num = glp_get_num_rows(lp);
   722       int col_num = glp_get_num_cols(lp);
   723 
   724       _primal_ray.resize(col_num + 1, 0.0);
   725 
   726       int index = glp_get_unbnd_ray(lp);
   727       if (index != 0) {
   728         // The primal ray is found in primal simplex second phase
   729         LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :
   730                       glp_get_col_stat(lp, index - row_num)) != GLP_BS,
   731                      "Wrong primal ray");
   732 
   733         bool negate = glp_get_obj_dir(lp) == GLP_MAX;
   734 
   735         if (index > row_num) {
   736           _primal_ray[index - row_num] = 1.0;
   737           if (glp_get_col_dual(lp, index - row_num) > 0) {
   738             negate = !negate;
   739           }
   740         } else {
   741           if (glp_get_row_dual(lp, index) > 0) {
   742             negate = !negate;
   743           }
   744         }
   745 
   746         std::vector<int> ray_indexes(row_num + 1);
   747         std::vector<Value> ray_values(row_num + 1);
   748         int ray_length = glp_eval_tab_col(lp, index, &ray_indexes.front(),
   749                                           &ray_values.front());
   750 
   751         for (int i = 1; i <= ray_length; ++i) {
   752           if (ray_indexes[i] > row_num) {
   753             _primal_ray[ray_indexes[i] - row_num] = ray_values[i];
   754           }
   755         }
   756 
   757         if (negate) {
   758           for (int i = 1; i <= col_num; ++i) {
   759             _primal_ray[i] = - _primal_ray[i];
   760           }
   761         }
   762       } else {
   763         for (int i = 1; i <= col_num; ++i) {
   764           _primal_ray[i] = glp_get_col_prim(lp, i);
   765         }
   766       }
   767     }
   768     return _primal_ray[i];
   769   }
   770 
   771   GlpkLp::Value GlpkLp::_getDualRay(int i) const {
   772     if (_dual_ray.empty()) {
   773       int row_num = glp_get_num_rows(lp);
   774 
   775       _dual_ray.resize(row_num + 1, 0.0);
   776 
   777       int index = glp_get_unbnd_ray(lp);
   778       if (index != 0) {
   779         // The dual ray is found in dual simplex second phase
   780         LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :
   781                       glp_get_col_stat(lp, index - row_num)) == GLP_BS,
   782 
   783                      "Wrong dual ray");
   784 
   785         int idx;
   786         bool negate = false;
   787 
   788         if (index > row_num) {
   789           idx = glp_get_col_bind(lp, index - row_num);
   790           if (glp_get_col_prim(lp, index - row_num) >
   791               glp_get_col_ub(lp, index - row_num)) {
   792             negate = true;
   793           }
   794         } else {
   795           idx = glp_get_row_bind(lp, index);
   796           if (glp_get_row_prim(lp, index) > glp_get_row_ub(lp, index)) {
   797             negate = true;
   798           }
   799         }
   800 
   801         _dual_ray[idx] = negate ?  - 1.0 : 1.0;
   802 
   803         glp_btran(lp, &_dual_ray.front());
   804       } else {
   805         double eps = 1e-7;
   806         // The dual ray is found in primal simplex first phase
   807         // We assume that the glpk minimizes the slack to get feasible solution
   808         for (int i = 1; i <= row_num; ++i) {
   809           int index = glp_get_bhead(lp, i);
   810           if (index <= row_num) {
   811             double res = glp_get_row_prim(lp, index);
   812             if (res > glp_get_row_ub(lp, index) + eps) {
   813               _dual_ray[i] = -1;
   814             } else if (res < glp_get_row_lb(lp, index) - eps) {
   815               _dual_ray[i] = 1;
   816             } else {
   817               _dual_ray[i] = 0;
   818             }
   819             _dual_ray[i] *= glp_get_rii(lp, index);
   820           } else {
   821             double res = glp_get_col_prim(lp, index - row_num);
   822             if (res > glp_get_col_ub(lp, index - row_num) + eps) {
   823               _dual_ray[i] = -1;
   824             } else if (res < glp_get_col_lb(lp, index - row_num) - eps) {
   825               _dual_ray[i] = 1;
   826             } else {
   827               _dual_ray[i] = 0;
   828             }
   829             _dual_ray[i] /= glp_get_sjj(lp, index - row_num);
   830           }
   831         }
   832 
   833         glp_btran(lp, &_dual_ray.front());
   834 
   835         for (int i = 1; i <= row_num; ++i) {
   836           _dual_ray[i] /= glp_get_rii(lp, i);
   837         }
   838       }
   839     }
   840     return _dual_ray[i];
   841   }
   842 
   843   GlpkLp::ProblemType GlpkLp::_getPrimalType() const {
   844     if (glp_get_status(lp) == GLP_OPT)
   845       return OPTIMAL;
   846     switch (glp_get_prim_stat(lp)) {
   847     case GLP_UNDEF:
   848       return UNDEFINED;
   849     case GLP_FEAS:
   850     case GLP_INFEAS:
   851       if (glp_get_dual_stat(lp) == GLP_NOFEAS) {
   852         return UNBOUNDED;
   853       } else {
   854         return UNDEFINED;
   855       }
   856     case GLP_NOFEAS:
   857       return INFEASIBLE;
   858     default:
   859       LEMON_ASSERT(false, "Wrong primal type");
   860       return  GlpkLp::ProblemType();
   861     }
   862   }
   863 
   864   GlpkLp::ProblemType GlpkLp::_getDualType() const {
   865     if (glp_get_status(lp) == GLP_OPT)
   866       return OPTIMAL;
   867     switch (glp_get_dual_stat(lp)) {
   868     case GLP_UNDEF:
   869       return UNDEFINED;
   870     case GLP_FEAS:
   871     case GLP_INFEAS:
   872       if (glp_get_prim_stat(lp) == GLP_NOFEAS) {
   873         return UNBOUNDED;
   874       } else {
   875         return UNDEFINED;
   876       }
   877     case GLP_NOFEAS:
   878       return INFEASIBLE;
   879     default:
   880       LEMON_ASSERT(false, "Wrong primal type");
   881       return  GlpkLp::ProblemType();
   882     }
   883   }
   884 
   885   void GlpkLp::presolver(bool presolve) {
   886     _presolve = presolve;
   887   }
   888 
   889   // GlpkMip members
   890 
   891   GlpkMip::GlpkMip()
   892     : LpBase(), MipSolver(), GlpkBase() {
   893   }
   894 
   895   GlpkMip::GlpkMip(const GlpkMip& other)
   896     : LpBase(), MipSolver(), GlpkBase(other) {
   897   }
   898 
   899   void GlpkMip::_setColType(int i, GlpkMip::ColTypes col_type) {
   900     switch (col_type) {
   901     case INTEGER:
   902       glp_set_col_kind(lp, i, GLP_IV);
   903       break;
   904     case REAL:
   905       glp_set_col_kind(lp, i, GLP_CV);
   906       break;
   907     }
   908   }
   909 
   910   GlpkMip::ColTypes GlpkMip::_getColType(int i) const {
   911     switch (glp_get_col_kind(lp, i)) {
   912     case GLP_IV:
   913     case GLP_BV:
   914       return INTEGER;
   915     default:
   916       return REAL;
   917     }
   918 
   919   }
   920 
   921   GlpkMip::SolveExitStatus GlpkMip::_solve() {
   922     glp_smcp smcp;
   923     glp_init_smcp(&smcp);
   924 
   925     smcp.msg_lev = _message_level;
   926     smcp.meth = GLP_DUAL;
   927 
   928     // If the basis is not valid we get an error return value.
   929     // In this case we can try to create a new basis.
   930     switch (glp_simplex(lp, &smcp)) {
   931     case 0:
   932       break;
   933     case GLP_EBADB:
   934     case GLP_ESING:
   935     case GLP_ECOND:
   936       glp_term_out(false);
   937       glp_adv_basis(lp, 0);
   938       glp_term_out(true);
   939       if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
   940       break;
   941     default:
   942       return UNSOLVED;
   943     }
   944 
   945     if (glp_get_status(lp) != GLP_OPT) return SOLVED;
   946 
   947     glp_iocp iocp;
   948     glp_init_iocp(&iocp);
   949 
   950     iocp.msg_lev = _message_level;
   951 
   952     if (glp_intopt(lp, &iocp) != 0) return UNSOLVED;
   953     return SOLVED;
   954   }
   955 
   956 
   957   GlpkMip::ProblemType GlpkMip::_getType() const {
   958     switch (glp_get_status(lp)) {
   959     case GLP_OPT:
   960       switch (glp_mip_status(lp)) {
   961       case GLP_UNDEF:
   962         return UNDEFINED;
   963       case GLP_NOFEAS:
   964         return INFEASIBLE;
   965       case GLP_FEAS:
   966         return FEASIBLE;
   967       case GLP_OPT:
   968         return OPTIMAL;
   969       default:
   970         LEMON_ASSERT(false, "Wrong problem type.");
   971         return GlpkMip::ProblemType();
   972       }
   973     case GLP_NOFEAS:
   974       return INFEASIBLE;
   975     case GLP_INFEAS:
   976     case GLP_FEAS:
   977       if (glp_get_dual_stat(lp) == GLP_NOFEAS) {
   978         return UNBOUNDED;
   979       } else {
   980         return UNDEFINED;
   981       }
   982     default:
   983       LEMON_ASSERT(false, "Wrong problem type.");
   984       return GlpkMip::ProblemType();
   985     }
   986   }
   987 
   988   GlpkMip::Value GlpkMip::_getSol(int i) const {
   989     return glp_mip_col_val(lp, i);
   990   }
   991 
   992   GlpkMip::Value GlpkMip::_getSolValue() const {
   993     return glp_mip_obj_val(lp);
   994   }
   995 
   996   GlpkMip* GlpkMip::newSolver() const { return new GlpkMip; }
   997   GlpkMip* GlpkMip::cloneSolver() const {return new GlpkMip(*this); }
   998 
   999   const char* GlpkMip::_solverName() const { return "GlpkMip"; }
  1000 
  1001 } //END OF NAMESPACE LEMON