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