inverse matrix in mexFunction

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

Answer 1

[Jan]

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.

[Jane Jean]

Thanks for the reply. My compiler did not recognize mwSignedIndex so I used size_t instead, which works perfectly for dgemm matrix multiplication. However Matlab still crashed after a couple of iterations...

[Titus Edelhofer]

I think Jan's right: I changed int to mwSignedIndex everywhere (esp. of course the sizeof(int) when claiming memory. No failures then ...

mwSignedIndex should be declared in mex.h

Titus

[Jane Jean]

I shall try again with mxSignedIndex. But how do I declare it in mwSignedIndex?

[Jane Jean]

I found the problem now. The code that I pasted below did not work as I forgot to change the 'int' in sizeof to 'size_t'. Thank you guys for pointing out the errors!!! ;) You saved my day!

[Jan]

mwSignedIndex is defined in tmwtypes.h, which is included by mex.h. So "#include mex.h" should be enough.

Answer 2

[Titus Edelhofer]

Hi Jane,

the code looks fine to me. Maybe it's the calling mex file causing the problem?

Titus

[Jane Jean]

Thanks for the reply. I've pasted the code in my mexFunction below.

Answer 3

[Jane Jean]

  Sigma_Psi_inv_t = mxMalloc(rowSigma_Psi*colSigma_Psi*sizeof(double));
  if(!Sigma_Psi_inv_t){
    mexErrMsgTxt("Memory allocation error in least-squares initialization.");
  }
  mexPrintf("entering inverseMatrix");
  inverseMatrix(rowSigma_Psi, arraySigma_Psi, Sigma_Psi_inv_t);

  outputRowLen = rowSigma_Psi;
  outputColLen = colSigma_Psi;
  plhs[0] = mxCreateDoubleMatrix(outputRowLen, outputColLen, mxREAL); 

  arrayOutput = mxGetPr(plhs[0]);
  for (i=0; i<outputRowLen; i++){
    for(j=0; j<outputColLen; j++){
      arrayOutput[i + j*outputRowLen] = Sigma_Psi_inv_t[i*outputColLen +j];
  
    }
  }
  
  mxFree(Sigma_Psi_inv_t);

[Jane Jean]

I still can't find the error of the code. I have tried mexPrintf but after a few iterations, Matlab just crashed without printing out any error message.

Is there anyway to efficiently debug to error?

Answer 4

[Jane Jean]

void inverseMatrix(size_t dim, double *matrix, double *invMatrix)
{
  // matrix and invMatrix are columnwise.
  int *IPIV, LWORK, INFO=0, i;
  double *WORK;
  mexPrintf("entered inverseMatrix");
  IPIV = (int*)mxMalloc((dim+1)*sizeof(int));
  LWORK = dim*dim;
  WORK = (double*)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;
}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{

  mxArray *xData;
  double *xValues, *outArray, *invMatrix;
  int i,j;
  int rowLen, colLen;
  size_t row;

  xData = prhs[0];

  xValues = mxGetPr(xData);
  rowLen = mxGetM(xData);
  colLen = mxGetN(xData);
  row = rowLen;

  plhs[0] = mxCreateDoubleMatrix(rowLen, colLen, mxREAL); 

  outArray = mxGetPr(plhs[0]);
  inverseMatrix(row, xValues, mxGetPr(plhs[0]));
  
  return;
}

to compile and run.

      for i=1:50
      blaslib = fullfile('C:\Program Files\MATLAB\R2011b\extern\lib\win64\microsoft\libmwblas.lib');
          lapacklib = fullfile('C:\Program Files\MATLAB\R2011b\extern\lib\win64\microsoft\libmwlapack.lib');
          mex('-largeArrayDims', 'test_id.c','mathFunction64.c',  blaslib, lapacklib)
      A=[0,1,2,3;4,5,6,7;8,9,10,11;12,13,14,15]
          [N] = test_id(A);
      end;