HDU 3549 Flow Problem【SAP】

Problem Description
Network flow is a well-known difficult problem for ACMers. Given a graph, your task is to find out the maximum flow for the weighted directed graph.
Input
The first line of input contains an integer T, denoting the number of test cases.
For each test case, the first line contains two integers N and M, denoting the number of vertexes and edges in the graph. (2 <= N <= 15, 0 <= M <= 1000)
Next M lines, each line contains three integers X, Y and C, there is an edge from X to Y and the capacity of it is C. (1 <= X, Y <= N, 1 <= C <= 1000)
Output
For each test cases, you should output the maximum flow from source 1 to sink N.
Sample Input
2 3 2 1 2 1 2 3 1 3 3 1 2 1 2 3 1 1 3 1
Sample Output
Case 1: 1
Case 2: 2
分析; SAP 加了GAP优化
邻接矩阵  
#include<stdio.h>
#include<string.h>
#include<stdlib.h>
#define INF 210000000
int n,m,s,u;
int c[17][17],r[17][17];
int gap[17];
int dis[17];
void init()
{
    int v,x,q[17],front=0,rear=0;
    memset(gap,0,sizeof(gap));
    memset(dis,0xff,sizeof(dis));
    q[rear++]=u;
    dis[u]=0;
    while(front<rear)
    {
        x=q[front++];
        gap[dis[x]]++;
        for(v=1;v<u;v++)
        {
            if(dis[v]==-1&&c[v][x]>0)
            {
                dis[v]=dis[x]+1;
                q[rear++]=v;
            }
        }
    }
}
int max_flow()
{
    init(); 
    int flow=0,top=s,i,j,k,pre[17],low[17];
    memset(low,0,sizeof(low));
    while(dis[s]<n)
    {
        int flag=0;
        low[s]=INF;
        for(i=1;i<=u;i++)
            if(c[top][i]>0&&dis[top]==dis[i]+1&&dis[i]>=0)
            {
                flag=1;
                break;
            }
        if(flag)
        {
            low[i]=c[top][i];
            if(low[i]>low[top])
                low[i]=low[top];
            pre[i]=top;
            top=i;
            if(top==u)
            {
                flow+=low[u];
                j=top;
                while(j!=s)
                {
                    k=pre[j];
                    c[k][j]-=low[u];
                    c[j][k]+=low[u];
                    j=k;
                }
                top=s;
                memset(low,0,sizeof(low));
            }
        }
        else 
        {
            int min=n-1;
            for(j=1;j<u;j++)
            {
                if(c[top][j]>0&&dis[j]+1<min&&dis[j]>=0)
                    min=dis[j]+1;
            }
            gap[dis[top]]--;
            if(gap[dis[top]]==0)
                break;
            gap[min]++;
            dis[top]=min;
            if(top!=s)
                top=pre[top];
        }
    }
    return flow;
}
int main()
{
    int a,b,w,i,ca,cases=1;
    scanf("%d",&ca);
    while(ca--)
    {
        scanf("%d%d",&n,&m);
        memset(c,0,sizeof(c));
        for(i=1;i<=m;i++)
        {
            scanf("%d%d%d",&a,&b,&w);
            c[a][b]+=w;
        }
        s=1;
        u=n;
        printf("Case %d: %d\n",cases++,max_flow());
    }
    return 0;
}
前向星 
#include<stdio.h>
#include<string.h>
#include<stdlib.h>
#define min(a,b)(a)<(b)?(a):(b)
const int INF=0x1f1f1f1f;
const int maxn=16;
const int maxm=100003;
struct node
{
    int c,next,to;
}e[maxm];
int tot;
int head[maxn];
void add(int s,int u,int flow)
{
    e[tot].to=u;
    e[tot].c=flow;
    e[tot].next=head[s];
    head[s]=tot++;
}
int max_flow(int st,int end,int n)
{
    int numh[maxn],h[maxn],curedge[maxn],pre[maxn];
    // numh:    用于GAP 优化的统计高度数量数组；
    // h:       距离标号数组
    // curedge: 当前弧数组
    // pre:     前驱数组
    int cur_flow,maxflow=0,u,tmp,neck,i;
    memset(h,0,sizeof(h));
    memset(numh,0,sizeof(numh));
    memset(pre,0xff,sizeof(pre));
    for(i=1;i<=n;i++)
        curedge[i]=head[i]; // 初始话当前弧为第一条邻接边
    numh[0]=n;
    u=st;
    while(h[st]<n)          //  当h[st]>=n 是，网络中肯定出现了 GAP 
    {
        if(u==end)
        {
            cur_flow=INF;
            for(i=st;i!=end;i=e[curedge[i]].to)
                if(cur_flow>e[curedge[i]].c)
                {
                    neck=i;
                    cur_flow=e[curedge[i]].c;
                }   // 增广成功，寻找‘瓶颈’ 边
            for(i=st;i!=end;i=e[curedge[i]].to)
            {
                tmp=curedge[i];
                e[tmp].c-=cur_flow;
                e[tmp^1].c+=cur_flow; // 修改路径上的边容量
            }
            maxflow+=cur_flow;
            u=neck;      // 下次增广从瓶颈开始
        }
        for(i=curedge[u];i!=-1;i=e[i].next)
            if(e[i].c&&h[u]==h[e[i].to]+1)
                break;   //  寻找可行弧
        if(i!=-1)
        {
            curedge[u]=i;
            pre[e[i].to]=u;
            u=e[i].to;
        }
        else 
        {
            if(--numh[h[u]]==0) break;   // GAP 优化
            curedge[u]=head[u];
            for(tmp=n,i=head[u];i!=-1;i=e[i].next)
                if(e[i].c)
                    tmp=min(tmp,h[e[i].to]);
            h[u]=tmp+1;
            ++numh[h[u]];
            if(u!=st)          // 重新标号并从当前点前驱重新增广
                u=pre[u];
        }
    }
    return maxflow;
}
int main()
{
    int flow,n,m,a,b,t,w,ca=1;
//    freopen("D:ce.txt","r",stdin);
    scanf("%d",&t);
    while(t--)
    {
        scanf("%d%d",&n,&m);
        tot=0;
        memset(head,0xff,sizeof(head));
        while(m--)
        {
            scanf("%d%d%d",&a,&b,&w);
            add(a,b,w);
            add(b,a,0);
            //add(b,a,w);
        }
        flow=max_flow(1,n,n);
        printf("Case %d: %d\n",ca++,flow);
    }
    return 0;
}

 

 SAP详解

关键概念与性质：

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

（1）d（t）= 0；

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

 

性质1：

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

证明：

 

[cpp] view plaincopy
  1. 令任意的从i结点到汇点t长度为k的路径，设为i= i1-i2-i3...ik+1=t  
  2. d(i(k)) <= d(i(k+1))+1 = d(t)+1 = 1  
  3. d(i(k-1)) <= d(i(k))+1 = 2  
  4. d(i(k-2)) <= d(i(k-1))+1 = 3  
  5.             :  
  6.             :  
  7. d(i(1)) <= d(i(2))+1 = k  
  8. 由此可见，d(i)是i到t的距离下限       

 

 

性质2：

允许弧（边）：如果对于边rij>0,且d(i)= d(j)+1，那么称为允许边

允许路：一条从源点到汇点的且有允许弧组成的路径

允许路是从源点到汇点的最短增广路径。

证明：

（1）因为rij>0所以必然是一条增广路

（2）假设路径p是一条允许路包含k条弧，那么由d(i) = d(j)+1 可知，d（s）= k；

又因为d（s）是点s到汇点的距离下限，所以距离下限为k，所以p便是一条最短路。

 

性质3：

在sap算法中距离标号始终是正确，有效的。并且每次的冲标号都会是距离标号严格递增

证明：略；

 

伪代码：

 

[cpp] view plaincopy
  1. //algorithm shortest augment path;  
  2. void sap()  
  3. {  
  4.     x = 0;  
  5.     obtain the exact distance labels d(i);  
  6.     i = s;  
  7.     while (d(s)<n)  
  8.     {  
  9.         if (i has an admissible arc)  
  10.         {  
  11.             advance(i);  
  12.             if (i == t)  
  13.             {  
  14.                 augment;  
  15.                 i = s;  
  16.             }  
  17.         }  
  18.         else  
  19.             retreat(i);  
  20.     }  
  21. }  
  22.   
  23. //procedure advance(i);  
  24. void advance(i)  
  25. {  
  26.     let(i,j)be an admissible arc in A(t);  
  27.     pre(j) = i;  
  28.     i = j;  
  29. }  
  30. //procedure retreat  
  31. void retreat()  
  32. {  
  33.     d(i) = min(d(j)+1):rij>0;  
  34.     if (i != s)  
  35.         i = pre(i);  
  36. }  

 

代码模板：

 

[c-sharp] view plaincopy
  1. #include <iostream>  
  2. #include <cstring>  
  3. #include <cstdlib>  
  4.   
  5. using namespace std;  
  6.   
  7. const int Max = 225;  
  8. const int oo = 210000000;  
  9.   
  10. int n,m,c[Max][Max],r[Max][Max],c1[Max][Max],source,sink；  
  11.   
  12. //c1是c的反向网络，用于dis数组的初始化  
  13. int dis[Max];  
  14.   
  15. void initialize()// 初始化dis数组  
  16. {  
  17.      int q[Max],head = 0, tail = 0;//BFS队列  
  18.      memset(dis,0,sizeof(dis));  
  19.      q[++head] = sink;  
  20.      while(tail < head)  
  21.      {  
  22.          int u = q[++tail], v;  
  23.          for(v = 0; v <= sink; v++)  
  24.          {  
  25.                if(!dis[v] && c1[u][v] > 0)  
  26.                {  
  27.                    dis[v] = dis[u] + 1;  
  28.                    q[++head] = v;  
  29.                }  
  30.          }  
  31.      }  
  32. }  
  33.   
  34. int maxflow_sap()  
  35. {  
  36.     initialize();  
  37.     int top = source, pre[Max], flow = 0, i, j, k, low[Max];  
  38.   
  39.     // top是当前增广路中最前面一个点。  
  40.     memset(low,0,sizeof(low));//low数组用于保存路径中的最小容量  
  41.     while(dis[source] < n)  
  42.     {  
  43.         bool flag = false;  
  44.         low[source] = oo;  
  45.         for(i = 0; i <= sink; i++)//找允许弧，根据允许弧的定义  
  46.         {  
  47.             if(r[top][i] > 0 && dis[top] == dis[i] + 1)  
  48.             {  
  49.                 flag = true;  
  50.                 break;  
  51.             }  
  52.         }  
  53.         if(flag)// 如果找到允许弧  
  54.         {  
  55.             low[i] = r[top][i];  
  56.             if(low[i] > low[top]) low[i] = low[top];//更新low  
  57.             pre[i] = top; top = i;  
  58.             if(top == sink)// 如果找到一条增广路了  
  59.             {  
  60.                 flow += low[sink];  
  61.                 j = top;  
  62.                 while(j != source)// 路径回溯更新残留网络  
  63.                 {  
  64.                     k = pre[j];  
  65.                     r[k][j] -= low[sink];  
  66.                     r[j][k] += low[sink];  
  67.                     j = k;  
  68.                 }  
  69.                 top = source;//从头开始再找最短路  
  70.                 memset(low,0,sizeof(low));  
  71.             }  
  72.         }  
  73.         else// 如果没有允许弧  
  74.         {  
  75.             int mindis = 10000000;  
  76.             for(j = 0; j <= sink; j++)//找和top相邻dis最小的点  
  77.             {  
  78.                 if(r[top][j] > 0 && mindis > dis[j] + 1)  
  79.                     mindis = dis[j] + 1;  
  80.             }  
  81.             dis[top] = mindis;//更新top的距离值  
  82.             if(top != source) top = pre[top];// 回溯找另外的路  
  83.         }  
  84.     }  
  85.     return(flow);  
  86. }  

 

运用gap优化：

即当标号中出现了不连续标号的情况时，即可以证明已经不存在新的增广流，此时的流量即为最大流。

简单证明下：

假设不存在标号为k的结点，那么这时候可以将所有的结点分成两部分，一部分为d（i）>k,另一部分为d（i）<k

如此分成两份，因为性质2可知，允许路为最短的增广路，又因为不存在从>k部分到<k部分的增广流，那么有最

大流最小割定理可知此时便是最大流。

优化代码：要注意在标号的时候不能直接把所有的初始为0，而应该为-1，否则会在gap优化的时候出现问题，不满足递增的性质，切记！

 

[c-sharp] view plaincopy
  1. #include <iostream>  
  2. #include <cstring>  
  3. #include <cstdlib>  
  4.   
  5. using namespace std;  
  6.   
  7. const int Max = 225;  
  8. const int oo = 210000000;  
  9.   
  10. int n,m,c[Max][Max],r[Max][Max],c1[Max][Max],source,sink;  
  11. int gap[Max];//用gap来记录当前标号为i的个数，用于gap优化；  
  12.   
  13. //c1是c的反向网络，用于dis数组的初始化  
  14. int dis[Max];  
  15.   
  16. void initialize()// 初始化dis数组  
  17. {  
  18.     memset(gap,0,sizeof(gap));  
  19.     int q[Max],head = 0, tail = 0;//BFS队列  
  20.     memset(dis,0xff,sizeof(dis));  
  21.     q[++head] = sink;  
  22.     dis[sink] = 0;  
  23.     gap[0]++;  
  24.     while(tail < head)  
  25.     {  
  26.         int u = q[++tail], v;  
  27.         for(v = 0; v <= sink; v++)  
  28.         {  
  29.             if(!dis[v] && c1[u][v] > 0)  
  30.             {  
  31.                 dis[v] = dis[u] + 1;  
  32.                 gap[dis[v]]++;  
  33.                 q[++head] = v;  
  34.             }  
  35.         }  
  36.      }  
  37. }  
  38.   
  39. int maxflow_sap()  
  40. {  
  41.     initialize();  
  42.     int top = source, pre[Max], flow = 0, i, j, k, low[Max];  
  43.   
  44.     // top是当前增广路中最前面一个点。  
  45.     memset(low,0,sizeof(low));//low数组用于保存路径中的最小容量  
  46.     while(dis[source] < n)  
  47.     {  
  48.         bool flag = false;  
  49.         low[source] = oo;  
  50.         for(i = 0; i <= sink; i++)//找允许弧，根据允许弧的定义  
  51.         {  
  52.             if(r[top][i] > 0 && dis[top] == dis[i] + 1&&dis[i]>=0)  
  53.             {  
  54.                 flag = true;  
  55.                 break;  
  56.             }  
  57.         }  
  58.         if(flag)// 如果找到允许弧  
  59.         {  
  60.             low[i] = r[top][i];  
  61.             if(low[i] > low[top])   
  62.                 low[i] = low[top];//更新low  
  63.             pre[i] = top; top = i;  
  64.             if(top == sink)// 如果找到一条增广路了  
  65.             {  
  66.                 flow += low[sink];  
  67.                 j = top;  
  68.                 while(j != source)// 路径回溯更新残留网络  
  69.                 {  
  70.                     k = pre[j];  
  71.                     r[k][j] -= low[sink];  
  72.                     r[j][k] += low[sink];  
  73.                     j = k;  
  74.                 }  
  75.                 top = source;//从头开始再找最短路  
  76.                 memset(low,0,sizeof(low));  
  77.             }  
  78.         }  
  79.         else// 如果没有允许弧  
  80.         {  
  81.             int mindis = 10000000;  
  82.             for(j = 0; j <= sink; j++)//找和top相邻dis最小的点  
  83.             {  
  84.                 if(r[top][j] > 0 && mindis > dis[j] + 1&&dis[j]>=0)  
  85.                     mindis = dis[j] + 1;  
  86.             }  
  87.             gap[dis[top]]--;  
  88.             if (gap[dis[top]] == 0)  
  89.                 break;  
  90.             gap[mindis]++;  
  91.             dis[top] = mindis;//更新top的距离值  
  92.             if(top != source) top = pre[top];// 回溯找另外的路  
  93.         }  
  94.     }  
  95.     return(flow);  
  96. }  

 

 

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

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

 

[cpp] view plaincopy
  1. #include <stdio.h>  
  2. #include <string.h>  
  3. #define INF 2100000000  
  4. #define MAXN 301  
  5.   
  6. int SAP(int map[][MAXN],int v_count,int s,int t)      //邻接矩阵，节点总数，始点，汇点  
  7. {  
  8.     int i;  
  9.     int cur_flow,max_flow,cur,min_label,temp;         //当前流，最大流，当前节点，最小标号，临时变量  
  10.     char flag;                                        //标志当前是否有可行流  
  11.     int cur_arc[MAXN],label[MAXN],neck[MAXN];         //当前弧，标号，瓶颈边的入点（姑且这么叫吧）  
  12.     int label_count[MAXN],back_up[MAXN],pre[MAXN];    //标号为i节点的数量，cur_flow的纪录,当前流路径中前驱  
  13.   
  14.     //初始化  
  15.     memset(label,0,MAXN*sizeof(int));  
  16.     memset(label_count,0,MAXN*sizeof(int));  
  17.   
  18.     memset(cur_arc,0,MAXN*sizeof(int));  
  19.     label_count[0]=v_count;                           //全部初始化为距离为0  
  20.   
  21.     neck[s]=s;  
  22.     max_flow=0;  
  23.     cur=s;  
  24.     cur_flow=INF;  
  25.   
  26.     //循环代替递归  
  27.     while(label[s]<v_count)  
  28.     {  
  29.         back_up[cur]=cur_flow;  
  30.         flag=0;  
  31.   
  32.         //选择允许路径（此处还可用邻接表优化）  
  33.         for(i=cur_arc[cur];i<v_count;i++)    //当前弧优化  
  34.         {  
  35.            if(map[cur][i]!=0&&label[cur]==label[i]+1)    //找到允许路径  
  36.            {  
  37.                flag=1;  
  38.                cur_arc[cur]=i;    //更新当前弧  
  39.                if(map[cur][i]<cur_flow)    //更新当前流  
  40.                {  
  41.                    cur_flow=map[cur][i];  
  42.                    neck[i]=cur;     //瓶颈为当前节点  
  43.                }  
  44.                else  
  45.                {  
  46.                    neck[i]=neck[cur];     //瓶颈相对前驱节点不变  
  47.                }  
  48.                pre[i]=cur;    //记录前驱  
  49.                cur=i;  
  50.                if(i==t)    //找到可行流  
  51.                {  
  52.                    max_flow+=cur_flow;    //更新最大流  
  53.   
  54.                    //修改残量网络  
  55.                    while(cur!=s)  
  56.                    {  
  57.                        if(map[pre[cur]][cur]!=INF)map[pre[cur]][cur]-=cur_flow;  
  58.                        back_up[cur] -= cur_flow;  
  59.                        if(map[cur][pre[cur]]!=INF)map[cur][pre[cur]]+=cur_flow;  
  60.                        cur=pre[cur];  
  61.                    }  
  62.   
  63.                    //优化，瓶颈之后的节点出栈  
  64.                    cur=neck[t];  
  65.                    cur_flow=back_up[cur];   
  66.                }  
  67.                break;  
  68.            }  
  69.         }  
  70.         if(flag)continue;  
  71.   
  72.         min_label=v_count-1;    //初始化min_label为节点总数-1  
  73.   
  74.         //找到相邻的标号最小的节点     
  75.         for(i=0;i<v_count;i++)  
  76.         {  
  77.             if(map[cur][i]!=0&&label[i]<min_label)  
  78.             {  
  79.                 min_label=label[i];  
  80.                 temp=i;  
  81.             }  
  82.         }  
  83.         cur_arc[cur]=temp;    //记录当前弧，下次从提供最小标号的节点开始搜索  
  84.         label_count[label[cur]]--;    //修改标号纪录  
  85.         if(label_count[label[cur]]==0)break;    //GAP优化  
  86.         label[cur]=min_label+1;    //修改当前节点标号  
  87.         label_count[label[cur]]++;     //修改标号记录  
  88.         if(cur!=s)  
  89.         {  
  90.            //从栈中弹出一个节点  
  91.            cur=pre[cur];  
  92.            cur_flow=back_up[cur];  
  93.         }  
  94.     }  
  95.     return(max_flow);  
  96. }  

 

posted @ 2012-04-11 19:03  'wind  阅读(343)  评论(0)  编辑  收藏  举报