Index: doc/groups.dox
===================================================================
--- doc/groups.dox	(revision 818)
+++ doc/groups.dox	(revision 879)
@@ -407,13 +407,11 @@
  - \ref NetworkSimplex Primal Network Simplex algorithm with various
    pivot strategies \ref dantzig63linearprog, \ref kellyoneill91netsimplex.
- - \ref CostScaling Push-Relabel and Augment-Relabel algorithms based on
-   cost scaling \ref goldberg90approximation, \ref goldberg97efficient,
+ - \ref CostScaling Cost Scaling algorithm based on push/augment and
+   relabel operations \ref goldberg90approximation, \ref goldberg97efficient,
    \ref bunnagel98efficient.
- - \ref CapacityScaling Successive Shortest %Path algorithm with optional
-   capacity scaling \ref edmondskarp72theoretical.
- - \ref CancelAndTighten The Cancel and Tighten algorithm
-   \ref goldberg89cyclecanceling.
- - \ref CycleCanceling Cycle-Canceling algorithms
-   \ref klein67primal, \ref goldberg89cyclecanceling.
+ - \ref CapacityScaling Capacity Scaling algorithm based on the successive
+   shortest path method \ref edmondskarp72theoretical.
+ - \ref CycleCanceling Cycle-Canceling algorithms, two of which are
+   strongly polynomial \ref klein67primal, \ref goldberg89cyclecanceling.
 
 In general NetworkSimplex is the most efficient implementation,
Index: lemon/Makefile.am
===================================================================
--- lemon/Makefile.am	(revision 863)
+++ lemon/Makefile.am	(revision 888)
@@ -63,4 +63,5 @@
 	lemon/binom_heap.h \
 	lemon/bucket_heap.h \
+	lemon/capacity_scaling.h \
 	lemon/cbc.h \
 	lemon/circulation.h \
@@ -69,7 +70,9 @@
 	lemon/concept_check.h \
 	lemon/connectivity.h \
+	lemon/core.h \
+	lemon/cost_scaling.h \
 	lemon/counter.h \
-	lemon/core.h \
 	lemon/cplex.h \
+	lemon/cycle_canceling.h \
 	lemon/dfs.h \
 	lemon/dijkstra.h \
Index: lemon/bellman_ford.h
===================================================================
--- lemon/bellman_ford.h	(revision 835)
+++ lemon/bellman_ford.h	(revision 870)
@@ -238,5 +238,7 @@
 	_dist = Traits::createDistMap(*_gr);
       }
-      _mask = new MaskMap(*_gr, false);
+      if(!_mask) {
+        _mask = new MaskMap(*_gr);
+      }
     }
     
@@ -404,4 +406,8 @@
 	  _process.push_back(it);
 	  _mask->set(it, true);
+	}
+      } else {
+	for (NodeIt it(*_gr); it != INVALID; ++it) {
+	  _mask->set(it, false);
 	}
       }
Index: lemon/capacity_scaling.h
===================================================================
--- lemon/capacity_scaling.h	(revision 887)
+++ lemon/capacity_scaling.h	(revision 887)
@@ -0,0 +1,951 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_CAPACITY_SCALING_H
+#define LEMON_CAPACITY_SCALING_H
+
+/// \ingroup min_cost_flow_algs
+///
+/// \file
+/// \brief Capacity Scaling algorithm for finding a minimum cost flow.
+
+#include <vector>
+#include <limits>
+#include <lemon/core.h>
+#include <lemon/bin_heap.h>
+
+namespace lemon {
+
+  /// \brief Default traits class of CapacityScaling algorithm.
+  ///
+  /// Default traits class of CapacityScaling algorithm.
+  /// \tparam GR Digraph type.
+  /// \tparam V The number type used for flow amounts, capacity bounds
+  /// and supply values. By default it is \c int.
+  /// \tparam C The number type used for costs and potentials.
+  /// By default it is the same as \c V.
+  template <typename GR, typename V = int, typename C = V>
+  struct CapacityScalingDefaultTraits
+  {
+    /// The type of the digraph
+    typedef GR Digraph;
+    /// The type of the flow amounts, capacity bounds and supply values
+    typedef V Value;
+    /// The type of the arc costs
+    typedef C Cost;
+
+    /// \brief The type of the heap used for internal Dijkstra computations.
+    ///
+    /// The type of the heap used for internal Dijkstra computations.
+    /// It must conform to the \ref lemon::concepts::Heap "Heap" concept,
+    /// its priority type must be \c Cost and its cross reference type
+    /// must be \ref RangeMap "RangeMap<int>".
+    typedef BinHeap<Cost, RangeMap<int> > Heap;
+  };
+
+  /// \addtogroup min_cost_flow_algs
+  /// @{
+
+  /// \brief Implementation of the Capacity Scaling algorithm for
+  /// finding a \ref min_cost_flow "minimum cost flow".
+  ///
+  /// \ref CapacityScaling implements the capacity scaling version
+  /// of the successive shortest path algorithm for finding a
+  /// \ref min_cost_flow "minimum cost flow" \ref amo93networkflows,
+  /// \ref edmondskarp72theoretical. It is an efficient dual
+  /// solution method.
+  ///
+  /// 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 V The number type used for flow amounts, capacity bounds
+  /// and supply values in the algorithm. By default it is \c int.
+  /// \tparam C The number type used for costs and potentials in the
+  /// algorithm. By default it is the same as \c V.
+  ///
+  /// \warning Both number types must be signed and all input data must
+  /// be integer.
+  /// \warning This algorithm does not support negative costs for such
+  /// arcs that have infinite upper bound.
+#ifdef DOXYGEN
+  template <typename GR, typename V, typename C, typename TR>
+#else
+  template < typename GR, typename V = int, typename C = V,
+             typename TR = CapacityScalingDefaultTraits<GR, V, C> >
+#endif
+  class CapacityScaling
+  {
+  public:
+
+    /// The type of the digraph
+    typedef typename TR::Digraph Digraph;
+    /// The type of the flow amounts, capacity bounds and supply values
+    typedef typename TR::Value Value;
+    /// The type of the arc costs
+    typedef typename TR::Cost Cost;
+
+    /// The type of the heap used for internal Dijkstra computations
+    typedef typename TR::Heap Heap;
+
+    /// The \ref CapacityScalingDefaultTraits "traits class" of the algorithm
+    typedef TR Traits;
+
+  public:
+
+    /// \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 digraph contains an arc of negative cost and infinite
+      /// upper bound. It means that the objective function is unbounded
+      /// on that arc, however, note that it could actually be bounded
+      /// over the feasible flows, but this algroithm cannot handle
+      /// these cases.
+      UNBOUNDED
+    };
+  
+  private:
+
+    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
+
+    typedef std::vector<int> IntVector;
+    typedef std::vector<char> BoolVector;
+    typedef std::vector<Value> ValueVector;
+    typedef std::vector<Cost> CostVector;
+
+  private:
+
+    // Data related to the underlying digraph
+    const GR &_graph;
+    int _node_num;
+    int _arc_num;
+    int _res_arc_num;
+    int _root;
+
+    // Parameters of the problem
+    bool _have_lower;
+    Value _sum_supply;
+
+    // Data structures for storing the digraph
+    IntNodeMap _node_id;
+    IntArcMap _arc_idf;
+    IntArcMap _arc_idb;
+    IntVector _first_out;
+    BoolVector _forward;
+    IntVector _source;
+    IntVector _target;
+    IntVector _reverse;
+
+    // Node and arc data
+    ValueVector _lower;
+    ValueVector _upper;
+    CostVector _cost;
+    ValueVector _supply;
+
+    ValueVector _res_cap;
+    CostVector _pi;
+    ValueVector _excess;
+    IntVector _excess_nodes;
+    IntVector _deficit_nodes;
+
+    Value _delta;
+    int _factor;
+    IntVector _pred;
+
+  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:
+
+    // Special implementation of the Dijkstra algorithm for finding
+    // shortest paths in the residual network of the digraph with
+    // respect to the reduced arc costs and modifying the node
+    // potentials according to the found distance labels.
+    class ResidualDijkstra
+    {
+    private:
+
+      int _node_num;
+      bool _geq;
+      const IntVector &_first_out;
+      const IntVector &_target;
+      const CostVector &_cost;
+      const ValueVector &_res_cap;
+      const ValueVector &_excess;
+      CostVector &_pi;
+      IntVector &_pred;
+      
+      IntVector _proc_nodes;
+      CostVector _dist;
+      
+    public:
+
+      ResidualDijkstra(CapacityScaling& cs) :
+        _node_num(cs._node_num), _geq(cs._sum_supply < 0),
+        _first_out(cs._first_out), _target(cs._target), _cost(cs._cost),
+        _res_cap(cs._res_cap), _excess(cs._excess), _pi(cs._pi),
+        _pred(cs._pred), _dist(cs._node_num)
+      {}
+
+      int run(int s, Value delta = 1) {
+        RangeMap<int> heap_cross_ref(_node_num, Heap::PRE_HEAP);
+        Heap heap(heap_cross_ref);
+        heap.push(s, 0);
+        _pred[s] = -1;
+        _proc_nodes.clear();
+
+        // Process nodes
+        while (!heap.empty() && _excess[heap.top()] > -delta) {
+          int u = heap.top(), v;
+          Cost d = heap.prio() + _pi[u], dn;
+          _dist[u] = heap.prio();
+          _proc_nodes.push_back(u);
+          heap.pop();
+
+          // Traverse outgoing residual arcs
+          int last_out = _geq ? _first_out[u+1] : _first_out[u+1] - 1;
+          for (int a = _first_out[u]; a != last_out; ++a) {
+            if (_res_cap[a] < delta) continue;
+            v = _target[a];
+            switch (heap.state(v)) {
+              case Heap::PRE_HEAP:
+                heap.push(v, d + _cost[a] - _pi[v]);
+                _pred[v] = a;
+                break;
+              case Heap::IN_HEAP:
+                dn = d + _cost[a] - _pi[v];
+                if (dn < heap[v]) {
+                  heap.decrease(v, dn);
+                  _pred[v] = a;
+                }
+                break;
+              case Heap::POST_HEAP:
+                break;
+            }
+          }
+        }
+        if (heap.empty()) return -1;
+
+        // Update potentials of processed nodes
+        int t = heap.top();
+        Cost dt = heap.prio();
+        for (int i = 0; i < int(_proc_nodes.size()); ++i) {
+          _pi[_proc_nodes[i]] += _dist[_proc_nodes[i]] - dt;
+        }
+
+        return t;
+      }
+
+    }; //class ResidualDijkstra
+
+  public:
+
+    /// \name Named Template Parameters
+    /// @{
+
+    template <typename T>
+    struct SetHeapTraits : public Traits {
+      typedef T Heap;
+    };
+
+    /// \brief \ref named-templ-param "Named parameter" for setting
+    /// \c Heap type.
+    ///
+    /// \ref named-templ-param "Named parameter" for setting \c Heap
+    /// type, which is used for internal Dijkstra computations.
+    /// It must conform to the \ref lemon::concepts::Heap "Heap" concept,
+    /// its priority type must be \c Cost and its cross reference type
+    /// must be \ref RangeMap "RangeMap<int>".
+    template <typename T>
+    struct SetHeap
+      : public CapacityScaling<GR, V, C, SetHeapTraits<T> > {
+      typedef  CapacityScaling<GR, V, C, SetHeapTraits<T> > Create;
+    };
+
+    /// @}
+
+  public:
+
+    /// \brief Constructor.
+    ///
+    /// The constructor of the class.
+    ///
+    /// \param graph The digraph the algorithm runs on.
+    CapacityScaling(const GR& graph) :
+      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
+      INF(std::numeric_limits<Value>::has_infinity ?
+          std::numeric_limits<Value>::infinity() :
+          std::numeric_limits<Value>::max())
+    {
+      // Check the number types
+      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
+        "The flow type of CapacityScaling must be signed");
+      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
+        "The cost type of CapacityScaling must be signed");
+
+      // Resize vectors
+      _node_num = countNodes(_graph);
+      _arc_num = countArcs(_graph);
+      _res_arc_num = 2 * (_arc_num + _node_num);
+      _root = _node_num;
+      ++_node_num;
+
+      _first_out.resize(_node_num + 1);
+      _forward.resize(_res_arc_num);
+      _source.resize(_res_arc_num);
+      _target.resize(_res_arc_num);
+      _reverse.resize(_res_arc_num);
+
+      _lower.resize(_res_arc_num);
+      _upper.resize(_res_arc_num);
+      _cost.resize(_res_arc_num);
+      _supply.resize(_node_num);
+      
+      _res_cap.resize(_res_arc_num);
+      _pi.resize(_node_num);
+      _excess.resize(_node_num);
+      _pred.resize(_node_num);
+
+      // Copy the graph
+      int i = 0, j = 0, k = 2 * _arc_num + _node_num - 1;
+      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
+        _node_id[n] = i;
+      }
+      i = 0;
+      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
+        _first_out[i] = j;
+        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
+          _arc_idf[a] = j;
+          _forward[j] = true;
+          _source[j] = i;
+          _target[j] = _node_id[_graph.runningNode(a)];
+        }
+        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
+          _arc_idb[a] = j;
+          _forward[j] = false;
+          _source[j] = i;
+          _target[j] = _node_id[_graph.runningNode(a)];
+        }
+        _forward[j] = false;
+        _source[j] = i;
+        _target[j] = _root;
+        _reverse[j] = k;
+        _forward[k] = true;
+        _source[k] = _root;
+        _target[k] = i;
+        _reverse[k] = j;
+        ++j; ++k;
+      }
+      _first_out[i] = j;
+      _first_out[_node_num] = k;
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        int fi = _arc_idf[a];
+        int bi = _arc_idb[a];
+        _reverse[fi] = bi;
+        _reverse[bi] = fi;
+      }
+      
+      // Reset parameters
+      reset();
+    }
+
+    /// \name Parameters
+    /// The parameters of the algorithm can be specified using these
+    /// functions.
+
+    /// @{
+
+    /// \brief Set the lower bounds on the arcs.
+    ///
+    /// This function sets the lower bounds on the 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 Value type
+    /// of the algorithm.
+    ///
+    /// \return <tt>(*this)</tt>
+    template <typename LowerMap>
+    CapacityScaling& lowerMap(const LowerMap& map) {
+      _have_lower = true;
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        _lower[_arc_idf[a]] = map[a];
+        _lower[_arc_idb[a]] = map[a];
+      }
+      return *this;
+    }
+
+    /// \brief Set the upper bounds (capacities) on the arcs.
+    ///
+    /// This function sets the upper bounds (capacities) on the 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).
+    ///
+    /// \param map An arc map storing the upper bounds.
+    /// Its \c Value type must be convertible to the \c Value type
+    /// of the algorithm.
+    ///
+    /// \return <tt>(*this)</tt>
+    template<typename UpperMap>
+    CapacityScaling& upperMap(const UpperMap& map) {
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        _upper[_arc_idf[a]] = map[a];
+      }
+      return *this;
+    }
+
+    /// \brief Set the costs of the arcs.
+    ///
+    /// This function sets the costs of the arcs.
+    /// If it is not used before calling \ref run(), the costs
+    /// will be set to \c 1 on all arcs.
+    ///
+    /// \param map An arc map storing the costs.
+    /// Its \c Value type must be convertible to the \c Cost type
+    /// of the algorithm.
+    ///
+    /// \return <tt>(*this)</tt>
+    template<typename CostMap>
+    CapacityScaling& costMap(const CostMap& map) {
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        _cost[_arc_idf[a]] =  map[a];
+        _cost[_arc_idb[a]] = -map[a];
+      }
+      return *this;
+    }
+
+    /// \brief Set the supply values of the nodes.
+    ///
+    /// This function sets the supply values of the nodes.
+    /// If neither this function nor \ref stSupply() is used before
+    /// calling \ref run(), the supply of each node will be set to zero.
+    ///
+    /// \param map A node map storing the supply values.
+    /// Its \c Value type must be convertible to the \c Value type
+    /// of the algorithm.
+    ///
+    /// \return <tt>(*this)</tt>
+    template<typename SupplyMap>
+    CapacityScaling& supplyMap(const SupplyMap& map) {
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        _supply[_node_id[n]] = map[n];
+      }
+      return *this;
+    }
+
+    /// \brief Set single source and target nodes and a supply value.
+    ///
+    /// This function sets a single source node and a single target node
+    /// and the required flow value.
+    /// If neither this function nor \ref supplyMap() is used before
+    /// calling \ref run(), the supply of each node will be set to zero.
+    ///
+    /// 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.
+    /// \param k The required amount of flow from node \c s to node \c t
+    /// (i.e. the supply of \c s and the demand of \c t).
+    ///
+    /// \return <tt>(*this)</tt>
+    CapacityScaling& 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;
+    }
+    
+    /// @}
+
+    /// \name Execution control
+    /// The algorithm can be executed using \ref run().
+
+    /// @{
+
+    /// \brief Run the algorithm.
+    ///
+    /// This function runs the algorithm.
+    /// The paramters can be specified using functions \ref lowerMap(),
+    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
+    /// For example,
+    /// \code
+    ///   CapacityScaling<ListDigraph> cs(graph);
+    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
+    ///     .supplyMap(sup).run();
+    /// \endcode
+    ///
+    /// This function can be called more than once. All the parameters
+    /// that have been given are kept for the next call, unless
+    /// \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 factor The capacity scaling factor. It must be larger than
+    /// one to use scaling. If it is less or equal to one, then scaling
+    /// will be disabled.
+    ///
+    /// \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 digraph contains an arc of negative cost
+    /// and infinite upper bound. It means that the objective function
+    /// is unbounded on that arc, however, note that it could actually be
+    /// bounded over the feasible flows, but this algroithm cannot handle
+    /// these cases.
+    ///
+    /// \see ProblemType
+    ProblemType run(int factor = 4) {
+      _factor = factor;
+      ProblemType pt = init();
+      if (pt != OPTIMAL) return pt;
+      return start();
+    }
+
+    /// \brief Reset all the parameters that have been given before.
+    ///
+    /// This function resets all the paramaters that have been given
+    /// before using functions \ref lowerMap(), \ref upperMap(),
+    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
+    ///
+    /// 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,
+    /// \code
+    ///   CapacityScaling<ListDigraph> cs(graph);
+    ///
+    ///   // First run
+    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
+    ///     .supplyMap(sup).run();
+    ///
+    ///   // Run again with modified cost map (reset() is not called,
+    ///   // so only the cost map have to be set again)
+    ///   cost[e] += 100;
+    ///   cs.costMap(cost).run();
+    ///
+    ///   // Run again from scratch using reset()
+    ///   // (the lower bounds will be set to zero on all arcs)
+    ///   cs.reset();
+    ///   cs.upperMap(capacity).costMap(cost)
+    ///     .supplyMap(sup).run();
+    /// \endcode
+    ///
+    /// \return <tt>(*this)</tt>
+    CapacityScaling& reset() {
+      for (int i = 0; i != _node_num; ++i) {
+        _supply[i] = 0;
+      }
+      for (int j = 0; j != _res_arc_num; ++j) {
+        _lower[j] = 0;
+        _upper[j] = INF;
+        _cost[j] = _forward[j] ? 1 : -1;
+      }
+      _have_lower = false;
+      return *this;
+    }
+
+    /// @}
+
+    /// \name Query Functions
+    /// The results of the algorithm can be obtained using these
+    /// functions.\n
+    /// The \ref run() function must be called before using them.
+
+    /// @{
+
+    /// \brief Return the total cost of the found flow.
+    ///
+    /// This function returns the total cost of the found flow.
+    /// Its complexity is O(e).
+    ///
+    /// \note The return type of the function can be specified as a
+    /// template parameter. For example,
+    /// \code
+    ///   cs.totalCost<double>();
+    /// \endcode
+    /// It is useful if the total cost cannot be stored in the \c Cost
+    /// type of the algorithm, which is the default return type of the
+    /// function.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    template <typename Number>
+    Number totalCost() const {
+      Number c = 0;
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        int i = _arc_idb[a];
+        c += static_cast<Number>(_res_cap[i]) *
+             (-static_cast<Number>(_cost[i]));
+      }
+      return c;
+    }
+
+#ifndef DOXYGEN
+    Cost totalCost() const {
+      return totalCost<Cost>();
+    }
+#endif
+
+    /// \brief Return the flow on the given arc.
+    ///
+    /// This function returns the flow on the given arc.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    Value flow(const Arc& a) const {
+      return _res_cap[_arc_idb[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.
+    template <typename FlowMap>
+    void flowMap(FlowMap &map) const {
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        map.set(a, _res_cap[_arc_idb[a]]);
+      }
+    }
+
+    /// \brief Return the potential (dual value) of the given node.
+    ///
+    /// This function returns the potential (dual value) of the
+    /// given node.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    Cost potential(const Node& n) const {
+      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.
+    template <typename PotentialMap>
+    void potentialMap(PotentialMap &map) const {
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        map.set(n, _pi[_node_id[n]]);
+      }
+    }
+
+    /// @}
+
+  private:
+
+    // Initialize the algorithm
+    ProblemType init() {
+      if (_node_num <= 1) return INFEASIBLE;
+
+      // Check the sum of supply values
+      _sum_supply = 0;
+      for (int i = 0; i != _root; ++i) {
+        _sum_supply += _supply[i];
+      }
+      if (_sum_supply > 0) return INFEASIBLE;
+      
+      // Initialize vectors
+      for (int i = 0; i != _root; ++i) {
+        _pi[i] = 0;
+        _excess[i] = _supply[i];
+      }
+
+      // Remove non-zero lower bounds
+      const Value MAX = std::numeric_limits<Value>::max();
+      int last_out;
+      if (_have_lower) {
+        for (int i = 0; i != _root; ++i) {
+          last_out = _first_out[i+1];
+          for (int j = _first_out[i]; j != last_out; ++j) {
+            if (_forward[j]) {
+              Value c = _lower[j];
+              if (c >= 0) {
+                _res_cap[j] = _upper[j] < MAX ? _upper[j] - c : INF;
+              } else {
+                _res_cap[j] = _upper[j] < MAX + c ? _upper[j] - c : INF;
+              }
+              _excess[i] -= c;
+              _excess[_target[j]] += c;
+            } else {
+              _res_cap[j] = 0;
+            }
+          }
+        }
+      } else {
+        for (int j = 0; j != _res_arc_num; ++j) {
+          _res_cap[j] = _forward[j] ? _upper[j] : 0;
+        }
+      }
+
+      // Handle negative costs
+      for (int i = 0; i != _root; ++i) {
+        last_out = _first_out[i+1] - 1;
+        for (int j = _first_out[i]; j != last_out; ++j) {
+          Value rc = _res_cap[j];
+          if (_cost[j] < 0 && rc > 0) {
+            if (rc >= MAX) return UNBOUNDED;
+            _excess[i] -= rc;
+            _excess[_target[j]] += rc;
+            _res_cap[j] = 0;
+            _res_cap[_reverse[j]] += rc;
+          }
+        }
+      }
+      
+      // Handle GEQ supply type
+      if (_sum_supply < 0) {
+        _pi[_root] = 0;
+        _excess[_root] = -_sum_supply;
+        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
+          int ra = _reverse[a];
+          _res_cap[a] = -_sum_supply + 1;
+          _res_cap[ra] = 0;
+          _cost[a] = 0;
+          _cost[ra] = 0;
+        }
+      } else {
+        _pi[_root] = 0;
+        _excess[_root] = 0;
+        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
+          int ra = _reverse[a];
+          _res_cap[a] = 1;
+          _res_cap[ra] = 0;
+          _cost[a] = 0;
+          _cost[ra] = 0;
+        }
+      }
+
+      // Initialize delta value
+      if (_factor > 1) {
+        // With scaling
+        Value max_sup = 0, max_dem = 0;
+        for (int i = 0; i != _node_num; ++i) {
+          Value ex = _excess[i];
+          if ( ex > max_sup) max_sup =  ex;
+          if (-ex > max_dem) max_dem = -ex;
+        }
+        Value max_cap = 0;
+        for (int j = 0; j != _res_arc_num; ++j) {
+          if (_res_cap[j] > max_cap) max_cap = _res_cap[j];
+        }
+        max_sup = std::min(std::min(max_sup, max_dem), max_cap);
+        for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2) ;
+      } else {
+        // Without scaling
+        _delta = 1;
+      }
+
+      return OPTIMAL;
+    }
+
+    ProblemType start() {
+      // Execute the algorithm
+      ProblemType pt;
+      if (_delta > 1)
+        pt = startWithScaling();
+      else
+        pt = startWithoutScaling();
+
+      // Handle non-zero lower bounds
+      if (_have_lower) {
+        int limit = _first_out[_root];
+        for (int j = 0; j != limit; ++j) {
+          if (!_forward[j]) _res_cap[j] += _lower[j];
+        }
+      }
+
+      // Shift potentials if necessary
+      Cost pr = _pi[_root];
+      if (_sum_supply < 0 || pr > 0) {
+        for (int i = 0; i != _node_num; ++i) {
+          _pi[i] -= pr;
+        }        
+      }
+      
+      return pt;
+    }
+
+    // Execute the capacity scaling algorithm
+    ProblemType startWithScaling() {
+      // Perform capacity scaling phases
+      int s, t;
+      ResidualDijkstra _dijkstra(*this);
+      while (true) {
+        // Saturate all arcs not satisfying the optimality condition
+        int last_out;
+        for (int u = 0; u != _node_num; ++u) {
+          last_out = _sum_supply < 0 ?
+            _first_out[u+1] : _first_out[u+1] - 1;
+          for (int a = _first_out[u]; a != last_out; ++a) {
+            int v = _target[a];
+            Cost c = _cost[a] + _pi[u] - _pi[v];
+            Value rc = _res_cap[a];
+            if (c < 0 && rc >= _delta) {
+              _excess[u] -= rc;
+              _excess[v] += rc;
+              _res_cap[a] = 0;
+              _res_cap[_reverse[a]] += rc;
+            }
+          }
+        }
+
+        // Find excess nodes and deficit nodes
+        _excess_nodes.clear();
+        _deficit_nodes.clear();
+        for (int u = 0; u != _node_num; ++u) {
+          Value ex = _excess[u];
+          if (ex >=  _delta) _excess_nodes.push_back(u);
+          if (ex <= -_delta) _deficit_nodes.push_back(u);
+        }
+        int next_node = 0, next_def_node = 0;
+
+        // Find augmenting shortest paths
+        while (next_node < int(_excess_nodes.size())) {
+          // Check deficit nodes
+          if (_delta > 1) {
+            bool delta_deficit = false;
+            for ( ; next_def_node < int(_deficit_nodes.size());
+                    ++next_def_node ) {
+              if (_excess[_deficit_nodes[next_def_node]] <= -_delta) {
+                delta_deficit = true;
+                break;
+              }
+            }
+            if (!delta_deficit) break;
+          }
+
+          // Run Dijkstra in the residual network
+          s = _excess_nodes[next_node];
+          if ((t = _dijkstra.run(s, _delta)) == -1) {
+            if (_delta > 1) {
+              ++next_node;
+              continue;
+            }
+            return INFEASIBLE;
+          }
+
+          // Augment along a shortest path from s to t
+          Value d = std::min(_excess[s], -_excess[t]);
+          int u = t;
+          int a;
+          if (d > _delta) {
+            while ((a = _pred[u]) != -1) {
+              if (_res_cap[a] < d) d = _res_cap[a];
+              u = _source[a];
+            }
+          }
+          u = t;
+          while ((a = _pred[u]) != -1) {
+            _res_cap[a] -= d;
+            _res_cap[_reverse[a]] += d;
+            u = _source[a];
+          }
+          _excess[s] -= d;
+          _excess[t] += d;
+
+          if (_excess[s] < _delta) ++next_node;
+        }
+
+        if (_delta == 1) break;
+        _delta = _delta <= _factor ? 1 : _delta / _factor;
+      }
+
+      return OPTIMAL;
+    }
+
+    // Execute the successive shortest path algorithm
+    ProblemType startWithoutScaling() {
+      // Find excess nodes
+      _excess_nodes.clear();
+      for (int i = 0; i != _node_num; ++i) {
+        if (_excess[i] > 0) _excess_nodes.push_back(i);
+      }
+      if (_excess_nodes.size() == 0) return OPTIMAL;
+      int next_node = 0;
+
+      // Find shortest paths
+      int s, t;
+      ResidualDijkstra _dijkstra(*this);
+      while ( _excess[_excess_nodes[next_node]] > 0 ||
+              ++next_node < int(_excess_nodes.size()) )
+      {
+        // Run Dijkstra in the residual network
+        s = _excess_nodes[next_node];
+        if ((t = _dijkstra.run(s)) == -1) return INFEASIBLE;
+
+        // Augment along a shortest path from s to t
+        Value d = std::min(_excess[s], -_excess[t]);
+        int u = t;
+        int a;
+        if (d > 1) {
+          while ((a = _pred[u]) != -1) {
+            if (_res_cap[a] < d) d = _res_cap[a];
+            u = _source[a];
+          }
+        }
+        u = t;
+        while ((a = _pred[u]) != -1) {
+          _res_cap[a] -= d;
+          _res_cap[_reverse[a]] += d;
+          u = _source[a];
+        }
+        _excess[s] -= d;
+        _excess[t] += d;
+      }
+
+      return OPTIMAL;
+    }
+
+  }; //class CapacityScaling
+
+  ///@}
+
+} //namespace lemon
+
+#endif //LEMON_CAPACITY_SCALING_H
Index: lemon/concepts/heap.h
===================================================================
--- lemon/concepts/heap.h	(revision 757)
+++ lemon/concepts/heap.h	(revision 883)
@@ -89,5 +89,9 @@
       /// handle the cross references. The assigned value must be
       /// \c PRE_HEAP (<tt>-1</tt>) for each item.
+#ifdef DOXYGEN
       explicit Heap(ItemIntMap &map) {}
+#else
+      explicit Heap(ItemIntMap&) {}
+#endif      
 
       /// \brief Constructor.
@@ -99,5 +103,9 @@
       /// \c PRE_HEAP (<tt>-1</tt>) for each item.
       /// \param comp The function object used for comparing the priorities.
+#ifdef DOXYGEN
       explicit Heap(ItemIntMap &map, const CMP &comp) {}
+#else
+      explicit Heap(ItemIntMap&, const CMP&) {}
+#endif      
 
       /// \brief The number of items stored in the heap.
@@ -127,5 +135,9 @@
       /// \param p The priority of the item.
       /// \pre \e i must not be stored in the heap.
+#ifdef DOXYGEN
       void push(const Item &i, const Prio &p) {}
+#else
+      void push(const Item&, const Prio&) {}
+#endif      
 
       /// \brief Return the item having minimum priority.
@@ -133,5 +145,5 @@
       /// This function returns the item having minimum priority.
       /// \pre The heap must be non-empty.
-      Item top() const {}
+      Item top() const { return Item(); }
 
       /// \brief The minimum priority.
@@ -139,5 +151,5 @@
       /// This function returns the minimum priority.
       /// \pre The heap must be non-empty.
-      Prio prio() const {}
+      Prio prio() const { return Prio(); }
 
       /// \brief Remove the item having minimum priority.
@@ -153,5 +165,9 @@
       /// \param i The item to delete.
       /// \pre \e i must be in the heap.
+#ifdef DOXYGEN
       void erase(const Item &i) {}
+#else
+      void erase(const Item&) {}
+#endif      
 
       /// \brief The priority of the given item.
@@ -160,5 +176,9 @@
       /// \param i The item.
       /// \pre \e i must be in the heap.
+#ifdef DOXYGEN
       Prio operator[](const Item &i) const {}
+#else
+      Prio operator[](const Item&) const { return Prio(); }
+#endif      
 
       /// \brief Set the priority of an item or insert it, if it is
@@ -171,5 +191,9 @@
       /// \param i The item.
       /// \param p The priority.
+#ifdef DOXYGEN
       void set(const Item &i, const Prio &p) {}
+#else
+      void set(const Item&, const Prio&) {}
+#endif      
 
       /// \brief Decrease the priority of an item to the given value.
@@ -179,5 +203,9 @@
       /// \param p The priority.
       /// \pre \e i must be stored in the heap with priority at least \e p.
+#ifdef DOXYGEN
       void decrease(const Item &i, const Prio &p) {}
+#else
+      void decrease(const Item&, const Prio&) {}
+#endif      
 
       /// \brief Increase the priority of an item to the given value.
@@ -187,5 +215,9 @@
       /// \param p The priority.
       /// \pre \e i must be stored in the heap with priority at most \e p.
+#ifdef DOXYGEN
       void increase(const Item &i, const Prio &p) {}
+#else
+      void increase(const Item&, const Prio&) {}
+#endif      
 
       /// \brief Return the state of an item.
@@ -197,5 +229,9 @@
       /// to the heap again.
       /// \param i The item.
+#ifdef DOXYGEN
       State state(const Item &i) const {}
+#else
+      State state(const Item&) const { return PRE_HEAP; }
+#endif      
 
       /// \brief Set the state of an item in the heap.
@@ -206,5 +242,9 @@
       /// \param i The item.
       /// \param st The state. It should not be \c IN_HEAP.
+#ifdef DOXYGEN
       void state(const Item& i, State st) {}
+#else
+      void state(const Item&, State) {}
+#endif      
 
 
Index: lemon/cost_scaling.h
===================================================================
--- lemon/cost_scaling.h	(revision 887)
+++ lemon/cost_scaling.h	(revision 887)
@@ -0,0 +1,1170 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_COST_SCALING_H
+#define LEMON_COST_SCALING_H
+
+/// \ingroup min_cost_flow_algs
+/// \file
+/// \brief Cost scaling algorithm for finding a minimum cost flow.
+
+#include <vector>
+#include <deque>
+#include <limits>
+
+#include <lemon/core.h>
+#include <lemon/maps.h>
+#include <lemon/math.h>
+#include <lemon/static_graph.h>
+#include <lemon/circulation.h>
+#include <lemon/bellman_ford.h>
+
+namespace lemon {
+
+  /// \brief Default traits class of CostScaling algorithm.
+  ///
+  /// Default traits class of CostScaling algorithm.
+  /// \tparam GR Digraph type.
+  /// \tparam V The number type used for flow amounts, capacity bounds
+  /// and supply values. By default it is \c int.
+  /// \tparam C The number type used for costs and potentials.
+  /// By default it is the same as \c V.
+#ifdef DOXYGEN
+  template <typename GR, typename V = int, typename C = V>
+#else
+  template < typename GR, typename V = int, typename C = V,
+             bool integer = std::numeric_limits<C>::is_integer >
+#endif
+  struct CostScalingDefaultTraits
+  {
+    /// The type of the digraph
+    typedef GR Digraph;
+    /// The type of the flow amounts, capacity bounds and supply values
+    typedef V Value;
+    /// The type of the arc costs
+    typedef C Cost;
+
+    /// \brief The large cost type used for internal computations
+    ///
+    /// The large cost type used for internal computations.
+    /// It is \c long \c long if the \c Cost type is integer,
+    /// otherwise it is \c double.
+    /// \c Cost must be convertible to \c LargeCost.
+    typedef double LargeCost;
+  };
+
+  // Default traits class for integer cost types
+  template <typename GR, typename V, typename C>
+  struct CostScalingDefaultTraits<GR, V, C, true>
+  {
+    typedef GR Digraph;
+    typedef V Value;
+    typedef C Cost;
+#ifdef LEMON_HAVE_LONG_LONG
+    typedef long long LargeCost;
+#else
+    typedef long LargeCost;
+#endif
+  };
+
+
+  /// \addtogroup min_cost_flow_algs
+  /// @{
+
+  /// \brief Implementation of the Cost Scaling algorithm for
+  /// finding a \ref min_cost_flow "minimum cost flow".
+  ///
+  /// \ref CostScaling implements a cost scaling algorithm that performs
+  /// push/augment and relabel operations for finding a \ref min_cost_flow
+  /// "minimum cost flow" \ref amo93networkflows, \ref goldberg90approximation,
+  /// \ref goldberg97efficient, \ref bunnagel98efficient. 
+  /// It is a highly efficient primal-dual solution method, which
+  /// can be viewed as the generalization of the \ref Preflow
+  /// "preflow push-relabel" algorithm for the maximum flow problem.
+  ///
+  /// 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 V The number type used for flow amounts, capacity bounds
+  /// and supply values in the algorithm. By default it is \c int.
+  /// \tparam C The number type used for costs and potentials in the
+  /// algorithm. By default it is the same as \c V.
+  ///
+  /// \warning Both number types must be signed and all input data must
+  /// be integer.
+  /// \warning This algorithm does not support negative costs for such
+  /// arcs that have infinite upper bound.
+  ///
+  /// \note %CostScaling provides three different internal methods,
+  /// from which the most efficient one is used by default.
+  /// For more information, see \ref Method.
+#ifdef DOXYGEN
+  template <typename GR, typename V, typename C, typename TR>
+#else
+  template < typename GR, typename V = int, typename C = V,
+             typename TR = CostScalingDefaultTraits<GR, V, C> >
+#endif
+  class CostScaling
+  {
+  public:
+
+    /// The type of the digraph
+    typedef typename TR::Digraph Digraph;
+    /// The type of the flow amounts, capacity bounds and supply values
+    typedef typename TR::Value Value;
+    /// The type of the arc costs
+    typedef typename TR::Cost Cost;
+
+    /// \brief The large cost type
+    ///
+    /// The large cost type used for internal computations.
+    /// Using the \ref CostScalingDefaultTraits "default traits class",
+    /// it is \c long \c long if the \c Cost type is integer,
+    /// otherwise it is \c double.
+    typedef typename TR::LargeCost LargeCost;
+
+    /// The \ref CostScalingDefaultTraits "traits class" of the algorithm
+    typedef TR Traits;
+
+  public:
+
+    /// \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 digraph contains an arc of negative cost and infinite
+      /// upper bound. It means that the objective function is unbounded
+      /// on that arc, however, note that it could actually be bounded
+      /// over the feasible flows, but this algroithm cannot handle
+      /// these cases.
+      UNBOUNDED
+    };
+
+    /// \brief Constants for selecting the internal method.
+    ///
+    /// Enum type containing constants for selecting the internal method
+    /// for the \ref run() function.
+    ///
+    /// \ref CostScaling provides three internal methods that differ mainly
+    /// in their base operations, which are used in conjunction with the
+    /// relabel operation.
+    /// By default, the so called \ref PARTIAL_AUGMENT
+    /// "Partial Augment-Relabel" method is used, which proved to be
+    /// the most efficient and the most robust on various test inputs.
+    /// However, the other methods can be selected using the \ref run()
+    /// function with the proper parameter.
+    enum Method {
+      /// Local push operations are used, i.e. flow is moved only on one
+      /// admissible arc at once.
+      PUSH,
+      /// Augment operations are used, i.e. flow is moved on admissible
+      /// paths from a node with excess to a node with deficit.
+      AUGMENT,
+      /// Partial augment operations are used, i.e. flow is moved on 
+      /// admissible paths started from a node with excess, but the
+      /// lengths of these paths are limited. This method can be viewed
+      /// as a combined version of the previous two operations.
+      PARTIAL_AUGMENT
+    };
+
+  private:
+
+    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
+
+    typedef std::vector<int> IntVector;
+    typedef std::vector<char> BoolVector;
+    typedef std::vector<Value> ValueVector;
+    typedef std::vector<Cost> CostVector;
+    typedef std::vector<LargeCost> LargeCostVector;
+
+  private:
+  
+    template <typename KT, typename VT>
+    class StaticVectorMap {
+    public:
+      typedef KT Key;
+      typedef VT Value;
+      
+      StaticVectorMap(std::vector<Value>& v) : _v(v) {}
+      
+      const Value& operator[](const Key& key) const {
+        return _v[StaticDigraph::id(key)];
+      }
+
+      Value& operator[](const Key& key) {
+        return _v[StaticDigraph::id(key)];
+      }
+      
+      void set(const Key& key, const Value& val) {
+        _v[StaticDigraph::id(key)] = val;
+      }
+
+    private:
+      std::vector<Value>& _v;
+    };
+
+    typedef StaticVectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap;
+    typedef StaticVectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap;
+
+  private:
+
+    // Data related to the underlying digraph
+    const GR &_graph;
+    int _node_num;
+    int _arc_num;
+    int _res_node_num;
+    int _res_arc_num;
+    int _root;
+
+    // Parameters of the problem
+    bool _have_lower;
+    Value _sum_supply;
+
+    // Data structures for storing the digraph
+    IntNodeMap _node_id;
+    IntArcMap _arc_idf;
+    IntArcMap _arc_idb;
+    IntVector _first_out;
+    BoolVector _forward;
+    IntVector _source;
+    IntVector _target;
+    IntVector _reverse;
+
+    // Node and arc data
+    ValueVector _lower;
+    ValueVector _upper;
+    CostVector _scost;
+    ValueVector _supply;
+
+    ValueVector _res_cap;
+    LargeCostVector _cost;
+    LargeCostVector _pi;
+    ValueVector _excess;
+    IntVector _next_out;
+    std::deque<int> _active_nodes;
+
+    // Data for scaling
+    LargeCost _epsilon;
+    int _alpha;
+
+    // Data for a StaticDigraph structure
+    typedef std::pair<int, int> IntPair;
+    StaticDigraph _sgr;
+    std::vector<IntPair> _arc_vec;
+    std::vector<LargeCost> _cost_vec;
+    LargeCostArcMap _cost_map;
+    LargeCostNodeMap _pi_map;
+  
+  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;
+
+  public:
+
+    /// \name Named Template Parameters
+    /// @{
+
+    template <typename T>
+    struct SetLargeCostTraits : public Traits {
+      typedef T LargeCost;
+    };
+
+    /// \brief \ref named-templ-param "Named parameter" for setting
+    /// \c LargeCost type.
+    ///
+    /// \ref named-templ-param "Named parameter" for setting \c LargeCost
+    /// type, which is used for internal computations in the algorithm.
+    /// \c Cost must be convertible to \c LargeCost.
+    template <typename T>
+    struct SetLargeCost
+      : public CostScaling<GR, V, C, SetLargeCostTraits<T> > {
+      typedef  CostScaling<GR, V, C, SetLargeCostTraits<T> > Create;
+    };
+
+    /// @}
+
+  public:
+
+    /// \brief Constructor.
+    ///
+    /// The constructor of the class.
+    ///
+    /// \param graph The digraph the algorithm runs on.
+    CostScaling(const GR& graph) :
+      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
+      _cost_map(_cost_vec), _pi_map(_pi),
+      INF(std::numeric_limits<Value>::has_infinity ?
+          std::numeric_limits<Value>::infinity() :
+          std::numeric_limits<Value>::max())
+    {
+      // Check the number types
+      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
+        "The flow type of CostScaling must be signed");
+      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
+        "The cost type of CostScaling must be signed");
+
+      // Resize vectors
+      _node_num = countNodes(_graph);
+      _arc_num = countArcs(_graph);
+      _res_node_num = _node_num + 1;
+      _res_arc_num = 2 * (_arc_num + _node_num);
+      _root = _node_num;
+
+      _first_out.resize(_res_node_num + 1);
+      _forward.resize(_res_arc_num);
+      _source.resize(_res_arc_num);
+      _target.resize(_res_arc_num);
+      _reverse.resize(_res_arc_num);
+
+      _lower.resize(_res_arc_num);
+      _upper.resize(_res_arc_num);
+      _scost.resize(_res_arc_num);
+      _supply.resize(_res_node_num);
+      
+      _res_cap.resize(_res_arc_num);
+      _cost.resize(_res_arc_num);
+      _pi.resize(_res_node_num);
+      _excess.resize(_res_node_num);
+      _next_out.resize(_res_node_num);
+
+      _arc_vec.reserve(_res_arc_num);
+      _cost_vec.reserve(_res_arc_num);
+
+      // Copy the graph
+      int i = 0, j = 0, k = 2 * _arc_num + _node_num;
+      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
+        _node_id[n] = i;
+      }
+      i = 0;
+      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
+        _first_out[i] = j;
+        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
+          _arc_idf[a] = j;
+          _forward[j] = true;
+          _source[j] = i;
+          _target[j] = _node_id[_graph.runningNode(a)];
+        }
+        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
+          _arc_idb[a] = j;
+          _forward[j] = false;
+          _source[j] = i;
+          _target[j] = _node_id[_graph.runningNode(a)];
+        }
+        _forward[j] = false;
+        _source[j] = i;
+        _target[j] = _root;
+        _reverse[j] = k;
+        _forward[k] = true;
+        _source[k] = _root;
+        _target[k] = i;
+        _reverse[k] = j;
+        ++j; ++k;
+      }
+      _first_out[i] = j;
+      _first_out[_res_node_num] = k;
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        int fi = _arc_idf[a];
+        int bi = _arc_idb[a];
+        _reverse[fi] = bi;
+        _reverse[bi] = fi;
+      }
+      
+      // Reset parameters
+      reset();
+    }
+
+    /// \name Parameters
+    /// The parameters of the algorithm can be specified using these
+    /// functions.
+
+    /// @{
+
+    /// \brief Set the lower bounds on the arcs.
+    ///
+    /// This function sets the lower bounds on the 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 Value type
+    /// of the algorithm.
+    ///
+    /// \return <tt>(*this)</tt>
+    template <typename LowerMap>
+    CostScaling& lowerMap(const LowerMap& map) {
+      _have_lower = true;
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        _lower[_arc_idf[a]] = map[a];
+        _lower[_arc_idb[a]] = map[a];
+      }
+      return *this;
+    }
+
+    /// \brief Set the upper bounds (capacities) on the arcs.
+    ///
+    /// This function sets the upper bounds (capacities) on the 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).
+    ///
+    /// \param map An arc map storing the upper bounds.
+    /// Its \c Value type must be convertible to the \c Value type
+    /// of the algorithm.
+    ///
+    /// \return <tt>(*this)</tt>
+    template<typename UpperMap>
+    CostScaling& upperMap(const UpperMap& map) {
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        _upper[_arc_idf[a]] = map[a];
+      }
+      return *this;
+    }
+
+    /// \brief Set the costs of the arcs.
+    ///
+    /// This function sets the costs of the arcs.
+    /// If it is not used before calling \ref run(), the costs
+    /// will be set to \c 1 on all arcs.
+    ///
+    /// \param map An arc map storing the costs.
+    /// Its \c Value type must be convertible to the \c Cost type
+    /// of the algorithm.
+    ///
+    /// \return <tt>(*this)</tt>
+    template<typename CostMap>
+    CostScaling& costMap(const CostMap& map) {
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        _scost[_arc_idf[a]] =  map[a];
+        _scost[_arc_idb[a]] = -map[a];
+      }
+      return *this;
+    }
+
+    /// \brief Set the supply values of the nodes.
+    ///
+    /// This function sets the supply values of the nodes.
+    /// If neither this function nor \ref stSupply() is used before
+    /// calling \ref run(), the supply of each node will be set to zero.
+    ///
+    /// \param map A node map storing the supply values.
+    /// Its \c Value type must be convertible to the \c Value type
+    /// of the algorithm.
+    ///
+    /// \return <tt>(*this)</tt>
+    template<typename SupplyMap>
+    CostScaling& supplyMap(const SupplyMap& map) {
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        _supply[_node_id[n]] = map[n];
+      }
+      return *this;
+    }
+
+    /// \brief Set single source and target nodes and a supply value.
+    ///
+    /// This function sets a single source node and a single target node
+    /// and the required flow value.
+    /// If neither this function nor \ref supplyMap() is used before
+    /// calling \ref run(), the supply of each node will be set to zero.
+    ///
+    /// 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.
+    /// \param k The required amount of flow from node \c s to node \c t
+    /// (i.e. the supply of \c s and the demand of \c t).
+    ///
+    /// \return <tt>(*this)</tt>
+    CostScaling& stSupply(const Node& s, const Node& t, Value k) {
+      for (int i = 0; i != _res_node_num; ++i) {
+        _supply[i] = 0;
+      }
+      _supply[_node_id[s]] =  k;
+      _supply[_node_id[t]] = -k;
+      return *this;
+    }
+    
+    /// @}
+
+    /// \name Execution control
+    /// The algorithm can be executed using \ref run().
+
+    /// @{
+
+    /// \brief Run the algorithm.
+    ///
+    /// This function runs the algorithm.
+    /// The paramters can be specified using functions \ref lowerMap(),
+    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
+    /// For example,
+    /// \code
+    ///   CostScaling<ListDigraph> cs(graph);
+    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
+    ///     .supplyMap(sup).run();
+    /// \endcode
+    ///
+    /// This function can be called more than once. All the parameters
+    /// that have been given are kept for the next call, unless
+    /// \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 method The internal method that will be used in the
+    /// algorithm. For more information, see \ref Method.
+    /// \param factor The cost scaling factor. It must be larger than one.
+    ///
+    /// \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 digraph contains an arc of negative cost
+    /// and infinite upper bound. It means that the objective function
+    /// is unbounded on that arc, however, note that it could actually be
+    /// bounded over the feasible flows, but this algroithm cannot handle
+    /// these cases.
+    ///
+    /// \see ProblemType, Method
+    ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 8) {
+      _alpha = factor;
+      ProblemType pt = init();
+      if (pt != OPTIMAL) return pt;
+      start(method);
+      return OPTIMAL;
+    }
+
+    /// \brief Reset all the parameters that have been given before.
+    ///
+    /// This function resets all the paramaters that have been given
+    /// before using functions \ref lowerMap(), \ref upperMap(),
+    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
+    ///
+    /// 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,
+    /// \code
+    ///   CostScaling<ListDigraph> cs(graph);
+    ///
+    ///   // First run
+    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
+    ///     .supplyMap(sup).run();
+    ///
+    ///   // Run again with modified cost map (reset() is not called,
+    ///   // so only the cost map have to be set again)
+    ///   cost[e] += 100;
+    ///   cs.costMap(cost).run();
+    ///
+    ///   // Run again from scratch using reset()
+    ///   // (the lower bounds will be set to zero on all arcs)
+    ///   cs.reset();
+    ///   cs.upperMap(capacity).costMap(cost)
+    ///     .supplyMap(sup).run();
+    /// \endcode
+    ///
+    /// \return <tt>(*this)</tt>
+    CostScaling& reset() {
+      for (int i = 0; i != _res_node_num; ++i) {
+        _supply[i] = 0;
+      }
+      int limit = _first_out[_root];
+      for (int j = 0; j != limit; ++j) {
+        _lower[j] = 0;
+        _upper[j] = INF;
+        _scost[j] = _forward[j] ? 1 : -1;
+      }
+      for (int j = limit; j != _res_arc_num; ++j) {
+        _lower[j] = 0;
+        _upper[j] = INF;
+        _scost[j] = 0;
+        _scost[_reverse[j]] = 0;
+      }      
+      _have_lower = false;
+      return *this;
+    }
+
+    /// @}
+
+    /// \name Query Functions
+    /// The results of the algorithm can be obtained using these
+    /// functions.\n
+    /// The \ref run() function must be called before using them.
+
+    /// @{
+
+    /// \brief Return the total cost of the found flow.
+    ///
+    /// This function returns the total cost of the found flow.
+    /// Its complexity is O(e).
+    ///
+    /// \note The return type of the function can be specified as a
+    /// template parameter. For example,
+    /// \code
+    ///   cs.totalCost<double>();
+    /// \endcode
+    /// It is useful if the total cost cannot be stored in the \c Cost
+    /// type of the algorithm, which is the default return type of the
+    /// function.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    template <typename Number>
+    Number totalCost() const {
+      Number c = 0;
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        int i = _arc_idb[a];
+        c += static_cast<Number>(_res_cap[i]) *
+             (-static_cast<Number>(_scost[i]));
+      }
+      return c;
+    }
+
+#ifndef DOXYGEN
+    Cost totalCost() const {
+      return totalCost<Cost>();
+    }
+#endif
+
+    /// \brief Return the flow on the given arc.
+    ///
+    /// This function returns the flow on the given arc.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    Value flow(const Arc& a) const {
+      return _res_cap[_arc_idb[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.
+    template <typename FlowMap>
+    void flowMap(FlowMap &map) const {
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        map.set(a, _res_cap[_arc_idb[a]]);
+      }
+    }
+
+    /// \brief Return the potential (dual value) of the given node.
+    ///
+    /// This function returns the potential (dual value) of the
+    /// given node.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    Cost potential(const Node& n) const {
+      return static_cast<Cost>(_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.
+    template <typename PotentialMap>
+    void potentialMap(PotentialMap &map) const {
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        map.set(n, static_cast<Cost>(_pi[_node_id[n]]));
+      }
+    }
+
+    /// @}
+
+  private:
+
+    // Initialize the algorithm
+    ProblemType init() {
+      if (_res_node_num <= 1) return INFEASIBLE;
+
+      // Check the sum of supply values
+      _sum_supply = 0;
+      for (int i = 0; i != _root; ++i) {
+        _sum_supply += _supply[i];
+      }
+      if (_sum_supply > 0) return INFEASIBLE;
+      
+
+      // Initialize vectors
+      for (int i = 0; i != _res_node_num; ++i) {
+        _pi[i] = 0;
+        _excess[i] = _supply[i];
+      }
+      
+      // Remove infinite upper bounds and check negative arcs
+      const Value MAX = std::numeric_limits<Value>::max();
+      int last_out;
+      if (_have_lower) {
+        for (int i = 0; i != _root; ++i) {
+          last_out = _first_out[i+1];
+          for (int j = _first_out[i]; j != last_out; ++j) {
+            if (_forward[j]) {
+              Value c = _scost[j] < 0 ? _upper[j] : _lower[j];
+              if (c >= MAX) return UNBOUNDED;
+              _excess[i] -= c;
+              _excess[_target[j]] += c;
+            }
+          }
+        }
+      } else {
+        for (int i = 0; i != _root; ++i) {
+          last_out = _first_out[i+1];
+          for (int j = _first_out[i]; j != last_out; ++j) {
+            if (_forward[j] && _scost[j] < 0) {
+              Value c = _upper[j];
+              if (c >= MAX) return UNBOUNDED;
+              _excess[i] -= c;
+              _excess[_target[j]] += c;
+            }
+          }
+        }
+      }
+      Value ex, max_cap = 0;
+      for (int i = 0; i != _res_node_num; ++i) {
+        ex = _excess[i];
+        _excess[i] = 0;
+        if (ex < 0) max_cap -= ex;
+      }
+      for (int j = 0; j != _res_arc_num; ++j) {
+        if (_upper[j] >= MAX) _upper[j] = max_cap;
+      }
+
+      // Initialize the large cost vector and the epsilon parameter
+      _epsilon = 0;
+      LargeCost lc;
+      for (int i = 0; i != _root; ++i) {
+        last_out = _first_out[i+1];
+        for (int j = _first_out[i]; j != last_out; ++j) {
+          lc = static_cast<LargeCost>(_scost[j]) * _res_node_num * _alpha;
+          _cost[j] = lc;
+          if (lc > _epsilon) _epsilon = lc;
+        }
+      }
+      _epsilon /= _alpha;
+
+      // Initialize maps for Circulation and remove non-zero lower bounds
+      ConstMap<Arc, Value> low(0);
+      typedef typename Digraph::template ArcMap<Value> ValueArcMap;
+      typedef typename Digraph::template NodeMap<Value> ValueNodeMap;
+      ValueArcMap cap(_graph), flow(_graph);
+      ValueNodeMap sup(_graph);
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        sup[n] = _supply[_node_id[n]];
+      }
+      if (_have_lower) {
+        for (ArcIt a(_graph); a != INVALID; ++a) {
+          int j = _arc_idf[a];
+          Value c = _lower[j];
+          cap[a] = _upper[j] - c;
+          sup[_graph.source(a)] -= c;
+          sup[_graph.target(a)] += c;
+        }
+      } else {
+        for (ArcIt a(_graph); a != INVALID; ++a) {
+          cap[a] = _upper[_arc_idf[a]];
+        }
+      }
+
+      // Find a feasible flow using Circulation
+      Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
+        circ(_graph, low, cap, sup);
+      if (!circ.flowMap(flow).run()) return INFEASIBLE;
+
+      // Set residual capacities and handle GEQ supply type
+      if (_sum_supply < 0) {
+        for (ArcIt a(_graph); a != INVALID; ++a) {
+          Value fa = flow[a];
+          _res_cap[_arc_idf[a]] = cap[a] - fa;
+          _res_cap[_arc_idb[a]] = fa;
+          sup[_graph.source(a)] -= fa;
+          sup[_graph.target(a)] += fa;
+        }
+        for (NodeIt n(_graph); n != INVALID; ++n) {
+          _excess[_node_id[n]] = sup[n];
+        }
+        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
+          int u = _target[a];
+          int ra = _reverse[a];
+          _res_cap[a] = -_sum_supply + 1;
+          _res_cap[ra] = -_excess[u];
+          _cost[a] = 0;
+          _cost[ra] = 0;
+          _excess[u] = 0;
+        }
+      } else {
+        for (ArcIt a(_graph); a != INVALID; ++a) {
+          Value fa = flow[a];
+          _res_cap[_arc_idf[a]] = cap[a] - fa;
+          _res_cap[_arc_idb[a]] = fa;
+        }
+        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
+          int ra = _reverse[a];
+          _res_cap[a] = 1;
+          _res_cap[ra] = 0;
+          _cost[a] = 0;
+          _cost[ra] = 0;
+        }
+      }
+      
+      return OPTIMAL;
+    }
+
+    // Execute the algorithm and transform the results
+    void start(Method method) {
+      // Maximum path length for partial augment
+      const int MAX_PATH_LENGTH = 4;
+      
+      // Execute the algorithm
+      switch (method) {
+        case PUSH:
+          startPush();
+          break;
+        case AUGMENT:
+          startAugment();
+          break;
+        case PARTIAL_AUGMENT:
+          startAugment(MAX_PATH_LENGTH);
+          break;
+      }
+
+      // Compute node potentials for the original costs
+      _arc_vec.clear();
+      _cost_vec.clear();
+      for (int j = 0; j != _res_arc_num; ++j) {
+        if (_res_cap[j] > 0) {
+          _arc_vec.push_back(IntPair(_source[j], _target[j]));
+          _cost_vec.push_back(_scost[j]);
+        }
+      }
+      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
+
+      typename BellmanFord<StaticDigraph, LargeCostArcMap>
+        ::template SetDistMap<LargeCostNodeMap>::Create bf(_sgr, _cost_map);
+      bf.distMap(_pi_map);
+      bf.init(0);
+      bf.start();
+
+      // Handle non-zero lower bounds
+      if (_have_lower) {
+        int limit = _first_out[_root];
+        for (int j = 0; j != limit; ++j) {
+          if (!_forward[j]) _res_cap[j] += _lower[j];
+        }
+      }
+    }
+
+    /// Execute the algorithm performing augment and relabel operations
+    void startAugment(int max_length = std::numeric_limits<int>::max()) {
+      // Paramters for heuristics
+      const int BF_HEURISTIC_EPSILON_BOUND = 1000;
+      const int BF_HEURISTIC_BOUND_FACTOR  = 3;
+
+      // Perform cost scaling phases
+      IntVector pred_arc(_res_node_num);
+      std::vector<int> path_nodes;
+      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
+                                        1 : _epsilon / _alpha )
+      {
+        // "Early Termination" heuristic: use Bellman-Ford algorithm
+        // to check if the current flow is optimal
+        if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
+          _arc_vec.clear();
+          _cost_vec.clear();
+          for (int j = 0; j != _res_arc_num; ++j) {
+            if (_res_cap[j] > 0) {
+              _arc_vec.push_back(IntPair(_source[j], _target[j]));
+              _cost_vec.push_back(_cost[j] + 1);
+            }
+          }
+          _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
+
+          BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
+          bf.init(0);
+          bool done = false;
+          int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
+          for (int i = 0; i < K && !done; ++i)
+            done = bf.processNextWeakRound();
+          if (done) break;
+        }
+
+        // Saturate arcs not satisfying the optimality condition
+        for (int a = 0; a != _res_arc_num; ++a) {
+          if (_res_cap[a] > 0 &&
+              _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
+            Value delta = _res_cap[a];
+            _excess[_source[a]] -= delta;
+            _excess[_target[a]] += delta;
+            _res_cap[a] = 0;
+            _res_cap[_reverse[a]] += delta;
+          }
+        }
+        
+        // Find active nodes (i.e. nodes with positive excess)
+        for (int u = 0; u != _res_node_num; ++u) {
+          if (_excess[u] > 0) _active_nodes.push_back(u);
+        }
+
+        // Initialize the next arcs
+        for (int u = 0; u != _res_node_num; ++u) {
+          _next_out[u] = _first_out[u];
+        }
+
+        // Perform partial augment and relabel operations
+        while (true) {
+          // Select an active node (FIFO selection)
+          while (_active_nodes.size() > 0 &&
+                 _excess[_active_nodes.front()] <= 0) {
+            _active_nodes.pop_front();
+          }
+          if (_active_nodes.size() == 0) break;
+          int start = _active_nodes.front();
+          path_nodes.clear();
+          path_nodes.push_back(start);
+
+          // Find an augmenting path from the start node
+          int tip = start;
+          while (_excess[tip] >= 0 &&
+                 int(path_nodes.size()) <= max_length) {
+            int u;
+            LargeCost min_red_cost, rc;
+            int last_out = _sum_supply < 0 ?
+              _first_out[tip+1] : _first_out[tip+1] - 1;
+            for (int a = _next_out[tip]; a != last_out; ++a) {
+              if (_res_cap[a] > 0 &&
+                  _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
+                u = _target[a];
+                pred_arc[u] = a;
+                _next_out[tip] = a;
+                tip = u;
+                path_nodes.push_back(tip);
+                goto next_step;
+              }
+            }
+
+            // Relabel tip node
+            min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
+            for (int a = _first_out[tip]; a != last_out; ++a) {
+              rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
+              if (_res_cap[a] > 0 && rc < min_red_cost) {
+                min_red_cost = rc;
+              }
+            }
+            _pi[tip] -= min_red_cost + _epsilon;
+
+            // Reset the next arc of tip
+            _next_out[tip] = _first_out[tip];
+
+            // Step back
+            if (tip != start) {
+              path_nodes.pop_back();
+              tip = path_nodes.back();
+            }
+
+          next_step: ;
+          }
+
+          // Augment along the found path (as much flow as possible)
+          Value delta;
+          int u, v = path_nodes.front(), pa;
+          for (int i = 1; i < int(path_nodes.size()); ++i) {
+            u = v;
+            v = path_nodes[i];
+            pa = pred_arc[v];
+            delta = std::min(_res_cap[pa], _excess[u]);
+            _res_cap[pa] -= delta;
+            _res_cap[_reverse[pa]] += delta;
+            _excess[u] -= delta;
+            _excess[v] += delta;
+            if (_excess[v] > 0 && _excess[v] <= delta)
+              _active_nodes.push_back(v);
+          }
+        }
+      }
+    }
+
+    /// Execute the algorithm performing push and relabel operations
+    void startPush() {
+      // Paramters for heuristics
+      const int BF_HEURISTIC_EPSILON_BOUND = 1000;
+      const int BF_HEURISTIC_BOUND_FACTOR  = 3;
+
+      // Perform cost scaling phases
+      BoolVector hyper(_res_node_num, false);
+      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
+                                        1 : _epsilon / _alpha )
+      {
+        // "Early Termination" heuristic: use Bellman-Ford algorithm
+        // to check if the current flow is optimal
+        if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
+          _arc_vec.clear();
+          _cost_vec.clear();
+          for (int j = 0; j != _res_arc_num; ++j) {
+            if (_res_cap[j] > 0) {
+              _arc_vec.push_back(IntPair(_source[j], _target[j]));
+              _cost_vec.push_back(_cost[j] + 1);
+            }
+          }
+          _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
+
+          BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
+          bf.init(0);
+          bool done = false;
+          int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
+          for (int i = 0; i < K && !done; ++i)
+            done = bf.processNextWeakRound();
+          if (done) break;
+        }
+
+        // Saturate arcs not satisfying the optimality condition
+        for (int a = 0; a != _res_arc_num; ++a) {
+          if (_res_cap[a] > 0 &&
+              _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
+            Value delta = _res_cap[a];
+            _excess[_source[a]] -= delta;
+            _excess[_target[a]] += delta;
+            _res_cap[a] = 0;
+            _res_cap[_reverse[a]] += delta;
+          }
+        }
+
+        // Find active nodes (i.e. nodes with positive excess)
+        for (int u = 0; u != _res_node_num; ++u) {
+          if (_excess[u] > 0) _active_nodes.push_back(u);
+        }
+
+        // Initialize the next arcs
+        for (int u = 0; u != _res_node_num; ++u) {
+          _next_out[u] = _first_out[u];
+        }
+
+        // Perform push and relabel operations
+        while (_active_nodes.size() > 0) {
+          LargeCost min_red_cost, rc;
+          Value delta;
+          int n, t, a, last_out = _res_arc_num;
+
+          // Select an active node (FIFO selection)
+        next_node:
+          n = _active_nodes.front();
+          last_out = _sum_supply < 0 ?
+            _first_out[n+1] : _first_out[n+1] - 1;
+
+          // Perform push operations if there are admissible arcs
+          if (_excess[n] > 0) {
+            for (a = _next_out[n]; a != last_out; ++a) {
+              if (_res_cap[a] > 0 &&
+                  _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
+                delta = std::min(_res_cap[a], _excess[n]);
+                t = _target[a];
+
+                // Push-look-ahead heuristic
+                Value ahead = -_excess[t];
+                int last_out_t = _sum_supply < 0 ?
+                  _first_out[t+1] : _first_out[t+1] - 1;
+                for (int ta = _next_out[t]; ta != last_out_t; ++ta) {
+                  if (_res_cap[ta] > 0 && 
+                      _cost[ta] + _pi[_source[ta]] - _pi[_target[ta]] < 0)
+                    ahead += _res_cap[ta];
+                  if (ahead >= delta) break;
+                }
+                if (ahead < 0) ahead = 0;
+
+                // Push flow along the arc
+                if (ahead < delta) {
+                  _res_cap[a] -= ahead;
+                  _res_cap[_reverse[a]] += ahead;
+                  _excess[n] -= ahead;
+                  _excess[t] += ahead;
+                  _active_nodes.push_front(t);
+                  hyper[t] = true;
+                  _next_out[n] = a;
+                  goto next_node;
+                } else {
+                  _res_cap[a] -= delta;
+                  _res_cap[_reverse[a]] += delta;
+                  _excess[n] -= delta;
+                  _excess[t] += delta;
+                  if (_excess[t] > 0 && _excess[t] <= delta)
+                    _active_nodes.push_back(t);
+                }
+
+                if (_excess[n] == 0) {
+                  _next_out[n] = a;
+                  goto remove_nodes;
+                }
+              }
+            }
+            _next_out[n] = a;
+          }
+
+          // Relabel the node if it is still active (or hyper)
+          if (_excess[n] > 0 || hyper[n]) {
+            min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
+            for (int a = _first_out[n]; a != last_out; ++a) {
+              rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
+              if (_res_cap[a] > 0 && rc < min_red_cost) {
+                min_red_cost = rc;
+              }
+            }
+            _pi[n] -= min_red_cost + _epsilon;
+            hyper[n] = false;
+
+            // Reset the next arc
+            _next_out[n] = _first_out[n];
+          }
+        
+          // Remove nodes that are not active nor hyper
+        remove_nodes:
+          while ( _active_nodes.size() > 0 &&
+                  _excess[_active_nodes.front()] <= 0 &&
+                  !hyper[_active_nodes.front()] ) {
+            _active_nodes.pop_front();
+          }
+        }
+      }
+    }
+
+  }; //class CostScaling
+
+  ///@}
+
+} //namespace lemon
+
+#endif //LEMON_COST_SCALING_H
Index: lemon/cycle_canceling.h
===================================================================
--- lemon/cycle_canceling.h	(revision 886)
+++ lemon/cycle_canceling.h	(revision 886)
@@ -0,0 +1,1140 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_CYCLE_CANCELING_H
+#define LEMON_CYCLE_CANCELING_H
+
+/// \ingroup min_cost_flow_algs
+/// \file
+/// \brief Cycle-canceling algorithms for finding a minimum cost flow.
+
+#include <vector>
+#include <limits>
+
+#include <lemon/core.h>
+#include <lemon/maps.h>
+#include <lemon/path.h>
+#include <lemon/math.h>
+#include <lemon/static_graph.h>
+#include <lemon/adaptors.h>
+#include <lemon/circulation.h>
+#include <lemon/bellman_ford.h>
+#include <lemon/howard.h>
+
+namespace lemon {
+
+  /// \addtogroup min_cost_flow_algs
+  /// @{
+
+  /// \brief Implementation of cycle-canceling algorithms for
+  /// finding a \ref min_cost_flow "minimum cost flow".
+  ///
+  /// \ref CycleCanceling implements three different cycle-canceling
+  /// algorithms for finding a \ref min_cost_flow "minimum cost flow"
+  /// \ref amo93networkflows, \ref klein67primal,
+  /// \ref goldberg89cyclecanceling.
+  /// The most efficent one (both theoretically and practically)
+  /// is the \ref CANCEL_AND_TIGHTEN "Cancel and Tighten" algorithm,
+  /// thus it is the default method.
+  /// It is strongly polynomial, but in practice, it is typically much
+  /// slower than the scaling algorithms and NetworkSimplex.
+  ///
+  /// 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 V The number type used for flow amounts, capacity bounds
+  /// and supply values in the algorithm. By default, it is \c int.
+  /// \tparam C The number type used for costs and potentials in the
+  /// algorithm. By default, it is the same as \c V.
+  ///
+  /// \warning Both number types must be signed and all input data must
+  /// be integer.
+  /// \warning This algorithm does not support negative costs for such
+  /// arcs that have infinite upper bound.
+  ///
+  /// \note For more information about the three available methods,
+  /// see \ref Method.
+#ifdef DOXYGEN
+  template <typename GR, typename V, typename C>
+#else
+  template <typename GR, typename V = int, typename C = V>
+#endif
+  class CycleCanceling
+  {
+  public:
+
+    /// The type of the digraph
+    typedef GR Digraph;
+    /// The type of the flow amounts, capacity bounds and supply values
+    typedef V Value;
+    /// The type of the arc costs
+    typedef C Cost;
+
+  public:
+
+    /// \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 digraph contains an arc of negative cost and infinite
+      /// upper bound. It means that the objective function is unbounded
+      /// on that arc, however, note that it could actually be bounded
+      /// over the feasible flows, but this algroithm cannot handle
+      /// these cases.
+      UNBOUNDED
+    };
+
+    /// \brief Constants for selecting the used method.
+    ///
+    /// Enum type containing constants for selecting the used method
+    /// for the \ref run() function.
+    ///
+    /// \ref CycleCanceling provides three different cycle-canceling
+    /// methods. By default, \ref CANCEL_AND_TIGHTEN "Cancel and Tighten"
+    /// is used, which proved to be the most efficient and the most robust
+    /// on various test inputs.
+    /// However, the other methods can be selected using the \ref run()
+    /// function with the proper parameter.
+    enum Method {
+      /// A simple cycle-canceling method, which uses the
+      /// \ref BellmanFord "Bellman-Ford" algorithm with limited iteration
+      /// number for detecting negative cycles in the residual network.
+      SIMPLE_CYCLE_CANCELING,
+      /// The "Minimum Mean Cycle-Canceling" algorithm, which is a
+      /// well-known strongly polynomial method
+      /// \ref goldberg89cyclecanceling. It improves along a
+      /// \ref min_mean_cycle "minimum mean cycle" in each iteration.
+      /// Its running time complexity is O(n<sup>2</sup>m<sup>3</sup>log(n)).
+      MINIMUM_MEAN_CYCLE_CANCELING,
+      /// The "Cancel And Tighten" algorithm, which can be viewed as an
+      /// improved version of the previous method
+      /// \ref goldberg89cyclecanceling.
+      /// It is faster both in theory and in practice, its running time
+      /// complexity is O(n<sup>2</sup>m<sup>2</sup>log(n)).
+      CANCEL_AND_TIGHTEN
+    };
+
+  private:
+
+    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
+    
+    typedef std::vector<int> IntVector;
+    typedef std::vector<char> CharVector;
+    typedef std::vector<double> DoubleVector;
+    typedef std::vector<Value> ValueVector;
+    typedef std::vector<Cost> CostVector;
+
+  private:
+  
+    template <typename KT, typename VT>
+    class StaticVectorMap {
+    public:
+      typedef KT Key;
+      typedef VT Value;
+      
+      StaticVectorMap(std::vector<Value>& v) : _v(v) {}
+      
+      const Value& operator[](const Key& key) const {
+        return _v[StaticDigraph::id(key)];
+      }
+
+      Value& operator[](const Key& key) {
+        return _v[StaticDigraph::id(key)];
+      }
+      
+      void set(const Key& key, const Value& val) {
+        _v[StaticDigraph::id(key)] = val;
+      }
+
+    private:
+      std::vector<Value>& _v;
+    };
+
+    typedef StaticVectorMap<StaticDigraph::Node, Cost> CostNodeMap;
+    typedef StaticVectorMap<StaticDigraph::Arc, Cost> CostArcMap;
+
+  private:
+
+
+    // Data related to the underlying digraph
+    const GR &_graph;
+    int _node_num;
+    int _arc_num;
+    int _res_node_num;
+    int _res_arc_num;
+    int _root;
+
+    // Parameters of the problem
+    bool _have_lower;
+    Value _sum_supply;
+
+    // Data structures for storing the digraph
+    IntNodeMap _node_id;
+    IntArcMap _arc_idf;
+    IntArcMap _arc_idb;
+    IntVector _first_out;
+    CharVector _forward;
+    IntVector _source;
+    IntVector _target;
+    IntVector _reverse;
+
+    // Node and arc data
+    ValueVector _lower;
+    ValueVector _upper;
+    CostVector _cost;
+    ValueVector _supply;
+
+    ValueVector _res_cap;
+    CostVector _pi;
+
+    // Data for a StaticDigraph structure
+    typedef std::pair<int, int> IntPair;
+    StaticDigraph _sgr;
+    std::vector<IntPair> _arc_vec;
+    std::vector<Cost> _cost_vec;
+    IntVector _id_vec;
+    CostArcMap _cost_map;
+    CostNodeMap _pi_map;
+  
+  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;
+
+  public:
+
+    /// \brief Constructor.
+    ///
+    /// The constructor of the class.
+    ///
+    /// \param graph The digraph the algorithm runs on.
+    CycleCanceling(const GR& graph) :
+      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
+      _cost_map(_cost_vec), _pi_map(_pi),
+      INF(std::numeric_limits<Value>::has_infinity ?
+          std::numeric_limits<Value>::infinity() :
+          std::numeric_limits<Value>::max())
+    {
+      // Check the number types
+      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
+        "The flow type of CycleCanceling must be signed");
+      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
+        "The cost type of CycleCanceling must be signed");
+
+      // Resize vectors
+      _node_num = countNodes(_graph);
+      _arc_num = countArcs(_graph);
+      _res_node_num = _node_num + 1;
+      _res_arc_num = 2 * (_arc_num + _node_num);
+      _root = _node_num;
+
+      _first_out.resize(_res_node_num + 1);
+      _forward.resize(_res_arc_num);
+      _source.resize(_res_arc_num);
+      _target.resize(_res_arc_num);
+      _reverse.resize(_res_arc_num);
+
+      _lower.resize(_res_arc_num);
+      _upper.resize(_res_arc_num);
+      _cost.resize(_res_arc_num);
+      _supply.resize(_res_node_num);
+      
+      _res_cap.resize(_res_arc_num);
+      _pi.resize(_res_node_num);
+
+      _arc_vec.reserve(_res_arc_num);
+      _cost_vec.reserve(_res_arc_num);
+      _id_vec.reserve(_res_arc_num);
+
+      // Copy the graph
+      int i = 0, j = 0, k = 2 * _arc_num + _node_num;
+      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
+        _node_id[n] = i;
+      }
+      i = 0;
+      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
+        _first_out[i] = j;
+        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
+          _arc_idf[a] = j;
+          _forward[j] = true;
+          _source[j] = i;
+          _target[j] = _node_id[_graph.runningNode(a)];
+        }
+        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
+          _arc_idb[a] = j;
+          _forward[j] = false;
+          _source[j] = i;
+          _target[j] = _node_id[_graph.runningNode(a)];
+        }
+        _forward[j] = false;
+        _source[j] = i;
+        _target[j] = _root;
+        _reverse[j] = k;
+        _forward[k] = true;
+        _source[k] = _root;
+        _target[k] = i;
+        _reverse[k] = j;
+        ++j; ++k;
+      }
+      _first_out[i] = j;
+      _first_out[_res_node_num] = k;
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        int fi = _arc_idf[a];
+        int bi = _arc_idb[a];
+        _reverse[fi] = bi;
+        _reverse[bi] = fi;
+      }
+      
+      // Reset parameters
+      reset();
+    }
+
+    /// \name Parameters
+    /// The parameters of the algorithm can be specified using these
+    /// functions.
+
+    /// @{
+
+    /// \brief Set the lower bounds on the arcs.
+    ///
+    /// This function sets the lower bounds on the 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 Value type
+    /// of the algorithm.
+    ///
+    /// \return <tt>(*this)</tt>
+    template <typename LowerMap>
+    CycleCanceling& lowerMap(const LowerMap& map) {
+      _have_lower = true;
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        _lower[_arc_idf[a]] = map[a];
+        _lower[_arc_idb[a]] = map[a];
+      }
+      return *this;
+    }
+
+    /// \brief Set the upper bounds (capacities) on the arcs.
+    ///
+    /// This function sets the upper bounds (capacities) on the 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).
+    ///
+    /// \param map An arc map storing the upper bounds.
+    /// Its \c Value type must be convertible to the \c Value type
+    /// of the algorithm.
+    ///
+    /// \return <tt>(*this)</tt>
+    template<typename UpperMap>
+    CycleCanceling& upperMap(const UpperMap& map) {
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        _upper[_arc_idf[a]] = map[a];
+      }
+      return *this;
+    }
+
+    /// \brief Set the costs of the arcs.
+    ///
+    /// This function sets the costs of the arcs.
+    /// If it is not used before calling \ref run(), the costs
+    /// will be set to \c 1 on all arcs.
+    ///
+    /// \param map An arc map storing the costs.
+    /// Its \c Value type must be convertible to the \c Cost type
+    /// of the algorithm.
+    ///
+    /// \return <tt>(*this)</tt>
+    template<typename CostMap>
+    CycleCanceling& costMap(const CostMap& map) {
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        _cost[_arc_idf[a]] =  map[a];
+        _cost[_arc_idb[a]] = -map[a];
+      }
+      return *this;
+    }
+
+    /// \brief Set the supply values of the nodes.
+    ///
+    /// This function sets the supply values of the nodes.
+    /// If neither this function nor \ref stSupply() is used before
+    /// calling \ref run(), the supply of each node will be set to zero.
+    ///
+    /// \param map A node map storing the supply values.
+    /// Its \c Value type must be convertible to the \c Value type
+    /// of the algorithm.
+    ///
+    /// \return <tt>(*this)</tt>
+    template<typename SupplyMap>
+    CycleCanceling& supplyMap(const SupplyMap& map) {
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        _supply[_node_id[n]] = map[n];
+      }
+      return *this;
+    }
+
+    /// \brief Set single source and target nodes and a supply value.
+    ///
+    /// This function sets a single source node and a single target node
+    /// and the required flow value.
+    /// If neither this function nor \ref supplyMap() is used before
+    /// calling \ref run(), the supply of each node will be set to zero.
+    ///
+    /// 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.
+    /// \param k The required amount of flow from node \c s to node \c t
+    /// (i.e. the supply of \c s and the demand of \c t).
+    ///
+    /// \return <tt>(*this)</tt>
+    CycleCanceling& stSupply(const Node& s, const Node& t, Value k) {
+      for (int i = 0; i != _res_node_num; ++i) {
+        _supply[i] = 0;
+      }
+      _supply[_node_id[s]] =  k;
+      _supply[_node_id[t]] = -k;
+      return *this;
+    }
+    
+    /// @}
+
+    /// \name Execution control
+    /// The algorithm can be executed using \ref run().
+
+    /// @{
+
+    /// \brief Run the algorithm.
+    ///
+    /// This function runs the algorithm.
+    /// The paramters can be specified using functions \ref lowerMap(),
+    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
+    /// For example,
+    /// \code
+    ///   CycleCanceling<ListDigraph> cc(graph);
+    ///   cc.lowerMap(lower).upperMap(upper).costMap(cost)
+    ///     .supplyMap(sup).run();
+    /// \endcode
+    ///
+    /// This function can be called more than once. All the parameters
+    /// that have been given are kept for the next call, unless
+    /// \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 method The cycle-canceling method that will be used.
+    /// For more information, see \ref Method.
+    ///
+    /// \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 digraph contains an arc of negative cost
+    /// and infinite upper bound. It means that the objective function
+    /// is unbounded on that arc, however, note that it could actually be
+    /// bounded over the feasible flows, but this algroithm cannot handle
+    /// these cases.
+    ///
+    /// \see ProblemType, Method
+    ProblemType run(Method method = CANCEL_AND_TIGHTEN) {
+      ProblemType pt = init();
+      if (pt != OPTIMAL) return pt;
+      start(method);
+      return OPTIMAL;
+    }
+
+    /// \brief Reset all the parameters that have been given before.
+    ///
+    /// This function resets all the paramaters that have been given
+    /// before using functions \ref lowerMap(), \ref upperMap(),
+    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
+    ///
+    /// 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,
+    /// \code
+    ///   CycleCanceling<ListDigraph> cs(graph);
+    ///
+    ///   // First run
+    ///   cc.lowerMap(lower).upperMap(upper).costMap(cost)
+    ///     .supplyMap(sup).run();
+    ///
+    ///   // Run again with modified cost map (reset() is not called,
+    ///   // so only the cost map have to be set again)
+    ///   cost[e] += 100;
+    ///   cc.costMap(cost).run();
+    ///
+    ///   // Run again from scratch using reset()
+    ///   // (the lower bounds will be set to zero on all arcs)
+    ///   cc.reset();
+    ///   cc.upperMap(capacity).costMap(cost)
+    ///     .supplyMap(sup).run();
+    /// \endcode
+    ///
+    /// \return <tt>(*this)</tt>
+    CycleCanceling& reset() {
+      for (int i = 0; i != _res_node_num; ++i) {
+        _supply[i] = 0;
+      }
+      int limit = _first_out[_root];
+      for (int j = 0; j != limit; ++j) {
+        _lower[j] = 0;
+        _upper[j] = INF;
+        _cost[j] = _forward[j] ? 1 : -1;
+      }
+      for (int j = limit; j != _res_arc_num; ++j) {
+        _lower[j] = 0;
+        _upper[j] = INF;
+        _cost[j] = 0;
+        _cost[_reverse[j]] = 0;
+      }      
+      _have_lower = false;
+      return *this;
+    }
+
+    /// @}
+
+    /// \name Query Functions
+    /// The results of the algorithm can be obtained using these
+    /// functions.\n
+    /// The \ref run() function must be called before using them.
+
+    /// @{
+
+    /// \brief Return the total cost of the found flow.
+    ///
+    /// This function returns the total cost of the found flow.
+    /// Its complexity is O(e).
+    ///
+    /// \note The return type of the function can be specified as a
+    /// template parameter. For example,
+    /// \code
+    ///   cc.totalCost<double>();
+    /// \endcode
+    /// It is useful if the total cost cannot be stored in the \c Cost
+    /// type of the algorithm, which is the default return type of the
+    /// function.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    template <typename Number>
+    Number totalCost() const {
+      Number c = 0;
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        int i = _arc_idb[a];
+        c += static_cast<Number>(_res_cap[i]) *
+             (-static_cast<Number>(_cost[i]));
+      }
+      return c;
+    }
+
+#ifndef DOXYGEN
+    Cost totalCost() const {
+      return totalCost<Cost>();
+    }
+#endif
+
+    /// \brief Return the flow on the given arc.
+    ///
+    /// This function returns the flow on the given arc.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    Value flow(const Arc& a) const {
+      return _res_cap[_arc_idb[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.
+    template <typename FlowMap>
+    void flowMap(FlowMap &map) const {
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        map.set(a, _res_cap[_arc_idb[a]]);
+      }
+    }
+
+    /// \brief Return the potential (dual value) of the given node.
+    ///
+    /// This function returns the potential (dual value) of the
+    /// given node.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    Cost potential(const Node& n) const {
+      return static_cast<Cost>(_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.
+    template <typename PotentialMap>
+    void potentialMap(PotentialMap &map) const {
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        map.set(n, static_cast<Cost>(_pi[_node_id[n]]));
+      }
+    }
+
+    /// @}
+
+  private:
+
+    // Initialize the algorithm
+    ProblemType init() {
+      if (_res_node_num <= 1) return INFEASIBLE;
+
+      // Check the sum of supply values
+      _sum_supply = 0;
+      for (int i = 0; i != _root; ++i) {
+        _sum_supply += _supply[i];
+      }
+      if (_sum_supply > 0) return INFEASIBLE;
+      
+
+      // Initialize vectors
+      for (int i = 0; i != _res_node_num; ++i) {
+        _pi[i] = 0;
+      }
+      ValueVector excess(_supply);
+      
+      // Remove infinite upper bounds and check negative arcs
+      const Value MAX = std::numeric_limits<Value>::max();
+      int last_out;
+      if (_have_lower) {
+        for (int i = 0; i != _root; ++i) {
+          last_out = _first_out[i+1];
+          for (int j = _first_out[i]; j != last_out; ++j) {
+            if (_forward[j]) {
+              Value c = _cost[j] < 0 ? _upper[j] : _lower[j];
+              if (c >= MAX) return UNBOUNDED;
+              excess[i] -= c;
+              excess[_target[j]] += c;
+            }
+          }
+        }
+      } else {
+        for (int i = 0; i != _root; ++i) {
+          last_out = _first_out[i+1];
+          for (int j = _first_out[i]; j != last_out; ++j) {
+            if (_forward[j] && _cost[j] < 0) {
+              Value c = _upper[j];
+              if (c >= MAX) return UNBOUNDED;
+              excess[i] -= c;
+              excess[_target[j]] += c;
+            }
+          }
+        }
+      }
+      Value ex, max_cap = 0;
+      for (int i = 0; i != _res_node_num; ++i) {
+        ex = excess[i];
+        if (ex < 0) max_cap -= ex;
+      }
+      for (int j = 0; j != _res_arc_num; ++j) {
+        if (_upper[j] >= MAX) _upper[j] = max_cap;
+      }
+
+      // Initialize maps for Circulation and remove non-zero lower bounds
+      ConstMap<Arc, Value> low(0);
+      typedef typename Digraph::template ArcMap<Value> ValueArcMap;
+      typedef typename Digraph::template NodeMap<Value> ValueNodeMap;
+      ValueArcMap cap(_graph), flow(_graph);
+      ValueNodeMap sup(_graph);
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        sup[n] = _supply[_node_id[n]];
+      }
+      if (_have_lower) {
+        for (ArcIt a(_graph); a != INVALID; ++a) {
+          int j = _arc_idf[a];
+          Value c = _lower[j];
+          cap[a] = _upper[j] - c;
+          sup[_graph.source(a)] -= c;
+          sup[_graph.target(a)] += c;
+        }
+      } else {
+        for (ArcIt a(_graph); a != INVALID; ++a) {
+          cap[a] = _upper[_arc_idf[a]];
+        }
+      }
+
+      // Find a feasible flow using Circulation
+      Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
+        circ(_graph, low, cap, sup);
+      if (!circ.flowMap(flow).run()) return INFEASIBLE;
+
+      // Set residual capacities and handle GEQ supply type
+      if (_sum_supply < 0) {
+        for (ArcIt a(_graph); a != INVALID; ++a) {
+          Value fa = flow[a];
+          _res_cap[_arc_idf[a]] = cap[a] - fa;
+          _res_cap[_arc_idb[a]] = fa;
+          sup[_graph.source(a)] -= fa;
+          sup[_graph.target(a)] += fa;
+        }
+        for (NodeIt n(_graph); n != INVALID; ++n) {
+          excess[_node_id[n]] = sup[n];
+        }
+        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
+          int u = _target[a];
+          int ra = _reverse[a];
+          _res_cap[a] = -_sum_supply + 1;
+          _res_cap[ra] = -excess[u];
+          _cost[a] = 0;
+          _cost[ra] = 0;
+        }
+      } else {
+        for (ArcIt a(_graph); a != INVALID; ++a) {
+          Value fa = flow[a];
+          _res_cap[_arc_idf[a]] = cap[a] - fa;
+          _res_cap[_arc_idb[a]] = fa;
+        }
+        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
+          int ra = _reverse[a];
+          _res_cap[a] = 1;
+          _res_cap[ra] = 0;
+          _cost[a] = 0;
+          _cost[ra] = 0;
+        }
+      }
+      
+      return OPTIMAL;
+    }
+    
+    // Build a StaticDigraph structure containing the current
+    // residual network
+    void buildResidualNetwork() {
+      _arc_vec.clear();
+      _cost_vec.clear();
+      _id_vec.clear();
+      for (int j = 0; j != _res_arc_num; ++j) {
+        if (_res_cap[j] > 0) {
+          _arc_vec.push_back(IntPair(_source[j], _target[j]));
+          _cost_vec.push_back(_cost[j]);
+          _id_vec.push_back(j);
+        }
+      }
+      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
+    }
+
+    // Execute the algorithm and transform the results
+    void start(Method method) {
+      // Execute the algorithm
+      switch (method) {
+        case SIMPLE_CYCLE_CANCELING:
+          startSimpleCycleCanceling();
+          break;
+        case MINIMUM_MEAN_CYCLE_CANCELING:
+          startMinMeanCycleCanceling();
+          break;
+        case CANCEL_AND_TIGHTEN:
+          startCancelAndTighten();
+          break;
+      }
+
+      // Compute node potentials
+      if (method != SIMPLE_CYCLE_CANCELING) {
+        buildResidualNetwork();
+        typename BellmanFord<StaticDigraph, CostArcMap>
+          ::template SetDistMap<CostNodeMap>::Create bf(_sgr, _cost_map);
+        bf.distMap(_pi_map);
+        bf.init(0);
+        bf.start();
+      }
+
+      // Handle non-zero lower bounds
+      if (_have_lower) {
+        int limit = _first_out[_root];
+        for (int j = 0; j != limit; ++j) {
+          if (!_forward[j]) _res_cap[j] += _lower[j];
+        }
+      }
+    }
+
+    // Execute the "Simple Cycle Canceling" method
+    void startSimpleCycleCanceling() {
+      // Constants for computing the iteration limits
+      const int BF_FIRST_LIMIT  = 2;
+      const double BF_LIMIT_FACTOR = 1.5;
+      
+      typedef StaticVectorMap<StaticDigraph::Arc, Value> FilterMap;
+      typedef FilterArcs<StaticDigraph, FilterMap> ResDigraph;
+      typedef StaticVectorMap<StaticDigraph::Node, StaticDigraph::Arc> PredMap;
+      typedef typename BellmanFord<ResDigraph, CostArcMap>
+        ::template SetDistMap<CostNodeMap>
+        ::template SetPredMap<PredMap>::Create BF;
+      
+      // Build the residual network
+      _arc_vec.clear();
+      _cost_vec.clear();
+      for (int j = 0; j != _res_arc_num; ++j) {
+        _arc_vec.push_back(IntPair(_source[j], _target[j]));
+        _cost_vec.push_back(_cost[j]);
+      }
+      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
+
+      FilterMap filter_map(_res_cap);
+      ResDigraph rgr(_sgr, filter_map);
+      std::vector<int> cycle;
+      std::vector<StaticDigraph::Arc> pred(_res_arc_num);
+      PredMap pred_map(pred);
+      BF bf(rgr, _cost_map);
+      bf.distMap(_pi_map).predMap(pred_map);
+
+      int length_bound = BF_FIRST_LIMIT;
+      bool optimal = false;
+      while (!optimal) {
+        bf.init(0);
+        int iter_num = 0;
+        bool cycle_found = false;
+        while (!cycle_found) {
+          // Perform some iterations of the Bellman-Ford algorithm
+          int curr_iter_num = iter_num + length_bound <= _node_num ?
+            length_bound : _node_num - iter_num;
+          iter_num += curr_iter_num;
+          int real_iter_num = curr_iter_num;
+          for (int i = 0; i < curr_iter_num; ++i) {
+            if (bf.processNextWeakRound()) {
+              real_iter_num = i;
+              break;
+            }
+          }
+          if (real_iter_num < curr_iter_num) {
+            // Optimal flow is found
+            optimal = true;
+            break;
+          } else {
+            // Search for node disjoint negative cycles
+            std::vector<int> state(_res_node_num, 0);
+            int id = 0;
+            for (int u = 0; u != _res_node_num; ++u) {
+              if (state[u] != 0) continue;
+              ++id;
+              int v = u;
+              for (; v != -1 && state[v] == 0; v = pred[v] == INVALID ?
+                   -1 : rgr.id(rgr.source(pred[v]))) {
+                state[v] = id;
+              }
+              if (v != -1 && state[v] == id) {
+                // A negative cycle is found
+                cycle_found = true;
+                cycle.clear();
+                StaticDigraph::Arc a = pred[v];
+                Value d, delta = _res_cap[rgr.id(a)];
+                cycle.push_back(rgr.id(a));
+                while (rgr.id(rgr.source(a)) != v) {
+                  a = pred_map[rgr.source(a)];
+                  d = _res_cap[rgr.id(a)];
+                  if (d < delta) delta = d;
+                  cycle.push_back(rgr.id(a));
+                }
+
+                // Augment along the cycle
+                for (int i = 0; i < int(cycle.size()); ++i) {
+                  int j = cycle[i];
+                  _res_cap[j] -= delta;
+                  _res_cap[_reverse[j]] += delta;
+                }
+              }
+            }
+          }
+
+          // Increase iteration limit if no cycle is found
+          if (!cycle_found) {
+            length_bound = static_cast<int>(length_bound * BF_LIMIT_FACTOR);
+          }
+        }
+      }
+    }
+
+    // Execute the "Minimum Mean Cycle Canceling" method
+    void startMinMeanCycleCanceling() {
+      typedef SimplePath<StaticDigraph> SPath;
+      typedef typename SPath::ArcIt SPathArcIt;
+      typedef typename Howard<StaticDigraph, CostArcMap>
+        ::template SetPath<SPath>::Create MMC;
+      
+      SPath cycle;
+      MMC mmc(_sgr, _cost_map);
+      mmc.cycle(cycle);
+      buildResidualNetwork();
+      while (mmc.findMinMean() && mmc.cycleLength() < 0) {
+        // Find the cycle
+        mmc.findCycle();
+
+        // Compute delta value
+        Value delta = INF;
+        for (SPathArcIt a(cycle); a != INVALID; ++a) {
+          Value d = _res_cap[_id_vec[_sgr.id(a)]];
+          if (d < delta) delta = d;
+        }
+
+        // Augment along the cycle
+        for (SPathArcIt a(cycle); a != INVALID; ++a) {
+          int j = _id_vec[_sgr.id(a)];
+          _res_cap[j] -= delta;
+          _res_cap[_reverse[j]] += delta;
+        }
+
+        // Rebuild the residual network        
+        buildResidualNetwork();
+      }
+    }
+
+    // Execute the "Cancel And Tighten" method
+    void startCancelAndTighten() {
+      // Constants for the min mean cycle computations
+      const double LIMIT_FACTOR = 1.0;
+      const int MIN_LIMIT = 5;
+
+      // Contruct auxiliary data vectors
+      DoubleVector pi(_res_node_num, 0.0);
+      IntVector level(_res_node_num);
+      CharVector reached(_res_node_num);
+      CharVector processed(_res_node_num);
+      IntVector pred_node(_res_node_num);
+      IntVector pred_arc(_res_node_num);
+      std::vector<int> stack(_res_node_num);
+      std::vector<int> proc_vector(_res_node_num);
+
+      // Initialize epsilon
+      double epsilon = 0;
+      for (int a = 0; a != _res_arc_num; ++a) {
+        if (_res_cap[a] > 0 && -_cost[a] > epsilon)
+          epsilon = -_cost[a];
+      }
+
+      // Start phases
+      Tolerance<double> tol;
+      tol.epsilon(1e-6);
+      int limit = int(LIMIT_FACTOR * std::sqrt(double(_res_node_num)));
+      if (limit < MIN_LIMIT) limit = MIN_LIMIT;
+      int iter = limit;
+      while (epsilon * _res_node_num >= 1) {
+        // Find and cancel cycles in the admissible network using DFS
+        for (int u = 0; u != _res_node_num; ++u) {
+          reached[u] = false;
+          processed[u] = false;
+        }
+        int stack_head = -1;
+        int proc_head = -1;
+        for (int start = 0; start != _res_node_num; ++start) {
+          if (reached[start]) continue;
+
+          // New start node
+          reached[start] = true;
+          pred_arc[start] = -1;
+          pred_node[start] = -1;
+
+          // Find the first admissible outgoing arc
+          double p = pi[start];
+          int a = _first_out[start];
+          int last_out = _first_out[start+1];
+          for (; a != last_out && (_res_cap[a] == 0 ||
+               !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
+          if (a == last_out) {
+            processed[start] = true;
+            proc_vector[++proc_head] = start;
+            continue;
+          }
+          stack[++stack_head] = a;
+
+          while (stack_head >= 0) {
+            int sa = stack[stack_head];
+            int u = _source[sa];
+            int v = _target[sa];
+
+            if (!reached[v]) {
+              // A new node is reached
+              reached[v] = true;
+              pred_node[v] = u;
+              pred_arc[v] = sa;
+              p = pi[v];
+              a = _first_out[v];
+              last_out = _first_out[v+1];
+              for (; a != last_out && (_res_cap[a] == 0 ||
+                   !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
+              stack[++stack_head] = a == last_out ? -1 : a;
+            } else {
+              if (!processed[v]) {
+                // A cycle is found
+                int n, w = u;
+                Value d, delta = _res_cap[sa];
+                for (n = u; n != v; n = pred_node[n]) {
+                  d = _res_cap[pred_arc[n]];
+                  if (d <= delta) {
+                    delta = d;
+                    w = pred_node[n];
+                  }
+                }
+
+                // Augment along the cycle
+                _res_cap[sa] -= delta;
+                _res_cap[_reverse[sa]] += delta;
+                for (n = u; n != v; n = pred_node[n]) {
+                  int pa = pred_arc[n];
+                  _res_cap[pa] -= delta;
+                  _res_cap[_reverse[pa]] += delta;
+                }
+                for (n = u; stack_head > 0 && n != w; n = pred_node[n]) {
+                  --stack_head;
+                  reached[n] = false;
+                }
+                u = w;
+              }
+              v = u;
+
+              // Find the next admissible outgoing arc
+              p = pi[v];
+              a = stack[stack_head] + 1;
+              last_out = _first_out[v+1];
+              for (; a != last_out && (_res_cap[a] == 0 ||
+                   !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
+              stack[stack_head] = a == last_out ? -1 : a;
+            }
+
+            while (stack_head >= 0 && stack[stack_head] == -1) {
+              processed[v] = true;
+              proc_vector[++proc_head] = v;
+              if (--stack_head >= 0) {
+                // Find the next admissible outgoing arc
+                v = _source[stack[stack_head]];
+                p = pi[v];
+                a = stack[stack_head] + 1;
+                last_out = _first_out[v+1];
+                for (; a != last_out && (_res_cap[a] == 0 ||
+                     !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
+                stack[stack_head] = a == last_out ? -1 : a;
+              }
+            }
+          }
+        }
+
+        // Tighten potentials and epsilon
+        if (--iter > 0) {
+          for (int u = 0; u != _res_node_num; ++u) {
+            level[u] = 0;
+          }
+          for (int i = proc_head; i > 0; --i) {
+            int u = proc_vector[i];
+            double p = pi[u];
+            int l = level[u] + 1;
+            int last_out = _first_out[u+1];
+            for (int a = _first_out[u]; a != last_out; ++a) {
+              int v = _target[a];
+              if (_res_cap[a] > 0 && tol.negative(_cost[a] + p - pi[v]) &&
+                  l > level[v]) level[v] = l;
+            }
+          }
+
+          // Modify potentials
+          double q = std::numeric_limits<double>::max();
+          for (int u = 0; u != _res_node_num; ++u) {
+            int lu = level[u];
+            double p, pu = pi[u];
+            int last_out = _first_out[u+1];
+            for (int a = _first_out[u]; a != last_out; ++a) {
+              if (_res_cap[a] == 0) continue;
+              int v = _target[a];
+              int ld = lu - level[v];
+              if (ld > 0) {
+                p = (_cost[a] + pu - pi[v] + epsilon) / (ld + 1);
+                if (p < q) q = p;
+              }
+            }
+          }
+          for (int u = 0; u != _res_node_num; ++u) {
+            pi[u] -= q * level[u];
+          }
+
+          // Modify epsilon
+          epsilon = 0;
+          for (int u = 0; u != _res_node_num; ++u) {
+            double curr, pu = pi[u];
+            int last_out = _first_out[u+1];
+            for (int a = _first_out[u]; a != last_out; ++a) {
+              if (_res_cap[a] == 0) continue;
+              curr = _cost[a] + pu - pi[_target[a]];
+              if (-curr > epsilon) epsilon = -curr;
+            }
+          }
+        } else {
+          typedef Howard<StaticDigraph, CostArcMap> MMC;
+          typedef typename BellmanFord<StaticDigraph, CostArcMap>
+            ::template SetDistMap<CostNodeMap>::Create BF;
+
+          // Set epsilon to the minimum cycle mean
+          buildResidualNetwork();
+          MMC mmc(_sgr, _cost_map);
+          mmc.findMinMean();
+          epsilon = -mmc.cycleMean();
+          Cost cycle_cost = mmc.cycleLength();
+          int cycle_size = mmc.cycleArcNum();
+          
+          // Compute feasible potentials for the current epsilon
+          for (int i = 0; i != int(_cost_vec.size()); ++i) {
+            _cost_vec[i] = cycle_size * _cost_vec[i] - cycle_cost;
+          }
+          BF bf(_sgr, _cost_map);
+          bf.distMap(_pi_map);
+          bf.init(0);
+          bf.start();
+          for (int u = 0; u != _res_node_num; ++u) {
+            pi[u] = static_cast<double>(_pi[u]) / cycle_size;
+          }
+        
+          iter = limit;
+        }
+      }
+    }
+
+  }; //class CycleCanceling
+
+  ///@}
+
+} //namespace lemon
+
+#endif //LEMON_CYCLE_CANCELING_H
Index: lemon/network_simplex.h
===================================================================
--- lemon/network_simplex.h	(revision 835)
+++ lemon/network_simplex.h	(revision 878)
@@ -44,11 +44,11 @@
   /// \ref amo93networkflows, \ref dantzig63linearprog,
   /// \ref kellyoneill91netsimplex.
-  /// This algorithm is a specialized version of the linear programming
-  /// simplex method directly for the minimum cost flow problem.
-  /// It is one of the most efficient solution methods.
+  /// This algorithm is a highly efficient specialized version of the
+  /// linear programming simplex method directly for the minimum cost
+  /// flow problem.
   ///
-  /// In general this class is the fastest implementation available
-  /// in LEMON for the minimum cost flow problem.
-  /// Moreover it supports both directions of the supply/demand inequality
+  /// In general, %NetworkSimplex is the fastest implementation available
+  /// in LEMON for this problem.
+  /// Moreover, it supports both directions of the supply/demand inequality
   /// constraints. For more information, see \ref SupplyType.
   ///
@@ -59,10 +59,10 @@
   ///
   /// \tparam GR The digraph type the algorithm runs on.
-  /// \tparam V The value type used for flow amounts, capacity bounds
+  /// \tparam V The number 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
+  /// \tparam C The number type used for costs and potentials in the
   /// algorithm. By default, it is the same as \c V.
   ///
-  /// \warning Both value types must be signed and all input data must
+  /// \warning Both number types must be signed and all input data must
   /// be integer.
   ///
@@ -127,5 +127,5 @@
     /// 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.
+    /// test inputs.
     /// However, another pivot rule can be selected using the \ref run()
     /// function with the proper parameter.
@@ -165,5 +165,5 @@
 
     typedef std::vector<int> IntVector;
-    typedef std::vector<bool> BoolVector;
+    typedef std::vector<char> CharVector;
     typedef std::vector<Value> ValueVector;
     typedef std::vector<Cost> CostVector;
@@ -213,6 +213,6 @@
     IntVector _last_succ;
     IntVector _dirty_revs;
-    BoolVector _forward;
-    IntVector _state;
+    CharVector _forward;
+    CharVector _state;
     int _root;
 
@@ -222,4 +222,6 @@
     int stem, par_stem, new_stem;
     Value delta;
+    
+    const Value MAX;
 
   public:
@@ -243,5 +245,5 @@
       const IntVector  &_target;
       const CostVector &_cost;
-      const IntVector  &_state;
+      const CharVector &_state;
       const CostVector &_pi;
       int &_in_arc;
@@ -295,5 +297,5 @@
       const IntVector  &_target;
       const CostVector &_cost;
-      const IntVector  &_state;
+      const CharVector &_state;
       const CostVector &_pi;
       int &_in_arc;
@@ -334,5 +336,5 @@
       const IntVector  &_target;
       const CostVector &_cost;
-      const IntVector  &_state;
+      const CharVector &_state;
       const CostVector &_pi;
       int &_in_arc;
@@ -407,5 +409,5 @@
       const IntVector  &_target;
       const CostVector &_cost;
-      const IntVector  &_state;
+      const CharVector &_state;
       const CostVector &_pi;
       int &_in_arc;
@@ -510,5 +512,5 @@
       const IntVector  &_target;
       const CostVector &_cost;
-      const IntVector  &_state;
+      const CharVector &_state;
       const CostVector &_pi;
       int &_in_arc;
@@ -632,9 +634,9 @@
     NetworkSimplex(const GR& graph, bool arc_mixing = false) :
       _graph(graph), _node_id(graph), _arc_id(graph),
+      MAX(std::numeric_limits<Value>::max()),
       INF(std::numeric_limits<Value>::has_infinity ?
-          std::numeric_limits<Value>::infinity() :
-          std::numeric_limits<Value>::max())
+          std::numeric_limits<Value>::infinity() : MAX)
     {
-      // Check the value types
+      // Check the number types
       LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
         "The flow type of NetworkSimplex must be signed");
@@ -728,5 +730,5 @@
     /// 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).
+    /// unbounded from above).
     ///
     /// \param map An arc map storing the upper bounds.
@@ -1021,7 +1023,7 @@
           Value c = _lower[i];
           if (c >= 0) {
-            _cap[i] = _upper[i] < INF ? _upper[i] - c : INF;
+            _cap[i] = _upper[i] < MAX ? _upper[i] - c : INF;
           } else {
-            _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF;
+            _cap[i] = _upper[i] < MAX + c ? _upper[i] - c : INF;
           }
           _supply[_source[i]] -= c;
@@ -1215,5 +1217,5 @@
         e = _pred[u];
         d = _forward[u] ?
-          _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]);
+          _flow[e] : (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]);
         if (d < delta) {
           delta = d;
@@ -1226,5 +1228,5 @@
         e = _pred[u];
         d = _forward[u] ? 
-          (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e];
+          (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]) : _flow[e];
         if (d <= delta) {
           delta = d;
@@ -1425,5 +1427,5 @@
         findJoinNode();
         bool change = findLeavingArc();
-        if (delta >= INF) return UNBOUNDED;
+        if (delta >= MAX) return UNBOUNDED;
         changeFlow(change);
         if (change) {
Index: test/min_cost_flow_test.cc
===================================================================
--- test/min_cost_flow_test.cc	(revision 716)
+++ test/min_cost_flow_test.cc	(revision 885)
@@ -25,6 +25,10 @@
 
 #include <lemon/network_simplex.h>
+#include <lemon/capacity_scaling.h>
+#include <lemon/cost_scaling.h>
+#include <lemon/cycle_canceling.h>
 
 #include <lemon/concepts/digraph.h>
+#include <lemon/concepts/heap.h>
 #include <lemon/concept_check.h>
 
@@ -33,4 +37,5 @@
 using namespace lemon;
 
+// Test networks
 char test_lgf[] =
   "@nodes\n"
@@ -77,4 +82,56 @@
   "target 12\n";
 
+char test_neg1_lgf[] =
+  "@nodes\n"
+  "label   sup\n"
+  "    1   100\n"
+  "    2     0\n"
+  "    3     0\n"
+  "    4  -100\n"
+  "    5     0\n"
+  "    6     0\n"
+  "    7     0\n"
+  "@arcs\n"
+  "      cost   low1   low2\n"
+  "1 2    100      0      0\n"
+  "1 3     30      0      0\n"
+  "2 4     20      0      0\n"
+  "3 4     80      0      0\n"
+  "3 2     50      0      0\n"
+  "5 3     10      0      0\n"
+  "5 6     80      0   1000\n"
+  "6 7     30      0  -1000\n"
+  "7 5   -120      0      0\n";
+  
+char test_neg2_lgf[] =
+  "@nodes\n"
+  "label   sup\n"
+  "    1   100\n"
+  "    2  -300\n"
+  "@arcs\n"
+  "      cost\n"
+  "1 2     -1\n";
+
+
+// Test data
+typedef ListDigraph Digraph;
+DIGRAPH_TYPEDEFS(ListDigraph);
+
+Digraph gr;
+Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);
+Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);
+ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
+Node v, w;
+
+Digraph neg1_gr;
+Digraph::ArcMap<int> neg1_c(neg1_gr), neg1_l1(neg1_gr), neg1_l2(neg1_gr);
+ConstMap<Arc, int> neg1_u1(std::numeric_limits<int>::max()), neg1_u2(5000);
+Digraph::NodeMap<int> neg1_s(neg1_gr);
+
+Digraph neg2_gr;
+Digraph::ArcMap<int> neg2_c(neg2_gr);
+ConstMap<Arc, int> neg2_l(0), neg2_u(1000);
+Digraph::NodeMap<int> neg2_s(neg2_gr);
+
 
 enum SupplyType {
@@ -83,4 +140,5 @@
   LEQ
 };
+
 
 // Check the interface of an MCF algorithm
@@ -269,28 +327,97 @@
 }
 
+template < typename MCF, typename Param >
+void runMcfGeqTests( Param param,
+                     const std::string &test_str = "",
+                     bool full_neg_cost_support = false )
+{
+  MCF mcf1(gr), mcf2(neg1_gr), mcf3(neg2_gr);
+  
+  // Basic tests
+  mcf1.upperMap(u).costMap(c).supplyMap(s1);
+  checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s1,
+           mcf1.OPTIMAL, true,     5240, test_str + "-1");
+  mcf1.stSupply(v, w, 27);
+  checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s2,
+           mcf1.OPTIMAL, true,     7620, test_str + "-2");
+  mcf1.lowerMap(l2).supplyMap(s1);
+  checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s1,
+           mcf1.OPTIMAL, true,     5970, test_str + "-3");
+  mcf1.stSupply(v, w, 27);
+  checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s2,
+           mcf1.OPTIMAL, true,     8010, test_str + "-4");
+  mcf1.reset().supplyMap(s1);
+  checkMcf(mcf1, mcf1.run(param), gr, l1, cu, cc, s1,
+           mcf1.OPTIMAL, true,       74, test_str + "-5");
+  mcf1.lowerMap(l2).stSupply(v, w, 27);
+  checkMcf(mcf1, mcf1.run(param), gr, l2, cu, cc, s2,
+           mcf1.OPTIMAL, true,       94, test_str + "-6");
+  mcf1.reset();
+  checkMcf(mcf1, mcf1.run(param), gr, l1, cu, cc, s3,
+           mcf1.OPTIMAL, true,        0, test_str + "-7");
+  mcf1.lowerMap(l2).upperMap(u);
+  checkMcf(mcf1, mcf1.run(param), gr, l2, u, cc, s3,
+           mcf1.INFEASIBLE, false,    0, test_str + "-8");
+  mcf1.lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4);
+  checkMcf(mcf1, mcf1.run(param), gr, l3, u, c, s4,
+           mcf1.OPTIMAL, true,     6360, test_str + "-9");
+
+  // Tests for the GEQ form
+  mcf1.reset().upperMap(u).costMap(c).supplyMap(s5);
+  checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s5,
+           mcf1.OPTIMAL, true,     3530, test_str + "-10", GEQ);
+  mcf1.lowerMap(l2);
+  checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s5,
+           mcf1.OPTIMAL, true,     4540, test_str + "-11", GEQ);
+  mcf1.supplyMap(s6);
+  checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s6,
+           mcf1.INFEASIBLE, false,    0, test_str + "-12", GEQ);
+
+  // Tests with negative costs
+  mcf2.lowerMap(neg1_l1).costMap(neg1_c).supplyMap(neg1_s);
+  checkMcf(mcf2, mcf2.run(param), neg1_gr, neg1_l1, neg1_u1, neg1_c, neg1_s,
+           mcf2.UNBOUNDED, false,     0, test_str + "-13");
+  mcf2.upperMap(neg1_u2);
+  checkMcf(mcf2, mcf2.run(param), neg1_gr, neg1_l1, neg1_u2, neg1_c, neg1_s,
+           mcf2.OPTIMAL, true,   -40000, test_str + "-14");
+  mcf2.reset().lowerMap(neg1_l2).costMap(neg1_c).supplyMap(neg1_s);
+  checkMcf(mcf2, mcf2.run(param), neg1_gr, neg1_l2, neg1_u1, neg1_c, neg1_s,
+           mcf2.UNBOUNDED, false,     0, test_str + "-15");
+
+  mcf3.costMap(neg2_c).supplyMap(neg2_s);
+  if (full_neg_cost_support) {
+    checkMcf(mcf3, mcf3.run(param), neg2_gr, neg2_l, neg2_u, neg2_c, neg2_s,
+             mcf3.OPTIMAL, true,   -300, test_str + "-16", GEQ);
+  } else {
+    checkMcf(mcf3, mcf3.run(param), neg2_gr, neg2_l, neg2_u, neg2_c, neg2_s,
+             mcf3.UNBOUNDED, false,   0, test_str + "-17", GEQ);
+  }
+  mcf3.upperMap(neg2_u);
+  checkMcf(mcf3, mcf3.run(param), neg2_gr, neg2_l, neg2_u, neg2_c, neg2_s,
+           mcf3.OPTIMAL, true,     -300, test_str + "-18", GEQ);
+}
+
+template < typename MCF, typename Param >
+void runMcfLeqTests( Param param,
+                     const std::string &test_str = "" )
+{
+  // Tests for the LEQ form
+  MCF mcf1(gr);
+  mcf1.supplyType(mcf1.LEQ);
+  mcf1.upperMap(u).costMap(c).supplyMap(s6);
+  checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s6,
+           mcf1.OPTIMAL, true,   5080, test_str + "-19", LEQ);
+  mcf1.lowerMap(l2);
+  checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s6,
+           mcf1.OPTIMAL, true,   5930, test_str + "-20", LEQ);
+  mcf1.supplyMap(s5);
+  checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s5,
+           mcf1.INFEASIBLE, false,  0, test_str + "-21", LEQ);
+}
+
+
 int main()
 {
-  // Check the interfaces
-  {
-    typedef concepts::Digraph GR;
-    checkConcept< McfClassConcept<GR, int, int>,
-                  NetworkSimplex<GR> >();
-    checkConcept< McfClassConcept<GR, double, double>,
-                  NetworkSimplex<GR, double> >();
-    checkConcept< McfClassConcept<GR, int, double>,
-                  NetworkSimplex<GR, int, double> >();
-  }
-
-  // Run various MCF tests
-  typedef ListDigraph Digraph;
-  DIGRAPH_TYPEDEFS(ListDigraph);
-
-  // Read the test digraph
-  Digraph gr;
-  Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);
-  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);
-  ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
-  Node v, w;
-
+  // Read the test networks
   std::istringstream input(test_lgf);
   DigraphReader<Digraph>(gr, input)
@@ -310,139 +437,104 @@
     .run();
   
-  // Build test digraphs with negative costs
-  Digraph neg_gr;
-  Node n1 = neg_gr.addNode();
-  Node n2 = neg_gr.addNode();
-  Node n3 = neg_gr.addNode();
-  Node n4 = neg_gr.addNode();
-  Node n5 = neg_gr.addNode();
-  Node n6 = neg_gr.addNode();
-  Node n7 = neg_gr.addNode();
-  
-  Arc a1 = neg_gr.addArc(n1, n2);
-  Arc a2 = neg_gr.addArc(n1, n3);
-  Arc a3 = neg_gr.addArc(n2, n4);
-  Arc a4 = neg_gr.addArc(n3, n4);
-  Arc a5 = neg_gr.addArc(n3, n2);
-  Arc a6 = neg_gr.addArc(n5, n3);
-  Arc a7 = neg_gr.addArc(n5, n6);
-  Arc a8 = neg_gr.addArc(n6, n7);
-  Arc a9 = neg_gr.addArc(n7, n5);
-  
-  Digraph::ArcMap<int> neg_c(neg_gr), neg_l1(neg_gr, 0), neg_l2(neg_gr, 0);
-  ConstMap<Arc, int> neg_u1(std::numeric_limits<int>::max()), neg_u2(5000);
-  Digraph::NodeMap<int> neg_s(neg_gr, 0);
-  
-  neg_l2[a7] =  1000;
-  neg_l2[a8] = -1000;
-  
-  neg_s[n1] =  100;
-  neg_s[n4] = -100;
-  
-  neg_c[a1] =  100;
-  neg_c[a2] =   30;
-  neg_c[a3] =   20;
-  neg_c[a4] =   80;
-  neg_c[a5] =   50;
-  neg_c[a6] =   10;
-  neg_c[a7] =   80;
-  neg_c[a8] =   30;
-  neg_c[a9] = -120;
-
-  Digraph negs_gr;
-  Digraph::NodeMap<int> negs_s(negs_gr);
-  Digraph::ArcMap<int> negs_c(negs_gr);
-  ConstMap<Arc, int> negs_l(0), negs_u(1000);
-  n1 = negs_gr.addNode();
-  n2 = negs_gr.addNode();
-  negs_s[n1] = 100;
-  negs_s[n2] = -300;
-  negs_c[negs_gr.addArc(n1, n2)] = -1;
-
-
-  // A. Test NetworkSimplex with the default pivot rule
-  {
-    NetworkSimplex<Digraph> mcf(gr);
-
-    // Check the equality form
-    mcf.upperMap(u).costMap(c);
-    checkMcf(mcf, mcf.supplyMap(s1).run(),
-             gr, l1, u, c, s1, mcf.OPTIMAL, true,   5240, "#A1");
-    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
-             gr, l1, u, c, s2, mcf.OPTIMAL, true,   7620, "#A2");
-    mcf.lowerMap(l2);
-    checkMcf(mcf, mcf.supplyMap(s1).run(),
-             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#A3");
-    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
-             gr, l2, u, c, s2, mcf.OPTIMAL, true,   8010, "#A4");
-    mcf.reset();
-    checkMcf(mcf, mcf.supplyMap(s1).run(),
-             gr, l1, cu, cc, s1, mcf.OPTIMAL, true,   74, "#A5");
-    checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
-             gr, l2, cu, cc, s2, mcf.OPTIMAL, true,   94, "#A6");
-    mcf.reset();
-    checkMcf(mcf, mcf.run(),
-             gr, l1, cu, cc, s3, mcf.OPTIMAL, true,    0, "#A7");
-    checkMcf(mcf, mcf.lowerMap(l2).upperMap(u).run(),
-             gr, l2, u, cc, s3, mcf.INFEASIBLE, false, 0, "#A8");
-    mcf.reset().lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4);
-    checkMcf(mcf, mcf.run(),
-             gr, l3, u, c, s4, mcf.OPTIMAL, true,   6360, "#A9");
-
-    // Check the GEQ form
-    mcf.reset().upperMap(u).costMap(c).supplyMap(s5);
-    checkMcf(mcf, mcf.run(),
-             gr, l1, u, c, s5, mcf.OPTIMAL, true,   3530, "#A10", GEQ);
-    mcf.supplyType(mcf.GEQ);
-    checkMcf(mcf, mcf.lowerMap(l2).run(),
-             gr, l2, u, c, s5, mcf.OPTIMAL, true,   4540, "#A11", GEQ);
-    mcf.supplyMap(s6);
-    checkMcf(mcf, mcf.run(),
-             gr, l2, u, c, s6, mcf.INFEASIBLE, false,  0, "#A12", GEQ);
-
-    // Check the LEQ form
-    mcf.reset().supplyType(mcf.LEQ);
-    mcf.upperMap(u).costMap(c).supplyMap(s6);
-    checkMcf(mcf, mcf.run(),
-             gr, l1, u, c, s6, mcf.OPTIMAL, true,   5080, "#A13", LEQ);
-    checkMcf(mcf, mcf.lowerMap(l2).run(),
-             gr, l2, u, c, s6, mcf.OPTIMAL, true,   5930, "#A14", LEQ);
-    mcf.supplyMap(s5);
-    checkMcf(mcf, mcf.run(),
-             gr, l2, u, c, s5, mcf.INFEASIBLE, false,  0, "#A15", LEQ);
-
-    // Check negative costs
-    NetworkSimplex<Digraph> neg_mcf(neg_gr);
-    neg_mcf.lowerMap(neg_l1).costMap(neg_c).supplyMap(neg_s);
-    checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l1, neg_u1,
-      neg_c, neg_s, neg_mcf.UNBOUNDED, false,    0, "#A16");
-    neg_mcf.upperMap(neg_u2);
-    checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l1, neg_u2,
-      neg_c, neg_s, neg_mcf.OPTIMAL, true,  -40000, "#A17");
-    neg_mcf.reset().lowerMap(neg_l2).costMap(neg_c).supplyMap(neg_s);
-    checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l2, neg_u1,
-      neg_c, neg_s, neg_mcf.UNBOUNDED, false,    0, "#A18");
-      
-    NetworkSimplex<Digraph> negs_mcf(negs_gr);
-    negs_mcf.costMap(negs_c).supplyMap(negs_s);
-    checkMcf(negs_mcf, negs_mcf.run(), negs_gr, negs_l, negs_u,
-      negs_c, negs_s, negs_mcf.OPTIMAL, true, -300, "#A19", GEQ);
-  }
-
-  // B. Test NetworkSimplex with each pivot rule
-  {
-    NetworkSimplex<Digraph> mcf(gr);
-    mcf.supplyMap(s1).costMap(c).upperMap(u).lowerMap(l2);
-
-    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
-             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B1");
-    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
-             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B2");
-    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
-             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B3");
-    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
-             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B4");
-    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
-             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B5");
+  std::istringstream neg_inp1(test_neg1_lgf);
+  DigraphReader<Digraph>(neg1_gr, neg_inp1)
+    .arcMap("cost", neg1_c)
+    .arcMap("low1", neg1_l1)
+    .arcMap("low2", neg1_l2)
+    .nodeMap("sup", neg1_s)
+    .run();
+  
+  std::istringstream neg_inp2(test_neg2_lgf);
+  DigraphReader<Digraph>(neg2_gr, neg_inp2)
+    .arcMap("cost", neg2_c)
+    .nodeMap("sup", neg2_s)
+    .run();
+  
+  // Check the interface of NetworkSimplex
+  {
+    typedef concepts::Digraph GR;
+    checkConcept< McfClassConcept<GR, int, int>,
+                  NetworkSimplex<GR> >();
+    checkConcept< McfClassConcept<GR, double, double>,
+                  NetworkSimplex<GR, double> >();
+    checkConcept< McfClassConcept<GR, int, double>,
+                  NetworkSimplex<GR, int, double> >();
+  }
+
+  // Check the interface of CapacityScaling
+  {
+    typedef concepts::Digraph GR;
+    checkConcept< McfClassConcept<GR, int, int>,
+                  CapacityScaling<GR> >();
+    checkConcept< McfClassConcept<GR, double, double>,
+                  CapacityScaling<GR, double> >();
+    checkConcept< McfClassConcept<GR, int, double>,
+                  CapacityScaling<GR, int, double> >();
+    typedef CapacityScaling<GR>::
+      SetHeap<concepts::Heap<int, RangeMap<int> > >::Create CAS;
+    checkConcept< McfClassConcept<GR, int, int>, CAS >();
+  }
+
+  // Check the interface of CostScaling
+  {
+    typedef concepts::Digraph GR;
+    checkConcept< McfClassConcept<GR, int, int>,
+                  CostScaling<GR> >();
+    checkConcept< McfClassConcept<GR, double, double>,
+                  CostScaling<GR, double> >();
+    checkConcept< McfClassConcept<GR, int, double>,
+                  CostScaling<GR, int, double> >();
+    typedef CostScaling<GR>::
+      SetLargeCost<double>::Create COS;
+    checkConcept< McfClassConcept<GR, int, int>, COS >();
+  }
+
+  // Check the interface of CycleCanceling
+  {
+    typedef concepts::Digraph GR;
+    checkConcept< McfClassConcept<GR, int, int>,
+                  CycleCanceling<GR> >();
+    checkConcept< McfClassConcept<GR, double, double>,
+                  CycleCanceling<GR, double> >();
+    checkConcept< McfClassConcept<GR, int, double>,
+                  CycleCanceling<GR, int, double> >();
+  }
+
+  // Test NetworkSimplex
+  { 
+    typedef NetworkSimplex<Digraph> MCF;
+    runMcfGeqTests<MCF>(MCF::FIRST_ELIGIBLE, "NS-FE", true);
+    runMcfLeqTests<MCF>(MCF::FIRST_ELIGIBLE, "NS-FE");
+    runMcfGeqTests<MCF>(MCF::BEST_ELIGIBLE,  "NS-BE", true);
+    runMcfLeqTests<MCF>(MCF::BEST_ELIGIBLE,  "NS-BE");
+    runMcfGeqTests<MCF>(MCF::BLOCK_SEARCH,   "NS-BS", true);
+    runMcfLeqTests<MCF>(MCF::BLOCK_SEARCH,   "NS-BS");
+    runMcfGeqTests<MCF>(MCF::CANDIDATE_LIST, "NS-CL", true);
+    runMcfLeqTests<MCF>(MCF::CANDIDATE_LIST, "NS-CL");
+    runMcfGeqTests<MCF>(MCF::ALTERING_LIST,  "NS-AL", true);
+    runMcfLeqTests<MCF>(MCF::ALTERING_LIST,  "NS-AL");
+  }
+  
+  // Test CapacityScaling
+  {
+    typedef CapacityScaling<Digraph> MCF;
+    runMcfGeqTests<MCF>(0, "SSP");
+    runMcfGeqTests<MCF>(2, "CAS");
+  }
+
+  // Test CostScaling
+  {
+    typedef CostScaling<Digraph> MCF;
+    runMcfGeqTests<MCF>(MCF::PUSH, "COS-PR");
+    runMcfGeqTests<MCF>(MCF::AUGMENT, "COS-AR");
+    runMcfGeqTests<MCF>(MCF::PARTIAL_AUGMENT, "COS-PAR");
+  }
+
+  // Test CycleCanceling
+  {
+    typedef CycleCanceling<Digraph> MCF;
+    runMcfGeqTests<MCF>(MCF::SIMPLE_CYCLE_CANCELING, "SCC");
+    runMcfGeqTests<MCF>(MCF::MINIMUM_MEAN_CYCLE_CANCELING, "MMCC");
+    runMcfGeqTests<MCF>(MCF::CANCEL_AND_TIGHTEN, "CAT");
   }
 
