src/work/marci/lp/lp_solver_base.h
author marci
Thu, 17 Feb 2005 15:14:13 +0000
changeset 1153 4b0468de3a31
parent 1152 1765ff9fefa1
child 1262 61f989e3e525
permissions -rw-r--r--
if you have a nuclear power plant and wanna compute small magic squares, then let's do it
     1 // -*- c++ -*-
     2 #ifndef LEMON_LP_SOLVER_BASE_H
     3 #define LEMON_LP_SOLVER_BASE_H
     4 
     5 ///\ingroup misc
     6 ///\file
     7 
     8 // #include <stdio.h>
     9 #include <stdlib.h>
    10 #include <iostream>
    11 #include <map>
    12 #include <limits>
    13 // #include <stdio>
    14 //#include <stdlib>
    15 extern "C" {
    16 #include "glpk.h"
    17 }
    18 
    19 #include <iostream>
    20 #include <vector>
    21 #include <string>
    22 #include <list>
    23 #include <memory>
    24 #include <utility>
    25 
    26 #include <lemon/invalid.h>
    27 #include <expression.h>
    28 //#include <stp.h>
    29 //#include <lemon/max_flow.h>
    30 //#include <augmenting_flow.h>
    31 //#include <iter_map.h>
    32 
    33 using std::cout;
    34 using std::cin;
    35 using std::endl;
    36 
    37 namespace lemon {
    38   
    39   /// \addtogroup misc
    40   /// @{
    41 
    42   /// \brief A partitioned vector with iterable classes.
    43   ///
    44   /// This class implements a container in which the data is stored in an 
    45   /// stl vector, the range is partitioned into sets and each set is 
    46   /// doubly linked in a list. 
    47   /// That is, each class is iterable by lemon iterators, and any member of 
    48   /// the vector can bo moved to an other class.
    49   template <typename T>
    50   class IterablePartition {
    51   protected:
    52     struct Node {
    53       T data;
    54       int prev; //invalid az -1
    55       int next; 
    56     };
    57     std::vector<Node> nodes;
    58     struct Tip {
    59       int first;
    60       int last;
    61     };
    62     std::vector<Tip> tips;
    63   public:
    64     /// The classes are indexed by integers from \c 0 to \c classNum()-1.
    65     int classNum() const { return tips.size(); }
    66     /// This lemon style iterator iterates through a class. 
    67     class Class;
    68     /// Constructor. The number of classes is to be given which is fixed 
    69     /// over the life of the container. 
    70     /// The partition classes are indexed from 0 to class_num-1. 
    71     IterablePartition(int class_num) { 
    72       for (int i=0; i<class_num; ++i) {
    73 	Tip t;
    74 	t.first=t.last=-1;
    75 	tips.push_back(t);
    76       }
    77     }
    78   protected:
    79     void befuz(Class it, int class_id) {
    80       if (tips[class_id].first==-1) {
    81 	if (tips[class_id].last==-1) {
    82 	  nodes[it.i].prev=nodes[it.i].next=-1;
    83 	  tips[class_id].first=tips[class_id].last=it.i;
    84 	}
    85       } else {
    86 	nodes[it.i].prev=tips[class_id].last;
    87 	nodes[it.i].next=-1;
    88 	nodes[tips[class_id].last].next=it.i;
    89 	tips[class_id].last=it.i;
    90       }
    91     }
    92     void kifuz(Class it, int class_id) {
    93       if (tips[class_id].first==it.i) {
    94 	if (tips[class_id].last==it.i) {
    95 	  tips[class_id].first=tips[class_id].last=-1;
    96 	} else {
    97 	  tips[class_id].first=nodes[it.i].next;
    98 	  nodes[nodes[it.i].next].prev=-1;
    99 	}
   100       } else {
   101 	if (tips[class_id].last==it.i) {
   102 	  tips[class_id].last=nodes[it.i].prev;
   103 	  nodes[nodes[it.i].prev].next=-1;
   104 	} else {
   105 	  nodes[nodes[it.i].next].prev=nodes[it.i].prev;
   106 	  nodes[nodes[it.i].prev].next=nodes[it.i].next;
   107 	}
   108       }
   109     }
   110   public:
   111     /// A new element with data \c t is pushed into the vector and into class 
   112     /// \c class_id.
   113     Class push_back(const T& t, int class_id) { 
   114       Node n;
   115       n.data=t;
   116       nodes.push_back(n);
   117       int i=nodes.size()-1;
   118       befuz(i, class_id);
   119       return i;
   120     }
   121     /// A member is moved to an other class.
   122     void set(Class it, int old_class_id, int new_class_id) {
   123       kifuz(it.i, old_class_id);
   124       befuz(it.i, new_class_id);
   125     }
   126     /// Returns the data pointed by \c it.
   127     T& operator[](Class it) { return nodes[it.i].data; }
   128     /// Returns the data pointed by \c it.
   129     const T& operator[](Class it) const { return nodes[it.i].data; }
   130     ///.
   131     class Class {
   132       friend class IterablePartition;
   133     protected:
   134       int i;
   135     public:
   136       /// Default constructor.
   137       Class() { }
   138       /// This constructor constructs an iterator which points
   139       /// to the member of th container indexed by the integer _i.
   140       Class(const int& _i) : i(_i) { }
   141       /// Invalid constructor.
   142       Class(const Invalid&) : i(-1) { }
   143       friend bool operator<(const Class& x, const Class& y);
   144       friend std::ostream& operator<<(std::ostream& os, 
   145 				      const Class& it);
   146       bool operator==(const Class& node) const {return i == node.i;}
   147       bool operator!=(const Class& node) const {return i != node.i;}
   148     };
   149     friend bool operator<(const Class& x, const Class& y) {
   150       return (x.i < y.i);
   151     }
   152     friend std::ostream& operator<<(std::ostream& os, 
   153 				    const Class& it) {
   154       os << it.i;
   155       return os;
   156     }
   157     /// First member of class \c class_id.
   158     Class& first(Class& it, int class_id) const {
   159       it.i=tips[class_id].first;
   160       return it;
   161     }
   162     /// Next member.
   163     Class& next(Class& it) const {
   164       it.i=nodes[it.i].next;
   165       return it;
   166     }
   167     /// True iff the iterator is valid.
   168     bool valid(const Class& it) const { return it.i!=-1; }
   169 
   170     class ClassIt : public Class {
   171       const IterablePartition* iterable_partition;
   172     public:
   173       ClassIt() { }
   174       ClassIt(Invalid i) : Class(i) { }
   175       ClassIt(const IterablePartition& _iterable_partition, 
   176 	      const int& i) : iterable_partition(&_iterable_partition) {
   177         _iterable_partition.first(*this, i);
   178       }
   179       ClassIt(const IterablePartition& _iterable_partition, 
   180 	      const Class& _class) : 
   181 	Class(_class), iterable_partition(&_iterable_partition) { }
   182       ClassIt& operator++() {
   183         iterable_partition->next(*this);
   184         return *this;
   185       }
   186     };
   187 
   188   };
   189 
   190 
   191   /*! \e
   192     \todo kellenene uj iterable structure bele, mert ez nem az igazi
   193     \todo A[x,y]-t cserel. Jobboldal, baloldal csere.
   194     \todo LEKERDEZESEK!!!
   195     \todo DOKSI!!!! Doxygen group!!!
   196     The aim of this class is to give a general surface to different 
   197     solvers, i.e. it makes possible to write algorithms using LP's, 
   198     in which the solver can be changed to an other one easily.
   199     \nosubgrouping
   200   */
   201   template <typename _Value>
   202   class LPSolverBase {
   203     
   204     /*! @name Uncategorized functions and types (public members)
   205     */
   206     //@{
   207   public:
   208 
   209     //UNCATEGORIZED
   210 
   211     /// \e
   212     typedef IterablePartition<int> Rows;
   213     /// \e
   214     typedef IterablePartition<int> Cols;
   215     /// \e
   216     typedef _Value Value;
   217     /// \e
   218     typedef Rows::Class Row;
   219     /// \e
   220     typedef Cols::Class Col;
   221   public:
   222     /// \e
   223     IterablePartition<int> row_iter_map;
   224     /// \e
   225     IterablePartition<int> col_iter_map;
   226     /// \e
   227     std::vector<Row> int_row_map;
   228     /// \e
   229     std::vector<Col> int_col_map;
   230     /// \e
   231     const int VALID_CLASS;
   232     /// \e
   233     const int INVALID_CLASS;
   234     /// \e 
   235     static const _Value INF;
   236   public:
   237     /// \e
   238     LPSolverBase() : row_iter_map(2), 
   239 		     col_iter_map(2), 
   240 		     VALID_CLASS(0), INVALID_CLASS(1) { }
   241     /// \e
   242     virtual ~LPSolverBase() { }
   243     //@}
   244 
   245     /*! @name Medium level interface (public members)
   246       These functions appear in the low level and also in the high level 
   247       interfaces thus these each of these functions have to be implemented 
   248       only once in the different interfaces.
   249       This means that these functions have to be reimplemented for all of the 
   250       different lp solvers. These are basic functions, and have the same 
   251       parameter lists in the low and high level interfaces. 
   252     */
   253     //@{
   254   public:
   255 
   256     //UNCATEGORIZED FUNCTIONS
   257 
   258     /// \e
   259     virtual void setMinimize() = 0;
   260     /// \e
   261     virtual void setMaximize() = 0;
   262 
   263     //SOLVER FUNCTIONS
   264 
   265     /// \e
   266     virtual void solveSimplex() = 0;
   267     /// \e
   268     virtual void solvePrimalSimplex() = 0;
   269     /// \e
   270     virtual void solveDualSimplex() = 0;
   271 
   272     //SOLUTION RETRIEVING
   273 
   274     /// \e
   275     virtual _Value getObjVal() = 0;
   276 
   277     //OTHER FUNCTIONS
   278 
   279     /// \e
   280     virtual int rowNum() const = 0;
   281     /// \e
   282     virtual int colNum() const = 0;
   283     /// \e
   284     virtual int warmUp() = 0;
   285     /// \e
   286     virtual void printWarmUpStatus(int i) = 0;
   287     /// \e
   288     virtual int getPrimalStatus() = 0;
   289     /// \e
   290     virtual void printPrimalStatus(int i) = 0;
   291     /// \e
   292     virtual int getDualStatus() = 0;
   293     /// \e
   294     virtual void printDualStatus(int i) = 0;
   295     /// Returns the status of the slack variable assigned to row \c row.
   296     virtual int getRowStat(const Row& row) = 0;
   297     /// \e
   298     virtual void printRowStatus(int i) = 0;
   299     /// Returns the status of the variable assigned to column \c col.
   300     virtual int getColStat(const Col& col) = 0;
   301     /// \e
   302     virtual void printColStatus(int i) = 0;
   303 
   304     //@}
   305 
   306     /*! @name Low level interface (protected members)
   307       Problem manipulating functions in the low level interface
   308     */
   309     //@{
   310   protected:
   311 
   312     //MATRIX MANIPULATING FUNCTIONS
   313 
   314     /// \e
   315     virtual int _addCol() = 0;
   316     /// \e
   317     virtual int _addRow() = 0;
   318     /// \e
   319     virtual void _eraseCol(int i) = 0;
   320     /// \e
   321     virtual void _eraseRow(int i) = 0;
   322     /// \e
   323     virtual void _setRowCoeffs(int i, 
   324 			       const std::vector<std::pair<int, _Value> >& coeffs) = 0;
   325     /// \e
   326     /// This routine modifies \c coeffs only by the \c push_back method.
   327     virtual void _getRowCoeffs(int i, 
   328 			       std::vector<std::pair<int, _Value> >& coeffs) = 0;
   329     /// \e
   330     virtual void _setColCoeffs(int i, 
   331 			       const std::vector<std::pair<int, _Value> >& coeffs) = 0;
   332     /// \e
   333     /// This routine modifies \c coeffs only by the \c push_back method.
   334     virtual void _getColCoeffs(int i, 
   335 			       std::vector<std::pair<int, _Value> >& coeffs) = 0;
   336     /// \e
   337     virtual void _setCoeff(int col, int row, _Value value) = 0;
   338     /// \e
   339     virtual _Value _getCoeff(int col, int row) = 0;
   340     //  public:
   341     //    /// \e
   342     //    enum Bound { FREE, LOWER, UPPER, DOUBLE, FIXED };
   343   protected:
   344     /// \e
   345     /// The lower bound of a variable (column) have to be given by an 
   346     /// extended number of type _Value, i.e. a finite number of type 
   347     /// _Value or -INF.
   348     virtual void _setColLowerBound(int i, _Value value) = 0;
   349     /// \e
   350     /// The lower bound of a variable (column) is an 
   351     /// extended number of type _Value, i.e. a finite number of type 
   352     /// _Value or -INF.
   353     virtual _Value _getColLowerBound(int i) = 0;
   354     /// \e
   355     /// The upper bound of a variable (column) have to be given by an 
   356     /// extended number of type _Value, i.e. a finite number of type 
   357     /// _Value or INF.
   358     virtual void _setColUpperBound(int i, _Value value) = 0;
   359     /// \e
   360     /// The upper bound of a variable (column) is an 
   361     /// extended number of type _Value, i.e. a finite number of type 
   362     /// _Value or INF.
   363     virtual _Value _getColUpperBound(int i) = 0;
   364     /// \e
   365     /// The lower bound of a linear expression (row) have to be given by an 
   366     /// extended number of type _Value, i.e. a finite number of type 
   367     /// _Value or -INF.
   368     virtual void _setRowLowerBound(int i, _Value value) = 0;
   369     /// \e
   370     /// The lower bound of a linear expression (row) is an 
   371     /// extended number of type _Value, i.e. a finite number of type 
   372     /// _Value or -INF.
   373     virtual _Value _getRowLowerBound(int i) = 0;
   374     /// \e
   375     /// The upper bound of a linear expression (row) have to be given by an 
   376     /// extended number of type _Value, i.e. a finite number of type 
   377     /// _Value or INF.
   378     virtual void _setRowUpperBound(int i, _Value value) = 0;
   379     /// \e
   380     /// The upper bound of a linear expression (row) is an 
   381     /// extended number of type _Value, i.e. a finite number of type 
   382     /// _Value or INF.
   383     virtual _Value _getRowUpperBound(int i) = 0;
   384     /// \e
   385     virtual void _setObjCoeff(int i, _Value obj_coef) = 0;
   386     /// \e
   387     virtual _Value _getObjCoeff(int i) = 0;
   388     
   389     //SOLUTION RETRIEVING
   390 
   391     /// \e
   392     virtual _Value _getPrimal(int i) = 0;
   393     //@}
   394     
   395     /*! @name High level interface (public members)
   396       Problem manipulating functions in the high level interface
   397     */
   398     //@{
   399   public:
   400 
   401     //MATRIX MANIPULATING FUNCTIONS
   402 
   403     /// \e
   404     Col addCol() {
   405       int i=_addCol();  
   406       Col col;
   407       col_iter_map.first(col, INVALID_CLASS);
   408       if (col_iter_map.valid(col)) { //van hasznalhato hely
   409 	col_iter_map.set(col, INVALID_CLASS, VALID_CLASS);
   410 	col_iter_map[col]=i;
   411       } else { //a cucc vegere kell inzertalni mert nincs szabad hely
   412 	col=col_iter_map.push_back(i, VALID_CLASS);
   413       }
   414       int_col_map.push_back(col);
   415       return col;
   416     }
   417     /// \e
   418     Row addRow() {
   419       int i=_addRow();
   420       Row row;
   421       row_iter_map.first(row, INVALID_CLASS);
   422       if (row_iter_map.valid(row)) { //van hasznalhato hely
   423 	row_iter_map.set(row, INVALID_CLASS, VALID_CLASS);
   424 	row_iter_map[row]=i;
   425       } else { //a cucc vegere kell inzertalni mert nincs szabad hely
   426 	row=row_iter_map.push_back(i, VALID_CLASS);
   427       }
   428       int_row_map.push_back(row);
   429       return row;
   430     }
   431     /// \e
   432     void eraseCol(const Col& col) {
   433       col_iter_map.set(col, VALID_CLASS, INVALID_CLASS);
   434       int cols[2];
   435       cols[1]=col_iter_map[col];
   436       _eraseCol(cols[1]);
   437       col_iter_map[col]=0; //glpk specifikus, de kell ez??
   438       Col it;
   439       for (col_iter_map.first(it, VALID_CLASS); 
   440 	   col_iter_map.valid(it); col_iter_map.next(it)) {
   441 	if (col_iter_map[it]>cols[1]) --col_iter_map[it];
   442       }
   443       int_col_map.erase(int_col_map.begin()+cols[1]);
   444     }
   445     /// \e
   446     void eraseRow(const Row& row) {
   447       row_iter_map.set(row, VALID_CLASS, INVALID_CLASS);
   448       int rows[2];
   449       rows[1]=row_iter_map[row];
   450       _eraseRow(rows[1]);
   451       row_iter_map[row]=0; //glpk specifikus, de kell ez??
   452       Row it;
   453       for (row_iter_map.first(it, VALID_CLASS); 
   454 	   row_iter_map.valid(it); row_iter_map.next(it)) {
   455 	if (row_iter_map[it]>rows[1]) --row_iter_map[it];
   456       }
   457       int_row_map.erase(int_row_map.begin()+rows[1]);
   458     }
   459     /// \e
   460     void setCoeff(Col col, Row row, _Value value) {
   461       _setCoeff(col_iter_map[col], row_iter_map[row], value);
   462     }
   463     /// \e
   464     _Value getCoeff(Col col, Row row) {
   465       return _getCoeff(col_iter_map[col], row_iter_map[row], value);
   466     }
   467     /// \e
   468     void setColLowerBound(Col col, _Value lo) {
   469       _setColLowerBound(col_iter_map[col], lo);
   470     }
   471     /// \e
   472     _Value getColLowerBound(Col col) {
   473       return _getColLowerBound(col_iter_map[col]);
   474     }
   475     /// \e
   476     void setColUpperBound(Col col, _Value up) {
   477       _setColUpperBound(col_iter_map[col], up);
   478     }
   479     /// \e
   480     _Value getColUpperBound(Col col) {      
   481       return _getColUpperBound(col_iter_map[col]);
   482     }
   483     /// \e
   484     void setRowLowerBound(Row row, _Value lo) {
   485       _setRowLowerBound(row_iter_map[row], lo);
   486     }
   487     /// \e
   488     _Value getRowLowerBound(Row row) {
   489       return _getRowLowerBound(row_iter_map[row]);
   490     }
   491     /// \e
   492     void setRowUpperBound(Row row, _Value up) {
   493       _setRowUpperBound(row_iter_map[row], up);
   494     }
   495     /// \e
   496     _Value getRowUpperBound(Row row) {      
   497       return _getRowUpperBound(row_iter_map[row]);
   498     }
   499     /// \e
   500     void setObjCoeff(const Col& col, _Value obj_coef) {
   501       _setObjCoeff(col_iter_map[col], obj_coef);
   502     }
   503     /// \e
   504     _Value getObjCoeff(const Col& col) {
   505       return _getObjCoeff(col_iter_map[col]);
   506     }
   507 
   508     //SOLUTION RETRIEVING FUNCTIONS
   509 
   510     /// \e
   511     _Value getPrimal(const Col& col) {
   512       return _getPrimal(col_iter_map[col]);
   513     }    
   514 
   515     //@}
   516 
   517     /*! @name User friend interface
   518       Problem manipulating functions in the user friend interface
   519     */
   520     //@{
   521 
   522     //EXPRESSION TYPES
   523 
   524     /// \e
   525     typedef Expr<Col, _Value> Expression;
   526     /// \e
   527     typedef Expr<Row, _Value> DualExpression;
   528     /// \e
   529     typedef Constr<Col, _Value> Constraint;
   530 
   531     //MATRIX MANIPULATING FUNCTIONS
   532 
   533     /// \e
   534     void setRowCoeffs(Row row, const Expression& expr) {
   535       std::vector<std::pair<int, _Value> > row_coeffs;
   536       for(typename Expression::Data::const_iterator i=expr.data.begin(); 
   537 	  i!=expr.data.end(); ++i) {
   538 	row_coeffs.push_back(std::make_pair
   539 			     (col_iter_map[(*i).first], (*i).second));
   540       }
   541       _setRowCoeffs(row_iter_map[row], row_coeffs);
   542     }
   543     /// \e 
   544     void setRow(Row row, const Constraint& constr) {
   545       setRowCoeffs(row, constr.expr);
   546       setRowLowerBound(row, constr.lo);
   547       setRowUpperBound(row, constr.up);
   548     }
   549     /// \e 
   550     Row addRow(const Constraint& constr) {
   551       Row row=addRow();
   552       setRowCoeffs(row, constr.expr);
   553       setRowLowerBound(row, constr.lo);
   554       setRowUpperBound(row, constr.up);
   555       return row;
   556     }
   557     /// \e
   558     /// This routine modifies \c expr by only adding to it.
   559     void getRowCoeffs(Row row, Expression& expr) {
   560       std::vector<std::pair<int, _Value> > row_coeffs;
   561       _getRowCoeffs(row_iter_map[row], row_coeffs);
   562       for(typename std::vector<std::pair<int, _Value> >::const_iterator 
   563  	    i=row_coeffs.begin(); i!=row_coeffs.end(); ++i) {
   564  	expr+= (*i).second*int_col_map[(*i).first];
   565       }
   566     }
   567     /// \e
   568     void setColCoeffs(Col col, const DualExpression& expr) {
   569       std::vector<std::pair<int, _Value> > col_coeffs;
   570       for(typename DualExpression::Data::const_iterator i=expr.data.begin(); 
   571 	  i!=expr.data.end(); ++i) {
   572 	col_coeffs.push_back(std::make_pair
   573 			     (row_iter_map[(*i).first], (*i).second));
   574       }
   575       _setColCoeffs(col_iter_map[col], col_coeffs);
   576     }
   577     /// \e
   578     /// This routine modifies \c expr by only adding to it.
   579     void getColCoeffs(Col col, DualExpression& expr) {
   580       std::vector<std::pair<int, _Value> > col_coeffs;
   581       _getColCoeffs(col_iter_map[col], col_coeffs);
   582       for(typename std::vector<std::pair<int, _Value> >::const_iterator 
   583  	    i=col_coeffs.begin(); i!=col_coeffs.end(); ++i) {
   584  	expr+= (*i).second*int_row_map[(*i).first];
   585       }
   586     }
   587     /// \e
   588     void setObjCoeffs(const Expression& expr) {
   589       // writing zero everywhere
   590       for(Cols::ClassIt it(col_iter_map, VALID_CLASS); it!=INVALID; ++it)
   591 	setObjCoeff(it, 0.0);
   592       // writing the data needed
   593       for(typename Expression::Data::const_iterator i=expr.data.begin(); 
   594 	  i!=expr.data.end(); ++i) {
   595 	setObjCoeff((*i).first, (*i).second);
   596       }
   597     }
   598     /// \e
   599     /// This routine modifies \c expr by only adding to it.
   600     void getObjCoeffs(Expression& expr) {
   601       for(Cols::ClassIt it(col_iter_map, VALID_CLASS); it!=INVALID; ++it)
   602 	expr+=getObjCoeff(it)*it;
   603     }
   604     //@}
   605 
   606 
   607     /*! @name MIP functions and types (public members)
   608     */
   609     //@{
   610   public:
   611     /// \e
   612     virtual void solveBandB() = 0;
   613     /// \e
   614     virtual void setLP() = 0;
   615     /// \e
   616     virtual void setMIP() = 0;
   617   protected:
   618    /// \e
   619     virtual void _setColCont(int i) = 0;
   620     /// \e
   621     virtual void _setColInt(int i) = 0;
   622     /// \e
   623     virtual _Value _getMIPPrimal(int i) = 0;
   624   public:
   625     /// \e
   626     void setColCont(Col col) {
   627       _setColCont(col_iter_map[col]);
   628     }
   629     /// \e
   630     void setColInt(Col col) {
   631       _setColInt(col_iter_map[col]);
   632     }
   633     /// \e
   634     _Value getMIPPrimal(Col col) {
   635       return _getMIPPrimal(col_iter_map[col]);
   636     }
   637     //@}
   638   };
   639   
   640   template <typename _Value>
   641   const _Value LPSolverBase<_Value>::INF=std::numeric_limits<_Value>::infinity();
   642 
   643 
   644   /// \brief Wrapper for GLPK solver
   645   /// 
   646   /// This class implements a lemon wrapper for GLPK.
   647   class LPGLPK : public LPSolverBase<double> {
   648   public:
   649     typedef LPSolverBase<double> Parent;
   650 
   651   public:
   652     /// \e
   653     LPX* lp;
   654 
   655   public:
   656     /// \e
   657     LPGLPK() : Parent(), 
   658 			lp(lpx_create_prob()) {
   659       int_row_map.push_back(Row());
   660       int_col_map.push_back(Col());
   661       lpx_set_int_parm(lp, LPX_K_DUAL, 1);
   662     }
   663     /// \e
   664     ~LPGLPK() {
   665       lpx_delete_prob(lp);
   666     }
   667 
   668     //MATRIX INDEPEDENT MANIPULATING FUNCTIONS
   669 
   670     /// \e
   671     void setMinimize() { 
   672       lpx_set_obj_dir(lp, LPX_MIN);
   673     }
   674     /// \e
   675     void setMaximize() { 
   676       lpx_set_obj_dir(lp, LPX_MAX);
   677     }
   678 
   679     //LOW LEVEL INTERFACE, MATRIX MANIPULATING FUNCTIONS
   680 
   681   protected:
   682     /// \e
   683     int _addCol() { 
   684       int i=lpx_add_cols(lp, 1);
   685       _setColLowerBound(i, -INF);
   686       _setColUpperBound(i, INF);
   687       return i;
   688     }
   689     /// \e
   690     int _addRow() { 
   691       int i=lpx_add_rows(lp, 1);
   692       return i;
   693     }
   694     /// \e
   695     virtual void _setRowCoeffs(int i, 
   696 			       const std::vector<std::pair<int, double> >& coeffs) {
   697       int mem_length=1+colNum();
   698       int* indices = new int[mem_length];
   699       double* doubles = new double[mem_length];
   700       int length=0;
   701       for (std::vector<std::pair<int, double> >::
   702 	     const_iterator it=coeffs.begin(); it!=coeffs.end(); ++it) {
   703 	++length;
   704 	indices[length]=it->first;
   705 	doubles[length]=it->second;
   706       }
   707       lpx_set_mat_row(lp, i, length, indices, doubles);
   708       delete [] indices;
   709       delete [] doubles;
   710     }
   711     /// \e
   712     virtual void _getRowCoeffs(int i, 
   713 			       std::vector<std::pair<int, double> >& coeffs) {
   714       int mem_length=1+colNum();
   715       int* indices = new int[mem_length];
   716       double* doubles = new double[mem_length];
   717       int length=lpx_get_mat_row(lp, i, indices, doubles);
   718       for (int i=1; i<=length; ++i) {
   719 	coeffs.push_back(std::make_pair(indices[i], doubles[i]));
   720       }
   721       delete [] indices;
   722       delete [] doubles;
   723     }
   724     /// \e
   725     virtual void _setColCoeffs(int i, 
   726 			       const std::vector<std::pair<int, double> >& coeffs) {
   727       int mem_length=1+rowNum();
   728       int* indices = new int[mem_length];
   729       double* doubles = new double[mem_length];
   730       int length=0;
   731       for (std::vector<std::pair<int, double> >::
   732 	     const_iterator it=coeffs.begin(); it!=coeffs.end(); ++it) {
   733 	++length;
   734 	indices[length]=it->first;
   735 	doubles[length]=it->second;
   736       }
   737       lpx_set_mat_col(lp, i, length, indices, doubles);
   738       delete [] indices;
   739       delete [] doubles;
   740     }
   741     /// \e
   742     virtual void _getColCoeffs(int i, 
   743 			       std::vector<std::pair<int, double> >& coeffs) {
   744       int mem_length=1+rowNum();
   745       int* indices = new int[mem_length];
   746       double* doubles = new double[mem_length];
   747       int length=lpx_get_mat_col(lp, i, indices, doubles);
   748       for (int i=1; i<=length; ++i) {
   749 	coeffs.push_back(std::make_pair(indices[i], doubles[i]));
   750       }
   751       delete [] indices;
   752       delete [] doubles;
   753     }
   754     /// \e
   755     virtual void _eraseCol(int i) {
   756       int cols[2];
   757       cols[1]=i;
   758       lpx_del_cols(lp, 1, cols);
   759     }
   760     virtual void _eraseRow(int i) {
   761       int rows[2];
   762       rows[1]=i;
   763       lpx_del_rows(lp, 1, rows);
   764     }
   765     void _setCoeff(int col, int row, double value) {
   766       /// FIXME not yet implemented
   767     }
   768     double _getCoeff(int col, int row) {
   769       /// FIXME not yet implemented
   770       return 0.0;
   771     }
   772     virtual void _setColLowerBound(int i, double lo) {
   773       if (lo==INF) {
   774 	//FIXME error
   775       }
   776       int b=lpx_get_col_type(lp, i);
   777       double up=lpx_get_col_ub(lp, i);	
   778       if (lo==-INF) {
   779 	switch (b) {
   780 	case LPX_FR:
   781 	case LPX_LO:
   782 	  lpx_set_col_bnds(lp, i, LPX_FR, lo, up);
   783 	  break;
   784 	case LPX_UP:
   785 	  break;
   786 	case LPX_DB:
   787 	case LPX_FX:
   788 	  lpx_set_col_bnds(lp, i, LPX_UP, lo, up);
   789 	  break;
   790 	default: ;
   791 	  //FIXME error
   792 	}
   793       } else {
   794 	switch (b) {
   795 	case LPX_FR:
   796 	case LPX_LO:
   797 	  lpx_set_col_bnds(lp, i, LPX_LO, lo, up);
   798 	  break;
   799 	case LPX_UP:	  
   800 	case LPX_DB:
   801 	case LPX_FX:
   802 	  if (lo==up) 
   803 	    lpx_set_col_bnds(lp, i, LPX_FX, lo, up);
   804 	  else 
   805 	    lpx_set_col_bnds(lp, i, LPX_DB, lo, up);
   806 	  break;
   807 	default: ;
   808 	  //FIXME error
   809 	}
   810       }
   811     }
   812     virtual double _getColLowerBound(int i) {
   813       int b=lpx_get_col_type(lp, i);
   814       switch (b) {
   815       case LPX_FR:
   816 	return -INF;
   817       case LPX_LO:
   818 	return lpx_get_col_lb(lp, i);
   819       case LPX_UP:
   820 	return -INF;
   821       case LPX_DB:
   822       case LPX_FX:
   823 	return lpx_get_col_lb(lp, i);
   824       default: ;
   825 	//FIXME error
   826 	return 0.0;
   827       }
   828     }
   829     virtual void _setColUpperBound(int i, double up) {
   830       if (up==-INF) {
   831 	//FIXME error
   832       }
   833       int b=lpx_get_col_type(lp, i);
   834       double lo=lpx_get_col_lb(lp, i);
   835       if (up==INF) {
   836 	switch (b) {
   837 	case LPX_FR:
   838 	case LPX_LO:
   839 	  break;
   840 	case LPX_UP:
   841 	  lpx_set_col_bnds(lp, i, LPX_FR, lo, up);
   842 	  break;
   843 	case LPX_DB:
   844 	case LPX_FX:
   845 	  lpx_set_col_bnds(lp, i, LPX_LO, lo, up);
   846 	  break;
   847 	default: ;
   848 	  //FIXME error
   849 	}
   850       } else {
   851 	switch (b) {
   852 	case LPX_FR:
   853 	  lpx_set_col_bnds(lp, i, LPX_UP, lo, up);
   854 	case LPX_LO:
   855 	  if (lo==up) 
   856 	    lpx_set_col_bnds(lp, i, LPX_FX, lo, up);
   857 	  else
   858 	    lpx_set_col_bnds(lp, i, LPX_DB, lo, up);
   859 	  break;
   860 	case LPX_UP:
   861 	  lpx_set_col_bnds(lp, i, LPX_UP, lo, up);
   862 	  break;
   863 	case LPX_DB:
   864 	case LPX_FX:
   865 	  if (lo==up) 
   866 	    lpx_set_col_bnds(lp, i, LPX_FX, lo, up);
   867 	  else 
   868 	    lpx_set_col_bnds(lp, i, LPX_DB, lo, up);
   869 	  break;
   870 	default: ;
   871 	  //FIXME error
   872 	}
   873       }
   874     }
   875     virtual double _getColUpperBound(int i) {
   876       int b=lpx_get_col_type(lp, i);
   877       switch (b) {
   878       case LPX_FR:
   879       case LPX_LO:
   880 	return INF;
   881       case LPX_UP:
   882       case LPX_DB:
   883       case LPX_FX:
   884 	return lpx_get_col_ub(lp, i);
   885       default: ;
   886 	//FIXME error
   887 	return 0.0;
   888       }
   889     }
   890     virtual void _setRowLowerBound(int i, double lo) {
   891       if (lo==INF) {
   892 	//FIXME error
   893       }
   894       int b=lpx_get_row_type(lp, i);
   895       double up=lpx_get_row_ub(lp, i);	
   896       if (lo==-INF) {
   897 	switch (b) {
   898 	case LPX_FR:
   899 	case LPX_LO:
   900 	  lpx_set_row_bnds(lp, i, LPX_FR, lo, up);
   901 	  break;
   902 	case LPX_UP:
   903 	  break;
   904 	case LPX_DB:
   905 	case LPX_FX:
   906 	  lpx_set_row_bnds(lp, i, LPX_UP, lo, up);
   907 	  break;
   908 	default: ;
   909 	  //FIXME error
   910 	}
   911       } else {
   912 	switch (b) {
   913 	case LPX_FR:
   914 	case LPX_LO:
   915 	  lpx_set_row_bnds(lp, i, LPX_LO, lo, up);
   916 	  break;
   917 	case LPX_UP:	  
   918 	case LPX_DB:
   919 	case LPX_FX:
   920 	  if (lo==up) 
   921 	    lpx_set_row_bnds(lp, i, LPX_FX, lo, up);
   922 	  else 
   923 	    lpx_set_row_bnds(lp, i, LPX_DB, lo, up);
   924 	  break;
   925 	default: ;
   926 	  //FIXME error
   927 	}
   928       }
   929     }
   930     virtual double _getRowLowerBound(int i) {
   931       int b=lpx_get_row_type(lp, i);
   932       switch (b) {
   933       case LPX_FR:
   934 	return -INF;
   935       case LPX_LO:
   936 	return lpx_get_row_lb(lp, i);
   937       case LPX_UP:
   938 	return -INF;
   939       case LPX_DB:
   940       case LPX_FX:
   941 	return lpx_get_row_lb(lp, i);
   942       default: ;
   943 	//FIXME error
   944 	return 0.0;
   945       }
   946     }
   947     virtual void _setRowUpperBound(int i, double up) {
   948       if (up==-INF) {
   949 	//FIXME error
   950       }
   951       int b=lpx_get_row_type(lp, i);
   952       double lo=lpx_get_row_lb(lp, i);
   953       if (up==INF) {
   954 	switch (b) {
   955 	case LPX_FR:
   956 	case LPX_LO:
   957 	  break;
   958 	case LPX_UP:
   959 	  lpx_set_row_bnds(lp, i, LPX_FR, lo, up);
   960 	  break;
   961 	case LPX_DB:
   962 	case LPX_FX:
   963 	  lpx_set_row_bnds(lp, i, LPX_LO, lo, up);
   964 	  break;
   965 	default: ;
   966 	  //FIXME error
   967 	}
   968       } else {
   969 	switch (b) {
   970 	case LPX_FR:
   971 	  lpx_set_row_bnds(lp, i, LPX_UP, lo, up);
   972 	case LPX_LO:
   973 	  if (lo==up) 
   974 	    lpx_set_row_bnds(lp, i, LPX_FX, lo, up);
   975 	  else
   976 	    lpx_set_row_bnds(lp, i, LPX_DB, lo, up);
   977 	  break;
   978 	case LPX_UP:
   979 	  lpx_set_row_bnds(lp, i, LPX_UP, lo, up);
   980 	  break;
   981 	case LPX_DB:
   982 	case LPX_FX:
   983 	  if (lo==up) 
   984 	    lpx_set_row_bnds(lp, i, LPX_FX, lo, up);
   985 	  else 
   986 	    lpx_set_row_bnds(lp, i, LPX_DB, lo, up);
   987 	  break;
   988 	default: ;
   989 	  //FIXME error
   990 	}
   991       }
   992     }
   993     virtual double _getRowUpperBound(int i) {
   994       int b=lpx_get_row_type(lp, i);
   995       switch (b) {
   996       case LPX_FR:
   997       case LPX_LO:
   998 	return INF;
   999       case LPX_UP:
  1000       case LPX_DB:
  1001       case LPX_FX:
  1002 	return lpx_get_row_ub(lp, i);
  1003       default: ;
  1004 	//FIXME error
  1005 	return 0.0;
  1006       }
  1007     }
  1008     /// \e
  1009     virtual double _getObjCoeff(int i) { 
  1010       return lpx_get_obj_coef(lp, i);
  1011     }
  1012     /// \e
  1013     virtual void _setObjCoeff(int i, double obj_coef) { 
  1014       lpx_set_obj_coef(lp, i, obj_coef);
  1015     }
  1016   public:
  1017     /// \e
  1018     void solveSimplex() { lpx_simplex(lp); }
  1019     /// \e
  1020     void solvePrimalSimplex() { lpx_simplex(lp); }
  1021     /// \e
  1022     void solveDualSimplex() { lpx_simplex(lp); }
  1023   protected:
  1024     virtual double _getPrimal(int i) {
  1025       return lpx_get_col_prim(lp, i);
  1026     }
  1027   public:
  1028     /// \e
  1029     double getObjVal() { return lpx_get_obj_val(lp); }
  1030     /// \e
  1031     int rowNum() const { return lpx_get_num_rows(lp); }
  1032     /// \e
  1033     int colNum() const { return lpx_get_num_cols(lp); }
  1034     /// \e
  1035     int warmUp() { return lpx_warm_up(lp); }
  1036     /// \e
  1037     void printWarmUpStatus(int i) {
  1038       switch (i) {
  1039       case LPX_E_OK: cout << "LPX_E_OK" << endl; break;
  1040       case LPX_E_EMPTY: cout << "LPX_E_EMPTY" << endl; break;	
  1041       case LPX_E_BADB: cout << "LPX_E_BADB" << endl; break;
  1042       case LPX_E_SING: cout << "LPX_E_SING" << endl; break;
  1043       }
  1044     }
  1045     /// \e
  1046     int getPrimalStatus() { return lpx_get_prim_stat(lp); }
  1047     /// \e
  1048     void printPrimalStatus(int i) {
  1049       switch (i) {
  1050       case LPX_P_UNDEF: cout << "LPX_P_UNDEF" << endl; break;
  1051       case LPX_P_FEAS: cout << "LPX_P_FEAS" << endl; break;	
  1052       case LPX_P_INFEAS: cout << "LPX_P_INFEAS" << endl; break;
  1053       case LPX_P_NOFEAS: cout << "LPX_P_NOFEAS" << endl; break;
  1054       }
  1055     }
  1056     /// \e
  1057     int getDualStatus() { return lpx_get_dual_stat(lp); }
  1058     /// \e
  1059     void printDualStatus(int i) {
  1060       switch (i) {
  1061       case LPX_D_UNDEF: cout << "LPX_D_UNDEF" << endl; break;
  1062       case LPX_D_FEAS: cout << "LPX_D_FEAS" << endl; break;	
  1063       case LPX_D_INFEAS: cout << "LPX_D_INFEAS" << endl; break;
  1064       case LPX_D_NOFEAS: cout << "LPX_D_NOFEAS" << endl; break;
  1065       }
  1066     }
  1067     /// Returns the status of the slack variable assigned to row \c row.
  1068     int getRowStat(const Row& row) { 
  1069       return lpx_get_row_stat(lp, row_iter_map[row]); 
  1070     }
  1071     /// \e
  1072     void printRowStatus(int i) {
  1073       switch (i) {
  1074       case LPX_BS: cout << "LPX_BS" << endl; break;
  1075       case LPX_NL: cout << "LPX_NL" << endl; break;	
  1076       case LPX_NU: cout << "LPX_NU" << endl; break;
  1077       case LPX_NF: cout << "LPX_NF" << endl; break;
  1078       case LPX_NS: cout << "LPX_NS" << endl; break;
  1079       }
  1080     }
  1081     /// Returns the status of the variable assigned to column \c col.
  1082     int getColStat(const Col& col) { 
  1083       return lpx_get_col_stat(lp, col_iter_map[col]); 
  1084     }
  1085     /// \e
  1086     void printColStatus(int i) {
  1087       switch (i) {
  1088       case LPX_BS: cout << "LPX_BS" << endl; break;
  1089       case LPX_NL: cout << "LPX_NL" << endl; break;	
  1090       case LPX_NU: cout << "LPX_NU" << endl; break;
  1091       case LPX_NF: cout << "LPX_NF" << endl; break;
  1092       case LPX_NS: cout << "LPX_NS" << endl; break;
  1093       }
  1094     }
  1095 
  1096     // MIP
  1097     /// \e
  1098     void solveBandB() { lpx_integer(lp); }
  1099     /// \e
  1100     void setLP() { lpx_set_class(lp, LPX_LP); }
  1101     /// \e
  1102     void setMIP() { lpx_set_class(lp, LPX_MIP); }
  1103   protected:
  1104     /// \e
  1105     void _setColCont(int i) { lpx_set_col_kind(lp, i, LPX_CV); }
  1106     /// \e
  1107     void _setColInt(int i) { lpx_set_col_kind(lp, i, LPX_IV); }
  1108     /// \e
  1109     double _getMIPPrimal(int i) { return lpx_mip_col_val(lp, i); }
  1110   };
  1111   
  1112   /// @}
  1113 
  1114 } //namespace lemon
  1115 
  1116 #endif //LEMON_LP_SOLVER_BASE_H