function [ a, det ] = i4mat_rref ( m, n, a )

%*****************************************************************************80
%
%% i4mat_rref() computes the reduced row echelon form of an I4MAT.
%
%  Discussion:
%
%    If a matrix A contains only integer entries, then when it is transformed
%    to row reduced echelon form, it is likely that many entries will no longer
%    be integers, due to the elimination process.
%
%    In some cases, tiny arithmetic errors in this elimination process can
%    result in spurious, tiny nonzero values which can invalidate the
%    calculation, particular if the elimination is being done in an effort
%    to determine the rank of the matrix.  These serious errors can easily
%    occur in very small matrices, such as of size 7x10.
%
%    If we, instead, insist on using only integer operations on an integer
%    matrix, we can guarantee that tiny roundoff errors will not cause
%    such problems.  On the other hand, as the elimination process proceeds,
%    we may instead calculate integer matrix entries of increasingly
%    large, and then ultimately meaningless magnitude.  I imagine this is 
%    likely to happen for moderate size matrices of order 50x50, say, but
%    this is a huge improvement over the unreliability of the real
%    arithmetic case.
%
%
%    Thus, we define "integer row reduced echelon form" (IRREF):
%
%
%    A matrix is in integer row reduced echelon form if:
%
%    * The leading nonzero in each row is positive.
%
%    * Each row has no common factor greater than 1.
%
%    * The leading nonzero in each row occurs in a column to
%      the right of the leading nonzero in the previous row.
%
%    * Rows which are entirely zero occur last.
%
%    * When a row contains a leading nonzero in column J, then column J
%      is otherwise entirely zero.
%
%  Example:
%
%    Input matrix:
%
%     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
%
%    Output matrix:
%
%     1    3    0    0    2    0    0
%     0    0    0    1    2    0    0
%     0    0    0    0    0    3    1
%     0    0    0    0    0    0    0
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    25 October 2020
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Charles Cullen,
%    An Introduction to Numerical Linear Algebra,
%    PWS Publishing Company, 1994,
%    ISBN: 978-0534936903,
%    LC: QA185.D37.C85.
%
%  Input:
%
%    integer M, N, the number of rows and columns.
%
%    integer A(M,N), the matrix to be analyzed. 
%
%  Output:
%
%    integer A(M,N), the IRREF of the matrix.
%
%    integer DET, the pseudo-determinant.
%
  if ( ~ i4mat_is_integer ( a ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_REF - Fatal error!' );
    fprintf ( 1, '  Input matrix A is not integral.\n' );
    error ( 'I4MAT_REF - Fatal error!' );
  end

  lead = 1;
  det = 1;
 
  for r = 1 : m

    if ( n < lead )
      break
    end
%
%  Start I at row R, and search for nonzero pivot entry A(I,LEAD).
%
    i = r;

    while ( a(i,lead) == 0.0 )

      i = i + 1;
%
%  If reach last row, reset I to R, and increment LEAD.
%
      if ( m < i )
        i = r;
        lead = lead + 1;
%
%  If reach last column, we can find no more pivots.
%
        if ( n < lead )
          lead = -1;
          break
        end
      end

    end

    if ( lead < 0 )
      break
    end
%
%  Move pivot I into row R.
%
    if ( i ~= r )
      temp     = a(i,1:n);
      a(i,1:n) = a(r,1:n);
      a(r,1:n) = temp;
    end
%
%  Ensure pivot is positive.
%
    if ( a(r,lead) < 0 )
      a(r,1:n) = - a(r,1:n);
      det = - det;
    end
%
%  Update the pseudo-determinant.
%
    det = det * a(r,lead);
%
%  Remove any common factor from row R.
%
    [ a(r,1:n), ifact ] = i4vec_red ( n, a(r,1:n), 1 );
%
%  Use a multiple of A(R,LEAD) to eliminate A(1:M,LEAD).
%
    for i = 1 : m

      if ( i ~= r )

        a(i,1:n) = a(r,lead) * a(i,1:n) - a(i,lead) * a(r,1:n);

        [ a(i,1:n), ifact ] = i4vec_red ( n, a(i,1:n), 1 );

      end

    end

    lead = lead + 1;

  end

  return
end

