function [ Q, R ] = housefact ( A )
%
%  function [ Q, R ] = housefact ( A )
%
%  HOUSEFACT computes the QR factorization of a matrix using Householder reflectors.
%
%  Discussion:
%
%    The QR factorization of a matrix A has the form
%
%      A = Q * R
%
%    where Q is orthogonal, and R is upper triangular.
%
%  Modified:
%
%    29 February 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A(N,N), the matrix whose QR factorization is desired.
%
%    Output, real Q(N,N), R(N,N), the QR factors of A.
%
[n,n] = size ( A );
Q = eye ( n );
R = A;

for k = 1 : n-1
%
%  Extract the diagonal and subdiagonal elements of column K of R.
%
  d = zeros ( n, 1 );
  d(k:n) = R(k:n, k );
%
%  Compute the appropriate vector V.
%
  alpha = - sign ( d(k) ) * norm ( d(k:n) );

  if ( alpha ~= 0.0 ) 

    v = zeros ( n, 1 );
    v(k) = sqrt ( 0.5 * ( 1.0 - d(k) / alpha ) );
    v(k+1:n) = - 0.5 * d(k+1:n) / ( alpha * v(k) );  
%
%  Define the reflector matrix H(v).
%
    H = reflector ( v );
%
%  Push H' into the back of Q, and H into the front of R.
%
    Q = Q * H';
    R = H * R;

  end

end
