RhodeCode

    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:

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

        }

        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

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

        }

      }

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

    }

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