How to speed up (or avoid) data copy from C matrix to mxArray?

There is a function in my Matlab code that is being called many many times, and takes most of the computation time. I've rewritten that function in C, and it's been dramatically sped up, but it is still not enough for my purposes. Although my function is more complicated than this, let me give you a very simplified example (this is not my actual function, but a very simplified example):

#include <<stdio.h>>
#include <<math.h>>
#include "mex.h"
#include "matrix.h"
__inline void test_func( const double *__restrict__ im, double *__restrict__ im1, 
        register const unsigned int si, register const unsigned int sj, register const unsigned int sk ) {

      /* Initial test function. */
      register unsigned int ii ;                  /* Counters */
      for( ii=0; ii<si*sj*sk; ii++ ){
          *( im1 + ii ) = 2.* *( im + ii ) ;
      }
      return ;
  }

void mexFunction( const int nlhs, mxArray *plhs[], const int nrhs, const mxArray *prhs[] )
{
    /* Gateway routine used to interface the C code and the Matlab code. */
    /* Declare variables. */
    register const mwSize num_dims = mxGetNumberOfDimensions( prhs[0] ) ;
    register const mwSize *__restrict__ dims = mxGetDimensions( prhs[0] ) ;
    mxArray *__restrict__ im ;
    mxArray *__restrict__ im1 ;
    double im1_st[ dims[0]*dims[1]*dims[2] ] ;

      /* Check for proper number of input and output arguments. */
      if (nrhs != 1) {
          mexErrMsgTxt("One input argument required.");
      }
      if (nlhs != 1) {
          mexErrMsgTxt("One output argument required.");
      }
      /* Check data type of input argument. */
      if (!(mxIsDouble(prhs[0]))) {
          mexErrMsgTxt("Input array must be of type double.");
      }
      /* Get the data. */
      im = mxGetPr( prhs[0] );
      /* Allocate the space for the return argument */
      plhs[0] = mxCreateNumericArray( num_dims, dims, mxDOUBLE_CLASS, mxREAL) ;
      im1 = mxGetPr( plhs[0] ) ;
      /* Call the C subroutine. */
      test_func( im, im1, dims[0], dims[1], dims[2] ) ;
  }

If instead of calling test_func with the mxArray im1 I use the C array im1_st as

test_func( im, im1_st, dims[0], dims[1], dims[2] ) ;

the computation time gets dramatically reduced (takes only half time in this example, only one tenth in my real function). However, when I copy im1_st back to im1

memcpy( im1, im1_st, dims[0]*dims[1]*dims[2] * sizeof(double) ) ;

the computation time degrades back again. This makes me think that the way I am copying the C array into the Matlab mxArray is very inefficient.

My question is, is there a faster way to output the results of my C function back into Matlab? Is there a way I can replace the pointers of the mxArray im1 with the pointers of the C array im1_st (everything I've read seems to indicate this is not possible)?

I would really appreciate any idea on how to fix this issue.

Answer 1

[James Tursa]

No. You cannot mix C/C++ native memory with mxArray's. That will mess up the MATLAB memory manager and lead to an eventual crash for sure. You have to use MATLAB allocated memory (either from one of the mxCreateEtc functions or from mxMalloc and friends).

That being said, since you have the source code, why can't you just rewrite things to avoid the data copy? Clearly you could do so in your example. Why can't you do it in your real case? The im1_st in your example is a local variable (i.e., created from stack memory) and there is no way you can attach this to an mxArray. But as you have shown, simply getting the data pointer from an mxArray works just fine. So why not do this in your real code?

Also, you don't really show your two cases that you used for timing (one using im1 and the other using im1_st), so I am not sure you are comparing apples to apples. In your example above you create both arrays, which doesn't seem to apply. Can you show the two separate mex codes you used for timing?

Finally, a couple of comments on your checks:

      if (nrhs != 1) {
          mexErrMsgTxt("One input argument required.");
      }

It doesn't really do any good to make the above check where you have it in your code since you already accessed prhs[0] in prior code. You should move all of your prhs[0] stuff to below this check.

      if (nlhs != 1) {
          mexErrMsgTxt("One output argument required.");
      }

In general, the above check is too restrictive. It requires the user to assign the output of the mex routine to a variable. Usually it is OK to allow the user to simply let the output go into ans. So I would advise changing that check to:

      if (nlhs > 1) {
          mexErrMsgTxt("Too many outputs.");
      }

Even if nlhs == 0, MATLAB always allocates room in the output for plhs[0], so it is safe to assign an output to plhs[0] in this case (and the output will go into ans).

Answer 2

[Manu Gomez]

Hi James, thank you for your quick answer. Let me show you the simplified codes I used for timing.

The C function that does all the calculations on a 3D array is "test_func"

__inline void test_func( const double *__restrict__ im, double *__restrict__ im1, 
        register const unsigned int si, register const unsigned int sj, register const unsigned int sk ) {

      /* Initial test function. */

      register unsigned int ii ;                  /* Counters */

      for( ii=0; ii<si*sj*sk; ii++ ){

          *( im1 + ii ) = 2.* *( im + ii ) ;

      }

      return ;

}

Now, the original mex function

#include stdio.h

#include math.h

#include "mex.h"

#include "matrix.h"

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

    /* Gateway routine used to interface the C code and the Matlab code. */

    /* Declare variables. */
    register mwSize num_dims ;
    register mwSize *__restrict__ dims ;
    mxArray *__restrict__ im ;
    mxArray *__restrict__ im1 ;

    /* Check for proper number of input and output arguments. */
    if (nrhs != 1) {
        mexErrMsgTxt("One input argument required.");
    }

    if (nlhs > 1) {
        mexErrMsgTxt("One output argument required.");
    }

    /* Check data type of input argument. */
    if (!(mxIsDouble(prhs[0]))) {
        mexErrMsgTxt("Input array must be of type double.");
    }

    /* Get the number of elements in the input argument. */
    num_dims = mxGetNumberOfDimensions( prhs[0] );
    dims = mxGetDimensions( prhs[0] ) ;

    /* Get the data. */
    im = mxGetPr( prhs[0] );

    /* Allocate the space for the return argument */
    plhs[0] = mxCreateNumericArray( num_dims, dims, mxDOUBLE_CLASS, mxREAL) ;
    im1 = mxGetPr( plhs[0] ) ;

    /* Call the C subroutine. */
    test_func( im, im1, dims[0], dims[1], dims[2] ) ;

}

As you see, here I just make some checks on the input variables and call test_func. Everything is operated on mxArrays. I call this function with a 128x128x32 array, and I measure a time t0, that is pretty low, but I need to lower it further. For this reason, I try to operate with C arrays on the stack, not on the heap, and use the following mex function:

#include stdio.h

#include math.h

#include "mex.h"

#include "matrix.h"

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

    /* Gateway routine used to interface the C code and the Matlab code. */

    /* Declare variables. */
    register const mwSize num_dims = mxGetNumberOfDimensions( prhs[0] ) ;
    register const mwSize *__restrict__ dims = mxGetDimensions( prhs[0] ) ;
    double *__restrict__ im ;
    /* double *__restrict__ im1 ; */
    double im1_st[ dims[0]*dims[1]*dims[2] ] ;

    /* Check for proper number of input and output arguments. */
    if (nrhs != 1) {
        mexErrMsgTxt("One input argument required.");
    }

    if (nlhs != 1) {
        mexErrMsgTxt("One output argument required.");
    }

    /* Check data type of input argument. */
    if (!(mxIsDouble(prhs[0]))) {
        mexErrMsgTxt("Input array must be of type double.");
    }

    /* Get the data. */
    im = mxGetPr( prhs[0] );

    /* Allocate the space for the return argument */
    plhs[0] = mxCreateNumericArray( num_dims, dims, mxDOUBLE_CLASS, mxREAL) ;
    /* im1 = mxGetPr( plhs[0] ) ; */

    /* Call the C subroutine. */
    test_func( im, im1_st, dims[0], dims[1], dims[2] ) ;

}

Note that, here, I am operating on the array im1_st, and I am not copying the results back to Matlab space. The reason is that I wanted to know how much time is actually spent on the calculations. I find t1 that ranges from t0/10 to t0/2.

The third try is to copy the data back to Matlab space:

#include stdio.h

#include math.h

#include "mex.h"

#include "matrix.h"

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

    /* Gateway routine used to interface the C code and the Matlab code. */

    /* Declare variables. */
    register const mwSize num_dims = mxGetNumberOfDimensions( prhs[0] ) ;
    register const mwSize *__restrict__ dims = mxGetDimensions( prhs[0] ) ;
    double *__restrict__ im ;
    /* double *__restrict__ im1 ; */
    double im1_st[ dims[0]*dims[1]*dims[2] ] ;

    /* Check for proper number of input and output arguments. */
    if (nrhs != 1) {
        mexErrMsgTxt("One input argument required.");
    }

    if (nlhs != 1) {
        mexErrMsgTxt("One output argument required.");
    }

    /* Check data type of input argument. */
    if (!(mxIsDouble(prhs[0]))) {
        mexErrMsgTxt("Input array must be of type double.");
    }

    /* Get the data. */
    im = mxGetPr( prhs[0] );

    /* Allocate the space for the return argument */
    plhs[0] = mxCreateNumericArray( num_dims, dims, mxDOUBLE_CLASS, mxREAL) ;
    /* im1 = mxGetPr( plhs[0] ) ; */

    /* Call the C subroutine. */
    test_func( im, im1_st, dims[0], dims[1], dims[2] ) ;
    memcpy( mxGetPr( plhs[0] ), im1_st, dims[0]*dims[1]*dims[2] * sizeof(double) ) ;

}

This takes a time t2 very similar to t0. So I need to either operate on a mxArray, or to operate on a C array and copy the data back to the mxArray. Both approaches take a similar computation time, higher than what I need.

Given that I cannot just return a C array, is there any way I can speed the process up? Can Matlab's mxArrays be allocated in the Stack rather than the heap? Can the copy of the C array to the mxArray be sped up?

Thanks for your help.

[James Tursa]

Typo, this:

    mxArray *__restrict__ im ;
    mxArray *__restrict__ im1 ;

should be this:

    double *__restrict__ im ;
    double *__restrict__ im1 ;

"Can Matlab's mxArrays be allocated in the Stack rather than the heap?"

No, mxArrays must be allocated using official API functions, which use the heap (necessary in order to return the mxArray back to the caller).

"Can the copy of the C array to the mxArray be sped up?"

Copying data could be done in parallel, of course. But you would think that the memcpy library function, if written well, would have accounted for this already if the data amount was large enough to justify it. But for a 128x128x32 array, the overhead of multi-threading will not be worth it.

One thing you haven't apparently tried is simply doing the allocation uninitialized which is faster. Try this (caution: writing this on the fly so untested):

    #include <stdio.h>
    #include <math.h>
    #include "mex.h"
    #include "matrix.h"  /* Note, you don't really need this since mex.h includes matrix.h */

    void mexFunction( const int nlhs, mxArray *plhs[], const int nrhs, const mxArray *prhs[] )
    {
        /* Gateway routine used to interface the C code and the Matlab code. */
        /* Declare variables. */
        mwSize num_dims ;
        const mwSize *__restrict__ dims ;
        double *__restrict__ im ;
        double *__restrict__ im1 ;

        /* Check for proper number of input and output arguments. */
        if (nrhs != 1) {
            mexErrMsgTxt("One input argument required.");
        }
        if (nlhs > 1) {
            mexErrMsgTxt("Too many outputs.");
        }
        /* Check data type of input argument. */
        if (!(mxIsDouble(prhs[0])) || mxIsSparse(prhs[0])) {
            mexErrMsgTxt("Input array must be of type full double.");
        }
        /* Get the number of elements in the input argument. */
        num_dims = mxGetNumberOfDimensions( prhs[0] );
        dims = mxGetDimensions( prhs[0] ) ;

        /* Get the data. */
        im = mxGetPr( prhs[0] );

        /* Allocate the space for the return argument */
        plhs[0] = mxCreateUninitNumericArray( num_dims, dims, mxDOUBLE_CLASS, mxREAL) ;
        im1 = mxGetPr( plhs[0] ) ;

        /* Call the C subroutine. */
        test_func( im, im1, dims[0], dims[1], dims[2] ) ;
    }

Also maybe change this loop:

      register unsigned int ii ;                  /* Counters */
      for( ii=0; ii<si*sj*sk; ii++ ){
          *( im1 + ii ) = 2.* *( im + ii ) ;
      }

to something like this instead:

      register unsigned int ii = si*sj*sk ;         /* Counters */
      while( ii-- ){
          *im1++ = 2.* *im++ ;
      }

Note that mxCreateUninitNumericArray has been in the API library for a long time, but only added to the official API interface relatively recently. If you have an older version of MATLAB that does not have this in the official API interface you will have to manually include the prototype:

mxArray *mxCreateUninitNumericArray(size_t ndim, size_t *dims, 
  mxClassID classid, mxComplexity ComplexFlag);

Answer 3

[Manu Gomez]

Also, it is my understanding that, when calling a C function on an mxArray, the function is acting on a copy of the data, not on the actual data. Is there any way I can act on the actual data of the mxArray and avoid all the internal duplications that matlab is doing?

[James Tursa]

"... when calling a C function on an mxArray ..."

You need to be more specific here. Most API functions work on the actual mxArray, not a copy. E.g., mxGetPr always returns the address of the actual data pointer of the mxArray, not a pointer to a copy of the data. mxGetCell and mxGetField always return a pointer to the actual cell or field element, not a pointer to a copy. Etc.

API functions that return deep copies of stuff would be in the minority, e.g., mxGetProperty and mexGetVariable always get copies. And calling mexCallMATLAB will return a deep copy even if the code involved would normally return a shared data copy at the m-file level.