choleskymatrix <- function ( a )

#*****************************************************************************80
#
## choleskymatrix computes the Cholesky L L' 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:
#
#    15 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.
#
{
  count.rows <- nrow ( a )
  count.cols <- ncol ( a )

  L <- diag ( 0, count.rows )

  for ( i in 1 : count.rows )
  {
    for ( k in 1 : i )
    {
      p.sum <- 0

      for ( j in 1 : k )
      {
        p.sum <- p.sum + L[j,i] * L[j,k]
      }

      if ( i == k )
      {
        L[k,i] <- sqrt ( A[i,i] - p.sum )
      }
      else
      {
        L[k,i] <- ( A[k,i] - p.sum ) / L[k,k]
      }
    }
  }
  return ( t(L) )
}