diff -r 58c330ad0b5c -r 6be1f9bd2ac0 lemon/network_simplex.h
--- a/lemon/network_simplex.h	Sun Oct 04 10:15:32 2009 +0200
+++ b/lemon/network_simplex.h	Wed Dec 09 11:14:06 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.
@@ -48,7 +50,7 @@
   /// 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
-  /// constraints. For more information see \ref SupplyType.
+  /// constraints. For more information, see \ref SupplyType.
   ///
   /// Most of the parameters of the problem (except for the digraph)
   /// can be given using separate functions, and the algorithm can be
@@ -57,16 +59,16 @@
   ///
   /// \tparam GR The digraph type the algorithm runs on.
   /// \tparam V The value type used for flow amounts, capacity bounds
-  /// and supply values in the algorithm. By default it is \c int.
+  /// and supply values in the algorithm. By default, it is \c int.
   /// \tparam C The value type used for costs and potentials in the
-  /// algorithm. By default it is the same as \c V.
+  /// algorithm. By default, it is the same as \c V.
   ///
   /// \warning Both value types must be signed and all input data must
   /// be integer.
   ///
   /// \note %NetworkSimplex provides five different pivot rule
   /// implementations, from which the most efficient one is used
-  /// by default. For more information see \ref PivotRule.
+  /// by default. For more information, see \ref PivotRule.
   template <typename GR, typename V = int, typename C = V>
   class NetworkSimplex
   {
@@ -122,35 +124,35 @@
     /// \ref NetworkSimplex provides five different pivot rule
     /// implementations that significantly affect the running time
     /// of the algorithm.
-    /// By default \ref BLOCK_SEARCH "Block Search" is used, which
+    /// By default, \ref BLOCK_SEARCH "Block Search" is used, which
     /// proved to be the most efficient and the most robust on various
     /// test inputs according to our benchmark tests.
-    /// However another pivot rule can be selected using the \ref run()
+    /// However, another pivot rule can be selected using the \ref run()
     /// function with the proper parameter.
     enum PivotRule {
 
-      /// The First Eligible pivot rule.
+      /// The \e First \e Eligible pivot rule.
       /// The next eligible arc is selected in a wraparound fashion
       /// in every iteration.
       FIRST_ELIGIBLE,
 
-      /// The Best Eligible pivot rule.
+      /// The \e Best \e Eligible pivot rule.
       /// The best eligible arc is selected in every iteration.
       BEST_ELIGIBLE,
 
-      /// The Block Search pivot rule.
+      /// The \e Block \e Search pivot rule.
       /// A specified number of arcs are examined in every iteration
       /// in a wraparound fashion and the best eligible arc is selected
       /// from this block.
       BLOCK_SEARCH,
 
-      /// The Candidate List pivot rule.
+      /// The \e Candidate \e List pivot rule.
       /// In a major iteration a candidate list is built from eligible arcs
       /// in a wraparound fashion and in the following minor iterations
       /// the best eligible arc is selected from this list.
       CANDIDATE_LIST,
 
-      /// The Altering Candidate List pivot rule.
+      /// The \e Altering \e Candidate \e List pivot rule.
       /// It is a modified version of the Candidate List method.
       /// It keeps only the several best eligible arcs from the former
       /// candidate list and extends this list in every iteration.
@@ -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
@@ -816,7 +812,7 @@
     /// If it is not used before calling \ref run(), the \ref GEQ supply
     /// type will be used.
     ///
-    /// For more information see \ref SupplyType.
+    /// For more information, see \ref SupplyType.
     ///
     /// \return <tt>(*this)</tt>
     NetworkSimplex& supplyType(SupplyType supply_type) {
@@ -848,11 +844,11 @@
     /// 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
+    /// However, the underlying digraph must not be modified after this
     /// class have been constructed, since it copies and extends the graph.
     ///
     /// \param pivot_rule The pivot rule that will be used during the
-    /// algorithm. For more information see \ref PivotRule.
+    /// algorithm. For more information, see \ref PivotRule.
     ///
     /// \return \c INFEASIBLE if no feasible flow exists,
     /// \n \c OPTIMAL if the problem has optimal solution
@@ -877,7 +873,7 @@
     /// 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
+    /// However, the underlying digraph must not be modified after this
     /// class have been constructed, since it copies and extends the graph.
     ///
     /// For example,