gravatar
alpar (Alpar Juttner)
alpar@cs.elte.hu
Fix newSolver()/cloneSolver() API in LP tools + doc improvements (#230) - More logical structure for newSolver()/cloneSolver() - Fix compilation problem with gcc-3.3 - Doc improvements
0 13 0
default
13 files changed with 118 insertions and 78 deletions:
↑ Collapse diff ↑
Ignore white space 1536 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 <lemon/clp.h>
20 20
#include <coin/ClpSimplex.hpp>
21 21

	
22 22
namespace lemon {
23 23

	
24 24
  ClpLp::ClpLp() {
25 25
    _prob = new ClpSimplex();
26 26
    _init_temporals();
27 27
    messageLevel(MESSAGE_NO_OUTPUT);
28 28
  }
29 29

	
30 30
  ClpLp::ClpLp(const ClpLp& other) {
31 31
    _prob = new ClpSimplex(*other._prob);
32 32
    rows = other.rows;
33 33
    cols = other.cols;
34 34
    _init_temporals();
35 35
    messageLevel(MESSAGE_NO_OUTPUT);
36 36
  }
37 37

	
38 38
  ClpLp::~ClpLp() {
39 39
    delete _prob;
40 40
    _clear_temporals();
41 41
  }
42 42

	
43 43
  void ClpLp::_init_temporals() {
44 44
    _primal_ray = 0;
45 45
    _dual_ray = 0;
46 46
  }
47 47

	
48 48
  void ClpLp::_clear_temporals() {
49 49
    if (_primal_ray) {
50 50
      delete[] _primal_ray;
51 51
      _primal_ray = 0;
52 52
    }
53 53
    if (_dual_ray) {
54 54
      delete[] _dual_ray;
55 55
      _dual_ray = 0;
56 56
    }
57 57
  }
58 58

	
59
  ClpLp* ClpLp::_newSolver() const {
59
  ClpLp* ClpLp::newSolver() const {
60 60
    ClpLp* newlp = new ClpLp;
61 61
    return newlp;
62 62
  }
63 63

	
64
  ClpLp* ClpLp::_cloneSolver() const {
64
  ClpLp* ClpLp::cloneSolver() const {
65 65
    ClpLp* copylp = new ClpLp(*this);
66 66
    return copylp;
67 67
  }
68 68

	
69 69
  const char* ClpLp::_solverName() const { return "ClpLp"; }
70 70

	
71 71
  int ClpLp::_addCol() {
72 72
    _prob->addColumn(0, 0, 0, -COIN_DBL_MAX, COIN_DBL_MAX, 0.0);
73 73
    return _prob->numberColumns() - 1;
74 74
  }
75 75

	
76 76
  int ClpLp::_addRow() {
77 77
    _prob->addRow(0, 0, 0, -COIN_DBL_MAX, COIN_DBL_MAX);
78 78
    return _prob->numberRows() - 1;
79 79
  }
80 80

	
81 81

	
82 82
  void ClpLp::_eraseCol(int c) {
83 83
    _col_names_ref.erase(_prob->getColumnName(c));
84 84
    _prob->deleteColumns(1, &c);
85 85
  }
86 86

	
87 87
  void ClpLp::_eraseRow(int r) {
88 88
    _row_names_ref.erase(_prob->getRowName(r));
89 89
    _prob->deleteRows(1, &r);
90 90
  }
91 91

	
92 92
  void ClpLp::_eraseColId(int i) {
93 93
    cols.eraseIndex(i);
94 94
    cols.shiftIndices(i);
95 95
  }
96 96

	
97 97
  void ClpLp::_eraseRowId(int i) {
98 98
    rows.eraseIndex(i);
99 99
    rows.shiftIndices(i);
100 100
  }
101 101

	
102 102
  void ClpLp::_getColName(int c, std::string& name) const {
103 103
    name = _prob->getColumnName(c);
104 104
  }
105 105

	
106 106
  void ClpLp::_setColName(int c, const std::string& name) {
107 107
    _prob->setColumnName(c, const_cast<std::string&>(name));
108 108
    _col_names_ref[name] = c;
109 109
  }
110 110

	
111 111
  int ClpLp::_colByName(const std::string& name) const {
112 112
    std::map<std::string, int>::const_iterator it = _col_names_ref.find(name);
113 113
    return it != _col_names_ref.end() ? it->second : -1;
114 114
  }
115 115

	
116 116
  void ClpLp::_getRowName(int r, std::string& name) const {
117 117
    name = _prob->getRowName(r);
118 118
  }
119 119

	
120 120
  void ClpLp::_setRowName(int r, const std::string& name) {
121 121
    _prob->setRowName(r, const_cast<std::string&>(name));
122 122
    _row_names_ref[name] = r;
123 123
  }
124 124

	
125 125
  int ClpLp::_rowByName(const std::string& name) const {
126 126
    std::map<std::string, int>::const_iterator it = _row_names_ref.find(name);
127 127
    return it != _row_names_ref.end() ? it->second : -1;
128 128
  }
129 129

	
130 130

	
131 131
  void ClpLp::_setRowCoeffs(int ix, ExprIterator b, ExprIterator e) {
132 132
    std::map<int, Value> coeffs;
133 133

	
134 134
    int n = _prob->clpMatrix()->getNumCols();
135 135

	
136 136
    const int* indices = _prob->clpMatrix()->getIndices();
137 137
    const double* elements = _prob->clpMatrix()->getElements();
138 138

	
139 139
    for (int i = 0; i < n; ++i) {
140 140
      CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[i];
141 141
      CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[i];
142 142

	
143 143
      const int* it = std::lower_bound(indices + begin, indices + end, ix);
144 144
      if (it != indices + end && *it == ix && elements[it - indices] != 0.0) {
145 145
        coeffs[i] = 0.0;
146 146
      }
147 147
    }
148 148

	
149 149
    for (ExprIterator it = b; it != e; ++it) {
150 150
      coeffs[it->first] = it->second;
151 151
    }
152 152

	
153 153
    for (std::map<int, Value>::iterator it = coeffs.begin();
154 154
         it != coeffs.end(); ++it) {
155 155
      _prob->modifyCoefficient(ix, it->first, it->second);
156 156
    }
157 157
  }
158 158

	
159 159
  void ClpLp::_getRowCoeffs(int ix, InsertIterator b) const {
160 160
    int n = _prob->clpMatrix()->getNumCols();
161 161

	
162 162
    const int* indices = _prob->clpMatrix()->getIndices();
163 163
    const double* elements = _prob->clpMatrix()->getElements();
164 164

	
165 165
    for (int i = 0; i < n; ++i) {
166 166
      CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[i];
167 167
      CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[i];
168 168

	
169 169
      const int* it = std::lower_bound(indices + begin, indices + end, ix);
170 170
      if (it != indices + end && *it == ix) {
171 171
        *b = std::make_pair(i, elements[it - indices]);
172 172
      }
173 173
    }
174 174
  }
175 175

	
176 176
  void ClpLp::_setColCoeffs(int ix, ExprIterator b, ExprIterator e) {
177 177
    std::map<int, Value> coeffs;
178 178

	
179 179
    CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[ix];
180 180
    CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[ix];
181 181

	
182 182
    const int* indices = _prob->clpMatrix()->getIndices();
183 183
    const double* elements = _prob->clpMatrix()->getElements();
184 184

	
185 185
    for (CoinBigIndex i = begin; i != end; ++i) {
186 186
      if (elements[i] != 0.0) {
187 187
        coeffs[indices[i]] = 0.0;
188 188
      }
189 189
    }
190 190
    for (ExprIterator it = b; it != e; ++it) {
191 191
      coeffs[it->first] = it->second;
192 192
    }
193 193
    for (std::map<int, Value>::iterator it = coeffs.begin();
194 194
         it != coeffs.end(); ++it) {
195 195
      _prob->modifyCoefficient(it->first, ix, it->second);
196 196
    }
197 197
  }
198 198

	
199 199
  void ClpLp::_getColCoeffs(int ix, InsertIterator b) const {
200 200
    CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[ix];
201 201
    CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[ix];
202 202

	
203 203
    const int* indices = _prob->clpMatrix()->getIndices();
204 204
    const double* elements = _prob->clpMatrix()->getElements();
205 205

	
206 206
    for (CoinBigIndex i = begin; i != end; ++i) {
207 207
      *b = std::make_pair(indices[i], elements[i]);
208 208
      ++b;
209 209
    }
210 210
  }
211 211

	
212 212
  void ClpLp::_setCoeff(int ix, int jx, Value value) {
213 213
    _prob->modifyCoefficient(ix, jx, value);
214 214
  }
215 215

	
216 216
  ClpLp::Value ClpLp::_getCoeff(int ix, int jx) const {
217 217
    CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[ix];
218 218
    CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[ix];
219 219

	
220 220
    const int* indices = _prob->clpMatrix()->getIndices();
221 221
    const double* elements = _prob->clpMatrix()->getElements();
222 222

	
223 223
    const int* it = std::lower_bound(indices + begin, indices + end, jx);
224 224
    if (it != indices + end && *it == jx) {
225 225
      return elements[it - indices];
226 226
    } else {
227 227
      return 0.0;
228 228
    }
229 229
  }
230 230

	
231 231
  void ClpLp::_setColLowerBound(int i, Value lo) {
232 232
    _prob->setColumnLower(i, lo == - INF ? - COIN_DBL_MAX : lo);
233 233
  }
234 234

	
235 235
  ClpLp::Value ClpLp::_getColLowerBound(int i) const {
236 236
    double val = _prob->getColLower()[i];
237 237
    return val == - COIN_DBL_MAX ? - INF : val;
238 238
  }
239 239

	
240 240
  void ClpLp::_setColUpperBound(int i, Value up) {
241 241
    _prob->setColumnUpper(i, up == INF ? COIN_DBL_MAX : up);
242 242
  }
243 243

	
244 244
  ClpLp::Value ClpLp::_getColUpperBound(int i) const {
245 245
    double val = _prob->getColUpper()[i];
246 246
    return val == COIN_DBL_MAX ? INF : val;
247 247
  }
248 248

	
249 249
  void ClpLp::_setRowLowerBound(int i, Value lo) {
250 250
    _prob->setRowLower(i, lo == - INF ? - COIN_DBL_MAX : lo);
251 251
  }
252 252

	
253 253
  ClpLp::Value ClpLp::_getRowLowerBound(int i) const {
254 254
    double val = _prob->getRowLower()[i];
255 255
    return val == - COIN_DBL_MAX ? - INF : val;
256 256
  }
257 257

	
258 258
  void ClpLp::_setRowUpperBound(int i, Value up) {
259 259
    _prob->setRowUpper(i, up == INF ? COIN_DBL_MAX : up);
260 260
  }
261 261

	
262 262
  ClpLp::Value ClpLp::_getRowUpperBound(int i) const {
263 263
    double val = _prob->getRowUpper()[i];
264 264
    return val == COIN_DBL_MAX ? INF : val;
265 265
  }
266 266

	
267 267
  void ClpLp::_setObjCoeffs(ExprIterator b, ExprIterator e) {
268 268
    int num = _prob->clpMatrix()->getNumCols();
269 269
    for (int i = 0; i < num; ++i) {
270 270
      _prob->setObjectiveCoefficient(i, 0.0);
271 271
    }
272 272
    for (ExprIterator it = b; it != e; ++it) {
273 273
      _prob->setObjectiveCoefficient(it->first, it->second);
274 274
    }
275 275
  }
276 276

	
277 277
  void ClpLp::_getObjCoeffs(InsertIterator b) const {
278 278
    int num = _prob->clpMatrix()->getNumCols();
279 279
    for (int i = 0; i < num; ++i) {
280 280
      Value coef = _prob->getObjCoefficients()[i];
281 281
      if (coef != 0.0) {
282 282
        *b = std::make_pair(i, coef);
283 283
        ++b;
284 284
      }
285 285
    }
286 286
  }
287 287

	
288 288
  void ClpLp::_setObjCoeff(int i, Value obj_coef) {
289 289
    _prob->setObjectiveCoefficient(i, obj_coef);
290 290
  }
291 291

	
292 292
  ClpLp::Value ClpLp::_getObjCoeff(int i) const {
293 293
    return _prob->getObjCoefficients()[i];
294 294
  }
295 295

	
296 296
  ClpLp::SolveExitStatus ClpLp::_solve() {
297 297
    return _prob->primal() >= 0 ? SOLVED : UNSOLVED;
298 298
  }
299 299

	
300 300
  ClpLp::SolveExitStatus ClpLp::solvePrimal() {
301 301
    return _prob->primal() >= 0 ? SOLVED : UNSOLVED;
302 302
  }
303 303

	
304 304
  ClpLp::SolveExitStatus ClpLp::solveDual() {
305 305
    return _prob->dual() >= 0 ? SOLVED : UNSOLVED;
306 306
  }
307 307

	
308 308
  ClpLp::SolveExitStatus ClpLp::solveBarrier() {
309 309
    return _prob->barrier() >= 0 ? SOLVED : UNSOLVED;
310 310
  }
311 311

	
312 312
  ClpLp::Value ClpLp::_getPrimal(int i) const {
313 313
    return _prob->primalColumnSolution()[i];
314 314
  }
315 315
  ClpLp::Value ClpLp::_getPrimalValue() const {
316 316
    return _prob->objectiveValue();
317 317
  }
318 318

	
319 319
  ClpLp::Value ClpLp::_getDual(int i) const {
320 320
    return _prob->dualRowSolution()[i];
321 321
  }
322 322

	
323 323
  ClpLp::Value ClpLp::_getPrimalRay(int i) const {
324 324
    if (!_primal_ray) {
325 325
      _primal_ray = _prob->unboundedRay();
326 326
      LEMON_ASSERT(_primal_ray != 0, "Primal ray is not provided");
327 327
    }
328 328
    return _primal_ray[i];
329 329
  }
330 330

	
331 331
  ClpLp::Value ClpLp::_getDualRay(int i) const {
332 332
    if (!_dual_ray) {
333 333
      _dual_ray = _prob->infeasibilityRay();
334 334
      LEMON_ASSERT(_dual_ray != 0, "Dual ray is not provided");
335 335
    }
336 336
    return _dual_ray[i];
337 337
  }
338 338

	
339 339
  ClpLp::VarStatus ClpLp::_getColStatus(int i) const {
340 340
    switch (_prob->getColumnStatus(i)) {
341 341
    case ClpSimplex::basic:
342 342
      return BASIC;
343 343
    case ClpSimplex::isFree:
344 344
      return FREE;
345 345
    case ClpSimplex::atUpperBound:
346 346
      return UPPER;
347 347
    case ClpSimplex::atLowerBound:
348 348
      return LOWER;
349 349
    case ClpSimplex::isFixed:
350 350
      return FIXED;
351 351
    case ClpSimplex::superBasic:
352 352
      return FREE;
353 353
    default:
354 354
      LEMON_ASSERT(false, "Wrong column status");
355 355
      return VarStatus();
356 356
    }
357 357
  }
358 358

	
359 359
  ClpLp::VarStatus ClpLp::_getRowStatus(int i) const {
360 360
    switch (_prob->getColumnStatus(i)) {
361 361
    case ClpSimplex::basic:
362 362
      return BASIC;
363 363
    case ClpSimplex::isFree:
364 364
      return FREE;
365 365
    case ClpSimplex::atUpperBound:
366 366
      return UPPER;
367 367
    case ClpSimplex::atLowerBound:
368 368
      return LOWER;
369 369
    case ClpSimplex::isFixed:
370 370
      return FIXED;
371 371
    case ClpSimplex::superBasic:
372 372
      return FREE;
373 373
    default:
374 374
      LEMON_ASSERT(false, "Wrong row status");
375 375
      return VarStatus();
376 376
    }
377 377
  }
378 378

	
379 379

	
380 380
  ClpLp::ProblemType ClpLp::_getPrimalType() const {
381 381
    if (_prob->isProvenOptimal()) {
382 382
      return OPTIMAL;
383 383
    } else if (_prob->isProvenPrimalInfeasible()) {
384 384
      return INFEASIBLE;
385 385
    } else if (_prob->isProvenDualInfeasible()) {
386 386
      return UNBOUNDED;
387 387
    } else {
388 388
      return UNDEFINED;
389 389
    }
390 390
  }
391 391

	
392 392
  ClpLp::ProblemType ClpLp::_getDualType() const {
393 393
    if (_prob->isProvenOptimal()) {
394 394
      return OPTIMAL;
395 395
    } else if (_prob->isProvenDualInfeasible()) {
396 396
      return INFEASIBLE;
397 397
    } else if (_prob->isProvenPrimalInfeasible()) {
398 398
      return INFEASIBLE;
399 399
    } else {
400 400
      return UNDEFINED;
401 401
    }
402 402
  }
403 403

	
404 404
  void ClpLp::_setSense(ClpLp::Sense sense) {
405 405
    switch (sense) {
406 406
    case MIN:
407 407
      _prob->setOptimizationDirection(1);
408 408
      break;
409 409
    case MAX:
410 410
      _prob->setOptimizationDirection(-1);
411 411
      break;
412 412
    }
413 413
  }
414 414

	
415 415
  ClpLp::Sense ClpLp::_getSense() const {
416 416
    double dir = _prob->optimizationDirection();
417 417
    if (dir > 0.0) {
418 418
      return MIN;
419 419
    } else {
420 420
      return MAX;
421 421
    }
422 422
  }
423 423

	
424 424
  void ClpLp::_clear() {
425 425
    delete _prob;
426 426
    _prob = new ClpSimplex();
427 427
    rows.clear();
428 428
    cols.clear();
429 429
    _col_names_ref.clear();
430 430
    _clear_temporals();
431 431
  }
432 432

	
433 433
  void ClpLp::messageLevel(MessageLevel m) {
434 434
    _prob->setLogLevel(static_cast<int>(m));
435 435
  }
436 436

	
437 437
} //END OF NAMESPACE LEMON
Ignore white space 1536 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_CLP_H
20 20
#define LEMON_CLP_H
21 21

	
22 22
///\file
23 23
///\brief Header of the LEMON-CLP lp solver interface.
24 24

	
25 25
#include <vector>
26 26
#include <string>
27 27

	
28 28
#include <lemon/lp_base.h>
29 29

	
30 30
class ClpSimplex;
31 31

	
32 32
namespace lemon {
33 33

	
34 34
  /// \ingroup lp_group
35 35
  ///
36 36
  /// \brief Interface for the CLP solver
37 37
  ///
38 38
  /// This class implements an interface for the Clp LP solver.  The
39 39
  /// Clp library is an object oriented lp solver library developed at
40 40
  /// the IBM. The CLP is part of the COIN-OR package and it can be
41 41
  /// used with Common Public License.
42 42
  class ClpLp : public LpSolver {
43 43
  protected:
44 44

	
45 45
    ClpSimplex* _prob;
46 46

	
47 47
    std::map<std::string, int> _col_names_ref;
48 48
    std::map<std::string, int> _row_names_ref;
49 49

	
50 50
  public:
51 51

	
52 52
    /// \e
53 53
    ClpLp();
54 54
    /// \e
55 55
    ClpLp(const ClpLp&);
56 56
    /// \e
57 57
    ~ClpLp();
58 58

	
59
    /// \e
60
    virtual ClpLp* newSolver() const;
61
    /// \e
62
    virtual ClpLp* cloneSolver() const;
63

	
59 64
  protected:
60 65

	
61 66
    mutable double* _primal_ray;
62 67
    mutable double* _dual_ray;
63 68

	
64 69
    void _init_temporals();
65 70
    void _clear_temporals();
66 71

	
67 72
  protected:
68 73

	
69
    virtual ClpLp* _newSolver() const;
70
    virtual ClpLp* _cloneSolver() const;
71

	
72 74
    virtual const char* _solverName() const;
73 75

	
74 76
    virtual int _addCol();
75 77
    virtual int _addRow();
76 78

	
77 79
    virtual void _eraseCol(int i);
78 80
    virtual void _eraseRow(int i);
79 81

	
80 82
    virtual void _eraseColId(int i);
81 83
    virtual void _eraseRowId(int i);
82 84

	
83 85
    virtual void _getColName(int col, std::string& name) const;
84 86
    virtual void _setColName(int col, const std::string& name);
85 87
    virtual int _colByName(const std::string& name) const;
86 88

	
87 89
    virtual void _getRowName(int row, std::string& name) const;
88 90
    virtual void _setRowName(int row, const std::string& name);
89 91
    virtual int _rowByName(const std::string& name) const;
90 92

	
91 93
    virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e);
92 94
    virtual void _getRowCoeffs(int i, InsertIterator b) const;
93 95

	
94 96
    virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e);
95 97
    virtual void _getColCoeffs(int i, InsertIterator b) const;
96 98

	
97 99
    virtual void _setCoeff(int row, int col, Value value);
98 100
    virtual Value _getCoeff(int row, int col) const;
99 101

	
100 102
    virtual void _setColLowerBound(int i, Value value);
101 103
    virtual Value _getColLowerBound(int i) const;
102 104
    virtual void _setColUpperBound(int i, Value value);
103 105
    virtual Value _getColUpperBound(int i) const;
104 106

	
105 107
    virtual void _setRowLowerBound(int i, Value value);
106 108
    virtual Value _getRowLowerBound(int i) const;
107 109
    virtual void _setRowUpperBound(int i, Value value);
108 110
    virtual Value _getRowUpperBound(int i) const;
109 111

	
110 112
    virtual void _setObjCoeffs(ExprIterator, ExprIterator);
111 113
    virtual void _getObjCoeffs(InsertIterator) const;
112 114

	
113 115
    virtual void _setObjCoeff(int i, Value obj_coef);
114 116
    virtual Value _getObjCoeff(int i) const;
115 117

	
116 118
    virtual void _setSense(Sense sense);
117 119
    virtual Sense _getSense() const;
118 120

	
119 121
    virtual SolveExitStatus _solve();
120 122

	
121 123
    virtual Value _getPrimal(int i) const;
122 124
    virtual Value _getDual(int i) const;
123 125

	
124 126
    virtual Value _getPrimalValue() const;
125 127

	
126 128
    virtual Value _getPrimalRay(int i) const;
127 129
    virtual Value _getDualRay(int i) const;
128 130

	
129 131
    virtual VarStatus _getColStatus(int i) const;
130 132
    virtual VarStatus _getRowStatus(int i) const;
131 133

	
132 134
    virtual ProblemType _getPrimalType() const;
133 135
    virtual ProblemType _getDualType() const;
134 136

	
135 137
    virtual void _clear();
136 138

	
137 139
  public:
138 140

	
139 141
    ///Solves LP with primal simplex method.
140 142
    SolveExitStatus solvePrimal();
141 143

	
142 144
    ///Solves LP with dual simplex method.
143 145
    SolveExitStatus solveDual();
144 146

	
145 147
    ///Solves LP with barrier method.
146 148
    SolveExitStatus solveBarrier();
147 149

	
148 150
    ///Returns the constraint identifier understood by CLP.
149 151
    int clpRow(Row r) const { return rows(id(r)); }
150 152

	
151 153
    ///Returns the variable identifier understood by CLP.
152 154
    int clpCol(Col c) const { return cols(id(c)); }
153 155

	
154 156
    ///Enum for \c messageLevel() parameter
155 157
    enum MessageLevel {
156 158
      /// no output (default value)
157 159
      MESSAGE_NO_OUTPUT = 0,
158 160
      /// print final solution
159 161
      MESSAGE_FINAL_SOLUTION = 1,
160 162
      /// print factorization
161 163
      MESSAGE_FACTORIZATION = 2,
162 164
      /// normal output
163 165
      MESSAGE_NORMAL_OUTPUT = 3,
164 166
      /// verbose output
165 167
      MESSAGE_VERBOSE_OUTPUT = 4
166 168
    };
167 169
    ///Set the verbosity of the messages
168 170

	
169 171
    ///Set the verbosity of the messages
170 172
    ///
171 173
    ///\param m is the level of the messages output by the solver routines.
172 174
    void messageLevel(MessageLevel m);
173 175

	
174 176
  };
175 177

	
176 178
} //END OF NAMESPACE LEMON
177 179

	
178 180
#endif //LEMON_CLP_H
179 181

	
Ignore white space 1536 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 <iostream>
20 20
#include <vector>
21 21
#include <cstring>
22 22

	
23 23
#include <lemon/cplex.h>
24 24

	
25 25
extern "C" {
26 26
#include <ilcplex/cplex.h>
27 27
}
28 28

	
29 29

	
30 30
///\file
31 31
///\brief Implementation of the LEMON-CPLEX lp solver interface.
32 32
namespace lemon {
33 33

	
34 34
  CplexEnv::LicenseError::LicenseError(int status) {
35 35
    if (!CPXgeterrorstring(0, status, _message)) {
36 36
      std::strcpy(_message, "Cplex unknown error");
37 37
    }
38 38
  }
39 39

	
40 40
  CplexEnv::CplexEnv() {
41 41
    int status;
42 42
    _cnt = new int;
43 43
    _env = CPXopenCPLEX(&status);
44 44
    if (_env == 0) {
45 45
      delete _cnt;
46 46
      _cnt = 0;
47 47
      throw LicenseError(status);
48 48
    }
49 49
  }
50 50

	
51 51
  CplexEnv::CplexEnv(const CplexEnv& other) {
52 52
    _env = other._env;
53 53
    _cnt = other._cnt;
54 54
    ++(*_cnt);
55 55
  }
56 56

	
57 57
  CplexEnv& CplexEnv::operator=(const CplexEnv& other) {
58 58
    _env = other._env;
59 59
    _cnt = other._cnt;
60 60
    ++(*_cnt);
61 61
    return *this;
62 62
  }
63 63

	
64 64
  CplexEnv::~CplexEnv() {
65 65
    --(*_cnt);
66 66
    if (*_cnt == 0) {
67 67
      delete _cnt;
68 68
      CPXcloseCPLEX(&_env);
69 69
    }
70 70
  }
71 71

	
72 72
  CplexBase::CplexBase() : LpBase() {
73 73
    int status;
74 74
    _prob = CPXcreateprob(cplexEnv(), &status, "Cplex problem");
75 75
  }
76 76

	
77 77
  CplexBase::CplexBase(const CplexEnv& env)
78 78
    : LpBase(), _env(env) {
79 79
    int status;
80 80
    _prob = CPXcreateprob(cplexEnv(), &status, "Cplex problem");
81 81
  }
82 82

	
83 83
  CplexBase::CplexBase(const CplexBase& cplex)
84 84
    : LpBase() {
85 85
    int status;
86 86
    _prob = CPXcloneprob(cplexEnv(), cplex._prob, &status);
87 87
    rows = cplex.rows;
88 88
    cols = cplex.cols;
89 89
  }
90 90

	
91 91
  CplexBase::~CplexBase() {
92 92
    CPXfreeprob(cplexEnv(),&_prob);
93 93
  }
94 94

	
95 95
  int CplexBase::_addCol() {
96 96
    int i = CPXgetnumcols(cplexEnv(), _prob);
97 97
    double lb = -INF, ub = INF;
98 98
    CPXnewcols(cplexEnv(), _prob, 1, 0, &lb, &ub, 0, 0);
99 99
    return i;
100 100
  }
101 101

	
102 102

	
103 103
  int CplexBase::_addRow() {
104 104
    int i = CPXgetnumrows(cplexEnv(), _prob);
105 105
    const double ub = INF;
106 106
    const char s = 'L';
107 107
    CPXnewrows(cplexEnv(), _prob, 1, &ub, &s, 0, 0);
108 108
    return i;
109 109
  }
110 110

	
111 111

	
112 112
  void CplexBase::_eraseCol(int i) {
113 113
    CPXdelcols(cplexEnv(), _prob, i, i);
114 114
  }
115 115

	
116 116
  void CplexBase::_eraseRow(int i) {
117 117
    CPXdelrows(cplexEnv(), _prob, i, i);
118 118
  }
119 119

	
120 120
  void CplexBase::_eraseColId(int i) {
121 121
    cols.eraseIndex(i);
122 122
    cols.shiftIndices(i);
123 123
  }
124 124
  void CplexBase::_eraseRowId(int i) {
125 125
    rows.eraseIndex(i);
126 126
    rows.shiftIndices(i);
127 127
  }
128 128

	
129 129
  void CplexBase::_getColName(int col, std::string &name) const {
130 130
    int size;
131 131
    CPXgetcolname(cplexEnv(), _prob, 0, 0, 0, &size, col, col);
132 132
    if (size == 0) {
133 133
      name.clear();
134 134
      return;
135 135
    }
136 136

	
137 137
    size *= -1;
138 138
    std::vector<char> buf(size);
139 139
    char *cname;
140 140
    int tmp;
141 141
    CPXgetcolname(cplexEnv(), _prob, &cname, &buf.front(), size,
142 142
                  &tmp, col, col);
143 143
    name = cname;
144 144
  }
145 145

	
146 146
  void CplexBase::_setColName(int col, const std::string &name) {
147 147
    char *cname;
148 148
    cname = const_cast<char*>(name.c_str());
149 149
    CPXchgcolname(cplexEnv(), _prob, 1, &col, &cname);
150 150
  }
151 151

	
152 152
  int CplexBase::_colByName(const std::string& name) const {
153 153
    int index;
154 154
    if (CPXgetcolindex(cplexEnv(), _prob,
155 155
                       const_cast<char*>(name.c_str()), &index) == 0) {
156 156
      return index;
157 157
    }
158 158
    return -1;
159 159
  }
160 160

	
161 161
  void CplexBase::_getRowName(int row, std::string &name) const {
162 162
    int size;
163 163
    CPXgetrowname(cplexEnv(), _prob, 0, 0, 0, &size, row, row);
164 164
    if (size == 0) {
165 165
      name.clear();
166 166
      return;
167 167
    }
168 168

	
169 169
    size *= -1;
170 170
    std::vector<char> buf(size);
171 171
    char *cname;
172 172
    int tmp;
173 173
    CPXgetrowname(cplexEnv(), _prob, &cname, &buf.front(), size,
174 174
                  &tmp, row, row);
175 175
    name = cname;
176 176
  }
177 177

	
178 178
  void CplexBase::_setRowName(int row, const std::string &name) {
179 179
    char *cname;
180 180
    cname = const_cast<char*>(name.c_str());
181 181
    CPXchgrowname(cplexEnv(), _prob, 1, &row, &cname);
182 182
  }
183 183

	
184 184
  int CplexBase::_rowByName(const std::string& name) const {
185 185
    int index;
186 186
    if (CPXgetrowindex(cplexEnv(), _prob,
187 187
                       const_cast<char*>(name.c_str()), &index) == 0) {
188 188
      return index;
189 189
    }
190 190
    return -1;
191 191
  }
192 192

	
193 193
  void CplexBase::_setRowCoeffs(int i, ExprIterator b,
194 194
                                      ExprIterator e)
195 195
  {
196 196
    std::vector<int> indices;
197 197
    std::vector<int> rowlist;
198 198
    std::vector<Value> values;
199 199

	
200 200
    for(ExprIterator it=b; it!=e; ++it) {
201 201
      indices.push_back(it->first);
202 202
      values.push_back(it->second);
203 203
      rowlist.push_back(i);
204 204
    }
205 205

	
206 206
    CPXchgcoeflist(cplexEnv(), _prob, values.size(),
207 207
                   &rowlist.front(), &indices.front(), &values.front());
208 208
  }
209 209

	
210 210
  void CplexBase::_getRowCoeffs(int i, InsertIterator b) const {
211 211
    int tmp1, tmp2, tmp3, length;
212 212
    CPXgetrows(cplexEnv(), _prob, &tmp1, &tmp2, 0, 0, 0, &length, i, i);
213 213

	
214 214
    length = -length;
215 215
    std::vector<int> indices(length);
216 216
    std::vector<double> values(length);
217 217

	
218 218
    CPXgetrows(cplexEnv(), _prob, &tmp1, &tmp2,
219 219
               &indices.front(), &values.front(),
220 220
               length, &tmp3, i, i);
221 221

	
222 222
    for (int i = 0; i < length; ++i) {
223 223
      *b = std::make_pair(indices[i], values[i]);
224 224
      ++b;
225 225
    }
226 226
  }
227 227

	
228 228
  void CplexBase::_setColCoeffs(int i, ExprIterator b, ExprIterator e) {
229 229
    std::vector<int> indices;
230 230
    std::vector<int> collist;
231 231
    std::vector<Value> values;
232 232

	
233 233
    for(ExprIterator it=b; it!=e; ++it) {
234 234
      indices.push_back(it->first);
235 235
      values.push_back(it->second);
236 236
      collist.push_back(i);
237 237
    }
238 238

	
239 239
    CPXchgcoeflist(cplexEnv(), _prob, values.size(),
240 240
                   &indices.front(), &collist.front(), &values.front());
241 241
  }
242 242

	
243 243
  void CplexBase::_getColCoeffs(int i, InsertIterator b) const {
244 244

	
245 245
    int tmp1, tmp2, tmp3, length;
246 246
    CPXgetcols(cplexEnv(), _prob, &tmp1, &tmp2, 0, 0, 0, &length, i, i);
247 247

	
248 248
    length = -length;
249 249
    std::vector<int> indices(length);
250 250
    std::vector<double> values(length);
251 251

	
252 252
    CPXgetcols(cplexEnv(), _prob, &tmp1, &tmp2,
253 253
               &indices.front(), &values.front(),
254 254
               length, &tmp3, i, i);
255 255

	
256 256
    for (int i = 0; i < length; ++i) {
257 257
      *b = std::make_pair(indices[i], values[i]);
258 258
      ++b;
259 259
    }
260 260

	
261 261
  }
262 262

	
263 263
  void CplexBase::_setCoeff(int row, int col, Value value) {
264 264
    CPXchgcoef(cplexEnv(), _prob, row, col, value);
265 265
  }
266 266

	
267 267
  CplexBase::Value CplexBase::_getCoeff(int row, int col) const {
268 268
    CplexBase::Value value;
269 269
    CPXgetcoef(cplexEnv(), _prob, row, col, &value);
270 270
    return value;
271 271
  }
272 272

	
273 273
  void CplexBase::_setColLowerBound(int i, Value value) {
274 274
    const char s = 'L';
275 275
    CPXchgbds(cplexEnv(), _prob, 1, &i, &s, &value);
276 276
  }
277 277

	
278 278
  CplexBase::Value CplexBase::_getColLowerBound(int i) const {
279 279
    CplexBase::Value res;
280 280
    CPXgetlb(cplexEnv(), _prob, &res, i, i);
281 281
    return res <= -CPX_INFBOUND ? -INF : res;
282 282
  }
283 283

	
284 284
  void CplexBase::_setColUpperBound(int i, Value value)
285 285
  {
286 286
    const char s = 'U';
287 287
    CPXchgbds(cplexEnv(), _prob, 1, &i, &s, &value);
288 288
  }
289 289

	
290 290
  CplexBase::Value CplexBase::_getColUpperBound(int i) const {
291 291
    CplexBase::Value res;
292 292
    CPXgetub(cplexEnv(), _prob, &res, i, i);
293 293
    return res >= CPX_INFBOUND ? INF : res;
294 294
  }
295 295

	
296 296
  CplexBase::Value CplexBase::_getRowLowerBound(int i) const {
297 297
    char s;
298 298
    CPXgetsense(cplexEnv(), _prob, &s, i, i);
299 299
    CplexBase::Value res;
300 300

	
301 301
    switch (s) {
302 302
    case 'G':
303 303
    case 'R':
304 304
    case 'E':
305 305
      CPXgetrhs(cplexEnv(), _prob, &res, i, i);
306 306
      return res <= -CPX_INFBOUND ? -INF : res;
307 307
    default:
308 308
      return -INF;
309 309
    }
310 310
  }
311 311

	
312 312
  CplexBase::Value CplexBase::_getRowUpperBound(int i) const {
313 313
    char s;
314 314
    CPXgetsense(cplexEnv(), _prob, &s, i, i);
315 315
    CplexBase::Value res;
316 316

	
317 317
    switch (s) {
318 318
    case 'L':
319 319
    case 'E':
320 320
      CPXgetrhs(cplexEnv(), _prob, &res, i, i);
321 321
      return res >= CPX_INFBOUND ? INF : res;
322 322
    case 'R':
323 323
      CPXgetrhs(cplexEnv(), _prob, &res, i, i);
324 324
      {
325 325
        double rng;
326 326
        CPXgetrngval(cplexEnv(), _prob, &rng, i, i);
327 327
        res += rng;
328 328
      }
329 329
      return res >= CPX_INFBOUND ? INF : res;
330 330
    default:
331 331
      return INF;
332 332
    }
333 333
  }
334 334

	
335 335
  //This is easier to implement
336 336
  void CplexBase::_set_row_bounds(int i, Value lb, Value ub) {
337 337
    if (lb == -INF) {
338 338
      const char s = 'L';
339 339
      CPXchgsense(cplexEnv(), _prob, 1, &i, &s);
340 340
      CPXchgrhs(cplexEnv(), _prob, 1, &i, &ub);
341 341
    } else if (ub == INF) {
342 342
      const char s = 'G';
343 343
      CPXchgsense(cplexEnv(), _prob, 1, &i, &s);
344 344
      CPXchgrhs(cplexEnv(), _prob, 1, &i, &lb);
345 345
    } else if (lb == ub){
346 346
      const char s = 'E';
347 347
      CPXchgsense(cplexEnv(), _prob, 1, &i, &s);
348 348
      CPXchgrhs(cplexEnv(), _prob, 1, &i, &lb);
349 349
    } else {
350 350
      const char s = 'R';
351 351
      CPXchgsense(cplexEnv(), _prob, 1, &i, &s);
352 352
      CPXchgrhs(cplexEnv(), _prob, 1, &i, &lb);
353 353
      double len = ub - lb;
354 354
      CPXchgrngval(cplexEnv(), _prob, 1, &i, &len);
355 355
    }
356 356
  }
357 357

	
358 358
  void CplexBase::_setRowLowerBound(int i, Value lb)
359 359
  {
360 360
    LEMON_ASSERT(lb != INF, "Invalid bound");
361 361
    _set_row_bounds(i, lb, CplexBase::_getRowUpperBound(i));
362 362
  }
363 363

	
364 364
  void CplexBase::_setRowUpperBound(int i, Value ub)
365 365
  {
366 366

	
367 367
    LEMON_ASSERT(ub != -INF, "Invalid bound");
368 368
    _set_row_bounds(i, CplexBase::_getRowLowerBound(i), ub);
369 369
  }
370 370

	
371 371
  void CplexBase::_setObjCoeffs(ExprIterator b, ExprIterator e)
372 372
  {
373 373
    std::vector<int> indices;
374 374
    std::vector<Value> values;
375 375
    for(ExprIterator it=b; it!=e; ++it) {
376 376
      indices.push_back(it->first);
377 377
      values.push_back(it->second);
378 378
    }
379 379
    CPXchgobj(cplexEnv(), _prob, values.size(),
380 380
              &indices.front(), &values.front());
381 381

	
382 382
  }
383 383

	
384 384
  void CplexBase::_getObjCoeffs(InsertIterator b) const
385 385
  {
386 386
    int num = CPXgetnumcols(cplexEnv(), _prob);
387 387
    std::vector<Value> x(num);
388 388

	
389 389
    CPXgetobj(cplexEnv(), _prob, &x.front(), 0, num - 1);
390 390
    for (int i = 0; i < num; ++i) {
391 391
      if (x[i] != 0.0) {
392 392
        *b = std::make_pair(i, x[i]);
393 393
        ++b;
394 394
      }
395 395
    }
396 396
  }
397 397

	
398 398
  void CplexBase::_setObjCoeff(int i, Value obj_coef)
399 399
  {
400 400
    CPXchgobj(cplexEnv(), _prob, 1, &i, &obj_coef);
401 401
  }
402 402

	
403 403
  CplexBase::Value CplexBase::_getObjCoeff(int i) const
404 404
  {
405 405
    Value x;
406 406
    CPXgetobj(cplexEnv(), _prob, &x, i, i);
407 407
    return x;
408 408
  }
409 409

	
410 410
  void CplexBase::_setSense(CplexBase::Sense sense) {
411 411
    switch (sense) {
412 412
    case MIN:
413 413
      CPXchgobjsen(cplexEnv(), _prob, CPX_MIN);
414 414
      break;
415 415
    case MAX:
416 416
      CPXchgobjsen(cplexEnv(), _prob, CPX_MAX);
417 417
      break;
418 418
    }
419 419
  }
420 420

	
421 421
  CplexBase::Sense CplexBase::_getSense() const {
422 422
    switch (CPXgetobjsen(cplexEnv(), _prob)) {
423 423
    case CPX_MIN:
424 424
      return MIN;
425 425
    case CPX_MAX:
426 426
      return MAX;
427 427
    default:
428 428
      LEMON_ASSERT(false, "Invalid sense");
429 429
      return CplexBase::Sense();
430 430
    }
431 431
  }
432 432

	
433 433
  void CplexBase::_clear() {
434 434
    CPXfreeprob(cplexEnv(),&_prob);
435 435
    int status;
436 436
    _prob = CPXcreateprob(cplexEnv(), &status, "Cplex problem");
437 437
    rows.clear();
438 438
    cols.clear();
439 439
  }
440 440

	
441 441
  // CplexLp members
442 442

	
443 443
  CplexLp::CplexLp()
444 444
    : LpBase(), CplexBase(), LpSolver() {}
445 445

	
446 446
  CplexLp::CplexLp(const CplexEnv& env)
447 447
    : LpBase(), CplexBase(env), LpSolver() {}
448 448

	
449 449
  CplexLp::CplexLp(const CplexLp& other)
450 450
    : LpBase(), CplexBase(other), LpSolver() {}
451 451

	
452 452
  CplexLp::~CplexLp() {}
453 453

	
454
  CplexLp* CplexLp::_newSolver() const { return new CplexLp; }
455
  CplexLp* CplexLp::_cloneSolver() const {return new CplexLp(*this); }
454
  CplexLp* CplexLp::newSolver() const { return new CplexLp; }
455
  CplexLp* CplexLp::cloneSolver() const {return new CplexLp(*this); }
456 456

	
457 457
  const char* CplexLp::_solverName() const { return "CplexLp"; }
458 458

	
459 459
  void CplexLp::_clear_temporals() {
460 460
    _col_status.clear();
461 461
    _row_status.clear();
462 462
    _primal_ray.clear();
463 463
    _dual_ray.clear();
464 464
  }
465 465

	
466 466
  // The routine returns zero unless an error occurred during the
467 467
  // optimization. Examples of errors include exhausting available
468 468
  // memory (CPXERR_NO_MEMORY) or encountering invalid data in the
469 469
  // CPLEX problem object (CPXERR_NO_PROBLEM). Exceeding a
470 470
  // user-specified CPLEX limit, or proving the model infeasible or
471 471
  // unbounded, are not considered errors. Note that a zero return
472 472
  // value does not necessarily mean that a solution exists. Use query
473 473
  // routines CPXsolninfo, CPXgetstat, and CPXsolution to obtain
474 474
  // further information about the status of the optimization.
475 475
  CplexLp::SolveExitStatus CplexLp::convertStatus(int status) {
476 476
#if CPX_VERSION >= 800
477 477
    if (status == 0) {
478 478
      switch (CPXgetstat(cplexEnv(), _prob)) {
479 479
      case CPX_STAT_OPTIMAL:
480 480
      case CPX_STAT_INFEASIBLE:
481 481
      case CPX_STAT_UNBOUNDED:
482 482
        return SOLVED;
483 483
      default:
484 484
        return UNSOLVED;
485 485
      }
486 486
    } else {
487 487
      return UNSOLVED;
488 488
    }
489 489
#else
490 490
    if (status == 0) {
491 491
      //We want to exclude some cases
492 492
      switch (CPXgetstat(cplexEnv(), _prob)) {
493 493
      case CPX_OBJ_LIM:
494 494
      case CPX_IT_LIM_FEAS:
495 495
      case CPX_IT_LIM_INFEAS:
496 496
      case CPX_TIME_LIM_FEAS:
497 497
      case CPX_TIME_LIM_INFEAS:
498 498
        return UNSOLVED;
499 499
      default:
500 500
        return SOLVED;
501 501
      }
502 502
    } else {
503 503
      return UNSOLVED;
504 504
    }
505 505
#endif
506 506
  }
507 507

	
508 508
  CplexLp::SolveExitStatus CplexLp::_solve() {
509 509
    _clear_temporals();
510 510
    return convertStatus(CPXlpopt(cplexEnv(), _prob));
511 511
  }
512 512

	
513 513
  CplexLp::SolveExitStatus CplexLp::solvePrimal() {
514 514
    _clear_temporals();
515 515
    return convertStatus(CPXprimopt(cplexEnv(), _prob));
516 516
  }
517 517

	
518 518
  CplexLp::SolveExitStatus CplexLp::solveDual() {
519 519
    _clear_temporals();
520 520
    return convertStatus(CPXdualopt(cplexEnv(), _prob));
521 521
  }
522 522

	
523 523
  CplexLp::SolveExitStatus CplexLp::solveBarrier() {
524 524
    _clear_temporals();
525 525
    return convertStatus(CPXbaropt(cplexEnv(), _prob));
526 526
  }
527 527

	
528 528
  CplexLp::Value CplexLp::_getPrimal(int i) const {
529 529
    Value x;
530 530
    CPXgetx(cplexEnv(), _prob, &x, i, i);
531 531
    return x;
532 532
  }
533 533

	
534 534
  CplexLp::Value CplexLp::_getDual(int i) const {
535 535
    Value y;
536 536
    CPXgetpi(cplexEnv(), _prob, &y, i, i);
537 537
    return y;
538 538
  }
539 539

	
540 540
  CplexLp::Value CplexLp::_getPrimalValue() const {
541 541
    Value objval;
542 542
    CPXgetobjval(cplexEnv(), _prob, &objval);
543 543
    return objval;
544 544
  }
545 545

	
546 546
  CplexLp::VarStatus CplexLp::_getColStatus(int i) const {
547 547
    if (_col_status.empty()) {
548 548
      _col_status.resize(CPXgetnumcols(cplexEnv(), _prob));
549 549
      CPXgetbase(cplexEnv(), _prob, &_col_status.front(), 0);
550 550
    }
551 551
    switch (_col_status[i]) {
552 552
    case CPX_BASIC:
553 553
      return BASIC;
554 554
    case CPX_FREE_SUPER:
555 555
      return FREE;
556 556
    case CPX_AT_LOWER:
557 557
      return LOWER;
558 558
    case CPX_AT_UPPER:
559 559
      return UPPER;
560 560
    default:
561 561
      LEMON_ASSERT(false, "Wrong column status");
562 562
      return CplexLp::VarStatus();
563 563
    }
564 564
  }
565 565

	
566 566
  CplexLp::VarStatus CplexLp::_getRowStatus(int i) const {
567 567
    if (_row_status.empty()) {
568 568
      _row_status.resize(CPXgetnumrows(cplexEnv(), _prob));
569 569
      CPXgetbase(cplexEnv(), _prob, 0, &_row_status.front());
570 570
    }
571 571
    switch (_row_status[i]) {
572 572
    case CPX_BASIC:
573 573
      return BASIC;
574 574
    case CPX_AT_LOWER:
575 575
      {
576 576
        char s;
577 577
        CPXgetsense(cplexEnv(), _prob, &s, i, i);
578 578
        return s != 'L' ? LOWER : UPPER;
579 579
      }
580 580
    case CPX_AT_UPPER:
581 581
      return UPPER;
582 582
    default:
583 583
      LEMON_ASSERT(false, "Wrong row status");
584 584
      return CplexLp::VarStatus();
585 585
    }
586 586
  }
587 587

	
588 588
  CplexLp::Value CplexLp::_getPrimalRay(int i) const {
589 589
    if (_primal_ray.empty()) {
590 590
      _primal_ray.resize(CPXgetnumcols(cplexEnv(), _prob));
591 591
      CPXgetray(cplexEnv(), _prob, &_primal_ray.front());
592 592
    }
593 593
    return _primal_ray[i];
594 594
  }
595 595

	
596 596
  CplexLp::Value CplexLp::_getDualRay(int i) const {
597 597
    if (_dual_ray.empty()) {
598 598

	
599 599
    }
600 600
    return _dual_ray[i];
601 601
  }
602 602

	
603 603
  //7.5-os cplex statusai (Vigyazat: a 9.0-asei masok!)
604 604
  // This table lists the statuses, returned by the CPXgetstat()
605 605
  // routine, for solutions to LP problems or mixed integer problems. If
606 606
  // no solution exists, the return value is zero.
607 607

	
608 608
  // For Simplex, Barrier
609 609
  // 1          CPX_OPTIMAL
610 610
  //          Optimal solution found
611 611
  // 2          CPX_INFEASIBLE
612 612
  //          Problem infeasible
613 613
  // 3    CPX_UNBOUNDED
614 614
  //          Problem unbounded
615 615
  // 4          CPX_OBJ_LIM
616 616
  //          Objective limit exceeded in Phase II
617 617
  // 5          CPX_IT_LIM_FEAS
618 618
  //          Iteration limit exceeded in Phase II
619 619
  // 6          CPX_IT_LIM_INFEAS
620 620
  //          Iteration limit exceeded in Phase I
621 621
  // 7          CPX_TIME_LIM_FEAS
622 622
  //          Time limit exceeded in Phase II
623 623
  // 8          CPX_TIME_LIM_INFEAS
624 624
  //          Time limit exceeded in Phase I
625 625
  // 9          CPX_NUM_BEST_FEAS
626 626
  //          Problem non-optimal, singularities in Phase II
627 627
  // 10         CPX_NUM_BEST_INFEAS
628 628
  //          Problem non-optimal, singularities in Phase I
629 629
  // 11         CPX_OPTIMAL_INFEAS
630 630
  //          Optimal solution found, unscaled infeasibilities
631 631
  // 12         CPX_ABORT_FEAS
632 632
  //          Aborted in Phase II
633 633
  // 13         CPX_ABORT_INFEAS
634 634
  //          Aborted in Phase I
635 635
  // 14          CPX_ABORT_DUAL_INFEAS
636 636
  //          Aborted in barrier, dual infeasible
637 637
  // 15          CPX_ABORT_PRIM_INFEAS
638 638
  //          Aborted in barrier, primal infeasible
639 639
  // 16          CPX_ABORT_PRIM_DUAL_INFEAS
640 640
  //          Aborted in barrier, primal and dual infeasible
641 641
  // 17          CPX_ABORT_PRIM_DUAL_FEAS
642 642
  //          Aborted in barrier, primal and dual feasible
643 643
  // 18          CPX_ABORT_CROSSOVER
644 644
  //          Aborted in crossover
645 645
  // 19          CPX_INForUNBD
646 646
  //          Infeasible or unbounded
647 647
  // 20   CPX_PIVOT
648 648
  //       User pivot used
649 649
  //
650 650
  //     Ezeket hova tegyem:
651 651
  // ??case CPX_ABORT_DUAL_INFEAS
652 652
  // ??case CPX_ABORT_CROSSOVER
653 653
  // ??case CPX_INForUNBD
654 654
  // ??case CPX_PIVOT
655 655

	
656 656
  //Some more interesting stuff:
657 657

	
658 658
  // CPX_PARAM_PROBMETHOD  1062  int  LPMETHOD
659 659
  // 0 Automatic
660 660
  // 1 Primal Simplex
661 661
  // 2 Dual Simplex
662 662
  // 3 Network Simplex
663 663
  // 4 Standard Barrier
664 664
  // Default: 0
665 665
  // Description: Method for linear optimization.
666 666
  // Determines which algorithm is used when CPXlpopt() (or "optimize"
667 667
  // in the Interactive Optimizer) is called. Currently the behavior of
668 668
  // the "Automatic" setting is that CPLEX simply invokes the dual
669 669
  // simplex method, but this capability may be expanded in the future
670 670
  // so that CPLEX chooses the method based on problem characteristics
671 671
#if CPX_VERSION < 900
672 672
  void statusSwitch(CPXENVptr cplexEnv(),int& stat){
673 673
    int lpmethod;
674 674
    CPXgetintparam (cplexEnv(),CPX_PARAM_PROBMETHOD,&lpmethod);
675 675
    if (lpmethod==2){
676 676
      if (stat==CPX_UNBOUNDED){
677 677
        stat=CPX_INFEASIBLE;
678 678
      }
679 679
      else{
680 680
        if (stat==CPX_INFEASIBLE)
681 681
          stat=CPX_UNBOUNDED;
682 682
      }
683 683
    }
684 684
  }
685 685
#else
686 686
  void statusSwitch(CPXENVptr,int&){}
687 687
#endif
688 688

	
689 689
  CplexLp::ProblemType CplexLp::_getPrimalType() const {
690 690
    // Unboundedness not treated well: the following is from cplex 9.0 doc
691 691
    // About Unboundedness
692 692

	
693 693
    // The treatment of models that are unbounded involves a few
694 694
    // subtleties. Specifically, a declaration of unboundedness means that
695 695
    // ILOG CPLEX has determined that the model has an unbounded
696 696
    // ray. Given any feasible solution x with objective z, a multiple of
697 697
    // the unbounded ray can be added to x to give a feasible solution
698 698
    // with objective z-1 (or z+1 for maximization models). Thus, if a
699 699
    // feasible solution exists, then the optimal objective is
700 700
    // unbounded. Note that ILOG CPLEX has not necessarily concluded that
701 701
    // a feasible solution exists. Users can call the routine CPXsolninfo
702 702
    // to determine whether ILOG CPLEX has also concluded that the model
703 703
    // has a feasible solution.
704 704

	
705 705
    int stat = CPXgetstat(cplexEnv(), _prob);
706 706
#if CPX_VERSION >= 800
707 707
    switch (stat)
708 708
      {
709 709
      case CPX_STAT_OPTIMAL:
710 710
        return OPTIMAL;
711 711
      case CPX_STAT_UNBOUNDED:
712 712
        return UNBOUNDED;
713 713
      case CPX_STAT_INFEASIBLE:
714 714
        return INFEASIBLE;
715 715
      default:
716 716
        return UNDEFINED;
717 717
      }
718 718
#else
719 719
    statusSwitch(cplexEnv(),stat);
720 720
    //CPXgetstat(cplexEnv(), _prob);
721 721
    //printf("A primal status: %d, CPX_OPTIMAL=%d \n",stat,CPX_OPTIMAL);
722 722
    switch (stat) {
723 723
    case 0:
724 724
      return UNDEFINED; //Undefined
725 725
    case CPX_OPTIMAL://Optimal
726 726
      return OPTIMAL;
727 727
    case CPX_UNBOUNDED://Unbounded
728 728
      return INFEASIBLE;//In case of dual simplex
729 729
      //return UNBOUNDED;
730 730
    case CPX_INFEASIBLE://Infeasible
731 731
      //    case CPX_IT_LIM_INFEAS:
732 732
      //     case CPX_TIME_LIM_INFEAS:
733 733
      //     case CPX_NUM_BEST_INFEAS:
734 734
      //     case CPX_OPTIMAL_INFEAS:
735 735
      //     case CPX_ABORT_INFEAS:
736 736
      //     case CPX_ABORT_PRIM_INFEAS:
737 737
      //     case CPX_ABORT_PRIM_DUAL_INFEAS:
738 738
      return UNBOUNDED;//In case of dual simplex
739 739
      //return INFEASIBLE;
740 740
      //     case CPX_OBJ_LIM:
741 741
      //     case CPX_IT_LIM_FEAS:
742 742
      //     case CPX_TIME_LIM_FEAS:
743 743
      //     case CPX_NUM_BEST_FEAS:
744 744
      //     case CPX_ABORT_FEAS:
745 745
      //     case CPX_ABORT_PRIM_DUAL_FEAS:
746 746
      //       return FEASIBLE;
747 747
    default:
748 748
      return UNDEFINED; //Everything else comes here
749 749
      //FIXME error
750 750
    }
751 751
#endif
752 752
  }
753 753

	
754 754
  //9.0-as cplex verzio statusai
755 755
  // CPX_STAT_ABORT_DUAL_OBJ_LIM
756 756
  // CPX_STAT_ABORT_IT_LIM
757 757
  // CPX_STAT_ABORT_OBJ_LIM
758 758
  // CPX_STAT_ABORT_PRIM_OBJ_LIM
759 759
  // CPX_STAT_ABORT_TIME_LIM
760 760
  // CPX_STAT_ABORT_USER
761 761
  // CPX_STAT_FEASIBLE_RELAXED
762 762
  // CPX_STAT_INFEASIBLE
763 763
  // CPX_STAT_INForUNBD
764 764
  // CPX_STAT_NUM_BEST
765 765
  // CPX_STAT_OPTIMAL
766 766
  // CPX_STAT_OPTIMAL_FACE_UNBOUNDED
767 767
  // CPX_STAT_OPTIMAL_INFEAS
768 768
  // CPX_STAT_OPTIMAL_RELAXED
769 769
  // CPX_STAT_UNBOUNDED
770 770

	
771 771
  CplexLp::ProblemType CplexLp::_getDualType() const {
772 772
    int stat = CPXgetstat(cplexEnv(), _prob);
773 773
#if CPX_VERSION >= 800
774 774
    switch (stat) {
775 775
    case CPX_STAT_OPTIMAL:
776 776
      return OPTIMAL;
777 777
    case CPX_STAT_UNBOUNDED:
778 778
      return INFEASIBLE;
779 779
    default:
780 780
      return UNDEFINED;
781 781
    }
782 782
#else
783 783
    statusSwitch(cplexEnv(),stat);
784 784
    switch (stat) {
785 785
    case 0:
786 786
      return UNDEFINED; //Undefined
787 787
    case CPX_OPTIMAL://Optimal
788 788
      return OPTIMAL;
789 789
    case CPX_UNBOUNDED:
790 790
      return INFEASIBLE;
791 791
    default:
792 792
      return UNDEFINED; //Everything else comes here
793 793
      //FIXME error
794 794
    }
795 795
#endif
796 796
  }
797 797

	
798 798
  // CplexMip members
799 799

	
800 800
  CplexMip::CplexMip()
801 801
    : LpBase(), CplexBase(), MipSolver() {
802 802

	
803 803
#if CPX_VERSION < 800
804 804
    CPXchgprobtype(cplexEnv(),  _prob, CPXPROB_MIP);
805 805
#else
806 806
    CPXchgprobtype(cplexEnv(),  _prob, CPXPROB_MILP);
807 807
#endif
808 808
  }
809 809

	
810 810
  CplexMip::CplexMip(const CplexEnv& env)
811 811
    : LpBase(), CplexBase(env), MipSolver() {
812 812

	
813 813
#if CPX_VERSION < 800
814 814
    CPXchgprobtype(cplexEnv(),  _prob, CPXPROB_MIP);
815 815
#else
816 816
    CPXchgprobtype(cplexEnv(),  _prob, CPXPROB_MILP);
817 817
#endif
818 818

	
819 819
  }
820 820

	
821 821
  CplexMip::CplexMip(const CplexMip& other)
822 822
    : LpBase(), CplexBase(other), MipSolver() {}
823 823

	
824 824
  CplexMip::~CplexMip() {}
825 825

	
826
  CplexMip* CplexMip::_newSolver() const { return new CplexMip; }
827
  CplexMip* CplexMip::_cloneSolver() const {return new CplexMip(*this); }
826
  CplexMip* CplexMip::newSolver() const { return new CplexMip; }
827
  CplexMip* CplexMip::cloneSolver() const {return new CplexMip(*this); }
828 828

	
829 829
  const char* CplexMip::_solverName() const { return "CplexMip"; }
830 830

	
831 831
  void CplexMip::_setColType(int i, CplexMip::ColTypes col_type) {
832 832

	
833 833
    // Note If a variable is to be changed to binary, a call to CPXchgbds
834 834
    // should also be made to change the bounds to 0 and 1.
835 835

	
836 836
    switch (col_type){
837 837
    case INTEGER: {
838 838
      const char t = 'I';
839 839
      CPXchgctype (cplexEnv(), _prob, 1, &i, &t);
840 840
    } break;
841 841
    case REAL: {
842 842
      const char t = 'C';
843 843
      CPXchgctype (cplexEnv(), _prob, 1, &i, &t);
844 844
    } break;
845 845
    default:
846 846
      break;
847 847
    }
848 848
  }
849 849

	
850 850
  CplexMip::ColTypes CplexMip::_getColType(int i) const {
851 851
    char t;
852 852
    CPXgetctype (cplexEnv(), _prob, &t, i, i);
853 853
    switch (t) {
854 854
    case 'I':
855 855
      return INTEGER;
856 856
    case 'C':
857 857
      return REAL;
858 858
    default:
859 859
      LEMON_ASSERT(false, "Invalid column type");
860 860
      return ColTypes();
861 861
    }
862 862

	
863 863
  }
864 864

	
865 865
  CplexMip::SolveExitStatus CplexMip::_solve() {
866 866
    int status;
867 867
    status = CPXmipopt (cplexEnv(), _prob);
868 868
    if (status==0)
869 869
      return SOLVED;
870 870
    else
871 871
      return UNSOLVED;
872 872

	
873 873
  }
874 874

	
875 875

	
876 876
  CplexMip::ProblemType CplexMip::_getType() const {
877 877

	
878 878
    int stat = CPXgetstat(cplexEnv(), _prob);
879 879

	
880 880
    //Fortunately, MIP statuses did not change for cplex 8.0
881 881
    switch (stat) {
882 882
    case CPXMIP_OPTIMAL:
883 883
      // Optimal integer solution has been found.
884 884
    case CPXMIP_OPTIMAL_TOL:
885 885
      // Optimal soluton with the tolerance defined by epgap or epagap has
886 886
      // been found.
887 887
      return OPTIMAL;
888 888
      //This also exists in later issues
889 889
      //    case CPXMIP_UNBOUNDED:
890 890
      //return UNBOUNDED;
891 891
      case CPXMIP_INFEASIBLE:
892 892
        return INFEASIBLE;
893 893
    default:
894 894
      return UNDEFINED;
895 895
    }
896 896
    //Unboundedness not treated well: the following is from cplex 9.0 doc
897 897
    // About Unboundedness
898 898

	
899 899
    // The treatment of models that are unbounded involves a few
900 900
    // subtleties. Specifically, a declaration of unboundedness means that
901 901
    // ILOG CPLEX has determined that the model has an unbounded
902 902
    // ray. Given any feasible solution x with objective z, a multiple of
903 903
    // the unbounded ray can be added to x to give a feasible solution
904 904
    // with objective z-1 (or z+1 for maximization models). Thus, if a
905 905
    // feasible solution exists, then the optimal objective is
906 906
    // unbounded. Note that ILOG CPLEX has not necessarily concluded that
907 907
    // a feasible solution exists. Users can call the routine CPXsolninfo
908 908
    // to determine whether ILOG CPLEX has also concluded that the model
909 909
    // has a feasible solution.
910 910
  }
911 911

	
912 912
  CplexMip::Value CplexMip::_getSol(int i) const {
913 913
    Value x;
914 914
    CPXgetmipx(cplexEnv(), _prob, &x, i, i);
915 915
    return x;
916 916
  }
917 917

	
918 918
  CplexMip::Value CplexMip::_getSolValue() const {
919 919
    Value objval;
920 920
    CPXgetmipobjval(cplexEnv(), _prob, &objval);
921 921
    return objval;
922 922
  }
923 923

	
924 924
} //namespace lemon
925 925

	
Ignore white space 1536 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_CPLEX_H
20 20
#define LEMON_CPLEX_H
21 21

	
22 22
///\file
23 23
///\brief Header of the LEMON-CPLEX lp solver interface.
24 24

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

	
27 27
struct cpxenv;
28 28
struct cpxlp;
29 29

	
30 30
namespace lemon {
31 31

	
32 32
  /// \brief Reference counted wrapper around cpxenv pointer
33 33
  ///
34 34
  /// The cplex uses environment object which is responsible for
35 35
  /// checking the proper license usage. This class provides a simple
36 36
  /// interface for share the environment object between different
37 37
  /// problems.
38 38
  class CplexEnv {
39 39
    friend class CplexBase;
40 40
  private:
41 41
    cpxenv* _env;
42 42
    mutable int* _cnt;
43 43

	
44 44
  public:
45 45

	
46 46
    /// \brief This exception is thrown when the license check is not
47 47
    /// sufficient
48 48
    class LicenseError : public Exception {
49 49
      friend class CplexEnv;
50 50
    private:
51 51

	
52 52
      LicenseError(int status);
53 53
      char _message[510];
54 54

	
55 55
    public:
56 56

	
57 57
      /// The short error message
58 58
      virtual const char* what() const throw() {
59 59
        return _message;
60 60
      }
61 61
    };
62 62

	
63 63
    /// Constructor
64 64
    CplexEnv();
65 65
    /// Shallow copy constructor
66 66
    CplexEnv(const CplexEnv&);
67 67
    /// Shallow assignement
68 68
    CplexEnv& operator=(const CplexEnv&);
69 69
    /// Destructor
70 70
    virtual ~CplexEnv();
71 71

	
72 72
  protected:
73 73

	
74 74
    cpxenv* cplexEnv() { return _env; }
75 75
    const cpxenv* cplexEnv() const { return _env; }
76 76
  };
77 77

	
78 78
  /// \brief Base interface for the CPLEX LP and MIP solver
79 79
  ///
80 80
  /// This class implements the common interface of the CPLEX LP and
81 81
  /// MIP solvers.  
82 82
  /// \ingroup lp_group
83 83
  class CplexBase : virtual public LpBase {
84 84
  protected:
85 85

	
86 86
    CplexEnv _env;
87 87
    cpxlp* _prob;
88 88

	
89 89
    CplexBase();
90 90
    CplexBase(const CplexEnv&);
91 91
    CplexBase(const CplexBase &);
92 92
    virtual ~CplexBase();
93 93

	
94 94
    virtual int _addCol();
95 95
    virtual int _addRow();
96 96

	
97 97
    virtual void _eraseCol(int i);
98 98
    virtual void _eraseRow(int i);
99 99

	
100 100
    virtual void _eraseColId(int i);
101 101
    virtual void _eraseRowId(int i);
102 102

	
103 103
    virtual void _getColName(int col, std::string& name) const;
104 104
    virtual void _setColName(int col, const std::string& name);
105 105
    virtual int _colByName(const std::string& name) const;
106 106

	
107 107
    virtual void _getRowName(int row, std::string& name) const;
108 108
    virtual void _setRowName(int row, const std::string& name);
109 109
    virtual int _rowByName(const std::string& name) const;
110 110

	
111 111
    virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e);
112 112
    virtual void _getRowCoeffs(int i, InsertIterator b) const;
113 113

	
114 114
    virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e);
115 115
    virtual void _getColCoeffs(int i, InsertIterator b) const;
116 116

	
117 117
    virtual void _setCoeff(int row, int col, Value value);
118 118
    virtual Value _getCoeff(int row, int col) const;
119 119

	
120 120
    virtual void _setColLowerBound(int i, Value value);
121 121
    virtual Value _getColLowerBound(int i) const;
122 122

	
123 123
    virtual void _setColUpperBound(int i, Value value);
124 124
    virtual Value _getColUpperBound(int i) const;
125 125

	
126 126
  private:
127 127
    void _set_row_bounds(int i, Value lb, Value ub);
128 128
  protected:
129 129

	
130 130
    virtual void _setRowLowerBound(int i, Value value);
131 131
    virtual Value _getRowLowerBound(int i) const;
132 132

	
133 133
    virtual void _setRowUpperBound(int i, Value value);
134 134
    virtual Value _getRowUpperBound(int i) const;
135 135

	
136 136
    virtual void _setObjCoeffs(ExprIterator b, ExprIterator e);
137 137
    virtual void _getObjCoeffs(InsertIterator b) const;
138 138

	
139 139
    virtual void _setObjCoeff(int i, Value obj_coef);
140 140
    virtual Value _getObjCoeff(int i) const;
141 141

	
142 142
    virtual void _setSense(Sense sense);
143 143
    virtual Sense _getSense() const;
144 144

	
145 145
    virtual void _clear();
146 146

	
147 147
  public:
148 148

	
149 149
    /// Returns the used \c CplexEnv instance
150 150
    const CplexEnv& env() const { return _env; }
151 151
    ///
152 152
    const cpxenv* cplexEnv() const { return _env.cplexEnv(); }
153 153

	
154 154
    cpxlp* cplexLp() { return _prob; }
155 155
    const cpxlp* cplexLp() const { return _prob; }
156 156

	
157 157
  };
158 158

	
159 159
  /// \brief Interface for the CPLEX LP solver
160 160
  ///
161 161
  /// This class implements an interface for the CPLEX LP solver.
162 162
  ///\ingroup lp_group
163
  class CplexLp : public CplexBase, public LpSolver {
163
  class CplexLp : public LpSolver, public CplexBase {
164 164
  public:
165 165
    /// \e
166 166
    CplexLp();
167 167
    /// \e
168 168
    CplexLp(const CplexEnv&);
169 169
    /// \e
170 170
    CplexLp(const CplexLp&);
171 171
    /// \e
172 172
    virtual ~CplexLp();
173 173

	
174
    /// \e
175
    virtual CplexLp* cloneSolver() const;
176
    /// \e
177
    virtual CplexLp* newSolver() const;
178

	
174 179
  private:
175 180

	
176 181
    // these values cannot retrieved element by element
177 182
    mutable std::vector<int> _col_status;
178 183
    mutable std::vector<int> _row_status;
179 184

	
180 185
    mutable std::vector<Value> _primal_ray;
181 186
    mutable std::vector<Value> _dual_ray;
182 187

	
183 188
    void _clear_temporals();
184 189

	
185 190
    SolveExitStatus convertStatus(int status);
186 191

	
187 192
  protected:
188 193

	
189
    virtual CplexLp* _cloneSolver() const;
190
    virtual CplexLp* _newSolver() const;
191

	
192 194
    virtual const char* _solverName() const;
193 195

	
194 196
    virtual SolveExitStatus _solve();
195 197
    virtual Value _getPrimal(int i) const;
196 198
    virtual Value _getDual(int i) const;
197 199
    virtual Value _getPrimalValue() const;
198 200

	
199 201
    virtual VarStatus _getColStatus(int i) const;
200 202
    virtual VarStatus _getRowStatus(int i) const;
201 203

	
202 204
    virtual Value _getPrimalRay(int i) const;
203 205
    virtual Value _getDualRay(int i) const;
204 206

	
205 207
    virtual ProblemType _getPrimalType() const;
206 208
    virtual ProblemType _getDualType() const;
207 209

	
208 210
  public:
209 211

	
210 212
    /// Solve with primal simplex method
211 213
    SolveExitStatus solvePrimal();
212 214

	
213 215
    /// Solve with dual simplex method
214 216
    SolveExitStatus solveDual();
215 217

	
216 218
    /// Solve with barrier method
217 219
    SolveExitStatus solveBarrier();
218 220

	
219 221
  };
220 222

	
221 223
  /// \brief Interface for the CPLEX MIP solver
222 224
  ///
223 225
  /// This class implements an interface for the CPLEX MIP solver.
224 226
  ///\ingroup lp_group
225
  class CplexMip : public CplexBase, public MipSolver {
227
  class CplexMip : public MipSolver, public CplexBase {
226 228
  public:
227 229
    /// \e
228 230
    CplexMip();
229 231
    /// \e
230 232
    CplexMip(const CplexEnv&);
231 233
    /// \e
232 234
    CplexMip(const CplexMip&);
233 235
    /// \e
234 236
    virtual ~CplexMip();
235 237

	
236 238
  protected:
237 239

	
238 240
    virtual CplexMip* _cloneSolver() const;
239 241
    virtual CplexMip* _newSolver() const;
240 242

	
241 243
    virtual const char* _solverName() const;
242 244

	
243 245
    virtual ColTypes _getColType(int col) const;
244 246
    virtual void _setColType(int col, ColTypes col_type);
245 247

	
246 248
    virtual SolveExitStatus _solve();
247 249
    virtual ProblemType _getType() const;
248 250
    virtual Value _getSol(int i) const;
249 251
    virtual Value _getSolValue() const;
250 252

	
251 253
  };
252 254

	
253 255
} //END OF NAMESPACE LEMON
254 256

	
255 257
#endif //LEMON_CPLEX_H
256 258

	
Ignore white space 1536 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
  // GlpkLp members
526 526

	
527 527
  GlpkLp::GlpkLp()
528 528
    : LpBase(), GlpkBase(), LpSolver() {
529 529
    messageLevel(MESSAGE_NO_OUTPUT);
530 530
  }
531 531

	
532 532
  GlpkLp::GlpkLp(const GlpkLp& other)
533 533
    : LpBase(other), GlpkBase(other), LpSolver(other) {
534 534
    messageLevel(MESSAGE_NO_OUTPUT);
535 535
  }
536 536

	
537
  GlpkLp* GlpkLp::_newSolver() const { return new GlpkLp; }
538
  GlpkLp* GlpkLp::_cloneSolver() const { return new GlpkLp(*this); }
537
  GlpkLp* GlpkLp::newSolver() const { return new GlpkLp; }
538
  GlpkLp* GlpkLp::cloneSolver() const { return new GlpkLp(*this); }
539 539

	
540 540
  const char* GlpkLp::_solverName() const { return "GlpkLp"; }
541 541

	
542 542
  void GlpkLp::_clear_temporals() {
543 543
    _primal_ray.clear();
544 544
    _dual_ray.clear();
545 545
  }
546 546

	
547 547
  GlpkLp::SolveExitStatus GlpkLp::_solve() {
548 548
    return solvePrimal();
549 549
  }
550 550

	
551 551
  GlpkLp::SolveExitStatus GlpkLp::solvePrimal() {
552 552
    _clear_temporals();
553 553

	
554 554
    glp_smcp smcp;
555 555
    glp_init_smcp(&smcp);
556 556

	
557 557
    switch (_message_level) {
558 558
    case MESSAGE_NO_OUTPUT:
559 559
      smcp.msg_lev = GLP_MSG_OFF;
560 560
      break;
561 561
    case MESSAGE_ERROR_MESSAGE:
562 562
      smcp.msg_lev = GLP_MSG_ERR;
563 563
      break;
564 564
    case MESSAGE_NORMAL_OUTPUT:
565 565
      smcp.msg_lev = GLP_MSG_ON;
566 566
      break;
567 567
    case MESSAGE_FULL_OUTPUT:
568 568
      smcp.msg_lev = GLP_MSG_ALL;
569 569
      break;
570 570
    }
571 571

	
572 572
    if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
573 573
    return SOLVED;
574 574
  }
575 575

	
576 576
  GlpkLp::SolveExitStatus GlpkLp::solveDual() {
577 577
    _clear_temporals();
578 578

	
579 579
    glp_smcp smcp;
580 580
    glp_init_smcp(&smcp);
581 581

	
582 582
    switch (_message_level) {
583 583
    case MESSAGE_NO_OUTPUT:
584 584
      smcp.msg_lev = GLP_MSG_OFF;
585 585
      break;
586 586
    case MESSAGE_ERROR_MESSAGE:
587 587
      smcp.msg_lev = GLP_MSG_ERR;
588 588
      break;
589 589
    case MESSAGE_NORMAL_OUTPUT:
590 590
      smcp.msg_lev = GLP_MSG_ON;
591 591
      break;
592 592
    case MESSAGE_FULL_OUTPUT:
593 593
      smcp.msg_lev = GLP_MSG_ALL;
594 594
      break;
595 595
    }
596 596
    smcp.meth = GLP_DUAL;
597 597

	
598 598
    if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
599 599
    return SOLVED;
600 600
  }
601 601

	
602 602
  GlpkLp::Value GlpkLp::_getPrimal(int i) const {
603 603
    return glp_get_col_prim(lp, i);
604 604
  }
605 605

	
606 606
  GlpkLp::Value GlpkLp::_getDual(int i) const {
607 607
    return glp_get_row_dual(lp, i);
608 608
  }
609 609

	
610 610
  GlpkLp::Value GlpkLp::_getPrimalValue() const {
611 611
    return glp_get_obj_val(lp);
612 612
  }
613 613

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

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

	
650 650
  GlpkLp::Value GlpkLp::_getPrimalRay(int i) const {
651 651
    if (_primal_ray.empty()) {
652 652
      int row_num = glp_get_num_rows(lp);
653 653
      int col_num = glp_get_num_cols(lp);
654 654

	
655 655
      _primal_ray.resize(col_num + 1, 0.0);
656 656

	
657 657
      int index = glp_get_unbnd_ray(lp);
658 658
      if (index != 0) {
659 659
        // The primal ray is found in primal simplex second phase
660 660
        LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :
661 661
                      glp_get_col_stat(lp, index - row_num)) != GLP_BS,
662 662
                     "Wrong primal ray");
663 663

	
664 664
        bool negate = glp_get_obj_dir(lp) == GLP_MAX;
665 665

	
666 666
        if (index > row_num) {
667 667
          _primal_ray[index - row_num] = 1.0;
668 668
          if (glp_get_col_dual(lp, index - row_num) > 0) {
669 669
            negate = !negate;
670 670
          }
671 671
        } else {
672 672
          if (glp_get_row_dual(lp, index) > 0) {
673 673
            negate = !negate;
674 674
          }
675 675
        }
676 676

	
677 677
        std::vector<int> ray_indexes(row_num + 1);
678 678
        std::vector<Value> ray_values(row_num + 1);
679 679
        int ray_length = glp_eval_tab_col(lp, index, &ray_indexes.front(),
680 680
                                          &ray_values.front());
681 681

	
682 682
        for (int i = 1; i <= ray_length; ++i) {
683 683
          if (ray_indexes[i] > row_num) {
684 684
            _primal_ray[ray_indexes[i] - row_num] = ray_values[i];
685 685
          }
686 686
        }
687 687

	
688 688
        if (negate) {
689 689
          for (int i = 1; i <= col_num; ++i) {
690 690
            _primal_ray[i] = - _primal_ray[i];
691 691
          }
692 692
        }
693 693
      } else {
694 694
        for (int i = 1; i <= col_num; ++i) {
695 695
          _primal_ray[i] = glp_get_col_prim(lp, i);
696 696
        }
697 697
      }
698 698
    }
699 699
    return _primal_ray[i];
700 700
  }
701 701

	
702 702
  GlpkLp::Value GlpkLp::_getDualRay(int i) const {
703 703
    if (_dual_ray.empty()) {
704 704
      int row_num = glp_get_num_rows(lp);
705 705

	
706 706
      _dual_ray.resize(row_num + 1, 0.0);
707 707

	
708 708
      int index = glp_get_unbnd_ray(lp);
709 709
      if (index != 0) {
710 710
        // The dual ray is found in dual simplex second phase
711 711
        LEMON_ASSERT((index <= row_num ? glp_get_row_stat(lp, index) :
712 712
                      glp_get_col_stat(lp, index - row_num)) == GLP_BS,
713 713

	
714 714
                     "Wrong dual ray");
715 715

	
716 716
        int idx;
717 717
        bool negate = false;
718 718

	
719 719
        if (index > row_num) {
720 720
          idx = glp_get_col_bind(lp, index - row_num);
721 721
          if (glp_get_col_prim(lp, index - row_num) >
722 722
              glp_get_col_ub(lp, index - row_num)) {
723 723
            negate = true;
724 724
          }
725 725
        } else {
726 726
          idx = glp_get_row_bind(lp, index);
727 727
          if (glp_get_row_prim(lp, index) > glp_get_row_ub(lp, index)) {
728 728
            negate = true;
729 729
          }
730 730
        }
731 731

	
732 732
        _dual_ray[idx] = negate ?  - 1.0 : 1.0;
733 733

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

	
764 764
        glp_btran(lp, &_dual_ray.front());
765 765

	
766 766
        for (int i = 1; i <= row_num; ++i) {
767 767
          _dual_ray[i] /= glp_get_rii(lp, i);
768 768
        }
769 769
      }
770 770
    }
771 771
    return _dual_ray[i];
772 772
  }
773 773

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

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

	
816 816
  void GlpkLp::presolver(bool b) {
817 817
    lpx_set_int_parm(lp, LPX_K_PRESOL, b ? 1 : 0);
818 818
  }
819 819

	
820 820
  void GlpkLp::messageLevel(MessageLevel m) {
821 821
    _message_level = m;
822 822
  }
823 823

	
824 824
  // GlpkMip members
825 825

	
826 826
  GlpkMip::GlpkMip()
827 827
    : LpBase(), GlpkBase(), MipSolver() {
828 828
    messageLevel(MESSAGE_NO_OUTPUT);
829 829
  }
830 830

	
831 831
  GlpkMip::GlpkMip(const GlpkMip& other)
832 832
    : LpBase(), GlpkBase(other), MipSolver() {
833 833
    messageLevel(MESSAGE_NO_OUTPUT);
834 834
  }
835 835

	
836 836
  void GlpkMip::_setColType(int i, GlpkMip::ColTypes col_type) {
837 837
    switch (col_type) {
838 838
    case INTEGER:
839 839
      glp_set_col_kind(lp, i, GLP_IV);
840 840
      break;
841 841
    case REAL:
842 842
      glp_set_col_kind(lp, i, GLP_CV);
843 843
      break;
844 844
    }
845 845
  }
846 846

	
847 847
  GlpkMip::ColTypes GlpkMip::_getColType(int i) const {
848 848
    switch (glp_get_col_kind(lp, i)) {
849 849
    case GLP_IV:
850 850
    case GLP_BV:
851 851
      return INTEGER;
852 852
    default:
853 853
      return REAL;
854 854
    }
855 855

	
856 856
  }
857 857

	
858 858
  GlpkMip::SolveExitStatus GlpkMip::_solve() {
859 859
    glp_smcp smcp;
860 860
    glp_init_smcp(&smcp);
861 861

	
862 862
    switch (_message_level) {
863 863
    case MESSAGE_NO_OUTPUT:
864 864
      smcp.msg_lev = GLP_MSG_OFF;
865 865
      break;
866 866
    case MESSAGE_ERROR_MESSAGE:
867 867
      smcp.msg_lev = GLP_MSG_ERR;
868 868
      break;
869 869
    case MESSAGE_NORMAL_OUTPUT:
870 870
      smcp.msg_lev = GLP_MSG_ON;
871 871
      break;
872 872
    case MESSAGE_FULL_OUTPUT:
873 873
      smcp.msg_lev = GLP_MSG_ALL;
874 874
      break;
875 875
    }
876 876
    smcp.meth = GLP_DUAL;
877 877

	
878 878
    if (glp_simplex(lp, &smcp) != 0) return UNSOLVED;
879 879
    if (glp_get_status(lp) != GLP_OPT) return SOLVED;
880 880

	
881 881
    glp_iocp iocp;
882 882
    glp_init_iocp(&iocp);
883 883

	
884 884
    switch (_message_level) {
885 885
    case MESSAGE_NO_OUTPUT:
886 886
      iocp.msg_lev = GLP_MSG_OFF;
887 887
      break;
888 888
    case MESSAGE_ERROR_MESSAGE:
889 889
      iocp.msg_lev = GLP_MSG_ERR;
890 890
      break;
891 891
    case MESSAGE_NORMAL_OUTPUT:
892 892
      iocp.msg_lev = GLP_MSG_ON;
893 893
      break;
894 894
    case MESSAGE_FULL_OUTPUT:
895 895
      iocp.msg_lev = GLP_MSG_ALL;
896 896
      break;
897 897
    }
898 898

	
899 899
    if (glp_intopt(lp, &iocp) != 0) return UNSOLVED;
900 900
    return SOLVED;
901 901
  }
902 902

	
903 903

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

	
935 935
  GlpkMip::Value GlpkMip::_getSol(int i) const {
936 936
    return glp_mip_col_val(lp, i);
937 937
  }
938 938

	
939 939
  GlpkMip::Value GlpkMip::_getSolValue() const {
940 940
    return glp_mip_obj_val(lp);
941 941
  }
942 942

	
943
  GlpkMip* GlpkMip::_newSolver() const { return new GlpkMip; }
944
  GlpkMip* GlpkMip::_cloneSolver() const {return new GlpkMip(*this); }
943
  GlpkMip* GlpkMip::newSolver() const { return new GlpkMip; }
944
  GlpkMip* GlpkMip::cloneSolver() const {return new GlpkMip(*this); }
945 945

	
946 946
  const char* GlpkMip::_solverName() const { return "GlpkMip"; }
947 947

	
948 948
  void GlpkMip::messageLevel(MessageLevel m) {
949 949
    _message_level = m;
950 950
  }
951 951

	
952 952
} //END OF NAMESPACE LEMON
Ignore white space 1536 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 103
  public:
104 104

	
105 105
    ///Pointer to the underlying GLPK data structure.
106 106
    LPX *lpx() {return lp;}
107 107
    ///Const pointer to the underlying GLPK data structure.
108 108
    const LPX *lpx() const {return lp;}
109 109

	
110 110
    ///Returns the constraint identifier understood by GLPK.
111 111
    int lpxRow(Row r) const { return rows(id(r)); }
112 112

	
113 113
    ///Returns the variable identifier understood by GLPK.
114 114
    int lpxCol(Col c) const { return cols(id(c)); }
115 115

	
116 116
  };
117 117

	
118 118
  /// \brief Interface for the GLPK LP solver
119 119
  ///
120 120
  /// This class implements an interface for the GLPK LP solver.
121 121
  ///\ingroup lp_group
122
  class GlpkLp : public GlpkBase, public LpSolver {
122
  class GlpkLp : public LpSolver, public GlpkBase {
123 123
  public:
124 124

	
125 125
    ///\e
126 126
    GlpkLp();
127 127
    ///\e
128 128
    GlpkLp(const GlpkLp&);
129 129

	
130
    ///\e
131
    virtual GlpkLp* cloneSolver() const;
132
    ///\e
133
    virtual GlpkLp* newSolver() const;
134

	
130 135
  private:
131 136

	
132 137
    mutable std::vector<double> _primal_ray;
133 138
    mutable std::vector<double> _dual_ray;
134 139

	
135 140
    void _clear_temporals();
136 141

	
137 142
  protected:
138 143

	
139
    virtual GlpkLp* _cloneSolver() const;
140
    virtual GlpkLp* _newSolver() const;
141

	
142 144
    virtual const char* _solverName() const;
143 145

	
144 146
    virtual SolveExitStatus _solve();
145 147
    virtual Value _getPrimal(int i) const;
146 148
    virtual Value _getDual(int i) const;
147 149

	
148 150
    virtual Value _getPrimalValue() const;
149 151

	
150 152
    virtual VarStatus _getColStatus(int i) const;
151 153
    virtual VarStatus _getRowStatus(int i) const;
152 154

	
153 155
    virtual Value _getPrimalRay(int i) const;
154 156
    virtual Value _getDualRay(int i) const;
155 157

	
156 158
    ///\todo It should be clarified
157 159
    ///
158 160
    virtual ProblemType _getPrimalType() const;
159 161
    virtual ProblemType _getDualType() const;
160 162

	
161 163
  public:
162 164

	
163 165
    ///Solve with primal simplex
164 166
    SolveExitStatus solvePrimal();
165 167

	
166 168
    ///Solve with dual simplex
167 169
    SolveExitStatus solveDual();
168 170

	
169 171
    ///Turns on or off the presolver
170 172

	
171 173
    ///Turns on (\c b is \c true) or off (\c b is \c false) the presolver
172 174
    ///
173 175
    ///The presolver is off by default.
174 176
    void presolver(bool b);
175 177

	
176 178
    ///Enum for \c messageLevel() parameter
177 179
    enum MessageLevel {
178 180
      /// no output (default value)
179 181
      MESSAGE_NO_OUTPUT = 0,
180 182
      /// error messages only
181 183
      MESSAGE_ERROR_MESSAGE = 1,
182 184
      /// normal output
183 185
      MESSAGE_NORMAL_OUTPUT = 2,
184 186
      /// full output (includes informational messages)
185 187
      MESSAGE_FULL_OUTPUT = 3
186 188
    };
187 189

	
188 190
  private:
189 191

	
190 192
    MessageLevel _message_level;
191 193

	
192 194
  public:
193 195

	
194 196
    ///Set the verbosity of the messages
195 197

	
196 198
    ///Set the verbosity of the messages
197 199
    ///
198 200
    ///\param m is the level of the messages output by the solver routines.
199 201
    void messageLevel(MessageLevel m);
200 202
  };
201 203

	
202 204
  /// \brief Interface for the GLPK MIP solver
203 205
  ///
204 206
  /// This class implements an interface for the GLPK MIP solver.
205 207
  ///\ingroup lp_group
206
  class GlpkMip : public GlpkBase, public MipSolver {
208
  class GlpkMip : public MipSolver, public GlpkBase {
207 209
  public:
208 210

	
209 211
    ///\e
210 212
    GlpkMip();
211 213
    ///\e
212 214
    GlpkMip(const GlpkMip&);
213 215

	
216
    virtual GlpkMip* cloneSolver() const;
217
    virtual GlpkMip* newSolver() const;
218

	
214 219
  protected:
215 220

	
216
    virtual GlpkMip* _cloneSolver() const;
217
    virtual GlpkMip* _newSolver() const;
218

	
219 221
    virtual const char* _solverName() const;
220 222

	
221 223
    virtual ColTypes _getColType(int col) const;
222 224
    virtual void _setColType(int col, ColTypes col_type);
223 225

	
224 226
    virtual SolveExitStatus _solve();
225 227
    virtual ProblemType _getType() const;
226 228
    virtual Value _getSol(int i) const;
227 229
    virtual Value _getSolValue() const;
228 230

	
229 231
    ///Enum for \c messageLevel() parameter
230 232
    enum MessageLevel {
231 233
      /// no output (default value)
232 234
      MESSAGE_NO_OUTPUT = 0,
233 235
      /// error messages only
234 236
      MESSAGE_ERROR_MESSAGE = 1,
235 237
      /// normal output
236 238
      MESSAGE_NORMAL_OUTPUT = 2,
237 239
      /// full output (includes informational messages)
238 240
      MESSAGE_FULL_OUTPUT = 3
239 241
    };
240 242

	
241 243
  private:
242 244

	
243 245
    MessageLevel _message_level;
244 246

	
245 247
  public:
246 248

	
247 249
    ///Set the verbosity of the messages
248 250

	
249 251
    ///Set the verbosity of the messages
250 252
    ///
251 253
    ///\param m is the level of the messages output by the solver routines.
252 254
    void messageLevel(MessageLevel m);
253 255
  };
254 256

	
255 257

	
256 258
} //END OF NAMESPACE LEMON
257 259

	
258 260
#endif //LEMON_GLPK_H
259 261

	
Ignore white space 1536 line context
... ...
@@ -153,1928 +153,1917 @@
153 153
      {
154 154
        _solver->cols.firstItem(_id);
155 155
      }
156 156
      /// Invalid constructor \& conversion
157 157
      
158 158
      /// Initialize the iterator to be invalid.
159 159
      /// \sa Invalid for more details.
160 160
      ColIt(const Invalid&) : Col(INVALID) {}
161 161
      /// Next column
162 162
      
163 163
      /// Assign the iterator to the next column.
164 164
      ///
165 165
      ColIt &operator++()
166 166
      {
167 167
        _solver->cols.nextItem(_id);
168 168
        return *this;
169 169
      }
170 170
    };
171 171

	
172 172
    /// \brief Returns the ID of the column.
173 173
    static int id(const Col& col) { return col._id; }
174 174
    /// \brief Returns the column with the given ID.
175 175
    ///
176 176
    /// \pre The argument should be a valid column ID in the LP problem.
177 177
    static Col colFromId(int id) { return Col(id); }
178 178

	
179 179
    ///Refer to a row of the LP.
180 180

	
181 181
    ///This type is used to refer to a row of the LP.
182 182
    ///
183 183
    ///Its value remains valid and correct even after the addition or erase of
184 184
    ///other rows.
185 185
    ///
186 186
    ///\note This class is similar to other Item types in LEMON, like
187 187
    ///Node and Arc types in digraph.
188 188
    class Row {
189 189
      friend class LpBase;
190 190
    protected:
191 191
      int _id;
192 192
      explicit Row(int id) : _id(id) {}
193 193
    public:
194 194
      typedef Value ExprValue;
195 195
      typedef True LpRow;
196 196
      /// Default constructor
197 197
      
198 198
      /// \warning The default constructor sets the Row to an
199 199
      /// undefined value.
200 200
      Row() {}
201 201
      /// Invalid constructor \& conversion.
202 202
      
203 203
      /// This constructor initializes the Row to be invalid.
204 204
      /// \sa Invalid for more details.      
205 205
      Row(const Invalid&) : _id(-1) {}
206 206
      /// Equality operator
207 207

	
208 208
      /// Two \ref Row "Row"s are equal if and only if they point to
209 209
      /// the same LP row or both are invalid.
210 210
      bool operator==(Row r) const  {return _id == r._id;}
211 211
      /// Inequality operator
212 212
      
213 213
      /// \sa operator==(Row r)
214 214
      ///
215 215
      bool operator!=(Row r) const  {return _id != r._id;}
216 216
      /// Artificial ordering operator.
217 217

	
218 218
      /// To allow the use of this object in std::map or similar
219 219
      /// associative container we require this.
220 220
      ///
221 221
      /// \note This operator only have to define some strict ordering of
222 222
      /// the items; this order has nothing to do with the iteration
223 223
      /// ordering of the items.
224 224
      bool operator<(Row r) const  {return _id < r._id;}
225 225
    };
226 226

	
227 227
    ///Iterator for iterate over the rows of an LP problem
228 228

	
229 229
    /// Its usage is quite simple, for example you can count the number
230 230
    /// of rows in an LP \c lp:
231 231
    ///\code
232 232
    /// int count=0;
233 233
    /// for (LpBase::RowIt c(lp); c!=INVALID; ++c) ++count;
234 234
    ///\endcode
235 235
    class RowIt : public Row {
236 236
      const LpBase *_solver;
237 237
    public:
238 238
      /// Default constructor
239 239
      
240 240
      /// \warning The default constructor sets the iterator
241 241
      /// to an undefined value.
242 242
      RowIt() {}
243 243
      /// Sets the iterator to the first Row
244 244
      
245 245
      /// Sets the iterator to the first Row.
246 246
      ///
247 247
      RowIt(const LpBase &solver) : _solver(&solver)
248 248
      {
249 249
        _solver->rows.firstItem(_id);
250 250
      }
251 251
      /// Invalid constructor \& conversion
252 252
      
253 253
      /// Initialize the iterator to be invalid.
254 254
      /// \sa Invalid for more details.
255 255
      RowIt(const Invalid&) : Row(INVALID) {}
256 256
      /// Next row
257 257
      
258 258
      /// Assign the iterator to the next row.
259 259
      ///
260 260
      RowIt &operator++()
261 261
      {
262 262
        _solver->rows.nextItem(_id);
263 263
        return *this;
264 264
      }
265 265
    };
266 266

	
267 267
    /// \brief Returns the ID of the row.
268 268
    static int id(const Row& row) { return row._id; }
269 269
    /// \brief Returns the row with the given ID.
270 270
    ///
271 271
    /// \pre The argument should be a valid row ID in the LP problem.
272 272
    static Row rowFromId(int id) { return Row(id); }
273 273

	
274 274
  public:
275 275

	
276 276
    ///Linear expression of variables and a constant component
277 277

	
278 278
    ///This data structure stores a linear expression of the variables
279 279
    ///(\ref Col "Col"s) and also has a constant component.
280 280
    ///
281 281
    ///There are several ways to access and modify the contents of this
282 282
    ///container.
283 283
    ///\code
284 284
    ///e[v]=5;
285 285
    ///e[v]+=12;
286 286
    ///e.erase(v);
287 287
    ///\endcode
288 288
    ///or you can also iterate through its elements.
289 289
    ///\code
290 290
    ///double s=0;
291 291
    ///for(LpBase::Expr::ConstCoeffIt i(e);i!=INVALID;++i)
292 292
    ///  s+=*i * primal(i);
293 293
    ///\endcode
294 294
    ///(This code computes the primal value of the expression).
295 295
    ///- Numbers (<tt>double</tt>'s)
296 296
    ///and variables (\ref Col "Col"s) directly convert to an
297 297
    ///\ref Expr and the usual linear operations are defined, so
298 298
    ///\code
299 299
    ///v+w
300 300
    ///2*v-3.12*(v-w/2)+2
301 301
    ///v*2.1+(3*v+(v*12+w+6)*3)/2
302 302
    ///\endcode
303 303
    ///are valid expressions.
304 304
    ///The usual assignment operations are also defined.
305 305
    ///\code
306 306
    ///e=v+w;
307 307
    ///e+=2*v-3.12*(v-w/2)+2;
308 308
    ///e*=3.4;
309 309
    ///e/=5;
310 310
    ///\endcode
311 311
    ///- The constant member can be set and read by dereference
312 312
    ///  operator (unary *)
313 313
    ///
314 314
    ///\code
315 315
    ///*e=12;
316 316
    ///double c=*e;
317 317
    ///\endcode
318 318
    ///
319 319
    ///\sa Constr
320 320
    class Expr {
321 321
      friend class LpBase;
322 322
    public:
323 323
      /// The key type of the expression
324 324
      typedef LpBase::Col Key;
325 325
      /// The value type of the expression
326 326
      typedef LpBase::Value Value;
327 327

	
328 328
    protected:
329 329
      Value const_comp;
330 330
      std::map<int, Value> comps;
331 331

	
332 332
    public:
333 333
      typedef True SolverExpr;
334 334
      /// Default constructor
335 335
      
336 336
      /// Construct an empty expression, the coefficients and
337 337
      /// the constant component are initialized to zero.
338 338
      Expr() : const_comp(0) {}
339 339
      /// Construct an expression from a column
340 340

	
341 341
      /// Construct an expression, which has a term with \c c variable
342 342
      /// and 1.0 coefficient.
343 343
      Expr(const Col &c) : const_comp(0) {
344 344
        typedef std::map<int, Value>::value_type pair_type;
345 345
        comps.insert(pair_type(id(c), 1));
346 346
      }
347 347
      /// Construct an expression from a constant
348 348

	
349 349
      /// Construct an expression, which's constant component is \c v.
350 350
      ///
351 351
      Expr(const Value &v) : const_comp(v) {}
352 352
      /// Returns the coefficient of the column
353 353
      Value operator[](const Col& c) const {
354 354
        std::map<int, Value>::const_iterator it=comps.find(id(c));
355 355
        if (it != comps.end()) {
356 356
          return it->second;
357 357
        } else {
358 358
          return 0;
359 359
        }
360 360
      }
361 361
      /// Returns the coefficient of the column
362 362
      Value& operator[](const Col& c) {
363 363
        return comps[id(c)];
364 364
      }
365 365
      /// Sets the coefficient of the column
366 366
      void set(const Col &c, const Value &v) {
367 367
        if (v != 0.0) {
368 368
          typedef std::map<int, Value>::value_type pair_type;
369 369
          comps.insert(pair_type(id(c), v));
370 370
        } else {
371 371
          comps.erase(id(c));
372 372
        }
373 373
      }
374 374
      /// Returns the constant component of the expression
375 375
      Value& operator*() { return const_comp; }
376 376
      /// Returns the constant component of the expression
377 377
      const Value& operator*() const { return const_comp; }
378 378
      /// \brief Removes the coefficients which's absolute value does
379 379
      /// not exceed \c epsilon. It also sets to zero the constant
380 380
      /// component, if it does not exceed epsilon in absolute value.
381 381
      void simplify(Value epsilon = 0.0) {
382 382
        std::map<int, Value>::iterator it=comps.begin();
383 383
        while (it != comps.end()) {
384 384
          std::map<int, Value>::iterator jt=it;
385 385
          ++jt;
386 386
          if (std::fabs((*it).second) <= epsilon) comps.erase(it);
387 387
          it=jt;
388 388
        }
389 389
        if (std::fabs(const_comp) <= epsilon) const_comp = 0;
390 390
      }
391 391

	
392 392
      void simplify(Value epsilon = 0.0) const {
393 393
        const_cast<Expr*>(this)->simplify(epsilon);
394 394
      }
395 395

	
396 396
      ///Sets all coefficients and the constant component to 0.
397 397
      void clear() {
398 398
        comps.clear();
399 399
        const_comp=0;
400 400
      }
401 401

	
402 402
      ///Compound assignment
403 403
      Expr &operator+=(const Expr &e) {
404 404
        for (std::map<int, Value>::const_iterator it=e.comps.begin();
405 405
             it!=e.comps.end(); ++it)
406 406
          comps[it->first]+=it->second;
407 407
        const_comp+=e.const_comp;
408 408
        return *this;
409 409
      }
410 410
      ///Compound assignment
411 411
      Expr &operator-=(const Expr &e) {
412 412
        for (std::map<int, Value>::const_iterator it=e.comps.begin();
413 413
             it!=e.comps.end(); ++it)
414 414
          comps[it->first]-=it->second;
415 415
        const_comp-=e.const_comp;
416 416
        return *this;
417 417
      }
418 418
      ///Multiply with a constant
419 419
      Expr &operator*=(const Value &v) {
420 420
        for (std::map<int, Value>::iterator it=comps.begin();
421 421
             it!=comps.end(); ++it)
422 422
          it->second*=v;
423 423
        const_comp*=v;
424 424
        return *this;
425 425
      }
426 426
      ///Division with a constant
427 427
      Expr &operator/=(const Value &c) {
428 428
        for (std::map<int, Value>::iterator it=comps.begin();
429 429
             it!=comps.end(); ++it)
430 430
          it->second/=c;
431 431
        const_comp/=c;
432 432
        return *this;
433 433
      }
434 434

	
435 435
      ///Iterator over the expression
436 436
      
437 437
      ///The iterator iterates over the terms of the expression. 
438 438
      /// 
439 439
      ///\code
440 440
      ///double s=0;
441 441
      ///for(LpBase::Expr::CoeffIt i(e);i!=INVALID;++i)
442 442
      ///  s+= *i * primal(i);
443 443
      ///\endcode
444 444
      class CoeffIt {
445 445
      private:
446 446

	
447 447
        std::map<int, Value>::iterator _it, _end;
448 448

	
449 449
      public:
450 450

	
451 451
        /// Sets the iterator to the first term
452 452
        
453 453
        /// Sets the iterator to the first term of the expression.
454 454
        ///
455 455
        CoeffIt(Expr& e)
456 456
          : _it(e.comps.begin()), _end(e.comps.end()){}
457 457

	
458 458
        /// Convert the iterator to the column of the term
459 459
        operator Col() const {
460 460
          return colFromId(_it->first);
461 461
        }
462 462

	
463 463
        /// Returns the coefficient of the term
464 464
        Value& operator*() { return _it->second; }
465 465

	
466 466
        /// Returns the coefficient of the term
467 467
        const Value& operator*() const { return _it->second; }
468 468
        /// Next term
469 469
        
470 470
        /// Assign the iterator to the next term.
471 471
        ///
472 472
        CoeffIt& operator++() { ++_it; return *this; }
473 473

	
474 474
        /// Equality operator
475 475
        bool operator==(Invalid) const { return _it == _end; }
476 476
        /// Inequality operator
477 477
        bool operator!=(Invalid) const { return _it != _end; }
478 478
      };
479 479

	
480 480
      /// Const iterator over the expression
481 481
      
482 482
      ///The iterator iterates over the terms of the expression. 
483 483
      /// 
484 484
      ///\code
485 485
      ///double s=0;
486 486
      ///for(LpBase::Expr::ConstCoeffIt i(e);i!=INVALID;++i)
487 487
      ///  s+=*i * primal(i);
488 488
      ///\endcode
489 489
      class ConstCoeffIt {
490 490
      private:
491 491

	
492 492
        std::map<int, Value>::const_iterator _it, _end;
493 493

	
494 494
      public:
495 495

	
496 496
        /// Sets the iterator to the first term
497 497
        
498 498
        /// Sets the iterator to the first term of the expression.
499 499
        ///
500 500
        ConstCoeffIt(const Expr& e)
501 501
          : _it(e.comps.begin()), _end(e.comps.end()){}
502 502

	
503 503
        /// Convert the iterator to the column of the term
504 504
        operator Col() const {
505 505
          return colFromId(_it->first);
506 506
        }
507 507

	
508 508
        /// Returns the coefficient of the term
509 509
        const Value& operator*() const { return _it->second; }
510 510

	
511 511
        /// Next term
512 512
        
513 513
        /// Assign the iterator to the next term.
514 514
        ///
515 515
        ConstCoeffIt& operator++() { ++_it; return *this; }
516 516

	
517 517
        /// Equality operator
518 518
        bool operator==(Invalid) const { return _it == _end; }
519 519
        /// Inequality operator
520 520
        bool operator!=(Invalid) const { return _it != _end; }
521 521
      };
522 522

	
523 523
    };
524 524

	
525 525
    ///Linear constraint
526 526

	
527 527
    ///This data stucture represents a linear constraint in the LP.
528 528
    ///Basically it is a linear expression with a lower or an upper bound
529 529
    ///(or both). These parts of the constraint can be obtained by the member
530 530
    ///functions \ref expr(), \ref lowerBound() and \ref upperBound(),
531 531
    ///respectively.
532 532
    ///There are two ways to construct a constraint.
533 533
    ///- You can set the linear expression and the bounds directly
534 534
    ///  by the functions above.
535 535
    ///- The operators <tt>\<=</tt>, <tt>==</tt> and  <tt>\>=</tt>
536 536
    ///  are defined between expressions, or even between constraints whenever
537 537
    ///  it makes sense. Therefore if \c e and \c f are linear expressions and
538 538
    ///  \c s and \c t are numbers, then the followings are valid expressions
539 539
    ///  and thus they can be used directly e.g. in \ref addRow() whenever
540 540
    ///  it makes sense.
541 541
    ///\code
542 542
    ///  e<=s
543 543
    ///  e<=f
544 544
    ///  e==f
545 545
    ///  s<=e<=t
546 546
    ///  e>=t
547 547
    ///\endcode
548 548
    ///\warning The validity of a constraint is checked only at run
549 549
    ///time, so e.g. \ref addRow(<tt>x[1]\<=x[2]<=5</tt>) will
550 550
    ///compile, but will fail an assertion.
551 551
    class Constr
552 552
    {
553 553
    public:
554 554
      typedef LpBase::Expr Expr;
555 555
      typedef Expr::Key Key;
556 556
      typedef Expr::Value Value;
557 557

	
558 558
    protected:
559 559
      Expr _expr;
560 560
      Value _lb,_ub;
561 561
    public:
562 562
      ///\e
563 563
      Constr() : _expr(), _lb(NaN), _ub(NaN) {}
564 564
      ///\e
565 565
      Constr(Value lb, const Expr &e, Value ub) :
566 566
        _expr(e), _lb(lb), _ub(ub) {}
567 567
      Constr(const Expr &e) :
568 568
        _expr(e), _lb(NaN), _ub(NaN) {}
569 569
      ///\e
570 570
      void clear()
571 571
      {
572 572
        _expr.clear();
573 573
        _lb=_ub=NaN;
574 574
      }
575 575

	
576 576
      ///Reference to the linear expression
577 577
      Expr &expr() { return _expr; }
578 578
      ///Cont reference to the linear expression
579 579
      const Expr &expr() const { return _expr; }
580 580
      ///Reference to the lower bound.
581 581

	
582 582
      ///\return
583 583
      ///- \ref INF "INF": the constraint is lower unbounded.
584 584
      ///- \ref NaN "NaN": lower bound has not been set.
585 585
      ///- finite number: the lower bound
586 586
      Value &lowerBound() { return _lb; }
587 587
      ///The const version of \ref lowerBound()
588 588
      const Value &lowerBound() const { return _lb; }
589 589
      ///Reference to the upper bound.
590 590

	
591 591
      ///\return
592 592
      ///- \ref INF "INF": the constraint is upper unbounded.
593 593
      ///- \ref NaN "NaN": upper bound has not been set.
594 594
      ///- finite number: the upper bound
595 595
      Value &upperBound() { return _ub; }
596 596
      ///The const version of \ref upperBound()
597 597
      const Value &upperBound() const { return _ub; }
598 598
      ///Is the constraint lower bounded?
599 599
      bool lowerBounded() const {
600 600
        return _lb != -INF && !isNaN(_lb);
601 601
      }
602 602
      ///Is the constraint upper bounded?
603 603
      bool upperBounded() const {
604 604
        return _ub != INF && !isNaN(_ub);
605 605
      }
606 606

	
607 607
    };
608 608

	
609 609
    ///Linear expression of rows
610 610

	
611 611
    ///This data structure represents a column of the matrix,
612 612
    ///thas is it strores a linear expression of the dual variables
613 613
    ///(\ref Row "Row"s).
614 614
    ///
615 615
    ///There are several ways to access and modify the contents of this
616 616
    ///container.
617 617
    ///\code
618 618
    ///e[v]=5;
619 619
    ///e[v]+=12;
620 620
    ///e.erase(v);
621 621
    ///\endcode
622 622
    ///or you can also iterate through its elements.
623 623
    ///\code
624 624
    ///double s=0;
625 625
    ///for(LpBase::DualExpr::ConstCoeffIt i(e);i!=INVALID;++i)
626 626
    ///  s+=*i;
627 627
    ///\endcode
628 628
    ///(This code computes the sum of all coefficients).
629 629
    ///- Numbers (<tt>double</tt>'s)
630 630
    ///and variables (\ref Row "Row"s) directly convert to an
631 631
    ///\ref DualExpr and the usual linear operations are defined, so
632 632
    ///\code
633 633
    ///v+w
634 634
    ///2*v-3.12*(v-w/2)
635 635
    ///v*2.1+(3*v+(v*12+w)*3)/2
636 636
    ///\endcode
637 637
    ///are valid \ref DualExpr dual expressions.
638 638
    ///The usual assignment operations are also defined.
639 639
    ///\code
640 640
    ///e=v+w;
641 641
    ///e+=2*v-3.12*(v-w/2);
642 642
    ///e*=3.4;
643 643
    ///e/=5;
644 644
    ///\endcode
645 645
    ///
646 646
    ///\sa Expr
647 647
    class DualExpr {
648 648
      friend class LpBase;
649 649
    public:
650 650
      /// The key type of the expression
651 651
      typedef LpBase::Row Key;
652 652
      /// The value type of the expression
653 653
      typedef LpBase::Value Value;
654 654

	
655 655
    protected:
656 656
      std::map<int, Value> comps;
657 657

	
658 658
    public:
659 659
      typedef True SolverExpr;
660 660
      /// Default constructor
661 661
      
662 662
      /// Construct an empty expression, the coefficients are
663 663
      /// initialized to zero.
664 664
      DualExpr() {}
665 665
      /// Construct an expression from a row
666 666

	
667 667
      /// Construct an expression, which has a term with \c r dual
668 668
      /// variable and 1.0 coefficient.
669 669
      DualExpr(const Row &r) {
670 670
        typedef std::map<int, Value>::value_type pair_type;
671 671
        comps.insert(pair_type(id(r), 1));
672 672
      }
673 673
      /// Returns the coefficient of the row
674 674
      Value operator[](const Row& r) const {
675 675
        std::map<int, Value>::const_iterator it = comps.find(id(r));
676 676
        if (it != comps.end()) {
677 677
          return it->second;
678 678
        } else {
679 679
          return 0;
680 680
        }
681 681
      }
682 682
      /// Returns the coefficient of the row
683 683
      Value& operator[](const Row& r) {
684 684
        return comps[id(r)];
685 685
      }
686 686
      /// Sets the coefficient of the row
687 687
      void set(const Row &r, const Value &v) {
688 688
        if (v != 0.0) {
689 689
          typedef std::map<int, Value>::value_type pair_type;
690 690
          comps.insert(pair_type(id(r), v));
691 691
        } else {
692 692
          comps.erase(id(r));
693 693
        }
694 694
      }
695 695
      /// \brief Removes the coefficients which's absolute value does
696 696
      /// not exceed \c epsilon. 
697 697
      void simplify(Value epsilon = 0.0) {
698 698
        std::map<int, Value>::iterator it=comps.begin();
699 699
        while (it != comps.end()) {
700 700
          std::map<int, Value>::iterator jt=it;
701 701
          ++jt;
702 702
          if (std::fabs((*it).second) <= epsilon) comps.erase(it);
703 703
          it=jt;
704 704
        }
705 705
      }
706 706

	
707 707
      void simplify(Value epsilon = 0.0) const {
708 708
        const_cast<DualExpr*>(this)->simplify(epsilon);
709 709
      }
710 710

	
711 711
      ///Sets all coefficients to 0.
712 712
      void clear() {
713 713
        comps.clear();
714 714
      }
715 715
      ///Compound assignment
716 716
      DualExpr &operator+=(const DualExpr &e) {
717 717
        for (std::map<int, Value>::const_iterator it=e.comps.begin();
718 718
             it!=e.comps.end(); ++it)
719 719
          comps[it->first]+=it->second;
720 720
        return *this;
721 721
      }
722 722
      ///Compound assignment
723 723
      DualExpr &operator-=(const DualExpr &e) {
724 724
        for (std::map<int, Value>::const_iterator it=e.comps.begin();
725 725
             it!=e.comps.end(); ++it)
726 726
          comps[it->first]-=it->second;
727 727
        return *this;
728 728
      }
729 729
      ///Multiply with a constant
730 730
      DualExpr &operator*=(const Value &v) {
731 731
        for (std::map<int, Value>::iterator it=comps.begin();
732 732
             it!=comps.end(); ++it)
733 733
          it->second*=v;
734 734
        return *this;
735 735
      }
736 736
      ///Division with a constant
737 737
      DualExpr &operator/=(const Value &v) {
738 738
        for (std::map<int, Value>::iterator it=comps.begin();
739 739
             it!=comps.end(); ++it)
740 740
          it->second/=v;
741 741
        return *this;
742 742
      }
743 743

	
744 744
      ///Iterator over the expression
745 745
      
746 746
      ///The iterator iterates over the terms of the expression. 
747 747
      /// 
748 748
      ///\code
749 749
      ///double s=0;
750 750
      ///for(LpBase::DualExpr::CoeffIt i(e);i!=INVALID;++i)
751 751
      ///  s+= *i * dual(i);
752 752
      ///\endcode
753 753
      class CoeffIt {
754 754
      private:
755 755

	
756 756
        std::map<int, Value>::iterator _it, _end;
757 757

	
758 758
      public:
759 759

	
760 760
        /// Sets the iterator to the first term
761 761
        
762 762
        /// Sets the iterator to the first term of the expression.
763 763
        ///
764 764
        CoeffIt(DualExpr& e)
765 765
          : _it(e.comps.begin()), _end(e.comps.end()){}
766 766

	
767 767
        /// Convert the iterator to the row of the term
768 768
        operator Row() const {
769 769
          return rowFromId(_it->first);
770 770
        }
771 771

	
772 772
        /// Returns the coefficient of the term
773 773
        Value& operator*() { return _it->second; }
774 774

	
775 775
        /// Returns the coefficient of the term
776 776
        const Value& operator*() const { return _it->second; }
777 777

	
778 778
        /// Next term
779 779
        
780 780
        /// Assign the iterator to the next term.
781 781
        ///
782 782
        CoeffIt& operator++() { ++_it; return *this; }
783 783

	
784 784
        /// Equality operator
785 785
        bool operator==(Invalid) const { return _it == _end; }
786 786
        /// Inequality operator
787 787
        bool operator!=(Invalid) const { return _it != _end; }
788 788
      };
789 789

	
790 790
      ///Iterator over the expression
791 791
      
792 792
      ///The iterator iterates over the terms of the expression. 
793 793
      /// 
794 794
      ///\code
795 795
      ///double s=0;
796 796
      ///for(LpBase::DualExpr::ConstCoeffIt i(e);i!=INVALID;++i)
797 797
      ///  s+= *i * dual(i);
798 798
      ///\endcode
799 799
      class ConstCoeffIt {
800 800
      private:
801 801

	
802 802
        std::map<int, Value>::const_iterator _it, _end;
803 803

	
804 804
      public:
805 805

	
806 806
        /// Sets the iterator to the first term
807 807
        
808 808
        /// Sets the iterator to the first term of the expression.
809 809
        ///
810 810
        ConstCoeffIt(const DualExpr& e)
811 811
          : _it(e.comps.begin()), _end(e.comps.end()){}
812 812

	
813 813
        /// Convert the iterator to the row of the term
814 814
        operator Row() const {
815 815
          return rowFromId(_it->first);
816 816
        }
817 817

	
818 818
        /// Returns the coefficient of the term
819 819
        const Value& operator*() const { return _it->second; }
820 820

	
821 821
        /// Next term
822 822
        
823 823
        /// Assign the iterator to the next term.
824 824
        ///
825 825
        ConstCoeffIt& operator++() { ++_it; return *this; }
826 826

	
827 827
        /// Equality operator
828 828
        bool operator==(Invalid) const { return _it == _end; }
829 829
        /// Inequality operator
830 830
        bool operator!=(Invalid) const { return _it != _end; }
831 831
      };
832 832
    };
833 833

	
834 834

	
835 835
  protected:
836 836

	
837 837
    class InsertIterator {
838 838
    private:
839 839

	
840 840
      std::map<int, Value>& _host;
841 841
      const _solver_bits::VarIndex& _index;
842 842

	
843 843
    public:
844 844

	
845 845
      typedef std::output_iterator_tag iterator_category;
846 846
      typedef void difference_type;
847 847
      typedef void value_type;
848 848
      typedef void reference;
849 849
      typedef void pointer;
850 850

	
851 851
      InsertIterator(std::map<int, Value>& host,
852 852
                   const _solver_bits::VarIndex& index)
853 853
        : _host(host), _index(index) {}
854 854

	
855 855
      InsertIterator& operator=(const std::pair<int, Value>& value) {
856 856
        typedef std::map<int, Value>::value_type pair_type;
857 857
        _host.insert(pair_type(_index[value.first], value.second));
858 858
        return *this;
859 859
      }
860 860

	
861 861
      InsertIterator& operator*() { return *this; }
862 862
      InsertIterator& operator++() { return *this; }
863 863
      InsertIterator operator++(int) { return *this; }
864 864

	
865 865
    };
866 866

	
867 867
    class ExprIterator {
868 868
    private:
869 869
      std::map<int, Value>::const_iterator _host_it;
870 870
      const _solver_bits::VarIndex& _index;
871 871
    public:
872 872

	
873 873
      typedef std::bidirectional_iterator_tag iterator_category;
874 874
      typedef std::ptrdiff_t difference_type;
875 875
      typedef const std::pair<int, Value> value_type;
876 876
      typedef value_type reference;
877 877

	
878 878
      class pointer {
879 879
      public:
880 880
        pointer(value_type& _value) : value(_value) {}
881 881
        value_type* operator->() { return &value; }
882 882
      private:
883 883
        value_type value;
884 884
      };
885 885

	
886 886
      ExprIterator(const std::map<int, Value>::const_iterator& host_it,
887 887
                   const _solver_bits::VarIndex& index)
888 888
        : _host_it(host_it), _index(index) {}
889 889

	
890 890
      reference operator*() {
891 891
        return std::make_pair(_index(_host_it->first), _host_it->second);
892 892
      }
893 893

	
894 894
      pointer operator->() {
895 895
        return pointer(operator*());
896 896
      }
897 897

	
898 898
      ExprIterator& operator++() { ++_host_it; return *this; }
899 899
      ExprIterator operator++(int) {
900 900
        ExprIterator tmp(*this); ++_host_it; return tmp;
901 901
      }
902 902

	
903 903
      ExprIterator& operator--() { --_host_it; return *this; }
904 904
      ExprIterator operator--(int) {
905 905
        ExprIterator tmp(*this); --_host_it; return tmp;
906 906
      }
907 907

	
908 908
      bool operator==(const ExprIterator& it) const {
909 909
        return _host_it == it._host_it;
910 910
      }
911 911

	
912 912
      bool operator!=(const ExprIterator& it) const {
913 913
        return _host_it != it._host_it;
914 914
      }
915 915

	
916 916
    };
917 917

	
918 918
  protected:
919 919

	
920 920
    //Abstract virtual functions
921
    virtual LpBase* _newSolver() const = 0;
922
    virtual LpBase* _cloneSolver() const = 0;
923 921

	
924 922
    virtual int _addColId(int col) { return cols.addIndex(col); }
925 923
    virtual int _addRowId(int row) { return rows.addIndex(row); }
926 924

	
927 925
    virtual void _eraseColId(int col) { cols.eraseIndex(col); }
928 926
    virtual void _eraseRowId(int row) { rows.eraseIndex(row); }
929 927

	
930 928
    virtual int _addCol() = 0;
931 929
    virtual int _addRow() = 0;
932 930

	
933 931
    virtual void _eraseCol(int col) = 0;
934 932
    virtual void _eraseRow(int row) = 0;
935 933

	
936 934
    virtual void _getColName(int col, std::string& name) const = 0;
937 935
    virtual void _setColName(int col, const std::string& name) = 0;
938 936
    virtual int _colByName(const std::string& name) const = 0;
939 937

	
940 938
    virtual void _getRowName(int row, std::string& name) const = 0;
941 939
    virtual void _setRowName(int row, const std::string& name) = 0;
942 940
    virtual int _rowByName(const std::string& name) const = 0;
943 941

	
944 942
    virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e) = 0;
945 943
    virtual void _getRowCoeffs(int i, InsertIterator b) const = 0;
946 944

	
947 945
    virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e) = 0;
948 946
    virtual void _getColCoeffs(int i, InsertIterator b) const = 0;
949 947

	
950 948
    virtual void _setCoeff(int row, int col, Value value) = 0;
951 949
    virtual Value _getCoeff(int row, int col) const = 0;
952 950

	
953 951
    virtual void _setColLowerBound(int i, Value value) = 0;
954 952
    virtual Value _getColLowerBound(int i) const = 0;
955 953

	
956 954
    virtual void _setColUpperBound(int i, Value value) = 0;
957 955
    virtual Value _getColUpperBound(int i) const = 0;
958 956

	
959 957
    virtual void _setRowLowerBound(int i, Value value) = 0;
960 958
    virtual Value _getRowLowerBound(int i) const = 0;
961 959

	
962 960
    virtual void _setRowUpperBound(int i, Value value) = 0;
963 961
    virtual Value _getRowUpperBound(int i) const = 0;
964 962

	
965 963
    virtual void _setObjCoeffs(ExprIterator b, ExprIterator e) = 0;
966 964
    virtual void _getObjCoeffs(InsertIterator b) const = 0;
967 965

	
968 966
    virtual void _setObjCoeff(int i, Value obj_coef) = 0;
969 967
    virtual Value _getObjCoeff(int i) const = 0;
970 968

	
971 969
    virtual void _setSense(Sense) = 0;
972 970
    virtual Sense _getSense() const = 0;
973 971

	
974 972
    virtual void _clear() = 0;
975 973

	
976 974
    virtual const char* _solverName() const = 0;
977 975

	
978 976
    //Own protected stuff
979 977

	
980 978
    //Constant component of the objective function
981 979
    Value obj_const_comp;
982 980

	
983 981
    LpBase() : rows(), cols(), obj_const_comp(0) {}
984 982

	
985 983
  public:
986 984

	
987 985
    /// Virtual destructor
988 986
    virtual ~LpBase() {}
989 987

	
990
    ///Creates a new LP problem
991
    LpBase* newSolver() {return _newSolver();}
992
    ///Makes a copy of the LP problem
993
    LpBase* cloneSolver() {return _cloneSolver();}
994

	
995 988
    ///Gives back the name of the solver.
996 989
    const char* solverName() const {return _solverName();}
997 990

	
998 991
    ///\name Build up and modify the LP
999 992

	
1000 993
    ///@{
1001 994

	
1002 995
    ///Add a new empty column (i.e a new variable) to the LP
1003 996
    Col addCol() { Col c; c._id = _addColId(_addCol()); return c;}
1004 997

	
1005 998
    ///\brief Adds several new columns (i.e variables) at once
1006 999
    ///
1007 1000
    ///This magic function takes a container as its argument and fills
1008 1001
    ///its elements with new columns (i.e. variables)
1009 1002
    ///\param t can be
1010 1003
    ///- a standard STL compatible iterable container with
1011 1004
    ///\ref Col as its \c values_type like
1012 1005
    ///\code
1013 1006
    ///std::vector<LpBase::Col>
1014 1007
    ///std::list<LpBase::Col>
1015 1008
    ///\endcode
1016 1009
    ///- a standard STL compatible iterable container with
1017 1010
    ///\ref Col as its \c mapped_type like
1018 1011
    ///\code
1019 1012
    ///std::map<AnyType,LpBase::Col>
1020 1013
    ///\endcode
1021 1014
    ///- an iterable lemon \ref concepts::WriteMap "write map" like
1022 1015
    ///\code
1023 1016
    ///ListGraph::NodeMap<LpBase::Col>
1024 1017
    ///ListGraph::ArcMap<LpBase::Col>
1025 1018
    ///\endcode
1026 1019
    ///\return The number of the created column.
1027 1020
#ifdef DOXYGEN
1028 1021
    template<class T>
1029 1022
    int addColSet(T &t) { return 0;}
1030 1023
#else
1031 1024
    template<class T>
1032 1025
    typename enable_if<typename T::value_type::LpCol,int>::type
1033 1026
    addColSet(T &t,dummy<0> = 0) {
1034 1027
      int s=0;
1035 1028
      for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addCol();s++;}
1036 1029
      return s;
1037 1030
    }
1038 1031
    template<class T>
1039 1032
    typename enable_if<typename T::value_type::second_type::LpCol,
1040 1033
                       int>::type
1041 1034
    addColSet(T &t,dummy<1> = 1) {
1042 1035
      int s=0;
1043 1036
      for(typename T::iterator i=t.begin();i!=t.end();++i) {
1044 1037
        i->second=addCol();
1045 1038
        s++;
1046 1039
      }
1047 1040
      return s;
1048 1041
    }
1049 1042
    template<class T>
1050 1043
    typename enable_if<typename T::MapIt::Value::LpCol,
1051 1044
                       int>::type
1052 1045
    addColSet(T &t,dummy<2> = 2) {
1053 1046
      int s=0;
1054 1047
      for(typename T::MapIt i(t); i!=INVALID; ++i)
1055 1048
        {
1056 1049
          i.set(addCol());
1057 1050
          s++;
1058 1051
        }
1059 1052
      return s;
1060 1053
    }
1061 1054
#endif
1062 1055

	
1063 1056
    ///Set a column (i.e a dual constraint) of the LP
1064 1057

	
1065 1058
    ///\param c is the column to be modified
1066 1059
    ///\param e is a dual linear expression (see \ref DualExpr)
1067 1060
    ///a better one.
1068 1061
    void col(Col c, const DualExpr &e) {
1069 1062
      e.simplify();
1070 1063
      _setColCoeffs(cols(id(c)), ExprIterator(e.comps.begin(), rows),
1071 1064
                    ExprIterator(e.comps.end(), rows));
1072 1065
    }
1073 1066

	
1074 1067
    ///Get a column (i.e a dual constraint) of the LP
1075 1068

	
1076 1069
    ///\param c is the column to get
1077 1070
    ///\return the dual expression associated to the column
1078 1071
    DualExpr col(Col c) const {
1079 1072
      DualExpr e;
1080 1073
      _getColCoeffs(cols(id(c)), InsertIterator(e.comps, rows));
1081 1074
      return e;
1082 1075
    }
1083 1076

	
1084 1077
    ///Add a new column to the LP
1085 1078

	
1086 1079
    ///\param e is a dual linear expression (see \ref DualExpr)
1087 1080
    ///\param o is the corresponding component of the objective
1088 1081
    ///function. It is 0 by default.
1089 1082
    ///\return The created column.
1090 1083
    Col addCol(const DualExpr &e, Value o = 0) {
1091 1084
      Col c=addCol();
1092 1085
      col(c,e);
1093 1086
      objCoeff(c,o);
1094 1087
      return c;
1095 1088
    }
1096 1089

	
1097 1090
    ///Add a new empty row (i.e a new constraint) to the LP
1098 1091

	
1099 1092
    ///This function adds a new empty row (i.e a new constraint) to the LP.
1100 1093
    ///\return The created row
1101 1094
    Row addRow() { Row r; r._id = _addRowId(_addRow()); return r;}
1102 1095

	
1103 1096
    ///\brief Add several new rows (i.e constraints) at once
1104 1097
    ///
1105 1098
    ///This magic function takes a container as its argument and fills
1106 1099
    ///its elements with new row (i.e. variables)
1107 1100
    ///\param t can be
1108 1101
    ///- a standard STL compatible iterable container with
1109 1102
    ///\ref Row as its \c values_type like
1110 1103
    ///\code
1111 1104
    ///std::vector<LpBase::Row>
1112 1105
    ///std::list<LpBase::Row>
1113 1106
    ///\endcode
1114 1107
    ///- a standard STL compatible iterable container with
1115 1108
    ///\ref Row as its \c mapped_type like
1116 1109
    ///\code
1117 1110
    ///std::map<AnyType,LpBase::Row>
1118 1111
    ///\endcode
1119 1112
    ///- an iterable lemon \ref concepts::WriteMap "write map" like
1120 1113
    ///\code
1121 1114
    ///ListGraph::NodeMap<LpBase::Row>
1122 1115
    ///ListGraph::ArcMap<LpBase::Row>
1123 1116
    ///\endcode
1124 1117
    ///\return The number of rows created.
1125 1118
#ifdef DOXYGEN
1126 1119
    template<class T>
1127 1120
    int addRowSet(T &t) { return 0;}
1128 1121
#else
1129 1122
    template<class T>
1130 1123
    typename enable_if<typename T::value_type::LpRow,int>::type
1131 1124
    addRowSet(T &t, dummy<0> = 0) {
1132 1125
      int s=0;
1133 1126
      for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addRow();s++;}
1134 1127
      return s;
1135 1128
    }
1136 1129
    template<class T>
1137 1130
    typename enable_if<typename T::value_type::second_type::LpRow, int>::type
1138 1131
    addRowSet(T &t, dummy<1> = 1) {
1139 1132
      int s=0;
1140 1133
      for(typename T::iterator i=t.begin();i!=t.end();++i) {
1141 1134
        i->second=addRow();
1142 1135
        s++;
1143 1136
      }
1144 1137
      return s;
1145 1138
    }
1146 1139
    template<class T>
1147 1140
    typename enable_if<typename T::MapIt::Value::LpRow, int>::type
1148 1141
    addRowSet(T &t, dummy<2> = 2) {
1149 1142
      int s=0;
1150 1143
      for(typename T::MapIt i(t); i!=INVALID; ++i)
1151 1144
        {
1152 1145
          i.set(addRow());
1153 1146
          s++;
1154 1147
        }
1155 1148
      return s;
1156 1149
    }
1157 1150
#endif
1158 1151

	
1159 1152
    ///Set a row (i.e a constraint) of the LP
1160 1153

	
1161 1154
    ///\param r is the row to be modified
1162 1155
    ///\param l is lower bound (-\ref INF means no bound)
1163 1156
    ///\param e is a linear expression (see \ref Expr)
1164 1157
    ///\param u is the upper bound (\ref INF means no bound)
1165 1158
    void row(Row r, Value l, const Expr &e, Value u) {
1166 1159
      e.simplify();
1167 1160
      _setRowCoeffs(rows(id(r)), ExprIterator(e.comps.begin(), cols),
1168 1161
                    ExprIterator(e.comps.end(), cols));
1169 1162
      _setRowLowerBound(rows(id(r)),l - *e);
1170 1163
      _setRowUpperBound(rows(id(r)),u - *e);
1171 1164
    }
1172 1165

	
1173 1166
    ///Set a row (i.e a constraint) of the LP
1174 1167

	
1175 1168
    ///\param r is the row to be modified
1176 1169
    ///\param c is a linear expression (see \ref Constr)
1177 1170
    void row(Row r, const Constr &c) {
1178 1171
      row(r, c.lowerBounded()?c.lowerBound():-INF,
1179 1172
          c.expr(), c.upperBounded()?c.upperBound():INF);
1180 1173
    }
1181 1174

	
1182 1175

	
1183 1176
    ///Get a row (i.e a constraint) of the LP
1184 1177

	
1185 1178
    ///\param r is the row to get
1186 1179
    ///\return the expression associated to the row
1187 1180
    Expr row(Row r) const {
1188 1181
      Expr e;
1189 1182
      _getRowCoeffs(rows(id(r)), InsertIterator(e.comps, cols));
1190 1183
      return e;
1191 1184
    }
1192 1185

	
1193 1186
    ///Add a new row (i.e a new constraint) to the LP
1194 1187

	
1195 1188
    ///\param l is the lower bound (-\ref INF means no bound)
1196 1189
    ///\param e is a linear expression (see \ref Expr)
1197 1190
    ///\param u is the upper bound (\ref INF means no bound)
1198 1191
    ///\return The created row.
1199 1192
    Row addRow(Value l,const Expr &e, Value u) {
1200 1193
      Row r=addRow();
1201 1194
      row(r,l,e,u);
1202 1195
      return r;
1203 1196
    }
1204 1197

	
1205 1198
    ///Add a new row (i.e a new constraint) to the LP
1206 1199

	
1207 1200
    ///\param c is a linear expression (see \ref Constr)
1208 1201
    ///\return The created row.
1209 1202
    Row addRow(const Constr &c) {
1210 1203
      Row r=addRow();
1211 1204
      row(r,c);
1212 1205
      return r;
1213 1206
    }
1214 1207
    ///Erase a column (i.e a variable) from the LP
1215 1208

	
1216 1209
    ///\param c is the column to be deleted
1217 1210
    void erase(Col c) {
1218 1211
      _eraseCol(cols(id(c)));
1219 1212
      _eraseColId(cols(id(c)));
1220 1213
    }
1221 1214
    ///Erase a row (i.e a constraint) from the LP
1222 1215

	
1223 1216
    ///\param r is the row to be deleted
1224 1217
    void erase(Row r) {
1225 1218
      _eraseRow(rows(id(r)));
1226 1219
      _eraseRowId(rows(id(r)));
1227 1220
    }
1228 1221

	
1229 1222
    /// Get the name of a column
1230 1223

	
1231 1224
    ///\param c is the coresponding column
1232 1225
    ///\return The name of the colunm
1233 1226
    std::string colName(Col c) const {
1234 1227
      std::string name;
1235 1228
      _getColName(cols(id(c)), name);
1236 1229
      return name;
1237 1230
    }
1238 1231

	
1239 1232
    /// Set the name of a column
1240 1233

	
1241 1234
    ///\param c is the coresponding column
1242 1235
    ///\param name The name to be given
1243 1236
    void colName(Col c, const std::string& name) {
1244 1237
      _setColName(cols(id(c)), name);
1245 1238
    }
1246 1239

	
1247 1240
    /// Get the column by its name
1248 1241

	
1249 1242
    ///\param name The name of the column
1250 1243
    ///\return the proper column or \c INVALID
1251 1244
    Col colByName(const std::string& name) const {
1252 1245
      int k = _colByName(name);
1253 1246
      return k != -1 ? Col(cols[k]) : Col(INVALID);
1254 1247
    }
1255 1248

	
1256 1249
    /// Get the name of a row
1257 1250

	
1258 1251
    ///\param r is the coresponding row
1259 1252
    ///\return The name of the row
1260 1253
    std::string rowName(Row r) const {
1261 1254
      std::string name;
1262 1255
      _getRowName(rows(id(r)), name);
1263 1256
      return name;
1264 1257
    }
1265 1258

	
1266 1259
    /// Set the name of a row
1267 1260

	
1268 1261
    ///\param r is the coresponding row
1269 1262
    ///\param name The name to be given
1270 1263
    void rowName(Row r, const std::string& name) {
1271 1264
      _setRowName(rows(id(r)), name);
1272 1265
    }
1273 1266

	
1274 1267
    /// Get the row by its name
1275 1268

	
1276 1269
    ///\param name The name of the row
1277 1270
    ///\return the proper row or \c INVALID
1278 1271
    Row rowByName(const std::string& name) const {
1279 1272
      int k = _rowByName(name);
1280 1273
      return k != -1 ? Row(rows[k]) : Row(INVALID);
1281 1274
    }
1282 1275

	
1283 1276
    /// Set an element of the coefficient matrix of the LP
1284 1277

	
1285 1278
    ///\param r is the row of the element to be modified
1286 1279
    ///\param c is the column of the element to be modified
1287 1280
    ///\param val is the new value of the coefficient
1288 1281
    void coeff(Row r, Col c, Value val) {
1289 1282
      _setCoeff(rows(id(r)),cols(id(c)), val);
1290 1283
    }
1291 1284

	
1292 1285
    /// Get an element of the coefficient matrix of the LP
1293 1286

	
1294 1287
    ///\param r is the row of the element
1295 1288
    ///\param c is the column of the element
1296 1289
    ///\return the corresponding coefficient
1297 1290
    Value coeff(Row r, Col c) const {
1298 1291
      return _getCoeff(rows(id(r)),cols(id(c)));
1299 1292
    }
1300 1293

	
1301 1294
    /// Set the lower bound of a column (i.e a variable)
1302 1295

	
1303 1296
    /// The lower bound of a variable (column) has to be given by an
1304 1297
    /// extended number of type Value, i.e. a finite number of type
1305 1298
    /// Value or -\ref INF.
1306 1299
    void colLowerBound(Col c, Value value) {
1307 1300
      _setColLowerBound(cols(id(c)),value);
1308 1301
    }
1309 1302

	
1310 1303
    /// Get the lower bound of a column (i.e a variable)
1311 1304

	
1312 1305
    /// This function returns the lower bound for column (variable) \c c
1313 1306
    /// (this might be -\ref INF as well).
1314 1307
    ///\return The lower bound for column \c c
1315 1308
    Value colLowerBound(Col c) const {
1316 1309
      return _getColLowerBound(cols(id(c)));
1317 1310
    }
1318 1311

	
1319 1312
    ///\brief Set the lower bound of  several columns
1320 1313
    ///(i.e variables) at once
1321 1314
    ///
1322 1315
    ///This magic function takes a container as its argument
1323 1316
    ///and applies the function on all of its elements.
1324 1317
    ///The lower bound of a variable (column) has to be given by an
1325 1318
    ///extended number of type Value, i.e. a finite number of type
1326 1319
    ///Value or -\ref INF.
1327 1320
#ifdef DOXYGEN
1328 1321
    template<class T>
1329 1322
    void colLowerBound(T &t, Value value) { return 0;}
1330 1323
#else
1331 1324
    template<class T>
1332 1325
    typename enable_if<typename T::value_type::LpCol,void>::type
1333 1326
    colLowerBound(T &t, Value value,dummy<0> = 0) {
1334 1327
      for(typename T::iterator i=t.begin();i!=t.end();++i) {
1335 1328
        colLowerBound(*i, value);
1336 1329
      }
1337 1330
    }
1338 1331
    template<class T>
1339 1332
    typename enable_if<typename T::value_type::second_type::LpCol,
1340 1333
                       void>::type
1341 1334
    colLowerBound(T &t, Value value,dummy<1> = 1) {
1342 1335
      for(typename T::iterator i=t.begin();i!=t.end();++i) {
1343 1336
        colLowerBound(i->second, value);
1344 1337
      }
1345 1338
    }
1346 1339
    template<class T>
1347 1340
    typename enable_if<typename T::MapIt::Value::LpCol,
1348 1341
                       void>::type
1349 1342
    colLowerBound(T &t, Value value,dummy<2> = 2) {
1350 1343
      for(typename T::MapIt i(t); i!=INVALID; ++i){
1351 1344
        colLowerBound(*i, value);
1352 1345
      }
1353 1346
    }
1354 1347
#endif
1355 1348

	
1356 1349
    /// Set the upper bound of a column (i.e a variable)
1357 1350

	
1358 1351
    /// The upper bound of a variable (column) has to be given by an
1359 1352
    /// extended number of type Value, i.e. a finite number of type
1360 1353
    /// Value or \ref INF.
1361 1354
    void colUpperBound(Col c, Value value) {
1362 1355
      _setColUpperBound(cols(id(c)),value);
1363 1356
    };
1364 1357

	
1365 1358
    /// Get the upper bound of a column (i.e a variable)
1366 1359

	
1367 1360
    /// This function returns the upper bound for column (variable) \c c
1368 1361
    /// (this might be \ref INF as well).
1369 1362
    /// \return The upper bound for column \c c
1370 1363
    Value colUpperBound(Col c) const {
1371 1364
      return _getColUpperBound(cols(id(c)));
1372 1365
    }
1373 1366

	
1374 1367
    ///\brief Set the upper bound of  several columns
1375 1368
    ///(i.e variables) at once
1376 1369
    ///
1377 1370
    ///This magic function takes a container as its argument
1378 1371
    ///and applies the function on all of its elements.
1379 1372
    ///The upper bound of a variable (column) has to be given by an
1380 1373
    ///extended number of type Value, i.e. a finite number of type
1381 1374
    ///Value or \ref INF.
1382 1375
#ifdef DOXYGEN
1383 1376
    template<class T>
1384 1377
    void colUpperBound(T &t, Value value) { return 0;}
1385 1378
#else
1386 1379
    template<class T1>
1387 1380
    typename enable_if<typename T1::value_type::LpCol,void>::type
1388 1381
    colUpperBound(T1 &t, Value value,dummy<0> = 0) {
1389 1382
      for(typename T1::iterator i=t.begin();i!=t.end();++i) {
1390 1383
        colUpperBound(*i, value);
1391 1384
      }
1392 1385
    }
1393 1386
    template<class T1>
1394 1387
    typename enable_if<typename T1::value_type::second_type::LpCol,
1395 1388
                       void>::type
1396 1389
    colUpperBound(T1 &t, Value value,dummy<1> = 1) {
1397 1390
      for(typename T1::iterator i=t.begin();i!=t.end();++i) {
1398 1391
        colUpperBound(i->second, value);
1399 1392
      }
1400 1393
    }
1401 1394
    template<class T1>
1402 1395
    typename enable_if<typename T1::MapIt::Value::LpCol,
1403 1396
                       void>::type
1404 1397
    colUpperBound(T1 &t, Value value,dummy<2> = 2) {
1405 1398
      for(typename T1::MapIt i(t); i!=INVALID; ++i){
1406 1399
        colUpperBound(*i, value);
1407 1400
      }
1408 1401
    }
1409 1402
#endif
1410 1403

	
1411 1404
    /// Set the lower and the upper bounds of a column (i.e a variable)
1412 1405

	
1413 1406
    /// The lower and the upper bounds of
1414 1407
    /// a variable (column) have to be given by an
1415 1408
    /// extended number of type Value, i.e. a finite number of type
1416 1409
    /// Value, -\ref INF or \ref INF.
1417 1410
    void colBounds(Col c, Value lower, Value upper) {
1418 1411
      _setColLowerBound(cols(id(c)),lower);
1419 1412
      _setColUpperBound(cols(id(c)),upper);
1420 1413
    }
1421 1414

	
1422 1415
    ///\brief Set the lower and the upper bound of several columns
1423 1416
    ///(i.e variables) at once
1424 1417
    ///
1425 1418
    ///This magic function takes a container as its argument
1426 1419
    ///and applies the function on all of its elements.
1427 1420
    /// The lower and the upper bounds of
1428 1421
    /// a variable (column) have to be given by an
1429 1422
    /// extended number of type Value, i.e. a finite number of type
1430 1423
    /// Value, -\ref INF or \ref INF.
1431 1424
#ifdef DOXYGEN
1432 1425
    template<class T>
1433 1426
    void colBounds(T &t, Value lower, Value upper) { return 0;}
1434 1427
#else
1435 1428
    template<class T2>
1436 1429
    typename enable_if<typename T2::value_type::LpCol,void>::type
1437 1430
    colBounds(T2 &t, Value lower, Value upper,dummy<0> = 0) {
1438 1431
      for(typename T2::iterator i=t.begin();i!=t.end();++i) {
1439 1432
        colBounds(*i, lower, upper);
1440 1433
      }
1441 1434
    }
1442 1435
    template<class T2>
1443 1436
    typename enable_if<typename T2::value_type::second_type::LpCol, void>::type
1444 1437
    colBounds(T2 &t, Value lower, Value upper,dummy<1> = 1) {
1445 1438
      for(typename T2::iterator i=t.begin();i!=t.end();++i) {
1446 1439
        colBounds(i->second, lower, upper);
1447 1440
      }
1448 1441
    }
1449 1442
    template<class T2>
1450 1443
    typename enable_if<typename T2::MapIt::Value::LpCol, void>::type
1451 1444
    colBounds(T2 &t, Value lower, Value upper,dummy<2> = 2) {
1452 1445
      for(typename T2::MapIt i(t); i!=INVALID; ++i){
1453 1446
        colBounds(*i, lower, upper);
1454 1447
      }
1455 1448
    }
1456 1449
#endif
1457 1450

	
1458 1451
    /// Set the lower bound of a row (i.e a constraint)
1459 1452

	
1460 1453
    /// The lower bound of a constraint (row) has to be given by an
1461 1454
    /// extended number of type Value, i.e. a finite number of type
1462 1455
    /// Value or -\ref INF.
1463 1456
    void rowLowerBound(Row r, Value value) {
1464 1457
      _setRowLowerBound(rows(id(r)),value);
1465 1458
    }
1466 1459

	
1467 1460
    /// Get the lower bound of a row (i.e a constraint)
1468 1461

	
1469 1462
    /// This function returns the lower bound for row (constraint) \c c
1470 1463
    /// (this might be -\ref INF as well).
1471 1464
    ///\return The lower bound for row \c r
1472 1465
    Value rowLowerBound(Row r) const {
1473 1466
      return _getRowLowerBound(rows(id(r)));
1474 1467
    }
1475 1468

	
1476 1469
    /// Set the upper bound of a row (i.e a constraint)
1477 1470

	
1478 1471
    /// The upper bound of a constraint (row) has to be given by an
1479 1472
    /// extended number of type Value, i.e. a finite number of type
1480 1473
    /// Value or -\ref INF.
1481 1474
    void rowUpperBound(Row r, Value value) {
1482 1475
      _setRowUpperBound(rows(id(r)),value);
1483 1476
    }
1484 1477

	
1485 1478
    /// Get the upper bound of a row (i.e a constraint)
1486 1479

	
1487 1480
    /// This function returns the upper bound for row (constraint) \c c
1488 1481
    /// (this might be -\ref INF as well).
1489 1482
    ///\return The upper bound for row \c r
1490 1483
    Value rowUpperBound(Row r) const {
1491 1484
      return _getRowUpperBound(rows(id(r)));
1492 1485
    }
1493 1486

	
1494 1487
    ///Set an element of the objective function
1495 1488
    void objCoeff(Col c, Value v) {_setObjCoeff(cols(id(c)),v); };
1496 1489

	
1497 1490
    ///Get an element of the objective function
1498 1491
    Value objCoeff(Col c) const { return _getObjCoeff(cols(id(c))); };
1499 1492

	
1500 1493
    ///Set the objective function
1501 1494

	
1502 1495
    ///\param e is a linear expression of type \ref Expr.
1503 1496
    ///
1504 1497
    void obj(const Expr& e) {
1505 1498
      _setObjCoeffs(ExprIterator(e.comps.begin(), cols),
1506 1499
                    ExprIterator(e.comps.end(), cols));
1507 1500
      obj_const_comp = *e;
1508 1501
    }
1509 1502

	
1510 1503
    ///Get the objective function
1511 1504

	
1512 1505
    ///\return the objective function as a linear expression of type
1513 1506
    ///Expr.
1514 1507
    Expr obj() const {
1515 1508
      Expr e;
1516 1509
      _getObjCoeffs(InsertIterator(e.comps, cols));
1517 1510
      *e = obj_const_comp;
1518 1511
      return e;
1519 1512
    }
1520 1513

	
1521 1514

	
1522 1515
    ///Set the direction of optimization
1523 1516
    void sense(Sense sense) { _setSense(sense); }
1524 1517

	
1525 1518
    ///Query the direction of the optimization
1526 1519
    Sense sense() const {return _getSense(); }
1527 1520

	
1528 1521
    ///Set the sense to maximization
1529 1522
    void max() { _setSense(MAX); }
1530 1523

	
1531 1524
    ///Set the sense to maximization
1532 1525
    void min() { _setSense(MIN); }
1533 1526

	
1534 1527
    ///Clears the problem
1535 1528
    void clear() { _clear(); }
1536 1529

	
1537 1530
    ///@}
1538 1531

	
1539 1532
  };
1540 1533

	
1541 1534
  /// Addition
1542 1535

	
1543 1536
  ///\relates LpBase::Expr
1544 1537
  ///
1545 1538
  inline LpBase::Expr operator+(const LpBase::Expr &a, const LpBase::Expr &b) {
1546 1539
    LpBase::Expr tmp(a);
1547 1540
    tmp+=b;
1548 1541
    return tmp;
1549 1542
  }
1550 1543
  ///Substraction
1551 1544

	
1552 1545
  ///\relates LpBase::Expr
1553 1546
  ///
1554 1547
  inline LpBase::Expr operator-(const LpBase::Expr &a, const LpBase::Expr &b) {
1555 1548
    LpBase::Expr tmp(a);
1556 1549
    tmp-=b;
1557 1550
    return tmp;
1558 1551
  }
1559 1552
  ///Multiply with constant
1560 1553

	
1561 1554
  ///\relates LpBase::Expr
1562 1555
  ///
1563 1556
  inline LpBase::Expr operator*(const LpBase::Expr &a, const LpBase::Value &b) {
1564 1557
    LpBase::Expr tmp(a);
1565 1558
    tmp*=b;
1566 1559
    return tmp;
1567 1560
  }
1568 1561

	
1569 1562
  ///Multiply with constant
1570 1563

	
1571 1564
  ///\relates LpBase::Expr
1572 1565
  ///
1573 1566
  inline LpBase::Expr operator*(const LpBase::Value &a, const LpBase::Expr &b) {
1574 1567
    LpBase::Expr tmp(b);
1575 1568
    tmp*=a;
1576 1569
    return tmp;
1577 1570
  }
1578 1571
  ///Divide with constant
1579 1572

	
1580 1573
  ///\relates LpBase::Expr
1581 1574
  ///
1582 1575
  inline LpBase::Expr operator/(const LpBase::Expr &a, const LpBase::Value &b) {
1583 1576
    LpBase::Expr tmp(a);
1584 1577
    tmp/=b;
1585 1578
    return tmp;
1586 1579
  }
1587 1580

	
1588 1581
  ///Create constraint
1589 1582

	
1590 1583
  ///\relates LpBase::Constr
1591 1584
  ///
1592 1585
  inline LpBase::Constr operator<=(const LpBase::Expr &e,
1593 1586
                                   const LpBase::Expr &f) {
1594 1587
    return LpBase::Constr(0, f - e, LpBase::INF);
1595 1588
  }
1596 1589

	
1597 1590
  ///Create constraint
1598 1591

	
1599 1592
  ///\relates LpBase::Constr
1600 1593
  ///
1601 1594
  inline LpBase::Constr operator<=(const LpBase::Value &e,
1602 1595
                                   const LpBase::Expr &f) {
1603 1596
    return LpBase::Constr(e, f, LpBase::NaN);
1604 1597
  }
1605 1598

	
1606 1599
  ///Create constraint
1607 1600

	
1608 1601
  ///\relates LpBase::Constr
1609 1602
  ///
1610 1603
  inline LpBase::Constr operator<=(const LpBase::Expr &e,
1611 1604
                                   const LpBase::Value &f) {
1612 1605
    return LpBase::Constr(- LpBase::INF, e, f);
1613 1606
  }
1614 1607

	
1615 1608
  ///Create constraint
1616 1609

	
1617 1610
  ///\relates LpBase::Constr
1618 1611
  ///
1619 1612
  inline LpBase::Constr operator>=(const LpBase::Expr &e,
1620 1613
                                   const LpBase::Expr &f) {
1621 1614
    return LpBase::Constr(0, e - f, LpBase::INF);
1622 1615
  }
1623 1616

	
1624 1617

	
1625 1618
  ///Create constraint
1626 1619

	
1627 1620
  ///\relates LpBase::Constr
1628 1621
  ///
1629 1622
  inline LpBase::Constr operator>=(const LpBase::Value &e,
1630 1623
                                   const LpBase::Expr &f) {
1631 1624
    return LpBase::Constr(LpBase::NaN, f, e);
1632 1625
  }
1633 1626

	
1634 1627

	
1635 1628
  ///Create constraint
1636 1629

	
1637 1630
  ///\relates LpBase::Constr
1638 1631
  ///
1639 1632
  inline LpBase::Constr operator>=(const LpBase::Expr &e,
1640 1633
                                   const LpBase::Value &f) {
1641 1634
    return LpBase::Constr(f, e, LpBase::INF);
1642 1635
  }
1643 1636

	
1644 1637
  ///Create constraint
1645 1638

	
1646 1639
  ///\relates LpBase::Constr
1647 1640
  ///
1648 1641
  inline LpBase::Constr operator==(const LpBase::Expr &e,
1649 1642
                                   const LpBase::Value &f) {
1650 1643
    return LpBase::Constr(f, e, f);
1651 1644
  }
1652 1645

	
1653 1646
  ///Create constraint
1654 1647

	
1655 1648
  ///\relates LpBase::Constr
1656 1649
  ///
1657 1650
  inline LpBase::Constr operator==(const LpBase::Expr &e,
1658 1651
                                   const LpBase::Expr &f) {
1659 1652
    return LpBase::Constr(0, f - e, 0);
1660 1653
  }
1661 1654

	
1662 1655
  ///Create constraint
1663 1656

	
1664 1657
  ///\relates LpBase::Constr
1665 1658
  ///
1666 1659
  inline LpBase::Constr operator<=(const LpBase::Value &n,
1667 1660
                                   const LpBase::Constr &c) {
1668 1661
    LpBase::Constr tmp(c);
1669 1662
    LEMON_ASSERT(isNaN(tmp.lowerBound()), "Wrong LP constraint");
1670 1663
    tmp.lowerBound()=n;
1671 1664
    return tmp;
1672 1665
  }
1673 1666
  ///Create constraint
1674 1667

	
1675 1668
  ///\relates LpBase::Constr
1676 1669
  ///
1677 1670
  inline LpBase::Constr operator<=(const LpBase::Constr &c,
1678 1671
                                   const LpBase::Value &n)
1679 1672
  {
1680 1673
    LpBase::Constr tmp(c);
1681 1674
    LEMON_ASSERT(isNaN(tmp.upperBound()), "Wrong LP constraint");
1682 1675
    tmp.upperBound()=n;
1683 1676
    return tmp;
1684 1677
  }
1685 1678

	
1686 1679
  ///Create constraint
1687 1680

	
1688 1681
  ///\relates LpBase::Constr
1689 1682
  ///
1690 1683
  inline LpBase::Constr operator>=(const LpBase::Value &n,
1691 1684
                                   const LpBase::Constr &c) {
1692 1685
    LpBase::Constr tmp(c);
1693 1686
    LEMON_ASSERT(isNaN(tmp.upperBound()), "Wrong LP constraint");
1694 1687
    tmp.upperBound()=n;
1695 1688
    return tmp;
1696 1689
  }
1697 1690
  ///Create constraint
1698 1691

	
1699 1692
  ///\relates LpBase::Constr
1700 1693
  ///
1701 1694
  inline LpBase::Constr operator>=(const LpBase::Constr &c,
1702 1695
                                   const LpBase::Value &n)
1703 1696
  {
1704 1697
    LpBase::Constr tmp(c);
1705 1698
    LEMON_ASSERT(isNaN(tmp.lowerBound()), "Wrong LP constraint");
1706 1699
    tmp.lowerBound()=n;
1707 1700
    return tmp;
1708 1701
  }
1709 1702

	
1710 1703
  ///Addition
1711 1704

	
1712 1705
  ///\relates LpBase::DualExpr
1713 1706
  ///
1714 1707
  inline LpBase::DualExpr operator+(const LpBase::DualExpr &a,
1715 1708
                                    const LpBase::DualExpr &b) {
1716 1709
    LpBase::DualExpr tmp(a);
1717 1710
    tmp+=b;
1718 1711
    return tmp;
1719 1712
  }
1720 1713
  ///Substraction
1721 1714

	
1722 1715
  ///\relates LpBase::DualExpr
1723 1716
  ///
1724 1717
  inline LpBase::DualExpr operator-(const LpBase::DualExpr &a,
1725 1718
                                    const LpBase::DualExpr &b) {
1726 1719
    LpBase::DualExpr tmp(a);
1727 1720
    tmp-=b;
1728 1721
    return tmp;
1729 1722
  }
1730 1723
  ///Multiply with constant
1731 1724

	
1732 1725
  ///\relates LpBase::DualExpr
1733 1726
  ///
1734 1727
  inline LpBase::DualExpr operator*(const LpBase::DualExpr &a,
1735 1728
                                    const LpBase::Value &b) {
1736 1729
    LpBase::DualExpr tmp(a);
1737 1730
    tmp*=b;
1738 1731
    return tmp;
1739 1732
  }
1740 1733

	
1741 1734
  ///Multiply with constant
1742 1735

	
1743 1736
  ///\relates LpBase::DualExpr
1744 1737
  ///
1745 1738
  inline LpBase::DualExpr operator*(const LpBase::Value &a,
1746 1739
                                    const LpBase::DualExpr &b) {
1747 1740
    LpBase::DualExpr tmp(b);
1748 1741
    tmp*=a;
1749 1742
    return tmp;
1750 1743
  }
1751 1744
  ///Divide with constant
1752 1745

	
1753 1746
  ///\relates LpBase::DualExpr
1754 1747
  ///
1755 1748
  inline LpBase::DualExpr operator/(const LpBase::DualExpr &a,
1756 1749
                                    const LpBase::Value &b) {
1757 1750
    LpBase::DualExpr tmp(a);
1758 1751
    tmp/=b;
1759 1752
    return tmp;
1760 1753
  }
1761 1754

	
1762 1755
  /// \ingroup lp_group
1763 1756
  ///
1764 1757
  /// \brief Common base class for LP solvers
1765 1758
  ///
1766 1759
  /// This class is an abstract base class for LP solvers. This class
1767 1760
  /// provides a full interface for set and modify an LP problem,
1768 1761
  /// solve it and retrieve the solution. You can use one of the
1769 1762
  /// descendants as a concrete implementation, or the \c Lp
1770 1763
  /// default LP solver. However, if you would like to handle LP
1771 1764
  /// solvers as reference or pointer in a generic way, you can use
1772 1765
  /// this class directly.
1773 1766
  class LpSolver : virtual public LpBase {
1774 1767
  public:
1775 1768

	
1776 1769
    /// The problem types for primal and dual problems
1777 1770
    enum ProblemType {
1778 1771
      ///Feasible solution hasn't been found (but may exist).
1779 1772
      UNDEFINED = 0,
1780 1773
      ///The problem has no feasible solution
1781 1774
      INFEASIBLE = 1,
1782 1775
      ///Feasible solution found
1783 1776
      FEASIBLE = 2,
1784 1777
      ///Optimal solution exists and found
1785 1778
      OPTIMAL = 3,
1786 1779
      ///The cost function is unbounded
1787 1780
      UNBOUNDED = 4
1788 1781
    };
1789 1782

	
1790 1783
    ///The basis status of variables
1791 1784
    enum VarStatus {
1792 1785
      /// The variable is in the basis
1793 1786
      BASIC, 
1794 1787
      /// The variable is free, but not basic
1795 1788
      FREE,
1796 1789
      /// The variable has active lower bound 
1797 1790
      LOWER,
1798 1791
      /// The variable has active upper bound
1799 1792
      UPPER,
1800 1793
      /// The variable is non-basic and fixed
1801 1794
      FIXED
1802 1795
    };
1803 1796

	
1804 1797
  protected:
1805 1798

	
1806 1799
    virtual SolveExitStatus _solve() = 0;
1807 1800

	
1808 1801
    virtual Value _getPrimal(int i) const = 0;
1809 1802
    virtual Value _getDual(int i) const = 0;
1810 1803

	
1811 1804
    virtual Value _getPrimalRay(int i) const = 0;
1812 1805
    virtual Value _getDualRay(int i) const = 0;
1813 1806

	
1814 1807
    virtual Value _getPrimalValue() const = 0;
1815 1808

	
1816 1809
    virtual VarStatus _getColStatus(int i) const = 0;
1817 1810
    virtual VarStatus _getRowStatus(int i) const = 0;
1818 1811

	
1819 1812
    virtual ProblemType _getPrimalType() const = 0;
1820 1813
    virtual ProblemType _getDualType() const = 0;
1821 1814

	
1822 1815
  public:
1823 1816

	
1817
    ///Allocate a new LP problem instance
1818
    virtual LpSolver* newSolver() const = 0;
1819
    ///Make a copy of the LP problem
1820
    virtual LpSolver* cloneSolver() const = 0;
1821

	
1824 1822
    ///\name Solve the LP
1825 1823

	
1826 1824
    ///@{
1827 1825

	
1828 1826
    ///\e Solve the LP problem at hand
1829 1827
    ///
1830 1828
    ///\return The result of the optimization procedure. Possible
1831 1829
    ///values and their meanings can be found in the documentation of
1832 1830
    ///\ref SolveExitStatus.
1833 1831
    SolveExitStatus solve() { return _solve(); }
1834 1832

	
1835 1833
    ///@}
1836 1834

	
1837 1835
    ///\name Obtain the solution
1838 1836

	
1839 1837
    ///@{
1840 1838

	
1841 1839
    /// The type of the primal problem
1842 1840
    ProblemType primalType() const {
1843 1841
      return _getPrimalType();
1844 1842
    }
1845 1843

	
1846 1844
    /// The type of the dual problem
1847 1845
    ProblemType dualType() const {
1848 1846
      return _getDualType();
1849 1847
    }
1850 1848

	
1851 1849
    /// Return the primal value of the column
1852 1850

	
1853 1851
    /// Return the primal value of the column.
1854 1852
    /// \pre The problem is solved.
1855 1853
    Value primal(Col c) const { return _getPrimal(cols(id(c))); }
1856 1854

	
1857 1855
    /// Return the primal value of the expression
1858 1856

	
1859 1857
    /// Return the primal value of the expression, i.e. the dot
1860 1858
    /// product of the primal solution and the expression.
1861 1859
    /// \pre The problem is solved.
1862 1860
    Value primal(const Expr& e) const {
1863 1861
      double res = *e;
1864 1862
      for (Expr::ConstCoeffIt c(e); c != INVALID; ++c) {
1865 1863
        res += *c * primal(c);
1866 1864
      }
1867 1865
      return res;
1868 1866
    }
1869 1867
    /// Returns a component of the primal ray
1870 1868
    
1871 1869
    /// The primal ray is solution of the modified primal problem,
1872 1870
    /// where we change each finite bound to 0, and we looking for a
1873 1871
    /// negative objective value in case of minimization, and positive
1874 1872
    /// objective value for maximization. If there is such solution,
1875 1873
    /// that proofs the unsolvability of the dual problem, and if a
1876 1874
    /// feasible primal solution exists, then the unboundness of
1877 1875
    /// primal problem.
1878 1876
    ///
1879 1877
    /// \pre The problem is solved and the dual problem is infeasible.
1880 1878
    /// \note Some solvers does not provide primal ray calculation
1881 1879
    /// functions.
1882 1880
    Value primalRay(Col c) const { return _getPrimalRay(cols(id(c))); }
1883 1881

	
1884 1882
    /// Return the dual value of the row
1885 1883

	
1886 1884
    /// Return the dual value of the row.
1887 1885
    /// \pre The problem is solved.
1888 1886
    Value dual(Row r) const { return _getDual(rows(id(r))); }
1889 1887

	
1890 1888
    /// Return the dual value of the dual expression
1891 1889

	
1892 1890
    /// Return the dual value of the dual expression, i.e. the dot
1893 1891
    /// product of the dual solution and the dual expression.
1894 1892
    /// \pre The problem is solved.
1895 1893
    Value dual(const DualExpr& e) const {
1896 1894
      double res = 0.0;
1897 1895
      for (DualExpr::ConstCoeffIt r(e); r != INVALID; ++r) {
1898 1896
        res += *r * dual(r);
1899 1897
      }
1900 1898
      return res;
1901 1899
    }
1902 1900

	
1903 1901
    /// Returns a component of the dual ray
1904 1902
    
1905 1903
    /// The dual ray is solution of the modified primal problem, where
1906 1904
    /// we change each finite bound to 0 (i.e. the objective function
1907 1905
    /// coefficients in the primal problem), and we looking for a
1908 1906
    /// ositive objective value. If there is such solution, that
1909 1907
    /// proofs the unsolvability of the primal problem, and if a
1910 1908
    /// feasible dual solution exists, then the unboundness of
1911 1909
    /// dual problem.
1912 1910
    ///
1913 1911
    /// \pre The problem is solved and the primal problem is infeasible.
1914 1912
    /// \note Some solvers does not provide dual ray calculation
1915 1913
    /// functions.
1916 1914
    Value dualRay(Row r) const { return _getDualRay(rows(id(r))); }
1917 1915

	
1918 1916
    /// Return the basis status of the column
1919 1917

	
1920 1918
    /// \see VarStatus
1921 1919
    VarStatus colStatus(Col c) const { return _getColStatus(cols(id(c))); }
1922 1920

	
1923 1921
    /// Return the basis status of the row
1924 1922

	
1925 1923
    /// \see VarStatus
1926 1924
    VarStatus rowStatus(Row r) const { return _getRowStatus(rows(id(r))); }
1927 1925

	
1928 1926
    ///The value of the objective function
1929 1927

	
1930 1928
    ///\return
1931 1929
    ///- \ref INF or -\ref INF means either infeasibility or unboundedness
1932 1930
    /// of the primal problem, depending on whether we minimize or maximize.
1933 1931
    ///- \ref NaN if no primal solution is found.
1934 1932
    ///- The (finite) objective value if an optimal solution is found.
1935 1933
    Value primal() const { return _getPrimalValue()+obj_const_comp;}
1936 1934
    ///@}
1937 1935

	
1938
    LpSolver* newSolver() {return _newSolver();}
1939
    LpSolver* cloneSolver() {return _cloneSolver();}
1940

	
1941 1936
  protected:
1942 1937

	
1943
    virtual LpSolver* _newSolver() const = 0;
1944
    virtual LpSolver* _cloneSolver() const = 0;
1945 1938
  };
1946 1939

	
1947 1940

	
1948 1941
  /// \ingroup lp_group
1949 1942
  ///
1950 1943
  /// \brief Common base class for MIP solvers
1951 1944
  ///
1952 1945
  /// This class is an abstract base class for MIP solvers. This class
1953 1946
  /// provides a full interface for set and modify an MIP problem,
1954 1947
  /// solve it and retrieve the solution. You can use one of the
1955 1948
  /// descendants as a concrete implementation, or the \c Lp
1956 1949
  /// default MIP solver. However, if you would like to handle MIP
1957 1950
  /// solvers as reference or pointer in a generic way, you can use
1958 1951
  /// this class directly.
1959 1952
  class MipSolver : virtual public LpBase {
1960 1953
  public:
1961 1954

	
1962 1955
    /// The problem types for MIP problems
1963 1956
    enum ProblemType {
1964 1957
      ///Feasible solution hasn't been found (but may exist).
1965 1958
      UNDEFINED = 0,
1966 1959
      ///The problem has no feasible solution
1967 1960
      INFEASIBLE = 1,
1968 1961
      ///Feasible solution found
1969 1962
      FEASIBLE = 2,
1970 1963
      ///Optimal solution exists and found
1971 1964
      OPTIMAL = 3,
1972 1965
      ///The cost function is unbounded
1973 1966
      ///
1974 1967
      ///The Mip or at least the relaxed problem is unbounded
1975 1968
      UNBOUNDED = 4
1976 1969
    };
1977 1970

	
1971
    ///Allocate a new MIP problem instance
1972
    virtual MipSolver* newSolver() const = 0;
1973
    ///Make a copy of the MIP problem
1974
    virtual MipSolver* cloneSolver() const = 0;
1975

	
1978 1976
    ///\name Solve the MIP
1979 1977

	
1980 1978
    ///@{
1981 1979

	
1982 1980
    /// Solve the MIP problem at hand
1983 1981
    ///
1984 1982
    ///\return The result of the optimization procedure. Possible
1985 1983
    ///values and their meanings can be found in the documentation of
1986 1984
    ///\ref SolveExitStatus.
1987 1985
    SolveExitStatus solve() { return _solve(); }
1988 1986

	
1989 1987
    ///@}
1990 1988

	
1991 1989
    ///\name Setting column type
1992 1990
    ///@{
1993 1991

	
1994 1992
    ///Possible variable (column) types (e.g. real, integer, binary etc.)
1995 1993
    enum ColTypes {
1996 1994
      ///Continuous variable (default)
1997 1995
      REAL = 0,
1998 1996
      ///Integer variable
1999 1997
      INTEGER = 1
2000 1998
    };
2001 1999

	
2002 2000
    ///Sets the type of the given column to the given type
2003 2001

	
2004 2002
    ///Sets the type of the given column to the given type.
2005 2003
    ///
2006 2004
    void colType(Col c, ColTypes col_type) {
2007 2005
      _setColType(cols(id(c)),col_type);
2008 2006
    }
2009 2007

	
2010 2008
    ///Gives back the type of the column.
2011 2009

	
2012 2010
    ///Gives back the type of the column.
2013 2011
    ///
2014 2012
    ColTypes colType(Col c) const {
2015 2013
      return _getColType(cols(id(c)));
2016 2014
    }
2017 2015
    ///@}
2018 2016

	
2019 2017
    ///\name Obtain the solution
2020 2018

	
2021 2019
    ///@{
2022 2020

	
2023 2021
    /// The type of the MIP problem
2024 2022
    ProblemType type() const {
2025 2023
      return _getType();
2026 2024
    }
2027 2025

	
2028 2026
    /// Return the value of the row in the solution
2029 2027

	
2030 2028
    ///  Return the value of the row in the solution.
2031 2029
    /// \pre The problem is solved.
2032 2030
    Value sol(Col c) const { return _getSol(cols(id(c))); }
2033 2031

	
2034 2032
    /// Return the value of the expression in the solution
2035 2033

	
2036 2034
    /// Return the value of the expression in the solution, i.e. the
2037 2035
    /// dot product of the solution and the expression.
2038 2036
    /// \pre The problem is solved.
2039 2037
    Value sol(const Expr& e) const {
2040 2038
      double res = *e;
2041 2039
      for (Expr::ConstCoeffIt c(e); c != INVALID; ++c) {
2042 2040
        res += *c * sol(c);
2043 2041
      }
2044 2042
      return res;
2045 2043
    }
2046 2044
    ///The value of the objective function
2047 2045
    
2048 2046
    ///\return
2049 2047
    ///- \ref INF or -\ref INF means either infeasibility or unboundedness
2050 2048
    /// of the problem, depending on whether we minimize or maximize.
2051 2049
    ///- \ref NaN if no primal solution is found.
2052 2050
    ///- The (finite) objective value if an optimal solution is found.
2053 2051
    Value solValue() const { return _getSolValue()+obj_const_comp;}
2054 2052
    ///@}
2055 2053

	
2056 2054
  protected:
2057 2055

	
2058 2056
    virtual SolveExitStatus _solve() = 0;
2059 2057
    virtual ColTypes _getColType(int col) const = 0;
2060 2058
    virtual void _setColType(int col, ColTypes col_type) = 0;
2061 2059
    virtual ProblemType _getType() const = 0;
2062 2060
    virtual Value _getSol(int i) const = 0;
2063 2061
    virtual Value _getSolValue() const = 0;
2064 2062

	
2065
  public:
2066

	
2067
    MipSolver* newSolver() {return _newSolver();}
2068
    MipSolver* cloneSolver() {return _cloneSolver();}
2069

	
2070
  protected:
2071

	
2072
    virtual MipSolver* _newSolver() const = 0;
2073
    virtual MipSolver* _cloneSolver() const = 0;
2074 2063
  };
2075 2064

	
2076 2065

	
2077 2066

	
2078 2067
} //namespace lemon
2079 2068

	
2080 2069
#endif //LEMON_LP_BASE_H
Ignore white space 1536 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 <lemon/lp_skeleton.h>
20 20

	
21 21
///\file
22 22
///\brief A skeleton file to implement LP solver interfaces
23 23
namespace lemon {
24 24

	
25 25
  int SkeletonSolverBase::_addCol()
26 26
  {
27 27
    return ++col_num;
28 28
  }
29 29

	
30 30
  int SkeletonSolverBase::_addRow()
31 31
  {
32 32
    return ++row_num;
33 33
  }
34 34

	
35 35
  void SkeletonSolverBase::_eraseCol(int) {}
36 36
  void SkeletonSolverBase::_eraseRow(int) {}
37 37

	
38 38
  void SkeletonSolverBase::_getColName(int, std::string &) const {}
39 39
  void SkeletonSolverBase::_setColName(int, const std::string &) {}
40 40
  int SkeletonSolverBase::_colByName(const std::string&) const { return -1; }
41 41

	
42 42
  void SkeletonSolverBase::_getRowName(int, std::string &) const {}
43 43
  void SkeletonSolverBase::_setRowName(int, const std::string &) {}
44 44
  int SkeletonSolverBase::_rowByName(const std::string&) const { return -1; }
45 45

	
46 46
  void SkeletonSolverBase::_setRowCoeffs(int, ExprIterator, ExprIterator) {}
47 47
  void SkeletonSolverBase::_getRowCoeffs(int, InsertIterator) const {}
48 48

	
49 49
  void SkeletonSolverBase::_setColCoeffs(int, ExprIterator, ExprIterator) {}
50 50
  void SkeletonSolverBase::_getColCoeffs(int, InsertIterator) const {}
51 51

	
52 52
  void SkeletonSolverBase::_setCoeff(int, int, Value) {}
53 53
  SkeletonSolverBase::Value SkeletonSolverBase::_getCoeff(int, int) const
54 54
  { return 0; }
55 55

	
56 56
  void SkeletonSolverBase::_setColLowerBound(int, Value) {}
57 57
  SkeletonSolverBase::Value SkeletonSolverBase::_getColLowerBound(int) const
58 58
  {  return 0; }
59 59

	
60 60
  void SkeletonSolverBase::_setColUpperBound(int, Value) {}
61 61
  SkeletonSolverBase::Value SkeletonSolverBase::_getColUpperBound(int) const
62 62
  {  return 0; }
63 63

	
64 64
  void SkeletonSolverBase::_setRowLowerBound(int, Value) {}
65 65
  SkeletonSolverBase::Value SkeletonSolverBase::_getRowLowerBound(int) const
66 66
  {  return 0; }
67 67

	
68 68
  void SkeletonSolverBase::_setRowUpperBound(int, Value) {}
69 69
  SkeletonSolverBase::Value SkeletonSolverBase::_getRowUpperBound(int) const
70 70
  {  return 0; }
71 71

	
72 72
  void SkeletonSolverBase::_setObjCoeffs(ExprIterator, ExprIterator) {}
73 73
  void SkeletonSolverBase::_getObjCoeffs(InsertIterator) const {};
74 74

	
75 75
  void SkeletonSolverBase::_setObjCoeff(int, Value) {}
76 76
  SkeletonSolverBase::Value SkeletonSolverBase::_getObjCoeff(int) const
77 77
  {  return 0; }
78 78

	
79 79
  void SkeletonSolverBase::_setSense(Sense) {}
80 80
  SkeletonSolverBase::Sense SkeletonSolverBase::_getSense() const
81 81
  { return MIN; }
82 82

	
83 83
  void SkeletonSolverBase::_clear() {
84 84
    row_num = col_num = 0;
85 85
  }
86 86

	
87 87
  LpSkeleton::SolveExitStatus LpSkeleton::_solve() { return SOLVED; }
88 88

	
89 89
  LpSkeleton::Value LpSkeleton::_getPrimal(int) const { return 0; }
90 90
  LpSkeleton::Value LpSkeleton::_getDual(int) const { return 0; }
91 91
  LpSkeleton::Value LpSkeleton::_getPrimalValue() const { return 0; }
92 92

	
93 93
  LpSkeleton::Value LpSkeleton::_getPrimalRay(int) const { return 0; }
94 94
  LpSkeleton::Value LpSkeleton::_getDualRay(int) const { return 0; }
95 95

	
96 96
  LpSkeleton::ProblemType LpSkeleton::_getPrimalType() const
97 97
  { return UNDEFINED; }
98 98

	
99 99
  LpSkeleton::ProblemType LpSkeleton::_getDualType() const
100 100
  { return UNDEFINED; }
101 101

	
102 102
  LpSkeleton::VarStatus LpSkeleton::_getColStatus(int) const
103 103
  { return BASIC; }
104 104

	
105 105
  LpSkeleton::VarStatus LpSkeleton::_getRowStatus(int) const
106 106
  { return BASIC; }
107 107

	
108
  LpSkeleton* LpSkeleton::_newSolver() const
108
  LpSkeleton* LpSkeleton::newSolver() const
109 109
  { return static_cast<LpSkeleton*>(0); }
110 110

	
111
  LpSkeleton* LpSkeleton::_cloneSolver() const
111
  LpSkeleton* LpSkeleton::cloneSolver() const
112 112
  { return static_cast<LpSkeleton*>(0); }
113 113

	
114 114
  const char* LpSkeleton::_solverName() const { return "LpSkeleton"; }
115 115

	
116 116
  MipSkeleton::SolveExitStatus MipSkeleton::_solve()
117 117
  { return SOLVED; }
118 118

	
119 119
  MipSkeleton::Value MipSkeleton::_getSol(int) const { return 0; }
120 120
  MipSkeleton::Value MipSkeleton::_getSolValue() const { return 0; }
121 121

	
122 122
  MipSkeleton::ProblemType MipSkeleton::_getType() const
123 123
  { return UNDEFINED; }
124 124

	
125
  MipSkeleton* MipSkeleton::_newSolver() const
125
  MipSkeleton* MipSkeleton::newSolver() const
126 126
  { return static_cast<MipSkeleton*>(0); }
127 127

	
128
  MipSkeleton* MipSkeleton::_cloneSolver() const
128
  MipSkeleton* MipSkeleton::cloneSolver() const
129 129
  { return static_cast<MipSkeleton*>(0); }
130 130

	
131 131
  const char* MipSkeleton::_solverName() const { return "MipSkeleton"; }
132 132

	
133 133
} //namespace lemon
134 134

	
Ignore white space 1536 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_LP_SKELETON_H
20 20
#define LEMON_LP_SKELETON_H
21 21

	
22 22
#include <lemon/lp_base.h>
23 23

	
24 24
///\file
25
///\brief A skeleton file to implement LP solver interfaces
25
///\brief Skeleton file to implement LP/MIP solver interfaces
26
///  
27
///The classes in this file do nothing, but they can serve as skeletons when
28
///implementing an interface to new solvers.
26 29
namespace lemon {
27 30

	
28
  ///A skeleton class to implement LP solver interfaces
31
  ///A skeleton class to implement LP/MIP solver base interface
32
  
33
  ///This class does nothing, but it can serve as a skeleton when
34
  ///implementing an interface to new solvers.
29 35
  class SkeletonSolverBase : public virtual LpBase {
30 36
    int col_num,row_num;
31 37

	
32 38
  protected:
33 39

	
34 40
    SkeletonSolverBase()
35 41
      : col_num(-1), row_num(-1) {}
36 42

	
37 43
    /// \e
38 44
    virtual int _addCol();
39 45
    /// \e
40 46
    virtual int _addRow();
41 47
    /// \e
42 48
    virtual void _eraseCol(int i);
43 49
    /// \e
44 50
    virtual void _eraseRow(int i);
45 51

	
46 52
    /// \e
47 53
    virtual void _getColName(int col, std::string& name) const;
48 54
    /// \e
49 55
    virtual void _setColName(int col, const std::string& name);
50 56
    /// \e
51 57
    virtual int _colByName(const std::string& name) const;
52 58

	
53 59
    /// \e
54 60
    virtual void _getRowName(int row, std::string& name) const;
55 61
    /// \e
56 62
    virtual void _setRowName(int row, const std::string& name);
57 63
    /// \e
58 64
    virtual int _rowByName(const std::string& name) const;
59 65

	
60 66
    /// \e
61 67
    virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e);
62 68
    /// \e
63 69
    virtual void _getRowCoeffs(int i, InsertIterator b) const;
64 70
    /// \e
65 71
    virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e);
66 72
    /// \e
67 73
    virtual void _getColCoeffs(int i, InsertIterator b) const;
68 74

	
69 75
    /// Set one element of the coefficient matrix
70 76
    virtual void _setCoeff(int row, int col, Value value);
71 77

	
72 78
    /// Get one element of the coefficient matrix
73 79
    virtual Value _getCoeff(int row, int col) const;
74 80

	
75 81
    /// The lower bound of a variable (column) have to be given by an
76 82
    /// extended number of type Value, i.e. a finite number of type
77 83
    /// Value or -\ref INF.
78 84
    virtual void _setColLowerBound(int i, Value value);
79 85
    /// \e
80 86

	
81 87
    /// The lower bound of a variable (column) is an
82 88
    /// extended number of type Value, i.e. a finite number of type
83 89
    /// Value or -\ref INF.
84 90
    virtual Value _getColLowerBound(int i) const;
85 91

	
86 92
    /// The upper bound of a variable (column) have to be given by an
87 93
    /// extended number of type Value, i.e. a finite number of type
88 94
    /// Value or \ref INF.
89 95
    virtual void _setColUpperBound(int i, Value value);
90 96
    /// \e
91 97

	
92 98
    /// The upper bound of a variable (column) is an
93 99
    /// extended number of type Value, i.e. a finite number of type
94 100
    /// Value or \ref INF.
95 101
    virtual Value _getColUpperBound(int i) const;
96 102

	
97 103
    /// The lower bound of a constraint (row) have to be given by an
98 104
    /// extended number of type Value, i.e. a finite number of type
99 105
    /// Value or -\ref INF.
100 106
    virtual void _setRowLowerBound(int i, Value value);
101 107
    /// \e
102 108

	
103 109
    /// The lower bound of a constraint (row) is an
104 110
    /// extended number of type Value, i.e. a finite number of type
105 111
    /// Value or -\ref INF.
106 112
    virtual Value _getRowLowerBound(int i) const;
107 113

	
108 114
    /// The upper bound of a constraint (row) have to be given by an
109 115
    /// extended number of type Value, i.e. a finite number of type
110 116
    /// Value or \ref INF.
111 117
    virtual void _setRowUpperBound(int i, Value value);
112 118
    /// \e
113 119

	
114 120
    /// The upper bound of a constraint (row) is an
115 121
    /// extended number of type Value, i.e. a finite number of type
116 122
    /// Value or \ref INF.
117 123
    virtual Value _getRowUpperBound(int i) const;
118 124

	
119 125
    /// \e
120 126
    virtual void _setObjCoeffs(ExprIterator b, ExprIterator e);
121 127
    /// \e
122 128
    virtual void _getObjCoeffs(InsertIterator b) const;
123 129

	
124 130
    /// \e
125 131
    virtual void _setObjCoeff(int i, Value obj_coef);
126 132
    /// \e
127 133
    virtual Value _getObjCoeff(int i) const;
128 134

	
129 135
    ///\e
130 136
    virtual void _setSense(Sense);
131 137
    ///\e
132 138
    virtual Sense _getSense() const;
133 139

	
134 140
    ///\e
135 141
    virtual void _clear();
136 142

	
137 143
  };
138 144

	
139
  /// \brief Interface for a skeleton LP solver
145
  /// \brief Skeleton class for an LP solver interface
140 146
  ///
141
  /// This class implements an interface for a skeleton LP solver.
147
  ///This class does nothing, but it can serve as a skeleton when
148
  ///implementing an interface to new solvers.
149

	
142 150
  ///\ingroup lp_group
143
  class LpSkeleton : public SkeletonSolverBase, public LpSolver {
151
  class LpSkeleton : public LpSolver, public SkeletonSolverBase {
144 152
  public:
145
    LpSkeleton() : SkeletonSolverBase(), LpSolver() {}
146

	
153
    ///\e
154
    LpSkeleton() : LpSolver(), SkeletonSolverBase() {}
155
    ///\e
156
    virtual LpSkeleton* newSolver() const;
157
    ///\e
158
    virtual LpSkeleton* cloneSolver() const;
147 159
  protected:
148 160

	
149 161
    ///\e
150 162
    virtual SolveExitStatus _solve();
151 163

	
152 164
    ///\e
153 165
    virtual Value _getPrimal(int i) const;
154 166
    ///\e
155 167
    virtual Value _getDual(int i) const;
156 168

	
157 169
    ///\e
158 170
    virtual Value _getPrimalValue() const;
159 171

	
160 172
    ///\e
161 173
    virtual Value _getPrimalRay(int i) const;
162 174
    ///\e
163 175
    virtual Value _getDualRay(int i) const;
164 176

	
165 177
    ///\e
166 178
    virtual ProblemType _getPrimalType() const;
167 179
    ///\e
168 180
    virtual ProblemType _getDualType() const;
169 181

	
170 182
    ///\e
171 183
    virtual VarStatus _getColStatus(int i) const;
172 184
    ///\e
173 185
    virtual VarStatus _getRowStatus(int i) const;
174 186

	
175 187
    ///\e
176
    virtual LpSkeleton* _newSolver() const;
177
    ///\e
178
    virtual LpSkeleton* _cloneSolver() const;
179
    ///\e
180 188
    virtual const char* _solverName() const;
181 189

	
182 190
  };
183 191

	
184
  /// \brief Interface for a skeleton MIP solver
192
  /// \brief Skeleton class for a MIP solver interface
185 193
  ///
186
  /// This class implements an interface for a skeleton MIP solver.
194
  ///This class does nothing, but it can serve as a skeleton when
195
  ///implementing an interface to new solvers.
187 196
  ///\ingroup lp_group
188
  class MipSkeleton : public SkeletonSolverBase, public MipSolver {
197
  class MipSkeleton : public MipSolver, public SkeletonSolverBase {
189 198
  public:
190
    MipSkeleton() : SkeletonSolverBase(), MipSolver() {}
199
    ///\e
200
    MipSkeleton() : MipSolver(), SkeletonSolverBase() {}
201
    ///\e
202
    virtual MipSkeleton* newSolver() const;
203
    ///\e
204
    virtual MipSkeleton* cloneSolver() const;
191 205

	
192 206
  protected:
193 207
    ///\e
194 208

	
195 209
    ///\bug Wrong interface
196 210
    ///
197 211
    virtual SolveExitStatus _solve();
198 212

	
199 213
    ///\e
200 214

	
201 215
    ///\bug Wrong interface
202 216
    ///
203 217
    virtual Value _getSol(int i) const;
204 218

	
205 219
    ///\e
206 220

	
207 221
    ///\bug Wrong interface
208 222
    ///
209 223
    virtual Value _getSolValue() const;
210 224

	
211 225
    ///\e
212 226

	
213 227
    ///\bug Wrong interface
214 228
    ///
215 229
    virtual ProblemType _getType() const;
216 230

	
217 231
    ///\e
218
    virtual MipSkeleton* _newSolver() const;
219

	
220
    ///\e
221
    virtual MipSkeleton* _cloneSolver() const;
222
    ///\e
223 232
    virtual const char* _solverName() const;
224

	
225 233
  };
226 234

	
227 235
} //namespace lemon
228 236

	
229 237
#endif
Ignore white space 1536 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 <iostream>
20 20
#include <lemon/soplex.h>
21 21

	
22 22
#include <soplex.h>
23 23

	
24 24

	
25 25
///\file
26 26
///\brief Implementation of the LEMON-SOPLEX lp solver interface.
27 27
namespace lemon {
28 28

	
29 29
  SoplexLp::SoplexLp() {
30 30
    soplex = new soplex::SoPlex;
31 31
  }
32 32

	
33 33
  SoplexLp::~SoplexLp() {
34 34
    delete soplex;
35 35
  }
36 36

	
37 37
  SoplexLp::SoplexLp(const SoplexLp& lp) {
38 38
    rows = lp.rows;
39 39
    cols = lp.cols;
40 40

	
41 41
    soplex = new soplex::SoPlex;
42 42
    (*static_cast<soplex::SPxLP*>(soplex)) = *(lp.soplex);
43 43

	
44 44
    _col_names = lp._col_names;
45 45
    _col_names_ref = lp._col_names_ref;
46 46

	
47 47
    _row_names = lp._row_names;
48 48
    _row_names_ref = lp._row_names_ref;
49 49

	
50 50
  }
51 51

	
52 52
  void SoplexLp::_clear_temporals() {
53 53
    _primal_values.clear();
54 54
    _dual_values.clear();
55 55
  }
56 56

	
57
  SoplexLp* SoplexLp::_newSolver() const {
57
  SoplexLp* SoplexLp::newSolver() const {
58 58
    SoplexLp* newlp = new SoplexLp();
59 59
    return newlp;
60 60
  }
61 61

	
62
  SoplexLp* SoplexLp::_cloneSolver() const {
62
  SoplexLp* SoplexLp::cloneSolver() const {
63 63
    SoplexLp* newlp = new SoplexLp(*this);
64 64
    return newlp;
65 65
  }
66 66

	
67 67
  const char* SoplexLp::_solverName() const { return "SoplexLp"; }
68 68

	
69 69
  int SoplexLp::_addCol() {
70 70
    soplex::LPCol c;
71 71
    c.setLower(-soplex::infinity);
72 72
    c.setUpper(soplex::infinity);
73 73
    soplex->addCol(c);
74 74

	
75 75
    _col_names.push_back(std::string());
76 76

	
77 77
    return soplex->nCols() - 1;
78 78
  }
79 79

	
80 80
  int SoplexLp::_addRow() {
81 81
    soplex::LPRow r;
82 82
    r.setLhs(-soplex::infinity);
83 83
    r.setRhs(soplex::infinity);
84 84
    soplex->addRow(r);
85 85

	
86 86
    _row_names.push_back(std::string());
87 87

	
88 88
    return soplex->nRows() - 1;
89 89
  }
90 90

	
91 91

	
92 92
  void SoplexLp::_eraseCol(int i) {
93 93
    soplex->removeCol(i);
94 94
    _col_names_ref.erase(_col_names[i]);
95 95
    _col_names[i] = _col_names.back();
96 96
    _col_names_ref[_col_names.back()] = i;
97 97
    _col_names.pop_back();
98 98
  }
99 99

	
100 100
  void SoplexLp::_eraseRow(int i) {
101 101
    soplex->removeRow(i);
102 102
    _row_names_ref.erase(_row_names[i]);
103 103
    _row_names[i] = _row_names.back();
104 104
    _row_names_ref[_row_names.back()] = i;
105 105
    _row_names.pop_back();
106 106
  }
107 107

	
108 108
  void SoplexLp::_eraseColId(int i) {
109 109
    cols.eraseIndex(i);
110 110
    cols.relocateIndex(i, cols.maxIndex());
111 111
  }
112 112
  void SoplexLp::_eraseRowId(int i) {
113 113
    rows.eraseIndex(i);
114 114
    rows.relocateIndex(i, rows.maxIndex());
115 115
  }
116 116

	
117 117
  void SoplexLp::_getColName(int c, std::string &name) const {
118 118
    name = _col_names[c];
119 119
  }
120 120

	
121 121
  void SoplexLp::_setColName(int c, const std::string &name) {
122 122
    _col_names_ref.erase(_col_names[c]);
123 123
    _col_names[c] = name;
124 124
    if (!name.empty()) {
125 125
      _col_names_ref.insert(std::make_pair(name, c));
126 126
    }
127 127
  }
128 128

	
129 129
  int SoplexLp::_colByName(const std::string& name) const {
130 130
    std::map<std::string, int>::const_iterator it =
131 131
      _col_names_ref.find(name);
132 132
    if (it != _col_names_ref.end()) {
133 133
      return it->second;
134 134
    } else {
135 135
      return -1;
136 136
    }
137 137
  }
138 138

	
139 139
  void SoplexLp::_getRowName(int r, std::string &name) const {
140 140
    name = _row_names[r];
141 141
  }
142 142

	
143 143
  void SoplexLp::_setRowName(int r, const std::string &name) {
144 144
    _row_names_ref.erase(_row_names[r]);
145 145
    _row_names[r] = name;
146 146
    if (!name.empty()) {
147 147
      _row_names_ref.insert(std::make_pair(name, r));
148 148
    }
149 149
  }
150 150

	
151 151
  int SoplexLp::_rowByName(const std::string& name) const {
152 152
    std::map<std::string, int>::const_iterator it =
153 153
      _row_names_ref.find(name);
154 154
    if (it != _row_names_ref.end()) {
155 155
      return it->second;
156 156
    } else {
157 157
      return -1;
158 158
    }
159 159
  }
160 160

	
161 161

	
162 162
  void SoplexLp::_setRowCoeffs(int i, ExprIterator b, ExprIterator e) {
163 163
    for (int j = 0; j < soplex->nCols(); ++j) {
164 164
      soplex->changeElement(i, j, 0.0);
165 165
    }
166 166
    for(ExprIterator it = b; it != e; ++it) {
167 167
      soplex->changeElement(i, it->first, it->second);
168 168
    }
169 169
  }
170 170

	
171 171
  void SoplexLp::_getRowCoeffs(int i, InsertIterator b) const {
172 172
    const soplex::SVector& vec = soplex->rowVector(i);
173 173
    for (int k = 0; k < vec.size(); ++k) {
174 174
      *b = std::make_pair(vec.index(k), vec.value(k));
175 175
      ++b;
176 176
    }
177 177
  }
178 178

	
179 179
  void SoplexLp::_setColCoeffs(int j, ExprIterator b, ExprIterator e) {
180 180
    for (int i = 0; i < soplex->nRows(); ++i) {
181 181
      soplex->changeElement(i, j, 0.0);
182 182
    }
183 183
    for(ExprIterator it = b; it != e; ++it) {
184 184
      soplex->changeElement(it->first, j, it->second);
185 185
    }
186 186
  }
187 187

	
188 188
  void SoplexLp::_getColCoeffs(int i, InsertIterator b) const {
189 189
    const soplex::SVector& vec = soplex->colVector(i);
190 190
    for (int k = 0; k < vec.size(); ++k) {
191 191
      *b = std::make_pair(vec.index(k), vec.value(k));
192 192
      ++b;
193 193
    }
194 194
  }
195 195

	
196 196
  void SoplexLp::_setCoeff(int i, int j, Value value) {
197 197
    soplex->changeElement(i, j, value);
198 198
  }
199 199

	
200 200
  SoplexLp::Value SoplexLp::_getCoeff(int i, int j) const {
201 201
    return soplex->rowVector(i)[j];
202 202
  }
203 203

	
204 204
  void SoplexLp::_setColLowerBound(int i, Value value) {
205 205
    LEMON_ASSERT(value != INF, "Invalid bound");
206 206
    soplex->changeLower(i, value != -INF ? value : -soplex::infinity);
207 207
  }
208 208

	
209 209
  SoplexLp::Value SoplexLp::_getColLowerBound(int i) const {
210 210
    double value = soplex->lower(i);
211 211
    return value != -soplex::infinity ? value : -INF;
212 212
  }
213 213

	
214 214
  void SoplexLp::_setColUpperBound(int i, Value value) {
215 215
    LEMON_ASSERT(value != -INF, "Invalid bound");
216 216
    soplex->changeUpper(i, value != INF ? value : soplex::infinity);
217 217
  }
218 218

	
219 219
  SoplexLp::Value SoplexLp::_getColUpperBound(int i) const {
220 220
    double value = soplex->upper(i);
221 221
    return value != soplex::infinity ? value : INF;
222 222
  }
223 223

	
224 224
  void SoplexLp::_setRowLowerBound(int i, Value lb) {
225 225
    LEMON_ASSERT(lb != INF, "Invalid bound");
226 226
    soplex->changeRange(i, lb != -INF ? lb : -soplex::infinity, soplex->rhs(i));
227 227
  }
228 228

	
229 229
  SoplexLp::Value SoplexLp::_getRowLowerBound(int i) const {
230 230
    double res = soplex->lhs(i);
231 231
    return res == -soplex::infinity ? -INF : res;
232 232
  }
233 233

	
234 234
  void SoplexLp::_setRowUpperBound(int i, Value ub) {
235 235
    LEMON_ASSERT(ub != -INF, "Invalid bound");
236 236
    soplex->changeRange(i, soplex->lhs(i), ub != INF ? ub : soplex::infinity);
237 237
  }
238 238

	
239 239
  SoplexLp::Value SoplexLp::_getRowUpperBound(int i) const {
240 240
    double res = soplex->rhs(i);
241 241
    return res == soplex::infinity ? INF : res;
242 242
  }
243 243

	
244 244
  void SoplexLp::_setObjCoeffs(ExprIterator b, ExprIterator e) {
245 245
    for (int j = 0; j < soplex->nCols(); ++j) {
246 246
      soplex->changeObj(j, 0.0);
247 247
    }
248 248
    for (ExprIterator it = b; it != e; ++it) {
249 249
      soplex->changeObj(it->first, it->second);
250 250
    }
251 251
  }
252 252

	
253 253
  void SoplexLp::_getObjCoeffs(InsertIterator b) const {
254 254
    for (int j = 0; j < soplex->nCols(); ++j) {
255 255
      Value coef = soplex->obj(j);
256 256
      if (coef != 0.0) {
257 257
        *b = std::make_pair(j, coef);
258 258
        ++b;
259 259
      }
260 260
    }
261 261
  }
262 262

	
263 263
  void SoplexLp::_setObjCoeff(int i, Value obj_coef) {
264 264
    soplex->changeObj(i, obj_coef);
265 265
  }
266 266

	
267 267
  SoplexLp::Value SoplexLp::_getObjCoeff(int i) const {
268 268
    return soplex->obj(i);
269 269
  }
270 270

	
271 271
  SoplexLp::SolveExitStatus SoplexLp::_solve() {
272 272

	
273 273
    _clear_temporals();
274 274

	
275 275
    soplex::SPxSolver::Status status = soplex->solve();
276 276

	
277 277
    switch (status) {
278 278
    case soplex::SPxSolver::OPTIMAL:
279 279
    case soplex::SPxSolver::INFEASIBLE:
280 280
    case soplex::SPxSolver::UNBOUNDED:
281 281
      return SOLVED;
282 282
    default:
283 283
      return UNSOLVED;
284 284
    }
285 285
  }
286 286

	
287 287
  SoplexLp::Value SoplexLp::_getPrimal(int i) const {
288 288
    if (_primal_values.empty()) {
289 289
      _primal_values.resize(soplex->nCols());
290 290
      soplex::Vector pv(_primal_values.size(), &_primal_values.front());
291 291
      soplex->getPrimal(pv);
292 292
    }
293 293
    return _primal_values[i];
294 294
  }
295 295

	
296 296
  SoplexLp::Value SoplexLp::_getDual(int i) const {
297 297
    if (_dual_values.empty()) {
298 298
      _dual_values.resize(soplex->nRows());
299 299
      soplex::Vector dv(_dual_values.size(), &_dual_values.front());
300 300
      soplex->getDual(dv);
301 301
    }
302 302
    return _dual_values[i];
303 303
  }
304 304

	
305 305
  SoplexLp::Value SoplexLp::_getPrimalValue() const {
306 306
    return soplex->objValue();
307 307
  }
308 308

	
309 309
  SoplexLp::VarStatus SoplexLp::_getColStatus(int i) const {
310 310
    switch (soplex->getBasisColStatus(i)) {
311 311
    case soplex::SPxSolver::BASIC:
312 312
      return BASIC;
313 313
    case soplex::SPxSolver::ON_UPPER:
314 314
      return UPPER;
315 315
    case soplex::SPxSolver::ON_LOWER:
316 316
      return LOWER;
317 317
    case soplex::SPxSolver::FIXED:
318 318
      return FIXED;
319 319
    case soplex::SPxSolver::ZERO:
320 320
      return FREE;
321 321
    default:
322 322
      LEMON_ASSERT(false, "Wrong column status");
323 323
      return VarStatus();
324 324
    }
325 325
  }
326 326

	
327 327
  SoplexLp::VarStatus SoplexLp::_getRowStatus(int i) const {
328 328
    switch (soplex->getBasisRowStatus(i)) {
329 329
    case soplex::SPxSolver::BASIC:
330 330
      return BASIC;
331 331
    case soplex::SPxSolver::ON_UPPER:
332 332
      return UPPER;
333 333
    case soplex::SPxSolver::ON_LOWER:
334 334
      return LOWER;
335 335
    case soplex::SPxSolver::FIXED:
336 336
      return FIXED;
337 337
    case soplex::SPxSolver::ZERO:
338 338
      return FREE;
339 339
    default:
340 340
      LEMON_ASSERT(false, "Wrong row status");
341 341
      return VarStatus();
342 342
    }
343 343
  }
344 344

	
345 345
  SoplexLp::Value SoplexLp::_getPrimalRay(int i) const {
346 346
    if (_primal_ray.empty()) {
347 347
      _primal_ray.resize(soplex->nCols());
348 348
      soplex::Vector pv(_primal_ray.size(), &_primal_ray.front());
349 349
      soplex->getDualfarkas(pv);
350 350
    }
351 351
    return _primal_ray[i];
352 352
  }
353 353

	
354 354
  SoplexLp::Value SoplexLp::_getDualRay(int i) const {
355 355
    if (_dual_ray.empty()) {
356 356
      _dual_ray.resize(soplex->nRows());
357 357
      soplex::Vector dv(_dual_ray.size(), &_dual_ray.front());
358 358
      soplex->getDualfarkas(dv);
359 359
    }
360 360
    return _dual_ray[i];
361 361
  }
362 362

	
363 363
  SoplexLp::ProblemType SoplexLp::_getPrimalType() const {
364 364
    switch (soplex->status()) {
365 365
    case soplex::SPxSolver::OPTIMAL:
366 366
      return OPTIMAL;
367 367
    case soplex::SPxSolver::UNBOUNDED:
368 368
      return UNBOUNDED;
369 369
    case soplex::SPxSolver::INFEASIBLE:
370 370
      return INFEASIBLE;
371 371
    default:
372 372
      return UNDEFINED;
373 373
    }
374 374
  }
375 375

	
376 376
  SoplexLp::ProblemType SoplexLp::_getDualType() const {
377 377
    switch (soplex->status()) {
378 378
    case soplex::SPxSolver::OPTIMAL:
379 379
      return OPTIMAL;
380 380
    case soplex::SPxSolver::UNBOUNDED:
381 381
      return UNBOUNDED;
382 382
    case soplex::SPxSolver::INFEASIBLE:
383 383
      return INFEASIBLE;
384 384
    default:
385 385
      return UNDEFINED;
386 386
    }
387 387
  }
388 388

	
389 389
  void SoplexLp::_setSense(Sense sense) {
390 390
    switch (sense) {
391 391
    case MIN:
392 392
      soplex->changeSense(soplex::SPxSolver::MINIMIZE);
393 393
      break;
394 394
    case MAX:
395 395
      soplex->changeSense(soplex::SPxSolver::MAXIMIZE);
396 396
    }
397 397
  }
398 398

	
399 399
  SoplexLp::Sense SoplexLp::_getSense() const {
400 400
    switch (soplex->spxSense()) {
401 401
    case soplex::SPxSolver::MAXIMIZE:
402 402
      return MAX;
403 403
    case soplex::SPxSolver::MINIMIZE:
404 404
      return MIN;
405 405
    default:
406 406
      LEMON_ASSERT(false, "Wrong sense.");
407 407
      return SoplexLp::Sense();
408 408
    }
409 409
  }
410 410

	
411 411
  void SoplexLp::_clear() {
412 412
    soplex->clear();
413 413
    _col_names.clear();
414 414
    _col_names_ref.clear();
415 415
    _row_names.clear();
416 416
    _row_names_ref.clear();
417 417
    cols.clear();
418 418
    rows.clear();
419 419
    _clear_temporals();
420 420
  }
421 421

	
422 422
} //namespace lemon
423 423

	
Ignore white space 1536 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_SOPLEX_H
20 20
#define LEMON_SOPLEX_H
21 21

	
22 22
///\file
23 23
///\brief Header of the LEMON-SOPLEX lp solver interface.
24 24

	
25 25
#include <vector>
26 26
#include <string>
27 27

	
28 28
#include <lemon/lp_base.h>
29 29

	
30 30
// Forward declaration
31 31
namespace soplex {
32 32
  class SoPlex;
33 33
}
34 34

	
35 35
namespace lemon {
36 36

	
37 37
  /// \ingroup lp_group
38 38
  ///
39 39
  /// \brief Interface for the SOPLEX solver
40 40
  ///
41 41
  /// This class implements an interface for the SoPlex LP solver.
42 42
  /// The SoPlex library is an object oriented lp solver library
43 43
  /// developed at the Konrad-Zuse-Zentrum f�r Informationstechnik
44 44
  /// Berlin (ZIB). You can find detailed information about it at the
45 45
  /// <tt>http://soplex.zib.de</tt> address.
46 46
  class SoplexLp : public LpSolver {
47 47
  private:
48 48

	
49 49
    soplex::SoPlex* soplex;
50 50

	
51 51
    std::vector<std::string> _col_names;
52 52
    std::map<std::string, int> _col_names_ref;
53 53

	
54 54
    std::vector<std::string> _row_names;
55 55
    std::map<std::string, int> _row_names_ref;
56 56

	
57 57
  private:
58 58

	
59 59
    // these values cannot be retrieved element by element
60 60
    mutable std::vector<Value> _primal_values;
61 61
    mutable std::vector<Value> _dual_values;
62 62

	
63 63
    mutable std::vector<Value> _primal_ray;
64 64
    mutable std::vector<Value> _dual_ray;
65 65

	
66 66
    void _clear_temporals();
67 67

	
68 68
  public:
69 69

	
70 70
    /// \e
71 71
    SoplexLp();
72 72
    /// \e
73 73
    SoplexLp(const SoplexLp&);
74 74
    /// \e
75 75
    ~SoplexLp();
76
    /// \e
77
    virtual SoplexLp* newSolver() const;
78
    /// \e
79
    virtual SoplexLp* cloneSolver() const;
76 80

	
77 81
  protected:
78 82

	
79
    virtual SoplexLp* _newSolver() const;
80
    virtual SoplexLp* _cloneSolver() const;
81

	
82 83
    virtual const char* _solverName() const;
83 84

	
84 85
    virtual int _addCol();
85 86
    virtual int _addRow();
86 87

	
87 88
    virtual void _eraseCol(int i);
88 89
    virtual void _eraseRow(int i);
89 90

	
90 91
    virtual void _eraseColId(int i);
91 92
    virtual void _eraseRowId(int i);
92 93

	
93 94
    virtual void _getColName(int col, std::string& name) const;
94 95
    virtual void _setColName(int col, const std::string& name);
95 96
    virtual int _colByName(const std::string& name) const;
96 97

	
97 98
    virtual void _getRowName(int row, std::string& name) const;
98 99
    virtual void _setRowName(int row, const std::string& name);
99 100
    virtual int _rowByName(const std::string& name) const;
100 101

	
101 102
    virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e);
102 103
    virtual void _getRowCoeffs(int i, InsertIterator b) const;
103 104

	
104 105
    virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e);
105 106
    virtual void _getColCoeffs(int i, InsertIterator b) const;
106 107

	
107 108
    virtual void _setCoeff(int row, int col, Value value);
108 109
    virtual Value _getCoeff(int row, int col) const;
109 110

	
110 111
    virtual void _setColLowerBound(int i, Value value);
111 112
    virtual Value _getColLowerBound(int i) const;
112 113
    virtual void _setColUpperBound(int i, Value value);
113 114
    virtual Value _getColUpperBound(int i) const;
114 115

	
115 116
    virtual void _setRowLowerBound(int i, Value value);
116 117
    virtual Value _getRowLowerBound(int i) const;
117 118
    virtual void _setRowUpperBound(int i, Value value);
118 119
    virtual Value _getRowUpperBound(int i) const;
119 120

	
120 121
    virtual void _setObjCoeffs(ExprIterator b, ExprIterator e);
121 122
    virtual void _getObjCoeffs(InsertIterator b) const;
122 123

	
123 124
    virtual void _setObjCoeff(int i, Value obj_coef);
124 125
    virtual Value _getObjCoeff(int i) const;
125 126

	
126 127
    virtual void _setSense(Sense sense);
127 128
    virtual Sense _getSense() const;
128 129

	
129 130
    virtual SolveExitStatus _solve();
130 131
    virtual Value _getPrimal(int i) const;
131 132
    virtual Value _getDual(int i) const;
132 133

	
133 134
    virtual Value _getPrimalValue() const;
134 135

	
135 136
    virtual Value _getPrimalRay(int i) const;
136 137
    virtual Value _getDualRay(int i) const;
137 138

	
138 139
    virtual VarStatus _getColStatus(int i) const;
139 140
    virtual VarStatus _getRowStatus(int i) const;
140 141

	
141 142
    virtual ProblemType _getPrimalType() const;
142 143
    virtual ProblemType _getDualType() const;
143 144

	
144 145
    virtual void _clear();
145 146

	
146 147
  };
147 148

	
148 149
} //END OF NAMESPACE LEMON
149 150

	
150 151
#endif //LEMON_SOPLEX_H
151 152

	
Ignore white space 1536 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
    //Test for clone/new
201
    LP* lpnew = lp.newSolver();
202
    LP* lpclone = lp.cloneSolver();
203
    delete lpnew;
204
    delete lpclone;
200 205

	
201 206
  }
202 207

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

	
208 213
    e[p1]=2;
209 214
    e[p1]+=2;
210 215
    e[p1]-=2;
211 216

	
212 217
    e=p1;
213 218
    e=f;
214 219

	
215 220
    e+=p1;
216 221
    e+=f;
217 222

	
218 223
    e-=p1;
219 224
    e-=f;
220 225

	
221 226
    e*=2;
222 227
    e*=2.2;
223 228
    e/=2;
224 229
    e/=2.2;
225 230

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

	
236 241
}
237 242

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

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

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

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

	
255 261
void aTest(LpSolver & lp)
256 262
{
257 263
  typedef LpSolver LP;
258 264

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

	
261 267
  typedef LpSolver::Row Row;
262 268
  typedef LpSolver::Col Col;
263 269

	
264 270

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

	
268 274

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

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

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

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

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

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

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

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

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

	
329 335
  delete clp;
330 336

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

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

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

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

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

	
356 362
}
357 363

	
364
template<class LP>
365
void cloneTest()
366
{
367
  //Test for clone/new
368
  
369
  LP* lp = new LP();
370
  LP* lpnew = lp->newSolver();
371
  LP* lpclone = lp->cloneSolver();
372
  delete lp;
373
  delete lpnew;
374
  delete lpclone;
375
}
376

	
358 377
int main()
359 378
{
360 379
  LpSkeleton lp_skel;
361 380
  lpTest(lp_skel);
362 381

	
363 382
#ifdef HAVE_GLPK
364 383
  {
365 384
    GlpkLp lp_glpk1,lp_glpk2;
366 385
    lpTest(lp_glpk1);
367 386
    aTest(lp_glpk2);
387
    cloneTest<GlpkLp>();
368 388
  }
369 389
#endif
370 390

	
371 391
#ifdef HAVE_CPLEX
372 392
  try {
373 393
    CplexLp lp_cplex1,lp_cplex2;
374 394
    lpTest(lp_cplex1);
375 395
    aTest(lp_cplex2);
376 396
  } catch (CplexEnv::LicenseError& error) {
377 397
#ifdef LEMON_FORCE_CPLEX_CHECK
378 398
    check(false, error.what());
379 399
#else
380 400
    std::cerr << error.what() << std::endl;
381 401
    std::cerr << "Cplex license check failed, lp check skipped" << std::endl;
382 402
#endif
383 403
  }
404
    cloneTest<CplexLp>();
384 405
#endif
385 406

	
386 407
#ifdef HAVE_SOPLEX
387 408
  {
388 409
    SoplexLp lp_soplex1,lp_soplex2;
389 410
    lpTest(lp_soplex1);
390 411
    aTest(lp_soplex2);
412
    cloneTest<SoplexLp>();
391 413
  }
392 414
#endif
393 415

	
394 416
#ifdef HAVE_CLP
395 417
  {
396 418
    ClpLp lp_clp1,lp_clp2;
397 419
    lpTest(lp_clp1);
398 420
    aTest(lp_clp2);
421
    cloneTest<ClpLp>();
399 422
  }
400 423
#endif
401 424

	
402 425
  return 0;
403 426
}
Ignore white space 1536 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
template<class MIP>
110
void cloneTest()
111
{
112
  
113
  MIP* mip = new MIP();
114
  MIP* mipnew = mip->newSolver();
115
  MIP* mipclone = mip->cloneSolver();
116
  delete mip;
117
  delete mipnew;
118
  delete mipclone;
119
}
109 120

	
110 121
int main()
111 122
{
112 123

	
113 124
#ifdef HAVE_GLPK
114 125
  {
115 126
    GlpkMip mip1;
116 127
    aTest(mip1);
128
    cloneTest<GlpkMip>();
117 129
  }
118 130
#endif
119 131

	
120 132
#ifdef HAVE_CPLEX
121 133
  try {
122 134
    CplexMip mip2;
123 135
    aTest(mip2);
124 136
  } catch (CplexEnv::LicenseError& error) {
125 137
#ifdef LEMON_FORCE_CPLEX_CHECK
126 138
    check(false, error.what());
127 139
#else
128 140
    std::cerr << error.what() << std::endl;
129 141
    std::cerr << "Cplex license check failed, lp check skipped" << std::endl;
130 142
#endif
131 143
  }
144
  cloneTest<CplexMip>();
132 145
#endif
133 146

	
134 147
  return 0;
135 148

	
136 149
}
0 comments (0 inline)