gravatar
deba@inf.elte.hu
deba@inf.elte.hu
Automatic GLPK env deallocation (#213)
0 4 0
default
4 files changed with 14 insertions and 10 deletions:
↑ Collapse diff ↑
Ignore white space 6 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2008
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
///\file
20 20
///\brief Implementation of the LEMON GLPK LP and MIP solver interface.
21 21

	
22 22
#include <lemon/glpk.h>
23 23
#include <glpk.h>
24 24

	
25 25
#include <lemon/assert.h>
26 26

	
27 27
namespace lemon {
28 28

	
29 29
  // GlpkBase members
30 30

	
31 31
  GlpkBase::GlpkBase() : LpBase() {
32 32
    lp = glp_create_prob();
33 33
    glp_create_index(lp);
34 34
  }
35 35

	
36 36
  GlpkBase::GlpkBase(const GlpkBase &other) : LpBase() {
37 37
    lp = glp_create_prob();
38 38
    glp_copy_prob(lp, other.lp, GLP_ON);
39 39
    glp_create_index(lp);
40 40
    rows = other.rows;
41 41
    cols = other.cols;
42 42
  }
43 43

	
44 44
  GlpkBase::~GlpkBase() {
45 45
    glp_delete_prob(lp);
46 46
  }
47 47

	
48 48
  int GlpkBase::_addCol() {
49 49
    int i = glp_add_cols(lp, 1);
50 50
    glp_set_col_bnds(lp, i, GLP_FR, 0.0, 0.0);
51 51
    return i;
52 52
  }
53 53

	
54 54
  int GlpkBase::_addRow() {
55 55
    int i = glp_add_rows(lp, 1);
56 56
    glp_set_row_bnds(lp, i, GLP_FR, 0.0, 0.0);
57 57
    return i;
58 58
  }
59 59

	
60 60
  void GlpkBase::_eraseCol(int i) {
61 61
    int ca[2];
62 62
    ca[1] = i;
63 63
    glp_del_cols(lp, 1, ca);
64 64
  }
65 65

	
66 66
  void GlpkBase::_eraseRow(int i) {
67 67
    int ra[2];
68 68
    ra[1] = i;
69 69
    glp_del_rows(lp, 1, ra);
70 70
  }
71 71

	
72 72
  void GlpkBase::_eraseColId(int i) {
73 73
    cols.eraseIndex(i);
74 74
    cols.shiftIndices(i);
75 75
  }
76 76

	
77 77
  void GlpkBase::_eraseRowId(int i) {
78 78
    rows.eraseIndex(i);
79 79
    rows.shiftIndices(i);
80 80
  }
81 81

	
82 82
  void GlpkBase::_getColName(int c, std::string& name) const {
83 83
    const char *str = glp_get_col_name(lp, c);
84 84
    if (str) name = str;
85 85
    else name.clear();
86 86
  }
87 87

	
88 88
  void GlpkBase::_setColName(int c, const std::string & name) {
89 89
    glp_set_col_name(lp, c, const_cast<char*>(name.c_str()));
90 90

	
91 91
  }
92 92

	
93 93
  int GlpkBase::_colByName(const std::string& name) const {
94 94
    int k = glp_find_col(lp, const_cast<char*>(name.c_str()));
95 95
    return k > 0 ? k : -1;
96 96
  }
97 97

	
98 98
  void GlpkBase::_getRowName(int r, std::string& name) const {
99 99
    const char *str = glp_get_row_name(lp, r);
100 100
    if (str) name = str;
101 101
    else name.clear();
102 102
  }
103 103

	
104 104
  void GlpkBase::_setRowName(int r, const std::string & name) {
105 105
    glp_set_row_name(lp, r, const_cast<char*>(name.c_str()));
106 106

	
107 107
  }
108 108

	
109 109
  int GlpkBase::_rowByName(const std::string& name) const {
110 110
    int k = glp_find_row(lp, const_cast<char*>(name.c_str()));
111 111
    return k > 0 ? k : -1;
112 112
  }
113 113

	
114 114
  void GlpkBase::_setRowCoeffs(int i, ExprIterator b, ExprIterator e) {
115 115
    std::vector<int> indexes;
116 116
    std::vector<Value> values;
117 117

	
118 118
    indexes.push_back(0);
119 119
    values.push_back(0);
120 120

	
121 121
    for(ExprIterator it = b; it != e; ++it) {
122 122
      indexes.push_back(it->first);
123 123
      values.push_back(it->second);
124 124
    }
125 125

	
126 126
    glp_set_mat_row(lp, i, values.size() - 1,
127 127
                    &indexes.front(), &values.front());
128 128
  }
129 129

	
130 130
  void GlpkBase::_getRowCoeffs(int ix, InsertIterator b) const {
131 131
    int length = glp_get_mat_row(lp, ix, 0, 0);
132 132

	
133 133
    std::vector<int> indexes(length + 1);
134 134
    std::vector<Value> values(length + 1);
135 135

	
136 136
    glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
137 137

	
138 138
    for (int i = 1; i <= length; ++i) {
139 139
      *b = std::make_pair(indexes[i], values[i]);
140 140
      ++b;
141 141
    }
142 142
  }
143 143

	
144 144
  void GlpkBase::_setColCoeffs(int ix, ExprIterator b,
145 145
                                     ExprIterator e) {
146 146

	
147 147
    std::vector<int> indexes;
148 148
    std::vector<Value> values;
149 149

	
150 150
    indexes.push_back(0);
151 151
    values.push_back(0);
152 152

	
153 153
    for(ExprIterator it = b; it != e; ++it) {
154 154
      indexes.push_back(it->first);
155 155
      values.push_back(it->second);
156 156
    }
157 157

	
158 158
    glp_set_mat_col(lp, ix, values.size() - 1,
159 159
                    &indexes.front(), &values.front());
160 160
  }
161 161

	
162 162
  void GlpkBase::_getColCoeffs(int ix, InsertIterator b) const {
163 163
    int length = glp_get_mat_col(lp, ix, 0, 0);
164 164

	
165 165
    std::vector<int> indexes(length + 1);
166 166
    std::vector<Value> values(length + 1);
167 167

	
168 168
    glp_get_mat_col(lp, ix, &indexes.front(), &values.front());
169 169

	
170 170
    for (int i = 1; i  <= length; ++i) {
171 171
      *b = std::make_pair(indexes[i], values[i]);
172 172
      ++b;
173 173
    }
174 174
  }
175 175

	
176 176
  void GlpkBase::_setCoeff(int ix, int jx, Value value) {
177 177

	
178 178
    if (glp_get_num_cols(lp) < glp_get_num_rows(lp)) {
179 179

	
180 180
      int length = glp_get_mat_row(lp, ix, 0, 0);
181 181

	
182 182
      std::vector<int> indexes(length + 2);
183 183
      std::vector<Value> values(length + 2);
184 184

	
185 185
      glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
186 186

	
187 187
      //The following code does not suppose that the elements of the
188 188
      //array indexes are sorted
189 189
      bool found = false;
190 190
      for (int i = 1; i  <= length; ++i) {
191 191
        if (indexes[i] == jx) {
192 192
          found = true;
193 193
          values[i] = value;
194 194
          break;
195 195
        }
196 196
      }
197 197
      if (!found) {
198 198
        ++length;
199 199
        indexes[length] = jx;
200 200
        values[length] = value;
201 201
      }
202 202

	
203 203
      glp_set_mat_row(lp, ix, length, &indexes.front(), &values.front());
204 204

	
205 205
    } else {
206 206

	
207 207
      int length = glp_get_mat_col(lp, jx, 0, 0);
208 208

	
209 209
      std::vector<int> indexes(length + 2);
210 210
      std::vector<Value> values(length + 2);
211 211

	
212 212
      glp_get_mat_col(lp, jx, &indexes.front(), &values.front());
213 213

	
214 214
      //The following code does not suppose that the elements of the
215 215
      //array indexes are sorted
216 216
      bool found = false;
217 217
      for (int i = 1; i <= length; ++i) {
218 218
        if (indexes[i] == ix) {
219 219
          found = true;
220 220
          values[i] = value;
221 221
          break;
222 222
        }
223 223
      }
224 224
      if (!found) {
225 225
        ++length;
226 226
        indexes[length] = ix;
227 227
        values[length] = value;
228 228
      }
229 229

	
230 230
      glp_set_mat_col(lp, jx, length, &indexes.front(), &values.front());
231 231
    }
232 232

	
233 233
  }
234 234

	
235 235
  GlpkBase::Value GlpkBase::_getCoeff(int ix, int jx) const {
236 236

	
237 237
    int length = glp_get_mat_row(lp, ix, 0, 0);
238 238

	
239 239
    std::vector<int> indexes(length + 1);
240 240
    std::vector<Value> values(length + 1);
241 241

	
242 242
    glp_get_mat_row(lp, ix, &indexes.front(), &values.front());
243 243

	
244 244
    for (int i = 1; i  <= length; ++i) {
245 245
      if (indexes[i] == jx) {
246 246
        return values[i];
247 247
      }
248 248
    }
249 249

	
250 250
    return 0;
251 251
  }
252 252

	
253 253
  void GlpkBase::_setColLowerBound(int i, Value lo) {
254 254
    LEMON_ASSERT(lo != INF, "Invalid bound");
255 255

	
256 256
    int b = glp_get_col_type(lp, i);
257 257
    double up = glp_get_col_ub(lp, i);
258 258
    if (lo == -INF) {
259 259
      switch (b) {
260 260
      case GLP_FR:
261 261
      case GLP_LO:
262 262
        glp_set_col_bnds(lp, i, GLP_FR, lo, up);
263 263
        break;
264 264
      case GLP_UP:
265 265
        break;
266 266
      case GLP_DB:
267 267
      case GLP_FX:
268 268
        glp_set_col_bnds(lp, i, GLP_UP, lo, up);
269 269
        break;
270 270
      default:
271 271
        break;
272 272
      }
273 273
    } else {
274 274
      switch (b) {
275 275
      case GLP_FR:
276 276
      case GLP_LO:
277 277
        glp_set_col_bnds(lp, i, GLP_LO, lo, up);
278 278
        break;
279 279
      case GLP_UP:
280 280
      case GLP_DB:
281 281
      case GLP_FX:
282 282
        if (lo == up)
283 283
          glp_set_col_bnds(lp, i, GLP_FX, lo, up);
284 284
        else
285 285
          glp_set_col_bnds(lp, i, GLP_DB, lo, up);
286 286
        break;
287 287
      default:
288 288
        break;
289 289
      }
290 290
    }
291 291
  }
292 292

	
293 293
  GlpkBase::Value GlpkBase::_getColLowerBound(int i) const {
294 294
    int b = glp_get_col_type(lp, i);
295 295
    switch (b) {
296 296
    case GLP_LO:
297 297
    case GLP_DB:
298 298
    case GLP_FX:
299 299
      return glp_get_col_lb(lp, i);
300 300
    default:
301 301
      return -INF;
302 302
    }
303 303
  }
304 304

	
305 305
  void GlpkBase::_setColUpperBound(int i, Value up) {
306 306
    LEMON_ASSERT(up != -INF, "Invalid bound");
307 307

	
308 308
    int b = glp_get_col_type(lp, i);
309 309
    double lo = glp_get_col_lb(lp, i);
310 310
    if (up == INF) {
311 311
      switch (b) {
312 312
      case GLP_FR:
313 313
      case GLP_LO:
314 314
        break;
315 315
      case GLP_UP:
316 316
        glp_set_col_bnds(lp, i, GLP_FR, lo, up);
317 317
        break;
318 318
      case GLP_DB:
319 319
      case GLP_FX:
320 320
        glp_set_col_bnds(lp, i, GLP_LO, lo, up);
321 321
        break;
322 322
      default:
323 323
        break;
324 324
      }
325 325
    } else {
326 326
      switch (b) {
327 327
      case GLP_FR:
328 328
        glp_set_col_bnds(lp, i, GLP_UP, lo, up);
329 329
        break;
330 330
      case GLP_UP:
331 331
        glp_set_col_bnds(lp, i, GLP_UP, lo, up);
332 332
        break;
333 333
      case GLP_LO:
334 334
      case GLP_DB:
335 335
      case GLP_FX:
336 336
        if (lo == up)
337 337
          glp_set_col_bnds(lp, i, GLP_FX, lo, up);
338 338
        else
339 339
          glp_set_col_bnds(lp, i, GLP_DB, lo, up);
340 340
        break;
341 341
      default:
342 342
        break;
343 343
      }
344 344
    }
345 345

	
346 346
  }
347 347

	
348 348
  GlpkBase::Value GlpkBase::_getColUpperBound(int i) const {
349 349
    int b = glp_get_col_type(lp, i);
350 350
      switch (b) {
351 351
      case GLP_UP:
352 352
      case GLP_DB:
353 353
      case GLP_FX:
354 354
        return glp_get_col_ub(lp, i);
355 355
      default:
356 356
        return INF;
357 357
      }
358 358
  }
359 359

	
360 360
  void GlpkBase::_setRowLowerBound(int i, Value lo) {
361 361
    LEMON_ASSERT(lo != INF, "Invalid bound");
362 362

	
363 363
    int b = glp_get_row_type(lp, i);
364 364
    double up = glp_get_row_ub(lp, i);
365 365
    if (lo == -INF) {
366 366
      switch (b) {
367 367
      case GLP_FR:
368 368
      case GLP_LO:
369 369
        glp_set_row_bnds(lp, i, GLP_FR, lo, up);
370 370
        break;
371 371
      case GLP_UP:
372 372
        break;
373 373
      case GLP_DB:
374 374
      case GLP_FX:
375 375
        glp_set_row_bnds(lp, i, GLP_UP, lo, up);
376 376
        break;
377 377
      default:
378 378
        break;
379 379
      }
380 380
    } else {
381 381
      switch (b) {
382 382
      case GLP_FR:
383 383
      case GLP_LO:
384 384
        glp_set_row_bnds(lp, i, GLP_LO, lo, up);
385 385
        break;
386 386
      case GLP_UP:
387 387
      case GLP_DB:
388 388
      case GLP_FX:
389 389
        if (lo == up)
390 390
          glp_set_row_bnds(lp, i, GLP_FX, lo, up);
391 391
        else
392 392
          glp_set_row_bnds(lp, i, GLP_DB, lo, up);
393 393
        break;
394 394
      default:
395 395
        break;
396 396
      }
397 397
    }
398 398

	
399 399
  }
400 400

	
401 401
  GlpkBase::Value GlpkBase::_getRowLowerBound(int i) const {
402 402
    int b = glp_get_row_type(lp, i);
403 403
    switch (b) {
404 404
    case GLP_LO:
405 405
    case GLP_DB:
406 406
    case GLP_FX:
407 407
      return glp_get_row_lb(lp, i);
408 408
    default:
409 409
      return -INF;
410 410
    }
411 411
  }
412 412

	
413 413
  void GlpkBase::_setRowUpperBound(int i, Value up) {
414 414
    LEMON_ASSERT(up != -INF, "Invalid bound");
415 415

	
416 416
    int b = glp_get_row_type(lp, i);
417 417
    double lo = glp_get_row_lb(lp, i);
418 418
    if (up == INF) {
419 419
      switch (b) {
420 420
      case GLP_FR:
421 421
      case GLP_LO:
422 422
        break;
423 423
      case GLP_UP:
424 424
        glp_set_row_bnds(lp, i, GLP_FR, lo, up);
425 425
        break;
426 426
      case GLP_DB:
427 427
      case GLP_FX:
428 428
        glp_set_row_bnds(lp, i, GLP_LO, lo, up);
429 429
        break;
430 430
      default:
431 431
        break;
432 432
      }
433 433
    } else {
434 434
      switch (b) {
435 435
      case GLP_FR:
436 436
        glp_set_row_bnds(lp, i, GLP_UP, lo, up);
437 437
        break;
438 438
      case GLP_UP:
439 439
        glp_set_row_bnds(lp, i, GLP_UP, lo, up);
440 440
        break;
441 441
      case GLP_LO:
442 442
      case GLP_DB:
443 443
      case GLP_FX:
444 444
        if (lo == up)
445 445
          glp_set_row_bnds(lp, i, GLP_FX, lo, up);
446 446
        else
447 447
          glp_set_row_bnds(lp, i, GLP_DB, lo, up);
448 448
        break;
449 449
      default:
450 450
        break;
451 451
      }
452 452
    }
453 453
  }
454 454

	
455 455
  GlpkBase::Value GlpkBase::_getRowUpperBound(int i) const {
456 456
    int b = glp_get_row_type(lp, i);
457 457
    switch (b) {
458 458
    case GLP_UP:
459 459
    case GLP_DB:
460 460
    case GLP_FX:
461 461
      return glp_get_row_ub(lp, i);
462 462
    default:
463 463
      return INF;
464 464
    }
465 465
  }
466 466

	
467 467
  void GlpkBase::_setObjCoeffs(ExprIterator b, ExprIterator e) {
468 468
    for (int i = 1; i <= glp_get_num_cols(lp); ++i) {
469 469
      glp_set_obj_coef(lp, i, 0.0);
470 470
    }
471 471
    for (ExprIterator it = b; it != e; ++it) {
472 472
      glp_set_obj_coef(lp, it->first, it->second);
473 473
    }
474 474
  }
475 475

	
476 476
  void GlpkBase::_getObjCoeffs(InsertIterator b) const {
477 477
    for (int i = 1; i <= glp_get_num_cols(lp); ++i) {
478 478
      Value val = glp_get_obj_coef(lp, i);
479 479
      if (val != 0.0) {
480 480
        *b = std::make_pair(i, val);
481 481
        ++b;
482 482
      }
483 483
    }
484 484
  }
485 485

	
486 486
  void GlpkBase::_setObjCoeff(int i, Value obj_coef) {
487 487
    //i = 0 means the constant term (shift)
488 488
    glp_set_obj_coef(lp, i, obj_coef);
489 489
  }
490 490

	
491 491
  GlpkBase::Value GlpkBase::_getObjCoeff(int i) const {
492 492
    //i = 0 means the constant term (shift)
493 493
    return glp_get_obj_coef(lp, i);
494 494
  }
495 495

	
496 496
  void GlpkBase::_setSense(GlpkBase::Sense sense) {
497 497
    switch (sense) {
498 498
    case MIN:
499 499
      glp_set_obj_dir(lp, GLP_MIN);
500 500
      break;
501 501
    case MAX:
502 502
      glp_set_obj_dir(lp, GLP_MAX);
503 503
      break;
504 504
    }
505 505
  }
506 506

	
507 507
  GlpkBase::Sense GlpkBase::_getSense() const {
508 508
    switch(glp_get_obj_dir(lp)) {
509 509
    case GLP_MIN:
510 510
      return MIN;
511 511
    case GLP_MAX:
512 512
      return MAX;
513 513
    default:
514 514
      LEMON_ASSERT(false, "Wrong sense");
515 515
      return GlpkBase::Sense();
516 516
    }
517 517
  }
518 518

	
519 519
  void GlpkBase::_clear() {
520 520
    glp_erase_prob(lp);
521 521
    rows.clear();
522 522
    cols.clear();
523 523
  }
524 524

	
525 525
  void GlpkBase::freeEnv() {
526 526
    glp_free_env();
527 527
  }
528 528

	
529
  GlpkBase::FreeEnvHelper GlpkBase::freeEnvHelper;
530

	
529 531
  // GlpkLp members
530 532

	
531 533
  GlpkLp::GlpkLp()
532 534
    : LpBase(), GlpkBase(), LpSolver() {
533 535
    messageLevel(MESSAGE_NO_OUTPUT);
534 536
  }
535 537

	
536 538
  GlpkLp::GlpkLp(const GlpkLp& other)
537 539
    : LpBase(other), GlpkBase(other), LpSolver(other) {
538 540
    messageLevel(MESSAGE_NO_OUTPUT);
539 541
  }
540 542

	
541 543
  GlpkLp* GlpkLp::_newSolver() const { return new GlpkLp; }
542 544
  GlpkLp* GlpkLp::_cloneSolver() const { return new GlpkLp(*this); }
543 545

	
544 546
  const char* GlpkLp::_solverName() const { return "GlpkLp"; }
545 547

	
546 548
  void GlpkLp::_clear_temporals() {
547 549
    _primal_ray.clear();
548 550
    _dual_ray.clear();
549 551
  }
550 552

	
551 553
  GlpkLp::SolveExitStatus GlpkLp::_solve() {
552 554
    return solvePrimal();
553 555
  }
554 556

	
555 557
  GlpkLp::SolveExitStatus GlpkLp::solvePrimal() {
556 558
    _clear_temporals();
557 559

	
558 560
    glp_smcp smcp;
559 561
    glp_init_smcp(&smcp);
560 562

	
561 563
    switch (_message_level) {
562 564
    case MESSAGE_NO_OUTPUT:
563 565
      smcp.msg_lev = GLP_MSG_OFF;
564 566
      break;
565 567
    case MESSAGE_ERROR_MESSAGE:
566 568
      smcp.msg_lev = GLP_MSG_ERR;
567 569
      break;
568 570
    case MESSAGE_NORMAL_OUTPUT:
569 571
      smcp.msg_lev = GLP_MSG_ON;
570 572
      break;
571 573
    case MESSAGE_FULL_OUTPUT:
572 574
      smcp.msg_lev = GLP_MSG_ALL;
573 575
      break;
574 576
    }
575 577

	
576 578
    if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
577 579
    return SOLVED;
578 580
  }
579 581

	
580 582
  GlpkLp::SolveExitStatus GlpkLp::solveDual() {
581 583
    _clear_temporals();
582 584

	
583 585
    glp_smcp smcp;
584 586
    glp_init_smcp(&smcp);
585 587

	
586 588
    switch (_message_level) {
587 589
    case MESSAGE_NO_OUTPUT:
588 590
      smcp.msg_lev = GLP_MSG_OFF;
589 591
      break;
590 592
    case MESSAGE_ERROR_MESSAGE:
591 593
      smcp.msg_lev = GLP_MSG_ERR;
592 594
      break;
593 595
    case MESSAGE_NORMAL_OUTPUT:
594 596
      smcp.msg_lev = GLP_MSG_ON;
595 597
      break;
596 598
    case MESSAGE_FULL_OUTPUT:
597 599
      smcp.msg_lev = GLP_MSG_ALL;
598 600
      break;
599 601
    }
600 602
    smcp.meth = GLP_DUAL;
601 603

	
602 604
    if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
603 605
    return SOLVED;
604 606
  }
605 607

	
606 608
  GlpkLp::Value GlpkLp::_getPrimal(int i) const {
607 609
    return glp_get_col_prim(lp, i);
608 610
  }
609 611

	
610 612
  GlpkLp::Value GlpkLp::_getDual(int i) const {
611 613
    return glp_get_row_dual(lp, i);
612 614
  }
613 615

	
614 616
  GlpkLp::Value GlpkLp::_getPrimalValue() const {
615 617
    return glp_get_obj_val(lp);
616 618
  }
617 619

	
618 620
  GlpkLp::VarStatus GlpkLp::_getColStatus(int i) const {
619 621
    switch (glp_get_col_stat(lp, i)) {
620 622
    case GLP_BS:
621 623
      return BASIC;
622 624
    case GLP_UP:
623 625
      return UPPER;
624 626
    case GLP_LO:
625 627
      return LOWER;
626 628
    case GLP_NF:
627 629
      return FREE;
628 630
    case GLP_NS:
629 631
      return FIXED;
630 632
    default:
631 633
      LEMON_ASSERT(false, "Wrong column status");
632 634
      return GlpkLp::VarStatus();
633 635
    }
634 636
  }
635 637

	
636 638
  GlpkLp::VarStatus GlpkLp::_getRowStatus(int i) const {
637 639
    switch (glp_get_row_stat(lp, i)) {
638 640
    case GLP_BS:
639 641
      return BASIC;
640 642
    case GLP_UP:
641 643
      return UPPER;
642 644
    case GLP_LO:
643 645
      return LOWER;
644 646
    case GLP_NF:
645 647
      return FREE;
646 648
    case GLP_NS:
647 649
      return FIXED;
648 650
    default:
649 651
      LEMON_ASSERT(false, "Wrong row status");
650 652
      return GlpkLp::VarStatus();
651 653
    }
652 654
  }
653 655

	
654 656
  GlpkLp::Value GlpkLp::_getPrimalRay(int i) const {
655 657
    if (_primal_ray.empty()) {
656 658
      int row_num = glp_get_num_rows(lp);
657 659
      int col_num = glp_get_num_cols(lp);
658 660

	
659 661
      _primal_ray.resize(col_num + 1, 0.0);
660 662

	
661 663
      int index = glp_get_unbnd_ray(lp);
662 664
      if (index != 0) {
663 665
        // The primal ray is found in primal simplex second phase
664 666
        LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :
665 667
                      glp_get_col_stat(lp, index - row_num)) != GLP_BS,
666 668
                     "Wrong primal ray");
667 669

	
668 670
        bool negate = glp_get_obj_dir(lp) == GLP_MAX;
669 671

	
670 672
        if (index > row_num) {
671 673
          _primal_ray[index - row_num] = 1.0;
672 674
          if (glp_get_col_dual(lp, index - row_num) > 0) {
673 675
            negate = !negate;
674 676
          }
675 677
        } else {
676 678
          if (glp_get_row_dual(lp, index) > 0) {
677 679
            negate = !negate;
678 680
          }
679 681
        }
680 682

	
681 683
        std::vector<int> ray_indexes(row_num + 1);
682 684
        std::vector<Value> ray_values(row_num + 1);
683 685
        int ray_length = glp_eval_tab_col(lp, index, &ray_indexes.front(),
684 686
                                          &ray_values.front());
685 687

	
686 688
        for (int i = 1; i <= ray_length; ++i) {
687 689
          if (ray_indexes[i] > row_num) {
688 690
            _primal_ray[ray_indexes[i] - row_num] = ray_values[i];
689 691
          }
690 692
        }
691 693

	
692 694
        if (negate) {
693 695
          for (int i = 1; i <= col_num; ++i) {
694 696
            _primal_ray[i] = - _primal_ray[i];
695 697
          }
696 698
        }
697 699
      } else {
698 700
        for (int i = 1; i <= col_num; ++i) {
699 701
          _primal_ray[i] = glp_get_col_prim(lp, i);
700 702
        }
701 703
      }
702 704
    }
703 705
    return _primal_ray[i];
704 706
  }
705 707

	
706 708
  GlpkLp::Value GlpkLp::_getDualRay(int i) const {
707 709
    if (_dual_ray.empty()) {
708 710
      int row_num = glp_get_num_rows(lp);
709 711

	
710 712
      _dual_ray.resize(row_num + 1, 0.0);
711 713

	
712 714
      int index = glp_get_unbnd_ray(lp);
713 715
      if (index != 0) {
714 716
        // The dual ray is found in dual simplex second phase
715 717
        LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :
716 718
                      glp_get_col_stat(lp, index - row_num)) == GLP_BS,
717 719

	
718 720
                     "Wrong dual ray");
719 721

	
720 722
        int idx;
721 723
        bool negate = false;
722 724

	
723 725
        if (index > row_num) {
724 726
          idx = glp_get_col_bind(lp, index - row_num);
725 727
          if (glp_get_col_prim(lp, index - row_num) >
726 728
              glp_get_col_ub(lp, index - row_num)) {
727 729
            negate = true;
728 730
          }
729 731
        } else {
730 732
          idx = glp_get_row_bind(lp, index);
731 733
          if (glp_get_row_prim(lp, index) > glp_get_row_ub(lp, index)) {
732 734
            negate = true;
733 735
          }
734 736
        }
735 737

	
736 738
        _dual_ray[idx] = negate ?  - 1.0 : 1.0;
737 739

	
738 740
        glp_btran(lp, &_dual_ray.front());
739 741
      } else {
740 742
        double eps = 1e-7;
741 743
        // The dual ray is found in primal simplex first phase
742 744
        // We assume that the glpk minimizes the slack to get feasible solution
743 745
        for (int i = 1; i <= row_num; ++i) {
744 746
          int index = glp_get_bhead(lp, i);
745 747
          if (index <= row_num) {
746 748
            double res = glp_get_row_prim(lp, index);
747 749
            if (res > glp_get_row_ub(lp, index) + eps) {
748 750
              _dual_ray[i] = -1;
749 751
            } else if (res < glp_get_row_lb(lp, index) - eps) {
750 752
              _dual_ray[i] = 1;
751 753
            } else {
752 754
              _dual_ray[i] = 0;
753 755
            }
754 756
            _dual_ray[i] *= glp_get_rii(lp, index);
755 757
          } else {
756 758
            double res = glp_get_col_prim(lp, index - row_num);
757 759
            if (res > glp_get_col_ub(lp, index - row_num) + eps) {
758 760
              _dual_ray[i] = -1;
759 761
            } else if (res < glp_get_col_lb(lp, index - row_num) - eps) {
760 762
              _dual_ray[i] = 1;
761 763
            } else {
762 764
              _dual_ray[i] = 0;
763 765
            }
764 766
            _dual_ray[i] /= glp_get_sjj(lp, index - row_num);
765 767
          }
766 768
        }
767 769

	
768 770
        glp_btran(lp, &_dual_ray.front());
769 771

	
770 772
        for (int i = 1; i <= row_num; ++i) {
771 773
          _dual_ray[i] /= glp_get_rii(lp, i);
772 774
        }
773 775
      }
774 776
    }
775 777
    return _dual_ray[i];
776 778
  }
777 779

	
778 780
  GlpkLp::ProblemType GlpkLp::_getPrimalType() const {
779 781
    if (glp_get_status(lp) == GLP_OPT)
780 782
      return OPTIMAL;
781 783
    switch (glp_get_prim_stat(lp)) {
782 784
    case GLP_UNDEF:
783 785
      return UNDEFINED;
784 786
    case GLP_FEAS:
785 787
    case GLP_INFEAS:
786 788
      if (glp_get_dual_stat(lp) == GLP_NOFEAS) {
787 789
        return UNBOUNDED;
788 790
      } else {
789 791
        return UNDEFINED;
790 792
      }
791 793
    case GLP_NOFEAS:
792 794
      return INFEASIBLE;
793 795
    default:
794 796
      LEMON_ASSERT(false, "Wrong primal type");
795 797
      return  GlpkLp::ProblemType();
796 798
    }
797 799
  }
798 800

	
799 801
  GlpkLp::ProblemType GlpkLp::_getDualType() const {
800 802
    if (glp_get_status(lp) == GLP_OPT)
801 803
      return OPTIMAL;
802 804
    switch (glp_get_dual_stat(lp)) {
803 805
    case GLP_UNDEF:
804 806
      return UNDEFINED;
805 807
    case GLP_FEAS:
806 808
    case GLP_INFEAS:
807 809
      if (glp_get_prim_stat(lp) == GLP_NOFEAS) {
808 810
        return UNBOUNDED;
809 811
      } else {
810 812
        return UNDEFINED;
811 813
      }
812 814
    case GLP_NOFEAS:
813 815
      return INFEASIBLE;
814 816
    default:
815 817
      LEMON_ASSERT(false, "Wrong primal type");
816 818
      return  GlpkLp::ProblemType();
817 819
    }
818 820
  }
819 821

	
820 822
  void GlpkLp::presolver(bool b) {
821 823
    lpx_set_int_parm(lp, LPX_K_PRESOL, b ? 1 : 0);
822 824
  }
823 825

	
824 826
  void GlpkLp::messageLevel(MessageLevel m) {
825 827
    _message_level = m;
826 828
  }
827 829

	
828 830
  // GlpkMip members
829 831

	
830 832
  GlpkMip::GlpkMip()
831 833
    : LpBase(), GlpkBase(), MipSolver() {
832 834
    messageLevel(MESSAGE_NO_OUTPUT);
833 835
  }
834 836

	
835 837
  GlpkMip::GlpkMip(const GlpkMip& other)
836 838
    : LpBase(), GlpkBase(other), MipSolver() {
837 839
    messageLevel(MESSAGE_NO_OUTPUT);
838 840
  }
839 841

	
840 842
  void GlpkMip::_setColType(int i, GlpkMip::ColTypes col_type) {
841 843
    switch (col_type) {
842 844
    case INTEGER:
843 845
      glp_set_col_kind(lp, i, GLP_IV);
844 846
      break;
845 847
    case REAL:
846 848
      glp_set_col_kind(lp, i, GLP_CV);
847 849
      break;
848 850
    }
849 851
  }
850 852

	
851 853
  GlpkMip::ColTypes GlpkMip::_getColType(int i) const {
852 854
    switch (glp_get_col_kind(lp, i)) {
853 855
    case GLP_IV:
854 856
    case GLP_BV:
855 857
      return INTEGER;
856 858
    default:
857 859
      return REAL;
858 860
    }
859 861

	
860 862
  }
861 863

	
862 864
  GlpkMip::SolveExitStatus GlpkMip::_solve() {
863 865
    glp_smcp smcp;
864 866
    glp_init_smcp(&smcp);
865 867

	
866 868
    switch (_message_level) {
867 869
    case MESSAGE_NO_OUTPUT:
868 870
      smcp.msg_lev = GLP_MSG_OFF;
869 871
      break;
870 872
    case MESSAGE_ERROR_MESSAGE:
871 873
      smcp.msg_lev = GLP_MSG_ERR;
872 874
      break;
873 875
    case MESSAGE_NORMAL_OUTPUT:
874 876
      smcp.msg_lev = GLP_MSG_ON;
875 877
      break;
876 878
    case MESSAGE_FULL_OUTPUT:
877 879
      smcp.msg_lev = GLP_MSG_ALL;
878 880
      break;
879 881
    }
880 882
    smcp.meth = GLP_DUAL;
881 883

	
882 884
    if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
883 885
    if (glp_get_status(lp) != GLP_OPT) return SOLVED;
884 886

	
885 887
    glp_iocp iocp;
886 888
    glp_init_iocp(&iocp);
887 889

	
888 890
    switch (_message_level) {
889 891
    case MESSAGE_NO_OUTPUT:
890 892
      iocp.msg_lev = GLP_MSG_OFF;
891 893
      break;
892 894
    case MESSAGE_ERROR_MESSAGE:
893 895
      iocp.msg_lev = GLP_MSG_ERR;
894 896
      break;
895 897
    case MESSAGE_NORMAL_OUTPUT:
896 898
      iocp.msg_lev = GLP_MSG_ON;
897 899
      break;
898 900
    case MESSAGE_FULL_OUTPUT:
899 901
      iocp.msg_lev = GLP_MSG_ALL;
900 902
      break;
901 903
    }
902 904

	
903 905
    if (glp_intopt(lp, &iocp) != 0) return UNSOLVED;
904 906
    return SOLVED;
905 907
  }
906 908

	
907 909

	
908 910
  GlpkMip::ProblemType GlpkMip::_getType() const {
909 911
    switch (glp_get_status(lp)) {
910 912
    case GLP_OPT:
911 913
      switch (glp_mip_status(lp)) {
912 914
      case GLP_UNDEF:
913 915
        return UNDEFINED;
914 916
      case GLP_NOFEAS:
915 917
        return INFEASIBLE;
916 918
      case GLP_FEAS:
917 919
        return FEASIBLE;
918 920
      case GLP_OPT:
919 921
        return OPTIMAL;
920 922
      default:
921 923
        LEMON_ASSERT(false, "Wrong problem type.");
922 924
        return GlpkMip::ProblemType();
923 925
      }
924 926
    case GLP_NOFEAS:
925 927
      return INFEASIBLE;
926 928
    case GLP_INFEAS:
927 929
    case GLP_FEAS:
928 930
      if (glp_get_dual_stat(lp) == GLP_NOFEAS) {
929 931
        return UNBOUNDED;
930 932
      } else {
931 933
        return UNDEFINED;
932 934
      }
933 935
    default:
934 936
      LEMON_ASSERT(false, "Wrong problem type.");
935 937
      return GlpkMip::ProblemType();
936 938
    }
937 939
  }
938 940

	
939 941
  GlpkMip::Value GlpkMip::_getSol(int i) const {
940 942
    return glp_mip_col_val(lp, i);
941 943
  }
942 944

	
943 945
  GlpkMip::Value GlpkMip::_getSolValue() const {
944 946
    return glp_mip_obj_val(lp);
945 947
  }
946 948

	
947 949
  GlpkMip* GlpkMip::_newSolver() const { return new GlpkMip; }
948 950
  GlpkMip* GlpkMip::_cloneSolver() const {return new GlpkMip(*this); }
949 951

	
950 952
  const char* GlpkMip::_solverName() const { return "GlpkMip"; }
951 953

	
952 954
  void GlpkMip::messageLevel(MessageLevel m) {
953 955
    _message_level = m;
954 956
  }
955 957

	
956 958
} //END OF NAMESPACE LEMON
Ignore white space 6 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2008
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#ifndef LEMON_GLPK_H
20 20
#define LEMON_GLPK_H
21 21

	
22 22
///\file
23 23
///\brief Header of the LEMON-GLPK lp solver interface.
24 24
///\ingroup lp_group
25 25

	
26 26
#include <lemon/lp_base.h>
27 27

	
28 28
// forward declaration
29 29
#ifndef _GLP_PROB
30 30
#define _GLP_PROB
31 31
typedef struct { double _prob; } glp_prob;
32 32
/* LP/MIP problem object */
33 33
#endif
34 34

	
35 35
namespace lemon {
36 36

	
37 37

	
38 38
  /// \brief Base interface for the GLPK LP and MIP solver
39 39
  ///
40 40
  /// This class implements the common interface of the GLPK LP and MIP solver.
41 41
  /// \ingroup lp_group
42 42
  class GlpkBase : virtual public LpBase {
43 43
  protected:
44 44

	
45 45
    typedef glp_prob LPX;
46 46
    glp_prob* lp;
47 47

	
48 48
    GlpkBase();
49 49
    GlpkBase(const GlpkBase&);
50 50
    virtual ~GlpkBase();
51 51

	
52 52
  protected:
53 53

	
54 54
    virtual int _addCol();
55 55
    virtual int _addRow();
56 56

	
57 57
    virtual void _eraseCol(int i);
58 58
    virtual void _eraseRow(int i);
59 59

	
60 60
    virtual void _eraseColId(int i);
61 61
    virtual void _eraseRowId(int i);
62 62

	
63 63
    virtual void _getColName(int col, std::string& name) const;
64 64
    virtual void _setColName(int col, const std::string& name);
65 65
    virtual int _colByName(const std::string& name) const;
66 66

	
67 67
    virtual void _getRowName(int row, std::string& name) const;
68 68
    virtual void _setRowName(int row, const std::string& name);
69 69
    virtual int _rowByName(const std::string& name) const;
70 70

	
71 71
    virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e);
72 72
    virtual void _getRowCoeffs(int i, InsertIterator b) const;
73 73

	
74 74
    virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e);
75 75
    virtual void _getColCoeffs(int i, InsertIterator b) const;
76 76

	
77 77
    virtual void _setCoeff(int row, int col, Value value);
78 78
    virtual Value _getCoeff(int row, int col) const;
79 79

	
80 80
    virtual void _setColLowerBound(int i, Value value);
81 81
    virtual Value _getColLowerBound(int i) const;
82 82

	
83 83
    virtual void _setColUpperBound(int i, Value value);
84 84
    virtual Value _getColUpperBound(int i) const;
85 85

	
86 86
    virtual void _setRowLowerBound(int i, Value value);
87 87
    virtual Value _getRowLowerBound(int i) const;
88 88

	
89 89
    virtual void _setRowUpperBound(int i, Value value);
90 90
    virtual Value _getRowUpperBound(int i) const;
91 91

	
92 92
    virtual void _setObjCoeffs(ExprIterator b, ExprIterator e);
93 93
    virtual void _getObjCoeffs(InsertIterator b) const;
94 94

	
95 95
    virtual void _setObjCoeff(int i, Value obj_coef);
96 96
    virtual Value _getObjCoeff(int i) const;
97 97

	
98 98
    virtual void _setSense(Sense);
99 99
    virtual Sense _getSense() const;
100 100

	
101 101
    virtual void _clear();
102 102

	
103
  private:
104

	
105
    static void freeEnv();
106

	
107
    struct FreeEnvHelper {
108
      ~FreeEnvHelper() {
109
        freeEnv();
110
      }
111
    };
112
    
113
    static FreeEnvHelper freeEnvHelper;
114
    
103 115
  public:
104 116

	
105
    /// \brief Deallocates the globally allocated memory of GLPK.
106

	
107
    /// Deallocates the globally allocated memory of GLPK.  \note
108
    /// Usually, it do not have to be called, because the GLPK use
109
    /// only a small amount of global memory, and it is deallocated
110
    /// automatically at the end of program.
111
    static void freeEnv();
112

	
113 117
    ///Pointer to the underlying GLPK data structure.
114 118
    LPX *lpx() {return lp;}
115 119
    ///Const pointer to the underlying GLPK data structure.
116 120
    const LPX *lpx() const {return lp;}
117 121

	
118 122
    ///Returns the constraint identifier understood by GLPK.
119 123
    int lpxRow(Row r) const { return rows(id(r)); }
120 124

	
121 125
    ///Returns the variable identifier understood by GLPK.
122 126
    int lpxCol(Col c) const { return cols(id(c)); }
123 127

	
124 128
  };
125 129

	
126 130
  /// \brief Interface for the GLPK LP solver
127 131
  ///
128 132
  /// This class implements an interface for the GLPK LP solver.
129 133
  ///\ingroup lp_group
130 134
  class GlpkLp : public GlpkBase, public LpSolver {
131 135
  public:
132 136

	
133 137
    ///\e
134 138
    GlpkLp();
135 139
    ///\e
136 140
    GlpkLp(const GlpkLp&);
137 141

	
138 142
  private:
139 143

	
140 144
    mutable std::vector<double> _primal_ray;
141 145
    mutable std::vector<double> _dual_ray;
142 146

	
143 147
    void _clear_temporals();
144 148

	
145 149
  protected:
146 150

	
147 151
    virtual GlpkLp* _cloneSolver() const;
148 152
    virtual GlpkLp* _newSolver() const;
149 153

	
150 154
    virtual const char* _solverName() const;
151 155

	
152 156
    virtual SolveExitStatus _solve();
153 157
    virtual Value _getPrimal(int i) const;
154 158
    virtual Value _getDual(int i) const;
155 159

	
156 160
    virtual Value _getPrimalValue() const;
157 161

	
158 162
    virtual VarStatus _getColStatus(int i) const;
159 163
    virtual VarStatus _getRowStatus(int i) const;
160 164

	
161 165
    virtual Value _getPrimalRay(int i) const;
162 166
    virtual Value _getDualRay(int i) const;
163 167

	
164 168
    ///\todo It should be clarified
165 169
    ///
166 170
    virtual ProblemType _getPrimalType() const;
167 171
    virtual ProblemType _getDualType() const;
168 172

	
169 173
  public:
170 174

	
171 175
    ///Solve with primal simplex
172 176
    SolveExitStatus solvePrimal();
173 177

	
174 178
    ///Solve with dual simplex
175 179
    SolveExitStatus solveDual();
176 180

	
177 181
    ///Turns on or off the presolver
178 182

	
179 183
    ///Turns on (\c b is \c true) or off (\c b is \c false) the presolver
180 184
    ///
181 185
    ///The presolver is off by default.
182 186
    void presolver(bool b);
183 187

	
184 188
    ///Enum for \c messageLevel() parameter
185 189
    enum MessageLevel {
186 190
      /// no output (default value)
187 191
      MESSAGE_NO_OUTPUT = 0,
188 192
      /// error messages only
189 193
      MESSAGE_ERROR_MESSAGE = 1,
190 194
      /// normal output
191 195
      MESSAGE_NORMAL_OUTPUT = 2,
192 196
      /// full output (includes informational messages)
193 197
      MESSAGE_FULL_OUTPUT = 3
194 198
    };
195 199

	
196 200
  private:
197 201

	
198 202
    MessageLevel _message_level;
199 203

	
200 204
  public:
201 205

	
202 206
    ///Set the verbosity of the messages
203 207

	
204 208
    ///Set the verbosity of the messages
205 209
    ///
206 210
    ///\param m is the level of the messages output by the solver routines.
207 211
    void messageLevel(MessageLevel m);
208 212
  };
209 213

	
210 214
  /// \brief Interface for the GLPK MIP solver
211 215
  ///
212 216
  /// This class implements an interface for the GLPK MIP solver.
213 217
  ///\ingroup lp_group
214 218
  class GlpkMip : public GlpkBase, public MipSolver {
215 219
  public:
216 220

	
217 221
    ///\e
218 222
    GlpkMip();
219 223
    ///\e
220 224
    GlpkMip(const GlpkMip&);
221 225

	
222 226
  protected:
223 227

	
224 228
    virtual GlpkMip* _cloneSolver() const;
225 229
    virtual GlpkMip* _newSolver() const;
226 230

	
227 231
    virtual const char* _solverName() const;
228 232

	
229 233
    virtual ColTypes _getColType(int col) const;
230 234
    virtual void _setColType(int col, ColTypes col_type);
231 235

	
232 236
    virtual SolveExitStatus _solve();
233 237
    virtual ProblemType _getType() const;
234 238
    virtual Value _getSol(int i) const;
235 239
    virtual Value _getSolValue() const;
236 240

	
237 241
    ///Enum for \c messageLevel() parameter
238 242
    enum MessageLevel {
239 243
      /// no output (default value)
240 244
      MESSAGE_NO_OUTPUT = 0,
241 245
      /// error messages only
242 246
      MESSAGE_ERROR_MESSAGE = 1,
243 247
      /// normal output
244 248
      MESSAGE_NORMAL_OUTPUT = 2,
245 249
      /// full output (includes informational messages)
246 250
      MESSAGE_FULL_OUTPUT = 3
247 251
    };
248 252

	
249 253
  private:
250 254

	
251 255
    MessageLevel _message_level;
252 256

	
253 257
  public:
254 258

	
255 259
    ///Set the verbosity of the messages
256 260

	
257 261
    ///Set the verbosity of the messages
258 262
    ///
259 263
    ///\param m is the level of the messages output by the solver routines.
260 264
    void messageLevel(MessageLevel m);
261 265
  };
262 266

	
263 267

	
264 268
} //END OF NAMESPACE LEMON
265 269

	
266 270
#endif //LEMON_GLPK_H
267 271

	
Ignore white space 49152 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2008
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#include <sstream>
20 20
#include <lemon/lp_skeleton.h>
21 21
#include "test_tools.h"
22 22
#include <lemon/tolerance.h>
23 23

	
24 24
#ifdef HAVE_CONFIG_H
25 25
#include <lemon/config.h>
26 26
#endif
27 27

	
28 28
#ifdef HAVE_GLPK
29 29
#include <lemon/glpk.h>
30 30
#endif
31 31

	
32 32
#ifdef HAVE_CPLEX
33 33
#include <lemon/cplex.h>
34 34
#endif
35 35

	
36 36
#ifdef HAVE_SOPLEX
37 37
#include <lemon/soplex.h>
38 38
#endif
39 39

	
40 40
#ifdef HAVE_CLP
41 41
#include <lemon/clp.h>
42 42
#endif
43 43

	
44 44
using namespace lemon;
45 45

	
46 46
void lpTest(LpSolver& lp)
47 47
{
48 48

	
49 49
  typedef LpSolver LP;
50 50

	
51 51
  std::vector<LP::Col> x(10);
52 52
  //  for(int i=0;i<10;i++) x.push_back(lp.addCol());
53 53
  lp.addColSet(x);
54 54
  lp.colLowerBound(x,1);
55 55
  lp.colUpperBound(x,1);
56 56
  lp.colBounds(x,1,2);
57 57

	
58 58
  std::vector<LP::Col> y(10);
59 59
  lp.addColSet(y);
60 60

	
61 61
  lp.colLowerBound(y,1);
62 62
  lp.colUpperBound(y,1);
63 63
  lp.colBounds(y,1,2);
64 64

	
65 65
  std::map<int,LP::Col> z;
66 66

	
67 67
  z.insert(std::make_pair(12,INVALID));
68 68
  z.insert(std::make_pair(2,INVALID));
69 69
  z.insert(std::make_pair(7,INVALID));
70 70
  z.insert(std::make_pair(5,INVALID));
71 71

	
72 72
  lp.addColSet(z);
73 73

	
74 74
  lp.colLowerBound(z,1);
75 75
  lp.colUpperBound(z,1);
76 76
  lp.colBounds(z,1,2);
77 77

	
78 78
  {
79 79
    LP::Expr e,f,g;
80 80
    LP::Col p1,p2,p3,p4,p5;
81 81
    LP::Constr c;
82 82

	
83 83
    p1=lp.addCol();
84 84
    p2=lp.addCol();
85 85
    p3=lp.addCol();
86 86
    p4=lp.addCol();
87 87
    p5=lp.addCol();
88 88

	
89 89
    e[p1]=2;
90 90
    *e=12;
91 91
    e[p1]+=2;
92 92
    *e+=12;
93 93
    e[p1]-=2;
94 94
    *e-=12;
95 95

	
96 96
    e=2;
97 97
    e=2.2;
98 98
    e=p1;
99 99
    e=f;
100 100

	
101 101
    e+=2;
102 102
    e+=2.2;
103 103
    e+=p1;
104 104
    e+=f;
105 105

	
106 106
    e-=2;
107 107
    e-=2.2;
108 108
    e-=p1;
109 109
    e-=f;
110 110

	
111 111
    e*=2;
112 112
    e*=2.2;
113 113
    e/=2;
114 114
    e/=2.2;
115 115

	
116 116
    e=((p1+p2)+(p1-p2)+(p1+12)+(12+p1)+(p1-12)+(12-p1)+
117 117
       (f+12)+(12+f)+(p1+f)+(f+p1)+(f+g)+
118 118
       (f-12)+(12-f)+(p1-f)+(f-p1)+(f-g)+
119 119
       2.2*f+f*2.2+f/2.2+
120 120
       2*f+f*2+f/2+
121 121
       2.2*p1+p1*2.2+p1/2.2+
122 122
       2*p1+p1*2+p1/2
123 123
       );
124 124

	
125 125

	
126 126
    c = (e  <= f  );
127 127
    c = (e  <= 2.2);
128 128
    c = (e  <= 2  );
129 129
    c = (e  <= p1 );
130 130
    c = (2.2<= f  );
131 131
    c = (2  <= f  );
132 132
    c = (p1 <= f  );
133 133
    c = (p1 <= p2 );
134 134
    c = (p1 <= 2.2);
135 135
    c = (p1 <= 2  );
136 136
    c = (2.2<= p2 );
137 137
    c = (2  <= p2 );
138 138

	
139 139
    c = (e  >= f  );
140 140
    c = (e  >= 2.2);
141 141
    c = (e  >= 2  );
142 142
    c = (e  >= p1 );
143 143
    c = (2.2>= f  );
144 144
    c = (2  >= f  );
145 145
    c = (p1 >= f  );
146 146
    c = (p1 >= p2 );
147 147
    c = (p1 >= 2.2);
148 148
    c = (p1 >= 2  );
149 149
    c = (2.2>= p2 );
150 150
    c = (2  >= p2 );
151 151

	
152 152
    c = (e  == f  );
153 153
    c = (e  == 2.2);
154 154
    c = (e  == 2  );
155 155
    c = (e  == p1 );
156 156
    c = (2.2== f  );
157 157
    c = (2  == f  );
158 158
    c = (p1 == f  );
159 159
    //c = (p1 == p2 );
160 160
    c = (p1 == 2.2);
161 161
    c = (p1 == 2  );
162 162
    c = (2.2== p2 );
163 163
    c = (2  == p2 );
164 164

	
165 165
    c = ((2 <= e) <= 3);
166 166
    c = ((2 <= p1) <= 3);
167 167

	
168 168
    c = ((2 >= e) >= 3);
169 169
    c = ((2 >= p1) >= 3);
170 170

	
171 171
    e[x[3]]=2;
172 172
    e[x[3]]=4;
173 173
    e[x[3]]=1;
174 174
    *e=12;
175 175

	
176 176
    lp.addRow(-LP::INF,e,23);
177 177
    lp.addRow(-LP::INF,3.0*(x[1]+x[2]/2)-x[3],23);
178 178
    lp.addRow(-LP::INF,3.0*(x[1]+x[2]*2-5*x[3]+12-x[4]/3)+2*x[4]-4,23);
179 179

	
180 180
    lp.addRow(x[1]+x[3]<=x[5]-3);
181 181
    lp.addRow((-7<=x[1]+x[3]-12)<=3);
182 182
    lp.addRow(x[1]<=x[5]);
183 183

	
184 184
    std::ostringstream buf;
185 185

	
186 186

	
187 187
    e=((p1+p2)+(p1-0.99*p2));
188 188
    //e.prettyPrint(std::cout);
189 189
    //(e<=2).prettyPrint(std::cout);
190 190
    double tolerance=0.001;
191 191
    e.simplify(tolerance);
192 192
    buf << "Coeff. of p2 should be 0.01";
193 193
    check(e[p2]>0, buf.str());
194 194

	
195 195
    tolerance=0.02;
196 196
    e.simplify(tolerance);
197 197
    buf << "Coeff. of p2 should be 0";
198 198
    check(const_cast<const LpSolver::Expr&>(e)[p2]==0, buf.str());
199 199

	
200 200

	
201 201
  }
202 202

	
203 203
  {
204 204
    LP::DualExpr e,f,g;
205 205
    LP::Row p1 = INVALID, p2 = INVALID, p3 = INVALID,
206 206
      p4 = INVALID, p5 = INVALID;
207 207

	
208 208
    e[p1]=2;
209 209
    e[p1]+=2;
210 210
    e[p1]-=2;
211 211

	
212 212
    e=p1;
213 213
    e=f;
214 214

	
215 215
    e+=p1;
216 216
    e+=f;
217 217

	
218 218
    e-=p1;
219 219
    e-=f;
220 220

	
221 221
    e*=2;
222 222
    e*=2.2;
223 223
    e/=2;
224 224
    e/=2.2;
225 225

	
226 226
    e=((p1+p2)+(p1-p2)+
227 227
       (p1+f)+(f+p1)+(f+g)+
228 228
       (p1-f)+(f-p1)+(f-g)+
229 229
       2.2*f+f*2.2+f/2.2+
230 230
       2*f+f*2+f/2+
231 231
       2.2*p1+p1*2.2+p1/2.2+
232 232
       2*p1+p1*2+p1/2
233 233
       );
234 234
  }
235 235

	
236 236
}
237 237

	
238 238
void solveAndCheck(LpSolver& lp, LpSolver::ProblemType stat,
239 239
                   double exp_opt) {
240 240
  using std::string;
241 241
  lp.solve();
242 242

	
243 243
  std::ostringstream buf;
244 244
  buf << "PrimalType should be: " << int(stat) << int(lp.primalType());
245 245

	
246 246
  check(lp.primalType()==stat, buf.str());
247 247

	
248 248
  if (stat ==  LpSolver::OPTIMAL) {
249 249
    std::ostringstream sbuf;
250 250
    sbuf << "Wrong optimal value: the right optimum is " << exp_opt;
251 251
    check(std::abs(lp.primal()-exp_opt) < 1e-3, sbuf.str());
252 252
  }
253 253
}
254 254

	
255 255
void aTest(LpSolver & lp)
256 256
{
257 257
  typedef LpSolver LP;
258 258

	
259 259
 //The following example is very simple
260 260

	
261 261
  typedef LpSolver::Row Row;
262 262
  typedef LpSolver::Col Col;
263 263

	
264 264

	
265 265
  Col x1 = lp.addCol();
266 266
  Col x2 = lp.addCol();
267 267

	
268 268

	
269 269
  //Constraints
270 270
  Row upright=lp.addRow(x1+2*x2 <=1);
271 271
  lp.addRow(x1+x2 >=-1);
272 272
  lp.addRow(x1-x2 <=1);
273 273
  lp.addRow(x1-x2 >=-1);
274 274
  //Nonnegativity of the variables
275 275
  lp.colLowerBound(x1, 0);
276 276
  lp.colLowerBound(x2, 0);
277 277
  //Objective function
278 278
  lp.obj(x1+x2);
279 279

	
280 280
  lp.sense(lp.MAX);
281 281

	
282 282
  //Testing the problem retrieving routines
283 283
  check(lp.objCoeff(x1)==1,"First term should be 1 in the obj function!");
284 284
  check(lp.sense() == lp.MAX,"This is a maximization!");
285 285
  check(lp.coeff(upright,x1)==1,"The coefficient in question is 1!");
286 286
  check(lp.colLowerBound(x1)==0,
287 287
        "The lower bound for variable x1 should be 0.");
288 288
  check(lp.colUpperBound(x1)==LpSolver::INF,
289 289
        "The upper bound for variable x1 should be infty.");
290 290
  check(lp.rowLowerBound(upright) == -LpSolver::INF,
291 291
        "The lower bound for the first row should be -infty.");
292 292
  check(lp.rowUpperBound(upright)==1,
293 293
        "The upper bound for the first row should be 1.");
294 294
  LpSolver::Expr e = lp.row(upright);
295 295
  check(e[x1] == 1, "The first coefficient should 1.");
296 296
  check(e[x2] == 2, "The second coefficient should 1.");
297 297

	
298 298
  lp.row(upright, x1+x2 <=1);
299 299
  e = lp.row(upright);
300 300
  check(e[x1] == 1, "The first coefficient should 1.");
301 301
  check(e[x2] == 1, "The second coefficient should 1.");
302 302

	
303 303
  LpSolver::DualExpr de = lp.col(x1);
304 304
  check(  de[upright] == 1, "The first coefficient should 1.");
305 305

	
306 306
  LpSolver* clp = lp.cloneSolver();
307 307

	
308 308
  //Testing the problem retrieving routines
309 309
  check(clp->objCoeff(x1)==1,"First term should be 1 in the obj function!");
310 310
  check(clp->sense() == clp->MAX,"This is a maximization!");
311 311
  check(clp->coeff(upright,x1)==1,"The coefficient in question is 1!");
312 312
  //  std::cout<<lp.colLowerBound(x1)<<std::endl;
313 313
  check(clp->colLowerBound(x1)==0,
314 314
        "The lower bound for variable x1 should be 0.");
315 315
  check(clp->colUpperBound(x1)==LpSolver::INF,
316 316
        "The upper bound for variable x1 should be infty.");
317 317

	
318 318
  check(lp.rowLowerBound(upright)==-LpSolver::INF,
319 319
        "The lower bound for the first row should be -infty.");
320 320
  check(lp.rowUpperBound(upright)==1,
321 321
        "The upper bound for the first row should be 1.");
322 322
  e = clp->row(upright);
323 323
  check(e[x1] == 1, "The first coefficient should 1.");
324 324
  check(e[x2] == 1, "The second coefficient should 1.");
325 325

	
326 326
  de = clp->col(x1);
327 327
  check(de[upright] == 1, "The first coefficient should 1.");
328 328

	
329 329
  delete clp;
330 330

	
331 331
  //Maximization of x1+x2
332 332
  //over the triangle with vertices (0,0) (0,1) (1,0)
333 333
  double expected_opt=1;
334 334
  solveAndCheck(lp, LpSolver::OPTIMAL, expected_opt);
335 335

	
336 336
  //Minimization
337 337
  lp.sense(lp.MIN);
338 338
  expected_opt=0;
339 339
  solveAndCheck(lp, LpSolver::OPTIMAL, expected_opt);
340 340

	
341 341
  //Vertex (-1,0) instead of (0,0)
342 342
  lp.colLowerBound(x1, -LpSolver::INF);
343 343
  expected_opt=-1;
344 344
  solveAndCheck(lp, LpSolver::OPTIMAL, expected_opt);
345 345

	
346 346
  //Erase one constraint and return to maximization
347 347
  lp.erase(upright);
348 348
  lp.sense(lp.MAX);
349 349
  expected_opt=LpSolver::INF;
350 350
  solveAndCheck(lp, LpSolver::UNBOUNDED, expected_opt);
351 351

	
352 352
  //Infeasibilty
353 353
  lp.addRow(x1+x2 <=-2);
354 354
  solveAndCheck(lp, LpSolver::INFEASIBLE, expected_opt);
355 355

	
356 356
}
357 357

	
358 358
int main()
359 359
{
360 360
  LpSkeleton lp_skel;
361 361
  lpTest(lp_skel);
362 362

	
363 363
#ifdef HAVE_GLPK
364 364
  {
365 365
    GlpkLp lp_glpk1,lp_glpk2;
366 366
    lpTest(lp_glpk1);
367 367
    aTest(lp_glpk2);
368 368
  }
369
  GlpkLp::freeEnv();
370 369
#endif
371 370

	
372 371
#ifdef HAVE_CPLEX
373 372
  try {
374 373
    CplexLp lp_cplex1,lp_cplex2;
375 374
    lpTest(lp_cplex1);
376 375
    aTest(lp_cplex2);
377 376
  } catch (CplexEnv::LicenseError& error) {
378 377
#ifdef LEMON_FORCE_CPLEX_CHECK
379 378
    check(false, error.what());
380 379
#else
381 380
    std::cerr << error.what() << std::endl;
382 381
    std::cerr << "Cplex license check failed, lp check skipped" << std::endl;
383 382
#endif
384 383
  }
385 384
#endif
386 385

	
387 386
#ifdef HAVE_SOPLEX
388 387
  {
389 388
    SoplexLp lp_soplex1,lp_soplex2;
390 389
    lpTest(lp_soplex1);
391 390
    aTest(lp_soplex2);
392 391
  }
393 392
#endif
394 393

	
395 394
#ifdef HAVE_CLP
396 395
  {
397 396
    ClpLp lp_clp1,lp_clp2;
398 397
    lpTest(lp_clp1);
399 398
    aTest(lp_clp2);
400 399
  }
401 400
#endif
402 401

	
403 402
  return 0;
404 403
}
Ignore white space 6 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2008
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#include "test_tools.h"
20 20

	
21 21

	
22 22
#ifdef HAVE_CONFIG_H
23 23
#include <lemon/config.h>
24 24
#endif
25 25

	
26 26
#ifdef HAVE_CPLEX
27 27
#include <lemon/cplex.h>
28 28
#endif
29 29

	
30 30
#ifdef HAVE_GLPK
31 31
#include <lemon/glpk.h>
32 32
#endif
33 33

	
34 34

	
35 35
using namespace lemon;
36 36

	
37 37
void solveAndCheck(MipSolver& mip, MipSolver::ProblemType stat,
38 38
                   double exp_opt) {
39 39
  using std::string;
40 40

	
41 41
  mip.solve();
42 42
  //int decimal,sign;
43 43
  std::ostringstream buf;
44 44
  buf << "Type should be: " << int(stat)<<" and it is "<<int(mip.type());
45 45

	
46 46

	
47 47
  //  itoa(stat,buf1, 10);
48 48
  check(mip.type()==stat, buf.str());
49 49

	
50 50
  if (stat ==  MipSolver::OPTIMAL) {
51 51
    std::ostringstream sbuf;
52 52
    buf << "Wrong optimal value: the right optimum is " << exp_opt;
53 53
    check(std::abs(mip.solValue()-exp_opt) < 1e-3, sbuf.str());
54 54
    //+ecvt(exp_opt,2)
55 55
  }
56 56
}
57 57

	
58 58
void aTest(MipSolver& mip)
59 59
{
60 60
 //The following example is very simple
61 61

	
62 62

	
63 63
  typedef MipSolver::Row Row;
64 64
  typedef MipSolver::Col Col;
65 65

	
66 66

	
67 67

	
68 68
  Col x1 = mip.addCol();
69 69
  Col x2 = mip.addCol();
70 70

	
71 71

	
72 72
  //Objective function
73 73
  mip.obj(x1);
74 74

	
75 75
  mip.max();
76 76

	
77 77

	
78 78
  //Unconstrained optimization
79 79
  mip.solve();
80 80
  //Check it out!
81 81

	
82 82
  //Constraints
83 83
  mip.addRow(2*x1+x2 <=2);
84 84
  mip.addRow(x1-2*x2 <=0);
85 85

	
86 86
  //Nonnegativity of the variable x1
87 87
  mip.colLowerBound(x1, 0);
88 88

	
89 89
  //Maximization of x1
90 90
  //over the triangle with vertices (0,0),(4/5,2/5),(0,2)
91 91
  double expected_opt=4.0/5.0;
92 92
  solveAndCheck(mip, MipSolver::OPTIMAL, expected_opt);
93 93

	
94 94
  //Restrict x2 to integer
95 95
  mip.colType(x2,MipSolver::INTEGER);
96 96
  expected_opt=1.0/2.0;
97 97
  solveAndCheck(mip, MipSolver::OPTIMAL, expected_opt);
98 98

	
99 99

	
100 100
  //Restrict both to integer
101 101
  mip.colType(x1,MipSolver::INTEGER);
102 102
  expected_opt=0;
103 103
  solveAndCheck(mip, MipSolver::OPTIMAL, expected_opt);
104 104

	
105 105

	
106 106

	
107 107
}
108 108

	
109 109

	
110 110
int main()
111 111
{
112 112

	
113 113
#ifdef HAVE_GLPK
114 114
  {
115 115
    GlpkMip mip1;
116 116
    aTest(mip1);
117 117
  }
118
  GlpkLp::freeEnv();
119 118
#endif
120 119

	
121 120
#ifdef HAVE_CPLEX
122 121
  try {
123 122
    CplexMip mip2;
124 123
    aTest(mip2);
125 124
  } catch (CplexEnv::LicenseError& error) {
126 125
#ifdef LEMON_FORCE_CPLEX_CHECK
127 126
    check(false, error.what());
128 127
#else
129 128
    std::cerr << error.what() << std::endl;
130 129
    std::cerr << "Cplex license check failed, lp check skipped" << std::endl;
131 130
#endif
132 131
  }
133 132
#endif
134 133

	
135 134
  return 0;
136 135

	
137 136
}
0 comments (0 inline)