Index: src/work/athos/lp/expression.h
===================================================================
--- src/work/athos/lp/expression.h	(revision 1240)
+++ src/work/athos/lp/expression.h	(revision 1240)
@@ -0,0 +1,197 @@
+// -*- c++ -*-
+#ifndef LEMON_EXPRESSION_H
+#define LEMON_EXPRESSION_H
+
+#include <iostream>
+#include <map>
+#include <limits>
+
+namespace lemon {
+
+  /*! \brief Linear expression
+
+    \c Expr<_Col,_Value> implements a class of linear expressions with the 
+    operations of addition and multiplication with scalar. 
+
+    \author Marton Makai
+   */
+  template <typename _Col, typename _Value>
+  class Expr;
+
+  template <typename _Col, typename _Value>
+  class Expr {
+//  protected:
+  public:
+    typedef 
+    typename std::map<_Col, _Value> Data; 
+    Data data;
+  public:
+    void simplify() {
+      for (typename Data::iterator i=data.begin(); 
+	   i!=data.end(); ++i) {
+	if ((*i).second==0) data.erase(i);
+      }
+    }
+    Expr() { }
+    Expr(_Col _col) { 
+      data.insert(std::make_pair(_col, 1));
+    }
+    Expr& operator*=(_Value _value) {
+      for (typename Data::iterator i=data.begin(); 
+	   i!=data.end(); ++i) {
+	(*i).second *= _value;
+      }
+      simplify();
+      return *this;
+    }
+    Expr& operator+=(const Expr<_Col, _Value>& expr) {
+      for (typename Data::const_iterator j=expr.data.begin(); 
+	   j!=expr.data.end(); ++j) {
+	typename Data::iterator i=data.find((*j).first);
+	if (i==data.end()) {
+	  data.insert(std::make_pair((*j).first, (*j).second));
+	} else {
+	  (*i).second+=(*j).second;
+	}
+      }
+      simplify();
+      return *this;
+    }
+    Expr& operator-=(const Expr<_Col, _Value>& expr) {
+      for (typename Data::const_iterator j=expr.data.begin(); 
+	   j!=expr.data.end(); ++j) {
+	typename Data::iterator i=data.find((*j).first);
+	if (i==data.end()) {
+	  data.insert(std::make_pair((*j).first, -(*j).second));
+	} else {
+	  (*i).second+=-(*j).second;
+	}
+      }
+      simplify();
+      return *this;
+    }
+    template <typename _C, typename _V> 
+    friend std::ostream& operator<<(std::ostream& os, 
+				    const Expr<_C, _V>& expr);
+  };
+
+  template <typename _Col, typename _Value>
+  Expr<_Col, _Value> operator*(_Value _value, _Col _col) {
+    Expr<_Col, _Value> tmp(_col);
+    tmp*=_value;
+    tmp.simplify();
+    return tmp;
+  }
+
+  template <typename _Col, typename _Value>
+  Expr<_Col, _Value> operator*(_Value _value, 
+			       const Expr<_Col, _Value>& expr) {
+    Expr<_Col, _Value> tmp(expr);
+    tmp*=_value;
+    tmp.simplify();
+    return tmp;
+  }
+
+  template <typename _Col, typename _Value>
+  Expr<_Col, _Value> operator+(const Expr<_Col, _Value>& expr1, 
+			       const Expr<_Col, _Value>& expr2) {
+    Expr<_Col, _Value> tmp(expr1);
+    tmp+=expr2;
+    tmp.simplify();
+    return tmp;
+  }
+
+  template <typename _Col, typename _Value>
+  Expr<_Col, _Value> operator-(const Expr<_Col, _Value>& expr1, 
+			       const Expr<_Col, _Value>& expr2) {
+    Expr<_Col, _Value> tmp(expr1);
+    tmp-=expr2;
+    tmp.simplify();
+    return tmp;
+  }
+
+  template <typename _Col, typename _Value>
+  std::ostream& operator<<(std::ostream& os, 
+			   const Expr<_Col, _Value>& expr) {
+    for (typename Expr<_Col, _Value>::Data::const_iterator i=
+	   expr.data.begin(); 
+	 i!=expr.data.end(); ++i) {
+      os << (*i).second << "*" << (*i).first << " ";
+    }
+    return os;
+  }
+
+  template <typename _Col, typename _Value>
+  class LConstr {
+    //  protected:
+  public:
+    Expr<_Col, _Value> expr;
+    _Value lo;
+  public:
+    LConstr(const Expr<_Col, _Value>& _expr, _Value _lo) : 
+      expr(_expr), lo(_lo) { }
+  };
+  
+  template <typename _Col, typename _Value>
+  LConstr<_Col, _Value> 
+  operator<=(_Value lo, const Expr<_Col, _Value>& expr) {
+    return LConstr<_Col, _Value>(expr, lo);
+  }
+
+  template <typename _Col, typename _Value>
+  class UConstr {
+    //  protected:
+  public:
+    Expr<_Col, _Value> expr;
+    _Value up;
+  public:
+    UConstr(const Expr<_Col, _Value>& _expr, _Value _up) : 
+      expr(_expr), up(_up) { }
+  };
+
+  template <typename _Col, typename _Value>
+  UConstr<_Col, _Value> 
+  operator<=(const Expr<_Col, _Value>& expr, _Value up) {
+    return UConstr<_Col, _Value>(expr, up);
+  }
+
+  template <typename _Col, typename _Value>
+  class Constr {
+    //  protected:
+  public:
+    Expr<_Col, _Value> expr;
+    _Value lo, up;
+  public:
+    Constr(const Expr<_Col, _Value>& _expr, _Value _lo, _Value _up) : 
+      expr(_expr), lo(_lo), up(_up) { }
+    Constr(const LConstr<_Col, _Value>& _lconstr) : 
+      expr(_lconstr.expr), 
+      lo(_lconstr.lo), 
+      up(std::numeric_limits<_Value>::infinity()) { }
+    Constr(const UConstr<_Col, _Value>& _uconstr) : 
+      expr(_uconstr.expr), 
+      lo(-std::numeric_limits<_Value>::infinity()), 
+      up(_uconstr.up) { }
+  };
+
+  template <typename _Col, typename _Value>
+  Constr<_Col, _Value> 
+  operator<=(const LConstr<_Col, _Value>& lconstr, _Value up) {
+    return Constr<_Col, _Value>(lconstr.expr, lconstr.lo, up);
+  }
+
+  template <typename _Col, typename _Value>
+  Constr<_Col, _Value> 
+  operator<=(_Value lo, const UConstr<_Col, _Value>& uconstr) {
+    return Constr<_Col, _Value>(uconstr.expr, lo, uconstr.up);
+  }
+
+  template <typename _Col, typename _Value>
+  Constr<_Col, _Value> 
+  operator==(const Expr<_Col, _Value>& expr, _Value value) {
+    return Constr<_Col, _Value>(expr, value, value);
+  }
+  
+} //namespace lemon
+
+#endif //LEMON_EXPRESSION_H
Index: src/work/athos/lp/expression_test.cc
===================================================================
--- src/work/athos/lp/expression_test.cc	(revision 1240)
+++ src/work/athos/lp/expression_test.cc	(revision 1240)
@@ -0,0 +1,23 @@
+#include <expression.h>
+#include <iostream>
+#include <string>
+
+using std::cout;
+using std::endl;
+using std::string;
+using namespace lemon;
+
+int main() {
+  Expr<string, double> b;
+  cout << b << endl;
+  Expr<string, double> c("f");
+  cout << c << endl;
+  Expr<string, double> d=8.0*string("g");
+  cout << d << endl;
+  c*=5;
+  cout << c << endl;
+  Expr<string, double> e=c;
+  e+=8.9*9.0*string("l");
+  cout << e << endl;
+  cout << c+d << endl;
+}
Index: src/work/athos/lp/lp_solver_base.h
===================================================================
--- src/work/athos/lp/lp_solver_base.h	(revision 1240)
+++ src/work/athos/lp/lp_solver_base.h	(revision 1240)
@@ -0,0 +1,1116 @@
+// -*- c++ -*-
+#ifndef LEMON_LP_SOLVER_BASE_H
+#define LEMON_LP_SOLVER_BASE_H
+
+///\ingroup misc
+///\file
+
+// #include <stdio.h>
+#include <stdlib.h>
+#include <iostream>
+#include <map>
+#include <limits>
+// #include <stdio>
+//#include <stdlib>
+extern "C" {
+#include "glpk.h"
+}
+
+#include <iostream>
+#include <vector>
+#include <string>
+#include <list>
+#include <memory>
+#include <utility>
+
+#include <lemon/invalid.h>
+#include <expression.h>
+//#include <stp.h>
+//#include <lemon/max_flow.h>
+//#include <augmenting_flow.h>
+//#include <iter_map.h>
+
+using std::cout;
+using std::cin;
+using std::endl;
+
+namespace lemon {
+  
+  /// \addtogroup misc
+  /// @{
+
+  /// \brief A partitioned vector with iterable classes.
+  ///
+  /// This class implements a container in which the data is stored in an 
+  /// stl vector, the range is partitioned into sets and each set is 
+  /// doubly linked in a list. 
+  /// That is, each class is iterable by lemon iterators, and any member of 
+  /// the vector can bo moved to an other class.
+  template <typename T>
+  class IterablePartition {
+  protected:
+    struct Node {
+      T data;
+      int prev; //invalid az -1
+      int next; 
+    };
+    std::vector<Node> nodes;
+    struct Tip {
+      int first;
+      int last;
+    };
+    std::vector<Tip> tips;
+  public:
+    /// The classes are indexed by integers from \c 0 to \c classNum()-1.
+    int classNum() const { return tips.size(); }
+    /// This lemon style iterator iterates through a class. 
+    class Class;
+    /// Constructor. The number of classes is to be given which is fixed 
+    /// over the life of the container. 
+    /// The partition classes are indexed from 0 to class_num-1. 
+    IterablePartition(int class_num) { 
+      for (int i=0; i<class_num; ++i) {
+	Tip t;
+	t.first=t.last=-1;
+	tips.push_back(t);
+      }
+    }
+  protected:
+    void befuz(Class it, int class_id) {
+      if (tips[class_id].first==-1) {
+	if (tips[class_id].last==-1) {
+	  nodes[it.i].prev=nodes[it.i].next=-1;
+	  tips[class_id].first=tips[class_id].last=it.i;
+	}
+      } else {
+	nodes[it.i].prev=tips[class_id].last;
+	nodes[it.i].next=-1;
+	nodes[tips[class_id].last].next=it.i;
+	tips[class_id].last=it.i;
+      }
+    }
+    void kifuz(Class it, int class_id) {
+      if (tips[class_id].first==it.i) {
+	if (tips[class_id].last==it.i) {
+	  tips[class_id].first=tips[class_id].last=-1;
+	} else {
+	  tips[class_id].first=nodes[it.i].next;
+	  nodes[nodes[it.i].next].prev=-1;
+	}
+      } else {
+	if (tips[class_id].last==it.i) {
+	  tips[class_id].last=nodes[it.i].prev;
+	  nodes[nodes[it.i].prev].next=-1;
+	} else {
+	  nodes[nodes[it.i].next].prev=nodes[it.i].prev;
+	  nodes[nodes[it.i].prev].next=nodes[it.i].next;
+	}
+      }
+    }
+  public:
+    /// A new element with data \c t is pushed into the vector and into class 
+    /// \c class_id.
+    Class push_back(const T& t, int class_id) { 
+      Node n;
+      n.data=t;
+      nodes.push_back(n);
+      int i=nodes.size()-1;
+      befuz(i, class_id);
+      return i;
+    }
+    /// A member is moved to an other class.
+    void set(Class it, int old_class_id, int new_class_id) {
+      kifuz(it.i, old_class_id);
+      befuz(it.i, new_class_id);
+    }
+    /// Returns the data pointed by \c it.
+    T& operator[](Class it) { return nodes[it.i].data; }
+    /// Returns the data pointed by \c it.
+    const T& operator[](Class it) const { return nodes[it.i].data; }
+    ///.
+    class Class {
+      friend class IterablePartition;
+    protected:
+      int i;
+    public:
+      /// Default constructor.
+      Class() { }
+      /// This constructor constructs an iterator which points
+      /// to the member of th container indexed by the integer _i.
+      Class(const int& _i) : i(_i) { }
+      /// Invalid constructor.
+      Class(const Invalid&) : i(-1) { }
+      friend bool operator<(const Class& x, const Class& y);
+      friend std::ostream& operator<<(std::ostream& os, 
+				      const Class& it);
+      bool operator==(const Class& node) const {return i == node.i;}
+      bool operator!=(const Class& node) const {return i != node.i;}
+    };
+    friend bool operator<(const Class& x, const Class& y) {
+      return (x.i < y.i);
+    }
+    friend std::ostream& operator<<(std::ostream& os, 
+				    const Class& it) {
+      os << it.i;
+      return os;
+    }
+    /// First member of class \c class_id.
+    Class& first(Class& it, int class_id) const {
+      it.i=tips[class_id].first;
+      return it;
+    }
+    /// Next member.
+    Class& next(Class& it) const {
+      it.i=nodes[it.i].next;
+      return it;
+    }
+    /// True iff the iterator is valid.
+    bool valid(const Class& it) const { return it.i!=-1; }
+
+    class ClassIt : public Class {
+      const IterablePartition* iterable_partition;
+    public:
+      ClassIt() { }
+      ClassIt(Invalid i) : Class(i) { }
+      ClassIt(const IterablePartition& _iterable_partition, 
+	      const int& i) : iterable_partition(&_iterable_partition) {
+        _iterable_partition.first(*this, i);
+      }
+      ClassIt(const IterablePartition& _iterable_partition, 
+	      const Class& _class) : 
+	Class(_class), iterable_partition(&_iterable_partition) { }
+      ClassIt& operator++() {
+        iterable_partition->next(*this);
+        return *this;
+      }
+    };
+
+  };
+
+
+  /*! \e
+    \todo kellenene uj iterable structure bele, mert ez nem az igazi
+    \todo A[x,y]-t cserel. Jobboldal, baloldal csere.
+    \todo LEKERDEZESEK!!!
+    \todo DOKSI!!!! Doxygen group!!!
+    The aim of this class is to give a general surface to different 
+    solvers, i.e. it makes possible to write algorithms using LP's, 
+    in which the solver can be changed to an other one easily.
+    \nosubgrouping
+  */
+  template <typename _Value>
+  class LPSolverBase {
+    
+    /*! @name Uncategorized functions and types (public members)
+    */
+    //@{
+  public:
+
+    //UNCATEGORIZED
+
+    /// \e
+    typedef IterablePartition<int> Rows;
+    /// \e
+    typedef IterablePartition<int> Cols;
+    /// \e
+    typedef _Value Value;
+    /// \e
+    typedef Rows::Class Row;
+    /// \e
+    typedef Cols::Class Col;
+  public:
+    /// \e
+    IterablePartition<int> row_iter_map;
+    /// \e
+    IterablePartition<int> col_iter_map;
+    /// \e
+    std::vector<Row> int_row_map;
+    /// \e
+    std::vector<Col> int_col_map;
+    /// \e
+    const int VALID_CLASS;
+    /// \e
+    const int INVALID_CLASS;
+    /// \e 
+    static const _Value INF;
+  public:
+    /// \e
+    LPSolverBase() : row_iter_map(2), 
+		     col_iter_map(2), 
+		     VALID_CLASS(0), INVALID_CLASS(1) { }
+    /// \e
+    virtual ~LPSolverBase() { }
+    //@}
+
+    /*! @name Medium level interface (public members)
+      These functions appear in the low level and also in the high level 
+      interfaces thus these each of these functions have to be implemented 
+      only once in the different interfaces.
+      This means that these functions have to be reimplemented for all of the 
+      different lp solvers. These are basic functions, and have the same 
+      parameter lists in the low and high level interfaces. 
+    */
+    //@{
+  public:
+
+    //UNCATEGORIZED FUNCTIONS
+
+    /// \e
+    virtual void setMinimize() = 0;
+    /// \e
+    virtual void setMaximize() = 0;
+
+    //SOLVER FUNCTIONS
+
+    /// \e
+    virtual void solveSimplex() = 0;
+    /// \e
+    virtual void solvePrimalSimplex() = 0;
+    /// \e
+    virtual void solveDualSimplex() = 0;
+
+    //SOLUTION RETRIEVING
+
+    /// \e
+    virtual _Value getObjVal() = 0;
+
+    //OTHER FUNCTIONS
+
+    /// \e
+    virtual int rowNum() const = 0;
+    /// \e
+    virtual int colNum() const = 0;
+    /// \e
+    virtual int warmUp() = 0;
+    /// \e
+    virtual void printWarmUpStatus(int i) = 0;
+    /// \e
+    virtual int getPrimalStatus() = 0;
+    /// \e
+    virtual void printPrimalStatus(int i) = 0;
+    /// \e
+    virtual int getDualStatus() = 0;
+    /// \e
+    virtual void printDualStatus(int i) = 0;
+    /// Returns the status of the slack variable assigned to row \c row.
+    virtual int getRowStat(const Row& row) = 0;
+    /// \e
+    virtual void printRowStatus(int i) = 0;
+    /// Returns the status of the variable assigned to column \c col.
+    virtual int getColStat(const Col& col) = 0;
+    /// \e
+    virtual void printColStatus(int i) = 0;
+
+    //@}
+
+    /*! @name Low level interface (protected members)
+      Problem manipulating functions in the low level interface
+    */
+    //@{
+  protected:
+
+    //MATRIX MANIPULATING FUNCTIONS
+
+    /// \e
+    virtual int _addCol() = 0;
+    /// \e
+    virtual int _addRow() = 0;
+    /// \e
+    virtual void _eraseCol(int i) = 0;
+    /// \e
+    virtual void _eraseRow(int i) = 0;
+    /// \e
+    virtual void _setRowCoeffs(int i, 
+			       const std::vector<std::pair<int, _Value> >& coeffs) = 0;
+    /// \e
+    /// This routine modifies \c coeffs only by the \c push_back method.
+    virtual void _getRowCoeffs(int i, 
+			       std::vector<std::pair<int, _Value> >& coeffs) = 0;
+    /// \e
+    virtual void _setColCoeffs(int i, 
+			       const std::vector<std::pair<int, _Value> >& coeffs) = 0;
+    /// \e
+    /// This routine modifies \c coeffs only by the \c push_back method.
+    virtual void _getColCoeffs(int i, 
+			       std::vector<std::pair<int, _Value> >& coeffs) = 0;
+    /// \e
+    virtual void _setCoeff(int col, int row, _Value value) = 0;
+    /// \e
+    virtual _Value _getCoeff(int col, int row) = 0;
+    //  public:
+    //    /// \e
+    //    enum Bound { FREE, LOWER, UPPER, DOUBLE, FIXED };
+  protected:
+    /// \e
+    /// The lower bound of a variable (column) have to be given by an 
+    /// extended number of type _Value, i.e. a finite number of type 
+    /// _Value or -INF.
+    virtual void _setColLowerBound(int i, _Value value) = 0;
+    /// \e
+    /// The lower bound of a variable (column) is an 
+    /// extended number of type _Value, i.e. a finite number of type 
+    /// _Value or -INF.
+    virtual _Value _getColLowerBound(int i) = 0;
+    /// \e
+    /// The upper bound of a variable (column) have to be given by an 
+    /// extended number of type _Value, i.e. a finite number of type 
+    /// _Value or INF.
+    virtual void _setColUpperBound(int i, _Value value) = 0;
+    /// \e
+    /// The upper bound of a variable (column) is an 
+    /// extended number of type _Value, i.e. a finite number of type 
+    /// _Value or INF.
+    virtual _Value _getColUpperBound(int i) = 0;
+    /// \e
+    /// The lower bound of a linear expression (row) have to be given by an 
+    /// extended number of type _Value, i.e. a finite number of type 
+    /// _Value or -INF.
+    virtual void _setRowLowerBound(int i, _Value value) = 0;
+    /// \e
+    /// The lower bound of a linear expression (row) is an 
+    /// extended number of type _Value, i.e. a finite number of type 
+    /// _Value or -INF.
+    virtual _Value _getRowLowerBound(int i) = 0;
+    /// \e
+    /// The upper bound of a linear expression (row) have to be given by an 
+    /// extended number of type _Value, i.e. a finite number of type 
+    /// _Value or INF.
+    virtual void _setRowUpperBound(int i, _Value value) = 0;
+    /// \e
+    /// The upper bound of a linear expression (row) is an 
+    /// extended number of type _Value, i.e. a finite number of type 
+    /// _Value or INF.
+    virtual _Value _getRowUpperBound(int i) = 0;
+    /// \e
+    virtual void _setObjCoeff(int i, _Value obj_coef) = 0;
+    /// \e
+    virtual _Value _getObjCoeff(int i) = 0;
+    
+    //SOLUTION RETRIEVING
+
+    /// \e
+    virtual _Value _getPrimal(int i) = 0;
+    //@}
+    
+    /*! @name High level interface (public members)
+      Problem manipulating functions in the high level interface
+    */
+    //@{
+  public:
+
+    //MATRIX MANIPULATING FUNCTIONS
+
+    /// \e
+    Col addCol() {
+      int i=_addCol();  
+      Col col;
+      col_iter_map.first(col, INVALID_CLASS);
+      if (col_iter_map.valid(col)) { //van hasznalhato hely
+	col_iter_map.set(col, INVALID_CLASS, VALID_CLASS);
+	col_iter_map[col]=i;
+      } else { //a cucc vegere kell inzertalni mert nincs szabad hely
+	col=col_iter_map.push_back(i, VALID_CLASS);
+      }
+      int_col_map.push_back(col);
+      return col;
+    }
+    /// \e
+    Row addRow() {
+      int i=_addRow();
+      Row row;
+      row_iter_map.first(row, INVALID_CLASS);
+      if (row_iter_map.valid(row)) { //van hasznalhato hely
+	row_iter_map.set(row, INVALID_CLASS, VALID_CLASS);
+	row_iter_map[row]=i;
+      } else { //a cucc vegere kell inzertalni mert nincs szabad hely
+	row=row_iter_map.push_back(i, VALID_CLASS);
+      }
+      int_row_map.push_back(row);
+      return row;
+    }
+    /// \e
+    void eraseCol(const Col& col) {
+      col_iter_map.set(col, VALID_CLASS, INVALID_CLASS);
+      int cols[2];
+      cols[1]=col_iter_map[col];
+      _eraseCol(cols[1]);
+      col_iter_map[col]=0; //glpk specifikus, de kell ez??
+      Col it;
+      for (col_iter_map.first(it, VALID_CLASS); 
+	   col_iter_map.valid(it); col_iter_map.next(it)) {
+	if (col_iter_map[it]>cols[1]) --col_iter_map[it];
+      }
+      int_col_map.erase(int_col_map.begin()+cols[1]);
+    }
+    /// \e
+    void eraseRow(const Row& row) {
+      row_iter_map.set(row, VALID_CLASS, INVALID_CLASS);
+      int rows[2];
+      rows[1]=row_iter_map[row];
+      _eraseRow(rows[1]);
+      row_iter_map[row]=0; //glpk specifikus, de kell ez??
+      Row it;
+      for (row_iter_map.first(it, VALID_CLASS); 
+	   row_iter_map.valid(it); row_iter_map.next(it)) {
+	if (row_iter_map[it]>rows[1]) --row_iter_map[it];
+      }
+      int_row_map.erase(int_row_map.begin()+rows[1]);
+    }
+    /// \e
+    void setCoeff(Col col, Row row, _Value value) {
+      _setCoeff(col_iter_map[col], row_iter_map[row], value);
+    }
+    /// \e
+    _Value getCoeff(Col col, Row row) {
+      return _getCoeff(col_iter_map[col], row_iter_map[row], value);
+    }
+    /// \e
+    void setColLowerBound(Col col, _Value lo) {
+      _setColLowerBound(col_iter_map[col], lo);
+    }
+    /// \e
+    _Value getColLowerBound(Col col) {
+      return _getColLowerBound(col_iter_map[col]);
+    }
+    /// \e
+    void setColUpperBound(Col col, _Value up) {
+      _setColUpperBound(col_iter_map[col], up);
+    }
+    /// \e
+    _Value getColUpperBound(Col col) {      
+      return _getColUpperBound(col_iter_map[col]);
+    }
+    /// \e
+    void setRowLowerBound(Row row, _Value lo) {
+      _setRowLowerBound(row_iter_map[row], lo);
+    }
+    /// \e
+    _Value getRowLowerBound(Row row) {
+      return _getRowLowerBound(row_iter_map[row]);
+    }
+    /// \e
+    void setRowUpperBound(Row row, _Value up) {
+      _setRowUpperBound(row_iter_map[row], up);
+    }
+    /// \e
+    _Value getRowUpperBound(Row row) {      
+      return _getRowUpperBound(row_iter_map[row]);
+    }
+    /// \e
+    void setObjCoeff(const Col& col, _Value obj_coef) {
+      _setObjCoeff(col_iter_map[col], obj_coef);
+    }
+    /// \e
+    _Value getObjCoeff(const Col& col) {
+      return _getObjCoeff(col_iter_map[col]);
+    }
+
+    //SOLUTION RETRIEVING FUNCTIONS
+
+    /// \e
+    _Value getPrimal(const Col& col) {
+      return _getPrimal(col_iter_map[col]);
+    }    
+
+    //@}
+
+    /*! @name User friend interface
+      Problem manipulating functions in the user friend interface
+    */
+    //@{
+
+    //EXPRESSION TYPES
+
+    /// \e
+    typedef Expr<Col, _Value> Expression;
+    /// \e
+    typedef Expr<Row, _Value> DualExpression;
+    /// \e
+    typedef Constr<Col, _Value> Constraint;
+
+    //MATRIX MANIPULATING FUNCTIONS
+
+    /// \e
+    void setRowCoeffs(Row row, const Expression& expr) {
+      std::vector<std::pair<int, _Value> > row_coeffs;
+      for(typename Expression::Data::const_iterator i=expr.data.begin(); 
+	  i!=expr.data.end(); ++i) {
+	row_coeffs.push_back(std::make_pair
+			     (col_iter_map[(*i).first], (*i).second));
+      }
+      _setRowCoeffs(row_iter_map[row], row_coeffs);
+    }
+    /// \e 
+    void setRow(Row row, const Constraint& constr) {
+      setRowCoeffs(row, constr.expr);
+      setRowLowerBound(row, constr.lo);
+      setRowUpperBound(row, constr.up);
+    }
+    /// \e 
+    Row addRow(const Constraint& constr) {
+      Row row=addRow();
+      setRowCoeffs(row, constr.expr);
+      setRowLowerBound(row, constr.lo);
+      setRowUpperBound(row, constr.up);
+      return row;
+    }
+    /// \e
+    /// This routine modifies \c expr by only adding to it.
+    void getRowCoeffs(Row row, Expression& expr) {
+      std::vector<std::pair<int, _Value> > row_coeffs;
+      _getRowCoeffs(row_iter_map[row], row_coeffs);
+      for(typename std::vector<std::pair<int, _Value> >::const_iterator 
+ 	    i=row_coeffs.begin(); i!=row_coeffs.end(); ++i) {
+ 	expr+= (*i).second*int_col_map[(*i).first];
+      }
+    }
+    /// \e
+    void setColCoeffs(Col col, const DualExpression& expr) {
+      std::vector<std::pair<int, _Value> > col_coeffs;
+      for(typename DualExpression::Data::const_iterator i=expr.data.begin(); 
+	  i!=expr.data.end(); ++i) {
+	col_coeffs.push_back(std::make_pair
+			     (row_iter_map[(*i).first], (*i).second));
+      }
+      _setColCoeffs(col_iter_map[col], col_coeffs);
+    }
+    /// \e
+    /// This routine modifies \c expr by only adding to it.
+    void getColCoeffs(Col col, DualExpression& expr) {
+      std::vector<std::pair<int, _Value> > col_coeffs;
+      _getColCoeffs(col_iter_map[col], col_coeffs);
+      for(typename std::vector<std::pair<int, _Value> >::const_iterator 
+ 	    i=col_coeffs.begin(); i!=col_coeffs.end(); ++i) {
+ 	expr+= (*i).second*int_row_map[(*i).first];
+      }
+    }
+    /// \e
+    void setObjCoeffs(const Expression& expr) {
+      // writing zero everywhere
+      for(Cols::ClassIt it(col_iter_map, VALID_CLASS); it!=INVALID; ++it)
+	setObjCoeff(it, 0.0);
+      // writing the data needed
+      for(typename Expression::Data::const_iterator i=expr.data.begin(); 
+	  i!=expr.data.end(); ++i) {
+	setObjCoeff((*i).first, (*i).second);
+      }
+    }
+    /// \e
+    /// This routine modifies \c expr by only adding to it.
+    void getObjCoeffs(Expression& expr) {
+      for(Cols::ClassIt it(col_iter_map, VALID_CLASS); it!=INVALID; ++it)
+	expr+=getObjCoeff(it)*it;
+    }
+    //@}
+
+
+    /*! @name MIP functions and types (public members)
+    */
+    //@{
+  public:
+    /// \e
+    virtual void solveBandB() = 0;
+    /// \e
+    virtual void setLP() = 0;
+    /// \e
+    virtual void setMIP() = 0;
+  protected:
+   /// \e
+    virtual void _setColCont(int i) = 0;
+    /// \e
+    virtual void _setColInt(int i) = 0;
+    /// \e
+    virtual _Value _getMIPPrimal(int i) = 0;
+  public:
+    /// \e
+    void setColCont(Col col) {
+      _setColCont(col_iter_map[col]);
+    }
+    /// \e
+    void setColInt(Col col) {
+      _setColInt(col_iter_map[col]);
+    }
+    /// \e
+    _Value getMIPPrimal(Col col) {
+      return _getMIPPrimal(col_iter_map[col]);
+    }
+    //@}
+  };
+  
+  template <typename _Value>
+  const _Value LPSolverBase<_Value>::INF=std::numeric_limits<_Value>::infinity();
+
+
+  /// \brief Wrapper for GLPK solver
+  /// 
+  /// This class implements a lemon wrapper for GLPK.
+  class LPGLPK : public LPSolverBase<double> {
+  public:
+    typedef LPSolverBase<double> Parent;
+
+  public:
+    /// \e
+    LPX* lp;
+
+  public:
+    /// \e
+    LPGLPK() : Parent(), 
+			lp(lpx_create_prob()) {
+      int_row_map.push_back(Row());
+      int_col_map.push_back(Col());
+      lpx_set_int_parm(lp, LPX_K_DUAL, 1);
+    }
+    /// \e
+    ~LPGLPK() {
+      lpx_delete_prob(lp);
+    }
+
+    //MATRIX INDEPEDENT MANIPULATING FUNCTIONS
+
+    /// \e
+    void setMinimize() { 
+      lpx_set_obj_dir(lp, LPX_MIN);
+    }
+    /// \e
+    void setMaximize() { 
+      lpx_set_obj_dir(lp, LPX_MAX);
+    }
+
+    //LOW LEVEL INTERFACE, MATRIX MANIPULATING FUNCTIONS
+
+  protected:
+    /// \e
+    int _addCol() { 
+      int i=lpx_add_cols(lp, 1);
+      _setColLowerBound(i, -INF);
+      _setColUpperBound(i, INF);
+      return i;
+    }
+    /// \e
+    int _addRow() { 
+      int i=lpx_add_rows(lp, 1);
+      return i;
+    }
+    /// \e
+    virtual void _setRowCoeffs(int i, 
+			       const std::vector<std::pair<int, double> >& coeffs) {
+      int mem_length=1+colNum();
+      int* indices = new int[mem_length];
+      double* doubles = new double[mem_length];
+      int length=0;
+      for (std::vector<std::pair<int, double> >::
+	     const_iterator it=coeffs.begin(); it!=coeffs.end(); ++it) {
+	++length;
+	indices[length]=it->first;
+	doubles[length]=it->second;
+      }
+      lpx_set_mat_row(lp, i, length, indices, doubles);
+      delete [] indices;
+      delete [] doubles;
+    }
+    /// \e
+    virtual void _getRowCoeffs(int i, 
+			       std::vector<std::pair<int, double> >& coeffs) {
+      int mem_length=1+colNum();
+      int* indices = new int[mem_length];
+      double* doubles = new double[mem_length];
+      int length=lpx_get_mat_row(lp, i, indices, doubles);
+      for (int i=1; i<=length; ++i) {
+	coeffs.push_back(std::make_pair(indices[i], doubles[i]));
+      }
+      delete [] indices;
+      delete [] doubles;
+    }
+    /// \e
+    virtual void _setColCoeffs(int i, 
+			       const std::vector<std::pair<int, double> >& coeffs) {
+      int mem_length=1+rowNum();
+      int* indices = new int[mem_length];
+      double* doubles = new double[mem_length];
+      int length=0;
+      for (std::vector<std::pair<int, double> >::
+	     const_iterator it=coeffs.begin(); it!=coeffs.end(); ++it) {
+	++length;
+	indices[length]=it->first;
+	doubles[length]=it->second;
+      }
+      lpx_set_mat_col(lp, i, length, indices, doubles);
+      delete [] indices;
+      delete [] doubles;
+    }
+    /// \e
+    virtual void _getColCoeffs(int i, 
+			       std::vector<std::pair<int, double> >& coeffs) {
+      int mem_length=1+rowNum();
+      int* indices = new int[mem_length];
+      double* doubles = new double[mem_length];
+      int length=lpx_get_mat_col(lp, i, indices, doubles);
+      for (int i=1; i<=length; ++i) {
+	coeffs.push_back(std::make_pair(indices[i], doubles[i]));
+      }
+      delete [] indices;
+      delete [] doubles;
+    }
+    /// \e
+    virtual void _eraseCol(int i) {
+      int cols[2];
+      cols[1]=i;
+      lpx_del_cols(lp, 1, cols);
+    }
+    virtual void _eraseRow(int i) {
+      int rows[2];
+      rows[1]=i;
+      lpx_del_rows(lp, 1, rows);
+    }
+    void _setCoeff(int col, int row, double value) {
+      /// FIXME not yet implemented
+    }
+    double _getCoeff(int col, int row) {
+      /// FIXME not yet implemented
+      return 0.0;
+    }
+    virtual void _setColLowerBound(int i, double lo) {
+      if (lo==INF) {
+	//FIXME error
+      }
+      int b=lpx_get_col_type(lp, i);
+      double up=lpx_get_col_ub(lp, i);	
+      if (lo==-INF) {
+	switch (b) {
+	case LPX_FR:
+	case LPX_LO:
+	  lpx_set_col_bnds(lp, i, LPX_FR, lo, up);
+	  break;
+	case LPX_UP:
+	  break;
+	case LPX_DB:
+	case LPX_FX:
+	  lpx_set_col_bnds(lp, i, LPX_UP, lo, up);
+	  break;
+	default: ;
+	  //FIXME error
+	}
+      } else {
+	switch (b) {
+	case LPX_FR:
+	case LPX_LO:
+	  lpx_set_col_bnds(lp, i, LPX_LO, lo, up);
+	  break;
+	case LPX_UP:	  
+	case LPX_DB:
+	case LPX_FX:
+	  if (lo==up) 
+	    lpx_set_col_bnds(lp, i, LPX_FX, lo, up);
+	  else 
+	    lpx_set_col_bnds(lp, i, LPX_DB, lo, up);
+	  break;
+	default: ;
+	  //FIXME error
+	}
+      }
+    }
+    virtual double _getColLowerBound(int i) {
+      int b=lpx_get_col_type(lp, i);
+      switch (b) {
+      case LPX_FR:
+	return -INF;
+      case LPX_LO:
+	return lpx_get_col_lb(lp, i);
+      case LPX_UP:
+	return -INF;
+      case LPX_DB:
+      case LPX_FX:
+	return lpx_get_col_lb(lp, i);
+      default: ;
+	//FIXME error
+	return 0.0;
+      }
+    }
+    virtual void _setColUpperBound(int i, double up) {
+      if (up==-INF) {
+	//FIXME error
+      }
+      int b=lpx_get_col_type(lp, i);
+      double lo=lpx_get_col_lb(lp, i);
+      if (up==INF) {
+	switch (b) {
+	case LPX_FR:
+	case LPX_LO:
+	  break;
+	case LPX_UP:
+	  lpx_set_col_bnds(lp, i, LPX_FR, lo, up);
+	  break;
+	case LPX_DB:
+	case LPX_FX:
+	  lpx_set_col_bnds(lp, i, LPX_LO, lo, up);
+	  break;
+	default: ;
+	  //FIXME error
+	}
+      } else {
+	switch (b) {
+	case LPX_FR:
+	  lpx_set_col_bnds(lp, i, LPX_UP, lo, up);
+	case LPX_LO:
+	  if (lo==up) 
+	    lpx_set_col_bnds(lp, i, LPX_FX, lo, up);
+	  else
+	    lpx_set_col_bnds(lp, i, LPX_DB, lo, up);
+	  break;
+	case LPX_UP:
+	  lpx_set_col_bnds(lp, i, LPX_UP, lo, up);
+	  break;
+	case LPX_DB:
+	case LPX_FX:
+	  if (lo==up) 
+	    lpx_set_col_bnds(lp, i, LPX_FX, lo, up);
+	  else 
+	    lpx_set_col_bnds(lp, i, LPX_DB, lo, up);
+	  break;
+	default: ;
+	  //FIXME error
+	}
+      }
+    }
+    virtual double _getColUpperBound(int i) {
+      int b=lpx_get_col_type(lp, i);
+      switch (b) {
+      case LPX_FR:
+      case LPX_LO:
+	return INF;
+      case LPX_UP:
+      case LPX_DB:
+      case LPX_FX:
+	return lpx_get_col_ub(lp, i);
+      default: ;
+	//FIXME error
+	return 0.0;
+      }
+    }
+    virtual void _setRowLowerBound(int i, double lo) {
+      if (lo==INF) {
+	//FIXME error
+      }
+      int b=lpx_get_row_type(lp, i);
+      double up=lpx_get_row_ub(lp, i);	
+      if (lo==-INF) {
+	switch (b) {
+	case LPX_FR:
+	case LPX_LO:
+	  lpx_set_row_bnds(lp, i, LPX_FR, lo, up);
+	  break;
+	case LPX_UP:
+	  break;
+	case LPX_DB:
+	case LPX_FX:
+	  lpx_set_row_bnds(lp, i, LPX_UP, lo, up);
+	  break;
+	default: ;
+	  //FIXME error
+	}
+      } else {
+	switch (b) {
+	case LPX_FR:
+	case LPX_LO:
+	  lpx_set_row_bnds(lp, i, LPX_LO, lo, up);
+	  break;
+	case LPX_UP:	  
+	case LPX_DB:
+	case LPX_FX:
+	  if (lo==up) 
+	    lpx_set_row_bnds(lp, i, LPX_FX, lo, up);
+	  else 
+	    lpx_set_row_bnds(lp, i, LPX_DB, lo, up);
+	  break;
+	default: ;
+	  //FIXME error
+	}
+      }
+    }
+    virtual double _getRowLowerBound(int i) {
+      int b=lpx_get_row_type(lp, i);
+      switch (b) {
+      case LPX_FR:
+	return -INF;
+      case LPX_LO:
+	return lpx_get_row_lb(lp, i);
+      case LPX_UP:
+	return -INF;
+      case LPX_DB:
+      case LPX_FX:
+	return lpx_get_row_lb(lp, i);
+      default: ;
+	//FIXME error
+	return 0.0;
+      }
+    }
+    virtual void _setRowUpperBound(int i, double up) {
+      if (up==-INF) {
+	//FIXME error
+      }
+      int b=lpx_get_row_type(lp, i);
+      double lo=lpx_get_row_lb(lp, i);
+      if (up==INF) {
+	switch (b) {
+	case LPX_FR:
+	case LPX_LO:
+	  break;
+	case LPX_UP:
+	  lpx_set_row_bnds(lp, i, LPX_FR, lo, up);
+	  break;
+	case LPX_DB:
+	case LPX_FX:
+	  lpx_set_row_bnds(lp, i, LPX_LO, lo, up);
+	  break;
+	default: ;
+	  //FIXME error
+	}
+      } else {
+	switch (b) {
+	case LPX_FR:
+	  lpx_set_row_bnds(lp, i, LPX_UP, lo, up);
+	case LPX_LO:
+	  if (lo==up) 
+	    lpx_set_row_bnds(lp, i, LPX_FX, lo, up);
+	  else
+	    lpx_set_row_bnds(lp, i, LPX_DB, lo, up);
+	  break;
+	case LPX_UP:
+	  lpx_set_row_bnds(lp, i, LPX_UP, lo, up);
+	  break;
+	case LPX_DB:
+	case LPX_FX:
+	  if (lo==up) 
+	    lpx_set_row_bnds(lp, i, LPX_FX, lo, up);
+	  else 
+	    lpx_set_row_bnds(lp, i, LPX_DB, lo, up);
+	  break;
+	default: ;
+	  //FIXME error
+	}
+      }
+    }
+    virtual double _getRowUpperBound(int i) {
+      int b=lpx_get_row_type(lp, i);
+      switch (b) {
+      case LPX_FR:
+      case LPX_LO:
+	return INF;
+      case LPX_UP:
+      case LPX_DB:
+      case LPX_FX:
+	return lpx_get_row_ub(lp, i);
+      default: ;
+	//FIXME error
+	return 0.0;
+      }
+    }
+    /// \e
+    virtual double _getObjCoeff(int i) { 
+      return lpx_get_obj_coef(lp, i);
+    }
+    /// \e
+    virtual void _setObjCoeff(int i, double obj_coef) { 
+      lpx_set_obj_coef(lp, i, obj_coef);
+    }
+  public:
+    /// \e
+    void solveSimplex() { lpx_simplex(lp); }
+    /// \e
+    void solvePrimalSimplex() { lpx_simplex(lp); }
+    /// \e
+    void solveDualSimplex() { lpx_simplex(lp); }
+  protected:
+    virtual double _getPrimal(int i) {
+      return lpx_get_col_prim(lp, i);
+    }
+  public:
+    /// \e
+    double getObjVal() { return lpx_get_obj_val(lp); }
+    /// \e
+    int rowNum() const { return lpx_get_num_rows(lp); }
+    /// \e
+    int colNum() const { return lpx_get_num_cols(lp); }
+    /// \e
+    int warmUp() { return lpx_warm_up(lp); }
+    /// \e
+    void printWarmUpStatus(int i) {
+      switch (i) {
+      case LPX_E_OK: cout << "LPX_E_OK" << endl; break;
+      case LPX_E_EMPTY: cout << "LPX_E_EMPTY" << endl; break;	
+      case LPX_E_BADB: cout << "LPX_E_BADB" << endl; break;
+      case LPX_E_SING: cout << "LPX_E_SING" << endl; break;
+      }
+    }
+    /// \e
+    int getPrimalStatus() { return lpx_get_prim_stat(lp); }
+    /// \e
+    void printPrimalStatus(int i) {
+      switch (i) {
+      case LPX_P_UNDEF: cout << "LPX_P_UNDEF" << endl; break;
+      case LPX_P_FEAS: cout << "LPX_P_FEAS" << endl; break;	
+      case LPX_P_INFEAS: cout << "LPX_P_INFEAS" << endl; break;
+      case LPX_P_NOFEAS: cout << "LPX_P_NOFEAS" << endl; break;
+      }
+    }
+    /// \e
+    int getDualStatus() { return lpx_get_dual_stat(lp); }
+    /// \e
+    void printDualStatus(int i) {
+      switch (i) {
+      case LPX_D_UNDEF: cout << "LPX_D_UNDEF" << endl; break;
+      case LPX_D_FEAS: cout << "LPX_D_FEAS" << endl; break;	
+      case LPX_D_INFEAS: cout << "LPX_D_INFEAS" << endl; break;
+      case LPX_D_NOFEAS: cout << "LPX_D_NOFEAS" << endl; break;
+      }
+    }
+    /// Returns the status of the slack variable assigned to row \c row.
+    int getRowStat(const Row& row) { 
+      return lpx_get_row_stat(lp, row_iter_map[row]); 
+    }
+    /// \e
+    void printRowStatus(int i) {
+      switch (i) {
+      case LPX_BS: cout << "LPX_BS" << endl; break;
+      case LPX_NL: cout << "LPX_NL" << endl; break;	
+      case LPX_NU: cout << "LPX_NU" << endl; break;
+      case LPX_NF: cout << "LPX_NF" << endl; break;
+      case LPX_NS: cout << "LPX_NS" << endl; break;
+      }
+    }
+    /// Returns the status of the variable assigned to column \c col.
+    int getColStat(const Col& col) { 
+      return lpx_get_col_stat(lp, col_iter_map[col]); 
+    }
+    /// \e
+    void printColStatus(int i) {
+      switch (i) {
+      case LPX_BS: cout << "LPX_BS" << endl; break;
+      case LPX_NL: cout << "LPX_NL" << endl; break;	
+      case LPX_NU: cout << "LPX_NU" << endl; break;
+      case LPX_NF: cout << "LPX_NF" << endl; break;
+      case LPX_NS: cout << "LPX_NS" << endl; break;
+      }
+    }
+
+    // MIP
+    /// \e
+    void solveBandB() { lpx_integer(lp); }
+    /// \e
+    void setLP() { lpx_set_class(lp, LPX_LP); }
+    /// \e
+    void setMIP() { lpx_set_class(lp, LPX_MIP); }
+  protected:
+    /// \e
+    void _setColCont(int i) { lpx_set_col_kind(lp, i, LPX_CV); }
+    /// \e
+    void _setColInt(int i) { lpx_set_col_kind(lp, i, LPX_IV); }
+    /// \e
+    double _getMIPPrimal(int i) { return lpx_mip_col_val(lp, i); }
+  };
+  
+  /// @}
+
+} //namespace lemon
+
+#endif //LEMON_LP_SOLVER_BASE_H
Index: src/work/athos/lp/lp_solver_wrapper.h
===================================================================
--- src/work/athos/lp/lp_solver_wrapper.h	(revision 1240)
+++ src/work/athos/lp/lp_solver_wrapper.h	(revision 1240)
@@ -0,0 +1,431 @@
+// -*- c++ -*-
+#ifndef LEMON_LP_SOLVER_WRAPPER_H
+#define LEMON_LP_SOLVER_WRAPPER_H
+
+///\ingroup misc
+///\file
+///\brief Dijkstra algorithm.
+
+// #include <stdio.h>
+#include <stdlib.h>
+// #include <stdio>
+//#include <stdlib>
+extern "C" {
+#include "glpk.h"
+}
+
+#include <iostream>
+#include <vector>
+#include <string>
+#include <list>
+#include <memory>
+#include <utility>
+
+//#include <sage_graph.h>
+//#include <lemon/list_graph.h>
+//#include <lemon/graph_wrapper.h>
+#include <lemon/invalid.h>
+//#include <bfs_dfs.h>
+//#include <stp.h>
+//#include <lemon/max_flow.h>
+//#include <augmenting_flow.h>
+//#include <iter_map.h>
+
+using std::cout;
+using std::cin;
+using std::endl;
+
+namespace lemon {
+
+  
+  /// \addtogroup misc
+  /// @{
+
+  /// \brief A partitioned vector with iterable classes.
+  ///
+  /// This class implements a container in which the data is stored in an 
+  /// stl vector, the range is partitioned into sets and each set is 
+  /// doubly linked in a list. 
+  /// That is, each class is iterable by lemon iterators, and any member of 
+  /// the vector can bo moved to an other class.
+  template <typename T>
+  class IterablePartition {
+  protected:
+    struct Node {
+      T data;
+      int prev; //invalid az -1
+      int next; 
+    };
+    std::vector<Node> nodes;
+    struct Tip {
+      int first;
+      int last;
+    };
+    std::vector<Tip> tips;
+  public:
+    /// The classes are indexed by integers from \c 0 to \c classNum()-1.
+    int classNum() const { return tips.size(); }
+    /// This lemon style iterator iterates through a class. 
+    class ClassIt;
+    /// Constructor. The number of classes is to be given which is fixed 
+    /// over the life of the container. 
+    /// The partition classes are indexed from 0 to class_num-1. 
+    IterablePartition(int class_num) { 
+      for (int i=0; i<class_num; ++i) {
+	Tip t;
+	t.first=t.last=-1;
+	tips.push_back(t);
+      }
+    }
+  protected:
+    void befuz(ClassIt it, int class_id) {
+      if (tips[class_id].first==-1) {
+	if (tips[class_id].last==-1) {
+	  nodes[it.i].prev=nodes[it.i].next=-1;
+	  tips[class_id].first=tips[class_id].last=it.i;
+	}
+      } else {
+	nodes[it.i].prev=tips[class_id].last;
+	nodes[it.i].next=-1;
+	nodes[tips[class_id].last].next=it.i;
+	tips[class_id].last=it.i;
+      }
+    }
+    void kifuz(ClassIt it, int class_id) {
+      if (tips[class_id].first==it.i) {
+	if (tips[class_id].last==it.i) {
+	  tips[class_id].first=tips[class_id].last=-1;
+	} else {
+	  tips[class_id].first=nodes[it.i].next;
+	  nodes[nodes[it.i].next].prev=-1;
+	}
+      } else {
+	if (tips[class_id].last==it.i) {
+	  tips[class_id].last=nodes[it.i].prev;
+	  nodes[nodes[it.i].prev].next=-1;
+	} else {
+	  nodes[nodes[it.i].next].prev=nodes[it.i].prev;
+	  nodes[nodes[it.i].prev].next=nodes[it.i].next;
+	}
+      }
+    }
+  public:
+    /// A new element with data \c t is pushed into the vector and into class 
+    /// \c class_id.
+    ClassIt push_back(const T& t, int class_id) { 
+      Node n;
+      n.data=t;
+      nodes.push_back(n);
+      int i=nodes.size()-1;
+      befuz(i, class_id);
+      return i;
+    }
+    /// A member is moved to an other class.
+    void set(ClassIt it, int old_class_id, int new_class_id) {
+      kifuz(it.i, old_class_id);
+      befuz(it.i, new_class_id);
+    }
+    /// Returns the data pointed by \c it.
+    T& operator[](ClassIt it) { return nodes[it.i].data; }
+    /// Returns the data pointed by \c it.
+    const T& operator[](ClassIt it) const { return nodes[it.i].data; }
+    ///.
+    class ClassIt {
+      friend class IterablePartition;
+    protected:
+      int i;
+    public:
+      /// Default constructor.
+      ClassIt() { }
+      /// This constructor constructs an iterator which points
+      /// to the member of th container indexed by the integer _i.
+      ClassIt(const int& _i) : i(_i) { }
+      /// Invalid constructor.
+      ClassIt(const Invalid&) : i(-1) { }
+    };
+    /// First member of class \c class_id.
+    ClassIt& first(ClassIt& it, int class_id) const {
+      it.i=tips[class_id].first;
+      return it;
+    }
+    /// Next member.
+    ClassIt& next(ClassIt& it) const {
+      it.i=nodes[it.i].next;
+      return it;
+    }
+    /// True iff the iterator is valid.
+    bool valid(const ClassIt& it) const { return it.i!=-1; }
+  };
+  
+  /// \brief Wrappers for LP solvers
+  /// 
+  /// This class implements a lemon wrapper for glpk.
+  /// Later other LP-solvers will be wrapped into lemon.
+  /// The aim of this class is to give a general surface to different 
+  /// solvers, i.e. it makes possible to write algorithms using LP's, 
+  /// in which the solver can be changed to an other one easily.
+  class LPSolverWrapper {
+  public:
+
+//   class Row {
+//   protected:
+//     int i;
+//   public:
+//     Row() { }
+//     Row(const Invalid&) : i(0) { }
+//     Row(const int& _i) : i(_i) { }
+//     operator int() const { return i; }
+//   };
+//   class RowIt : public Row {
+//   public:
+//     RowIt(const Row& row) : Row(row) { }
+//   };
+
+//   class Col {
+//   protected:
+//     int i;
+//   public:
+//     Col() { }
+//     Col(const Invalid&) : i(0) { }
+//     Col(const int& _i) : i(_i) { }
+//     operator int() const { return i; }
+//   };
+//   class ColIt : public Col {
+//     ColIt(const Col& col) : Col(col) { }
+//   };
+
+  public:
+    ///.
+    LPX* lp;
+    ///.
+    typedef IterablePartition<int>::ClassIt RowIt;
+    ///.
+    IterablePartition<int> row_iter_map;
+    ///.
+    typedef IterablePartition<int>::ClassIt ColIt;
+    ///.
+    IterablePartition<int> col_iter_map;
+    //std::vector<int> row_id_to_lp_row_id;
+    //std::vector<int> col_id_to_lp_col_id;
+    ///.
+    const int VALID_ID;
+    ///.
+    const int INVALID_ID;
+
+  public:
+    ///.
+    LPSolverWrapper() : lp(lpx_create_prob()), 
+			row_iter_map(2), 
+			col_iter_map(2), 
+			//row_id_to_lp_row_id(), col_id_to_lp_col_id(), 
+			VALID_ID(0), INVALID_ID(1) {
+      lpx_set_int_parm(lp, LPX_K_DUAL, 1);
+    }
+    ///.
+    ~LPSolverWrapper() {
+      lpx_delete_prob(lp);
+    }
+    ///.
+    void setMinimize() { 
+      lpx_set_obj_dir(lp, LPX_MIN);
+    }
+    ///.
+    void setMaximize() { 
+      lpx_set_obj_dir(lp, LPX_MAX);
+    }
+    ///.
+    ColIt addCol() {
+      int i=lpx_add_cols(lp, 1);  
+      ColIt col_it;
+      col_iter_map.first(col_it, INVALID_ID);
+      if (col_iter_map.valid(col_it)) { //van hasznalhato hely
+	col_iter_map.set(col_it, INVALID_ID, VALID_ID);
+	col_iter_map[col_it]=i;
+	//col_id_to_lp_col_id[col_iter_map[col_it]]=i;
+      } else { //a cucc vegere kell inzertalni mert nincs szabad hely
+	//col_id_to_lp_col_id.push_back(i);
+	//int j=col_id_to_lp_col_id.size()-1;
+	col_it=col_iter_map.push_back(i, VALID_ID);
+      }
+//    edge_index_map.set(e, i);
+//    lpx_set_col_bnds(lp, i, LPX_DB, 0.0, 1.0);
+//    lpx_set_obj_coef(lp, i, cost[e]);    
+      return col_it;
+    }
+    ///.
+    RowIt addRow() {
+      int i=lpx_add_rows(lp, 1);  
+      RowIt row_it;
+      row_iter_map.first(row_it, INVALID_ID);
+      if (row_iter_map.valid(row_it)) { //van hasznalhato hely
+	row_iter_map.set(row_it, INVALID_ID, VALID_ID);
+	row_iter_map[row_it]=i;
+      } else { //a cucc vegere kell inzertalni mert nincs szabad hely
+	row_it=row_iter_map.push_back(i, VALID_ID);
+      }
+      return row_it;
+    }
+    //pair<RowIt, double>-bol kell megadni egy std range-et
+    ///.
+    template <typename Begin, typename End>
+    void setColCoeffs(const ColIt& col_it, 
+		      Begin begin, End end) {
+      int mem_length=1+lpx_get_num_rows(lp);
+      int* indices = new int[mem_length];
+      double* doubles = new double[mem_length];
+      int length=0;
+      for ( ; begin!=end; ++begin) {
+	++length;
+	indices[length]=row_iter_map[begin->first];
+	doubles[length]=begin->second;
+      }
+      lpx_set_mat_col(lp, col_iter_map[col_it], length, indices, doubles);
+      delete [] indices;
+      delete [] doubles;
+    }
+    //pair<ColIt, double>-bol kell megadni egy std range-et
+    ///.
+    template <typename Begin, typename End>
+    void setRowCoeffs(const RowIt& row_it, 
+		      Begin begin, End end) {
+      int mem_length=1+lpx_get_num_cols(lp);
+      int* indices = new int[mem_length];
+      double* doubles = new double[mem_length];
+      int length=0;
+      for ( ; begin!=end; ++begin) {
+	++length;
+	indices[length]=col_iter_map[begin->first];
+	doubles[length]=begin->second;
+      }
+      lpx_set_mat_row(lp, row_iter_map[row_it], length, indices, doubles);
+      delete [] indices;
+      delete [] doubles;
+    }
+    ///.
+    void eraseCol(const ColIt& col_it) {
+      col_iter_map.set(col_it, VALID_ID, INVALID_ID);
+      int cols[2];
+      cols[1]=col_iter_map[col_it];
+      lpx_del_cols(lp, 1, cols);
+      col_iter_map[col_it]=0; //glpk specifikus
+      ColIt it;
+      for (col_iter_map.first(it, VALID_ID); 
+	   col_iter_map.valid(it); col_iter_map.next(it)) {
+	if (col_iter_map[it]>cols[1]) --col_iter_map[it];
+      }
+    }
+    ///.
+    void eraseRow(const RowIt& row_it) {
+      row_iter_map.set(row_it, VALID_ID, INVALID_ID);
+      int rows[2];
+      rows[1]=row_iter_map[row_it];
+      lpx_del_rows(lp, 1, rows);
+      row_iter_map[row_it]=0; //glpk specifikus
+      RowIt it;
+      for (row_iter_map.first(it, VALID_ID); 
+	   row_iter_map.valid(it); row_iter_map.next(it)) {
+	if (row_iter_map[it]>rows[1]) --row_iter_map[it];
+      }
+    }
+    ///.
+    void setColBounds(const ColIt& col_it, int bound_type, 
+		      double lo, double up) {
+      lpx_set_col_bnds(lp, col_iter_map[col_it], bound_type, lo, up);
+    }
+    ///.
+    double getObjCoef(const ColIt& col_it) { 
+      return lpx_get_obj_coef(lp, col_iter_map[col_it]);
+    }
+    ///.
+    void setRowBounds(const RowIt& row_it, int bound_type, 
+		      double lo, double up) {
+      lpx_set_row_bnds(lp, row_iter_map[row_it], bound_type, lo, up);
+    }
+    ///.
+    void setObjCoef(const ColIt& col_it, double obj_coef) { 
+      lpx_set_obj_coef(lp, col_iter_map[col_it], obj_coef);
+    }
+    ///.
+    void solveSimplex() { lpx_simplex(lp); }
+    ///.
+    void solvePrimalSimplex() { lpx_simplex(lp); }
+    ///.
+    void solveDualSimplex() { lpx_simplex(lp); }
+    ///.
+    double getPrimal(const ColIt& col_it) {
+      return lpx_get_col_prim(lp, col_iter_map[col_it]);
+    }
+    ///.
+    double getObjVal() { return lpx_get_obj_val(lp); }
+    ///.
+    int rowNum() const { return lpx_get_num_rows(lp); }
+    ///.
+    int colNum() const { return lpx_get_num_cols(lp); }
+    ///.
+    int warmUp() { return lpx_warm_up(lp); }
+    ///.
+    void printWarmUpStatus(int i) {
+      switch (i) {
+	case LPX_E_OK: cout << "LPX_E_OK" << endl; break;
+	case LPX_E_EMPTY: cout << "LPX_E_EMPTY" << endl; break;	
+	case LPX_E_BADB: cout << "LPX_E_BADB" << endl; break;
+	case LPX_E_SING: cout << "LPX_E_SING" << endl; break;
+      }
+    }
+    ///.
+    int getPrimalStatus() { return lpx_get_prim_stat(lp); }
+    ///.
+    void printPrimalStatus(int i) {
+      switch (i) {
+	case LPX_P_UNDEF: cout << "LPX_P_UNDEF" << endl; break;
+	case LPX_P_FEAS: cout << "LPX_P_FEAS" << endl; break;	
+	case LPX_P_INFEAS: cout << "LPX_P_INFEAS" << endl; break;
+	case LPX_P_NOFEAS: cout << "LPX_P_NOFEAS" << endl; break;
+      }
+    }
+    ///.
+    int getDualStatus() { return lpx_get_dual_stat(lp); }
+    ///.
+    void printDualStatus(int i) {
+      switch (i) {
+	case LPX_D_UNDEF: cout << "LPX_D_UNDEF" << endl; break;
+	case LPX_D_FEAS: cout << "LPX_D_FEAS" << endl; break;	
+	case LPX_D_INFEAS: cout << "LPX_D_INFEAS" << endl; break;
+	case LPX_D_NOFEAS: cout << "LPX_D_NOFEAS" << endl; break;
+      }
+    }
+    /// Returns the status of the slack variable assigned to row \c row_it.
+    int getRowStat(const RowIt& row_it) { 
+      return lpx_get_row_stat(lp, row_iter_map[row_it]); 
+    }
+    ///.
+    void printRowStatus(int i) {
+      switch (i) {
+	case LPX_BS: cout << "LPX_BS" << endl; break;
+	case LPX_NL: cout << "LPX_NL" << endl; break;	
+	case LPX_NU: cout << "LPX_NU" << endl; break;
+	case LPX_NF: cout << "LPX_NF" << endl; break;
+	case LPX_NS: cout << "LPX_NS" << endl; break;
+      }
+    }
+    /// Returns the status of the variable assigned to column \c col_it.
+    int getColStat(const ColIt& col_it) { 
+      return lpx_get_col_stat(lp, col_iter_map[col_it]); 
+    }
+    ///.
+    void printColStatus(int i) {
+      switch (i) {
+	case LPX_BS: cout << "LPX_BS" << endl; break;
+	case LPX_NL: cout << "LPX_NL" << endl; break;	
+	case LPX_NU: cout << "LPX_NU" << endl; break;
+	case LPX_NF: cout << "LPX_NF" << endl; break;
+	case LPX_NS: cout << "LPX_NS" << endl; break;
+      }
+    }
+  };
+  
+  /// @}
+
+} //namespace lemon
+
+#endif //LEMON_LP_SOLVER_WRAPPER_H
Index: src/work/athos/lp/magic_square.cc
===================================================================
--- src/work/athos/lp/magic_square.cc	(revision 1240)
+++ src/work/athos/lp/magic_square.cc	(revision 1240)
@@ -0,0 +1,90 @@
+// -*- c++ -*-
+#include <iostream>
+#include <fstream>
+
+#include <lemon/time_measure.h>
+#include <lp_solver_base.h>
+
+using std::cout;
+using std::endl;
+using namespace lemon;
+
+/*
+  On an 1537Mhz PC, the run times with 
+  glpk are the following.
+  for n=3,4, some secondes
+  for n=5, 25 hours
+ */
+
+int main(int, char **) {
+  const int n=4;
+  const double row_sum=(1.0+n*n)*n/2;
+  Timer ts;
+  ts.reset();
+  typedef LPGLPK LPSolver;
+  typedef LPSolver::Col Col;
+  LPSolver lp;
+  typedef std::map<std::pair<int, int>, Col> Coords;
+  Coords x;
+  // we create a new variable for each entry 
+  // of the magic square
+  for (int i=1; i<=n; ++i) {
+    for (int j=1; j<=n; ++j) { 
+      Col col=lp.addCol();
+      x[std::make_pair(i,j)]=col;
+      lp.setColLowerBound(col, 1.0);
+      lp.setColUpperBound(col, double(n*n));
+    }
+  }
+  LPSolver::Expression expr3, expr4;
+  for (int i=1; i<=n; ++i) {
+    LPSolver::Expression expr1, expr2;
+    for (int j=1; j<=n; ++j) {
+      expr1+=x[std::make_pair(i, j)];
+      expr2+=x[std::make_pair(j, i)];
+    }
+    // sum of rows and columns
+    lp.addRow(expr1==row_sum);
+    lp.addRow(expr2==row_sum);
+    expr3+=x[std::make_pair(i, i)];
+    expr4+=x[std::make_pair(i, (n+1)-i)];
+  }
+  // sum of the diagonal entries
+  lp.addRow(expr3==row_sum);
+  lp.addRow(expr4==row_sum);
+  lp.solveSimplex();
+  cout << "elapsed time: " << ts << endl;
+  for (int i=1; i<=n; ++i) {
+    for (int j=1; j<=n; ++j) { 
+      cout << "x("<<i<<","<<j<<")="<<lp.getPrimal(x[std::make_pair(i,j)]) 
+	   << endl;
+    }
+  }
+  // we make new binary variables for each pair of 
+  // entries of the square to achieve that 
+  // the values of different entries are different
+  lp.setMIP();
+  for (Coords::const_iterator it=x.begin(); it!=x.end(); ++it) {
+    Coords::const_iterator jt=it; ++jt;
+    for(; jt!=x.end(); ++jt) {
+      Col col1=(*it).second;
+      Col col2=(*jt).second;
+      Col col=lp.addCol();
+      lp.setColLowerBound(col, 0.0);
+      lp.setColUpperBound(col, 1.0);
+      lp.addRow(double(-n*n+1.0)<=1.0*col2-1.0*col1-double(n*n)*col<=-1.0);
+      lp.setColInt(col);
+    }
+  }
+  cout << "elapsed time: " << ts << endl;
+  lp.solveSimplex();
+  // let's solve the integer problem
+  lp.solveBandB();
+  cout << "elapsed time: " << ts << endl;
+  for (int i=1; i<=n; ++i) {
+    for (int j=1; j<=n; ++j) { 
+      cout << "x("<<i<<","<<j<<")="<<lp.getMIPPrimal(x[std::make_pair(i,j)]) 
+	   << endl;
+    }
+  }
+}
Index: src/work/athos/lp/makefile
===================================================================
--- src/work/athos/lp/makefile	(revision 1240)
+++ src/work/athos/lp/makefile	(revision 1240)
@@ -0,0 +1,73 @@
+#INCLUDEDIRS ?= -I.. -I. -I./{marci,jacint,alpar,klao,akos}
+#GLPKROOT = /usr/local/glpk-4.4
+INCLUDEDIRS ?= -I/home/marci/boost -I/usr/local/cplex/cplex75/include -I../../{marci,alpar,klao,akos,athos} -I. -I../../.. -I../.. -I..# -I$(GLPKROOT)/include
+#INCLUDEDIRS ?= -I../.. -I../.. -I../../{marci,jacint,alpar,klao,akos} -I/usr/local/glpk-4.4/include
+CXXFLAGS = -g -O2 -W -Wall $(INCLUDEDIRS) -ansi -pedantic
+LDFLAGS  =  -lglpk#-lcplex -lm -lpthread -lilocplex -L/usr/local/cplex/cplex75/lib/i86_linux2_glibc2.2_gcc3.0/static_mt# -L$(GLPKROOT)/lib
+
+BINARIES = magic_square max_flow_expression expression_test max_flow_by_lp# sample sample2 sample11 sample15
+
+#include ../makefile
+
+# Hat, ez elismerem, hogy nagyon ronda, de mukodik minden altalam
+# ismert rendszeren :-)  (Misi)
+ifdef GCCVER
+CXX := g++-$(GCCVER)
+else
+CXX := $(shell type -p g++-3.3 || type -p g++-3.2 || type -p g++-3.0 || type -p g++-3 || echo g++)
+endif
+
+ifdef DEBUG
+CXXFLAGS += -DDEBUG
+endif
+
+CC := $(CXX)
+
+
+all: $(BINARIES)
+
+################
+# Minden binarishoz egy sor, felsorolva, hogy mely object file-okbol
+# all elo.
+# Kiveve ha siman file.cc -> file  esetrol van szo, amikor is nem kell
+# irni semmit.
+
+#proba: proba.o seged.o
+
+################
+
+
+# .depend dep depend:
+# 	-$(CXX) $(CXXFLAGS) -M $(BINARIES:=.cc) > .depend
+
+#makefile: .depend
+#sinclude .depend
+#moert nem megy az eredeti /usr/bin/ld-vel?
+
+# %: %.o
+# 	$(CXX) -o $@ $< $(LDFLAGS)
+
+# %.o: %.cc
+# 	$(CXX) $(CXXFLAGS) -c $<
+
+%: %.cc
+	$(CXX) $(CXXFLAGS) -o $@ $< $(LDFLAGS)
+
+sample11prof: sample11prof.o
+	 $(CXX) -pg -o sample11prof sample11prof.o -L$(GLPKROOT)/lib -lglpk
+sample11prof.o: sample11.cc
+	$(CXX) -pg $(CXXFLAGS) -c -o sample11prof.o sample11.cc
+
+# sample.o: sample.cc
+# 	$(CXX) $(CXXFLAGS) -c -o sample.o sample.cc
+
+# sample2: sample2.o
+# 	$(CXX) -o sample2 sample2.o -L/usr/local/glpk-4.4/lib -lglpk
+# sample2.o: sample2.cc
+# 	$(CXX) $(CXXFLAGS) -c -o sample2.o sample2.cc
+
+
+clean:
+	$(RM) *.o $(BINARIES) .depend
+
+.PHONY: all clean dep depend
Index: src/work/athos/lp/max_flow_by_lp.cc
===================================================================
--- src/work/athos/lp/max_flow_by_lp.cc	(revision 1240)
+++ src/work/athos/lp/max_flow_by_lp.cc	(revision 1240)
@@ -0,0 +1,186 @@
+// -*- c++ -*-
+#include <iostream>
+#include <fstream>
+
+#include <lemon/smart_graph.h>
+#include <lemon/list_graph.h>
+#include <lemon/dimacs.h>
+#include <lemon/time_measure.h>
+//#include <graph_wrapper.h>
+#include <lemon/preflow.h>
+#include <augmenting_flow.h>
+//#include <preflow_res.h>
+//#include <lp_solver_wrapper_2.h>
+#include <min_cost_gen_flow.h>
+
+// Use a DIMACS max flow file as stdin.
+// max_flow_demo < dimacs_max_flow_file
+
+using namespace lemon;
+
+int main(int, char **) {
+
+  typedef ListGraph MutableGraph;
+  typedef ListGraph Graph;
+  typedef Graph::Node Node;
+  typedef Graph::Edge Edge;
+  typedef Graph::EdgeIt EdgeIt;
+
+  Graph g;
+
+  Node s, t;
+  Graph::EdgeMap<int> cap(g);
+  //readDimacsMaxFlow(std::cin, g, s, t, cap);
+  readDimacs(std::cin, g, cap, s, t);
+  Timer ts;
+  Graph::EdgeMap<int> flow(g); //0 flow
+  Preflow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> > 
+    max_flow_test(g, s, t, cap, flow);
+  AugmentingFlow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> > 
+    augmenting_flow_test(g, s, t, cap, flow);
+  
+  Graph::NodeMap<bool> cut(g);
+
+  {
+    std::cout << "preflow ..." << std::endl;
+    ts.reset();
+    max_flow_test.run();
+    std::cout << "elapsed time: " << ts << std::endl;
+    std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
+    max_flow_test.minCut(cut);
+
+    for (EdgeIt e(g); e!=INVALID; ++e) {
+      if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) 
+	std::cout << "Slackness does not hold!" << std::endl;
+      if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0) 
+	std::cout << "Slackness does not hold!" << std::endl;
+    }
+  }
+
+//   {
+//     std::cout << "preflow ..." << std::endl;
+//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
+//     ts.reset();
+//     max_flow_test.preflow(Preflow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >::GEN_FLOW);
+//     std::cout << "elapsed time: " << ts << std::endl;
+//     std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
+
+//     FOR_EACH_LOC(Graph::EdgeIt, e, g) {
+//       if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) 
+// 	std::cout << "Slackness does not hold!" << std::endl;
+//       if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0) 
+// 	std::cout << "Slackness does not hold!" << std::endl;
+//     }
+//   }
+
+//   {
+//     std::cout << "wrapped preflow ..." << std::endl;
+//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
+//     ts.reset();
+//     pre_flow_res.run();
+//     std::cout << "elapsed time: " << ts << std::endl;
+//     std::cout << "flow value: "<< pre_flow_test.flowValue() << std::endl;
+//   }
+
+  {
+    std::cout << "physical blocking flow augmentation ..." << std::endl;
+    for (EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
+    ts.reset();
+    int i=0;
+    while (augmenting_flow_test.augmentOnBlockingFlow<MutableGraph>()) { ++i; }
+    std::cout << "elapsed time: " << ts << std::endl;
+    std::cout << "number of augmentation phases: " << i << std::endl; 
+    std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
+
+    for (EdgeIt e(g); e!=INVALID; ++e) {
+      if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) 
+	std::cout << "Slackness does not hold!" << std::endl;
+      if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0) 
+	std::cout << "Slackness does not hold!" << std::endl;
+    }
+  }
+
+//   {
+//     std::cout << "faster physical blocking flow augmentation ..." << std::endl;
+//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
+//     ts.reset();
+//     int i=0;
+//     while (max_flow_test.augmentOnBlockingFlow1<MutableGraph>()) { ++i; }
+//     std::cout << "elapsed time: " << ts << std::endl;
+//     std::cout << "number of augmentation phases: " << i << std::endl; 
+//     std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
+//   }
+
+  {
+    std::cout << "on-the-fly blocking flow augmentation ..." << std::endl;
+    for (EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
+    ts.reset();
+    int i=0;
+    while (augmenting_flow_test.augmentOnBlockingFlow2()) { ++i; }
+    std::cout << "elapsed time: " << ts << std::endl;
+    std::cout << "number of augmentation phases: " << i << std::endl; 
+    std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
+
+    for (EdgeIt e(g); e!=INVALID; ++e) {
+      if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) 
+	std::cout << "Slackness does not hold!" << std::endl;
+      if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0) 
+	std::cout << "Slackness does not hold!" << std::endl;
+    }
+  }
+
+//   {
+//     std::cout << "on-the-fly shortest path augmentation ..." << std::endl;
+//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
+//     ts.reset();
+//     int i=0;
+//     while (augmenting_flow_test.augmentOnShortestPath()) { ++i; }
+//     std::cout << "elapsed time: " << ts << std::endl;
+//     std::cout << "number of augmentation phases: " << i << std::endl; 
+//     std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
+
+//     FOR_EACH_LOC(Graph::EdgeIt, e, g) {
+//       if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) 
+// 	std::cout << "Slackness does not hold!" << std::endl;
+//       if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0) 
+// 	std::cout << "Slackness does not hold!" << std::endl;
+//     }
+//   }
+
+//   {
+//     std::cout << "on-the-fly shortest path augmentation ..." << std::endl;
+//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
+//     ts.reset();
+//     int i=0;
+//     while (augmenting_flow_test.augmentOnShortestPath2()) { ++i; }
+//     std::cout << "elapsed time: " << ts << std::endl;
+//     std::cout << "number of augmentation phases: " << i << std::endl; 
+//     std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
+
+//     FOR_EACH_LOC(Graph::EdgeIt, e, g) {
+//       if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) 
+// 	std::cout << "Slackness does not hold!" << std::endl;
+//       if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0) 
+// 	std::cout << "Slackness does not hold!" << std::endl;
+//     }
+//   }
+
+  ts.reset();
+
+  Edge e=g.addEdge(t, s);
+  Graph::EdgeMap<int> cost(g, 0);
+  cost.set(e, -1);
+  cap.set(e, 10000);
+  typedef ConstMap<Node, int> Excess;
+  Excess excess(0);
+  typedef ConstMap<Edge, int> LCap;
+  LCap lcap(0);
+
+  MinCostGenFlow<Graph, int, Excess, LCap> 
+    min_cost(g, excess, lcap, cap, flow, cost); 
+  min_cost.feasible();
+  min_cost.runByLP();
+
+  std::cout << "elapsed time: " << ts << std::endl;
+  std::cout << "flow value: "<< flow[e] << std::endl;
+}
Index: src/work/athos/lp/max_flow_expression.cc
===================================================================
--- src/work/athos/lp/max_flow_expression.cc	(revision 1240)
+++ src/work/athos/lp/max_flow_expression.cc	(revision 1240)
@@ -0,0 +1,109 @@
+// -*- c++ -*-
+#include <iostream>
+#include <fstream>
+
+#include <lemon/graph_utils.h>
+#include <lemon/smart_graph.h>
+#include <lemon/list_graph.h>
+#include <lemon/dimacs.h>
+#include <lemon/time_measure.h>
+#include <lp_solver_base.h>
+
+using std::cout;
+using std::endl;
+using namespace lemon;
+
+template<typename Edge, typename EdgeIndexMap> 
+class PrimalMap {
+protected:
+  LPGLPK* lp;
+  EdgeIndexMap* edge_index_map;
+public:
+  PrimalMap(LPGLPK& _lp, EdgeIndexMap& _edge_index_map) : 
+    lp(&_lp), edge_index_map(&_edge_index_map) { }
+  double operator[](Edge e) const { 
+    return lp->getPrimal((*edge_index_map)[e]);
+  }
+};
+
+// Use a DIMACS max flow file as stdin.
+// max_flow_expresion < dimacs_max_flow_file
+
+int main(int, char **) {
+
+  typedef ListGraph Graph;
+  typedef Graph::Node Node;
+  typedef Graph::Edge Edge;
+  typedef Graph::EdgeIt EdgeIt;
+
+  Graph g;
+  
+  Node s, t;
+  Graph::EdgeMap<int> cap(g);
+  readDimacs(std::cin, g, cap, s, t);
+  Timer ts;
+  
+  typedef LPGLPK LPSolver;
+  LPSolver lp;
+  lp.setMaximize();
+  typedef LPSolver::Col Col;
+  typedef LPSolver::Row Row;
+  typedef Graph::EdgeMap<Col> EdgeIndexMap;
+  typedef Graph::NodeMap<Row> NodeIndexMap;
+  EdgeIndexMap edge_index_map(g);
+  NodeIndexMap node_index_map(g);
+  PrimalMap<Edge, EdgeIndexMap> flow(lp, edge_index_map);
+
+  // nonnegativity of flow and capacity function
+  for (Graph::EdgeIt e(g); e!=INVALID; ++e) {
+    Col col=lp.addCol();
+    edge_index_map.set(e, col);
+    // interesting property in GLPK:
+    // if you change the order of the following two lines, the 
+    // two runs of GLPK are extremely different
+      lp.setColLowerBound(col, 0);
+      lp.setColUpperBound(col, cap[e]);
+  }
+  
+  for (Graph::NodeIt n(g); n!=INVALID; ++n) {
+    LPSolver::Expression expr;
+    for (Graph::OutEdgeIt e(g, n); e!=INVALID; ++e)
+      expr+=edge_index_map[e];
+    for (Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
+      expr-=edge_index_map[e];
+    // cost function
+    if (n==s) {
+      lp.setObjCoeffs(expr);      
+    }
+    // flow conservation constraints
+    if ((n!=s) && (n!=t)) {
+      node_index_map.set(n, lp.addRow(expr == 0.0));
+    }
+  }
+  lp.solveSimplex();
+  cout << "elapsed time: " << ts << endl;
+//   cout << "rows:" << endl;
+//   for (
+//        LPSolver::Rows::ClassIt i(lp.row_iter_map, 0);
+//        i!=INVALID;
+//        ++i) { 
+//     cout << i << " ";
+//   }
+//   cout << endl;
+//   cout << "cols:" << endl;
+//   for (
+//        LPSolver::Cols::ClassIt i(lp.col_iter_map, 0);
+//        i!=INVALID;
+//        ++i) { 
+//     cout << i << " ";
+//   }
+//   cout << endl;
+  lp.setMIP();
+  cout << "elapsed time: " << ts << endl;
+  for (LPSolver::Cols::ClassIt it(lp.col_iter_map ,1); it!=INVALID; ++it) {
+    lp.setColInt(it);
+  }
+  cout << "elapsed time: " << ts << endl;
+  lp.solveBandB();
+  cout << "elapsed time: " << ts << endl;
+}
Index: src/work/athos/lp/min_cost_gen_flow.h
===================================================================
--- src/work/athos/lp/min_cost_gen_flow.h	(revision 1240)
+++ src/work/athos/lp/min_cost_gen_flow.h	(revision 1240)
@@ -0,0 +1,268 @@
+// -*- c++ -*-
+#ifndef LEMON_MIN_COST_GEN_FLOW_H
+#define LEMON_MIN_COST_GEN_FLOW_H
+#include <iostream>
+//#include <fstream>
+
+#include <lemon/smart_graph.h>
+#include <lemon/list_graph.h>
+//#include <lemon/dimacs.h>
+//#include <lemon/time_measure.h>
+//#include <graph_wrapper.h>
+#include <lemon/preflow.h>
+#include <lemon/min_cost_flow.h>
+//#include <augmenting_flow.h>
+//#include <preflow_res.h>
+#include <work/marci/merge_node_graph_wrapper.h>
+#include <work/marci/lp/lp_solver_wrapper_3.h>
+
+namespace lemon {
+
+  template<typename Edge, typename EdgeIndexMap> 
+  class PrimalMap {
+  protected:
+    LPGLPK* lp;
+    EdgeIndexMap* edge_index_map;
+  public:
+    PrimalMap(LPGLPK& _lp, EdgeIndexMap& _edge_index_map) : 
+      lp(&_lp), edge_index_map(&_edge_index_map) { }
+    double operator[](Edge e) const { 
+      return lp->getPrimal((*edge_index_map)[e]);
+    }
+  };
+
+  // excess: rho-delta egyelore csak =0-ra.
+  template <typename Graph, typename Num,
+	    typename Excess=typename Graph::template NodeMap<Num>, 
+	    typename LCapMap=typename Graph::template EdgeMap<Num>,
+	    typename CapMap=typename Graph::template EdgeMap<Num>,
+            typename FlowMap=typename Graph::template EdgeMap<Num>,
+	    typename CostMap=typename Graph::template EdgeMap<Num> >
+  class MinCostGenFlow {
+  protected:
+    const Graph& g;
+    const Excess& excess;
+    const LCapMap& lcapacity;
+    const CapMap& capacity;
+    FlowMap& flow;
+    const CostMap& cost;
+  public:
+    MinCostGenFlow(const Graph& _g, const Excess& _excess, 
+		   const LCapMap& _lcapacity, const CapMap& _capacity, 
+		   FlowMap& _flow, 
+		   const CostMap& _cost) :
+      g(_g), excess(_excess), lcapacity(_lcapacity),
+      capacity(_capacity), flow(_flow), cost(_cost) { }
+    bool feasible() {
+      //      std::cout << "making new vertices..." << std::endl; 
+      typedef ListGraph Graph2;
+      Graph2 g2;
+      typedef MergeEdgeGraphWrapper<const Graph, Graph2> GW;
+      //      std::cout << "merging..." << std::endl; 
+      GW gw(g, g2);
+      typename GW::Node s(INVALID, g2.addNode(), true);
+      typename GW::Node t(INVALID, g2.addNode(), true);
+      typedef SmartGraph Graph3;
+      //      std::cout << "making extender graph..." << std::endl; 
+      typedef NewEdgeSetGraphWrapper2<GW, Graph3> GWW;
+//       {
+// 	checkConcept<StaticGraph, GWW>();   
+//       }
+      GWW gww(gw);
+      typedef AugmentingGraphWrapper<GW, GWW> GWWW;
+      GWWW gwww(gw, gww);
+
+      //      std::cout << "making new edges..." << std::endl; 
+      typename GWWW::template EdgeMap<Num> translated_cap(gwww);
+
+      for (typename GW::EdgeIt e(gw); e!=INVALID; ++e) {
+	translated_cap.set(typename GWWW::Edge(e,INVALID,false), 
+			   capacity[e]-lcapacity[e]);
+	//	cout << "t_cap " << gw.id(e) << " " 
+	//	     << translated_cap[typename GWWW::Edge(e,INVALID,false)] << endl;
+      }
+
+      Num expected=0;
+
+      //      std::cout << "making new edges 2..." << std::endl; 
+      for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
+	Num a=0;
+	for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
+	  a+=lcapacity[e];
+	for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) 
+	  a-=lcapacity[e];
+	if (excess[n]>a) {
+	  typename GWW::Edge e=
+	    gww.addEdge(typename GW::Node(n,INVALID,false), t);
+	  translated_cap.set(typename GWWW::Edge(INVALID, e, true), 
+			     excess[n]-a);
+	  //	  std::cout << g.id(n) << "->t " << excess[n]-a << std::endl;
+	}
+	if (excess[n]<a) {
+	  typename GWW::Edge e=
+	    gww.addEdge(s, typename GW::Node(n,INVALID,false));
+	  translated_cap.set(typename GWWW::Edge(INVALID, e, true), 
+			     a-excess[n]);
+	  expected+=a-excess[n];
+	  //	  std::cout << "s->" << g.id(n) << " "<< a-excess[n] <<std:: endl;
+	}
+      }
+
+      //      std::cout << "preflow..." << std::endl; 
+      typename GWWW::template EdgeMap<Num> translated_flow(gwww, 0);
+      Preflow<GWWW, Num> preflow(gwww, s, t, 
+				 translated_cap, translated_flow);
+      preflow.run();
+      //      std::cout << "fv: " << preflow.flowValue() << std::endl; 
+      //      std::cout << "expected: " << expected << std::endl; 
+
+      for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
+	typename GW::Edge ew(e, INVALID, false);
+	typename GWWW::Edge ewww(ew, INVALID, false);
+	flow.set(e, translated_flow[ewww]+lcapacity[e]);
+      }
+      return (preflow.flowValue()>=expected);
+    }
+    // for nonnegative costs
+    bool run() {
+      //      std::cout << "making new vertices..." << std::endl; 
+      typedef ListGraph Graph2;
+      Graph2 g2;
+      typedef MergeEdgeGraphWrapper<const Graph, Graph2> GW;
+      //      std::cout << "merging..." << std::endl; 
+      GW gw(g, g2);
+      typename GW::Node s(INVALID, g2.addNode(), true);
+      typename GW::Node t(INVALID, g2.addNode(), true);
+      typedef SmartGraph Graph3;
+      //      std::cout << "making extender graph..." << std::endl; 
+      typedef NewEdgeSetGraphWrapper2<GW, Graph3> GWW;
+//       {
+// 	checkConcept<StaticGraph, GWW>();   
+//       }
+      GWW gww(gw);
+      typedef AugmentingGraphWrapper<GW, GWW> GWWW;
+      GWWW gwww(gw, gww);
+
+      //      std::cout << "making new edges..." << std::endl; 
+      typename GWWW::template EdgeMap<Num> translated_cap(gwww);
+
+      for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
+	typename GW::Edge ew(e, INVALID, false);
+	typename GWWW::Edge ewww(ew, INVALID, false);
+	translated_cap.set(ewww, capacity[e]-lcapacity[e]);
+	//	cout << "t_cap " << g.id(e) << " " 
+	//	     << translated_cap[ewww] << endl;
+      }
+
+      Num expected=0;
+
+      //      std::cout << "making new edges 2..." << std::endl; 
+      for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
+	//	std::cout << "node: " << g.id(n) << std::endl;
+	Num a=0;
+	for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e) {
+	  a+=lcapacity[e];
+	  //	  std::cout << "bee: " << g.id(e) << " " << lcapacity[e] << std::endl;
+	}
+	for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) {
+	  a-=lcapacity[e];
+	  //	  std::cout << "kie: " << g.id(e) << " " << lcapacity[e] << std::endl;
+	}
+	//	std::cout << "excess " << g.id(n) << ": " << a << std::endl;
+	if (0>a) {
+	  typename GWW::Edge e=
+	    gww.addEdge(typename GW::Node(n,INVALID,false), t);
+	  translated_cap.set(typename GWWW::Edge(INVALID, e, true), 
+			     -a);
+	  //	  std::cout << g.id(n) << "->t " << -a << std::endl;
+	}
+	if (0<a) {
+	  typename GWW::Edge e=
+	    gww.addEdge(s, typename GW::Node(n,INVALID,false));
+	  translated_cap.set(typename GWWW::Edge(INVALID, e, true), 
+			     a);
+	  expected+=a;
+	  //	  std::cout << "s->" << g.id(n) << " "<< a <<std:: endl;
+	}
+      }
+
+      //      std::cout << "preflow..." << std::endl; 
+      typename GWWW::template EdgeMap<Num> translated_cost(gwww, 0);
+      for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
+	translated_cost.set(typename GWWW::Edge(
+        typename GW::Edge(e, INVALID, false), INVALID, false), cost[e]);
+      }
+      //      typename GWWW::template EdgeMap<Num> translated_flow(gwww, 0);
+      MinCostFlow<GWWW, typename GWWW::template EdgeMap<Num>, 
+      typename GWWW::template EdgeMap<Num> > 
+      min_cost_flow(gwww, translated_cost, translated_cap, 
+		    s, t);
+      while (min_cost_flow.augment()) { }
+      std::cout << "fv: " << min_cost_flow.flowValue() << std::endl; 
+      std::cout << "expected: " << expected << std::endl; 
+
+      for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
+	typename GW::Edge ew(e, INVALID, false);
+	typename GWWW::Edge ewww(ew, INVALID, false);
+	//	std::cout << g.id(e) << " " << flow[e] << std::endl;
+	flow.set(e, lcapacity[e]+
+		 min_cost_flow.getFlow()[ewww]);
+      }
+      return (min_cost_flow.flowValue()>=expected);
+    }
+    void runByLP() {
+      typedef LPGLPK LPSolver;
+      LPSolver lp;
+      lp.setMinimize();
+      typedef LPSolver::ColIt ColIt;
+      typedef LPSolver::RowIt RowIt;
+      typedef typename Graph::template EdgeMap<ColIt> EdgeIndexMap;
+      EdgeIndexMap edge_index_map(g);
+      PrimalMap<typename Graph::Edge, EdgeIndexMap> lp_flow(lp, edge_index_map);
+      for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
+	ColIt col_it=lp.addCol();
+	edge_index_map.set(e, col_it);
+	if (lcapacity[e]==capacity[e])
+	  lp.setColBounds(col_it, LPSolver::FIXED, lcapacity[e], capacity[e]);
+	else 
+	  lp.setColBounds(col_it, LPSolver::DOUBLE, lcapacity[e], capacity[e]);
+	lp.setObjCoef(col_it, cost[e]);
+      }
+      LPSolver::ColIt col_it;
+      for (lp.col_iter_map.first(col_it, lp.VALID_CLASS); 
+	   lp.col_iter_map.valid(col_it); 
+	   lp.col_iter_map.next(col_it)) {
+//	std::cout << "ize " << lp.col_iter_map[col_it] << std::endl;
+      }
+      for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
+	typename Graph::template EdgeMap<Num> coeffs(g, 0);
+	for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
+	coeffs.set(e, coeffs[e]+1);
+	for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) 
+	coeffs.set(e, coeffs[e]-1);
+	RowIt row_it=lp.addRow();
+	typename std::vector< std::pair<ColIt, double> > row;
+	//std::cout << "node:" <<g.id(n)<<std::endl;
+	for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
+	  if (coeffs[e]!=0) {
+	    //std::cout << " edge:" <<g.id(e)<<" "<<coeffs[e];
+	    row.push_back(std::make_pair(edge_index_map[e], coeffs[e]));
+	  }
+	}
+	//std::cout << std::endl;
+	//std::cout << " " << g.id(n) << " " << row.size() << std::endl;
+	lp.setRowCoeffs(row_it, row.begin(), row.end());
+	lp.setRowBounds(row_it, LPSolver::FIXED, 0.0, 0.0);
+      }
+      lp.solveSimplex();
+      //std::cout << lp.colNum() << std::endl;
+      //std::cout << lp.rowNum() << std::endl;
+      //std::cout << "flow value: "<< lp.getObjVal() << std::endl;
+      for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) 
+      flow.set(e, lp_flow[e]);
+    }
+  };
+
+} // namespace lemon
+
+#endif //LEMON_MIN_COST_GEN_FLOW_H