Published on *OpenMx* (http://openmx.psyc.virginia.edu)

By *neale*

Created *03/20/2012 - 10:48*

Tue, 03/20/2012 - 10:48 — neale [1]

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);

}

}

**Links:**

[1] http://openmx.psyc.virginia.edu/users/neale

[2] http://openmx.psyc.virginia.edu/thread/2014

[3] http://openmx.psyc.virginia.edu/thread/1012

[4] http://openmx.psyc.virginia.edu/forums/openmx-developer-forums/new-feature-proposals