function [rU1, rU2, rU3, rLambda, resultF] = hopm( A, nTerms )
% HOPM( F, NTERMS )   Higher Order Power Method for NTERMS
%
% From the paper "On the Best Rank-1 and Rank-(R1,R2...,Rn)
% Approximation of Higer-Order Tensors".
% This function implements the Rank-1 approximation.
%
% Author: Greg Coombe
% Date: August 4, 2003
%
N = 3; % The matrices are (I1 x I2 x I3), for now

% Initialize the result
resultF = zeros(size(A));

% Power iteration for nTerms
for term=1:nTerms
    term
    tic;
	F = A - resultF;
    
	% Generate the initial values
	%tF = flatten( F, 1 );
	%[U1, s, v] = svd( tF );
	%U1 = U1(:,1);
	U1 = rand(size(F,1),1);
	
	%tF = flatten( F, 2 );
	%[U2, s, v] = svd( tF );
	%U2 = U2(:,1);
	U2 = rand(size(F,2),1);
	
	%tF = flatten( F, 3 );
	%[U3, s, v] = svd( tF ); % takes too long!!!
	U3 = rand(size(F,3),1);
	
	% Setup the convergence criteria
	thd = 1e-7;
	old_lambda = 1e5 * [1 1 1];
	converged = false;
	
	while ( ~converged )
        tU1 = squeeze( tmul( tmul(F, U2', 2), U3', 3) );
        lambda(1) = norm(tU1);
        U1 = tU1 / lambda(1);
        
        tU2 = squeeze( tmul( tmul(F, U1', 1), U3', 3));
        lambda(2) = norm(tU2);
        U2 = tU2' / lambda(2);
	
        tU3 = squeeze( tmul( tmul(F, U1', 1), U2', 2));
        lambda(3) = norm(tU3);
        U3 = tU3 / lambda(3);
        
        converged = ( abs(norm(lambda - old_lambda)) < thd );
        old_lambda = lambda;
	end	
    
    % Now reconstruct the matrix, and figure out the remaining
    % residual (to be approximated on the next pass )
	for i=1:size(U1,1),
        for j=1:size(U2,1),
            for k=1:size(U3,1),
                val = lambda(1) * U1(i)*U2(j)*U3(k);
                resultF(i,j,k) = resultF(i,j,k) + val;
            end
        end
	end
    
    % Return the vectors
    rU1(:,term) = U1;
    rU2(:,term) = U2;
    rU3(:,term) = U3;
    rLambda(:,term) = lambda';
    
    iteration_time = toc
end % num_terms