% Unpreconditioned Lanczos algorithm :
%               A*V_k = V_k*T_k
%       [alfa1  beta2			  ]
%       [beta2  alfa2  beta3		  ]
% T_k = [       beta3  alfa3 beta4	  ]
%       [ ...				  ]
%       [ 		   beta_k   alfa_k]
%
function [lanmin,lanmax]=Lanczos(A)

eps_cmp=0.01;
disp(['This routine will estimate the larges and smallest eigenvalue'])
disp(['of a matrix A within a relative tolerance ' num2str(eps_cmp) ',']) 
disp(['i.e., abs(eig(k)-eig(k-1))<=' num2str(eps_cmp)])

lanmax=100000;
lanmin=0;
s=size(A,1);
d(:,1)=rand(s,1); d(:,1)=d(:,1)/norm(d(:,1)); % normalized random
g       =  A*d(:,1);
alfa(1) = (g'*d(:,1));
beta(1) = 0;
g       = g - alfa(1)*d(:,1);
d(:,2)  = g/norm(g);
it=2; T(1,1)=alfa(1);
% ----------------------
while it<s,
     h        = A*d(:,it);
     alfa(it) = h'*d(:,it);
     beta(it) = h'*d(:,it-1);
     g        = h - alfa(it)*d(:,it) - beta(it)*d(:,it-1);
     T(it,it)=alfa(it);T(it,it-1)=beta(it);T(it-1,it)=beta(it);
     if it>1,
        clear e  %,size(T), T,wait
        e=sort(eig(T)); %min(e),max(e),eps_cmp,wait 
% - - - - - - - - - - - - - - check for the largest and the smallest
        if (abs(lanmax-max(e))<eps_cmp)&(abs(lanmin-min(e))<eps_cmp),
           return
        end
        if abs(lanmax-max(e))>eps_cmp, lanmax = max(e);  end
        if abs(lanmin-min(e))>eps_cmp, lanmin = min(e);  end
     end
     d(:,it+1)=g/norm(g);
     it=it+1;
end