Symmetric Matrix Eigenvalue Problems - Numerical Methods in Engineering with MATLAB - page 334

P = eye(n); % Initialize rotation matrix
fori=1:maxRot %BeginJacobirotations
[Amax,k,L] = maxElem(A);
if Amax < tol;
eVals = diag(A); eVecs = P;
return
end
[A,P] = rotate(A,P,k,L);
end
error('Too many Jacobi rotations')
function [Amax,k,L] = maxElem(A)
% Finds Amax = A(k,L) (largest off-diag. elem. of A).
n = size(A,1);
Amax = 0;
fori=1:n-1
forj=i+1:n
if abs(A(i,j)) >= Amax
Amax = abs(A(i,j));
k=i;L=j;
end
end
end
function [A,P] = rotate(A,P,k,L)
% zeros A(k,L) by a Jacobi rotation and updates
% transformation matrix P.
n = size(A,1);
diff = A(L,L) - A(k,k);
if abs(A(k,L)) < abs(diff)*1.0e-36
t = A(k,L);
else
phi = diff/(2*A(k,L));
t = 1/(abs(phi) + sqrt(phiˆ2 + 1));
ifphi<0;t=-t;end;
end
c=1/sqrt(tˆ2+1);s=t*c;
tau=s/(1+c);
temp = A(k,L); A(k,L) = 0;
A(k,k) = A(k,k) - t*temp;