function test_lu ( )

%*****************************************************************************80
%
%% test_lu() tests the LU factors.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    27 October 2021
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'test_lu():\n' );
  fprintf ( 1, '  A = a test matrix of order M by N\n' );
  fprintf ( 1, '  L, U are the LU factors.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  |A| = Frobenius norm of A.\n' );
  fprintf ( 1, '  |A-LU| = Frobenius norm of A-L*U.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Title                    M     N      ' );
  fprintf ( 1, '|A|            |A-LU|\n' );
  fprintf ( 1, '\n' );
%
%  bodewig
%
  title = 'bodewig';
  m = 4;
  n = 4;
  a = bodewig_matrix ( );
  [ l, u ] = bodewig_lu ( );
  error_frobenius = r8mat_is_lu ( m, n, a, l, u );
  norm_a_frobenius = norm ( a, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  borderband
%
  title = 'borderband';
  m = 5;
  n = 5;
  a = borderband_matrix ( n );
  [ l, u ] = borderband_lu ( n );
  error_frobenius = r8mat_is_lu ( m, n, a, l, u );
  norm_a_frobenius = norm ( a, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  dif2
%
  title = 'dif2';
  m = 5;
  n = 5;
  a = dif2_matrix ( n );
  [ l, u ] = dif2_lu ( n );
  error_frobenius = r8mat_is_lu ( m, n, a, l, u );
  norm_a_frobenius = norm ( a, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  fibonacci2
%
  title = 'fibonacci2';
  m = 5;
  n = 5;
  a = fibonacci2_matrix ( n );
  [ l, u ] = fibonacci2_lu ( n );
  error_frobenius = r8mat_is_lu ( m, n, a, l, u );
  norm_a_frobenius = norm ( a, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  gfpp
%
  title = 'gfpp';
  m = 5;
  n = 5;
  alpha = 1.0;
  a = gfpp_matrix ( n, alpha );
  [ l, u ] = gfpp_lu ( n, alpha );
  error_frobenius = r8mat_is_lu ( m, n, a, l, u );
  norm_a_frobenius = norm ( a, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  givens
%
  title = 'givens';
  m = 5;
  n = 5;
  a = givens_matrix ( m, n );
  [ l, u ] = givens_lu ( n );
  error_frobenius = r8mat_is_lu ( m, n, a, l, u );
  norm_a_frobenius = norm ( a, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  golub
%
  title = 'golub';
  m = 5;
  n = 5;
  key = 123456789;
  a = golub_matrix ( n, key );
  [ l, u ] = golub_lu ( n, key );
  error_frobenius = r8mat_is_lu ( m, n, a, l, u );
  norm_a_frobenius = norm ( a, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  kms
%
  title = 'kms';
  m = 5;
  n = 5;
  alpha = 0.5;
  a = kms_matrix ( alpha, m, n );
  [ l, u ] = kms_lu ( alpha, n );
  error_frobenius = r8mat_is_lu ( m, n, a, l, u );
  norm_a_frobenius = norm ( a, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  lehmer
%
  title = 'lehmer';
  m = 5;
  n = 5;
  a = lehmer_matrix ( m, n );
  [ l, u ] = lehmer_lu ( n );
  error_frobenius = r8mat_is_lu ( m, n, a, l, u );
  norm_a_frobenius = norm ( a, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  minij
%
  title = 'minij';
  m = 5;
  n = 5;
  a = minij_matrix ( m, n );
  [ l, u ] = minij_lu ( n );
  error_frobenius = r8mat_is_lu ( m, n, a, l, u );
  norm_a_frobenius = norm ( a, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  moler1
%
  title = 'moler1';
  m = 5;
  n = 5;
  alpha = 2.0;
  a = moler1_matrix ( alpha, m, n );
  [ l, u ] = moler1_lu ( alpha, n );
  error_frobenius = r8mat_is_lu ( m, n, a, l, u );
  norm_a_frobenius = norm ( a, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  moler3
%
  title = 'moler3';
  m = 5;
  n = 5;
  a = moler3_matrix ( m, n );
  [ l, u ] = moler3_lu ( n );
  error_frobenius = r8mat_is_lu ( m, n, a, l, u );
  norm_a_frobenius = norm ( a, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  oto
%
  title = 'oto';
  m = 5;
  n = 5;
  a = oto_matrix ( m, n );
  [ l, u ] = oto_lu ( n );
  error_frobenius = r8mat_is_lu ( m, n, a, l, u );
  norm_a_frobenius = norm ( a, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  pascal2
%
  title = 'pascal2';
  m = 5;
  n = 5;
  a = pascal2_matrix ( n );
  [ l, u ] = pascal2_lu ( n );
  error_frobenius = r8mat_is_lu ( m, n, a, l, u );
  norm_a_frobenius = norm ( a, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );
%
%  vand2
%
  title = 'vand2';
  m = 5;
  n = 5;
  x = 0.0 + 5.0 * randn ( n );
  a = vand2_matrix ( n, x );
  [ l, u ] = vand2_lu ( n, x );
  error_frobenius = r8mat_is_lu ( m, n, a, l, u );
  norm_a_frobenius = norm ( a, 'fro' );
  fprintf ( 1, '  %-20s  %4d  %4d  %14g  %14g\n', ...
    title, m, n, norm_a_frobenius, error_frobenius );

  return
end
