function eros_test ( )

%*****************************************************************************80
%
%% eros_test() tests eros().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    08 January 2020
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../eros' );

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'eros_test():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  Test eros().\n' );

  gauss_test1 ( );
  gauss_test2 ( );
  gauss_det_test ( );
  gauss_inverse_test ( );
  gauss_plu_test ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'eros_test():\n' );
  fprintf ( 1, '  Normal end of execution\n' );
  timestamp ( );

  rmpath ( '../eros' );

  return
end
function gauss_test1 ( )

%*****************************************************************************80
%
%% gauss_test1() tests gauss() with a random 3x3 system.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    01 January 2020
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'gauss_test1():\n' );
  fprintf ( 1, '  gauss() solves a 3x3 linear system with one right hand side.\n' );

  A = rand ( 3, 3 );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix A:\n' );
  fprintf ( 1, '\n' );
  disp ( A );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Exact solution x:\n' );
  fprintf ( 1, '\n' );
  x = [ 2.0; 1.0; 3.0 ];
  disp ( x );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Right hand side b:\n' );
  fprintf ( 1, '\n' );
  b = A * x;
  disp ( b );

  wait = 1;
  gauss ( A, b, wait );

  return
end
function gauss_test2 ( )

%*****************************************************************************80
%
%% gauss_test2() tests gauss with a random 3x3 system and two right hand sides.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    01 January 2020
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'gauss_test2():\n' );
  fprintf ( 1, '  gauss solves a random 3x3 linear system\n' );
  fprintf ( 1, '  with two right hand sides.\n' );

  A = rand ( 3, 3 );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix A:\n' );
  fprintf ( 1, '\n' );
  disp ( A );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Exact solution x:\n' );
  fprintf ( 1, '\n' );
  x = [ 2.0, 10.0; ...
        1.0,  1.5; ...
        3.0, -2.0 ];
  disp ( x );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Right hand side b:\n' );
  fprintf ( 1, '\n' );
  b = A * x;
  disp ( b );

  wait = 1;
  gauss ( A, b, wait );

  return
end
function gauss_det_test ( )

%*****************************************************************************80
%
%% gauss_det_test() tests gauss_det().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    08 January 2020
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'gauss_det_test():\n' );
  fprintf ( 1, '  gauss_det() uses Gauss elimination to find the determinant\n' );
  fprintf ( 1, '  of a matrix.\n' );

  A = [ 1, 2, 3; 
        4, 5, 6;
        7, 8, 9 ];
  value = gauss_det ( A );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix A:\n' );
  fprintf ( 1, '\n' );
  disp ( A );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Computed determinant = %g:\n', value );
  fprintf ( 1, '  MATLAB det(A) =        %g:\n', det ( A ) );

  A = golub_matrix ( 5 );
  value = gauss_det ( A );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix A:\n' );
  fprintf ( 1, '\n' );
  disp ( A );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Computed determinant = %g:\n', value );
  fprintf ( 1, '  MATLAB det(A) =        %g:\n', det ( A ) );

  return
end
function gauss_inverse_test ( )

%*****************************************************************************80
%
%% gauss_inverse_test() tests gauss_inverse().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    08 January 2020
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'gauss_inverse_test():\n' );
  fprintf ( 1, '  gauss_inverse() uses Gauss elimination to find the\n' );
  fprintf ( 1, '  inverse of a matrix.\n' );

  A = [ 1, 2, 3; 
        4, 5, 8;
        7, 8, 9 ];
  B = gauss_inverse ( A );
  I = eye ( size ( A ) );
  AB = A * B;
  r = norm ( AB - I );
  B2 = inv ( A );
  e = norm ( B2 - B );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix A:\n' );
  fprintf ( 1, '\n' );
  disp ( A );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Computed inverse B:\n' );
  fprintf ( 1, '\n' );
  disp ( B );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  MATLAB inverse B2 = inv(A):\n' );
  fprintf ( 1, '\n' );
  disp ( B2 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Residual norm ||A*B-I|| = %g:\n', r );
  fprintf ( 1, '  Error norm    ||B2-B|| =  %g:\n', e );

  A = golub_matrix ( 5 );
  B = gauss_inverse ( A );
  I = eye ( size ( A ) );
  AB = A * B;
  r = norm ( AB - I );
  B2 = inv ( A );
  e = norm ( B2 - B );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix A:\n' );
  fprintf ( 1, '\n' );
  disp ( A );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Computed inverse B:\n' );
  fprintf ( 1, '\n' );
  disp ( B );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  MATLAB inverse B2 = inv(A):\n' );
  fprintf ( 1, '\n' );
  disp ( B2 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Residual norm ||A*B-I|| = %g:\n', r );
  fprintf ( 1, '  Error norm    ||B2-B|| =  %g:\n', e );

  return
end
function gauss_plu_test ( )

%*****************************************************************************80
%
%% gauss_plu_test() tests gauss_plu().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    08 January 2020
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'gauss_plu_test():\n' );
  fprintf ( 1, '  gauss_plu() uses Gauss elimination to find the\n' );
  fprintf ( 1, '  PLU factors of a matrix.\n' );

  A = [ 1, 2, 3; 
        4, 5, 8;
        7, 8, 9 ];
  [ P, L, U ] = gauss_plu ( A );
  PLU = P' * L * U;
  r = norm ( A - PLU );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix A:\n' );
  fprintf ( 1, '\n' );
  disp ( A );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Permutation matrix P:\n' );
  fprintf ( 1, '\n' );
  disp ( P );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Unit lower triangular matrix L:\n' );
  fprintf ( 1, '\n' );
  disp ( L );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Upper triangular matrix U:\n' );
  fprintf ( 1, '\n' );
  disp ( U );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Residual norm ||A-P''*L*U|| = %g:\n', r );

  A = golub_matrix ( 5 );
  [ P, L, U ] = gauss_plu ( A );
  PLU = P' * L * U;
  r = norm ( A - PLU );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix A:\n' );
  fprintf ( 1, '\n' );
  disp ( A );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Permutation matrix P:\n' );
  fprintf ( 1, '\n' );
  disp ( P );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Unit lower triangular matrix L:\n' );
  fprintf ( 1, '\n' );
  disp ( L );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Upper triangular matrix U:\n' );
  fprintf ( 1, '\n' );
  disp ( U );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Residual norm ||A-P''*L*U|| = %g:\n', r );

  return
end
function A = golub_matrix ( n )

%*****************************************************************************80
%
%% golub_matrix() returns the golub matrix.
%
%  Discussion:
%
%    These matrices are the product of random unit lower and unit upper 
%    triangular matrices.
%
%    These matrices tend to be badly conditioned.
%
%  Properties:
%
%    A can be LU factored without pivoting.
%
%    det(A) = 1.
%
%    For values of n greater than 10, the determinant cannot may not be 
%    reliably computed.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    07 January 2020
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Cleve Moler,
%    Numerical Computing with MATLAB,
%    SIAM, 2004,
%    ISBN13: 978-0-898716-60-3,
%    LC: QA297.M625
%
%  Input:
%
%    integer N, the order of the matrix.
%
%  Output:
%
%    real A(N,N), the matrix.
%
  s = 10.0;

  L = eye ( n, n );
  for i = 1 : n
    for j = 1 : i - 1
      L(i,j) = s * randn ( );
    end
  end

  U = eye ( n, n );
  for i = 1 : n
    for j = i + 1 : n
      U(i,j) = s * randn ( );
    end
  end

  A = L * U;

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% timestamp() prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

