row major order versus column major order

["E. Robert Tisdale" <edwin at netwood dot net>]

The only difference between multidimensional Fortran and C arrays
is the order in which the subscripts appear.
Fortran multidimensional array subscripts appear in order
from least to most significant and C multidimensional array subscripts
appear in reverse order from most to least significant so that

        A(i1, i2, ..., iN) == A[iN]...[i2][i1]

which means that multidimensional Fortran and C arrays
appear to be transposes of each other.  Both Fortran and C programmers
interpret one dimensional arrays of numbers as vector objects and
interpret two dimensional arrays of numbers as matrix objects and
both Fortran and C programmers interpret the first and second subscripts
of a two dimensional array as row and column indices respectively
so two dimensional Fortran and C arrays are said to be stored
in column major order and row major order respectively
even if both Fortran and C programs refer to exactly the same object.

Fortran programmers naturally tend to think of one dimensional arrays
as column vectors and two dimensional arrays as collections of column vectors.
C programmers naturally tend to think of one dimensional arrays
as row vectors and two dimensional arrays as collections of row vectors.
Fortran programs can call C vector and matrix functions and
C programs can call Fortran vector and matrix subroutines
but the operation may be interpreted differently.
Suppose, for example, that a Fortran subroutine implements
a matrix-vector multiplication

        Ax

where A is an m by n matrix and x is an n element column vector.
C programmers interpret the Fortran subroutine as an implementation
of a vector-matrix multiplication

        x^TA^T = vM

where v = x^T is an n element row vector and M = A^T is an n by m matrix
but both x and v reference the same vector object in memory
and both A and M reference the same matrix object in memory.

A Fortran subroutine which computes an LDU decomposition

        PA = L(DU)

appears to the C program as

        A^TP^T = (DU)^TL^T

which the C programmer might interpret as

        AP = (LD)U

Similarly, a C function which computes and LDU decomposition

        PA = L(DU)

appears to the Fortran programmer as

        AP = (LD)U

A really good linear algebra subroutine library
would support both decompositions
but Fortran favors the most natural implementation of AP = (LD)U
and C favors the most natural implementation of PA = L(DU)
because they minimize cache misses due to striding.