Classic Mx had a \chol function to obtain L for symmetric positive (semi-)definite A in A = L * L'. I figure that the BLAS function
dpbtrf could be used: SUBROUTINE DPBTRF( UPLO, N, KD, AB, LDAB, INFO )
and that it wouldn't be too hard to add this to the list of omxAlgebraFunctions.c - but that it would be harder for me than for say Tim Brick to do this. I made a start but have been hauled off to do other frantically urgent things. I suspect we need to tell the algebra cruncher routine the priority of this operation also.
void omxMatrixCholesky(omxMatrix** matList, int numArgs, omxMatrix* result)
{
if(OMX_DEBUG_ALGEBRA) { Rprintf("ALGEBRA: Matrix Cholesky Decomposition.\n");}
omxMatrix* inMat = matList[0];
int ipiv[inMat->rows], lwork;
lwork = 4 * inMat->rows * inMat->cols;
double work[lwork];
int l = 0;
if(rows != cols) {
char *errstr = Calloc(250, char);
sprintf(errstr, "Cholesky decomposition of non-square matrix cannot even be attempted\n");
omxRaiseError(result->currentState, -1, errstr);
Free(errstr);
}
omxCopyMatrix(result, inMat);
F77_CALL(dpbtrf)('L', &(result->cols), result->data, &(result->leading), ipiv, &l);
if(l != 0) {
char *errstr = Calloc(250, char);
sprintf(errstr, "Attempted to Choleskify non-invertable matrix.");
omxRaiseError(result->currentState, -1, errstr);
Free(errstr);
}
}