COIN-OR::LEMON - Graph Library

source: lemon/lemon/glpk.cc @ 1337:4add05447ca0

Last change on this file since 1337:4add05447ca0 was 1336:0759d974de81, checked in by Gabor Gevay <ggab90@…>, 10 years ago

STL style iterators (#325)

For

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