Index: lemon/network_simplex.h
===================================================================
--- lemon/network_simplex.h	(revision 2629)
+++ lemon/network_simplex.h	(revision 2630)
@@ -27,4 +27,5 @@
 #include <vector>
 #include <limits>
+#include <algorithm>
 
 #include <lemon/graph_adaptor.h>
@@ -138,5 +139,5 @@
       Cost operator[](const Edge &e) const {
         return _cost_map[e] + _pot_map[_gr.source(e)]
-                           - _pot_map[_gr.target(e)];
+                            - _pot_map[_gr.target(e)];
       }
 
@@ -169,5 +170,5 @@
 
       /// Find next entering edge
-      inline bool findEnteringEdge() {
+      bool findEnteringEdge() {
         Edge e;
         for (int i = _next_edge; i < int(_edges.size()); ++i) {
@@ -213,5 +214,5 @@
 
       /// Find next entering edge
-      inline bool findEnteringEdge() {
+      bool findEnteringEdge() {
         Cost min = 0;
         Edge e;
@@ -260,5 +261,5 @@
 
       /// Find next entering edge
-      inline bool findEnteringEdge() {
+      bool findEnteringEdge() {
         Cost curr, min = 0;
         Edge e;
@@ -337,9 +338,9 @@
 
       /// Find next entering edge
-      inline bool findEnteringEdge() {
+      bool findEnteringEdge() {
         Cost min, curr;
         if (_curr_length > 0 && _minor_count < _minor_limit) {
-          // Minor iteration: selecting the best eligible edge from
-          // the current candidate list
+          // Minor iteration: select the best eligible edge from the
+          // current candidate list
           ++_minor_count;
           Edge e;
@@ -359,5 +360,5 @@
         }
 
-        // Major iteration: building a new candidate list
+        // Major iteration: build a new candidate list
         Edge e;
         min = 0;
@@ -424,5 +425,5 @@
         SortFunc(const SCostMap &map) : _map(map) {}
         bool operator()(const Edge &e1, const Edge &e2) {
-	  return _map[e1] < _map[e2];
+          return _map[e1] > _map[e2];
         }
       };
@@ -438,8 +439,8 @@
       {
         // The main parameters of the pivot rule
-        const double BLOCK_SIZE_FACTOR = 1.0;
+        const double BLOCK_SIZE_FACTOR = 1.5;
         const int MIN_BLOCK_SIZE = 10;
         const double HEAD_LENGTH_FACTOR = 0.1;
-        const int MIN_HEAD_LENGTH = 5;
+        const int MIN_HEAD_LENGTH = 3;
 
         _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edges.size())),
@@ -452,15 +453,15 @@
 
       /// Find next entering edge
-      inline bool findEnteringEdge() {
-        // Checking the current candidate list
+      bool findEnteringEdge() {
+        // Check the current candidate list
         Edge e;
-        for (int idx = 0; idx < _curr_length; ++idx) {
-          e = _candidates[idx];
+        for (int ix = 0; ix < _curr_length; ++ix) {
+          e = _candidates[ix];
           if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) >= 0) {
-            _candidates[idx--] = _candidates[--_curr_length];
-          }
-        }
-
-        // Extending the list
+            _candidates[ix--] = _candidates[--_curr_length];
+          }
+        }
+
+        // Extend the list
         int cnt = _block_size;
         int last_edge = 0;
@@ -495,22 +496,13 @@
         _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);
-
+        // Make heap of the candidate list (approximating a partial sort)
+        make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
+                   _sort_func );
+
+        // Pop the first element of the heap
         _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;
-        }
-
+        pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
+                  _sort_func );
+        _curr_length = std::min(_head_length, _curr_length - 1);
         return true;
       }
@@ -559,5 +551,5 @@
     // The state edge map
     IntEdgeMap _state;
-    // The root node of the starting spanning tree
+    // The artificial root node of the spanning tree structure
     Node _root;
 
@@ -576,13 +568,13 @@
     EdgeRefMap _edge_ref;
 
-    // The entering edge of the current pivot iteration.
+    // The entering edge of the current pivot iteration
     Edge _in_edge;
 
-    // Temporary nodes used in the current pivot iteration.
+    // Temporary nodes used in the current pivot iteration
     Node join, u_in, v_in, u_out, v_out;
     Node right, first, second, last;
     Node stem, par_stem, new_stem;
     // The maximum augment amount along the found cycle in the current
-    // pivot iteration.
+    // pivot iteration
     Capacity delta;
 
@@ -838,6 +830,6 @@
     ///
     /// 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.
+    /// pivot rule proved to be 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.
@@ -912,10 +904,9 @@
   private:
 
-    /// \brief Extend the underlying graph and initialize all the
-    /// node and edge maps.
+    // Extend the underlying graph and initialize all the node and edge maps
     bool init() {
       if (!_valid_supply) return false;
 
-      // Initializing result maps
+      // Initialize result maps
       if (!_flow_result) {
         _flow_result = new FlowMap(_graph_ref);
@@ -927,5 +918,5 @@
       }
 
-      // Initializing the edge vector and the edge maps
+      // Initialize the edge vector and the edge maps
       _edges.reserve(countEdges(_graph));
       for (EdgeIt e(_graph); e != INVALID; ++e) {
@@ -935,5 +926,5 @@
       }
 
-      // Adding an artificial root node to the graph
+      // Add an artificial root node to the graph
       _root = _graph.addNode();
       _parent[_root] = INVALID;
@@ -943,6 +934,6 @@
       _potential[_root] = 0;
 
-      // Adding artificial edges to the graph and initializing the node
-      // maps of the spanning tree data structure
+      // Add artificial edges to the graph and initialize the node maps
+      // of the spanning tree data structure
       Node last = _root;
       Edge e;
@@ -975,6 +966,6 @@
     }
 
-    /// Find the join node.
-    inline Node findJoinNode() {
+    // Find the join node
+    void findJoinNode() {
       Node u = _graph.source(_in_edge);
       Node v = _graph.target(_in_edge);
@@ -987,12 +978,11 @@
         else v = _parent[v];
       }
-      return u;
-    }
-
-    /// \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
+      join = u;
+    }
+
+    // Find the leaving edge of the cycle and returns true if the 
+    // leaving edge is not the same as the entering edge
+    bool findLeavingEdge() {
+      // Initialize first and second nodes according to the direction
       // of the cycle
       if (_state[_in_edge] == STATE_LOWER) {
@@ -1008,6 +998,5 @@
       Edge e;
 
-      // Searching the cycle along the path form the first node to the
-      // root node
+      // Search the cycle along the path form the first node to the root
       for (Node u = first; u != join; u = _parent[u]) {
         e = _pred_edge[u];
@@ -1021,6 +1010,5 @@
         }
       }
-      // Searching the cycle along the path form the second node to the
-      // root node
+      // Search the cycle along the path form the second node to the root
       for (Node u = second; u != join; u = _parent[u]) {
         e = _pred_edge[u];
@@ -1037,7 +1025,7 @@
     }
 
-    /// Change \c flow and \c state edge maps.
-    inline void changeFlows(bool change) {
-      // Augmenting along the cycle
+    // Change _flow and _state edge maps
+    void changeFlows(bool change) {
+      // Augment along the cycle
       if (delta > 0) {
         Capacity val = _state[_in_edge] * delta;
@@ -1050,5 +1038,5 @@
         }
       }
-      // Updating the state of the entering and leaving edges
+      // Update the state of the entering and leaving edges
       if (change) {
         _state[_in_edge] = STATE_TREE;
@@ -1060,10 +1048,10 @@
     }
 
-    /// Update \c thread and \c parent node maps.
-    inline void updateThreadParent() {
+    // Update _thread and _parent node maps
+    void updateThreadParent() {
       Node u;
       v_out = _parent[u_out];
 
-      // Handling the case when join and v_out coincide
+      // Handle the case when join and v_out coincide
       bool par_first = false;
       if (join == v_out) {
@@ -1076,6 +1064,6 @@
       }
 
-      // Finding the last successor of u_in (u) and the node after it
-      // (right) according to the thread index
+      // Find the last successor of u_in (u) and the node after it (right)
+      // according to the thread index
       for (u = u_in; _depth[_thread[u]] > _depth[u_in]; u = _thread[u]) ;
       right = _thread[u];
@@ -1086,5 +1074,5 @@
       else last = _thread[v_in];
 
-      // Updating stem nodes
+      // Update stem nodes
       _thread[v_in] = stem = u_in;
       par_stem = v_in;
@@ -1092,17 +1080,16 @@
         _thread[u] = new_stem = _parent[stem];
 
-        // Finding the node just before the stem node (u) according to
+        // Find the node just before the stem node (u) according to
         // the original thread index
         for (u = new_stem; _thread[u] != stem; u = _thread[u]) ;
         _thread[u] = right;
 
-        // Changing the parent node of stem and shifting stem and
-        // par_stem nodes
+        // Change the parent node of stem and shift stem and par_stem nodes
         _parent[stem] = par_stem;
         par_stem = stem;
         stem = new_stem;
 
-        // Finding the last successor of stem (u) and the node after it
-        // (right) according to the thread index
+        // Find the last successor of stem (u) and the node after it (right)
+        // according to the thread index
         for (u = stem; _depth[_thread[u]] > _depth[stem]; u = _thread[u]) ;
         right = _thread[u];
@@ -1119,6 +1106,6 @@
     }
 
-    /// Update \c pred_edge and \c forward node maps.
-    inline void updatePredEdge() {
+    // Update _pred_edge and _forward node maps.
+    void updatePredEdge() {
       Node u = u_out, v;
       while (u != u_in) {
@@ -1132,6 +1119,6 @@
     }
 
-    /// Update \c depth and \c potential node maps.
-    inline void updateDepthPotential() {
+    // Update _depth and _potential node maps
+    void updateDepthPotential() {
       _depth[u_in] = _depth[v_in] + 1;
       Cost sigma = _forward[u_in] ?
@@ -1146,7 +1133,7 @@
     }
 
-    /// Execute the algorithm.
+    // Execute the algorithm
     bool start(PivotRuleEnum pivot_rule) {
-      // Selecting the pivot rule implementation
+      // Select the pivot rule implementation
       switch (pivot_rule) {
         case FIRST_ELIGIBLE_PIVOT:
@@ -1168,7 +1155,7 @@
       PivotRuleImplementation pivot(*this, _edges);
 
-      // Executing the network simplex algorithm
+      // Execute the network simplex algorithm
       while (pivot.findEnteringEdge()) {
-        join = findJoinNode();
+        findJoinNode();
         bool change = findLeavingEdge();
         changeFlows(change);
@@ -1180,6 +1167,5 @@
       }
 
-      // Checking if the flow amount equals zero on all the artificial
-      // edges
+      // Check if the flow amount equals zero on all the artificial edges
       for (InEdgeIt e(_graph, _root); e != INVALID; ++e)
         if (_flow[e] > 0) return false;
@@ -1187,5 +1173,5 @@
         if (_flow[e] > 0) return false;
 
-      // Copying flow values to _flow_result
+      // Copy flow values to _flow_result
       if (_lower) {
         for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
@@ -1195,5 +1181,5 @@
           (*_flow_result)[e] = _flow[_edge_ref[e]];
       }
-      // Copying potential values to _potential_result
+      // Copy potential values to _potential_result
       for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
         (*_potential_result)[n] = _potential[_node_ref[n]];
