function [Q,T] = SkewReduce(A)
% A is an n-by-n skew-symmetric matrix
% Q is orthogonal so that Q'*A*Q = T is tridiagonal

% Note that if C is skew-symmetric then
%   (I - 2vv')C(I - 2vv') = C - (2v'C) - (2Cv)v' + 4v(v'Cv)v'
%                         = C + vw' - wv'
% where w = C*(2*v). (v'*Cv = 0)
[n,n] = size(A);
Q = eye(n,n);
T = zeros(n,n);
for k=1:n-2
    [v,alfa] = House(A(k+1:n,k));
    w = A(k+1:n,k+1:n)*(2*v);
    % This update is not quite optimal but good enough....
    A(k+1:n,k+1:n) = A(k+1:n,k+1:n) + v*w' - w*v';
    Q(:,k+1:n)     = Q(:,k+1:n) - (Q(:,k+1:n)*v)*(2*v)';
    T(k+1,k)       =  alfa;
    T(k,k+1)       = -alfa;
end
T(n,n-1) = A(n,n-1);
T(n-1,n) = -T(n,n-1);
