% capacity(M,z,d) expects a network's incidence matrix M, the vector of arc
% capacities z, and a sequence of demand vectors d, i.e., d may be a matrix of
% column-wise concatenated demand vectors (or a single column vector).
%
% The function returns the value the max-flow (sum over all flows in case of
% multiple flows), the solution vector x, and the block-wise expanded incidence
% matrix N. The product N*x gives the injection / reception rates for each flow.

function [N x val] = capacity(M,z,d)

% read dimensions
[n m] = size(M);
[n flows] = size(d);

% n = number of nodes
% m = number of arcs
% f = number of flows

N = [];
D = [];
A = [];
for i=1:flows
	N = blkdiag(N,M);
	D = blkdiag(D,d(:,i));
	A = [A eye(m,m)];
end

[na ma] = size(N);

lb = zeros(ma,1);
b = zeros(na,1);

B = N - D*inv(D'*D)*D'*N;

F = inv(D'*D)*D'*N;
f = ones(1,flows)*F;

%     x = linprog(f, A,b,Aeq,beq,lb,ub)
[x val] = linprog(-f,A,z,B,b,lb,[]);
val = -val
x
