function q=SVDGolubReinsch(A)
% SVDGolubReinsch singular values by the Golub Reinsch algorithm
%   q=SVDGolubReinsch(A) computes the singular values of A using the
%   Golub-Reinsch Algorithm and a machine independent criterium
%   following G. W. Stewart.

[m,n]=size(A); 
if n>m, A=A'; [m,n]=size(A);end
[q,e,A]=Bidiagonalize(A);
k=n;                        
while k~=1,
  splitting=false;
  l=k; 
  while ~splitting & l>1,
    t=abs(q(l-1))+abs(q(l));
    if t+abs(e(l))==t;                % splitting: e(l) is small
      splitting=true; e(l)=0;
    else 
      t=abs(e(l))+abs(e(l-1));
      if t+abs(q(l-1))==t,            % q(l-1) is small, cancellation
        splitting=true; q(l-1)=0;
        bulge=e(l); e(l)=0;           % introduce zero row
        for kk=l:k-1  
          [c s]=GivensRotation1(bulge,q(kk));
          q(kk)=s*bulge+c*q(kk);  
          bulge=-s*e(kk+1);
          e(kk+1)=c*e(kk+1); 
        end
        [c s]=GivensRotation1(bulge, q(k));
        q(k)=s*bulge+c*q(k);
      else
        l=l-1;                        % check elements on row before 
      end 
    end
  end 
  if splitting & l==k,                % convergence 
    if q(k)<0, q(k)=-q(k); end        % adapt sign, singular values >= 0
    k=k-1;
  else
    [e,q]=QRStep(l,k,e,q);        
  end
end 
if q(1)<0, q(1)=-q(1); end
[q]=sort(q); q=q(n:-1:1)';            % sort singular values