题解 P3381 【【模板】最小费用最大流】
panda_2134
2018-02-28 22:49:47
感谢 @WuHongxun 神牛介绍原始对偶算法
按照[TopCoder上的资料](https://www.topcoder.com/community/data-science/data-science-tutorials/minimum-cost-flow-part-two-algorithms/)，楼下的原始对偶并不是很规范呢，因为只有SPFA= =
正规的原始对偶是先用Bellman-Ford算法来构造距离顶标，再进行Reduce操作。可以证明，在Reduce操作后的每条边费用一定非负。
那么就可以愉快的Dijkstra啦，我们每次进行“Dijkstra-Reduce-多路增广”的循环，直到无法增广为止。
附送一个各类费用流算法性能测试表格：其中的“Binary Heap with decrease_key”是手写的，其他的都是STL / pb_ds。
| Algorithm | Accepted | Time |
| ---------------------------------------- | -------- | ------ |
| Dijkstra+Pairing Heap+Primal Dual(O2) | Yes | 820ms |
| Dijkstra+std::priority_queue+Primal Dual(O2) | Yes | 832ms |
| Dijkstra+Binary Heap with decrease_key<br />+Primal Dual(O2) | Yes | 888ms |
| Dijkstra+Pairing Heap+Primal Dual | Yes | 1236ms |
| Dijkstra+Binary Heap with decrease_key+Primal Dual | Yes | 1528ms |
| SPFA+Primal Dual(O2) | Yes | 1548ms |
| SPFA+Edmond Karp | Yes | 1596ms |
| SPFA+Primal Dual | Yes | 2184ms |
| Dijkstra+std::priority_queue+Primal Dual | No | 3036ms |
| SPFA+SLF+Primal Dual(O2) | No | 3204ms |
| SPFA+SLF+Primal Dual | No | 4740ms |
Dijkstra + STL堆版本的原始对偶在开O2的情况下是很优秀的，当然，如果用pb_ds的配对堆，开不开O2速度都很不错。手写二叉堆性能并不比STL好多少。
**UPD-20180321:原来的代码有点问题，把改过的发上来**
```cpp
// luogu-judger-enable-o2
#include <bits/stdc++.h>
#include <bits/extc++.h>
using namespace std;
typedef long long int64;
typedef pair<int64, int> HeapNode;
typedef __gnu_pbds::priority_queue<HeapNode, greater<HeapNode>, __gnu_pbds::pairing_heap_tag>
pairing_heap;
struct Edge {
int u, v;
int64 flow, cap, cost;
int next;
};
const int MAXN = 5000, MAXM = 50000;
const int64 LLINF = 0x3f3f3f3f3f3f3f3fLL;
int e_ptr = 1, s, t, head[MAXN+10]; Edge E[(MAXM+10)<<1];
int64 dist[MAXN+10], MaxFlow, MinCost, delta;
int inq[MAXN+10], done[MAXN+10], vis[MAXN+10];
pairing_heap pq; pairing_heap :: point_iterator it[MAXN+10];
void AddEdge(int u, int v, int64 cap, int64 cost) {
E[++e_ptr] = (Edge) { u, v, 0, cap, cost, head[u] }; head[u] = e_ptr;
E[++e_ptr] = (Edge) { v, u, 0, 0, -cost, head[v] }; head[v] = e_ptr;
}
void Reduce() {
for(int i = 2; i <= e_ptr; i++)
E[i].cost += dist[E[i].v] - dist[E[i].u];
delta += dist[s];
}
bool BellmanFord() {
queue<int> Q;
memset(dist, 0x3f, sizeof(dist));
dist[t] = 0; Q.push(t); inq[t] = true;
while(!Q.empty()) {
int u = Q.front(); Q.pop(); inq[u] = false;
for(int j=head[u]; j; j=E[j].next) {
int v = E[j].v; int64 f = E[j^1].flow, c = E[j^1].cap, len = E[j^1].cost;
if(f < c && dist[v] > dist[u] + len) {
dist[v] = dist[u] + len;
if(!inq[v]) {
inq[v] = true;
Q.push(v);
}
}
}
}
return dist[s] != LLINF;
}
bool Dijkstra() {
memset(dist, 0x3f, sizeof(dist));
memset(it, 0, sizeof(it));
dist[t] = 0; it[t] = pq.push(make_pair(dist[t], t));
while(!pq.empty()) {
HeapNode p = pq.top(); pq.pop();
int u = p.second;
for(int j=head[u]; j; j=E[j].next) {
int v = E[j].v; int64 f = E[j^1].flow, c = E[j^1].cap, len = E[j^1].cost;
if(f < c && dist[v] > dist[u] + len) {
dist[v] = dist[u] + len;
if(it[v] == NULL)
it[v] = pq.push(make_pair(dist[v], v));
else
pq.modify(it[v], make_pair(dist[v], v));
}
}
}
return dist[s] != LLINF;
}
int64 DFS(int u, int64 flow) {
if(u == t || flow == 0) return flow;
vis[u] = true; // differ from dinic
int64 res = flow;
for(int j=head[u]; j; j=E[j].next) {
int v = E[j].v; int64 f = E[j].flow, c = E[j].cap, len = E[j].cost;
if(!vis[v] && f < c && len == 0) { // not `dist[v] == dist[u]` ! they do not equal !
int64 tmp = DFS(v, min(res, c-f)); // len = 0 <=> on the shortest path
E[j].flow += tmp;
E[j^1].flow -= tmp;
res -= tmp;
}
}
return flow - res;
}
void Augment() {
int64 CurFlow = 0;
while(memset(vis, 0, sizeof(vis)),
(CurFlow = DFS(s, LLINF))) {
MaxFlow += CurFlow;
MinCost += CurFlow * delta;
}
}
void PrimalDual() {
if(!BellmanFord()) return;
Reduce(); Augment();
while(Dijkstra()) {
Reduce(); Augment();
}
}
int main() {
int n, m, u, v, cap, cost;
scanf("%d%d%d%d", &n, &m, &s, &t);
for(int i=1; i<=m; i++) {
scanf("%d%d%d%d", &u, &v, &cap, &cost);
AddEdge(u, v, cap, cost);
}
PrimalDual();
printf("%lld %lld", MaxFlow, MinCost);
return 0;
}
```