function cake_division=cake_division(T);

[M,N]=size(T);

f=(T(1,:))'; f(N+1:M*N)=0; f=-f;

num_c=2*N*M+2*N+2*(M-1); % number of inequality constraints; 
                         % note that an equation constraint 
                         % is equivalent to two inequality
                         % constraints
                      
A=zeros(num_c,N*M); b=zeros(num_c,1);

for i=1:N*M,
    A(i,i)=1; b(i)=1;
end;
for i=1:N*M,
    A(N*M+i,i)=-1; b(N*M+i)=0;
end;
for i=1:N,
    A(2*N*M+i,i:N:N*M)=1; b(2*N*M+i)=1;
end;
for i=1:N,
    A(2*N*M+N+i,i:N:N*M)=-1; b(2*N*M+N+i)=-1;
end;
for i=2:M,
    A(2*N*M+2*N+i-1,1:N)=T(1,1:N);
    A(2*N*M+2*N+i-1,(i-1)*N+1:i*N)=-T(i,1:N); 
    b(2*N*M+2*N+i-1)=0;
end;
for i=2:M,
    A(2*N*M+2*N+M+i-2,1:N)=-T(1,1:N);
    A(2*N*M+2*N+M+i-2,(i-1)*N+1:i*N)=T(i,1:N); 
    b(2*N*M+2*N+M+i-2)=0;
end;

options=optimset('display','off');
X=linprog(f,A,b,[],[],[],[],[],options);

for i=1:M,
    D(i,:)=X((i-1)*N+1:i*N);
end;

cake_division=D;

