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