function AINV = rref_inverse ( A )

%*****************************************************************************80
%
%% rref_inverse() uses reduced row echelon form (RREF) to compute an inverse.
%
%  Discussion:
%
%    The procedure will fail if A is not square, or not invertible.
%
%    This is simply a demonstration of how RREF can be used to compute
%    the inverse.  But:
%    * there are better ways to compute the inverse, including
%      B = inv ( A )
%    * the inverse matrix is usually not the appropriate tool for solving
%      linear systems.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    05 December 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real A(N,N), a square, invertible matrix.
%
%  Output:
%
%    real AINV(N,N), the inverse of A, computed using the rref_compute().
%
  [ m, n ] = size ( A );

  if ( m ~= n )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'rref_inverse(): Fatal error!\n' );
    fprintf ( 1, '  Input matrix is not square.\n' );
    error ( 'rref_inverse(): Fatal error!\n' );
  end
%
%  First do an RREF on A alone.
%  The second argument is the independent columns found in the input matrix.
%
  [ ARREF, a_cols ] = rref_compute ( A );

  if ( length ( a_cols ) < n )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'rref_inverse(): Warning!\n' );
    fprintf ( 1, '  The input matrix seems to be singular.\n' );
    fprintf ( 1, '  The inverse could not be computed.\n' );
    error ( '  The input matrix seems to be singular.\n' );
  end
%
%  If A has N independent columns, append the identity and repeat the RREF.
%
  AI = [ A, eye(m,n) ];

  [ AIRREF, a_cols ] = rref_compute ( AI );
%
%  The last N columns of AIRREF are the inverse.
%
  AINV = AIRREF(1:m,n+1:2*n);

  return
end

