# HG changeset patch
# User deba
# Date 1195333091 0
# Node ID 57143c09dc20a4b6d65d0c58f709331be66dfb72
# Parent  26983135fd6d2842904bc7b2a486ab8ae57bd5b7
Redesign the maximum flow algorithms

Redesigned interface
Preflow changed to use elevator
Edmonds-Karp does not use the ResGraphAdaptor
Goldberg-Tarjan algorithm (Preflow with Dynamic Trees)
Dinitz-Sleator-Tarjan (Blocking flow with Dynamic Tree)

diff -r 26983135fd6d -r 57143c09dc20 benchmark/hcube.cc
--- a/benchmark/hcube.cc	Wed Nov 14 17:53:08 2007 +0000
+++ b/benchmark/hcube.cc	Sat Nov 17 20:58:11 2007 +0000
@@ -112,9 +112,9 @@
   {
    Graph::EdgeMap<int> flow(G);
    
-    Preflow<Graph,int> MF(G,nodes[0],nodes[1<<dim-1],map,flow);
+    Preflow<Graph> MF(G,map,nodes[0],nodes[1<<dim-1]);
     for(int i=0;i<10;i++)
-      MF.run(MF.NO_FLOW);
+      MF.run();
   }
   PrintTime("PREFLOW",T);
 
diff -r 26983135fd6d -r 57143c09dc20 demo/disjoint_paths_demo.cc
--- a/demo/disjoint_paths_demo.cc	Wed Nov 14 17:53:08 2007 +0000
+++ b/demo/disjoint_paths_demo.cc	Sat Nov 17 20:58:11 2007 +0000
@@ -63,7 +63,7 @@
 
   Graph::EdgeMap<int> flow(graph);
 
-  Preflow<Graph, int, Capacity> preflow(graph, source, target, capacity, flow); 
+  Preflow<Graph, Capacity> preflow(graph, capacity, source, target); 
   
   preflow.run();
 
@@ -72,7 +72,7 @@
   graphToEps(graph, "edge_disjoint_paths.eps").
     title("edge disjoint paths").scaleToA4().
     copyright("(C) 2003-2007 LEMON Project").drawArrows().
-    edgeColors(composeMap(functorMap(color), flow)).
+    edgeColors(composeMap(functorMap(color), preflow.flowMap())).
     coords(coords).run();
 
 
@@ -87,9 +87,9 @@
 
   SGraph::EdgeMap<int> sflow(sgraph);
 
-  Preflow<SGraph, int, SCapacity> spreflow(sgraph, SGraph::outNode(source), 
-                                           SGraph::inNode(target), 
-                                           scapacity, sflow);
+  Preflow<SGraph, SCapacity> spreflow(sgraph, scapacity, 
+				      SGraph::outNode(source), 
+				      SGraph::inNode(target));
 
   spreflow.run();
 
@@ -99,7 +99,7 @@
   graphToEps(sgraph, "node_disjoint_paths.eps").
     title("node disjoint paths").scaleToA4().
     copyright("(C) 2003-2007 LEMON Project").drawArrows().
-    edgeColors(composeMap(functorMap(color), sflow)).
+    edgeColors(composeMap(functorMap(color), spreflow.flowMap())).
     coords(SGraph::combinedNodeMap(coords,
 				   shiftMap(coords,
 					    dim2::Point<double>(5, 0)))).
diff -r 26983135fd6d -r 57143c09dc20 demo/sub_graph_adaptor_demo.cc
--- a/demo/sub_graph_adaptor_demo.cc	Wed Nov 14 17:53:08 2007 +0000
+++ b/demo/sub_graph_adaptor_demo.cc	Sat Nov 17 20:58:11 2007 +0000
@@ -109,10 +109,8 @@
   SubGW gw(g, tight_edge_filter);
 
   ConstMap<Edge, int> const_1_map(1);
-  Graph::EdgeMap<int> flow(g, 0);
   // Max flow between s and t in the graph of tight edges.
-  Preflow<SubGW, int, ConstMap<Edge, int>, Graph::EdgeMap<int> > 
-    preflow(gw, s, t, const_1_map, flow);
+  Preflow<SubGW, ConstMap<Edge, int> > preflow(gw, const_1_map, s, t);
   preflow.run();
 
   cout << "maximum number of edge-disjoint shortest paths: " 
@@ -120,7 +118,7 @@
   cout << "edges of the maximum number of edge-disjoint shortest s-t paths: " 
        << endl;
   for(EdgeIt e(g); e!=INVALID; ++e) 
-    if (flow[e])
+    if (preflow.flow(e) != 0)
       cout << " " << g.id(e) << ", "
 	   << g.id(g.source(e)) << "--" 
 	   << length[e] << "->" << g.id(g.target(e)) << endl;
diff -r 26983135fd6d -r 57143c09dc20 doc/groups.dox
--- a/doc/groups.dox	Wed Nov 14 17:53:08 2007 +0000
+++ b/doc/groups.dox	Sat Nov 17 20:58:11 2007 +0000
@@ -223,8 +223,27 @@
 This group describes the algorithms for finding maximum flows and
 feasible circulations.
 
-\image html flow.png
-\image latex flow.eps "Graph flow" width=\textwidth
+The maximum flow problem is to find a flow between a single-source and
+single-target that is maximum. Formally, there is \f$G=(V,A)\f$
+directed graph, an \f$c_a:A\rightarrow\mathbf{R}^+_0\f$ capacity
+function and given \f$s, t \in V\f$ source and target node. The
+maximum flow is the solution of the next optimization problem:
+
+\f[ 0 \le f_a \le c_a \f]
+\f[ \sum_{v\in\delta^{-}(u)}f_{vu}=\sum_{v\in\delta^{+}(u)}f_{uv} \quad u \in V \setminus \{s,t\}\f]
+\f[ \max \sum_{v\in\delta^{+}(s)}f_{uv} - \sum_{v\in\delta^{-}(s)}f_{vu}\f]
+
+The lemon contains several algorithms for solve maximum flow problems:
+- \ref lemon::EdmondsKarp "Edmonds-Karp" 
+- \ref lemon::Preflow "Goldberg's Preflow algorithm"
+- \ref lemon::DinitzSleatorTarjan "Dinitz's blocking flow algorithm with dynamic tree"
+- \ref lemon::GoldbergTarjan "Preflow algorithm with dynamic trees"
+
+In most cases the \ref lemon::Preflow "preflow" algorithm provides the
+fastest method to compute the maximum flow. All impelementations
+provides functions for query the minimum cut, which is the dual linear
+programming probelm of the maximum flow.
+
 */
 
 /**
diff -r 26983135fd6d -r 57143c09dc20 doc/images/flow.eps
--- a/doc/images/flow.eps	Wed Nov 14 17:53:08 2007 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,847 +0,0 @@
-%!PS-Adobe-3.0 EPSF-3.0
-%%BoundingBox: 13 59 827 531
-%%HiResBoundingBox: 13.1911 59.6331 826.078 530.254
-%%Creator: Karbon14 EPS Exportfilter 0.5
-%%CreationDate: (04/18/06 08:36:16)
-%%For: (Balazs Dezso) ()
-%%Title: ()
-
-/N {newpath} def
-/C {closepath} def
-/m {moveto} def
-/c {curveto} def
-/l {lineto} def
-/s {stroke} def
-/f {fill} def
-/w {setlinewidth} def
-/d {setdash} def
-/r {setrgbcolor} def
-/S {gsave} def
-/R {grestore} def
-
-N
-694.242 342.699 m
-748.648 255.656 753.472 247.937 799.086 174.957 c
-[] 0 d 0 0 0.5 r 3.58125 w s
-
-N
-807.312 161.797 m
-794.531 172.109 l
-803.64 177.805 l
-807.312 161.797 l
-C
-S 0 0 0.5 r f R
-
-N
-678.324 199.43 m
-736.031 179.062 744.519 176.066 787.746 160.812 c
-[] 0 d 0.6 0 0.2 r 3.58125 w s
-
-N
-802.382 155.644 m
-785.957 155.746 l
-789.531 165.879 l
-802.382 155.644 l
-C
-S 0.6 0 0.2 r f R
-
-N
-654.445 135.75 m
-723.707 142.676 732.703 143.578 786.316 148.937 c
-[] 0 d 0.6 0 0.2 r 3.58125 w s
-
-N
-801.761 150.484 m
-786.851 143.594 l
-785.781 154.285 l
-801.761 150.484 l
-C
-S 0.6 0 0.2 r f R
-
-N
-606.687 72.0741 m
-700.828 108.281 709.32 111.547 788.007 141.812 c
-[] 0 d 0.8 0 0.1 r 3.58125 w s
-
-N
-802.496 147.387 m
-789.937 136.797 l
-786.082 146.828 l
-802.496 147.387 l
-C
-S 0.8 0 0.1 r f R
-
-N
-694.242 342.699 m
-643.015 305.133 635.746 299.805 596.992 271.383 c
-[] 0 d 1 0 0 r 3.58125 w s
-
-N
-584.476 262.207 m
-593.816 275.715 l
-600.168 267.051 l
-584.476 262.207 l
-C
-S 1 0 0 r f R
-
-N
-582.808 358.621 m
-627.601 352.223 636.382 350.965 667.058 346.586 c
-[] 0 d 0.4 0 0.3 r 3.58125 w s
-
-N
-682.425 344.39 m
-666.3 341.266 l
-667.82 351.902 l
-682.425 344.39 l
-C
-S 0.4 0 0.3 r f R
-
-N
-511.171 517.816 m
-595.562 437.09 602.148 430.793 674.398 361.683 c
-[] 0 d 0.6 1 0.2 r 3.58125 w s
-
-N
-685.617 350.953 m
-670.687 357.801 l
-678.113 365.566 l
-685.617 350.953 l
-C
-S 0.6 1 0.2 r f R
-
-N
-495.25 143.711 m
-577.074 168.613 585.757 171.258 652.054 191.433 c
-[] 0 d 1 0 0 r 3.58125 w s
-
-N
-666.902 195.953 m
-653.617 186.293 l
-650.488 196.574 l
-666.902 195.953 l
-C
-S 1 0 0 r f R
-
-N
-574.847 255.148 m
-616.957 232.473 624.793 228.254 654.144 212.449 c
-[] 0 d 0.6 0 0.2 r 3.58125 w s
-
-N
-667.812 205.09 m
-651.597 207.719 l
-656.691 217.18 l
-667.812 205.09 l
-C
-S 0.6 0 0.2 r f R
-
-N
-495.25 143.711 m
-564.468 140.25 573.496 139.797 627.019 137.125 c
-[] 0 d 0.6 0 0.2 r 3.58125 w s
-
-N
-642.519 136.348 m
-626.75 131.758 l
-627.285 142.488 l
-642.519 136.348 l
-C
-S 0.6 0 0.2 r f R
-
-N
-391.773 72.0741 m
-489.199 72.0741 498.293 72.0741 579.226 72.0741 c
-[] 0 d 0.8 0 0.1 r 3.58125 w s
-
-N
-594.746 72.0741 m
-579.226 66.6991 l
-579.226 77.4451 l
-594.746 72.0741 l
-C
-S 0.8 0 0.1 r f R
-
-N
-495.25 143.711 m
-528.867 190.777 534.089 198.086 558.886 232.801 c
-[] 0 d 1 0 0 r 3.58125 w s
-
-N
-567.91 245.433 m
-563.257 229.68 l
-554.515 235.926 l
-567.91 245.433 l
-C
-S 1 0 0 r f R
-
-N
-431.574 271.062 m
-458.687 216.84 462.711 208.793 482.968 168.273 c
-[] 0 d 1 0 0 r 3.58125 w s
-
-N
-489.91 154.39 m
-478.164 165.871 l
-487.777 170.676 l
-489.91 154.39 l
-C
-S 1 0 0 r f R
-
-N
-391.773 72.0741 m
-434.636 101.75 441.992 106.84 472.671 128.082 c
-[] 0 d 0.6 0 0.2 r 3.58125 w s
-
-N
-485.433 136.914 m
-475.73 123.664 l
-469.613 132.496 l
-485.433 136.914 l
-C
-S 0.6 0 0.2 r f R
-
-N
-574.847 255.148 m
-577.964 295.68 578.64 304.457 580.703 331.238 c
-[] 0 d 1 0 0 r 3.58125 w s
-
-N
-581.89 346.715 m
-586.058 330.828 l
-575.343 331.652 l
-581.89 346.715 l
-C
-S 1 0 0 r f R
-
-N
-431.574 271.062 m
-492.73 264.269 501.675 263.273 547.554 258.18 c
-[] 0 d 0.6 0 0.2 r 3.58125 w s
-
-N
-562.98 256.465 m
-546.961 252.84 l
-548.148 263.516 l
-562.98 256.465 l
-C
-S 0.6 0 0.2 r f R
-
-N
-511.171 517.816 m
-542.777 447.578 546.492 439.32 571.539 383.664 c
-[] 0 d 0.6 1 0.2 r 3.58125 w s
-
-N
-577.906 369.508 m
-566.64 381.457 l
-576.437 385.867 l
-577.906 369.508 l
-C
-S 0.6 1 0.2 r f R
-
-N
-391.773 374.539 m
-477.171 367.422 486.218 366.668 555.441 360.902 c
-[] 0 d 0.8 1 0.1 r 3.58125 w s
-
-N
-570.91 359.613 m
-554.996 355.547 l
-555.886 366.254 l
-570.91 359.613 l
-C
-S 0.8 1 0.1 r f R
-
-N
-431.574 509.855 m
-459.371 512.637 467.753 513.473 483.843 515.082 c
-[] 0 d 0.2 0 0.4 r 3.58125 w s
-
-N
-499.289 516.625 m
-484.379 509.734 l
-483.312 520.43 l
-499.289 516.625 l
-C
-S 0.2 0 0.4 r f R
-
-N
-391.773 374.539 m
-444.953 438.351 450.761 445.324 493.589 496.719 c
-[] 0 d 1 0 0 r 3.58125 w s
-
-N
-503.527 508.645 m
-497.718 493.277 l
-489.461 500.16 l
-503.527 508.645 l
-C
-S 1 0 0 r f R
-
-N
-391.773 374.539 m
-408.687 432.043 411.226 440.676 423.824 483.508 c
-[] 0 d 1 0 0 r 3.58125 w s
-
-N
-428.203 498.402 m
-428.98 481.992 l
-418.668 485.027 l
-428.203 498.402 l
-C
-S 1 0 0 r f R
-
-N
-168.906 454.137 m
-290.617 479.953 299.535 481.848 404.711 504.156 c
-[] 0 d 0.2 0 0.4 r 3.58125 w s
-
-N
-419.894 507.379 m
-405.824 498.902 l
-403.593 509.414 l
-419.894 507.379 l
-C
-S 0.2 0 0.4 r f R
-
-N
-391.773 374.539 m
-407.699 333.133 410.879 324.863 421.714 296.695 c
-[] 0 d 0.6 1 0.2 r 3.58125 w s
-
-N
-427.285 282.207 m
-416.699 294.766 l
-426.73 298.621 l
-427.285 282.207 l
-C
-S 0.6 1 0.2 r f R
-
-N
-168.906 454.137 m
-270.98 417.683 279.554 414.621 365.918 383.777 c
-[] 0 d 0.4 0 0.3 r 3.58125 w s
-
-N
-380.535 378.555 m
-364.109 378.715 l
-367.722 388.836 l
-380.535 378.555 l
-C
-S 0.4 0 0.3 r f R
-
-N
-431.574 271.062 m
-392.125 238.195 385.261 232.473 357.156 209.051 c
-[] 0 d 1 1 0 r 3.58125 w s
-
-N
-345.23 199.113 m
-353.714 213.176 l
-360.597 204.922 l
-345.23 199.113 l
-C
-S 1 1 0 r f R
-
-N
-336.058 191.469 m
-359.39 141.473 363.183 133.348 380.164 96.9571 c
-[] 0 d 0.8 1 0.1 r 3.58125 w s
-
-N
-386.726 82.8941 m
-375.296 94.6871 l
-385.031 99.2301 l
-386.726 82.8941 l
-C
-S 0.8 1 0.1 r f R
-
-N
-208.703 143.711 m
-290.812 111.582 299.269 108.273 366.203 82.0821 c
-[] 0 d 0.6 1 0.2 r 3.58125 w s
-
-N
-380.66 76.4261 m
-364.246 77.0781 l
-368.164 87.0861 l
-380.66 76.4261 l
-C
-S 0.6 1 0.2 r f R
-
-N
-336.058 191.469 m
-297.203 262.117 292.836 270.055 261.738 326.598 c
-[] 0 d 0.8 0 0.1 r 3.58125 w s
-
-N
-254.257 340.199 m
-266.445 329.187 l
-257.031 324.008 l
-254.257 340.199 l
-C
-S 0.8 0 0.1 r f R
-
-N
-208.703 143.711 m
-262.414 163.851 270.82 167.004 310.347 181.828 c
-[] 0 d 0.6 0 0.2 r 3.58125 w s
-
-N
-324.879 187.277 m
-312.234 176.797 l
-308.461 186.859 l
-324.879 187.277 l
-C
-S 0.6 0 0.2 r f R
-
-N
-248.503 350.66 m
-215.242 393.902 209.777 401.004 185.648 432.371 c
-[] 0 d 0.4 0 0.3 r 3.58125 w s
-
-N
-176.183 444.672 m
-189.906 435.644 l
-181.39 429.094 l
-176.183 444.672 l
-C
-S 0.4 0 0.3 r f R
-
-N
-25.6286 271.062 m
-127.703 307.519 136.281 310.582 222.64 341.426 c
-[] 0 d 0.6 0 0.2 r 3.58125 w s
-
-N
-237.257 346.644 m
-224.449 336.363 l
-220.836 346.484 l
-237.257 346.644 l
-C
-S 0.6 0 0.2 r f R
-
-N
-25.6286 271.062 m
-91.1325 354.758 96.7419 361.93 151.98 432.512 c
-[] 0 d 0.2 0 0.4 r 3.58125 w s
-
-N
-161.546 444.734 m
-156.211 429.199 l
-147.75 435.824 l
-161.546 444.734 l
-C
-S 0.2 0 0.4 r f R
-
-N
-25.6286 271.062 m
-108.961 213.098 116.433 207.902 186.16 159.394 c
-[] 0 d 0.2 0 0.4 r 3.58125 w s
-
-N
-198.902 150.531 m
-183.093 154.984 l
-189.23 163.805 l
-198.902 150.531 l
-C
-S 0.2 0 0.4 r f R
-
-N
-37.57 271.062 m
-37.57 271.062 l
-37.57 277.656 32.2224 283.004 25.6286 283.004 c
-19.0348 283.004 13.6911 277.656 13.6911 271.062 c
-13.6911 264.469 19.0348 259.129 25.6286 259.129 c
-32.2224 259.129 37.57 264.469 37.57 271.062 c
-37.57 271.062 l
-C
-S 0 0 0 r f R
-
-N
-36.4841 271.062 m
-36.4841 271.062 l
-36.4841 277.055 31.6208 281.918 25.6286 281.918 c
-19.6364 281.918 14.777 277.055 14.777 271.062 c
-14.777 265.07 19.6364 260.215 25.6286 260.215 c
-31.6208 260.215 36.4841 265.07 36.4841 271.062 c
-36.4841 271.062 l
-C
-S 0 0 1 r f R
-
-N
-220.644 143.711 m
-220.644 143.711 l
-220.644 150.305 215.296 155.652 208.703 155.652 c
-202.109 155.652 196.765 150.305 196.765 143.711 c
-196.765 137.117 202.109 131.773 208.703 131.773 c
-215.296 131.773 220.644 137.117 220.644 143.711 c
-220.644 143.711 l
-C
-S 0 0 0 r f R
-
-N
-219.558 143.711 m
-219.558 143.711 l
-219.558 149.703 214.695 154.566 208.703 154.566 c
-202.711 154.566 197.851 149.703 197.851 143.711 c
-197.851 137.719 202.711 132.855 208.703 132.855 c
-214.695 132.855 219.558 137.719 219.558 143.711 c
-219.558 143.711 l
-C
-S 0.5 0.5 1 r f R
-
-N
-180.843 454.137 m
-180.843 454.137 l
-180.843 460.73 175.5 466.078 168.906 466.078 c
-162.312 466.078 156.964 460.73 156.964 454.137 c
-156.964 447.543 162.312 442.199 168.906 442.199 c
-175.5 442.199 180.843 447.543 180.843 454.137 c
-180.843 454.137 l
-C
-S 0 0 0 r f R
-
-N
-179.757 454.137 m
-179.757 454.137 l
-179.757 460.129 174.898 464.992 168.906 464.992 c
-162.914 464.992 158.05 460.129 158.05 454.137 c
-158.05 448.144 162.914 443.281 168.906 443.281 c
-174.898 443.281 179.757 448.144 179.757 454.137 c
-179.757 454.137 l
-C
-S 0.5 0.5 1 r f R
-
-N
-260.441 350.66 m
-260.441 350.66 l
-260.441 357.254 255.097 362.601 248.503 362.601 c
-241.91 362.601 236.562 357.254 236.562 350.66 c
-236.562 344.066 241.91 338.723 248.503 338.723 c
-255.097 338.723 260.441 344.066 260.441 350.66 c
-260.441 350.66 l
-C
-S 0 0 0 r f R
-
-N
-259.359 350.66 m
-259.359 350.66 l
-259.359 356.652 254.496 361.516 248.503 361.516 c
-242.511 361.516 237.648 356.652 237.648 350.66 c
-237.648 344.668 242.511 339.805 248.503 339.805 c
-254.496 339.805 259.359 344.668 259.359 350.66 c
-259.359 350.66 l
-C
-S 0.5 0.5 1 r f R
-
-N
-348 191.469 m
-348 191.469 l
-348 198.062 342.652 203.41 336.058 203.41 c
-329.464 203.41 324.121 198.062 324.121 191.469 c
-324.121 184.875 329.464 179.531 336.058 179.531 c
-342.652 179.531 348 184.875 348 191.469 c
-348 191.469 l
-C
-S 0 0 0 r f R
-
-N
-346.914 191.469 m
-346.914 191.469 l
-346.914 197.461 342.05 202.324 336.058 202.324 c
-330.066 202.324 325.207 197.461 325.207 191.469 c
-325.207 185.476 330.066 180.617 336.058 180.617 c
-342.05 180.617 346.914 185.476 346.914 191.469 c
-346.914 191.469 l
-C
-S 0.5 0.5 1 r f R
-
-N
-403.714 72.0741 m
-403.714 72.0741 l
-403.714 78.6681 398.367 84.0121 391.773 84.0121 c
-385.179 84.0121 379.839 78.6681 379.839 72.0741 c
-379.839 65.4801 385.179 60.1331 391.773 60.1331 c
-398.367 60.1331 403.714 65.4801 403.714 72.0741 c
-403.714 72.0741 l
-C
-S 0 0 0 r f R
-
-N
-402.629 72.0741 m
-402.629 72.0741 l
-402.629 78.0661 397.765 82.9301 391.773 82.9301 c
-385.781 82.9301 380.925 78.0661 380.925 72.0741 c
-380.925 66.0821 385.781 61.2191 391.773 61.2191 c
-397.765 61.2191 402.629 66.0821 402.629 72.0741 c
-402.629 72.0741 l
-C
-S 1 0.5 1 r f R
-
-N
-443.511 271.062 m
-443.511 271.062 l
-443.511 277.656 438.168 283.004 431.574 283.004 c
-424.98 283.004 419.632 277.656 419.632 271.062 c
-419.632 264.469 424.98 259.129 431.574 259.129 c
-438.168 259.129 443.511 264.469 443.511 271.062 c
-443.511 271.062 l
-C
-S 0 0 0 r f R
-
-N
-442.425 271.062 m
-442.425 271.062 l
-442.425 277.055 437.566 281.918 431.574 281.918 c
-425.582 281.918 420.718 277.055 420.718 271.062 c
-420.718 265.07 425.582 260.215 431.574 260.215 c
-437.566 260.215 442.425 265.07 442.425 271.062 c
-442.425 271.062 l
-C
-S 1 0.5 1 r f R
-
-N
-403.714 374.539 m
-403.714 374.539 l
-403.714 381.133 398.367 386.48 391.773 386.48 c
-385.179 386.48 379.839 381.133 379.839 374.539 c
-379.839 367.945 385.179 362.601 391.773 362.601 c
-398.367 362.601 403.714 367.945 403.714 374.539 c
-403.714 374.539 l
-C
-S 0 0 0 r f R
-
-N
-402.629 374.539 m
-402.629 374.539 l
-402.629 380.531 397.765 385.394 391.773 385.394 c
-385.781 385.394 380.925 380.531 380.925 374.539 c
-380.925 368.547 385.781 363.687 391.773 363.687 c
-397.765 363.687 402.629 368.547 402.629 374.539 c
-402.629 374.539 l
-C
-S 0.5 0.5 1 r f R
-
-N
-443.511 509.855 m
-443.511 509.855 l
-443.511 516.449 438.168 521.797 431.574 521.797 c
-424.98 521.797 419.632 516.449 419.632 509.855 c
-419.632 503.262 424.98 497.918 431.574 497.918 c
-438.168 497.918 443.511 503.262 443.511 509.855 c
-443.511 509.855 l
-C
-S 0 0 0 r f R
-
-N
-442.425 509.855 m
-442.425 509.855 l
-442.425 515.848 437.566 520.711 431.574 520.711 c
-425.582 520.711 420.718 515.848 420.718 509.855 c
-420.718 503.863 425.582 499 431.574 499 c
-437.566 499 442.425 503.863 442.425 509.855 c
-442.425 509.855 l
-C
-S 0.5 0.5 1 r f R
-
-N
-523.109 517.816 m
-523.109 517.816 l
-523.109 524.41 517.765 529.754 511.171 529.754 c
-504.578 529.754 499.23 524.41 499.23 517.816 c
-499.23 511.223 504.578 505.875 511.171 505.875 c
-517.765 505.875 523.109 511.223 523.109 517.816 c
-523.109 517.816 l
-C
-S 0 0 0 r f R
-
-N
-522.023 517.816 m
-522.023 517.816 l
-522.023 523.809 517.164 528.668 511.171 528.668 c
-505.179 528.668 500.316 523.809 500.316 517.816 c
-500.316 511.824 505.179 506.961 511.171 506.961 c
-517.164 506.961 522.023 511.824 522.023 517.816 c
-522.023 517.816 l
-C
-S 0.5 0.5 1 r f R
-
-N
-594.746 358.621 m
-594.746 358.621 l
-594.746 365.215 589.402 370.558 582.808 370.558 c
-576.214 370.558 570.867 365.215 570.867 358.621 c
-570.867 352.027 576.214 346.68 582.808 346.68 c
-589.402 346.68 594.746 352.027 594.746 358.621 c
-594.746 358.621 l
-C
-S 0 0 0 r f R
-
-N
-593.664 358.621 m
-593.664 358.621 l
-593.664 364.613 588.8 369.473 582.808 369.473 c
-576.816 369.473 571.953 364.613 571.953 358.621 c
-571.953 352.629 576.816 347.766 582.808 347.766 c
-588.8 347.766 593.664 352.629 593.664 358.621 c
-593.664 358.621 l
-C
-S 1 0.5 1 r f R
-
-N
-586.789 255.148 m
-586.789 255.148 l
-586.789 261.742 581.441 267.082 574.847 267.082 c
-568.253 267.082 562.91 261.742 562.91 255.148 c
-562.91 248.555 568.253 243.207 574.847 243.207 c
-581.441 243.207 586.789 248.555 586.789 255.148 c
-586.789 255.148 l
-C
-S 0 0 0 r f R
-
-N
-585.703 255.148 m
-585.703 255.148 l
-585.703 261.14 580.839 265.996 574.847 265.996 c
-568.855 265.996 563.992 261.14 563.992 255.148 c
-563.992 249.156 568.855 244.293 574.847 244.293 c
-580.839 244.293 585.703 249.156 585.703 255.148 c
-585.703 255.148 l
-C
-S 1 0.5 1 r f R
-
-N
-507.191 143.711 m
-507.191 143.711 l
-507.191 150.305 501.843 155.652 495.25 155.652 c
-488.656 155.652 483.312 150.305 483.312 143.711 c
-483.312 137.117 488.656 131.773 495.25 131.773 c
-501.843 131.773 507.191 137.117 507.191 143.711 c
-507.191 143.711 l
-C
-S 0 0 0 r f R
-
-N
-506.105 143.711 m
-506.105 143.711 l
-506.105 149.703 501.242 154.566 495.25 154.566 c
-489.257 154.566 484.398 149.703 484.398 143.711 c
-484.398 137.719 489.257 132.855 495.25 132.855 c
-501.242 132.855 506.105 137.719 506.105 143.711 c
-506.105 143.711 l
-C
-S 1 0.5 1 r f R
-
-N
-618.629 72.0741 m
-618.629 72.0741 l
-618.629 78.6681 613.281 84.0121 606.687 84.0121 c
-600.093 84.0121 594.746 78.6681 594.746 72.0741 c
-594.746 65.4801 600.093 60.1331 606.687 60.1331 c
-613.281 60.1331 618.629 65.4801 618.629 72.0741 c
-618.629 72.0741 l
-C
-S 0 0 0 r f R
-
-N
-617.543 72.0741 m
-617.543 72.0741 l
-617.543 78.0661 612.679 82.9301 606.687 82.9301 c
-600.695 82.9301 595.832 78.0661 595.832 72.0741 c
-595.832 66.0821 600.695 61.2191 606.687 61.2191 c
-612.679 61.2191 617.543 66.0821 617.543 72.0741 c
-617.543 72.0741 l
-C
-S 1 0.5 1 r f R
-
-N
-666.386 135.75 m
-666.386 135.75 l
-666.386 142.344 661.039 147.691 654.445 147.691 c
-647.851 147.691 642.507 142.344 642.507 135.75 c
-642.507 129.156 647.851 123.812 654.445 123.812 c
-661.039 123.812 666.386 129.156 666.386 135.75 c
-666.386 135.75 l
-C
-S 0 0 0 r f R
-
-N
-665.3 135.75 m
-665.3 135.75 l
-665.3 141.742 660.437 146.605 654.445 146.605 c
-648.453 146.605 643.589 141.742 643.589 135.75 c
-643.589 129.758 648.453 124.898 654.445 124.898 c
-660.437 124.898 665.3 129.758 665.3 135.75 c
-665.3 135.75 l
-C
-S 1 0.5 1 r f R
-
-N
-690.265 199.43 m
-690.265 199.43 l
-690.265 206.023 684.918 211.371 678.324 211.371 c
-671.73 211.371 666.386 206.023 666.386 199.43 c
-666.386 192.836 671.73 187.488 678.324 187.488 c
-684.918 187.488 690.265 192.836 690.265 199.43 c
-690.265 199.43 l
-C
-S 0 0 0 r f R
-
-N
-689.179 199.43 m
-689.179 199.43 l
-689.179 205.422 684.316 210.285 678.324 210.285 c
-672.332 210.285 667.472 205.422 667.472 199.43 c
-667.472 193.437 672.332 188.574 678.324 188.574 c
-684.316 188.574 689.179 193.437 689.179 199.43 c
-689.179 199.43 l
-C
-S 1 0.5 1 r f R
-
-N
-706.183 342.699 m
-706.183 342.699 l
-706.183 349.293 700.836 354.64 694.242 354.64 c
-687.648 354.64 682.304 349.293 682.304 342.699 c
-682.304 336.105 687.648 330.762 694.242 330.762 c
-700.836 330.762 706.183 336.105 706.183 342.699 c
-706.183 342.699 l
-C
-S 0 0 0 r f R
-
-N
-705.097 342.699 m
-705.097 342.699 l
-705.097 348.691 700.234 353.555 694.242 353.555 c
-688.25 353.555 683.39 348.691 683.39 342.699 c
-683.39 336.707 688.25 331.848 694.242 331.848 c
-700.234 331.848 705.097 336.707 705.097 342.699 c
-705.097 342.699 l
-C
-S 1 0.5 1 r f R
-
-N
-825.578 151.672 m
-825.578 151.672 l
-825.578 158.266 820.234 163.609 813.64 163.609 c
-807.046 163.609 801.699 158.266 801.699 151.672 c
-801.699 145.078 807.046 139.73 813.64 139.73 c
-820.234 139.73 825.578 145.078 825.578 151.672 c
-825.578 151.672 l
-C
-S 0 0 0 r f R
-
-N
-824.496 151.672 m
-824.496 151.672 l
-824.496 157.664 819.632 162.523 813.64 162.523 c
-807.648 162.523 802.785 157.664 802.785 151.672 c
-802.785 145.68 807.648 140.816 813.64 140.816 c
-819.632 140.816 824.496 145.68 824.496 151.672 c
-824.496 151.672 l
-C
-S 1 0 1 r f R
-
-%%EOF
diff -r 26983135fd6d -r 57143c09dc20 doc/images/flow.png
Binary file doc/images/flow.png has changed
diff -r 26983135fd6d -r 57143c09dc20 lemon/Makefile.am
--- a/lemon/Makefile.am	Wed Nov 14 17:53:08 2007 +0000
+++ b/lemon/Makefile.am	Sat Nov 17 20:58:11 2007 +0000
@@ -52,9 +52,11 @@
 	lemon/dag_shortest_path.h \
 	lemon/dfs.h \
 	lemon/dijkstra.h \
+	lemon/dinitz_sleator_tarjan.h \
 	lemon/dist_log.h \
 	lemon/dim2.h \
 	lemon/dimacs.h \
+	lemon/dynamic_tree.h \
 	lemon/edge_set.h \
 	lemon/edmonds_karp.h \
 	lemon/elevator.h \
@@ -71,6 +73,7 @@
 	lemon/graph_utils.h \
 	lemon/graph_writer.h \
 	lemon/grid_ugraph.h \
+	lemon/goldberg_tarjan.h \
 	lemon/hao_orlin.h \
 	lemon/hypercube_graph.h \
 	lemon/iterable_maps.h \
diff -r 26983135fd6d -r 57143c09dc20 lemon/dinitz_sleator_tarjan.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/dinitz_sleator_tarjan.h	Sat Nov 17 20:58:11 2007 +0000
@@ -0,0 +1,762 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2007
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_DINITZ_SLEATOR_TARJAN_H
+#define LEMON_DINITZ_SLEATOR_TARJAN_H
+
+/// \file 
+/// \ingroup max_flow 
+/// \brief Implementation the dynamic tree data structure of Sleator
+/// and Tarjan.
+
+#include <lemon/time_measure.h>
+#include <queue>
+#include <lemon/graph_utils.h>
+#include <lemon/tolerance.h>
+#include <lemon/graph_adaptor.h>
+#include <lemon/bfs.h>
+#include <lemon/edge_set.h>
+#include <lemon/dynamic_tree.h>
+
+#include <vector>
+#include <limits>
+#include <fstream>
+
+
+namespace lemon {
+
+  /// \brief Default traits class of DinitzSleatorTarjan class.
+  ///
+  /// Default traits class of DinitzSleatorTarjan class.
+  /// \param _Graph Graph type.
+  /// \param _CapacityMap Type of capacity map.
+  template <typename _Graph, typename _CapacityMap>
+  struct DinitzSleatorTarjanDefaultTraits {
+
+    /// \brief The graph type the algorithm runs on. 
+    typedef _Graph Graph;
+
+    /// \brief The type of the map that stores the edge capacities.
+    ///
+    /// The type of the map that stores the edge capacities.
+    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
+    typedef _CapacityMap CapacityMap;
+
+    /// \brief The type of the length of the edges.
+    typedef typename CapacityMap::Value Value;
+
+    /// \brief The map type that stores the flow values.
+    ///
+    /// The map type that stores the flow values. 
+    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
+    typedef typename Graph::template EdgeMap<Value> FlowMap;
+
+    /// \brief Instantiates a FlowMap.
+    ///
+    /// This function instantiates a \ref FlowMap. 
+    /// \param graph The graph, to which we would like to define the flow map.
+    static FlowMap* createFlowMap(const Graph& graph) {
+      return new FlowMap(graph);
+    }
+
+    /// \brief The tolerance used by the algorithm
+    ///
+    /// The tolerance used by the algorithm to handle inexact computation.
+    typedef Tolerance<Value> Tolerance;
+
+  };
+
+  /// \ingroup max_flow
+  ///
+  /// \brief Dinitz-Sleator-Tarjan algorithms class.
+  ///
+  /// This class provides an implementation of the \e
+  /// Dinitz-Sleator-Tarjan \e algorithm producing a flow of maximum
+  /// value in a directed graph. The DinitzSleatorTarjan algorithm is
+  /// the fastest known max flow algorithms wich using blocking
+  /// flow. It is an improvement of the Dinitz algorithm by using the
+  /// \ref DynamicTree "dynamic tree" data structure of Sleator and
+  /// Tarjan.
+  ///
+  /// This blocking flow algorithms builds a layered graph according
+  /// to \ref Bfs "breadth-first search" distance from the target node
+  /// in the reversed residual graph. The layered graph contains each
+  /// residual edge which steps one level down. After that the
+  /// algorithm constructs a blocking flow on the layered graph and
+  /// augments the overall flow with it. The number of the levels in
+  /// the layered graph is strictly increasing in each augmenting
+  /// phase therefore the number of the augmentings is at most
+  /// \f$n-1\f$.  The length of each phase is at most
+  /// \f$O(m\log(n))\f$, that the overall time complexity is
+  /// \f$O(nm\log(n))\f$.
+  ///
+  /// \param _Graph The directed graph type the algorithm runs on.
+  /// \param _CapacityMap The capacity map type.
+  /// \param _Traits Traits class to set various data types used by
+  /// the algorithm.  The default traits class is \ref
+  /// DinitzSleatorTarjanDefaultTraits.  See \ref
+  /// DinitzSleatorTarjanDefaultTraits for the documentation of a
+  /// Dinitz-Sleator-Tarjan traits class.
+  ///
+  /// \author Tamas Hamori and Balazs Dezso
+#ifdef DOXYGEN
+  template <typename _Graph, typename _CapacityMap, typename _Traits>
+#else
+  template <typename _Graph, 
+	    typename _CapacityMap = typename _Graph::template EdgeMap<int>,
+	    typename _Traits = 
+	    DinitzSleatorTarjanDefaultTraits<_Graph, _CapacityMap> >
+#endif
+  class DinitzSleatorTarjan {
+  public:
+
+    typedef _Traits Traits;
+    typedef typename Traits::Graph Graph;
+    typedef typename Traits::CapacityMap CapacityMap;
+    typedef typename Traits::Value Value; 
+
+    typedef typename Traits::FlowMap FlowMap;
+    typedef typename Traits::Tolerance Tolerance;
+
+
+  private:
+
+    GRAPH_TYPEDEFS(typename Graph);
+
+
+    typedef typename Graph::template NodeMap<int> LevelMap;
+    typedef typename Graph::template NodeMap<int> IntNodeMap;
+    typedef typename Graph::template NodeMap<Edge> EdgeNodeMap;
+    typedef DynamicTree<Value, IntNodeMap, Tolerance, false> DynTree;
+
+  private:
+    
+    const Graph& _graph;
+    const CapacityMap* _capacity;
+
+    Node _source, _target;
+
+    FlowMap* _flow;
+    bool _local_flow;
+
+    IntNodeMap* _level;
+    EdgeNodeMap* _dt_edges;
+    
+    IntNodeMap* _dt_index;
+    DynTree* _dt;
+
+    Tolerance _tolerance;
+    
+    Value _flow_value;
+    Value _max_value;
+
+
+  public:
+
+    typedef DinitzSleatorTarjan Create;
+
+    ///\name Named template parameters
+
+    ///@{
+
+    template <typename _FlowMap>
+    struct DefFlowMapTraits : public Traits {
+      typedef _FlowMap FlowMap;
+      static FlowMap *createFlowMap(const Graph&) {
+	throw UninitializedParameter();
+      }
+    };
+
+    /// \brief \ref named-templ-param "Named parameter" for setting
+    /// FlowMap type
+    ///
+    /// \ref named-templ-param "Named parameter" for setting FlowMap
+    /// type
+    template <typename _FlowMap>
+    struct DefFlowMap 
+      : public DinitzSleatorTarjan<Graph, CapacityMap, 
+			      DefFlowMapTraits<_FlowMap> > {
+      typedef DinitzSleatorTarjan<Graph, CapacityMap, 
+			     DefFlowMapTraits<_FlowMap> > Create;
+    };
+
+    template <typename _Elevator>
+    struct DefElevatorTraits : public Traits {
+      typedef _Elevator Elevator;
+      static Elevator *createElevator(const Graph&, int) {
+	throw UninitializedParameter();
+      }
+    };
+
+    /// @}
+
+    /// \brief \ref Exception for the case when the source equals the target.
+    ///
+    /// \ref Exception for the case when the source equals the target.
+    ///
+    class InvalidArgument : public lemon::LogicError {
+    public:
+      virtual const char* what() const throw() {
+	return "lemon::DinitzSleatorTarjan::InvalidArgument";
+      }
+    };
+
+    /// \brief The constructor of the class.
+    ///
+    /// The constructor of the class. 
+    /// \param graph The directed graph the algorithm runs on. 
+    /// \param capacity The capacity of the edges. 
+    /// \param source The source node.
+    /// \param target The target node.
+    DinitzSleatorTarjan(const Graph& graph, const CapacityMap& capacity,
+			Node source, Node target)
+      : _graph(graph), _capacity(&capacity),
+	_source(source), _target(target),
+	_flow(0), _local_flow(false),
+	_level(0), _dt_edges(0),
+	_dt_index(0), _dt(0),
+	_tolerance(), _flow_value(), _max_value()
+    {
+      if (_source == _target) {
+	throw InvalidArgument();
+      }
+    }
+
+    /// \brief Destrcutor.
+    ///
+    /// Destructor.
+    ~DinitzSleatorTarjan() {
+      destroyStructures();
+    }
+
+    /// \brief Sets the capacity map.
+    ///
+    /// Sets the capacity map.
+    /// \return \c (*this)
+    DinitzSleatorTarjan& capacityMap(const CapacityMap& map) {
+      _capacity = &map;
+      return *this;
+    }
+
+    /// \brief Sets the flow map.
+    ///
+    /// Sets the flow map.
+    /// \return \c (*this)
+    DinitzSleatorTarjan& flowMap(FlowMap& map) {
+      if (_local_flow) {
+	delete _flow;
+	_local_flow = false;
+      }
+      _flow = &map;
+      return *this;
+    }
+
+    /// \brief Returns the flow map.
+    ///
+    /// \return The flow map.
+    const FlowMap& flowMap() {
+      return *_flow;
+    }
+
+    /// \brief Sets the source node.
+    ///
+    /// Sets the source node.
+    /// \return \c (*this)
+    DinitzSleatorTarjan& source(const Node& node) {
+      _source = node;
+      return *this;
+    }
+
+    /// \brief Sets the target node.
+    ///
+    /// Sets the target node.
+    /// \return \c (*this)
+    DinitzSleatorTarjan& target(const Node& node) {
+      _target = node;
+      return *this;
+    }
+
+    /// \brief Sets the tolerance used by algorithm.
+    ///
+    /// Sets the tolerance used by algorithm.
+    DinitzSleatorTarjan& tolerance(const Tolerance& tolerance) const {
+      _tolerance = tolerance;
+      if (_dt) {
+	_dt.tolerance(_tolerance);
+      }
+      return *this;
+    } 
+
+    /// \brief Returns the tolerance used by algorithm.
+    ///
+    /// Returns the tolerance used by algorithm.
+    const Tolerance& tolerance() const {
+      return tolerance;
+    } 
+
+  private:
+        
+    void createStructures() {
+      if (!_flow) {
+	_flow = Traits::createFlowMap(_graph);
+	_local_flow = true;
+      }
+      if (!_level) {
+	_level = new LevelMap(_graph);
+      }
+      if (!_dt_index && !_dt) {
+	_dt_index = new IntNodeMap(_graph);
+	_dt = new DynTree(*_dt_index, _tolerance);
+      }
+      if (!_dt_edges) {
+	_dt_edges = new EdgeNodeMap(_graph);
+      }
+      _max_value = _dt->maxValue();
+    }
+
+    void destroyStructures() {
+      if (_local_flow) {
+	delete _flow;
+      }
+      if (_level) {
+	delete _level;
+      }
+      if (_dt) {
+	delete _dt;
+      }
+      if (_dt_index) {
+	delete _dt_index;
+      }
+      if (_dt_edges) {
+	delete _dt_edges;
+      }
+    }
+
+    bool createLayeredGraph() {
+
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+	_level->set(n, -2);
+      }
+      
+      int level = 0;
+
+      std::vector<Node> queue;
+      queue.push_back(_target);
+      _level->set(_target, level);
+      
+      while ((*_level)[_source] == -2 && !queue.empty()) {
+	std::vector<Node> nqueue;
+	++level;
+	
+	for (int i = 0; i < int(queue.size()); ++i) {
+	  Node n = queue[i];
+	  
+	  for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+	    Node v = _graph.target(e);
+	    if ((*_level)[v] != -2) continue;
+	    Value rem = (*_flow)[e];
+	    if (!_tolerance.positive(rem)) continue;
+	    _level->set(v, level);
+	    nqueue.push_back(v);
+	  }
+
+	  for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+	    Node v = _graph.source(e);
+	    if ((*_level)[v] != -2) continue;
+	    Value rem = (*_capacity)[e] - (*_flow)[e];
+	    if (!_tolerance.positive(rem)) continue;
+	    _level->set(v, level);
+	    nqueue.push_back(v);
+	  }
+
+	}
+	queue.swap(nqueue);
+      }
+      
+      return (*_level)[_source] != -2;
+    }
+
+    void initEdges() {
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+	_graph.firstOut((*_dt_edges)[n], n);
+      }
+    }
+        
+    
+    void augmentPath() {
+      Value rem;
+      Node n = _dt->findCost(_source, rem);
+      _flow_value += rem;
+      _dt->addCost(_source, - rem);
+
+      _dt->cut(n);
+      _dt->addCost(n, _max_value);
+
+      Edge e = (*_dt_edges)[n];
+      if (_graph.source(e) == n) {
+	_flow->set(e, (*_capacity)[e]);
+	
+	_graph.nextOut(e);
+	if (e == INVALID) {
+	  _graph.firstIn(e, n);
+	}
+      } else {
+	_flow->set(e, 0);
+	_graph.nextIn(e);
+      }
+      _dt_edges->set(n, e);
+
+    }
+
+    bool advance(Node n) {
+      Edge e = (*_dt_edges)[n];
+      if (e == INVALID) return false;
+
+      Node u;
+      Value rem;      
+      if (_graph.source(e) == n) {
+	u = _graph.target(e);
+	while ((*_level)[n] != (*_level)[u] + 1 || 
+	       !_tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
+	  _graph.nextOut(e);
+	  if (e == INVALID) break;
+	  u = _graph.target(e);
+	}
+	if (e != INVALID) {
+	  rem = (*_capacity)[e] - (*_flow)[e];
+	} else {
+	  _graph.firstIn(e, n);
+	  if (e == INVALID) {
+	    _dt_edges->set(n, INVALID);
+	    return false;
+	  }
+	  u = _graph.source(e);
+	  while ((*_level)[n] != (*_level)[u] + 1 ||
+		 !_tolerance.positive((*_flow)[e])) {
+	    _graph.nextIn(e);
+	    if (e == INVALID) {
+	      _dt_edges->set(n, INVALID);
+	      return false;
+	    }
+	    u = _graph.source(e);
+	  }
+	  rem = (*_flow)[e];
+	}
+      } else {
+	u = _graph.source(e);
+	while ((*_level)[n] != (*_level)[u] + 1 ||
+	       !_tolerance.positive((*_flow)[e])) {
+	  _graph.nextIn(e);
+	  if (e == INVALID) {
+	    _dt_edges->set(n, INVALID);
+	    return false;
+	  }
+	  u = _graph.source(e);
+	}
+	rem = (*_flow)[e];
+      }
+
+      _dt->addCost(n, - std::numeric_limits<Value>::max());
+      _dt->addCost(n, rem);
+      _dt->link(n, u);
+      _dt_edges->set(n, e);
+      return true;
+    }
+
+    void retreat(Node n) {
+      _level->set(n, -1);
+      
+      for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+	Node u = _graph.target(e);
+	if ((*_dt_edges)[u] == e && _dt->findRoot(u) == n) {
+	  Value rem;
+	  _dt->findCost(u, rem);
+	  _flow->set(e, rem);
+	  _dt->cut(u);
+	  _dt->addCost(u, - rem);
+	  _dt->addCost(u, _max_value);
+	}
+      }
+      for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+	Node u = _graph.source(e);
+	if ((*_dt_edges)[u] == e && _dt->findRoot(u) == n) {
+	  Value rem;
+	  _dt->findCost(u, rem);
+	  _flow->set(e, (*_capacity)[e] - rem);
+	  _dt->cut(u);
+	  _dt->addCost(u, - rem);
+	  _dt->addCost(u, _max_value);
+	}
+      }
+    }
+
+    void extractTrees() {
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+	
+	Node w = _dt->findRoot(n);
+      
+	while (w != n) {
+      
+	  Value rem;      
+	  Node u = _dt->findCost(n, rem);
+
+	  _dt->cut(u);
+	  _dt->addCost(u, - rem);
+	  _dt->addCost(u, _max_value);
+	  
+	  Edge e = (*_dt_edges)[u];
+	  _dt_edges->set(u, INVALID);
+	  
+	  if (u == _graph.source(e)) {
+	    _flow->set(e, (*_capacity)[e] - rem);
+	  } else {
+	    _flow->set(e, rem);
+	  }
+	  
+	  w = _dt->findRoot(n);
+	}      
+      }
+    }
+
+
+  public:
+    
+    /// \name Execution control The simplest way to execute the
+    /// algorithm is to use the \c run() member functions.
+    /// \n
+    /// If you need more control on initial solution or
+    /// execution then you have to call one \ref init() function and then
+    /// the start() or multiple times the \c augment() member function.  
+    
+    ///@{
+
+    /// \brief Initializes the algorithm
+    /// 
+    /// It sets the flow to empty flow.
+    void init() {
+      createStructures();
+
+      _dt->clear();
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        _dt->makeTree(n);
+        _dt->addCost(n, _max_value);
+      }
+
+      for (EdgeIt it(_graph); it != INVALID; ++it) {
+        _flow->set(it, 0);
+      }
+      _flow_value = 0;
+    }
+    
+    /// \brief Initializes the algorithm
+    /// 
+    /// Initializes the flow to the \c flowMap. The \c flowMap should
+    /// contain a feasible flow, ie. in each node excluding the source
+    /// and the target the incoming flow should be equal to the
+    /// outgoing flow.
+    template <typename FlowMap>
+    void flowInit(const FlowMap& flowMap) {
+      createStructures();
+
+      _dt->clear();
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        _dt->makeTree(n);
+        _dt->addCost(n, _max_value);
+      }
+
+      for (EdgeIt e(_graph); e != INVALID; ++e) {
+	_flow->set(e, flowMap[e]);
+      }
+      _flow_value = 0;
+      for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
+        _flow_value += (*_flow)[jt];
+      }
+      for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
+        _flow_value -= (*_flow)[jt];
+      }
+    }
+
+    /// \brief Initializes the algorithm
+    /// 
+    /// Initializes the flow to the \c flowMap. The \c flowMap should
+    /// contain a feasible flow, ie. in each node excluding the source
+    /// and the target the incoming flow should be equal to the
+    /// outgoing flow.  
+    /// \return %False when the given flowMap does not contain
+    /// feasible flow.
+    template <typename FlowMap>
+    bool checkedFlowInit(const FlowMap& flowMap) {
+      createStructures();
+
+      _dt->clear();
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        _dt->makeTree(n);
+        _dt->addCost(n, _max_value);
+      }
+
+      for (EdgeIt e(_graph); e != INVALID; ++e) {
+	_flow->set(e, flowMap[e]);
+      }
+      for (NodeIt it(_graph); it != INVALID; ++it) {
+        if (it == _source || it == _target) continue;
+        Value outFlow = 0;
+        for (OutEdgeIt jt(_graph, it); jt != INVALID; ++jt) {
+          outFlow += (*_flow)[jt];
+        }
+        Value inFlow = 0;
+        for (InEdgeIt jt(_graph, it); jt != INVALID; ++jt) {
+          inFlow += (*_flow)[jt];
+        }
+        if (_tolerance.different(outFlow, inFlow)) {
+          return false;
+        }
+      }
+      for (EdgeIt it(_graph); it != INVALID; ++it) {
+        if (_tolerance.less((*_flow)[it], 0)) return false;
+        if (_tolerance.less((*_capacity)[it], (*_flow)[it])) return false;
+      }
+      _flow_value = 0;
+      for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
+        _flow_value += (*_flow)[jt];
+      }
+      for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
+        _flow_value -= (*_flow)[jt];
+      }
+      return true;
+    }
+
+    /// \brief Executes the algorithm
+    ///
+    /// It runs augmenting phases by adding blocking flow until the
+    /// optimal solution is reached.
+    void start() {
+      while (augment());
+    }
+
+    /// \brief Augments the flow with a blocking flow on a layered
+    /// graph.
+    /// 
+    /// This function builds a layered graph and then find a blocking
+    /// flow on this graph. The number of the levels in the layered
+    /// graph is strictly increasing in each augmenting phase
+    /// therefore the number of the augmentings is at most \f$ n-1
+    /// \f$.  The length of each phase is at most \f$ O(m \log(n))
+    /// \f$, that the overall time complexity is \f$ O(nm \log(n)) \f$.
+    /// \return %False when there is not residual path between the
+    /// source and the target so the current flow is a feasible and
+    /// optimal solution.
+    bool augment() {
+      Node n;
+
+      if (createLayeredGraph()) {
+	
+	Timer bf_timer;
+	initEdges();
+
+	n = _dt->findRoot(_source);
+	while (true) {
+	  Edge e;
+	  if (n == _target) {
+	    augmentPath();
+	  } else if (!advance(n)) {
+	    if (n != _source) {
+	      retreat(n);
+	    } else {
+	      break;
+	    }
+	  }
+	  n = _dt->findRoot(_source);
+	}     
+	extractTrees();
+
+	return true;
+      } else {
+	return false;
+      }
+    }
+    
+    /// \brief runs the algorithm.
+    /// 
+    /// It is just a shorthand for:
+    ///
+    ///\code 
+    /// ek.init();
+    /// ek.start();
+    ///\endcode
+    void run() {
+      init();
+      start();
+    }
+
+    /// @}
+
+    /// \name Query Functions
+    /// The result of the %Dijkstra algorithm can be obtained using these
+    /// functions.\n
+    /// Before the use of these functions,
+    /// either run() or start() must be called.
+    
+    ///@{
+
+    /// \brief Returns the value of the maximum flow.
+    ///
+    /// Returns the value of the maximum flow by returning the excess
+    /// of the target node \c t. This value equals to the value of
+    /// the maximum flow already after the first phase.
+    Value flowValue() const {
+      return _flow_value;
+    }
+
+
+    /// \brief Returns the flow on the edge.
+    ///
+    /// Sets the \c flowMap to the flow on the edges. This method can
+    /// be called after the second phase of algorithm.
+    Value flow(const Edge& edge) const {
+      return (*_flow)[edge];
+    }
+
+    /// \brief Returns true when the node is on the source side of minimum cut.
+    ///
+
+    /// Returns true when the node is on the source side of minimum
+    /// cut. This method can be called both after running \ref
+    /// startFirstPhase() and \ref startSecondPhase().
+    bool minCut(const Node& node) const {
+      return (*_level)[node] == -2;
+    }
+
+    /// \brief Returns a minimum value cut.
+    ///
+    /// Sets \c cut to the characteristic vector of a minimum value cut
+    /// It simply calls the minMinCut member.
+    /// \retval cut Write node bool map. 
+    template <typename CutMap>
+    void minCutMap(CutMap& cutMap) const {
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+	cutMap.set(n, (*_level)[n] == -2);
+      }
+      cutMap.set(_source, true);
+    }    
+
+    /// @}
+
+  };
+}
+
+#endif
diff -r 26983135fd6d -r 57143c09dc20 lemon/dynamic_tree.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/dynamic_tree.h	Sat Nov 17 20:58:11 2007 +0000
@@ -0,0 +1,1076 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2007
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+#ifndef LEMON_DYNAMIC_TREE_H
+#define LEMON_DYNAMIC_TREE_H
+
+/// \ingroup auxdata
+/// \file
+/// \brief The dynamic tree data structure of Sleator and Tarjan.
+///
+
+#include <vector>
+#include <limits>
+#include <lemon/tolerance.h>
+
+namespace lemon {
+
+  /// \ingroup auxdata
+  ///
+  /// \brief The dynamic tree data structure of Sleator and Tarjan.
+  ///
+  /// This class provides an implementation of the dynamic tree data
+  /// structure for maintaining a set of node-disjoint rooted
+  /// trees. Each item has an associated value, and the item with
+  /// minimum value can be find in \f$O(\log(n)\f$ on the path from a
+  /// node to the its root, and the items on such path can be
+  /// increased or decreased globally with a certain value in the same
+  /// running time. We regard a tree edge as directed toward the root,
+  /// that is from child to parent. Its structure can be modified by
+  /// two basic operations: \e link(v,w) adds an edge between a root v
+  /// and a node w in a different component; \e cut(v) removes the
+  /// edge between v and its parent.
+  /// 
+  /// \param _Value The value type of the items.  
+  /// \param _ItemIntMap Converts item type of node to integer.
+  /// \param _Tolerance The tolerance class to handle computation
+  /// problems.
+  /// \param _enableSize If true then the data structre manatain the
+  /// size of each tree. The feature is used in \ref GoldbergTarjan
+  /// algorithm. The default value is true.
+  ///
+  /// \author Hamori Tamas
+#ifdef DOXYGEN
+  template <typename _Value, typename _ItemIntMap, 
+	    typename _Tolerance, bool _enableSize>
+#else
+  template <typename _Value, typename _ItemIntMap, 
+	    typename _Tolerance = lemon::Tolerance<_Value>,
+	    bool _enableSize = true>
+#endif
+  class DynamicTree {
+  public:
+    /// \brief The integer map on the items.
+    typedef _ItemIntMap ItemIntMap;
+    /// \brief The item type of nodes.
+    typedef typename ItemIntMap::Key Item;
+    /// \brief The value type of the algorithms.
+    typedef _Value Value;
+    /// \brief The tolerance used by the algorithm.
+    typedef _Tolerance Tolerance;
+
+  private:
+  
+    class ItemData;
+
+    std::vector<ItemData> _data;
+    ItemIntMap &_iim;
+    Value _max_value;
+    Tolerance _tolerance;
+
+  public:
+    /// \brief The constructor of the class.
+    ///
+    /// \param iim The integer map on the items. 
+    /// \param tolerance Tolerance class.
+    DynamicTree(ItemIntMap &iim, const Tolerance& tolerance = Tolerance())
+      : _iim(iim), _max_value(std::numeric_limits<Value>::max()), 
+	_tolerance(tolerance) {}
+  
+    ~DynamicTree() {}
+
+    /// \brief Clears the data structure
+    ///
+    /// Clears the data structure
+    void clear() {
+      _data.clear();
+    }
+
+    /// \brief Sets the tolerance used by algorithm.
+    ///
+    /// Sets the tolerance used by algorithm.
+    void tolerance(const Tolerance& tolerance) const {
+      _tolerance = tolerance;
+      return *this;
+    } 
+  
+    /// \brief Returns the tolerance used by algorithm.
+    ///
+    /// Returns the tolerance used by algorithm.
+    const Tolerance& tolerance() const {
+      return tolerance;
+    } 
+  
+    /// \brief Create a new tree containing a single node with cost zero.
+    void makeTree(const Item &item) {
+      _data[makePath(item)].successor = -1;
+    }
+    
+    /// \brief Return the root of the tree containing node with itemtype
+    /// \e item.
+    Item findRoot(const Item &item) {
+      return _data[findTail(expose(_iim[item]))].id;
+    }
+    
+    /// \brief Return the the value of nodes in the tree containing
+    /// node with itemtype \e item.
+    int findSize(const Item &item) {
+      return _data[expose(_iim[item])].size;
+    }
+    
+    /// \brief Return the minimum cost containing node.
+    /// 
+    /// Return into \e d the minimum cost on the tree path from \e item
+    /// to findRoot(item).  Return the last item (closest to its root)
+    /// on this path of the minimum cost.
+    Item findCost(const Item &item, Value& d){
+      return _data[findPathCost(expose(_iim[item]),d)].id;
+    }
+    
+    /// \brief Add \e x value to the cost of every node on the path from
+    /// \e item to findRoot(item).
+    void addCost(const Item &item, Value x){
+      addPathCost(expose(_iim[item]), x);
+    }
+    
+    /// \brief Combine the trees containing nodes \e item1 and \e item2
+    /// by adding an edge from \e item1 \e item2.
+    /// 
+    /// This operation assumes that \e item1 is root and \e item2 is in
+    /// a different tree.
+    void link(const Item &item1, const Item &item2){
+      int v = _iim[item1];
+      int w = _iim[item2];
+      int p = expose(w);
+      join(-1, expose(v), p);
+      _data[v].successor = -1;
+      _data[v].size += _data[p].size;
+
+    }    
+    
+    /// \brief Divide the tree containing node \e item into two trees by
+    /// deleting the edge out of \e item.
+    /// 
+    /// This operation assumes that \e item is not a tree root.
+    void cut(const Item &item) {
+      int v = _iim[item];
+      int p, q;
+      expose(v);
+      split(p, v, q);
+      if (p != -1) {
+	_data[p].successor = v;
+      }
+      _data[v].size -= _data[q].size;
+      if (q != -1) {
+	_data[q].successor = _data[v].successor;
+      }
+      _data[v].successor = -1;
+    }
+
+    ///\brief 
+    Item parent(const Item &item){
+      return _data[_iim[item].p].id;
+    }
+
+    ///\brief Return the upper bound of the costs.
+    Value maxValue() const {
+      return _max_value;
+    }
+    
+  private:
+
+    int makePath(const Item &item) {
+      _iim.set(item, _data.size());
+      ItemData v(item);
+      _data.push_back(v);
+      return _iim[item];
+    }
+
+    int findPath(int v){
+      splay(v);
+      return v;
+    }
+    
+    int findPathCost(int p, Value &d){
+      while ((_data[p].right != -1 && 
+	      !_tolerance.less(0, _data[_data[p].right].dmin)) || 
+	     (_data[p].left != -1 && _tolerance.less(0, _data[p].dcost))) {
+	if (_data[p].right != -1 && 
+	    !_tolerance.less(0, _data[_data[p].right].dmin)) {
+	  p = _data[p].right;
+	} else if (_data[p].left != -1 && 
+		   !_tolerance.less(0, _data[_data[p].left].dmin)){
+	  p = _data[p].left;
+	}
+      }
+      splay(p);
+      d = _data[p].dmin;
+      return p; 
+    }
+
+    int findTail(int p){
+      while (_data[p].right != -1) {
+	p = _data[p].right;
+      }
+      splay(p);
+      return p;
+    }
+    
+    void addPathCost(int p, Value x){
+      if (!_tolerance.less(x, _max_value)) {
+	_data[p].dmin = x;_data[p].dcost = x;
+      } else if (!_tolerance.less(-x, _max_value)) {
+	_data[p].dmin = 0;
+	_data[p].dcost = 0;
+      } else {
+	_data[p].dmin += x;
+      }
+    }
+
+    void join(int p, int v, int q) {
+      Value min = _max_value;
+      Value pmin = _max_value;
+      Value vmin = _data[v].dmin;
+      Value qmin = _max_value;
+      if (p != -1){
+	pmin = _data[p].dmin;
+      }
+      if (q != -1){
+	qmin = _data[q].dmin;
+      }
+        
+      if (_tolerance.less(vmin, qmin)) {
+	if (_tolerance.less(vmin,pmin)) {
+	  min = vmin;
+	} else {
+	  min = pmin;
+	}
+      } else if (_tolerance.less(qmin,pmin)) {
+	min = qmin;
+      } else {
+	min = pmin;
+      }
+
+      if (p != -1){
+	_data[p].parent = v;
+	_data[p].dmin -= min;
+      }
+      if (q!=-1){
+	_data[q].parent = v;
+	if (_tolerance.less(_data[q].dmin,_max_value)) {
+	  _data[q].dmin -= min;
+	}
+      }
+      _data[v].left = p;
+      _data[v].right = q;
+      if (_tolerance.less(min,_max_value)) {
+	_data[v].dcost = _data[v].dmin - min;
+      }
+      _data[v].dmin = min;
+    }
+
+    void split(int &p, int v, int &q){
+      splay(v);
+      p = -1;
+      if (_data[v].left != -1){
+	p = _data[v].left;
+	_data[p].dmin += _data[v].dmin;
+	_data[p].parent = -1;
+	_data[v].left = -1;
+      }
+      q = -1;
+      if (_data[v].right != -1) {
+	q=_data[v].right;
+	if (_tolerance.less(_data[q].dmin, _max_value)) {
+	  _data[q].dmin += _data[v].dmin;
+	}
+	_data[q].parent = -1;
+	_data[v].right = -1;
+      } 
+      if (_tolerance.less(_data[v].dcost, _max_value)) {
+	_data[v].dmin += _data[v].dcost;
+	_data[v].dcost = 0;
+      } else {
+	_data[v].dmin = _data[v].dcost;
+      }
+    }
+ 
+    int expose(int v) {
+      int p, q, r, w;
+      p = -1;
+      while (v != -1) {
+	w = _data[findPath(v)].successor;
+	split(q, v, r);
+	if (q != -1) {
+	  _data[q].successor = v;
+	}
+	join(p, v, r);
+	p = v;
+	v = w;
+      }
+      _data[p].successor = -1;
+      return p;
+    }
+
+    void splay(int v) {
+      while (_data[v].parent != -1) {
+	if (v == _data[_data[v].parent].left) {
+	  if (_data[_data[v].parent].parent == -1) {
+	    zig(v);
+	  } else {
+	    if (_data[v].parent == _data[_data[_data[v].parent].parent].left) {
+	      zig(_data[v].parent);
+	      zig(v);
+	    } else {
+	      zig(v);
+	      zag(v);
+	    }
+	  }
+	} else {
+	  if (_data[_data[v].parent].parent == -1) {
+	    zag(v);
+	  } else {
+	    if (_data[v].parent == _data[_data[_data[v].parent].parent].left) {
+	      zag(v);
+	      zig(v);
+	    } else {
+	      zag(_data[v].parent);
+	      zag(v);
+	    }
+	  }
+	}
+      }
+    }
+
+
+    void zig(int v) {
+      Value min = _data[_data[v].parent].dmin;
+      int a = _data[v].parent;
+        
+      Value aa = _data[a].dcost;
+      if (_tolerance.less(aa, _max_value)) { 
+	aa+= min;
+      }
+
+
+      int b = v;
+      Value ab = min + _data[b].dmin;
+      Value bb = _data[b].dcost;
+      if (_tolerance.less(bb, _max_value)) { 
+	bb+= ab;
+      }
+
+      int c = -1;
+      Value cc = _max_value;
+      if (_data[a].right != -1) {
+	c = _data[a].right;
+	cc = _data[c].dmin;
+	if (_tolerance.less(cc, _max_value)) {
+	  cc+=min;
+	}
+      }
+
+      int d = -1;
+      Value dd = _max_value;
+      if (_data[v].left != -1){
+	d = _data[v].left;
+	dd = ab + _data[d].dmin;
+      }
+
+      int e = -1;
+      Value ee = _max_value;
+      if (_data[v].right != -1) {
+	e = _data[v].right;
+	ee = ab + _data[e].dmin;
+      }
+
+      Value min2;
+      if (_tolerance.less(0, _data[b].dmin) || 
+	  (e != -1 && !_tolerance.less(0, _data[e].dmin))) {
+	min2 = min;
+      } else {
+	if (_tolerance.less(aa, cc)) {
+	  if (_tolerance.less(aa, ee)) {
+	    min2 = aa;
+	  } else {
+	    min2 = ee;
+	  }
+	} else if (_tolerance.less(cc, ee)) {
+	  min2 = cc;
+	} else {
+	  min2 = ee;
+	}
+      }
+        
+      _data[a].dcost = aa;
+      if (_tolerance.less(aa, _max_value)) { 
+	_data[a].dcost -= min2;
+      }
+      _data[a].dmin = min2;
+      if (_tolerance.less(min2,_max_value)) { 
+	_data[a].dmin -= min; 
+      }
+      _data[a].size -= _data[b].size;
+      _data[b].dcost = bb;
+      if (_tolerance.less(bb, _max_value)) { 
+	_data[b].dcost -= min;
+      }
+      _data[b].dmin = min;
+      _data[b].size += _data[a].size;
+      if (c != -1) {
+	_data[c].dmin = cc;
+	if (_tolerance.less(cc, _max_value)) {
+	  _data[c].dmin -= min2;
+	}
+      }
+      if (d != -1) {
+	_data[d].dmin = dd - min;
+	_data[a].size += _data[d].size;
+	_data[b].size -= _data[d].size;
+      }
+      if (e != -1) {
+	_data[e].dmin = ee - min2;
+      }
+        
+      int w = _data[v].parent;
+      _data[v].successor = _data[w].successor;
+      _data[w].successor = -1;
+      _data[v].parent = _data[w].parent;
+      _data[w].parent = v;
+      _data[w].left = _data[v].right;
+      _data[v].right = w;
+      if (_data[v].parent != -1){
+	if (_data[_data[v].parent].right == w) {
+	  _data[_data[v].parent].right = v;
+	} else {
+	  _data[_data[v].parent].left = v;
+	}
+      }
+      if (_data[w].left != -1){
+	_data[_data[w].left].parent = w;
+      }
+    }
+
+
+    void zag(int v) {
+
+      Value min = _data[_data[v].parent].dmin;
+
+      int a = _data[v].parent;
+      Value aa = _data[a].dcost;
+      if (_tolerance.less(aa, _max_value)) { 
+	aa += min;
+      }
+        
+      int b = v;
+      Value ab = min + _data[b].dmin;
+      Value bb = _data[b].dcost;
+      if (_tolerance.less(bb, _max_value)) {
+	bb += ab;
+      }
+
+      int c = -1;
+      Value cc = _max_value;
+      if (_data[a].left != -1){
+	c = _data[a].left;
+	cc = min + _data[c].dmin;
+      }
+
+      int d = -1;
+      Value dd = _max_value;
+      if (_data[v].right!=-1) {
+	d = _data[v].right;
+	dd = _data[d].dmin;
+	if (_tolerance.less(dd, _max_value)) {
+	  dd += ab;
+	}
+      }
+
+      int e = -1;
+      Value ee = _max_value;
+      if (_data[v].left != -1){
+	e = _data[v].left;
+	ee = ab + _data[e].dmin;
+      }
+
+      Value min2;
+      if (_tolerance.less(0, _data[b].dmin) || 
+	  (e != -1 && !_tolerance.less(0, _data[e].dmin))) {
+	min2 = min;
+      } else {
+	if (_tolerance.less(aa, cc)) {
+	  if (_tolerance.less(aa, ee)) {
+	    min2 = aa;
+	  } else {
+	    min2 = ee;
+	  }
+	} else if (_tolerance.less(cc, ee)) {
+	  min2 = cc;
+	} else {
+	  min2 = ee;
+	}
+      }
+      _data[a].dcost = aa;
+      if (_tolerance.less(aa, _max_value)) { 
+	_data[a].dcost -= min2;
+      }
+      _data[a].dmin = min2;
+      if (_tolerance.less(min2, _max_value)) {
+	_data[a].dmin -= min;
+      }
+      _data[a].size -= _data[b].size;
+      _data[b].dcost = bb;
+      if (_tolerance.less(bb, _max_value)) { 
+	_data[b].dcost -= min;
+      }
+      _data[b].dmin = min;
+      _data[b].size += _data[a].size;
+      if (c != -1) {
+	_data[c].dmin = cc - min2;
+      }
+      if (d != -1) {
+	_data[d].dmin = dd;
+	_data[a].size += _data[d].size;
+	_data[b].size -= _data[d].size;
+	if (_tolerance.less(dd, _max_value)) {
+	  _data[d].dmin -= min;
+	}
+      }
+      if (e != -1) {
+	_data[e].dmin = ee - min2;
+      }
+        
+      int w = _data[v].parent;
+      _data[v].successor = _data[w].successor;
+      _data[w].successor = -1;
+      _data[v].parent = _data[w].parent;
+      _data[w].parent = v;
+      _data[w].right = _data[v].left;
+      _data[v].left = w;
+      if (_data[v].parent != -1){
+	if (_data[_data[v].parent].left == w) {
+	  _data[_data[v].parent].left = v;
+	} else {
+	  _data[_data[v].parent].right = v;
+	}
+      }
+      if (_data[w].right != -1){
+	_data[_data[w].right].parent = w;
+      }
+    }
+
+  private:
+
+    class ItemData {
+    public:
+      Item id;
+      int size;
+      int successor;
+      int parent;
+      int left;
+      int right;
+      Value dmin;
+      Value dcost;
+        
+    public:
+      ItemData(const Item &item)
+	: id(item), size(1), successor(), parent(-1), 
+	  left(-1), right(-1), dmin(0), dcost(0) {}
+    };
+     
+  };
+
+  template <typename _Value, typename _ItemIntMap, typename _Tolerance>
+  class DynamicTree<_Value, _ItemIntMap, _Tolerance, false> {
+  public:
+    typedef _ItemIntMap ItemIntMap;
+    typedef typename ItemIntMap::Key Item;
+    typedef _Value Value;
+    typedef _Tolerance Tolerance;
+
+  private:
+  
+    class ItemData;
+
+    std::vector<ItemData> _data;
+    ItemIntMap &_iim;
+    Value _max_value;
+    Tolerance _tolerance;
+
+  public:
+    DynamicTree(ItemIntMap &iim, const Tolerance& tolerance = Tolerance())
+      : _iim(iim), _max_value(std::numeric_limits<Value>::max()), 
+	_tolerance(tolerance) {}
+  
+    ~DynamicTree() {}
+
+    void clear() {
+      _data.clear();
+    }
+
+    void tolerance(const Tolerance& tolerance) const {
+      _tolerance = tolerance;
+      return *this;
+    } 
+  
+    const Tolerance& tolerance() const {
+      return tolerance;
+    } 
+  
+    void makeTree(const Item &item) {
+      _data[makePath(item)].successor = -1;
+    }
+    
+    Item findRoot(const Item &item) {
+      return _data[findTail(expose(_iim[item]))].id;
+    }
+    
+    Item findCost(const Item &item, Value& d){
+      return _data[findPathCost(expose(_iim[item]),d)].id;
+    }
+    
+    void addCost(const Item &item, Value x){
+      addPathCost(expose(_iim[item]), x);
+    }
+    
+    void link(const Item &item1, const Item &item2){
+      int v = _iim[item1];
+      int w = _iim[item2];
+      int p = expose(w);
+      join(-1, expose(v), p);
+      _data[v].successor = -1;
+    }    
+    
+    void cut(const Item &item) {
+      int v = _iim[item];
+      int p, q;
+      expose(v);
+      split(p, v, q);
+      if (p != -1) {
+	_data[p].successor = v;
+      }
+      if (q != -1) {
+	_data[q].successor = _data[v].successor;
+      }
+      _data[v].successor = -1;
+    }
+
+    Item parent(const Item &item){
+      return _data[_iim[item].p].id;
+    }
+
+    Value maxValue() const {
+      return _max_value;
+    }
+    
+  private:
+
+    int makePath(const Item &item) {
+      _iim.set(item, _data.size());
+      ItemData v(item);
+      _data.push_back(v);
+      return _iim[item];
+    }
+
+    int findPath(int v){
+      splay(v);
+      return v;
+    }
+    
+    int findPathCost(int p, Value &d){
+      while ((_data[p].right != -1 && 
+	      !_tolerance.less(0, _data[_data[p].right].dmin)) || 
+	     (_data[p].left != -1 && _tolerance.less(0, _data[p].dcost))) {
+	if (_data[p].right != -1 && 
+	    !_tolerance.less(0, _data[_data[p].right].dmin)) {
+	  p = _data[p].right;
+	} else if (_data[p].left != -1 && 
+		   !_tolerance.less(0, _data[_data[p].left].dmin)){
+	  p = _data[p].left;
+	}
+      }
+      splay(p);
+      d = _data[p].dmin;
+      return p; 
+    }
+
+    int findTail(int p){
+      while (_data[p].right != -1) {
+	p = _data[p].right;
+      }
+      splay(p);
+      return p;
+    }
+    
+    void addPathCost(int p, Value x){
+      if (!_tolerance.less(x, _max_value)) {
+	_data[p].dmin = x;_data[p].dcost = x;
+      } else if (!_tolerance.less(-x, _max_value)) {
+	_data[p].dmin = 0;
+	_data[p].dcost = 0;
+      } else {
+	_data[p].dmin += x;
+      }
+    }
+
+    void join(int p, int v, int q) {
+      Value min = _max_value;
+      Value pmin = _max_value;
+      Value vmin = _data[v].dmin;
+      Value qmin = _max_value;
+      if (p != -1){
+	pmin = _data[p].dmin;
+      }
+      if (q != -1){
+	qmin = _data[q].dmin;
+      }
+        
+      if (_tolerance.less(vmin, qmin)) {
+	if (_tolerance.less(vmin,pmin)) {
+	  min = vmin;
+	} else {
+	  min = pmin;
+	}
+      } else if (_tolerance.less(qmin,pmin)) {
+	min = qmin;
+      } else {
+	min = pmin;
+      }
+
+      if (p != -1){
+	_data[p].parent = v;
+	_data[p].dmin -= min;
+      }
+      if (q!=-1){
+	_data[q].parent = v;
+	if (_tolerance.less(_data[q].dmin,_max_value)) {
+	  _data[q].dmin -= min;
+	}
+      }
+      _data[v].left = p;
+      _data[v].right = q;
+      if (_tolerance.less(min,_max_value)) {
+	_data[v].dcost = _data[v].dmin - min;
+      }
+      _data[v].dmin = min;
+    }
+
+    void split(int &p, int v, int &q){
+      splay(v);
+      p = -1;
+      if (_data[v].left != -1){
+	p = _data[v].left;
+	_data[p].dmin += _data[v].dmin;
+	_data[p].parent = -1;
+	_data[v].left = -1;
+      }
+      q = -1;
+      if (_data[v].right != -1) {
+	q=_data[v].right;
+	if (_tolerance.less(_data[q].dmin, _max_value)) {
+	  _data[q].dmin += _data[v].dmin;
+	}
+	_data[q].parent = -1;
+	_data[v].right = -1;
+      } 
+      if (_tolerance.less(_data[v].dcost, _max_value)) {
+	_data[v].dmin += _data[v].dcost;
+	_data[v].dcost = 0;
+      } else {
+	_data[v].dmin = _data[v].dcost;
+      }
+    }
+ 
+    int expose(int v) {
+      int p, q, r, w;
+      p = -1;
+      while (v != -1) {
+	w = _data[findPath(v)].successor;
+	split(q, v, r);
+	if (q != -1) {
+	  _data[q].successor = v;
+	}
+	join(p, v, r);
+	p = v;
+	v = w;
+      }
+      _data[p].successor = -1;
+      return p;
+    }
+
+    void splay(int v) {
+      while (_data[v].parent != -1) {
+	if (v == _data[_data[v].parent].left) {
+	  if (_data[_data[v].parent].parent == -1) {
+	    zig(v);
+	  } else {
+	    if (_data[v].parent == _data[_data[_data[v].parent].parent].left) {
+	      zig(_data[v].parent);
+	      zig(v);
+	    } else {
+	      zig(v);
+	      zag(v);
+	    }
+	  }
+	} else {
+	  if (_data[_data[v].parent].parent == -1) {
+	    zag(v);
+	  } else {
+	    if (_data[v].parent == _data[_data[_data[v].parent].parent].left) {
+	      zag(v);
+	      zig(v);
+	    } else {
+	      zag(_data[v].parent);
+	      zag(v);
+	    }
+	  }
+	}
+      }
+    }
+
+
+    void zig(int v) {
+      Value min = _data[_data[v].parent].dmin;
+      int a = _data[v].parent;
+        
+      Value aa = _data[a].dcost;
+      if (_tolerance.less(aa, _max_value)) { 
+	aa+= min;
+      }
+
+
+      int b = v;
+      Value ab = min + _data[b].dmin;
+      Value bb = _data[b].dcost;
+      if (_tolerance.less(bb, _max_value)) { 
+	bb+= ab;
+      }
+
+      int c = -1;
+      Value cc = _max_value;
+      if (_data[a].right != -1) {
+	c = _data[a].right;
+	cc = _data[c].dmin;
+	if (_tolerance.less(cc, _max_value)) {
+	  cc+=min;
+	}
+      }
+
+      int d = -1;
+      Value dd = _max_value;
+      if (_data[v].left != -1){
+	d = _data[v].left;
+	dd = ab + _data[d].dmin;
+      }
+
+      int e = -1;
+      Value ee = _max_value;
+      if (_data[v].right != -1) {
+	e = _data[v].right;
+	ee = ab + _data[e].dmin;
+      }
+
+      Value min2;
+      if (_tolerance.less(0, _data[b].dmin) || 
+	  (e != -1 && !_tolerance.less(0, _data[e].dmin))) {
+	min2 = min;
+      } else {
+	if (_tolerance.less(aa, cc)) {
+	  if (_tolerance.less(aa, ee)) {
+	    min2 = aa;
+	  } else {
+	    min2 = ee;
+	  }
+	} else if (_tolerance.less(cc, ee)) {
+	  min2 = cc;
+	} else {
+	  min2 = ee;
+	}
+      }
+        
+      _data[a].dcost = aa;
+      if (_tolerance.less(aa, _max_value)) { 
+	_data[a].dcost -= min2;
+      }
+      _data[a].dmin = min2;
+      if (_tolerance.less(min2,_max_value)) { 
+	_data[a].dmin -= min; 
+      }
+      _data[b].dcost = bb;
+      if (_tolerance.less(bb, _max_value)) { 
+	_data[b].dcost -= min;
+      }
+      _data[b].dmin = min;
+      if (c != -1) {
+	_data[c].dmin = cc;
+	if (_tolerance.less(cc, _max_value)) {
+	  _data[c].dmin -= min2;
+	}
+      }
+      if (d != -1) {
+	_data[d].dmin = dd - min;
+      }
+      if (e != -1) {
+	_data[e].dmin = ee - min2;
+      }
+        
+      int w = _data[v].parent;
+      _data[v].successor = _data[w].successor;
+      _data[w].successor = -1;
+      _data[v].parent = _data[w].parent;
+      _data[w].parent = v;
+      _data[w].left = _data[v].right;
+      _data[v].right = w;
+      if (_data[v].parent != -1){
+	if (_data[_data[v].parent].right == w) {
+	  _data[_data[v].parent].right = v;
+	} else {
+	  _data[_data[v].parent].left = v;
+	}
+      }
+      if (_data[w].left != -1){
+	_data[_data[w].left].parent = w;
+      }
+    }
+
+
+    void zag(int v) {
+
+      Value min = _data[_data[v].parent].dmin;
+
+      int a = _data[v].parent;
+      Value aa = _data[a].dcost;
+      if (_tolerance.less(aa, _max_value)) { 
+	aa += min;
+      }
+        
+      int b = v;
+      Value ab = min + _data[b].dmin;
+      Value bb = _data[b].dcost;
+      if (_tolerance.less(bb, _max_value)) {
+	bb += ab;
+      }
+
+      int c = -1;
+      Value cc = _max_value;
+      if (_data[a].left != -1){
+	c = _data[a].left;
+	cc = min + _data[c].dmin;
+      }
+
+      int d = -1;
+      Value dd = _max_value;
+      if (_data[v].right!=-1) {
+	d = _data[v].right;
+	dd = _data[d].dmin;
+	if (_tolerance.less(dd, _max_value)) {
+	  dd += ab;
+	}
+      }
+
+      int e = -1;
+      Value ee = _max_value;
+      if (_data[v].left != -1){
+	e = _data[v].left;
+	ee = ab + _data[e].dmin;
+      }
+
+      Value min2;
+      if (_tolerance.less(0, _data[b].dmin) || 
+	  (e != -1 && !_tolerance.less(0, _data[e].dmin))) {
+	min2 = min;
+      } else {
+	if (_tolerance.less(aa, cc)) {
+	  if (_tolerance.less(aa, ee)) {
+	    min2 = aa;
+	  } else {
+	    min2 = ee;
+	  }
+	} else if (_tolerance.less(cc, ee)) {
+	  min2 = cc;
+	} else {
+	  min2 = ee;
+	}
+      }
+      _data[a].dcost = aa;
+      if (_tolerance.less(aa, _max_value)) { 
+	_data[a].dcost -= min2;
+      }
+      _data[a].dmin = min2;
+      if (_tolerance.less(min2, _max_value)) {
+	_data[a].dmin -= min;
+      }
+      _data[b].dcost = bb;
+      if (_tolerance.less(bb, _max_value)) { 
+	_data[b].dcost -= min;
+      }
+      _data[b].dmin = min;
+      if (c != -1) {
+	_data[c].dmin = cc - min2;
+      }
+      if (d != -1) {
+	_data[d].dmin = dd;
+	if (_tolerance.less(dd, _max_value)) {
+	  _data[d].dmin -= min;
+	}
+      }
+      if (e != -1) {
+	_data[e].dmin = ee - min2;
+      }
+        
+      int w = _data[v].parent;
+      _data[v].successor = _data[w].successor;
+      _data[w].successor = -1;
+      _data[v].parent = _data[w].parent;
+      _data[w].parent = v;
+      _data[w].right = _data[v].left;
+      _data[v].left = w;
+      if (_data[v].parent != -1){
+	if (_data[_data[v].parent].left == w) {
+	  _data[_data[v].parent].left = v;
+	} else {
+	  _data[_data[v].parent].right = v;
+	}
+      }
+      if (_data[w].right != -1){
+	_data[_data[w].right].parent = w;
+      }
+    }
+
+  private:
+
+    class ItemData {
+    public:
+      Item id;
+      int successor;
+      int parent;
+      int left;
+      int right;
+      Value dmin;
+      Value dcost;
+        
+    public:
+      ItemData(const Item &item)
+	: id(item), successor(), parent(-1), 
+	  left(-1), right(-1), dmin(0), dcost(0) {}
+    };
+     
+  };
+
+}
+
+#endif
diff -r 26983135fd6d -r 57143c09dc20 lemon/edmonds_karp.h
--- a/lemon/edmonds_karp.h	Wed Nov 14 17:53:08 2007 +0000
+++ b/lemon/edmonds_karp.h	Sat Nov 17 20:58:11 2007 +0000
@@ -23,47 +23,95 @@
 /// \ingroup max_flow
 /// \brief Implementation of the Edmonds-Karp algorithm.
 
-#include <lemon/graph_adaptor.h>
 #include <lemon/tolerance.h>
-#include <lemon/bfs.h>
+#include <vector>
 
 namespace lemon {
 
+  /// \brief Default traits class of EdmondsKarp class.
+  ///
+  /// Default traits class of EdmondsKarp class.
+  /// \param _Graph Graph type.
+  /// \param _CapacityMap Type of capacity map.
+  template <typename _Graph, typename _CapacityMap>
+  struct EdmondsKarpDefaultTraits {
+
+    /// \brief The graph type the algorithm runs on. 
+    typedef _Graph Graph;
+
+    /// \brief The type of the map that stores the edge capacities.
+    ///
+    /// The type of the map that stores the edge capacities.
+    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
+    typedef _CapacityMap CapacityMap;
+
+    /// \brief The type of the length of the edges.
+    typedef typename CapacityMap::Value Value;
+
+    /// \brief The map type that stores the flow values.
+    ///
+    /// The map type that stores the flow values. 
+    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
+    typedef typename Graph::template EdgeMap<Value> FlowMap;
+
+    /// \brief Instantiates a FlowMap.
+    ///
+    /// This function instantiates a \ref FlowMap. 
+    /// \param graph The graph, to which we would like to define the flow map.
+    static FlowMap* createFlowMap(const Graph& graph) {
+      return new FlowMap(graph);
+    }
+
+    /// \brief The tolerance used by the algorithm
+    ///
+    /// The tolerance used by the algorithm to handle inexact computation.
+    typedef Tolerance<Value> Tolerance;
+
+  };
+
   /// \ingroup max_flow
+  ///
   /// \brief Edmonds-Karp algorithms class.
   ///
   /// This class provides an implementation of the \e Edmonds-Karp \e
   /// algorithm producing a flow of maximum value in a directed
-  /// graph. The Edmonds-Karp algorithm is slower than the Preflow algorithm
-  /// but it has an advantage of the step-by-step execution control with
-  /// feasible flow solutions. The \e source node, the \e target node, the \e
-  /// capacity of the edges and the \e starting \e flow value of the
-  /// edges should be passed to the algorithm through the
-  /// constructor.
+  /// graphs. The Edmonds-Karp algorithm is slower than the Preflow
+  /// algorithm but it has an advantage of the step-by-step execution
+  /// control with feasible flow solutions. The \e source node, the \e
+  /// target node, the \e capacity of the edges and the \e starting \e
+  /// flow value of the edges should be passed to the algorithm
+  /// through the constructor.
   ///
-  /// The time complexity of the algorithm is \f$ O(n * e^2) \f$ in
+  /// The time complexity of the algorithm is \f$ O(nm^2) \f$ in
   /// worst case.  Always try the preflow algorithm instead of this if
   /// you just want to compute the optimal flow.
   ///
   /// \param _Graph The directed graph type the algorithm runs on.
-  /// \param _Number The number type of the capacities and the flow values.
   /// \param _CapacityMap The capacity map type.
-  /// \param _FlowMap The flow map type.
-  /// \param _Tolerance The tolerance class to handle computation problems.
+  /// \param _Traits Traits class to set various data types used by
+  /// the algorithm.  The default traits class is \ref
+  /// EdmondsKarpDefaultTraits.  See \ref EdmondsKarpDefaultTraits for the
+  /// documentation of a Edmonds-Karp traits class. 
   ///
   /// \author Balazs Dezso 
 #ifdef DOXYGEN
-  template <typename _Graph, typename _Number,
-	    typename _CapacityMap, typename _FlowMap, typename _Tolerance>
-#else
-  template <typename _Graph, typename _Number,
-	    typename _CapacityMap = typename _Graph::template EdgeMap<_Number>,
-            typename _FlowMap = typename _Graph::template EdgeMap<_Number>,
-	    typename _Tolerance = Tolerance<_Number> >
+  template <typename _Graph, typename _CapacityMap, typename _Traits>
+#else 
+  template <typename _Graph,
+	    typename _CapacityMap = typename _Graph::template EdgeMap<int>,
+            typename _Traits = EdmondsKarpDefaultTraits<_Graph, _CapacityMap> >
 #endif
   class EdmondsKarp {
   public:
 
+    typedef _Traits Traits;
+    typedef typename Traits::Graph Graph;
+    typedef typename Traits::CapacityMap CapacityMap;
+    typedef typename Traits::Value Value; 
+
+    typedef typename Traits::FlowMap FlowMap;
+    typedef typename Traits::Tolerance Tolerance;
+
     /// \brief \ref Exception for the case when the source equals the target.
     ///
     /// \ref Exception for the case when the source equals the target.
@@ -76,110 +124,239 @@
     };
 
 
-    /// \brief The graph type the algorithm runs on. 
-    typedef _Graph Graph;
-    /// \brief The value type of the algorithms.
-    typedef _Number Number;
-    /// \brief The capacity map on the edges.
-    typedef _CapacityMap CapacityMap;
-    /// \brief The flow map on the edges.
-    typedef _FlowMap FlowMap;
-    /// \brief The tolerance used by the algorithm.
-    typedef _Tolerance Tolerance;
-
-    typedef ResGraphAdaptor<Graph, Number, CapacityMap, 
-                            FlowMap, Tolerance> ResGraph;
-
   private:
 
-    typedef typename Graph::Node Node;
-    typedef typename Graph::Edge Edge;
+    GRAPH_TYPEDEFS(typename Graph);
+    typedef typename Graph::template NodeMap<Edge> PredMap;
     
-    typedef typename Graph::NodeIt NodeIt;
-    typedef typename Graph::EdgeIt EdgeIt;
-    typedef typename Graph::InEdgeIt InEdgeIt;
-    typedef typename Graph::OutEdgeIt OutEdgeIt;
+    const Graph& _graph;
+    const CapacityMap* _capacity;
+
+    Node _source, _target;
+
+    FlowMap* _flow;
+    bool _local_flow;
+
+    PredMap* _pred;
+    std::vector<Node> _queue;
+    
+    Tolerance _tolerance;
+    Value _flow_value;
+
+    void createStructures() {
+      if (!_flow) {
+	_flow = Traits::createFlowMap(_graph);
+	_local_flow = true;
+      }
+      if (!_pred) {
+	_pred = new PredMap(_graph);
+      }
+      _queue.resize(countNodes(_graph));
+    }
+
+    void destroyStructures() {
+      if (_local_flow) {
+	delete _flow;
+      }
+      if (_pred) {
+	delete _pred;
+      }
+    }
     
   public:
 
+    ///\name Named template parameters
+
+    ///@{
+
+    template <typename _FlowMap>
+    struct DefFlowMapTraits : public Traits {
+      typedef _FlowMap FlowMap;
+      static FlowMap *createFlowMap(const Graph&) {
+	throw UninitializedParameter();
+      }
+    };
+
+    /// \brief \ref named-templ-param "Named parameter" for setting
+    /// FlowMap type
+    ///
+    /// \ref named-templ-param "Named parameter" for setting FlowMap
+    /// type
+    template <typename _FlowMap>
+    struct DefFlowMap 
+      : public EdmondsKarp<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> > {
+      typedef EdmondsKarp<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> > 
+      Create;
+    };
+
+
+    /// @}
+
     /// \brief The constructor of the class.
     ///
     /// The constructor of the class. 
     /// \param graph The directed graph the algorithm runs on. 
+    /// \param capacity The capacity of the edges. 
     /// \param source The source node.
     /// \param target The target node.
-    /// \param capacity The capacity of the edges. 
-    /// \param flow The flow of the edges. 
-    /// \param tolerance Tolerance class.
-    EdmondsKarp(const Graph& graph, Node source, Node target,
-                const CapacityMap& capacity, FlowMap& flow, 
-                const Tolerance& tolerance = Tolerance())
-      : _graph(graph), _capacity(capacity), _flow(flow), 
-        _tolerance(tolerance), _resgraph(graph, capacity, flow, tolerance), 
-        _source(source), _target(target) 
+    EdmondsKarp(const Graph& graph, const CapacityMap& capacity,
+		Node source, Node target)
+      : _graph(graph), _capacity(&capacity), _source(source), _target(target),
+	_flow(0), _local_flow(false), _pred(0), _tolerance(), _flow_value()
     {
       if (_source == _target) {
         throw InvalidArgument();
       }
     }
 
+    /// \brief Destrcutor.
+    ///
+    /// Destructor.
+    ~EdmondsKarp() {
+      destroyStructures();
+    }
+
+    /// \brief Sets the capacity map.
+    ///
+    /// Sets the capacity map.
+    /// \return \c (*this)
+    EdmondsKarp& capacityMap(const CapacityMap& map) {
+      _capacity = &map;
+      return *this;
+    }
+
+    /// \brief Sets the flow map.
+    ///
+    /// Sets the flow map.
+    /// \return \c (*this)
+    EdmondsKarp& flowMap(FlowMap& map) {
+      if (_local_flow) {
+	delete _flow;
+	_local_flow = false;
+      }
+      _flow = &map;
+      return *this;
+    }
+
+    /// \brief Returns the flow map.
+    ///
+    /// \return The flow map.
+    const FlowMap& flowMap() {
+      return *_flow;
+    }
+
+    /// \brief Sets the source node.
+    ///
+    /// Sets the source node.
+    /// \return \c (*this)
+    EdmondsKarp& source(const Node& node) {
+      _source = node;
+      return *this;
+    }
+
+    /// \brief Sets the target node.
+    ///
+    /// Sets the target node.
+    /// \return \c (*this)
+    EdmondsKarp& target(const Node& node) {
+      _target = node;
+      return *this;
+    }
+
+    /// \brief Sets the tolerance used by algorithm.
+    ///
+    /// Sets the tolerance used by algorithm.
+    EdmondsKarp& tolerance(const Tolerance& tolerance) const {
+      _tolerance = tolerance;
+      return *this;
+    } 
+
+    /// \brief Returns the tolerance used by algorithm.
+    ///
+    /// Returns the tolerance used by algorithm.
+    const Tolerance& tolerance() const {
+      return tolerance;
+    } 
+
+    /// \name Execution control The simplest way to execute the
+    /// algorithm is to use the \c run() member functions.
+    /// \n
+    /// If you need more control on initial solution or
+    /// execution then you have to call one \ref init() function and then
+    /// the start() or multiple times the \c augment() member function.  
+    
+    ///@{
+
     /// \brief Initializes the algorithm
     /// 
     /// It sets the flow to empty flow.
     void init() {
+      createStructures();
       for (EdgeIt it(_graph); it != INVALID; ++it) {
-        _flow.set(it, 0);
+        _flow->set(it, 0);
       }
-      _value = 0;
+      _flow_value = 0;
     }
     
     /// \brief Initializes the algorithm
     /// 
-    /// If the flow map initially flow this let the flow map
-    /// unchanged but the flow value will be set by the flow
-    /// on the outedges from the source.
-    void flowInit() {
-      _value = 0;
+    /// Initializes the flow to the \c flowMap. The \c flowMap should
+    /// contain a feasible flow, ie. in each node excluding the source
+    /// and the target the incoming flow should be equal to the
+    /// outgoing flow.
+    template <typename FlowMap>
+    void flowInit(const FlowMap& flowMap) {
+      createStructures();
+      for (EdgeIt e(_graph); e != INVALID; ++e) {
+	_flow->set(e, flowMap[e]);
+      }
+      _flow_value = 0;
       for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
-        _value += _flow[jt];
+        _flow_value += (*_flow)[jt];
       }
       for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
-        _value -= _flow[jt];
+        _flow_value -= (*_flow)[jt];
       }
     }
 
     /// \brief Initializes the algorithm
     /// 
-    /// If the flow map initially flow this let the flow map
-    /// unchanged but the flow value will be set by the flow
-    /// on the outedges from the source. It also checks that
-    /// the flow map really contains a flow.
-    /// \return %True when the flow map really a flow.
-    bool checkedFlowInit() {
-      _value = 0;
-      for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
-        _value += _flow[jt];
-      }
-      for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
-        _value -= _flow[jt];
+    /// Initializes the flow to the \c flowMap. The \c flowMap should
+    /// contain a feasible flow, ie. in each node excluding the source
+    /// and the target the incoming flow should be equal to the
+    /// outgoing flow.  
+    /// \return %False when the given flowMap does not contain
+    /// feasible flow.
+    template <typename FlowMap>
+    bool checkedFlowInit(const FlowMap& flowMap) {
+      createStructures();
+      for (EdgeIt e(_graph); e != INVALID; ++e) {
+	_flow->set(e, flowMap[e]);
       }
       for (NodeIt it(_graph); it != INVALID; ++it) {
         if (it == _source || it == _target) continue;
-        Number outFlow = 0;
+        Value outFlow = 0;
         for (OutEdgeIt jt(_graph, it); jt != INVALID; ++jt) {
-          outFlow += _flow[jt];
+          outFlow += (*_flow)[jt];
         }
-        Number inFlow = 0;
+        Value inFlow = 0;
         for (InEdgeIt jt(_graph, it); jt != INVALID; ++jt) {
-          inFlow += _flow[jt];
+          inFlow += (*_flow)[jt];
         }
         if (_tolerance.different(outFlow, inFlow)) {
           return false;
         }
       }
       for (EdgeIt it(_graph); it != INVALID; ++it) {
-        if (_tolerance.less(_flow[it], 0)) return false;
-        if (_tolerance.less(_capacity[it], _flow[it])) return false;
+        if (_tolerance.less((*_flow)[it], 0)) return false;
+        if (_tolerance.less((*_capacity)[it], (*_flow)[it])) return false;
+      }
+      _flow_value = 0;
+      for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
+        _flow_value += (*_flow)[jt];
+      }
+      for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
+        _flow_value -= (*_flow)[jt];
       }
       return true;
     }
@@ -195,32 +372,76 @@
     /// \return %False when the augmenting is not success so the
     /// current flow is a feasible and optimal solution.
     bool augment() {
-      typename Bfs<ResGraph>
-      ::template DefDistMap<NullMap<Node, int> >
-      ::Create bfs(_resgraph);
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+	_pred->set(n, INVALID);
+      }
+      
+      int first = 0, last = 1;
+      
+      _queue[0] = _source;
+      _pred->set(_source, OutEdgeIt(_graph, _source));
 
-      NullMap<Node, int> distMap;
-      bfs.distMap(distMap);
-      
-      bfs.init();
-      bfs.addSource(_source);
-      bfs.start(_target);
+      while (first != last && (*_pred)[_target] == INVALID) {
+	Node n = _queue[first++];
+	
+	for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+	  Value rem = (*_capacity)[e] - (*_flow)[e];
+	  Node t = _graph.target(e);
+	  if (_tolerance.positive(rem) && (*_pred)[t] == INVALID) {
+	    _pred->set(t, e);
+	    _queue[last++] = t;
+	  }
+	}
+	for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+	  Value rem = (*_flow)[e];
+	  Node t = _graph.source(e);
+	  if (_tolerance.positive(rem) && (*_pred)[t] == INVALID) {
+	    _pred->set(t, e);
+	    _queue[last++] = t;
+	  }
+	}
+      }
 
-      if (!bfs.reached(_target)) {
-        return false;
+      if ((*_pred)[_target] != INVALID) {
+	Node n = _target;
+	Edge e = (*_pred)[n];
+
+	Value prem = (*_capacity)[e] - (*_flow)[e];
+	n = _graph.source(e);
+	while (n != _source) {
+	  e = (*_pred)[n];
+	  if (_graph.target(e) == n) {
+	    Value rem = (*_capacity)[e] - (*_flow)[e];
+	    if (rem < prem) prem = rem;
+	    n = _graph.source(e);
+	  } else {
+	    Value rem = (*_flow)[e];
+	    if (rem < prem) prem = rem;
+	    n = _graph.target(e);   
+	  } 
+	}
+
+	n = _target;
+	e = (*_pred)[n];
+
+	_flow->set(e, (*_flow)[e] + prem);
+	n = _graph.source(e);
+	while (n != _source) {
+	  e = (*_pred)[n];
+	  if (_graph.target(e) == n) {
+	    _flow->set(e, (*_flow)[e] + prem);
+	    n = _graph.source(e);
+	  } else {
+	    _flow->set(e, (*_flow)[e] - prem);
+	    n = _graph.target(e);   
+	  } 
+	}
+
+	_flow_value += prem;	
+	return true;
+      } else {
+	return false;
       }
-      Number min = _resgraph.rescap(bfs.predEdge(_target));
-      for (Node it = bfs.predNode(_target); it != _source; 
-           it = bfs.predNode(it)) {
-        if (min > _resgraph.rescap(bfs.predEdge(it))) {
-          min = _resgraph.rescap(bfs.predEdge(it));
-        }
-      } 
-      for (Node it = _target; it != _source; it = bfs.predNode(it)) {
-        _resgraph.augment(bfs.predEdge(it), min);
-      }
-      _value += min;
-      return true;
     }
 
     /// \brief Executes the algorithm
@@ -230,13 +451,6 @@
       while (augment()) {}
     }
 
-    /// \brief Gives back the current flow value.
-    ///
-    /// Gives back the current flow _value.
-    Number flowValue() const {
-      return _value;
-    }
-
     /// \brief runs the algorithm.
     /// 
     /// It is just a shorthand for:
@@ -250,119 +464,59 @@
       start();
     }
 
+    /// @}
+
+    /// \name Query Functions
+    /// The result of the %Dijkstra algorithm can be obtained using these
+    /// functions.\n
+    /// Before the use of these functions,
+    /// either run() or start() must be called.
+    
+    ///@{
+
+    /// \brief Returns the value of the maximum flow.
+    ///
+    /// Returns the value of the maximum flow by returning the excess
+    /// of the target node \c t. This value equals to the value of
+    /// the maximum flow already after the first phase.
+    Value flowValue() const {
+      return _flow_value;
+    }
+
+
+    /// \brief Returns the flow on the edge.
+    ///
+    /// Sets the \c flowMap to the flow on the edges. This method can
+    /// be called after the second phase of algorithm.
+    Value flow(const Edge& edge) const {
+      return (*_flow)[edge];
+    }
+
+    /// \brief Returns true when the node is on the source side of minimum cut.
+    ///
+
+    /// Returns true when the node is on the source side of minimum
+    /// cut. This method can be called both after running \ref
+    /// startFirstPhase() and \ref startSecondPhase().
+    bool minCut(const Node& node) const {
+      return (*_pred)[node] != INVALID;
+    }
+
     /// \brief Returns a minimum value cut.
     ///
     /// Sets \c cut to the characteristic vector of a minimum value cut
     /// It simply calls the minMinCut member.
     /// \retval cut Write node bool map. 
     template <typename CutMap>
-    void minCut(CutMap& cut) const {
-      minMinCut(cut);
-    }
+    void minCutMap(CutMap& cutMap) const {
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+	cutMap.set(n, (*_pred)[n] != INVALID);
+      }
+      cutMap.set(_source, true);
+    }    
 
-    /// \brief Returns the inclusionwise minimum of the minimum value cuts.
-    ///
-    /// Sets \c cut to the characteristic vector of the minimum value cut
-    /// which is inclusionwise minimum. It is computed by processing a
-    /// bfs from the source node \c source in the residual graph.  
-    /// \retval cut Write node bool map. 
-    template <typename CutMap>
-    void minMinCut(CutMap& cut) const {
+    /// @}
 
-      typename Bfs<ResGraph>
-      ::template DefDistMap<NullMap<Node, int> >
-      ::template DefProcessedMap<CutMap>
-      ::Create bfs(_resgraph);
-
-      NullMap<Node, int> distMap;
-      bfs.distMap(distMap);
-
-      bfs.processedMap(cut);
-     
-      bfs.run(_source);
-    }
-
-    /// \brief Returns the inclusionwise minimum of the minimum value cuts.
-    ///
-    /// Sets \c cut to the characteristic vector of the minimum value cut
-    /// which is inclusionwise minimum. It is computed by processing a
-    /// bfs from the source node \c source in the residual graph.  
-    /// \retval cut Write node bool map. 
-    template <typename CutMap>
-    void maxMinCut(CutMap& cut) const {
-
-      typedef RevGraphAdaptor<const ResGraph> RevGraph;
-
-      RevGraph revgraph(_resgraph);
-
-      typename Bfs<RevGraph>
-      ::template DefDistMap<NullMap<Node, int> >
-      ::template DefPredMap<NullMap<Node, Edge> >
-      ::template DefProcessedMap<NotWriteMap<CutMap> >
-      ::Create bfs(revgraph);
-
-      NullMap<Node, int> distMap;
-      bfs.distMap(distMap);
-
-      NullMap<Node, Edge> predMap;
-      bfs.predMap(predMap);
-
-      NotWriteMap<CutMap> notcut(cut);
-      bfs.processedMap(notcut);
-     
-      bfs.run(_target);
-    }
-
-    /// \brief Returns the source node.
-    ///
-    /// Returns the source node.
-    /// 
-    Node source() const { 
-      return _source;
-    }
-
-    /// \brief Returns the target node.
-    ///
-    /// Returns the target node.
-    /// 
-    Node target() const { 
-      return _target;
-    }
-
-    /// \brief Returns a reference to capacity map.
-    ///
-    /// Returns a reference to capacity map.
-    /// 
-    const CapacityMap &capacityMap() const { 
-      return *_capacity;
-    }
-     
-    /// \brief Returns a reference to flow map.
-    ///
-    /// Returns a reference to flow map.
-    /// 
-    const FlowMap &flowMap() const { 
-      return *_flow;
-    }
-
-    /// \brief Returns the tolerance used by algorithm.
-    ///
-    /// Returns the tolerance used by algorithm.
-    const Tolerance& tolerance() const {
-      return tolerance;
-    } 
-    
-  private:
-    
-    const Graph& _graph;
-    const CapacityMap& _capacity;
-    FlowMap& _flow;
-    Tolerance _tolerance;
-    
-    ResGraph _resgraph;
-    Node _source, _target;
-    Number _value;
-    
   };
 
 }
diff -r 26983135fd6d -r 57143c09dc20 lemon/goldberg_tarjan.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/goldberg_tarjan.h	Sat Nov 17 20:58:11 2007 +0000
@@ -0,0 +1,1020 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2007
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_GOLDBERG_TARJAN_H
+#define LEMON_GOLDBERG_TARJAN_H
+
+#include <vector>
+#include <queue>
+
+#include <lemon/error.h>
+#include <lemon/bits/invalid.h>
+#include <lemon/tolerance.h>
+#include <lemon/maps.h>
+#include <lemon/graph_utils.h>
+#include <lemon/dynamic_tree.h>
+#include <limits>
+
+/// \file
+/// \ingroup max_flow
+/// \brief Implementation of the preflow algorithm.
+
+namespace lemon {
+
+  /// \brief Default traits class of GoldbergTarjan class.
+  ///
+  /// Default traits class of GoldbergTarjan class.
+  /// \param _Graph Graph type.
+  /// \param _CapacityMap Type of capacity map.
+  template <typename _Graph, typename _CapacityMap>
+  struct GoldbergTarjanDefaultTraits {
+
+    /// \brief The graph type the algorithm runs on. 
+    typedef _Graph Graph;
+
+    /// \brief The type of the map that stores the edge capacities.
+    ///
+    /// The type of the map that stores the edge capacities.
+    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
+    typedef _CapacityMap CapacityMap;
+
+    /// \brief The type of the length of the edges.
+    typedef typename CapacityMap::Value Value;
+
+    /// \brief The map type that stores the flow values.
+    ///
+    /// The map type that stores the flow values. 
+    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
+    typedef typename Graph::template EdgeMap<Value> FlowMap;
+
+    /// \brief Instantiates a FlowMap.
+    ///
+    /// This function instantiates a \ref FlowMap. 
+    /// \param graph The graph, to which we would like to define the flow map.
+    static FlowMap* createFlowMap(const Graph& graph) {
+      return new FlowMap(graph);
+    }
+
+    /// \brief The eleavator type used by GoldbergTarjan algorithm.
+    /// 
+    /// The elevator type used by GoldbergTarjan algorithm.
+    ///
+    /// \sa Elevator
+    /// \sa LinkedElevator
+    typedef LinkedElevator<Graph, typename Graph::Node> Elevator;
+    
+    /// \brief Instantiates an Elevator.
+    ///
+    /// This function instantiates a \ref Elevator. 
+    /// \param graph The graph, to which we would like to define the elevator.
+    /// \param max_level The maximum level of the elevator.
+    static Elevator* createElevator(const Graph& graph, int max_level) {
+      return new Elevator(graph, max_level);
+    }
+
+    /// \brief The tolerance used by the algorithm
+    ///
+    /// The tolerance used by the algorithm to handle inexact computation.
+    typedef Tolerance<Value> Tolerance;
+
+  };
+
+  /// \ingroup max_flow
+  /// \brief Goldberg-Tarjan algorithms class.
+  ///
+  /// This class provides an implementation of the \e GoldbergTarjan
+  /// \e algorithm producing a flow of maximum value in a directed
+  /// graph. The GoldbergTarjan algorithm is a theoretical improvement
+  /// of the Goldberg's \ref Preflow "preflow" algorithm by using the \ref
+  /// DynamicTree "dynamic tree" data structure of Sleator and Tarjan.
+  /// 
+  /// The original preflow algorithm with \e "highest label" or \e
+  /// FIFO heuristic has \f$O(n^3)\f$ time complexity. The current
+  /// algorithm improved this complexity to
+  /// \f$O(nm\log(\frac{n^2}{m}))\f$. The algorithm builds limited
+  /// size dynamic trees on the residual graph, which can be used to
+  /// make multilevel push operations and gives a better bound on the
+  /// number of non-saturating pushes.
+  ///
+  /// \param Graph The directed graph type the algorithm runs on.
+  /// \param CapacityMap The capacity map type.
+  /// \param _Traits Traits class to set various data types used by
+  /// the algorithm.  The default traits class is \ref
+  /// GoldbergTarjanDefaultTraits.  See \ref
+  /// GoldbergTarjanDefaultTraits for the documentation of a
+  /// Goldberg-Tarjan traits class.
+  ///
+  /// \author Tamas Hamori and Balazs Dezso
+#ifdef DOXYGEN
+  template <typename _Graph, typename _CapacityMap, typename _Traits> 
+#else
+  template <typename _Graph,
+	    typename _CapacityMap = typename _Graph::template EdgeMap<int>,
+            typename _Traits = 
+	    GoldbergTarjanDefaultTraits<_Graph, _CapacityMap> >
+#endif
+  class GoldbergTarjan {
+  public:
+
+    typedef _Traits Traits;
+    typedef typename Traits::Graph Graph;
+    typedef typename Traits::CapacityMap CapacityMap;
+    typedef typename Traits::Value Value; 
+
+    typedef typename Traits::FlowMap FlowMap;
+    typedef typename Traits::Elevator Elevator;
+    typedef typename Traits::Tolerance Tolerance;
+    
+  protected:
+
+    GRAPH_TYPEDEFS(typename Graph);
+
+    typedef typename Graph::template NodeMap<Node> NodeNodeMap;
+    typedef typename Graph::template NodeMap<int> IntNodeMap;
+
+    typedef typename Graph::template NodeMap<Edge> EdgeNodeMap;
+    typedef typename Graph::template EdgeMap<Edge> EdgeEdgeMap;
+
+    typedef typename std::vector<Node> VecNode;
+ 
+    typedef DynamicTree<Value,IntNodeMap,Tolerance> DynTree;
+
+    const Graph& _graph;
+    const CapacityMap* _capacity;
+    int _node_num; //the number of nodes of G
+
+    Node _source;
+    Node _target;
+
+    FlowMap* _flow;
+    bool _local_flow;
+
+    Elevator* _level;
+    bool _local_level;
+
+    typedef typename Graph::template NodeMap<Value> ExcessMap;
+    ExcessMap* _excess;
+    
+    Tolerance _tolerance;
+
+    bool _phase;
+
+    // constant for treesize
+    static const double _tree_bound = 2;
+    int _max_tree_size;
+    
+    //tags for the dynamic tree
+    DynTree *_dt; 
+    //datastructure of dyanamic tree
+    IntNodeMap *_dt_index;
+    //datastrucure for solution of the communication between the two class
+    EdgeNodeMap *_dt_edges; 
+    //nodeMap for storing the outgoing edge from the node in the tree  
+
+    //max of the type Value
+    const Value _max_value;
+            
+  public:
+
+    typedef GoldbergTarjan Create;
+
+    ///\name Named template parameters
+
+    ///@{
+
+    template <typename _FlowMap>
+    struct DefFlowMapTraits : public Traits {
+      typedef _FlowMap FlowMap;
+      static FlowMap *createFlowMap(const Graph&) {
+	throw UninitializedParameter();
+      }
+    };
+
+    /// \brief \ref named-templ-param "Named parameter" for setting
+    /// FlowMap type
+    ///
+    /// \ref named-templ-param "Named parameter" for setting FlowMap
+    /// type
+    template <typename _FlowMap>
+    struct DefFlowMap 
+      : public GoldbergTarjan<Graph, CapacityMap, 
+			      DefFlowMapTraits<_FlowMap> > {
+      typedef GoldbergTarjan<Graph, CapacityMap, 
+			     DefFlowMapTraits<_FlowMap> > Create;
+    };
+
+    template <typename _Elevator>
+    struct DefElevatorTraits : public Traits {
+      typedef _Elevator Elevator;
+      static Elevator *createElevator(const Graph&, int) {
+	throw UninitializedParameter();
+      }
+    };
+
+    /// \brief \ref named-templ-param "Named parameter" for setting
+    /// Elevator type
+    ///
+    /// \ref named-templ-param "Named parameter" for setting Elevator
+    /// type
+    template <typename _Elevator>
+    struct DefElevator 
+      : public GoldbergTarjan<Graph, CapacityMap, 
+			      DefElevatorTraits<_Elevator> > {
+      typedef GoldbergTarjan<Graph, CapacityMap, 
+			     DefElevatorTraits<_Elevator> > Create;
+    };
+
+    template <typename _Elevator>
+    struct DefStandardElevatorTraits : public Traits {
+      typedef _Elevator Elevator;
+      static Elevator *createElevator(const Graph& graph, int max_level) {
+	return new Elevator(graph, max_level);
+      }
+    };
+
+    /// \brief \ref named-templ-param "Named parameter" for setting
+    /// Elevator type
+    ///
+    /// \ref named-templ-param "Named parameter" for setting Elevator
+    /// type. The Elevator should be standard constructor interface, ie.
+    /// the graph and the maximum level should be passed to it.
+    template <typename _Elevator>
+    struct DefStandardElevator 
+      : public GoldbergTarjan<Graph, CapacityMap, 
+			      DefStandardElevatorTraits<_Elevator> > {
+      typedef GoldbergTarjan<Graph, CapacityMap, 
+			     DefStandardElevatorTraits<_Elevator> > Create;
+    };    
+
+
+    ///\ref Exception for the case when s=t.
+
+    ///\ref Exception for the case when the source equals the target.
+    class InvalidArgument : public lemon::LogicError {
+    public:
+      virtual const char* what() const throw() {
+	return "lemon::GoldbergTarjan::InvalidArgument";
+      }
+    };
+
+  public: 
+
+    /// \brief The constructor of the class.
+    ///
+    /// The constructor of the class. 
+    /// \param graph The directed graph the algorithm runs on. 
+    /// \param capacity The capacity of the edges. 
+    /// \param source The source node.
+    /// \param target The target node.
+    /// Except the graph, all of these parameters can be reset by
+    /// calling \ref source, \ref target, \ref capacityMap and \ref
+    /// flowMap, resp.
+    GoldbergTarjan(const Graph& graph, const CapacityMap& capacity,
+		   Node source, Node target) 
+      : _graph(graph), _capacity(&capacity), 
+	_node_num(), _source(source), _target(target),
+	_flow(0), _local_flow(false),
+	_level(0), _local_level(false),
+	_excess(0), _tolerance(), 
+	_phase(), _max_tree_size(),
+	_dt(0), _dt_index(0), _dt_edges(0), 
+	_max_value(std::numeric_limits<Value>::max()) { 
+      if (_source == _target) throw InvalidArgument();
+    }
+
+    /// \brief Destrcutor.
+    ///
+    /// Destructor.
+    ~GoldbergTarjan() {
+      destroyStructures();
+    }
+    
+    /// \brief Sets the capacity map.
+    ///
+    /// Sets the capacity map.
+    /// \return \c (*this)
+    GoldbergTarjan& capacityMap(const CapacityMap& map) {
+      _capacity = &map;
+      return *this;
+    }
+
+    /// \brief Sets the flow map.
+    ///
+    /// Sets the flow map.
+    /// \return \c (*this)
+    GoldbergTarjan& flowMap(FlowMap& map) {
+      if (_local_flow) {
+	delete _flow;
+	_local_flow = false;
+      }
+      _flow = &map;
+      return *this;
+    }
+
+    /// \brief Returns the flow map.
+    ///
+    /// \return The flow map.
+    const FlowMap& flowMap() {
+      return *_flow;
+    }
+
+    /// \brief Sets the elevator.
+    ///
+    /// Sets the elevator.
+    /// \return \c (*this)
+    GoldbergTarjan& elevator(Elevator& elevator) {
+      if (_local_level) {
+	delete _level;
+	_local_level = false;
+      }
+      _level = &elevator;
+      return *this;
+    }
+
+    /// \brief Returns the elevator.
+    ///
+    /// \return The elevator.
+    const Elevator& elevator() {
+      return *_level;
+    }
+
+    /// \brief Sets the source node.
+    ///
+    /// Sets the source node.
+    /// \return \c (*this)
+    GoldbergTarjan& source(const Node& node) {
+      _source = node;
+      return *this;
+    }
+
+    /// \brief Sets the target node.
+    ///
+    /// Sets the target node.
+    /// \return \c (*this)
+    GoldbergTarjan& target(const Node& node) {
+      _target = node;
+      return *this;
+    }
+ 
+    /// \brief Sets the tolerance used by algorithm.
+    ///
+    /// Sets the tolerance used by algorithm.
+    GoldbergTarjan& tolerance(const Tolerance& tolerance) const {
+      _tolerance = tolerance;
+      if (_dt) {
+	_dt->tolerance(_tolerance);
+      }
+      return *this;
+    } 
+
+    /// \brief Returns the tolerance used by algorithm.
+    ///
+    /// Returns the tolerance used by algorithm.
+    const Tolerance& tolerance() const {
+      return tolerance;
+    } 
+
+         
+  private:
+    
+    void createStructures() {
+      _node_num = countNodes(_graph);
+
+      _max_tree_size = (double(_node_num) * double(_node_num)) / 
+	double(countEdges(_graph));
+
+      if (!_flow) {
+	_flow = Traits::createFlowMap(_graph);
+	_local_flow = true;
+      }
+      if (!_level) {
+	_level = Traits::createElevator(_graph, _node_num);
+	_local_level = true;
+      }
+      if (!_excess) {
+	_excess = new ExcessMap(_graph);
+      }
+      if (!_dt_index && !_dt) {
+	_dt_index = new IntNodeMap(_graph);
+	_dt = new DynTree(*_dt_index, _tolerance);
+      }
+      if (!_dt_edges) {
+	_dt_edges = new EdgeNodeMap(_graph);
+      }
+    }
+
+    void destroyStructures() {
+      if (_local_flow) {
+	delete _flow;
+      }
+      if (_local_level) {
+	delete _level;
+      }
+      if (_excess) {
+	delete _excess;
+      }
+      if (_dt) {
+	delete _dt;
+      }
+      if (_dt_index) {
+	delete _dt_index;
+      }
+      if (_dt_edges) {
+	delete _dt_edges;
+      }
+    }
+
+    bool sendOut(Node n, Value& excess) {
+
+      Node w = _dt->findRoot(n);
+      
+      while (w != n) {
+      
+	Value rem;      
+	Node u = _dt->findCost(n, rem);
+
+        if (_tolerance.positive(rem) && !_level->active(w) && w != _target) {
+	  _level->activate(w);
+        }
+
+        if (!_tolerance.less(rem, excess)) {
+	  _dt->addCost(n, - excess);
+	  _excess->set(w, (*_excess)[w] + excess);
+	  excess = 0;
+	  return true;
+	}
+	
+	
+        _dt->addCost(n, - rem);
+
+        excess -= rem;
+        _excess->set(w, (*_excess)[w] + rem);
+        
+	_dt->cut(u);
+	_dt->addCost(u, _max_value);
+	  
+	Edge e = (*_dt_edges)[u];
+	_dt_edges->set(u, INVALID);
+        
+	if (u == _graph.source(e)) {
+	  _flow->set(e, (*_capacity)[e]);
+	} else {
+	  _flow->set(e, 0);
+	}           
+
+        w = _dt->findRoot(n);
+      }      
+      return false;
+    }
+
+    bool sendIn(Node n, Value& excess) {
+
+      Node w = _dt->findRoot(n);
+      
+      while (w != n) {
+      
+	Value rem;      
+	Node u = _dt->findCost(n, rem);
+
+        if (_tolerance.positive(rem) && !_level->active(w) && w != _source) {
+	  _level->activate(w);
+        }
+
+        if (!_tolerance.less(rem, excess)) {
+	  _dt->addCost(n, - excess);
+	  _excess->set(w, (*_excess)[w] + excess);
+	  excess = 0;
+	  return true;
+	}
+	
+	
+        _dt->addCost(n, - rem);
+
+        excess -= rem;
+        _excess->set(w, (*_excess)[w] + rem);
+        
+	_dt->cut(u);
+	_dt->addCost(u, _max_value);
+	  
+	Edge e = (*_dt_edges)[u];
+	_dt_edges->set(u, INVALID);
+        
+	if (u == _graph.source(e)) {
+	  _flow->set(e, (*_capacity)[e]);
+	} else {
+	  _flow->set(e, 0);
+	}           
+
+        w = _dt->findRoot(n);
+      }      
+      return false;
+    }
+
+    void cutChildren(Node n) {
+    
+      for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+        
+        Node v = _graph.target(e);
+        
+        if ((*_dt_edges)[v] != INVALID && (*_dt_edges)[v] == e) {
+          _dt->cut(v);
+          _dt_edges->set(v, INVALID);
+	  Value rem;
+          _dt->findCost(v, rem);
+          _dt->addCost(v, - rem);
+          _dt->addCost(v, _max_value);
+          _flow->set(e, rem);
+
+        }      
+      }
+
+      for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+        
+        Node v = _graph.source(e);
+        
+        if ((*_dt_edges)[v] != INVALID && (*_dt_edges)[v] == e) {
+          _dt->cut(v);
+          _dt_edges->set(v, INVALID);
+	  Value rem;
+          _dt->findCost(v, rem);
+          _dt->addCost(v, - rem);
+          _dt->addCost(v, _max_value);
+          _flow->set(e, (*_capacity)[e] - rem);
+
+        }      
+      }
+    }
+
+    void extractTrees() {
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+	
+	Node w = _dt->findRoot(n);
+      
+	while (w != n) {
+      
+	  Value rem;      
+	  Node u = _dt->findCost(n, rem);
+
+	  _dt->cut(u);
+	  _dt->addCost(u, - rem);
+	  _dt->addCost(u, _max_value);
+	  
+	  Edge e = (*_dt_edges)[u];
+	  _dt_edges->set(u, INVALID);
+	  
+	  if (u == _graph.source(e)) {
+	    _flow->set(e, (*_capacity)[e] - rem);
+	  } else {
+	    _flow->set(e, rem);
+	  }
+	  
+	  w = _dt->findRoot(n);
+	}      
+      }
+    }
+
+  public:    
+
+    /// \name Execution control The simplest way to execute the
+    /// algorithm is to use one of the member functions called \c
+    /// run().  
+    /// \n
+    /// If you need more control on initial solution or
+    /// execution then you have to call one \ref init() function and then
+    /// the startFirstPhase() and if you need the startSecondPhase().  
+    
+    ///@{
+
+    /// \brief Initializes the internal data structures.
+    ///
+    /// Initializes the internal data structures.
+    ///
+    void init() {
+      createStructures();
+
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+	_excess->set(n, 0);
+      }
+
+      for (EdgeIt e(_graph); e != INVALID; ++e) {
+	_flow->set(e, 0);
+      }
+
+      _dt->clear();
+      for (NodeIt v(_graph); v != INVALID; ++v) {
+        (*_dt_edges)[v] = INVALID;
+        _dt->makeTree(v);
+        _dt->addCost(v, _max_value);
+      }
+
+      typename Graph::template NodeMap<bool> reached(_graph, false);
+
+      _level->initStart();
+      _level->initAddItem(_target);
+
+      std::vector<Node> queue;
+      reached.set(_source, true);
+
+      queue.push_back(_target);
+      reached.set(_target, true);
+      while (!queue.empty()) {
+	_level->initNewLevel();
+	std::vector<Node> nqueue;
+	for (int i = 0; i < int(queue.size()); ++i) {
+	  Node n = queue[i];
+	  for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+	    Node u = _graph.source(e);
+	    if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
+	      reached.set(u, true);
+	      _level->initAddItem(u);
+	      nqueue.push_back(u);
+	    }
+	  }
+	}
+	queue.swap(nqueue);
+      }
+      _level->initFinish();
+
+      for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) {
+	if (_tolerance.positive((*_capacity)[e])) {
+	  Node u = _graph.target(e);
+	  if ((*_level)[u] == _level->maxLevel()) continue;
+	  _flow->set(e, (*_capacity)[e]);
+	  _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
+	  if (u != _target && !_level->active(u)) {
+	    _level->activate(u);
+	  }
+	}
+      }
+    }
+
+    /// \brief Starts the first phase of the preflow algorithm.
+    ///
+    /// The preflow algorithm consists of two phases, this method runs
+    /// the first phase. After the first phase the maximum flow value
+    /// and a minimum value cut can already be computed, although a
+    /// maximum flow is not yet obtained. So after calling this method
+    /// \ref flowValue() returns the value of a maximum flow and \ref
+    /// minCut() returns a minimum cut.     
+    /// \pre One of the \ref init() functions should be called. 
+    void startFirstPhase() {
+      _phase = true;
+      Node n;
+
+      while ((n = _level->highestActive()) != INVALID) {
+	Value excess = (*_excess)[n];
+	int level = _level->highestActiveLevel();
+	int new_level = _level->maxLevel();
+
+	if (_dt->findRoot(n) != n) {
+	  if (sendOut(n, excess)) goto no_more_push;
+	}
+	
+	for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+	  Value rem = (*_capacity)[e] - (*_flow)[e];
+	  Node v = _graph.target(e);
+	  
+	  if (!_tolerance.positive(rem) && (*_dt_edges)[v] != e) continue;
+
+	  if ((*_level)[v] < level) {
+	    
+	    if (_dt->findSize(n) + _dt->findSize(v) < 
+		_tree_bound * _max_tree_size) {
+	      _dt->addCost(n, -_max_value);
+	      _dt->addCost(n, rem);
+	      _dt->link(n, v);
+	      _dt_edges->set(n, e);
+	      if (sendOut(n, excess)) goto no_more_push;
+	    } else {
+	      if (!_level->active(v) && v != _target) {
+		_level->activate(v);
+	      }
+	      if (!_tolerance.less(rem, excess)) {
+		_flow->set(e, (*_flow)[e] + excess);
+		_excess->set(v, (*_excess)[v] + excess);
+		excess = 0;		  
+		goto no_more_push;
+	      } else {
+		excess -= rem;
+		_excess->set(v, (*_excess)[v] + rem);
+		_flow->set(e, (*_capacity)[e]);
+	      }		
+	    }
+	  } else if (new_level > (*_level)[v]) {
+	    new_level = (*_level)[v];
+	  }
+	}
+
+	for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+	  Value rem = (*_flow)[e];
+	  Node v = _graph.source(e);
+	  if (!_tolerance.positive(rem) && (*_dt_edges)[v] != e) continue;
+
+	  if ((*_level)[v] < level) {
+	    
+	    if (_dt->findSize(n) + _dt->findSize(v) < 
+		_tree_bound * _max_tree_size) {
+	      _dt->addCost(n, - _max_value);
+	      _dt->addCost(n, rem);
+	      _dt->link(n, v);
+	      _dt_edges->set(n, e);
+	      if (sendOut(n, excess)) goto no_more_push;
+	    } else {
+	      if (!_level->active(v) && v != _target) {
+		_level->activate(v);
+	      }
+	      if (!_tolerance.less(rem, excess)) {
+		_flow->set(e, (*_flow)[e] - excess);
+		_excess->set(v, (*_excess)[v] + excess);
+		excess = 0;		  
+		goto no_more_push;
+	      } else {
+		excess -= rem;
+		_excess->set(v, (*_excess)[v] + rem);
+		_flow->set(e, 0);
+	      }		
+	    }
+	  } else if (new_level > (*_level)[v]) {
+	    new_level = (*_level)[v];
+	  }
+	}
+		
+      no_more_push:
+
+	_excess->set(n, excess);
+	
+	if (excess != 0) {
+	  cutChildren(n);
+	  if (new_level + 1 < _level->maxLevel()) {
+	    _level->liftHighestActive(new_level + 1);
+	  } else {
+	    _level->liftHighestActiveToTop();
+	  }
+	  if (_level->emptyLevel(level)) {
+	    _level->liftToTop(level);
+	  }
+	} else {
+	  _level->deactivate(n);
+	}	
+      }
+      extractTrees();
+    }
+
+    /// \brief Starts the second phase of the preflow algorithm.
+    ///
+    /// The preflow algorithm consists of two phases, this method runs
+    /// the second phase. After calling \ref init() and \ref
+    /// startFirstPhase() and then \ref startSecondPhase(), \ref
+    /// flowMap() return a maximum flow, \ref flowValue() returns the
+    /// value of a maximum flow, \ref minCut() returns a minimum cut
+    /// \pre The \ref init() and startFirstPhase() functions should be
+    /// called before.
+    void startSecondPhase() {
+      _phase = false;
+      
+      typename Graph::template NodeMap<bool> reached(_graph, false);
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+	reached.set(n, (*_level)[n] < _level->maxLevel());
+      }
+
+      _level->initStart();
+      _level->initAddItem(_source);
+ 
+      std::vector<Node> queue;
+      queue.push_back(_source);
+      reached.set(_source, true);
+
+      while (!queue.empty()) {
+	_level->initNewLevel();
+	std::vector<Node> nqueue;
+	for (int i = 0; i < int(queue.size()); ++i) {
+	  Node n = queue[i];
+	  for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+	    Node v = _graph.target(e);
+	    if (!reached[v] && _tolerance.positive((*_flow)[e])) {
+	      reached.set(v, true);
+	      _level->initAddItem(v);
+	      nqueue.push_back(v);
+	    }
+	  }
+	  for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+	    Node u = _graph.source(e);
+	    if (!reached[u] && 
+		_tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
+	      reached.set(u, true);
+	      _level->initAddItem(u);
+	      nqueue.push_back(u);
+	    }
+	  }
+	}
+	queue.swap(nqueue);
+      }
+      _level->initFinish();
+
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+	if ((*_excess)[n] > 0 && _target != n) {
+	  _level->activate(n);
+	}
+      }
+
+      Node n;
+
+      while ((n = _level->highestActive()) != INVALID) {
+	Value excess = (*_excess)[n];
+	int level = _level->highestActiveLevel();
+	int new_level = _level->maxLevel();
+
+	if (_dt->findRoot(n) != n) {
+	  if (sendIn(n, excess)) goto no_more_push;
+	}
+	
+	for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+	  Value rem = (*_capacity)[e] - (*_flow)[e];
+	  Node v = _graph.target(e);
+	  
+	  if (!_tolerance.positive(rem) && (*_dt_edges)[v] != e) continue;
+
+	  if ((*_level)[v] < level) {
+	    
+	    if (_dt->findSize(n) + _dt->findSize(v) < 
+		_tree_bound * _max_tree_size) {
+	      _dt->addCost(n, -_max_value);
+	      _dt->addCost(n, rem);
+	      _dt->link(n, v);
+	      _dt_edges->set(n, e);
+	      if (sendIn(n, excess)) goto no_more_push;
+	    } else {
+	      if (!_level->active(v) && v != _source) {
+		_level->activate(v);
+	      }
+	      if (!_tolerance.less(rem, excess)) {
+		_flow->set(e, (*_flow)[e] + excess);
+		_excess->set(v, (*_excess)[v] + excess);
+		excess = 0;		  
+		goto no_more_push;
+	      } else {
+		excess -= rem;
+		_excess->set(v, (*_excess)[v] + rem);
+		_flow->set(e, (*_capacity)[e]);
+	      }		
+	    }
+	  } else if (new_level > (*_level)[v]) {
+	    new_level = (*_level)[v];
+	  }
+	}
+
+	for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+	  Value rem = (*_flow)[e];
+	  Node v = _graph.source(e);
+	  if (!_tolerance.positive(rem) && (*_dt_edges)[v] != e) continue;
+
+	  if ((*_level)[v] < level) {
+	    
+	    if (_dt->findSize(n) + _dt->findSize(v) < 
+		_tree_bound * _max_tree_size) {
+	      _dt->addCost(n, - _max_value);
+	      _dt->addCost(n, rem);
+	      _dt->link(n, v);
+	      _dt_edges->set(n, e);
+	      if (sendIn(n, excess)) goto no_more_push;
+	    } else {
+	      if (!_level->active(v) && v != _source) {
+		_level->activate(v);
+	      }
+	      if (!_tolerance.less(rem, excess)) {
+		_flow->set(e, (*_flow)[e] - excess);
+		_excess->set(v, (*_excess)[v] + excess);
+		excess = 0;		  
+		goto no_more_push;
+	      } else {
+		excess -= rem;
+		_excess->set(v, (*_excess)[v] + rem);
+		_flow->set(e, 0);
+	      }		
+	    }
+	  } else if (new_level > (*_level)[v]) {
+	    new_level = (*_level)[v];
+	  }
+	}
+		
+      no_more_push:
+
+	_excess->set(n, excess);
+	
+	if (excess != 0) {
+	  cutChildren(n);
+	  if (new_level + 1 < _level->maxLevel()) {
+	    _level->liftHighestActive(new_level + 1);
+	  } else {
+	    _level->liftHighestActiveToTop();
+	  }
+	  if (_level->emptyLevel(level)) {
+	    _level->liftToTop(level);
+	  }
+	} else {
+	  _level->deactivate(n);
+	}	
+      }
+      extractTrees();
+    }
+
+    /// \brief Runs the Goldberg-Tarjan algorithm.  
+    ///
+    /// Runs the Goldberg-Tarjan algorithm.
+    /// \note pf.run() is just a shortcut of the following code.
+    /// \code
+    ///   pf.init();
+    ///   pf.startFirstPhase();
+    ///   pf.startSecondPhase();
+    /// \endcode
+    void run() {
+      init();
+      startFirstPhase();
+      startSecondPhase();
+    }
+
+    /// \brief Runs the Goldberg-Tarjan algorithm to compute the minimum cut.  
+    ///
+    /// Runs the Goldberg-Tarjan algorithm to compute the minimum cut.
+    /// \note pf.runMinCut() is just a shortcut of the following code.
+    /// \code
+    ///   pf.init();
+    ///   pf.startFirstPhase();
+    /// \endcode
+    void runMinCut() {
+      init();
+      startFirstPhase();
+    }
+
+    /// @}
+
+    /// \name Query Functions
+    /// The result of the %Dijkstra algorithm can be obtained using these
+    /// functions.\n
+    /// Before the use of these functions,
+    /// either run() or start() must be called.
+    
+    ///@{
+
+    /// \brief Returns the value of the maximum flow.
+    ///
+    /// Returns the value of the maximum flow by returning the excess
+    /// of the target node \c t. This value equals to the value of
+    /// the maximum flow already after the first phase.
+    Value flowValue() const {
+      return (*_excess)[_target];
+    }
+
+    /// \brief Returns true when the node is on the source side of minimum cut.
+    ///
+    /// Returns true when the node is on the source side of minimum
+    /// cut. This method can be called both after running \ref
+    /// startFirstPhase() and \ref startSecondPhase().
+    bool minCut(const Node& node) const {
+      return ((*_level)[node] == _level->maxLevel()) == _phase;
+    }
+ 
+    /// \brief Returns a minimum value cut.
+    ///
+    /// Sets the \c cutMap to the characteristic vector of a minimum value
+    /// cut. This method can be called both after running \ref
+    /// startFirstPhase() and \ref startSecondPhase(). The result after second
+    /// phase could be changed slightly if inexact computation is used.
+    /// \pre The \c cutMap should be a bool-valued node-map.
+    template <typename CutMap>
+    void minCutMap(CutMap& cutMap) const {
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+	cutMap.set(n, minCut(n));
+      }
+    }
+
+    /// \brief Returns the flow on the edge.
+    ///
+    /// Sets the \c flowMap to the flow on the edges. This method can
+    /// be called after the second phase of algorithm.
+    Value flow(const Edge& edge) const {
+      return (*_flow)[edge];
+    }
+
+    /// @}
+
+  }; 
+  
+} //namespace lemon
+
+#endif
diff -r 26983135fd6d -r 57143c09dc20 lemon/preflow.h
--- a/lemon/preflow.h	Wed Nov 14 17:53:08 2007 +0000
+++ b/lemon/preflow.h	Sat Nov 17 20:58:11 2007 +0000
@@ -19,897 +19,897 @@
 #ifndef LEMON_PREFLOW_H
 #define LEMON_PREFLOW_H
 
-#include <vector>
-#include <queue>
-
 #include <lemon/error.h>
 #include <lemon/bits/invalid.h>
 #include <lemon/tolerance.h>
 #include <lemon/maps.h>
 #include <lemon/graph_utils.h>
+#include <lemon/elevator.h>
 
 /// \file
 /// \ingroup max_flow
 /// \brief Implementation of the preflow algorithm.
 
-namespace lemon {
+namespace lemon { 
+  
+  /// \brief Default traits class of Preflow class.
+  ///
+  /// Default traits class of Preflow class.
+  /// \param _Graph Graph type.
+  /// \param _CapacityMap Type of capacity map.
+  template <typename _Graph, typename _CapacityMap>
+  struct PreflowDefaultTraits {
 
-  ///\ingroup max_flow
-  ///\brief %Preflow algorithms class.
+    /// \brief The graph type the algorithm runs on. 
+    typedef _Graph Graph;
+
+    /// \brief The type of the map that stores the edge capacities.
+    ///
+    /// The type of the map that stores the edge capacities.
+    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
+    typedef _CapacityMap CapacityMap;
+
+    /// \brief The type of the length of the edges.
+    typedef typename CapacityMap::Value Value;
+
+    /// \brief The map type that stores the flow values.
+    ///
+    /// The map type that stores the flow values. 
+    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
+    typedef typename Graph::template EdgeMap<Value> FlowMap;
+
+    /// \brief Instantiates a FlowMap.
+    ///
+    /// This function instantiates a \ref FlowMap. 
+    /// \param graph The graph, to which we would like to define the flow map.
+    static FlowMap* createFlowMap(const Graph& graph) {
+      return new FlowMap(graph);
+    }
+
+    /// \brief The eleavator type used by Preflow algorithm.
+    /// 
+    /// The elevator type used by Preflow algorithm.
+    ///
+    /// \sa Elevator
+    /// \sa LinkedElevator
+    typedef LinkedElevator<Graph, typename Graph::Node> Elevator;
+    
+    /// \brief Instantiates an Elevator.
+    ///
+    /// This function instantiates a \ref Elevator. 
+    /// \param graph The graph, to which we would like to define the elevator.
+    /// \param max_level The maximum level of the elevator.
+    static Elevator* createElevator(const Graph& graph, int max_level) {
+      return new Elevator(graph, max_level);
+    }
+
+    /// \brief The tolerance used by the algorithm
+    ///
+    /// The tolerance used by the algorithm to handle inexact computation.
+    typedef Tolerance<Value> Tolerance;
+
+  };
+  
+
+  /// \ingroup max_flow
   ///
-  ///This class provides an implementation of the \e preflow \e
-  ///algorithm producing a flow of maximum value in a directed
-  ///graph. The preflow algorithms are the fastest known max flow algorithms. 
-  ///The \e source node, the \e target node, the \e
-  ///capacity of the edges and the \e starting \e flow value of the
-  ///edges should be passed to the algorithm through the
-  ///constructor. It is possible to change these quantities using the
-  ///functions \ref source, \ref target, \ref capacityMap and \ref
-  ///flowMap.
+  /// \brief %Preflow algorithms class.
   ///
-  ///After running \ref lemon::Preflow::phase1() "phase1()"
-  ///or \ref lemon::Preflow::run() "run()", the maximal flow
-  ///value can be obtained by calling \ref flowValue(). The minimum
-  ///value cut can be written into a <tt>bool</tt> node map by
-  ///calling \ref minCut(). (\ref minMinCut() and \ref maxMinCut() writes
-  ///the inclusionwise minimum and maximum of the minimum value cuts,
-  ///resp.)
+  /// This class provides an implementation of the Goldberg's \e
+  /// preflow \e algorithm producing a flow of maximum value in a
+  /// directed graph. The preflow algorithms are the fastest known max
+  /// flow algorithms. The current implementation use a mixture of the
+  /// \e "highest label" and the \e "bound decrease" heuristics.
+  /// The worst case time complexity of the algorithm is \f$O(n^3)\f$.
   ///
-  ///\param Graph The directed graph type the algorithm runs on.
-  ///\param Num The number type of the capacities and the flow values.
-  ///\param CapacityMap The capacity map type.
-  ///\param FlowMap The flow map type.
-  ///\param Tol The tolerance type. 
+  /// The algorithm consists from two phases. After the first phase
+  /// the maximal flow value and the minimum cut can be obtained. The
+  /// second phase constructs the feasible maximum flow on each edge.
   ///
-  ///\author Jacint Szabo 
-  ///\todo Second template parameter is superfluous
-  template <typename Graph, typename Num,
-	    typename CapacityMap=typename Graph::template EdgeMap<Num>,
-            typename FlowMap=typename Graph::template EdgeMap<Num>,
-	    typename Tol=Tolerance<Num> >
+  /// \param _Graph The directed graph type the algorithm runs on.
+  /// \param _CapacityMap The flow map type.
+  /// \param _Traits Traits class to set various data types used by
+  /// the algorithm.  The default traits class is \ref
+  /// PreflowDefaultTraits.  See \ref PreflowDefaultTraits for the
+  /// documentation of a %Preflow traits class. 
+  ///
+  ///\author Jacint Szabo and Balazs Dezso
+#ifdef DOXYGEN
+  template <typename _Graph, typename _CapacityMap, typename _Traits>
+#else
+  template <typename _Graph, 
+	    typename _CapacityMap = typename _Graph::template EdgeMap<int>,
+	    typename _Traits = PreflowDefaultTraits<_Graph, _CapacityMap> >
+#endif
   class Preflow {
-  protected:
-    typedef typename Graph::Node Node;
-    typedef typename Graph::NodeIt NodeIt;
-    typedef typename Graph::EdgeIt EdgeIt;
-    typedef typename Graph::OutEdgeIt OutEdgeIt;
-    typedef typename Graph::InEdgeIt InEdgeIt;
+  public:
+    
+    typedef _Traits Traits;
+    typedef typename Traits::Graph Graph;
+    typedef typename Traits::CapacityMap CapacityMap;
+    typedef typename Traits::Value Value; 
 
-    typedef typename Graph::template NodeMap<Node> NNMap;
-    typedef typename std::vector<Node> VecNode;
+    typedef typename Traits::FlowMap FlowMap;
+    typedef typename Traits::Elevator Elevator;
+    typedef typename Traits::Tolerance Tolerance;
 
-    const Graph* _g;
-    Node _source;
-    Node _target;
+    /// \brief \ref Exception for uninitialized parameters.
+    ///
+    /// This error represents problems in the initialization
+    /// of the parameters of the algorithms.
+    class UninitializedParameter : public lemon::UninitializedParameter {
+    public:
+      virtual const char* what() const throw() {
+	return "lemon::Preflow::UninitializedParameter";
+      }
+    };
+
+  private:
+
+    GRAPH_TYPEDEFS(typename Graph);
+
+    const Graph& _graph;
     const CapacityMap* _capacity;
+
+    int _node_num;
+
+    Node _source, _target;
+
     FlowMap* _flow;
+    bool _local_flow;
 
-    Tol _surely;
-    
-    int _node_num;      //the number of nodes of G
-    
-    typename Graph::template NodeMap<int> level;  
-    typename Graph::template NodeMap<Num> excess;
+    Elevator* _level;
+    bool _local_level;
 
-    // constants used for heuristics
-    static const int H0=20;
-    static const int H1=1;
+    typedef typename Graph::template NodeMap<Value> ExcessMap;
+    ExcessMap* _excess;
+
+    Tolerance _tolerance;
+
+    bool _phase;
+
+    void createStructures() {
+      _node_num = countNodes(_graph);
+
+      if (!_flow) {
+	_flow = Traits::createFlowMap(_graph);
+	_local_flow = true;
+      }
+      if (!_level) {
+	_level = Traits::createElevator(_graph, _node_num);
+	_local_level = true;
+      }
+      if (!_excess) {
+	_excess = new ExcessMap(_graph);
+      }
+    }
+
+    void destroyStructures() {
+      if (_local_flow) {
+	delete _flow;
+      }
+      if (_local_level) {
+	delete _level;
+      }
+      if (_excess) {
+	delete _excess;
+      }
+    }
 
   public:
 
-    ///\ref Exception for the case when s=t.
+    typedef Preflow Create;
 
-    ///\ref Exception for the case when the source equals the target.
-    class InvalidArgument : public lemon::LogicError {
-    public:
-      virtual const char* what() const throw() {
-	return "lemon::Preflow::InvalidArgument";
+    ///\name Named template parameters
+
+    ///@{
+
+    template <typename _FlowMap>
+    struct DefFlowMapTraits : public Traits {
+      typedef _FlowMap FlowMap;
+      static FlowMap *createFlowMap(const Graph&) {
+	throw UninitializedParameter();
       }
     };
-    
-    
-    ///Indicates the property of the starting flow map.
-    
-    ///Indicates the property of the starting flow map.
+
+    /// \brief \ref named-templ-param "Named parameter" for setting
+    /// FlowMap type
     ///
-    enum FlowEnum{
-      ///indicates an unspecified edge map. \c flow will be 
-      ///set to the constant zero flow in the beginning of
-      ///the algorithm in this case.
-      NO_FLOW,
-      ///constant zero flow
-      ZERO_FLOW,
-      ///any flow, i.e. the sum of the in-flows equals to
-      ///the sum of the out-flows in every node except the \c source and
-      ///the \c target.
-      GEN_FLOW,
-      ///any preflow, i.e. the sum of the in-flows is at 
-      ///least the sum of the out-flows in every node except the \c source.
-      PRE_FLOW
+    /// \ref named-templ-param "Named parameter" for setting FlowMap
+    /// type
+    template <typename _FlowMap>
+    struct DefFlowMap 
+      : public Preflow<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> > {
+      typedef Preflow<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> > Create;
     };
 
-    ///Indicates the state of the preflow algorithm.
+    template <typename _Elevator>
+    struct DefElevatorTraits : public Traits {
+      typedef _Elevator Elevator;
+      static Elevator *createElevator(const Graph&, int) {
+	throw UninitializedParameter();
+      }
+    };
 
-    ///Indicates the state of the preflow algorithm.
+    /// \brief \ref named-templ-param "Named parameter" for setting
+    /// Elevator type
     ///
-    enum StatusEnum {
-      ///before running the algorithm or
-      ///at an unspecified state.
-      AFTER_NOTHING,
-      ///right after running \ref phase1()
-      AFTER_PREFLOW_PHASE_1,      
-      ///after running \ref phase2()
-      AFTER_PREFLOW_PHASE_2
+    /// \ref named-templ-param "Named parameter" for setting Elevator
+    /// type
+    template <typename _Elevator>
+    struct DefElevator 
+      : public Preflow<Graph, CapacityMap, DefElevatorTraits<_Elevator> > {
+      typedef Preflow<Graph, CapacityMap, DefElevatorTraits<_Elevator> > Create;
     };
-    
-  protected: 
-    FlowEnum flow_prop;
-    StatusEnum status; // Do not needle this flag only if necessary.
-    
-  public: 
-    ///The constructor of the class.
 
-    ///The constructor of the class. 
-    ///\param _gr The directed graph the algorithm runs on. 
-    ///\param _s The source node.
-    ///\param _t The target node.
-    ///\param _cap The capacity of the edges. 
-    ///\param _f The flow of the edges. 
-    ///\param _sr Tol class.
-    ///Except the graph, all of these parameters can be reset by
-    ///calling \ref source, \ref target, \ref capacityMap and \ref
-    ///flowMap, resp.
-    Preflow(const Graph& _gr, Node _s, Node _t, 
-            const CapacityMap& _cap, FlowMap& _f,
-            const Tol &_sr=Tol()) :
-	_g(&_gr), _source(_s), _target(_t), _capacity(&_cap),
-	_flow(&_f), _surely(_sr),
-	_node_num(countNodes(_gr)), level(_gr), excess(_gr,0), 
-	flow_prop(NO_FLOW), status(AFTER_NOTHING) { 
-	if ( _source==_target )
-	  throw InvalidArgument();
-    }
-    
-    ///Give a reference to the tolerance handler class
+    template <typename _Elevator>
+    struct DefStandardElevatorTraits : public Traits {
+      typedef _Elevator Elevator;
+      static Elevator *createElevator(const Graph& graph, int max_level) {
+	return new Elevator(graph, max_level);
+      }
+    };
 
-    ///Give a reference to the tolerance handler class
-    ///\sa Tolerance
-    Tol &tolerance() { return _surely; }
+    /// \brief \ref named-templ-param "Named parameter" for setting
+    /// Elevator type
+    ///
+    /// \ref named-templ-param "Named parameter" for setting Elevator
+    /// type. The Elevator should be standard constructor interface, ie.
+    /// the graph and the maximum level should be passed to it.
+    template <typename _Elevator>
+    struct DefStandardElevator 
+      : public Preflow<Graph, CapacityMap, 
+		       DefStandardElevatorTraits<_Elevator> > {
+      typedef Preflow<Graph, CapacityMap, 
+		      DefStandardElevatorTraits<_Elevator> > Create;
+    };    
 
-    ///Runs the preflow algorithm.  
+    /// @}
 
-    ///Runs the preflow algorithm.
+    /// \brief The constructor of the class.
     ///
-    void run() {
-      phase1(flow_prop);
-      phase2();
-    }
-    
-    ///Runs the preflow algorithm.  
-    
-    ///Runs the preflow algorithm. 
-    ///\pre The starting flow map must be
-    /// - a constant zero flow if \c fp is \c ZERO_FLOW,
-    /// - an arbitrary flow if \c fp is \c GEN_FLOW,
-    /// - an arbitrary preflow if \c fp is \c PRE_FLOW,
-    /// - any map if \c fp is NO_FLOW.
-    ///If the starting flow map is a flow or a preflow then 
-    ///the algorithm terminates faster.
-    void run(FlowEnum fp) {
-      flow_prop=fp;
-      run();
-    }
-      
-    ///Runs the first phase of the preflow algorithm.
+    /// The constructor of the class. 
+    /// \param graph The directed graph the algorithm runs on. 
+    /// \param capacity The capacity of the edges. 
+    /// \param source The source node.
+    /// \param target The target node.
+    Preflow(const Graph& graph, const CapacityMap& capacity, 
+	       Node source, Node target) 
+      : _graph(graph), _capacity(&capacity), 
+	_node_num(0), _source(source), _target(target), 
+	_flow(0), _local_flow(false),
+	_level(0), _local_level(false),
+	_excess(0), _tolerance(), _phase() {}
 
-    ///The preflow algorithm consists of two phases, this method runs
-    ///the first phase. After the first phase the maximum flow value
-    ///and a minimum value cut can already be computed, although a
-    ///maximum flow is not yet obtained. So after calling this method
-    ///\ref flowValue returns the value of a maximum flow and \ref
-    ///minCut returns a minimum cut.     
-    ///\warning \ref minMinCut and \ref maxMinCut do not give minimum
-    ///value cuts unless calling \ref phase2.  
-    ///\warning A real flow map (i.e. not \ref lemon::NullMap "NullMap")
-    ///is needed for this phase.
-    ///\pre The starting flow must be 
-    ///- a constant zero flow if \c fp is \c ZERO_FLOW, 
-    ///- an arbitary flow if \c fp is \c GEN_FLOW, 
-    ///- an arbitary preflow if \c fp is \c PRE_FLOW, 
-    ///- any map if \c fp is NO_FLOW.
-    void phase1(FlowEnum fp)
-    {
-      flow_prop=fp;
-      phase1();
+    /// \brief Destrcutor.
+    ///
+    /// Destructor.
+    ~Preflow() {
+      destroyStructures();
     }
 
+    /// \brief Sets the capacity map.
+    ///
+    /// Sets the capacity map.
+    /// \return \c (*this)
+    Preflow& capacityMap(const CapacityMap& map) {
+      _capacity = &map;
+      return *this;
+    }
+
+    /// \brief Sets the flow map.
+    ///
+    /// Sets the flow map.
+    /// \return \c (*this)
+    Preflow& flowMap(FlowMap& map) {
+      if (_local_flow) {
+	delete _flow;
+	_local_flow = false;
+      }
+      _flow = &map;
+      return *this;
+    }
+
+    /// \brief Returns the flow map.
+    ///
+    /// \return The flow map.
+    const FlowMap& flowMap() {
+      return *_flow;
+    }
+
+    /// \brief Sets the elevator.
+    ///
+    /// Sets the elevator.
+    /// \return \c (*this)
+    Preflow& elevator(Elevator& elevator) {
+      if (_local_level) {
+	delete _level;
+	_local_level = false;
+      }
+      _level = &elevator;
+      return *this;
+    }
+
+    /// \brief Returns the elevator.
+    ///
+    /// \return The elevator.
+    const Elevator& elevator() {
+      return *_level;
+    }
+
+    /// \brief Sets the source node.
+    ///
+    /// Sets the source node.
+    /// \return \c (*this)
+    Preflow& source(const Node& node) {
+      _source = node;
+      return *this;
+    }
+
+    /// \brief Sets the target node.
+    ///
+    /// Sets the target node.
+    /// \return \c (*this)
+    Preflow& target(const Node& node) {
+      _target = node;
+      return *this;
+    }
+ 
+    /// \brief Sets the tolerance used by algorithm.
+    ///
+    /// Sets the tolerance used by algorithm.
+    Preflow& tolerance(const Tolerance& tolerance) const {
+      _tolerance = tolerance;
+      return *this;
+    } 
+
+    /// \brief Returns the tolerance used by algorithm.
+    ///
+    /// Returns the tolerance used by algorithm.
+    const Tolerance& tolerance() const {
+      return tolerance;
+    } 
+
+    /// \name Execution control The simplest way to execute the
+    /// algorithm is to use one of the member functions called \c
+    /// run().  
+    /// \n
+    /// If you need more control on initial solution or
+    /// execution then you have to call one \ref init() function and then
+    /// the startFirstPhase() and if you need the startSecondPhase().  
     
-    ///Runs the first phase of the preflow algorithm.
+    ///@{
 
-    ///The preflow algorithm consists of two phases, this method runs
-    ///the first phase. After the first phase the maximum flow value
-    ///and a minimum value cut can already be computed, although a
-    ///maximum flow is not yet obtained. So after calling this method
-    ///\ref flowValue returns the value of a maximum flow and \ref
-    ///minCut returns a minimum cut.
-    ///\warning \ref minMinCut() and \ref maxMinCut() do not
-    ///give minimum value cuts unless calling \ref phase2().
-    ///\warning A real flow map (i.e. not \ref lemon::NullMap "NullMap")
-    ///is needed for this phase.
-    void phase1()
-    {
-      int heur0=int(H0*_node_num);  //time while running 'bound decrease'
-      int heur1=int(H1*_node_num);  //time while running 'highest label'
-      int heur=heur1;         //starting time interval (#of relabels)
-      int numrelabel=0;
+    /// \brief Initializes the internal data structures.
+    ///
+    /// Initializes the internal data structures.
+    ///
+    void init() {
+      createStructures();
 
-      bool what_heur=1;
-      //It is 0 in case 'bound decrease' and 1 in case 'highest label'
+      _phase = true;
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+	_excess->set(n, 0);
+      }
 
-      bool end=false;
-      //Needed for 'bound decrease', true means no active 
-      //nodes are above bound b.
+      for (EdgeIt e(_graph); e != INVALID; ++e) {
+	_flow->set(e, 0);
+      }
 
-      int k=_node_num-2;  //bound on the highest level under n containing a node
-      int b=k;    //bound on the highest level under n containing an active node
+      typename Graph::template NodeMap<bool> reached(_graph, false);
 
-      VecNode first(_node_num, INVALID);
-      NNMap next(*_g, INVALID);
+      _level->initStart();
+      _level->initAddItem(_target);
 
-      NNMap left(*_g, INVALID);
-      NNMap right(*_g, INVALID);
-      VecNode level_list(_node_num,INVALID);
-      //List of the nodes in level i<n, set to n.
+      std::vector<Node> queue;
+      reached.set(_source, true);
 
-      preflowPreproc(first, next, level_list, left, right);
-
-      //Push/relabel on the highest level active nodes.
-      while ( true ) {
-	if ( b == 0 ) {
-	  if ( !what_heur && !end && k > 0 ) {
-	    b=k;
-	    end=true;
-	  } else break;
-	}
-
-	if ( first[b]==INVALID ) --b;
-	else {
-	  end=false;
-	  Node w=first[b];
-	  first[b]=next[w];
-	  int newlevel=push(w, next, first);
-	  if ( excess[w] != 0 ) {
-            relabel(w, newlevel, first, next, level_list, 
-                    left, right, b, k, what_heur);
-          }
-
-	  ++numrelabel;
-	  if ( numrelabel >= heur ) {
-	    numrelabel=0;
-	    if ( what_heur ) {
-	      what_heur=0;
-	      heur=heur0;
-	      end=false;
-	    } else {
-	      what_heur=1;
-	      heur=heur1;
-	      b=k;
+      queue.push_back(_target);
+      reached.set(_target, true);
+      while (!queue.empty()) {
+	std::vector<Node> nqueue;
+	for (int i = 0; i < int(queue.size()); ++i) {
+	  Node n = queue[i];
+	  for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+	    Node u = _graph.source(e);
+	    if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
+	      reached.set(u, true);
+	      _level->initAddItem(u);
+	      nqueue.push_back(u);
 	    }
 	  }
 	}
+	queue.swap(nqueue);
       }
-      flow_prop=PRE_FLOW;
-      status=AFTER_PREFLOW_PHASE_1;
-    }
-    // Heuristics:
-    //   2 phase
-    //   gap
-    //   list 'level_list' on the nodes on level i implemented by hand
-    //   stack 'active' on the active nodes on level i      
-    //   runs heuristic 'highest label' for H1*n relabels
-    //   runs heuristic 'bound decrease' for H0*n relabels,
-    //        starts with 'highest label'
-    //   Parameters H0 and H1 are initialized to 20 and 1.
+      _level->initFinish();
 
-
-    ///Runs the second phase of the preflow algorithm.
-
-    ///The preflow algorithm consists of two phases, this method runs
-    ///the second phase. After calling \ref phase1() and then
-    ///\ref phase2(),
-    /// \ref flowMap() return a maximum flow, \ref flowValue
-    ///returns the value of a maximum flow, \ref minCut returns a
-    ///minimum cut, while the methods \ref minMinCut and \ref
-    ///maxMinCut return the inclusionwise minimum and maximum cuts of
-    ///minimum value, resp.  \pre \ref phase1 must be called before.
-    ///
-    /// \todo The inexact computation can cause positive excess on a set of 
-    /// unpushable nodes. We may have to watch the empty level in this case 
-    /// due to avoid the terrible long running time.
-    void phase2()
-    {
-
-      int k=_node_num-2;  //bound on the highest level under n containing a node
-      int b=k;    //bound on the highest level under n of an active node
-
-    
-      VecNode first(_node_num, INVALID);
-      NNMap next(*_g, INVALID); 
-      level.set(_source,0);
-      std::queue<Node> bfs_queue;
-      bfs_queue.push(_source);
-
-      while ( !bfs_queue.empty() ) {
-
-	Node v=bfs_queue.front();
-	bfs_queue.pop();
-	int l=level[v]+1;
-
-	for(InEdgeIt e(*_g,v); e!=INVALID; ++e) {
-	  if ( !_surely.positive((*_capacity)[e] - (*_flow)[e])) continue;
-	  Node u=_g->source(e);
-	  if ( level[u] >= _node_num ) {
-	    bfs_queue.push(u);
-	    level.set(u, l);
-	    if ( excess[u] != 0 ) {
-	      next.set(u,first[l]);
-	      first[l]=u;
-	    }
-	  }
-	}
-
-	for(OutEdgeIt e(*_g,v); e!=INVALID; ++e) {
-	  if ( !_surely.positive((*_flow)[e]) ) continue;
-	  Node u=_g->target(e);
-	  if ( level[u] >= _node_num ) {
-	    bfs_queue.push(u);
-	    level.set(u, l);
-	    if ( excess[u] != 0 ) {
-	      next.set(u,first[l]);
-	      first[l]=u;
-	    }
-	  }
-	}
-      }
-      b=_node_num-2;
-
-      while ( true ) {
-
-	if ( b == 0 ) break;
-	if ( first[b]==INVALID ) --b;
-	else {
-	  Node w=first[b];
-	  first[b]=next[w];
-	  int newlevel=push(w,next, first);
-	  
-          if ( newlevel == _node_num) {
-            excess.set(w, 0);
-	    level.set(w,_node_num);
-          }
-	  //relabel
-	  if ( excess[w] != 0 ) {
-	    level.set(w,++newlevel);
-	    next.set(w,first[newlevel]);
-	    first[newlevel]=w;
-	    b=newlevel;
-	  }
-	} 
-      } // while(true)
-      flow_prop=GEN_FLOW;
-      status=AFTER_PREFLOW_PHASE_2;
-    }
-
-    /// Returns the value of the maximum flow.
-
-    /// Returns the value of the maximum flow by returning the excess
-    /// of the target node \c t. This value equals to the value of
-    /// the maximum flow already after running \ref phase1.
-    Num flowValue() const {
-      return excess[_target];
-    }
-
-
-    ///Returns a minimum value cut.
-
-    ///Sets \c M to the characteristic vector of a minimum value
-    ///cut. This method can be called both after running \ref
-    ///phase1 and \ref phase2. It is much faster after
-    ///\ref phase1.  \pre M should be a bool-valued node-map. \pre
-    ///If \ref minCut() is called after \ref phase2() then M should
-    ///be initialized to false.
-    template<typename _CutMap>
-    void minCut(_CutMap& M) const {
-      switch ( status ) {
-	case AFTER_PREFLOW_PHASE_1:
-	for(NodeIt v(*_g); v!=INVALID; ++v) {
-	  if (level[v] < _node_num) {
-	    M.set(v, false);
-	  } else {
-	    M.set(v, true);
-	  }
-	}
-	break;
-	case AFTER_PREFLOW_PHASE_2:
-	minMinCut(M);
-	break;
-	case AFTER_NOTHING:
-	break;
-      }
-    }
-
-    ///Returns the inclusionwise minimum of the minimum value cuts.
-
-    ///Sets \c M to the characteristic vector of the minimum value cut
-    ///which is inclusionwise minimum. It is computed by processing a
-    ///bfs from the source node \c s in the residual graph.  \pre M
-    ///should be a node map of bools initialized to false.  \pre \ref
-    ///phase2 should already be run.
-    template<typename _CutMap>
-    void minMinCut(_CutMap& M) const {
-
-      std::queue<Node> queue;
-      M.set(_source,true);
-      queue.push(_source);
-      
-      while (!queue.empty()) {
-	Node w=queue.front();
-	queue.pop();
-	
-	for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
-	  Node v=_g->target(e);
-	  if (!M[v] && _surely.positive((*_capacity)[e] -(*_flow)[e]) ) {
-	    queue.push(v);
-	    M.set(v, true);
-	  }
-	}
-	
-	for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
-	  Node v=_g->source(e);
-	  if (!M[v] && _surely.positive((*_flow)[e]) ) {
-	    queue.push(v);
-	    M.set(v, true);
-	  }
-	}
-      }
-    }
-    
-    ///Returns the inclusionwise maximum of the minimum value cuts.
-
-    ///Sets \c M to the characteristic vector of the minimum value cut
-    ///which is inclusionwise maximum. It is computed by processing a
-    ///backward bfs from the target node \c t in the residual graph.
-    ///\pre \ref phase2() or run() should already be run.
-    template<typename _CutMap>
-    void maxMinCut(_CutMap& M) const {
-
-      for(NodeIt v(*_g) ; v!=INVALID; ++v) M.set(v, true);
-
-      std::queue<Node> queue;
-
-      M.set(_target,false);
-      queue.push(_target);
-
-      while (!queue.empty()) {
-        Node w=queue.front();
-	queue.pop();
-
-	for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
-	  Node v=_g->source(e);
-	  if (M[v] && _surely.positive((*_capacity)[e] - (*_flow)[e]) ) {
-	    queue.push(v);
-	    M.set(v, false);
-	  }
-	}
-
-	for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
-	  Node v=_g->target(e);
-	  if (M[v] && _surely.positive((*_flow)[e]) ) {
-	    queue.push(v);
-	    M.set(v, false);
+      for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) {
+	if (_tolerance.positive((*_capacity)[e])) {
+	  Node u = _graph.target(e);
+	  if ((*_level)[u] == _level->maxLevel()) continue;
+	  _flow->set(e, (*_capacity)[e]);
+	  _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
+	  if (u != _target && !_level->active(u)) {
+	    _level->activate(u);
 	  }
 	}
       }
     }
 
-    ///Sets the source node to \c _s.
+    /// \brief Initializes the internal data structures.
+    ///
+    /// Initializes the internal data structures and sets the initial
+    /// flow to the given \c flowMap. The \c flowMap should contain a
+    /// flow or at least a preflow, ie. in each node excluding the
+    /// target the incoming flow should greater or equal to the
+    /// outgoing flow.
+    /// \return %False when the given \c flowMap is not a preflow. 
+    template <typename FlowMap>
+    bool flowInit(const FlowMap& flowMap) {
+      createStructures();
 
-    ///Sets the source node to \c _s.
-    /// 
-    void source(Node _s) { 
-      _source=_s; 
-      if ( flow_prop != ZERO_FLOW ) flow_prop=NO_FLOW;
-      status=AFTER_NOTHING; 
-    }
+      for (EdgeIt e(_graph); e != INVALID; ++e) {
+	_flow->set(e, flowMap[e]);
+      }
 
-    ///Returns the source node.
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+	Value excess = 0;
+	for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+	  excess += (*_flow)[e];
+	}
+	for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+	  excess -= (*_flow)[e];
+	}
+	if (excess < 0 && n != _source) return false;
+	_excess->set(n, excess);
+      }
 
-    ///Returns the source node.
-    /// 
-    Node source() const { 
-      return _source;
-    }
+      typename Graph::template NodeMap<bool> reached(_graph, false);
 
-    ///Sets the target node to \c _t.
+      _level->initStart();
+      _level->initAddItem(_target);
 
-    ///Sets the target node to \c _t.
-    ///
-    void target(Node _t) { 
-      _target=_t; 
-      if ( flow_prop == GEN_FLOW ) flow_prop=PRE_FLOW;
-      status=AFTER_NOTHING; 
-    }
+      std::vector<Node> queue;
+      reached.set(_source, true);
 
-    ///Returns the target node.
-
-    ///Returns the target node.
-    /// 
-    Node target() const { 
-      return _target;
-    }
-
-    /// Sets the edge map of the capacities to _cap.
-
-    /// Sets the edge map of the capacities to _cap.
-    /// 
-    void capacityMap(const CapacityMap& _cap) { 
-      _capacity=&_cap; 
-      status=AFTER_NOTHING; 
-    }
-    /// Returns a reference to capacity map.
-
-    /// Returns a reference to capacity map.
-    /// 
-    const CapacityMap &capacityMap() const { 
-      return *_capacity;
-    }
-
-    /// Sets the edge map of the flows to _flow.
-
-    /// Sets the edge map of the flows to _flow.
-    /// 
-    void flowMap(FlowMap& _f) { 
-      _flow=&_f; 
-      flow_prop=NO_FLOW;
-      status=AFTER_NOTHING; 
-    }
-     
-    /// Returns a reference to flow map.
-
-    /// Returns a reference to flow map.
-    /// 
-    const FlowMap &flowMap() const { 
-      return *_flow;
-    }
-
-  private:
-
-    int push(Node w, NNMap& next, VecNode& first) {
-
-      int lev=level[w];
-      Num exc=excess[w];
-      int newlevel=_node_num;       //bound on the next level of w
-
-      for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
-	if ( !_surely.positive((*_capacity)[e] - (*_flow)[e])) continue;
-	Node v=_g->target(e);
-	
-	if( lev > level[v] ) { //Push is allowed now
-	  
-	  if ( excess[v] == 0 && v!=_target && v!=_source ) {
-	    next.set(v,first[level[v]]);
-	    first[level[v]]=v;
-	  }
-
-	  Num cap=(*_capacity)[e];
-	  Num flo=(*_flow)[e];
-	  Num remcap=cap-flo;
-	  
-	  if ( ! _surely.less(remcap, exc) ) { //A nonsaturating push.
-	    
-	    _flow->set(e, flo+exc);
-	    excess.set(v, excess[v]+exc);
-	    exc=0;
-	    break;
-
-	  } else { //A saturating push.
-	    _flow->set(e, cap);
-	    excess.set(v, excess[v]+remcap);
-	    exc-=remcap;
-	  }
-	} else if ( newlevel > level[v] ) newlevel = level[v];
-      } //for out edges wv
-
-      if ( exc != 0 ) {
-	for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
-	  
-	  if ( !_surely.positive((*_flow)[e]) ) continue;
-	  Node v=_g->source(e);
-	  
-	  if( lev > level[v] ) { //Push is allowed now
-
-	    if ( excess[v] == 0 && v!=_target && v!=_source ) {
-	      next.set(v,first[level[v]]);
-	      first[level[v]]=v;
-	    }
-
-	    Num flo=(*_flow)[e];
-
-	    if ( !_surely.less(flo, exc) ) { //A nonsaturating push.
-
-	      _flow->set(e, flo-exc);
-	      excess.set(v, excess[v]+exc);
-	      exc=0;
-	      break;
-	    } else {  //A saturating push.
-
-	      excess.set(v, excess[v]+flo);
-	      exc-=flo;
-	      _flow->set(e,0);
-	    }
-	  } else if ( newlevel > level[v] ) newlevel = level[v];
-	} //for in edges vw
-
-      } // if w still has excess after the out edge for cycle
-
-      excess.set(w, exc);
-      
-      return newlevel;
-    }
-    
-    
-    
-    void preflowPreproc(VecNode& first, NNMap& next, 
-			VecNode& level_list, NNMap& left, NNMap& right)
-    {
-      for(NodeIt v(*_g); v!=INVALID; ++v) level.set(v,_node_num);
-      std::queue<Node> bfs_queue;
-      
-      if ( flow_prop == GEN_FLOW || flow_prop == PRE_FLOW ) {
-	//Reverse_bfs from t in the residual graph,
-	//to find the starting level.
-	level.set(_target,0);
-	bfs_queue.push(_target);
-	
-	while ( !bfs_queue.empty() ) {
-	  
-	  Node v=bfs_queue.front();
-	  bfs_queue.pop();
-	  int l=level[v]+1;
-	  
-	  for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
-	    if ( !_surely.positive((*_capacity)[e] - (*_flow)[e] )) continue;
-	    Node w=_g->source(e);
-	    if ( level[w] == _node_num && w != _source ) {
-	      bfs_queue.push(w);
-	      Node z=level_list[l];
-	      if ( z!=INVALID ) left.set(z,w);
-	      right.set(w,z);
-	      level_list[l]=w;
-	      level.set(w, l);
+      queue.push_back(_target);
+      reached.set(_target, true);
+      while (!queue.empty()) {
+	_level->initNewLevel();
+	std::vector<Node> nqueue;
+	for (int i = 0; i < int(queue.size()); ++i) {
+	  Node n = queue[i];
+	  for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+	    Node u = _graph.source(e);
+	    if (!reached[u] && 
+		_tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
+	      reached.set(u, true);
+	      _level->initAddItem(u);
+	      nqueue.push_back(u);
 	    }
 	  }
-	  
-	  for(OutEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
-	    if ( !_surely.positive((*_flow)[e]) ) continue;
-	    Node w=_g->target(e);
-	    if ( level[w] == _node_num && w != _source ) {
-	      bfs_queue.push(w);
-	      Node z=level_list[l];
-	      if ( z!=INVALID ) left.set(z,w);
-	      right.set(w,z);
-	      level_list[l]=w;
-	      level.set(w, l);
-	    }
-	  }
-	} //while
-      } //if
-
-
-      switch (flow_prop) {
-	case NO_FLOW:  
-	for(EdgeIt e(*_g); e!=INVALID; ++e) _flow->set(e,0);
-	case ZERO_FLOW:
-	for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
-	
-	//Reverse_bfs from t, to find the starting level.
-	level.set(_target,0);
-	bfs_queue.push(_target);
-	
-	while ( !bfs_queue.empty() ) {
-	  
-	  Node v=bfs_queue.front();
-	  bfs_queue.pop();
-	  int l=level[v]+1;
-	  
-	  for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
-	    Node w=_g->source(e);
-	    if ( level[w] == _node_num && w != _source ) {
-	      bfs_queue.push(w);
-	      Node z=level_list[l];
-	      if ( z!=INVALID ) left.set(z,w);
-	      right.set(w,z);
-	      level_list[l]=w;
-	      level.set(w, l);
+	  for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+	    Node v = _graph.target(e);
+	    if (!reached[v] && _tolerance.positive((*_flow)[e])) {
+	      reached.set(v, true);
+	      _level->initAddItem(v);
+	      nqueue.push_back(v);
 	    }
 	  }
 	}
-	
-	//the starting flow
-	for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
-	  Num c=(*_capacity)[e];
-	  if ( !_surely.positive(c) ) continue;
-	  Node w=_g->target(e);
-	  if ( level[w] < _node_num ) {
-	    if ( excess[w] == 0 && w!=_target ) { //putting into the stack
-	      next.set(w,first[level[w]]);
-	      first[level[w]]=w;
-	    }
-	    _flow->set(e, c);
-	    excess.set(w, excess[w]+c);
+	queue.swap(nqueue);
+      }
+      _level->initFinish();
+
+      for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) {
+	Value rem = (*_capacity)[e] - (*_flow)[e]; 
+	if (_tolerance.positive(rem)) {
+	  Node u = _graph.target(e);
+	  if ((*_level)[u] == _level->maxLevel()) continue;
+	  _flow->set(e, (*_capacity)[e]);
+	  _excess->set(u, (*_excess)[u] + rem);
+	  if (u != _target && !_level->active(u)) {
+	    _level->activate(u);
 	  }
 	}
-	break;
+      }
+      for (InEdgeIt e(_graph, _source); e != INVALID; ++e) {
+	Value rem = (*_flow)[e]; 
+	if (_tolerance.positive(rem)) {
+	  Node v = _graph.source(e);
+	  if ((*_level)[v] == _level->maxLevel()) continue;
+	  _flow->set(e, 0);
+	  _excess->set(v, (*_excess)[v] + rem);
+	  if (v != _target && !_level->active(v)) {
+	    _level->activate(v);
+	  }
+	}
+      }
+      return true;
+    }
 
-	case GEN_FLOW:
-	for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
-	{
-	  Num exc=0;
-	  for(InEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc+=(*_flow)[e];
-	  for(OutEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc-=(*_flow)[e];
-          if (!_surely.positive(exc)) {
-            exc = 0;
-          }
-          excess.set(_target,exc);
+    /// \brief Starts the first phase of the preflow algorithm.
+    ///
+    /// The preflow algorithm consists of two phases, this method runs
+    /// the first phase. After the first phase the maximum flow value
+    /// and a minimum value cut can already be computed, although a
+    /// maximum flow is not yet obtained. So after calling this method
+    /// \ref flowValue() returns the value of a maximum flow and \ref
+    /// minCut() returns a minimum cut.     
+    /// \pre One of the \ref init() functions should be called. 
+    void startFirstPhase() {
+      _phase = true;
+
+      Node n = _level->highestActive();
+      int level = _level->highestActiveLevel();
+      while (n != INVALID) {
+	int num = _node_num;
+
+	while (num > 0 && n != INVALID) {
+	  Value excess = (*_excess)[n];
+	  int new_level = _level->maxLevel();
+
+	  for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+	    Value rem = (*_capacity)[e] - (*_flow)[e];
+	    if (!_tolerance.positive(rem)) continue;
+	    Node v = _graph.target(e);
+	    if ((*_level)[v] < level) {
+	      if (!_level->active(v) && v != _target) {
+		_level->activate(v);
+	      }
+	      if (!_tolerance.less(rem, excess)) {
+		_flow->set(e, (*_flow)[e] + excess);
+		_excess->set(v, (*_excess)[v] + excess);
+		excess = 0;
+		goto no_more_push_1;
+	      } else {
+		excess -= rem;
+		_excess->set(v, (*_excess)[v] + rem);
+		_flow->set(e, (*_capacity)[e]);
+	      }
+	    } else if (new_level > (*_level)[v]) {
+	      new_level = (*_level)[v];
+	    }
+	  }
+
+	  for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+	    Value rem = (*_flow)[e];
+	    if (!_tolerance.positive(rem)) continue;
+	    Node v = _graph.source(e);
+	    if ((*_level)[v] < level) {
+	      if (!_level->active(v) && v != _target) {
+		_level->activate(v);
+	      }
+	      if (!_tolerance.less(rem, excess)) {
+		_flow->set(e, (*_flow)[e] - excess);
+		_excess->set(v, (*_excess)[v] + excess);
+		excess = 0;
+		goto no_more_push_1;
+	      } else {
+		excess -= rem;
+		_excess->set(v, (*_excess)[v] + rem);
+		_flow->set(e, 0);
+	      }
+	    } else if (new_level > (*_level)[v]) {
+	      new_level = (*_level)[v];
+	    }
+	  }
+
+	no_more_push_1:
+
+	  _excess->set(n, excess);
+
+	  if (excess != 0) {
+	    if (new_level + 1 < _level->maxLevel()) {
+	      _level->liftHighestActive(new_level + 1);
+	    } else {
+	      _level->liftHighestActiveToTop();
+	    }
+	    if (_level->emptyLevel(level)) {
+	      _level->liftToTop(level);
+	    }
+	  } else {
+	    _level->deactivate(n);
+	  }
+	  
+	  n = _level->highestActive();
+	  level = _level->highestActiveLevel();
+	  --num;
+	}
+	
+	num = _node_num * 20;
+	while (num > 0 && n != INVALID) {
+	  Value excess = (*_excess)[n];
+	  int new_level = _level->maxLevel();
+
+	  for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+	    Value rem = (*_capacity)[e] - (*_flow)[e];
+	    if (!_tolerance.positive(rem)) continue;
+	    Node v = _graph.target(e);
+	    if ((*_level)[v] < level) {
+	      if (!_level->active(v) && v != _target) {
+		_level->activate(v);
+	      }
+	      if (!_tolerance.less(rem, excess)) {
+		_flow->set(e, (*_flow)[e] + excess);
+		_excess->set(v, (*_excess)[v] + excess);
+		excess = 0;
+		goto no_more_push_2;
+	      } else {
+		excess -= rem;
+		_excess->set(v, (*_excess)[v] + rem);
+		_flow->set(e, (*_capacity)[e]);
+	      }
+	    } else if (new_level > (*_level)[v]) {
+	      new_level = (*_level)[v];
+	    }
+	  }
+
+	  for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+	    Value rem = (*_flow)[e];
+	    if (!_tolerance.positive(rem)) continue;
+	    Node v = _graph.source(e);
+	    if ((*_level)[v] < level) {
+	      if (!_level->active(v) && v != _target) {
+		_level->activate(v);
+	      }
+	      if (!_tolerance.less(rem, excess)) {
+		_flow->set(e, (*_flow)[e] - excess);
+		_excess->set(v, (*_excess)[v] + excess);
+		excess = 0;
+		goto no_more_push_2;
+	      } else {
+		excess -= rem;
+		_excess->set(v, (*_excess)[v] + rem);
+		_flow->set(e, 0);
+	      }
+	    } else if (new_level > (*_level)[v]) {
+	      new_level = (*_level)[v];
+	    }
+	  }
+
+	no_more_push_2:
+
+	  _excess->set(n, excess);
+
+	  if (excess != 0) {
+	    if (new_level + 1 < _level->maxLevel()) {
+	      _level->liftActiveOn(level, new_level + 1);
+	    } else {
+	      _level->liftActiveToTop(level);
+	    }
+	    if (_level->emptyLevel(level)) {
+	      _level->liftToTop(level);
+	    }
+	  } else {
+	    _level->deactivate(n);
+	  }
+
+	  while (level >= 0 && _level->activeFree(level)) {
+	    --level;
+	  }
+	  if (level == -1) {
+	    n = _level->highestActive();
+	    level = _level->highestActiveLevel();
+	  } else {
+	    n = _level->activeOn(level);
+	  }
+	  --num;
+	}
+      }
+    }
+
+    /// \brief Starts the second phase of the preflow algorithm.
+    ///
+    /// The preflow algorithm consists of two phases, this method runs
+    /// the second phase. After calling \ref init() and \ref
+    /// startFirstPhase() and then \ref startSecondPhase(), \ref
+    /// flowMap() return a maximum flow, \ref flowValue() returns the
+    /// value of a maximum flow, \ref minCut() returns a minimum cut
+    /// \pre The \ref init() and startFirstPhase() functions should be
+    /// called before.
+    void startSecondPhase() {
+      _phase = false;
+
+      typename Graph::template NodeMap<bool> reached(_graph, false);
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+	reached.set(n, (*_level)[n] < _level->maxLevel());
+      }
+
+      _level->initStart();
+      _level->initAddItem(_source);
+ 
+      std::vector<Node> queue;
+      queue.push_back(_source);
+      reached.set(_source, true);
+
+      while (!queue.empty()) {
+	_level->initNewLevel();
+	std::vector<Node> nqueue;
+	for (int i = 0; i < int(queue.size()); ++i) {
+	  Node n = queue[i];
+	  for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+	    Node v = _graph.target(e);
+	    if (!reached[v] && _tolerance.positive((*_flow)[e])) {
+	      reached.set(v, true);
+	      _level->initAddItem(v);
+	      nqueue.push_back(v);
+	    }
+	  }
+	  for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+	    Node u = _graph.source(e);
+	    if (!reached[u] && 
+		_tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
+	      reached.set(u, true);
+	      _level->initAddItem(u);
+	      nqueue.push_back(u);
+	    }
+	  }
+	}
+	queue.swap(nqueue);
+      }
+      _level->initFinish();
+
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+	if ((*_excess)[n] > 0 && _target != n) {
+	  _level->activate(n);
+	}
+      }
+
+      Node n;
+      while ((n = _level->highestActive()) != INVALID) {
+	Value excess = (*_excess)[n];
+	int level = _level->highestActiveLevel();
+	int new_level = _level->maxLevel();
+
+	for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+	  Value rem = (*_capacity)[e] - (*_flow)[e];
+	  if (!_tolerance.positive(rem)) continue;
+	  Node v = _graph.target(e);
+	  if ((*_level)[v] < level) {
+	    if (!_level->active(v) && v != _source) {
+	      _level->activate(v);
+	    }
+	    if (!_tolerance.less(rem, excess)) {
+	      _flow->set(e, (*_flow)[e] + excess);
+	      _excess->set(v, (*_excess)[v] + excess);
+	      excess = 0;
+	      goto no_more_push;
+	    } else {
+	      excess -= rem;
+	      _excess->set(v, (*_excess)[v] + rem);
+	      _flow->set(e, (*_capacity)[e]);
+	    }
+	  } else if (new_level > (*_level)[v]) {
+	    new_level = (*_level)[v];
+	  }
 	}
 
-	//the starting flow
-	for(OutEdgeIt e(*_g,_source); e!=INVALID; ++e)	{
-	  Num rem=(*_capacity)[e]-(*_flow)[e];
-	  if ( !_surely.positive(rem) ) continue;
-	  Node w=_g->target(e);
-	  if ( level[w] < _node_num ) {
-	    if ( excess[w] == 0 && w!=_target ) { //putting into the stack
-	      next.set(w,first[level[w]]);
-	      first[level[w]]=w;
-	    }   
-	    _flow->set(e, (*_capacity)[e]);
-	    excess.set(w, excess[w]+rem);
+	for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+	  Value rem = (*_flow)[e];
+	  if (!_tolerance.positive(rem)) continue;
+	  Node v = _graph.source(e);
+	  if ((*_level)[v] < level) {
+	    if (!_level->active(v) && v != _source) {
+	      _level->activate(v);
+	    }
+	    if (!_tolerance.less(rem, excess)) {
+	      _flow->set(e, (*_flow)[e] - excess);
+	      _excess->set(v, (*_excess)[v] + excess);
+	      excess = 0;
+	      goto no_more_push;
+	    } else {
+	      excess -= rem;
+	      _excess->set(v, (*_excess)[v] + rem);
+	      _flow->set(e, 0);
+	    }
+	  } else if (new_level > (*_level)[v]) {
+	    new_level = (*_level)[v];
 	  }
 	}
-	
-	for(InEdgeIt e(*_g,_source); e!=INVALID; ++e) {
-	  if ( !_surely.positive((*_flow)[e]) ) continue;
-	  Node w=_g->source(e);
-	  if ( level[w] < _node_num ) {
-	    if ( excess[w] == 0 && w!=_target ) {
-	      next.set(w,first[level[w]]);
-	      first[level[w]]=w;
-	    }  
-	    excess.set(w, excess[w]+(*_flow)[e]);
-	    _flow->set(e, 0);
+
+      no_more_push:
+
+	_excess->set(n, excess);
+      
+	if (excess != 0) {
+	  if (new_level + 1 < _level->maxLevel()) {
+	    _level->liftHighestActive(new_level + 1);
+	  } else {
+	    // Calculation error 
+	    _level->liftHighestActiveToTop();
 	  }
-	}
-	break;
-
-	case PRE_FLOW:	
-	//the starting flow
-	for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
-	  Num rem=(*_capacity)[e]-(*_flow)[e];
-	  if ( !_surely.positive(rem) ) continue;
-	  Node w=_g->target(e);
-	  if ( level[w] < _node_num ) _flow->set(e, (*_capacity)[e]);
-	}
-	
-	for(InEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
-	  if ( !_surely.positive((*_flow)[e]) ) continue;
-	  Node w=_g->source(e);
-	  if ( level[w] < _node_num ) _flow->set(e, 0);
-	}
-	
-	//computing the excess
-	for(NodeIt w(*_g); w!=INVALID; ++w) {
-	  Num exc=0;
-	  for(InEdgeIt e(*_g,w); e!=INVALID; ++e) exc+=(*_flow)[e];
-	  for(OutEdgeIt e(*_g,w); e!=INVALID; ++e) exc-=(*_flow)[e];
-          if (!_surely.positive(exc)) {
-            exc = 0;
-          }
-	  excess.set(w,exc);
-	  
-	  //putting the active nodes into the stack
-	  int lev=level[w];
-	    if ( exc != 0 && lev < _node_num && Node(w) != _target ) {
-	      next.set(w,first[lev]);
-	      first[lev]=w;
-	    }
-	}
-	break;
-      } //switch
-    } //preflowPreproc
-
-
-    void relabel(Node w, int newlevel, VecNode& first, NNMap& next, 
-		 VecNode& level_list, NNMap& left,
-		 NNMap& right, int& b, int& k, bool what_heur )
-    {
-
-      int lev=level[w];
-
-      Node right_n=right[w];
-      Node left_n=left[w];
-
-      //unlacing starts
-      if ( right_n!=INVALID ) {
-	if ( left_n!=INVALID ) {
-	  right.set(left_n, right_n);
-	  left.set(right_n, left_n);
+	  if (_level->emptyLevel(level)) {
+	    // Calculation error 
+	    _level->liftToTop(level);
+	  }
 	} else {
-	  level_list[lev]=right_n;
-	  left.set(right_n, INVALID);
-	}
-      } else {
-	if ( left_n!=INVALID ) {
-	  right.set(left_n, INVALID);
-	} else {
-	  level_list[lev]=INVALID;
-	}
-      }
-      //unlacing ends
-
-      if ( level_list[lev]==INVALID ) {
-
-	//gapping starts
-	for (int i=lev; i!=k ; ) {
-	  Node v=level_list[++i];
-	  while ( v!=INVALID ) {
-	    level.set(v,_node_num);
-	    v=right[v];
-	  }
-	  level_list[i]=INVALID;
-	  if ( !what_heur ) first[i]=INVALID;
+	  _level->deactivate(n);
 	}
 
-	level.set(w,_node_num);
-	b=lev-1;
-	k=b;
-	//gapping ends
+      }
+    }
 
-      } else {
+    /// \brief Runs the preflow algorithm.  
+    ///
+    /// Runs the preflow algorithm.
+    /// \note pf.run() is just a shortcut of the following code.
+    /// \code
+    ///   pf.init();
+    ///   pf.startFirstPhase();
+    ///   pf.startSecondPhase();
+    /// \endcode
+    void run() {
+      init();
+      startFirstPhase();
+      startSecondPhase();
+    }
 
-	if ( newlevel == _node_num ) level.set(w,_node_num);
-	else {
-	  level.set(w,++newlevel);
-	  next.set(w,first[newlevel]);
-	  first[newlevel]=w;
-	  if ( what_heur ) b=newlevel;
-	  if ( k < newlevel ) ++k;      //now k=newlevel
-	  Node z=level_list[newlevel];
-	  if ( z!=INVALID ) left.set(z,w);
-	  right.set(w,z);
-	  left.set(w,INVALID);
-	  level_list[newlevel]=w;
-	}
+    /// \brief Runs the preflow algorithm to compute the minimum cut.  
+    ///
+    /// Runs the preflow algorithm to compute the minimum cut.
+    /// \note pf.runMinCut() is just a shortcut of the following code.
+    /// \code
+    ///   pf.init();
+    ///   pf.startFirstPhase();
+    /// \endcode
+    void runMinCut() {
+      init();
+      startFirstPhase();
+    }
+
+    /// @}
+
+    /// \name Query Functions
+    /// The result of the %Dijkstra algorithm can be obtained using these
+    /// functions.\n
+    /// Before the use of these functions,
+    /// either run() or start() must be called.
+    
+    ///@{
+
+    /// \brief Returns the value of the maximum flow.
+    ///
+    /// Returns the value of the maximum flow by returning the excess
+    /// of the target node \c t. This value equals to the value of
+    /// the maximum flow already after the first phase.
+    Value flowValue() const {
+      return (*_excess)[_target];
+    }
+
+    /// \brief Returns true when the node is on the source side of minimum cut.
+    ///
+    /// Returns true when the node is on the source side of minimum
+    /// cut. This method can be called both after running \ref
+    /// startFirstPhase() and \ref startSecondPhase().
+    bool minCut(const Node& node) const {
+      return ((*_level)[node] == _level->maxLevel()) == _phase;
+    }
+ 
+    /// \brief Returns a minimum value cut.
+    ///
+    /// Sets the \c cutMap to the characteristic vector of a minimum value
+    /// cut. This method can be called both after running \ref
+    /// startFirstPhase() and \ref startSecondPhase(). The result after second
+    /// phase could be changed slightly if inexact computation is used.
+    /// \pre The \c cutMap should be a bool-valued node-map.
+    template <typename CutMap>
+    void minCutMap(CutMap& cutMap) const {
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+	cutMap.set(n, minCut(n));
       }
-    } //relabel
+    }
 
-  }; 
+    /// \brief Returns the flow on the edge.
+    ///
+    /// Sets the \c flowMap to the flow on the edges. This method can
+    /// be called after the second phase of algorithm.
+    Value flow(const Edge& edge) const {
+      return (*_flow)[edge];
+    }
+    
+    /// @}    
+  };
+}
 
-  ///\ingroup max_flow
-  ///\brief Function type interface for Preflow algorithm.
-  ///
-  ///Function type interface for Preflow algorithm.
-  ///\sa Preflow
-  template<class GR, class CM, class FM>
-  Preflow<GR,typename CM::Value,CM,FM> preflow(const GR &g,
-			    typename GR::Node source,
-			    typename GR::Node target,
-			    const CM &cap,
-			    FM &flow
-			    )
-  {
-    return Preflow<GR,typename CM::Value,CM,FM>(g,source,target,cap,flow);
-  }
-
-} //namespace lemon
-
-#endif //LEMON_PREFLOW_H
+#endif
diff -r 26983135fd6d -r 57143c09dc20 test/preflow_test.cc
--- a/test/preflow_test.cc	Wed Nov 14 17:53:08 2007 +0000
+++ b/test/preflow_test.cc	Sat Nov 17 20:58:11 2007 +0000
@@ -28,7 +28,7 @@
 
 using namespace lemon;
 
-void check_Preflow() 
+void checkPreflow() 
 {
   typedef int VType;
   typedef concepts::Graph Graph;
@@ -37,35 +37,40 @@
   typedef Graph::Edge Edge;
   typedef concepts::ReadMap<Edge,VType> CapMap;
   typedef concepts::ReadWriteMap<Edge,VType> FlowMap;
-  typedef concepts::ReadWriteMap<Node,bool> CutMap;
+  typedef concepts::WriteMap<Node,bool> CutMap;
  
-  typedef Preflow<Graph, int, CapMap, FlowMap> PType;
-
   Graph g;
   Node n;
+  Edge e;
   CapMap cap;
   FlowMap flow;
   CutMap cut;
 
-  PType preflow_test(g,n,n,cap,flow);
+  Preflow<Graph, CapMap>::DefFlowMap<FlowMap>::Create preflow_test(g,cap,n,n);
 
+  preflow_test.capacityMap(cap);
+  flow = preflow_test.flowMap();
+  preflow_test.flowMap(flow);
+  preflow_test.source(n);
+  preflow_test.target(n);
+  
+  preflow_test.init();
+  preflow_test.flowInit(cap);
+  preflow_test.startFirstPhase();
+  preflow_test.startSecondPhase();
   preflow_test.run();
+  preflow_test.runMinCut();
+
   preflow_test.flowValue();
-  preflow_test.source(n);
-  preflow_test.flowMap(flow);
+  preflow_test.minCut(n);
+  preflow_test.minCutMap(cut);
+  preflow_test.flow(e);
 
-  preflow_test.phase1(PType::NO_FLOW);
-  preflow_test.minCut(cut);
-
-  preflow_test.phase2();
-  preflow_test.target(n);
-  preflow_test.capacityMap(cap);
-  preflow_test.minMinCut(cut);
-  preflow_test.maxMinCut(cut);
 }
 
-int cut_value ( SmartGraph& g, SmartGraph::NodeMap<bool>& cut, 
-		SmartGraph::EdgeMap<int>& cap) {
+int cutValue (const SmartGraph& g, 
+	      const SmartGraph::NodeMap<bool>& cut, 
+	      const SmartGraph::EdgeMap<int>& cap) {
   
   int c=0;
   for(SmartGraph::EdgeIt e(g); e!=INVALID; ++e) {
@@ -74,6 +79,29 @@
   return c;
 }
 
+bool checkFlow(const SmartGraph& g, 
+	       const SmartGraph::EdgeMap<int>& flow, 
+	       const SmartGraph::EdgeMap<int>& cap, 
+	       SmartGraph::Node s, SmartGraph::Node t) {
+
+  for (SmartGraph::EdgeIt e(g); e != INVALID; ++e) {
+    if (flow[e] < 0 || flow[e] > cap[e]) return false;
+  }
+
+  for (SmartGraph::NodeIt n(g); n != INVALID; ++n) {
+    if (n == s || n == t) continue;
+    int sum = 0;
+    for (SmartGraph::OutEdgeIt e(g, n); e != INVALID; ++e) {
+      sum += flow[e];
+    }
+    for (SmartGraph::InEdgeIt e(g, n); e != INVALID; ++e) {
+      sum -= flow[e];
+    }
+    if (sum != 0) return false;
+  }
+  return true;
+}
+
 int main() {
 
   typedef SmartGraph Graph;
@@ -85,7 +113,7 @@
   typedef Graph::EdgeMap<int> FlowMap;
   typedef Graph::NodeMap<bool> CutMap;
 
-  typedef Preflow<Graph, int> PType;
+  typedef Preflow<Graph, CapMap> PType;
 
   std::string f_name;
   if( getenv("srcdir") )
@@ -102,74 +130,50 @@
   CapMap cap(g);
   readDimacs(file, g, cap, s, t);
 
-  FlowMap flow(g,0);
+  PType preflow_test(g, cap, s, t);
+  preflow_test.run();
+    
+  check(checkFlow(g, preflow_test.flowMap(), cap, s, t),
+	"The flow is not feasible.");
 
- 
+  CutMap min_cut(g);
+  preflow_test.minCutMap(min_cut); 
+  int min_cut_value=cutValue(g,min_cut,cap);
+  
+  check(preflow_test.flowValue() == min_cut_value,
+	"The max flow value is not equal to the three min cut values.");
 
-  PType preflow_test(g, s, t, cap, flow);
-  preflow_test.run(PType::ZERO_FLOW);
-    
-  CutMap min_cut(g,false);
-  preflow_test.minCut(min_cut); 
-  int min_cut_value=cut_value(g,min_cut,cap);
-   
-  CutMap min_min_cut(g,false);
-  preflow_test.minMinCut(min_min_cut); 
-  int min_min_cut_value=cut_value(g,min_min_cut,cap);
-   
-  CutMap max_min_cut(g,false);
-  preflow_test.maxMinCut(max_min_cut); 
-  int max_min_cut_value=cut_value(g,max_min_cut,cap);
-
-  check(preflow_test.flowValue() == min_cut_value &&
-	min_cut_value == min_min_cut_value &&
-	min_min_cut_value == max_min_cut_value,
-	"The max flow value is not equal to the three min cut values.");
+  FlowMap flow(g);
+  flow = preflow_test.flowMap();
 
   int flow_value=preflow_test.flowValue();
 
+  for(EdgeIt e(g); e!=INVALID; ++e) cap[e]=2*cap[e]; 
+  preflow_test.flowInit(flow);
+  preflow_test.startFirstPhase();
 
-
-  for(EdgeIt e(g); e!=INVALID; ++e) cap[e]=2*cap[e]; 
-  preflow_test.capacityMap(cap);  
-
-  preflow_test.phase1(PType::PRE_FLOW);
-
-  CutMap min_cut1(g,false);
-  preflow_test.minCut(min_cut1); 
-  min_cut_value=cut_value(g,min_cut1,cap);
+  CutMap min_cut1(g);
+  preflow_test.minCutMap(min_cut1); 
+  min_cut_value=cutValue(g,min_cut1,cap);
    
   check(preflow_test.flowValue() == min_cut_value &&
 	min_cut_value == 2*flow_value,
 	"The max flow value or the min cut value is wrong.");
 
-  preflow_test.phase2();
+  preflow_test.startSecondPhase();
 
-  CutMap min_cut2(g,false);
-  preflow_test.minCut(min_cut2); 
-  min_cut_value=cut_value(g,min_cut2,cap);
+  check(checkFlow(g, preflow_test.flowMap(), cap, s, t),
+	"The flow is not feasible.");
+
+  CutMap min_cut2(g);
+  preflow_test.minCutMap(min_cut2); 
+  min_cut_value=cutValue(g,min_cut2,cap);
    
-  CutMap min_min_cut2(g,false);
-  preflow_test.minMinCut(min_min_cut2); 
-  min_min_cut_value=cut_value(g,min_min_cut2,cap);
- 
-  preflow_test.maxMinCut(max_min_cut); 
-  max_min_cut_value=cut_value(g,max_min_cut,cap);
-
   check(preflow_test.flowValue() == min_cut_value &&
-	min_cut_value == min_min_cut_value &&
-	min_min_cut_value == max_min_cut_value &&
 	min_cut_value == 2*flow_value,
 	"The max flow value or the three min cut values were not doubled");
 
 
-
-  EdgeIt e(g);
-  for( int i=1; i==10; ++i ) {
-    flow.set(e,0);
-    ++e;
-  }
-
   preflow_test.flowMap(flow); 
 
   NodeIt tmp1(g,s);
@@ -185,19 +189,13 @@
   
   preflow_test.run();
 
-  CutMap min_cut3(g,false);
-  preflow_test.minCut(min_cut3); 
-  min_cut_value=cut_value(g,min_cut3,cap);
+  CutMap min_cut3(g);
+  preflow_test.minCutMap(min_cut3); 
+  min_cut_value=cutValue(g,min_cut3,cap);
    
-  CutMap min_min_cut3(g,false);
-  preflow_test.minMinCut(min_min_cut3); 
-  min_min_cut_value=cut_value(g,min_min_cut3,cap);
-   
-  preflow_test.maxMinCut(max_min_cut); 
-  max_min_cut_value=cut_value(g,max_min_cut,cap);
 
-  check(preflow_test.flowValue() == min_cut_value &&
-	min_cut_value == min_min_cut_value &&
-	min_min_cut_value == max_min_cut_value,
+  check(preflow_test.flowValue() == min_cut_value,
 	"The max flow value or the three min cut values are incorrect.");
+  
+  return 0;
 }