想法
如何判定在当前流量下，一个费用流是否是最小费用流？
这个问题等同于，是否存在一种方案，在不改变总流量的情况下改变一些边的流量，最终减小费用。网络流的消圈定理就对此作出了解答。
消圈定理
可行流 (f) 是当前流量下的最小费用流 (Leftrightarrow) 残余网络不存在负环
把单点看成一个长度为 0 的环。
设 (f') 是一个流量不变，费用更小的可行流，那么 (f'-f) 的每个点都满足流入流量都等于流出流量。这是因为 (f,f') 的除源点汇点的其他点都分别满足这个条件，而总流量不变，所以源汇的流量也平衡。
这时 (f'-f) 是由一些环组成的。若不存在环，那么不可能满足流量平衡；若除了环还存在一些链，那么链头和链尾也就不满足要求。
由于 (w(f')<w(f)) ，那么 (w(f'-f)<0) ，所以这些环中一定存在一个负环。
应用
有了这个定理，我们就能判定当前流量下费用流是否是最小费用流了——若有负环就不是，否则是。
这也给我们提供了一个得到最小费用流的方法，即不断在负环上增广，直到找不到负环 。
直接做的复杂度上界为 (O(nm^2cw)) ，其中 (c,w) 分别为最大容量和最大费用，整个过程可能执行 (mcu) 次，一次最坏找到负环为 (O(nm)) 。
代码
POJ2175 学会了这个方法就很简单了。
用spfa来判断负环，用队列实现 600ms 左右，换成栈变成 300ms 左右。然而dfs TLE 了。讲道理，栈不就是dfs吗？！
#include<cstdio>
#include<cctype>
#include<climits>
#include<cstring>
#include<algorithm>
using namespace std;
inline char nchar() {
static const int bufl=1<<20;
static char buf[bufl],*a,*b;
return a==b && (b=(a=buf)+fread(buf,1,bufl,stdin),a==b)?EOF:*a++;
}
inline int read() {
int x=0,f=1;
char c=nchar();
for (;!isdigit(c);c=nchar()) if (c=='-') f=-1;
for (;isdigit(c);c=nchar()) x=x*10+c-'0';
return x*f;
}
const int maxn=205;
const int maxe=maxn*maxn*2;
int n,m,bt[maxn],ct[maxn];
struct D {
int x,y,d;
inline int operator - (const D &b) const {
return abs(x-b.x)+abs(y-b.y)+1;
}
} b[maxn],c[maxn];
namespace graph {
struct edge {
int u,v,c,w,nxt;
} e[maxe];
int h[maxn],tot=1,pre[maxn],d[maxn],cnt[maxn],que[maxe],ql,qr; // pre e
bool inq[maxn],vis[maxn];
inline int add(int x,int y,int c,int w) {
e[++tot]=(edge){x,y,c,w,h[x]};
return h[x]=tot;
}
void print(int x) {
for (int i=x;;i=e[pre[i]].u) {
const int &k=pre[i];
--e[k].c,++e[k^1].c;
if (e[pre[i]].u==x) break;
}
puts("SUBOPTIMAL");
for (int i=1;i<=n;++i) {
static int a[maxn];
for (int j=h[i],v=e[j].v;j;j=e[j].nxt,v=e[j].v) a[v-n]=e[j^1].c;
for (int j=1;j<=m;++j) printf("%d%c",a[j],"
"[j==m]);
}
}
inline int find(int x) {
for (int i=x;;i=e[pre[i]].u) if (!vis[i]) vis[i]=true; else return i;
}
inline void cancle(int s,int t) {
memset(d,0x3f,sizeof d);
d[t]=0;
que[ql=qr=1]=t;
while (ql<=qr) {
int x=que[qr--];
for (int i=h[x],v=e[i].v;i;i=e[i].nxt,v=e[i].v) if (e[i].c && d[v]>d[x]+e[i].w) {
d[v]=d[x]+e[i].w;
pre[v]=i;
if (inq[v]) continue;
inq[v]=true;
if (++cnt[que[++qr]=v]>t+1) {
int p=find(v);
print(p);
exit(0);
}
}
inq[x]=false;
}
}
}
using graph::add;
int main() {
#ifndef ONLINE_JUDGE
freopen("test.in","r",stdin);
#endif
n=read(),m=read();
int S=0,T=n+m+1;
for (int i=1;i<=n;++i) {
b[i].x=read(),b[i].y=read();
bt[i]=add(S,i,b[i].d=read(),0);
add(i,S,0,0);
}
for (int i=1;i<=m;++i) {
c[i].x=read(),c[i].y=read();
ct[i]=add(n+i,T,c[i].d=read(),0);
add(T,n+i,0,0);
}
for (int i=1;i<=n;++i) for (int j=1;j<=m;++j) {
int d=b[i]-c[j],x=read();
add(i,n+j,INT_MAX-x,d);
add(n+j,i,x,-d);
graph::e[bt[i]].c-=x,graph::e[bt[i]^1].c+=x;
graph::e[ct[j]].c-=x,graph::e[ct[j]^1].c+=x;
}
graph::cancle(S,T);
puts("OPTIMAL");
return 0;
}