Return double pointer 2D matrix in mex function

Question

I am working with 2D arrays defined with double pointers, e.g,

double** array;  
array = (double**) calloc(numRows, sizeof(double*));  
for (int i = 0; i < numRows; i++)  
{
   array[i] = (double*) calloc(numCols, sizeof(double));  
   /* then array[i][j] = values */
}

// code to return matlab array
plhs[0] = mxCreateDoubleMatrix(numRows, numCols, mxREAL);

// memory copy
// ?????????

for (i = 0; i < numRows; i++){  
       free(array[i]);  
    }     
 free(array);

I want to return array in matlab. An execution I have until now for the // memory copy part and I think it is fine, please correct me is:

stacked1D = mxGetPr(plhs[0]);
int n = 0;
for ( int r = 0; r < max_degree; r++ )
    for ( int c = 0; c < n_vars; c++ )
        stacked1D[n++] = stacked2D[r][c];

I am wondered if we can do it with a mem-copy function like this mxSetPr(OUT, *stacked2D); which is not working in this syntax.

Could you please give a hint-explanation or possible answer?

Answer 1

Row and column iterations should be reversed in your code, and what PaulMcKenzie suggested, although it's a good idea in principle, will not work with Mex matrices (they are laid out column-by-column, so with that method you would have to access your matrix with M[column][row] which is un-natural and confusing).

Alternatively, you could use a simple wrapper like the following:

template <class T, class I = unsigned>
struct MatrixWrapper
{
    T *data;
    I nrows, ncols;

    MatrixWrapper() 
        { clear(); }
    MatrixWrapper( T *data_, I nrows_, I ncols_ )
        { set(data_,nrows_,ncols_); }

    inline void clear()
        { data = NULL; nrows = ncols = 0; }

    inline void set( T *data_, I nrows_, I ncols_ )
        { data = data_; nrows = nrows_; ncols = ncols_; }

    inline T& operator() ( I r, I c ) const
        { return data[ r + nrows*c ]; }
};

and your Mex function would look like this:

// allocate a temporary matrix
double *M_data = new double[ nrows*ncols ];
MatrixWrapper<double> M( M_data, nrows, ncols );

// access values with M(i,j) as a "normal" matrix

// allocate a Mex output
plhs[0] = mxCreateDoubleMatrix( nrows, ncols, mxREAL );
MatrixWrapper<double> out0( mxGetPr(plhs[0]), nrows, ncols );

// copy to the Mex output
for ( unsigned c = 0; c < ncols; c++ )
for ( unsigned r = 0; r < nrows; r++ )
    out0(r,c) = M(r,c);

// free temporary allocation
M.clear();
delete[] M_data;