function S = singvalsdemo(A)
%SINGVALSDEMO Demonstrates the computation of the singular values
%of an m x n matrix A with m >= n.
%
%   S = singvalsdemo(A), where S is an n x 1 vector computes the singular
%   values of A in order from greatest to least, and outputs what
%   actions take place. First, reduce A to a bidiagonal matrix using
%   Householder reflections. Then apply the Demmel & Kahan zero-shift
%   QR downward sweep. The bidiagonal matrix has increasingly
%   smaller elements above the diagonal created by successive sweeps
%   of the algorithm.

disp('Start with the upper bidiagonal matrix:');
A = bidiag(A)

[m, n] = size(A);
k = min(m,n);

disp('Completed reduction to upper bidiagonal form.');
disp('Ready to begin the chase.');
A(abs(A)<1.e-13)=0;
spy(A);
pause;

while k >= 2
	% Convergence test.
	if abs(A(k-1,k)) <= 1.0e-4*(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 it makes A(i+1,i) non-zero.
			[c, s] = givensparms(A(i,i), A(i,i+1));
         if i > 1
            fprintf('Zero out A(%d,%d) and A(%d,%d), but leave a nonzero at A(%d,%d)\n',i,i+1,i-1,i+1,i+1,i);
         else
            fprintf('Zero out A(%d,%d) and leave a nonzero at A(%d,%d)\n',i,i+1,i+1,i);
         end
			% Apply the rotation by forming a post-product
			A(1:k,1:k) = givensmulpsvd(A(1:k,1:k),i,i+1,c,s);
         A
			A(abs(A)<1.e-5)=0;
			spy(A);
         pause;
         fprintf('\n\n');
			% 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+2) and A(i,i+1) 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);
         
         if i ~= k-1
            fprintf('Zero out A(%d,%d) and leave nonzeros at A(%d,%d) and A(%d,%d).\n',i+1,i,i,i+1,i,i+2);
         else
            fprintf('Zero out A(%d,%d)\n',i+1,i);
         end
         A
 			A(abs(A)<1.e-5)=0;
			spy(A);
			pause;
         fprintf('\n\n');
		end
	end
end

A = diag(diag(A));
spy(A);
pause;
S = sort(abs(diag(A)), 'descend');
