MATLAB: Inverse matrix in mexFunction

dgetrfdgetriinversemex

void inverseMatrix(int dim, double *matrix, double *invMatrix)
{
  // matrix and invMatrix are columnwise.
  int *IPIV, LWORK, INFO, i;
  double *WORK;
  mexPrintf("entered inverseMatrix");
  IPIV = mxMalloc((dim+1)*sizeof(int));
  LWORK = dim*dim;
  WORK = mxMalloc(LWORK*sizeof(double));
  for (i=0;i<dim*dim;i++){
    invMatrix[i] = matrix[i];
  }
  mexPrintf("before dgetrf");
  dgetrf(&dim, &dim, invMatrix, &dim, IPIV, &INFO);
  mexPrintf("before dgetri");
  dgetri(&dim, invMatrix, &dim, IPIV, WORK, &LWORK, &INFO);
  mxFree(IPIV);
  mxFree(WORK);
  return;
}
I am trying to use dgetrf and dgetri to inverse a matrix in C but Matlab crashes after successfully giving the correct answer 2 times (I did an interation to try the stability of the c code).

Best Answer

The copy of the data to the output would be faster using memcpy:
plhs[0] = mxCreateDoubleMatrix(outputRowLen, outputColLen, mxREAL); 
memcpy(mxGetPr(plhs[0]), Sigma_Psi_inv_t, 
       outputRowLen*outputColLen * sizeof(double));
But you can use the data of the output array directly:
plhs[0] = mxCreateDoubleMatrix(rowSigma_Psi, colSigma_Psi, mxREAL); 
mexPrintf("entering inverseMatrix");
inverseMatrix(rowSigma_Psi, arraySigma_Psi, mxGetPr(plhs[0]));
Do not use int for LAPACK calls in an 64 bit environment. Use mwSignedIndex instead, which is a 64 bit integer.
Related Question