function varargout = jacobisvd(A, tol, maxsweeps)
%JACOBISVD One-sided Jacobi method for computing
%the reduced SVD of an m x n matrix A, m >= n.
%
%   [U S V] = jacobisvd(A,tol,maxsweeps) returns an m x n orthogonal
%   matrix U, an n x n matrix diagonal matrix S containing
%   the singular values of A, and an n x n orthogonal
%   matrix V. numsweeps is the number of sweeps required.
%   If the error tolerance is not attained in maxsweeps,
%   numsweeps is set to -1
%
%   See Algorithm 4.1, p. 32 of the the 1989 article
%   "Jacobi's method is more accurate than QR"
%   by James Demmel and Kresimir Veselic. It can be found at
%   http://www.netlib.org/lapack/lawnspdf/lawn15.pdf

n=size(A,2);
if nargin == 1
   tol = 1.0e-10;
   maxsweeps = 10;
elseif nargin == 2
   maxsweeps = 10;
end

U=A;
V=eye(n);
singvals = zeros(n,1);
errormeasure = tol + 1;
numsweeps = 0;

while errormeasure >= tol && numsweeps <= maxsweeps
	numsweeps = numsweeps + 1;
	for i=1:n-1
		errormeasure = 0;
		for j = i+1:n
			% prepare to calculate the Jacobi rotation
			normcoli = norm(U(:,i));
			normcolj = norm(U(:,j));
			% swap columns if norm of column i of U is less than
			% the norm of column j. Suggested by Hestenes (1958)
			% must also swap the same columns in V.
			% the singular values will appear in decreasing order
			% in S
			if normcoli < normcolj
				tmp = U(:,i);
				U(:,i) = U(:,j);
				U(:,j) = tmp;
				tmp = V(:,i);
				V(:,i) = V(:,j);
				V(:,j) = tmp;
			end
			alpha=sum(U(:,i).^2);
			beta=sum(U(:,j).^2);
			gamma=sum(U(:,i).*U(:,j));
		if alpha*beta ~= 0
			errormeasure = max(errormeasure, abs(gamma)/sqrt(alpha*beta));
		end
			% compute Jacobi rotation that makes columns i and j of U
			% orthogonal and also zeros out A'A(i,j) and A'A(j,i).
			if gamma ~= 0
				zeta=(beta-alpha)/(2*gamma);
				if zeta >= 0
					t = 1/(abs(zeta)+sqrt(1+zeta^2));
				else
					t = -1/(abs(zeta)+sqrt(1+zeta^2));
				end
				c=1.0/sqrt(1+t^2);
				s=c*t;
			else
				c = 1;
				s = 0;
			end
			% update columns i and j of U
			t=U(:,i);
			U(:,i)=c*t-s*U(:,j);
			U(:,j)=s*t+c*U(:,j);
			% update matrix V of right singular vectors
			t=V(:,i);
			V(:,i)=c*t - s*V(:,j);
			V(:,j)=s*t + c*V(:,j);
		end
	end
end

% the singular values are the norms of the columns of U
% the left singular vectors are the normalized columns of U
for j=1:n
	singvals(j)=norm(U(:,j));
	if singvals(j) > eps
		U(:,j) = U(:,j)/singvals(j);
	end
end

S = diag(singvals);
if errormeasure >= tol
	numsweeps = -1;
end

if nargout == 4
   varargout{1} = U;
   varargout{2} = S;
   varargout{3} = V;
   varargout{4} = maxsweeps; 
elseif nargout == 3 || nargout == 2
   varargout{1} = U;
   varargout{2} = S;
   varargout{3} = V;
else
   varargout{1} = singvals;
end

end
