function [ A, a_cols ] = rref_compute ( A )

%*****************************************************************************80
%
%% rref_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:
%
%    27 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:
%
%    real A(M,N), the matrix to be analyzed. 
%
%  Output:
%
%    real A(M,N), the reduced row echelon form of the matrix.
%
%    integer a_cols(*): the columns for which a pivot entry was found.
%

%
%  Initialize an empty list of independent columns.
%
  a_cols = [];
%
%  Get the dimensions of the matrix.
%
  [ rows, cols ] = size ( A );
%
%  Set a tolerance for the size of a pivot value.
%
  tol = sqrt ( eps );
%
%  Start from the first row.
%
  row = 1;
%
%  Seek a pivot for each column.
%
  for col = 1 : cols
%
%  Exit if we have run out of rows to examine.
%
    if ( rows < row )
      break;
    end
%
%  Find row index P of maximum element in subvector A(R:M,C).
%
    [ pivot_value, pivot_row ] = max ( abs ( A(row:rows,col) ) );
%
%  Adjust value of P as an entry in full vector A(1:M,C).
%
    pivot_row = pivot_row + row - 1;
%
%  If A(row:rows,col) was essentially zero, there is no pivot for this column.
%
    if ( abs ( pivot_value ) <= tol )
      A(pivot_row,col) = 0.0;
      continue;
    end
%
%  A pivot was found for column col.
%
    a_cols = [ a_cols, col ];
%
%  Swap current row and pivot row
%
    temp = A(row,col:cols);
    A(row,col:cols) = A(pivot_row,col:cols);
    A(pivot_row,col:cols) = temp;
%
%  Normalize row R so that A(row,col) = 1.
%
    A(row,col:cols) = A(row,col:cols) / A(row,col);
%
%  Eliminate nonzeros in column C using multiples of row R, except for A(R,C).
%
    for i = 1 : rows
      if ( i ~= row )
        A(i,col:cols) = A(i,col:cols) - A(i,col) * A(row,col:cols);
      end
    end
%
%  Move to the next row.
%
    row = row + 1;

  end

  return
end

