inverse matrix in mexFunction

2012 4 4
4 个回答
回答已采纳
6 次查看（30 天）

0 个投票

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

请先登录，再回答此问题。

 采纳的回答

Jan
Jan 2012-4-4

0 个投票

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.

5 个评论
显示 3更早的评论 隐藏 3更早的评论

Jane Jean
Jane Jean 2012-4-4
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
Titus Edelhofer 2012-4-4
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
Jane Jean 2012-4-4
I shall try again with mxSignedIndex. But how do I declare it in mwSignedIndex?
Jane Jean
Jane Jean 2012-4-4
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
Jan 2012-4-4
mwSignedIndex is defined in tmwtypes.h, which is included by mex.h. So "#include mex.h" should be enough.

请先登录，再进行评论。

更多回答（3 个）

Titus Edelhofer
Titus Edelhofer 2012-4-4

0 个投票

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

1 个评论
显示 -1更早的评论 隐藏 -1更早的评论

Jane Jean
Jane Jean 2012-4-4
Thanks for the reply. I've pasted the code in my mexFunction below.

请先登录，再进行评论。

Jane Jean
Jane Jean 2012-4-4

0 个投票

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

1 个评论
显示 -1更早的评论 隐藏 -1更早的评论

Jane Jean
Jane Jean 2012-4-4
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?

请先登录，再进行评论。

Jane Jean
Jane Jean 2012-4-4

0 个投票

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;

类别

在 帮助中心 和 File Exchange 中查找有关 Write C Functions Callable from MATLAB (MEX Files) 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by