FUNCTION chol,S
;+
;L=chol(S) returns the lower triangular matrix such that
;S=LL'. S must be positive definite.
;-

R=S
choldc,R,D
;Find the elements below the diagonal
dim=size(S,/dim)
k=indgen(dim[0],dim[1])
k1=where(k mod dim[0] LT k/dim[0])
L=diag(D)
L[k1]=R[k1]
Return,L
END