function d=ImplicitQR(a,b,test)
% IMPLICITQR computes eigevalues of a symmetric tridiagonal matrix 
%   d=ImplicitQR(a,b) computes the eigenvalues d of the matrix
%   T=diag(a)+diag(b,1)+diag(b,-1) with the QR-algorithm with implicit
%   Wilkinson shift.
%   If test==1 we get testoutput for didactical reasons after each sweep

if nargin<3, test=0; end
n=length(a); I=1; 
while n>1;
  for k=I:n-1                               % Check for small b_i and 
    an=abs(a(k))+abs(a(k+1));               % possible deflation
    if an+b(k)==an, b(k)=0; I=k+1; end
  end
  if I==n;                                  % deflation
    n=n-1; I=1;
  elseif I+1==n  
    g=100*abs(b(I)); h=a(I+1)-a(I);
    if abs(h)+g==abs(h), t=b(I)/h;          % 2x2 block: annihilate b(I) 
    else                                    % explicitly by one rotation 
      theta=0.5*h/b(I);
      t=1/(abs(theta)+sqrt(1+theta^2));
      if theta<0, t=-t; end
    end
    a(I)=a(I)-b(I)*t;
    a(I+1)=a(I+1)+b(I)*t;
    b(I)=0;
    n=n-2; I=1;                             % deflation
  else  
    delta=(a(n)-a(n-1))/2;                  % QR-step from I to n 
    if delta==0,                            % with  Wilkinson shift
      sigma=a(n)-abs(b(n-1));
    else
      sigma=a(n)-b(n-1)^2/(delta+sign(delta)*sqrt(delta^2+b(n-1)^2));
    end
    [c s]=GivensRotation(a(I)-sigma, b(I));   % first transformation
    d=a(I)-a(I+1); z=s*(s*d+2*c*b(I));
    a(I)=a(I)-z;
    a(I+1)=a(I+1)+z;
    b(I)=c*s*d+b(I)*(c^2-s^2);
    bulge=-s*b(I+1);
    b(I+1)=c*b(I+1);
    for k= I+2:n                              % chasing the bulge 
      [c s]=GivensRotation(b(k-2), bulge);
      balt= b(k-1);
      d=a(k-1)-a(k); z=s*(s*d+2*c*b(k-1));
      a(k-1)=a(k-1)-z;
      a(k)=a(k)+z;
      b(k-2)=c*b(k-2)-s*bulge;
      b(k-1)=c*s*d+b(k-1)*(c^2-s^2);
      if k<n
        bulge=-s*b(k);
        b(k)=c*b(k);
      end
    end 
    if test                                   % test-output
      clc; ['working from I=' int2str(I) ' to n=' int2str(n)]
      T=diag(a)+diag(b,-1)+diag(b,1); T(1:n,1:n), pause
    end
  end 
end
d=sort(a);
