function sigma = BiSVD(a,b)
%%%
%%%  This code computes the smallest singular
%%%  values of an  upper bidiagonal matrix to high 
%%%  relative accuracy, using the dqds method without shift. 
%%%  We will not consider overflow/underflow in this code.
%%%
%%%  Written by Ming Gu for MA221, Fall 2007
%%%
ITER = 200;
n    = length(a);
A    = diag(a)+diag(b(1:end-1),1);
a    = a.^2;
b    = b.^2;
for k=1:ITER
    d = a(1);
%%%    disp([k,b(n-1),sqrt(a(n))]);
    for j=1:n-1
        a(j)=d+b(j);
        t   = a(j+1)/a(j);
        b(j)= b(j)*t;
        d   = d*t;
    end
    a(n) = d;
    sigma(k) = sqrt(d);
    if (k>1)
       if (b(n-1)== 0 | abs((sigma(end)-sigma(end-1))/sigma(end))<5e-16)
           break;
       end
    end
end
%%%disp([length(sigma), sigma(end),abs(sigma(end)-min(svd(A)))/sigma(end)]);