注明这是个转载...来自百度百科的SAP，只是想方便查找...
SAP算法
SAP算法（by dd_engi):求最大流有一种经典的算法，就是每次找增广路时用BFS找，保证找到的增广路是弧数最少的，也就是所谓的Edmonds-Karp算法。可以证明的是在使用最短路增广时增广过程不超过V*E次，每次BFS的时间都是O(E)，所以Edmonds-Karp的时间复杂度就是O(V*E^2)。
如果能让每次寻找增广路时的时间复杂度降下来，那么就能提高算法效率了，使用距离标号的最短增广路算法就是这样的。所谓距离标号，就是某个点到汇点的最少的弧的数量（另外一种距离标号是从源点到该点的最少的弧的数量，本质上没什么区别）。设点i的标号为D[i]，那么如果将满足D[i]=D[j]+1的弧(i,j)叫做允许弧，且增广时只走允许弧，那么就可以达到“怎么走都是最短路”的效果。每个点的初始标号可以在一开始用一次从汇点沿所有反向边的BFS求出，实践中可以初始设全部点的距离标号为0，问题就是如何在增广过程中维护这个距离标号。
维护距离标号的方法是这样的：当找增广路过程中发现某点出发没有允许弧时，将这个点的距离标号设为由它出发的所有弧的终点的距离标号的最小值加一。这种维护距离标号的方法的正确性我就不证了。由于距离标号的存在，由于“怎么走都是最短路”，所以就可以采用DFS找增广路，用一个栈保存当前路径的弧即可。当某个点的距离标号被改变时，栈中指向它的那条弧肯定已经不是允许弧了，所以就让它出栈，并继续用栈顶的弧的端点增广。为了使每次找增广路的时间变成均摊O(V)，还有一个重要的优化是对于每个点保存“当前弧”：初始时当前弧是邻接表的第一条弧；在邻接表中查找时从当前弧开始查找，找到了一条允许弧，就把这条弧设为当前弧；改变距离标号时，把当前弧重新设为邻接表的第一条弧，还有一种在常数上有所优化的写法是改变距离标号时把当前弧设为那条提供了最小标号的弧。当前弧的写法之所以正确就在于任何时候我们都能保证在邻接表中当前弧的前面肯定不存在允许弧。
还有一个常数优化是在每次找到路径并增广完毕之后不要将路径中所有的顶点退栈，而是只将瓶颈边以及之后的边退栈，这是借鉴了Dinic算法的思想。注意任何时候待增广的“当前点”都应该是栈顶的点的终点。这的确只是一个常数优化，由于当前边结构的存在，我们肯定可以在O(n)的时间内复原路径中瓶颈边之前的所有边。
代码：
var
flag:boolean;
jl,min,flow,aug,j,m,n,tmp,a,b,c,i:longint;
his,pre,dis,vh,di:array[0..1024] of longint;
map:array[1..1024,1..1024] of longint;
begin
assign(input,'ditch.in');reset(input);
assign(output,'ditch.out');rewrite(output);
readln(m,n);
for i:=1 to m do
begin
readln(a,b,c);
inc(map[a,b],c);
end;
vh[0]:=n;
for i:=1 to n do di[i]:=1;{当前弧初始为第一条弧}
i:=1;{从S开始搜}
aug:=maxlongint;
while dis[1]<n do
begin
his[i]:=aug;{标记，以便以后返回这个值}
flag:=false;{这个是是否有找到允许弧的标志}
for j:=di[i] to n do{从当前弧开始搜}
begin
if (map[i,j]>0)and(dis[j]+1=dis[i]) then{找到允许弧}
begin
flag:=true;
di[i]:=j;{标记当前弧}
if map[i,j]<aug then aug:=map[i,j];
pre[j]:=i;{记录前驱}
i:=j;
if i=n then{找到增广路}
begin
inc(flow,aug);
while i<>1 do
begin
tmp:=i;
i:=pre[i];
dec(map[i,tmp],aug);
inc(map[tmp,i],aug);
end;
aug:=maxlongint;
end;
break;{找到允许弧则退出查找}
end;
end;
if flag then continue;{找到允许弧}
min:=n-1;{没有允许弧了，需要重标号}
for j:=1 to n do
begin
if (map[i,j]>0)and(dis[j]<min) then
begin
jl:=j;
min:=dis[j];
end;
end;
di[i]:=jl;
dec(vh[dis[i]]);{GAP优化}
if vh[dis[i]]=0 then break;
dis[i]:=min+1;
inc(vh[dis[i]]);
if i<>1 then
begin
i:=pre[i];{返回上一层}
aug:=his[i];{知道之前记录这个值的用处了吧}
end;
end;
write(flow);
close(input);
close(output);
end.
优化：
1.邻接表优化：
如果顶点多的话，往往N^2存不下，这时候就要存边：
存每条边的出发点，终止点和价值，然后排序一下，再记录每个出发点的位置。以后要调用从出发点出发的边时候，只需要从记录的位置开始找即可（其实可以用链表）。优点是时间加快空间节省，缺点是编程复杂度将变大，所以在题目允许的情况下，建议使用邻接矩阵。
2.GAP优化：
如果一次重标号时，出现距离断层，则可以证明ST无可行流，此时则可以直接退出算法。
3.当前弧优化：
为了使每次找增广路的时间变成均摊O(V)，还有一个重要的优化是对于每个点保存“当前弧”：初始时当前弧是邻接表的第一条弧；在邻接表中查找时从当前弧开始查找，找到了一条允许弧，就把这条弧设为当前弧；改变距离标号时，把当前弧重新设为邻接表的第一条弧。
sap算法心得
网络最大流算法是网络流算法的基础，实现方法很多，但时间复杂度与编程复杂
度难于兼顾。一方面，诸如预流推进等高效算法难于编写调试，有时实际效果欠佳
（参见dd_engi的评测）；另一方面，基于增广路的算法大多时间效率不高。于是许多人选择了相对简单的Dinic算法。事实上，SAP算法更易于理解，时间效率更高，编程更简单，思路更清晰。
（名词解释）SAP（Shortest Augmenting Paths）: 最短增广路
算法基本框架：
· 定义距离标号为到汇点距离的下界。
· 在初始距离标号的基础上，不断沿可行弧找增广路增广，一般采用深度优先。可行弧定义为:{(i,j)|h[i]=h[j]+1}
· 遍历当前节点完成后，为了使下次再来的时候有路可走（当然要满足距离标号的性质：不超过真实距离），重标号当前节点为:min{h[j]|(i,j)}+1
· 重标号后当前节点处理完毕。当源点被重标号后，检查其距离标号，当大于顶点数时，图中已不存在可增广路，此时算法结束；否则再次从源点遍历。
· 理论复杂度为:O(n2m)
我的心得：
· 理论上初始标号要用反向BFS求得，实践中可以全部设为0，可以证明：这样做
不改变渐进时间复杂度。
· 理论上可以写出许多子程序并迭代实现，但判断琐碎，没有层次，实践中用递归简单明了，增加的常数复杂度比起SAP速度微乎其微却极大降低编程复杂度，易
于编写调试。
· ★GAP优化★ （性价比极高，推荐！详见程序“//GAP”部分）
注意到我们在某次增广后，最大流可能已经求出，因此算法做了许多无用功。可
以发现，距离标号是单调增的。这启示我们如果标号中存在“间隙”，则图中不
会再有可增广路，于是算法提前终止。实践中我们使用数组vh[i]记录标号为i的
顶点个数，若重标号使得vh中原标号项变为0，则停止算法。
附效率测试与比较：
算法高标推进SAP(GAP) Dinic SAP(NOGAP) BFS+EK
时间(s) 5.67 6.02 11.35 23.90 55.77
注:高标推进程序在180行左右。Dinic用的是dd_engi带BFS的标程，带GAP优化的SAP参见程序。
附源程序与注释（注意以下程序未用BFS，只用到递归，实现简单）：
program ditch; //以USACO的ditch为例
var
c:array[0..1000,0..1000] of cardinal; //邻接矩阵
h,vh:array[0..1000] of cardinal;
n,augc,rd1,rd2,m,i:cardinal; //augc为增广路容量
flow:cardinal=0;
found:boolean; //记录是否已达汇点
procedure aug(const m:cardinal);
var
i,augco,minh:cardinal;
begin
minh:=n-1;
augco:=augc;
if m=n then begin
found:=true;
inc(flow,augc);
exit;
end;
for i:=1 to n do
if c[m,i]>0 then begin
if h[i]+1=h[m] then begin
if c[m,i]<augc then augc:=c[m,i];
aug(i);
if h[1]>=n then exit; //GAP
if found then break;
augc:=augco;
end;
if h[i]<minh then minh:=h[i];
end;
if not found then begin //重标号
dec(vh[h[m]]); //GAP
if vh[h[m]]=0 then h[1]:=n; //GAP
h[m]:=minh+1;
inc(vh[h[m]]); //GAP
end else begin //修改残量
dec(c[m,i],augc);
inc(c[i,m],augc);
end;
end;
begin
assign(input,'ditch.in');
assign(output,'ditch.out');
reset(input);
rewrite(output);
readln(m,n);
fillchar(c,sizeof(c),0);
for i:=1 to m do begin
readln(rd1,rd2,augc);
inc(c[rd1,rd2],augc);
end;
fillchar(h,sizeof(h),0);
fillchar(vh,sizeof(vh),0);
vh[0]:=n;
while h[1]<n do begin
augc:=$FFFFFFFF;
found:=false;
aug(1);
end;
writeln(flow);
close(input);
close(output);
end.
Reference:
1.几种高级网络流算法的实测, dd_engi
2.Network Flows: Theory, Algorithms, and Applications, Ahuja, R.K.
3.http://www.gnocuil.cn/blog/2007/11/网络流的sap算法, gnocuil.
c++语言版本：（pku 1273）
#include<iostream>
using namespace std;
#define maxn 9999999
#define N 110
int n,augc,m,i,flow;
int h[N],vh[N];
int c[N][N];
bool found;
void aug(int m)
{
int i,augco,minh;
minh=n-1;
augco=augc;
if(n==m) { found=true; flow+=augc; return; }
for(i=1;i<=n;i++)
{
if(c[m][i]>0){
if(h[i]+1==h[m]){
if(c[m][i]<augc) augc=c[m][i];
aug(i);
if(h[1]>=n) return;//GAP
if(found) break;
augc=augco;
}
if(h[i]<minh) minh=h[i];
}
}
if(!found){
vh[h[m]]--;
if(vh[h[m]]==0) h[1]=n;
h[m]=minh+1;
vh[h[m]]++;
}
else { c[m][i]-=augc;c[i][m]+=augc; }
}
int main()
{
int x,y,z;
while(scanf("%d%d",&m,&n)!=EOF)
{
memset(c,0,sizeof(c));
for(i=1;i<=m;i++)
{
scanf("%d%d%d",&x,&y,&z);
c[x][y]+=z;
}
memset(h,0,sizeof(h));
memset(vh,0,sizeof(vh));
vh[0]=n;flow=0;
while(h[1]<n)
{
augc=maxn;
found=false;
aug(1);
}
printf("%d\n",flow);
}
return 0;
}
## bfs改进版 pascal ##
program profit;
const
name='profit';
type
toetype=^etype;
etype=record
t,c:longint;
p,next:toetype;
end;
var
e:array[-55000..5500] of toetype;
h:array[-55000..5500] of longint;
vh:array[0..70000] of longint;
n,m,i,rd,rd1,rd2:longint;
per:longint=0;
procedure addedge(const sn,tn,cn:longint);
var
newe:toetype;
begin
new(newe);
with newe^ do begin
t:=tn;
c:=cn;
next:=e[sn];
new(p);
with p^ do begin
t:=sn;
c:=0;
next:=e[tn];
p:=newe;
end;
e[tn]:=p;
end;
e[sn]:=newe;
end;
function aug(const u,t:longint):longint;
var
minh,l,d:longint;
ie:toetype;
begin
if u=0 then begin aug:=t;exit;end;
l:=t;
minh:=n+m+1;
ie:=e[u];
while ie<>nil do begin
if ie^.c>0 then begin
if h[ie^.t]+1=h[u] then begin
if l<ie^.c then d:=l else d:=ie^.c;
d:=aug(ie^.t,d);
dec(ie^.c,d);
inc(ie^.p^.c,d);
dec(l,d);
if h[n+1]>=n+m+2 then begin aug:=t-l;exit;end;
if l=0 then break;
end;
if h[ie^.t]<minh then minh:=h[ie^.t];
end;
ie:=ie^.next;
end;
if l=t then begin
dec(vh[h[u]]);
if vh[h[u]]=0 then h[n+1]:=n+m+2;
h[u]:=minh+1;
inc(vh[h[u]]);
end;
aug:=t-l;
end;
begin
assign(input,name+'.in');
assign(output,name+'.out');
reset(input);
rewrite(output);
readln(n,m);
fillchar(e,sizeof(e),0);
for i:=1 to n do begin
read(rd);
addedge(n+1,i,rd);
end;
readln;
for i:=1 to m do begin
readln(rd1,rd2,rd);
inc(per,rd);
addedge(rd1,-i,maxlongint);
addedge(rd2,-i,maxlongint);
addedge(-i,0,rd);
end;
fillchar(h,sizeof(h),0);
fillchar(vh,sizeof(vh),0);
vh[0]:=n+m+2;
while h[n+1]<n+m+2 do
dec(per,aug(n+1,maxlongint));
writeln(per);
close(input);
close(output);
end.
==这个是SAP==
program ditch;
Const Maxn = 3000;
Inf=MaxLongint;
Fin ='Ditch.in';
Fout = 'Ditch.out';
Var Cf,g:Array[0..Maxn,0..Maxn]Of Longint;
d,pre,list,Min,Cur,Numb:Array[0..Maxn]Of Longint;
n,m,MaxFlow:Longint;
Procedure Init;
Var i,a,b,c,Open,Closed:Longint;
Begin
Assign(Input,Fin);Reset(Input);
Readln(m,n);
For i:=1 To m Do Begin
Readln(a,b,c);
Inc(Cf[a,b],c);
Inc(G[a,0]);G[a,G[a,0]]:=b;
Inc(G[b,0]);G[b,G[b,0]]:=a;
End;
Fillchar(D,Sizeof(D),$7f);
D[n]:=0;List[1]:=n;Open:=1;Closed:=0;
While Open>Closed Do Begin
Inc(Closed);a:=List[Closed];
For i:=1 To G[a,0] Do
If D[a]+1<D[G[a,i]] Then Begin
Inc(Open);List[Open]:=G[a,i];D[G[a,i]]:=D[a]+1;
End;
End;
For i:=1 To n Do Inc(Numb[D[i]]);
Close(Input);
End;
Procedure Modify;
Var i:Longint;
Begin
i:=n;
While Pre[i]<>0 Do Begin
Inc(Cf[i,Pre[i]],Min[n]);
Dec(Cf[Pre[i],i],Min[n]);
i:=Pre[i];
End;
Inc(MaxFlow,Min[n]);
End;
Procedure Print;
Begin
Assign(Output,Fout);Rewrite(Output);
Writeln(MaxFlow);
Close(Output);
Halt;
End;
Procedure Main;
Var i,j,k:Longint;
Flag:Boolean;
Begin
i:=1;MaxFlow:=0;Min[1]:=Inf;
While D[1]<n Do Begin
Flag:=False;k:=Inf;
While Cur[i]<=G[i,0] Do Begin
If (Cf[i,G[i,cur[i]]] > 0)And(D[G[i,cur[i]]]+1=D[i]) Then Begin
Flag:=True;
Pre[G[i,cur[i]]]:=i;
If Cf[i,G[i,cur[i]]]<Min[i] Then
Min[G[i,cur[i]]]:=Cf[i,G[i,cur[i]]]
Else Min[G[i,cur[i]]]:=Min[i];
i:=G[i,cur[i]];
If i=n Then Begin
Modify;
i:=1;
End;
End
Else Inc(Cur[i]);
End;
If Not Flag Then Begin
Dec(Numb[D[i]]);
If Numb[D[i]]=0 Then Print;
For j:=1 To G[i,0] Do If (Cf[i,G[i,j]]>0) And (D[G[i,j]]<k) Then
k:=D[G[i,j]];
If k<Inf Then
D[i]:=k+1
Else
D[i]:=n;
Cur[i]:=1;
Inc(Numb[D[i]]);
If i<>1 Then i:=Pre[i];
End;
End;
End;
begin
Init;
Main;
Print;
end.
## c++ 语言版 sap ##
//sap 网络最大流算法
#include<iostream>
using namespace std;
const int N = 3000;
const int inf=0x7fffffff;
int cf[N][N],g[N][N];
int d[N],pre[N],list[N],minn[N],cur[N],numb[N];
int n,m,maxflow,i,j,x,y,z,open,closed;
void init()
{
memset(cf,0,sizeof(cf));
memset(g,0,sizeof(g));
memset(d,inf,sizeof(d));
memset(list,0,sizeof(list));
memset(pre,0,sizeof(pre));
memset(numb,0,sizeof(numb));
memset(minn,0,sizeof(minn));
memset(cur,0,sizeof(cur));
for(i=1;i<=m;i++)
{
scanf("%d%d%d",&x,&y,&z);
cf[x][y]+=z;
g[x][0]++;g[x][g[x][0]]=y;
g[y][0]++;g[y][g[y][0]]=x;
}
d[n]=0;list[1]=n;open=1;closed=0;
while(open>closed)
{
closed++;x=list[closed];
for(i=1;i<=g[x][0];i++)
if(d[x]+1<d[g[x][i]]){
open++;
list[open]=g[x][i];
d[g[x][i]]=d[x]+1;
}
}
for(i=1;i<=n;i++) numb[d[i]]++;
}
void modify()
{
int i=n;
while(pre[i]!=0)
{
cf[i][pre[i]]+=minn[n];
cf[pre[i]][i]-=minn[n];
i=pre[i];
}
maxflow+=minn[n];
}
void sap()
{
int i,j,k;bool flag;
i=1;maxflow=0;minn[1]=inf;
while(d[1]<n)
{
flag=false;k=inf;
while(cur[i]<=g[i][0]) {
if(cf[i][g[i][cur[i]]]>0 && d[g[i][cur[i]]]+1==d[i]){
flag=true;
pre[g[i][cur[i]]]=i;
if(cf[i][g[i][cur[i]]]<minn[i]) minn[g[i][cur[i]]]=cf[i][g[i][cur[i]]];
else minn[g[i][cur[i]]]=minn[i];
i=g[i][cur[i]];
if(i==n) { modify(); i=1; }
}
else cur[i]++;
}
if(!flag){
numb[d[i]]--;
if(numb[d[i]]==0) return;
for(j=1;j<=g[i][0];j++)
if(cf[i][g[i][j]]>0 && d[g[i][j]]<k) k=d[g[i][j]];
if(k<inf) d[i]=k+1; else d[i]=n;
cur[i]=1;
numb[d[i]]++;
if(i!=1) i=pre[i];
}
}
}
int main()
{
while(scanf("%d%d",&m,&n)!=EOF)
{
init();
sap();
printf("%d\n",maxflow);
}
return 0;
}
评论