warming up your workspace

Cache-friendly code, or how swapping two loop lines doubles your speed

Here is a performance fact that surprises people who have written code for years: you can take a correct, optimized nested loop, swap the order of its two lines, change nothing about the result, and make it run twice as slow. No extra work, no different arithmetic, the exact same number of additions. The only thing that changed is the order in which you touched memory, and that order is everything, because the gap between your CPU and your RAM is the dominant cost in most numerical code. Understanding it is the entry fee to high-performance computing, and you can measure it in one file.

The one idea

Your CPU does not read one number from memory at a time. When you ask for a value, the hardware fetches a whole cache line, typically 64 bytes, about eight double-precision numbers, into a small, fast memory close to the core. If your next access is one of those neighbors, it is already there and nearly free. If your next access is far away, the cache line was wasted and the CPU stalls waiting for RAM, which can take a hundred times longer than an add. So fast code is code that, once it touches a piece of memory, uses the neighbors that came along for free before moving on. That means walking memory in the order it is physically laid out.

The catch: how a 2D array is laid out depends on the language. Julia, Fortran, and MATLAB store matrices column-major, one whole column contiguous in memory, then the next. C, C++, and NumPy store them row-major. Get the loop order to match the layout and you fly; mismatch it and you crawl.

Build both

Julia is column-major, so consecutive elements down a column, A[1,j], A[2,j], A[3,j], sit next to each other in memory. The cache-friendly loop therefore varies the first index on the inside.

function sum_columnwise(A)      # inner loop walks DOWN a column: contiguous
    s = 0.0
    @inbounds for j in 1:size(A,2), i in 1:size(A,1)
        s += A[i,j]
    end
    s
end

function sum_rowwise(A)         # inner loop jumps ACROSS columns: strided
    s = 0.0
    @inbounds for i in 1:size(A,1), j in 1:size(A,2)
        s += A[i,j]
    end
    s
end

Look at the two functions. They are identical except for which index the inner loop varies. Both compute the sum of every element. Both do exactly n^2 additions. One is cache-friendly, one is not.

Proof: same answer, half the speed

Run them on a 6000 by 6000 matrix:

column-major order (cache-friendly): 34.5 ms
row-major order    (cache-hostile):  72.8 ms
slowdown from wrong loop order: 2.1x
results identical: true

Same sum, verified equal, more than twice the time. The row-wise version asks for A[1,1], then A[1,2], which is a whole column away in memory, so each access lands on a fresh cache line and uses one number out of the eight it dragged in. The column-wise version marches straight down memory, using all eight before fetching the next line. Nothing changed but the order of two for clauses.

Three details that matter:

  • The speed gap grows with the data. When the matrix fits in cache, the order barely matters. As it outgrows the cache, and real HPC data always does, the penalty for the wrong order climbs, because every wasted cache line now costs a full trip to RAM.
  • @inbounds here just removes bounds-checking so we measure memory behavior, not safety checks. The cache effect is entirely separate and would show up in C, Fortran, or NumPy the same way, with the opposite loop order for the row-major languages.
  • This is why library authors obsess over memory layout. A matrix multiply that respects cache lines and processes data in cache-sized blocks can be ten times faster than the naive triple loop, same operations, radically different access pattern.

Where this shows up

Every high-performance numerical library, BLAS, LAPACK, the kernels under NumPy and Julia, is built on this. The famous optimizations, loop tiling, blocking, choosing the right loop order in matrix multiply, are all in service of the cache. It is why transposing a matrix before a computation can speed it up despite the cost of the transpose, and why the same algorithm can be memory-bound on one machine and compute-bound on another. Before you reach for parallelism or a GPU, getting your access pattern to match your memory layout is often the single biggest win available.

If you want to build the blocked algorithms, measure cache misses directly, and write kernels that saturate memory bandwidth, that is the ground the hpc track on IWTLP covers in Julia, where the cache is a first-class concern.