Row/Col major order in mex function

I have a question regarding row/column order in a C Mex-Function.

I am trying to access elements in a matrix contiguously, for speed.

I wrote a mex function in C, which takes a matrix as input:

im1 = mxGetPr(prhs[0]);

Let's say that

im1 = [1 2 3;  
       4 5 6]

Now, C has row major order, so linear indices are obtained from the row and column coordinates (r,c) as follows:

index = c + nbr_of_col * r;

However, when I print the values with

for( int r = 0; r < nbr_of_row; r++ ){
   for( int c = 0; c < nbr_of_col; c++){
      mexPrintf("im1[%d][%d] =  %0.2f \t", r, c, im1[c + nCol*r]);
   }
   mexPrintf("\n")
}

I get:

1   4   2
5   3   6

Which shows that it is traversing the matrix columnwise (using the linear indexing convention of Matlab): 1,4,2,5,3,6.

My question is : I can use the indexing formula

index = r + nbr_of_row*c

to traverse the matrix row-wise. Which major order is faster? I thought Row-major was faster since it's in C, but mxGetPr linearizes matrices according to the MATLAB convention...

[dpb]

mxGetPr might just produce a 1D array that is linearized according to Matlab convention (column order), and so the fact that C is row order is irrelevant.

mxGetPr doesn't "produce" anything -- it simply returns a pointer to the beginning of the first element of the data array no matter its apparent dimension in Matlab.

Since array storage is allocated contiguously, from a memory standpoint it is indeed as if it were a 1D array as that is what contiguous memory is.

The order of elements is, as you note, "column major" because the data were allocated in Matlab and that's the way it operates; if you were to have created an array in C and then looked at it the order would have been reversed.

Answer 1

[James Tursa]

Other things being equal, the fastest way is to access the data in the order that it is physically stored in memory. This makes the most efficient use of the processor cache system. Bottom line for MATLAB arrays, access them in column order even if you are inside a C mex routine. And I would advise to simply use your own linear indexing formulas for this (as you have shown above) rather than trying to coax a 2D array syntax such as variable[i][j] in your C code.

[Nicolas]

Great! Really appreciate the quick answer, that's exactly what I wanted to know. I definitely won't try to use a 2D array in the C code. Thank you!

Answer 2

[Mike Hosea]

I have written numerical library routines for the better part of the last 25 years, and it didn't take me very long to realize that I just could not use C's multidimensional arrays for anything because the dimensions needed to be compile-time constants. My advice is just to forget about how C stores matrices, and forget about vague ideas of about this or that being faster "because it is C". Allocating large matrices on the stack isn't a great idea, anyway. What C provides for matrices is showing you its historical roots as a mid-level language for writing an OS. If the focus had been anywhere near in the same vicinity as matrix calculations, all that was needed was to look at FORTRAN for an example of something useful .

If you have an application that in fact involves using smallish matrices with dimensions fixed at compile time, then maybe use C arrays if you want. When writing code to work with any matrix size, however, if you want to work in C using matrices in row-major order, then malloc the memory and instead of A[i][j], store and reference using A[i*ncols + j]. If you want to work in C using matrices in column-major order, store and reference using A[i + j*nrows]. In most cases there is effectively no inherent ordering in C. It is whatever you want it to be. If you are interfacing with MATLAB it is likely that column-major order will be preferred, since if you work in row-major order you will have to transpose on the way in and on the way out to switch from one ordering to the other.