lumatrix <- function ( a )

#*****************************************************************************80
#
## lumatrix computes the Lower-Upper (LU) triangular factors 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:
#
#    14 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

  L <- diag ( count.cols )
  P <- diag ( count.cols )
  U <- A

  for ( row.curr in 1:count.rows )
  {
    if ( piv <= count.cols )
    {

      i <- row.curr

      while ( U[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 ( list ( P=P, L=L, U=U ) )
          }
        }
      }

      if ( i != row.curr )
      {
        U <- swaprows ( U, i, row.curr )
        P <- swaprows ( P, i, row.curr )
      }

      for ( j in row.curr : count.rows )
      {
        if ( j != row.curr )
        {
          k <- U[j,piv] / U[row.curr,piv]
          U <- replacerow ( U, j, row.curr, - k )
          L[j,piv] <- k
        }
      }
      piv <- piv + 1
    }
  }
  return ( list ( P=P, L=L, U=U ) )
}