function x = h_solve ( Q, R, b )
%
%  function x = h_solve ( Q, R, b )
%
%  H_SOLVE solves a linear system factored by H_FACTOR.
%
%  Discussion:
%
%    The Householder 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.
%
%    The solution procedure is not sophisticated.  It tries to handle
%    the case where a diagonal value of R is zero.
%
%    If the matrix A is square and has full rank, then X is a solution
%    in the usual sense.
%
%    If the matrix A has more columns than rows, then the last N-M
%    entries of X are simply forced to be zero.  Assuming the matrix
%    has full rank, then X is one of many solutions.
%
%    If the matrix A has more rows than columns, there is probably no
%    solution with zero residual.  The procedure returns a solution
%    in the least squares sense, which minimizes the L2 norm of
%    the residual.
%
%  Modified:
%
%    27 March 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real Q(M,M), R(M,N), the QR factors of A.
%
%    Input, real B(M), the right hand side of the linear system.
%
%    Output, real X(N), the QR solution of the linear system.
%
TOL = 0.000001;

[m,n] = size ( R );
%
%  If you don't zero out X this way, and M < N, then entries of X
%  won't be set.
%
x = zeros (n,1);
%
%  Compute the modified right hand side.
%
qtb = Q' * b;
%
%  Solve the upper triangular system.
%
for j = min ( m, n ) : -1 : 1

  if ( abs ( R(j,j) ) <= 0.0 ) 
    x(j) = 0.0;
  else
    x(j) = qtb(j) / R(j,j);
    for i = 1 : j-1
      qtb(i) = qtb(i) - R(i,j) * x(j);
    end
  end

end
