c++ - 欧氏距离使用内在指令

Question

对于一个研究项目，我需要计算很多欧几里得距离，其中必须选择某些尺寸而将其他尺寸丢弃。在程序的当前状态下，选定维度的数组包含100个ish元素，我计算出大约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。但是我不是内在指令的新手，而且还没有找到可行的解决方案。

我想知道你们是否可以提出一些建议，甚至代码建议，以提高此功能的运行速度。附带说明，我无权访问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;
}

Answer 1

基本建议:

不要将uint16_t用作循环计数器，除非您真的想强制编译器每次都将其截断为16位。至少使用unsigned，否则有时您可以通过使用uintptr_t(或更常见的是size_t)获得更好的汇编。仅从使用32位操作数大小的asm指令就可以在x86-64上免费实现从32位到指针宽度的零扩​​展，但是有时候编译器仍然做不到。

使用五个或更多单独的累加器而不是一个euclidean，因此可以运行多个sub / FMA指令，而不会阻塞将FMA集成到一个累加器中的循环依赖链的延迟。

FMA在Intel Haswell上具有5个周期的延迟，但每0.5个周期1个吞吐量。另请参阅latency vs throughput in intel intrinsics和我对Why does mulss take only 3 cycles on Haswell, different from Agner's instruction tables?的回答，以获取更高级的版本。

避免将args传递给全局变量。显然，您的n是一个编译时常量(很好)，但是selected_dimensions不是，不是吗？如果是这样，那么您在整个程序中只使用一组掩码，因此不要忘记下面有关压缩掩码的内容。

当使用全局变量将函数内联到调用程序中时会破坏编译器的优化，该调用程序会在调用之前设置全局变量。 (通常仅在设置全局变量和使用全局变量之间存在非内联函数调用，但这并不少见。)

更新:您的数组很小，只有约100个元素，因此仅展开2个元素可能很好，以减少启动/清理开销。乱序执行可能会在这短暂的迭代中隐藏FMA延迟，尤其是在不需要此函数调用的最终结果来确定下一个调用的输入参数的情况下。

总的函数调用开销很重要，而不仅仅是大型数组矢量化的效率。

As discussed in comments，剥离循环的第一次迭代，可以通过初始化euc1 = stuff(p1[0], p2[0]);而不是_mm256_setzero_ps()来避免第一次FMA。

将数组填充为零的完整 vector (甚至是2个 vector 的完整展开循环体)，可以完全避免标量清理循环，并使整个函数非常紧凑。

如果不能只是填充，您仍然可以通过加载一直到输入末尾的未对齐 vector 来避免标量清理，并对其进行掩码以避免重复计数。 (请参见this Q&A以获取基于未对齐计数生成蒙版的方法)。在编写输出数组的其他类型的问题中，可以重做重叠的元素。

您没有显示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中。但是甚至比这更好的是，使用具有全零或全一的 bool(boolean) _mm256_and_ps。这使值保持不变( 1 & x == x)或清零( 0 & x == 0，并且float 0.0的二进制表示形式为全零。)

如果您的面具在高速缓存中没有丢失，请完全拆开包装存放，以便可以加载它们。

如果您使用具有相同 p1和 p2的不同蒙版，则可以预先计算平方的 p1-p2，然后对此进行蒙版 add_ps缩减。 (但请注意，在Intel Skylake之前，FMA具有比ADD更好的吞吐量。H​​aswell/ Broadwell有2个FMA单元，但在专用单元上运行ADDPS时延较低(3c对5c)。只有一个vector-FP加法单元。 Skylake只需在FMA单元上以4个周期的延迟运行所有内容。)无论如何，这意味着将FMA用作 1.0 * x + y实际上可能是吞吐量的胜利。但是您可能还不错，因为您仍然需要分别加载掩码和 square(p1-p2)，因此每个FP添加2次加载，因此每个周期进行一次加载就可以满足加载吞吐量。除非您(或编译器)在前面剥离一些迭代并将这些迭代的浮点数据保存在跨多个不同 local_selected掩码的寄存器中。

更新:我是在假设数组大小为2-3百万而不是〜100的情况下编写的。 L1D缓存的配置文件无法决定是否花更多的CPU指令来减少缓存占用空间是否值得。如果您始终对300万个 call 使用相同的掩码，那么压缩它可能不值得。

您可以将掩码压缩为每个元素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的情况下会稍微影响瓶颈，因为对于三个 vector ALU操作，您有三个负载。 (使用微融合，在Intel CPU上总共只有4个融合域uops，因此在前端没有瓶颈)。而且这三个ALU操作系统可以在不同的端口上运行(对于Intel pre-Skylake上的端口5， vandps是1 uop。在SKL上，它可以在任何端口上运行)。

在port5上(在Haswell / Broadwell上)添加随机播放( pmovsx)可能存在的瓶颈。您可能想要使用 vpand进行屏蔽，因此，即使您要调整HSW / BDW，即使它们在整数AND和FP数学指令之间有额外的旁路等待时间，它也可以在任何端口上运行。有了足够的累加器，您就不会受到延迟的限制。 (Skylake根据其运行在哪个端口上而为 VANDPS提供了额外的旁路延迟)。

blendv比AND慢:始终至少2 oups。

对于大型阵列，可进一步压缩蒙版

如果您的数组大于L2高速缓存，并且您的mask数组具有与float数组一样多的元素，则很可能会成为加载带宽的瓶颈(至少在使用多个 vector 累加器展开之后)。这意味着花费更多的指令来解压缩掩码数据值得减少带宽需求。

我认为您的 mask 数据的理想格式是32个交错的 mask vector ，这使得在运行中“解压缩”非常便宜。使用移位将正确的掩码带入每个32位元素的高位，并与 vblendvps一起使用，以通过与0混合来有条件地将元素归零。 (或算术右移+ bool(boolean) 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个寄存器的使用。但是，由于每个掩码都依赖于前一个掩码，因此它对引起掩码延迟链延迟的资源冲突可能更敏感。

英特尔Haswell仅在port5上运行两个 vblendvps uops，因此您可以考虑使用 _mm256_srai_epi32 + _mm256_and_ps。但是Skylake可以在p015的任何一个上运行2 oups，因此blend在那里很好(尽管它确实占用了保存全零 vector 的 vector 寄存器)。

生成带有打包比较的交错格式的 mask ，然后将 _mm256_srli_epi32(cmp_result, 31)和或与要构建的 vector 相乘。然后将其左移1。重复32次。

如果阵列中的完整数据 vector 少于32个，则仍可以使用此格式。低位将不使用。或者您可以为每个 vector 设置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 vector )，然后使用 r0 = _mm512_maskz_sub_ps(p1,p2, local_selected[i]); 对减法进行零掩码。

如果实际上是限制加载端口uop吞吐量的瓶颈(每个时钟2个负载)，则可以尝试一次加载64位掩码数据，并使用mask-shift获取不同的低16位数据。但是，除非您的数据在L1D缓存中很热，否则这可能不会成为问题。

首先，使用“与掩码比较”即可轻松生成掩码数据，而无需交织。

理想情况下，您可以缓存阻止调用此代码的代码，以便可以在缓存中热数据时重用它们。例如从p1和p2的前64kiB中获得所需的所有组合，然后移至后面的元素，并在高速缓存中进行处理。