COIN-OR::LEMON - Graph Library

source: lemon-1.2/lemon/glpk.cc @ 955:8d281761dea4

Last change on this file since 955:8d281761dea4 was 955:8d281761dea4, checked in by Alpar Juttner <alpar@…>, 7 years ago

Fix buggy reinitialization in _solver_bits::VarIndex::clear() (#441)

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