function [d, NIT, flg] = QR(a,b);
%
% QR Algorithm on the symmetric 
% tridiagonal problem
%
% On input, a is diagonal, b is off-diagonal.
% On output, d is eigenvalues, NIT is # of iterations for each eigenvalue
% flg = 0 is success and ~=0 otherwise.
%
% Ming Gu for Math 128B, 2009
%
a = a(:);
b = b(:);
n = length(a);
d = NaN(n,1);
NIT = zeros(n,1);
tol = max(max(abs(a)),max(abs(b))) * eps / sqrt(n);
for k = n:-1:3
    nit = max(k * 10, 50);
    flg = 1;
    for iter = 1:nit
       t1 = (a(k)+a(k-1))/2;
       t2 = sqrt(((a(k)-a(k-1))/2)^2 + b(k)^2);
       sh = [t1-t2;t1+t2];
       [junk,ind]= min(abs(sh-a(k)));
       shift = sh(ind);
       [G, Y]= planerot([a(1)-shift;b(2)]);
       Mat2  = [G*[a(1) b(2); b(2) a(2)];[0 b(3)]]* G';
       for j = 2:k-2
           a(j-1) = Mat2(1,1);
           [G, Y] = planerot(Mat2(2:3,1));
           b(j)   = Y(1);
           Mat2   = [G*[Mat2(2:3,2),[Mat2(3,2);a(j+1)]]; [0 b(j+2)]] * G';
       end
       a(k-2) = Mat2(1,1);
       [G, Y] = planerot(Mat2(2:3,1));
       Mat2(1:2,:) = G * [Mat2(2:3,2),[Mat2(3,2); a(k)]] * G';
       a(k-1:k) = [Mat2(1,1);Mat2(2,2)];
       b(k-1:k) = [Y(1);Mat2(2,1)];
       if (abs(b(k))<=tol)
           flg = 0;
           NIT(k) = iter;
           break;
       end
    end
    if (flg == 1)
        return;
    end
end
d   = a;
flg = 0;
if (n==1)
    return;
end
d(1:2) = eig([a(1), b(2); b(2) a(2)]);
b(1:2)   = 0;
disp([abs(b), NIT]);
return