The cholesky-decomposition can be made "robust" - just keep track of zeros and negative signs in the diagonal.
Here is a code-snippet in my (proprietary,sorry) MatMate-language, which should be easily translatable into a C- or Basic - or something-else routine:
TMP = X // X is the symmetric matrix
L = Null(X) // Null-matrix of same size as X;
// shall become the cholesky-factor
sg = L[*,1] // vector which keeps track of the signs in the diagonal
c = 0 // index for the current column/row
// repeat the following up to the number of rows/colums of X
c=c+1
d = w(TMP[c,c]) // the value of the diagonal element TMP[c,c]
// if zero or lower/equal required precision finish procedure
s = sign(d)
sg[c] = s // store the sign of the diagonal element
L[*,c] = TMP[*,c]/sqrt(s*d) // fill the c'th column
TMP = TMP - s * L[*,c] * L[*,c]' // reduce tmp;
// the apostroph means transposition
Check result:
chk = L * mkdiag(sg) * L' // check: re-compose cholesky-factors
// with signed-diagonal
err = X - chk // should be zero
We can even include, that the algorithm terminates when d becomes zero: the cholesky-factor L reflects then exactly the rank of X