Index: lemon/network_simplex.h
===================================================================
--- lemon/network_simplex.h	(revision 665)
+++ lemon/network_simplex.h	(revision 687)
@@ -31,7 +31,4 @@
 #include <lemon/core.h>
 #include <lemon/math.h>
-#include <lemon/maps.h>
-#include <lemon/circulation.h>
-#include <lemon/adaptors.h>
 
 namespace lemon {
@@ -51,6 +48,11 @@
   /// In general this class is the fastest implementation available
   /// in LEMON for the minimum cost flow problem.
-  /// Moreover it supports both direction of the supply/demand inequality
-  /// constraints. For more information see \ref ProblemType.
+  /// Moreover it supports both directions of the supply/demand inequality
+  /// 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
+  /// executed using the \ref run() function. If some parameters are not
+  /// specified, then default values will be used.
   ///
   /// \tparam GR The digraph type the algorithm runs on.
@@ -89,65 +91,39 @@
   public:
 
-    /// \brief Enum type for selecting the pivot rule.
-    ///
-    /// Enum type for selecting the pivot rule for the \ref run()
-    /// function.
-    ///
-    /// \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
-    /// 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()
-    /// function with the proper parameter.
-    enum PivotRule {
-
-      /// The First Eligible pivot rule.
-      /// The next eligible arc is selected in a wraparound fashion
-      /// in every iteration.
-      FIRST_ELIGIBLE,
-
-      /// The Best Eligible pivot rule.
-      /// The best eligible arc is selected in every iteration.
-      BEST_ELIGIBLE,
-
-      /// The Block 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.
-      /// 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.
-      /// 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.
-      ALTERING_LIST
+    /// \brief Problem type constants for the \c run() function.
+    ///
+    /// Enum type containing the problem type constants that can be
+    /// returned by the \ref run() function of the algorithm.
+    enum ProblemType {
+      /// The problem has no feasible solution (flow).
+      INFEASIBLE,
+      /// The problem has optimal solution (i.e. it is feasible and
+      /// bounded), and the algorithm has found optimal flow and node
+      /// potentials (primal and dual solutions).
+      OPTIMAL,
+      /// The objective function of the problem is unbounded, i.e.
+      /// there is a directed cycle having negative total cost and
+      /// infinite upper bound.
+      UNBOUNDED
     };
     
-    /// \brief Enum type for selecting the problem type.
-    ///
-    /// Enum type for selecting the problem type, i.e. the direction of
-    /// the inequalities in the supply/demand constraints of the
-    /// \ref min_cost_flow "minimum cost flow problem".
-    ///
-    /// The default problem type is \c GEQ, since this form is supported
+    /// \brief Constants for selecting the type of the supply constraints.
+    ///
+    /// Enum type containing constants for selecting the supply type,
+    /// i.e. the direction of the inequalities in the supply/demand
+    /// constraints of the \ref min_cost_flow "minimum cost flow problem".
+    ///
+    /// The default supply type is \c GEQ, since this form is supported
     /// by other minimum cost flow algorithms and the \ref Circulation
-    /// algorithm as well.
-    /// The \c LEQ problem type can be selected using the \ref problemType()
+    /// algorithm, as well.
+    /// The \c LEQ problem type can be selected using the \ref supplyType()
     /// function.
     ///
-    /// Note that the equality form is a special case of both problem type.
-    enum ProblemType {
-
-      /// This option means that there are "<em>greater or equal</em>"
-      /// constraints in the defintion, i.e. the exact formulation of the
-      /// problem is the following.
+    /// Note that the equality form is a special case of both supply types.
+    enum SupplyType {
+
+      /// This option means that there are <em>"greater or equal"</em>
+      /// supply/demand constraints in the definition, i.e. the exact
+      /// formulation of the problem is the following.
       /**
           \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
@@ -165,7 +141,7 @@
       CARRY_SUPPLIES = GEQ,
 
-      /// This option means that there are "<em>less or equal</em>"
-      /// constraints in the defintion, i.e. the exact formulation of the
-      /// problem is the following.
+      /// This option means that there are <em>"less or equal"</em>
+      /// supply/demand constraints in the definition, i.e. the exact
+      /// formulation of the problem is the following.
       /**
           \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
@@ -183,5 +159,48 @@
       SATISFY_DEMANDS = LEQ
     };
-
+    
+    /// \brief Constants for selecting the pivot rule.
+    ///
+    /// Enum type containing constants for selecting the pivot rule for
+    /// the \ref run() function.
+    ///
+    /// \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
+    /// 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()
+    /// function with the proper parameter.
+    enum PivotRule {
+
+      /// The First Eligible pivot rule.
+      /// The next eligible arc is selected in a wraparound fashion
+      /// in every iteration.
+      FIRST_ELIGIBLE,
+
+      /// The Best Eligible pivot rule.
+      /// The best eligible arc is selected in every iteration.
+      BEST_ELIGIBLE,
+
+      /// The Block 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.
+      /// 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.
+      /// 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.
+      ALTERING_LIST
+    };
+    
   private:
 
@@ -221,5 +240,7 @@
     Node _psource, _ptarget;
     Flow _pstflow;
-    ProblemType _ptype;
+    SupplyType _stype;
+    
+    Flow _sum_supply;
 
     // Result maps
@@ -259,4 +280,13 @@
     int stem, par_stem, new_stem;
     Flow delta;
+
+  public:
+  
+    /// \brief Constant for infinite upper bounds (capacities).
+    ///
+    /// Constant for infinite upper bounds (capacities).
+    /// It is \c std::numeric_limits<Flow>::infinity() if available,
+    /// \c std::numeric_limits<Flow>::max() otherwise.
+    const Flow INF;
 
   private:
@@ -662,15 +692,17 @@
       _graph(graph),
       _plower(NULL), _pupper(NULL), _pcost(NULL),
-      _psupply(NULL), _pstsup(false), _ptype(GEQ),
+      _psupply(NULL), _pstsup(false), _stype(GEQ),
       _flow_map(NULL), _potential_map(NULL),
       _local_flow(false), _local_potential(false),
-      _node_id(graph)
+      _node_id(graph),
+      INF(std::numeric_limits<Flow>::has_infinity ?
+          std::numeric_limits<Flow>::infinity() :
+          std::numeric_limits<Flow>::max())
     {
-      LEMON_ASSERT(std::numeric_limits<Flow>::is_integer &&
-                   std::numeric_limits<Flow>::is_signed,
-        "The flow type of NetworkSimplex must be signed integer");
-      LEMON_ASSERT(std::numeric_limits<Cost>::is_integer &&
-                   std::numeric_limits<Cost>::is_signed,
-        "The cost type of NetworkSimplex must be signed integer");
+      // Check the value types
+      LEMON_ASSERT(std::numeric_limits<Flow>::is_signed,
+        "The flow type of NetworkSimplex must be signed");
+      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
+        "The cost type of NetworkSimplex must be signed");
     }
 
@@ -690,7 +722,6 @@
     ///
     /// This function sets the lower bounds on the arcs.
-    /// If neither this function nor \ref boundMaps() is used before
-    /// calling \ref run(), the lower bounds will be set to zero
-    /// on all arcs.
+    /// If it is not used before calling \ref run(), the lower bounds
+    /// will be set to zero on all arcs.
     ///
     /// \param map An arc map storing the lower bounds.
@@ -699,6 +730,6 @@
     ///
     /// \return <tt>(*this)</tt>
-    template <typename LOWER>
-    NetworkSimplex& lowerMap(const LOWER& map) {
+    template <typename LowerMap>
+    NetworkSimplex& lowerMap(const LowerMap& map) {
       delete _plower;
       _plower = new FlowArcMap(_graph);
@@ -712,8 +743,7 @@
     ///
     /// This function sets the upper bounds (capacities) on the arcs.
-    /// If none of the functions \ref upperMap(), \ref capacityMap()
-    /// and \ref boundMaps() is used before calling \ref run(),
-    /// the upper bounds (capacities) will be set to
-    /// \c std::numeric_limits<Flow>::max() on all arcs.
+    /// If it is not used before calling \ref run(), the upper bounds
+    /// will be set to \ref INF on all arcs (i.e. the flow value will be
+    /// unbounded from above on each arc).
     ///
     /// \param map An arc map storing the upper bounds.
@@ -722,6 +752,6 @@
     ///
     /// \return <tt>(*this)</tt>
-    template<typename UPPER>
-    NetworkSimplex& upperMap(const UPPER& map) {
+    template<typename UpperMap>
+    NetworkSimplex& upperMap(const UpperMap& map) {
       delete _pupper;
       _pupper = new FlowArcMap(_graph);
@@ -732,41 +762,4 @@
     }
 
-    /// \brief Set the upper bounds (capacities) on the arcs.
-    ///
-    /// This function sets the upper bounds (capacities) on the arcs.
-    /// It is just an alias for \ref upperMap().
-    ///
-    /// \return <tt>(*this)</tt>
-    template<typename CAP>
-    NetworkSimplex& capacityMap(const CAP& map) {
-      return upperMap(map);
-    }
-
-    /// \brief Set the lower and upper bounds on the arcs.
-    ///
-    /// This function sets the lower and upper bounds on the arcs.
-    /// If neither this function nor \ref lowerMap() is used before
-    /// calling \ref run(), the lower bounds will be set to zero
-    /// on all arcs.
-    /// If none of the functions \ref upperMap(), \ref capacityMap()
-    /// and \ref boundMaps() is used before calling \ref run(),
-    /// the upper bounds (capacities) will be set to
-    /// \c std::numeric_limits<Flow>::max() on all arcs.
-    ///
-    /// \param lower An arc map storing the lower bounds.
-    /// \param upper An arc map storing the upper bounds.
-    ///
-    /// The \c Value type of the maps must be convertible to the
-    /// \c Flow type of the algorithm.
-    ///
-    /// \note This function is just a shortcut of calling \ref lowerMap()
-    /// and \ref upperMap() separately.
-    ///
-    /// \return <tt>(*this)</tt>
-    template <typename LOWER, typename UPPER>
-    NetworkSimplex& boundMaps(const LOWER& lower, const UPPER& upper) {
-      return lowerMap(lower).upperMap(upper);
-    }
-
     /// \brief Set the costs of the arcs.
     ///
@@ -780,6 +773,6 @@
     ///
     /// \return <tt>(*this)</tt>
-    template<typename COST>
-    NetworkSimplex& costMap(const COST& map) {
+    template<typename CostMap>
+    NetworkSimplex& costMap(const CostMap& map) {
       delete _pcost;
       _pcost = new CostArcMap(_graph);
@@ -802,6 +795,6 @@
     ///
     /// \return <tt>(*this)</tt>
-    template<typename SUP>
-    NetworkSimplex& supplyMap(const SUP& map) {
+    template<typename SupplyMap>
+    NetworkSimplex& supplyMap(const SupplyMap& map) {
       delete _psupply;
       _pstsup = false;
@@ -820,4 +813,8 @@
     /// 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
+    /// assigned to \c t and all other nodes have zero supply value.
     ///
     /// \param s The source node.
@@ -837,15 +834,15 @@
     }
     
-    /// \brief Set the problem type.
-    ///
-    /// This function sets the problem type for the algorithm.
-    /// If it is not used before calling \ref run(), the \ref GEQ problem
+    /// \brief Set the type of the supply constraints.
+    ///
+    /// This function sets the type of the supply/demand constraints.
+    /// If it is not used before calling \ref run(), the \ref GEQ supply
     /// type will be used.
     ///
-    /// For more information see \ref ProblemType.
+    /// For more information see \ref SupplyType.
     ///
     /// \return <tt>(*this)</tt>
-    NetworkSimplex& problemType(ProblemType problem_type) {
-      _ptype = problem_type;
+    NetworkSimplex& supplyType(SupplyType supply_type) {
+      _stype = supply_type;
       return *this;
     }
@@ -897,11 +894,10 @@
     /// This function runs the algorithm.
     /// The paramters can be specified using functions \ref lowerMap(),
-    /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
-    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), 
-    /// \ref problemType(), \ref flowMap() and \ref potentialMap().
+    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), 
+    /// \ref supplyType(), \ref flowMap() and \ref potentialMap().
     /// For example,
     /// \code
     ///   NetworkSimplex<ListDigraph> ns(graph);
-    ///   ns.boundMaps(lower, upper).costMap(cost)
+    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
     ///     .supplyMap(sup).run();
     /// \endcode
@@ -915,7 +911,16 @@
     /// algorithm. For more information see \ref PivotRule.
     ///
-    /// \return \c true if a feasible flow can be found.
-    bool run(PivotRule pivot_rule = BLOCK_SEARCH) {
-      return init() && start(pivot_rule);
+    /// \return \c INFEASIBLE if no feasible flow exists,
+    /// \n \c OPTIMAL if the problem has optimal solution
+    /// (i.e. it is feasible and bounded), and the algorithm has found
+    /// optimal flow and node potentials (primal and dual solutions),
+    /// \n \c UNBOUNDED if the objective function of the problem is
+    /// unbounded, i.e. there is a directed cycle having negative total
+    /// cost and infinite upper bound.
+    ///
+    /// \see ProblemType, PivotRule
+    ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
+      if (!init()) return INFEASIBLE;
+      return start(pivot_rule);
     }
 
@@ -924,6 +929,5 @@
     /// This function resets all the paramaters that have been given
     /// before using functions \ref lowerMap(), \ref upperMap(),
-    /// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
-    /// \ref supplyMap(), \ref stSupply(), \ref problemType(), 
+    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(),
     /// \ref flowMap() and \ref potentialMap().
     ///
@@ -937,5 +941,5 @@
     ///
     ///   // First run
-    ///   ns.lowerMap(lower).capacityMap(cap).costMap(cost)
+    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
     ///     .supplyMap(sup).run();
     ///
@@ -948,5 +952,5 @@
     ///   // (the lower bounds will be set to zero on all arcs)
     ///   ns.reset();
-    ///   ns.capacityMap(cap).costMap(cost)
+    ///   ns.upperMap(capacity).costMap(cost)
     ///     .supplyMap(sup).run();
     /// \endcode
@@ -963,5 +967,5 @@
       _psupply = NULL;
       _pstsup = false;
-      _ptype = GEQ;
+      _stype = GEQ;
       if (_local_flow) delete _flow_map;
       if (_local_potential) delete _potential_map;
@@ -986,5 +990,5 @@
     ///
     /// This function returns the total cost of the found flow.
-    /// The complexity of the function is O(e).
+    /// Its complexity is O(e).
     ///
     /// \note The return type of the function can be specified as a
@@ -998,7 +1002,7 @@
     ///
     /// \pre \ref run() must be called before using this function.
-    template <typename Num>
-    Num totalCost() const {
-      Num c = 0;
+    template <typename Value>
+    Value totalCost() const {
+      Value c = 0;
       if (_pcost) {
         for (ArcIt e(_graph); e != INVALID; ++e)
@@ -1051,5 +1055,5 @@
     /// This function returns a const reference to a node map storing
     /// the found potentials, which form the dual solution of the
-    /// \ref min_cost_flow "minimum cost flow" problem.
+    /// \ref min_cost_flow "minimum cost flow problem".
     ///
     /// \pre \ref run() must be called before using this function.
@@ -1102,5 +1106,5 @@
       // Initialize node related data
       bool valid_supply = true;
-      Flow sum_supply = 0;
+      _sum_supply = 0;
       if (!_pstsup && !_psupply) {
         _pstsup = true;
@@ -1113,8 +1117,8 @@
           _node_id[n] = i;
           _supply[i] = (*_psupply)[n];
-          sum_supply += _supply[i];
-        }
-        valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
-                       (_ptype == LEQ && sum_supply >= 0);
+          _sum_supply += _supply[i];
+        }
+        valid_supply = (_stype == GEQ && _sum_supply <= 0) ||
+                       (_stype == LEQ && _sum_supply >= 0);
       } else {
         int i = 0;
@@ -1128,89 +1132,15 @@
       if (!valid_supply) return false;
 
-      // Infinite capacity value
-      Flow inf_cap =
-        std::numeric_limits<Flow>::has_infinity ?
-        std::numeric_limits<Flow>::infinity() :
-        std::numeric_limits<Flow>::max();
-
       // Initialize artifical cost
-      Cost art_cost;
+      Cost ART_COST;
       if (std::numeric_limits<Cost>::is_exact) {
-        art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
+        ART_COST = std::numeric_limits<Cost>::max() / 4 + 1;
       } else {
-        art_cost = std::numeric_limits<Cost>::min();
+        ART_COST = std::numeric_limits<Cost>::min();
         for (int i = 0; i != _arc_num; ++i) {
-          if (_cost[i] > art_cost) art_cost = _cost[i];
-        }
-        art_cost = (art_cost + 1) * _node_num;
-      }
-
-      // Run Circulation to check if a feasible solution exists
-      typedef ConstMap<Arc, Flow> ConstArcMap;
-      ConstArcMap zero_arc_map(0), inf_arc_map(inf_cap);
-      FlowNodeMap *csup = NULL;
-      bool local_csup = false;
-      if (_psupply) {
-        csup = _psupply;
-      } else {
-        csup = new FlowNodeMap(_graph, 0);
-        (*csup)[_psource] =  _pstflow;
-        (*csup)[_ptarget] = -_pstflow;
-        local_csup = true;
-      }
-      bool circ_result = false;
-      if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {
-        // GEQ problem type
-        if (_plower) {
-          if (_pupper) {
-            Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>
-              circ(_graph, *_plower, *_pupper, *csup);
-            circ_result = circ.run();
-          } else {
-            Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>
-              circ(_graph, *_plower, inf_arc_map, *csup);
-            circ_result = circ.run();
-          }
-        } else {
-          if (_pupper) {
-            Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>
-              circ(_graph, zero_arc_map, *_pupper, *csup);
-            circ_result = circ.run();
-          } else {
-            Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>
-              circ(_graph, zero_arc_map, inf_arc_map, *csup);
-            circ_result = circ.run();
-          }
-        }
-      } else {
-        // LEQ problem type
-        typedef ReverseDigraph<const GR> RevGraph;
-        typedef NegMap<FlowNodeMap> NegNodeMap;
-        RevGraph rgraph(_graph);
-        NegNodeMap neg_csup(*csup);
-        if (_plower) {
-          if (_pupper) {
-            Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>
-              circ(rgraph, *_plower, *_pupper, neg_csup);
-            circ_result = circ.run();
-          } else {
-            Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>
-              circ(rgraph, *_plower, inf_arc_map, neg_csup);
-            circ_result = circ.run();
-          }
-        } else {
-          if (_pupper) {
-            Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>
-              circ(rgraph, zero_arc_map, *_pupper, neg_csup);
-            circ_result = circ.run();
-          } else {
-            Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>
-              circ(rgraph, zero_arc_map, inf_arc_map, neg_csup);
-            circ_result = circ.run();
-          }
-        }
-      }
-      if (local_csup) delete csup;
-      if (!circ_result) return false;
+          if (_cost[i] > ART_COST) ART_COST = _cost[i];
+        }
+        ART_COST = (ART_COST + 1) * _node_num;
+      }
 
       // Set data for the artificial root node
@@ -1222,9 +1152,9 @@
       _succ_num[_root] = all_node_num;
       _last_succ[_root] = _root - 1;
-      _supply[_root] = -sum_supply;
-      if (sum_supply < 0) {
-        _pi[_root] = -art_cost;
+      _supply[_root] = -_sum_supply;
+      if (_sum_supply < 0) {
+        _pi[_root] = -ART_COST;
       } else {
-        _pi[_root] = art_cost;
+        _pi[_root] = ART_COST;
       }
 
@@ -1261,5 +1191,5 @@
         } else {
           for (int i = 0; i != _arc_num; ++i)
-            _cap[i] = inf_cap;
+            _cap[i] = INF;
         }
         if (_pcost) {
@@ -1276,6 +1206,15 @@
         for (int i = 0; i != _arc_num; ++i) {
           Flow c = (*_plower)[_arc_ref[i]];
-          if (c != 0) {
-            _cap[i] -= c;
+          if (c > 0) {
+            if (_cap[i] < INF) _cap[i] -= c;
+            _supply[_source[i]] -= c;
+            _supply[_target[i]] += c;
+          }
+          else if (c < 0) {
+            if (_cap[i] < INF + c) {
+              _cap[i] -= c;
+            } else {
+              _cap[i] = INF;
+            }
             _supply[_source[i]] -= c;
             _supply[_target[i]] += c;
@@ -1292,15 +1231,15 @@
         _parent[u] = _root;
         _pred[u] = e;
-        _cost[e] = art_cost;
-        _cap[e] = inf_cap;
+        _cost[e] = ART_COST;
+        _cap[e] = INF;
         _state[e] = STATE_TREE;
-        if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
+        if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {
           _flow[e] = _supply[u];
           _forward[u] = true;
-          _pi[u] = -art_cost + _pi[_root];
+          _pi[u] = -ART_COST + _pi[_root];
         } else {
           _flow[e] = -_supply[u];
           _forward[u] = false;
-          _pi[u] = art_cost + _pi[_root];
+          _pi[u] = ART_COST + _pi[_root];
         }
       }
@@ -1343,5 +1282,6 @@
       for (int u = first; u != join; u = _parent[u]) {
         e = _pred[u];
-        d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
+        d = _forward[u] ?
+          _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]);
         if (d < delta) {
           delta = d;
@@ -1353,5 +1293,6 @@
       for (int u = second; u != join; u = _parent[u]) {
         e = _pred[u];
-        d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
+        d = _forward[u] ? 
+          (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e];
         if (d <= delta) {
           delta = d;
@@ -1527,5 +1468,5 @@
 
     // Execute the algorithm
-    bool start(PivotRule pivot_rule) {
+    ProblemType start(PivotRule pivot_rule) {
       // Select the pivot rule implementation
       switch (pivot_rule) {
@@ -1541,9 +1482,9 @@
           return start<AlteringListPivotRule>();
       }
-      return false;
+      return INFEASIBLE; // avoid warning
     }
 
     template <typename PivotRuleImpl>
-    bool start() {
+    ProblemType start() {
       PivotRuleImpl pivot(*this);
 
@@ -1552,8 +1493,26 @@
         findJoinNode();
         bool change = findLeavingArc();
+        if (delta >= INF) return UNBOUNDED;
         changeFlow(change);
         if (change) {
           updateTreeStructure();
           updatePotential();
+        }
+      }
+      
+      // Check feasibility
+      if (_sum_supply < 0) {
+        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
+          if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE;
+        }
+      }
+      else if (_sum_supply > 0) {
+        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
+          if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE;
+        }
+      }
+      else {
+        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
+          if (_flow[e] != 0) return INFEASIBLE;
         }
       }
@@ -1575,5 +1534,5 @@
       }
 
-      return true;
+      return OPTIMAL;
     }
 
