function [U,S,V]= jacobi_svd(A)  

TOL=1.e-8;
n=size(A,1);
U=A;
V=eye(n);
converge=TOL+1;
while converge>TOL
  converge=0;
  for j=2:n
    for i=1:j-1

      alpha=U(:,i)'*U(:,i);  
      beta=U(:,j)'*U(:,j);   
      gamma=U(:,i)'*U(:,j);  
      converge=max(converge,abs(gamma)/sqrt(alpha*beta));
      zeta=(beta-alpha)/(2*gamma);
      t=sign(zeta)/(abs(zeta)+sqrt(1+zeta^2));
      c=1/((1+t^2)^(.5));
      s=c*t;
      
      t=U(:,i);
      U(:,i)=c*t-s*U(:,j);
      U(:,j)=s*t+c*U(:,j);
      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));
  U(:,j)=U(:,j)/singvals(j);
end
S=diag(singvals);