COIN-OR::LEMON - Graph Library

source: lemon-1.2/lemon/glpk.cc @ 565:7ab97e2a0c33

Last change on this file since 565:7ab97e2a0c33 was 565:7ab97e2a0c33, checked in by Balazs Dezso <deba@…>, 11 years ago

Fixing presolver and basis handling (#255)

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