% *****************
% Lanczos iteration
% *****************
% Vesa Apaja

clear all;close all;
N=3000;
e = rand(N,N);
A = spdiags(e,0,N,N);


B=zeros(N,N);
for i=1:N 
    for k=1:2
        j = 1+mod(i+k,N);
        B(i,j) = rand;
        B(j,i) = B(i,j);  
    end
end
A=sparse(A+B);


tic
% built-in result; ordinary diagonalization
all = eig(full(A));
% try built-in sparse matrix diagonalization; find only 6
% eigenvalues
%all = eigs(A)
toc
all = sort(all,'descend');

% sparsity pattern
spy(A);
title ('Sparsity pattern of the matrix')


tic
% Lanczos 
% =======
% start with vector b
b = ones(N,1);
% normalized b to vector q
q = b/norm(b);
% store values of q to Q
Q = q;

n    = 100;
alpha = zeros(n,1);
beta  = zeros(n,1);
for i=1:n
    v = A*q;
    alpha(i) = q'*v;    
    v = v - alpha(i)*q;
    if (i>1) v = v - beta(i-1)*Q(:,i-1); end
    beta(i) = norm(v);
    q  = v/beta(i);
    
    % append q to storage Q
    Q = [Q,q];
end
% Form a tridigonal matrix T
T = diag(alpha(1:n)) + diag(beta(1:n-1),1) ...
    + diag(beta(1:n-1),-1);
% Find eigenvalues D and eigenvectors V of matrix T
figure (2)
hold on;
xlabel('Iteration #')
ylabel('Ritz eigenvalues')
for i=1:n
    [V,D] = eig(T(1:i,1:i));
    ritz = sort(diag(D),'descend');
    % plot(i,ritz,'ok','MarkerSize',1);
end
toc
hold off;
figure(3)
plot(ritz,'o') 
hold on;
plot(all,'*')
legend('Lanczos','full')
hold off;

% Test unitarity of Q
err = zeros(n,1);
for i=1:n
    err(i) = norm(Q(:,1:i)'*Q(:,1:i)-eye(i,i));
end
figure(4)
semilogy(1:n,err,'*');
xlabel('Iteration')
title('Loss of orthogonality')

for i = 1:size(all)
    [min_diff,ind] = min(abs(ritz-all(i)));
    fprintf(['i = %5d     ritz = %20.10f full = %20.10f diff= %20g\n'],i,ritz(ind),all(i),ritz(ind)-all(i))
end

