With mxMalloc() and mxSetDoubles(), mex function crash matlab

I try to build a mex interface for a C function r8mat_expm1() that calculates matrix exponential,

double *r8mat_expm1 ( int n, double a[] );

Where n is the dimension and a is the matrix.

My mex function my_expm.c is as follows,

void mexFunction(int nlhs,mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    double *tmpResult;
	double *inputMat;
    int numRows;
    
    numRows = mxGetM(prhs[0]);
    inputMat = mxGetPr(prhs[0]); 
    
    tmpResult = mxMalloc(numRows*numRows*sizeof(*tmpResult));
	tmpResult = r8mat_expm1(numRows, inputMat); // calculate matrix exponential

    plhs[0] = mxCreateDoubleMatrix(numRows, numRows, mxREAL);
    mxSetDoubles(plhs[0], tmpResult);
    
    mxFree(tmpResult);
}

Which parse the input and call the computational routine r8mat_expm1().

After building on linux, I use

c=my_expm(2,[1,0;0,2]') % this is problematic, thanks Bruno Luong for pointing out

to run the compiled function. However, it just crashes matlab, and "Segmentation fault (core dumped) matlab" appears on linux terminal.

I suspect this is related to mxMalloc() and mxSetDoubles(). When I did not use these two functions, matlab will not crash. I chose to use mxSetDoubles() to fix an error, as suggested in this post.

I also tried another suggestion in that post, i.e. change r8mat_expm1() from

double *r8mat_expm1 ( int n, double a[] );

to

void r8mat_expm1 ( int n, double a[], doule output[] );

to avoid return all zeros error, but it still return all zeros.

So what is the problem behind this and how can I fix this crash error or get the right ouput instead of all zeros?

[Steven Lord]

Is your goal to write a function to compute the matrix exponential or is your goal to compute the matrix exponential? If the latter, consider simply calling the expm function in MATLAB. [Even if you want to or are required to write the function yourself (for a homework assignment, perhaps) you can check your results against expm.]

[Bruno Luong]

He did compare his code againts expm.

[Xingwang Yong]

The latter, but I want to do this in C language. I compiled it to see if it works by comparing its output to matlab built-in expm(). But later I found there is already a test case to demonstrate its effectiveness. So no need to compile it for comparison...

Answer 1

[James Tursa]

You have memory leaks in your current code. E.g.,

    tmpResult = mxMalloc(numRows*numRows*sizeof(*tmpResult)); // this memory gets leaked by the very next line
    tmpResult = r8mat_expm1(numRows, inputMat); // leaks memory from the mxMalloc( ) call
    plhs[0] = mxCreateDoubleMatrix(numRows, numRows, mxREAL); // this memory gets leaked by the very next line
    mxSetDoubles(plhs[0], tmpResult); // leaks memory from the mxCreateDoubleMatrix( ) call

This leaked memory is bad but does not cause the crash.

It appears you attach native C/C++ allocated memory to an mxArray, and this will cause a crash. You can't mix native C/C++ allocated memory in an mxArray. All memory in an mxArray must come from MATLAB API memory allocation functions.

You really need to post some details about the r8mat_expm1( ) function for us to be sure about what the real problem is. But from your description it sounds like it allocates the memory for the result using native C/C++ routines (e.g., malloc or calloc) and then returns the pointer to that, which you then attach to an mxArray ... resulting in a crash. If this is true, then you will either need to change the function to use MATLAB API allocation functions, or copy the results. To copy the results you would do something like this:

#include <string.h> // for memcpy( )
#include <stdlib.h> // for free( )
    double *tmpResult;
    double *inputMat;
    double *pr;
    int numRows;
    numRows = mxGetM(prhs[0]);
    inputMat = mxGetDoubles(prhs[0]);
    tmpResult = r8mat_expm1(numRows, inputMat); // calculate matrix exponential
    plhs[0] = mxCreateDoubleMatrix(numRows, numRows, mxREAL);
    pr = mxGetDoubles(plhs[0]);
    memcpy(pr,tmpResult,numRows*numRows*sizeof(double));
    free(tmpResult); // use free( ) for native C/C++ allocated memory

To avoid the data copy and use MATLAB API allocation functions, you would need to alter the source code of r8mat_expm1( ) and pass in the pointer instead of allocating memory inside the function. Then your mexFunction( ) code would look like this instead of the above:

    double *inputMat;
    double *pr;
    int numRows;
    numRows = mxGetM(prhs[0]);
    inputMat = mxGetDoubles(prhs[0]);
    plhs[0] = mxCreateDoubleMatrix(numRows, numRows, mxREAL);
    pr = mxGetDoubles(plhs[0]);
    r8mat_expm1(numRows, inputMat, pr); // calculate matrix exponential

And the signature for the r8mat_expm1( ) function would change to this:

    void r8mat_expm1(int, double *, double *);

where the 3rd input argument is used inside r8mat_expm1( ) for the output, and no memory allocation is done inside this function.

Another way to do this is to change the malloc( ) and calloc( ) inside r8mat_expm1( ) to be mxMalloc( ) and mxCalloc( ) instead, and then attach that result to the mxArray. If you want I can post that code, but frankly it is more complicated than what I have posted above so I would see no practical reason to code it this way.

[Xingwang Yong]

I used the 1st solution, i.e. memcpy() , it worked. Thank you very much, James. Your expertise really helped me a lot.

The source code of r8mat_expm1() can be found here and the sub-routines can be found here.

As for the details of r8mat_expm1(), you are correct, saying " it allocates the memory for the result using native C/C++ routines". The function r8mat_expm1() does not use malloc() itself, but it called a sub-routine, i.e. r8mat_identity_new(), which use malloc() to allocate memory for the final result.

For the 2nd solution, there was a time I modified only the signature of r8mat_expm1() as you suggested,

// the original code
double *r8mat_expm1 ( int n, double a[] )
{	 
  double *e;
  e = r8mat_identity_new(); 
  return e;
}

// modified code
void r8mat_expm1 ( int n, double a[], double e[] )
{	 
  //double *e;
  e = r8mat_identity_new(); 
  return;
}

In the mex function, I used

double *outMat;
r8mat_expm1(numRows, inputMat, outMat)

but the outMat are all zeros even when I debug in gdb. I can understand I would get wrong result in matlab because I mixed native C memory and matlab API memory. But I can not understand why outMat are all zeros even in native C memory. All thing I've done is changing the signature. Is this related to malloc() inside the sub-routine r8mat_identity_new() ?

For the 3rd solution, as you said, it is more complicated since r8mat_expm1() and its sub-routines extensively use malloc().

Surprisingly, the compiled my_expm() is 15x times slower than matlab built-in expm().

Again, thanks for your help!

[Bruno Luong]

You cannot do this

// modified code
void r8mat_expm1 ( int n, double a[], double e[] )
{	 
  //double *e;
  e = r8mat_identity_new(); 
  return;
}

The statement

e = r8mat_identity_new();

changes the pointer e to the newly allocated array by r8mat_identity_new. Your input pointer e is lost (and the INPUT array is not filled).

You should remove the MALLOC allocation of r8mat_identity_new pass your e as input and populate it with 1/0.

This is quite similar to the mistake you did with

tmpResult = mxMalloc(...);
tmpResult = r8mat_expm1(...);

Until you fully understand how really pointer works, you'll keep doing such mistakes.

[Xingwang Yong]

Why would my input pointer e lost? I thought

e = r8mat_identity_new();

is an assign statement that assgins the newly allocated array to my input pointer e. I can understand you saying "This is quite similar to the mistake you did with", because tmpResult is allocated twice, which caused memory leak. However, e is allocated only once here.

In the original code, the statement does the assignment job, the only difference is that e is declared inside the function.

It seems this assign stament behaves differently for (1) e is input parament , and (2) e is declared inside.

Did I miss anything here? Or is there any tutorial about this pointer-thing?

Many thanks!

[Bruno Luong]

After

e = r8mat_identity_new()

you erases the pointer declared as input

void r8mat_expm1 ( int n, double a[], double e[] )

This is the one you want to fill, since it was the data pointer of resulting mxArray lhs.

The statements

e = r8mat_identity_new()

allocates another pointer (with MALLOC inside the body) affects to "e" therefore the input pointer value is just lost. The new pointer has nothing to do with the input pointer value declared as INPUT argument of "double e[]."

What I explains here is nothing more than the basic (but fundamental) of C programming.

Each people learns C-programming (and pointers) according to they level, practice and speed. I can't recommend anything as specific tutorial for anyone.

[Xingwang Yong]

I think I fully understand the problem here. Thanks for your continuous help, Bruno!

[Xingwang Yong]

Finally, I tried the 3rd solution suggested by James, i.e. replace all the malloc with mxMalloc. It also worked and it is as fast as matlab built-in expm().

In the very first comment, I said "the compiled my_expm() is 15x times slower than matlab built-in expm().", which is not accurate because I forgot to remove the '-g' option when calling mex(). After removing '-g', the my_expm() modified using the memcpy() solution is 6~7 times slower than matlab built-in expm(). This is reasonable. I think copy will take some time.

[James Tursa]

I guess for completeness, here is how I would code the third method:

    double *tmpResult;
    double *inputMat;
    int numRows;
    numRows = mxGetM(prhs[0]);
    inputMat = mxGetDoubles(prhs[0]);
    tmpResult = r8mat_expm1(numRows, inputMat); // Changed to use mxMalloc and mxCalloc!!!
    plhs[0] = mxCreateDoubleMatrix(0, 0, mxREAL); // Data pointer starts out as NULL
    mxSetDoubles(plhs[0], tmpResult); // Set the data pointer
    mxSetM(plhs[0], numRows); // Set the row size
    mxSetN(plhs[0], numRows); // Set the column size

plhs[0] first gets allocated with 0x0 size. Data pointer is NULL and there is no data memory allocated. Then the pointer returned by r8mat_expm1( ) gets attached directly to plhs[0] ... but this is only possible IF the r8mat_expm1( ) function has been altered to use only MATLAB API memory allocation functions mxMalloc( ) and mxCalloc( ). Then set the size of plhs[0]. This method of coding avoids unnecessary memory allocation and memory leaks, but there are a few extra steps to make sure this happens properly. Note that there is no mxFree(tmpResult) listed above. This is not included in the code because the tmpResult pointer is now attached directly to plhs[0].

[Xingwang Yong]

I did exactly as you said above except that I did not set the size of plhs[0], i.e. I did not use

mxSetM(plhs[0], numRows); // Set the row size
mxSetN(plhs[0], numRows); // Set the column size

Is this used to speed up? I tried to append them to the last, but there is no perceivable speedup. Do I have to add these extra two lines?

[James Tursa]

I would suggest that you post the complete current code you are using for your mexFunction. Then we can comment on it for any potential problems or inefficiencies.

[Xingwang Yong]

I am using matrix exponential to solve differential equations,

where M is 5-by-1 vectors, A is a 5-by-5 matrix and A is changing with time. During a small time , A can be assumed to be constant, so the solution is

The complete code is too long, I'll post the most relative part,

M = M0;
for k = 1:length(time)
    A = calculate_A(time(k));
    M = my_expm(A*delta_t)*M;
end

So is it necessary to append mxSetM() to the last of mexFunction?

[James Tursa]

I meant the complete C code for mexFunction ...

[Xingwang Yong]

The code is similar to your version, except that I did not use mxSetM() at the last.

void mexFunction(int nlhs,mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    double *outputMat;
	double *inputMat;
    int numRows;
    
    numRows = mxGetM(prhs[0]);
    inputMat = mxGetPr(prhs[0]); 
    
    outputMat = mxMalloc(numRows*numRows*sizeof(*outputMat));
	outputMat = r8mat_expm1(numRows, inputMat); // calculate matrix exponential
    
    plhs[0] = mxCreateDoubleMatrix(numRows, numRows, mxREAL);
    mxSetDoubles(plhs[0], outputMat);
}

[James Tursa]

Assuming you have changed r8mat_expm1( ) to use mxMalloc or mxCalloc, this will work. However, this code has memory leaks as follows:

outputMat = mxMalloc(numRows*numRows*sizeof(*outputMat)); // (1)

Line (1) allocates memory and stores the pointer to it in outputMat. All fine and good.

outputMat = r8mat_expm1(numRows, inputMat); // calculate matrix exponential (2)

Line (2) allocates more memory and stores the pointer for this additional memory in the same variable outputMat as Line (1). This wipes out the pointer from Line (1) and you have now lost your only handle to the memory allocated from Line (1). You can now never free this memory. It is leaked.

plhs[0] = mxCreateDoubleMatrix(numRows, numRows, mxREAL); // (3)

Line (3) allocates memory for the plhs[0] mxArray, and the data pointer for this memory is stored inside the mxArray structure. All fine and good.

mxSetDoubles(plhs[0], outputMat); // (4)

Line (4) wipes out the data pointer from Line (3) and you have now lost your only handle to the data memory allocated from Line (3). You can now never free this memory. It is leaked.

For the Line (1) memory that you leaked, fortunately for you the MATLAB memory manager has a garbage collection facility when the mex function returns to the caller and this will get cleaned up for you automatically. But it is still bad programming practice to rely on this.

For the Line (3) data memory that you leaked, there is no safety net for you. This data pointer is not on the garbage collection list and will not get cleaned up by the MATLAB memory manager. It is leaked forever and can only be recovered by restarting MATLAB.

If you insist on this strategy rather than the code I posted, the proper way to do things to avoid the memory leaks would be:

//    outputMat = mxMalloc(numRows*numRows*sizeof(*outputMat));  DELETE THIS LINE COMPLETELY
	outputMat = r8mat_expm1(numRows, inputMat); // calculate matrix exponential
    
    plhs[0] = mxCreateDoubleMatrix(numRows, numRows, mxREAL);
    mxFree(mxGetDoubles(plhs[0])); // free the data memory currently in plhs[0]
    mxSetDoubles(plhs[0], outputMat); // set the data pointer in plhs[0]

This does extra unnecessary work than the method I posted, but at least it avoids memory leaks.

[Xingwang Yong]

Thanks, James, I got your two points.

outputMat = mxMalloc(numRows*numRows*sizeof(*outputMat)); // (1)
outputMat = r8mat_expm1(numRows, inputMat); // calculate matrix exponential (2)

This is totally wrong, as both you and Bruno has pointed out.

My misunderstading about pointers leads to another problem, in the original version, I wrote

 plhs[0] = mxCreateDoubleMatrix(numRows, numRows, mxREAL);
 mxSetDoubles(plhs[0], outputMat); // set the data pointer in plhs[0]

I thought this would "copy" the contents that outputMat points to into contents that plh[0] point to. But my thought is totally wrong. Actually, mxSetDoubles(plhs[0], outputMat) only set the pointer plhs[0] equal to the pointer outputMat.

By the way,

 mxFree(mxGetDoubles(plhs[0])); // free the data memory currently in plhs[0]

maybe this mxFree() should be replaced by mxDestroyArray(), because the doc of mxFree says "Do not use mxFree for an mxArray created by any other functions in the Matrix Library API. Use mxDestroyArray instead." I tried to do so, but the result is wrong.

Maybe an mxArray contains a lot of information, mxDestroyArray() destorys everything, including memory, array diemension, etc. But mxFree() only free the memory but leave the other information unchanged. Am I getting this right? I think what I say can explain why we use

plhs[0] = mxCreateDoubleMatrix(numRows, numRows, mxREAL);
mxFree(mxGetDoubles(plhs[0])); // free the data memory currently in plhs[0]
mxSetDoubles(plhs[0], outputMat); // set the data pointer in plhs[0]

i.e. first, create an array with all the information and then only free the memory.

[James Tursa]

mxDestroyArray( ) is only to be used on mxArray variables. I.e., those created with functions like mxCreateDoubleMatrix( ), mxCreateNumericArray( ), etc.

mxFree( ) is only to be used on "stand-alone" memory. I.e., allocated with functions such as mxMalloc( ), mxCalloc( ), and memory behind data pointers such as that returned with mxGetDoubles( ), mxGetComplexDoubles( ), etc.

mxFree(mxGetDoubles(plhs[0])) does not call mxFree( ) for an mxArray. Rather, it calls mxFree( ) for the data pointer of an mxArray. This is OK as long as you replace the data pointer in the very next line with something valid, which we did.

What would not be correct is the following:

mxFree(plhs[0]); // <-- WRONG!

[Xingwang Yong]

Yes, you are right.

I used

plhs[0] = mxCreateDoubleMatrix(numRows, numRows, mxREAL);
mxDestroyArray(plhs[0]);

This just deletes everything and I can not get the desired output.

[Bruno Luong]

Wait, why do you want to destroy/free with plhr[0] in your mex function? This supposes to be the result returned by the mex function, so you must leave it alone. What are trying to do?

I'm just afraid that you missundestand completely the dataflow of mex file.

[Xingwang Yong]

No, no, no, I just try it. I won't use it.

Answer 2

[Bruno Luong]

You must not call

mxFree(tmpResult);

Since the data is attached to the output.

[Xingwang Yong]

I deleted

mxFree(tmpResult);

and matlab still crashes with "Segmentation fault (core dumped) matlab"

[Bruno Luong]

Then you have other problem in the code r8mat_expm1.

[Xingwang Yong]

Yes, r8mat_expm1() called some other C sub-routines, which contains alloc() and free().

Do you have any idea how can I copy the output of 8mat_expm1(), which is double * , to the output of mex function?

[Xingwang Yong]

Yes, r8mat_expm1() called some other C sub-routines, which contains alloc() and free().

Do you have any idea how can I copy the output of 8mat_expm1(), which is double * , to the output of mex function?

Answer 3

[Bruno Luong]

You seem to call

mxGetM(prhs[0])

To get dimensions of the first parameter of calling

c=my_expm(2,[1,0;0,2]')

which is 2 and not the 2x2 matrix.

Who knows what else. The mex code seem to have quality issue, which is not tolerable and resulting crash.

[Xingwang Yong]

OH! Thanks, that is the problem!

Forgive me. In the 1st version, my mex function takes 2 inputs, and in the version I posted here, it takes only 1 input.

I forgot to change the command in matlab from

c=my_expm(2,[1,0;0,2]')

to

c=my_expm([1,0;0,2]')

Thank you very much! Everything works now!

[Xingwang Yong]

When I run the command for the first time,

c=my_expm([1,0;0,2])

everything is fine, however, when I run it again, it will crash matlab and without the previous "Segmentation fault (core dumped) matlab" message. It just says "killed" in linux terminal.

This is peculiar. Do you have any ideas?

[Bruno Luong]

i have no idea since you don't post anything related to your function

    tmpResult = mxMalloc(numRows*numRows*sizeof(*tmpResult));
	tmpResult = r8mat_expm1(numRows, inputMat)

The fact that you allocate the pointer, don't use it as input of the function, then return another pointer telling me that you don't know what are pointers and how to manipulate them.

As I said the code has ùany flaws (I just detect free of them with the small snip you have posted).

[Xingwang Yong]

Sorry for the inconvience. The full C code contains many lines and several sub-routines, I thought it would be concise if I posted a snippet only. The source code of r8mat_expm1() can be found here and its sub-routines can be found here.

I seldom use C, so actually I know little about pointers.

I still can not see why the code below does not work,

tmpResult = mxMalloc(numRows*numRows*sizeof(*tmpResult));
	tmpResult = r8mat_expm1(numRows, inputMat)

The signature of r8mat_expm1() is double *r8mat_expm1 ( int n, double a[] );

Why can't I assign the output of r8mat_expm1() to a pointer tmpResult? Do I have to pass tmpResult as input?

By the way,

double *tmpResult;
tmpResult = r8mat_expm1(numRows, inputMat);

if I do not use mxMalloc() for tmpResutl, this would give right result to tmpResult in mex fucntion although the result returned to matlab is wrong.

[Bruno Luong]

I think James has answered many question you raised already.

First, r8mat_expm1 use MALLOC to allocate resulting array "e", see the code of

double *r8mat_identity_new ( int n )
...
    a = ( double * ) malloc ( n * n * sizeof ( double ) );
...
    return a;

Your claim "The function r8mat_expm1() does not use malloc() for the result" is just plain wrong.

You are not allowed to assign a data array of an mxArray allocated with generic C MALLOC. You have to use mxMalloc.

The safest modification is to do exactly what he advises

#include <string.h> // for memcpy( )
#include <stdlib.h> // for free( )
    double *tmpResult;
    double *inputMat;
    double *pr;
    int numRows;
    numRows = mxGetM(prhs[0]);
    inputMat = mxGetDoubles(prhs[0]);
    tmpResult = r8mat_expm1(numRows, inputMat); // calculate matrix exponential
    plhs[0] = mxCreateDoubleMatrix(numRows, numRows, mxREAL);
    pr = mxGetDoubles(plhs[0]);
    memcpy(pr,tmpResult,numRows*numRows*sizeof(double));
    free(tmpResult); // use free( ) for native C/C++ allocated memory

Obviously you are not the owner of r8mat_expm1, If I was you I would avoid to change r8mat_expm1 in case the authors of those function updates the code in the futur.

"if I do not use mxMalloc() for tmpResutl, this would give right result to tmpResult in mex fucntion although the result returned to matlab is wrong."

The fact that you use mxMalloc before calling r8mat_expm1 is just make memory leaks (you allocate but then throw it away since the original pointer is lost after function returns). And beside it is still not allowed since r8mat_expm1 uses MALLOC.

[Xingwang Yong]

Thank you, Bruno. As you wrote this comment, I did my homework about pointers and digged up in the C source code and found the memory of final result is allocated by r8mat_identity_new(), same as you found.

"Your claim "The function r8mat_expm1() does not use malloc() for the result" is just plain wrong." Sorry for the mistake, at first I only searched the word 'malloc' inside r8mat_expm1() and did not go further into sub-routines. So I said in my previous comment "it did not call malloc() directly". I'll edit my comment under James's answer to avoid ambiguity.

"The fact that you use mxMalloc before calling mxGetDoubles", is this a typo error? I think mxGetDoubles should be replaced by r8mat_expm1 or mxSetDoubles. I did not use mxGetDoubles in my code.

Anyway, I've got your point about my wrong usage of pointers. Thank you.

[Bruno Luong]

Sorry, you are right I make a typo. Shoul read "The fact that you use mxMalloc before calling r8mat_expm1 "