function [rez]=tensorprod(x,y,varargin)
% [rez]=tensorprod(x,y,z,w, etc)
%
% TENSORPROD
% takes tensor products of multidimensional arrays by summing over
% dimensions of matching size.
%
% the trailing dimension(s) of x must match the leading dimensions of y
% the trailing dimension(s) of y must match the leading dimensions of z
% and so on.
%
% e.g., 
% x=rand(4,5,6,7); y=rand(6,7,8,9,10); z=rand(8,9,10,11,12,13); w=rand(13,14);
% tensorprod(x,y,z,w);
% 
% e.g., 
% size(tensorprod(rand(1,2,5,6),rand(5,6,10,11,12),rand(11,12,20)))
% is [1,2,10,20]
%
% Sam Roweis, June 1999

if(~isempty(varargin)) 
  rez=tensorprod(tensorprod(x,y),varargin{1},varargin{2:end});
  return;
end

sx=size(x); sy=size(y); kk=NaN;
for jj=1:min(length(sx),length(sy)) 
  if(all(sx((end-jj+1):end)==sy(1:jj))) kk=jj; end
end
if(isnan(kk)) rez=[]; return; end

srest=num2cell([sx(1:(end-kk)),sy((kk+1):end)]);
rez=reshape(reshape(x(:),prod(sx(1:(end-kk))),prod(sx((end-kk+1):end)))*...
        reshape(y(:),prod(sy(1:kk)),prod(sy((kk+1):end))),srest{:});
  
