 function  [V, T]=symlanz(A, v0, k)
%
%-----------------------------------------------------------
%  k-step (symmetric) Lanczos:    
%
%              A V(:,1:k) = V(:,1:k+1) T 
%-----------------------------------------------------------
%
%  Input:
%
%      A      (n-by-n) Hermitian
%      v0     n-vector, used to start Lanczos process (applied A)
%      k      # of Lanczos steps  requested
%      rtol   relative tolerance for regarding a vector a zero
%             after projection.  
%
%  Output:
%
%      V      n-by-(k+1)  Lanczos vectors
%      T      (k+1)-by-k 
%
%                                RCL 06/28/09
 [m,n]=size(A);
 if m ~= n,
    disp('symlanz: the input matrix is not square');
    return
 elseif max(max(abs(A-A'))) > 0.0
    disp('symlanz: the input matrix is not Hermitian');
    return
 end
%
 nrm2 = norm(v0);
 if nrm2 == 0.0,
    disp('symlanz: input v0 is a zero vector');
    return
 end
 
 V = zeros(n,k+1); T = zeros(k+1,k); 
 V(:,1)=v0/nrm2; 
%
 vh = A*V(:,1); nrm1 = norm(vh,1); 
 j=1;
 T(1,1) = V(:,1)' * vh;
 vh = vh-T(1,1)*V(:,1); 
 nrm2 = norm(vh);
 T(2,1)=nrm2; V(:,2) = vh/nrm2; 
 for j=2:k,
     T(j-1,j)=T(j,j-1);
     vh = A*V(:,j);  
     vh = vh - T(j-1,j)*V(:,j-1);
     T(j,j) = V(:,j)' * vh;
     vh = vh - T(j,j)*V(:,j);
     % if reorthogonalization
     %
     vh =vh-V(:,1:j)*( V(:,1:j)'*vh );
     %
     nrm2 = norm(vh);
     T(j+1,j) = nrm2; V(:,j+1) = vh/nrm2; 
 end
