function [ Q, R ] = h_factor ( A )
%
%  function [ Q, R ] = h_factor ( A )
%
%  H_FACTOR computes the QR factorization of a matrix using Householder matrices.
%
%  Discussion:
%
%    The QR factorization of an M by N matrix A has the form
%
%      A = Q * R
%
%    where Q is M by M and orthogonal, and R is M by N and upper triangular.
%
%  Modified:
%
%    21 March 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A(M,N), the matrix whose QR factorization is desired.
%
%    Output, real Q(M,M), R(M,N), the QR factors of A.
%
[m,n] = size ( A );
%
%  Initialize.
%
Q = eye ( m, m );
R = A;

khi = min ( m-1, n );

for k = 1 : khi
%
%  Get the Householder matrix that zeros out the K-th subdiagonal of R.
%
  H = householder ( R, k );
%
%  Push H' into the back of Q, and H into the front of R.
%
  Q = Q * H';
  R = H * R;

end
