function S = singvals(A,tol)
%singvals Computes the singular values of a m x n matrix.
%
%   S = singvals(A), where S is an n x 1 vector computes the singular
%   values of A in order from greatest to least. First, reduce A using
%   Householder reflections to a bidiagonal matrix. Then apply the
%   Demmel & Kahan zero-shift QR downward sweep. The bidiagonal matrix has
%   small elements above the diagonal created by successive sweeps of the
%   algorithm.

[m, n] = size(A);
% if A is symmetric, just compute the absolute value of
% its eigenvalues
if (n == m) && (norm(A-A','inf') < eps)
   S = abs(eigsymb(A));
   S = sort(S, 'descend');
   return;
end

if m < n
   A = A';
   tmp = m;
   m = n;
   n = tmp;
end
A = bidiag(A);
k = min(m,n);

if nargin == 1
   tol = 1.0e-8;
end

while k >= 2
	% convergence test.
	if (abs(A(k-1,k)) <= tol*(abs(A(k-1,k-1)) + abs(A(k,k))))% ||...
		A(k-1,k) = 0;
		% tolerance met. reduce the size of the submatrix.
		k = k-1;
	else
	  for i=1:k-1
			% Compute the Givens parameters for a rotation
			% that will zero-out A(i,i+1) and A(i-1,i+1),
         % but makes A(i+1,i) non-zero.
			[c, s] = givensparms(A(i,i), A(i,i+1));
			% Apply the rotation by forming a post-product
			A(1:k,1:k) = givensmulpsvd(A(1:k,1:k),i,i+1,c,s);
			% Compute the Givens parameters for a rotation
			% that will zero-out A(i+1,i) to correct the result
			% of the previous rotation. The rotation makes
			% A(i,i+1) and A(i,i+2) non-zero.
			[c, s] = givensparms(A(i,i),A(i+1,i));
			% Apply the rotation as a pre-product
			A(1:k,1:k) = givensmul(A(1:k,1:k),i,i+1,c,s);
     end
   end
end
S = abs(diag(A));
