  function [U,Sigma,V,sweeps] = ParJacSVD(A,tol)
% Parallel Jacobi algorithm for the svd.
%
% A is nxn. 
% U and V are nxn and orthogonal with the
% property that  U'*AV = Sigma and sqrt(off(Sigma))<=tol*norm(A,'fro').
% Initializations...
n = length(A);
m = n/2;
U = eye(n,n);
V = eye(n,n);
delta = tol*norm(A,'fro');
sqrtoffA = norm(A-diag(diag(A)),'fro'); 
sweeps = 0;
while sqrtoffA>delta 
    for p=1:n-1
       % Generate the pth rotation set
       if p==1
           ivals = 1:2:n;
           jvals = 2:2:n;
       else
           iold = ivals;
           jold = jvals;
           ivals = [iold(1) jold(1) iold(2:m-1)];
           jvals = [jold(2:m) iold(m)];
       end
       for q=1:m
            % Solve the qth subproblem of the current rotation set and
            % update...
            i = ivals(q);
            j = jvals(q);
            [c1,s1,c2,s2] = svd2(A([i j],[i j]));
            A([i j],:) = [c1 s1;-s1 c1]'*A([i j],:);
            U(:,[i j]) = U(:,[i j])*[c1 s1;-s1 c1];
            A(:,[i j]) = A(:,[i j])*[c2 s2;-s2 c2];
            V(:,[i j]) = V(:,[i j])*[c2 s2;-s2 c2]; 
       end
    end
    sqrtoffA = norm(A-diag(diag(A)),'fro'); 
    sweeps = sweeps+1;
end
Sigma = A;

  function [c1,s1,c2,s2] = svd2(B)
% B a 2x2 matrix and [c1 s1;-s1 c1]'B[c2 s2;-s2 c2] is diagonal
% Compute [c s;-s c] so [c s;-s c]'*B is symmetric.
% Equating the (1,2) and (2,1) entries we see that
%     cB(1,2)-sB(2,2) = sB(1,1)+cB(2,1)
% i.e., 
%     c(B(1,2)-B(2,1)) = s(B(1,1)+B(2,2))
alfa = B(1,2)-B(2,1);
beta = B(1,1)+B(2,2);
if alfa==0
    c = 1; s = 0;
elseif abs(alfa)>=abs(beta)
    ctn = beta/alfa;
    s = 1/sqrt(1+ctn^2);
    c = s*ctn;
else
    tan = alfa/beta;
    c = 1/sqrt(1+tan^2);
    s = tan*c;
end
[c2,s2] = symSchur2([c s;-s c]'*B,1,2);
% The left rotation is a product:
%        [c1 s1;-s1 c1] = [c s; -s c]*[c2 s2 -s c2]
c1 = c*c2-s*s2;
s1 = c*s2+s*c2;
          
  function [c,s] = symSchur2(A,p,q)
% function [c,s] = symSChur2(A,p,q)
% 2x2 symmetric SChur decomposition
% A is nxn and symmetric and 1<=p<q<=n
% c and s satisfy c^2+s^2=1 and if Q = [c s;-s c], then
% Q'*A([p q],[p q])Q is diagonal.
% GVL4: Algorithm 8.5.1
if A(p,q)~=0
    tau = (A(q,q)-A(p,p))/(2*A(p,q));
    if tau>=0
        t = 1/(tau+sqrt(1+tau^2));
    else
        t = 1/(tau-sqrt(1+tau^2));
    end
    c = 1/sqrt(1+t^2); s = t*c;
else
    c = 1; s = 0;
end
  
% -1 S not symmetric
% -1 for incorrect updates to V and/or U



