Intuition about memory layout for fast SIMD / data oriented design

Question

I have been watching data-oriented-design talks recently, but I never understood the reasoning behind their unanimously chosen memory layout.

Lets say we have a 3D animation to render, and in each frame we need to re-normalize our orientation vectors.

The "Scalar code"

They always show code that might look something like this:

let scene = [{"camera1", vec4{1, 1, 1, 1}}, ...]

for object in scene
    object.orientation = normalize(object.orientation)

So far so good... The memory at &scene might look roughly thus:

[string,X,Y,Z,W,string,X,Y,Z,W,string,X,Y,Z,W,...]

"SSE aware code"

Every talk then shows the improved, cookie-cutter, version:

let xs = [1, ...]
let ys = [1, ...]
let zs = [1, ...]
let ws = [1, ...]
let scene = [{"camera1", ptr_vec4{&xs[1], &ys[1], &zs[1], &ws[1]}}, ...]

for (o1, o2, o3, o4) in scene
    (o1, o2, o3, o4) = normalize_sse(o1, o2, o3, o4)

Which, due to it's memory layout, is not only more memory-efficient, but can also process our scene 4 objects at a time.
Memory at &xs, &ys, &zs, and &ws

[X,X,X,X,X,X,...]
[Y,Y,Y,Y,Y,Y,...]
[Z,Z,Z,Z,Z,Z,...]
[W,W,W,W,W,W,...]

But why 4 separate arrays?

If the __m128 (packed-4-singles) is the predominant type in engines,
    which i believe it is;
and if the type is 128 bits long,
    which it definitely is;
and if the cache line width / 128 = 4,
    which it almost always does;
and if x86_64 is only capable of writing a full cache line,
    which I am almost certain of
- why is the data not structured as follows instead?!

Memory at &packed_orientations:

[X,X,X,X,Y,Y,Y,Y,Z,Z,Z,Z,W,W,W,W,X,X,...]
 ^---------cache-line------------^

I have no benchmark to test this on, and I don't understand intrinsics enough to even try, but by my intuition, should this not be way faster? We would be saving 4x page loads and writes, simplifying allocations, and saving pointers, and the code would be simpler since instead of 4 pointers we can do pointer addition. Am I wrong?

Thank you! :)

Answer 1

The amount of data you need to get through your memory subsystem is identical no matter whether you do 4 separate arrays or your suggested interleaving. You therefore don't save page loads or writes (I don't see why the "separate arrays" case should read and write each page or cache line more than once).

You do disperse the memory transfers more - you might have 1 L1 cache miss every iteration in your case and 4 cache misses every 4th iteration in the "separate arrays" case. I don't know which one would be preferred.

Anyway, the main point is to not have unnecessary memory pushed through your caches that you don't interact with. In your example, having string values that are neither read nor written but still pushed through the caches needlessly takes up bandwidth.

Answer 2

One major downside to interleaving on the vector-width is that you need to change the layout to take advantage of wider vectors. (AVX, AVX512).

But yes, when you're purely manually vectorizing (with no loops that a compiler might auto-vectorize with its choice of vector width), this could maybe be worth it if all your (important) loops always use all struct members.

Otherwise Max's point applies: a loop that touches only x and y will waste bandwidth on the z and w members.

---

It won't be way faster, though; with a reasonable amount of loop unrolling, indexing 4 arrays or incrementing 4 pointers is barely worse than 1. HW prefetch on Intel CPUs can track one forward + 1 backward stream per 4k page, so 4 input streams is basically fine.

(But L2 is 4-way associative in Skylake, down from 8 in earlier, so more than 4 input streams all with the same alignment relative to a 4k page would cause conflict misses / defeat prefetching. So with more than 4 large / page-aligned arrays, an interleaved format could avoid that problem.)

For small arrays, where the whole interleaved thing fits in one 4k page, yes that's a potential advantage. Otherwise it's about the same total amount of pages touched and potential TLB misses, with one 4x as often instead of coming in groups of 4. This might well be better for TLB prefetch, if it can be doing one page-walk ahead of time instead of being swamped with multiple TLB misses coming at the same time.

---

Tweaking the SoA struct:

It may help to let the compiler know that the memory pointed to by each of the pointers is non-overlapping. Most C++ compilers (including all 4 major x86 compilers, gcc/clang/MSVC/ICC) support __restrict as a keyword with the same semantics as C99 restrict. Or for portability, use #ifdef / #define to define a restrict keyword as empty or __restrict or whatever, as appropriate for the compiler.

struct SoA_scene {
        size_t size;
        float *__restrict xs;
        float *__restrict ys;
        float *__restrict zs;
        float *__restrict ws;
};

This can definitely help with auto-vectorization, otherwise the compiler doesn't know that xs[i] = foo; doesn't change the value of ys[i+1] for the next iteration.

If you read those vars into local variables (so the compiler is sure that pointer assignments don't modify the pointer itself in the struct), you might declare them as float *__restrict xs = soa.xs; and so on.

The interleaved format inherently avoids this aliasing possibility.

Answer 3

One of the things not mentioned yet is that memory access has quite a bit of latency. And of course, when reading from 4 pointers, the result is available when the last value arrives. So even if 3 out of 4 values are in cache, the last value may need to come from memory and stall your whole operation.

That's why SSE doesn't even support this mode. All your values have to be contiguous in memory, and for quite a while they had to be aligned (so they couldn't cross a cache line boundary).

Importantly, that means your example (Structure of Arrays) does not work in SSE hardware. You can't use element [1] from 4 different vectors in a single operation. You can use elements [0] to [3] from a single vector.

Answer 4

I have implemented a simple benchmark for both methods.

Result: The striped layout is at best 10 % faster than the standard layout*. But with SSE4.1 we can do much better.

*When compiled with gcc -Ofast on an i5-7200U cpu.

The structure is slightly easier to work with, but much less versatile. It might, however, have a bit of an advantage in a real scenario, once the allocator is sufficiently busy.

Striped layout

Time 4624 ms

Memory usage summary: heap total: 713728, heap peak: 713728, stack peak: 2896
         total calls   total memory   failed calls
 malloc|          3         713728              0
realloc|          0              0              0  (nomove:0, dec:0, free:0)
 calloc|          0              0              0
   free|          1         640000

#include <chrono>
#include <cstdio>
#include <random>
#include <vector>
#include <xmmintrin.h>

/* -----------------------------------------------------------------------------
        Striped layout [X,X,X,X,y,y,y,y,Z,Z,Z,Z,w,w,w,w,X,X,X,X...]
----------------------------------------------------------------------------- */

using AoSoA_scene = std::vector<__m128>;

void print_scene(AoSoA_scene const &scene)
{
        // This is likely undefined behavior. Data might need to be stored
        // differently, but this is simpler to index.
        auto &&punned_data = reinterpret_cast<float const *>(scene.data());
        auto scene_size = std::size(scene);

        // Limit to 8 lines
        for(size_t j = 0lu; j < std::min(scene_size, 8lu); ++j) {
                for(size_t i = 0lu; i < 4lu; ++i) {
                        printf("%10.3e ", punned_data[j + 4lu * i]);
                }
                printf("\n");
        }
        if(scene_size > 8lu) {
                printf("(%lu more)...\n", scene_size - 8lu);
        }
        printf("\n");
}

void normalize(AoSoA_scene &scene)
{
        // Euclidean norm, SIMD 4 x 4D-vectors at a time.
        for(size_t i = 0lu; i < scene.size(); i += 4lu) {
                __m128 xs = scene[i + 0lu];
                __m128 ys = scene[i + 1lu];
                __m128 zs = scene[i + 2lu];
                __m128 ws = scene[i + 3lu];

                __m128 xxs = _mm_mul_ps(xs, xs);
                __m128 yys = _mm_mul_ps(ys, ys);
                __m128 zzs = _mm_mul_ps(zs, zs);
                __m128 wws = _mm_mul_ps(ws, ws);

                __m128 xx_yys = _mm_add_ps(xxs, yys);
                __m128 zz_wws = _mm_add_ps(zzs, wws);

                __m128 xx_yy_zz_wws = _mm_add_ps(xx_yys, zz_wws);

                __m128 norms = _mm_sqrt_ps(xx_yy_zz_wws);

                scene[i + 0lu] = _mm_div_ps(xs, norms);
                scene[i + 1lu] = _mm_div_ps(ys, norms);
                scene[i + 2lu] = _mm_div_ps(zs, norms);
                scene[i + 3lu] = _mm_div_ps(ws, norms);
        }
}

float randf()
{
        std::random_device random_device;
        std::default_random_engine random_engine{random_device()};
        std::uniform_real_distribution<float> distribution(-10.0f, 10.0f);
        return distribution(random_engine);
}

int main()
{
        // Scene description, e.g. cameras, or particles, or boids etc.
        // Has to be a multiple of 4!   -- No edge case handling.
        std::vector<__m128> scene(40'000);

        for(size_t i = 0lu; i < std::size(scene); ++i) {
                scene[i] = _mm_set_ps(randf(), randf(), randf(), randf());
        }

        // Print, normalize 100'000 times, print again

        // Compiler is hopefully not smart enough to realize
        // idempotence of normalization
        using std::chrono::steady_clock;
        using std::chrono::duration_cast;
        using std::chrono::milliseconds;
        // >:(

        print_scene(scene);
        printf("Working...\n");

        auto begin = steady_clock::now();
        for(int j = 0; j < 100'000; ++j) {
                normalize(scene);
        }
        auto end = steady_clock::now();
        auto duration = duration_cast<milliseconds>(end - begin);

        printf("Time %lu ms\n", duration.count());
        print_scene(scene);

        return 0;
}

SoA layout

Time 4982 ms

Memory usage summary: heap total: 713728, heap peak: 713728, stack peak: 2992
         total calls   total memory   failed calls
 malloc|          6         713728              0
realloc|          0              0              0  (nomove:0, dec:0, free:0)
 calloc|          0              0              0
   free|          4         640000

#include <chrono>
#include <cstdio>
#include <random>
#include <vector>
#include <xmmintrin.h>

/* -----------------------------------------------------------------------------
        SoA layout [X,X,X,X,...], [y,y,y,y,...], [Z,Z,Z,Z,...], ...
----------------------------------------------------------------------------- */

struct SoA_scene {
        size_t size;
        float *xs;
        float *ys;
        float *zs;
        float *ws;
};

void print_scene(SoA_scene const &scene)
{
        // This is likely undefined behavior. Data might need to be stored
        // differently, but this is simpler to index.

        // Limit to 8 lines
        for(size_t j = 0lu; j < std::min(scene.size, 8lu); ++j) {
                printf("%10.3e ", scene.xs[j]);
                printf("%10.3e ", scene.ys[j]);
                printf("%10.3e ", scene.zs[j]);
                printf("%10.3e ", scene.ws[j]);
                printf("\n");
        }
        if(scene.size > 8lu) {
                printf("(%lu more)...\n", scene.size - 8lu);
        }
        printf("\n");
}

void normalize(SoA_scene &scene)
{
        // Euclidean norm, SIMD 4 x 4D-vectors at a time.
        for(size_t i = 0lu; i < scene.size; i += 4lu) {
                __m128 xs = _mm_load_ps(&scene.xs[i]);
                __m128 ys = _mm_load_ps(&scene.ys[i]);
                __m128 zs = _mm_load_ps(&scene.zs[i]);
                __m128 ws = _mm_load_ps(&scene.ws[i]);

                __m128 xxs = _mm_mul_ps(xs, xs);
                __m128 yys = _mm_mul_ps(ys, ys);
                __m128 zzs = _mm_mul_ps(zs, zs);
                __m128 wws = _mm_mul_ps(ws, ws);

                __m128 xx_yys = _mm_add_ps(xxs, yys);
                __m128 zz_wws = _mm_add_ps(zzs, wws);

                __m128 xx_yy_zz_wws = _mm_add_ps(xx_yys, zz_wws);

                __m128 norms = _mm_sqrt_ps(xx_yy_zz_wws);

                __m128 normed_xs = _mm_div_ps(xs, norms);
                __m128 normed_ys = _mm_div_ps(ys, norms);
                __m128 normed_zs = _mm_div_ps(zs, norms);
                __m128 normed_ws = _mm_div_ps(ws, norms);

                _mm_store_ps(&scene.xs[i], normed_xs);
                _mm_store_ps(&scene.ys[i], normed_ys);
                _mm_store_ps(&scene.zs[i], normed_zs);
                _mm_store_ps(&scene.ws[i], normed_ws);
        }
}

float randf()
{
        std::random_device random_device;
        std::default_random_engine random_engine{random_device()};
        std::uniform_real_distribution<float> distribution(-10.0f, 10.0f);
        return distribution(random_engine);
}

int main()
{
        // Scene description, e.g. cameras, or particles, or boids etc.
        // Has to be a multiple of 4!   -- No edge case handling.
        auto scene_size = 40'000lu;
        std::vector<float> xs(scene_size);
        std::vector<float> ys(scene_size);
        std::vector<float> zs(scene_size);
        std::vector<float> ws(scene_size);

        for(size_t i = 0lu; i < scene_size; ++i) {
                xs[i] = randf();
                ys[i] = randf();
                zs[i] = randf();
                ws[i] = randf();
        }

        SoA_scene scene{
                scene_size,
                std::data(xs),
                std::data(ys),
                std::data(zs),
                std::data(ws)
        };
        // Print, normalize 100'000 times, print again

        // Compiler is hopefully not smart enough to realize
        // idempotence of normalization
        using std::chrono::steady_clock;
        using std::chrono::duration_cast;
        using std::chrono::milliseconds;
        // >:(

        print_scene(scene);
        printf("Working...\n");

        auto begin = steady_clock::now();
        for(int j = 0; j < 100'000; ++j) {
                normalize(scene);
        }
        auto end = steady_clock::now();
        auto duration = duration_cast<milliseconds>(end - begin);

        printf("Time %lu ms\n", duration.count());
        print_scene(scene);

        return 0;
}

AoS layout

Since SSE4.1 there seems to be a third option -- by far the simplest and fastest one.

Time 3074 ms

Memory usage summary: heap total: 746552, heap peak: 713736, stack peak: 2720
         total calls   total memory   failed calls
 malloc|          5         746552              0
realloc|          0              0              0  (nomove:0, dec:0, free:0)
 calloc|          0              0              0
   free|          2         672816
Histogram for block sizes:
    0-15              1  20% =========================
 1024-1039            1  20% =========================
32816-32831           1  20% =========================
   large              2  40% ==================================================

/* -----------------------------------------------------------------------------
        AoS layout [{X,y,Z,w},{X,y,Z,w},{X,y,Z,w},{X,y,Z,w},...]
----------------------------------------------------------------------------- */

using AoS_scene = std::vector<__m128>;

void print_scene(AoS_scene const &scene)
{
        // This is likely undefined behavior. Data might need to be stored
        // differently, but this is simpler to index.
        auto &&punned_data = reinterpret_cast<float const *>(scene.data());
        auto scene_size = std::size(scene);

        // Limit to 8 lines
        for(size_t j = 0lu; j < std::min(scene_size, 8lu); ++j) {
                for(size_t i = 0lu; i < 4lu; ++i) {
                        printf("%10.3e ", punned_data[j * 4lu + i]);
                }
                printf("\n");
        }
        if(scene_size > 8lu) {
                printf("(%lu more)...\n", scene_size - 8lu);
        }
        printf("\n");
}

void normalize(AoS_scene &scene)
{
        // Euclidean norm, SIMD 4 x 4D-vectors at a time.
        for(size_t i = 0lu; i < scene.size(); i += 4lu) {
                __m128 vec = scene[i];
                __m128 dot = _mm_dp_ps(vec, vec, 255);
                __m128 norms = _mm_sqrt_ps(dot);
                scene[i] = _mm_div_ps(vec, norms);
        }
}

float randf()
{
        std::random_device random_device;
        std::default_random_engine random_engine{random_device()};
        std::uniform_real_distribution<float> distribution(-10.0f, 10.0f);
        return distribution(random_engine);
}

int main()
{
        // Scene description, e.g. cameras, or particles, or boids etc.
        std::vector<__m128> scene(40'000);

        for(size_t i = 0lu; i < std::size(scene); ++i) {
                scene[i] = _mm_set_ps(randf(), randf(), randf(), randf());
        }

        // Print, normalize 100'000 times, print again

        // Compiler is hopefully not smart enough to realize
        // idempotence of normalization
        using std::chrono::steady_clock;
        using std::chrono::duration_cast;
        using std::chrono::milliseconds;
        // >:(

        print_scene(scene);
        printf("Working...\n");

        auto begin = steady_clock::now();
        for(int j = 0; j < 100'000; ++j) {
                normalize(scene);
                //break;
        }
        auto end = steady_clock::now();
        auto duration = duration_cast<milliseconds>(end - begin);

        printf("Time %lu ms\n", duration.count());
        print_scene(scene);

        return 0;
}