function [d,rot,Z,V]=Jacobi(A)
% JACOBI computes eigenvalues of a symmetric matrix
%   [d,rot,Z,V]=Jacobi(A) computes the eigenvalues d and eigenvectors
%   V of the symmetric matrix A.  d=Jacobi(A) computes only the
%   eigenvalues without accumulating the vectors.  Only the upper
%   triangular triu(A,1) of the matrix A is used.  rot counts the
%   rotations and Z stores the sum of the absolute values of the off
%   diagonal elements after each sweep.  Jacobi is a translation of
%   the ALGOL 60 Procedure of H. Rutishauser from the Handbook
%   Wilkinson-Reinsch: Linear Algebra.

Z=[];
n=max(size(A));
eivec=1; if nargout<4, eivec=0; end
if eivec, V=eye(n); end
d=diag(A); b=d; z=zeros(size(d));
rot=0; i =0;
sm=sum(sum(abs(triu(A,1))));
while (sm ~=0  )& (i<50)
  i=i+1;                                   % sweep 
  if  i<4, tresh=0.2*sm/n^2;
  else tresh=0;
  end
  for p=1:n-1
    for q=p+1:n
      g=100*abs(A(p,q));
      if (i>4)&(abs(d(p))+g==abs(d(p)))&(abs(d(q))+g==abs(d(q))) 
        A(p,q)=0;
      else
        if abs(A(p,q))>tresh
          h=d(q)-d(p);                     % rotate 
          if abs(h)+g==abs(h), t =A(p,q)/h;
          else		   
            theta=0.5*h/A(p,q);			
            t=1/(abs(theta)+sqrt(1+theta^2));
            if theta < 0, t =-t; end 	
          end                              % end of tan calculation	
          c=1/sqrt(1+ t^2); s=t*c; tau=s/(1+c); 			
          h=t*A(p,q);						
          z(p)=z(p)-h; z(q)=z(q)+h;			
          d(p)=d(p)-h; d(q)=d(q)+h;		
          A(p,q) =0; 				
          g=A(1:p-1,p); h=A(1:p-1,q);	
          A(1:p-1,p)=g -s*(h+g*tau);	
          A(1:p-1,q)=h +s*(g-h*tau);       % end of case 1<=j<p 
          g=A(p,p+1:q-1); h=A(p+1:q-1,q);	
          A(p,p+1:q-1)=g -s*(h'+g*tau);		
          A(p+1:q-1,q)=h +s*(g'-h*tau);    % end of case p<j<q
          g=A(p,q+1:n); h=A(q,q+1:n);						
          A(p,q+1:n)=g -s*(h+g*tau);						
          A(q,q+1:n)=h +s*(g-h*tau);       % of cases q<j<=n 
          if eivec 							
            g=V(:,p); h=V(:,q);		
            V(:,p)=g -s*(h+g*tau);					   
            V(:,q)=h +s*(g-h*tau);			
          end                                                             
          rot=rot + 1;
        end 
      end 
    end 
  end 
  b=b+z; d=b; z=zeros(size(d));
  sm=sum(sum(abs(triu(A,1))));
  Z=[Z,sm];
end  
