参考

[1]《矩阵计算六讲》，P141~143.
[2]《数值分析》. Timothy Sauer著，裴玉茹，马庚宇译. P486.
[3]对称矩阵的特征值（EVD）计算的雅各宾（Jacobi）方法.

代码

代码来源于[1]的伪代码，一些matlab不自带的函数参考[3]，测试用例是从[2]选取的，可能不够多。
测试用例
$A =USV^T\\ \begin{bmatrix} 0 & -\frac{1}{2}\\ 3 & 0\\ 0 & 0 \end{bmatrix} = \begin{bmatrix} 0 & -1 & 0\\ 1 & 0 & 0\\ 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} 3 & 0\\ 0 & \frac{1}{2}\\ 0 & 0 \end{bmatrix} \begin{bmatrix} 1 & 0\\ 0 & 1\\ \end{bmatrix}$

A = [0,-0.5;3,0;0,0];
tic;
[U,S,V] = myJacobiSVD(A,1e-14)
toc;

function [U,sigma,V] = myJacobiSVD(A,epsilon)
% [U,S,V] = myJacobiSVD(A)
% 计算A=USV
% 其中V表示理论的V'
% A可能不是方阵

%% 初始化
[m_1,n_1] = size(A);
n = max(m_1,n_1);
V = eye(n_1);
U = zeros(m_1);
rots=1;
sigma = zeros(n_1,1);
for i=1:n_1
    sigma(i) = vectorMul(A(:,i));
end

%% 迭代
while rots>=1
    rots=0;
    for p=1:n_1-1
        %%%%%%%%%%%%%%%%%%%%%%%%
        % 这个地方放加速代码
        %%%%%%%%%%%%%%%%%%%%%%%%
        
        for q=p+1:n_1
            beta = vectorMul(A(:,p),A(:,q));
            
            sig = sigma(p)*sigma(q);
            if sig>0 & abs(beta)>epsilon*sqrt(sig)
                rots=rots+1;
                [c,s,t] = jacobi(sigma(p),beta,sigma(q));
                sigma(p)=sigma(p)-beta*t;   sigma(q)=sigma(q)+beta*t;
                % 下面这两步实际上是givensMat的计算
                A = givensMatMulRight(c,s,A,p,q);
                V = givensMatMulRight(c,s,V,p,q);
            end
        end
    end
end

%% 重新排列
sort(sigma,'descend');%降序排列

%% 计算矩阵的秩
% 这儿最好还是可以判断一下矩阵的秩是否已知
tol = max(m_1,n_1)*eps(sigma(1));
r = sum(sigma > tol);

%% 计算Sigma和U
for k=1:r
    %sigma(k) = sqrt(sigma(k));
    sigma(k) = sqrt(sigma(k));
    for j=1:n_1
        U(j,k) = A(j,k)/sigma(k);
    end
end


end

下面这个其实直接调用matlab的就行

function [s] = vectorMul(v1,v2)
% function [s] = vectorMul(v)
% 计算s = v1'*v2


if nargin<2
    [m,~]=size(v1);
    s = 0;
    for i=1:m
        s=s+v1(i)^2;
    end
else
    [m,~]=size(v1);
    s=0;
    for i=1:m
        s=s+v1(i)*v2(i);
    end
end


end

其他用到但是没有直接给出的函数在[3]中给出了。
结果是准确的

U =

     0    -1     0
     1     0     0
     0     0     0


S =

    3.0000
    0.5000


V =

     1     0
     0     1

问题补充

在测试上一节的矩阵
$A = \begin{bmatrix} 10^{40} & 10^{29} & 10^{19}\\ 10^{29} & 10^{20} &10^{9}\\ 10^{19} & 10^{9} & 1 \end{bmatrix}$
的时候，发现获得的秩是不对的，本来应该是3，但是实际计算得到的却是1。
MATLAB计算的 $A^TA$ 为

1.00000000000000e+80    1.00000000000000e+69    1.00000000000000e+59
1.00000000000000e+69    1.00000000000000e+58    1.00000000000000e+48
1.00000000000000e+59    1.00000000000000e+48    1.00000000000000e+38

问题还是次数差别太大，超出了机器精度造成的，但是对这种情况确实不行，要考虑别的方法。

一是方阵先求特征值再开方，
二是计算 $B = \begin{bmatrix} 0& A^T\\ A & 0 \end{bmatrix}$ 的特征值。但是后者我还不甚明了。