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

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

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