MEX function doesn't initialise left-hand value

Question

I have to use a MEX function which is supposed to allocate a matrix and fill it (for example, make an identity matrix). Here is a code sample:

#include <stdlib.h>
#include "mex.h"

void sample(int n, double** T)
{
     (*T) = (double*)malloc(n*n*sizeof(double));  
     int i, j;
     if ((*T) == NULL)
     {
        return;
    }
    else
    {
        for (i = 0; i < n; i++)
        {
            for (j = 0; j < n; j++)
            {
                if (i != j)
                    (*T)[i*n + j] = 0;
                else 
                    (*T)[i*n + j] = 1;
             }
        }
    }
}

void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[] )
{
    mwSize n = *(mwSize*)mxGetData(prhs[0]);
    double *T;

    plhs[0] = mxCreateDoubleMatrix(1, (size_t)n*(size_t)n, mxREAL);

    T = mxGetPr(plhs[0]);

    sample(n, &T);  
}

I use it like this:

n = 5;
T = [];
T = sample(5, T);

This call returns a 1-by-0 empty matrix. So, it doesn't allocate a n-by-n matrix and doesn't fill it.

Answer 1

You have three very crucial errors and a few minor hiccups:

1. This line: mwSize n = *(mwSize*)mxGetData(prhs[0]); at the beginning of your MEX function.

   This makes absolutely no sense. We can clearly see that the input is a single integer. mxGetData actually gets a pointer to the actual data that is coming in as inputs into the MEX function, yet you are casting the pointer to mwSize then dereferencing it. I have absolutely no idea what this would give you, but it's definitely not right. What you can do instead is get the actual pointer to the data using mxGetPr where the pointer would provide you with data that conforms to a double precision real type. Because this would be an integer, you can cast this to int:

   int n = (int) *mxGetPr(prhs[0]);
   

2. Because of the correction with (1), you can now allocate your matrix without explicitly casting with size_t. This I wouldn't call an error, but it's more of a style correction:

   plhs[0] = mxCreateDoubleMatrix(1, n*n, mxREAL);
   

3. The most important error of all.

   Take note that your sample function inside your MEX code accepts a pointer to a double pointer in memory for the purposes of modifying the matrix that you created in your main MEX gateway (i.e. mexFunction). However, you first allocate the memory for the matrix to be written to the output, but are allocating the memory again inside your sample function. Because of this, you are actually creating another pointer to memory and writing the data to this pointer and not with the original pointer that you had passed to the function. Specifically, inside sample you have this code:

   (*T) = (double*)malloc(n*n*sizeof(double));
   

   Therefore, upon exit of this function the memory that you intended to modify has not been modified. The point is that allocation inside the function is not required as you have already allocated the memory of the matrix outside of the function. This can be removed from your code.

4. This is very important in terms of memory access. Remember that MATLAB is a column-major based language, which means that the columns of a matrix are laid out in a contiguous fashion instead of the rows like in languages such as C. Therefore, your method of accessing the matrix given that your matrix is of size m x n with m rows and n columns should be (*T)[j*m + i] where the index i accesses the rows and j accesses the columns. Right now you have it set to (*T)[i*n + j]. For now as your matrix shares the same rows and columns and that you are creating an identity matrix, either way of accessing is the same but be weary of rectangular sized matrices.

   Therefore, the corrected code for your sample function is thus:

   void sample(int n, double** T)
   {     
       int i, j;
   
       if (*T == NULL)
       {
            return;
       }
   
       for (i = 0; i < n; i++)
       {
           for (j = 0; j < n; j++)
           {
               if (i != j)
                   (*T)[j*m + i] = 0;
               else 
                   (*T)[j*m + i] = 1;
           }
       }
   }
   

   Once you run this code, it now produces the desired results. However, you have created a n*n row vector and not a n x n matrix, but I'm assuming you did that on purpose. To get this to run, I placed your source in a file called sample.c, compiled it then ran it:

   >> mex -O sample.c
   Building with 'gcc'.
   MEX completed successfully.
   >> n = 5;
   >> T = sample(n)
   
   T =
   
     Columns 1 through 19
   
    1     0     0     0     0     0     1     0     0     0     0     0     1     0     0     0     0     0     1
   
     Columns 20 through 25
   
    0     0     0     0     0     1
   

   If you wish to create a n x n matrix instead, change the call to mxCreateDoubleMatrix so that both the first and second arguments are set to n:

   plhs[0] = mxCreateDoubleMatrix(n, n, mxREAL);
   

5. This isn't an error but worth nitpicking about. Your function only expects 1 input so there's no need to submit T as the second parameter. Your code does not check for a second parameter anyway. In addition, you need to include the proper headers for your code to compile:

   #include <stdlib.h> /* For malloc */
   #include "mex.h" /* For the MEX library */
   

---

The final code that is corrected is thus:

#include <stdlib.h>
#include "mex.h"

void sample(int n, double** T)
{     
    int i, j;
    if (*T == NULL)
    {
        return;
    }

    for (i = 0; i < n; i++)
    {
        for (j = 0; j < n; j++)
        {
            if (i != j)
                (*T)[j*m + i] = 0;
            else 
                (*T)[j*m + i] = 1;
         }
    }
}

void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[] )
{
    int n = (int) *mxGetPr(prhs[0]);
    double *T;

    plhs[0] = mxCreateDoubleMatrix(1, n*n, mxREAL);
    /* Change to this line if you wish to have a n x n matrix */
    /* plhs[0] = mxCreateDoubleMatrix(n, n, mxREAL); */

    T = mxGetPr(plhs[0]);

    sample(n, &T);  
}

Just for the sake of argument, this is what happens if change mxCreateDoubleMatrix so that it outputs a n x n matrix:

>> mex -O sample.c
Building with 'gcc'.
MEX completed successfully.
>> n = 5;
>> T = sample(n)

T =

     1     0     0     0     0
     0     1     0     0     0
     0     0     1     0     0
     0     0     0     1     0
     0     0     0     0     1