function rref_test ( )

%*****************************************************************************80
%
%% RREF_TEST tests RREF.
%
%  Modified:
%
%    28 September 2018
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'RREF_TEST:\n' );
  fprintf ( 1, '  RREF computes the reduced row echelon form (RREF) of\n' );
  fprintf ( 1, '  a matrix.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The RREF is susceptible to roundoff error.\n' );
  fprintf ( 1, '  With that in mind, it can be used to determine:\n' );
  fprintf ( 1, '  * the rank of a square matrix;\n' );
  fprintf ( 1, '  * the column rank of a rectangular matrix;\n' );
  fprintf ( 1, '  * the independent columns of a rectangular matrix\n' );
  fprintf ( 1, '  * if a linear system is consistent;\n' );
  fprintf ( 1, '  * the inverse of a square invertible matrix.\n' );
%
%  Case 1, just get the RREF.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Test 1, just get the RREF.\n' );

  m = 4;
  n1 = 7;
  n2 = 0;
  n = n1 + n2;
  a = [ ...
     1,  3,  0,  2,  6,  3,  1; ...
    -2, -6,  0, -2, -8,  3,  1; ...
     3,  9,  0,  0,  6,  6,  2; ...
    -1, -3,  0,  1,  0,  9,  3 ];

  r8mat_print ( m, n, a, '  Matrix A:' );
  [ a_rref, a_cols ] = rref ( a );
  a_rank = length ( a_cols );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  A has rank %d\n', a_rank );
  r8mat_print ( m, n, a_rref, '  RREF of A:' );
%
%  Case 2, look for inconsistency.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Test 2, look for inconsistency.\n' );

  m = 7;
  n = 4;

  a = [ ...
    1, -2, 3, -1; ...
    3, -6, 9, -3; ...
    0,  0, 0,  0; ...
    2, -2, 0,  1; ...
    6, -8, 6,  0; ...
    3,  3, 6,  9; ...
    1,  1, 2,  3 ];

  r8mat_print ( m, n, a, '  Matrix A:' );
  [ a_rref, a_cols ] = rref ( a );
  a_rank = length ( a_cols );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  A has rank %d\n', a_rank );

  b1 = [ 1;  2; 3; 4;  5; 5; 7 ];
  r8vec_print ( m, b1, '  RHS b1:' );
  ab1 = [ a, b1 ];
  [ ab1_rref, ab1_cols ] = rref ( ab1 );
  ab1_rank = length ( ab1_cols );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  [A;b1] has rank %d\n', ab1_rank );
  if ( a_rank < ab1_rank )
    fprintf ( 1, '  Linear system A*x=b1 is INCONSISTENT.\n' );
  else
    fprintf ( 1, '  Linear system A*x=b1 is CONSISTENT.\n' );
    x1 = rref_solve ( a, b1 );
    r8vec_print ( n, x1, '  Solution of A*x1=b1:' );
    ax1 = a * x1;
    r8vec2_print ( m, b1, ax1, '  Compare b1 and A*x1:' );
  end

  b2 = [ 4; 12; 0; 2; 12; 9; 3 ];
  r8vec_print ( m, b2, '  RHS b2:' );
  ab2 = [ a, b2 ];
  [ ab2_rref, ab2_cols ] = rref ( ab2 );
  ab2_rank = length ( ab2_cols );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  [A;b2] has rank %d\n', ab2_rank );
  if ( a_rank < ab2_rank )
    fprintf ( 1, '  Linear system A*x2=b2 is INCONSISTENT.\n' );
  else
    fprintf ( 1, '  Linear system A*x2=b2 is CONSISTENT.\n' );
    x2 = rref_solve ( a, b2 );
    r8vec_print ( n, x2, '  Solution of A*x2=b2:' );
    ax2 = a * x2;
    r8vec2_print ( m, b2, ax2, '  Compare b2 and A*x2:' );
  end
%
%  Case 3, inverse
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Test 3, inverse.\n' );
  n = 4;

  a = [ ...
     5,  7,  6,  5; ...
     7, 10,  8,  7; ...
     6,  8, 10,  9; ...
     5,  7,  9, 10 ];
  r8mat_print ( n, n, a, '  Matrix A:' );

  a_inv = rref_inverse ( a );
  r8mat_print ( n, n, a_inv, '  Estimated A_inv:' );
  p = a_inv * a;
  r8mat_print ( n, n, p, '  A_inv * A:' );
%
%  Case 4, independent columns
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Test 4, independent columns.\n' );

  m = 7;
  n = 4;

  a = [ ...
        1, 1,  2, 3; ...
        2, 3,  4, 9; ...
        3, 0,  6, 0; ...
        4, 2,  8, 0; ...
        5, 6, 10, 6; ...
        6, 3, 12, 6; ...
        7, 1, 14, 2 ];

  r8mat_print ( m, n, a, '  Matrix A:' );

  a_columns = rref_columns ( a );
  [ m, n2 ] = size ( a_columns );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of independent columns is %d\n', n2 );

  r8mat_print ( m, n2, a_columns, '  Independent columns of A:' );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'RREF_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


