lemon/glpk.cc
author Peter Kovacs <kpeter@inf.elte.hu>
Tue, 24 Mar 2009 00:18:25 +0100
changeset 596 8c3112a66878
parent 461 08d495d48089
child 525 0fec6a017ead
child 528 9db62975c32b
permissions -rw-r--r--
Use XTI implementation instead of ATI in NetworkSimplex (#234)

XTI (eXtended Threaded Index) is an imporved version of the widely
known ATI (Augmented Threaded Index) method for storing and updating
the spanning tree structure in Network Simplex algorithms.

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