%mode-k matricization M back to m-th order, n-dimensional tensor A
%size of M must be n by n^(m-1) for n>=2 and m>=2
%k is optional unless k>=2, default k=1
function A=mat2t(M,k)
if nargin<2
    k=1;
end
s=size(M);
n=s(1);
m=round(log(s(2))/log(s(1)))+1;
N1=n^(m-1);
s=repmat(n,1,m);
A=zeros(s);
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)]);
        A(indA{:})=M(i,j);
    end
end



