[Numpy-discussion] rewriting NumPy code in C or C++ or similar

Sturla Molden sturla at molden.no
Tue Mar 8 09:25:12 EST 2011

Den 08.03.2011 05:05, skrev Dan Halbert:
> Thanks, that's a good suggestion. I have not written Fortran since 1971,
> but it's come a long way. I was a little worried about the row-major vs
> column-major issue, but perhaps that can be handled just by remembering
> to reverse the subscript order between C and Fortran.

In practice this is not a problem. Most numerical libraries for C assume 
Fortran-ordering, even OpenGL assumes Fortran-ordering. People program 
MEX files for Matlab in C all the time. Fortran-ordering is assumed in 
MEX files too.

In ANSI C, array bounds must be known at compile time, so a Fortran 
routine with the interface

     subroutine foobar( lda, A )
         integer lda
         double precision A(lda,*)
     end subroutine

will usually be written like

     void foobar( int lda, double A[]);

in C, ignoring different calling convention for lda.

Now we would index A(row,col) in Fortran and A[row + col*lda] in C. Is 
that too difficult to remember?

In ANSI C the issue actually only arises with small "array of arrays" 
having static shape, or convoluted contructs like "pointer to an array 
of pointers to arrays". Just avoid those and stay with 1D arrays in C -- 
do the 1D to 2D mapping in you mind.

In C99 arrays are allowed to have dynamic size, which mean we can do

    void foobar( int lda, double *pA )
       typedef double row_t [lda];
       vector_t *A = (vector_t*)((void*)&pA[0]);

Here we must index A[k][i] to match A(i,k) in Fortran. I still have not 
seen anyone use C99 like this, so I think it is merely theoretical.

Chances are if you know how to do this with C99, you also know how to 
get the subscripts right. If you are afraid to forget "to reverse the 
subscript order between C and Fortran", it just tells me you don't 
really know what you are doing when using C, and should probably use 
something else.

Why not Cython? It has "native support" for NumPy arrays.

Personally I prefer Fortran 95, but that is merely a matter of taste.


More information about the NumPy-Discussion mailing list