I found this code through a friend. However, I am not able to understand how this works. Can someone please explain?

Implementing Jacobi SVD:

function [Vj D] = JacobiEig(A)

Size = size(A,1);
E = eye(Size);
G=cell(1,2);    if true
      % code
    end
Vj = E; % Start with unit matrix
for Rotations=[1:Size^2*20] % Limit number of rotations
  % Find maximum off-diagonal element
  Max = 0;
  for r=1:Size-1
    for c=r+1:Size
      if abs(A(r,c))>Max % New Max found
        p = r; q = c;
        Max = abs(A(r,c));
      end
    end
  end

    % Compare Max with working precision
    if Max<eps
      break % A is diagonalized, stop now   
    end

    % Find sin and cos of rotation angle
    theta = (A(q,q)-A(p,p))/(2*A(p,q));
    t = 1/(abs(theta)+sqrt(theta^2+1));
    if theta<0
      t = -t;
    end
    c = 1/sqrt(t^2+1);
    s = t*c;

    % Build rotation matrix
    P = E;
    P(p,p) = c;
    P(q,q) = c;
    P(p,q) = s;
    P(q,p) = -s;

    % Do rotation
    A = P'*A*P;
    Vj = Vj*P;
  end
  D = diag(spdiags(A,0)); % Return diagonal
  G{1,1}=Vj;
  G{1,2}=D;
  end

Compression of Data:

function [K,m,L,n] = Compress(S,D)

G=cell(1,4);
m=size(D,1);
n=size(S,1);
i=1;
j=1;

for r=1:m
    for c=1:m
        if(D(r,c)~=0)
        K(i)=D(r,c);
        i=i+1;
        end;
    end
end

for r=1:n
    for c=1:n
        if(S(r,c)~=0)
        L(j)=S(r,c);
        j=j+1;
        end;
    end
end

G{1,1}=K;
G{1,2}=m;
G{1,3}=L;
G{1,4}=n;

end

Reconstruction of Data:

function [Mj,Ms] = Recon(K,m,L,n,S,Vj,V,U)

G=cell(1,8);

X=zeros(m,m);
X=diag(L);
Y=zeros(n,n);
Y=diag(K);

for r=1:m
    for c=1:m
        if(abs(X(r,c))<2e-64)
        X(r,c)=0;
        end;
    end
end
for r=1:n
    for c=1:n
        if(abs(Y(r,c))<2e-64)
        Y(r,c)=0;
        end;
    end
end

Ms=U*X*V';
Mj=U*Y*Vj';

G{1,7}=Mj;
G{1,8}=Ms;

End