function [A, b, c]=transport(a, d, C)
m=length(a)
n=length(d)
if sum(a)~=sum(d) return; end;
A=zeros(m+n, m*n);
c=zeros(1,m*n);
for j=1:m*n
    s=ceil(j/n);
    u=mod(j,n);
    if u==0, u=n; end;
    A(s,j)=1;
    A(m+u,j)=1;
    c(j)=C(s,u);
end;
for i=1:m
    b(i)=a(i);
end;
for j=1:n
    b(m+j)=d(j);
end;
