sap算法详解与模板

关键概念与性质：

距离函数（distance function），我们说一个距离函数是有效的当且仅当满足有效条件（valid function）

（1）d（t）= 0；

（2）d（i）<= d(j)+1(如果弧rij在残余网络G（x）中)；

 

性质1：

如果距离标号是有效的，那么d（i）便是残余网络中从点i到汇点的距离的下限；

证明：

//algorithm shortest augment path;
void sap()
{
x = 0;
obtain the exact distance labels d(i);
i = s;
while (d(s)<n)
{
if (i has an admissible arc)
{
advance(i);
if (i == t)
{
augment;
i = s;
}
}
else
retreat(i);
}
}

//procedure advance(i);
void advance(i)
{
let(i,j)be an admissible arc in A(t);
pre(j) = i;
i = j;
}
//procedure retreat
void retreat()
{
d(i) = min(d(j)+1):rij>0;
if (i != s)
i = pre(i);
}

代码模板：

#include <iostream>
#include <cstring>
#include <cstdlib>

using namespace std;

const int Max = 225;
const int oo = 210000000;

int n,m,c[Max][Max],r[Max][Max],c1[Max][Max],source,sink;
int gap[Max];//用gap来记录当前标号为i的个数，用于gap优化；

//c1是c的反向网络，用于dis数组的初始化
int dis[Max];

void initialize()// 初始化dis数组
{
memset(gap,0,sizeof(gap));
int q[Max],head = 0, tail = 0;//BFS队列
memset(dis,0xff,sizeof(dis));
q[++head] = sink;
dis[sink] = 0;
gap[0]++;
while(tail < head)
{
int u = q[++tail], v;
for(v = 0; v <= sink; v++)
{
if(!dis[v] && c1[u][v] > 0)
{
dis[v] = dis[u] + 1;
gap[dis[v]]++;
q[++head] = v;
}
}
}
}

int maxflow_sap()
{
initialize();
int top = source, pre[Max], flow = 0, i, j, k, low[Max];

// top是当前增广路中最前面一个点。
memset(low,0,sizeof(low));//low数组用于保存路径中的最小容量
while(dis[source] < n)
{
bool flag = false;
low[source] = oo;
for(i = 0; i <= sink; i++)//找允许弧，根据允许弧的定义
{
if(r[top][i] > 0 && dis[top] == dis[i] + 1&&dis[i]>=0)
{
flag = true;
break;
}
}
if(flag)// 如果找到允许弧
{
low[i] = r[top][i];
if(low[i] > low[top])
low[i] = low[top];//更新low
pre[i] = top; top = i;
if(top == sink)// 如果找到一条增广路了
{
flow += low[sink];
j = top;
while(j != source)// 路径回溯更新残留网络
{
k = pre[j];
r[k][j] -= low[sink];
r[j][k] += low[sink];
j = k;
}
top = source;//从头开始再找最短路
memset(low,0,sizeof(low));
}
}
else// 如果没有允许弧
{
int mindis = 10000000;
for(j = 0; j <= sink; j++)//找和top相邻dis最小的点
{
if(r[top][j] > 0 && mindis > dis[j] + 1&&dis[j]>=0)
mindis = dis[j] + 1;
}
gap[dis[top]]--;
if (gap[dis[top]] == 0)
break;
gap[mindis]++;
dis[top] = mindis;//更新top的距离值
if(top != source) top = pre[top];// 回溯找另外的路
}
}
return(flow);
}

 

注意：在运用sap的时候必须要时刻的保证标号的两个性质，因此，不能再初始标号的时候将全部初始为0层，因为有些点是不具有层数的，或者说是层数是无穷大的，不可达的。

网上找的比较清晰地模板sap，有各种优化

[cpp]
view plaincopyprint?

#include <stdio.h>  
#include <string.h>  
#define INF 2100000000  
#define MAXN 301  
  
int SAP(int map[][MAXN],int v_count,int s,int t)      //邻接矩阵，节点总数，始点，汇点  
{  
    int i;  
    int cur_flow,max_flow,cur,min_label,temp;         //当前流，最大流，当前节点，最小标号，临时变量  
    char flag;                                        //标志当前是否有可行流  
    int cur_arc[MAXN],label[MAXN],neck[MAXN];         //当前弧，标号，瓶颈边的入点（姑且这么叫吧）  
    int label_count[MAXN],back_up[MAXN],pre[MAXN];    //标号为i节点的数量，cur_flow的纪录,当前流路径中前驱  
  
    //初始化
  
    memset(label,0,MAXN*sizeof(int));  
    memset(label_count,0,MAXN*sizeof(int));  
  
    memset(cur_arc,0,MAXN*sizeof(int));  
    label_count[0]=v_count;                           //全部初始化为距离为0  
  
    neck[s]=s;  
    max_flow=0;  
    cur=s;  
    cur_flow=INF;  
  
    //循环代替递归
  
    while(label[s]<v_count)  
    {  
        back_up[cur]=cur_flow;  
        flag=0;  
  
        //选择允许路径（此处还可用邻接表优化）  
        for(i=cur_arc[cur];i<v_count;i++)    //当前弧优化  
        {  
           if(map[cur][i]!=0&&label[cur]==label[i]+1)    //找到允许路径  
           {  
               flag=1;  
               cur_arc[cur]=i;    //更新当前弧  
               if(map[cur][i]<cur_flow)    //更新当前流  
               {  
                   cur_flow=map[cur][i];  
                   neck[i]=cur;     //瓶颈为当前节点  
               }  
               else  
               {  
                   neck[i]=neck[cur];     //瓶颈相对前驱节点不变  
               }  
               pre[i]=cur;    //记录前驱  
               cur=i;  
               if(i==t)    //找到可行流  
               {  
                   max_flow+=cur_flow;    //更新最大流  
  
                   //修改残量网络  
                   while(cur!=s)  
                   {  
                       if(map[pre[cur]][cur]!=INF)map[pre[cur]][cur]-=cur_flow;  
                       back_up[cur] -= cur_flow;  
                       if(map[cur][pre[cur]]!=INF)map[cur][pre[cur]]+=cur_flow;  
                       cur=pre[cur];  
                   }  
  
                   //优化，瓶颈之后的节点出栈  
                   cur=neck[t];  
                   cur_flow=back_up[cur];   
               }  
               break;  
           }  
        }  
        if(flag)continue;  
  
        min_label=v_count-1;    //初始化min_label为节点总数-1  
  
        //找到相邻的标号最小的节点     
        for(i=0;i<v_count;i++)  
        {  
            if(map[cur][i]!=0&&label[i]<min_label)  
            {  
                min_label=label[i];  
                temp=i;  
            }  
        }  
        cur_arc[cur]=temp;    //记录当前弧，下次从提供最小标号的节点开始搜索  
        label_count[label[cur]]--;    //修改标号纪录  
        if(label_count[label[cur]]==0)break;    //GAP优化  
        label[cur]=min_label+1;    //修改当前节点标号  
        label_count[label[cur]]++;     //修改标号记录  
        if(cur!=s)  
        {  
           //从栈中弹出一个节点  
           cur=pre[cur];  
           cur_flow=back_up[cur];  
        }  
    }  
    return(max_flow);  
}