Index: lemon/Makefile.am
===================================================================
--- lemon/Makefile.am	(revision 953)
+++ lemon/Makefile.am	(revision 999)
@@ -91,4 +91,5 @@
 	lemon/graph_to_eps.h \
 	lemon/grid_graph.h \
+	lemon/grosso_locatelli_pullan_mc.h \
 	lemon/hartmann_orlin_mmc.h \
 	lemon/howard_mmc.h \
Index: lemon/grosso_locatelli_pullan_mc.h
===================================================================
--- lemon/grosso_locatelli_pullan_mc.h	(revision 999)
+++ lemon/grosso_locatelli_pullan_mc.h	(revision 999)
@@ -0,0 +1,680 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2010
+ * 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_GROSSO_LOCATELLI_PULLAN_MC_H
+#define LEMON_GROSSO_LOCATELLI_PULLAN_MC_H
+
+/// \ingroup approx_algs
+///
+/// \file
+/// \brief The iterated local search algorithm of Grosso, Locatelli, and Pullan
+/// for the maximum clique problem
+
+#include <vector>
+#include <limits>
+#include <lemon/core.h>
+#include <lemon/random.h>
+
+namespace lemon {
+
+  /// \addtogroup approx_algs
+  /// @{
+
+  /// \brief Implementation of the iterated local search algorithm of Grosso,
+  /// Locatelli, and Pullan for the maximum clique problem
+  ///
+  /// \ref GrossoLocatelliPullanMc implements the iterated local search
+  /// algorithm of Grosso, Locatelli, and Pullan for solving the \e maximum
+  /// \e clique \e problem \ref grosso08maxclique.
+  /// It is to find the largest complete subgraph (\e clique) in an
+  /// undirected graph, i.e., the largest set of nodes where each
+  /// pair of nodes is connected.
+  ///
+  /// This class provides a simple but highly efficient and robust heuristic
+  /// method that quickly finds a large clique, but not necessarily the
+  /// largest one.
+  ///
+  /// \tparam GR The undirected graph type the algorithm runs on.
+  ///
+  /// \note %GrossoLocatelliPullanMc provides three different node selection
+  /// rules, from which the most powerful one is used by default.
+  /// For more information, see \ref SelectionRule.
+  template <typename GR>
+  class GrossoLocatelliPullanMc
+  {
+  public:
+
+    /// \brief Constants for specifying the node selection rule.
+    ///
+    /// Enum type containing constants for specifying the node selection rule
+    /// for the \ref run() function.
+    ///
+    /// During the algorithm, nodes are selected for addition to the current
+    /// clique according to the applied rule.
+    /// In general, the PENALTY_BASED rule turned out to be the most powerful
+    /// and the most robust, thus it is the default option.
+    /// However, another selection rule can be specified using the \ref run()
+    /// function with the proper parameter.
+    enum SelectionRule {
+
+      /// A node is selected randomly without any evaluation at each step.
+      RANDOM,
+
+      /// A node of maximum degree is selected randomly at each step.
+      DEGREE_BASED,
+
+      /// A node of minimum penalty is selected randomly at each step.
+      /// The node penalties are updated adaptively after each stage of the
+      /// search process.
+      PENALTY_BASED
+    };
+
+  private:
+
+    TEMPLATE_GRAPH_TYPEDEFS(GR);
+
+    typedef std::vector<int> IntVector;
+    typedef std::vector<char> BoolVector;
+    typedef std::vector<BoolVector> BoolMatrix;
+    // Note: vector<char> is used instead of vector<bool> for efficiency reasons
+
+    const GR &_graph;
+    IntNodeMap _id;
+
+    // Internal matrix representation of the graph
+    BoolMatrix _gr;
+    int _n;
+
+    // The current clique
+    BoolVector _clique;
+    int _size;
+
+    // The best clique found so far
+    BoolVector _best_clique;
+    int _best_size;
+
+    // The "distances" of the nodes from the current clique.
+    // _delta[u] is the number of nodes in the clique that are
+    // not connected with u.
+    IntVector _delta;
+
+    // The current tabu set
+    BoolVector _tabu;
+
+    // Random number generator
+    Random _rnd;
+
+  private:
+
+    // Implementation of the RANDOM node selection rule.
+    class RandomSelectionRule
+    {
+    private:
+
+      // References to the algorithm instance
+      const BoolVector &_clique;
+      const IntVector  &_delta;
+      const BoolVector &_tabu;
+      Random &_rnd;
+
+      // Pivot rule data
+      int _n;
+
+    public:
+
+      // Constructor
+      RandomSelectionRule(GrossoLocatelliPullanMc &mc) :
+        _clique(mc._clique), _delta(mc._delta), _tabu(mc._tabu),
+        _rnd(mc._rnd), _n(mc._n)
+      {}
+
+      // Return a node index for a feasible add move or -1 if no one exists
+      int nextFeasibleAddNode() const {
+        int start_node = _rnd[_n];
+        for (int i = start_node; i != _n; i++) {
+          if (_delta[i] == 0 && !_tabu[i]) return i;
+        }
+        for (int i = 0; i != start_node; i++) {
+          if (_delta[i] == 0 && !_tabu[i]) return i;
+        }
+        return -1;
+      }
+
+      // Return a node index for a feasible swap move or -1 if no one exists
+      int nextFeasibleSwapNode() const {
+        int start_node = _rnd[_n];
+        for (int i = start_node; i != _n; i++) {
+          if (!_clique[i] && _delta[i] == 1 && !_tabu[i]) return i;
+        }
+        for (int i = 0; i != start_node; i++) {
+          if (!_clique[i] && _delta[i] == 1 && !_tabu[i]) return i;
+        }
+        return -1;
+      }
+
+      // Return a node index for an add move or -1 if no one exists
+      int nextAddNode() const {
+        int start_node = _rnd[_n];
+        for (int i = start_node; i != _n; i++) {
+          if (_delta[i] == 0) return i;
+        }
+        for (int i = 0; i != start_node; i++) {
+          if (_delta[i] == 0) return i;
+        }
+        return -1;
+      }
+
+      // Update internal data structures between stages (if necessary)
+      void update() {}
+
+    }; //class RandomSelectionRule
+
+
+    // Implementation of the DEGREE_BASED node selection rule.
+    class DegreeBasedSelectionRule
+    {
+    private:
+
+      // References to the algorithm instance
+      const BoolVector &_clique;
+      const IntVector  &_delta;
+      const BoolVector &_tabu;
+      Random &_rnd;
+
+      // Pivot rule data
+      int _n;
+      IntVector _deg;
+
+    public:
+
+      // Constructor
+      DegreeBasedSelectionRule(GrossoLocatelliPullanMc &mc) :
+        _clique(mc._clique), _delta(mc._delta), _tabu(mc._tabu),
+        _rnd(mc._rnd), _n(mc._n), _deg(_n)
+      {
+        for (int i = 0; i != _n; i++) {
+          int d = 0;
+          BoolVector &row = mc._gr[i];
+          for (int j = 0; j != _n; j++) {
+            if (row[j]) d++;
+          }
+          _deg[i] = d;
+        }
+      }
+
+      // Return a node index for a feasible add move or -1 if no one exists
+      int nextFeasibleAddNode() const {
+        int start_node = _rnd[_n];
+        int node = -1, max_deg = -1;
+        for (int i = start_node; i != _n; i++) {
+          if (_delta[i] == 0 && !_tabu[i] && _deg[i] > max_deg) {
+            node = i;
+            max_deg = _deg[i];
+          }
+        }
+        for (int i = 0; i != start_node; i++) {
+          if (_delta[i] == 0 && !_tabu[i] && _deg[i] > max_deg) {
+            node = i;
+            max_deg = _deg[i];
+          }
+        }
+        return node;
+      }
+
+      // Return a node index for a feasible swap move or -1 if no one exists
+      int nextFeasibleSwapNode() const {
+        int start_node = _rnd[_n];
+        int node = -1, max_deg = -1;
+        for (int i = start_node; i != _n; i++) {
+          if (!_clique[i] && _delta[i] == 1 && !_tabu[i] &&
+              _deg[i] > max_deg) {
+            node = i;
+            max_deg = _deg[i];
+          }
+        }
+        for (int i = 0; i != start_node; i++) {
+          if (!_clique[i] && _delta[i] == 1 && !_tabu[i] &&
+              _deg[i] > max_deg) {
+            node = i;
+            max_deg = _deg[i];
+          }
+        }
+        return node;
+      }
+
+      // Return a node index for an add move or -1 if no one exists
+      int nextAddNode() const {
+        int start_node = _rnd[_n];
+        int node = -1, max_deg = -1;
+        for (int i = start_node; i != _n; i++) {
+          if (_delta[i] == 0 && _deg[i] > max_deg) {
+            node = i;
+            max_deg = _deg[i];
+          }
+        }
+        for (int i = 0; i != start_node; i++) {
+          if (_delta[i] == 0 && _deg[i] > max_deg) {
+            node = i;
+            max_deg = _deg[i];
+          }
+        }
+        return node;
+      }
+
+      // Update internal data structures between stages (if necessary)
+      void update() {}
+
+    }; //class DegreeBasedSelectionRule
+
+
+    // Implementation of the PENALTY_BASED node selection rule.
+    class PenaltyBasedSelectionRule
+    {
+    private:
+
+      // References to the algorithm instance
+      const BoolVector &_clique;
+      const IntVector  &_delta;
+      const BoolVector &_tabu;
+      Random &_rnd;
+
+      // Pivot rule data
+      int _n;
+      IntVector _penalty;
+
+    public:
+
+      // Constructor
+      PenaltyBasedSelectionRule(GrossoLocatelliPullanMc &mc) :
+        _clique(mc._clique), _delta(mc._delta), _tabu(mc._tabu),
+        _rnd(mc._rnd), _n(mc._n), _penalty(_n, 0)
+      {}
+
+      // Return a node index for a feasible add move or -1 if no one exists
+      int nextFeasibleAddNode() const {
+        int start_node = _rnd[_n];
+        int node = -1, min_p = std::numeric_limits<int>::max();
+        for (int i = start_node; i != _n; i++) {
+          if (_delta[i] == 0 && !_tabu[i] && _penalty[i] < min_p) {
+            node = i;
+            min_p = _penalty[i];
+          }
+        }
+        for (int i = 0; i != start_node; i++) {
+          if (_delta[i] == 0 && !_tabu[i] && _penalty[i] < min_p) {
+            node = i;
+            min_p = _penalty[i];
+          }
+        }
+        return node;
+      }
+
+      // Return a node index for a feasible swap move or -1 if no one exists
+      int nextFeasibleSwapNode() const {
+        int start_node = _rnd[_n];
+        int node = -1, min_p = std::numeric_limits<int>::max();
+        for (int i = start_node; i != _n; i++) {
+          if (!_clique[i] && _delta[i] == 1 && !_tabu[i] &&
+              _penalty[i] < min_p) {
+            node = i;
+            min_p = _penalty[i];
+          }
+        }
+        for (int i = 0; i != start_node; i++) {
+          if (!_clique[i] && _delta[i] == 1 && !_tabu[i] &&
+              _penalty[i] < min_p) {
+            node = i;
+            min_p = _penalty[i];
+          }
+        }
+        return node;
+      }
+
+      // Return a node index for an add move or -1 if no one exists
+      int nextAddNode() const {
+        int start_node = _rnd[_n];
+        int node = -1, min_p = std::numeric_limits<int>::max();
+        for (int i = start_node; i != _n; i++) {
+          if (_delta[i] == 0 && _penalty[i] < min_p) {
+            node = i;
+            min_p = _penalty[i];
+          }
+        }
+        for (int i = 0; i != start_node; i++) {
+          if (_delta[i] == 0 && _penalty[i] < min_p) {
+            node = i;
+            min_p = _penalty[i];
+          }
+        }
+        return node;
+      }
+
+      // Update internal data structures between stages (if necessary)
+      void update() {}
+
+    }; //class PenaltyBasedSelectionRule
+
+  public:
+
+    /// \brief Constructor.
+    ///
+    /// Constructor.
+    /// The global \ref rnd "random number generator instance" is used
+    /// during the algorithm.
+    ///
+    /// \param graph The undirected graph the algorithm runs on.
+    GrossoLocatelliPullanMc(const GR& graph) :
+      _graph(graph), _id(_graph), _rnd(rnd)
+    {}
+
+    /// \brief Constructor with random seed.
+    ///
+    /// Constructor with random seed.
+    ///
+    /// \param graph The undirected graph the algorithm runs on.
+    /// \param seed Seed value for the internal random number generator
+    /// that is used during the algorithm.
+    GrossoLocatelliPullanMc(const GR& graph, int seed) :
+      _graph(graph), _id(_graph), _rnd(seed)
+    {}
+
+    /// \brief Constructor with random number generator.
+    ///
+    /// Constructor with random number generator.
+    ///
+    /// \param graph The undirected graph the algorithm runs on.
+    /// \param random A random number generator that is used during the
+    /// algorithm.
+    GrossoLocatelliPullanMc(const GR& graph, const Random& random) :
+      _graph(graph), _id(_graph), _rnd(random)
+    {}
+
+    /// \name Execution Control
+    /// @{
+
+    /// \brief Runs the algorithm.
+    ///
+    /// This function runs the algorithm.
+    ///
+    /// \param step_num The maximum number of node selections (steps)
+    /// during the search process.
+    /// This parameter controls the running time and the success of the
+    /// algorithm. For larger values, the algorithm runs slower but it more
+    /// likely finds larger cliques. For smaller values, the algorithm is
+    /// faster but probably gives worse results.
+    /// \param rule The node selection rule. For more information, see
+    /// \ref SelectionRule.
+    ///
+    /// \return The size of the found clique.
+    int run(int step_num = 100000,
+            SelectionRule rule = PENALTY_BASED)
+    {
+      init();
+      switch (rule) {
+        case RANDOM:
+          return start<RandomSelectionRule>(step_num);
+        case DEGREE_BASED:
+          return start<DegreeBasedSelectionRule>(step_num);
+        case PENALTY_BASED:
+          return start<PenaltyBasedSelectionRule>(step_num);
+      }
+      return 0; // avoid warning
+    }
+
+    /// @}
+
+    /// \name Query Functions
+    /// @{
+
+    /// \brief The size of the found clique
+    ///
+    /// This function returns the size of the found clique.
+    ///
+    /// \pre run() must be called before using this function.
+    int cliqueSize() const {
+      return _best_size;
+    }
+
+    /// \brief Gives back the found clique in a \c bool node map
+    ///
+    /// This function gives back the characteristic vector of the found
+    /// clique in the given node map.
+    /// It must be a \ref concepts::WriteMap "writable" node map with
+    /// \c bool (or convertible) value type.
+    ///
+    /// \pre run() must be called before using this function.
+    template <typename CliqueMap>
+    void cliqueMap(CliqueMap &map) const {
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        map[n] = static_cast<bool>(_best_clique[_id[n]]);
+      }
+    }
+
+    /// \brief Iterator to list the nodes of the found clique
+    ///
+    /// This iterator class lists the nodes of the found clique.
+    /// Before using it, you must allocate a GrossoLocatelliPullanMc instance
+    /// and call its \ref GrossoLocatelliPullanMc::run() "run()" method.
+    ///
+    /// The following example prints out the IDs of the nodes in the found
+    /// clique.
+    /// \code
+    ///   GrossoLocatelliPullanMc<Graph> mc(g);
+    ///   mc.run();
+    ///   for (GrossoLocatelliPullanMc<Graph>::CliqueNodeIt n(mc);
+    ///        n != INVALID; ++n)
+    ///   {
+    ///     std::cout << g.id(n) << std::endl;
+    ///   }
+    /// \endcode
+    class CliqueNodeIt
+    {
+    private:
+      NodeIt _it;
+      BoolNodeMap _map;
+
+    public:
+
+      /// Constructor
+
+      /// Constructor.
+      /// \param mc The algorithm instance.
+      CliqueNodeIt(const GrossoLocatelliPullanMc &mc)
+       : _map(mc._graph)
+      {
+        mc.cliqueMap(_map);
+        for (_it = NodeIt(mc._graph); _it != INVALID && !_map[_it]; ++_it) ;
+      }
+
+      /// Conversion to \c Node
+      operator Node() const { return _it; }
+
+      bool operator==(Invalid) const { return _it == INVALID; }
+      bool operator!=(Invalid) const { return _it != INVALID; }
+
+      /// Next node
+      CliqueNodeIt &operator++() {
+        for (++_it; _it != INVALID && !_map[_it]; ++_it) ;
+        return *this;
+      }
+
+      /// Postfix incrementation
+
+      /// Postfix incrementation.
+      ///
+      /// \warning This incrementation returns a \c Node, not a
+      /// \c CliqueNodeIt as one may expect.
+      typename GR::Node operator++(int) {
+        Node n=*this;
+        ++(*this);
+        return n;
+      }
+
+    };
+
+    /// @}
+
+  private:
+
+    // Adds a node to the current clique
+    void addCliqueNode(int u) {
+      if (_clique[u]) return;
+      _clique[u] = true;
+      _size++;
+      BoolVector &row = _gr[u];
+      for (int i = 0; i != _n; i++) {
+        if (!row[i]) _delta[i]++;
+      }
+    }
+
+    // Removes a node from the current clique
+    void delCliqueNode(int u) {
+      if (!_clique[u]) return;
+      _clique[u] = false;
+      _size--;
+      BoolVector &row = _gr[u];
+      for (int i = 0; i != _n; i++) {
+        if (!row[i]) _delta[i]--;
+      }
+    }
+
+    // Initialize data structures
+    void init() {
+      _n = countNodes(_graph);
+      int ui = 0;
+      for (NodeIt u(_graph); u != INVALID; ++u) {
+        _id[u] = ui++;
+      }
+      _gr.clear();
+      _gr.resize(_n, BoolVector(_n, false));
+      ui = 0;
+      for (NodeIt u(_graph); u != INVALID; ++u) {
+        for (IncEdgeIt e(_graph, u); e != INVALID; ++e) {
+          int vi = _id[_graph.runningNode(e)];
+          _gr[ui][vi] = true;
+          _gr[vi][ui] = true;
+        }
+        ++ui;
+      }
+
+      _clique.clear();
+      _clique.resize(_n, false);
+      _size = 0;
+      _best_clique.clear();
+      _best_clique.resize(_n, false);
+      _best_size = 0;
+      _delta.clear();
+      _delta.resize(_n, 0);
+      _tabu.clear();
+      _tabu.resize(_n, false);
+    }
+
+    // Executes the algorithm
+    template <typename SelectionRuleImpl>
+    int start(int max_select) {
+      // Options for the restart rule
+      const bool delta_based_restart = true;
+      const int restart_delta_limit = 4;
+
+      if (_n == 0) return 0;
+      if (_n == 1) {
+        _best_clique[0] = true;
+        _best_size = 1;
+        return _best_size;
+      }
+
+      // Iterated local search
+      SelectionRuleImpl sel_method(*this);
+      int select = 0;
+      IntVector restart_nodes;
+
+      while (select < max_select) {
+
+        // Perturbation/restart
+        if (delta_based_restart) {
+          restart_nodes.clear();
+          for (int i = 0; i != _n; i++) {
+            if (_delta[i] >= restart_delta_limit)
+              restart_nodes.push_back(i);
+          }
+        }
+        int rs_node = -1;
+        if (restart_nodes.size() > 0) {
+          rs_node = restart_nodes[_rnd[restart_nodes.size()]];
+        } else {
+          rs_node = _rnd[_n];
+        }
+        BoolVector &row = _gr[rs_node];
+        for (int i = 0; i != _n; i++) {
+          if (_clique[i] && !row[i]) delCliqueNode(i);
+        }
+        addCliqueNode(rs_node);
+
+        // Local search
+        _tabu.clear();
+        _tabu.resize(_n, false);
+        bool tabu_empty = true;
+        int max_swap = _size;
+        while (select < max_select) {
+          select++;
+          int u;
+          if ((u = sel_method.nextFeasibleAddNode()) != -1) {
+            // Feasible add move
+            addCliqueNode(u);
+            if (tabu_empty) max_swap = _size;
+          }
+          else if ((u = sel_method.nextFeasibleSwapNode()) != -1) {
+            // Feasible swap move
+            int v = -1;
+            BoolVector &row = _gr[u];
+            for (int i = 0; i != _n; i++) {
+              if (_clique[i] && !row[i]) {
+                v = i;
+                break;
+              }
+            }
+            addCliqueNode(u);
+            delCliqueNode(v);
+            _tabu[v] = true;
+            tabu_empty = false;
+            if (--max_swap <= 0) break;
+          }
+          else if ((u = sel_method.nextAddNode()) != -1) {
+            // Non-feasible add move
+            addCliqueNode(u);
+          }
+          else break;
+        }
+        if (_size > _best_size) {
+          _best_clique = _clique;
+          _best_size = _size;
+          if (_best_size == _n) return _best_size;
+        }
+        sel_method.update();
+      }
+
+      return _best_size;
+    }
+
+  }; //class GrossoLocatelliPullanMc
+
+  ///@}
+
+} //namespace lemon
+
+#endif //LEMON_GROSSO_LOCATELLI_PULLAN_MC_H
Index: lemon/network_simplex.h
===================================================================
--- lemon/network_simplex.h	(revision 991)
+++ lemon/network_simplex.h	(revision 978)
@@ -167,7 +167,6 @@
     typedef std::vector<Value> ValueVector;
     typedef std::vector<Cost> CostVector;
-    typedef std::vector<signed char> CharVector;
-    // Note: vector<signed char> is used instead of vector<ArcState> and 
-    // vector<ArcDirection> for efficiency reasons
+    typedef std::vector<char> BoolVector;
+    // Note: vector<char> is used instead of vector<bool> for efficiency reasons
 
     // State constants for arcs
@@ -178,9 +177,7 @@
     };
 
-    // Direction constants for tree arcs
-    enum ArcDirection {
-      DIR_DOWN = -1,
-      DIR_UP   =  1
-    };
+    typedef std::vector<signed char> StateVector;
+    // Note: vector<signed char> is used instead of vector<ArcState> for
+    // efficiency reasons
 
   private:
@@ -221,11 +218,13 @@
     IntVector _succ_num;
     IntVector _last_succ;
-    CharVector _pred_dir;
-    CharVector _state;
     IntVector _dirty_revs;
+    BoolVector _forward;
+    StateVector _state;
     int _root;
 
     // Temporary data used in the current pivot iteration
     int in_arc, join, u_in, v_in, u_out, v_out;
+    int first, second, right, last;
+    int stem, par_stem, new_stem;
     Value delta;
 
@@ -252,5 +251,5 @@
       const IntVector  &_target;
       const CostVector &_cost;
-      const CharVector &_state;
+      const StateVector &_state;
       const CostVector &_pi;
       int &_in_arc;
@@ -304,5 +303,5 @@
       const IntVector  &_target;
       const CostVector &_cost;
-      const CharVector &_state;
+      const StateVector &_state;
       const CostVector &_pi;
       int &_in_arc;
@@ -343,5 +342,5 @@
       const IntVector  &_target;
       const CostVector &_cost;
-      const CharVector &_state;
+      const StateVector &_state;
       const CostVector &_pi;
       int &_in_arc;
@@ -416,5 +415,5 @@
       const IntVector  &_target;
       const CostVector &_cost;
-      const CharVector &_state;
+      const StateVector &_state;
       const CostVector &_pi;
       int &_in_arc;
@@ -519,5 +518,5 @@
       const IntVector  &_target;
       const CostVector &_cost;
-      const CharVector &_state;
+      const StateVector &_state;
       const CostVector &_pi;
       int &_in_arc;
@@ -572,11 +571,9 @@
         // Check the current candidate list
         int e;
-        Cost c;
         for (int i = 0; i != _curr_length; ++i) {
           e = _candidates[i];
-          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
-          if (c < 0) {
-            _cand_cost[e] = c;
-          } else {
+          _cand_cost[e] = _state[e] *
+            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          if (_cand_cost[e] >= 0) {
             _candidates[i--] = _candidates[--_curr_length];
           }
@@ -588,7 +585,7 @@
 
         for (e = _next_arc; e != _search_arc_num; ++e) {
-          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
-          if (c < 0) {
-            _cand_cost[e] = c;
+          _cand_cost[e] = _state[e] *
+            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          if (_cand_cost[e] < 0) {
             _candidates[_curr_length++] = e;
           }
@@ -637,10 +634,9 @@
     ///
     /// \param graph The digraph the algorithm runs on.
-    /// \param arc_mixing Indicate if the arcs will be stored in a
+    /// \param arc_mixing Indicate if the arcs have to be stored in a
     /// mixed order in the internal data structure.
-    /// In general, it leads to similar performance as using the original
-    /// arc order, but it makes the algorithm more robust and in special
-    /// cases, even significantly faster. Therefore, it is enabled by default.
-    NetworkSimplex(const GR& graph, bool arc_mixing = true) :
+    /// 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),
       _arc_mixing(arc_mixing),
@@ -918,5 +914,5 @@
       _parent.resize(all_node_num);
       _pred.resize(all_node_num);
-      _pred_dir.resize(all_node_num);
+      _forward.resize(all_node_num);
       _thread.resize(all_node_num);
       _rev_thread.resize(all_node_num);
@@ -932,5 +928,5 @@
       if (_arc_mixing) {
         // Store the arcs in a mixed order
-        const int skip = std::max(_arc_num / _node_num, 3);
+        int k = std::max(int(std::sqrt(double(_arc_num))), 10);
         int i = 0, j = 0;
         for (ArcIt a(_graph); a != INVALID; ++a) {
@@ -938,5 +934,5 @@
           _source[i] = _node_id[_graph.source(a)];
           _target[i] = _node_id[_graph.target(a)];
-          if ((i += skip) >= _arc_num) i = ++j;
+          if ((i += k) >= _arc_num) i = ++j;
         }
       } else {
@@ -1121,5 +1117,5 @@
           _state[e] = STATE_TREE;
           if (_supply[u] >= 0) {
-            _pred_dir[u] = DIR_UP;
+            _forward[u] = true;
             _pi[u] = 0;
             _source[e] = u;
@@ -1128,5 +1124,5 @@
             _cost[e] = 0;
           } else {
-            _pred_dir[u] = DIR_DOWN;
+            _forward[u] = false;
             _pi[u] = ART_COST;
             _source[e] = _root;
@@ -1148,5 +1144,5 @@
           _last_succ[u] = u;
           if (_supply[u] >= 0) {
-            _pred_dir[u] = DIR_UP;
+            _forward[u] = true;
             _pi[u] = 0;
             _pred[u] = e;
@@ -1158,5 +1154,5 @@
             _state[e] = STATE_TREE;
           } else {
-            _pred_dir[u] = DIR_DOWN;
+            _forward[u] = false;
             _pi[u] = ART_COST;
             _pred[u] = f;
@@ -1189,5 +1185,5 @@
           _last_succ[u] = u;
           if (_supply[u] <= 0) {
-            _pred_dir[u] = DIR_DOWN;
+            _forward[u] = false;
             _pi[u] = 0;
             _pred[u] = e;
@@ -1199,5 +1195,5 @@
             _state[e] = STATE_TREE;
           } else {
-            _pred_dir[u] = DIR_UP;
+            _forward[u] = true;
             _pi[u] = -ART_COST;
             _pred[u] = f;
@@ -1242,5 +1238,4 @@
       // Initialize first and second nodes according to the direction
       // of the cycle
-      int first, second;
       if (_state[in_arc] == STATE_LOWER) {
         first  = _source[in_arc];
@@ -1252,15 +1247,12 @@
       delta = _cap[in_arc];
       int result = 0;
-      Value c, d;
+      Value d;
       int e;
 
-      // Search the cycle form the first node to the join node
+      // Search the cycle along the path form the first node to the root
       for (int u = first; u != join; u = _parent[u]) {
         e = _pred[u];
-        d = _flow[e];
-        if (_pred_dir[u] == DIR_DOWN) {
-          c = _cap[e];
-          d = c >= MAX ? INF : c - d;
-        }
+        d = _forward[u] ?
+          _flow[e] : (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]);
         if (d < delta) {
           delta = d;
@@ -1269,13 +1261,9 @@
         }
       }
-
-      // Search the cycle form the second node to the join node
+      // Search the cycle along the path form the second node to the root
       for (int u = second; u != join; u = _parent[u]) {
         e = _pred[u];
-        d = _flow[e];
-        if (_pred_dir[u] == DIR_UP) {
-          c = _cap[e];
-          d = c >= MAX ? INF : c - d;
-        }
+        d = _forward[u] ?
+          (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]) : _flow[e];
         if (d <= delta) {
           delta = d;
@@ -1302,8 +1290,8 @@
         _flow[in_arc] += val;
         for (int u = _source[in_arc]; u != join; u = _parent[u]) {
-          _flow[_pred[u]] -= _pred_dir[u] * val;
+          _flow[_pred[u]] += _forward[u] ? -val : val;
         }
         for (int u = _target[in_arc]; u != join; u = _parent[u]) {
-          _flow[_pred[u]] += _pred_dir[u] * val;
+          _flow[_pred[u]] += _forward[u] ? val : -val;
         }
       }
@@ -1320,4 +1308,5 @@
     // Update the tree structure
     void updateTreeStructure() {
+      int u, w;
       int old_rev_thread = _rev_thread[u_out];
       int old_succ_num = _succ_num[u_out];
@@ -1325,127 +1314,122 @@
       v_out = _parent[u_out];
 
-      // Check if u_in and u_out coincide
-      if (u_in == u_out) {
-        // Update _parent, _pred, _pred_dir
-        _parent[u_in] = v_in;
-        _pred[u_in] = in_arc;
-        _pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN;
-
-        // Update _thread and _rev_thread
-        if (_thread[v_in] != u_out) {
-          int after = _thread[old_last_succ];
-          _thread[old_rev_thread] = after;
-          _rev_thread[after] = old_rev_thread;
-          after = _thread[v_in];
-          _thread[v_in] = u_out;
-          _rev_thread[u_out] = v_in;
-          _thread[old_last_succ] = after;
-          _rev_thread[after] = old_last_succ;
-        }
+      u = _last_succ[u_in];  // the last successor of u_in
+      right = _thread[u];    // the node after it
+
+      // Handle the case when old_rev_thread equals to v_in
+      // (it also means that join and v_out coincide)
+      if (old_rev_thread == v_in) {
+        last = _thread[_last_succ[u_out]];
       } else {
-        // Handle the case when old_rev_thread equals to v_in
-        // (it also means that join and v_out coincide)
-        int thread_continue = old_rev_thread == v_in ?
-          _thread[old_last_succ] : _thread[v_in];
-
-        // Update _thread and _parent along the stem nodes (i.e. the nodes
-        // between u_in and u_out, whose parent have to be changed)
-        int stem = u_in;              // the current stem node
-        int par_stem = v_in;          // the new parent of stem
-        int next_stem;                // the next stem node
-        int last = _last_succ[u_in];  // the last successor of stem
-        int before, after = _thread[last];
-        _thread[v_in] = u_in;
-        _dirty_revs.clear();
-        _dirty_revs.push_back(v_in);
-        while (stem != u_out) {
-          // Insert the next stem node into the thread list
-          next_stem = _parent[stem];
-          _thread[last] = next_stem;
-          _dirty_revs.push_back(last);
-
-          // Remove the subtree of stem from the thread list
-          before = _rev_thread[stem];
-          _thread[before] = after;
-          _rev_thread[after] = before;
-
-          // Change the parent node and shift stem nodes
-          _parent[stem] = par_stem;
-          par_stem = stem;
-          stem = next_stem;
-
-          // Update last and after
-          last = _last_succ[stem] == _last_succ[par_stem] ?
-            _rev_thread[par_stem] : _last_succ[stem];
-          after = _thread[last];
-        }
-        _parent[u_out] = par_stem;
-        _thread[last] = thread_continue;
-        _rev_thread[thread_continue] = last;
-        _last_succ[u_out] = last;
-
-        // Remove the subtree of u_out from the thread list except for
-        // the case when old_rev_thread equals to v_in
-        if (old_rev_thread != v_in) {
-          _thread[old_rev_thread] = after;
-          _rev_thread[after] = old_rev_thread;
-        }
-
-        // Update _rev_thread using the new _thread values
-        for (int i = 0; i != int(_dirty_revs.size()); ++i) {
-          int u = _dirty_revs[i];
-          _rev_thread[_thread[u]] = u;
-        }
-
-        // Update _pred, _pred_dir, _last_succ and _succ_num for the
-        // stem nodes from u_out to u_in
-        int tmp_sc = 0, tmp_ls = _last_succ[u_out];
-        for (int u = u_out, p = _parent[u]; u != u_in; u = p, p = _parent[u]) {
-          _pred[u] = _pred[p];
-          _pred_dir[u] = -_pred_dir[p];
-          tmp_sc += _succ_num[u] - _succ_num[p];
-          _succ_num[u] = tmp_sc;
-          _last_succ[p] = tmp_ls;
-        }
-        _pred[u_in] = in_arc;
-        _pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN;
-        _succ_num[u_in] = old_succ_num;
+        last = _thread[v_in];
+      }
+
+      // Update _thread and _parent along the stem nodes (i.e. the nodes
+      // between u_in and u_out, whose parent have to be changed)
+      _thread[v_in] = stem = u_in;
+      _dirty_revs.clear();
+      _dirty_revs.push_back(v_in);
+      par_stem = v_in;
+      while (stem != u_out) {
+        // Insert the next stem node into the thread list
+        new_stem = _parent[stem];
+        _thread[u] = new_stem;
+        _dirty_revs.push_back(u);
+
+        // Remove the subtree of stem from the thread list
+        w = _rev_thread[stem];
+        _thread[w] = right;
+        _rev_thread[right] = w;
+
+        // Change the parent node and shift stem nodes
+        _parent[stem] = par_stem;
+        par_stem = stem;
+        stem = new_stem;
+
+        // Update u and right
+        u = _last_succ[stem] == _last_succ[par_stem] ?
+          _rev_thread[par_stem] : _last_succ[stem];
+        right = _thread[u];
+      }
+      _parent[u_out] = par_stem;
+      _thread[u] = last;
+      _rev_thread[last] = u;
+      _last_succ[u_out] = u;
+
+      // Remove the subtree of u_out from the thread list except for
+      // the case when old_rev_thread equals to v_in
+      // (it also means that join and v_out coincide)
+      if (old_rev_thread != v_in) {
+        _thread[old_rev_thread] = right;
+        _rev_thread[right] = old_rev_thread;
+      }
+
+      // Update _rev_thread using the new _thread values
+      for (int i = 0; i != int(_dirty_revs.size()); ++i) {
+        u = _dirty_revs[i];
+        _rev_thread[_thread[u]] = u;
+      }
+
+      // Update _pred, _forward, _last_succ and _succ_num for the
+      // stem nodes from u_out to u_in
+      int tmp_sc = 0, tmp_ls = _last_succ[u_out];
+      u = u_out;
+      while (u != u_in) {
+        w = _parent[u];
+        _pred[u] = _pred[w];
+        _forward[u] = !_forward[w];
+        tmp_sc += _succ_num[u] - _succ_num[w];
+        _succ_num[u] = tmp_sc;
+        _last_succ[w] = tmp_ls;
+        u = w;
+      }
+      _pred[u_in] = in_arc;
+      _forward[u_in] = (u_in == _source[in_arc]);
+      _succ_num[u_in] = old_succ_num;
+
+      // Set limits for updating _last_succ form v_in and v_out
+      // towards the root
+      int up_limit_in = -1;
+      int up_limit_out = -1;
+      if (_last_succ[join] == v_in) {
+        up_limit_out = join;
+      } else {
+        up_limit_in = join;
       }
 
       // Update _last_succ from v_in towards the root
-      int up_limit_out = _last_succ[join] == v_in ? join : -1;
-      int last_succ_out = _last_succ[u_out];
-      for (int u = v_in; u != -1 && _last_succ[u] == v_in; u = _parent[u]) {
-        _last_succ[u] = last_succ_out;
-      }
-
+      for (u = v_in; u != up_limit_in && _last_succ[u] == v_in;
+           u = _parent[u]) {
+        _last_succ[u] = _last_succ[u_out];
+      }
       // Update _last_succ from v_out towards the root
       if (join != old_rev_thread && v_in != old_rev_thread) {
-        for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
+        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
              u = _parent[u]) {
           _last_succ[u] = old_rev_thread;
         }
-      }
-      else if (last_succ_out != old_last_succ) {
-        for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
+      } else {
+        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
              u = _parent[u]) {
-          _last_succ[u] = last_succ_out;
+          _last_succ[u] = _last_succ[u_out];
         }
       }
 
       // Update _succ_num from v_in to join
-      for (int u = v_in; u != join; u = _parent[u]) {
+      for (u = v_in; u != join; u = _parent[u]) {
         _succ_num[u] += old_succ_num;
       }
       // Update _succ_num from v_out to join
-      for (int u = v_out; u != join; u = _parent[u]) {
+      for (u = v_out; u != join; u = _parent[u]) {
         _succ_num[u] -= old_succ_num;
       }
     }
 
-    // Update potentials in the subtree that has been moved
+    // Update potentials
     void updatePotential() {
-      Cost sigma = _pi[v_in] - _pi[u_in] -
-                   _pred_dir[u_in] * _cost[in_arc];
+      Cost sigma = _forward[u_in] ?
+        _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
+        _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
+      // Update potentials in the subtree, which has been moved
       int end = _thread[_last_succ[u_in]];
       for (int u = u_in; u != end; u = _thread[u]) {
