#! /usr/bin/env python3
#
def rref_compute_test ( ):

#*****************************************************************************80
#
## rref_compute_test() tests rref_compute().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    05 December 2024
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  print ( '' )
  print ( 'rref_compute_test():' )
  print ( '  rref_compute() computes the reduced row echelon form (RREF)' )
  print ( '  of a matrix.' )

  A = np.array ( [ \
    [  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 ] ] )

  print ( '' )
  print ( '  Matrix A:' )
  print ( A )

  A_RREF = rref_compute ( A )

  print ( '' )
  print ( '  rref_compute(A):' )
  print ( A_RREF )

  return

def 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:
#
#    Original Python version by ChatGPT.
#    This version by 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 B(M,N), the reduced row echelon form of the matrix.
#
#    integer a_cols(*): the columns for which a pivot entry was found.
#
  import numpy as np
#
#  Ensure the matrix elements are floating-point for division.
#
  A = A.astype ( float )
#
#  Get the array shape.
#
  rows, cols = A.shape
#
#  Set a tolerance for the size of a pivot value.
#
  tol = np.sqrt ( np.finfo(float).eps )
#
#  Start from the first row.
#
  row = 0
#
#  Seek a pivot column COL.
#
  for col in range ( cols ):
#
#  Exit if we have run out of rows to examine.
#
    if ( rows <= row ):
      break
#        
#  Find the pivot row.
#
    p = np.argmax ( np.abs ( A[row:rows, col] ) ) + row
#
#  Skip this column if all values are below the tolerance.
#
    if ( np.abs ( A[p, col] ) <= tol ):
      continue
#
#  ERO: Swap current row and pivot row.
#
    A[[row, p]] = A[[p, row]]
#  
#  ERO: Scale the pivot row to make the pivot element equal to 1.
#
    A[row] = A[row] / A[row, col]
#   
#  ERO: Add multiples of pivot row to eliminate columns above and below.
#
    for i in range ( rows ):
      if ( i != row ):
        A[i] = A[i] - A[i, col] * A[row]
#
#  Move to the next row.
#
    row = row + 1
    
  return A

if ( __name__ == "__main__" ):
  rref_compute_test ( )