# HG changeset patch # User Peter Kovacs # Date 1242122800 -7200 # Node ID 8b0df68370a49ce0073f49133d998301fd07a16e # Parent 4d3d1a2cd23dce4b6d484435ebd940aece722131 Fix the GEQ/LEQ handling in NetworkSimplex + improve doc (#291) - Fix the optimality conditions for the GEQ/LEQ form. - Fix the initialization of the algortihm. It ensures correct solutions and it is much faster for the inequality forms. - Fix the pivot rules to search all the arcs that have to be allowed to get in the basis. - Better block size for the Block Search pivot rule. - Improve documentation of the problem and move it to a separate page. diff -r 4d3d1a2cd23d -r 8b0df68370a4 doc/Makefile.am --- a/doc/Makefile.am Mon May 11 16:38:21 2009 +0100 +++ b/doc/Makefile.am Tue May 12 12:06:40 2009 +0200 @@ -8,6 +8,7 @@ doc/license.dox \ doc/mainpage.dox \ doc/migration.dox \ + doc/min_cost_flow.dox \ doc/named-param.dox \ doc/namespaces.dox \ doc/html \ diff -r 4d3d1a2cd23d -r 8b0df68370a4 doc/groups.dox --- a/doc/groups.dox Mon May 11 16:38:21 2009 +0100 +++ b/doc/groups.dox Tue May 12 12:06:40 2009 +0200 @@ -335,91 +335,16 @@ */ /** -@defgroup min_cost_flow Minimum Cost Flow Algorithms +@defgroup min_cost_flow_algs Minimum Cost Flow Algorithms @ingroup algs \brief Algorithms for finding minimum cost flows and circulations. This group contains the algorithms for finding minimum cost flows and -circulations. +circulations. For more information about this problem and its dual +solution see \ref min_cost_flow "Minimum Cost Flow Problem". -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 (lower and upper bounds) -and arc costs. -Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{Z}\f$, -\f$upper: A\rightarrow\mathbf{Z}\cup\{+\infty\}\f$ denote the lower and -upper bounds for the flow values on the arcs, for which -\f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$, -\f$cost: A\rightarrow\mathbf{Z}\f$ denotes the cost per unit flow -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}\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}\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)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{R}\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. + + +\section mcf_algs Algorithms + +LEMON contains several algorithms for solving this problem, for more +information see \ref min_cost_flow_algs "Minimum Cost Flow Algorithms". + +A feasible solution for this problem can be found using \ref Circulation. + + +\section mcf_dual Dual Solution + +The dual solution of the minimum cost flow problem is represented by +node potentials \f$\pi: V\rightarrow\mathbf{R}\f$. +An \f$f: A\rightarrow\mathbf{R}\f$ primal feasible solution is optimal +if and only if for some \f$\pi: V\rightarrow\mathbf{R}\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)"less or equal" (LEQ) supply/demand constraints, +instead of the "greater or equal" (GEQ) constraints. + +\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. +The equality form is also a special case of this form, of course. + +You could easily transform this case to the \ref mcf_def "GEQ form" +of the problem 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. + +Note that the optimality conditions for this supply constraint type are +slightly differ from the conditions that are discussed for the GEQ form, +namely the potentials have to be non-negative instead of non-positive. +An \f$f: A\rightarrow\mathbf{R}\f$ feasible solution of this problem +is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{R}\f$ +node potentials the following 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)=0\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$. + +*/ +} diff -r 4d3d1a2cd23d -r 8b0df68370a4 lemon/network_simplex.h --- a/lemon/network_simplex.h Mon May 11 16:38:21 2009 +0100 +++ b/lemon/network_simplex.h Tue May 12 12:06:40 2009 +0200 @@ -19,7 +19,7 @@ #ifndef LEMON_NETWORK_SIMPLEX_H #define LEMON_NETWORK_SIMPLEX_H -/// \ingroup min_cost_flow +/// \ingroup min_cost_flow_algs /// /// \file /// \brief Network Simplex algorithm for finding a minimum cost flow. @@ -33,7 +33,7 @@ namespace lemon { - /// \addtogroup min_cost_flow + /// \addtogroup min_cost_flow_algs /// @{ /// \brief Implementation of the primal Network Simplex algorithm @@ -102,50 +102,16 @@ /// 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 supplyType() - /// function. - /// - /// Note that the equality form is a special case of both supply types. + /// The default supply type is \c GEQ, the \c LEQ type can be + /// selected using \ref supplyType(). + /// The equality form is a special case of both supply types. enum SupplyType { - /// This option means that there are "greater or equal" - /// 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] - \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. + /// supply/demand constraints in the definition of the problem. GEQ, - /// It is just an alias for the \c GEQ option. - CARRY_SUPPLIES = GEQ, - /// This option means that there are "less or equal" - /// 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] - \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 + /// supply/demand constraints in the definition of the problem. + LEQ }; /// \brief Constants for selecting the pivot rule. @@ -215,6 +181,8 @@ const GR &_graph; int _node_num; int _arc_num; + int _all_arc_num; + int _search_arc_num; // Parameters of the problem bool _have_lower; @@ -277,7 +245,7 @@ const IntVector &_state; const CostVector &_pi; int &_in_arc; - int _arc_num; + int _search_arc_num; // Pivot rule data int _next_arc; @@ -288,13 +256,14 @@ FirstEligiblePivotRule(NetworkSimplex &ns) : _source(ns._source), _target(ns._target), _cost(ns._cost), _state(ns._state), _pi(ns._pi), - _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0) + _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), + _next_arc(0) {} // Find next entering arc bool findEnteringArc() { Cost c; - for (int e = _next_arc; e < _arc_num; ++e) { + for (int e = _next_arc; e < _search_arc_num; ++e) { c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); if (c < 0) { _in_arc = e; @@ -328,7 +297,7 @@ const IntVector &_state; const CostVector &_pi; int &_in_arc; - int _arc_num; + int _search_arc_num; public: @@ -336,13 +305,13 @@ BestEligiblePivotRule(NetworkSimplex &ns) : _source(ns._source), _target(ns._target), _cost(ns._cost), _state(ns._state), _pi(ns._pi), - _in_arc(ns.in_arc), _arc_num(ns._arc_num) + _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num) {} // Find next entering arc bool findEnteringArc() { Cost c, min = 0; - for (int e = 0; e < _arc_num; ++e) { + for (int e = 0; e < _search_arc_num; ++e) { c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); if (c < min) { min = c; @@ -367,7 +336,7 @@ const IntVector &_state; const CostVector &_pi; int &_in_arc; - int _arc_num; + int _search_arc_num; // Pivot rule data int _block_size; @@ -379,14 +348,15 @@ BlockSearchPivotRule(NetworkSimplex &ns) : _source(ns._source), _target(ns._target), _cost(ns._cost), _state(ns._state), _pi(ns._pi), - _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0) + _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), + _next_arc(0) { // The main parameters of the pivot rule - const double BLOCK_SIZE_FACTOR = 2.0; + const double BLOCK_SIZE_FACTOR = 0.5; const int MIN_BLOCK_SIZE = 10; _block_size = std::max( int(BLOCK_SIZE_FACTOR * - std::sqrt(double(_arc_num))), + std::sqrt(double(_search_arc_num))), MIN_BLOCK_SIZE ); } @@ -395,7 +365,7 @@ Cost c, min = 0; int cnt = _block_size; int e, min_arc = _next_arc; - for (e = _next_arc; e < _arc_num; ++e) { + for (e = _next_arc; e < _search_arc_num; ++e) { c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); if (c < min) { min = c; @@ -440,7 +410,7 @@ const IntVector &_state; const CostVector &_pi; int &_in_arc; - int _arc_num; + int _search_arc_num; // Pivot rule data IntVector _candidates; @@ -454,7 +424,8 @@ CandidateListPivotRule(NetworkSimplex &ns) : _source(ns._source), _target(ns._target), _cost(ns._cost), _state(ns._state), _pi(ns._pi), - _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0) + _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), + _next_arc(0) { // The main parameters of the pivot rule const double LIST_LENGTH_FACTOR = 1.0; @@ -463,7 +434,7 @@ const int MIN_MINOR_LIMIT = 3; _list_length = std::max( int(LIST_LENGTH_FACTOR * - std::sqrt(double(_arc_num))), + std::sqrt(double(_search_arc_num))), MIN_LIST_LENGTH ); _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length), MIN_MINOR_LIMIT ); @@ -500,7 +471,7 @@ // Major iteration: build a new candidate list min = 0; _curr_length = 0; - for (e = _next_arc; e < _arc_num; ++e) { + for (e = _next_arc; e < _search_arc_num; ++e) { c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); if (c < 0) { _candidates[_curr_length++] = e; @@ -546,7 +517,7 @@ const IntVector &_state; const CostVector &_pi; int &_in_arc; - int _arc_num; + int _search_arc_num; // Pivot rule data int _block_size, _head_length, _curr_length; @@ -574,8 +545,8 @@ AlteringListPivotRule(NetworkSimplex &ns) : _source(ns._source), _target(ns._target), _cost(ns._cost), _state(ns._state), _pi(ns._pi), - _in_arc(ns.in_arc), _arc_num(ns._arc_num), - _next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost) + _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), + _next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost) { // The main parameters of the pivot rule const double BLOCK_SIZE_FACTOR = 1.5; @@ -584,7 +555,7 @@ const int MIN_HEAD_LENGTH = 3; _block_size = std::max( int(BLOCK_SIZE_FACTOR * - std::sqrt(double(_arc_num))), + std::sqrt(double(_search_arc_num))), MIN_BLOCK_SIZE ); _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size), MIN_HEAD_LENGTH ); @@ -610,7 +581,7 @@ int last_arc = 0; int limit = _head_length; - for (int e = _next_arc; e < _arc_num; ++e) { + for (int e = _next_arc; e < _search_arc_num; ++e) { _cand_cost[e] = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); if (_cand_cost[e] < 0) { @@ -678,17 +649,17 @@ _node_num = countNodes(_graph); _arc_num = countArcs(_graph); int all_node_num = _node_num + 1; - int all_arc_num = _arc_num + _node_num; + int max_arc_num = _arc_num + 2 * _node_num; - _source.resize(all_arc_num); - _target.resize(all_arc_num); + _source.resize(max_arc_num); + _target.resize(max_arc_num); - _lower.resize(all_arc_num); - _upper.resize(all_arc_num); - _cap.resize(all_arc_num); - _cost.resize(all_arc_num); + _lower.resize(_arc_num); + _upper.resize(_arc_num); + _cap.resize(max_arc_num); + _cost.resize(max_arc_num); _supply.resize(all_node_num); - _flow.resize(all_arc_num); + _flow.resize(max_arc_num); _pi.resize(all_node_num); _parent.resize(all_node_num); @@ -698,7 +669,7 @@ _rev_thread.resize(all_node_num); _succ_num.resize(all_node_num); _last_succ.resize(all_node_num); - _state.resize(all_arc_num); + _state.resize(max_arc_num); // Copy the graph (store the arcs in a mixed order) int i = 0; @@ -1069,7 +1040,7 @@ // Initialize artifical cost Cost ART_COST; if (std::numeric_limits::is_exact) { - ART_COST = std::numeric_limits::max() / 4 + 1; + ART_COST = std::numeric_limits::max() / 2 + 1; } else { ART_COST = std::numeric_limits::min(); for (int i = 0; i != _arc_num; ++i) { @@ -1093,29 +1064,121 @@ _succ_num[_root] = _node_num + 1; _last_succ[_root] = _root - 1; _supply[_root] = -_sum_supply; - _pi[_root] = _sum_supply < 0 ? -ART_COST : ART_COST; + _pi[_root] = 0; // Add artificial arcs and initialize the spanning tree data structure - for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { - _parent[u] = _root; - _pred[u] = e; - _thread[u] = u + 1; - _rev_thread[u + 1] = u; - _succ_num[u] = 1; - _last_succ[u] = u; - _cost[e] = ART_COST; - _cap[e] = INF; - _state[e] = STATE_TREE; - if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) { - _flow[e] = _supply[u]; - _forward[u] = true; - _pi[u] = -ART_COST + _pi[_root]; - } else { - _flow[e] = -_supply[u]; - _forward[u] = false; - _pi[u] = ART_COST + _pi[_root]; + 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; + _rev_thread[u + 1] = u; + _succ_num[u] = 1; + _last_succ[u] = u; + _cap[e] = INF; + _state[e] = STATE_TREE; + if (_supply[u] >= 0) { + _forward[u] = true; + _pi[u] = 0; + _source[e] = u; + _target[e] = _root; + _flow[e] = _supply[u]; + _cost[e] = 0; + } else { + _forward[u] = false; + _pi[u] = ART_COST; + _source[e] = _root; + _target[e] = u; + _flow[e] = -_supply[u]; + _cost[e] = ART_COST; + } } } + else if (_sum_supply > 0) { + // LEQ supply constraints + _search_arc_num = _arc_num + _node_num; + int f = _arc_num + _node_num; + for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { + _parent[u] = _root; + _thread[u] = u + 1; + _rev_thread[u + 1] = u; + _succ_num[u] = 1; + _last_succ[u] = u; + if (_supply[u] >= 0) { + _forward[u] = true; + _pi[u] = 0; + _pred[u] = e; + _source[e] = u; + _target[e] = _root; + _cap[e] = INF; + _flow[e] = _supply[u]; + _cost[e] = 0; + _state[e] = STATE_TREE; + } else { + _forward[u] = false; + _pi[u] = ART_COST; + _pred[u] = f; + _source[f] = _root; + _target[f] = u; + _cap[f] = INF; + _flow[f] = -_supply[u]; + _cost[f] = ART_COST; + _state[f] = STATE_TREE; + _source[e] = u; + _target[e] = _root; + _cap[e] = INF; + _flow[e] = 0; + _cost[e] = 0; + _state[e] = STATE_LOWER; + ++f; + } + } + _all_arc_num = f; + } + else { + // GEQ supply constraints + _search_arc_num = _arc_num + _node_num; + int f = _arc_num + _node_num; + for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { + _parent[u] = _root; + _thread[u] = u + 1; + _rev_thread[u + 1] = u; + _succ_num[u] = 1; + _last_succ[u] = u; + if (_supply[u] <= 0) { + _forward[u] = false; + _pi[u] = 0; + _pred[u] = e; + _source[e] = _root; + _target[e] = u; + _cap[e] = INF; + _flow[e] = -_supply[u]; + _cost[e] = 0; + _state[e] = STATE_TREE; + } else { + _forward[u] = true; + _pi[u] = -ART_COST; + _pred[u] = f; + _source[f] = u; + _target[f] = _root; + _cap[f] = INF; + _flow[f] = _supply[u]; + _state[f] = STATE_TREE; + _cost[f] = ART_COST; + _source[e] = _root; + _target[e] = u; + _cap[e] = INF; + _flow[e] = 0; + _cost[e] = 0; + _state[e] = STATE_LOWER; + ++f; + } + } + _all_arc_num = f; + } return true; } @@ -1374,20 +1437,8 @@ } // 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; - } + 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 @@ -1401,6 +1452,30 @@ } } } + + // Shift potentials to meet the requirements of the GEQ/LEQ type + // optimality conditions + if (_sum_supply == 0) { + if (_stype == GEQ) { + Cost max_pot = std::numeric_limits::min(); + 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::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; }