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