optimization - N 个向量的内存布局

Question

我想知道我应该如何存储 N 个小维度的向量(比如 X、Y、Z)以提高效率。

出于缓存局部性的原因，我期望将向量一个接一个地打包 [N][3](主要行)会产生比布局 [3][N] 更好的结果(其中维度 X，Y 然后Z 在使用 OpenMP 进行矢量运算时连续布局。

但是，将每个向量乘以 3x3 矩阵，并使用英特尔 MKL BLAS，我发现布局 [3][N] 的速度是原来的两倍。

我猜缓存局部性被 SSE 指令用于非跨步数据这一事实所抵消。这让我想知道人们(例如在计算机图形学中)如何存储他们的向量，以及是否还有其他优缺点。

Answer 1

通常使用两种数据布局:结构数组 (AoS) 和数组结构 (SoA)。

服务质量:

struct
{
   float x;
   float y;
   float z;
} points[N];

SoA:

struct
{
   float x[N];
   float y[N];
   float z[N];
} points;

为了将 AoS 情况下的每个点与 3x3 矩阵 M 相乘，循环体如下所示:

r[i].x = M[0][0]*points[i].x +
         M[0][1]*points[i].y +
         M[0][2]*points[i].z;
// ditto for r[i].y and r[i].z

SSE 可以一次乘以 4 个 float (AVX 可以乘以 8 个 float )并且它还提供点积运算，但问题是在向量寄存器中加载 3 个 float 是非常低效的操作。可以添加一个额外的 float 字段来填充结构，但仍然会损失 1/4 的计算能力，因为两个向量操作数中的第 4 个 float 未被使用(或不包含有用信息)。您也可以不对这些点进行矢量化，例如一次处理 4 个点，因为一次加载 points[i].x 到 points[i+3].x 需要收集加载，目前还不支持在 x86 上(尽管一旦支持 AVX2 的 CPU 可用，这种情况就会改变)。

在 SoA 的情况下，内循环是:

r.x[i] = M[0][0]*points.x[i] +
         M[0][1]*points.y[i] +
         M[0][2]*points.z[i];
// ditto for r.y[i] and r.z[i]

基本上看起来是一样的，但是有一个很关键的区别。现在，编译器可以使用矢量指令并一次处理 4 个点(如果使用 AVX，甚至可以处理 8 个点)。例如。它可以展开循环并将其转换为以下等效向量:

<r.x[i], r.x[i+1], r.x[i+2], r.x[i+3]> =
    M[0][0]*<x[i], x[i+1], x[i+2], x[i+3]> +
    M[0][1]*<y[i], y[i+1], y[i+2], y[i+3]> +
    M[0][2]*<z[i], z[i+1], z[i+2], z[i+3]>

这里有三个向量加载、三个标量向量乘法、三个向量加法和一个向量存储。它们都 100% 地利用了 SSE 的矢量功能。唯一的问题是当点数不能被 4 整除时，但可以轻松填充数组，否则编译器可能会生成标量代码以串行执行剩余的迭代。无论哪种方式，如果您有很多点，只为剩余的 1 到 3 个点损失一些性能比在每个点不断地利用硬件要有益得多。

另一方面，如果您经常需要访问随机点的 (x,y,z) 元组，那么 SoA 实现将导致三个缓存行读取(如果数据没有适合缓存)，而 AoS 实现将需要一个或两个(使用填充可以避免需要两次加载的情况)。所以答案是 - 数据结构取决于哪种操作主导算法。