Fast algorithm for LU factorization If A is a symmetric matrix, is there a fast algorithm for LU factorization? I know this algorithm for non-symmetric matrix.
    For k = 1,..,n
       For i = k + 1,...,n
         mult := a_{ik}/a_{kk}
         a_{jk} := mult
         For j = k+1,...,n
             a_{ij} := a_{ij} - mult * a_{kj}
         Endfor
     Endfor
   Endfor

 A: See Wikipedia article, try using Cholesky decomposition. 
Added Later: The matrix has to be positive definite. So this is just a partial answer to the question.
A: 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 
