How Much Linear Memory Access Is Enough?

For basically any high-performance computation, memory layout and access pattern are critical. Common wisdom is that linear, contiguous memory performs best and should almost always be preferred. However, it should be intuitively clear that this has diminishing returns: processing a single 32 GB block vs processing two 16 GB blocks will not meaningfully differ in performance. Working with smaller blocks enables some interesting data structures, so I've set out to experimentally determine what block size is needed to effectively capture the full performance.

Findings

Setup and detailed analysis below, but my personal takeaway is:

• 1 MB blocks are enough for basically any workload of this kind
• 128 kB blocks suffice once you have at least ~1 cycle per processed byte
• 4 kB blocks are already enough once you're above roughly ~10 cycles per processed byte

(for raw data processing, not necessarily if there are other per-block costs)

This is the full results chart for my Ryzen 9 7950X3D, effectively showing the block sizes needed for peak performance across different workloads. The rest of this post will go over the setup and discuss a few isolated graphs.

Code and results are available here: github.com/solidean/bench-linear-access

Setup

The memory hierarchy of modern CPUs is famously complex, so I've tried to create an experimental setup where we can isolate and control the effects well enough to make our results generalize.

Our main question is: when we have to process a data set in blocks of linear memory, how large do the individual blocks have to be so that any penalty in jumping between blocks is amortized?

For that, our "processing kernel" takes a span<span<float const> const> (or float**) and returns an opaque uint64_t result (to make sure the optimizer cannot elide our computation). This is effectively the non-owning version of "vector of blocks".

For example, the scalar_stats kernel is

// Scalar stats kernel: computes running stats over all blocks.
// Returns a mini hash (no quality requirements) to prevent elision.
uint64_t kernel_scalar_stats(std::span<std::span<float const> const> data)
{
struct stats
{
float m0 = 0, m1 = 0, m2 = 0;
float min = +std::numeric_limits<float>::max();
float max = -std::numeric_limits<float>::max();
};
stats s;
for (auto block : data)
for (auto d : block)
{
s.m0 += 1;
s.m1 += d;
s.m2 += d * d;
if (d < s.min)
s.min = d;
if (d > s.max)
s.max = d;
}
auto b = [](float f) { return std::bit_cast<uint32_t>(f); };
return b(s.m0) ^ b(s.m1) ^ b(s.m2) ^ b(s.min) ^ b(s.max);
}

This represents a lightweight but scalar computation, reaching ~7 GB/s on my system (for large blocks). So roughly 3 cycles per float, or 0.75 cycles per byte.

The return value is simply some value dependent on all computation. We do a similar aggregation over all experiments and write the resulting "hash" into the .csv. That constraints the optimizer enough without going through volatile.

A single run is a single call of the kernel function on a set of blocks. Each block is contiguous and has a fixed "block size". The total size of all blocks that we execute on is the "working set". We can partition the same working set using different block sizes, which is the degree of freedom we have in practice.

Benchmarking on a modern system is notoriously noisy, so we need to repeat the same run multiple times. We chose to report medians to guard against outliers on either end. However, if the working set fits into a cache, we would measure warm-ish timings while our real use case might work on cold data.

To counteract this, we place the blocks in random order and random positions inside a sufficiently large "backing memory" (4 GB here).

We still clobber the cache before each run:

// folds a 256 MB block of memory into the seed and returns a dependent result
uint64_t clobber_cache(uint64_t seed)
{
static std::vector<uint64_t> clobber_data;
if (clobber_data.empty())
{
clobber_data.resize((256uLL << 20) / sizeof(uint64_t));
std::default_random_engine rng;
for (auto& d : clobber_data)
d = rng();
}
// slightly dependency heavy but very hard to optimize away
for (auto d : clobber_data)
seed = seed * d + 7;
return seed;
}

Clobbering is not part of the timing. Its purpose is simply to ensure no useful data remains in the cache hierarchy before each run.

The scalar_stats kernel across block sizes

The first experiment runs the above scalar_stats kernel across varying block sizes (32 byte to 2 MB) and working sets (1 MB to 64 MB).

While there is some noise in the measurements, this graph shows a few things:

• Our layout randomization and clobbering is effective. There is basically no effect of the cache hierarchy here and the results are effectively independent of working set size.
• Beyond 128 kB there is little benefit to increasing the block size.
• We converge to slightly beyond 7 GB/s here, which is ~3 cycles per float (or ~0.75 cycles per byte).

Each individual run starts with a randomized layout and a cold cache, so the curve shape cleanly captures the contiguity vs. jumping around in memory tradeoff.

The repeated variant of the previous experiment

We get one of these scenarios by relaxing the "coldness" a bit. The "repeated" scenario only clobbers and chooses a random layout at the start of each set of K runs. This models a scenario where we hit the data repeatedly and any data or TLB cache that fits can be reused.

Unsurprisingly, smaller working sets can be (partially) kept in the cache hierarchy, which leads to earlier peaking. The conservative 128 kB works for all sizes but you might get away with 4..16 kB blocks if your total data size is small.

Comparing different workloads

To represent the high-end SIMD processing, we have simd_sum, an AVX2 (or NEON) sum over all floats:

uint64_t kernel_simd_sum(std::span<std::span<float const> const> data)
{
__m256 acc = _mm256_setzero_ps();
for (auto block : data)
{
auto const* ptr = block.data();
auto const count = block.size();
for (size_t i = 0; i < count; i += 8)
acc = _mm256_add_ps(acc, _mm256_load_ps(ptr + i));
}
// reduce: bitcast to uint32 array and XOR
auto words = std::bit_cast<std::array<uint32_t, 8>>(acc);
uint32_t h = 0;
for (auto w : words)
h ^= w;
return uint64_t(h);
}

This reaches 50+ GB/s on my Ryzen, which is roughly 0.1 cycles/byte.