function [V,D,offVals] = JacobiClassical(A,tol)
% Classical Jacobi algorithm for symmetric Schur decompoistion.
%
% A is nxn and symmetric. V is nxn and orthogonal with the
% property that V'*AV = D and off(D)<=tol*norm(A,'fro').
%
% offVals is a vector that reports the value of off(A) after
% every update.

n = length(A);
V = eye(n,n);
delta = tol*norm(A,'fro');
offA = off(A);
offVals = offA;
while offA>delta
% Find the largest off diagonal entry...
p=1;q=2; maxAij = abs(A(1,2));
for i=1:n
for j=i+1:n
if abs(A(i,j))>maxAij
p=i;q=j;maxAij=abs(A(p,q));
end
end
end
% Solve the (p,q) subproblem...
[V2,D2] = symSchur2(A([p q],[p q]));
% Update: A <--- Vpq'*A*Vpq where Vpq is the identity except
% Vpq([p q],[p q]) == V2. A highly structured similarity transform.
% First do A<-- Vpq'*A...
A([p q],:) = V2'*A([p q],:);
% then A <-- A*Vpq
A(:,[p q]) = A(:,[p q])*V2;
% then V <-- V*Vpq.
V(:,[p q]) = V(:,[p q])*V2;
% We NEVER formed Vpq!
% Now tabulate the reduction in off(A)...
% Math fact: off(A_new)^2 = off(A_old)^2 - 2A(p,q)^2
offA = off(A);
offVals = [offVals;offA];
% Note that the search for the largest off-diagonal entry is O(n^2)
% while the actual update cost is O(n).
end
D = A;

function z = off(A)
% A is an nxn symmetric matrix and z is the
% Frobenius norm of its off-diagonal part.
z = norm(A-diag(diag(A)),'fro');
% i.e., z = sqrt( sum over all i~=j A(i,j)^2)