function L = cholesky_factor ( A )
%
%  function L = cholesky_factor ( A )
%
%  CHOLESKY_FACTOR computes the Cholesky factor of a matrix A.
%
%  Discussion:
%
%    The matrix A must be symmetric and positive semidefinite.
%
%    The Cholesky factorization is a lower triangular matrix L such that:
%
%      A = L * L'  
%
%  Modified:
%
%    21 March 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A(N,N), the matrix to be factored.
%
%    Output, real L(N,N), the Cholesky factor of A.
%
L = A;

[ n, n ] = size ( A );

for j = 1 : n

  L(1:j-1,j) = zeros ( j-1, 1 );

  for i = j : n

    sum = L(j,i) - L(j,1:j-1)' * L(i,1:j-1);

    if ( i == j ) 
      if ( sum <= 0.0 )
        fprintf ( '\n' );
        fprintf ( 'CHOLESKY_FACTOR - Fatal error!\n');
        fprintf ( '  Matrix is not positive definite.\n');
        return;
      else
        L(i,j) = sqrt ( sum );
      end
    else
      L(i,j) = sum / L(j,j);
    end

  end

end


