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

#*****************************************************************************80
#
## rref_test() tests rref().
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    14 December 2024
#
#  Author:
#
#    John Burkardt
#
  print ( '' )
  print ( 'rref_test():' )
  print ( '  Test rref(), reduced row echelon format functions.' )

  is_rref_test ( )
  pivot_row_test ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'rref_test():' )
  print ( '  Normal end of execution.' )

  return

def is_rref_test ( ):

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

  print ( '' )
  print ( 'is_rref_test():' )
  print ( '  is_rref() reports whether a matrix is in reduced row echelon format.' )
#
#  Zero rows must come last.
#
  A0 = np.array ( [    \
    [ 1, 0, 0, 9, 4 ], \
    [ 0, 0, 1, 0, 8 ], \
    [ 0, 0, 0, 0, 0 ], \
    [ 0, 0, 0, 1, 0 ]  \
  ] )
#
#  First nonzero must be to right of first nonzero in previous row.
#
  A1 = np.array ( [    \
    [ 1, 0, 0, 9, 4 ], \
    [ 0, 0, 0, 1, 0 ], \
    [ 0, 0, 1, 0, 8 ], \
    [ 0, 0, 0, 0, 0 ] \
  ] )
#
#  First nonzero must be a 1.
#
  A2 = np.array ( [    \
    [ 1, 0, 0, 9, 4 ], \
    [ 0, 1, 0, 2, 8 ], \
    [ 0, 0, 3, 0, 0 ], \
    [ 0, 0, 0, 0, 0 ] \
  ] )
#
#  First nonzero must only nonzero in its column.
#
  A3 = np.array ( [    \
    [ 1, 0, 3, 9, 4 ], \
    [ 0, 1, 0, 2, 8 ], \
    [ 0, 0, 1, 0, 0 ], \
    [ 0, 0, 0, 0, 0 ] \
  ] )
#
#  RREF example.
#
  A4 = np.array ( [    \
    [ 1, 0, 3, 0, 4 ], \
    [ 0, 1, 2, 0, 8 ], \
    [ 0, 0, 0, 1, 0 ], \
    [ 0, 0, 0, 0, 0 ] \
  ] )

  for A, A_name in [ [ A0, 'A0' ], [ A1, 'A1' ], [ A2, 'A2' ], [ A3, 'A3' ], [ A4, 'A4' ] ]:
    print ( '' )
    print ( '  ' + A_name + ':' )
    print ( A )
    print ( '  is_rref(' + A_name + ') = ', is_rref ( A ) )

  return

def is_rref ( A ):

#*****************************************************************************80
#
## is_rref() determines if a matrix is in reduced row echelon format.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    14 December 2024
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    A[m,n]: a numpy array defining a matrix.
#
#  Output:
#
#    is_rref: True if A is in reduced row echelon form.
#
  m, n = A.shape

  c = -1

  for r in range ( 0, m ):
#
#  Increment the first legal column for next pivot.
#
    c = c + 1
#
#  Search for pivot p in this row.
#  If none, set p = n.
#
    p = n
    for j in range ( 0, n ):
      if ( A[r,j] != 0.0 ):
        p = j
        break
#
#  If p == n, go to next row.
#
    if ( p == n ):
      continue
#
#  If p is too early, fail.
#
    if ( p < c ):
      return False
#
#  Accept p as new c.
#
    c = p
#
#  If A(r,c) is not 1, fail
#
    if ( A[r,c] != 1.0 ):
      return False
#
#  If A(r,c) is not the only nonzero in column c, fail
#
    for i in range ( 0, m ):
      if ( i != r ):
        if ( A[i,c] != 0.0 ):
          return False

  return True

def pivot_row ( A, row, col ):

#*****************************************************************************80
#
## pivot_row() seeks the maximum matrix entry in A[row:,col].
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    05 December 2024
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    A[m,n]: the matrix to be analyzed.
#
#    row, col: search for the pivot in entries A[row:,col].
#
#  Output:
#
#    p: if p < m, then A[p,col] contains the pivot entry.
#    if p == m, then A[row:,col] is entirely zero, and no pivot could be found.
#
  import numpy as np

  m = A.shape[0]
#
#  Find the index of the entry of maximum magnitude in column COL, 
#  in or below row ROW.
#
  p = np.argmax ( np.abs ( A[row:,col] ) ) + row
#
#  If the maximum magnitude entry is zero, return the special value P = M.
#
  if ( A[p,col] == 0 ):
    p = m

  return p

def pivot_row_test ( ):

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

  print ( '' )
  print ( 'pivot_row_test():' )

  A = np.array ( [    \
    [  1, 10, 0, 0, 0 ], \
    [  0,  2, 4, 0, 0 ], \
    [ -5,  0, 8, 0, 0 ], \
    [  0,  9, 7, 0, 6 ]  \
  ] )

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

  m = A.shape[0]

  for row, col in [ [ 0, 0 ], [ 0, 1 ], [ 1, 2 ], [ 3, 2 ], [ 0, 3 ], [ 1, 4 ] ]:
    p = pivot_row ( A, row, col )

    print ( '' )
    if ( p == m ):
      print ( '  pivot_row( A,', row, ',', col, ') = ', p, '= NO PIVOT' )
    else:
      print ( '  pivot_row( A,', row, ',', col, ') = ', p )
      print ( '  A[', p, ',', row, '] = ', A[p,col] )

  return

def timestamp ( ):

#*****************************************************************************80
#
## timestamp() prints the date as a timestamp.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    21 August 2019
#
#  Author:
#
#    John Burkardt
#
  import time

  t = time.time ( )
  print ( time.ctime ( t ) )

  return

if ( __name__ == '__main__' ):
  timestamp ( )
  rref_test ( )
  timestamp ( )