function [ a, a_rank ] = i4mat_rref2 ( m, n1, n2, a )

%*****************************************************************************80
%
%% I4MAT_RREF2 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.
%
%    Some algorithms augment the matrix A with additional columns.
%    These might represent right hand side vectors of the system A*x=b,
%    or they might be a copy of the identity matrix.  In either case,
%    the intention is that the reduction operations will be applied to
%    the elements of these additional columns, but the pivot operation
%    will not proceed beyond the last column of the original matrix A.
%    This function allows for this procedure by specifying that the
%    pivot matrix extends from column 1 to column N1, while there may
%    be additional columns labeled N1+1 through N1+N2.
%
%  Example:
%
%    M = 4, N1 = 7, N2 = 0
%
%    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 GNU LGPL license.
%
%  Modified:
%
%    19 September 2018
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Charles Cullen,
%    An Introduction to Numerical Linear Algebra,
%    PWS Publishing Company, 1994,
%    ISBN: 978-0534936903,
%    LC: QA185.D37.C85.
%
%  Parameters:
%
%    Input, integer M, N1, N2, the number of rows and columns.
%
%    Input, integer A(M,1:N1+N2), the matrix to be analyzed. 
%
%    Output, integer A(M,1:N1+N2), the IRREF of the matrix.
%
  n = n1 + n2;
  a_rank = 0;

  lead = 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;
%
%  Finish if we completed column N1.
%
        if ( n1 < lead )
          lead = -1;
          break
        end
      end

    end

    if ( lead < 0 )
      break
    end
%
%  Increase rank by 1.
%
    a_rank = a_rank + 1;
%
%  Move pivot I into row R.
%
    if ( i ~= r )
      a = i4mat_row_swap ( m, n, a, i, r );
    end
%
%  Ensure pivot is positive.
%
    if ( a(r,lead) < 0 )
      a(r,1:n) = - a(r,1:n);
    end
%
%  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

