lemon/glpk.cc
author Alpar Juttner <alpar@cs.elte.hu>
Wed, 25 Aug 2010 21:34:58 +0200
changeset 927 16d466589b54
parent 566 e7017ec2d5cd
child 746 e4554cd6b2bf
child 988 8d281761dea4
permissions -rw-r--r--
Backport the CMAKE related changesets from the main branch

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