  function [U,S,V] = jacSVD(A,tol)
% A is an n-by-n matrix. 
% U and V are orthogonal matrices with the property that
% if S = U'*A*V = S, then Off(S) <= tol*norm(A,'fro') where
% Off(S)^2 = sum of the squares of all the off diagonal entries

[n,n] = size(A);
U = eye(n,n);
V = eye(n,n);
c = norm(A,'fro');
while off(A) > tol*c
    for i=1:n-1
        for j=i+1:n
            % Solve the (i,j) problem. 
            % off(A)^2 will be reduced by A(i,j)^2 + A(j,i)^2.
            
            % Compute the SVD of [A(i,i) A(i,j); A(j,i) A(j,j)]
            [Utilde,Stilde,Vtilde] = svd(A([i,j],[i j]));
            
            % Update A, U, and V...
            A(:,[i j]) = A(:,[i j])*Vtilde;
            V(:,[i j]) = V(:,[i j])*Vtilde;
            A([i j],:) = Utilde'*A([i j],:);
            U(:,[i j]) = U(:,[i j])*Utilde;
        end
    end
end
S = A;

function s = off(A)
% A is an n-by-n matrix.
% s = sqrt(sum of the squares of all |A((i,j)| with i not equal to j)
   [n,n] = size(A);
   s = 0;
   for j=1:n
       s = sum(abs(A(1:j-1,j)).^2) + sum(abs(A(j+1:n,j)).^2);
   end
   s = sqrt(s);
