function ShowLanczos()
% Examines "Math Book" Lanczos method for large sparse symmetric eigenvalue
% problems.

% Generate a random sparse example...
n = 1000;
density = .01;
A = sprandsym(n,density);
% Compute the three smallest and the three largest eigenvalues of A
% and the corresponding eigenvectors...
[X,D] = eigs(A,6,'BE');
lambda = diag(D);
% Check these results...
clc
disp(' lambda ||Ax - lambda*x || ')
disp('------------------------------------------------')
for j=1:6
res = A*X(:,j) - lambda(j)*X(:,j);
fprintf('%20.15f %10.3e\n',D(j,j),norm(res))
end
% Now run "math book" Lanczos and compare...
% Random initial vector...
q1 = randn(n,1); q1 = q1/norm(q1);
% Perform kmax steps of "math book" Lanczos...
kMax = 200;
[QkMax,TkMax,rkMax] = Lanczos(A,q1,kMax);
% Check out the approximate eigenvalues...
for k=20:20:kMax
% Assess the approximation to the three smallest and three
% largest eigenvalues of A by looking at the three smallest/largest
% eigenvalues of Tk...
Tk = TkMax(1:k,1:k);
Dk = eig(Tk);
[z,idx] = sort(diag(Dk));
theta = [Dk(idx(1:3));Dk(idx(k-2:k))]';
Qk = QkMax(:,1:k);
fprintf('\n\n+++++++++++++++++++++++++++++++++++++++++++++++++++\n')
fprintf('\nk = %3d\n\n',k)
fprintf('|| Qk''Qk - eye(k,k) || = %10.3e\n',norm(Qk'*Qk-eye(k,k)));
fprintf('\nThe three smallest...\n\n')
fprintf('A: ')
fprintf('%20.15f',lambda(1:3))
fprintf('\nTk: ')
fprintf('%20.15f',theta(1:3))
fprintf('\n\nThe three largest...\n\n')
fprintf('A: ')
fprintf('%20.15f',lambda(4:6))
fprintf('\nTk: ')
fprintf('%20.15f',theta(4:6))
fprintf('\n')
end

function [Qk,Tk,rk] = Lanczos(A,q1,k)
% A*Qk = Qk*Tk + rk*ek'
alfa = zeros(k,1);
beta = zeros(k-1,1);
Qk = zeros(length(q1),k);
Qk(:,1) = q1;
for j=1:k
alfa(j) = Qk(:,j)'*A*Qk(:,j);
if j==1
rk = A*Qk(:,j) - alfa(j)*Qk(:,j);
else
rk = A*Qk(:,j) - alfa(j)*Qk(:,j) - beta(j-1)*Qk(:,j-1);
end
if j<k
beta(j) = norm(rk);
Qk(:,j+1) = rk/beta(j);
end
end
% Form Tk...
Tk = diag(beta(1:k-1),-1)+ diag(alfa)+ diag(beta(1:k-1),1);