关于快速 SIMD / 面向数据设计的内存布局的直觉
Intuition about memory layout for fast SIMD / data oriented design
我最近一直在看面向数据的设计讲座，但我一直不明白他们一致选择内存布局背后的原因。
假设我们有一个 3D 动画要渲染，并且在每一帧中我们需要重新标准化我们的方向向量。
"Scalar code"
他们总是显示可能看起来像这样的代码：
let scene = [{"camera1", vec4{1, 1, 1, 1}}, ...]
for object in scene
object.orientation = normalize(object.orientation)
到目前为止一切顺利... &scene 处的内存可能大致如下所示：
[string,X,Y,Z,W,string,X,Y,Z,W,string,X,Y,Z,W,...]
"SSE aware code"
Every talk then shows the improved,千篇一律，版本：
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)
由于它的内存布局，它不仅内存效率更高，而且还可以一次处理我们场景中的 4 个对象。
&xs、&ys、&zs 和 &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,...]
但为什么是 4 个单独的数组？
如果 __m128（packed-4-singles）是引擎中的主要类型，
我相信是这样；
如果类型是 128 位长，
这绝对是；
如果 缓存行宽度 / 128 = 4,
它几乎总是这样做；
如果 x86_64 只能写入完整的缓存行，
我几乎可以肯定
- 为什么数据的结构不是如下所示？！
&packed_orientations 处的内存：
[X,X,X,X,Y,Y,Y,Y,Z,Z,Z,Z,W,W,W,W,X,X,...]
^---------cache-line------------^
我没有基准来测试它，而且我对内在函数的理解还不足以尝试，但根据我的直觉，这不应该是 方式 更快？
我们将节省 4 倍的页面加载和写入、简化分配和保存指针，并且代码会更简单，因为我们可以进行指针加法而不是 4 个指针。我错了吗？
谢谢！ :)
无论您是使用 4 个单独的数组还是建议的交错，您需要通过内存子系统获取的数据量都是相同的。因此，您不会保存页面加载或写入（我不明白为什么 "separate arrays" 案例应该多次读取和写入每个页面或缓存行）。
你确实分散了内存传输——在你的情况下，你的每次迭代可能有 1 个 L1 缓存未命中，在 "separate arrays" 的情况下，每第 4 次迭代有 4 个缓存未命中。我不知道哪个更受欢迎。
无论如何，要点是不要将不必要的内存推送到您不与之交互的缓存中。在您的示例中，string 既不读取也不写入但仍通过缓存推送的值会不必要地占用带宽。
尚未提及的一件事是内存访问有相当多的延迟。当然，当从 4 个指针读取时，结果在 last 值到达时可用。因此，即使 4 个值中有 3 个在缓存中，最后一个值可能需要来自内存并停止整个操作。
这就是 SSE 甚至不支持这种模式的原因。您所有的值在内存中都必须是连续的，并且在很长一段时间内它们必须对齐（因此它们不能跨越缓存行边界）。
重要的是，这意味着您的示例（数组结构）在 SSE 硬件中不起作用。您不能在单个操作中使用来自 4 个不同向量的元素 [1]。您可以使用单个向量中的元素[0]到[3]。
我已经为这两种方法实施了一个简单的基准测试。
结果：条纹布局最多比标准布局快 10%*。但是有了 SSE4.1，我们可以做得更好。
*在 i5-7200U cpu 上使用 gcc -Ofast 编译时。
结构 使用起来稍微容易一些，但通用性差得多。但是，一旦分配器足够繁忙，它在实际场景中可能会有一些优势。
条纹布局
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布局
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 布局
从 SSE4.1 开始，似乎有第三种选择 -- 迄今为止最简单、最快的选择。
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;
}
矢量宽度交错的一个主要缺点是您需要更改布局以利用更宽的矢量。 (AVX, AVX512).
但是，是的，当您纯手动矢量化时（没有编译器可能会根据其选择的矢量宽度自动矢量化的循环），如果您所有（重要的）循环始终使用 all结构成员。
否则 Max 的观点适用：只涉及 x 和 y 的循环将浪费 z 和 w 成员的带宽。
它不会方式更快，但是；通过合理数量的循环展开，索引 4 个数组或递增 4 个指针几乎不比 1 差。Intel CPU 上的硬件预取可以跟踪每 4k 页一个前向 + 1 个后向流，所以 4 个输入流基本上没问题。
（但 L2 在 Skylake 中是 4 向关联的，低于之前的 8 向关联，因此超过 4 个输入流都具有相同的相对于 4k 页面的对齐方式会导致冲突未命中/预取失败。因此超过4 个大/页面对齐的数组，交错格式可以避免这个问题。）
对于小型阵列，其中整个交错的东西适合一个 4k 页，是的，这是一个潜在的优势。否则，触摸的页面总数和潜在的 TLB 未命中大致相同，一个 4 倍，而不是 4 个一组。这对于 TLB 预取可能更好，如果它可以提前进行一页遍历被同时出现的多个 TLB 未命中所淹没。
调整 SoA 结构：
让编译器知道每个指针指向的内存是非重叠的可能会有所帮助。大多数 C++ 编译器（包括所有 4 个主要的 x86 编译器，gcc/clang/MSVC/ICC）支持 __restrict 作为与 C99 restrict 具有相同语义的关键字。或者为了可移植性，使用 #ifdef / #define 将 restrict 关键字定义为空或 __restrict 或其他适合编译器的关键字。
struct SoA_scene {
size_t size;
float *__restrict xs;
float *__restrict ys;
float *__restrict zs;
float *__restrict ws;
};
这绝对有助于自动矢量化，否则编译器不知道 xs[i] = foo; 不会为下一次迭代更改 ys[i+1] 的值。
如果您将这些 var 读入局部变量（因此编译器确保指针赋值不会修改结构中的指针本身），您可以将 它们 声明为 float *__restrict xs = soa.xs;等等。
交错格式本质上避免了这种混叠的可能性。
我最近一直在看面向数据的设计讲座，但我一直不明白他们一致选择内存布局背后的原因。
假设我们有一个 3D 动画要渲染，并且在每一帧中我们需要重新标准化我们的方向向量。
"Scalar code"
他们总是显示可能看起来像这样的代码：
let scene = [{"camera1", vec4{1, 1, 1, 1}}, ...]
for object in scene
object.orientation = normalize(object.orientation)
到目前为止一切顺利... &scene 处的内存可能大致如下所示：
[string,X,Y,Z,W,string,X,Y,Z,W,string,X,Y,Z,W,...]
"SSE aware code"
Every talk then shows the improved,千篇一律，版本：
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)
由于它的内存布局，它不仅内存效率更高，而且还可以一次处理我们场景中的 4 个对象。
&xs、&ys、&zs 和 &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,...]
但为什么是 4 个单独的数组？
如果 __m128（packed-4-singles）是引擎中的主要类型，
我相信是这样；
如果类型是 128 位长，
这绝对是；
如果 缓存行宽度 / 128 = 4,
它几乎总是这样做；
如果 x86_64 只能写入完整的缓存行，
我几乎可以肯定
- 为什么数据的结构不是如下所示？！
&packed_orientations 处的内存：
[X,X,X,X,Y,Y,Y,Y,Z,Z,Z,Z,W,W,W,W,X,X,...]
^---------cache-line------------^
我没有基准来测试它，而且我对内在函数的理解还不足以尝试，但根据我的直觉，这不应该是 方式 更快？
我们将节省 4 倍的页面加载和写入、简化分配和保存指针，并且代码会更简单，因为我们可以进行指针加法而不是 4 个指针。我错了吗？
谢谢！ :)
无论您是使用 4 个单独的数组还是建议的交错，您需要通过内存子系统获取的数据量都是相同的。因此，您不会保存页面加载或写入（我不明白为什么 "separate arrays" 案例应该多次读取和写入每个页面或缓存行）。
你确实分散了内存传输——在你的情况下，你的每次迭代可能有 1 个 L1 缓存未命中，在 "separate arrays" 的情况下，每第 4 次迭代有 4 个缓存未命中。我不知道哪个更受欢迎。
无论如何，要点是不要将不必要的内存推送到您不与之交互的缓存中。在您的示例中，string 既不读取也不写入但仍通过缓存推送的值会不必要地占用带宽。
尚未提及的一件事是内存访问有相当多的延迟。当然，当从 4 个指针读取时，结果在 last 值到达时可用。因此，即使 4 个值中有 3 个在缓存中，最后一个值可能需要来自内存并停止整个操作。
这就是 SSE 甚至不支持这种模式的原因。您所有的值在内存中都必须是连续的，并且在很长一段时间内它们必须对齐（因此它们不能跨越缓存行边界）。
重要的是，这意味着您的示例（数组结构）在 SSE 硬件中不起作用。您不能在单个操作中使用来自 4 个不同向量的元素 [1]。您可以使用单个向量中的元素[0]到[3]。
我已经为这两种方法实施了一个简单的基准测试。
结果：条纹布局最多比标准布局快 10%*。但是有了 SSE4.1，我们可以做得更好。
*在 i5-7200U cpu 上使用 gcc -Ofast 编译时。
结构 使用起来稍微容易一些，但通用性差得多。但是，一旦分配器足够繁忙，它在实际场景中可能会有一些优势。
条纹布局
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布局
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 布局
从 SSE4.1 开始，似乎有第三种选择 -- 迄今为止最简单、最快的选择。
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;
}
矢量宽度交错的一个主要缺点是您需要更改布局以利用更宽的矢量。 (AVX, AVX512).
但是，是的，当您纯手动矢量化时（没有编译器可能会根据其选择的矢量宽度自动矢量化的循环），如果您所有（重要的）循环始终使用 all结构成员。
否则 Max 的观点适用：只涉及 x 和 y 的循环将浪费 z 和 w 成员的带宽。
它不会方式更快，但是；通过合理数量的循环展开，索引 4 个数组或递增 4 个指针几乎不比 1 差。Intel CPU 上的硬件预取可以跟踪每 4k 页一个前向 + 1 个后向流，所以 4 个输入流基本上没问题。
（但 L2 在 Skylake 中是 4 向关联的，低于之前的 8 向关联，因此超过 4 个输入流都具有相同的相对于 4k 页面的对齐方式会导致冲突未命中/预取失败。因此超过4 个大/页面对齐的数组，交错格式可以避免这个问题。）
对于小型阵列，其中整个交错的东西适合一个 4k 页，是的，这是一个潜在的优势。否则，触摸的页面总数和潜在的 TLB 未命中大致相同，一个 4 倍，而不是 4 个一组。这对于 TLB 预取可能更好，如果它可以提前进行一页遍历被同时出现的多个 TLB 未命中所淹没。
调整 SoA 结构：
让编译器知道每个指针指向的内存是非重叠的可能会有所帮助。大多数 C++ 编译器（包括所有 4 个主要的 x86 编译器，gcc/clang/MSVC/ICC）支持 __restrict 作为与 C99 restrict 具有相同语义的关键字。或者为了可移植性，使用 #ifdef / #define 将 restrict 关键字定义为空或 __restrict 或其他适合编译器的关键字。
struct SoA_scene {
size_t size;
float *__restrict xs;
float *__restrict ys;
float *__restrict zs;
float *__restrict ws;
};
这绝对有助于自动矢量化，否则编译器不知道 xs[i] = foo; 不会为下一次迭代更改 ys[i+1] 的值。
如果您将这些 var 读入局部变量（因此编译器确保指针赋值不会修改结构中的指针本身），您可以将 它们 声明为 float *__restrict xs = soa.xs;等等。
交错格式本质上避免了这种混叠的可能性。