%mode-k matricization of m-th order, n-dimensional tensor A
%k is optional unless k>=2, default k=1
function M=t2mat(A,k)
if nargin<2
    k=1;
end
s=size(A);
n=s(1);
m=length(s);
N1=n^(m-1);
M=zeros(n,N1);
ss=s(2:end);
ind=cell(1,m-1);
for j=1:N1
    [ind{:}]=ind2sub(ss,j);
    indj=[ind{:}];
    for i=1:n
        indA=num2cell([indj(1:k-1) i indj(k:end)]);
        M(i,j)=A(indA{:});
    end
end

