function [Q,T] = SkewReduce(A)
% A is skew-symmetric.
% Q is orthogonal and Q'*A*Q = T is tridiagonal

% Note that if M is skew-symmetric then
% (I - 2vv')M(I - 2vv') = M - 2vv'M - 2Mvv' + 4v(v'Mv)v'
% = M + vz' - zv'
% where z =2* M*v

[n,n] = size(A);
T = zeros(n,n);
for k=1:n-2
% Generate the kth Householder vector and store in A(k+1:n,k)...
[v,alpha] = House(A(k+1:n,k));
A(k+1:n,k) = v;
% Generate the new T-entry and update the rest of A...
T(k+1,k) = alpha;
T(k,k+1) = -alpha;
z = 2*(A(k+1:n,k+1:n)*v);
A(k+1:n,k+1:n) = A(k+1:n,k+1:n) + v*z' - z*v';
end
T(n,n-1) = A(n,n-1);
T(n-1,n) = A(n-1,n);
Q = eye(n,n);
% Q = Q1*Q2*...*Qn2
for k=n-2:-1:1
% Q = Qk*Q ...
v = A(k+1:n,k);
Q(k+1:n,k+1:n) = Q(k+1:n,k+1:n) - (2*v)*(v'*Q(k+1:n,k+1:n));
end

function [v,alpha] = House(x)
% x is an n-vector
% v is a unit vector so that if y = (I - 2*v*v')x then y(2:n) = 0
% and y(1) = alpha. Returning alpha makes the implementation of SkewReduce
% a little simpler.
v = x;
if v(1)< 0
alpha = norm(x);
else
alpha = -norm(x);
end
v(1) = v(1)-alpha;
v = v/norm(v);