function [ B ] = do_householder( u, A )
%DO_HOUSEHOULDER( u, A ) Apply a specified Householder transformation 
%   A is the (square!) input matrix.
%   u is the normal to the plane of the Householder reflection.
%   B (returned) is the transformed output.
 
% Apply the HH reflection
    normu=norm(u);
    if normu~=0;
        w = u'*A;
        B = A - 2.*(u*w)./((normu)^2);
    else
        B=A; % unchanged
    end
end

