
% Demonstration of Jacobi method for finding the eigenvalues and
% eigenvectors of a symmetric matrix.  H. Motteler, Feb 1999
%
% This example constructs and applies the full rotation matices
% explicitly, which makes it easy to see what is happening, but 
% is very inefficient.
% 
% (note: the eigenvalues and vectors may not appear in the same
% order as produced by the matlab "eig" routine)
%

% the symmetric matrix A from example 4.7 in the text
n = 3;
A0 = [ 1  0  2;
       0  2  1;
       2  1  1];

% generate a random symmetric matrix A0
n = 10
A0 = rand(n, n);
A0 = A0'*A0;

% j_eps is an upper bound on off-diagonal elements
j_eps = 1e-4;

A = A0;
V = eye(n);
jflag = 1;
k = 1;

% Loop on "sweeps" of the above-diagonal part of A, until
% all off-diagonal elements are bounded by j_eps

while jflag

  % jflag is set to one if we find an off-diag element > j_eps
  jflag = 0;

  % "sweep" loops
  for i = 1:n
    for j = n:-1:i+1

      % see if this element needs zeroing
      if abs(A(i,j)) > j_eps;
        a = A(i,i);
        b = A(i,j);
        d = A(j,j);

        % calculate sine and cosine of rotation angle
        [s, c] = getrot(a, b, d);

        % generate the rotation matrix 
        Jrot = eye(n);
        Jrot(i,i) =  c;
        Jrot(i,j) =  s;
        Jrot(j,i) = -s;
        Jrot(j,j) =  c;

        % apply the rotation to A
        A = Jrot' * A * Jrot ;

        % eigenvectors are the product of the rotations
        V = V * Jrot;
        jflag = 1;

      end
    end
  end
  k = k + 1; % sweep count
end

% show sweep count
k

% show the results
% diag(A)
% V

% compare to Matlab eigenvector routine
[V2,D2] = eig(A0);
% diag(D2)
% V2

norm(sort(diag(A)) - sort(diag(D2)))

