ISPA 算法

 4. Improved SAP 算法

 本次介绍的重头戏。通常的 SAP 类算法在寻找增广路时总要先进行 BFS，BFS 的最坏情况下复杂度为 O(E)，这样使得普通 SAP 类算法最坏情况下时间复杂度达到了 O(VE²)。为了避免这种情况，Ahuja 和 Orlin 在1987年提出了Improved SAP 算法，它充分利用了距离标号的作用，每次发现顶点无出弧时不是像 Dinic 算法那样到最后进行 BFS，而是就地对顶点距离重标号，这样相当于在遍历的同时顺便构建了新的分层网络，每轮寻找之间不必再插入全图的 BFS 操作，极大提高了运行效率。国内一般把这个算法称为 SAP...显然这是不准确的，毕竟从字面意思上来看 E-K 和 Dinic 都属于 SAP，我还是习惯称为 ISAP 或改进的 SAP 算法。

 与 Dinic 算法不同，ISAP 中的距离标号是每个顶点到达终点 t 的距离。同样也不需显式构造分层网络，只要保存每个顶点的距离标号即可。程序开始时用一个反向 BFS 初始化所有顶点的距离标号，之后从源点开始，进行如下三种操作：(1)当前顶点 i 为终点时增广 (2) 当前顶点有满足 dist[i] = dist[j] + 1 的出弧时前进 (3) 当前顶点无满足条件的出弧时重标号并回退一步。整个循环当源点 s 的距离标号 dist[s] >= n 时结束。对 i 点的重标号操作可概括为 dist[i] = 1 + min{dist[j] : (i,j)属于残量网络G_f}。具体算法描述如下：

algorithm Improved-Shortest-Augmenting-Path 1 f <-- 0 2 从终点 t 开始进行一遍反向 BFS 求得所有顶点的起始距离标号 d(i) 3 i <-- s 4 while d(s) < n do 5 if i = t then // 找到增广路 6 Augment 7 i <-- s // 从源点 s 开始下次寻找 8 if Gf 包含从 i 出发的一条允许弧 (i,j) 9 then Advance(i) 10 else Retreat(i) // 没有从 i 出发的允许弧则回退 11 return f procedure Advance(i) 1 设 (i,j) 为从 i 出发的一条允许弧 2 pi(j) <-- i // 保存一条反向路径，为回退时准备 3 i <-- j // 前进一步，使 j 成为当前结点 procedure Retreat(i) 1 d(i) <-- 1 + min{d(j):(i,j)属于残量网络Gf} // 称为重标号的操作 2 if i != s 3 then i <-- pi(i) // 回退一步 procedure Augment 1 pi 中记录为当前找到的增广路 P 2 delta <-- min{rij:(i,j)属于P} 3 沿路径 P 增广 delta 的流量 4 更新残量网络 Gf

算法中的允许弧是指在残量网络中满足 dist[i] = dist[j] + 1 的弧。Retreat 过程中若从 i 出发没有弧属于残量网络 G_f 则把顶点距离重标号为 n 。

 虽然 ISAP 算法时间复杂度与 Dinic 相同都是 O(V²E)，但在实际表现中要好得多。要提的一点是关于 ISAP 的一个所谓 GAP 优化。由于从 s 到 t 的一条最短路径的顶点距离标号单调递减，且相邻顶点标号差严格等于1，因此可以预见如果在当前网络中距离标号为 k (0 <= k < n) 的顶点数为 0，那么可以知道一定不存在一条从 s 到 t 的增广路径，此时可直接跳出主循环。在我的实测中，这个优化是绝对不能少的，一方面可以提高速度，另外可增强 ISAP 算法时间上的稳定性，不然某些情况下 ISAP 会出奇的费时，而且大大慢于 Dinic 算法。相对的，初始的一遍 BFS 却是可有可无，因为 ISAP 可在循环中自动建立起分层网络。实测加不加 BFS 运行时间差只有 5% 左右，代码量可节省 15~20 行。

下面是一个牛人的模板：

#include <stdio.h> #define N 2007 // Number of nodes #define oo 1000000000 // Infinity // Nodes, Arcs, the source node and the sink node int n, m, source, sink; // Matrixes for maintaining // Graph and Flow int G[N][N], F[N][N]; int pi[N]; // predecessor list int CurrentNode[N]; // Current edge for each node int queue[N]; // Queue for reverse BFS int d[N]; // Distance function int numbs[N]; // numbs[k] is the number of nodes i with d[i]==k // Reverse breadth-first search // to establish distance function d int rev_BFS() { int i, j, head(0), tail(0); // Initially, all d[i]=n for(i = 1; i <= n; i++) numbs[ d[i] = n ] ++; // Start from the sink numbs[n]--; d[sink] = 0; numbs[0]++; queue[ ++tail ] = sink; // While queue is not empty while( head != tail ) { i = queue[++head]; // Get the next node // Check all adjacent nodes for(j = 1; j <= n; j++) { // If it was reached before or there is no edge // then continue if(d[j] < n || G[j][i] == 0) continue; // j is reached first time // put it into queue queue[ ++tail ] = j; // Update distance function numbs[n]--; d[j]=d[i]+1; numbs[d[j]]++; } } return 0; } // Augmenting the flow using predecessor list pi[] int Augment() { int i, j, tmp, width(oo); // Find the capacity of the path for(i = sink, j = pi[i]; i != source; i = j, j = pi[j]) { tmp = G[j][i]; if(tmp < width) width = tmp; } // Augmentation itself for(i = sink, j = pi[i]; i != source; i = j, j = pi[j]) { G[j][i] -= width; F[j][i] += width; G[i][j] += width; F[i][j] -= width; } return width; } // Relabel and backtrack int Retreat(int &i) { int tmp; int j, mind(n-1); // Check all adjacent edges // to find nearest for(j=1; j <= n; j++) // If there is an arc // and j is "nearer" if(G[i][j] > 0 && d[j] < mind) mind = d[j]; tmp = d[i]; // Save previous distance // Relabel procedure itself numbs[d[i]]--; d[i] = 1 + mind; numbs[d[i]]++; // Backtrack, if possible (i is not a local variable! ) if( i != source ) i = pi[i]; // If numbs[ tmp ] is zero, algorithm will stop return numbs[ tmp ]; } // Main procedure int find_max_flow() { int flow(0), i, j; rev_BFS(); // Establish exact distance function // For each node current arc is the first arc for(i=1; i<=n; i++) CurrentNode[i] = 1; // Begin searching from the source i = source; // The main cycle (while the source is not "far" from the sink) for( ; d[source] < n ; ) { // Start searching an admissible arc from the current arc for(j = CurrentNode[i]; j <= n; j++) // If the arc exists in the residual network // and if it is an admissible if( G[i][j] > 0 && d[i] == d[j] + 1 ) // Then finish searhing break; // If the admissible arc is found if( j <= n ) { CurrentNode[i] = j; // Mark the arc as "current" pi[j] = i; // j is reachable from i i = j; // Go forward // If we found an augmenting path if( i == sink ) { flow += Augment(); // Augment the flow i = source; // Begin from the source again } } // If no an admissible arc found else { CurrentNode[i] = 1; // Current arc is the first arc again // If numbs[ d[i] ] == 0 then the flow is the maximal if( Retreat(i) == 0 ) break; } } // End of the main cycle // We return flow value return flow; } // The main function // Graph is represented in input as triples <from, to, capacity> // No comments here int main() { int i, p, q, r; scanf("%d %d %d %d", &n, &m, &source, &sink); for(i = 0; i < m; i++) { scanf("%d %d %d", &p, &q, &r); G[p][q] += r; } printf("%d", find_max_flow()); return 0; }

原文请参考这个牛的blog：http://www.topcoder.com/tc?module=Static&d1=tutorials&d2=maxFlowRevisited#1