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