Index: doc/groups.dox
===================================================================
--- doc/groups.dox	(revision 478)
+++ doc/groups.dox	(revision 656)
@@ -318,13 +318,13 @@
 The \e maximum \e flow \e problem is to find a flow of maximum value between
 a single source and a single target. Formally, there is a \f$G=(V,A)\f$
-digraph, a \f$cap:A\rightarrow\mathbf{R}^+_0\f$ capacity function and
+digraph, a \f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function and
 \f$s, t \in V\f$ source and target nodes.
-A maximum flow is an \f$f:A\rightarrow\mathbf{R}^+_0\f$ solution of the
+A maximum flow is an \f$f: A\rightarrow\mathbf{R}^+_0\f$ solution of the
 following optimization problem.
 
-\f[ \max\sum_{a\in\delta_{out}(s)}f(a) - \sum_{a\in\delta_{in}(s)}f(a) \f]
-\f[ \sum_{a\in\delta_{out}(v)} f(a) = \sum_{a\in\delta_{in}(v)} f(a)
-    \qquad \forall v\in V\setminus\{s,t\} \f]
-\f[ 0 \leq f(a) \leq cap(a) \qquad \forall a\in A \f]
+\f[ \max\sum_{sv\in A} f(sv) - \sum_{vs\in A} f(vs) \f]
+\f[ \sum_{uv\in A} f(uv) = \sum_{vu\in A} f(vu)
+    \quad \forall u\in V\setminus\{s,t\} \f]
+\f[ 0 \leq f(uv) \leq cap(uv) \quad \forall uv\in A \f]
 
 LEMON contains several algorithms for solving maximum flow problems:
@@ -346,33 +346,101 @@
 \brief Algorithms for finding minimum cost flows and circulations.
 
-This group describes the algorithms for finding minimum cost flows and
+This group contains the algorithms for finding minimum cost flows and
 circulations.
 
 The \e minimum \e cost \e flow \e problem is to find a feasible flow of
 minimum total cost from a set of supply nodes to a set of demand nodes
-in a network with capacity constraints and arc costs.
+in a network with capacity constraints (lower and upper bounds)
+and arc costs.
 Formally, let \f$G=(V,A)\f$ be a digraph,
 \f$lower, upper: A\rightarrow\mathbf{Z}^+_0\f$ denote the lower and
-upper bounds for the flow values on the arcs,
+upper bounds for the flow values on the arcs, for which
+\f$0 \leq lower(uv) \leq upper(uv)\f$ holds for all \f$uv\in A\f$.
 \f$cost: A\rightarrow\mathbf{Z}^+_0\f$ denotes the cost per unit flow
-on the arcs, and
-\f$supply: V\rightarrow\mathbf{Z}\f$ denotes the supply/demand values
-of the nodes.
-A minimum cost flow is an \f$f:A\rightarrow\mathbf{R}^+_0\f$ solution of
-the following optimization problem.
-
-\f[ \min\sum_{a\in A} f(a) cost(a) \f]
-\f[ \sum_{a\in\delta_{out}(v)} f(a) - \sum_{a\in\delta_{in}(v)} f(a) =
-    supply(v) \qquad \forall v\in V \f]
-\f[ lower(a) \leq f(a) \leq upper(a) \qquad \forall a\in A \f]
-
-LEMON contains several algorithms for solving minimum cost flow problems:
- - \ref CycleCanceling Cycle-canceling algorithms.
- - \ref CapacityScaling Successive shortest path algorithm with optional
+on the arcs, and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
+signed supply values of the nodes.
+If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
+supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
+\f$-sup(u)\f$ demand.
+A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}^+_0\f$ solution
+of the following optimization problem.
+
+\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
+\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
+    sup(u) \quad \forall u\in V \f]
+\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
+
+The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
+zero or negative in order to have a feasible solution (since the sum
+of the expressions on the left-hand side of the inequalities is zero).
+It means that the total demand must be greater or equal to the total
+supply and all the supplies have to be carried out from the supply nodes,
+but there could be demands that are not satisfied.
+If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand
+constraints have to be satisfied with equality, i.e. all demands
+have to be satisfied and all supplies have to be used.
+
+If you need the opposite inequalities in the supply/demand constraints
+(i.e. the total demand is less than the total supply and all the demands
+have to be satisfied while there could be supplies that are not used),
+then you could easily transform the problem to the above form by reversing
+the direction of the arcs and taking the negative of the supply values
+(e.g. using \ref ReverseDigraph and \ref NegMap adaptors).
+However \ref NetworkSimplex algorithm also supports this form directly
+for the sake of convenience.
+
+A feasible solution for this problem can be found using \ref Circulation.
+
+Note that the above formulation is actually more general than the usual
+definition of the minimum cost flow problem, in which strict equalities
+are required in the supply/demand contraints, i.e.
+
+\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) =
+    sup(u) \quad \forall u\in V. \f]
+
+However if the sum of the supply values is zero, then these two problems
+are equivalent. So if you need the equality form, you have to ensure this
+additional contraint for the algorithms.
+
+The dual solution of the minimum cost flow problem is represented by node 
+potentials \f$\pi: V\rightarrow\mathbf{Z}\f$.
+An \f$f: A\rightarrow\mathbf{Z}^+_0\f$ feasible solution of the problem
+is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$
+node potentials the following \e complementary \e slackness optimality
+conditions hold.
+
+ - For all \f$uv\in A\f$ arcs:
+   - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
+   - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
+   - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
+ - For all \f$u\in V\f$:
+   - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
+     then \f$\pi(u)=0\f$.
+ 
+Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc
+\f$uv\in A\f$ with respect to the node potentials \f$\pi\f$, i.e.
+\f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f]
+
+All algorithms provide dual solution (node potentials) as well
+if an optimal flow is found.
+
+LEMON contains several algorithms for solving minimum cost flow problems.
+ - \ref NetworkSimplex Primal Network Simplex algorithm with various
+   pivot strategies.
+ - \ref CostScaling Push-Relabel and Augment-Relabel algorithms based on
+   cost scaling.
+ - \ref CapacityScaling Successive Shortest %Path algorithm with optional
    capacity scaling.
- - \ref CostScaling Push-relabel and augment-relabel algorithms based on
-   cost scaling.
- - \ref NetworkSimplex Primal network simplex algorithm with various
-   pivot strategies.
+ - \ref CancelAndTighten The Cancel and Tighten algorithm.
+ - \ref CycleCanceling Cycle-Canceling algorithms.
+
+Most of these implementations support the general inequality form of the
+minimum cost flow problem, but CancelAndTighten and CycleCanceling
+only support the equality form due to the primal method they use.
+
+In general NetworkSimplex is the most efficient implementation,
+but in special cases other algorithms could be faster.
+For example, if the total supply and/or capacities are rather small,
+CapacityScaling is usually the fastest algorithm (without effective scaling).
 */
 
Index: lemon/network_simplex.h
===================================================================
--- lemon/network_simplex.h	(revision 655)
+++ lemon/network_simplex.h	(revision 656)
@@ -31,4 +31,7 @@
 #include <lemon/core.h>
 #include <lemon/math.h>
+#include <lemon/maps.h>
+#include <lemon/circulation.h>
+#include <lemon/adaptors.h>
 
 namespace lemon {
@@ -48,4 +51,6 @@
   /// 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.
   ///
   /// \tparam GR The digraph type the algorithm runs on.
@@ -59,5 +64,6 @@
   ///
   /// \note %NetworkSimplex provides five different pivot rule
-  /// implementations. For more information see \ref PivotRule.
+  /// implementations, from which the most efficient one is used
+  /// by default. For more information see \ref PivotRule.
   template <typename GR, typename F = int, typename C = F>
   class NetworkSimplex
@@ -69,8 +75,15 @@
     /// The cost type of the algorithm
     typedef C Cost;
+#ifdef DOXYGEN
+    /// The type of the flow map
+    typedef GR::ArcMap<Flow> FlowMap;
+    /// The type of the potential map
+    typedef GR::NodeMap<Cost> PotentialMap;
+#else
     /// The type of the flow map
     typedef typename GR::template ArcMap<Flow> FlowMap;
     /// The type of the potential map
     typedef typename GR::template NodeMap<Cost> PotentialMap;
+#endif
 
   public:
@@ -118,4 +131,56 @@
       ALTERING_LIST
     };
+    
+    /// \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
+    /// 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()
+    /// 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.
+      /**
+          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
+          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
+              sup(u) \quad \forall u\in V \f]
+          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
+      */
+      /// It means that the total demand must be greater or equal to the 
+      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
+      /// negative) and all the supplies have to be carried out from 
+      /// the supply nodes, but there could be demands that are not 
+      /// satisfied.
+      GEQ,
+      /// It is just an alias for the \c GEQ option.
+      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.
+      /**
+          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
+          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
+              sup(u) \quad \forall u\in V \f]
+          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
+      */
+      /// It means that the total demand must be less or equal to the 
+      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
+      /// positive) and all the demands have to be satisfied, but there
+      /// could be supplies that are not carried out from the supply
+      /// nodes.
+      LEQ,
+      /// It is just an alias for the \c LEQ option.
+      SATISFY_DEMANDS = LEQ
+    };
 
   private:
@@ -156,4 +221,5 @@
     Node _psource, _ptarget;
     Flow _pstflow;
+    ProblemType _ptype;
 
     // Result maps
@@ -587,5 +653,5 @@
     /// \brief Constructor.
     ///
-    /// Constructor.
+    /// The constructor of the class.
     ///
     /// \param graph The digraph the algorithm runs on.
@@ -593,5 +659,5 @@
       _graph(graph),
       _plower(NULL), _pupper(NULL), _pcost(NULL),
-      _psupply(NULL), _pstsup(false),
+      _psupply(NULL), _pstsup(false), _ptype(GEQ),
       _flow_map(NULL), _potential_map(NULL),
       _local_flow(false), _local_potential(false),
@@ -612,4 +678,10 @@
     }
 
+    /// \name Parameters
+    /// The parameters of the algorithm can be specified using these
+    /// functions.
+
+    /// @{
+
     /// \brief Set the lower bounds on the arcs.
     ///
@@ -761,4 +833,18 @@
       return *this;
     }
+    
+    /// \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
+    /// type will be used.
+    ///
+    /// For more information see \ref ProblemType.
+    ///
+    /// \return <tt>(*this)</tt>
+    NetworkSimplex& problemType(ProblemType problem_type) {
+      _ptype = problem_type;
+      return *this;
+    }
 
     /// \brief Set the flow map.
@@ -796,4 +882,6 @@
       return *this;
     }
+    
+    /// @}
 
     /// \name Execution Control
@@ -805,8 +893,9 @@
     ///
     /// This function runs the algorithm.
-    /// The paramters can be specified using \ref lowerMap(),
+    /// The paramters can be specified using functions \ref lowerMap(),
     /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
-    /// \ref costMap(), \ref supplyMap() and \ref stSupply()
-    /// functions. For example,
+    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), 
+    /// \ref problemType(), \ref flowMap() and \ref potentialMap().
+    /// For example,
     /// \code
     ///   NetworkSimplex<ListDigraph> ns(graph);
@@ -831,7 +920,8 @@
     ///
     /// This function resets all the paramaters that have been given
-    /// using \ref lowerMap(), \ref upperMap(), \ref capacityMap(),
-    /// \ref boundMaps(), \ref costMap(), \ref supplyMap() and
-    /// \ref stSupply() functions before.
+    /// before using functions \ref lowerMap(), \ref upperMap(),
+    /// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
+    /// \ref supplyMap(), \ref stSupply(), \ref problemType(), 
+    /// \ref flowMap() and \ref potentialMap().
     ///
     /// It is useful for multiple run() calls. If this function is not
@@ -870,4 +960,12 @@
       _psupply = NULL;
       _pstsup = false;
+      _ptype = GEQ;
+      if (_local_flow) delete _flow_map;
+      if (_local_potential) delete _potential_map;
+      _flow_map = NULL;
+      _potential_map = NULL;
+      _local_flow = false;
+      _local_potential = false;
+
       return *this;
     }
@@ -1001,4 +1099,5 @@
       // Initialize node related data
       bool valid_supply = true;
+      Flow sum_supply = 0;
       if (!_pstsup && !_psupply) {
         _pstsup = true;
@@ -1007,12 +1106,12 @@
       }
       if (_psupply) {
-        Flow sum = 0;
         int i = 0;
         for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
           _node_id[n] = i;
           _supply[i] = (*_psupply)[n];
-          sum += _supply[i];
-        }
-        valid_supply = (sum == 0);
+          sum_supply += _supply[i];
+        }
+        valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
+                       (_ptype == LEQ && sum_supply >= 0);
       } else {
         int i = 0;
@@ -1022,7 +1121,92 @@
         }
         _supply[_node_id[_psource]] =  _pstflow;
-        _supply[_node_id[_ptarget]]   = -_pstflow;
+        _supply[_node_id[_ptarget]] = -_pstflow;
       }
       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;
+      if (std::numeric_limits<Cost>::is_exact) {
+        art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
+      } else {
+        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;
+      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, ConstArcMap(inf_cap), *csup);
+            circ_result = circ.run();
+          }
+        } else {
+          if (_pupper) {
+            Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>
+              circ(_graph, ConstArcMap(0), *_pupper, *csup);
+            circ_result = circ.run();
+          } else {
+            Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>
+              circ(_graph, ConstArcMap(0), ConstArcMap(inf_cap), *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, ConstArcMap(inf_cap), neg_csup);
+            circ_result = circ.run();
+          }
+        } else {
+          if (_pupper) {
+            Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>
+              circ(rgraph, ConstArcMap(0), *_pupper, neg_csup);
+            circ_result = circ.run();
+          } else {
+            Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>
+              circ(rgraph, ConstArcMap(0), ConstArcMap(inf_cap), neg_csup);
+            circ_result = circ.run();
+          }
+        }
+      }
+      if (local_csup) delete csup;
+      if (!circ_result) return false;
 
       // Set data for the artificial root node
@@ -1034,6 +1218,10 @@
       _succ_num[_root] = all_node_num;
       _last_succ[_root] = _root - 1;
-      _supply[_root] = 0;
-      _pi[_root] = 0;
+      _supply[_root] = -sum_supply;
+      if (sum_supply < 0) {
+        _pi[_root] = -art_cost;
+      } else {
+        _pi[_root] = art_cost;
+      }
 
       // Store the arcs in a mixed order
@@ -1046,8 +1234,4 @@
 
       // Initialize arc maps
-      Flow inf_cap =
-        std::numeric_limits<Flow>::has_infinity ?
-        std::numeric_limits<Flow>::infinity() :
-        std::numeric_limits<Flow>::max();
       if (_pupper && _pcost) {
         for (int i = 0; i != _arc_num; ++i) {
@@ -1084,16 +1268,4 @@
       }
       
-      // Initialize artifical cost
-      Cost art_cost;
-      if (std::numeric_limits<Cost>::is_exact) {
-        art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
-      } else {
-        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;
-      }
-
       // Remove non-zero lower bounds
       if (_plower) {
@@ -1119,12 +1291,12 @@
         _cap[e] = inf_cap;
         _state[e] = STATE_TREE;
-        if (_supply[u] >= 0) {
+        if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
           _flow[e] = _supply[u];
           _forward[u] = true;
-          _pi[u] = -art_cost;
+          _pi[u] = -art_cost + _pi[_root];
         } else {
           _flow[e] = -_supply[u];
           _forward[u] = false;
-          _pi[u] = art_cost;
+          _pi[u] = art_cost + _pi[_root];
         }
       }
@@ -1383,9 +1555,4 @@
       }
 
-      // 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 (_plower) {
Index: test/min_cost_flow_test.cc
===================================================================
--- test/min_cost_flow_test.cc	(revision 654)
+++ test/min_cost_flow_test.cc	(revision 656)
@@ -34,17 +34,17 @@
 char test_lgf[] =
   "@nodes\n"
-  "label  sup1 sup2 sup3\n"
-  "    1    20   27    0\n"
-  "    2    -4    0    0\n"
-  "    3     0    0    0\n"
-  "    4     0    0    0\n"
-  "    5     9    0    0\n"
-  "    6    -6    0    0\n"
-  "    7     0    0    0\n"
-  "    8     0    0    0\n"
-  "    9     3    0    0\n"
-  "   10    -2    0    0\n"
-  "   11     0    0    0\n"
-  "   12   -20  -27    0\n"
+  "label  sup1 sup2 sup3 sup4 sup5\n"
+  "    1    20   27    0   20   30\n"
+  "    2    -4    0    0   -8   -3\n"
+  "    3     0    0    0    0    0\n"
+  "    4     0    0    0    0    0\n"
+  "    5     9    0    0    6   11\n"
+  "    6    -6    0    0   -5   -6\n"
+  "    7     0    0    0    0    0\n"
+  "    8     0    0    0    0    3\n"
+  "    9     3    0    0    0    0\n"
+  "   10    -2    0    0   -7   -2\n"
+  "   11     0    0    0  -10    0\n"
+  "   12   -20  -27    0  -30  -20\n"
   "\n"
   "@arcs\n"
@@ -77,4 +77,10 @@
 
 
+enum ProblemType {
+  EQ,
+  GEQ,
+  LEQ
+};
+
 // Check the interface of an MCF algorithm
 template <typename GR, typename Flow, typename Cost>
@@ -98,15 +104,17 @@
              .supplyMap(sup)
              .stSupply(n, n, k)
+             .flowMap(flow)
+             .potentialMap(pot)
              .run();
-
-      const typename MCF::FlowMap &fm = mcf.flowMap();
-      const typename MCF::PotentialMap &pm = mcf.potentialMap();
-
-      v = mcf.totalCost();
-      double x = mcf.template totalCost<double>();
-      v = mcf.flow(a);
-      v = mcf.potential(n);
-      mcf.flowMap(flow);
-      mcf.potentialMap(pot);
+      
+      const MCF& const_mcf = mcf;
+
+      const typename MCF::FlowMap &fm = const_mcf.flowMap();
+      const typename MCF::PotentialMap &pm = const_mcf.potentialMap();
+
+      v = const_mcf.totalCost();
+      double x = const_mcf.template totalCost<double>();
+      v = const_mcf.flow(a);
+      v = const_mcf.potential(n);
 
       ignore_unused_variable_warning(fm);
@@ -143,5 +151,6 @@
            typename SM, typename FM >
 bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
-                const SM& supply, const FM& flow )
+                const SM& supply, const FM& flow,
+                ProblemType type = EQ )
 {
   TEMPLATE_DIGRAPH_TYPEDEFS(GR);
@@ -157,5 +166,8 @@
     for (InArcIt e(gr, n); e != INVALID; ++e)
       sum -= flow[e];
-    if (sum != supply[n]) return false;
+    bool b = (type ==  EQ && sum == supply[n]) ||
+             (type == GEQ && sum >= supply[n]) ||
+             (type == LEQ && sum <= supply[n]);
+    if (!b) return false;
   }
 
@@ -166,7 +178,8 @@
 // using the "Complementary Slackness" optimality condition
 template < typename GR, typename LM, typename UM,
-           typename CM, typename FM, typename PM >
+           typename CM, typename SM, typename FM, typename PM >
 bool checkPotential( const GR& gr, const LM& lower, const UM& upper,
-                     const CM& cost, const FM& flow, const PM& pi )
+                     const CM& cost, const SM& supply, const FM& flow, 
+                     const PM& pi )
 {
   TEMPLATE_DIGRAPH_TYPEDEFS(GR);
@@ -180,4 +193,14 @@
           (red_cost < 0 && flow[e] == upper[e]);
   }
+  
+  for (NodeIt n(gr); opt && n != INVALID; ++n) {
+    typename SM::Value sum = 0;
+    for (OutArcIt e(gr, n); e != INVALID; ++e)
+      sum += flow[e];
+    for (InArcIt e(gr, n); e != INVALID; ++e)
+      sum -= flow[e];
+    opt = (sum == supply[n]) || (pi[n] == 0);
+  }
+  
   return opt;
 }
@@ -191,12 +214,13 @@
                const CM& cost, const SM& supply,
                bool result, typename CM::Value total,
-               const std::string &test_id = "" )
+               const std::string &test_id = "",
+               ProblemType type = EQ )
 {
   check(mcf_result == result, "Wrong result " + test_id);
   if (result) {
-    check(checkFlow(gr, lower, upper, supply, mcf.flowMap()),
+    check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type),
           "The flow is not feasible " + test_id);
     check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
-    check(checkPotential(gr, lower, upper, cost, mcf.flowMap(),
+    check(checkPotential(gr, lower, upper, cost, supply, mcf.flowMap(),
                          mcf.potentialMap()),
           "Wrong potentials " + test_id);
@@ -227,5 +251,5 @@
   Digraph gr;
   Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), u(gr);
-  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr);
+  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr);
   ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
   Node v, w;
@@ -240,4 +264,6 @@
     .nodeMap("sup2", s2)
     .nodeMap("sup3", s3)
+    .nodeMap("sup4", s4)
+    .nodeMap("sup5", s5)
     .node("source", v)
     .node("target", w)
@@ -248,4 +274,5 @@
     NetworkSimplex<Digraph> mcf(gr);
 
+    // Check the equality form
     mcf.upperMap(u).costMap(c);
     checkMcf(mcf, mcf.supplyMap(s1).run(),
@@ -268,4 +295,26 @@
     checkMcf(mcf, mcf.boundMaps(l2, u).run(),
              gr, l2, u, cc, s3, false,   0, "#A8");
+
+    // Check the GEQ form
+    mcf.reset().upperMap(u).costMap(c).supplyMap(s4);
+    checkMcf(mcf, mcf.run(),
+             gr, l1, u, c, s4, true,  3530, "#A9", GEQ);
+    mcf.problemType(mcf.GEQ);
+    checkMcf(mcf, mcf.lowerMap(l2).run(),
+             gr, l2, u, c, s4, true,  4540, "#A10", GEQ);
+    mcf.problemType(mcf.CARRY_SUPPLIES).supplyMap(s5);
+    checkMcf(mcf, mcf.run(),
+             gr, l2, u, c, s5, false,    0, "#A11", GEQ);
+
+    // Check the LEQ form
+    mcf.reset().problemType(mcf.LEQ);
+    mcf.upperMap(u).costMap(c).supplyMap(s5);
+    checkMcf(mcf, mcf.run(),
+             gr, l1, u, c, s5, true,  5080, "#A12", LEQ);
+    checkMcf(mcf, mcf.lowerMap(l2).run(),
+             gr, l2, u, c, s5, true,  5930, "#A13", LEQ);
+    mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4);
+    checkMcf(mcf, mcf.run(),
+             gr, l2, u, c, s4, false,    0, "#A14", LEQ);
   }
 
