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

[Ondrej Certik]

Hi Sturla,

On Tue, Mar 8, 2011 at 6:25 AM, Sturla Molden <sturla at molden.no> wrote:
> 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.

+1 to all that you wrote about Fortran. I am pretty much switching to
it from C/C++ for all my numerical work, and then I use Cython to call
it from Python, as well as cmake for the build system.

Ondrej