function gradient_descent_linear_test ( )

%*****************************************************************************80
%
%% gradient_descent_linear_test() tests gradient_descent_linear().
%
%  Discussion:
%
%    This code suggests how slowly a simple gradient descent method will
%    converge to the solution of a linear least squares problem.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    27 July 2022
%
%  Author:
%
%    John Burkardt.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'gradient_descent_linear_test():\n' );
  fprintf ( 1, '  gradient_descent_linear() approximates the solution of a\n' );
  fprintf ( 1, '  least squares problem, to find a solution x to A*x=b,\n' );
  fprintf ( 1, '  by minimizing ||Ax-b|| using gradient descent.\n' );

  A = [ ...
   1.47, 1.0; ...
   1.50, 1.0; ...
   1.52, 1.0; ...
   1.55, 1.0; ...
   1.57, 1.0; ...
   1.60, 1.0; ...
   1.63, 1.0; ...
   1.65, 1.0; ...
   1.68, 1.0; ...
   1.70, 1.0; ...
   1.73, 1.0; ...
   1.75, 1.0; ...
   1.78, 1.0; ...
   1.80, 1.0; ...
   1.83, 1.0 ];

  [ m, n ] = size ( A );

  b = [ ...
    52.21; ...
    53.12; ...
    54.48; ...
    55.84; ...
    57.20; ...
    58.57; ...
    59.93; ...
    61.29; ...
    63.11; ...
    64.47; ...
    66.28; ...
    68.10; ...
    69.92; ...
    72.19; ...
    74.46 ];

  alpha = 0.02;
  tol = 1.0e-6;
  it_max = 10000;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Learning rate = %g\n', alpha );
  fprintf ( 1, '  Stepsize tolerance = %g\n', tol );
  fprintf ( 1, '  Maximum iterations = %d\n', it_max );

  x_exact = [ 61.272; -39.062 ];

  [ x, it ] = gradient_descent_linear ( A, b, alpha, tol, it_max );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of iterations = %d\n', it );
  disp ( '  Estimated solution:' );
  disp ( x );
  disp ( '  Exact solution:' );
  disp ( x_exact );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Error = %g\n', norm ( x - x_exact ) );
  fprintf ( 1, '\n' )
  fprintf ( 1, '  Residual for estimated solution = %g\n', norm ( b - A * x ) );
  fprintf ( 1, '  Residual for exact solution =     %g\n', norm ( b - A * x_exact ) );

  return
end

