Index: lemon/network_simplex.h
===================================================================
--- lemon/network_simplex.h	(revision 2593)
+++ lemon/network_simplex.h	(revision 2619)
@@ -33,4 +33,6 @@
 #include <lemon/math.h>
 
+#define _DEBUG_
+
 namespace lemon {
 
@@ -38,9 +40,9 @@
   /// @{
 
-  /// \brief Implementation of the network simplex algorithm for
-  /// finding a minimum cost flow.
+  /// \brief Implementation of the primal network simplex algorithm
+  /// for finding a minimum cost flow.
   ///
-  /// \ref NetworkSimplex implements the network simplex algorithm for
-  /// finding a minimum cost flow.
+  /// \ref NetworkSimplex implements the primal network simplex algorithm
+  /// for finding a minimum cost flow.
   ///
   /// \tparam Graph The directed graph type the algorithm runs on.
@@ -56,14 +58,13 @@
   /// - \c CostMap::Value must be signed type.
   ///
-  /// \note \ref NetworkSimplex provides six different pivot rule
+  /// \note \ref NetworkSimplex provides five different pivot rule
   /// implementations that significantly affect the efficiency of the
   /// algorithm.
-  /// By default a combined pivot rule is used, which is the fastest
-  /// implementation according to our benchmark tests.
-  /// Another pivot rule can be selected using \ref run() function
-  /// with the proper parameter.
+  /// By default "Block Search" pivot rule is used, which proved to be
+  /// by far the most efficient according to our benchmark tests.
+  /// However another pivot rule can be selected using \ref run()
+  /// function with the proper parameter.
   ///
   /// \author Peter Kovacs
-
   template < typename Graph,
              typename LowerMap = typename Graph::template EdgeMap<int>,
@@ -90,7 +91,10 @@
     typedef typename SGraph::template NodeMap<Edge> EdgeNodeMap;
     typedef typename SGraph::template EdgeMap<int> IntEdgeMap;
+    typedef typename SGraph::template EdgeMap<bool> BoolEdgeMap;
 
     typedef typename Graph::template NodeMap<Node> NodeRefMap;
     typedef typename Graph::template EdgeMap<Edge> EdgeRefMap;
+
+    typedef std::vector<Edge> EdgeVector;
 
   public:
@@ -108,7 +112,6 @@
       BEST_ELIGIBLE_PIVOT,
       BLOCK_SEARCH_PIVOT,
-      LIMITED_SEARCH_PIVOT,
       CANDIDATE_LIST_PIVOT,
-      COMBINED_PIVOT
+      ALTERING_LIST_PIVOT
     };
 
@@ -149,30 +152,38 @@
     /// This class implements the "First Eligible" pivot rule
     /// for the \ref NetworkSimplex "network simplex" algorithm.
+    ///
+    /// For more information see \ref NetworkSimplex::run().
     class FirstEligiblePivotRule
     {
     private:
 
+      // References to the NetworkSimplex class
       NetworkSimplex &_ns;
-      EdgeIt _next_edge;
+      EdgeVector &_edges;
+
+      int _next_edge;
 
     public:
 
-      /// Constructor.
-      FirstEligiblePivotRule(NetworkSimplex &ns) :
-        _ns(ns), _next_edge(ns._graph) {}
-
-      /// Finds the next entering edge.
-      bool findEnteringEdge() {
-        for (EdgeIt e = _next_edge; e != INVALID; ++e) {
+      /// Constructor
+      FirstEligiblePivotRule(NetworkSimplex &ns, EdgeVector &edges) :
+        _ns(ns), _edges(edges), _next_edge(0) {}
+
+      /// Find next entering edge
+      inline bool findEnteringEdge() {
+        Edge e;
+        for (int i = _next_edge; i < int(_edges.size()); ++i) {
+          e = _edges[i];
           if (_ns._state[e] * _ns._red_cost[e] < 0) {
             _ns._in_edge = e;
-            _next_edge = ++e;
+            _next_edge = i + 1;
             return true;
           }
         }
-        for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
+        for (int i = 0; i < _next_edge; ++i) {
+          e = _edges[i];
           if (_ns._state[e] * _ns._red_cost[e] < 0) {
             _ns._in_edge = e;
-            _next_edge = ++e;
+            _next_edge = i + 1;
             return true;
           }
@@ -187,19 +198,26 @@
     /// This class implements the "Best Eligible" pivot rule
     /// for the \ref NetworkSimplex "network simplex" algorithm.
+    ///
+    /// For more information see \ref NetworkSimplex::run().
     class BestEligiblePivotRule
     {
     private:
 
+      // References to the NetworkSimplex class
       NetworkSimplex &_ns;
+      EdgeVector &_edges;
 
     public:
 
-      /// Constructor.
-      BestEligiblePivotRule(NetworkSimplex &ns) : _ns(ns) {}
-
-      /// Finds the next entering edge.
-      bool findEnteringEdge() {
+      /// Constructor
+      BestEligiblePivotRule(NetworkSimplex &ns, EdgeVector &edges) :
+        _ns(ns), _edges(edges) {}
+
+      /// Find next entering edge
+      inline bool findEnteringEdge() {
         Cost min = 0;
-        for (EdgeIt e(_ns._graph); e != INVALID; ++e) {
+        Edge e;
+        for (int i = 0; i < int(_edges.size()); ++i) {
+          e = _edges[i];
           if (_ns._state[e] * _ns._red_cost[e] < min) {
             min = _ns._state[e] * _ns._red_cost[e];
@@ -216,109 +234,67 @@
     /// This class implements the "Block Search" pivot rule
     /// for the \ref NetworkSimplex "network simplex" algorithm.
+    ///
+    /// For more information see \ref NetworkSimplex::run().
     class BlockSearchPivotRule
     {
     private:
 
+      // References to the NetworkSimplex class
       NetworkSimplex &_ns;
-      EdgeIt _next_edge, _min_edge;
+      EdgeVector &_edges;
+
       int _block_size;
-
-      static const int MIN_BLOCK_SIZE = 10;
+      int _next_edge, _min_edge;
 
     public:
 
-      /// Constructor.
-      BlockSearchPivotRule(NetworkSimplex &ns) :
-        _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph)
+      /// Constructor
+      BlockSearchPivotRule(NetworkSimplex &ns, EdgeVector &edges) :
+        _ns(ns), _edges(edges), _next_edge(0), _min_edge(0)
       {
-        _block_size = 2 * int(sqrt(countEdges(_ns._graph)));
-        if (_block_size < MIN_BLOCK_SIZE) _block_size = MIN_BLOCK_SIZE;
-      }
-
-      /// Finds the next entering edge.
-      bool findEnteringEdge() {
+        // The main parameters of the pivot rule
+        const double BLOCK_SIZE_FACTOR = 2.0;
+        const int MIN_BLOCK_SIZE = 10;
+
+        _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edges.size())),
+                                MIN_BLOCK_SIZE );
+      }
+
+      /// Find next entering edge
+      inline bool findEnteringEdge() {
         Cost curr, min = 0;
-        int cnt = 0;
-        for (EdgeIt e = _next_edge; e != INVALID; ++e) {
+        Edge e;
+        int cnt = _block_size;
+        int i;
+        for (i = _next_edge; i < int(_edges.size()); ++i) {
+          e = _edges[i];
           if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
             min = curr;
-            _min_edge = e;
+            _min_edge = i;
           }
-          if (++cnt == _block_size) {
+          if (--cnt == 0) {
             if (min < 0) break;
-            cnt = 0;
+            cnt = _block_size;
           }
         }
-        if (min == 0) {
-          for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
+        if (min == 0 || cnt > 0) {
+          for (i = 0; i < _next_edge; ++i) {
+            e = _edges[i];
             if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
               min = curr;
-              _min_edge = e;
+              _min_edge = i;
             }
-            if (++cnt == _block_size) {
+            if (--cnt == 0) {
               if (min < 0) break;
-              cnt = 0;
+              cnt = _block_size;
             }
           }
         }
-        _ns._in_edge = _min_edge;
-        _next_edge = ++_min_edge;
-        return min < 0;
+        if (min >= 0) return false;
+        _ns._in_edge = _edges[_min_edge];
+        _next_edge = i;
+        return true;
       }
     }; //class BlockSearchPivotRule
-
-    /// \brief Implementation of the "Limited Search" pivot rule for the
-    /// \ref NetworkSimplex "network simplex" algorithm.
-    ///
-    /// This class implements the "Limited Search" pivot rule
-    /// for the \ref NetworkSimplex "network simplex" algorithm.
-    class LimitedSearchPivotRule
-    {
-    private:
-
-      NetworkSimplex &_ns;
-      EdgeIt _next_edge, _min_edge;
-      int _sample_size;
-
-      static const int SAMPLE_SIZE_FACTOR = 15;
-      static const int MIN_SAMPLE_SIZE = 10;
-
-    public:
-
-      /// Constructor.
-      LimitedSearchPivotRule(NetworkSimplex &ns) :
-        _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph)
-      {
-        _sample_size = countEdges(_ns._graph) *
-                       SAMPLE_SIZE_FACTOR / 10000;
-        if (_sample_size < MIN_SAMPLE_SIZE)
-          _sample_size = MIN_SAMPLE_SIZE;
-      }
-
-      /// Finds the next entering edge.
-      bool findEnteringEdge() {
-        Cost curr, min = 0;
-        int cnt = 0;
-        for (EdgeIt e = _next_edge; e != INVALID; ++e) {
-          if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
-            min = curr;
-            _min_edge = e;
-          }
-          if (curr < 0 && ++cnt == _sample_size) break;
-        }
-        if (min == 0) {
-          for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
-            if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
-              min = curr;
-              _min_edge = e;
-            }
-            if (curr < 0 && ++cnt == _sample_size) break;
-          }
-        }
-        _ns._in_edge = _min_edge;
-        _next_edge = ++_min_edge;
-        return min < 0;
-      }
-    }; //class LimitedSearchPivotRule
 
     /// \brief Implementation of the "Candidate List" pivot rule for the
@@ -327,92 +303,219 @@
     /// This class implements the "Candidate List" pivot rule
     /// for the \ref NetworkSimplex "network simplex" algorithm.
+    ///
+    /// For more information see \ref NetworkSimplex::run().
     class CandidateListPivotRule
     {
     private:
 
+      // References to the NetworkSimplex class
       NetworkSimplex &_ns;
-
-      // The list of candidate edges.
-      std::vector<Edge> _candidates;
-      // The maximum length of the edge list.
-      int _list_length;
-      // The maximum number of minor iterations between two major
-      // itarations.
-      int _minor_limit;
-
-      int _minor_count;
-      EdgeIt _next_edge;
-
-      static const int LIST_LENGTH_FACTOR = 20;
-      static const int MINOR_LIMIT_FACTOR = 10;
-      static const int MIN_LIST_LENGTH = 10;
-      static const int MIN_MINOR_LIMIT = 2;
+      EdgeVector &_edges;
+
+      EdgeVector _candidates;
+      int _list_length, _minor_limit;
+      int _curr_length, _minor_count;
+      int _next_edge, _min_edge;
 
     public:
 
-      /// Constructor.
-      CandidateListPivotRule(NetworkSimplex &ns) :
-        _ns(ns), _next_edge(ns._graph)
+      /// Constructor
+      CandidateListPivotRule(NetworkSimplex &ns, EdgeVector &edges) :
+        _ns(ns), _edges(edges), _next_edge(0), _min_edge(0)
       {
-        int edge_num = countEdges(_ns._graph);
-        _minor_count = 0;
-        _list_length = edge_num * LIST_LENGTH_FACTOR / 10000;
-        if (_list_length < MIN_LIST_LENGTH)
-          _list_length = MIN_LIST_LENGTH;
-        _minor_limit = _list_length * MINOR_LIMIT_FACTOR / 100;
-        if (_minor_limit < MIN_MINOR_LIMIT)
-          _minor_limit = MIN_MINOR_LIMIT;
-      }
-
-      /// Finds the next entering edge.
-      bool findEnteringEdge() {
+        // The main parameters of the pivot rule
+        const double LIST_LENGTH_FACTOR = 1.0;
+        const int MIN_LIST_LENGTH = 10;
+        const double MINOR_LIMIT_FACTOR = 0.1;
+        const int MIN_MINOR_LIMIT = 3;
+
+        _list_length = std::max( int(LIST_LENGTH_FACTOR * sqrt(_edges.size())),
+                                 MIN_LIST_LENGTH );
+        _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
+                                 MIN_MINOR_LIMIT );
+        _curr_length = _minor_count = 0;
+        _candidates.resize(_list_length);
+      }
+
+      /// Find next entering edge
+      inline bool findEnteringEdge() {
         Cost min, curr;
-        if (_minor_count < _minor_limit && _candidates.size() > 0) {
-          // Minor iteration
+        if (_curr_length > 0 && _minor_count < _minor_limit) {
+          // Minor iteration: selecting the best eligible edge from
+          // the current candidate list
           ++_minor_count;
           Edge e;
           min = 0;
-          for (int i = 0; i < int(_candidates.size()); ++i) {
+          for (int i = 0; i < _curr_length; ++i) {
             e = _candidates[i];
-            if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
-              min = curr;
-              _ns._in_edge = e;
-            }
-          }
-          if (min < 0) return true;
-        }
-
-        // Major iteration
-        _candidates.clear();
-        EdgeIt e = _next_edge;
-        min = 0;
-        for ( ; e != INVALID; ++e) {
-          if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
-            _candidates.push_back(e);
+            curr = _ns._state[e] * _ns._red_cost[e];
             if (curr < min) {
               min = curr;
               _ns._in_edge = e;
             }
-            if (int(_candidates.size()) == _list_length) break;
+            if (curr >= 0) {
+              _candidates[i--] = _candidates[--_curr_length];
+            }
           }
-        }
-        if (int(_candidates.size()) < _list_length) {
-          for (e = EdgeIt(_ns._graph); e != _next_edge; ++e) {
+          if (min < 0) return true;
+        }
+
+        // Major iteration: building a new candidate list
+        Edge e;
+        min = 0;
+        _curr_length = 0;
+        int i;
+        for (i = _next_edge; i < int(_edges.size()); ++i) {
+          e = _edges[i];
+          if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
+            _candidates[_curr_length++] = e;
+            if (curr < min) {
+              min = curr;
+              _min_edge = i;
+            }
+            if (_curr_length == _list_length) break;
+          }
+        }
+        if (_curr_length < _list_length) {
+          for (i = 0; i < _next_edge; ++i) {
+            e = _edges[i];
             if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
-              _candidates.push_back(e);
+              _candidates[_curr_length++] = e;
               if (curr < min) {
                 min = curr;
-                _ns._in_edge = e;
+                _min_edge = i;
               }
-              if (int(_candidates.size()) == _list_length) break;
+              if (_curr_length == _list_length) break;
             }
           }
         }
-        if (_candidates.size() == 0) return false;
+        if (_curr_length == 0) return false;
         _minor_count = 1;
-        _next_edge = ++e;
+        _ns._in_edge = _edges[_min_edge];
+        _next_edge = i;
         return true;
       }
     }; //class CandidateListPivotRule
+
+    /// \brief Implementation of the "Altering Candidate List" pivot rule
+    /// for the \ref NetworkSimplex "network simplex" algorithm.
+    ///
+    /// This class implements the "Altering Candidate List" pivot rule
+    /// for the \ref NetworkSimplex "network simplex" algorithm.
+    ///
+    /// For more information see \ref NetworkSimplex::run().
+    class AlteringListPivotRule
+    {
+    private:
+
+      // References to the NetworkSimplex class
+      NetworkSimplex &_ns;
+      EdgeVector &_edges;
+
+      EdgeVector _candidates;
+      SCostMap _cand_cost;
+      int _block_size, _head_length, _curr_length;
+      int _next_edge;
+
+      // Functor class to compare edges during sort of the candidate list
+      class SortFunc
+      {
+      private:
+        const SCostMap &_map;
+      public:
+        SortFunc(const SCostMap &map) : _map(map) {}
+        bool operator()(const Edge &e1, const Edge &e2) {
+	  return _map[e1] < _map[e2];
+        }
+      };
+
+      SortFunc _sort_func;
+
+    public:
+
+      /// Constructor
+      AlteringListPivotRule(NetworkSimplex &ns, EdgeVector &edges) :
+        _ns(ns), _edges(edges), _cand_cost(_ns._graph),
+        _next_edge(0), _sort_func(_cand_cost)
+      {
+        // The main parameters of the pivot rule
+        const double BLOCK_SIZE_FACTOR = 1.0;
+        const int MIN_BLOCK_SIZE = 10;
+        const double HEAD_LENGTH_FACTOR = 0.1;
+        const int MIN_HEAD_LENGTH = 5;
+
+        _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edges.size())),
+                                MIN_BLOCK_SIZE );
+        _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
+                                 MIN_HEAD_LENGTH );
+        _candidates.resize(_head_length + _block_size);
+        _curr_length = 0;
+      }
+
+      /// Find next entering edge
+      inline bool findEnteringEdge() {
+        // Checking the current candidate list
+        Edge e;
+        for (int idx = 0; idx < _curr_length; ++idx) {
+          e = _candidates[idx];
+          if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) >= 0) {
+            _candidates[idx--] = _candidates[--_curr_length];
+          }
+        }
+
+        // Extending the list
+        int cnt = _block_size;
+        int last_edge = 0;
+        int limit = _head_length;
+        for (int i = _next_edge; i < int(_edges.size()); ++i) {
+          e = _edges[i];
+          if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) < 0) {
+            _candidates[_curr_length++] = e;
+            last_edge = i;
+          }
+          if (--cnt == 0) {
+            if (_curr_length > limit) break;
+            limit = 0;
+            cnt = _block_size;
+          }
+        }
+        if (_curr_length <= limit) {
+          for (int i = 0; i < _next_edge; ++i) {
+            e = _edges[i];
+            if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) < 0) {
+              _candidates[_curr_length++] = e;
+              last_edge = i;
+            }
+            if (--cnt == 0) {
+              if (_curr_length > limit) break;
+              limit = 0;
+              cnt = _block_size;
+            }
+          }
+        }
+        if (_curr_length == 0) return false;
+        _next_edge = last_edge + 1;
+
+        // Sorting the list partially
+        EdgeVector::iterator sort_end = _candidates.begin();
+        EdgeVector::iterator vector_end = _candidates.begin();
+        for (int idx = 0; idx < _curr_length; ++idx) {
+          ++vector_end;
+          if (idx <= _head_length) ++sort_end;
+        }
+        partial_sort(_candidates.begin(), sort_end, vector_end, _sort_func);
+
+        _ns._in_edge = _candidates[0];
+        if (_curr_length > _head_length) {
+          _candidates[0] = _candidates[_head_length - 1];
+          _curr_length = _head_length - 1;
+        } else {
+          _candidates[0] = _candidates[_curr_length - 1];
+          --_curr_length;
+        }
+
+        return true;
+      }
+    }; //class AlteringListPivotRule
 
   private:
@@ -424,7 +527,4 @@
       STATE_LOWER =  1
     };
-
-    // Constant for the combined pivot rule.
-    static const int COMBINED_PIVOT_MAX_DEG = 5;
 
   private:
@@ -466,4 +566,7 @@
     // The reduced cost map
     ReducedCostMap _red_cost;
+
+    // The non-artifical edges
+    EdgeVector _edges;
 
     // Members for handling the original graph
@@ -673,7 +776,7 @@
     }
 
-    /// \brief Sets the flow map.
-    ///
-    /// Sets the flow map.
+    /// \brief Set the flow map.
+    ///
+    /// Set the flow map.
     ///
     /// \return \c (*this)
@@ -687,7 +790,7 @@
     }
 
-    /// \brief Sets the potential map.
-    ///
-    /// Sets the potential map.
+    /// \brief Set the potential map.
+    ///
+    /// Set the potential map.
     ///
     /// \return \c (*this)
@@ -702,6 +805,4 @@
 
     /// \name Execution control
-    /// The only way to execute the algorithm is to call the run()
-    /// function.
 
     /// @{
@@ -727,22 +828,20 @@
     /// edge is selected from this block (\ref BlockSearchPivotRule).
     ///
-    /// - LIMITED_SEARCH_PIVOT A specified number of eligible edges are
-    /// examined in every iteration in a wraparound fashion and the best
-    /// one is selected from them (\ref LimitedSearchPivotRule).
-    ///
-    /// - CANDIDATE_LIST_PIVOT In major iterations a candidate list is
-    /// built from eligible edges and it is used for edge selection in
-    /// the following minor iterations (\ref CandidateListPivotRule).
-    ///
-    /// - COMBINED_PIVOT This is a combined version of the two fastest
-    /// pivot rules.
-    /// For rather sparse graphs \ref LimitedSearchPivotRule
-    /// "Limited Search" implementation is used, otherwise
-    /// \ref BlockSearchPivotRule "Block Search" pivot rule is used.
-    /// According to our benchmark tests this combined method is the
-    /// most efficient.
+    /// - CANDIDATE_LIST_PIVOT In a major iteration a candidate list is
+    /// built from eligible edges in a wraparound fashion and in the
+    /// following minor iterations the best eligible edge is selected
+    /// from this list (\ref CandidateListPivotRule).
+    ///
+    /// - ALTERING_LIST_PIVOT It is a modified version of the
+    /// "Candidate List" pivot rule. It keeps only the several best
+    /// eligible edges from the former candidate list and extends this
+    /// list in every iteration (\ref AlteringListPivotRule).
+    ///
+    /// According to our comprehensive benchmark tests the "Block Search"
+    /// pivot rule proved to be by far the fastest and the most robust
+    /// on various test inputs. Thus it is the default option.
     ///
     /// \return \c true if a feasible flow can be found.
-    bool run(PivotRuleEnum pivot_rule = COMBINED_PIVOT) {
+    bool run(PivotRuleEnum pivot_rule = BLOCK_SEARCH_PIVOT) {
       return init() && start(pivot_rule);
     }
@@ -751,14 +850,15 @@
 
     /// \name Query Functions
-    /// The result of the algorithm can be obtained using these
-    /// functions.
-    /// \n run() must be called before using them.
+    /// The results of the algorithm can be obtained using these
+    /// functions.\n
+    /// \ref lemon::NetworkSimplex::run() "run()" must be called before
+    /// using them.
 
     /// @{
 
-    /// \brief Returns a const reference to the edge map storing the
+    /// \brief Return a const reference to the edge map storing the
     /// found flow.
     ///
-    /// Returns a const reference to the edge map storing the found flow.
+    /// Return a const reference to the edge map storing the found flow.
     ///
     /// \pre \ref run() must be called before using this function.
@@ -767,8 +867,8 @@
     }
 
-    /// \brief Returns a const reference to the node map storing the
+    /// \brief Return a const reference to the node map storing the
     /// found potentials (the dual solution).
     ///
-    /// Returns a const reference to the node map storing the found
+    /// Return a const reference to the node map storing the found
     /// potentials (the dual solution).
     ///
@@ -778,7 +878,7 @@
     }
 
-    /// \brief Returns the flow on the given edge.
-    ///
-    /// Returns the flow on the given edge.
+    /// \brief Return the flow on the given edge.
+    ///
+    /// Return the flow on the given edge.
     ///
     /// \pre \ref run() must be called before using this function.
@@ -787,7 +887,7 @@
     }
 
-    /// \brief Returns the potential of the given node.
-    ///
-    /// Returns the potential of the given node.
+    /// \brief Return the potential of the given node.
+    ///
+    /// Return the potential of the given node.
     ///
     /// \pre \ref run() must be called before using this function.
@@ -796,7 +896,7 @@
     }
 
-    /// \brief Returns the total cost of the found flow.
-    ///
-    /// Returns the total cost of the found flow. The complexity of the
+    /// \brief Return the total cost of the found flow.
+    ///
+    /// Return the total cost of the found flow. The complexity of the
     /// function is \f$ O(e) \f$.
     ///
@@ -813,5 +913,5 @@
   private:
 
-    /// \brief Extends the underlaying graph and initializes all the
+    /// \brief Extend the underlying graph and initialize all the
     /// node and edge maps.
     bool init() {
@@ -828,6 +928,8 @@
       }
 
-      // Initializing state and flow maps
+      // Initializing the edge vector and the edge maps
+      _edges.reserve(countEdges(_graph));
       for (EdgeIt e(_graph); e != INVALID; ++e) {
+        _edges.push_back(e);
         _flow[e] = 0;
         _state[e] = STATE_LOWER;
@@ -874,6 +976,6 @@
     }
 
-    /// Finds the join node.
-    Node findJoinNode() {
+    /// Find the join node.
+    inline Node findJoinNode() {
       Node u = _graph.source(_in_edge);
       Node v = _graph.target(_in_edge);
@@ -889,7 +991,8 @@
     }
 
-    /// \brief Finds the leaving edge of the cycle. Returns \c true if
-    /// the leaving edge is not the same as the entering edge.
-    bool findLeavingEdge() {
+    /// \brief Find the leaving edge of the cycle.
+    /// \return \c true if the leaving edge is not the same as the
+    /// entering edge.
+    inline bool findLeavingEdge() {
       // Initializing first and second nodes according to the direction
       // of the cycle
@@ -935,6 +1038,6 @@
     }
 
-    /// Changes \c flow and \c state edge maps.
-    void changeFlows(bool change) {
+    /// Change \c flow and \c state edge maps.
+    inline void changeFlows(bool change) {
       // Augmenting along the cycle
       if (delta > 0) {
@@ -958,6 +1061,6 @@
     }
 
-    /// Updates \c thread and \c parent node maps.
-    void updateThreadParent() {
+    /// Update \c thread and \c parent node maps.
+    inline void updateThreadParent() {
       Node u;
       v_out = _parent[u_out];
@@ -1017,6 +1120,6 @@
     }
 
-    /// Updates \c pred_edge and \c forward node maps.
-    void updatePredEdge() {
+    /// Update \c pred_edge and \c forward node maps.
+    inline void updatePredEdge() {
       Node u = u_out, v;
       while (u != u_in) {
@@ -1030,6 +1133,6 @@
     }
 
-    /// Updates \c depth and \c potential node maps.
-    void updateDepthPotential() {
+    /// Update \c depth and \c potential node maps.
+    inline void updateDepthPotential() {
       _depth[u_in] = _depth[v_in] + 1;
       _potential[u_in] = _forward[u_in] ?
@@ -1050,6 +1153,7 @@
     }
 
-    /// Executes the algorithm.
+    /// Execute the algorithm.
     bool start(PivotRuleEnum pivot_rule) {
+      // Selecting the pivot rule implementation
       switch (pivot_rule) {
         case FIRST_ELIGIBLE_PIVOT:
@@ -1059,14 +1163,8 @@
         case BLOCK_SEARCH_PIVOT:
           return start<BlockSearchPivotRule>();
-        case LIMITED_SEARCH_PIVOT:
-          return start<LimitedSearchPivotRule>();
         case CANDIDATE_LIST_PIVOT:
           return start<CandidateListPivotRule>();
-        case COMBINED_PIVOT:
-          if ( countEdges(_graph) / countNodes(_graph) <=
-               COMBINED_PIVOT_MAX_DEG )
-            return start<LimitedSearchPivotRule>();
-          else
-            return start<BlockSearchPivotRule>();
+        case ALTERING_LIST_PIVOT:
+          return start<AlteringListPivotRule>();
       }
       return false;
@@ -1075,5 +1173,5 @@
     template<class PivotRuleImplementation>
     bool start() {
-      PivotRuleImplementation pivot(*this);
+      PivotRuleImplementation pivot(*this, _edges);
 
       // Executing the network simplex algorithm
