RhodeCode

      // Check the sum of supply values

      _sum_supply = 0;

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

        _sum_supply += _supply[i];

      }

      if ( !((_stype == GEQ && _sum_supply <= 0) ||

             (_stype == LEQ && _sum_supply >= 0)) ) return false;

      // Remove non-zero lower bounds

      if (_have_lower) {

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

          Value c = _lower[i];

          if (c >= 0) {

            _cap[i] = _upper[i] < INF ? _upper[i] - c : INF;

          } else {

            _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF;

          }

          _supply[_source[i]] -= c;

          _supply[_target[i]] += c;

        }

      } else {

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

          _cap[i] = _upper[i];

        }

      }

      // Initialize artifical cost

      Cost ART_COST;

      if (std::numeric_limits<Cost>::is_exact) {

        ART_COST = std::numeric_limits<Cost>::max() / 2 + 1;

      } else {

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

          if (_cost[i] > ART_COST) ART_COST = _cost[i];

        }

        ART_COST = (ART_COST + 1) * _node_num;

      }

      // Initialize arc maps

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

        _flow[i] = 0;

        _state[i] = STATE_LOWER;

      }

      // Set data for the artificial root node

      _root = _node_num;

      _parent[_root] = -1;

      _pred[_root] = -1;

      _thread[_root] = 0;

      _rev_thread[0] = _root;

      _succ_num[_root] = _node_num + 1;

      _last_succ[_root] = _root - 1;

      _supply[_root] = -_sum_supply;

      _pi[_root] = 0;

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

      if (_sum_supply == 0) {

        // EQ supply constraints

        _search_arc_num = _arc_num;

        _all_arc_num = _arc_num + _node_num;

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

          _parent[u] = _root;

          _pred[u] = e;

          _thread[u] = u + 1;

      while (pivot.findEnteringArc()) {

        findJoinNode();

        bool change = findLeavingArc();

        if (delta >= INF) return UNBOUNDED;

        changeFlow(change);

        if (change) {

          updateTreeStructure();

          updatePotential();

        }

      }

      // Check feasibility

      for (int e = _search_arc_num; e != _all_arc_num; ++e) {

        if (_flow[e] != 0) return INFEASIBLE;

      }

      // Transform the solution and the supply map to the original form

      if (_have_lower) {

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

          Value c = _lower[i];

          if (c != 0) {

            _flow[i] += c;

            _supply[_source[i]] += c;

            _supply[_target[i]] -= c;

          }

        }

      }

      // Shift potentials to meet the requirements of the GEQ/LEQ type

      // optimality conditions

      if (_sum_supply == 0) {

        if (_stype == GEQ) {

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

            if (_pi[i] > max_pot) max_pot = _pi[i];

          }

          if (max_pot > 0) {

            for (int i = 0; i != _node_num; ++i)

              _pi[i] -= max_pot;

          }

        } else {

          Cost min_pot = std::numeric_limits<Cost>::max();

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

            if (_pi[i] < min_pot) min_pot = _pi[i];

          }

          if (min_pot < 0) {

            for (int i = 0; i != _node_num; ++i)

              _pi[i] -= min_pot;

          }

        }

      }

      return OPTIMAL;

    }

  }; //class NetworkSimplex

  ///@}

} //namespace lemon

#endif //LEMON_NETWORK_SIMPLEX_H