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

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

  is_ref_test ( )
  pivot_row_test ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'ref_test():' )
  print ( '  Normal end of execution.' )

  return

def is_ref_test ( ):

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

  print ( '' )
  print ( 'is_ref_test():' )
  print ( '  is_ref() reports whether a matrix is in row echelon format.' )

  A0 = np.array ( [    \
    [ 5, 3, 0, 9, 4 ], \
    [ 0, 0, 2, 7, 8 ], \
    [ 0, 0, 0, 6, 1 ], \
    [ 0, 0, 0, 0, 0 ]  \
  ] )

  A1 = np.array ( [    \
    [ 5, 3, 0, 9, 4 ], \
    [ 0, 0, 2, 7, 8 ], \
    [ 0, 0, 0, 0, 0 ], \
    [ 0, 0, 0, 6, 1 ]
  ] )

  A2 = np.array ( [    \
    [ 5, 3, 0, 9, 4 ], \
    [ 0, 0, 2, 7, 8 ], \
    [ 0, 0, 0, 6, 1 ], \
    [ 0, 0, 0, 4, 0 ]  \
  ] )

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

  return

def is_ref ( A ):

#*****************************************************************************80
#
## is_ref() determines if a matrix is in row echelon format.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    04 December 2024
#
#  Author:
#
#    John Burkardt
#
  m, n = A.shape

  pold = -1

  for r in range ( 0, m ):
#
#  Search for pivot p in this row.
#  If none, p = n.
#
    p = n
    for j in range ( 0, n ):
      if ( A[r,j] != 0.0 ):
        p = j
        break
#
#  If there was a nonzero pivot,
#  if must occur later than the previous one.
#
    if ( p < n ):
      if ( p <= pold ):
        return False

    pold = p

  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

  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 ( )
  ref_test ( )
  timestamp ( )