function lambda = eigenvalue(T)

% function lambda = eigenvalue(T)
% QR iteration 
% eigenvalues lambda of a symmetric tridiagonal matrix 
% lower and main diagonal stored in T(:,1:2)

TOL = 10^(-10); MAX_IT = 100; 

n = size(T,1);
if n == 1; 
   lambda = T(1,2);
elseif n == 2; 
   lambda = shift(T);
else;
   T = [T; 0 0];
   for i = 1:MAX_IT;
      [t_min, k_min] = min(abs(T(2:n,1)));
      if t_min < TOL; 
         break;
      end;
      s = shift([T(n-1,2) T(n,1); T(n,1) T(n,2)]);
      v = [T(1,2) - s(1); T(2,1)]; 
      for j = 1:n-1;
         A = [T(j,:) T(j+1,1);v(2) T(j+1,:);0 0 T(j+2,1)];
         Q = g_par(v); 
         A(1:2,:) = Q*A(1:2,:);
         A(:,2:3) = A(:,2:3)*Q';
         T(j:j+1,:) = [A(1,1:2); A(2,2:3)];
         T(j+2,1) = A(3,3);
         v = [A(2,2); A(3,2)];
      end;
   end;
   lambda = [eigenvalue(T(1:k_min,:)); ...
             eigenvalue(T(k_min+1:n,:))];
end;

lambda = sort(lambda);

%%%%% local functions 

function s = shift(T)

% eigenvalues of [T(1,2) T(2,1); T(2,1) T(2,2)]

m = sum(T(:,2))/2;
r = sqrt(m^2 - T(1,2)*T(2,2) + T(2,1)^2);
s = [m-r; m+r]; 
if T(2,2) > T(1,2); 
   s = s([2 1]);
end;

%%%%%

function Q = g_par(v) 

% function Q = g_par(v) 
% Givens matrix based on the vector v 

if norm(v) == 0; 
   v = [1; 0];
end;

r = norm(v);
sigma = sign(v(1)); 
if sigma == 0; sigma = 1; end; 
c =  sigma*v(1)/r;
s = -sigma*v(2)/r;

Q = [c -s; s c];