使用内在指令的欧几里得距离
Euclidean distance using intrinsic instruction
对于一个研究项目，我需要计算很多欧几里得距离，其中必须选择某些维度而丢弃其他维度。在程序的当前状态下，选定维度的数组有 100 个左右的元素，我计算了大约 2-3 百万个距离。我当前的代码如下：
float compute_distance(const float* p1, const float* p2) const
{
__m256 euclidean = _mm256_setzero_ps();
const uint16_t n = nbr_dimensions;
const uint16_t aligend_n = n - n % 16;
const float* local_selected = selected_dimensions;
for (uint16_t i = 0; i < aligend_n; i += 16)
{
const __m256 r1 = _mm256_sub_ps(_mm256_load_ps(&p1[i]), _mm256_load_ps(&p2[i]));
euclidean = _mm256_fmadd_ps(_mm256_mul_ps(r1, r1), _mm256_load_ps(&local_selected[i]), euclidean);
const __m256 r2 = _mm256_sub_ps(_mm256_load_ps(&p1[i + 8]), _mm256_load_ps(&p2[i + 8]));
euclidean = _mm256_fmadd_ps(_mm256_mul_ps(r2, r2), _mm256_load_ps(&local_selected[i + 8]), euclidean);
}
float distance = hsum256_ps_avx(euclidean);
for (uint16_t i = aligend_n; i < n; ++i)
{
const float num = p1[i] - p2[i];
distance += num * num * local_selected[i];
}
return distance;
}
所选尺寸是预先确定的。因此，我可以预先计算 __m256 的数组以传递给 _mm256_blendv_ps，而不是在行 euclidean = _mm256_fmadd_ps(_mm256_mul_ps(r1, r1), _mm256_load_ps(&local_selected[i]), euclidean); 中乘以 0 或 1。但我是 intrinsic instructions 的新手，我还没有找到可行的解决方案。
我想知道你们是否可以提供一些建议，甚至代码建议，以提高此功能的 运行 速度。作为旁注，我无法访问 AVX-512 指令。
更新：
使用上述第一个解决方案，得到：
float compute_distance(const float* p1, const float* p2) const
{
const size_t n = nbr_dimensions;
const size_t aligend_n = n - n % 16;
const unsigned int* local_selected = selected_dimensions;
const __m256* local_masks = masks;
__m256 euc1 = _mm256_setzero_ps(), euc2 = _mm256_setzero_ps(),
euc3 = _mm256_setzero_ps(), euc4 = _mm256_setzero_ps();
const size_t n_max = aligend_n/8;
for (size_t i = 0; i < n_max; i += 4)
{
const __m256 r1 = _mm256_sub_ps(_mm256_load_ps(&p1[i * 8 + 0]), _mm256_load_ps(&p2[i * 8 + 0]));
const __m256 r1_1 = _mm256_and_ps(r1, local_masks[i + 0]);
euc1 = _mm256_fmadd_ps(r1_1, r1_1, euc1);
const __m256 r2 = _mm256_sub_ps(_mm256_load_ps(&p1[i * 8 + 8]), _mm256_load_ps(&p2[i * 8 + 8]));
const __m256 r2_1 = _mm256_and_ps(r2, local_masks[i + 1]);
euc2 = _mm256_fmadd_ps(r2_1, r2_1, euc2);
const __m256 r3 = _mm256_sub_ps(_mm256_load_ps(&p1[i * 8 + 16]), _mm256_load_ps(&p2[i * 8 + 16]));
const __m256 r3_1 = _mm256_and_ps(r3, local_masks[i + 2]);
euc3 = _mm256_fmadd_ps(r3_1, r3_1, euc3);
const __m256 r4 = _mm256_sub_ps(_mm256_load_ps(&p1[i * 8 + 24]), _mm256_load_ps(&p2[i * 8 + 24]));
const __m256 r4_1 = _mm256_and_ps(r4, local_masks[i + 3]);
euc4 = _mm256_fmadd_ps(r4_1, r4_1, euc4);
}
float distance = hsum256_ps_avx(_mm256_add_ps(_mm256_add_ps(euc1, euc2), _mm256_add_ps(euc3, euc4)));
for (size_t i = aligend_n; i < n; ++i)
{
const float num = p1[i] - p2[i];
distance += num * num * local_selected[i];
}
return distance;
}
基本建议：
不要使用 uint16_t 作为你的循环计数器，除非你真的想强制编译器每次都 t运行cate 到 16 位。至少使用 unsigned，或者有时使用 uintptr_t（或更常规地，size_t）可以获得更好的 asm。从 32 位到指针宽度的零扩展仅通过使用 32 位操作数大小的 asm 指令在 x86-64 上免费发生，但有时编译器仍然做得不好。
使用五个或更多单独的累加器而不是一个 euclidean，因此多个 sub / FMA 指令可以运行，而不会在将 FMA 放入一个累加器的循环携带依赖链的延迟上出现瓶颈。
FMA 在 Intel Haswell 上的延迟为 5 个周期，但吞吐量为每 0.5 个周期一个。另请参阅 latency vs throughput in intel intrinsics, and also my answer on 以获得更高级的版本。
避免通过全局变量传递参数。显然你的 n 是一个编译时常量（这很好），但 selected_dimensions 不是，对吗？如果是，那么你在整个程序中只使用了一组掩码，所以不要在意下面关于压缩掩码的内容。
使用全局变量可能会破坏编译器优化，因为它将您的函数内联到调用者中，调用者在调用它之前设置了全局变量。 （通常只有在设置全局变量和使用它之间存在非内联函数调用时才会发生，但这并不少见。）
更新：你的数组很小，只有约 100 个元素，所以只展开 2 个元素可能很好，以减少启动/清理开销。乱序执行可以在这么短的迭代次数内隐藏 FMA 延迟，特别是如果不需要此函数调用的最终结果来决定下一次调用的输入参数。
总的函数调用开销很重要，而不仅仅是大型数组的向量化效率有多高。
，剥离循环的第一次迭代可以让您通过初始化 euc1 = stuff(p1[0], p2[0]); 而不是 _mm256_setzero_ps().
来避免第一个 FMA
用零将数组填充为完整向量（甚至是 2 个向量的完整展开循环体）可以让您完全避免标量清理循环，并使整个函数非常紧凑。
如果你不能只是填充，你仍然可以通过加载一个一直到输入末尾的未对齐向量来避免标量清理，并屏蔽它以避免重复计算。 （请参阅 以了解基于未对齐计数生成掩码的方法）。在您编写输出数组的其他类型的问题中，重做重叠元素没问题。
您没有显示您的 hsum256_ps_avx 代码，但这是总延迟和函数吞吐量的相当一部分。确保针对吞吐量对其进行优化：例如避免 haddps / _mm_hadd_ps。请参阅我在 Fastest way to do horizontal float vector sum on x86 上的回答。
您的具体情况：
I could thus pre-compute an array of __m256 to pass to a _mm256_blendv_ps instead of multiplying by 0 or 1 in the FMA.
是的，那会更好，尤其是如果它可以让您将其他内容折叠到 FMAdd / FMSub 中。但更好的是，使用全零或全一的布尔值 _mm256_and_ps。这会使值保持不变 (1 & x == x) 或归零 (0 & x == 0，并且 float 0.0 的二进制表示全为零。)
如果您的面具在缓存中没有丢失，请将它们完全解压存储，以便可以加载它们。
如果您使用具有相同 p1 和 p2 的不同掩码，您可以预先计算 p1-p2 平方，然后只进行掩码 add_ps减少。 （但请注意，在 Intel pre-Skylake 上，FMA 的吞吐量优于 ADD。Haswell/Broadwell 有 2 个 FMA 单元，但 运行 ADDPS 在延迟较低的专用单元上（3c 与 5c）。只有一个vector-FP 添加单元。Skylake 只是 运行s FMA 单元上的所有内容，具有 4 个周期延迟。）无论如何，这意味着使用 FMA 作为 1.0 * x + y 实际上可以赢得吞吐量。但是您可能没问题，因为您仍然需要分别加载掩码和 square(p1-p2)，因此每个 FP 添加 2 次加载，因此每个周期一次跟上加载吞吐量。除非您（或编译器）在前面剥离一些迭代并将这些迭代的浮点数据保存在跨多个不同 local_selected 掩码的寄存器中。
update：我写这篇文章时假设数组大小是 2-3 百万，而不是 ~100。 L1D 缓存未命中的配置文件，以确定是否值得花费更多 CPU 指令来减少缓存占用空间。如果你总是对所有 300 万次调用使用相同的掩码，那么压缩可能不值得。
您可以将掩码压缩到每个元素 8 位并使用 pmovsx (_mm256_cvtepi8_epi32) 加载它们（对全一值进行符号扩展会产生更宽的全一，因为这就是 2 的补码 -1 作品）。不幸的是，将它用作负载很烦人；编译器有时无法将 _mm256_cvtepi8_epi32(_mm_cvtsi64x_si128(foo)) 优化为 vpmovsxbd ymm0, [mem]，而是实际使用单独的 vmovq 指令。
const uint64_t *local_selected = something; // packed to 1B per element
__m256 euc1 = _mm256_setzero_ps(), euc2 = _mm256_setzero_ps(),
euc3 = _mm256_setzero_ps(), euc4 = _mm256_setzero_ps();
for (i = 0 ; i < n ; i += 8*4) { // 8 floats * an unroll of 4
__m256 mask = _mm256_castsi256_ps( _mm256_cvtepi8_epi32(_mm_cvtsi64x_si128(local_selected[i*1 + 0])) );
// __m256 mask = _mm256_load_ps(local_selected[i*8 + 0]); // without packing
const __m256 r1 = _mm256_sub_ps(_mm256_load_ps(&p1[i*8 + 0]), _mm256_load_ps(&p2[i*8 + 0]));
r1 = _mm256_and_ps(r1, mask); // zero r1 or leave it untouched.
euc1 = _mm256_fmadd_ps(r1, r1, euc1); // euc1 += r1*r1
// ... same for r2 with local_selected[i + 1]
// and p1/p2[i*8 + 8]
// euc2 += (r2*r2) & mask2
// and again for euc3 (local_selected[i + 2], p1/p2[i*8 + 16]
// and again for euc3 (local_selected[i + 3], p1/p2[i*8 + 24]
}
euclidean = hsum (euc1+euc2+euc3+euc4);
我猜你在没有 pmovsx 的情况下在负载吞吐量上略有瓶颈，因为你有三个负载用于三个向量 ALU 操作。 （对于微融合，在 Intel CPU 上总共只有 4 个融合域 uops，因此它在前端没有瓶颈）。并且三个 ALU uops 可以在不同的端口上 运行（vandps 在 Intel pre-Skylake 上对于端口 5 是 1 uop。在 SKL 上它可以在任何端口上 运行）。
添加随机播放（pmovsx）可能会在端口 5（Haswell/Broadwell 上）出现瓶颈。您可能想要使用 vpand 作为掩码，这样它就可以在任何端口上 运行，如果您正在调整 HSW/BDW，即使它们在整数 AND 和FP 数学指令。有了足够的累加器，您就不会受到延迟限制。 （Skylake 有 VANDPS 的额外旁路延迟，具体取决于它碰巧 运行 在哪个端口上。
blendv 比 AND 慢：总是至少 2 微码。
为大型数组进一步压缩掩码
如果您的数组大于 L2 缓存，并且您的掩码数组具有与浮点数组一样多的元素，您很可能会在加载带宽上遇到瓶颈（至少在使用多个向量累加器展开后）。这意味着花更多的指令解包掩码数据是值得的，以减少那部分带宽需求。
我认为掩码数据的理想格式是交错的 32 个掩码向量，这使得动态 "unpack" 非常便宜。使用移位将正确的掩码带入每个 32 位元素的高位，并将其与 vblendvps 一起使用，通过与零混合来有条件地将元素置零。 （或者算术右移+布尔AND）
__m256i masks = _mm256_load_si256(...);
// this actually needs a cast to __m256, omitted for readability
r0 = _mm256_blendv_ps(_mm256_setzero_ps(), r0, masks);
...
__m256i mask1 = _mm256_slli_epi32(masks, 1);
r1 = _mm256_blendv_ps(_mm256_setzero_ps(), r1, mask1);
...
__m256i mask2 = _mm256_slli_epi32(masks, 2);
r2 = _mm256_blendv_ps(_mm256_setzero_ps(), r2, mask2);
...
// fully unrolling is overkill; you can set up for a loop back to r0 with
masks = _mm256_slli_epi32(masks, 4);
你也可以在每一步都做masks = _mm256_slli_epi32(masks, 1);，这可能会更好，因为它使用了 1 个寄存器。但它可能对导致掩码 dep 链延迟的资源冲突更敏感，因为每个掩码都依赖于前一个掩码。
Intel Haswell 运行仅在端口 5 上有两个 vblendvps 微指令，因此您可以考虑使用 _mm256_srai_epi32 + _mm256_and_ps。但是 Skylake 可以 运行 任何 p015 上的 2 微指令，所以混合在那里很好（尽管它确实占用了一个包含全零向量的向量寄存器）。
使用打包比较生成交错格式的掩码，然后 _mm256_srli_epi32(cmp_result, 31) 并将其或运算到您正在构建的向量中。然后左移1。重复32次。
如果数组中的完整数据向量少于 32 个，您仍然可以使用此格式。较低的位将不会被使用。或者您可以为每个向量设置 2 个或更多 selected_dimensions 个掩码。例如每个元素的前 16 位用于一个 selected_dimensions，而后 16 位用于另一个。你可以做类似
的事情
__m256i masks = _mm256_load_si256(dimensions[selector/2]);
masks = _mm256_sll_epi32(masks, 16 * (selector % 2));
// or maybe
if (selector % 2) {
masks = _mm256_slli_epi32(masks, 16);
}
AVX512:
AVX512可以直接使用位图掩码，所以效率更高一些。只需使用 const __mmask16 *local_selected = whatever; 声明一个 16 位掩码数组（用于 16 个浮点数的 512b 向量），并使用 r0 = _mm512_maskz_sub_ps(p1,p2, local_selected[i]); 对减法进行零掩码。
如果加载端口 uop 吞吐量确实存在瓶颈（每个时钟 2 个加载），您可以尝试一次加载 64 位掩码数据，并使用掩码移位获得不同的低 16 位数据。不过，除非您的数据在 L1D 缓存中很热，否则这可能不会成为问题。
首先使用 compare-into-mask 生成掩码数据非常容易，不需要交织。
理想情况下，您可以缓存阻止调用它的代码，这样您就可以在缓存中很热的时候重用数据。例如从 p1 和 p2 的前 64kiB 中获取所有你想要的组合，然后移动到后面的元素并在它们在缓存中很热的时候执行它们。
对于一个研究项目，我需要计算很多欧几里得距离，其中必须选择某些维度而丢弃其他维度。在程序的当前状态下，选定维度的数组有 100 个左右的元素，我计算了大约 2-3 百万个距离。我当前的代码如下：
float compute_distance(const float* p1, const float* p2) const
{
__m256 euclidean = _mm256_setzero_ps();
const uint16_t n = nbr_dimensions;
const uint16_t aligend_n = n - n % 16;
const float* local_selected = selected_dimensions;
for (uint16_t i = 0; i < aligend_n; i += 16)
{
const __m256 r1 = _mm256_sub_ps(_mm256_load_ps(&p1[i]), _mm256_load_ps(&p2[i]));
euclidean = _mm256_fmadd_ps(_mm256_mul_ps(r1, r1), _mm256_load_ps(&local_selected[i]), euclidean);
const __m256 r2 = _mm256_sub_ps(_mm256_load_ps(&p1[i + 8]), _mm256_load_ps(&p2[i + 8]));
euclidean = _mm256_fmadd_ps(_mm256_mul_ps(r2, r2), _mm256_load_ps(&local_selected[i + 8]), euclidean);
}
float distance = hsum256_ps_avx(euclidean);
for (uint16_t i = aligend_n; i < n; ++i)
{
const float num = p1[i] - p2[i];
distance += num * num * local_selected[i];
}
return distance;
}
所选尺寸是预先确定的。因此，我可以预先计算 __m256 的数组以传递给 _mm256_blendv_ps，而不是在行 euclidean = _mm256_fmadd_ps(_mm256_mul_ps(r1, r1), _mm256_load_ps(&local_selected[i]), euclidean); 中乘以 0 或 1。但我是 intrinsic instructions 的新手，我还没有找到可行的解决方案。
我想知道你们是否可以提供一些建议，甚至代码建议，以提高此功能的 运行 速度。作为旁注，我无法访问 AVX-512 指令。
更新： 使用上述第一个解决方案，得到：
float compute_distance(const float* p1, const float* p2) const
{
const size_t n = nbr_dimensions;
const size_t aligend_n = n - n % 16;
const unsigned int* local_selected = selected_dimensions;
const __m256* local_masks = masks;
__m256 euc1 = _mm256_setzero_ps(), euc2 = _mm256_setzero_ps(),
euc3 = _mm256_setzero_ps(), euc4 = _mm256_setzero_ps();
const size_t n_max = aligend_n/8;
for (size_t i = 0; i < n_max; i += 4)
{
const __m256 r1 = _mm256_sub_ps(_mm256_load_ps(&p1[i * 8 + 0]), _mm256_load_ps(&p2[i * 8 + 0]));
const __m256 r1_1 = _mm256_and_ps(r1, local_masks[i + 0]);
euc1 = _mm256_fmadd_ps(r1_1, r1_1, euc1);
const __m256 r2 = _mm256_sub_ps(_mm256_load_ps(&p1[i * 8 + 8]), _mm256_load_ps(&p2[i * 8 + 8]));
const __m256 r2_1 = _mm256_and_ps(r2, local_masks[i + 1]);
euc2 = _mm256_fmadd_ps(r2_1, r2_1, euc2);
const __m256 r3 = _mm256_sub_ps(_mm256_load_ps(&p1[i * 8 + 16]), _mm256_load_ps(&p2[i * 8 + 16]));
const __m256 r3_1 = _mm256_and_ps(r3, local_masks[i + 2]);
euc3 = _mm256_fmadd_ps(r3_1, r3_1, euc3);
const __m256 r4 = _mm256_sub_ps(_mm256_load_ps(&p1[i * 8 + 24]), _mm256_load_ps(&p2[i * 8 + 24]));
const __m256 r4_1 = _mm256_and_ps(r4, local_masks[i + 3]);
euc4 = _mm256_fmadd_ps(r4_1, r4_1, euc4);
}
float distance = hsum256_ps_avx(_mm256_add_ps(_mm256_add_ps(euc1, euc2), _mm256_add_ps(euc3, euc4)));
for (size_t i = aligend_n; i < n; ++i)
{
const float num = p1[i] - p2[i];
distance += num * num * local_selected[i];
}
return distance;
}
基本建议：
不要使用 uint16_t 作为你的循环计数器，除非你真的想强制编译器每次都 t运行cate 到 16 位。至少使用 unsigned，或者有时使用 uintptr_t（或更常规地，size_t）可以获得更好的 asm。从 32 位到指针宽度的零扩展仅通过使用 32 位操作数大小的 asm 指令在 x86-64 上免费发生，但有时编译器仍然做得不好。
使用五个或更多单独的累加器而不是一个 euclidean，因此多个 sub / FMA 指令可以运行，而不会在将 FMA 放入一个累加器的循环携带依赖链的延迟上出现瓶颈。
FMA 在 Intel Haswell 上的延迟为 5 个周期，但吞吐量为每 0.5 个周期一个。另请参阅 latency vs throughput in intel intrinsics, and also my answer on
避免通过全局变量传递参数。显然你的 n 是一个编译时常量（这很好），但 selected_dimensions 不是，对吗？如果是，那么你在整个程序中只使用了一组掩码，所以不要在意下面关于压缩掩码的内容。
使用全局变量可能会破坏编译器优化，因为它将您的函数内联到调用者中，调用者在调用它之前设置了全局变量。 （通常只有在设置全局变量和使用它之间存在非内联函数调用时才会发生，但这并不少见。）
更新：你的数组很小，只有约 100 个元素，所以只展开 2 个元素可能很好，以减少启动/清理开销。乱序执行可以在这么短的迭代次数内隐藏 FMA 延迟，特别是如果不需要此函数调用的最终结果来决定下一次调用的输入参数。
总的函数调用开销很重要，而不仅仅是大型数组的向量化效率有多高。
euc1 = stuff(p1[0], p2[0]); 而不是 _mm256_setzero_ps().
用零将数组填充为完整向量（甚至是 2 个向量的完整展开循环体）可以让您完全避免标量清理循环，并使整个函数非常紧凑。
如果你不能只是填充，你仍然可以通过加载一个一直到输入末尾的未对齐向量来避免标量清理，并屏蔽它以避免重复计算。 （请参阅
您没有显示您的 hsum256_ps_avx 代码，但这是总延迟和函数吞吐量的相当一部分。确保针对吞吐量对其进行优化：例如避免 haddps / _mm_hadd_ps。请参阅我在 Fastest way to do horizontal float vector sum on x86 上的回答。
您的具体情况：
I could thus pre-compute an array of __m256 to pass to a
_mm256_blendv_psinstead of multiplying by 0 or 1 in the FMA.
是的，那会更好，尤其是如果它可以让您将其他内容折叠到 FMAdd / FMSub 中。但更好的是，使用全零或全一的布尔值 _mm256_and_ps。这会使值保持不变 (1 & x == x) 或归零 (0 & x == 0，并且 float 0.0 的二进制表示全为零。)
如果您的面具在缓存中没有丢失，请将它们完全解压存储，以便可以加载它们。
如果您使用具有相同 p1 和 p2 的不同掩码，您可以预先计算 p1-p2 平方，然后只进行掩码 add_ps减少。 （但请注意，在 Intel pre-Skylake 上，FMA 的吞吐量优于 ADD。Haswell/Broadwell 有 2 个 FMA 单元，但 运行 ADDPS 在延迟较低的专用单元上（3c 与 5c）。只有一个vector-FP 添加单元。Skylake 只是 运行s FMA 单元上的所有内容，具有 4 个周期延迟。）无论如何，这意味着使用 FMA 作为 1.0 * x + y 实际上可以赢得吞吐量。但是您可能没问题，因为您仍然需要分别加载掩码和 square(p1-p2)，因此每个 FP 添加 2 次加载，因此每个周期一次跟上加载吞吐量。除非您（或编译器）在前面剥离一些迭代并将这些迭代的浮点数据保存在跨多个不同 local_selected 掩码的寄存器中。
update：我写这篇文章时假设数组大小是 2-3 百万，而不是 ~100。 L1D 缓存未命中的配置文件，以确定是否值得花费更多 CPU 指令来减少缓存占用空间。如果你总是对所有 300 万次调用使用相同的掩码，那么压缩可能不值得。
您可以将掩码压缩到每个元素 8 位并使用 pmovsx (_mm256_cvtepi8_epi32) 加载它们（对全一值进行符号扩展会产生更宽的全一，因为这就是 2 的补码 -1 作品）。不幸的是，将它用作负载很烦人；编译器有时无法将 _mm256_cvtepi8_epi32(_mm_cvtsi64x_si128(foo)) 优化为 vpmovsxbd ymm0, [mem]，而是实际使用单独的 vmovq 指令。
const uint64_t *local_selected = something; // packed to 1B per element
__m256 euc1 = _mm256_setzero_ps(), euc2 = _mm256_setzero_ps(),
euc3 = _mm256_setzero_ps(), euc4 = _mm256_setzero_ps();
for (i = 0 ; i < n ; i += 8*4) { // 8 floats * an unroll of 4
__m256 mask = _mm256_castsi256_ps( _mm256_cvtepi8_epi32(_mm_cvtsi64x_si128(local_selected[i*1 + 0])) );
// __m256 mask = _mm256_load_ps(local_selected[i*8 + 0]); // without packing
const __m256 r1 = _mm256_sub_ps(_mm256_load_ps(&p1[i*8 + 0]), _mm256_load_ps(&p2[i*8 + 0]));
r1 = _mm256_and_ps(r1, mask); // zero r1 or leave it untouched.
euc1 = _mm256_fmadd_ps(r1, r1, euc1); // euc1 += r1*r1
// ... same for r2 with local_selected[i + 1]
// and p1/p2[i*8 + 8]
// euc2 += (r2*r2) & mask2
// and again for euc3 (local_selected[i + 2], p1/p2[i*8 + 16]
// and again for euc3 (local_selected[i + 3], p1/p2[i*8 + 24]
}
euclidean = hsum (euc1+euc2+euc3+euc4);
我猜你在没有 pmovsx 的情况下在负载吞吐量上略有瓶颈，因为你有三个负载用于三个向量 ALU 操作。 （对于微融合，在 Intel CPU 上总共只有 4 个融合域 uops，因此它在前端没有瓶颈）。并且三个 ALU uops 可以在不同的端口上 运行（vandps 在 Intel pre-Skylake 上对于端口 5 是 1 uop。在 SKL 上它可以在任何端口上 运行）。
添加随机播放（pmovsx）可能会在端口 5（Haswell/Broadwell 上）出现瓶颈。您可能想要使用 vpand 作为掩码，这样它就可以在任何端口上 运行，如果您正在调整 HSW/BDW，即使它们在整数 AND 和FP 数学指令。有了足够的累加器，您就不会受到延迟限制。 （Skylake 有 VANDPS 的额外旁路延迟，具体取决于它碰巧 运行 在哪个端口上。
blendv 比 AND 慢：总是至少 2 微码。
为大型数组进一步压缩掩码
如果您的数组大于 L2 缓存，并且您的掩码数组具有与浮点数组一样多的元素，您很可能会在加载带宽上遇到瓶颈（至少在使用多个向量累加器展开后）。这意味着花更多的指令解包掩码数据是值得的，以减少那部分带宽需求。
我认为掩码数据的理想格式是交错的 32 个掩码向量，这使得动态 "unpack" 非常便宜。使用移位将正确的掩码带入每个 32 位元素的高位，并将其与 vblendvps 一起使用，通过与零混合来有条件地将元素置零。 （或者算术右移+布尔AND）
__m256i masks = _mm256_load_si256(...);
// this actually needs a cast to __m256, omitted for readability
r0 = _mm256_blendv_ps(_mm256_setzero_ps(), r0, masks);
...
__m256i mask1 = _mm256_slli_epi32(masks, 1);
r1 = _mm256_blendv_ps(_mm256_setzero_ps(), r1, mask1);
...
__m256i mask2 = _mm256_slli_epi32(masks, 2);
r2 = _mm256_blendv_ps(_mm256_setzero_ps(), r2, mask2);
...
// fully unrolling is overkill; you can set up for a loop back to r0 with
masks = _mm256_slli_epi32(masks, 4);
你也可以在每一步都做masks = _mm256_slli_epi32(masks, 1);，这可能会更好，因为它使用了 1 个寄存器。但它可能对导致掩码 dep 链延迟的资源冲突更敏感，因为每个掩码都依赖于前一个掩码。
Intel Haswell 运行仅在端口 5 上有两个 vblendvps 微指令，因此您可以考虑使用 _mm256_srai_epi32 + _mm256_and_ps。但是 Skylake 可以 运行 任何 p015 上的 2 微指令，所以混合在那里很好（尽管它确实占用了一个包含全零向量的向量寄存器）。
使用打包比较生成交错格式的掩码，然后 _mm256_srli_epi32(cmp_result, 31) 并将其或运算到您正在构建的向量中。然后左移1。重复32次。
如果数组中的完整数据向量少于 32 个，您仍然可以使用此格式。较低的位将不会被使用。或者您可以为每个向量设置 2 个或更多 selected_dimensions 个掩码。例如每个元素的前 16 位用于一个 selected_dimensions，而后 16 位用于另一个。你可以做类似
__m256i masks = _mm256_load_si256(dimensions[selector/2]);
masks = _mm256_sll_epi32(masks, 16 * (selector % 2));
// or maybe
if (selector % 2) {
masks = _mm256_slli_epi32(masks, 16);
}
AVX512:
AVX512可以直接使用位图掩码，所以效率更高一些。只需使用 const __mmask16 *local_selected = whatever; 声明一个 16 位掩码数组（用于 16 个浮点数的 512b 向量），并使用 r0 = _mm512_maskz_sub_ps(p1,p2, local_selected[i]); 对减法进行零掩码。
如果加载端口 uop 吞吐量确实存在瓶颈（每个时钟 2 个加载），您可以尝试一次加载 64 位掩码数据，并使用掩码移位获得不同的低 16 位数据。不过，除非您的数据在 L1D 缓存中很热，否则这可能不会成为问题。
首先使用 compare-into-mask 生成掩码数据非常容易，不需要交织。
理想情况下，您可以缓存阻止调用它的代码，这样您就可以在缓存中很热的时候重用数据。例如从 p1 和 p2 的前 64kiB 中获取所有你想要的组合，然后移动到后面的元素并在它们在缓存中很热的时候执行它们。