function P3Grade()
% Checks out SkewJacobi

% -1 if block off incorrect a little bit
% -3 if updates of the form Q = Q*J are O(n^3). Should only be O(n)

clc
disp('Small example test for correctness.')
n = 17;
A = randn(n,n);
A = A-A';
tol = 10^(-12);
MaxNumSweeps = 20;
[Q,D,offTraceYOU] = SkewJacobi(A,tol,MaxNumSweeps);
errQ= norm(Q'*Q-eye(n,n));
errA = norm(Q'*A*Q - D);
fprintf('n = %1d\n',n)
fprintf('|| Q''Q  - I || = %10.3e\n',errQ)
fprintf('|| Q''AQ - D || = %10.3e\n',errA)
[Q,DCVL,offTraceCVL] = SkewJacobiCVL(A,tol,MaxNumSweeps);
[p,q] = size(offTraceYOU);
if q>p
    offTraceYOU = offTraceYOU';
end
offTraceYOU = [offTraceYOU;zeros(21-length(offTraceYOU),1)];
offTraceCVL = [offTraceCVL;zeros(21-length(offTraceCVL),1)];
offTRACES = [offTraceYOU offTraceCVL]
TheEigs = [sort(imag(eig(D))) sort(imag(eig(DCVL)))]

disp('Large example test for efficiency.')
n = 201;
A = randn(n,n);
A = A-A';
tol = 10^(-12);
MaxNumSweeps = 20;
tic
[Q,D,offTraceYOU] = SkewJacobi(A,tol,MaxNumSweeps);
tYOU = toc;
tic
[Q,DCVL,offTraceCVL] = SkewJacobiCVL(A,tol,MaxNumSweeps);
tCVL = toc;
fprintf('\nn = %3d   tYOU/tCVL = %6.3f\n',n,tYOU/tCVL)



    


  function [Q,D,BlockOffTrace] = SkewJacobiCVL(A,tol,MaxNumSweeps)
% A is an nxn skew symmetric matrix and n = 2m+1.
% Q is nxn orthogonal.
% D is nxn and block diagonal with 1x1 and 2x2 diagonal blocks.
% BlockOff(D) <= tol*norm(A,'fro') where D = Q'*A*Q.
% BlockOffTrace is a column vector with the property that BlockOffTrace(k) is
% the value of BlockOff(A) after k-1 sweeps.

sweepCount = 0;
errBound = tol*norm(A,'fro');
BlockOffTrace = BlockOff(A);
[n,n] = size(A);
m = floor(n/2);
Q = eye(n,n);
while (sweepCount<=MaxNumSweeps && BlockOffTrace(end) >= errBound)
    for i=1:m
        for j=i+1:m+1
            if j<=m
               v = [2*i-1 2*i 2*j-1 2*j];
               [Qtilde,T] = schur(A(v,v));
            else
               v = [2*i-1 2*i 2*j-1];
               Qtilde = RealSchur3(A(v,v));
            end
            A(:,v) = A(:,v)*Qtilde;
            Q(:,v) = Q(:,v)*Qtilde;
            A(v,:) = Qtilde'*A(v,:); 
        end
    end
    sweepCount = sweepCount+1;
    BlockOffTrace = [BlockOffTrace;BlockOff(A)];
end
D = zeros(n,n);
for k=1:2:n-1
    D(k:k+1,k:k+1) = A(k:k+1,k:k+1);
end

   function alfa = BlockOff(A)
 [n,n] = size(A);
 m = floor(n/2);
 alfa = 0;
 for i=1:m
     iRan = (2*i-1):(2*i);
     for j=i+1:m+1
         if j<=m
             jRan = (2*j-1):(2*j);
         else
             jRan = (2*j-1):(2*j-1);
         end
         alfa = alfa + norm(A(iRan,jRan),'fro')^2;  
     end
 end
 alfa = sqrt(2*alfa);
    
            
            


   function Q = RealSchur3(A)
 % A is a 3x3 real skew-symmetric matrix.
 % Q is a 3x3 orthogonal matrix so that Q'*A*Q has the form
 %                   Q'*A*Q = [ 0  mu  0 ; -mu 0 0 ; 0 0 0]
 
 x = A(1,2); y = A(1,3); z = A(2,3);
 % Note that
 %              [ 0  x  y ] [  z ]   [ 0 ]
 %              [-x  0  z ] [ -y ] = [ 0 ]
 %              [-y -z  0 ] [  x ]   [ 0 ]
 
 [v,beta] = House([z;-y;x]);
 v = v/norm(v);
 Q = eye(3,3) - 2*v*v';
 
 % Since Qtilde*u = [*;0;0] it follows that the
 % first column of Q is a null vector for A. Thus
 % the first column of Q'*A*Q is all zeros. Since this matrix
 % is skew symmetric we must have the form
 %
 %                  [ 0  0   0 ]
 %      Q'*A*Q   =  [ 0  0  mu ]
 %                  [ 0 -mu  0 ]
 Q = Q(:,[ 2 3 1]);
 % Now we have
 %                  [ 0  mu   0 ]
 %      Q'*A*Q   =  [-mu  0   0 ]
 %                  [ 0   0   0 ]