function [d,V]=Jacobi2(A);
% JACOBI2 computes the eigen-decomposition of a symmetric matrix
%   [d,V]=Jacobi2(A) computes the eigen-decomposition of the 
%   symmetric matrix A as Jacobi1, but now using the symmetry

n=length(A); V=eye(size(A));
S=1;
while S~=0,                                            % sweep
  for p=1:n-1
    for q=p+1:n
      if A(p,q)~=0,                                    % rotate a_pq=0
        theta=(A(q,q)-A(p,p))/2/A(p,q);
        t=1/(abs(theta)+sqrt(1+theta^2));
        if theta<0, t=-t; end
        c=1/sqrt(1+t^2); s=c*t;
        h=A(1:p-1,p)*c-A(1:p-1,q)*s;                   % A=U^T*A
        A(1:p-1,q)=A(1:p-1,p)*s+A(1:p-1,q)*c;
        A(1:p-1,p)=h;
        h=A(p,p+1:q-1)*c-A(p+1:q-1,q)'*s;
        A(p+1:q-1,q)=A(p,p+1:q-1)'*s+A(p+1:q-1,q)*c;
        A(p,p+1:q-1)=h;
        h=A(p,q+1:n)*c-A(q,q+1:n)*s;                   % A=A*U
        A(q,q+1:n)=A(p,q+1:n)*s+A(q,q+1:n)*c;
        A(p,q+1:n)=h;
        h=A(p,q)*t;
        A(p,p)=A(p,p)-h;                               % update diagonal 
        A(q,q)=A(q,q)+h;
        A(p,q)=0;
        h=V(:,p)*c-V(:,q)*s;                           % V=V*U
        V(:,q)=V(:,p)*s+V(:,q)*c;
        V(:,p)=h;
      end
    end 
  end 
  A                                                    % test output 
  S=2*(norm(triu(A),'fro')^2-norm(diag(A))^2)
  Sneu=2*norm(triu(A,1),'fro')^2
end 
d=diag(A);
