%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% prob1c.m
%   A is gallery('wathen',20,20)
%   This file calls lanczos() to get Arnoldi factorization of A with 100 
%   iterations and compute approximated eigenpairs of A
% Sungwoo Park
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc;
A = gallery('wathen',20,20);

n = size(A,1);
tol = 10^-8;
v0 = spones(ones(n,1));
k = 100;

tic;
[T,V,r] = lanczos(A,v0,k,tol);
fprintf('||r|| = %f\n',norm(r))
L1 = eig(T);
t1 = toc;
tic;
L2 = eig(A);
t2=toc;

fprintf('[time to compute 100 approximated eigenvalues] = %f sec\n',t1);
fprintf('[time to compute full exact eigenvalues] = %f sec\n',t2);

i=1;
j=1;
err = [];
L2 = [-inf; L2; inf];
fprintf('\n[Result]\n');
fprintf('[ i]\t\teig of T |\tclosest eig of A |\trelative error\n');
while(i<=k && j<=n+1)
    if L1(i)>=L2(j) && L1(i)<=L2(j+1)
        if (L1(i)-L2(j))<(L2(j+1)-L1(i))
            r = (L1(i)-L2(j))/L2(j);
            fprintf('[%2d]\t%f\t%f\t\t%f\n',i,L1(i),L2(j),r);
        else
            r = (L2(j+1)-L1(i))/L2(j+1);
            fprintf('[%2d]\t%f\t%f\t\t%f\n',i,L1(i),L2(j+1),r);
        end
        err = [err r];
        i=i+1;
    else
        j=j+1;
    end
end
 
figure(1);
semilogy([1:k],err,'b.');
title('[Relative error between the i-th smallest eigenvalue of T and the closest eigenvalue of A]')
ylabel('Relative error in log scale');
xlabel('i');