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