%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% lanczos()
%   This function iterates Arnoldi factorization by k steps 
%   for a symmetric matrix
%   Input : 
%         A - Symmetric matrix to be factorized
%         v1 - initial vector for V
%         k - # of iterations
%         tol - tolerance for breakdown
%     Output
%         T - k by k trigoanl matrix 
%         V - n by k Arnoldi matrix
% Sungwoo Park
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [T,V,w] = lanczos(A,v1,k,tol)

% Initialization
v0=0;
v1 = v1/norm(v1);
V = v1;
beta = 0;
alpha = 0;
T = spalloc(k,k,3*k);

% Compute norm of A
if issparse(A)
    anorm =  normest(A);
else
    anorm = norm(A);
end

% Iterate k-steps Arnoldi
for j=1:k

    w = A*v1 - beta*v0;
    alpha= w'*v1;
    T(j,j) = alpha;
    w = w-alpha*v1;
    beta= norm(w);
    if j<k,
        T(j,j+1) = beta;
        T(j+1,j) = beta;
    end
    if j~=k && beta<anorm*tol,
        fprintf('[Break Down] Early termination with %d-th iteration\n',j);
        break;
    end
    v0 = v1;
    v1 = w/beta;
    V = [V v1];
end

end