function [U,S,V] = svd_J(A,tol)
% SVDJ Singular value decomposition using Jacobi algorithm.
%
if nargin == 1, tol = eps; end % Tolerancia por defecto
[M,N] = size(A); K = min(M,N); % K es el número de valores sigulares
On = 0;
for c=A, On=On + sum(abs(c).^2); end
On=On./N; % Suma coef.^2/N
Previous_Off = Inf; V = eye(N);
while true
R = 0; % Contar las rotationes
for r = 1:N - 1
for c = r + 1:N
% Calculate the three elements of the implicit matrix B that are
% needed to calculate a Jacobi rotation. Since B is Hermitian, the
% fourth element (b_cr) is not needed.
b_rr = sum(abs(A(:,r)).^2); % Real value.
b_cc = sum(abs(A(:,c)).^2); % Real value.
b_rc = A(:,r)' * A(:,c); % Same type as A.
% Calculate a Jacobi rotation (four elements of G). The two values
% that we calculate are a real value, C = cos(theta) and S, a value
% of the same type as A, such that |S| = sin(theta).
m = abs(b_rc);
if m ~= 0 % If the off-diagonal element is zero, we don't rotate.
tau = (b_cc - b_rr)/(2*m); % tau is real and will be zero if
% the two on-diagonal elements are
% equal. In this case G will be an
% identity matrix, and there is no
% point in further calculating it.
if tau ~= 0
R = R + 1; % Count the rotation we are about to perform.
t = sign(tau)./(abs(tau) + sqrt(1+tau.^ 2));
C = 1./sqrt(1 + t.^ 2);
S = (b_rc.* t.* C)./ m;
% Initialize the rotation matrix, which is the same size as the
% implicit matrix B.
% We have to create an identity matrix here of the same type as A,
% that is, quaternion if A is a quaternion, double if A is double.
% To do this we use a function handle (q.v.) constructed from the
% class type of A. This was done before the loop, since the type
% of A is invariant.
G = eye(N); G(r,r) = C; G(c,c) = C; G(r,c) = S; G(c,r) =-conj(S);
A = A * G;
V = V * G;
end
end
end
end
if R == 0, error('No rotations performed during sweep.'), end
% Calculate the sum of the off-diagonal elements of the matrix B.
B = A' * A;
Off = sum(sum(abs(triu(B,1))))/(N.^2); % Normalise by the matrix size!
if (Off/On) < tol, break; end % Off-diagonal sum is small enough to stop.
if Previous_Off < Off
warning('QTFM:information', ...
'Terminating sweeps: off diagonal sum increased on last sweep.')
break;
end;
Previous_Off = Off;
end
% Extract and sort the singular values. The vector T may be longer than the
% number of singular values (K) in cases where A is not square.
[T,IX] = sort(sqrt(abs(diag(B))),'descend');
if nargout == 0 || nargout == 1 % .. only the singular values are needed.
U = T(1:K);
end
if nargout == 3 % .. the singular vectors and singular values are needed.
A = A(:, IX); % Map the columns of A and V into the same order as the
V = V(:, IX); % singular values, using the sort indices in IX.
% Construct the left singular vectors. These are in A but we need
% to divide each column by the corresponding singular value. This
% calculation is done by replicating T to make a matrix which can
% then be divided into A element-by-element (vectorized division).
U = A./ repmat(T',M,1);
S = diag(T); % Construct a diagonal matrix of singular values from
% the vector T, because when there are three output
% parameters, S is required to be a matrix.
end
end