我有一个函数，它接受两个数组并用两个数组坐标的双线性组合填充第三个数组。如果您好奇的话，该函数是 8 维实向量空间上的克利福德乘积。
为了尽可能快地完成此操作，我预先计算了定义产品所需的“访问模式”的索引。每个访问模式都是一个无符号 32 位整数。这是一个访问模式示例：
00001100.00000100.00001000.00000001
假设输入数组是arr1和arr2，输出数组是ret，这个模式的含义如下：
ret[00001100] += arr1[00000100] * arr2[00001000] * (-1)^00000001;
(-1)^s 项实际上并未计算。相反，我只是根据最后 8 位是否为 0 进行分支。非零值表示 -= 而不是 +=。有些位被浪费了，但拥有 32 位模式似乎是值得的。
访问模式存储在一个大的排序数组中。写入被捆绑在一起，因此我只需要每个块访问 ret[i] 一次 - 据推测，这通过让 CPU 在写入内存之前从 arr1 和 arr2 积累所需的所有内容来提高缓存局部性。说得不同，而不是做：
ret[i] += arr1[j]*arr2[k]*(-1)^s;
ret[i] += arr1[j']*arr2[k']*(-1)^s';
...
我愿意：
double x = 0.0;
x += arr1[j]*arr2[k]*(-1)^s;
x += arr1[j']*arr2[k']*(-1)^s';
...
ret[i] += x;
在具有常数 i 的块上每次迭代一次。
在上面，访问模式是：(i,j,k,s), (i,j',k',s'),...当 i 改变时，我移动到下一个迭代。
我想利用 SIMD 和指令级并行性。这就是我认为的样子：
对于ILP，添加更多中间变量：
double x = 0.0;
double y = 0.0;
x += arr1[j]*arr2[k]*(-1)^s;
y += arr1[j']*arr2[k']*(-1)^s';
x += arr1[j'']*arr2[k'']*(-1)^s'';
y += arr1[j''']*arr2[k''']*(-1)^s''';
...
ret[i] += x + y;
我只是在每次迭代中采用更多的访问模式，并且我试图说服 CPU 并行修改 x 和 y。我认为，向其中添加 SIMD 需要在每次迭代中采用更多的访问模式，但这只是我所了解的。在我开始深入研究 SIMD 操作之前，请帮助引导我走上正确的道路。编写这样的函数时应该注意什么？如何才能使其尽可能快？
我认为你的算法将在内存加载事务上遇到瓶颈。大多数现代 CPU 每个周期只能执行 2 个负载事务。因为您是从随机位置加载的，所以每个加载的元素都需要 1 次加载。您的目标是保存这些负载交易。
在您写的注释中，您的访问模式数组有 130k 元素。这意味着要计算完整的内容，您需要从两个输入数组加载 260k 个元素，再加上访问模式表中的 130k 个整数。
但是，对于 2 个输入（每个输入 256 个元素），输入元素只有 64k 个唯一乘积。请注意，256×256 元素的 2D 数组需要 512 kb 内存，适合 CPU 缓存。
考虑将算法分为两步。
第一步，计算这 64k 个独特的乘积。为了通过 AVX 展开 4×16 块来有效地实现它，这样你只需要从第一个数组加载 64 个向量，从第二个数组加载 1k 个向量，并存储 16k 个向量。
然后在第二步中，对于输出向量的每个元素，您需要累积
±products[ i ]uint32_t这是一个如何在没有分支的情况下实现第二步来应用符号的示例。该代码未经测试。
与您的算法相比，此方法节省了大约 1/3 的内存负载。
// Increment or decrement acc by an element from source array
// Index to load the element is in bytes 1-2 of the packed integer
// The sign is in the lowest bit of that integer
inline void addSigned( __m128d& acc, const double* source, uint32_t packed )
{
// Load element from the temporary buffer
__m128d inc = _mm_load_sd( source + ( ( packed >> 8 ) & 0xFFFFu ) );
// Make FP64 scalar with either 0.0 or -0.0 depending on eth requested sign
__m128i iv = _mm_cvtsi32_si128( packed & 1 );
iv = _mm_slli_epi64( iv, 63 );
// Conditionally negate the increment, and accumulate
inc = _mm_xor_pd( inc, _mm_castsi128_pd( iv ) );
acc = _mm_add_sd( acc, inc );
}
// Reduction table
const std::array<uint32_t, 130000> s_table;
void computeStuff( double* result, const double* productsMatrix )
{
// Assuming the table is sparse, fill output vector with zeros
const __m256d zero = _mm256_setzero_pd();
for( size_t i = 0; i < 0x100; i += 4 )
_mm256_storeu_pd( result + i, zero );
uint32_t destElement = 0x100;
__m128d acc = _mm_setzero_pd();
for( uint32_t packed : s_table )
{
const uint32_t rdi = packed >> 24;
if( rdi != destElement )
{
if( destElement < 0x100 )
_mm_store_sd( result + destElement, acc );
acc = _mm_setzero_pd();
destElement = rdi;
}
addSigned( acc, productsMatrix, packed );
}
if( destElement < 0x100 )
_mm_store_sd( result + destElement, acc );
}
另一个想法，编写一个工具将访问模式表编译成相当高效的 AVX C++ 源代码。通过相当高效，我的意思是尽可能使用全矢量加载，并重用加载的元素来计算多个输出值。
不幸的是，如果您的访问模式确实是随机的，则生成的函数将太长，将成为指令获取和解码的瓶颈。然而，如果您的访问模式不是真正随机的，并且您可以滥用对称性、输入向量的循环排列或类似技巧将算法拆分为处理不同数据的同一函数的多个调用，则该方法可能会有所帮助。