Like others above, I recommend cholesky. I've found that the increased number of additions, subtractions and memory accesses means that LDLt is slower than cholesky.
There are in fact a number a number of variations on cholesky, and which one will be fastest depends on the representation you choose for your matrices. I generally use a fortran style representation, that is a matrix M is a double* M with M(i,j) being m[i+dim*j]; for this I reckon that an upper triangular cholesky is (a little) the fastest, that is one seeks upper triangular U with U'*U = M.
For fixed, small, dimension it is definitely worth considering writing a version that uses no loops. A relatively straightforward way to do this is to write a program to do it. As I recall, using a routine that deals with the general case as a template, it only took a morning to write a program that would write a specific fixed dimension version. The savings can be considerable. For example my general version takes 0.47 seconds to do a million 9x9 factorisations, while the loopless version takes 0.17 seconds -- these timings running single threaded on a 2.6GHz pc.
To show that this is not a major task, I've included the source of such a program below. It includes the general version of the factorisation as a comment. I've used this code in circumstances where the matrices are not close to singular, and I reckon it works ok there; however it may well be too crude for more delicate work.
/* ----------------------------------------------------------------
** to write fixed dimension ut cholesky routines
** ----------------------------------------------------------------
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <strings.h>
/* ----------------------------------------------------------------
*/
#if 0
static inline double vec_dot_1_1( int dim, const double* x, const double* y)
{
double d = 0.0;
while( --dim >= 0)
{ d += *x++ * *y++;
}
return d;
}
/* ----------------------------------------------------------------
** ut cholesky: solve U'*U = P for ut U in P (only ut of P accessed)
** ----------------------------------------------------------------
*/
int mat_ut_cholesky( int dim, double* P)
{
int i, j;
double d;
double* Ucoli;
for( Ucoli=P, i=0; i<dim; ++i, Ucoli+=dim)
{ /* U[i,i] = P[i,i] - Sum{ k<i | U[k,i]*U[k,i]} */
d = Ucoli[i] - vec_dot_1_1( i, Ucoli, Ucoli);
if ( d < 0.0)
{ return 0;
}
Ucoli[i] = sqrt( d);
d = 1.0/Ucoli[i];
for( j=i+1; j<dim; ++j)
{ /* U[i,j] = (P[i,j] - Sum{ k<i | U[k,i]*U[k,j]})/U[i,i] */
P[i+j*dim] = d*(P[i+j*dim] - vec_dot_1_1( i, Ucoli, P+j*dim));
}
}
return 1;
}
/* ----------------------------------------------------------------
*/
#endif
/* ----------------------------------------------------------------
**
** ----------------------------------------------------------------
*/
static void write_ut_inner_step( int dim, int i, int off)
{
int j, k, l;
printf( "\td = 1.0/P[%d];\n", i+off);
for( j=i+1; j<dim; ++j)
{ k = i+j*dim;
printf( "\tP[%d] = d * ", k);
if ( i)
{ printf( "(P[%d]", k);
printf( " - (P[%d]*P[%d]", off, j*dim);
for( l=1; l<i; ++l)
{ printf( " + P[%d]*P[%d]", l+off, l+j*dim);
}
printf( "));");
}
else
{ printf( "P[%d];", k);
}
printf( "\n");
}
}
static void write_dot( int n, int off)
{
int i;
printf( "P[%d]*P[%d]", off, off);
for( i=1; i<n; ++i)
{ printf( "+P[%d]*P[%d]", off+i, off+i);
}
}
static void write_ut_outer_step( int dim, int i, int off)
{
printf( "\td = P[%d]", off+i);
if ( i)
{ printf( " - (");
write_dot( i, off);
printf( ")");
}
printf( ";\n");
printf( "\tif ( d <= 0.0)\n");
printf( "\t{\treturn 0;\n");
printf( "\t}\n");
printf( "\tP[%d] = sqrt( d);\n", i+off);
if ( i < dim-1)
{ write_ut_inner_step( dim, i, off);
}
}
static void write_ut_chol( int dim)
{
int i;
int off=0;
printf( "int\tut_chol_%.2d( double* P)\n", dim);
printf( "{\n");
printf( "double\td;\n");
for( i=0; i<dim; ++i)
{ write_ut_outer_step( dim, i, off);
printf( "\n");
off += dim;
}
printf( "\treturn 1;\n");
printf( "}\n");
}
/* ----------------------------------------------------------------
*/
/* ----------------------------------------------------------------
**
** ----------------------------------------------------------------
*/
static int read_args( int* dim, int argc, char** argv)
{
while( argc)
{ if ( strcmp( *argv, "-h") == 0)
{ return 0;
}
else if (strcmp( *argv, "-d") == 0)
{ --argc; ++argv;
*dim = atoi( (--argc, *argv++));
}
else
{ break;
}
}
return 1;
}
int main( int argc, char** argv)
{
int dim = 9;
if( read_args( &dim, --argc, ++argv))
{ write_ut_chol( dim);
}
else
{ fprintf( stderr, "usage: wchol (-d dim)? -- writes to stdout\n");
}
return EXIT_SUCCESS;
}
/* ----------------------------------------------------------------
*/
