RhodeCode

    // Result maps

    FlowMap *_flow_map;

    PotentialMap *_potential_map;

    bool _local_flow;

    bool _local_potential;

    // The number of nodes and arcs in the original graph

    int _node_num;

    int _arc_num;

    // Data structures for storing the graph

    IntNodeMap _node_id;

    ArcVector _arc_ref;

    IntVector _source;

    IntVector _target;

    // Node and arc maps

    CapacityVector _cap;

    CostVector _cost;

    CostVector _supply;

    CapacityVector _flow;

    CostVector _pi;

    // Data for storing the spanning tree structure

    IntVector _parent;

    IntVector _pred;

    IntVector _thread;

    BoolVector _forward;

    IntVector _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;

    Capacity delta;

  private:

    /// \brief Implementation of the "First Eligible" pivot rule for the

    /// \ref NetworkSimplex "network simplex" algorithm.

    ///

    /// 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

        _flow_map = new FlowMap(_graph);

        _local_flow = true;

      }

      if (!_potential_map) {

        _potential_map = new PotentialMap(_graph);

        _local_potential = true;

      }

      // Initialize vectors

      _node_num = countNodes(_graph);

      _arc_num = countArcs(_graph);

      int all_node_num = _node_num + 1;

      int all_arc_num = _arc_num + _node_num;

      _arc_ref.resize(_arc_num);

      _source.resize(all_arc_num);

      _target.resize(all_arc_num);

      _cap.resize(all_arc_num);

      _cost.resize(all_arc_num);

      _supply.resize(all_node_num);

      _flow.resize(all_arc_num, 0);

      _pi.resize(all_node_num, 0);

      _parent.resize(all_node_num);

      _pred.resize(all_node_num);

      _forward.resize(all_node_num);

      _thread.resize(all_node_num);

      _state.resize(all_arc_num, STATE_LOWER);

      // Initialize node related data

      bool valid_supply = true;

      if (_orig_supply) {

        Supply sum = 0;

        int i = 0;

        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {

          _node_id[n] = i;

          _supply[i] = (*_orig_supply)[n];

          sum += _supply[i];

        }

        valid_supply = (sum == 0);

      } else {

        int i = 0;

        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {

          _node_id[n] = i;

          _supply[i] = 0;

        }

        _supply[_node_id[_orig_source]] =  _orig_flow_value;

        _supply[_node_id[_orig_target]] = -_orig_flow_value;

      }

      if (!valid_supply) return false;

      // Set data for the artificial root node

      _root = _node_num;

      _parent[_root] = -1;

      _pred[_root] = -1;

      _thread[_root] = 0;

      _supply[_root] = 0;

      _pi[_root] = 0;

      // Store the arcs in a mixed order

      int k = std::max(int(sqrt(_arc_num)), 10);

      int i = 0;

      for (ArcIt e(_graph); e != INVALID; ++e) {

        _arc_ref[i] = e;

        if ((i += k) >= _arc_num) i = (i % k) + 1;

      }

      // Initialize arc maps

      for (int i = 0; i != _arc_num; ++i) {

        Arc e = _arc_ref[i];

        _source[i] = _node_id[_graph.source(e)];

        _target[i] = _node_id[_graph.target(e)];

        _cost[i] = _orig_cost[e];

        _cap[i] = _orig_cap[e];

      }

      // Remove non-zero lower bounds

      if (_orig_lower) {

        for (int i = 0; i != _arc_num; ++i) {

          Capacity c = (*_orig_lower)[_arc_ref[i]];

          if (c != 0) {

            _cap[i] -= c;

            _supply[_source[i]] -= c;

            _supply[_target[i]] += c;

          }

        }

      }

      // Add artificial arcs and initialize the spanning tree data structure

      Cost max_cost = std::numeric_limits<Cost>::max() / 4;

      Capacity max_cap = std::numeric_limits<Capacity>::max();

      for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {

        _thread[u] = u + 1;

        _parent[u] = _root;

        _pred[u] = e;

        if (_supply[u] >= 0) {

          _flow[e] = _supply[u];

          _forward[u] = true;

          _pi[u] = -max_cost;

        } else {

          _flow[e] = -_supply[u];

          _forward[u] = false;

          _pi[u] = max_cost;

        }

        _cost[e] = max_cost;

        _cap[e] = max_cap;

        _state[e] = STATE_TREE;

      }

      return true;

    }

    // Find the join node

    void findJoinNode() {

      int u = _source[in_arc];

      int v = _target[in_arc];

      while (u != v) {

      }

      join = u;

    }

    // Find the leaving arc of the cycle and returns true if the

    // leaving arc is not the same as the entering arc

    bool findLeavingArc() {

      // Initialize first and second nodes according to the direction

      // of the cycle

      if (_state[in_arc] == STATE_LOWER) {

        first  = _source[in_arc];

        second = _target[in_arc];

      } else {

        first  = _target[in_arc];

        second = _source[in_arc];

      }

      delta = _cap[in_arc];

      int result = 0;

      Capacity d;

      int e;

      // 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];

    // Change _flow and _state vectors

    void changeFlow(bool change) {

      // Augment along the cycle

      if (delta > 0) {

        Capacity val = _state[in_arc] * delta;

        _flow[in_arc] += val;

        for (int u = _source[in_arc]; u != join; u = _parent[u]) {

          _flow[_pred[u]] += _forward[u] ? -val : val;

        }

        for (int u = _target[in_arc]; u != join; u = _parent[u]) {

          _flow[_pred[u]] += _forward[u] ? val : -val;

        }

      }

      // Update the state of the entering and leaving arcs

      if (change) {

        _state[in_arc] = STATE_TREE;

        _state[_pred[u_out]] =

          (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;

      } else {

        _state[in_arc] = -_state[in_arc];

      }

    }

      v_out = _parent[u_out];

      }

      _thread[v_in] = stem = u_in;

      par_stem = v_in;

      while (stem != u_out) {

        _parent[stem] = par_stem;

        par_stem = stem;

        stem = new_stem;

        right = _thread[u];

      }

      _parent[u_out] = par_stem;

      _thread[u] = last;

      } else {

      }

    }

      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]];

      }

    }

    // Execute the algorithm

    bool start(PivotRuleEnum pivot_rule) {

      // Select the pivot rule implementation

      switch (pivot_rule) {

        case FIRST_ELIGIBLE_PIVOT:

          return start<FirstEligiblePivotRule>();

        case BEST_ELIGIBLE_PIVOT:

          return start<BestEligiblePivotRule>();

        case BLOCK_SEARCH_PIVOT:

          return start<BlockSearchPivotRule>();

        case CANDIDATE_LIST_PIVOT:

          return start<CandidateListPivotRule>();

        case ALTERING_LIST_PIVOT:

          return start<AlteringListPivotRule>();

      }

      return false;

    }

    template<class PivotRuleImplementation>

    bool start() {

      PivotRuleImplementation pivot(*this);

      // Execute the network simplex algorithm

      while (pivot.findEnteringArc()) {

        findJoinNode();

        bool change = findLeavingArc();

        changeFlow(change);

        if (change) {

        }

      }

      // Check if the flow amount equals zero on all the artificial arcs

      for (int e = _arc_num; e != _arc_num + _node_num; ++e) {

        if (_flow[e] > 0) return false;

      }

      // Copy flow values to _flow_map

      if (_orig_lower) {

        for (int i = 0; i != _arc_num; ++i) {

          Arc e = _arc_ref[i];

          _flow_map->set(e, (*_orig_lower)[e] + _flow[i]);

        }

      } else {

        for (int i = 0; i != _arc_num; ++i) {

          _flow_map->set(_arc_ref[i], _flow[i]);

        }

      }

      // Copy potential values to _potential_map

      for (NodeIt n(_graph); n != INVALID; ++n) {

        _potential_map->set(n, _pi[_node_id[n]]);

      }