[模板] 最大流和费用流分别的两种做法
注意：以下所有说明均以帮助理解模板为目的，不保证正确性。
最大流
dinic
考虑每次找一条S到T的不满流的路径并进行增广，但需要解决转圈圈的问题
所以首先用bfs给你的网络分层，之后再按照层(每个点只能走到下一层的点)来dfs，尽可能地把流量都占满
当前弧优化：在一次dfs中，每个点已经被访问完的出边再被访问没有意义，所以可以记录下来已经做到了哪个边
LOJ101
1 #include<bits/stdc++.h> 2 #define pa pair<int,int> 3 #define CLR(a,x) memset(a,x,sizeof(a)) 4 #define MP make_pair 5 #define fi first 6 #define se second 7 using namespace std; 8 typedef long long ll; 9 typedef unsigned long long ull; 10 typedef unsigned int ui; 11 typedef long double ld; 12 const int maxn=105,maxm=5005; 13 14 inline char gc(){ 15 return getchar(); 16 static const int maxs=1<<16;static char buf[maxs],*p1=buf,*p2=buf; 17 return p1==p2&&(p2=(p1=buf)+fread(buf,1,maxs,stdin),p1==p2)?EOF:*p1++; 18 } 19 inline ll rd(){ 20 ll x=0;char c=gc();bool neg=0; 21 while(c<'0'||c>'9'){if(c=='-') neg=1;c=gc();} 22 while(c>='0'&&c<='9') x=(x<<1)+(x<<3)+c-'0',c=gc(); 23 return neg?(~x+1):x; 24 } 25 26 struct Edge{ 27 int b,l,ne; 28 }eg[maxm*2]; 29 int egh[maxn],ect=1,cur[maxn]; 30 int N,M,S,T; 31 int dep[maxn]; 32 33 inline void adeg(int a,int b,int l){ 34 eg[++ect].b=b,eg[ect].l=l,eg[ect].ne=egh[a],egh[a]=ect; 35 } 36 37 inline bool bfs(){ 38 static queue<int> q; 39 for(int i=1;i<=N;i++) dep[i]=1e9; 40 dep[S]=0;q.push(S); 41 while(!q.empty()){ 42 int p=q.front();q.pop(); 43 for(int i=egh[p];i;i=eg[i].ne){ 44 int b=eg[i].b;if(!eg[i].l||dep[b]<=dep[p]+1) continue; 45 dep[b]=dep[p]+1;q.push(b); 46 } 47 } 48 return dep[T]<1e9; 49 } 50 51 ll dinic(int x,ll y){ 52 if(x==T) return y; 53 ll tmp=y; 54 for(int &i=cur[x];i;i=eg[i].ne){ 55 int b=eg[i].b;if(!eg[i].l||dep[b]!=dep[x]+1) continue; 56 ll r=dinic(b,min(tmp,(ll)eg[i].l)); 57 tmp-=r,eg[i].l-=r,eg[i^1].l+=r; 58 if(!tmp) break; 59 }return y-tmp; 60 } 61 62 int main(){ 63 //freopen("","r",stdin); 64 N=rd(),M=rd(),S=rd(),T=rd(); 65 for(int i=1;i<=M;i++){ 66 int a=rd(),b=rd(),l=rd(); 67 adeg(a,b,l),adeg(b,a,0); 68 } 69 ll ans=0; 70 while(bfs()){ 71 memcpy(cur,egh,sizeof(cur)); 72 ans+=dinic(S,1e15); 73 } 74 printf("%lld\n",ans); 75 return 0; 76 }
预流推进(HLPP)
考虑对某个点，尽可能地把它目前有的流量都推出去，当然这会导致流量不平衡，也就是说，每个点都会有一个额外流(我记做exc)
假设S的exc是inf，那么当除S和T外所有点的exc都为0时，T的exc的最大值就是最大流量
我们考虑给每个点附上一个高度h，并钦定只能从h推向h-1，以及钦定S的高度是N，T的高度是0不会改变
于是就可以每次取出最高的有exc的高度，然后推流，之后如果它的exc还有剩，就把它的h改成它出边中最小的h+1
所以用一个优先队列维护一下就好了
优化1：先用一次T到S的bfs，给每个点附上一个初始的h
优化2：如果某个高度x没有点了，那么所有高度大于x，小于N+1的点都不会再有机会流向T，所以直接把它们的高度设成N+1，所以拿一个cnt记一下每个高度的点的个数即可。记得cnt要开两倍大小
LOJ101
1 #include<bits/stdc++.h> 2 #define pa pair<int,int> 3 #define CLR(a,x) memset(a,x,sizeof(a)) 4 #define MP make_pair 5 #define fi first 6 #define se second 7 using namespace std; 8 typedef long long ll; 9 typedef unsigned long long ull; 10 typedef unsigned int ui; 11 typedef long double ld; 12 const int maxn=105,maxm=5005; 13 14 inline char gc(){ 15 return getchar(); 16 static const int maxs=1<<16;static char buf[maxs],*p1=buf,*p2=buf; 17 return p1==p2&&(p2=(p1=buf)+fread(buf,1,maxs,stdin),p1==p2)?EOF:*p1++; 18 } 19 inline ll rd(){ 20 ll x=0;char c=gc();bool neg=0; 21 while(c<'0'||c>'9'){if(c=='-') neg=1;c=gc();} 22 while(c>='0'&&c<='9') x=(x<<1)+(x<<3)+c-'0',c=gc(); 23 return neg?(~x+1):x; 24 } 25 26 struct Edge{ 27 int b,l,ne; 28 }eg[maxm*2]; 29 int egh[maxn],ect=1; 30 int N,M,S,T; 31 int h[maxn],cnt[maxn*2]; 32 ll exc[maxn]; 33 bool flag[maxn]; 34 struct cmp{ 35 inline bool operator ()(int a,int b) const{return h[a]<h[b];} 36 }; 37 priority_queue<int,vector<int>,cmp> q; 38 39 inline void adeg(int a,int b,int l){ 40 eg[++ect].b=b,eg[ect].l=l,eg[ect].ne=egh[a],egh[a]=ect; 41 } 42 43 inline bool bfs(){ 44 for(int i=1;i<=N;i++) h[i]=1e9; 45 queue<int> q; 46 q.push(T);h[T]=0; 47 while(!q.empty()){ 48 int p=q.front();q.pop(); 49 for(int i=egh[p];i;i=eg[i].ne){ 50 int b=eg[i].b;if(!eg[i^1].l||h[b]<=h[p]+1) continue; 51 h[b]=h[p]+1,q.push(b); 52 } 53 }return h[S]<1e9; 54 } 55 56 inline ll HLPP(){ 57 if(!bfs()) return 0; 58 h[S]=N;flag[S]=flag[T]=1; 59 for(int i=1;i<=N;i++) if(h[i]<1e9) cnt[h[i]]++; 60 for(int i=egh[S];i;i=eg[i].ne){ 61 int b=eg[i].b,l=eg[i].l;if(!l) continue; 62 exc[b]+=l,eg[i].l-=l,eg[i^1].l+=l; 63 if(!flag[b]) q.push(b),flag[b]=1; 64 } 65 while(!q.empty()){ 66 int p=q.top();q.pop(); 67 flag[p]=0; 68 for(int i=egh[p];i;i=eg[i].ne){ 69 int b=eg[i].b;if(!eg[i].l||h[b]!=h[p]-1) continue; 70 int l=min(exc[p],(ll)eg[i].l); 71 exc[p]-=l,exc[b]+=l,eg[i].l-=l,eg[i^1].l+=l; 72 if(!flag[b]) q.push(b),flag[b]=1; 73 if(!exc[p]) break; 74 } 75 if(exc[p]){ 76 if(!--cnt[h[p]]){ 77 for(int i=1;i<=N;i++){ 78 if(!flag[i]&&h[i]>h[p]&&h[i]<N+1) h[i]=N+1; 79 } 80 } 81 int mi=1e9; 82 for(int i=egh[p];i;i=eg[i].ne){ 83 int b=eg[i].b;if(!eg[i].l) continue; 84 mi=min(mi,h[b]+1); 85 } 86 cnt[h[p]=mi]++; 87 q.push(p),flag[p]=1; 88 } 89 } 90 return exc[T]; 91 } 92 93 int main(){ 94 //freopen("","r",stdin); 95 N=rd(),M=rd(),S=rd(),T=rd(); 96 for(int i=1;i<=M;i++){ 97 int a=rd(),b=rd(),l=rd(); 98 adeg(a,b,l),adeg(b,a,0); 99 } 100 ll ans=HLPP(); 101 printf("%lld\n",ans); 102 return 0; 103 }
费用流
spfa
考虑使用spfa来找到一条从S到T的费用和最小的路径，然后增广这条路径。记一下每个点最短距离取的那个入边即可
LOJ102
1 #include<bits/stdc++.h> 2 #define pa pair<int,int> 3 #define CLR(a,x) memset(a,x,sizeof(a)) 4 #define MP make_pair 5 #define fi first 6 #define se second 7 using namespace std; 8 typedef long long ll; 9 typedef unsigned long long ull; 10 typedef unsigned int ui; 11 typedef long double ld; 12 const int maxn=405,maxm=15005,inf=(1u<<31)-1; 13 14 inline char gc(){ 15 return getchar(); 16 static const int maxs=1<<16;static char buf[maxs],*p1=buf,*p2=buf; 17 return p1==p2&&(p2=(p1=buf)+fread(buf,1,maxs,stdin),p1==p2)?EOF:*p1++; 18 } 19 inline ll rd(){ 20 ll x=0;char c=gc();bool neg=0; 21 while(c<'0'||c>'9'){if(c=='-') neg=1;c=gc();} 22 while(c>='0'&&c<='9') x=(x<<1)+(x<<3)+c-'0',c=gc(); 23 return neg?(~x+1):x; 24 } 25 26 struct Edge{ 27 int b,l,c,ne; 28 }eg[maxm*2]; 29 int N,M,egh[maxn],ect=1,S,T; 30 int dis[maxn],ine[maxn]; 31 bool flag[maxn]; 32 queue<int> q; 33 34 inline void adeg(int a,int b,int c,int l){ 35 eg[++ect]=Edge{b,l,c,egh[a]},egh[a]=ect; 36 } 37 38 inline bool spfa(){ 39 for(int i=1;i<=N;i++) dis[i]=inf,flag[i]=0; 40 q.push(S);dis[S]=0; 41 while(!q.empty()){ 42 int p=q.front();q.pop(); 43 flag[p]=0; 44 for(int i=egh[p];i;i=eg[i].ne){ 45 int b=eg[i].b;if(!eg[i].l||dis[b]<=dis[p]+eg[i].c) continue; 46 dis[b]=dis[p]+eg[i].c,ine[b]=i; 47 if(!flag[b]) q.push(b),flag[b]=1; 48 } 49 }return dis[T]<inf; 50 } 51 52 int main(){ 53 //freopen("","r",stdin); 54 N=rd(),M=rd(),S=1,T=N; 55 for(int i=1;i<=M;i++){ 56 int a=rd(),b=rd(),l=rd(),c=rd(); 57 adeg(a,b,c,l),adeg(b,a,-c,0); 58 } 59 int flw=0,cst=0; 60 while(spfa()){ 61 int mi=inf; 62 for(int x=T;x!=S;x=eg[ine[x]^1].b) mi=min(mi,eg[ine[x]].l); 63 for(int x=T;x!=S;x=eg[ine[x]^1].b) eg[ine[x]].l-=mi,eg[ine[x]^1].l+=mi; 64 flw+=mi,cst+=mi*dis[T]; 65 } 66 printf("%d %d\n",flw,cst); 67 return 0; 68 }
原始对偶
首先，我们刚刚的spfa实际上已经像dinic的bfs一样给图分了个层，所以我们可以仿照dinic的做法，按照这些层来选增广路（就是只能走最短路）。但要注意因为可能有环，已经走过的点就不能再走了，要打个标记，也正是因为这个标记，一次增广可能并不完全，所以我需要做多次增广直到找不到增广路。而且如果spfa是从S到T求的的话，要从T到S来增广，因为你只知道从S到x的距离，你再从S开始去增广的话，它有可能根本到不了T。
然而这还不够，我们想试着用更高效的方法来分层。
考虑dijkstra，然而图中有负边。但我们可以先进行一次spfa，之后每次都根据上次算出来的dis改变边权，使得它们（还有剩余流量的）全是非负的。
考虑连接i和j的边(i,j)，一定有dis[i]+cost[(i,j)]>=dis[j]，所以我们可以把(i,j)的边权变成dis[i]+cost[(i,j)]-dis[j]。之后我们会发现，S到T的路径的和就变成了原来的和加上dis[S]减去dis[T]。于是就可以dijkstra了
但要注意的是，每次dijkstra完都要更新边权，因为会出现新边（就是说 原来没有剩余流量的边现在有剩余流量了，而因为它原来没有剩余流量 所以它的边权并不一定是非负的）
LOJ102
1 #include<bits/stdc++.h> 2 #define pa pair<int,int> 3 #define CLR(a,x) memset(a,x,sizeof(a)) 4 #define MP make_pair 5 #define fi first 6 #define se second 7 using namespace std; 8 typedef long long ll; 9 typedef unsigned long long ull; 10 typedef unsigned int ui; 11 typedef long double ld; 12 const int maxn=405,maxm=15005,inf=(1u<<31)-1; 13 14 inline char gc(){ 15 return getchar(); 16 static const int maxs=1<<16;static char buf[maxs],*p1=buf,*p2=buf; 17 return p1==p2&&(p2=(p1=buf)+fread(buf,1,maxs,stdin),p1==p2)?EOF:*p1++; 18 } 19 inline ll rd(){ 20 ll x=0;char c=gc();bool neg=0; 21 while(c<'0'||c>'9'){if(c=='-') neg=1;c=gc();} 22 while(c>='0'&&c<='9') x=(x<<1)+(x<<3)+c-'0',c=gc(); 23 return neg?(~x+1):x; 24 } 25 26 struct Edge{ 27 int b,l,c,ne; 28 }eg[maxm*2]; 29 int N,M,egh[maxn],ect=1,S,T; 30 int dis[maxn],off; 31 bool flag[maxn]; 32 33 inline void adeg(int a,int b,int c,int l){ 34 eg[++ect]=Edge{b,l,c,egh[a]},egh[a]=ect; 35 } 36 37 inline bool spfa(){ 38 static queue<int> q; 39 for(int i=1;i<=N;i++) dis[i]=inf,flag[i]=0; 40 q.push(S);dis[S]=0; 41 while(!q.empty()){ 42 int p=q.front();q.pop(); 43 flag[p]=0; 44 for(int i=egh[p];i;i=eg[i].ne){ 45 int b=eg[i].b;if(!eg[i].l||dis[b]<=dis[p]+eg[i].c) continue; 46 dis[b]=dis[p]+eg[i].c; 47 if(!flag[b]) q.push(b),flag[b]=1; 48 } 49 } 50 for(int i=2;i<=ect;i++){ 51 eg[i].c+=dis[eg[i^1].b]-dis[eg[i].b]; 52 } 53 off+=dis[T]-dis[S]; 54 return dis[T]<inf; 55 } 56 57 inline bool dijk(){ 58 static priority_queue<pa,vector<pa>,greater<pa> > q; 59 for(int i=1;i<=N;i++) dis[i]=inf,flag[i]=0; 60 q.push(MP(0,S)),dis[S]=0; 61 while(!q.empty()){ 62 int p=q.top().se;q.pop(); 63 if(flag[p]) continue; 64 flag[p]=1; 65 for(int i=egh[p];i;i=eg[i].ne){ 66 int b=eg[i].b;if(!eg[i].l||dis[b]<=dis[p]+eg[i].c) continue; 67 dis[b]=dis[p]+eg[i].c; 68 q.push(MP(dis[b],b)); 69 } 70 } 71 for(int i=2;i<=ect;i++){ 72 eg[i].c+=dis[eg[i^1].b]-dis[eg[i].b]; 73 } 74 off+=dis[T]-dis[S]; 75 return dis[T]<inf; 76 } 77 78 int dinic(int x,int y){ 79 if(x==S) return y; 80 int tmp=y;flag[x]=1; 81 for(int i=egh[x];i;i=eg[i].ne){ 82 int b=eg[i].b;if(flag[b]||!eg[i^1].l||eg[i^1].c) continue; 83 int r=dinic(b,min(tmp,eg[i^1].l)); 84 tmp-=r,eg[i^1].l-=r,eg[i].l+=r; 85 if(!tmp) break; 86 }return y-tmp; 87 } 88 89 int main(){ 90 //freopen("","r",stdin); 91 N=rd(),M=rd(),S=1,T=N; 92 for(int i=1;i<=M;i++){ 93 int a=rd(),b=rd(),l=rd(),c=rd(); 94 adeg(a,b,c,l),adeg(b,a,-c,0); 95 } 96 int flw=0,cst=0; 97 spfa(); 98 while(dijk()){ 99 int a=-1; 100 while(a){ 101 CLR(flag,0); 102 a=dinic(T,inf); 103 flw+=a,cst+=a*off; 104 } 105 } 106 printf("%d %d\n",flw,cst); 107 return 0; 108 }
