Index: lemon/Makefile.am
===================================================================
--- lemon/Makefile.am	(revision 674)
+++ lemon/Makefile.am	(revision 686)
@@ -16,5 +16,6 @@
 	lemon/bits/windows.cc
 
-
+nodist_lemon_HEADERS = lemon/config.h	
+	
 lemon_libemon_la_CXXFLAGS = \
 	$(AM_CXXFLAGS) \
@@ -58,9 +59,9 @@
 	lemon/bfs.h \
 	lemon/bin_heap.h \
+	lemon/cbc.h \
 	lemon/circulation.h \
 	lemon/clp.h \
 	lemon/color.h \
 	lemon/concept_check.h \
-	lemon/config.h \
 	lemon/connectivity.h \
 	lemon/counter.h \
Index: lemon/circulation.h
===================================================================
--- lemon/circulation.h	(revision 669)
+++ lemon/circulation.h	(revision 688)
@@ -65,6 +65,6 @@
     typedef SM SupplyMap;
 
-    /// \brief The type of the flow values.
-    typedef typename SupplyMap::Value Flow;
+    /// \brief The type of the flow and supply values.
+    typedef typename SupplyMap::Value Value;
 
     /// \brief The type of the map that stores the flow values.
@@ -73,5 +73,5 @@
     /// It must conform to the \ref concepts::ReadWriteMap "ReadWriteMap"
     /// concept.
-    typedef typename Digraph::template ArcMap<Flow> FlowMap;
+    typedef typename Digraph::template ArcMap<Value> FlowMap;
 
     /// \brief Instantiates a FlowMap.
@@ -105,5 +105,5 @@
     ///
     /// The tolerance used by the algorithm to handle inexact computation.
-    typedef lemon::Tolerance<Flow> Tolerance;
+    typedef lemon::Tolerance<Value> Tolerance;
 
   };
@@ -188,6 +188,6 @@
     ///The type of the digraph the algorithm runs on.
     typedef typename Traits::Digraph Digraph;
-    ///The type of the flow values.
-    typedef typename Traits::Flow Flow;
+    ///The type of the flow and supply values.
+    typedef typename Traits::Value Value;
 
     ///The type of the lower bound map.
@@ -222,5 +222,5 @@
     bool _local_level;
 
-    typedef typename Digraph::template NodeMap<Flow> ExcessMap;
+    typedef typename Digraph::template NodeMap<Value> ExcessMap;
     ExcessMap* _excess;
 
@@ -531,5 +531,5 @@
           (*_excess)[_g.source(e)] -= (*_lo)[e];
         } else {
-          Flow fc = -(*_excess)[_g.target(e)];
+          Value fc = -(*_excess)[_g.target(e)];
           _flow->set(e, fc);
           (*_excess)[_g.target(e)] = 0;
@@ -564,9 +564,9 @@
         int actlevel=(*_level)[act];
         int mlevel=_node_num;
-        Flow exc=(*_excess)[act];
+        Value exc=(*_excess)[act];
 
         for(OutArcIt e(_g,act);e!=INVALID; ++e) {
           Node v = _g.target(e);
-          Flow fc=(*_up)[e]-(*_flow)[e];
+          Value fc=(*_up)[e]-(*_flow)[e];
           if(!_tol.positive(fc)) continue;
           if((*_level)[v]<actlevel) {
@@ -592,5 +592,5 @@
         for(InArcIt e(_g,act);e!=INVALID; ++e) {
           Node v = _g.source(e);
-          Flow fc=(*_flow)[e]-(*_lo)[e];
+          Value fc=(*_flow)[e]-(*_lo)[e];
           if(!_tol.positive(fc)) continue;
           if((*_level)[v]<actlevel) {
@@ -662,11 +662,11 @@
     ///@{
 
-    /// \brief Returns the flow on the given arc.
-    ///
-    /// Returns the flow on the given arc.
+    /// \brief Returns the flow value on the given arc.
+    ///
+    /// Returns the flow value on the given arc.
     ///
     /// \pre Either \ref run() or \ref init() must be called before
     /// using this function.
-    Flow flow(const Arc& arc) const {
+    Value flow(const Arc& arc) const {
       return (*_flow)[arc];
     }
@@ -751,5 +751,5 @@
       for(NodeIt n(_g);n!=INVALID;++n)
         {
-          Flow dif=-(*_supply)[n];
+          Value dif=-(*_supply)[n];
           for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
           for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
@@ -766,8 +766,8 @@
     bool checkBarrier() const
     {
-      Flow delta=0;
-      Flow inf_cap = std::numeric_limits<Flow>::has_infinity ?
-        std::numeric_limits<Flow>::infinity() :
-        std::numeric_limits<Flow>::max();
+      Value delta=0;
+      Value inf_cap = std::numeric_limits<Value>::has_infinity ?
+        std::numeric_limits<Value>::infinity() :
+        std::numeric_limits<Value>::max();
       for(NodeIt n(_g);n!=INVALID;++n)
         if(barrier(n))
Index: lemon/connectivity.h
===================================================================
--- lemon/connectivity.h	(revision 633)
+++ lemon/connectivity.h	(revision 694)
@@ -1371,4 +1371,5 @@
     typedef typename Graph::NodeIt NodeIt;
     typedef typename Graph::Arc Arc;
+    if (NodeIt(graph) == INVALID) return true;
     Dfs<Graph> dfs(graph);
     dfs.init();
@@ -1546,15 +1547,14 @@
   ///
   /// Returns true when there are not parallel edges in the graph.
-  template <typename Digraph>
-  bool parallelFree(const Digraph& digraph) {
-    typename Digraph::template NodeMap<bool> reached(digraph, false);
-    for (typename Digraph::NodeIt n(digraph); n != INVALID; ++n) {
-      for (typename Digraph::OutArcIt a(digraph, n); a != INVALID; ++a) {
-        if (reached[digraph.target(a)]) return false;
-        reached.set(digraph.target(a), true);
-      }
-      for (typename Digraph::OutArcIt a(digraph, n); a != INVALID; ++a) {
-        reached.set(digraph.target(a), false);
-      }
+  template <typename Graph>
+  bool parallelFree(const Graph& graph) {
+    typename Graph::template NodeMap<int> reached(graph, 0);
+    int cnt = 1;
+    for (typename Graph::NodeIt n(graph); n != INVALID; ++n) {
+      for (typename Graph::OutArcIt a(graph, n); a != INVALID; ++a) {
+        if (reached[graph.target(a)] == cnt) return false;
+        reached[graph.target(a)] = cnt;
+      }
+      ++cnt;
     }
     return true;
@@ -1566,17 +1566,15 @@
   /// Returns true when there are not loop edges and parallel edges in
   /// the graph.
-  template <typename Digraph>
-  bool simpleDigraph(const Digraph& digraph) {
-    typename Digraph::template NodeMap<bool> reached(digraph, false);
-    for (typename Digraph::NodeIt n(digraph); n != INVALID; ++n) {
-      reached.set(n, true);
-      for (typename Digraph::OutArcIt a(digraph, n); a != INVALID; ++a) {
-        if (reached[digraph.target(a)]) return false;
-        reached.set(digraph.target(a), true);
-      }
-      for (typename Digraph::OutArcIt a(digraph, n); a != INVALID; ++a) {
-        reached.set(digraph.target(a), false);
-      }
-      reached.set(n, false);
+  template <typename Graph>
+  bool simpleGraph(const Graph& graph) {
+    typename Graph::template NodeMap<int> reached(graph, 0);
+    int cnt = 1;
+    for (typename Graph::NodeIt n(graph); n != INVALID; ++n) {
+      reached[n] = cnt;
+      for (typename Graph::OutArcIt a(graph, n); a != INVALID; ++a) {
+        if (reached[graph.target(a)] == cnt) return false;
+        reached[graph.target(a)] = cnt;
+      }
+      ++cnt;
     }
     return true;
Index: lemon/core.h
===================================================================
--- lemon/core.h	(revision 664)
+++ lemon/core.h	(revision 686)
@@ -23,5 +23,5 @@
 #include <algorithm>
 
-#include <lemon/core.h>
+#include <lemon/config.h>
 #include <lemon/bits/enable_if.h>
 #include <lemon/bits/traits.h>
Index: lemon/network_simplex.h
===================================================================
--- lemon/network_simplex.h	(revision 665)
+++ lemon/network_simplex.h	(revision 690)
@@ -31,7 +31,4 @@
 #include <lemon/core.h>
 #include <lemon/math.h>
-#include <lemon/maps.h>
-#include <lemon/circulation.h>
-#include <lemon/adaptors.h>
 
 namespace lemon {
@@ -51,12 +48,17 @@
   /// In general this class is the fastest implementation available
   /// in LEMON for the minimum cost flow problem.
-  /// Moreover it supports both direction of the supply/demand inequality
-  /// constraints. For more information see \ref ProblemType.
+  /// Moreover it supports both directions of the supply/demand inequality
+  /// constraints. For more information see \ref SupplyType.
+  ///
+  /// Most of the parameters of the problem (except for the digraph)
+  /// can be given using separate functions, and the algorithm can be
+  /// executed using the \ref run() function. If some parameters are not
+  /// specified, then default values will be used.
   ///
   /// \tparam GR The digraph type the algorithm runs on.
-  /// \tparam F The value type used for flow amounts, capacity bounds
+  /// \tparam V The value type used for flow amounts, capacity bounds
   /// and supply values in the algorithm. By default it is \c int.
   /// \tparam C The value type used for costs and potentials in the
-  /// algorithm. By default it is the same as \c F.
+  /// algorithm. By default it is the same as \c V.
   ///
   /// \warning Both value types must be signed and all input data must
@@ -66,88 +68,51 @@
   /// implementations, from which the most efficient one is used
   /// by default. For more information see \ref PivotRule.
-  template <typename GR, typename F = int, typename C = F>
+  template <typename GR, typename V = int, typename C = V>
   class NetworkSimplex
   {
   public:
 
-    /// The flow type of the algorithm
-    typedef F Flow;
-    /// The cost type of the algorithm
+    /// The type of the flow amounts, capacity bounds and supply values
+    typedef V Value;
+    /// The type of the arc costs
     typedef C Cost;
-#ifdef DOXYGEN
-    /// The type of the flow map
-    typedef GR::ArcMap<Flow> FlowMap;
-    /// The type of the potential map
-    typedef GR::NodeMap<Cost> PotentialMap;
-#else
-    /// The type of the flow map
-    typedef typename GR::template ArcMap<Flow> FlowMap;
-    /// The type of the potential map
-    typedef typename GR::template NodeMap<Cost> PotentialMap;
-#endif
 
   public:
 
-    /// \brief Enum type for selecting the pivot rule.
-    ///
-    /// Enum type for selecting the pivot rule for the \ref run()
-    /// function.
-    ///
-    /// \ref NetworkSimplex provides five different pivot rule
-    /// implementations that significantly affect the running time
-    /// of the algorithm.
-    /// By default \ref BLOCK_SEARCH "Block Search" is used, which
-    /// proved to be the most efficient and the most robust on various
-    /// test inputs according to our benchmark tests.
-    /// However another pivot rule can be selected using the \ref run()
-    /// function with the proper parameter.
-    enum PivotRule {
-
-      /// The First Eligible pivot rule.
-      /// The next eligible arc is selected in a wraparound fashion
-      /// in every iteration.
-      FIRST_ELIGIBLE,
-
-      /// The Best Eligible pivot rule.
-      /// The best eligible arc is selected in every iteration.
-      BEST_ELIGIBLE,
-
-      /// The Block Search pivot rule.
-      /// A specified number of arcs are examined in every iteration
-      /// in a wraparound fashion and the best eligible arc is selected
-      /// from this block.
-      BLOCK_SEARCH,
-
-      /// The Candidate List pivot rule.
-      /// In a major iteration a candidate list is built from eligible arcs
-      /// in a wraparound fashion and in the following minor iterations
-      /// the best eligible arc is selected from this list.
-      CANDIDATE_LIST,
-
-      /// The Altering Candidate List pivot rule.
-      /// It is a modified version of the Candidate List method.
-      /// It keeps only the several best eligible arcs from the former
-      /// candidate list and extends this list in every iteration.
-      ALTERING_LIST
+    /// \brief Problem type constants for the \c run() function.
+    ///
+    /// Enum type containing the problem type constants that can be
+    /// returned by the \ref run() function of the algorithm.
+    enum ProblemType {
+      /// The problem has no feasible solution (flow).
+      INFEASIBLE,
+      /// The problem has optimal solution (i.e. it is feasible and
+      /// bounded), and the algorithm has found optimal flow and node
+      /// potentials (primal and dual solutions).
+      OPTIMAL,
+      /// The objective function of the problem is unbounded, i.e.
+      /// there is a directed cycle having negative total cost and
+      /// infinite upper bound.
+      UNBOUNDED
     };
     
-    /// \brief Enum type for selecting the problem type.
-    ///
-    /// Enum type for selecting the problem type, i.e. the direction of
-    /// the inequalities in the supply/demand constraints of the
-    /// \ref min_cost_flow "minimum cost flow problem".
-    ///
-    /// The default problem type is \c GEQ, since this form is supported
+    /// \brief Constants for selecting the type of the supply constraints.
+    ///
+    /// Enum type containing constants for selecting the supply type,
+    /// i.e. the direction of the inequalities in the supply/demand
+    /// constraints of the \ref min_cost_flow "minimum cost flow problem".
+    ///
+    /// The default supply type is \c GEQ, since this form is supported
     /// by other minimum cost flow algorithms and the \ref Circulation
-    /// algorithm as well.
-    /// The \c LEQ problem type can be selected using the \ref problemType()
+    /// algorithm, as well.
+    /// The \c LEQ problem type can be selected using the \ref supplyType()
     /// function.
     ///
-    /// Note that the equality form is a special case of both problem type.
-    enum ProblemType {
-
-      /// This option means that there are "<em>greater or equal</em>"
-      /// constraints in the defintion, i.e. the exact formulation of the
-      /// problem is the following.
+    /// Note that the equality form is a special case of both supply types.
+    enum SupplyType {
+
+      /// This option means that there are <em>"greater or equal"</em>
+      /// supply/demand constraints in the definition, i.e. the exact
+      /// formulation of the problem is the following.
       /**
           \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
@@ -165,7 +130,7 @@
       CARRY_SUPPLIES = GEQ,
 
-      /// This option means that there are "<em>less or equal</em>"
-      /// constraints in the defintion, i.e. the exact formulation of the
-      /// problem is the following.
+      /// This option means that there are <em>"less or equal"</em>
+      /// supply/demand constraints in the definition, i.e. the exact
+      /// formulation of the problem is the following.
       /**
           \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
@@ -183,12 +148,51 @@
       SATISFY_DEMANDS = LEQ
     };
-
+    
+    /// \brief Constants for selecting the pivot rule.
+    ///
+    /// Enum type containing constants for selecting the pivot rule for
+    /// the \ref run() function.
+    ///
+    /// \ref NetworkSimplex provides five different pivot rule
+    /// implementations that significantly affect the running time
+    /// of the algorithm.
+    /// By default \ref BLOCK_SEARCH "Block Search" is used, which
+    /// proved to be the most efficient and the most robust on various
+    /// test inputs according to our benchmark tests.
+    /// However another pivot rule can be selected using the \ref run()
+    /// function with the proper parameter.
+    enum PivotRule {
+
+      /// The First Eligible pivot rule.
+      /// The next eligible arc is selected in a wraparound fashion
+      /// in every iteration.
+      FIRST_ELIGIBLE,
+
+      /// The Best Eligible pivot rule.
+      /// The best eligible arc is selected in every iteration.
+      BEST_ELIGIBLE,
+
+      /// The Block Search pivot rule.
+      /// A specified number of arcs are examined in every iteration
+      /// in a wraparound fashion and the best eligible arc is selected
+      /// from this block.
+      BLOCK_SEARCH,
+
+      /// The Candidate List pivot rule.
+      /// In a major iteration a candidate list is built from eligible arcs
+      /// in a wraparound fashion and in the following minor iterations
+      /// the best eligible arc is selected from this list.
+      CANDIDATE_LIST,
+
+      /// The Altering Candidate List pivot rule.
+      /// It is a modified version of the Candidate List method.
+      /// It keeps only the several best eligible arcs from the former
+      /// candidate list and extends this list in every iteration.
+      ALTERING_LIST
+    };
+    
   private:
 
     TEMPLATE_DIGRAPH_TYPEDEFS(GR);
-
-    typedef typename GR::template ArcMap<Flow> FlowArcMap;
-    typedef typename GR::template ArcMap<Cost> CostArcMap;
-    typedef typename GR::template NodeMap<Flow> FlowNodeMap;
 
     typedef std::vector<Arc> ArcVector;
@@ -196,5 +200,5 @@
     typedef std::vector<int> IntVector;
     typedef std::vector<bool> BoolVector;
-    typedef std::vector<Flow> FlowVector;
+    typedef std::vector<Value> ValueVector;
     typedef std::vector<Cost> CostVector;
 
@@ -214,30 +218,21 @@
 
     // Parameters of the problem
-    FlowArcMap *_plower;
-    FlowArcMap *_pupper;
-    CostArcMap *_pcost;
-    FlowNodeMap *_psupply;
-    bool _pstsup;
-    Node _psource, _ptarget;
-    Flow _pstflow;
-    ProblemType _ptype;
-
-    // Result maps
-    FlowMap *_flow_map;
-    PotentialMap *_potential_map;
-    bool _local_flow;
-    bool _local_potential;
+    bool _have_lower;
+    SupplyType _stype;
+    Value _sum_supply;
 
     // Data structures for storing the digraph
     IntNodeMap _node_id;
-    ArcVector _arc_ref;
+    IntArcMap _arc_id;
     IntVector _source;
     IntVector _target;
 
     // Node and arc data
-    FlowVector _cap;
+    ValueVector _lower;
+    ValueVector _upper;
+    ValueVector _cap;
     CostVector _cost;
-    FlowVector _supply;
-    FlowVector _flow;
+    ValueVector _supply;
+    ValueVector _flow;
     CostVector _pi;
 
@@ -258,5 +253,14 @@
     int first, second, right, last;
     int stem, par_stem, new_stem;
-    Flow delta;
+    Value delta;
+
+  public:
+  
+    /// \brief Constant for infinite upper bounds (capacities).
+    ///
+    /// Constant for infinite upper bounds (capacities).
+    /// It is \c std::numeric_limits<Value>::infinity() if available,
+    /// \c std::numeric_limits<Value>::max() otherwise.
+    const Value INF;
 
   private:
@@ -660,23 +664,66 @@
     /// \param graph The digraph the algorithm runs on.
     NetworkSimplex(const GR& graph) :
-      _graph(graph),
-      _plower(NULL), _pupper(NULL), _pcost(NULL),
-      _psupply(NULL), _pstsup(false), _ptype(GEQ),
-      _flow_map(NULL), _potential_map(NULL),
-      _local_flow(false), _local_potential(false),
-      _node_id(graph)
+      _graph(graph), _node_id(graph), _arc_id(graph),
+      INF(std::numeric_limits<Value>::has_infinity ?
+          std::numeric_limits<Value>::infinity() :
+          std::numeric_limits<Value>::max())
     {
-      LEMON_ASSERT(std::numeric_limits<Flow>::is_integer &&
-                   std::numeric_limits<Flow>::is_signed,
-        "The flow type of NetworkSimplex must be signed integer");
-      LEMON_ASSERT(std::numeric_limits<Cost>::is_integer &&
-                   std::numeric_limits<Cost>::is_signed,
-        "The cost type of NetworkSimplex must be signed integer");
-    }
-
-    /// Destructor.
-    ~NetworkSimplex() {
-      if (_local_flow) delete _flow_map;
-      if (_local_potential) delete _potential_map;
+      // Check the value types
+      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
+        "The flow type of NetworkSimplex must be signed");
+      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
+        "The cost type of NetworkSimplex must be signed");
+        
+      // Resize vectors
+      _node_num = countNodes(_graph);
+      _arc_num = countArcs(_graph);
+      int all_node_num = _node_num + 1;
+      int all_arc_num = _arc_num + _node_num;
+
+      _source.resize(all_arc_num);
+      _target.resize(all_arc_num);
+
+      _lower.resize(all_arc_num);
+      _upper.resize(all_arc_num);
+      _cap.resize(all_arc_num);
+      _cost.resize(all_arc_num);
+      _supply.resize(all_node_num);
+      _flow.resize(all_arc_num);
+      _pi.resize(all_node_num);
+
+      _parent.resize(all_node_num);
+      _pred.resize(all_node_num);
+      _forward.resize(all_node_num);
+      _thread.resize(all_node_num);
+      _rev_thread.resize(all_node_num);
+      _succ_num.resize(all_node_num);
+      _last_succ.resize(all_node_num);
+      _state.resize(all_arc_num);
+
+      // Copy the graph (store the arcs in a mixed order)
+      int i = 0;
+      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
+        _node_id[n] = i;
+      }
+      int k = std::max(int(std::sqrt(double(_arc_num))), 10);
+      i = 0;
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        _arc_id[a] = i;
+        _source[i] = _node_id[_graph.source(a)];
+        _target[i] = _node_id[_graph.target(a)];
+        if ((i += k) >= _arc_num) i = (i % k) + 1;
+      }
+      
+      // Initialize maps
+      for (int i = 0; i != _node_num; ++i) {
+        _supply[i] = 0;
+      }
+      for (int i = 0; i != _arc_num; ++i) {
+        _lower[i] = 0;
+        _upper[i] = INF;
+        _cost[i] = 1;
+      }
+      _have_lower = false;
+      _stype = GEQ;
     }
 
@@ -690,19 +737,17 @@
     ///
     /// This function sets the lower bounds on the arcs.
-    /// If neither this function nor \ref boundMaps() is used before
-    /// calling \ref run(), the lower bounds will be set to zero
-    /// on all arcs.
+    /// If it is not used before calling \ref run(), the lower bounds
+    /// will be set to zero on all arcs.
     ///
     /// \param map An arc map storing the lower bounds.
-    /// Its \c Value type must be convertible to the \c Flow type
+    /// Its \c Value type must be convertible to the \c Value type
     /// of the algorithm.
     ///
     /// \return <tt>(*this)</tt>
-    template <typename LOWER>
-    NetworkSimplex& lowerMap(const LOWER& map) {
-      delete _plower;
-      _plower = new FlowArcMap(_graph);
+    template <typename LowerMap>
+    NetworkSimplex& lowerMap(const LowerMap& map) {
+      _have_lower = true;
       for (ArcIt a(_graph); a != INVALID; ++a) {
-        (*_plower)[a] = map[a];
+        _lower[_arc_id[a]] = map[a];
       }
       return *this;
@@ -712,59 +757,19 @@
     ///
     /// This function sets the upper bounds (capacities) on the arcs.
-    /// If none of the functions \ref upperMap(), \ref capacityMap()
-    /// and \ref boundMaps() is used before calling \ref run(),
-    /// the upper bounds (capacities) will be set to
-    /// \c std::numeric_limits<Flow>::max() on all arcs.
+    /// If it is not used before calling \ref run(), the upper bounds
+    /// will be set to \ref INF on all arcs (i.e. the flow value will be
+    /// unbounded from above on each arc).
     ///
     /// \param map An arc map storing the upper bounds.
-    /// Its \c Value type must be convertible to the \c Flow type
+    /// Its \c Value type must be convertible to the \c Value type
     /// of the algorithm.
     ///
     /// \return <tt>(*this)</tt>
-    template<typename UPPER>
-    NetworkSimplex& upperMap(const UPPER& map) {
-      delete _pupper;
-      _pupper = new FlowArcMap(_graph);
+    template<typename UpperMap>
+    NetworkSimplex& upperMap(const UpperMap& map) {
       for (ArcIt a(_graph); a != INVALID; ++a) {
-        (*_pupper)[a] = map[a];
+        _upper[_arc_id[a]] = map[a];
       }
       return *this;
-    }
-
-    /// \brief Set the upper bounds (capacities) on the arcs.
-    ///
-    /// This function sets the upper bounds (capacities) on the arcs.
-    /// It is just an alias for \ref upperMap().
-    ///
-    /// \return <tt>(*this)</tt>
-    template<typename CAP>
-    NetworkSimplex& capacityMap(const CAP& map) {
-      return upperMap(map);
-    }
-
-    /// \brief Set the lower and upper bounds on the arcs.
-    ///
-    /// This function sets the lower and upper bounds on the arcs.
-    /// If neither this function nor \ref lowerMap() is used before
-    /// calling \ref run(), the lower bounds will be set to zero
-    /// on all arcs.
-    /// If none of the functions \ref upperMap(), \ref capacityMap()
-    /// and \ref boundMaps() is used before calling \ref run(),
-    /// the upper bounds (capacities) will be set to
-    /// \c std::numeric_limits<Flow>::max() on all arcs.
-    ///
-    /// \param lower An arc map storing the lower bounds.
-    /// \param upper An arc map storing the upper bounds.
-    ///
-    /// The \c Value type of the maps must be convertible to the
-    /// \c Flow type of the algorithm.
-    ///
-    /// \note This function is just a shortcut of calling \ref lowerMap()
-    /// and \ref upperMap() separately.
-    ///
-    /// \return <tt>(*this)</tt>
-    template <typename LOWER, typename UPPER>
-    NetworkSimplex& boundMaps(const LOWER& lower, const UPPER& upper) {
-      return lowerMap(lower).upperMap(upper);
     }
 
@@ -780,10 +785,8 @@
     ///
     /// \return <tt>(*this)</tt>
-    template<typename COST>
-    NetworkSimplex& costMap(const COST& map) {
-      delete _pcost;
-      _pcost = new CostArcMap(_graph);
+    template<typename CostMap>
+    NetworkSimplex& costMap(const CostMap& map) {
       for (ArcIt a(_graph); a != INVALID; ++a) {
-        (*_pcost)[a] = map[a];
+        _cost[_arc_id[a]] = map[a];
       }
       return *this;
@@ -798,15 +801,12 @@
     ///
     /// \param map A node map storing the supply values.
-    /// Its \c Value type must be convertible to the \c Flow type
+    /// Its \c Value type must be convertible to the \c Value type
     /// of the algorithm.
     ///
     /// \return <tt>(*this)</tt>
-    template<typename SUP>
-    NetworkSimplex& supplyMap(const SUP& map) {
-      delete _psupply;
-      _pstsup = false;
-      _psupply = new FlowNodeMap(_graph);
+    template<typename SupplyMap>
+    NetworkSimplex& supplyMap(const SupplyMap& map) {
       for (NodeIt n(_graph); n != INVALID; ++n) {
-        (*_psupply)[n] = map[n];
+        _supply[_node_id[n]] = map[n];
       }
       return *this;
@@ -821,4 +821,8 @@
     /// (It makes sense only if non-zero lower bounds are given.)
     ///
+    /// Using this function has the same effect as using \ref supplyMap()
+    /// with such a map in which \c k is assigned to \c s, \c -k is
+    /// assigned to \c t and all other nodes have zero supply value.
+    ///
     /// \param s The source node.
     /// \param t The target node.
@@ -827,63 +831,27 @@
     ///
     /// \return <tt>(*this)</tt>
-    NetworkSimplex& stSupply(const Node& s, const Node& t, Flow k) {
-      delete _psupply;
-      _psupply = NULL;
-      _pstsup = true;
-      _psource = s;
-      _ptarget = t;
-      _pstflow = k;
+    NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) {
+      for (int i = 0; i != _node_num; ++i) {
+        _supply[i] = 0;
+      }
+      _supply[_node_id[s]] =  k;
+      _supply[_node_id[t]] = -k;
       return *this;
     }
     
-    /// \brief Set the problem type.
-    ///
-    /// This function sets the problem type for the algorithm.
-    /// If it is not used before calling \ref run(), the \ref GEQ problem
+    /// \brief Set the type of the supply constraints.
+    ///
+    /// This function sets the type of the supply/demand constraints.
+    /// If it is not used before calling \ref run(), the \ref GEQ supply
     /// type will be used.
     ///
-    /// For more information see \ref ProblemType.
+    /// For more information see \ref SupplyType.
     ///
     /// \return <tt>(*this)</tt>
-    NetworkSimplex& problemType(ProblemType problem_type) {
-      _ptype = problem_type;
+    NetworkSimplex& supplyType(SupplyType supply_type) {
+      _stype = supply_type;
       return *this;
     }
 
-    /// \brief Set the flow map.
-    ///
-    /// This function sets the flow map.
-    /// If it is not used before calling \ref run(), an instance will
-    /// be allocated automatically. The destructor deallocates this
-    /// automatically allocated map, of course.
-    ///
-    /// \return <tt>(*this)</tt>
-    NetworkSimplex& flowMap(FlowMap& map) {
-      if (_local_flow) {
-        delete _flow_map;
-        _local_flow = false;
-      }
-      _flow_map = &map;
-      return *this;
-    }
-
-    /// \brief Set the potential map.
-    ///
-    /// This function sets the potential map, which is used for storing
-    /// the dual solution.
-    /// If it is not used before calling \ref run(), an instance will
-    /// be allocated automatically. The destructor deallocates this
-    /// automatically allocated map, of course.
-    ///
-    /// \return <tt>(*this)</tt>
-    NetworkSimplex& potentialMap(PotentialMap& map) {
-      if (_local_potential) {
-        delete _potential_map;
-        _local_potential = false;
-      }
-      _potential_map = &map;
-      return *this;
-    }
-    
     /// @}
 
@@ -897,11 +865,10 @@
     /// This function runs the algorithm.
     /// The paramters can be specified using functions \ref lowerMap(),
-    /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
-    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), 
-    /// \ref problemType(), \ref flowMap() and \ref potentialMap().
+    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), 
+    /// \ref supplyType().
     /// For example,
     /// \code
     ///   NetworkSimplex<ListDigraph> ns(graph);
-    ///   ns.boundMaps(lower, upper).costMap(cost)
+    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
     ///     .supplyMap(sup).run();
     /// \endcode
@@ -911,11 +878,22 @@
     /// \ref reset() is called, thus only the modified parameters
     /// have to be set again. See \ref reset() for examples.
+    /// However the underlying digraph must not be modified after this
+    /// class have been constructed, since it copies and extends the graph.
     ///
     /// \param pivot_rule The pivot rule that will be used during the
     /// algorithm. For more information see \ref PivotRule.
     ///
-    /// \return \c true if a feasible flow can be found.
-    bool run(PivotRule pivot_rule = BLOCK_SEARCH) {
-      return init() && start(pivot_rule);
+    /// \return \c INFEASIBLE if no feasible flow exists,
+    /// \n \c OPTIMAL if the problem has optimal solution
+    /// (i.e. it is feasible and bounded), and the algorithm has found
+    /// optimal flow and node potentials (primal and dual solutions),
+    /// \n \c UNBOUNDED if the objective function of the problem is
+    /// unbounded, i.e. there is a directed cycle having negative total
+    /// cost and infinite upper bound.
+    ///
+    /// \see ProblemType, PivotRule
+    ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
+      if (!init()) return INFEASIBLE;
+      return start(pivot_rule);
     }
 
@@ -924,11 +902,11 @@
     /// This function resets all the paramaters that have been given
     /// before using functions \ref lowerMap(), \ref upperMap(),
-    /// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
-    /// \ref supplyMap(), \ref stSupply(), \ref problemType(), 
-    /// \ref flowMap() and \ref potentialMap().
+    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType().
     ///
     /// It is useful for multiple run() calls. If this function is not
     /// used, all the parameters given before are kept for the next
     /// \ref run() call.
+    /// However the underlying digraph must not be modified after this
+    /// class have been constructed, since it copies and extends the graph.
     ///
     /// For example,
@@ -937,5 +915,5 @@
     ///
     ///   // First run
-    ///   ns.lowerMap(lower).capacityMap(cap).costMap(cost)
+    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
     ///     .supplyMap(sup).run();
     ///
@@ -948,5 +926,5 @@
     ///   // (the lower bounds will be set to zero on all arcs)
     ///   ns.reset();
-    ///   ns.capacityMap(cap).costMap(cost)
+    ///   ns.upperMap(capacity).costMap(cost)
     ///     .supplyMap(sup).run();
     /// \endcode
@@ -954,21 +932,14 @@
     /// \return <tt>(*this)</tt>
     NetworkSimplex& reset() {
-      delete _plower;
-      delete _pupper;
-      delete _pcost;
-      delete _psupply;
-      _plower = NULL;
-      _pupper = NULL;
-      _pcost = NULL;
-      _psupply = NULL;
-      _pstsup = false;
-      _ptype = GEQ;
-      if (_local_flow) delete _flow_map;
-      if (_local_potential) delete _potential_map;
-      _flow_map = NULL;
-      _potential_map = NULL;
-      _local_flow = false;
-      _local_potential = false;
-
+      for (int i = 0; i != _node_num; ++i) {
+        _supply[i] = 0;
+      }
+      for (int i = 0; i != _arc_num; ++i) {
+        _lower[i] = 0;
+        _upper[i] = INF;
+        _cost[i] = 1;
+      }
+      _have_lower = false;
+      _stype = GEQ;
       return *this;
     }
@@ -986,5 +957,5 @@
     ///
     /// This function returns the total cost of the found flow.
-    /// The complexity of the function is O(e).
+    /// Its complexity is O(e).
     ///
     /// \note The return type of the function can be specified as a
@@ -998,13 +969,10 @@
     ///
     /// \pre \ref run() must be called before using this function.
-    template <typename Num>
-    Num totalCost() const {
-      Num c = 0;
-      if (_pcost) {
-        for (ArcIt e(_graph); e != INVALID; ++e)
-          c += (*_flow_map)[e] * (*_pcost)[e];
-      } else {
-        for (ArcIt e(_graph); e != INVALID; ++e)
-          c += (*_flow_map)[e];
+    template <typename Number>
+    Number totalCost() const {
+      Number c = 0;
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        int i = _arc_id[a];
+        c += Number(_flow[i]) * Number(_cost[i]);
       }
       return c;
@@ -1022,16 +990,20 @@
     ///
     /// \pre \ref run() must be called before using this function.
-    Flow flow(const Arc& a) const {
-      return (*_flow_map)[a];
-    }
-
-    /// \brief Return a const reference to the flow map.
-    ///
-    /// This function returns a const reference to an arc map storing
-    /// the found flow.
+    Value flow(const Arc& a) const {
+      return _flow[_arc_id[a]];
+    }
+
+    /// \brief Return the flow map (the primal solution).
+    ///
+    /// This function copies the flow value on each arc into the given
+    /// map. The \c Value type of the algorithm must be convertible to
+    /// the \c Value type of the map.
     ///
     /// \pre \ref run() must be called before using this function.
-    const FlowMap& flowMap() const {
-      return *_flow_map;
+    template <typename FlowMap>
+    void flowMap(FlowMap &map) const {
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        map.set(a, _flow[_arc_id[a]]);
+      }
     }
 
@@ -1043,17 +1015,20 @@
     /// \pre \ref run() must be called before using this function.
     Cost potential(const Node& n) const {
-      return (*_potential_map)[n];
-    }
-
-    /// \brief Return a const reference to the potential map
-    /// (the dual solution).
-    ///
-    /// This function returns a const reference to a node map storing
-    /// the found potentials, which form the dual solution of the
-    /// \ref min_cost_flow "minimum cost flow" problem.
+      return _pi[_node_id[n]];
+    }
+
+    /// \brief Return the potential map (the dual solution).
+    ///
+    /// This function copies the potential (dual value) of each node
+    /// into the given map.
+    /// The \c Cost type of the algorithm must be convertible to the
+    /// \c Value type of the map.
     ///
     /// \pre \ref run() must be called before using this function.
-    const PotentialMap& potentialMap() const {
-      return *_potential_map;
+    template <typename PotentialMap>
+    void potentialMap(PotentialMap &map) const {
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        map.set(n, _pi[_node_id[n]]);
+      }
     }
 
@@ -1064,154 +1039,50 @@
     // Initialize internal data structures
     bool init() {
-      // Initialize result maps
-      if (!_flow_map) {
-        _flow_map = new FlowMap(_graph);
-        _local_flow = true;
-      }
-      if (!_potential_map) {
-        _potential_map = new PotentialMap(_graph);
-        _local_potential = true;
-      }
-
-      // Initialize vectors
-      _node_num = countNodes(_graph);
-      _arc_num = countArcs(_graph);
-      int all_node_num = _node_num + 1;
-      int all_arc_num = _arc_num + _node_num;
       if (_node_num == 0) return false;
 
-      _arc_ref.resize(_arc_num);
-      _source.resize(all_arc_num);
-      _target.resize(all_arc_num);
-
-      _cap.resize(all_arc_num);
-      _cost.resize(all_arc_num);
-      _supply.resize(all_node_num);
-      _flow.resize(all_arc_num);
-      _pi.resize(all_node_num);
-
-      _parent.resize(all_node_num);
-      _pred.resize(all_node_num);
-      _forward.resize(all_node_num);
-      _thread.resize(all_node_num);
-      _rev_thread.resize(all_node_num);
-      _succ_num.resize(all_node_num);
-      _last_succ.resize(all_node_num);
-      _state.resize(all_arc_num);
-
-      // Initialize node related data
-      bool valid_supply = true;
-      Flow sum_supply = 0;
-      if (!_pstsup && !_psupply) {
-        _pstsup = true;
-        _psource = _ptarget = NodeIt(_graph);
-        _pstflow = 0;
-      }
-      if (_psupply) {
-        int i = 0;
-        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
-          _node_id[n] = i;
-          _supply[i] = (*_psupply)[n];
-          sum_supply += _supply[i];
-        }
-        valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
-                       (_ptype == LEQ && sum_supply >= 0);
+      // Check the sum of supply values
+      _sum_supply = 0;
+      for (int i = 0; i != _node_num; ++i) {
+        _sum_supply += _supply[i];
+      }
+      if ( !((_stype == GEQ && _sum_supply <= 0) ||
+             (_stype == LEQ && _sum_supply >= 0)) ) return false;
+
+      // Remove non-zero lower bounds
+      if (_have_lower) {
+        for (int i = 0; i != _arc_num; ++i) {
+          Value c = _lower[i];
+          if (c >= 0) {
+            _cap[i] = _upper[i] < INF ? _upper[i] - c : INF;
+          } else {
+            _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF;
+          }
+          _supply[_source[i]] -= c;
+          _supply[_target[i]] += c;
+        }
       } else {
-        int i = 0;
-        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
-          _node_id[n] = i;
-          _supply[i] = 0;
-        }
-        _supply[_node_id[_psource]] =  _pstflow;
-        _supply[_node_id[_ptarget]] = -_pstflow;
-      }
-      if (!valid_supply) return false;
-
-      // Infinite capacity value
-      Flow inf_cap =
-        std::numeric_limits<Flow>::has_infinity ?
-        std::numeric_limits<Flow>::infinity() :
-        std::numeric_limits<Flow>::max();
+        for (int i = 0; i != _arc_num; ++i) {
+          _cap[i] = _upper[i];
+        }
+      }
 
       // Initialize artifical cost
-      Cost art_cost;
+      Cost ART_COST;
       if (std::numeric_limits<Cost>::is_exact) {
-        art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
+        ART_COST = std::numeric_limits<Cost>::max() / 4 + 1;
       } else {
-        art_cost = std::numeric_limits<Cost>::min();
+        ART_COST = std::numeric_limits<Cost>::min();
         for (int i = 0; i != _arc_num; ++i) {
-          if (_cost[i] > art_cost) art_cost = _cost[i];
-        }
-        art_cost = (art_cost + 1) * _node_num;
-      }
-
-      // Run Circulation to check if a feasible solution exists
-      typedef ConstMap<Arc, Flow> ConstArcMap;
-      ConstArcMap zero_arc_map(0), inf_arc_map(inf_cap);
-      FlowNodeMap *csup = NULL;
-      bool local_csup = false;
-      if (_psupply) {
-        csup = _psupply;
-      } else {
-        csup = new FlowNodeMap(_graph, 0);
-        (*csup)[_psource] =  _pstflow;
-        (*csup)[_ptarget] = -_pstflow;
-        local_csup = true;
-      }
-      bool circ_result = false;
-      if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {
-        // GEQ problem type
-        if (_plower) {
-          if (_pupper) {
-            Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>
-              circ(_graph, *_plower, *_pupper, *csup);
-            circ_result = circ.run();
-          } else {
-            Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>
-              circ(_graph, *_plower, inf_arc_map, *csup);
-            circ_result = circ.run();
-          }
-        } else {
-          if (_pupper) {
-            Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>
-              circ(_graph, zero_arc_map, *_pupper, *csup);
-            circ_result = circ.run();
-          } else {
-            Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>
-              circ(_graph, zero_arc_map, inf_arc_map, *csup);
-            circ_result = circ.run();
-          }
-        }
-      } else {
-        // LEQ problem type
-        typedef ReverseDigraph<const GR> RevGraph;
-        typedef NegMap<FlowNodeMap> NegNodeMap;
-        RevGraph rgraph(_graph);
-        NegNodeMap neg_csup(*csup);
-        if (_plower) {
-          if (_pupper) {
-            Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>
-              circ(rgraph, *_plower, *_pupper, neg_csup);
-            circ_result = circ.run();
-          } else {
-            Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>
-              circ(rgraph, *_plower, inf_arc_map, neg_csup);
-            circ_result = circ.run();
-          }
-        } else {
-          if (_pupper) {
-            Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>
-              circ(rgraph, zero_arc_map, *_pupper, neg_csup);
-            circ_result = circ.run();
-          } else {
-            Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>
-              circ(rgraph, zero_arc_map, inf_arc_map, neg_csup);
-            circ_result = circ.run();
-          }
-        }
-      }
-      if (local_csup) delete csup;
-      if (!circ_result) return false;
-
+          if (_cost[i] > ART_COST) ART_COST = _cost[i];
+        }
+        ART_COST = (ART_COST + 1) * _node_num;
+      }
+
+      // Initialize arc maps
+      for (int i = 0; i != _arc_num; ++i) {
+        _flow[i] = 0;
+        _state[i] = STATE_LOWER;
+      }
+      
       // Set data for the artificial root node
       _root = _node_num;
@@ -1220,87 +1091,28 @@
       _thread[_root] = 0;
       _rev_thread[0] = _root;
-      _succ_num[_root] = all_node_num;
+      _succ_num[_root] = _node_num + 1;
       _last_succ[_root] = _root - 1;
-      _supply[_root] = -sum_supply;
-      if (sum_supply < 0) {
-        _pi[_root] = -art_cost;
-      } else {
-        _pi[_root] = art_cost;
-      }
-
-      // Store the arcs in a mixed order
-      int k = std::max(int(std::sqrt(double(_arc_num))), 10);
-      int i = 0;
-      for (ArcIt e(_graph); e != INVALID; ++e) {
-        _arc_ref[i] = e;
-        if ((i += k) >= _arc_num) i = (i % k) + 1;
-      }
-
-      // Initialize arc maps
-      if (_pupper && _pcost) {
-        for (int i = 0; i != _arc_num; ++i) {
-          Arc e = _arc_ref[i];
-          _source[i] = _node_id[_graph.source(e)];
-          _target[i] = _node_id[_graph.target(e)];
-          _cap[i] = (*_pupper)[e];
-          _cost[i] = (*_pcost)[e];
-          _flow[i] = 0;
-          _state[i] = STATE_LOWER;
-        }
-      } else {
-        for (int i = 0; i != _arc_num; ++i) {
-          Arc e = _arc_ref[i];
-          _source[i] = _node_id[_graph.source(e)];
-          _target[i] = _node_id[_graph.target(e)];
-          _flow[i] = 0;
-          _state[i] = STATE_LOWER;
-        }
-        if (_pupper) {
-          for (int i = 0; i != _arc_num; ++i)
-            _cap[i] = (*_pupper)[_arc_ref[i]];
-        } else {
-          for (int i = 0; i != _arc_num; ++i)
-            _cap[i] = inf_cap;
-        }
-        if (_pcost) {
-          for (int i = 0; i != _arc_num; ++i)
-            _cost[i] = (*_pcost)[_arc_ref[i]];
-        } else {
-          for (int i = 0; i != _arc_num; ++i)
-            _cost[i] = 1;
-        }
-      }
-      
-      // Remove non-zero lower bounds
-      if (_plower) {
-        for (int i = 0; i != _arc_num; ++i) {
-          Flow c = (*_plower)[_arc_ref[i]];
-          if (c != 0) {
-            _cap[i] -= c;
-            _supply[_source[i]] -= c;
-            _supply[_target[i]] += c;
-          }
-        }
-      }
+      _supply[_root] = -_sum_supply;
+      _pi[_root] = _sum_supply < 0 ? -ART_COST : ART_COST;
 
       // Add artificial arcs and initialize the spanning tree data structure
       for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
+        _parent[u] = _root;
+        _pred[u] = e;
         _thread[u] = u + 1;
         _rev_thread[u + 1] = u;
         _succ_num[u] = 1;
         _last_succ[u] = u;
-        _parent[u] = _root;
-        _pred[u] = e;
-        _cost[e] = art_cost;
-        _cap[e] = inf_cap;
+        _cost[e] = ART_COST;
+        _cap[e] = INF;
         _state[e] = STATE_TREE;
-        if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
+        if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {
           _flow[e] = _supply[u];
           _forward[u] = true;
-          _pi[u] = -art_cost + _pi[_root];
+          _pi[u] = -ART_COST + _pi[_root];
         } else {
           _flow[e] = -_supply[u];
           _forward[u] = false;
-          _pi[u] = art_cost + _pi[_root];
+          _pi[u] = ART_COST + _pi[_root];
         }
       }
@@ -1337,5 +1149,5 @@
       delta = _cap[in_arc];
       int result = 0;
-      Flow d;
+      Value d;
       int e;
 
@@ -1343,5 +1155,6 @@
       for (int u = first; u != join; u = _parent[u]) {
         e = _pred[u];
-        d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
+        d = _forward[u] ?
+          _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]);
         if (d < delta) {
           delta = d;
@@ -1353,5 +1166,6 @@
       for (int u = second; u != join; u = _parent[u]) {
         e = _pred[u];
-        d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
+        d = _forward[u] ? 
+          (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e];
         if (d <= delta) {
           delta = d;
@@ -1375,5 +1189,5 @@
       // Augment along the cycle
       if (delta > 0) {
-        Flow val = _state[in_arc] * delta;
+        Value val = _state[in_arc] * delta;
         _flow[in_arc] += val;
         for (int u = _source[in_arc]; u != join; u = _parent[u]) {
@@ -1527,5 +1341,5 @@
 
     // Execute the algorithm
-    bool start(PivotRule pivot_rule) {
+    ProblemType start(PivotRule pivot_rule) {
       // Select the pivot rule implementation
       switch (pivot_rule) {
@@ -1541,9 +1355,9 @@
           return start<AlteringListPivotRule>();
       }
-      return false;
+      return INFEASIBLE; // avoid warning
     }
 
     template <typename PivotRuleImpl>
-    bool start() {
+    ProblemType start() {
       PivotRuleImpl pivot(*this);
 
@@ -1552,4 +1366,5 @@
         findJoinNode();
         bool change = findLeavingArc();
+        if (delta >= INF) return UNBOUNDED;
         changeFlow(change);
         if (change) {
@@ -1558,22 +1373,35 @@
         }
       }
-
-      // Copy flow values to _flow_map
-      if (_plower) {
+      
+      // Check feasibility
+      if (_sum_supply < 0) {
+        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
+          if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE;
+        }
+      }
+      else if (_sum_supply > 0) {
+        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
+          if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE;
+        }
+      }
+      else {
+        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
+          if (_flow[e] != 0) return INFEASIBLE;
+        }
+      }
+
+      // Transform the solution and the supply map to the original form
+      if (_have_lower) {
         for (int i = 0; i != _arc_num; ++i) {
-          Arc e = _arc_ref[i];
-          _flow_map->set(e, (*_plower)[e] + _flow[i]);
-        }
-      } else {
-        for (int i = 0; i != _arc_num; ++i) {
-          _flow_map->set(_arc_ref[i], _flow[i]);
-        }
-      }
-      // Copy potential values to _potential_map
-      for (NodeIt n(_graph); n != INVALID; ++n) {
-        _potential_map->set(n, _pi[_node_id[n]]);
-      }
-
-      return true;
+          Value c = _lower[i];
+          if (c != 0) {
+            _flow[i] += c;
+            _supply[_source[i]] += c;
+            _supply[_target[i]] -= c;
+          }
+        }
+      }
+
+      return OPTIMAL;
     }
 
Index: lemon/preflow.h
===================================================================
--- lemon/preflow.h	(revision 658)
+++ lemon/preflow.h	(revision 688)
@@ -47,5 +47,5 @@
 
     /// \brief The type of the flow values.
-    typedef typename CapacityMap::Value Flow;
+    typedef typename CapacityMap::Value Value;
 
     /// \brief The type of the map that stores the flow values.
@@ -53,5 +53,5 @@
     /// The type of the map that stores the flow values.
     /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
-    typedef typename Digraph::template ArcMap<Flow> FlowMap;
+    typedef typename Digraph::template ArcMap<Value> FlowMap;
 
     /// \brief Instantiates a FlowMap.
@@ -85,5 +85,5 @@
     ///
     /// The tolerance used by the algorithm to handle inexact computation.
-    typedef lemon::Tolerance<Flow> Tolerance;
+    typedef lemon::Tolerance<Value> Tolerance;
 
   };
@@ -126,5 +126,5 @@
     typedef typename Traits::CapacityMap CapacityMap;
     ///The type of the flow values.
-    typedef typename Traits::Flow Flow;
+    typedef typename Traits::Value Value;
 
     ///The type of the flow map.
@@ -152,5 +152,5 @@
     bool _local_level;
 
-    typedef typename Digraph::template NodeMap<Flow> ExcessMap;
+    typedef typename Digraph::template NodeMap<Value> ExcessMap;
     ExcessMap* _excess;
 
@@ -471,5 +471,5 @@
 
       for (NodeIt n(_graph); n != INVALID; ++n) {
-        Flow excess = 0;
+        Value excess = 0;
         for (InArcIt e(_graph, n); e != INVALID; ++e) {
           excess += (*_flow)[e];
@@ -520,5 +520,5 @@
 
       for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
-        Flow rem = (*_capacity)[e] - (*_flow)[e];
+        Value rem = (*_capacity)[e] - (*_flow)[e];
         if (_tolerance.positive(rem)) {
           Node u = _graph.target(e);
@@ -532,5 +532,5 @@
       }
       for (InArcIt e(_graph, _source); e != INVALID; ++e) {
-        Flow rem = (*_flow)[e];
+        Value rem = (*_flow)[e];
         if (_tolerance.positive(rem)) {
           Node v = _graph.source(e);
@@ -565,9 +565,9 @@
 
         while (num > 0 && n != INVALID) {
-          Flow excess = (*_excess)[n];
+          Value excess = (*_excess)[n];
           int new_level = _level->maxLevel();
 
           for (OutArcIt e(_graph, n); e != INVALID; ++e) {
-            Flow rem = (*_capacity)[e] - (*_flow)[e];
+            Value rem = (*_capacity)[e] - (*_flow)[e];
             if (!_tolerance.positive(rem)) continue;
             Node v = _graph.target(e);
@@ -592,5 +592,5 @@
 
           for (InArcIt e(_graph, n); e != INVALID; ++e) {
-            Flow rem = (*_flow)[e];
+            Value rem = (*_flow)[e];
             if (!_tolerance.positive(rem)) continue;
             Node v = _graph.source(e);
@@ -638,9 +638,9 @@
         num = _node_num * 20;
         while (num > 0 && n != INVALID) {
-          Flow excess = (*_excess)[n];
+          Value excess = (*_excess)[n];
           int new_level = _level->maxLevel();
 
           for (OutArcIt e(_graph, n); e != INVALID; ++e) {
-            Flow rem = (*_capacity)[e] - (*_flow)[e];
+            Value rem = (*_capacity)[e] - (*_flow)[e];
             if (!_tolerance.positive(rem)) continue;
             Node v = _graph.target(e);
@@ -665,5 +665,5 @@
 
           for (InArcIt e(_graph, n); e != INVALID; ++e) {
-            Flow rem = (*_flow)[e];
+            Value rem = (*_flow)[e];
             if (!_tolerance.positive(rem)) continue;
             Node v = _graph.source(e);
@@ -779,10 +779,10 @@
       Node n;
       while ((n = _level->highestActive()) != INVALID) {
-        Flow excess = (*_excess)[n];
+        Value excess = (*_excess)[n];
         int level = _level->highestActiveLevel();
         int new_level = _level->maxLevel();
 
         for (OutArcIt e(_graph, n); e != INVALID; ++e) {
-          Flow rem = (*_capacity)[e] - (*_flow)[e];
+          Value rem = (*_capacity)[e] - (*_flow)[e];
           if (!_tolerance.positive(rem)) continue;
           Node v = _graph.target(e);
@@ -807,5 +807,5 @@
 
         for (InArcIt e(_graph, n); e != INVALID; ++e) {
-          Flow rem = (*_flow)[e];
+          Value rem = (*_flow)[e];
           if (!_tolerance.positive(rem)) continue;
           Node v = _graph.source(e);
@@ -898,16 +898,16 @@
     /// \pre Either \ref run() or \ref init() must be called before
     /// using this function.
-    Flow flowValue() const {
+    Value flowValue() const {
       return (*_excess)[_target];
     }
 
-    /// \brief Returns the flow on the given arc.
-    ///
-    /// Returns the flow on the given arc. This method can
+    /// \brief Returns the flow value on the given arc.
+    ///
+    /// Returns the flow value on the given arc. This method can
     /// be called after the second phase of the algorithm.
     ///
     /// \pre Either \ref run() or \ref init() must be called before
     /// using this function.
-    Flow flow(const Arc& arc) const {
+    Value flow(const Arc& arc) const {
       return (*_flow)[arc];
     }
