double chol_check(double x)
/*
 * returns  sqrt(x) for x > 0,  0 for x <= 0
 */
{
    static string msg =
	(string)[7, "chol: warning, matrix may not be positive definite"];

    if( x <= 0.0) {
	print(msg);
	return 0.0;
    } else
	return sqrt(x);
}

matrix chol( matrix A)
/*
 * Cholesky decomposition of symmetric non-negative definite matrix A
 * returns lower triangular G such that  A = G * G'
 */
{
    int i, k, n = rows(A);
    matrix G = zeros(n,n);

    G(1,1) = chol_check( A(1,1));
    if( n > 1)
	G(2:n,1) = A(2:n,1) / G(1,1);

    for( k = 2; k <= n; ++k) {

	G(k,k) = chol_check( A(k,k) -  G(k,1:k-1) * G(k,1:k-1)' );
	if( G(k,k) != 0.0  &&  k < n)
	    G(k+1:n,k) = ( A(k+1:n,k) - G(k+1:n,1:k-1) * G(k,1:k-1)' ) / G(k,k);
    }

    return G;
}
