_maxflow.c

数据类型和算法库LEDA 数据类型和算法库LEDA

/*******************************************************************************
+
+  LEDA  3.0
+
+
+  _maxflow.c
+
+
+  Copyright (c) 1992  by  Max-Planck-Institut fuer Informatik
+  Im Stadtwald, 6600 Saarbruecken, FRG     
+  All rights reserved.
+ 
*******************************************************************************/


/******************************************************************************\
*                                                                              *
*  MAX_FLOW                                                                    *
*                                                                              *
*  Max Flow Algorithm, using the Goldberg-Tarjan algorithm, Journal ACM 35,    *
*  Oct 1988, p921-940 .                                                        *
*                                                                              *
*  Applies push steps to vertices with positive preflow, based on              *
*  distance labels.  Uses FIFO rule (implemented by FIFO queue)                *
*  for selecting a vertex with +ve preflow. )                                  *
*                                                                              *
*                                                                              *
*  Implemented by                                                              *
*                                                                              *
*  Joseph Cheriyan & Stefan Naeher                                             *
*                                                                              *
*  Fachbereich Informatik                                                      *
*  Universitaet des Saarlandes                                                 *
*  6600 Saarbruecken                                                           *
*                                                                              *
\******************************************************************************/


#include <LEDA/graph_alg.h>
#include <LEDA/queue.h>
#include <math.h>


#if defined(REAL_NUMBERS)
typedef double num_type;
typedef node_array<double> num_node_array;
typedef edge_array<double> num_edge_array;
#else
typedef int num_type;
typedef node_array<int> num_node_array;
typedef edge_array<int> num_edge_array;
#endif



static void gt_bfs(node s_vert, int dist_s_vert, 
                   const graph& G, node s, 
                   const num_edge_array& flow,
                   const num_edge_array& cap_minus_flow,
                         node_array<int>& dist);


num_type MAX_FLOW(graph& G, node s, node t, const num_edge_array& cap, 
                                                  num_edge_array& flow)
{

/* Computes a maximum flow from source s to sink t, using
   Goldberg-Tarjan algorithm [Journal ACM 35 (Oct 1988) p921-940].
   ( Applies push steps to vertices with positive preflow, based on 
     distance labels.  Uses FIFO rule (implemented by FIFO queue) 
     for selecting a vertex with +ve preflow. )

   input:   cap[e]  gives capacity of edge e 
 
   output:  flow[e] flow on edge e
            returns total flow into sink

   node_arrays used by the algorithm:
   dist[v]   = lower bound on shortest distance from v to t
   excess[v] = (total incoming preflow[v]) - (total outgoing preflow[v])

   Implemented by  Joseph Cheriyan & Stefan Naeher.
*/

node_array<int> dist(G,0);
num_node_array  excess(G,0);
num_edge_array  cap_minus_flow(G,0);

list<node> U;     // FIFO Queue used for selecting vertex with +ve preflow
node     v,w;
edge     e;

int N = G.number_of_nodes();

/*
 heuristic: parameter for heuristic suggested by Goldberg to speed up algorithm 
            compute exact distance labels after every "heuristic" relabel steps

            experiments indicate that sqrt(|E|) is a good choice 
*/

int heuristic = int(sqrt(G.number_of_edges()));


int limit_heur   = heuristic;
int num_relabels = 0;

if (s == t) error_handler(1,"MAXFLOW: source == sink");
 
forall_edges(e,G)
{ flow[e] = 0;
  cap_minus_flow[e] = cap[e];
 }
 
forall_adj_edges(e,s)
if (source(e) == s)
  { v = target(e);
    flow[e] = cap[e];
    cap_minus_flow[e] = 0;
    excess[v] += flow[e];
    if (v!=s && v!=t) U.append(v);
  }
 
dist[s] = N;

G.make_undirected();

while (! U.empty() )      /* main loop */
{ 
  v = U.pop();

  forall_adj_edges(e,v)
  { 
    num_type ev = excess[v];

    if(ev <= 0) break;
  
    // push flow along each edge e = (v,w) in the residual graph with
    // dist[v] == dist[w] + 1
  
    if (v == source(e))
      { w = target(e);
        int dv = dist[v];
        if (dv < N && dv == (dist[w]+1))
        { num_type delta = Min(cap_minus_flow[e],ev);
          // pushes towards s that increase flow[e] are not allowed
          if (delta > 0)
          { if (excess[w] == 0 && w != t) U.append(w);
            flow[e] += delta;
            cap_minus_flow[e] -= delta;
            excess[w] += delta;
            ev -= delta;
           }
         }
       }
    else
      { w = source(e);

        if (dist[v] == (dist[w]+1))
        { num_type delta = Min(flow[e],ev);
          if (delta > 0)
          { if (excess[w] == 0 && w != t && w != s) U.append(w);
            flow[e] -= delta;
            cap_minus_flow[e] += delta;
            excess[w] += delta;
            ev -= delta;
           }
         }
       }

     excess[v] = ev;

    }  /* forall_adj_edges(e,v) */

   G.init_adj_iterator(v);
  
   if (excess[v] > 0)
   {
     // relabel vertex v (i.e. update dist[v]) because all
     // admissible edges in the residual graph have been saturated 
 
     num_relabels++;
   
     if (heuristic && (num_relabels == limit_heur))
       { 
         // heuristic suggested by Goldberg to reduce number of relabels:
         // periodically compute exact dist[] labels by two backward bfs 
         // one starting at t and another starting at s
 
         limit_heur += heuristic;
 
        node x;
        forall_nodes(x,G) dist[x] = MAXINT;

        gt_bfs(t,0,G,s,flow,cap_minus_flow,dist);
        gt_bfs(s,N,G,s,flow,cap_minus_flow,dist);
       }
     else
       { 
         int min_dist = MAXINT;
 
         forall_adj_edges(e,v)
         {
           if (source(e) == v && cap_minus_flow[e] > 0 && dist[v] < N)
             //pushes towards s that increase flow[e] are not allowed
             min_dist = Min(min_dist,dist[target(e)]);
 
           if (target(e) == v && flow[e] > 0)
             min_dist = Min(min_dist,dist[source(e)]);
          }
 
         if (min_dist != MAXINT) min_dist++;
         dist[v] = min_dist;
        }
 
     U.push(v);
      
    } // if (excess[v] > 0)
 
 }  // while (!U.empty())  
       
   

/*

//code to verify that flow is legal

num_type tot_s_flow = 0;
forall_adj_edges(e,s)  
   if (source(e) == s)   tot_s_flow += cap[e];
if (tot_s_flow != (excess[s] + excess[t]))
 { cout << "gt_max_flow: total out cap s != excess[s] + excess[t]\n";
   forall_nodes(v,G)
     if ((excess[v] != 0) && (v != t) && (v != s))
      cout << "gt_max_flow: v!=s,t has excess[v]=" << excess[v] << "\n";
  }
*/


G.make_directed();

return excess[t];  // value of maximum flow from s to t
  
}




//procedure to perform backward bfs starting at vertex s_vert with 
//distance label dist_s_vert

static
void gt_bfs(node s_vert, int dist_s_vert,
             const graph& G, node s, 
             const num_edge_array& flow,
             const num_edge_array& cap_minus_flow, 
                   node_array<int>& dist)
{ node x,y;
  edge e;
  queue<node> bfs_Q;
  num_type c;

  dist[s_vert] = dist_s_vert;
  bfs_Q.append(s_vert);

  while (! bfs_Q.empty())
  { 
    x = bfs_Q.pop();

    int d = dist[x]+1;

    forall_adj_edges(e,x)
    {
      if (x == source(e)) 
        { y = target(e); 
          c = flow[e]; 
         }
      else 
        { y = source(e); 
          c = cap_minus_flow[e]; 
         }

      if (c > 0 && dist[y] == MAXINT && (x != target(e) || s_vert != s))
        // while doing bfs with s_vert==s, 
        // only edges with flow[e]>0 are reachable from s
      { 
        dist[y] = d;
        bfs_Q.append(y);

        //next statement is not really needed, but is added for
        //protection, to ensure that s does not get relabelled by
        //bfs starting at t; this may lead to decrease in d[v]!
        if (y == s) dist[y] = G.number_of_nodes();
       }
     }  //forall_adj_edges(e,x)
   }  //while bfs_Q not empty
}