% Aufgabe 27
clear all

A = gallery('poisson',30);
n = size(A,1);

disp(['Size of A: ', num2str(n)])

q1 = rand(n,1); %(1:n)';
q1 = q1/norm(q1);

ak = A*q1;                     
alpha(1) = q1'*ak;                 

r = ak - alpha(1)*q1;              
beta(1)=norm(r); 

% Compute Eigenvalues of T_1
e_max(1) = eigSturm(alpha,beta,1,1e-14);
e_min(1) = e_max(1);

k = 2;
tic
while k<=100 %abs(beta(end))>1e-4 && k <= n 
  % Lanczos iteration
  q0 = q1;
  q1 = r/beta(end);  
  ak = A * q1;               
  alpha(k) = q1'*ak;    

  % Compute Eigenvalues of T_k
  e_max(k) = eigSturm(alpha,beta,1, 1e-15);
  e_min(k) = eigSturm(alpha,beta,k, 1e-15);
 
  r = ak - alpha(k) *q1 - beta(k-1)*q0;                                 
  beta(k) = norm(r);  
  
  if abs(e_max(k)-e_max(k-1)) < 1e-15 * abs(e_max(k)) && ...
     abs(e_min(k)-e_min(k-1)) < 1e-15 * abs(e_min(k))
     break
  end
  k= k+1;  
end
disp(['Time for Lanczos: ', num2str(toc)])

% figure(1)
% tic
% ew = eig(A);
% disp(['Time for eig:    ', num2str(toc)])

% plot([1,k],[max(ew), max(ew)],'k--')
% hold on 
% plot([1,k],[min(ew), min(ew)],'r--')
plot(1:k,e_max,'-o',1:k,e_min,'-o')
% hold off
% 
% figure(2)
% semilogy(1:k,abs(max(ew)-e_max),'--b', 1:k,abs(min(ew)-e_min),'-r')