function x = gs_solve ( Q, R, b )
%
%  function x = gs_solve ( Q, R, b )
%
%  GS_SOLVE solves a linear system factored by GS_FACTOR.
%
%  Discussion:
%
%    The Gram-Schmidt QR factorization of an M by N matrix A has the form
%
%      A = Q * R
%
%    where Q is M by N with orthogonal columns, and R is N 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,N), R(N,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.00001;

[m,n] = size ( R );
%
%  Do this so X is a column vector.
%
x = zeros(n,1);
%
%  Compute the modified right hand side.
%
qtb = Q' * b;
%
%  Solve the upper triangular system.
%
for j = 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
