function P2Grade()
% Checks out TwoExpos

% A small example to compare to naive approach
n = 5;
% Generate a random skew-symmetric matrix
S = randn(n,n);
S = S-S';
% Stilde is S with the (1,n) and(n,1) entries zeroed...
Stilde = S; Stilde(1,n) = 0; Stilde(n,1) = 0;
clc
[F,Ftilde] = TwoExpos(S)
I = eye(n,n);
F0 = (I-S/2)\(I+S/2);
F0tilde = (I-Stilde/2)\(I+Stilde/2);
fprintf('Error in F = %10.3e\n',norm(F-F0))
fprintf('Error in Ftilde = %10.3e\n',norm(Ftilde-F0tilde))

% A larger example to see if Ftilde is basically "free" given that F
% has been computed via a single nxn LU factorization
n = 1000; nRepeat = 1;
% Generate a random skew-symmetric matrix
S = randn(n,n);
S = S-S';
% How long does your code take?
tic
for k=1:nRepeat
[F,Ftilde] = TwoExpos(S);
end
tTwoExposYOU = toc;
% How long does my code take?
tic
for k=1:nRepeat
[F,Ftilde] = TwoExposCVL(S);
end
tTwoExposCVL = toc;
% How long does it take to get just F?
tic
for k=1:nRepeat
I = eye(n,n);
[L,U,piv] = lu(I-S/2,'vector');
F0 = (I+S/2);
F0 = U\(L\F0(piv,:));
end
tF = toc;
fprintf('\n\nn=%1d tTwoExposYOU/tF = %5.2f\n',n,tTwoExposYOU/tF)
fprintf('n=%1d tTwoExposCVL/tF = %5.2f\n',n,tTwoExposCVL/tF)

% Up to -3 for serious unnecessary O(n^3) overheads

function [F,Ftilde] = TwoExposCVL(S)
% S is an nxn skew-symmetric matrix
% Returns F = R11(S) and Ftilde = R11(Stilde) where Stilde is
% obtained from S by setting the (1,n) and (n,1) entries to zero.

[n,n] = size(S);
% Compute F = R11(S). This is the solution to MF =N where
M = eye(n,n)-S/2;
N = eye(n,n)+S/2;
% Use the LU factorization...
[L,U,piv] = lu(M,'vector');
F = U\(L\N(piv,:));

% We now compute Ftilde, which is the solution to Mtilde*Ftilde = Ntilde
% where
% Mtilde = M + WY'
% Ntilde = N - WY'
% W = (S(1,n)/2)*[e1 en]
% Y = [en -e1]
% eye(n,n) = [e1,...,en]
%
% This follows from the fact that Stilde = S - S(1,n)*(e1*en'-en*e1'.
%
% Sherman-Morrison-Woodbury says that
%
% inv(Mtilde) = inv(M) - inv(M)*W*inv(I+Y'*inv(A)*W)*Y'*inv(M)
%
% and so with the definitions
% Wtilde = inv(M)*W n-by-2
% Qtilde = I + Y'*Wtilde 2-by-2
% Ytilde' = inv(Qtilde)*Y' 2-by-n
%
% we see that
% Ftilde = inv(Mtilde)*(N - WY')
% = (inv(M) - inv(M)*W*inv(I+Y'*inv(A)*W)*Y'*inv(M))*(N - W*Y')
% = (inv(M) - Wtilde*Ytilde'*inv(M))*(N - W*Y')
% = (I - Wtilde*Ytilde')*(F - Wtilde*Y')

% Use the LU of M to compute Wtilde and Ytilde...
W = zeros(n,2); W(1,1) = S(1,n)/2; W(n,2) = S(1,n)/2;
Y = zeros(n,2); Y(n,1)=1; Y(1,2) = -1;
Wtilde = U\(L\W(piv,:));
Qtilde = [1+Wtilde(n,1) Wtilde(n,2) ; -Wtilde(1,1) 1-Wtilde(1,2)];
Ytilde = (Qtilde\Y')';

% Did not exploit stucture in W and Y. That is because we don't care about O(n^2)
% savings given that we have spent O(n^3) factoring M.

% Now compute Ftilde = (I - Wtilde*Ytilde')*(F - Wtilde*Y')...
Ftilde = F - Wtilde*Y';
Ftilde = Ftilde - Wtilde*(Ytilde'*Ftilde);