function varargout = sltensor_svd(T, n)
%SLTENSOR_SVD Performs a Higher-Order SVD on a Tensor
%
% $ Syntax $
%   - [S, U1, U2, ...] = sltensor_svd(T)
%   - [S, U1, U2, ...] = sltensor_svd(T, n)
%   - [S, Us] = sltensor_svd(...)
%
% $ Description $
%   - [S, U1, U2, ...] = sltensor_svd(T) Performs a Higher Order Singular Value
%   Decomposition to a Tensor T. In the output arguments, S is the core
%   tensor, U1, U2, ... are the singular vector matrices of mode 1, 2,...
%
%   - [S, U1, U2, ...] = sltensor_svd(T, n) Here n specifies the order of the
%   tensor T, n should be not less than ndims(T). If n > ndims(T) is just
%   regarded as a tensor of dimensions d1xd2x...x1. 
%
%   - [S, Us] = sltensor_svd(...) where all mode matrices are returned to Us,
%   which is an 1 x n cell array, with each cell containing a mode matrix.
% 
% $ History $
%   - Created by Dahua Lin on Dec 17th, 2005
%

%% parse and verify the arguments
slchknargs(nargin, 1);
n0 = ndims(T);
if nargin == 1
    n = n0;
else
    if n < n0;
        error('sltoolbox:invalidarg', ...
            'The order %d is too small', n);
    end
end
if nargout > 2 && nargout ~= n+1
    error('sltoolbox:invalidnargout', ...
        'The number of outputs is not valid');
end

%% compute

Us = cell(n, 1);
S = T;
for i = 1 : n0
    M = sltensor_unfold(T, i);
    [d1, d2] = size(M);
    if d1 < d2
        [V, D, U] = svd(M', 0);
    else
        [U, D, V] = svd(M, 0);
    end
    slignorevars(D, V);
    clear D V;
    
    Us{i} = U;
    S = sltensor_multiply(S, U', i);
end
if n > n0
    for i = n0+1 : n
        Us{i} = 1;
    end
end

%% output
varargout{1} = S;
if nargout == 2
    varargout{2} = Us;
else
    varargout(2:n+1) = Us;
end


    






    
    