function Mar8Lec()

n = 6;
A = randn(n,n);
A = A+A';
[V,D,offVals] = JacobiCyclic(A,.00001);

function [V,D,offVals] = JacobiCyclic(A,tol)
% function [V,D,offVals] = JacobiCyclic(A,tol)
% Cyclic 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.
% GVL4: Algorithm 8.5.3

n = length(A);
V = eye(n,n);
delta = tol*norm(A,'fro');
offA = norm(A-diag(diag(A)),'fro');
offVals = offA;
while offA>delta
    for p=1:n-1
        for q = p+1:n
           [c,s] = symSchur2(A,p,q);
           A([p q],:) = [c s;-s c]'*A([p q],:);
           A(:,[p q]) = A(:,[p q])*[c s;-s c];
           V(:,[p q]) = V(:,[p q])*[c s;-s c]; 
           offA = norm(A-diag(diag(A)),'fro'); 
           offVals = [offVals;offA];
           clc
           format short
           A(p,q) = 0;
           A(q,p) = 0;
           fprintf('p = %1d  q = %1d\n',p,q)
           A = A
           pause(.31)
        end
    end
end
D = A;

  function [c,s] = symSchur2(A,p,q)
% function [c,s] = symSChur2(A,p,q)
% 2x2 symmetric SChur decomposition
% A is nxn and symmetric and 1<=p<q<=n
% c and s satisfy c^2+s^2=1 and if Q = [c s;-s c], then
% Q'*A([p q],[p q])Q is diagonal.
% GVL4: Algorithm 8.5.1
if A(p,q)~=0
    tau = (A(q,q)-A(p,p))/(2*A(p,q));
    if tau>=0
        t = 1/(tau+sqrt(1+tau^2));
    else
        t = 1/(tau-sqrt(1+tau^2));
    end
    c = 1/sqrt(1+t^2); s = t*c;
else
    c = 1; s = 0;
end
  


