[Numpy-discussion] rewriting NumPy code in C or C++ or similar
Sturla Molden
sturla@molden...
Tue Mar 8 08:25:12 CST 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.
Sturla
More information about the NumPy-Discussion
mailing list