function [alpha,beta,Q,v]=Lanczos(A,b,k)
% LANCZOS Lanczos algorithm
%   [alpha,beta,Q,v]=Lanczos(A,b,k) applies k<=n steps of the Lanczos
%   algorithm to the symmetric matrix A starting with the vector b.
%   Computes Q (theoretically) orthogonal and a symmetric tridiagonal 
%   matrix given by the diagonal in the vector alpha, and the super- and
%   sub-diagonal in the vector beta.

Q(:,1)=b/norm(b); 
if k==1, beta=[]; end
for j=1:k
  v=A*Q(:,j);
  alpha(j)=Q(:,j)'*v;
  v=v-alpha(j)*Q(:,j);
  if j>1
    v=v-beta(j-1)*Q(:,j-1);
  end
  if j<k
    beta(j)=norm(v);
    Q(:,j+1)=v/beta(j);
  end
end
