单边Jacobi求解SVD和双边Jacobi求解SVD方法的Matlab代码

单边Jacobi求解SVD方法的Matlab代码：

function [U,S,V] = svd_jacobi(A)
% SVD-Jacobi
% U*S*V' = A
% A: m*n matrix
% U: m*m matrix
% S: m*n matrix
% V: n*n matrix
% Jacobi algorithm, EISPACK subroutine TRED2, TQLI
% Reference: J.H.Wilkinson, C.Reinsch, Handbook for Automatic Computation, Vol.II, Linear Algebra, Springer, 1971.
[m,n] = size(A);
if m < n
    A = A';
    [m,n] = size(A);
end
U = eye(m);
S = A';
V = eye(n);
maxiter = 30;
eps = 1e-15;
for iter = 1 : maxiter
    for p = 1 : n - 1
        for q = p + 1 : n
            if abs(S(p,q)) > eps
                theta = (S(q,q) - S(p,p)) / (2 * S(p,q));
                t = sign(theta) / (abs(theta) + sqrt(1 + theta^2));
                c = 1 / sqrt(1 + t^2);
                s = c * t;
                G = eye(n);
                G(p,p) = c;
                G(q,q) = c;
                G(p,q) = -s;
                G(q,p) = s;
                S = G' * S * G;
                V = V * G;
            end
        end
    end
    for p = 1 : m - 1
        for q = p + 1 : m
            if abs(S(p,q)) > eps
                theta = (S(q,q) - S(p,p)) / (2 * S(p,q));
                t = sign(theta) / (abs(theta) + sqrt(1 + theta^2));
                c = 1 / sqrt(1 + t^2);
                s = c * t;
                G = eye(m);
                G(p,p) = c;
                G(q,q) = c;
                G(p,q) = -s;
                G(q,p) = s;
                S = G' * S * G;
                U = U * G;
            end
        end
    end
    if max(max(triu(abs(S),1)))) < eps
        break;
    end
end
if m > n
    U = U(:,1:n);
    S = S(1:n,:);
else
    V = V(:,1:m);
    S = S(:,1:m);
end
end

双边Jacobi求解SVD方法的Matlab代码：

function [U,S,V] = svd_jacobi2(A)
% SVD-Jacobi
% U*S*V' = A
% A: m*n matrix
% U: m*m matrix
% S: m*n matrix
% V: n*n matrix
% Jacobi algorithm, EISPACK subroutine TRED2, TQLI
% Reference: J.H.Wilkinson, C.Reinsch, Handbook for Automatic Computation, Vol.II, Linear Algebra, Springer, 1971.
[m,n] = size(A);
if m < n
    A = A';
    [m,n] = size(A);
end
U = eye(m);
S = A';
V = eye(n);
maxiter = 30;
eps = 1e-15;
for iter = 1 : maxiter
    for p = 1 : n - 1
        for q = p + 1 : n
            if abs(S(p,q)) > eps
                theta = (S(q,q) - S(p,p)) / (2 * S(p,q));
                t = sign(theta) / (abs(theta) + sqrt(1 + theta^2));
                c = 1 / sqrt(1 + t^2);
                s = c * t;
                G = eye(n);
                G(p,p) = c;
                G(q,q) = c;
                G(p,q) = -s;
                G(q,p) = s;
                S = G' * S * G;
                V = V * G;
            end
        end
    end
    for p = 1 : m - 1
        for q = p + 1 : m
            if abs(S(p,q)) > eps
                theta = (S(q,q) - S(p,p)) / (2 * S(p,q));
                t = sign(theta) / (abs(theta) + sqrt(1 + theta^2));
                c = 1 / sqrt(1 + t^2);
                s = c * t;
                G = eye(m);
                G(p,p) = c;
                G(q,q) = c;
                G(p,q) = -s;
                G(q,p) = s;
                S = G' * S * G;
                U = U * G;
            end
        end
    end
    if max(max(triu(abs(S),1)))) < eps
        break;
    end
end
if m > n
    U = U(:,1:n);
    S = S(1:n,:);
else
    V = V(:,1:m);
    S = S(:,1:m);
end
end