function x = svd_solve ( m, n, a, b )

%*****************************************************************************80
%
%% SVD_SOLVE solves a linear system in the least squares sense.
%
%  Discussion:
%
%    The vector X returned by this routine should always minimize the 
%    Euclidean norm of the residual ||A*x-b||.
%
%    If the matrix A does not have full column rank, then there are multiple
%    vectors that attain the minimum residual.  In that case, the vector
%    X returned by this routine is the unique such minimizer that has the 
%    the minimum possible Euclidean norm, that is, ||A*x-b|| and ||x||
%    are both minimized.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 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.
%
%    Input, real A(M,N), the matrix.
%
%    Input, real B(M), the right hand side.
%
%    Output, real X(N), the least squares solution.
%

%
%  Get the SVD.
%
  a_copy = a;
  lda = m;
  ldu = m;
  ldv = n;
  job = 11;

  [ a_copy, sdiag, e, u, v, info ] = dsvdc ( a_copy, lda, m, n, ldu, ldv, job );

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SVD_SOLVE - Fatal error!\n' );
    fprintf ( 1, '  The SVD could not be calculated.\n' );
    fprintf ( 1, '  LINPACK routine DSVDC returned a nonzero\n' );
    fprintf ( 1, '  value of the error flag, INFO = %d\n', info );
    error ( 'SVD_SOLVE - Fatal error!' );
  end

  ub = u' * b;

  sub = zeros ( n, 1 );
%
%  For singular problems, there may be tiny but nonzero singular values
%  that should be ignored.  This is a reasonable attempt to avoid such 
%  problems, although in general, the user might wish to control the tolerance.
%
  smax = max ( sdiag );
  if ( smax <= eps )
    smax = 1.0
  end

  stol = eps * smax;

  for i = 1 : n
    if ( i <= m )
      if ( stol <= sdiag(i) )
        sub(i) = ub(i) / sdiag(i);
      end
    end
  end

  x = v * sub;

  return
end
