% Runs Lanczos on a Hermitian matrix, 
% selective orthogonalization as in Parlett et al Fortran cod
% Incoming parameters:
% n matrix order
% nev number of eigenvalues sought, default=1
% jmax maximum size of basis, max number of matrix vector multiplies
% jtestint frequency of testing for convergence
% vstart starting vector default randn(n,1)
%
cpustart=cputime;
reps=10*sqrt(eps);
eps1=eps;
nconv=0;
lmbc=0;
r=vstart;
v=[];
alpha=[];
beta=[];
wjps=[];
beta(1)=norm(r);
sqn=1;
bijn=1;
fprintf('   j cev  cputime    betaji     lambda        squotient \n')
for j=1:jmax,
  % Basic recursion
  v(:,j)=r/beta(j);
  r=A*v(:,j);
  if j>1, r=r-v(:,j-1)*beta(j); end;
  alpha(j)=r'*v(:,j);
  r=r-v(:,j)*alpha(j);
  beta(j+1)=norm(r);
  % Reorthogonalize
  h1=v'*r;
  r=r-v*h1;
  nh1=norm(h1);
  wjps(j)=nh1;
  if nh1>reps,
     r=r-v*(v'*r);
  end;
  
    % Compute Ritz value and test for convergence or stagnation
    Tj=diag(alpha)+diag(beta(2:j),1)+diag(beta(2:j),-1);
    [s d]=eig(Tj);
    thetas=diag(d);
    bndv=abs(beta(j+1)*s(j,:));
    % Take Ritz vector closest to start vector
    [s1c it]=max(abs(s(1,:)));
    bndt=bndv(it);
    lambdat=d(it,it);
    bijo=bijn;
    bijn=bndt;
    sqo=sqn;
    sqn=abs(s(j,it)/s(1,it));
    fprintf('%4.0f %3.0f %8.2e %8.2e %15.8e %8.2e \n',j,nconv,cputime-cpustart,bndt,lambdat,sqn)
  if beta(j+1)<epert | bndt<tolconv , break; end;
end;
lmb=d(it,it);
xv=v*s(:,it);



  
