diff -r c160bf9f18ef -r 580af8cf2f6a lemon/network_simplex.h
--- a/lemon/network_simplex.h	Thu Nov 05 10:01:02 2009 +0100
+++ b/lemon/network_simplex.h	Thu Nov 05 10:23:16 2009 +0100
@@ -40,7 +40,9 @@
   /// for finding a \ref min_cost_flow "minimum cost flow".
   ///
   /// \ref NetworkSimplex implements the primal Network Simplex algorithm
-  /// for finding a \ref min_cost_flow "minimum cost flow".
+  /// for finding a \ref min_cost_flow "minimum cost flow"
+  /// \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.
@@ -161,8 +163,6 @@
 
     TEMPLATE_DIGRAPH_TYPEDEFS(GR);
 
-    typedef std::vector<Arc> ArcVector;
-    typedef std::vector<Node> NodeVector;
     typedef std::vector<int> IntVector;
     typedef std::vector<bool> BoolVector;
     typedef std::vector<Value> ValueVector;
@@ -364,33 +364,32 @@
       bool findEnteringArc() {
         Cost c, min = 0;
         int cnt = _block_size;
-        int e, min_arc = _next_arc;
+        int e;
         for (e = _next_arc; e < _search_arc_num; ++e) {
           c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
           if (c < min) {
             min = c;
-            min_arc = e;
+            _in_arc = e;
           }
           if (--cnt == 0) {
-            if (min < 0) break;
+            if (min < 0) goto search_end;
             cnt = _block_size;
           }
         }
-        if (min == 0 || cnt > 0) {
-          for (e = 0; e < _next_arc; ++e) {
-            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
-            if (c < min) {
-              min = c;
-              min_arc = e;
-            }
-            if (--cnt == 0) {
-              if (min < 0) break;
-              cnt = _block_size;
-            }
+        for (e = 0; e < _next_arc; ++e) {
+          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          if (c < min) {
+            min = c;
+            _in_arc = e;
+          }
+          if (--cnt == 0) {
+            if (min < 0) goto search_end;
+            cnt = _block_size;
           }
         }
         if (min >= 0) return false;
-        _in_arc = min_arc;
+
+      search_end:
         _next_arc = e;
         return true;
       }
@@ -428,7 +427,7 @@
         _next_arc(0)
       {
         // The main parameters of the pivot rule
-        const double LIST_LENGTH_FACTOR = 1.0;
+        const double LIST_LENGTH_FACTOR = 0.25;
         const int MIN_LIST_LENGTH = 10;
         const double MINOR_LIMIT_FACTOR = 0.1;
         const int MIN_MINOR_LIMIT = 3;
@@ -445,7 +444,7 @@
       /// Find next entering arc
       bool findEnteringArc() {
         Cost min, c;
-        int e, min_arc = _next_arc;
+        int e;
         if (_curr_length > 0 && _minor_count < _minor_limit) {
           // Minor iteration: select the best eligible arc from the
           // current candidate list
@@ -456,16 +455,13 @@
             c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
             if (c < min) {
               min = c;
-              min_arc = e;
+              _in_arc = e;
             }
-            if (c >= 0) {
+            else if (c >= 0) {
               _candidates[i--] = _candidates[--_curr_length];
             }
           }
-          if (min < 0) {
-            _in_arc = min_arc;
-            return true;
-          }
+          if (min < 0) return true;
         }
 
         // Major iteration: build a new candidate list
@@ -477,27 +473,26 @@
             _candidates[_curr_length++] = e;
             if (c < min) {
               min = c;
-              min_arc = e;
+              _in_arc = e;
             }
-            if (_curr_length == _list_length) break;
+            if (_curr_length == _list_length) goto search_end;
           }
         }
-        if (_curr_length < _list_length) {
-          for (e = 0; e < _next_arc; ++e) {
-            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
-            if (c < 0) {
-              _candidates[_curr_length++] = e;
-              if (c < min) {
-                min = c;
-                min_arc = e;
-              }
-              if (_curr_length == _list_length) break;
+        for (e = 0; e < _next_arc; ++e) {
+          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          if (c < 0) {
+            _candidates[_curr_length++] = e;
+            if (c < min) {
+              min = c;
+              _in_arc = e;
             }
+            if (_curr_length == _list_length) goto search_end;
           }
         }
         if (_curr_length == 0) return false;
+      
+      search_end:        
         _minor_count = 1;
-        _in_arc = min_arc;
         _next_arc = e;
         return true;
       }
@@ -549,7 +544,7 @@
         _next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost)
       {
         // The main parameters of the pivot rule
-        const double BLOCK_SIZE_FACTOR = 1.5;
+        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 = 3;
@@ -578,39 +573,35 @@
 
         // Extend the list
         int cnt = _block_size;
-        int last_arc = 0;
         int limit = _head_length;
 
-        for (int e = _next_arc; e < _search_arc_num; ++e) {
+        for (e = _next_arc; e < _search_arc_num; ++e) {
           _cand_cost[e] = _state[e] *
             (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
           if (_cand_cost[e] < 0) {
             _candidates[_curr_length++] = e;
-            last_arc = e;
           }
           if (--cnt == 0) {
-            if (_curr_length > limit) break;
+            if (_curr_length > limit) goto search_end;
             limit = 0;
             cnt = _block_size;
           }
         }
-        if (_curr_length <= limit) {
-          for (int e = 0; e < _next_arc; ++e) {
-            _cand_cost[e] = _state[e] *
-              (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
-            if (_cand_cost[e] < 0) {
-              _candidates[_curr_length++] = e;
-              last_arc = e;
-            }
-            if (--cnt == 0) {
-              if (_curr_length > limit) break;
-              limit = 0;
-              cnt = _block_size;
-            }
+        for (e = 0; e < _next_arc; ++e) {
+          _cand_cost[e] = _state[e] *
+            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          if (_cand_cost[e] < 0) {
+            _candidates[_curr_length++] = e;
+          }
+          if (--cnt == 0) {
+            if (_curr_length > limit) goto search_end;
+            limit = 0;
+            cnt = _block_size;
           }
         }
         if (_curr_length == 0) return false;
-        _next_arc = last_arc + 1;
+        
+      search_end:
 
         // Make heap of the candidate list (approximating a partial sort)
         make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
@@ -618,6 +609,7 @@
 
         // Pop the first element of the heap
         _in_arc = _candidates[0];
+        _next_arc = e;
         pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
                   _sort_func );
         _curr_length = std::min(_head_length, _curr_length - 1);
@@ -633,7 +625,11 @@
     /// The constructor of the class.
     ///
     /// \param graph The digraph the algorithm runs on.
-    NetworkSimplex(const GR& graph) :
+    /// \param arc_mixing Indicate if the arcs have to be stored in a
+    /// mixed order in the internal data structure. 
+    /// In special cases, it could lead to better overall performance,
+    /// but it is usually slower. Therefore it is disabled by default.
+    NetworkSimplex(const GR& graph, bool arc_mixing = false) :
       _graph(graph), _node_id(graph), _arc_id(graph),
       INF(std::numeric_limits<Value>::has_infinity ?
           std::numeric_limits<Value>::infinity() :
@@ -671,31 +667,33 @@
       _last_succ.resize(all_node_num);
       _state.resize(max_arc_num);
 
-      // Copy the graph (store the arcs in a mixed order)
+      // Copy the graph
       int i = 0;
       for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
         _node_id[n] = i;
       }
-      int k = std::max(int(std::sqrt(double(_arc_num))), 10);
-      i = 0;
-      for (ArcIt a(_graph); a != INVALID; ++a) {
-        _arc_id[a] = i;
-        _source[i] = _node_id[_graph.source(a)];
-        _target[i] = _node_id[_graph.target(a)];
-        if ((i += k) >= _arc_num) i = (i % k) + 1;
+      if (arc_mixing) {
+        // Store the arcs in a mixed order
+        int k = std::max(int(std::sqrt(double(_arc_num))), 10);
+        int i = 0, j = 0;
+        for (ArcIt a(_graph); a != INVALID; ++a) {
+          _arc_id[a] = i;
+          _source[i] = _node_id[_graph.source(a)];
+          _target[i] = _node_id[_graph.target(a)];
+          if ((i += k) >= _arc_num) i = ++j;
+        }
+      } else {
+        // Store the arcs in the original order
+        int i = 0;
+        for (ArcIt a(_graph); a != INVALID; ++a, ++i) {
+          _arc_id[a] = i;
+          _source[i] = _node_id[_graph.source(a)];
+          _target[i] = _node_id[_graph.target(a)];
+        }
       }
       
-      // Initialize maps
-      for (int i = 0; i != _node_num; ++i) {
-        _supply[i] = 0;
-      }
-      for (int i = 0; i != _arc_num; ++i) {
-        _lower[i] = 0;
-        _upper[i] = INF;
-        _cost[i] = 1;
-      }
-      _have_lower = false;
-      _stype = GEQ;
+      // Reset parameters
+      reset();
     }
 
     /// \name Parameters
@@ -768,7 +766,6 @@
     /// 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.
-    /// (It makes sense only if non-zero lower bounds are given.)
     ///
     /// \param map A node map storing the supply values.
     /// Its \c Value type must be convertible to the \c Value type
@@ -789,7 +786,6 @@
     /// 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.
-    /// (It makes sense only if non-zero lower bounds are given.)
     ///
     /// Using this function has the same effect as using \ref supplyMap()
     /// with such a map in which \c k is assigned to \c s, \c -k is