function i4mat_rref2_test ( )

%*****************************************************************************80
%
%% i4mat_rref2_test() tests i4mat_rref2().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    02 February 2019
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../i4mat_rref2' );

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'i4mat_rref2_test():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  i4mat_rref2() computes the reduced row echelon form of\n' );
  fprintf ( 1, '  an augmented integer matrix.\n' );
%
%  Case 1, just get the IRREF.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Test 1, just get the IRREF.\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 ];

  i4mat_print ( m, n, a, '  Matrix A:' );
  [ a_rref, a_rank ] = i4mat_rref2 ( m, n1, n2, a );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  A has rank %d\n', a_rank );
  i4mat_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;
  n1 = 4;

  a1 = [ ...
    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 ];

  i4mat_print ( m, n1, a1, '  Matrix A:' );

  b1 = [ 1;  2; 3; 4;  5; 5; 7 ];

  i4vec_print ( m, b1, '  RHS b1:' );

  b2 = [ 4; 12; 0; 2; 12; 9; 3 ];
  i4vec_print ( m, b2, '  RHS b2:' );

  a = a1;
  n1 = 4;
  n2 = 0;
  [ a_rref, a_rank0 ] = i4mat_rref2 ( m, n1, n2, a );
  fprintf ( 1, '  A has rank %d\n', a_rank );
  i4mat_print ( m, n1, a_rref, '  RREF of A' );

  a = [ a1, b1 ];
  n1 = 4;
  n2 = 1;
  [ a_rref, a_rank1 ] = i4mat_rref2 ( m, n1+n2, 0, a );
  fprintf ( 1, '  A|b1 has rank %d\n', a_rank1 );
  i4mat_print ( m, n1+n2, a_rref, '  RREF of A|b1' );
  if ( a_rank < a_rank1 )
    fprintf ( 1, '  Linear system A*x=b1 is INCONSISTENT.\n' );
  else
    fprintf ( 1, '  Linear system A*x=b1 is CONSISTENT.\n' );
    b_rref = a_rref(1:m,n1+1);
    i4vec_print ( m, b_rref, '  Right hand side of A_rref*x_rref=b1_rref:' );
    x_rref = i4mat_rref2_solve ( m, n1, n2, a_rref );
    r8vec_print ( n1, x_rref, '  Solution of A*x=b1:' );
  end

  a = [ a1, b2 ];
  n1 = 4;
  n2 = 1;
  [ a_rref, a_rank2 ] = i4mat_rref2 ( m, n1+n2, 0, a );
  fprintf ( 1, '  A|b2 has rank %d\n', a_rank2 );
  i4mat_print ( m, n1+n2, a_rref, '  RREF of A|b2' );
  if ( a_rank < a_rank2 )
    fprintf ( 1, '  Linear system A*x=b1 is INCONSISTENT.\n' );
  else
    fprintf ( 1, '  Linear system A*x=b1 is CONSISTENT.\n' );
    b_rref = a_rref(1:m,n1+n2);
    i4vec_print ( m, b_rref, '  Right hand side of A_rref*x_rref=b1_rref:' );
    x_rref = i4mat_rref2_solve ( m, n1, n2, a_rref );
    r8vec_print ( n1, x_rref, '  Solution of A*x=b2:' );
  end
%
%  Case 3, inverse
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Test 3, inverse.\n' );
  m = 4;
  n1 = 4;
  n2 = 4;

  a1 = [ ...
     5,  7,  6,  5; ...
     7, 10,  8,  7; ...
     6,  8, 10,  9; ...
     5,  7,  9, 10 ];

  a2 = eye ( m, n2 );
  a = [ a1, a2 ];
  [ a_rref, a_rank ] = i4mat_rref2 ( m, n1, n2, a );
  a_inv = a_rref(1:m,n1+1:n1+n2);
  p = a_inv * a;
  r8mat_print ( m, n1, p, '  A_inv * A:' );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'i4mat_rref2_test():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../i4mat_rref2' );

  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

