refmatrix <- function ( a )

#*****************************************************************************80
#
## refmatrix computes the Row Echelon Form (REF) of a matrix.
#
#  Licensing:
#
#    Copyright 2016 James P. Howard, II
#
#    The computer code and data files on this web page are distributed under
#    https://opensource.org/licenses/BSD-2-Clause, the BSD-2-Clause license.
#
#  Modified:
#
#    10 February 2020
#
#  Author:
#
#    Original R code by James Howard;
#    Modifications by John Burkardt.
#
#  Reference:
#
#    James Howard,
#    Computational Methods for Numerical Analysis with R,
#    CRC Press, 2017,
#    ISBN13: 978-1-4987-2363-3.
#
{
  source ( "/home/burkardt/public_html/r_src/eros/replacerow.R" )
  source ( "/home/burkardt/public_html/r_src/eros/swaprows.R" )

  count.rows <- nrow ( a )
  count.cols <- ncol ( a )
  piv <- 1

  for ( row.curr in 1:count.rows )
  {
    if ( piv <= count.cols )
    {
      i <- row.curr
      while ( a[i,piv] == 0.0 && i < count.rows )
      {
        i <- i + 1
        if ( i > count.rows )
        {
          i <- row.curr
          piv <- piv + 1
          if ( piv > count.cols )
          {
            return ( a )
          }
        }
      }
      if ( i != row.curr )
      {
        a <- swaprows ( a, i, row.curr )
      }
      for ( j in row.curr : count.rows )
      {
        if ( j != row.curr )
        {
          k <- a[j,piv] / a[row.curr,piv]
          a <- replacerow ( a, j, row.curr, - k )
        }
      }
      piv <- piv + 1
    }
  }
  return ( a )
}