function [alpha,beta,delta,eta,gamma,rho,W,V] = lanczos(A,step,r,l,n)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  File name: lanczos.m
%  Purpose:   Lanczos process 
%  Inputs:    A:     the matrix to be tridiagonalized
%             r:     the right initial vector
%             l:     the left initial vector
%             n:     dimension of A, r and l 
%             step:  the number of Lanczos steps, step <= n 
%  Outputs:   alpha,beta,eta,gamma,rho, parameters of tridiagonal matrices  
%             delta: biorthogonal property 
%             W, V   left and right Lanczos vectors  
%
%  Inputs and outputs satisfy the following governing equations: 
%
%    T = tridiag(kmax,rho(2:kmax),alpha(1:kmax),beta(2:kmax));
%    Tt = tridiag(kmax,eta(2:kmax),alpha(1:kmax),gamma(2:kmax));
%    ej = zeros(kmax,1);
%    ej(kmax) = 1;
%    res_r = norm(A*V(:,1:kmax)-V(:,1:kmax)*T-V(:,kmax+1)*ej'*rho(kmax+1),1)
%    res_l = norm(A.'*W(:,1:kmax)-W(:,1:kmax)*Tt-W(:,kmax+1)*ej'*eta(kmax+1),1)
%    D = diag(delta(1:kmax));
%    biorth_WV = norm(W(:,1:kmax).'*V(:,1:kmax) - D,1 )
%    Tdiff = norm( D*full(T)/D - full(Tt).', 1)
%
%  Note: the implementation has full rebiorthogonalization 
%
%  Written by Z. Bai, Jan 27  1997 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

j=1;

rho(j)=norm(r,2);
eta(j)=norm(l,2);
V(1:n,j)=r/rho(j);
W(1:n,j)=l/eta(j);

delta(j)=W(1:n,j).'*V(1:n,j);

temp = A*V(1:n,j);

alpha(j)=(W(1:n,j).'*temp)/delta(j);

v=temp-alpha(j)*V(1:n,j);
rho(j+1)=norm(v,2);
V(1:n,j+1)=v/rho(j+1);

tempt = A.'*W(1:n,j); 
w=tempt-alpha(j)*W(1:n,j);
eta(j+1)=norm(w,2);
W(1:n,j+1)=w/eta(j+1);

% Lanczos loop

for j=2:step,

   delta(j)=W(1:n,j).'*V(1:n,j);

   temp = A*V(1:n,j);

   alpha(j)=( W(1:n,j).'*temp ) / delta(j);

   beta(j)=delta(j)/delta(j-1)*eta(j);
   gamma(j)=delta(j)/delta(j-1)*rho(j);

   %construct next pair of Lanczos vectors 

   if j<n, 

      v = temp - alpha(j)*V(1:n,j) - V(1:n,j-1)*beta(j);

      % do full re-orthogonalization  (option)
      for k = 1:j,
          v = v - V(:,k)*(W(:,k).'*v)/delta(k);
      end 

      rho(j+1)=norm(v,2);
      V(1:n,j+1)=v/rho(j+1);

      tempt = A.'*W(1:n,j);

      w = tempt-alpha(j)*W(1:n,j)-gamma(j)*W(1:n,j-1);

      % do full re-orthogonalization  (option)
      for k = 1:j,
          w = w - W(:,k)*(V(:,k).'*w)/delta(k);
      end

      eta(j+1)=norm(w,2);
      W(1:n,j+1)=w/eta(j+1);

   end

end