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:
%
%    28 September 2018
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real A(N,N), a square, invertible matrix.
%
%  Output:
%
%    real AINV(N,N), the inverse of A, as computed using MATLAB's
%    rref() function for reduced row echelon form.
%    If A is not square, or not invertible, then an empty matrix is returned.
%
  [ 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' );
    ainv = [];
    return;
  end
%
%  First do an RREF on A alone.
%  The second argument K is essentially the number of independent columns
%  found in the input matrix.
%
  [ arref, k ] = rref ( a );

  if ( length ( k ) < 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' );
    ainv = [];
    return
  end
%
%  If A has N independent columns, append the identity and repeat the RREF.
%
  ai = [ a, eye(m,n) ];

  [ airref, k ] = rref ( ai );
%
%  The last N columns of AIRREF are the inverse.
%
  ainv = airref(1:m,n+1:2*n);

  return
end

