Row-Major Layout and the Address Arithmetic Behind Indexing

A two-dimensional array in C is not a grid. int a[4][8] declares an array of four objects, each of which is an array of eight int. Every property that matters about layout follows from that single sentence: the elements of one row are adjacent, the rows are adjacent to each other, and the address of a[i][j] is a linear function of i and j. Fortran and the BLAS libraries chose the mirror image of this rule, and the mismatch between the two conventions is a recurring source of wrong results when numerical code crosses a language boundary.

The subscript is one multiply and one add

The element index of a[i][j] in a row-major array with C columns is i * C + j. Multiplying by sizeof(int) converts that to a byte offset. On the common 32-bit and 64-bit machines where int is four bytes, the scaling is a shift by two, so the whole address computation is a shift, a multiply by the row length, and an add. The multiply by C is the only term that depends on the array’s shape; it is why a compiler that knows C at compile time can strength-reduce a row walk into a pointer increment.

static inline size_t elem_index(size_t i, size_t j, size_t cols) {
    return i * cols + j;          /* element offset, not bytes */
}

static inline int *at(int *base, size_t i, size_t j, size_t cols) {
    return base + elem_index(i, j, cols);   /* scales by sizeof(int) */
}

Nothing in the C standard contains the phrase “row-major”. The layout is forced instead by the type rules. a[i] is equivalent to *(a + i), and the pointer arithmetic on a advances in units of sizeof(int[8]), which is 32 bytes for the example above. Dereferencing that pointer yields an array, which decays to a pointer to its first element, and [j] then advances by sizeof(int). Because pointer arithmetic is defined in terms of whole objects laid out consecutively, the second subscript must vary fastest in memory. The same reasoning applies at any nesting depth: the rightmost subscript always has stride one element.

The declared type also fixes the stride between rows. For int a[4][8] the row stride is 32 bytes regardless of how many rows are actually used, so a function that accepts int (*a)[8] and a function that accepts a flat int * with a separate cols argument describe different layouts unless cols happens to be 8.

The compiler’s knowledge of the row length is what makes the difference. A parameter declared as int (*a)[8] carries the row length in its type, so a[i][j] compiles to a multiply by the constant 8 and a shift by 2. A parameter declared as int *a with a runtime cols carries nothing, so every access computes i * cols + j with a real multiply, and the compiler cannot hoist the multiply out of an inner loop unless it can prove cols does not change. The two-argument form is also what allows the compiler to reason about aliasing: with int (*a)[8] it knows the extent of each row, while with a bare pointer it knows only that some number of ints follow. The restrict qualifier, added in C99, is the explicit promise that two pointer parameters do not overlap, and it is the mechanism that recovers the aliasing information the type no longer carries.

The transposed convention and where it breaks

Fortran stores the leftmost subscript contiguously, so a(i,j) sits at (j-1)*rows + (i-1) in one-based terms. This is not an arbitrary historical choice: it matches the way numerical analysts wrote matrix formulas, and the reference BLAS interface inherited it. Every BLAS routine that takes a matrix takes a leading dimension and assumes column-major storage. There is no layout parameter in the classic interface.

Passing a row-major buffer to a column-major routine is the classic interoperation bug. For a square matrix the routine silently computes with the transpose, which for symmetric operations such as a matrix-vector product may produce a plausible-looking but wrong vector. For a non-square matrix the leading-dimension arithmetic walks outside the allocation. Two workarounds exist: store the data column-major in the first place, or call a wrapper that transposes on entry. CBLAS, the C binding to BLAS, exposes a CblasRowMajor/CblasColMajor parameter, but that parameter describes the caller’s buffer; the underlying Fortran routine still receives a column-major view, and the wrapper may copy.

Traversal order and the cache line

The cache line is 64 bytes on the x86-64 and AArch64 machines that dominate server and desktop hardware. A line therefore holds sixteen consecutive 32-bit integers. For a 1024×1024 int matrix, one row is 4096 bytes, exactly 64 lines, and the whole matrix is 4 MiB.

A loop that fixes i and varies j consumes one line per sixteen elements: the first access to a line misses, the next fifteen hit. A loop that fixes j and varies i touches a different line on every iteration, because consecutive i values are 4096 bytes apart. Each of those accesses maps to a distinct line, and the column’s lines are spread across the entire 4 MiB, which exceeds the L1 and L2 capacities of typical cores. The result is one cache miss per element instead of one per sixteen, a factor of sixteen in lines fetched, plus a translation lookaside buffer pressure problem discussed in Streaming Through Large Arrays.

Hardware prefetchers mitigate but do not remove this. A stream prefetcher detects a constant stride and issues requests ahead of the demand loads. A 4096-byte stride is constant, so the prefetcher can follow it, but the requests it generates still consume memory bandwidth and still miss in cache; the latency is hidden, the traffic is not.

Blocking converts the bad traversal into a sequence of good ones. Processing a TILE × TILE square at a time means the working set of the inner pair of loops is TILE rows of TILE elements from each array the loop touches. A three-array elementwise kernel at TILE = 64 with 32-bit elements holds three such tiles, 48 KiB, which is larger than the L1 data cache on many current cores but well inside L2; the block still converts each line’s use from once per pass into TILE uses while the block is resident. Sizing TILE by the cache level being targeted is the tuning step: TILE of 16 with 32-bit elements gives three 1 KiB tiles that fit in L1, at the cost of more loop and bounds overhead. The same total number of elements is visited either way, but each line is reused many times before eviction instead of once.

#define TILE 64

void add_tiled(size_t n, const float *restrict a,
               const float *restrict b, float *restrict c) {
    for (size_t ii = 0; ii < n; ii += TILE)
        for (size_t jj = 0; jj < n; jj += TILE)
            for (size_t i = ii; i < ii + TILE && i < n; i++)
                for (size_t j = jj; j < jj + TILE && j < n; j++)
                    c[i * n + j] = a[i * n + j] + b[i * n + j];
}

Tiling costs address arithmetic: the index expression i * n + j is recomputed rather than incremented, and the loop bounds carry the partial-tile predicates. For a memory-bound kernel the trade is favourable by a wide margin; for a kernel that already fits in L1 the tiling overhead is pure loss.

What the loop order changes for a transposed walk

A loop that must visit a row-major array column by column cannot be tiled without changing what it computes, because the tiles are square only if the access pattern permits. Two options remain. The first is to transpose the array once into a second buffer, paying 2·n·m element moves and n·m extra storage, after which both the original and the transposed traversal are sequential in one of the two buffers. The second is to accept the strided access and rely on the prefetcher, which follows the constant stride but cannot reduce the traffic: each element still costs a line.

The choice is a space-for-time trade with a clear threshold. Transposing pays n·m element moves once; the strided walk pays a line fetch per element on every pass. For a single pass over the array the transpose is pure loss. For a kernel that makes many passes — an iterative solver, for example — the transpose is amortized over the passes and wins. The deciding quantity is the number of passes, and it is usually known from the algorithm rather than measured.

Where the index computation goes wrong

Three failure modes recur in code that computes offsets by hand.

Overflow in the element count. i * cols + j is computed in size_t, which wraps modulo 2^64 on a 64-bit platform. If the product exceeds the address space, the wrapped value points somewhere else entirely and the access is undefined rather than trapped. On a platform with 32-bit size_t, an index above 2^32 − 1 cannot be represented at all, so the product wraps and the access addresses the wrong element long before the address space is exhausted. The allocation size has the same problem: malloc(n * sizeof(T)) computes the product in size_t and passes the wrapped result, so an overflow produces a small allocation and a later out-of-range write. calloc takes the element count and the element size as separate arguments, and the C standard says an implementation should detect a multiplication overflow there; the common C libraries do check it and return a null pointer instead of a short block, which is the reason to prefer calloc over malloc plus memset when the element count comes from untrusted input.

Signed indices. a[i - 1] with i of type size_t and i == 0 wraps to SIZE_MAX and produces a pointer far outside the object. The resulting access is undefined behaviour, not a bounds error; it may read valid memory belonging to another object, which is why the bug survives testing. Compilers exploit the undefined-ness: the assumption that a pointer derived from an array stays within it licenses optimizations that remove the check that would have caught the mistake.

Pointer arithmetic past the end. Forming the address one past the last element is legal; forming the address two past it, or one before the first, is not, even if the pointer is never dereferenced. A loop written as for (p = a + n; p != a; --p) use(p[-1]) stays inside the rule. A loop that computes a + n + 1 as a sentinel does not.

Choosing the convention deliberately

Layout is a property of the producer and the consumer, and the two must agree. A program that only ever traverses its own arrays can pick either convention as long as it is consistent. A program that hands buffers to an external library inherits that library’s convention, and the cost of disagreeing ranges from a copy to silent corruption. When the array is one-dimensional, the question disappears: the row-major and column-major rules coincide, and the layout reduces to the flat sequence discussed in What a Slice Descriptor Holds.

For arrays of records rather than scalars, the same address arithmetic applies with the element size replaced by the struct’s stride, including its padding; the consequences for a scan that reads one field are covered in Array of Structs or Struct of Arrays.