function [ a2, b2, incon, freedom_num, freedom ] = i4mat_rref_system ( m, n, a, b )

%*****************************************************************************80
%
%% i4mat_rref_system() sets up an augmented IRREF linear system.
%
%  Discussion:
%
%    An MxN integer linear system A * X = B is considered.
% 
%    The matrix A and right hand side B are assumed to have been converted
%    to integer row-reduced echelon form (IRREF).
%
%    To create, if possible, a solvable NxN system, this function removes
%    trailing zero rows, and inserts where necessary, rows of the identity
%    matrix in A and zeros in B, corresponding to undetermined degrees of 
%    freedom, producing the NxN system:
%
%      A2 * X = B2
%
%    The function also indicates whether the initial system was inconsistent,
%    and identifies those rows of A2 that correspond to degrees of freedom.
%    
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    09 September 2018
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer M, N, the number of rows and columns of the IRREF matrix A.
%
%    integer A(M,N), the IRREF matrix to be analyzed. 
%
%    integer B(M), the IRREF right hand side.
%
%  Output:
%
%    integer A2(N,N), the modified IRREF matrix.
%
%    integer B2(N), the modified IRREF right hand side.
%
%    logical INCON, is TRUE if the system A*X=B is inconsistent.
%
%    integer FREEDOM_NUM, the number of degrees of freedom.
%    If FREEDOM_NUM < 0, then there are no degrees of freedom and the
%    system is overdetermined.
%
%    integer FREEDOM(FREEDOM_NUM), the indices of the degrees
%    of freedom, presuming 0 <= FREEDOM_NUM.
%

%
%  Determine 0 <= M2 <= M, the location of the last nonzero row in A.
%  If any zero row of A has a nonzero B, then the equations are inconsistent.
%
  m2 = m;
  incon = false;

  while ( 0 < m2 )

    if ( any ( a(m2,1:n) ~= 0 ) )
      break
    end

    if ( b(m2) ~= 0 )
      incon = true;
    end

    m2 = m2 - 1;

  end
%
%  Copy rows 1 through M2 of A and B into A2 and B2.
%
  a2 = zeros ( n, n );
  b2 = zeros ( n, 1 );
  a2(1:m2,1:n) = a(1:m2,1:n);
  b2(1:m2,1)   = b(1:m2,1);
%
%  Count the indeterminate variables.
%
  freedom_num = n - m2;
%
%  If pivot in column J is missing,
%  modify matrix and right hand side.
%  Add J to list of indeterminate variables.
%
  freedom = [];

  if ( 0 < freedom_num )

    for j = 1 : n
      if ( m2 < j )
        a2(j,1:n) = 0;
        a2(j,j) = 1;
        b2(j) = 0;
        freedom = [ freedom; j ];
        m2 = m2 + 1;
      elseif ( a2(j,j) == 0 )
        a2(j+1:m2+1,1:n) = a2(j:m2,1:n);
        a2(j,1:n) = 0;
        a2(j,j) = 1;
        b2(j+1:m2+1,1) = b2(j:m2,1);
        b2(j) = 0;
        freedom = [ freedom; j ];
        m2 = m2 + 1;
      end
    end

  end

  return
end

