I think this does the trick:
function B = MatsDiag(A)
[ l c m ] = size(A)
B = sparse( ...
reshape(repmat(reshape(1:m*l,l,[]),c,1),[],1), ...
kron(1:m*c,ones(1,l))', ...
reshape(A(:,:),[],1), ...
m*l,m*c);
end
where I assume that A(:,:,i) is the i-th matrix with size l times c. I shortly tested it with a small number of parameter choices and it seems to work, but I guess there quite some room for optimisation left.
Regards, Laurent