N个向量的内存布局

Question

我想知道如何高效地存储N个小维度向量（比如X，Y，Z）。

为了缓存本地性的原因，我原本希望将向量按顺序打包[N] [3]（行主要）而不是以维度X、Y、Z顺序排列[3] [N]。我认为这样在使用OpenMP进行向量操作时会得到更好的结果。

然而，在使用Intel MKL BLAS计算每个向量与3×3矩阵相乘时，我发现布局[3][N]速度快了两倍。

我猜缓存本地性被SSE指令对非跨步数据的作用所抵消。这让我想知道人们（比如计算机图形学领域）是如何存储他们的向量的，是否有其他优缺点。

- pch

Answer 1

有两种常见的数据布局：结构数组（AoS）和数组结构（SoA）。

AoS：

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个浮点数（AVX可以一次性乘以8个浮点数），并且它还提供了点积运算，但问题在于将3个浮点数加载到向量寄存器中的效率非常低。可以向结构添加一个额外的float字段来填充结构，但仍然会损失1/4的计算能力，因为两个向量操作数中的第四个浮点数未使用（或不包含有用信息）。您也不能对点进行矢量化处理，例如一次处理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]>

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

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

- Hristo Iliev