function [ x, flag ] = normal_solve ( m, n, a, b )

%*****************************************************************************80
%
%% NORMAL_SOLVE solves a linear system using the normal equations.
%
%  Discussion:
%
%    Given a presumably rectangular MxN system of equations A*x=b, this routine
%    sets up the NxN system A'*A*x=A'b.  Assuming N <= M, and that A has full
%    column rank, the system will be solvable, and the vector x that is returned
%    will minimize the Euclidean norm of the residual.
%
%    One drawback to this approach is that the condition number of the linear
%    system A'*A is effectively the square of the condition number of A, 
%    meaning that there is a substantial loss of accuracy.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 August 2016
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    David Kahaner, Cleve Moler, Steven Nash,
%    Numerical Methods and Software,
%    Prentice Hall, 1989,
%    ISBN: 0-13-627258-4,
%    LC: TA345.K34.
%
%  Parameters:
%
%    Input, integer M, the number of rows of A.
%
%    Input, integer N, the number of columns of A.
%    It must be the case that N <= M.
%
%    Input, real A(M,N), the matrix.
%    The matrix must have full column rank.
%
%    Input, real B(M), the right hand side.
%
%    Output, real X(N), the least squares solution.
%
%    Output, integer FLAG,
%    0, no error was detected.
%    1, an error occurred.
%
  flag = 0;
  x = zeros ( n, 1 );

  if ( m < n )
    flag = 1;
    return
  end

  ata = a' * a;

  atb = a' * b;

  [ ata_c, flag ] = r8mat_cholesky_factor ( n, ata );

  if ( flag ~= 0 )
    return
  end

  x = r8mat_cholesky_solve ( n, ata_c, atb );

  return
end
 
