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
