费用流 Dijkstra 原始对偶方法(primal-dual method)

简单叙述用Dijkstra求费用流

Dijkstra不能求有负权边的最短路。
类似于Johnson算法，我们也可以设计一个势函数，以满足在与原图等价的新图中的边权非负。

但是这个算法并不能处理有负圈的情况（可能需要消圈算法）。

对网络(G)中的每一个点设置一个势函数(h(u))，在任意残留网络G'的任意边((u, v))都需要满足(w_{u, v} + h(u) - h(v) ge 0)。
令图G的对偶图（不知道能不能这么说）为(G')，其对应的边((u, v))的权值为(w'_{u, v} = w_{u, v} + h(u) - h(v))
对于原图中的任意一条路径((u_1, u_2, cdots, u_k))，它在(G)中的权值为(w_{u_1, u_2} + w_{u_2, u_3} + cdots + w_{u_{k - 1}, u_k})，在(G')中的权值为
(w'_{u_1, u_2} + w'_{u_2, u_3} + cdots + w'_{u_{k - 1}, u_k})
(= w_{u_1, u_2} + w_{u_2, u_3} + cdots + w_{u_{k - 1}, u_k} + h(u_1) - h(u_2) + h(u_2) - h(u_3) + cdots + h(u_{k - 1}) - h(u_k))
(= w_{u_1, u_2} + w_{u_2, u_3} + cdots + w_{u_{k - 1}, u_k} + h(u_1) - h(u_k))

所以，我们在(G')求出的路径都可以对应到(G)上，令(dist_{u, v})为图(G)点(u)到点(v)的最短路径，(dist'_{u, v})为图(G')点(u)到点(v)的最短路径，显然有(dist_{u, v} = dist'_{u, v} - h(u) + h(v))
所以我们只需要求(G')的最短路径，就能对应回原图的最短路径。

若网络(G)初始边权非负，我们可令(h(u) = 0)
否则我们令(h(u) = dist_{s, u})，这个可以用Bellman-Ford算法解决。（这与Johnson算法是一模一样的）

我们考虑怎么维护势函数(h(u))。
令网路(G)上(dist_{s, u} = d_u)，网路(G')上(dist'_{s, u} = d'_u)
令残余网路(G')上新的势函数为(h'(u))。
对于残余网络上的一条边((u, v))，有两种可能：

((u, v) in G)，那么有(d_u + w_{u, v} + h(u) - h(v) ge d_v)
移项得(w_{u, v} + (h(u) + d_u) - (h(v) + d_v) ge 0)

((u, v) in G)，那么((v, u) in G的增广路)，就有(d_v + w_{v, u} + h(v) - h(u) = d_u)
移项得(-w_{v, u} + (h(u) + d_u) - (h(v) - d_v) = 0)
由(w_{u, v} = -w_{v, u})可得(w_{u, v} + (h(u) + d_u) - (h(v) - d_v) = 0)
也就有(w_{u, v} + (h(u) + d_u) - (h(v) - d_v) ge 0)

所以我们不妨令(h'(u) = h(u) + d_u)，这就维护好了势函数(h(u))

luogu P3381 【模板】最小费用最大流

单路增广

#include <bits/stdc++.h>

using namespace std;

typedef pair<int, int> PII;

#define fi first
#define se second
#define mp make_pair

const int N = 5005, M = 50005;
const int INF = ~0u >> 1;

int n, m, s, t;
int tot = -1, head[N];
struct Edge {
    int p, nxt, c, w;
    Edge(int p = 0, int nxt = 0, int c = 0, int w = 0) : p(p), nxt(nxt), c(c), w(w) {}
} edge[M * 2];
inline void Add_Edge(int u, int v, int c, int w) {
    edge[++tot] = Edge(v, head[u], c, w);
    head[u] = tot;
    return;
}

int h[N], d[N], dc[N], pr[N];
priority_queue<PII> pq;

bool Dijkstra(int s, int t, int &mf, int &mc) {
    fill(d + 1, d + n + 1, INF);
    dc[s] = INF;
    d[s] = 0;
    pr[s] = -1;
    pq.push(mp(0, s));
    while (!pq.empty()) {
        PII cur = pq.top();
        pq.pop();
        int u = cur.se;
        if (-cur.fi > d[u]) continue;
        for (int i = head[u]; ~i; i = edge[i].nxt) {
            int v = edge[i].p, c = edge[i].c, w = edge[i].w + h[u] - h[v];
            if (!c) continue;
            if (d[v] > d[u] + w) {
                d[v] = d[u] + w;
                pr[v] = i;
                dc[v] = min(dc[u], c);
                pq.push(mp(-d[v], v));
            }
        }
    }
    if (d[t] == INF) return 0;
    for (int i = 1; i <= n; ++i)
        if (d[i] < INF) h[i] += d[i];
    int c = dc[t];
    mf += c;
    mc += c * h[t];
    for (int x = t; ~pr[x]; x = edge[pr[x] ^ 1].p) {
        edge[pr[x]].c -= c;
        edge[pr[x] ^ 1].c += c;
    }
    return 1;
}

int main() {
    scanf("%d%d%d%d", &n, &m, &s, &t);
    memset(head, -1, sizeof(head));
    for (int i = 1; i <= m; ++i) {
        int u, v, c, w;
        scanf("%d%d%d%d", &u, &v, &c, &w);
        Add_Edge(u, v, c, w);
        Add_Edge(v, u, 0, -w);
    }
    int mf = 0, mc = 0;
    while (Dijkstra(s, t, mf, mc));
    printf("%d %d
", mf, mc);
    return 0;
}

多路增广（复杂度不变）

#include <bits/stdc++.h>

using namespace std;

typedef pair<int, int> PII;

#define fi first
#define se second
#define mp make_pair

const int N = 5005, M = 50005;
const int INF = ~0u >> 1;

int n, m, s, t;
int tot = -1, head[N], cur[N];
struct Edge {
    int p, nxt, c, w;
    Edge(int p = 0, int nxt = 0, int c = 0, int w = 0) : p(p), nxt(nxt), c(c), w(w) {}
} edge[M * 2];
inline void Add_Edge(int u, int v, int c, int w) {
    edge[++tot] = Edge(v, head[u], c, w);
    head[u] = tot;
    return;
}

int h[N], d[N];
priority_queue<PII> pq;

bool Dijkstra(int s, int t) {
    fill(d + 1, d + n + 1, INF);
    pq.push(mp(d[s] = 0, s));
    while (!pq.empty()) {
        PII cur = pq.top();
        pq.pop();
        int u = cur.se;
        if (-cur.fi > d[u]) continue;
        for (int i = head[u]; ~i; i = edge[i].nxt) {
            int v = edge[i].p, c = edge[i].c, w = edge[i].w + h[u] - h[v];
            if (c && d[v] > d[u] + w) pq.push(mp(-(d[v] = d[u] + w), v));
        }
    }
    return d[t] < INF;
}

int v[N];

int DFS(int u, int c, int t) {
    if (u == t) return c;
    int r = c;
    v[u] = 1;
    for (int &i = cur[u]; ~i && r; i = edge[i].nxt) {
        int v = edge[i].p, c = edge[i].c, w = edge[i].w + h[u] - h[v];
        if (!::v[v] && c && d[u] + w == d[v]) {
            int x = DFS(v, min(r, c), t);
            r -= x;
            edge[i].c -= x;
            edge[i ^ 1].c += x;
        }
    }
    v[u] = 0;
    return c - r;
}

int main() {
    scanf("%d%d%d%d", &n, &m, &s, &t);
    memset(head, -1, sizeof(head));
    for (int i = 1; i <= m; ++i) {
        int u, v, c, w;
        scanf("%d%d%d%d", &u, &v, &c, &w);
        Add_Edge(u, v, c, w);
        Add_Edge(v, u, 0, -w);
    }
    int mf = 0, mc = 0;
    while (Dijkstra(s, t)) {
        memcpy(cur, head, sizeof(cur));
        int c = DFS(s, INF, t);
        for (int i = 1; i <= n; ++i)
            if (d[i] < INF) h[i] += d[i];
        mf += c;
        mc += c * h[t];
    }
    printf("%d %d
", mf, mc);
    return 0;
}