function [ A, a_cols ] = rref_compute ( A )

%*****************************************************************************80
%
%% rref2_compute() computes the reduced row echelon form of a matrix.
%
%  Discussion:
%
%    A rectangular matrix is in row reduced echelon form if:
%
%    * The leading nonzero entry in each row has the value 1.
%
%    * All entries are zero above and below the leading nonzero entry 
%      in each row.
%
%    * 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.
%
%  Example:
%
%    M = 4, N = 7
%
%    Matrix A:
%
%     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
%
%    RREF(A):
%
%     1    3    0    0    2    0    0
%     0    0    0    1    2    0    0
%     0    0    0    0    0    1   1/3
%     0    0    0    0    0    0    0
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    05 December 2024
%
%  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 A(M,N), the matrix to be analyzed. 
%
%  Output:
%
%    integer A(M,N), the reduced row echelon form of the matrix.
%
%    integer a_cols(*): the columns for which a pivot entry was found.
%
  [ m, n ] = size ( A );
%
%  Initialize rank to 0.
%
  a_cols = [];
%
%  Search for a pivot for each column, starting with column C = 1.
%
  c = 1;
%
%  Start the search for a pivot in row R = 1.
%
  r = 1;
%
%  Here, we use a tolerance TOL instead of checking for an exact zero.
%
  tol = sqrt ( eps );

  while ( true )
%
%  Find row index P of maximum element in subvector A(R:M,C).
%
    [ pval, p ] = max ( abs ( A(r:m,c) ) );
%
%  Adjust value of P as an entry in full vector A(1:M,C).
%
    p = p + r - 1;
%
%  If A(r:m,c) was all zero, there is no pivot for this column.
%
    if ( abs ( pval ) <= tol )

      A(p,c) = 0.0;
%
%  Increment the column index for the next search.
%
      c = c + 1;
%
%  If there are no more columns to search, terminate the algorithm.
%
      if ( n < c )
        break
      end
%
%  Restart loop with the next column index.
%
      continue;

    end
%
%  A pivot was found.  Increase rank by 1.
%
    a_cols = [ a_cols, c ];
%
%  Swap rows R and P.
%
    temp = A(r,c:n);
    A(r,c:n) = A(p,c:n);
    A(p,c:n) = temp;
%
%  Normalize row R so that A(r,c) = 1.
%
    A(r,c:n) = A(r,c:n) / A(r,c);
%
%  Eliminate nonzeros in column C using multiples of row R, except for A(R,C).
%
    for i = 1 : m
      if ( i ~= r )
        A(i,c:n) = A(i,c:n) - A(i,c) * A(r,c:n);
      end
    end
%
%  Increment C, but terminate if we reached the value C=N.
%
    if ( c < n )
      c = c + 1;
    else
      break;
    end
%
%  Increment R, but terminate if we reached the value R=M.
%
    if ( r < m )
      r = r + 1;
    else
      break;
    end

  end

  return
end

