function [v, beta] = house(x)
%HOUSE Householder matrix.
%   [V, BETA] = GALLERY('HOUSE',X) takes X, an N-element column vector,
%   and returns V and BETA such that EYE(N) - BETA*V*V'
%   is a Householder matrix. A Householder matrix H satisfies the
%   relationship
%      H*X = -sign(X(1))*norm(X)*e1,
%   where e1 is the first column of EYE(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, third edition,
%       Johns Hopkins University Press, Baltimore, Maryland, 1996, Sec. 5.1.
%   [2] N. J. Higham, Accuracy and Stability of Numerical Algorithms,
%       Society for Industrial and Applied Mathematics, Philadelphia,
%       1996, Chap. 18.
%   [3] G. W. Stewart, Introduction to Matrix Computations, Academic Press,
%       New York, 1973, pp. 231-234, 262.
%   [4] J. H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University
%       Press, 1965, pp. 48-50.
%
%   Nicholas J. Higham, Dec 1999.
%   Copyright 1984-2000 The MathWorks, Inc. 
%   $Revision: 1.9 $  $Date: 2000/06/01 16:46:23 $

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

s = norm(x) * mysign(x(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.
