Port NetworkSimplex from SVN -r3520 (#234)
authorPeter Kovacs <kpeter@inf.elte.hu>
Tue, 24 Feb 2009 09:46:02 +0100
changeset 593e8349c6f12ca
parent 522 6a17a722b50e
child 594 a79ef774fae1
Port NetworkSimplex from SVN -r3520 (#234)
lemon/Makefile.am
lemon/network_simplex.h
test/CMakeLists.txt
test/Makefile.am
test/min_cost_flow_test.cc
     1.1 --- a/lemon/Makefile.am	Mon Feb 23 18:01:14 2009 +0000
     1.2 +++ b/lemon/Makefile.am	Tue Feb 24 09:46:02 2009 +0100
     1.3 @@ -85,6 +85,7 @@
     1.4  	lemon/max_matching.h \
     1.5  	lemon/min_cost_arborescence.h \
     1.6  	lemon/nauty_reader.h \
     1.7 +	lemon/network_simplex.h \
     1.8  	lemon/path.h \
     1.9  	lemon/preflow.h \
    1.10  	lemon/radix_sort.h \
     2.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     2.2 +++ b/lemon/network_simplex.h	Tue Feb 24 09:46:02 2009 +0100
     2.3 @@ -0,0 +1,1191 @@
     2.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2.5 + *
     2.6 + * This file is a part of LEMON, a generic C++ optimization library.
     2.7 + *
     2.8 + * Copyright (C) 2003-2009
     2.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    2.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
    2.11 + *
    2.12 + * Permission to use, modify and distribute this software is granted
    2.13 + * provided that this copyright notice appears in all copies. For
    2.14 + * precise terms see the accompanying LICENSE file.
    2.15 + *
    2.16 + * This software is provided "AS IS" with no warranty of any kind,
    2.17 + * express or implied, and with no claim as to its suitability for any
    2.18 + * purpose.
    2.19 + *
    2.20 + */
    2.21 +
    2.22 +#ifndef LEMON_NETWORK_SIMPLEX_H
    2.23 +#define LEMON_NETWORK_SIMPLEX_H
    2.24 +
    2.25 +/// \ingroup min_cost_flow
    2.26 +///
    2.27 +/// \file
    2.28 +/// \brief Network simplex algorithm for finding a minimum cost flow.
    2.29 +
    2.30 +#include <vector>
    2.31 +#include <limits>
    2.32 +#include <algorithm>
    2.33 +
    2.34 +#include <lemon/math.h>
    2.35 +
    2.36 +namespace lemon {
    2.37 +
    2.38 +  /// \addtogroup min_cost_flow
    2.39 +  /// @{
    2.40 +
    2.41 +  /// \brief Implementation of the primal network simplex algorithm
    2.42 +  /// for finding a \ref min_cost_flow "minimum cost flow".
    2.43 +  ///
    2.44 +  /// \ref NetworkSimplex implements the primal network simplex algorithm
    2.45 +  /// for finding a \ref min_cost_flow "minimum cost flow".
    2.46 +  ///
    2.47 +  /// \tparam Digraph The digraph type the algorithm runs on.
    2.48 +  /// \tparam LowerMap The type of the lower bound map.
    2.49 +  /// \tparam CapacityMap The type of the capacity (upper bound) map.
    2.50 +  /// \tparam CostMap The type of the cost (length) map.
    2.51 +  /// \tparam SupplyMap The type of the supply map.
    2.52 +  ///
    2.53 +  /// \warning
    2.54 +  /// - Arc capacities and costs should be \e non-negative \e integers.
    2.55 +  /// - Supply values should be \e signed \e integers.
    2.56 +  /// - The value types of the maps should be convertible to each other.
    2.57 +  /// - \c CostMap::Value must be signed type.
    2.58 +  ///
    2.59 +  /// \note \ref NetworkSimplex provides five different pivot rule
    2.60 +  /// implementations that significantly affect the efficiency of the
    2.61 +  /// algorithm.
    2.62 +  /// By default "Block Search" pivot rule is used, which proved to be
    2.63 +  /// by far the most efficient according to our benchmark tests.
    2.64 +  /// However another pivot rule can be selected using \ref run()
    2.65 +  /// function with the proper parameter.
    2.66 +#ifdef DOXYGEN
    2.67 +  template < typename Digraph,
    2.68 +             typename LowerMap,
    2.69 +             typename CapacityMap,
    2.70 +             typename CostMap,
    2.71 +             typename SupplyMap >
    2.72 +
    2.73 +#else
    2.74 +  template < typename Digraph,
    2.75 +             typename LowerMap = typename Digraph::template ArcMap<int>,
    2.76 +             typename CapacityMap = typename Digraph::template ArcMap<int>,
    2.77 +             typename CostMap = typename Digraph::template ArcMap<int>,
    2.78 +             typename SupplyMap = typename Digraph::template NodeMap<int> >
    2.79 +#endif
    2.80 +  class NetworkSimplex
    2.81 +  {
    2.82 +    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
    2.83 +
    2.84 +    typedef typename CapacityMap::Value Capacity;
    2.85 +    typedef typename CostMap::Value Cost;
    2.86 +    typedef typename SupplyMap::Value Supply;
    2.87 +
    2.88 +    typedef std::vector<Arc> ArcVector;
    2.89 +    typedef std::vector<Node> NodeVector;
    2.90 +    typedef std::vector<int> IntVector;
    2.91 +    typedef std::vector<bool> BoolVector;
    2.92 +    typedef std::vector<Capacity> CapacityVector;
    2.93 +    typedef std::vector<Cost> CostVector;
    2.94 +    typedef std::vector<Supply> SupplyVector;
    2.95 +
    2.96 +  public:
    2.97 +
    2.98 +    /// The type of the flow map
    2.99 +    typedef typename Digraph::template ArcMap<Capacity> FlowMap;
   2.100 +    /// The type of the potential map
   2.101 +    typedef typename Digraph::template NodeMap<Cost> PotentialMap;
   2.102 +
   2.103 +  public:
   2.104 +
   2.105 +    /// Enum type for selecting the pivot rule used by \ref run()
   2.106 +    enum PivotRuleEnum {
   2.107 +      FIRST_ELIGIBLE_PIVOT,
   2.108 +      BEST_ELIGIBLE_PIVOT,
   2.109 +      BLOCK_SEARCH_PIVOT,
   2.110 +      CANDIDATE_LIST_PIVOT,
   2.111 +      ALTERING_LIST_PIVOT
   2.112 +    };
   2.113 +
   2.114 +  private:
   2.115 +
   2.116 +    // State constants for arcs
   2.117 +    enum ArcStateEnum {
   2.118 +      STATE_UPPER = -1,
   2.119 +      STATE_TREE  =  0,
   2.120 +      STATE_LOWER =  1
   2.121 +    };
   2.122 +
   2.123 +  private:
   2.124 +
   2.125 +    // References for the original data
   2.126 +    const Digraph &_orig_graph;
   2.127 +    const LowerMap *_orig_lower;
   2.128 +    const CapacityMap &_orig_cap;
   2.129 +    const CostMap &_orig_cost;
   2.130 +    const SupplyMap *_orig_supply;
   2.131 +    Node _orig_source;
   2.132 +    Node _orig_target;
   2.133 +    Capacity _orig_flow_value;
   2.134 +
   2.135 +    // Result maps
   2.136 +    FlowMap *_flow_result;
   2.137 +    PotentialMap *_potential_result;
   2.138 +    bool _local_flow;
   2.139 +    bool _local_potential;
   2.140 +
   2.141 +    // Data structures for storing the graph
   2.142 +    ArcVector _arc;
   2.143 +    NodeVector _node;
   2.144 +    IntNodeMap _node_id;
   2.145 +    IntVector _source;
   2.146 +    IntVector _target;
   2.147 +
   2.148 +    // The number of nodes and arcs in the original graph
   2.149 +    int _node_num;
   2.150 +    int _arc_num;
   2.151 +
   2.152 +    // Node and arc maps
   2.153 +    CapacityVector _cap;
   2.154 +    CostVector _cost;
   2.155 +    CostVector _supply;
   2.156 +    CapacityVector _flow;
   2.157 +    CostVector _pi;
   2.158 +
   2.159 +    // Node and arc maps for the spanning tree structure
   2.160 +    IntVector _depth;
   2.161 +    IntVector _parent;
   2.162 +    IntVector _pred;
   2.163 +    IntVector _thread;
   2.164 +    BoolVector _forward;
   2.165 +    IntVector _state;
   2.166 +
   2.167 +    // The root node
   2.168 +    int _root;
   2.169 +
   2.170 +    // The entering arc in the current pivot iteration
   2.171 +    int _in_arc;
   2.172 +
   2.173 +    // Temporary data used in the current pivot iteration
   2.174 +    int join, u_in, v_in, u_out, v_out;
   2.175 +    int right, first, second, last;
   2.176 +    int stem, par_stem, new_stem;
   2.177 +    Capacity delta;
   2.178 +
   2.179 +  private:
   2.180 +
   2.181 +    /// \brief Implementation of the "First Eligible" pivot rule for the
   2.182 +    /// \ref NetworkSimplex "network simplex" algorithm.
   2.183 +    ///
   2.184 +    /// This class implements the "First Eligible" pivot rule
   2.185 +    /// for the \ref NetworkSimplex "network simplex" algorithm.
   2.186 +    ///
   2.187 +    /// For more information see \ref NetworkSimplex::run().
   2.188 +    class FirstEligiblePivotRule
   2.189 +    {
   2.190 +    private:
   2.191 +
   2.192 +      // References to the NetworkSimplex class
   2.193 +      const ArcVector &_arc;
   2.194 +      const IntVector  &_source;
   2.195 +      const IntVector  &_target;
   2.196 +      const CostVector &_cost;
   2.197 +      const IntVector  &_state;
   2.198 +      const CostVector &_pi;
   2.199 +      int &_in_arc;
   2.200 +      int _arc_num;
   2.201 +
   2.202 +      // Pivot rule data
   2.203 +      int _next_arc;
   2.204 +
   2.205 +    public:
   2.206 +
   2.207 +      /// Constructor
   2.208 +      FirstEligiblePivotRule(NetworkSimplex &ns) :
   2.209 +        _arc(ns._arc), _source(ns._source), _target(ns._target),
   2.210 +        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
   2.211 +        _in_arc(ns._in_arc), _arc_num(ns._arc_num), _next_arc(0)
   2.212 +      {}
   2.213 +
   2.214 +      /// Find next entering arc
   2.215 +      bool findEnteringArc() {
   2.216 +        Cost c;
   2.217 +        for (int e = _next_arc; e < _arc_num; ++e) {
   2.218 +          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   2.219 +          if (c < 0) {
   2.220 +            _in_arc = e;
   2.221 +            _next_arc = e + 1;
   2.222 +            return true;
   2.223 +          }
   2.224 +        }
   2.225 +        for (int e = 0; e < _next_arc; ++e) {
   2.226 +          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   2.227 +          if (c < 0) {
   2.228 +            _in_arc = e;
   2.229 +            _next_arc = e + 1;
   2.230 +            return true;
   2.231 +          }
   2.232 +        }
   2.233 +        return false;
   2.234 +      }
   2.235 +
   2.236 +    }; //class FirstEligiblePivotRule
   2.237 +
   2.238 +
   2.239 +    /// \brief Implementation of the "Best Eligible" pivot rule for the
   2.240 +    /// \ref NetworkSimplex "network simplex" algorithm.
   2.241 +    ///
   2.242 +    /// This class implements the "Best Eligible" pivot rule
   2.243 +    /// for the \ref NetworkSimplex "network simplex" algorithm.
   2.244 +    ///
   2.245 +    /// For more information see \ref NetworkSimplex::run().
   2.246 +    class BestEligiblePivotRule
   2.247 +    {
   2.248 +    private:
   2.249 +
   2.250 +      // References to the NetworkSimplex class
   2.251 +      const ArcVector &_arc;
   2.252 +      const IntVector  &_source;
   2.253 +      const IntVector  &_target;
   2.254 +      const CostVector &_cost;
   2.255 +      const IntVector  &_state;
   2.256 +      const CostVector &_pi;
   2.257 +      int &_in_arc;
   2.258 +      int _arc_num;
   2.259 +
   2.260 +    public:
   2.261 +
   2.262 +      /// Constructor
   2.263 +      BestEligiblePivotRule(NetworkSimplex &ns) :
   2.264 +        _arc(ns._arc), _source(ns._source), _target(ns._target),
   2.265 +        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
   2.266 +        _in_arc(ns._in_arc), _arc_num(ns._arc_num)
   2.267 +      {}
   2.268 +
   2.269 +      /// Find next entering arc
   2.270 +      bool findEnteringArc() {
   2.271 +        Cost c, min = 0;
   2.272 +        for (int e = 0; e < _arc_num; ++e) {
   2.273 +          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   2.274 +          if (c < min) {
   2.275 +            min = c;
   2.276 +            _in_arc = e;
   2.277 +          }
   2.278 +        }
   2.279 +        return min < 0;
   2.280 +      }
   2.281 +
   2.282 +    }; //class BestEligiblePivotRule
   2.283 +
   2.284 +
   2.285 +    /// \brief Implementation of the "Block Search" pivot rule for the
   2.286 +    /// \ref NetworkSimplex "network simplex" algorithm.
   2.287 +    ///
   2.288 +    /// This class implements the "Block Search" pivot rule
   2.289 +    /// for the \ref NetworkSimplex "network simplex" algorithm.
   2.290 +    ///
   2.291 +    /// For more information see \ref NetworkSimplex::run().
   2.292 +    class BlockSearchPivotRule
   2.293 +    {
   2.294 +    private:
   2.295 +
   2.296 +      // References to the NetworkSimplex class
   2.297 +      const ArcVector &_arc;
   2.298 +      const IntVector  &_source;
   2.299 +      const IntVector  &_target;
   2.300 +      const CostVector &_cost;
   2.301 +      const IntVector  &_state;
   2.302 +      const CostVector &_pi;
   2.303 +      int &_in_arc;
   2.304 +      int _arc_num;
   2.305 +
   2.306 +      // Pivot rule data
   2.307 +      int _block_size;
   2.308 +      int _next_arc;
   2.309 +
   2.310 +    public:
   2.311 +
   2.312 +      /// Constructor
   2.313 +      BlockSearchPivotRule(NetworkSimplex &ns) :
   2.314 +        _arc(ns._arc), _source(ns._source), _target(ns._target),
   2.315 +        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
   2.316 +        _in_arc(ns._in_arc), _arc_num(ns._arc_num), _next_arc(0)
   2.317 +      {
   2.318 +        // The main parameters of the pivot rule
   2.319 +        const double BLOCK_SIZE_FACTOR = 2.0;
   2.320 +        const int MIN_BLOCK_SIZE = 10;
   2.321 +
   2.322 +        _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_arc_num)),
   2.323 +                                MIN_BLOCK_SIZE );
   2.324 +      }
   2.325 +
   2.326 +      /// Find next entering arc
   2.327 +      bool findEnteringArc() {
   2.328 +        Cost c, min = 0;
   2.329 +        int cnt = _block_size;
   2.330 +        int e, min_arc = _next_arc;
   2.331 +        for (e = _next_arc; e < _arc_num; ++e) {
   2.332 +          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   2.333 +          if (c < min) {
   2.334 +            min = c;
   2.335 +            min_arc = e;
   2.336 +          }
   2.337 +          if (--cnt == 0) {
   2.338 +            if (min < 0) break;
   2.339 +            cnt = _block_size;
   2.340 +          }
   2.341 +        }
   2.342 +        if (min == 0 || cnt > 0) {
   2.343 +          for (e = 0; e < _next_arc; ++e) {
   2.344 +            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   2.345 +            if (c < min) {
   2.346 +              min = c;
   2.347 +              min_arc = e;
   2.348 +            }
   2.349 +            if (--cnt == 0) {
   2.350 +              if (min < 0) break;
   2.351 +              cnt = _block_size;
   2.352 +            }
   2.353 +          }
   2.354 +        }
   2.355 +        if (min >= 0) return false;
   2.356 +        _in_arc = min_arc;
   2.357 +        _next_arc = e;
   2.358 +        return true;
   2.359 +      }
   2.360 +
   2.361 +    }; //class BlockSearchPivotRule
   2.362 +
   2.363 +
   2.364 +    /// \brief Implementation of the "Candidate List" pivot rule for the
   2.365 +    /// \ref NetworkSimplex "network simplex" algorithm.
   2.366 +    ///
   2.367 +    /// This class implements the "Candidate List" pivot rule
   2.368 +    /// for the \ref NetworkSimplex "network simplex" algorithm.
   2.369 +    ///
   2.370 +    /// For more information see \ref NetworkSimplex::run().
   2.371 +    class CandidateListPivotRule
   2.372 +    {
   2.373 +    private:
   2.374 +
   2.375 +      // References to the NetworkSimplex class
   2.376 +      const ArcVector &_arc;
   2.377 +      const IntVector  &_source;
   2.378 +      const IntVector  &_target;
   2.379 +      const CostVector &_cost;
   2.380 +      const IntVector  &_state;
   2.381 +      const CostVector &_pi;
   2.382 +      int &_in_arc;
   2.383 +      int _arc_num;
   2.384 +
   2.385 +      // Pivot rule data
   2.386 +      IntVector _candidates;
   2.387 +      int _list_length, _minor_limit;
   2.388 +      int _curr_length, _minor_count;
   2.389 +      int _next_arc;
   2.390 +
   2.391 +    public:
   2.392 +
   2.393 +      /// Constructor
   2.394 +      CandidateListPivotRule(NetworkSimplex &ns) :
   2.395 +        _arc(ns._arc), _source(ns._source), _target(ns._target),
   2.396 +        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
   2.397 +        _in_arc(ns._in_arc), _arc_num(ns._arc_num), _next_arc(0)
   2.398 +      {
   2.399 +        // The main parameters of the pivot rule
   2.400 +        const double LIST_LENGTH_FACTOR = 1.0;
   2.401 +        const int MIN_LIST_LENGTH = 10;
   2.402 +        const double MINOR_LIMIT_FACTOR = 0.1;
   2.403 +        const int MIN_MINOR_LIMIT = 3;
   2.404 +
   2.405 +        _list_length = std::max( int(LIST_LENGTH_FACTOR * sqrt(_arc_num)),
   2.406 +                                 MIN_LIST_LENGTH );
   2.407 +        _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
   2.408 +                                 MIN_MINOR_LIMIT );
   2.409 +        _curr_length = _minor_count = 0;
   2.410 +        _candidates.resize(_list_length);
   2.411 +      }
   2.412 +
   2.413 +      /// Find next entering arc
   2.414 +      bool findEnteringArc() {
   2.415 +        Cost min, c;
   2.416 +        int e, min_arc = _next_arc;
   2.417 +        if (_curr_length > 0 && _minor_count < _minor_limit) {
   2.418 +          // Minor iteration: select the best eligible arc from the
   2.419 +          // current candidate list
   2.420 +          ++_minor_count;
   2.421 +          min = 0;
   2.422 +          for (int i = 0; i < _curr_length; ++i) {
   2.423 +            e = _candidates[i];
   2.424 +            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   2.425 +            if (c < min) {
   2.426 +              min = c;
   2.427 +              min_arc = e;
   2.428 +            }
   2.429 +            if (c >= 0) {
   2.430 +              _candidates[i--] = _candidates[--_curr_length];
   2.431 +            }
   2.432 +          }
   2.433 +          if (min < 0) {
   2.434 +            _in_arc = min_arc;
   2.435 +            return true;
   2.436 +          }
   2.437 +        }
   2.438 +
   2.439 +        // Major iteration: build a new candidate list
   2.440 +        min = 0;
   2.441 +        _curr_length = 0;
   2.442 +        for (e = _next_arc; e < _arc_num; ++e) {
   2.443 +          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   2.444 +          if (c < 0) {
   2.445 +            _candidates[_curr_length++] = e;
   2.446 +            if (c < min) {
   2.447 +              min = c;
   2.448 +              min_arc = e;
   2.449 +            }
   2.450 +            if (_curr_length == _list_length) break;
   2.451 +          }
   2.452 +        }
   2.453 +        if (_curr_length < _list_length) {
   2.454 +          for (e = 0; e < _next_arc; ++e) {
   2.455 +            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   2.456 +            if (c < 0) {
   2.457 +              _candidates[_curr_length++] = e;
   2.458 +              if (c < min) {
   2.459 +                min = c;
   2.460 +                min_arc = e;
   2.461 +              }
   2.462 +              if (_curr_length == _list_length) break;
   2.463 +            }
   2.464 +          }
   2.465 +        }
   2.466 +        if (_curr_length == 0) return false;
   2.467 +        _minor_count = 1;
   2.468 +        _in_arc = min_arc;
   2.469 +        _next_arc = e;
   2.470 +        return true;
   2.471 +      }
   2.472 +
   2.473 +    }; //class CandidateListPivotRule
   2.474 +
   2.475 +
   2.476 +    /// \brief Implementation of the "Altering Candidate List" pivot rule
   2.477 +    /// for the \ref NetworkSimplex "network simplex" algorithm.
   2.478 +    ///
   2.479 +    /// This class implements the "Altering Candidate List" pivot rule
   2.480 +    /// for the \ref NetworkSimplex "network simplex" algorithm.
   2.481 +    ///
   2.482 +    /// For more information see \ref NetworkSimplex::run().
   2.483 +    class AlteringListPivotRule
   2.484 +    {
   2.485 +    private:
   2.486 +
   2.487 +      // References to the NetworkSimplex class
   2.488 +      const ArcVector &_arc;
   2.489 +      const IntVector  &_source;
   2.490 +      const IntVector  &_target;
   2.491 +      const CostVector &_cost;
   2.492 +      const IntVector  &_state;
   2.493 +      const CostVector &_pi;
   2.494 +      int &_in_arc;
   2.495 +      int _arc_num;
   2.496 +
   2.497 +      // Pivot rule data
   2.498 +      int _block_size, _head_length, _curr_length;
   2.499 +      int _next_arc;
   2.500 +      IntVector _candidates;
   2.501 +      CostVector _cand_cost;
   2.502 +
   2.503 +      // Functor class to compare arcs during sort of the candidate list
   2.504 +      class SortFunc
   2.505 +      {
   2.506 +      private:
   2.507 +        const CostVector &_map;
   2.508 +      public:
   2.509 +        SortFunc(const CostVector &map) : _map(map) {}
   2.510 +        bool operator()(int left, int right) {
   2.511 +          return _map[left] > _map[right];
   2.512 +        }
   2.513 +      };
   2.514 +
   2.515 +      SortFunc _sort_func;
   2.516 +
   2.517 +    public:
   2.518 +
   2.519 +      /// Constructor
   2.520 +      AlteringListPivotRule(NetworkSimplex &ns) :
   2.521 +        _arc(ns._arc), _source(ns._source), _target(ns._target),
   2.522 +        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
   2.523 +        _in_arc(ns._in_arc), _arc_num(ns._arc_num),
   2.524 +        _next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost)
   2.525 +      {
   2.526 +        // The main parameters of the pivot rule
   2.527 +        const double BLOCK_SIZE_FACTOR = 1.5;
   2.528 +        const int MIN_BLOCK_SIZE = 10;
   2.529 +        const double HEAD_LENGTH_FACTOR = 0.1;
   2.530 +        const int MIN_HEAD_LENGTH = 3;
   2.531 +
   2.532 +        _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_arc_num)),
   2.533 +                                MIN_BLOCK_SIZE );
   2.534 +        _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
   2.535 +                                 MIN_HEAD_LENGTH );
   2.536 +        _candidates.resize(_head_length + _block_size);
   2.537 +        _curr_length = 0;
   2.538 +      }
   2.539 +
   2.540 +      /// Find next entering arc
   2.541 +      bool findEnteringArc() {
   2.542 +        // Check the current candidate list
   2.543 +        int e;
   2.544 +        for (int i = 0; i < _curr_length; ++i) {
   2.545 +          e = _candidates[i];
   2.546 +          _cand_cost[e] = _state[e] *
   2.547 +            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   2.548 +          if (_cand_cost[e] >= 0) {
   2.549 +            _candidates[i--] = _candidates[--_curr_length];
   2.550 +          }
   2.551 +        }
   2.552 +
   2.553 +        // Extend the list
   2.554 +        int cnt = _block_size;
   2.555 +        int last_edge = 0;
   2.556 +        int limit = _head_length;
   2.557 +
   2.558 +        for (int e = _next_arc; e < _arc_num; ++e) {
   2.559 +          _cand_cost[e] = _state[e] *
   2.560 +            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   2.561 +          if (_cand_cost[e] < 0) {
   2.562 +            _candidates[_curr_length++] = e;
   2.563 +            last_edge = e;
   2.564 +          }
   2.565 +          if (--cnt == 0) {
   2.566 +            if (_curr_length > limit) break;
   2.567 +            limit = 0;
   2.568 +            cnt = _block_size;
   2.569 +          }
   2.570 +        }
   2.571 +        if (_curr_length <= limit) {
   2.572 +          for (int e = 0; e < _next_arc; ++e) {
   2.573 +            _cand_cost[e] = _state[e] *
   2.574 +              (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
   2.575 +            if (_cand_cost[e] < 0) {
   2.576 +              _candidates[_curr_length++] = e;
   2.577 +              last_edge = e;
   2.578 +            }
   2.579 +            if (--cnt == 0) {
   2.580 +              if (_curr_length > limit) break;
   2.581 +              limit = 0;
   2.582 +              cnt = _block_size;
   2.583 +            }
   2.584 +          }
   2.585 +        }
   2.586 +        if (_curr_length == 0) return false;
   2.587 +        _next_arc = last_edge + 1;
   2.588 +
   2.589 +        // Make heap of the candidate list (approximating a partial sort)
   2.590 +        make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
   2.591 +                   _sort_func );
   2.592 +
   2.593 +        // Pop the first element of the heap
   2.594 +        _in_arc = _candidates[0];
   2.595 +        pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
   2.596 +                  _sort_func );
   2.597 +        _curr_length = std::min(_head_length, _curr_length - 1);
   2.598 +        return true;
   2.599 +      }
   2.600 +
   2.601 +    }; //class AlteringListPivotRule
   2.602 +
   2.603 +  public:
   2.604 +
   2.605 +    /// \brief General constructor (with lower bounds).
   2.606 +    ///
   2.607 +    /// General constructor (with lower bounds).
   2.608 +    ///
   2.609 +    /// \param digraph The digraph the algorithm runs on.
   2.610 +    /// \param lower The lower bounds of the arcs.
   2.611 +    /// \param capacity The capacities (upper bounds) of the arcs.
   2.612 +    /// \param cost The cost (length) values of the arcs.
   2.613 +    /// \param supply The supply values of the nodes (signed).
   2.614 +    NetworkSimplex( const Digraph &digraph,
   2.615 +                    const LowerMap &lower,
   2.616 +                    const CapacityMap &capacity,
   2.617 +                    const CostMap &cost,
   2.618 +                    const SupplyMap &supply ) :
   2.619 +      _orig_graph(digraph), _orig_lower(&lower), _orig_cap(capacity),
   2.620 +      _orig_cost(cost), _orig_supply(&supply),
   2.621 +      _flow_result(NULL), _potential_result(NULL),
   2.622 +      _local_flow(false), _local_potential(false),
   2.623 +      _node_id(digraph)
   2.624 +    {}
   2.625 +
   2.626 +    /// \brief General constructor (without lower bounds).
   2.627 +    ///
   2.628 +    /// General constructor (without lower bounds).
   2.629 +    ///
   2.630 +    /// \param digraph The digraph the algorithm runs on.
   2.631 +    /// \param capacity The capacities (upper bounds) of the arcs.
   2.632 +    /// \param cost The cost (length) values of the arcs.
   2.633 +    /// \param supply The supply values of the nodes (signed).
   2.634 +    NetworkSimplex( const Digraph &digraph,
   2.635 +                    const CapacityMap &capacity,
   2.636 +                    const CostMap &cost,
   2.637 +                    const SupplyMap &supply ) :
   2.638 +      _orig_graph(digraph), _orig_lower(NULL), _orig_cap(capacity),
   2.639 +      _orig_cost(cost), _orig_supply(&supply),
   2.640 +      _flow_result(NULL), _potential_result(NULL),
   2.641 +      _local_flow(false), _local_potential(false),
   2.642 +      _node_id(digraph)
   2.643 +    {}
   2.644 +
   2.645 +    /// \brief Simple constructor (with lower bounds).
   2.646 +    ///
   2.647 +    /// Simple constructor (with lower bounds).
   2.648 +    ///
   2.649 +    /// \param digraph The digraph the algorithm runs on.
   2.650 +    /// \param lower The lower bounds of the arcs.
   2.651 +    /// \param capacity The capacities (upper bounds) of the arcs.
   2.652 +    /// \param cost The cost (length) values of the arcs.
   2.653 +    /// \param s The source node.
   2.654 +    /// \param t The target node.
   2.655 +    /// \param flow_value The required amount of flow from node \c s
   2.656 +    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
   2.657 +    NetworkSimplex( const Digraph &digraph,
   2.658 +                    const LowerMap &lower,
   2.659 +                    const CapacityMap &capacity,
   2.660 +                    const CostMap &cost,
   2.661 +                    Node s, Node t,
   2.662 +                    Capacity flow_value ) :
   2.663 +      _orig_graph(digraph), _orig_lower(&lower), _orig_cap(capacity),
   2.664 +      _orig_cost(cost), _orig_supply(NULL),
   2.665 +      _orig_source(s), _orig_target(t), _orig_flow_value(flow_value),
   2.666 +      _flow_result(NULL), _potential_result(NULL),
   2.667 +      _local_flow(false), _local_potential(false),
   2.668 +      _node_id(digraph)
   2.669 +    {}
   2.670 +
   2.671 +    /// \brief Simple constructor (without lower bounds).
   2.672 +    ///
   2.673 +    /// Simple constructor (without lower bounds).
   2.674 +    ///
   2.675 +    /// \param digraph The digraph the algorithm runs on.
   2.676 +    /// \param capacity The capacities (upper bounds) of the arcs.
   2.677 +    /// \param cost The cost (length) values of the arcs.
   2.678 +    /// \param s The source node.
   2.679 +    /// \param t The target node.
   2.680 +    /// \param flow_value The required amount of flow from node \c s
   2.681 +    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
   2.682 +    NetworkSimplex( const Digraph &digraph,
   2.683 +                    const CapacityMap &capacity,
   2.684 +                    const CostMap &cost,
   2.685 +                    Node s, Node t,
   2.686 +                    Capacity flow_value ) :
   2.687 +      _orig_graph(digraph), _orig_lower(NULL), _orig_cap(capacity),
   2.688 +      _orig_cost(cost), _orig_supply(NULL),
   2.689 +      _orig_source(s), _orig_target(t), _orig_flow_value(flow_value),
   2.690 +      _flow_result(NULL), _potential_result(NULL),
   2.691 +      _local_flow(false), _local_potential(false),
   2.692 +      _node_id(digraph)
   2.693 +    {}
   2.694 +
   2.695 +    /// Destructor.
   2.696 +    ~NetworkSimplex() {
   2.697 +      if (_local_flow) delete _flow_result;
   2.698 +      if (_local_potential) delete _potential_result;
   2.699 +    }
   2.700 +
   2.701 +    /// \brief Set the flow map.
   2.702 +    ///
   2.703 +    /// This function sets the flow map.
   2.704 +    ///
   2.705 +    /// \return <tt>(*this)</tt>
   2.706 +    NetworkSimplex& flowMap(FlowMap &map) {
   2.707 +      if (_local_flow) {
   2.708 +        delete _flow_result;
   2.709 +        _local_flow = false;
   2.710 +      }
   2.711 +      _flow_result = &map;
   2.712 +      return *this;
   2.713 +    }
   2.714 +
   2.715 +    /// \brief Set the potential map.
   2.716 +    ///
   2.717 +    /// This function sets the potential map.
   2.718 +    ///
   2.719 +    /// \return <tt>(*this)</tt>
   2.720 +    NetworkSimplex& potentialMap(PotentialMap &map) {
   2.721 +      if (_local_potential) {
   2.722 +        delete _potential_result;
   2.723 +        _local_potential = false;
   2.724 +      }
   2.725 +      _potential_result = &map;
   2.726 +      return *this;
   2.727 +    }
   2.728 +
   2.729 +    /// \name Execution control
   2.730 +    /// The algorithm can be executed using the
   2.731 +    /// \ref lemon::NetworkSimplex::run() "run()" function.
   2.732 +    /// @{
   2.733 +
   2.734 +    /// \brief Run the algorithm.
   2.735 +    ///
   2.736 +    /// This function runs the algorithm.
   2.737 +    ///
   2.738 +    /// \param pivot_rule The pivot rule that is used during the
   2.739 +    /// algorithm.
   2.740 +    ///
   2.741 +    /// The available pivot rules:
   2.742 +    ///
   2.743 +    /// - FIRST_ELIGIBLE_PIVOT The next eligible arc is selected in
   2.744 +    /// a wraparound fashion in every iteration
   2.745 +    /// (\ref FirstEligiblePivotRule).
   2.746 +    ///
   2.747 +    /// - BEST_ELIGIBLE_PIVOT The best eligible arc is selected in
   2.748 +    /// every iteration (\ref BestEligiblePivotRule).
   2.749 +    ///
   2.750 +    /// - BLOCK_SEARCH_PIVOT A specified number of arcs are examined in
   2.751 +    /// every iteration in a wraparound fashion and the best eligible
   2.752 +    /// arc is selected from this block (\ref BlockSearchPivotRule).
   2.753 +    ///
   2.754 +    /// - CANDIDATE_LIST_PIVOT In a major iteration a candidate list is
   2.755 +    /// built from eligible arcs in a wraparound fashion and in the
   2.756 +    /// following minor iterations the best eligible arc is selected
   2.757 +    /// from this list (\ref CandidateListPivotRule).
   2.758 +    ///
   2.759 +    /// - ALTERING_LIST_PIVOT It is a modified version of the
   2.760 +    /// "Candidate List" pivot rule. It keeps only the several best
   2.761 +    /// eligible arcs from the former candidate list and extends this
   2.762 +    /// list in every iteration (\ref AlteringListPivotRule).
   2.763 +    ///
   2.764 +    /// According to our comprehensive benchmark tests the "Block Search"
   2.765 +    /// pivot rule proved to be the fastest and the most robust on
   2.766 +    /// various test inputs. Thus it is the default option.
   2.767 +    ///
   2.768 +    /// \return \c true if a feasible flow can be found.
   2.769 +    bool run(PivotRuleEnum pivot_rule = BLOCK_SEARCH_PIVOT) {
   2.770 +      return init() && start(pivot_rule);
   2.771 +    }
   2.772 +
   2.773 +    /// @}
   2.774 +
   2.775 +    /// \name Query Functions
   2.776 +    /// The results of the algorithm can be obtained using these
   2.777 +    /// functions.\n
   2.778 +    /// \ref lemon::NetworkSimplex::run() "run()" must be called before
   2.779 +    /// using them.
   2.780 +    /// @{
   2.781 +
   2.782 +    /// \brief Return a const reference to the flow map.
   2.783 +    ///
   2.784 +    /// This function returns a const reference to an arc map storing
   2.785 +    /// the found flow.
   2.786 +    ///
   2.787 +    /// \pre \ref run() must be called before using this function.
   2.788 +    const FlowMap& flowMap() const {
   2.789 +      return *_flow_result;
   2.790 +    }
   2.791 +
   2.792 +    /// \brief Return a const reference to the potential map
   2.793 +    /// (the dual solution).
   2.794 +    ///
   2.795 +    /// This function returns a const reference to a node map storing
   2.796 +    /// the found potentials (the dual solution).
   2.797 +    ///
   2.798 +    /// \pre \ref run() must be called before using this function.
   2.799 +    const PotentialMap& potentialMap() const {
   2.800 +      return *_potential_result;
   2.801 +    }
   2.802 +
   2.803 +    /// \brief Return the flow on the given arc.
   2.804 +    ///
   2.805 +    /// This function returns the flow on the given arc.
   2.806 +    ///
   2.807 +    /// \pre \ref run() must be called before using this function.
   2.808 +    Capacity flow(const Arc& arc) const {
   2.809 +      return (*_flow_result)[arc];
   2.810 +    }
   2.811 +
   2.812 +    /// \brief Return the potential of the given node.
   2.813 +    ///
   2.814 +    /// This function returns the potential of the given node.
   2.815 +    ///
   2.816 +    /// \pre \ref run() must be called before using this function.
   2.817 +    Cost potential(const Node& node) const {
   2.818 +      return (*_potential_result)[node];
   2.819 +    }
   2.820 +
   2.821 +    /// \brief Return the total cost of the found flow.
   2.822 +    ///
   2.823 +    /// This function returns the total cost of the found flow.
   2.824 +    /// The complexity of the function is \f$ O(e) \f$.
   2.825 +    ///
   2.826 +    /// \pre \ref run() must be called before using this function.
   2.827 +    Cost totalCost() const {
   2.828 +      Cost c = 0;
   2.829 +      for (ArcIt e(_orig_graph); e != INVALID; ++e)
   2.830 +        c += (*_flow_result)[e] * _orig_cost[e];
   2.831 +      return c;
   2.832 +    }
   2.833 +
   2.834 +    /// @}
   2.835 +
   2.836 +  private:
   2.837 +
   2.838 +    // Initialize internal data structures
   2.839 +    bool init() {
   2.840 +      // Initialize result maps
   2.841 +      if (!_flow_result) {
   2.842 +        _flow_result = new FlowMap(_orig_graph);
   2.843 +        _local_flow = true;
   2.844 +      }
   2.845 +      if (!_potential_result) {
   2.846 +        _potential_result = new PotentialMap(_orig_graph);
   2.847 +        _local_potential = true;
   2.848 +      }
   2.849 +
   2.850 +      // Initialize vectors
   2.851 +      _node_num = countNodes(_orig_graph);
   2.852 +      _arc_num = countArcs(_orig_graph);
   2.853 +      int all_node_num = _node_num + 1;
   2.854 +      int all_edge_num = _arc_num + _node_num;
   2.855 +
   2.856 +      _arc.resize(_arc_num);
   2.857 +      _node.reserve(_node_num);
   2.858 +      _source.resize(all_edge_num);
   2.859 +      _target.resize(all_edge_num);
   2.860 +
   2.861 +      _cap.resize(all_edge_num);
   2.862 +      _cost.resize(all_edge_num);
   2.863 +      _supply.resize(all_node_num);
   2.864 +      _flow.resize(all_edge_num, 0);
   2.865 +      _pi.resize(all_node_num, 0);
   2.866 +
   2.867 +      _depth.resize(all_node_num);
   2.868 +      _parent.resize(all_node_num);
   2.869 +      _pred.resize(all_node_num);
   2.870 +      _thread.resize(all_node_num);
   2.871 +      _forward.resize(all_node_num);
   2.872 +      _state.resize(all_edge_num, STATE_LOWER);
   2.873 +
   2.874 +      // Initialize node related data
   2.875 +      bool valid_supply = true;
   2.876 +      if (_orig_supply) {
   2.877 +        Supply sum = 0;
   2.878 +        int i = 0;
   2.879 +        for (NodeIt n(_orig_graph); n != INVALID; ++n, ++i) {
   2.880 +          _node.push_back(n);
   2.881 +          _node_id[n] = i;
   2.882 +          _supply[i] = (*_orig_supply)[n];
   2.883 +          sum += _supply[i];
   2.884 +        }
   2.885 +        valid_supply = (sum == 0);
   2.886 +      } else {
   2.887 +        int i = 0;
   2.888 +        for (NodeIt n(_orig_graph); n != INVALID; ++n, ++i) {
   2.889 +          _node.push_back(n);
   2.890 +          _node_id[n] = i;
   2.891 +          _supply[i] = 0;
   2.892 +        }
   2.893 +        _supply[_node_id[_orig_source]] =  _orig_flow_value;
   2.894 +        _supply[_node_id[_orig_target]] = -_orig_flow_value;
   2.895 +      }
   2.896 +      if (!valid_supply) return false;
   2.897 +
   2.898 +      // Set data for the artificial root node
   2.899 +      _root = _node_num;
   2.900 +      _depth[_root] = 0;
   2.901 +      _parent[_root] = -1;
   2.902 +      _pred[_root] = -1;
   2.903 +      _thread[_root] = 0;
   2.904 +      _supply[_root] = 0;
   2.905 +      _pi[_root] = 0;
   2.906 +
   2.907 +      // Store the arcs in a mixed order
   2.908 +      int k = std::max(int(sqrt(_arc_num)), 10);
   2.909 +      int i = 0;
   2.910 +      for (ArcIt e(_orig_graph); e != INVALID; ++e) {
   2.911 +        _arc[i] = e;
   2.912 +        if ((i += k) >= _arc_num) i = (i % k) + 1;
   2.913 +      }
   2.914 +
   2.915 +      // Initialize arc maps
   2.916 +      for (int i = 0; i != _arc_num; ++i) {
   2.917 +        Arc e = _arc[i];
   2.918 +        _source[i] = _node_id[_orig_graph.source(e)];
   2.919 +        _target[i] = _node_id[_orig_graph.target(e)];
   2.920 +        _cost[i] = _orig_cost[e];
   2.921 +        _cap[i] = _orig_cap[e];
   2.922 +      }
   2.923 +
   2.924 +      // Remove non-zero lower bounds
   2.925 +      if (_orig_lower) {
   2.926 +        for (int i = 0; i != _arc_num; ++i) {
   2.927 +          Capacity c = (*_orig_lower)[_arc[i]];
   2.928 +          if (c != 0) {
   2.929 +            _cap[i] -= c;
   2.930 +            _supply[_source[i]] -= c;
   2.931 +            _supply[_target[i]] += c;
   2.932 +          }
   2.933 +        }
   2.934 +      }
   2.935 +
   2.936 +      // Add artificial arcs and initialize the spanning tree data structure
   2.937 +      Cost max_cost = std::numeric_limits<Cost>::max() / 4;
   2.938 +      Capacity max_cap = std::numeric_limits<Capacity>::max();
   2.939 +      for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
   2.940 +        _thread[u] = u + 1;
   2.941 +        _depth[u] = 1;
   2.942 +        _parent[u] = _root;
   2.943 +        _pred[u] = e;
   2.944 +        if (_supply[u] >= 0) {
   2.945 +          _flow[e] = _supply[u];
   2.946 +          _forward[u] = true;
   2.947 +          _pi[u] = -max_cost;
   2.948 +        } else {
   2.949 +          _flow[e] = -_supply[u];
   2.950 +          _forward[u] = false;
   2.951 +          _pi[u] = max_cost;
   2.952 +        }
   2.953 +        _cost[e] = max_cost;
   2.954 +        _cap[e] = max_cap;
   2.955 +        _state[e] = STATE_TREE;
   2.956 +      }
   2.957 +
   2.958 +      return true;
   2.959 +    }
   2.960 +
   2.961 +    // Find the join node
   2.962 +    void findJoinNode() {
   2.963 +      int u = _source[_in_arc];
   2.964 +      int v = _target[_in_arc];
   2.965 +      while (_depth[u] > _depth[v]) u = _parent[u];
   2.966 +      while (_depth[v] > _depth[u]) v = _parent[v];
   2.967 +      while (u != v) {
   2.968 +        u = _parent[u];
   2.969 +        v = _parent[v];
   2.970 +      }
   2.971 +      join = u;
   2.972 +    }
   2.973 +
   2.974 +    // Find the leaving arc of the cycle and returns true if the
   2.975 +    // leaving arc is not the same as the entering arc
   2.976 +    bool findLeavingArc() {
   2.977 +      // Initialize first and second nodes according to the direction
   2.978 +      // of the cycle
   2.979 +      if (_state[_in_arc] == STATE_LOWER) {
   2.980 +        first  = _source[_in_arc];
   2.981 +        second = _target[_in_arc];
   2.982 +      } else {
   2.983 +        first  = _target[_in_arc];
   2.984 +        second = _source[_in_arc];
   2.985 +      }
   2.986 +      delta = _cap[_in_arc];
   2.987 +      int result = 0;
   2.988 +      Capacity d;
   2.989 +      int e;
   2.990 +
   2.991 +      // Search the cycle along the path form the first node to the root
   2.992 +      for (int u = first; u != join; u = _parent[u]) {
   2.993 +        e = _pred[u];
   2.994 +        d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
   2.995 +        if (d < delta) {
   2.996 +          delta = d;
   2.997 +          u_out = u;
   2.998 +          result = 1;
   2.999 +        }
  2.1000 +      }
  2.1001 +      // Search the cycle along the path form the second node to the root
  2.1002 +      for (int u = second; u != join; u = _parent[u]) {
  2.1003 +        e = _pred[u];
  2.1004 +        d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
  2.1005 +        if (d <= delta) {
  2.1006 +          delta = d;
  2.1007 +          u_out = u;
  2.1008 +          result = 2;
  2.1009 +        }
  2.1010 +      }
  2.1011 +
  2.1012 +      if (result == 1) {
  2.1013 +        u_in = first;
  2.1014 +        v_in = second;
  2.1015 +      } else {
  2.1016 +        u_in = second;
  2.1017 +        v_in = first;
  2.1018 +      }
  2.1019 +      return result != 0;
  2.1020 +    }
  2.1021 +
  2.1022 +    // Change _flow and _state vectors
  2.1023 +    void changeFlow(bool change) {
  2.1024 +      // Augment along the cycle
  2.1025 +      if (delta > 0) {
  2.1026 +        Capacity val = _state[_in_arc] * delta;
  2.1027 +        _flow[_in_arc] += val;
  2.1028 +        for (int u = _source[_in_arc]; u != join; u = _parent[u]) {
  2.1029 +          _flow[_pred[u]] += _forward[u] ? -val : val;
  2.1030 +        }
  2.1031 +        for (int u = _target[_in_arc]; u != join; u = _parent[u]) {
  2.1032 +          _flow[_pred[u]] += _forward[u] ? val : -val;
  2.1033 +        }
  2.1034 +      }
  2.1035 +      // Update the state of the entering and leaving arcs
  2.1036 +      if (change) {
  2.1037 +        _state[_in_arc] = STATE_TREE;
  2.1038 +        _state[_pred[u_out]] =
  2.1039 +          (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
  2.1040 +      } else {
  2.1041 +        _state[_in_arc] = -_state[_in_arc];
  2.1042 +      }
  2.1043 +    }
  2.1044 +
  2.1045 +    // Update _thread and _parent vectors
  2.1046 +    void updateThreadParent() {
  2.1047 +      int u;
  2.1048 +      v_out = _parent[u_out];
  2.1049 +
  2.1050 +      // Handle the case when join and v_out coincide
  2.1051 +      bool par_first = false;
  2.1052 +      if (join == v_out) {
  2.1053 +        for (u = join; u != u_in && u != v_in; u = _thread[u]) ;
  2.1054 +        if (u == v_in) {
  2.1055 +          par_first = true;
  2.1056 +          while (_thread[u] != u_out) u = _thread[u];
  2.1057 +          first = u;
  2.1058 +        }
  2.1059 +      }
  2.1060 +
  2.1061 +      // Find the last successor of u_in (u) and the node after it (right)
  2.1062 +      // according to the thread index
  2.1063 +      for (u = u_in; _depth[_thread[u]] > _depth[u_in]; u = _thread[u]) ;
  2.1064 +      right = _thread[u];
  2.1065 +      if (_thread[v_in] == u_out) {
  2.1066 +        for (last = u; _depth[last] > _depth[u_out]; last = _thread[last]) ;
  2.1067 +        if (last == u_out) last = _thread[last];
  2.1068 +      }
  2.1069 +      else last = _thread[v_in];
  2.1070 +
  2.1071 +      // Update stem nodes
  2.1072 +      _thread[v_in] = stem = u_in;
  2.1073 +      par_stem = v_in;
  2.1074 +      while (stem != u_out) {
  2.1075 +        _thread[u] = new_stem = _parent[stem];
  2.1076 +
  2.1077 +        // Find the node just before the stem node (u) according to
  2.1078 +        // the original thread index
  2.1079 +        for (u = new_stem; _thread[u] != stem; u = _thread[u]) ;
  2.1080 +        _thread[u] = right;
  2.1081 +
  2.1082 +        // Change the parent node of stem and shift stem and par_stem nodes
  2.1083 +        _parent[stem] = par_stem;
  2.1084 +        par_stem = stem;
  2.1085 +        stem = new_stem;
  2.1086 +
  2.1087 +        // Find the last successor of stem (u) and the node after it (right)
  2.1088 +        // according to the thread index
  2.1089 +        for (u = stem; _depth[_thread[u]] > _depth[stem]; u = _thread[u]) ;
  2.1090 +        right = _thread[u];
  2.1091 +      }
  2.1092 +      _parent[u_out] = par_stem;
  2.1093 +      _thread[u] = last;
  2.1094 +
  2.1095 +      if (join == v_out && par_first) {
  2.1096 +        if (first != v_in) _thread[first] = right;
  2.1097 +      } else {
  2.1098 +        for (u = v_out; _thread[u] != u_out; u = _thread[u]) ;
  2.1099 +        _thread[u] = right;
  2.1100 +      }
  2.1101 +    }
  2.1102 +
  2.1103 +    // Update _pred and _forward vectors
  2.1104 +    void updatePredArc() {
  2.1105 +      int u = u_out, v;
  2.1106 +      while (u != u_in) {
  2.1107 +        v = _parent[u];
  2.1108 +        _pred[u] = _pred[v];
  2.1109 +        _forward[u] = !_forward[v];
  2.1110 +        u = v;
  2.1111 +      }
  2.1112 +      _pred[u_in] = _in_arc;
  2.1113 +      _forward[u_in] = (u_in == _source[_in_arc]);
  2.1114 +    }
  2.1115 +
  2.1116 +    // Update _depth and _potential vectors
  2.1117 +    void updateDepthPotential() {
  2.1118 +      _depth[u_in] = _depth[v_in] + 1;
  2.1119 +      Cost sigma = _forward[u_in] ?
  2.1120 +        _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
  2.1121 +        _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
  2.1122 +      _pi[u_in] += sigma;
  2.1123 +      for(int u = _thread[u_in]; _parent[u] != -1; u = _thread[u]) {
  2.1124 +        _depth[u] = _depth[_parent[u]] + 1;
  2.1125 +        if (_depth[u] <= _depth[u_in]) break;
  2.1126 +        _pi[u] += sigma;
  2.1127 +      }
  2.1128 +    }
  2.1129 +
  2.1130 +    // Execute the algorithm
  2.1131 +    bool start(PivotRuleEnum pivot_rule) {
  2.1132 +      // Select the pivot rule implementation
  2.1133 +      switch (pivot_rule) {
  2.1134 +        case FIRST_ELIGIBLE_PIVOT:
  2.1135 +          return start<FirstEligiblePivotRule>();
  2.1136 +        case BEST_ELIGIBLE_PIVOT:
  2.1137 +          return start<BestEligiblePivotRule>();
  2.1138 +        case BLOCK_SEARCH_PIVOT:
  2.1139 +          return start<BlockSearchPivotRule>();
  2.1140 +        case CANDIDATE_LIST_PIVOT:
  2.1141 +          return start<CandidateListPivotRule>();
  2.1142 +        case ALTERING_LIST_PIVOT:
  2.1143 +          return start<AlteringListPivotRule>();
  2.1144 +      }
  2.1145 +      return false;
  2.1146 +    }
  2.1147 +
  2.1148 +    template<class PivotRuleImplementation>
  2.1149 +    bool start() {
  2.1150 +      PivotRuleImplementation pivot(*this);
  2.1151 +
  2.1152 +      // Execute the network simplex algorithm
  2.1153 +      while (pivot.findEnteringArc()) {
  2.1154 +        findJoinNode();
  2.1155 +        bool change = findLeavingArc();
  2.1156 +        changeFlow(change);
  2.1157 +        if (change) {
  2.1158 +          updateThreadParent();
  2.1159 +          updatePredArc();
  2.1160 +          updateDepthPotential();
  2.1161 +        }
  2.1162 +      }
  2.1163 +
  2.1164 +      // Check if the flow amount equals zero on all the artificial arcs
  2.1165 +      for (int e = _arc_num; e != _arc_num + _node_num; ++e) {
  2.1166 +        if (_flow[e] > 0) return false;
  2.1167 +      }
  2.1168 +
  2.1169 +      // Copy flow values to _flow_result
  2.1170 +      if (_orig_lower) {
  2.1171 +        for (int i = 0; i != _arc_num; ++i) {
  2.1172 +          Arc e = _arc[i];
  2.1173 +          (*_flow_result)[e] = (*_orig_lower)[e] + _flow[i];
  2.1174 +        }
  2.1175 +      } else {
  2.1176 +        for (int i = 0; i != _arc_num; ++i) {
  2.1177 +          (*_flow_result)[_arc[i]] = _flow[i];
  2.1178 +        }
  2.1179 +      }
  2.1180 +      // Copy potential values to _potential_result
  2.1181 +      for (int i = 0; i != _node_num; ++i) {
  2.1182 +        (*_potential_result)[_node[i]] = _pi[i];
  2.1183 +      }
  2.1184 +
  2.1185 +      return true;
  2.1186 +    }
  2.1187 +
  2.1188 +  }; //class NetworkSimplex
  2.1189 +
  2.1190 +  ///@}
  2.1191 +
  2.1192 +} //namespace lemon
  2.1193 +
  2.1194 +#endif //LEMON_NETWORK_SIMPLEX_H
     3.1 --- a/test/CMakeLists.txt	Mon Feb 23 18:01:14 2009 +0000
     3.2 +++ b/test/CMakeLists.txt	Tue Feb 24 09:46:02 2009 +0100
     3.3 @@ -30,6 +30,7 @@
     3.4    maps_test
     3.5    max_matching_test
     3.6    min_cost_arborescence_test
     3.7 +  min_cost_flow_test
     3.8    path_test
     3.9    preflow_test
    3.10    radix_sort_test
     4.1 --- a/test/Makefile.am	Mon Feb 23 18:01:14 2009 +0000
     4.2 +++ b/test/Makefile.am	Tue Feb 24 09:46:02 2009 +0100
     4.3 @@ -26,6 +26,7 @@
     4.4  	test/maps_test \
     4.5  	test/max_matching_test \
     4.6  	test/min_cost_arborescence_test \
     4.7 +	test/min_cost_flow_test \
     4.8  	test/path_test \
     4.9  	test/preflow_test \
    4.10  	test/radix_sort_test \
    4.11 @@ -68,6 +69,7 @@
    4.12  test_mip_test_SOURCES = test/mip_test.cc
    4.13  test_max_matching_test_SOURCES = test/max_matching_test.cc
    4.14  test_min_cost_arborescence_test_SOURCES = test/min_cost_arborescence_test.cc
    4.15 +test_min_cost_flow_test_SOURCES = test/min_cost_flow_test.cc
    4.16  test_path_test_SOURCES = test/path_test.cc
    4.17  test_preflow_test_SOURCES = test/preflow_test.cc
    4.18  test_radix_sort_test_SOURCES = test/radix_sort_test.cc
     5.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     5.2 +++ b/test/min_cost_flow_test.cc	Tue Feb 24 09:46:02 2009 +0100
     5.3 @@ -0,0 +1,455 @@
     5.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
     5.5 + *
     5.6 + * This file is a part of LEMON, a generic C++ optimization library.
     5.7 + *
     5.8 + * Copyright (C) 2003-2009
     5.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    5.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
    5.11 + *
    5.12 + * Permission to use, modify and distribute this software is granted
    5.13 + * provided that this copyright notice appears in all copies. For
    5.14 + * precise terms see the accompanying LICENSE file.
    5.15 + *
    5.16 + * This software is provided "AS IS" with no warranty of any kind,
    5.17 + * express or implied, and with no claim as to its suitability for any
    5.18 + * purpose.
    5.19 + *
    5.20 + */
    5.21 +
    5.22 +#include <iostream>
    5.23 +#include <fstream>
    5.24 +
    5.25 +#include <lemon/list_graph.h>
    5.26 +#include <lemon/smart_graph.h>
    5.27 +#include <lemon/lgf_reader.h>
    5.28 +
    5.29 +//#include <lemon/cycle_canceling.h>
    5.30 +//#include <lemon/capacity_scaling.h>
    5.31 +//#include <lemon/cost_scaling.h>
    5.32 +#include <lemon/network_simplex.h>
    5.33 +//#include <lemon/min_cost_flow.h>
    5.34 +//#include <lemon/min_cost_max_flow.h>
    5.35 +
    5.36 +#include <lemon/concepts/digraph.h>
    5.37 +#include <lemon/concept_check.h>
    5.38 +
    5.39 +#include "test_tools.h"
    5.40 +
    5.41 +using namespace lemon;
    5.42 +
    5.43 +char test_lgf[] =
    5.44 +  "@nodes\n"
    5.45 +  "label  sup1 sup2 sup3\n"
    5.46 +  "    1    20   27    0\n"
    5.47 +  "    2    -4    0    0\n"
    5.48 +  "    3     0    0    0\n"
    5.49 +  "    4     0    0    0\n"
    5.50 +  "    5     9    0    0\n"
    5.51 +  "    6    -6    0    0\n"
    5.52 +  "    7     0    0    0\n"
    5.53 +  "    8     0    0    0\n"
    5.54 +  "    9     3    0    0\n"
    5.55 +  "   10    -2    0    0\n"
    5.56 +  "   11     0    0    0\n"
    5.57 +  "   12   -20  -27    0\n"
    5.58 +  "\n"
    5.59 +  "@arcs\n"
    5.60 +  "       cost  cap low1 low2\n"
    5.61 +  " 1  2    70   11    0    8\n"
    5.62 +  " 1  3   150    3    0    1\n"
    5.63 +  " 1  4    80   15    0    2\n"
    5.64 +  " 2  8    80   12    0    0\n"
    5.65 +  " 3  5   140    5    0    3\n"
    5.66 +  " 4  6    60   10    0    1\n"
    5.67 +  " 4  7    80    2    0    0\n"
    5.68 +  " 4  8   110    3    0    0\n"
    5.69 +  " 5  7    60   14    0    0\n"
    5.70 +  " 5 11   120   12    0    0\n"
    5.71 +  " 6  3     0    3    0    0\n"
    5.72 +  " 6  9   140    4    0    0\n"
    5.73 +  " 6 10    90    8    0    0\n"
    5.74 +  " 7  1    30    5    0    0\n"
    5.75 +  " 8 12    60   16    0    4\n"
    5.76 +  " 9 12    50    6    0    0\n"
    5.77 +  "10 12    70   13    0    5\n"
    5.78 +  "10  2   100    7    0    0\n"
    5.79 +  "10  7    60   10    0    0\n"
    5.80 +  "11 10    20   14    0    6\n"
    5.81 +  "12 11    30   10    0    0\n"
    5.82 +  "\n"
    5.83 +  "@attributes\n"
    5.84 +  "source 1\n"
    5.85 +  "target 12\n";
    5.86 +
    5.87 +
    5.88 +// Check the interface of an MCF algorithm
    5.89 +template <typename GR, typename Value>
    5.90 +class McfClassConcept
    5.91 +{
    5.92 +public:
    5.93 +
    5.94 +  template <typename MCF>
    5.95 +  struct Constraints {
    5.96 +    void constraints() {
    5.97 +      checkConcept<concepts::Digraph, GR>();
    5.98 +
    5.99 +      MCF mcf_test1(g, lower, upper, cost, sup);
   5.100 +      MCF mcf_test2(g, upper, cost, sup);
   5.101 +      MCF mcf_test3(g, lower, upper, cost, n, n, k);
   5.102 +      MCF mcf_test4(g, upper, cost, n, n, k);
   5.103 +
   5.104 +      // TODO: This part should be enabled and the next part
   5.105 +      // should be removed if map copying is supported
   5.106 +/*
   5.107 +      flow = mcf_test1.flowMap();
   5.108 +      mcf_test1.flowMap(flow);
   5.109 +
   5.110 +      pot = mcf_test1.potentialMap();
   5.111 +      mcf_test1.potentialMap(pot);
   5.112 +*/
   5.113 +/**/
   5.114 +      const typename MCF::FlowMap &fm =
   5.115 +        mcf_test1.flowMap();
   5.116 +      mcf_test1.flowMap(flow);
   5.117 +      const typename MCF::PotentialMap &pm =
   5.118 +        mcf_test1.potentialMap();
   5.119 +      mcf_test1.potentialMap(pot);
   5.120 +      ignore_unused_variable_warning(fm);
   5.121 +      ignore_unused_variable_warning(pm);
   5.122 +/**/
   5.123 +
   5.124 +      mcf_test1.run();
   5.125 +
   5.126 +      v = mcf_test1.totalCost();
   5.127 +      v = mcf_test1.flow(a);
   5.128 +      v = mcf_test1.potential(n);
   5.129 +    }
   5.130 +
   5.131 +    typedef typename GR::Node Node;
   5.132 +    typedef typename GR::Arc Arc;
   5.133 +    typedef concepts::ReadMap<Node, Value> NM;
   5.134 +    typedef concepts::ReadMap<Arc, Value> AM;
   5.135 +
   5.136 +    const GR &g;
   5.137 +    const AM &lower;
   5.138 +    const AM &upper;
   5.139 +    const AM &cost;
   5.140 +    const NM &sup;
   5.141 +    const Node &n;
   5.142 +    const Arc &a;
   5.143 +    const Value &k;
   5.144 +    Value v;
   5.145 +
   5.146 +    typename MCF::FlowMap &flow;
   5.147 +    typename MCF::PotentialMap &pot;
   5.148 +  };
   5.149 +
   5.150 +};
   5.151 +
   5.152 +
   5.153 +// Check the feasibility of the given flow (primal soluiton)
   5.154 +template < typename GR, typename LM, typename UM,
   5.155 +           typename SM, typename FM >
   5.156 +bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
   5.157 +                const SM& supply, const FM& flow )
   5.158 +{
   5.159 +  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   5.160 +
   5.161 +  for (ArcIt e(gr); e != INVALID; ++e) {
   5.162 +    if (flow[e] < lower[e] || flow[e] > upper[e]) return false;
   5.163 +  }
   5.164 +
   5.165 +  for (NodeIt n(gr); n != INVALID; ++n) {
   5.166 +    typename SM::Value sum = 0;
   5.167 +    for (OutArcIt e(gr, n); e != INVALID; ++e)
   5.168 +      sum += flow[e];
   5.169 +    for (InArcIt e(gr, n); e != INVALID; ++e)
   5.170 +      sum -= flow[e];
   5.171 +    if (sum != supply[n]) return false;
   5.172 +  }
   5.173 +
   5.174 +  return true;
   5.175 +}
   5.176 +
   5.177 +// Check the feasibility of the given potentials (dual soluiton)
   5.178 +// using the Complementary Slackness optimality condition
   5.179 +template < typename GR, typename LM, typename UM,
   5.180 +           typename CM, typename FM, typename PM >
   5.181 +bool checkPotential( const GR& gr, const LM& lower, const UM& upper,
   5.182 +                     const CM& cost, const FM& flow, const PM& pi )
   5.183 +{
   5.184 +  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   5.185 +
   5.186 +  bool opt = true;
   5.187 +  for (ArcIt e(gr); opt && e != INVALID; ++e) {
   5.188 +    typename CM::Value red_cost =
   5.189 +      cost[e] + pi[gr.source(e)] - pi[gr.target(e)];
   5.190 +    opt = red_cost == 0 ||
   5.191 +          (red_cost > 0 && flow[e] == lower[e]) ||
   5.192 +          (red_cost < 0 && flow[e] == upper[e]);
   5.193 +  }
   5.194 +  return opt;
   5.195 +}
   5.196 +
   5.197 +// Run a minimum cost flow algorithm and check the results
   5.198 +template < typename MCF, typename GR,
   5.199 +           typename LM, typename UM,
   5.200 +           typename CM, typename SM >
   5.201 +void checkMcf( const MCF& mcf, bool mcf_result,
   5.202 +               const GR& gr, const LM& lower, const UM& upper,
   5.203 +               const CM& cost, const SM& supply,
   5.204 +               bool result, typename CM::Value total,
   5.205 +               const std::string &test_id = "" )
   5.206 +{
   5.207 +  check(mcf_result == result, "Wrong result " + test_id);
   5.208 +  if (result) {
   5.209 +    check(checkFlow(gr, lower, upper, supply, mcf.flowMap()),
   5.210 +          "The flow is not feasible " + test_id);
   5.211 +    check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
   5.212 +    check(checkPotential(gr, lower, upper, cost, mcf.flowMap(),
   5.213 +                         mcf.potentialMap()),
   5.214 +          "Wrong potentials " + test_id);
   5.215 +  }
   5.216 +}
   5.217 +
   5.218 +int main()
   5.219 +{
   5.220 +  // Check the interfaces
   5.221 +  {
   5.222 +    typedef int Value;
   5.223 +    // This typedef should be enabled if the standard maps are
   5.224 +    // reference maps in the graph concepts
   5.225 +    //typedef concepts::Digraph GR;
   5.226 +    typedef ListDigraph GR;
   5.227 +    typedef concepts::ReadMap<GR::Node, Value> NM;
   5.228 +    typedef concepts::ReadMap<GR::Arc, Value> AM;
   5.229 +
   5.230 +    //checkConcept< McfClassConcept<GR, Value>,
   5.231 +    //              CycleCanceling<GR, AM, AM, AM, NM> >();
   5.232 +    //checkConcept< McfClassConcept<GR, Value>,
   5.233 +    //              CapacityScaling<GR, AM, AM, AM, NM> >();
   5.234 +    //checkConcept< McfClassConcept<GR, Value>,
   5.235 +    //              CostScaling<GR, AM, AM, AM, NM> >();
   5.236 +    checkConcept< McfClassConcept<GR, Value>,
   5.237 +                  NetworkSimplex<GR, AM, AM, AM, NM> >();
   5.238 +    //checkConcept< MinCostFlow<GR, Value>,
   5.239 +    //              NetworkSimplex<GR, AM, AM, AM, NM> >();
   5.240 +  }
   5.241 +
   5.242 +  // Run various MCF tests
   5.243 +  typedef ListDigraph Digraph;
   5.244 +  DIGRAPH_TYPEDEFS(ListDigraph);
   5.245 +
   5.246 +  // Read the test digraph
   5.247 +  Digraph gr;
   5.248 +  Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), u(gr);
   5.249 +  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr);
   5.250 +  Node v, w;
   5.251 +
   5.252 +  std::istringstream input(test_lgf);
   5.253 +  DigraphReader<Digraph>(gr, input)
   5.254 +    .arcMap("cost", c)
   5.255 +    .arcMap("cap", u)
   5.256 +    .arcMap("low1", l1)
   5.257 +    .arcMap("low2", l2)
   5.258 +    .nodeMap("sup1", s1)
   5.259 +    .nodeMap("sup2", s2)
   5.260 +    .nodeMap("sup3", s3)
   5.261 +    .node("source", v)
   5.262 +    .node("target", w)
   5.263 +    .run();
   5.264 +
   5.265 +/*
   5.266 +  // A. Test CapacityScaling with scaling
   5.267 +  {
   5.268 +    CapacityScaling<Digraph> mcf1(gr, u, c, s1);
   5.269 +    CapacityScaling<Digraph> mcf2(gr, u, c, v, w, 27);
   5.270 +    CapacityScaling<Digraph> mcf3(gr, u, c, s3);
   5.271 +    CapacityScaling<Digraph> mcf4(gr, l2, u, c, s1);
   5.272 +    CapacityScaling<Digraph> mcf5(gr, l2, u, c, v, w, 27);
   5.273 +    CapacityScaling<Digraph> mcf6(gr, l2, u, c, s3);
   5.274 +
   5.275 +    checkMcf(mcf1, mcf1.run(), gr, l1, u, c, s1, true,  5240, "#A1");
   5.276 +    checkMcf(mcf2, mcf2.run(), gr, l1, u, c, s2, true,  7620, "#A2");
   5.277 +    checkMcf(mcf3, mcf3.run(), gr, l1, u, c, s3, true,     0, "#A3");
   5.278 +    checkMcf(mcf4, mcf4.run(), gr, l2, u, c, s1, true,  5970, "#A4");
   5.279 +    checkMcf(mcf5, mcf5.run(), gr, l2, u, c, s2, true,  8010, "#A5");
   5.280 +    checkMcf(mcf6, mcf6.run(), gr, l2, u, c, s3, false,    0, "#A6");
   5.281 +  }
   5.282 +
   5.283 +  // B. Test CapacityScaling without scaling
   5.284 +  {
   5.285 +    CapacityScaling<Digraph> mcf1(gr, u, c, s1);
   5.286 +    CapacityScaling<Digraph> mcf2(gr, u, c, v, w, 27);
   5.287 +    CapacityScaling<Digraph> mcf3(gr, u, c, s3);
   5.288 +    CapacityScaling<Digraph> mcf4(gr, l2, u, c, s1);
   5.289 +    CapacityScaling<Digraph> mcf5(gr, l2, u, c, v, w, 27);
   5.290 +    CapacityScaling<Digraph> mcf6(gr, l2, u, c, s3);
   5.291 +
   5.292 +    checkMcf(mcf1, mcf1.run(false), gr, l1, u, c, s1, true,  5240, "#B1");
   5.293 +    checkMcf(mcf2, mcf2.run(false), gr, l1, u, c, s2, true,  7620, "#B2");
   5.294 +    checkMcf(mcf3, mcf3.run(false), gr, l1, u, c, s3, true,     0, "#B3");
   5.295 +    checkMcf(mcf4, mcf4.run(false), gr, l2, u, c, s1, true,  5970, "#B4");
   5.296 +    checkMcf(mcf5, mcf5.run(false), gr, l2, u, c, s2, true,  8010, "#B5");
   5.297 +    checkMcf(mcf6, mcf6.run(false), gr, l2, u, c, s3, false,    0, "#B6");
   5.298 +  }
   5.299 +
   5.300 +  // C. Test CostScaling using partial augment-relabel method
   5.301 +  {
   5.302 +    CostScaling<Digraph> mcf1(gr, u, c, s1);
   5.303 +    CostScaling<Digraph> mcf2(gr, u, c, v, w, 27);
   5.304 +    CostScaling<Digraph> mcf3(gr, u, c, s3);
   5.305 +    CostScaling<Digraph> mcf4(gr, l2, u, c, s1);
   5.306 +    CostScaling<Digraph> mcf5(gr, l2, u, c, v, w, 27);
   5.307 +    CostScaling<Digraph> mcf6(gr, l2, u, c, s3);
   5.308 +
   5.309 +    checkMcf(mcf1, mcf1.run(), gr, l1, u, c, s1, true,  5240, "#C1");
   5.310 +    checkMcf(mcf2, mcf2.run(), gr, l1, u, c, s2, true,  7620, "#C2");
   5.311 +    checkMcf(mcf3, mcf3.run(), gr, l1, u, c, s3, true,     0, "#C3");
   5.312 +    checkMcf(mcf4, mcf4.run(), gr, l2, u, c, s1, true,  5970, "#C4");
   5.313 +    checkMcf(mcf5, mcf5.run(), gr, l2, u, c, s2, true,  8010, "#C5");
   5.314 +    checkMcf(mcf6, mcf6.run(), gr, l2, u, c, s3, false,    0, "#C6");
   5.315 +  }
   5.316 +
   5.317 +  // D. Test CostScaling using push-relabel method
   5.318 +  {
   5.319 +    CostScaling<Digraph> mcf1(gr, u, c, s1);
   5.320 +    CostScaling<Digraph> mcf2(gr, u, c, v, w, 27);
   5.321 +    CostScaling<Digraph> mcf3(gr, u, c, s3);
   5.322 +    CostScaling<Digraph> mcf4(gr, l2, u, c, s1);
   5.323 +    CostScaling<Digraph> mcf5(gr, l2, u, c, v, w, 27);
   5.324 +    CostScaling<Digraph> mcf6(gr, l2, u, c, s3);
   5.325 +
   5.326 +    checkMcf(mcf1, mcf1.run(false), gr, l1, u, c, s1, true,  5240, "#D1");
   5.327 +    checkMcf(mcf2, mcf2.run(false), gr, l1, u, c, s2, true,  7620, "#D2");
   5.328 +    checkMcf(mcf3, mcf3.run(false), gr, l1, u, c, s3, true,     0, "#D3");
   5.329 +    checkMcf(mcf4, mcf4.run(false), gr, l2, u, c, s1, true,  5970, "#D4");
   5.330 +    checkMcf(mcf5, mcf5.run(false), gr, l2, u, c, s2, true,  8010, "#D5");
   5.331 +    checkMcf(mcf6, mcf6.run(false), gr, l2, u, c, s3, false,    0, "#D6");
   5.332 +  }
   5.333 +*/
   5.334 +
   5.335 +  // E. Test NetworkSimplex with FIRST_ELIGIBLE_PIVOT
   5.336 +  {
   5.337 +    NetworkSimplex<Digraph>::PivotRuleEnum pr =
   5.338 +      NetworkSimplex<Digraph>::FIRST_ELIGIBLE_PIVOT;
   5.339 +    NetworkSimplex<Digraph> mcf1(gr, u, c, s1);
   5.340 +    NetworkSimplex<Digraph> mcf2(gr, u, c, v, w, 27);
   5.341 +    NetworkSimplex<Digraph> mcf3(gr, u, c, s3);
   5.342 +    NetworkSimplex<Digraph> mcf4(gr, l2, u, c, s1);
   5.343 +    NetworkSimplex<Digraph> mcf5(gr, l2, u, c, v, w, 27);
   5.344 +    NetworkSimplex<Digraph> mcf6(gr, l2, u, c, s3);
   5.345 +
   5.346 +    checkMcf(mcf1, mcf1.run(pr), gr, l1, u, c, s1, true,  5240, "#E1");
   5.347 +    checkMcf(mcf2, mcf2.run(pr), gr, l1, u, c, s2, true,  7620, "#E2");
   5.348 +    checkMcf(mcf3, mcf3.run(pr), gr, l1, u, c, s3, true,     0, "#E3");
   5.349 +    checkMcf(mcf4, mcf4.run(pr), gr, l2, u, c, s1, true,  5970, "#E4");
   5.350 +    checkMcf(mcf5, mcf5.run(pr), gr, l2, u, c, s2, true,  8010, "#E5");
   5.351 +    checkMcf(mcf6, mcf6.run(pr), gr, l2, u, c, s3, false,    0, "#E6");
   5.352 +  }
   5.353 +
   5.354 +  // F. Test NetworkSimplex with BEST_ELIGIBLE_PIVOT
   5.355 +  {
   5.356 +    NetworkSimplex<Digraph>::PivotRuleEnum pr =
   5.357 +      NetworkSimplex<Digraph>::BEST_ELIGIBLE_PIVOT;
   5.358 +    NetworkSimplex<Digraph> mcf1(gr, u, c, s1);
   5.359 +    NetworkSimplex<Digraph> mcf2(gr, u, c, v, w, 27);
   5.360 +    NetworkSimplex<Digraph> mcf3(gr, u, c, s3);
   5.361 +    NetworkSimplex<Digraph> mcf4(gr, l2, u, c, s1);
   5.362 +    NetworkSimplex<Digraph> mcf5(gr, l2, u, c, v, w, 27);
   5.363 +    NetworkSimplex<Digraph> mcf6(gr, l2, u, c, s3);
   5.364 +
   5.365 +    checkMcf(mcf1, mcf1.run(pr), gr, l1, u, c, s1, true,  5240, "#F1");
   5.366 +    checkMcf(mcf2, mcf2.run(pr), gr, l1, u, c, s2, true,  7620, "#F2");
   5.367 +    checkMcf(mcf3, mcf3.run(pr), gr, l1, u, c, s3, true,     0, "#F3");
   5.368 +    checkMcf(mcf4, mcf4.run(pr), gr, l2, u, c, s1, true,  5970, "#F4");
   5.369 +    checkMcf(mcf5, mcf5.run(pr), gr, l2, u, c, s2, true,  8010, "#F5");
   5.370 +    checkMcf(mcf6, mcf6.run(pr), gr, l2, u, c, s3, false,    0, "#F6");
   5.371 +  }
   5.372 +
   5.373 +  // G. Test NetworkSimplex with BLOCK_SEARCH_PIVOT
   5.374 +  {
   5.375 +    NetworkSimplex<Digraph>::PivotRuleEnum pr =
   5.376 +      NetworkSimplex<Digraph>::BLOCK_SEARCH_PIVOT;
   5.377 +    NetworkSimplex<Digraph> mcf1(gr, u, c, s1);
   5.378 +    NetworkSimplex<Digraph> mcf2(gr, u, c, v, w, 27);
   5.379 +    NetworkSimplex<Digraph> mcf3(gr, u, c, s3);
   5.380 +    NetworkSimplex<Digraph> mcf4(gr, l2, u, c, s1);
   5.381 +    NetworkSimplex<Digraph> mcf5(gr, l2, u, c, v, w, 27);
   5.382 +    NetworkSimplex<Digraph> mcf6(gr, l2, u, c, s3);
   5.383 +
   5.384 +    checkMcf(mcf1, mcf1.run(pr), gr, l1, u, c, s1, true,  5240, "#G1");
   5.385 +    checkMcf(mcf2, mcf2.run(pr), gr, l1, u, c, s2, true,  7620, "#G2");
   5.386 +    checkMcf(mcf3, mcf3.run(pr), gr, l1, u, c, s3, true,     0, "#G3");
   5.387 +    checkMcf(mcf4, mcf4.run(pr), gr, l2, u, c, s1, true,  5970, "#G4");
   5.388 +    checkMcf(mcf5, mcf5.run(pr), gr, l2, u, c, s2, true,  8010, "#G5");
   5.389 +    checkMcf(mcf6, mcf6.run(pr), gr, l2, u, c, s3, false,    0, "#G6");
   5.390 +  }
   5.391 +
   5.392 +  // H. Test NetworkSimplex with CANDIDATE_LIST_PIVOT
   5.393 +  {
   5.394 +    NetworkSimplex<Digraph>::PivotRuleEnum pr =
   5.395 +      NetworkSimplex<Digraph>::CANDIDATE_LIST_PIVOT;
   5.396 +    NetworkSimplex<Digraph> mcf1(gr, u, c, s1);
   5.397 +    NetworkSimplex<Digraph> mcf2(gr, u, c, v, w, 27);
   5.398 +    NetworkSimplex<Digraph> mcf3(gr, u, c, s3);
   5.399 +    NetworkSimplex<Digraph> mcf4(gr, l2, u, c, s1);
   5.400 +    NetworkSimplex<Digraph> mcf5(gr, l2, u, c, v, w, 27);
   5.401 +    NetworkSimplex<Digraph> mcf6(gr, l2, u, c, s3);
   5.402 +
   5.403 +    checkMcf(mcf1, mcf1.run(pr), gr, l1, u, c, s1, true,  5240, "#H1");
   5.404 +    checkMcf(mcf2, mcf2.run(pr), gr, l1, u, c, s2, true,  7620, "#H2");
   5.405 +    checkMcf(mcf3, mcf3.run(pr), gr, l1, u, c, s3, true,     0, "#H3");
   5.406 +    checkMcf(mcf4, mcf4.run(pr), gr, l2, u, c, s1, true,  5970, "#H4");
   5.407 +    checkMcf(mcf5, mcf5.run(pr), gr, l2, u, c, s2, true,  8010, "#H5");
   5.408 +    checkMcf(mcf6, mcf6.run(pr), gr, l2, u, c, s3, false,    0, "#H6");
   5.409 +  }
   5.410 +
   5.411 +  // I. Test NetworkSimplex with ALTERING_LIST_PIVOT
   5.412 +  {
   5.413 +    NetworkSimplex<Digraph>::PivotRuleEnum pr =
   5.414 +      NetworkSimplex<Digraph>::ALTERING_LIST_PIVOT;
   5.415 +    NetworkSimplex<Digraph> mcf1(gr, u, c, s1);
   5.416 +    NetworkSimplex<Digraph> mcf2(gr, u, c, v, w, 27);
   5.417 +    NetworkSimplex<Digraph> mcf3(gr, u, c, s3);
   5.418 +    NetworkSimplex<Digraph> mcf4(gr, l2, u, c, s1);
   5.419 +    NetworkSimplex<Digraph> mcf5(gr, l2, u, c, v, w, 27);
   5.420 +    NetworkSimplex<Digraph> mcf6(gr, l2, u, c, s3);
   5.421 +
   5.422 +    checkMcf(mcf1, mcf1.run(pr), gr, l1, u, c, s1, true,  5240, "#I1");
   5.423 +    checkMcf(mcf2, mcf2.run(pr), gr, l1, u, c, s2, true,  7620, "#I2");
   5.424 +    checkMcf(mcf3, mcf3.run(pr), gr, l1, u, c, s3, true,     0, "#I3");
   5.425 +    checkMcf(mcf4, mcf4.run(pr), gr, l2, u, c, s1, true,  5970, "#I4");
   5.426 +    checkMcf(mcf5, mcf5.run(pr), gr, l2, u, c, s2, true,  8010, "#I5");
   5.427 +    checkMcf(mcf6, mcf6.run(pr), gr, l2, u, c, s3, false,    0, "#I6");
   5.428 +  }
   5.429 +
   5.430 +/*
   5.431 +  // J. Test MinCostFlow
   5.432 +  {
   5.433 +    MinCostFlow<Digraph> mcf1(gr, u, c, s1);
   5.434 +    MinCostFlow<Digraph> mcf2(gr, u, c, v, w, 27);
   5.435 +    MinCostFlow<Digraph> mcf3(gr, u, c, s3);
   5.436 +    MinCostFlow<Digraph> mcf4(gr, l2, u, c, s1);
   5.437 +    MinCostFlow<Digraph> mcf5(gr, l2, u, c, v, w, 27);
   5.438 +    MinCostFlow<Digraph> mcf6(gr, l2, u, c, s3);
   5.439 +
   5.440 +    checkMcf(mcf1, mcf1.run(), gr, l1, u, c, s1, true,  5240, "#J1");
   5.441 +    checkMcf(mcf2, mcf2.run(), gr, l1, u, c, s2, true,  7620, "#J2");
   5.442 +    checkMcf(mcf3, mcf3.run(), gr, l1, u, c, s3, true,     0, "#J3");
   5.443 +    checkMcf(mcf4, mcf4.run(), gr, l2, u, c, s1, true,  5970, "#J4");
   5.444 +    checkMcf(mcf5, mcf5.run(), gr, l2, u, c, s2, true,  8010, "#J5");
   5.445 +    checkMcf(mcf6, mcf6.run(), gr, l2, u, c, s3, false,    0, "#J6");
   5.446 +  }
   5.447 +*/
   5.448 +/*
   5.449 +  // K. Test MinCostMaxFlow
   5.450 +  {
   5.451 +    MinCostMaxFlow<Digraph> mcmf(gr, u, c, v, w);
   5.452 +    mcmf.run();
   5.453 +    checkMcf(mcmf, true, gr, l1, u, c, s3, true, 7620, "#K1");
   5.454 +  }
   5.455 +*/
   5.456 +
   5.457 +  return 0;
   5.458 +}