lemon/capacity_scaling.h
author ladanyi
Sun, 02 Mar 2008 22:55:27 +0000
changeset 2592 f1fb0c31f952
parent 2588 4d3bc1d04c1d
child 2620 8f41a3129746
permissions -rw-r--r--
Revert to long long int since currently I don't know a better solution.
     1 /* -*- C++ -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library
     4  *
     5  * Copyright (C) 2003-2008
     6  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7  * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8  *
     9  * Permission to use, modify and distribute this software is granted
    10  * provided that this copyright notice appears in all copies. For
    11  * precise terms see the accompanying LICENSE file.
    12  *
    13  * This software is provided "AS IS" with no warranty of any kind,
    14  * express or implied, and with no claim as to its suitability for any
    15  * purpose.
    16  *
    17  */
    18 
    19 #ifndef LEMON_CAPACITY_SCALING_H
    20 #define LEMON_CAPACITY_SCALING_H
    21 
    22 /// \ingroup min_cost_flow
    23 ///
    24 /// \file
    25 /// \brief Capacity scaling algorithm for finding a minimum cost flow.
    26 
    27 #include <vector>
    28 #include <lemon/bin_heap.h>
    29 
    30 namespace lemon {
    31 
    32   /// \addtogroup min_cost_flow
    33   /// @{
    34 
    35   /// \brief Implementation of the capacity scaling algorithm for
    36   /// finding a minimum cost flow.
    37   ///
    38   /// \ref CapacityScaling implements the capacity scaling version
    39   /// of the successive shortest path algorithm for finding a minimum
    40   /// cost flow.
    41   ///
    42   /// \tparam Graph The directed graph type the algorithm runs on.
    43   /// \tparam LowerMap The type of the lower bound map.
    44   /// \tparam CapacityMap The type of the capacity (upper bound) map.
    45   /// \tparam CostMap The type of the cost (length) map.
    46   /// \tparam SupplyMap The type of the supply map.
    47   ///
    48   /// \warning
    49   /// - Edge capacities and costs should be \e non-negative \e integers.
    50   /// - Supply values should be \e signed \e integers.
    51   /// - The value types of the maps should be convertible to each other.
    52   /// - \c CostMap::Value must be signed type.
    53   ///
    54   /// \author Peter Kovacs
    55 
    56   template < typename Graph,
    57              typename LowerMap = typename Graph::template EdgeMap<int>,
    58              typename CapacityMap = typename Graph::template EdgeMap<int>,
    59              typename CostMap = typename Graph::template EdgeMap<int>,
    60              typename SupplyMap = typename Graph::template NodeMap<int> >
    61   class CapacityScaling
    62   {
    63     GRAPH_TYPEDEFS(typename Graph);
    64 
    65     typedef typename CapacityMap::Value Capacity;
    66     typedef typename CostMap::Value Cost;
    67     typedef typename SupplyMap::Value Supply;
    68     typedef typename Graph::template EdgeMap<Capacity> CapacityEdgeMap;
    69     typedef typename Graph::template NodeMap<Supply> SupplyNodeMap;
    70     typedef typename Graph::template NodeMap<Edge> PredMap;
    71 
    72   public:
    73 
    74     /// The type of the flow map.
    75     typedef typename Graph::template EdgeMap<Capacity> FlowMap;
    76     /// The type of the potential map.
    77     typedef typename Graph::template NodeMap<Cost> PotentialMap;
    78 
    79   private:
    80 
    81     /// \brief Special implementation of the \ref Dijkstra algorithm
    82     /// for finding shortest paths in the residual network.
    83     ///
    84     /// \ref ResidualDijkstra is a special implementation of the
    85     /// \ref Dijkstra algorithm for finding shortest paths in the
    86     /// residual network of the graph with respect to the reduced edge
    87     /// costs and modifying the node potentials according to the
    88     /// distance of the nodes.
    89     class ResidualDijkstra
    90     {
    91       typedef typename Graph::template NodeMap<int> HeapCrossRef;
    92       typedef BinHeap<Cost, HeapCrossRef> Heap;
    93 
    94     private:
    95 
    96       // The directed graph the algorithm runs on
    97       const Graph &_graph;
    98 
    99       // The main maps
   100       const FlowMap &_flow;
   101       const CapacityEdgeMap &_res_cap;
   102       const CostMap &_cost;
   103       const SupplyNodeMap &_excess;
   104       PotentialMap &_potential;
   105 
   106       // The distance map
   107       PotentialMap _dist;
   108       // The pred edge map
   109       PredMap &_pred;
   110       // The processed (i.e. permanently labeled) nodes
   111       std::vector<Node> _proc_nodes;
   112 
   113     public:
   114 
   115       /// Constructor.
   116       ResidualDijkstra( const Graph &graph,
   117                         const FlowMap &flow,
   118                         const CapacityEdgeMap &res_cap,
   119                         const CostMap &cost,
   120                         const SupplyMap &excess,
   121                         PotentialMap &potential,
   122                         PredMap &pred ) :
   123         _graph(graph), _flow(flow), _res_cap(res_cap), _cost(cost),
   124         _excess(excess), _potential(potential), _dist(graph),
   125         _pred(pred)
   126       {}
   127 
   128       /// Runs the algorithm from the given source node.
   129       Node run(Node s, Capacity delta = 1) {
   130         HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);
   131         Heap heap(heap_cross_ref);
   132         heap.push(s, 0);
   133         _pred[s] = INVALID;
   134         _proc_nodes.clear();
   135 
   136         // Processing nodes
   137         while (!heap.empty() && _excess[heap.top()] > -delta) {
   138           Node u = heap.top(), v;
   139           Cost d = heap.prio() + _potential[u], nd;
   140           _dist[u] = heap.prio();
   141           heap.pop();
   142           _proc_nodes.push_back(u);
   143 
   144           // Traversing outgoing edges
   145           for (OutEdgeIt e(_graph, u); e != INVALID; ++e) {
   146             if (_res_cap[e] >= delta) {
   147               v = _graph.target(e);
   148               switch(heap.state(v)) {
   149               case Heap::PRE_HEAP:
   150                 heap.push(v, d + _cost[e] - _potential[v]);
   151                 _pred[v] = e;
   152                 break;
   153               case Heap::IN_HEAP:
   154                 nd = d + _cost[e] - _potential[v];
   155                 if (nd < heap[v]) {
   156                   heap.decrease(v, nd);
   157                   _pred[v] = e;
   158                 }
   159                 break;
   160               case Heap::POST_HEAP:
   161                 break;
   162               }
   163             }
   164           }
   165 
   166           // Traversing incoming edges
   167           for (InEdgeIt e(_graph, u); e != INVALID; ++e) {
   168             if (_flow[e] >= delta) {
   169               v = _graph.source(e);
   170               switch(heap.state(v)) {
   171               case Heap::PRE_HEAP:
   172                 heap.push(v, d - _cost[e] - _potential[v]);
   173                 _pred[v] = e;
   174                 break;
   175               case Heap::IN_HEAP:
   176                 nd = d - _cost[e] - _potential[v];
   177                 if (nd < heap[v]) {
   178                   heap.decrease(v, nd);
   179                   _pred[v] = e;
   180                 }
   181                 break;
   182               case Heap::POST_HEAP:
   183                 break;
   184               }
   185             }
   186           }
   187         }
   188         if (heap.empty()) return INVALID;
   189 
   190         // Updating potentials of processed nodes
   191         Node t = heap.top();
   192         Cost t_dist = heap.prio();
   193         for (int i = 0; i < int(_proc_nodes.size()); ++i)
   194           _potential[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist;
   195 
   196         return t;
   197       }
   198 
   199     }; //class ResidualDijkstra
   200 
   201   private:
   202 
   203     // The directed graph the algorithm runs on
   204     const Graph &_graph;
   205     // The original lower bound map
   206     const LowerMap *_lower;
   207     // The modified capacity map
   208     CapacityEdgeMap _capacity;
   209     // The original cost map
   210     const CostMap &_cost;
   211     // The modified supply map
   212     SupplyNodeMap _supply;
   213     bool _valid_supply;
   214 
   215     // Edge map of the current flow
   216     FlowMap *_flow;
   217     bool _local_flow;
   218     // Node map of the current potentials
   219     PotentialMap *_potential;
   220     bool _local_potential;
   221 
   222     // The residual capacity map
   223     CapacityEdgeMap _res_cap;
   224     // The excess map
   225     SupplyNodeMap _excess;
   226     // The excess nodes (i.e. nodes with positive excess)
   227     std::vector<Node> _excess_nodes;
   228     // The deficit nodes (i.e. nodes with negative excess)
   229     std::vector<Node> _deficit_nodes;
   230 
   231     // The delta parameter used for capacity scaling
   232     Capacity _delta;
   233     // The maximum number of phases
   234     int _phase_num;
   235 
   236     // The pred edge map
   237     PredMap _pred;
   238     // Implementation of the Dijkstra algorithm for finding augmenting
   239     // shortest paths in the residual network
   240     ResidualDijkstra *_dijkstra;
   241 
   242   public:
   243 
   244     /// \brief General constructor (with lower bounds).
   245     ///
   246     /// General constructor (with lower bounds).
   247     ///
   248     /// \param graph The directed graph the algorithm runs on.
   249     /// \param lower The lower bounds of the edges.
   250     /// \param capacity The capacities (upper bounds) of the edges.
   251     /// \param cost The cost (length) values of the edges.
   252     /// \param supply The supply values of the nodes (signed).
   253     CapacityScaling( const Graph &graph,
   254                      const LowerMap &lower,
   255                      const CapacityMap &capacity,
   256                      const CostMap &cost,
   257                      const SupplyMap &supply ) :
   258       _graph(graph), _lower(&lower), _capacity(graph), _cost(cost),
   259       _supply(graph), _flow(0), _local_flow(false),
   260       _potential(0), _local_potential(false),
   261       _res_cap(graph), _excess(graph), _pred(graph)
   262     {
   263       // Removing non-zero lower bounds
   264       _capacity = subMap(capacity, lower);
   265       _res_cap = _capacity;
   266       Supply sum = 0;
   267       for (NodeIt n(_graph); n != INVALID; ++n) {
   268         Supply s = supply[n];
   269         for (InEdgeIt e(_graph, n); e != INVALID; ++e)
   270           s += lower[e];
   271         for (OutEdgeIt e(_graph, n); e != INVALID; ++e)
   272           s -= lower[e];
   273         _supply[n] = s;
   274         sum += s;
   275       }
   276       _valid_supply = sum == 0;
   277     }
   278 
   279     /// \brief General constructor (without lower bounds).
   280     ///
   281     /// General constructor (without lower bounds).
   282     ///
   283     /// \param graph The directed graph the algorithm runs on.
   284     /// \param capacity The capacities (upper bounds) of the edges.
   285     /// \param cost The cost (length) values of the edges.
   286     /// \param supply The supply values of the nodes (signed).
   287     CapacityScaling( const Graph &graph,
   288                      const CapacityMap &capacity,
   289                      const CostMap &cost,
   290                      const SupplyMap &supply ) :
   291       _graph(graph), _lower(NULL), _capacity(capacity), _cost(cost),
   292       _supply(supply), _flow(0), _local_flow(false),
   293       _potential(0), _local_potential(false),
   294       _res_cap(capacity), _excess(graph), _pred(graph)
   295     {
   296       // Checking the sum of supply values
   297       Supply sum = 0;
   298       for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
   299       _valid_supply = sum == 0;
   300     }
   301 
   302     /// \brief Simple constructor (with lower bounds).
   303     ///
   304     /// Simple constructor (with lower bounds).
   305     ///
   306     /// \param graph The directed graph the algorithm runs on.
   307     /// \param lower The lower bounds of the edges.
   308     /// \param capacity The capacities (upper bounds) of the edges.
   309     /// \param cost The cost (length) values of the edges.
   310     /// \param s The source node.
   311     /// \param t The target node.
   312     /// \param flow_value The required amount of flow from node \c s
   313     /// to node \c t (i.e. the supply of \c s and the demand of \c t).
   314     CapacityScaling( const Graph &graph,
   315                      const LowerMap &lower,
   316                      const CapacityMap &capacity,
   317                      const CostMap &cost,
   318                      Node s, Node t,
   319                      Supply flow_value ) :
   320       _graph(graph), _lower(&lower), _capacity(graph), _cost(cost),
   321       _supply(graph), _flow(0), _local_flow(false),
   322       _potential(0), _local_potential(false),
   323       _res_cap(graph), _excess(graph), _pred(graph)
   324     {
   325       // Removing non-zero lower bounds
   326       _capacity = subMap(capacity, lower);
   327       _res_cap = _capacity;
   328       for (NodeIt n(_graph); n != INVALID; ++n) {
   329         Supply sum = 0;
   330         if (n == s) sum =  flow_value;
   331         if (n == t) sum = -flow_value;
   332         for (InEdgeIt e(_graph, n); e != INVALID; ++e)
   333           sum += lower[e];
   334         for (OutEdgeIt e(_graph, n); e != INVALID; ++e)
   335           sum -= lower[e];
   336         _supply[n] = sum;
   337       }
   338       _valid_supply = true;
   339     }
   340 
   341     /// \brief Simple constructor (without lower bounds).
   342     ///
   343     /// Simple constructor (without lower bounds).
   344     ///
   345     /// \param graph The directed graph the algorithm runs on.
   346     /// \param capacity The capacities (upper bounds) of the edges.
   347     /// \param cost The cost (length) values of the edges.
   348     /// \param s The source node.
   349     /// \param t The target node.
   350     /// \param flow_value The required amount of flow from node \c s
   351     /// to node \c t (i.e. the supply of \c s and the demand of \c t).
   352     CapacityScaling( const Graph &graph,
   353                      const CapacityMap &capacity,
   354                      const CostMap &cost,
   355                      Node s, Node t,
   356                      Supply flow_value ) :
   357       _graph(graph), _lower(NULL), _capacity(capacity), _cost(cost),
   358       _supply(graph, 0), _flow(0), _local_flow(false),
   359       _potential(0), _local_potential(false),
   360       _res_cap(capacity), _excess(graph), _pred(graph)
   361     {
   362       _supply[s] =  flow_value;
   363       _supply[t] = -flow_value;
   364       _valid_supply = true;
   365     }
   366 
   367     /// Destructor.
   368     ~CapacityScaling() {
   369       if (_local_flow) delete _flow;
   370       if (_local_potential) delete _potential;
   371       delete _dijkstra;
   372     }
   373 
   374     /// \brief Sets the flow map.
   375     ///
   376     /// Sets the flow map.
   377     ///
   378     /// \return \c (*this)
   379     CapacityScaling& flowMap(FlowMap &map) {
   380       if (_local_flow) {
   381         delete _flow;
   382         _local_flow = false;
   383       }
   384       _flow = &map;
   385       return *this;
   386     }
   387 
   388     /// \brief Sets the potential map.
   389     ///
   390     /// Sets the potential map.
   391     ///
   392     /// \return \c (*this)
   393     CapacityScaling& potentialMap(PotentialMap &map) {
   394       if (_local_potential) {
   395         delete _potential;
   396         _local_potential = false;
   397       }
   398       _potential = &map;
   399       return *this;
   400     }
   401 
   402     /// \name Execution control
   403     /// The only way to execute the algorithm is to call the run()
   404     /// function.
   405 
   406     /// @{
   407 
   408     /// \brief Runs the algorithm.
   409     ///
   410     /// Runs the algorithm.
   411     ///
   412     /// \param scaling Enable or disable capacity scaling.
   413     /// If the maximum edge capacity and/or the amount of total supply
   414     /// is rather small, the algorithm could be slightly faster without
   415     /// scaling.
   416     ///
   417     /// \return \c true if a feasible flow can be found.
   418     bool run(bool scaling = true) {
   419       return init(scaling) && start();
   420     }
   421 
   422     /// @}
   423 
   424     /// \name Query Functions
   425     /// The result of the algorithm can be obtained using these
   426     /// functions.
   427     /// \n run() must be called before using them.
   428 
   429     /// @{
   430 
   431     /// \brief Returns a const reference to the edge map storing the
   432     /// found flow.
   433     ///
   434     /// Returns a const reference to the edge map storing the found flow.
   435     ///
   436     /// \pre \ref run() must be called before using this function.
   437     const FlowMap& flowMap() const {
   438       return *_flow;
   439     }
   440 
   441     /// \brief Returns a const reference to the node map storing the
   442     /// found potentials (the dual solution).
   443     ///
   444     /// Returns a const reference to the node map storing the found
   445     /// potentials (the dual solution).
   446     ///
   447     /// \pre \ref run() must be called before using this function.
   448     const PotentialMap& potentialMap() const {
   449       return *_potential;
   450     }
   451 
   452     /// \brief Returns the flow on the given edge.
   453     ///
   454     /// Returns the flow on the given edge.
   455     ///
   456     /// \pre \ref run() must be called before using this function.
   457     Capacity flow(const Edge& edge) const {
   458       return (*_flow)[edge];
   459     }
   460 
   461     /// \brief Returns the potential of the given node.
   462     ///
   463     /// Returns the potential of the given node.
   464     ///
   465     /// \pre \ref run() must be called before using this function.
   466     Cost potential(const Node& node) const {
   467       return (*_potential)[node];
   468     }
   469 
   470     /// \brief Returns the total cost of the found flow.
   471     ///
   472     /// Returns the total cost of the found flow. The complexity of the
   473     /// function is \f$ O(e) \f$.
   474     ///
   475     /// \pre \ref run() must be called before using this function.
   476     Cost totalCost() const {
   477       Cost c = 0;
   478       for (EdgeIt e(_graph); e != INVALID; ++e)
   479         c += (*_flow)[e] * _cost[e];
   480       return c;
   481     }
   482 
   483     /// @}
   484 
   485   private:
   486 
   487     /// Initializes the algorithm.
   488     bool init(bool scaling) {
   489       if (!_valid_supply) return false;
   490 
   491       // Initializing maps
   492       if (!_flow) {
   493         _flow = new FlowMap(_graph);
   494         _local_flow = true;
   495       }
   496       if (!_potential) {
   497         _potential = new PotentialMap(_graph);
   498         _local_potential = true;
   499       }
   500       for (EdgeIt e(_graph); e != INVALID; ++e) (*_flow)[e] = 0;
   501       for (NodeIt n(_graph); n != INVALID; ++n) (*_potential)[n] = 0;
   502       _excess = _supply;
   503 
   504       _dijkstra = new ResidualDijkstra( _graph, *_flow, _res_cap, _cost,
   505                                         _excess, *_potential, _pred );
   506 
   507       // Initializing delta value
   508       if (scaling) {
   509         // With scaling
   510         Supply max_sup = 0, max_dem = 0;
   511         for (NodeIt n(_graph); n != INVALID; ++n) {
   512           if ( _supply[n] > max_sup) max_sup =  _supply[n];
   513           if (-_supply[n] > max_dem) max_dem = -_supply[n];
   514         }
   515         Capacity max_cap = 0;
   516         for (EdgeIt e(_graph); e != INVALID; ++e) {
   517           if (_capacity[e] > max_cap) max_cap = _capacity[e];
   518         }
   519         max_sup = std::min(std::min(max_sup, max_dem), max_cap);
   520         _phase_num = 0;
   521         for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2)
   522           ++_phase_num;
   523       } else {
   524         // Without scaling
   525         _delta = 1;
   526       }
   527 
   528       return true;
   529     }
   530 
   531     bool start() {
   532       if (_delta > 1)
   533         return startWithScaling();
   534       else
   535         return startWithoutScaling();
   536     }
   537 
   538     /// Executes the capacity scaling algorithm.
   539     bool startWithScaling() {
   540       // Processing capacity scaling phases
   541       Node s, t;
   542       int phase_cnt = 0;
   543       int factor = 4;
   544       while (true) {
   545         // Saturating all edges not satisfying the optimality condition
   546         for (EdgeIt e(_graph); e != INVALID; ++e) {
   547           Node u = _graph.source(e), v = _graph.target(e);
   548           Cost c = _cost[e] + (*_potential)[u] - (*_potential)[v];
   549           if (c < 0 && _res_cap[e] >= _delta) {
   550             _excess[u] -= _res_cap[e];
   551             _excess[v] += _res_cap[e];
   552             (*_flow)[e] = _capacity[e];
   553             _res_cap[e] = 0;
   554           }
   555           else if (c > 0 && (*_flow)[e] >= _delta) {
   556             _excess[u] += (*_flow)[e];
   557             _excess[v] -= (*_flow)[e];
   558             (*_flow)[e] = 0;
   559             _res_cap[e] = _capacity[e];
   560           }
   561         }
   562 
   563         // Finding excess nodes and deficit nodes
   564         _excess_nodes.clear();
   565         _deficit_nodes.clear();
   566         for (NodeIt n(_graph); n != INVALID; ++n) {
   567           if (_excess[n] >=  _delta) _excess_nodes.push_back(n);
   568           if (_excess[n] <= -_delta) _deficit_nodes.push_back(n);
   569         }
   570         int next_node = 0;
   571 
   572         // Finding augmenting shortest paths
   573         while (next_node < int(_excess_nodes.size())) {
   574           // Checking deficit nodes
   575           if (_delta > 1) {
   576             bool delta_deficit = false;
   577             for (int i = 0; i < int(_deficit_nodes.size()); ++i) {
   578               if (_excess[_deficit_nodes[i]] <= -_delta) {
   579                 delta_deficit = true;
   580                 break;
   581               }
   582             }
   583             if (!delta_deficit) break;
   584           }
   585 
   586           // Running Dijkstra
   587           s = _excess_nodes[next_node];
   588           if ((t = _dijkstra->run(s, _delta)) == INVALID) {
   589             if (_delta > 1) {
   590               ++next_node;
   591               continue;
   592             }
   593             return false;
   594           }
   595 
   596           // Augmenting along a shortest path from s to t.
   597           Capacity d = std::min(_excess[s], -_excess[t]);
   598           Node u = t;
   599           Edge e;
   600           if (d > _delta) {
   601             while ((e = _pred[u]) != INVALID) {
   602               Capacity rc;
   603               if (u == _graph.target(e)) {
   604                 rc = _res_cap[e];
   605                 u = _graph.source(e);
   606               } else {
   607                 rc = (*_flow)[e];
   608                 u = _graph.target(e);
   609               }
   610               if (rc < d) d = rc;
   611             }
   612           }
   613           u = t;
   614           while ((e = _pred[u]) != INVALID) {
   615             if (u == _graph.target(e)) {
   616               (*_flow)[e] += d;
   617               _res_cap[e] -= d;
   618               u = _graph.source(e);
   619             } else {
   620               (*_flow)[e] -= d;
   621               _res_cap[e] += d;
   622               u = _graph.target(e);
   623             }
   624           }
   625           _excess[s] -= d;
   626           _excess[t] += d;
   627 
   628           if (_excess[s] < _delta) ++next_node;
   629         }
   630 
   631         if (_delta == 1) break;
   632         if (++phase_cnt > _phase_num / 4) factor = 2;
   633         _delta = _delta <= factor ? 1 : _delta / factor;
   634       }
   635 
   636       // Handling non-zero lower bounds
   637       if (_lower) {
   638         for (EdgeIt e(_graph); e != INVALID; ++e)
   639           (*_flow)[e] += (*_lower)[e];
   640       }
   641       return true;
   642     }
   643 
   644     /// Executes the successive shortest path algorithm.
   645     bool startWithoutScaling() {
   646       // Finding excess nodes
   647       for (NodeIt n(_graph); n != INVALID; ++n)
   648         if (_excess[n] > 0) _excess_nodes.push_back(n);
   649       if (_excess_nodes.size() == 0) return true;
   650       int next_node = 0;
   651 
   652       // Finding shortest paths
   653       Node s, t;
   654       while ( _excess[_excess_nodes[next_node]] > 0 ||
   655               ++next_node < int(_excess_nodes.size()) )
   656       {
   657         // Running Dijkstra
   658         s = _excess_nodes[next_node];
   659         if ((t = _dijkstra->run(s)) == INVALID) return false;
   660 
   661         // Augmenting along a shortest path from s to t
   662         Capacity d = std::min(_excess[s], -_excess[t]);
   663         Node u = t;
   664         Edge e;
   665         if (d > 1) {
   666           while ((e = _pred[u]) != INVALID) {
   667             Capacity rc;
   668             if (u == _graph.target(e)) {
   669               rc = _res_cap[e];
   670               u = _graph.source(e);
   671             } else {
   672               rc = (*_flow)[e];
   673               u = _graph.target(e);
   674             }
   675             if (rc < d) d = rc;
   676           }
   677         }
   678         u = t;
   679         while ((e = _pred[u]) != INVALID) {
   680           if (u == _graph.target(e)) {
   681             (*_flow)[e] += d;
   682             _res_cap[e] -= d;
   683             u = _graph.source(e);
   684           } else {
   685             (*_flow)[e] -= d;
   686             _res_cap[e] += d;
   687             u = _graph.target(e);
   688           }
   689         }
   690         _excess[s] -= d;
   691         _excess[t] += d;
   692       }
   693 
   694       // Handling non-zero lower bounds
   695       if (_lower) {
   696         for (EdgeIt e(_graph); e != INVALID; ++e)
   697           (*_flow)[e] += (*_lower)[e];
   698       }
   699       return true;
   700     }
   701 
   702   }; //class CapacityScaling
   703 
   704   ///@}
   705 
   706 } //namespace lemon
   707 
   708 #endif //LEMON_CAPACITY_SCALING_H