[Numpy-discussion] Cython numerical syntax revisited

[David Cournapeau]

Sturla Molden wrote:
> The justification for Fortran ordering is in the mathematics. Say we have
> a set of linear equations
>
>    A * X = B
>
> and are going to solve for X, using some magical subroutine 'solve'. The
> most efficient way to store these arrays becomes the Fortran ordering.
> That is,
>
>    call solve(A, X, B)
>
> will be mathematically equivalent to the loop
>
>    do i = i, n
>       call solve(A, X(:,i), B)
>    end do
>
> All the arrays in the call to solve are still kept contigous! This would
> not be the case with C ordering, which is an important reason that C sucks
> so much for numerical computing. To write efficient linear algebra in C,
> we have to store matrices in memory transposed to how they appear in
> mathematical equations. In fact, Matlab uses Fortran ordering because of
> this.
>   

I don't think that's true: I am pretty sure matlab follows fortran
ordering because matlab was born as a "shell" around BLAS, LINPACK and co.

I don't understand your argument about Row vs Column matters: which one
is best depends on your linear algebra equations. You give an example
where Fortran is better, but I can find example where C order will be
more appropriate. Most of the time, for anything non trivial, which one
is best depends on the dimensions of the problem (Kalman filtering in
high dimensional spaces for example), because some parts of the
equations are better handled in a row-order fashion, and some other
parts in a column order fashion.

I don't know whether this is true, but I have read several times that
the column order in Fortran is historical and due to some  specificities
of the early IBM - but I have of course no idea what the hardware in
1954 looks like from a programming point of view :) This does not
prevent it from being a happy accident with regard to performances, though.

cheers,

David