Fix the GEQ/LEQ handling in NetworkSimplex + improve doc (#291)
authorPeter Kovacs <kpeter@inf.elte.hu>
Tue, 12 May 2009 12:06:40 +0200
changeset 6638b0df68370a4
parent 662 4d3d1a2cd23d
child 664 cc61d09f053b
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.
doc/Makefile.am
doc/groups.dox
doc/min_cost_flow.dox
lemon/network_simplex.h
     1.1 --- a/doc/Makefile.am	Mon May 11 16:38:21 2009 +0100
     1.2 +++ b/doc/Makefile.am	Tue May 12 12:06:40 2009 +0200
     1.3 @@ -8,6 +8,7 @@
     1.4  	doc/license.dox \
     1.5  	doc/mainpage.dox \
     1.6  	doc/migration.dox \
     1.7 +	doc/min_cost_flow.dox \
     1.8  	doc/named-param.dox \
     1.9  	doc/namespaces.dox \
    1.10  	doc/html \
     2.1 --- a/doc/groups.dox	Mon May 11 16:38:21 2009 +0100
     2.2 +++ b/doc/groups.dox	Tue May 12 12:06:40 2009 +0200
     2.3 @@ -335,91 +335,16 @@
     2.4  */
     2.5  
     2.6  /**
     2.7 -@defgroup min_cost_flow Minimum Cost Flow Algorithms
     2.8 +@defgroup min_cost_flow_algs Minimum Cost Flow Algorithms
     2.9  @ingroup algs
    2.10  
    2.11  \brief Algorithms for finding minimum cost flows and circulations.
    2.12  
    2.13  This group contains the algorithms for finding minimum cost flows and
    2.14 -circulations.
    2.15 +circulations. For more information about this problem and its dual
    2.16 +solution see \ref min_cost_flow "Minimum Cost Flow Problem".
    2.17  
    2.18 -The \e minimum \e cost \e flow \e problem is to find a feasible flow of
    2.19 -minimum total cost from a set of supply nodes to a set of demand nodes
    2.20 -in a network with capacity constraints (lower and upper bounds)
    2.21 -and arc costs.
    2.22 -Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{Z}\f$,
    2.23 -\f$upper: A\rightarrow\mathbf{Z}\cup\{+\infty\}\f$ denote the lower and
    2.24 -upper bounds for the flow values on the arcs, for which
    2.25 -\f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$,
    2.26 -\f$cost: A\rightarrow\mathbf{Z}\f$ denotes the cost per unit flow
    2.27 -on the arcs and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
    2.28 -signed supply values of the nodes.
    2.29 -If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
    2.30 -supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
    2.31 -\f$-sup(u)\f$ demand.
    2.32 -A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}\f$ solution
    2.33 -of the following optimization problem.
    2.34 -
    2.35 -\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
    2.36 -\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
    2.37 -    sup(u) \quad \forall u\in V \f]
    2.38 -\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
    2.39 -
    2.40 -The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
    2.41 -zero or negative in order to have a feasible solution (since the sum
    2.42 -of the expressions on the left-hand side of the inequalities is zero).
    2.43 -It means that the total demand must be greater or equal to the total
    2.44 -supply and all the supplies have to be carried out from the supply nodes,
    2.45 -but there could be demands that are not satisfied.
    2.46 -If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand
    2.47 -constraints have to be satisfied with equality, i.e. all demands
    2.48 -have to be satisfied and all supplies have to be used.
    2.49 -
    2.50 -If you need the opposite inequalities in the supply/demand constraints
    2.51 -(i.e. the total demand is less than the total supply and all the demands
    2.52 -have to be satisfied while there could be supplies that are not used),
    2.53 -then you could easily transform the problem to the above form by reversing
    2.54 -the direction of the arcs and taking the negative of the supply values
    2.55 -(e.g. using \ref ReverseDigraph and \ref NegMap adaptors).
    2.56 -However \ref NetworkSimplex algorithm also supports this form directly
    2.57 -for the sake of convenience.
    2.58 -
    2.59 -A feasible solution for this problem can be found using \ref Circulation.
    2.60 -
    2.61 -Note that the above formulation is actually more general than the usual
    2.62 -definition of the minimum cost flow problem, in which strict equalities
    2.63 -are required in the supply/demand contraints, i.e.
    2.64 -
    2.65 -\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) =
    2.66 -    sup(u) \quad \forall u\in V. \f]
    2.67 -
    2.68 -However if the sum of the supply values is zero, then these two problems
    2.69 -are equivalent. So if you need the equality form, you have to ensure this
    2.70 -additional contraint for the algorithms.
    2.71 -
    2.72 -The dual solution of the minimum cost flow problem is represented by node 
    2.73 -potentials \f$\pi: V\rightarrow\mathbf{Z}\f$.
    2.74 -An \f$f: A\rightarrow\mathbf{Z}\f$ feasible solution of the problem
    2.75 -is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$
    2.76 -node potentials the following \e complementary \e slackness optimality
    2.77 -conditions hold.
    2.78 -
    2.79 - - For all \f$uv\in A\f$ arcs:
    2.80 -   - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
    2.81 -   - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
    2.82 -   - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
    2.83 - - For all \f$u\in V\f$ nodes:
    2.84 -   - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
    2.85 -     then \f$\pi(u)=0\f$.
    2.86 - 
    2.87 -Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc
    2.88 -\f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e.
    2.89 -\f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f]
    2.90 -
    2.91 -All algorithms provide dual solution (node potentials) as well,
    2.92 -if an optimal flow is found.
    2.93 -
    2.94 -LEMON contains several algorithms for solving minimum cost flow problems.
    2.95 +LEMON contains several algorithms for this problem.
    2.96   - \ref NetworkSimplex Primal Network Simplex algorithm with various
    2.97     pivot strategies.
    2.98   - \ref CostScaling Push-Relabel and Augment-Relabel algorithms based on
    2.99 @@ -429,10 +354,6 @@
   2.100   - \ref CancelAndTighten The Cancel and Tighten algorithm.
   2.101   - \ref CycleCanceling Cycle-Canceling algorithms.
   2.102  
   2.103 -Most of these implementations support the general inequality form of the
   2.104 -minimum cost flow problem, but CancelAndTighten and CycleCanceling
   2.105 -only support the equality form due to the primal method they use.
   2.106 -
   2.107  In general NetworkSimplex is the most efficient implementation,
   2.108  but in special cases other algorithms could be faster.
   2.109  For example, if the total supply and/or capacities are rather small,
     3.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     3.2 +++ b/doc/min_cost_flow.dox	Tue May 12 12:06:40 2009 +0200
     3.3 @@ -0,0 +1,153 @@
     3.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
     3.5 + *
     3.6 + * This file is a part of LEMON, a generic C++ optimization library.
     3.7 + *
     3.8 + * Copyright (C) 2003-2009
     3.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    3.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
    3.11 + *
    3.12 + * Permission to use, modify and distribute this software is granted
    3.13 + * provided that this copyright notice appears in all copies. For
    3.14 + * precise terms see the accompanying LICENSE file.
    3.15 + *
    3.16 + * This software is provided "AS IS" with no warranty of any kind,
    3.17 + * express or implied, and with no claim as to its suitability for any
    3.18 + * purpose.
    3.19 + *
    3.20 + */
    3.21 +
    3.22 +namespace lemon {
    3.23 +
    3.24 +/**
    3.25 +\page min_cost_flow Minimum Cost Flow Problem
    3.26 +
    3.27 +\section mcf_def Definition (GEQ form)
    3.28 +
    3.29 +The \e minimum \e cost \e flow \e problem is to find a feasible flow of
    3.30 +minimum total cost from a set of supply nodes to a set of demand nodes
    3.31 +in a network with capacity constraints (lower and upper bounds)
    3.32 +and arc costs.
    3.33 +
    3.34 +Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{R}\f$,
    3.35 +\f$upper: A\rightarrow\mathbf{R}\cup\{+\infty\}\f$ denote the lower and
    3.36 +upper bounds for the flow values on the arcs, for which
    3.37 +\f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$,
    3.38 +\f$cost: A\rightarrow\mathbf{R}\f$ denotes the cost per unit flow
    3.39 +on the arcs and \f$sup: V\rightarrow\mathbf{R}\f$ denotes the
    3.40 +signed supply values of the nodes.
    3.41 +If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
    3.42 +supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
    3.43 +\f$-sup(u)\f$ demand.
    3.44 +A minimum cost flow is an \f$f: A\rightarrow\mathbf{R}\f$ solution
    3.45 +of the following optimization problem.
    3.46 +
    3.47 +\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
    3.48 +\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
    3.49 +    sup(u) \quad \forall u\in V \f]
    3.50 +\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
    3.51 +
    3.52 +The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
    3.53 +zero or negative in order to have a feasible solution (since the sum
    3.54 +of the expressions on the left-hand side of the inequalities is zero).
    3.55 +It means that the total demand must be greater or equal to the total
    3.56 +supply and all the supplies have to be carried out from the supply nodes,
    3.57 +but there could be demands that are not satisfied.
    3.58 +If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand
    3.59 +constraints have to be satisfied with equality, i.e. all demands
    3.60 +have to be satisfied and all supplies have to be used.
    3.61 +
    3.62 +
    3.63 +\section mcf_algs Algorithms
    3.64 +
    3.65 +LEMON contains several algorithms for solving this problem, for more
    3.66 +information see \ref min_cost_flow_algs "Minimum Cost Flow Algorithms".
    3.67 +
    3.68 +A feasible solution for this problem can be found using \ref Circulation.
    3.69 +
    3.70 +
    3.71 +\section mcf_dual Dual Solution
    3.72 +
    3.73 +The dual solution of the minimum cost flow problem is represented by
    3.74 +node potentials \f$\pi: V\rightarrow\mathbf{R}\f$.
    3.75 +An \f$f: A\rightarrow\mathbf{R}\f$ primal feasible solution is optimal
    3.76 +if and only if for some \f$\pi: V\rightarrow\mathbf{R}\f$ node potentials
    3.77 +the following \e complementary \e slackness optimality conditions hold.
    3.78 +
    3.79 + - For all \f$uv\in A\f$ arcs:
    3.80 +   - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
    3.81 +   - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
    3.82 +   - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
    3.83 + - For all \f$u\in V\f$ nodes:
    3.84 +   - \f$\pi(u)<=0\f$;
    3.85 +   - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
    3.86 +     then \f$\pi(u)=0\f$.
    3.87 + 
    3.88 +Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc
    3.89 +\f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e.
    3.90 +\f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f]
    3.91 +
    3.92 +All algorithms provide dual solution (node potentials), as well,
    3.93 +if an optimal flow is found.
    3.94 +
    3.95 +
    3.96 +\section mcf_eq Equality Form
    3.97 +
    3.98 +The above \ref mcf_def "definition" is actually more general than the
    3.99 +usual formulation of the minimum cost flow problem, in which strict
   3.100 +equalities are required in the supply/demand contraints.
   3.101 +
   3.102 +\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
   3.103 +\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) =
   3.104 +    sup(u) \quad \forall u\in V \f]
   3.105 +\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
   3.106 +
   3.107 +However if the sum of the supply values is zero, then these two problems
   3.108 +are equivalent.
   3.109 +The \ref min_cost_flow_algs "algorithms" in LEMON support the general
   3.110 +form, so if you need the equality form, you have to ensure this additional
   3.111 +contraint manually.
   3.112 +
   3.113 +
   3.114 +\section mcf_leq Opposite Inequalites (LEQ Form)
   3.115 +
   3.116 +Another possible definition of the minimum cost flow problem is
   3.117 +when there are <em>"less or equal"</em> (LEQ) supply/demand constraints,
   3.118 +instead of the <em>"greater or equal"</em> (GEQ) constraints.
   3.119 +
   3.120 +\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
   3.121 +\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
   3.122 +    sup(u) \quad \forall u\in V \f]
   3.123 +\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
   3.124 +
   3.125 +It means that the total demand must be less or equal to the 
   3.126 +total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
   3.127 +positive) and all the demands have to be satisfied, but there
   3.128 +could be supplies that are not carried out from the supply
   3.129 +nodes.
   3.130 +The equality form is also a special case of this form, of course.
   3.131 +
   3.132 +You could easily transform this case to the \ref mcf_def "GEQ form"
   3.133 +of the problem by reversing the direction of the arcs and taking the
   3.134 +negative of the supply values (e.g. using \ref ReverseDigraph and
   3.135 +\ref NegMap adaptors).
   3.136 +However \ref NetworkSimplex algorithm also supports this form directly
   3.137 +for the sake of convenience.
   3.138 +
   3.139 +Note that the optimality conditions for this supply constraint type are
   3.140 +slightly differ from the conditions that are discussed for the GEQ form,
   3.141 +namely the potentials have to be non-negative instead of non-positive.
   3.142 +An \f$f: A\rightarrow\mathbf{R}\f$ feasible solution of this problem
   3.143 +is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{R}\f$
   3.144 +node potentials the following conditions hold.
   3.145 +
   3.146 + - For all \f$uv\in A\f$ arcs:
   3.147 +   - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
   3.148 +   - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
   3.149 +   - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
   3.150 + - For all \f$u\in V\f$ nodes:
   3.151 +   - \f$\pi(u)>=0\f$;
   3.152 +   - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
   3.153 +     then \f$\pi(u)=0\f$.
   3.154 +
   3.155 +*/
   3.156 +}
     4.1 --- a/lemon/network_simplex.h	Mon May 11 16:38:21 2009 +0100
     4.2 +++ b/lemon/network_simplex.h	Tue May 12 12:06:40 2009 +0200
     4.3 @@ -19,7 +19,7 @@
     4.4  #ifndef LEMON_NETWORK_SIMPLEX_H
     4.5  #define LEMON_NETWORK_SIMPLEX_H
     4.6  
     4.7 -/// \ingroup min_cost_flow
     4.8 +/// \ingroup min_cost_flow_algs
     4.9  ///
    4.10  /// \file
    4.11  /// \brief Network Simplex algorithm for finding a minimum cost flow.
    4.12 @@ -33,7 +33,7 @@
    4.13  
    4.14  namespace lemon {
    4.15  
    4.16 -  /// \addtogroup min_cost_flow
    4.17 +  /// \addtogroup min_cost_flow_algs
    4.18    /// @{
    4.19  
    4.20    /// \brief Implementation of the primal Network Simplex algorithm
    4.21 @@ -102,50 +102,16 @@
    4.22      /// i.e. the direction of the inequalities in the supply/demand
    4.23      /// constraints of the \ref min_cost_flow "minimum cost flow problem".
    4.24      ///
    4.25 -    /// The default supply type is \c GEQ, since this form is supported
    4.26 -    /// by other minimum cost flow algorithms and the \ref Circulation
    4.27 -    /// algorithm, as well.
    4.28 -    /// The \c LEQ problem type can be selected using the \ref supplyType()
    4.29 -    /// function.
    4.30 -    ///
    4.31 -    /// Note that the equality form is a special case of both supply types.
    4.32 +    /// The default supply type is \c GEQ, the \c LEQ type can be
    4.33 +    /// selected using \ref supplyType().
    4.34 +    /// The equality form is a special case of both supply types.
    4.35      enum SupplyType {
    4.36 -
    4.37        /// This option means that there are <em>"greater or equal"</em>
    4.38 -      /// supply/demand constraints in the definition, i.e. the exact
    4.39 -      /// formulation of the problem is the following.
    4.40 -      /**
    4.41 -          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
    4.42 -          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
    4.43 -              sup(u) \quad \forall u\in V \f]
    4.44 -          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
    4.45 -      */
    4.46 -      /// It means that the total demand must be greater or equal to the 
    4.47 -      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
    4.48 -      /// negative) and all the supplies have to be carried out from 
    4.49 -      /// the supply nodes, but there could be demands that are not 
    4.50 -      /// satisfied.
    4.51 +      /// supply/demand constraints in the definition of the problem.
    4.52        GEQ,
    4.53 -      /// It is just an alias for the \c GEQ option.
    4.54 -      CARRY_SUPPLIES = GEQ,
    4.55 -
    4.56        /// This option means that there are <em>"less or equal"</em>
    4.57 -      /// supply/demand constraints in the definition, i.e. the exact
    4.58 -      /// formulation of the problem is the following.
    4.59 -      /**
    4.60 -          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
    4.61 -          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
    4.62 -              sup(u) \quad \forall u\in V \f]
    4.63 -          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
    4.64 -      */
    4.65 -      /// It means that the total demand must be less or equal to the 
    4.66 -      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
    4.67 -      /// positive) and all the demands have to be satisfied, but there
    4.68 -      /// could be supplies that are not carried out from the supply
    4.69 -      /// nodes.
    4.70 -      LEQ,
    4.71 -      /// It is just an alias for the \c LEQ option.
    4.72 -      SATISFY_DEMANDS = LEQ
    4.73 +      /// supply/demand constraints in the definition of the problem.
    4.74 +      LEQ
    4.75      };
    4.76      
    4.77      /// \brief Constants for selecting the pivot rule.
    4.78 @@ -215,6 +181,8 @@
    4.79      const GR &_graph;
    4.80      int _node_num;
    4.81      int _arc_num;
    4.82 +    int _all_arc_num;
    4.83 +    int _search_arc_num;
    4.84  
    4.85      // Parameters of the problem
    4.86      bool _have_lower;
    4.87 @@ -277,7 +245,7 @@
    4.88        const IntVector  &_state;
    4.89        const CostVector &_pi;
    4.90        int &_in_arc;
    4.91 -      int _arc_num;
    4.92 +      int _search_arc_num;
    4.93  
    4.94        // Pivot rule data
    4.95        int _next_arc;
    4.96 @@ -288,13 +256,14 @@
    4.97        FirstEligiblePivotRule(NetworkSimplex &ns) :
    4.98          _source(ns._source), _target(ns._target),
    4.99          _cost(ns._cost), _state(ns._state), _pi(ns._pi),
   4.100 -        _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
   4.101 +        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
   4.102 +        _next_arc(0)
   4.103        {}
   4.104  
   4.105        // Find next entering arc
   4.106        bool findEnteringArc() {
   4.107          Cost c;
   4.108 -        for (int e = _next_arc; e < _arc_num; ++e) {
   4.109 +        for (int e = _next_arc; e < _search_arc_num; ++e) {
   4.110            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   4.111            if (c < 0) {
   4.112              _in_arc = e;
   4.113 @@ -328,7 +297,7 @@
   4.114        const IntVector  &_state;
   4.115        const CostVector &_pi;
   4.116        int &_in_arc;
   4.117 -      int _arc_num;
   4.118 +      int _search_arc_num;
   4.119  
   4.120      public:
   4.121  
   4.122 @@ -336,13 +305,13 @@
   4.123        BestEligiblePivotRule(NetworkSimplex &ns) :
   4.124          _source(ns._source), _target(ns._target),
   4.125          _cost(ns._cost), _state(ns._state), _pi(ns._pi),
   4.126 -        _in_arc(ns.in_arc), _arc_num(ns._arc_num)
   4.127 +        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num)
   4.128        {}
   4.129  
   4.130        // Find next entering arc
   4.131        bool findEnteringArc() {
   4.132          Cost c, min = 0;
   4.133 -        for (int e = 0; e < _arc_num; ++e) {
   4.134 +        for (int e = 0; e < _search_arc_num; ++e) {
   4.135            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   4.136            if (c < min) {
   4.137              min = c;
   4.138 @@ -367,7 +336,7 @@
   4.139        const IntVector  &_state;
   4.140        const CostVector &_pi;
   4.141        int &_in_arc;
   4.142 -      int _arc_num;
   4.143 +      int _search_arc_num;
   4.144  
   4.145        // Pivot rule data
   4.146        int _block_size;
   4.147 @@ -379,14 +348,15 @@
   4.148        BlockSearchPivotRule(NetworkSimplex &ns) :
   4.149          _source(ns._source), _target(ns._target),
   4.150          _cost(ns._cost), _state(ns._state), _pi(ns._pi),
   4.151 -        _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
   4.152 +        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
   4.153 +        _next_arc(0)
   4.154        {
   4.155          // The main parameters of the pivot rule
   4.156 -        const double BLOCK_SIZE_FACTOR = 2.0;
   4.157 +        const double BLOCK_SIZE_FACTOR = 0.5;
   4.158          const int MIN_BLOCK_SIZE = 10;
   4.159  
   4.160          _block_size = std::max( int(BLOCK_SIZE_FACTOR *
   4.161 -                                    std::sqrt(double(_arc_num))),
   4.162 +                                    std::sqrt(double(_search_arc_num))),
   4.163                                  MIN_BLOCK_SIZE );
   4.164        }
   4.165  
   4.166 @@ -395,7 +365,7 @@
   4.167          Cost c, min = 0;
   4.168          int cnt = _block_size;
   4.169          int e, min_arc = _next_arc;
   4.170 -        for (e = _next_arc; e < _arc_num; ++e) {
   4.171 +        for (e = _next_arc; e < _search_arc_num; ++e) {
   4.172            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   4.173            if (c < min) {
   4.174              min = c;
   4.175 @@ -440,7 +410,7 @@
   4.176        const IntVector  &_state;
   4.177        const CostVector &_pi;
   4.178        int &_in_arc;
   4.179 -      int _arc_num;
   4.180 +      int _search_arc_num;
   4.181  
   4.182        // Pivot rule data
   4.183        IntVector _candidates;
   4.184 @@ -454,7 +424,8 @@
   4.185        CandidateListPivotRule(NetworkSimplex &ns) :
   4.186          _source(ns._source), _target(ns._target),
   4.187          _cost(ns._cost), _state(ns._state), _pi(ns._pi),
   4.188 -        _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
   4.189 +        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
   4.190 +        _next_arc(0)
   4.191        {
   4.192          // The main parameters of the pivot rule
   4.193          const double LIST_LENGTH_FACTOR = 1.0;
   4.194 @@ -463,7 +434,7 @@
   4.195          const int MIN_MINOR_LIMIT = 3;
   4.196  
   4.197          _list_length = std::max( int(LIST_LENGTH_FACTOR *
   4.198 -                                     std::sqrt(double(_arc_num))),
   4.199 +                                     std::sqrt(double(_search_arc_num))),
   4.200                                   MIN_LIST_LENGTH );
   4.201          _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
   4.202                                   MIN_MINOR_LIMIT );
   4.203 @@ -500,7 +471,7 @@
   4.204          // Major iteration: build a new candidate list
   4.205          min = 0;
   4.206          _curr_length = 0;
   4.207 -        for (e = _next_arc; e < _arc_num; ++e) {
   4.208 +        for (e = _next_arc; e < _search_arc_num; ++e) {
   4.209            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   4.210            if (c < 0) {
   4.211              _candidates[_curr_length++] = e;
   4.212 @@ -546,7 +517,7 @@
   4.213        const IntVector  &_state;
   4.214        const CostVector &_pi;
   4.215        int &_in_arc;
   4.216 -      int _arc_num;
   4.217 +      int _search_arc_num;
   4.218  
   4.219        // Pivot rule data
   4.220        int _block_size, _head_length, _curr_length;
   4.221 @@ -574,8 +545,8 @@
   4.222        AlteringListPivotRule(NetworkSimplex &ns) :
   4.223          _source(ns._source), _target(ns._target),
   4.224          _cost(ns._cost), _state(ns._state), _pi(ns._pi),
   4.225 -        _in_arc(ns.in_arc), _arc_num(ns._arc_num),
   4.226 -        _next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost)
   4.227 +        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
   4.228 +        _next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost)
   4.229        {
   4.230          // The main parameters of the pivot rule
   4.231          const double BLOCK_SIZE_FACTOR = 1.5;
   4.232 @@ -584,7 +555,7 @@
   4.233          const int MIN_HEAD_LENGTH = 3;
   4.234  
   4.235          _block_size = std::max( int(BLOCK_SIZE_FACTOR *
   4.236 -                                    std::sqrt(double(_arc_num))),
   4.237 +                                    std::sqrt(double(_search_arc_num))),
   4.238                                  MIN_BLOCK_SIZE );
   4.239          _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
   4.240                                   MIN_HEAD_LENGTH );
   4.241 @@ -610,7 +581,7 @@
   4.242          int last_arc = 0;
   4.243          int limit = _head_length;
   4.244  
   4.245 -        for (int e = _next_arc; e < _arc_num; ++e) {
   4.246 +        for (int e = _next_arc; e < _search_arc_num; ++e) {
   4.247            _cand_cost[e] = _state[e] *
   4.248              (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   4.249            if (_cand_cost[e] < 0) {
   4.250 @@ -678,17 +649,17 @@
   4.251        _node_num = countNodes(_graph);
   4.252        _arc_num = countArcs(_graph);
   4.253        int all_node_num = _node_num + 1;
   4.254 -      int all_arc_num = _arc_num + _node_num;
   4.255 +      int max_arc_num = _arc_num + 2 * _node_num;
   4.256  
   4.257 -      _source.resize(all_arc_num);
   4.258 -      _target.resize(all_arc_num);
   4.259 +      _source.resize(max_arc_num);
   4.260 +      _target.resize(max_arc_num);
   4.261  
   4.262 -      _lower.resize(all_arc_num);
   4.263 -      _upper.resize(all_arc_num);
   4.264 -      _cap.resize(all_arc_num);
   4.265 -      _cost.resize(all_arc_num);
   4.266 +      _lower.resize(_arc_num);
   4.267 +      _upper.resize(_arc_num);
   4.268 +      _cap.resize(max_arc_num);
   4.269 +      _cost.resize(max_arc_num);
   4.270        _supply.resize(all_node_num);
   4.271 -      _flow.resize(all_arc_num);
   4.272 +      _flow.resize(max_arc_num);
   4.273        _pi.resize(all_node_num);
   4.274  
   4.275        _parent.resize(all_node_num);
   4.276 @@ -698,7 +669,7 @@
   4.277        _rev_thread.resize(all_node_num);
   4.278        _succ_num.resize(all_node_num);
   4.279        _last_succ.resize(all_node_num);
   4.280 -      _state.resize(all_arc_num);
   4.281 +      _state.resize(max_arc_num);
   4.282  
   4.283        // Copy the graph (store the arcs in a mixed order)
   4.284        int i = 0;
   4.285 @@ -1069,7 +1040,7 @@
   4.286        // Initialize artifical cost
   4.287        Cost ART_COST;
   4.288        if (std::numeric_limits<Cost>::is_exact) {
   4.289 -        ART_COST = std::numeric_limits<Cost>::max() / 4 + 1;
   4.290 +        ART_COST = std::numeric_limits<Cost>::max() / 2 + 1;
   4.291        } else {
   4.292          ART_COST = std::numeric_limits<Cost>::min();
   4.293          for (int i = 0; i != _arc_num; ++i) {
   4.294 @@ -1093,29 +1064,121 @@
   4.295        _succ_num[_root] = _node_num + 1;
   4.296        _last_succ[_root] = _root - 1;
   4.297        _supply[_root] = -_sum_supply;
   4.298 -      _pi[_root] = _sum_supply < 0 ? -ART_COST : ART_COST;
   4.299 +      _pi[_root] = 0;
   4.300  
   4.301        // Add artificial arcs and initialize the spanning tree data structure
   4.302 -      for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
   4.303 -        _parent[u] = _root;
   4.304 -        _pred[u] = e;
   4.305 -        _thread[u] = u + 1;
   4.306 -        _rev_thread[u + 1] = u;
   4.307 -        _succ_num[u] = 1;
   4.308 -        _last_succ[u] = u;
   4.309 -        _cost[e] = ART_COST;
   4.310 -        _cap[e] = INF;
   4.311 -        _state[e] = STATE_TREE;
   4.312 -        if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {
   4.313 -          _flow[e] = _supply[u];
   4.314 -          _forward[u] = true;
   4.315 -          _pi[u] = -ART_COST + _pi[_root];
   4.316 -        } else {
   4.317 -          _flow[e] = -_supply[u];
   4.318 -          _forward[u] = false;
   4.319 -          _pi[u] = ART_COST + _pi[_root];
   4.320 +      if (_sum_supply == 0) {
   4.321 +        // EQ supply constraints
   4.322 +        _search_arc_num = _arc_num;
   4.323 +        _all_arc_num = _arc_num + _node_num;
   4.324 +        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
   4.325 +          _parent[u] = _root;
   4.326 +          _pred[u] = e;
   4.327 +          _thread[u] = u + 1;
   4.328 +          _rev_thread[u + 1] = u;
   4.329 +          _succ_num[u] = 1;
   4.330 +          _last_succ[u] = u;
   4.331 +          _cap[e] = INF;
   4.332 +          _state[e] = STATE_TREE;
   4.333 +          if (_supply[u] >= 0) {
   4.334 +            _forward[u] = true;
   4.335 +            _pi[u] = 0;
   4.336 +            _source[e] = u;
   4.337 +            _target[e] = _root;
   4.338 +            _flow[e] = _supply[u];
   4.339 +            _cost[e] = 0;
   4.340 +          } else {
   4.341 +            _forward[u] = false;
   4.342 +            _pi[u] = ART_COST;
   4.343 +            _source[e] = _root;
   4.344 +            _target[e] = u;
   4.345 +            _flow[e] = -_supply[u];
   4.346 +            _cost[e] = ART_COST;
   4.347 +          }
   4.348          }
   4.349        }
   4.350 +      else if (_sum_supply > 0) {
   4.351 +        // LEQ supply constraints
   4.352 +        _search_arc_num = _arc_num + _node_num;
   4.353 +        int f = _arc_num + _node_num;
   4.354 +        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
   4.355 +          _parent[u] = _root;
   4.356 +          _thread[u] = u + 1;
   4.357 +          _rev_thread[u + 1] = u;
   4.358 +          _succ_num[u] = 1;
   4.359 +          _last_succ[u] = u;
   4.360 +          if (_supply[u] >= 0) {
   4.361 +            _forward[u] = true;
   4.362 +            _pi[u] = 0;
   4.363 +            _pred[u] = e;
   4.364 +            _source[e] = u;
   4.365 +            _target[e] = _root;
   4.366 +            _cap[e] = INF;
   4.367 +            _flow[e] = _supply[u];
   4.368 +            _cost[e] = 0;
   4.369 +            _state[e] = STATE_TREE;
   4.370 +          } else {
   4.371 +            _forward[u] = false;
   4.372 +            _pi[u] = ART_COST;
   4.373 +            _pred[u] = f;
   4.374 +            _source[f] = _root;
   4.375 +            _target[f] = u;
   4.376 +            _cap[f] = INF;
   4.377 +            _flow[f] = -_supply[u];
   4.378 +            _cost[f] = ART_COST;
   4.379 +            _state[f] = STATE_TREE;
   4.380 +            _source[e] = u;
   4.381 +            _target[e] = _root;
   4.382 +            _cap[e] = INF;
   4.383 +            _flow[e] = 0;
   4.384 +            _cost[e] = 0;
   4.385 +            _state[e] = STATE_LOWER;
   4.386 +            ++f;
   4.387 +          }
   4.388 +        }
   4.389 +        _all_arc_num = f;
   4.390 +      }
   4.391 +      else {
   4.392 +        // GEQ supply constraints
   4.393 +        _search_arc_num = _arc_num + _node_num;
   4.394 +        int f = _arc_num + _node_num;
   4.395 +        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
   4.396 +          _parent[u] = _root;
   4.397 +          _thread[u] = u + 1;
   4.398 +          _rev_thread[u + 1] = u;
   4.399 +          _succ_num[u] = 1;
   4.400 +          _last_succ[u] = u;
   4.401 +          if (_supply[u] <= 0) {
   4.402 +            _forward[u] = false;
   4.403 +            _pi[u] = 0;
   4.404 +            _pred[u] = e;
   4.405 +            _source[e] = _root;
   4.406 +            _target[e] = u;
   4.407 +            _cap[e] = INF;
   4.408 +            _flow[e] = -_supply[u];
   4.409 +            _cost[e] = 0;
   4.410 +            _state[e] = STATE_TREE;
   4.411 +          } else {
   4.412 +            _forward[u] = true;
   4.413 +            _pi[u] = -ART_COST;
   4.414 +            _pred[u] = f;
   4.415 +            _source[f] = u;
   4.416 +            _target[f] = _root;
   4.417 +            _cap[f] = INF;
   4.418 +            _flow[f] = _supply[u];
   4.419 +            _state[f] = STATE_TREE;
   4.420 +            _cost[f] = ART_COST;
   4.421 +            _source[e] = _root;
   4.422 +            _target[e] = u;
   4.423 +            _cap[e] = INF;
   4.424 +            _flow[e] = 0;
   4.425 +            _cost[e] = 0;
   4.426 +            _state[e] = STATE_LOWER;
   4.427 +            ++f;
   4.428 +          }
   4.429 +        }
   4.430 +        _all_arc_num = f;
   4.431 +      }
   4.432  
   4.433        return true;
   4.434      }
   4.435 @@ -1374,20 +1437,8 @@
   4.436        }
   4.437        
   4.438        // Check feasibility
   4.439 -      if (_sum_supply < 0) {
   4.440 -        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
   4.441 -          if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE;
   4.442 -        }
   4.443 -      }
   4.444 -      else if (_sum_supply > 0) {
   4.445 -        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
   4.446 -          if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE;
   4.447 -        }
   4.448 -      }
   4.449 -      else {
   4.450 -        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
   4.451 -          if (_flow[e] != 0) return INFEASIBLE;
   4.452 -        }
   4.453 +      for (int e = _search_arc_num; e != _all_arc_num; ++e) {
   4.454 +        if (_flow[e] != 0) return INFEASIBLE;
   4.455        }
   4.456  
   4.457        // Transform the solution and the supply map to the original form
   4.458 @@ -1401,6 +1452,30 @@
   4.459            }
   4.460          }
   4.461        }
   4.462 +      
   4.463 +      // Shift potentials to meet the requirements of the GEQ/LEQ type
   4.464 +      // optimality conditions
   4.465 +      if (_sum_supply == 0) {
   4.466 +        if (_stype == GEQ) {
   4.467 +          Cost max_pot = std::numeric_limits<Cost>::min();
   4.468 +          for (int i = 0; i != _node_num; ++i) {
   4.469 +            if (_pi[i] > max_pot) max_pot = _pi[i];
   4.470 +          }
   4.471 +          if (max_pot > 0) {
   4.472 +            for (int i = 0; i != _node_num; ++i)
   4.473 +              _pi[i] -= max_pot;
   4.474 +          }
   4.475 +        } else {
   4.476 +          Cost min_pot = std::numeric_limits<Cost>::max();
   4.477 +          for (int i = 0; i != _node_num; ++i) {
   4.478 +            if (_pi[i] < min_pot) min_pot = _pi[i];
   4.479 +          }
   4.480 +          if (min_pot < 0) {
   4.481 +            for (int i = 0; i != _node_num; ++i)
   4.482 +              _pi[i] -= min_pot;
   4.483 +          }
   4.484 +        }
   4.485 +      }
   4.486  
   4.487        return OPTIMAL;
   4.488      }