function P2()
% Tests SkewSchur3

I = eye(3,3);

clc
disp('Example 1')
A = [ 0 1 2; -1 0 3; -2 -3 0]
[v, alfa] = SkewSchur3(A);
Q = I - 2*v*v';
alfa = alfa
v = v'
Atilde = Q'*A*Q


disp('Example 2')
A = [ 0 1 0; -1 0 3; 0 -3 0]
[v, alfa] = SkewSchur3(A);
Q = I - 2*v*v';
alfa = alfa
v = v'
Atilde = Q'*A*Q

disp('Example 3')
A = [ 0 1 2; -1 0 0; -2 0 0]
[v, alfa] = SkewSchur3(A);
Q = I - 2*v*v';
alfa = alfa
v = v'
Atilde = Q'*A*Q

disp('Example 4')
A = [ 0 0 2; 0 0 3; -2 -3 0]
[v, alfa] = SkewSchur3(A);
Q = I - 2*v*v';
alfa = alfa
v = v'
Atilde = Q'*A*Q

disp('Example 5')
A = [ 0 0 0; 0 0 0; 0 0 0]
[v, alfa] = SkewSchur3(A);
Q = I - 2*v*v';
alfa = alfa
v = v'
Atilde = Q'*A*Q


function [v,alfa] = SkewSchur3CVL(A)
% A is a 3-by-3 skew-symmetric matrix
% v is a unit vector with the property that if Q = I - 2vv', and
% Q'*A*Q = D, then Q'*A*Q = [0 0 0; 0 0 alfa; 0 -alfa 0]

% Compute null vector for A...

z = [A(2,3);-A(1,3); A(1,2)];
if norm(z) == 0
    alfa = 0; v = [1;0;0];
else
    v = House(z);
    Atilde = A - 2*(A*v)*v';   % Atilde = A*Q
    Atilde = Atilde - (2*v)*(v'*A);  % Atilde = Q'*(A*Q)
    alfa = Atilde(3,2)
end

function v = House(x)
n = length(x);
s = norm(x,2);
if s==0
    v = [1;zeros(n-1,1)];
else
    v = x;
    if x(1) >=0
        v(1) = v(1)+s;
    else
        v(1) = v(1)-s;
    end
    v = v/norm(v);
end
    