function [v, beta] = house(x)
%HOUSE  Householder matrix.
%   [V, BETA] = GALLERY('HOUSE',X) takes X, a scalar or N-element 
%   column vector, and returns V and BETA such that EYE(N,N) - BETA*V*V'
%   is a Householder matrix. A Householder matrix H satisfies the 
%   following relationship: 
%      H*X = -sign(X(1))*norm(X)*e1, 
%   where e1 is the first column of EYE(N,N). 
%
%   If X is complex, then 
%       sign(X) = exp(i*arg(X)) 
%   which equals X./abs(X) when X ~= 0. 
%   If X = 0, then V = 0 and BETA = 1.

%   Theory:
%   Hx = y: (I - beta*v*v')x = -s*e1.
%   Must have |s| = norm(x), v = x+s*e1, and
%   x'y = x'Hx =(x'Hx)' real => arg(s) = arg(x(1)).
%   So take s = sign(x(1))*norm(x) (which avoids cancellation).
%   v'v = (x(1)+s)^2 + x(2)^2 + ... + x(n)^2
%       = 2*norm(x)*(norm(x) + |x(1)|).
%
%   References:
%   [1] G.H. Golub and C.F. Van Loan, Matrix Computations, second edition,
%   Johns Hopkins University Press, Baltimore, Maryland, 1989, pp.38-43.
%   [2] G.W. Stewart, Introduction to Matrix Computations, Academic Press,
%   New York, 1973, pp. 231-234, 262.
%   [3] J.H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University
%   Press, 1965, pp. 48-50.
%
%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 1.6 $  $Date: 1997/11/21 23:29:36 $

[m, n] = size(x);
if n > 1, error('Argument must be a column vector.'), end

s = norm(x) * (sign(x(1)) + (x(1)==0));    % Modification for sign(0)=1.
v = x;
if s == 0, beta = 1; return, end           % Quit if x is the zero vector.
v(1) = v(1) + s;
beta = 1/(s'*v(1));                        % NB the conjugated s.

% beta = 1/(abs(s)*(abs(s)+abs(x(1)) would guarantee beta real.
% But beta as above can be non-real (due to rounding) only when x is complex.
