function [T v]=Lanczos(A,disp,MaxIter,vstart)

%  [T V]=Lanczos(A,disp,vstart)
%
% Use Paige's symmetric Lanzcos Algorithm (given as Alg 1.1 p. 8 in Meurant and in Cullum &
% Willoughby p.48 and 103) 
% to compute the factorization the tridiagonal matrix T similar to A:
%     V'AV = T
% and V is orthonormal Lanczos Vectors 


% Initializing
n=size(A,1);

if ~exist('disp','var')
  disp=1;
end
disp=1*disp;

if ~exist('MaxIter','var')|isempty(MaxIter);
    MaxIter=n;
   
end

m=MaxIter;



%m=2*n;   %<-- do this if I expect badly behaved problems
vhat=zeros(n,m);                           
v=zeros(n,m);
u=zeros(n,m);

alpha=zeros(m,1);
eta=zeros(m,1);

% Initialize
if ~exist('vstart','var')
  vhat(:,1)=randn(n,1);       %r=b-A*x;                           % residual, steepest descent
else
  vhat(:,1)=vstart;
end
eta(1)=norm(vhat(:,1));
v(:,1)=vhat(:,1)/eta(1);
alpha(1)=v(:,1)'*A*v(:,1);
vhat(:,2)=(A - diag(ones(n,1)*alpha(1)))*v(:,1);

i=1;
lognvtol=-8;
logdtol=-8;

while 1             % CG stopping tolerance - see 6.3, p 136
    
  i=i+1;
  
  eta(i)=sqrt(vhat(:,i)'*vhat(:,i));
  
  if (eta(i)==0)|(log10(eta(i))< logdtol)   
       if disp
         fprintf('\nforced exit: log10(eta(i)) is less than tolerance = %d',logdtol)
       end
       %keyboard
       break
  end
  
  v(:,i)=vhat(:,i)/eta(i);
  u(:,i-1)=A*v(:,i)-eta(i)*v(:,i-1);
  alpha(i)=v(:,i)'*u(:,i-1);
  vhat(:,i+1)=u(:,i-1) - v(:,i)*alpha(i);
  
  if (alpha(i)==0)|(log10(alpha(i))< logdtol)   
       if disp
         fprintf('\nforced exit: log10(alpha(i)) is less than tolerance = %d',logdtol)
       end
       %keyboard
       break
  end
  
  
  
  if (log10(eta(i))<lognvtol) | (i>=m)   
     if disp; fprintf('\nforced exit: At iter=%d (n=%d)',i,n); end
     
     %T=spdiags([[eta(2:i); 0] alpha [0; eta(2:i)]],-1:1,i,i);    
          
     %keyboard
     break
  end
  
end

%T=GetTriDiag(alpha,eta);
T=spdiags([[eta(2:i); 0] alpha(1:i) [0; eta(2:i)]],-1:1,i,i);     
          

if disp
  if alpha(i)&disp
         fprintf('\nExit Lanczos, log10||v||=%.2e\tlog10(d)=%f at iter=%d (n=%d)\n',log10(eta(i)),log10(alpha(i)),i,n)
  elseif disp
         fprintf('\nExit Lanczos, ||r||=%.2e\td=%f at iter=%d (n=%d)\n',eta(i),alpha(i),i,n)
  end
  
  
end

%u(i,i+1)=spdiags([ones(n,1),[0;alpha]],0:1,n,n);