Memory Layout & Roofline

Same matrix, two layouts, 6× bandwidth difference

A 2-D matrix doesn't actually live as 2-D in memory — RAM is one big 1-D array. The choice of how to flatten it changes which threads can read it efficiently.

Row-major (C order) — addresses march along a row: M[i][j] and M[i][j+1] are adjacent in memory. Default in C/C++, Python/NumPy, PyTorch, MLX, Metal.
Column-major (Fortran order) — addresses march down a column: M[i][j] and M[i+1][j] are adjacent. Default in Fortran, MATLAB, R, and (because BLAS was specified in 1979) cuBLAS.

This matters because GPUs read memory in coalesced chunks: when 32 neighboring threads in a warp ask for 32 contiguous addresses, the hardware fetches them in one transaction. If those threads stride non-contiguously, the hardware needs many transactions to satisfy them — bandwidth tanks.

Concrete number from a 4096×4096 FP32 sum on RTX 4090:

Row-major, threads scan rows (coalesced): 610 GB/s
Column-major, threads scan rows (strided): 105 GB/s

That's a 6× swing caused entirely by memory access pattern.

Roofline: combining bandwidth + compute into one ceiling

The roofline model says: achievable performance is bounded by

min(peak FLOPs, peak bandwidth × arithmetic intensity).

Plot the kernel's intensity on the X axis and its achieved FLOPs on the Y axis; you're either under the slanted bandwidth ceiling (memory-bound, low intensity) or under the flat compute ceiling (compute-bound, high intensity). Knowing which side of the roofline knee you're on tells you which optimization to chase.

Code

Row-major vs column-major in NumPy — same data, different strides·python

import numpy as np

# Row-major (default)
A = np.zeros((4, 5), order='C')
print(A.strides)   # (40, 8) — step 40 bytes between rows, 8 between cols

# Column-major
B = np.zeros((4, 5), order='F')
print(B.strides)   # (8, 32) — step 8 bytes between rows, 32 between cols

# When you summing rows in a row-major matrix, the inner loop
# walks contiguous memory (good).
# When you sum columns in a row-major matrix, the inner loop
# strides by row-stride bytes per step (bad).
# GPUs amplify this difference because warps coalesce 32-wide.

# A.sum(axis=1) on row-major  — fast
# A.sum(axis=0) on row-major  — slower
# (Reverse for column-major)

Roofline arithmetic — what ceiling am I under?·python

# RTX 4090 specs (FP32 non-Tensor-Core)
PEAK_FLOPS_FP32 = 82e12       # 82 TFLOP/s
PEAK_BW         = 1.0e12      # 1.0 TB/s GDDR6X

def roofline(intensity):
    return min(PEAK_FLOPS_FP32, PEAK_BW * intensity)

# A small AXPY: 0.08 FLOPs/byte
print(roofline(0.08))      # ~ 80 GFLOP/s — bandwidth-bound

# A 4096³ GEMM: ~1365 FLOPs/byte
print(roofline(1365))      # capped at 82 TFLOP/s — compute-bound

# Knee of the roofline (where bandwidth = compute):
knee = PEAK_FLOPS_FP32 / PEAK_BW
print(f'Knee at intensity = {knee:.1f} FLOPs/byte')   # ~ 82

Exercise

In a Python REPL, time A.sum(axis=0) vs A.sum(axis=1) for a 4096×4096 row-major NumPy matrix. You should see the axis=1 (along rows) version is meaningfully faster despite doing the same number of FLOPs. Then create the column-major copy B = np.asfortranarray(A) and re-time — the bias flips. This is the same effect that makes coalesced GPU access 6× faster than strided.

Same matrix, two layouts, 6× bandwidth difference

Roofline: combining bandwidth + compute into one ceiling

Code

External links

Exercise

Progress

Comments 0